! source file: /raid21/jmuglia/iron/clean/updates/setmom.F subroutine setmom (is, ie, js, je) !======================================================================= ! set up everything which must be done only once per run !======================================================================= implicit none character (120) :: fname, new_file_name integer ie, is, je, js, i, ioun, j, jrow, n, m, indp, ip, k, kr integer jq, mask, kz, nv, ke, ks, iou, isle integer ib(10), ic(10) logical error, vmixset, hmixset, exists real snapint, snapls, snaple, snapde, trajint, cksum, checksum real ocnp, boxat, boxau, dvolt, dvolu, sum, pctocn, diffa, amix real runstep, dtatms, ahbkg, spnep, senep, c1e3 include "size.h" include "param.h" include "pconst.h" include "stdunits.h" include "accel.h" include "calendar.h" logical annlev, annlevobc, initpt include "coord.h" include "cpolar.h" include "cprnts.h" include "cregin.h" include "csbc.h" include "ctmb.h" include "diag.h" include "docnam.h" include "emode.h" include "grdvar.h" include "hmixc.h" include "index.h" include "iounit.h" include "isleperim.h" include "isopyc.h" include "levind.h" include "mw.h" include "scalar.h" include "stab.h" include "state.h" include "switch.h" include "tmngr.h" include "vmixc.h" integer kmz(imt,jmt) include "tidal_kv.h" real tmpij(imtm2,jmtm2) gravrho0r = grav*rho0r zetar = 1./500.e2 ! inverse vertical decay scale, zeta = 500m ogamma = 0.2*rho0r*zetar ! Osborn, 1980's mixing efficiency, gamma ! scaled by 1./(zeta*rho0) ! error in tracer conservation generated by wide_open_mw = .true. ! if (jmw .eq. jmt) then ! wide_open_mw = .true. ! else wide_open_mw = .false. ! endif ! stability diagnostic call stabi error = .false. vmixset = .false. hmixset = .false. visc_cbu_limit = 1.0e6 diff_cbt_limit = 1.0e6 ! user specified tracer names are placed into "trname" here. do m=1,nt trname(m) = '**unknown***' enddo trname(1) = 'potentl_temp' trname(2) = 'salinity ' !----------------------------------------------------------------------- ! get i/o units needed for prognostic variables !----------------------------------------------------------------------- call getunitnumber (kflds) call getunitnumber (latdisk(1)) call getunitnumber (latdisk(2)) !----------------------------------------------------------------------- ! read in 2d array of tidal energy dissipation rate ! similar to that of Jayne and St. Laurent, GRL 2001 ! except here we use E0 = W/(A*Nlev) ! no need to convert MOM native cgs ! units lesson: W = kg m^2 /s^3, therefore, W/m^2 = kg/s^3 ! therefore multiply edr by N0*1e3 to get g/s^3, where N0 is some ! reference value of N. ! in tidal.nc the energy flux already includes the Nb, so it's only ! multiplied by 1e3 to convert W/m^2 in g/s^3 !----------------------------------------------------------------------- edrm2(:,:,:)=0. edrs2(:,:,:)=0. edrk1(:,:,:)=0. edro1(:,:,:)=0. c edr(:,:) = 0. ib(:) = 1 ic(:) = 1 ic(1) = imt ic(2) = jmt ic(3) = 1 fname = new_file_name ("O_tidenrg_egbert.nc") inquire (file=trim(fname), exist=exists) if (exists) then c1e3 = 1.e3 call openfile (trim(fname), iou) do k=1,km ib(3) = k c print*, "read tidal_egbert, k=",k call getvara ('M2', iou, imt*jmt, ib, ic & ,edrm2(1:imt,k,1:jmt), c1e3, c0) call getvara ('S2', iou, imt*jmt, ib, ic & ,edrs2(1:imt,k,1:jmt), c1e3, c0) call getvara ('K1', iou, imt*jmt, ib, ic & ,edrk1(1:imt,k,1:jmt), c1e3, c0) call getvara ('O1', iou, imt*jmt, ib, ic & ,edro1(1:imt,k,1:jmt), c1e3, c0) c print*, edrm2(1:imt,k,1:jmt) enddo call closefile (iou) c call getvara ('tidenrg', iou, imtm2*jmtm2, ib, ic, tmpij c &, c1e3, c0) c edr(2:imtm1,2:jmtm1) = tmpij(1:imtm2,1:jmtm2) else print*, "Warning => Can not find ",trim(fname) endif !----------------------------------------------------------------------- ! compute density coefficients based on depth of grid points !----------------------------------------------------------------------- call eqstate (zt, km, ro0, to, so, c, tmink, tmaxk, smink, smaxk) cksum = checksum(c, km, 9) write (stdout &,'(6x,"Checksum for density coefficients =",e14.7)') cksum !----------------------------------------------------------------------- ! if the MW is not fully opened, then set time level indicators in ! the MW ("tau-1" "tau" "tau+1") to constant values. !----------------------------------------------------------------------- if (.not. wide_open_mw) then taum1 = -1 tau = 0 taup1 = +1 endif !----------------------------------------------------------------------- ! set prognostic quantities to either initial conditions or restart !----------------------------------------------------------------------- ! initialize two dimensional fields on disk ! initialize stream function fields in memory do n=1,2 do jrow=1,jmt do i=1,imt psi(i,jrow,n) = c0 ptd(i,jrow) = c0 enddo enddo enddo ! initialize stream function guess fields on disk ! block`s 1 & 2 are for the stream function guess field on disk do n=1,nkflds call oput (kflds, nwds, n, ptd(1,1)) enddo ! initialize all latitude rows call rowi ! initialize time step counter = 0 itt = 0 irstdy = 0 msrsdy = 0 ! initialize all prognostic quantities from the restart if (.not. init) then fname = new_file_name ("restart_mom.nc") inquire (file=trim(fname), exist=exists) if (exists) call mom_rest_in (fname, is, ie, js, je) endif ! construct depth arrays associated with "u" cells call depth_u (kmt, imt, jmt, zw, km, kmu, h, hr) ! compute a topography checksum cksum = 0.0 ocnp = 0 do jrow=2,jmt-1 do i=2,imt-1 cksum = cksum + i*jrow*kmt(i,jrow) ocnp = ocnp + float(kmt(i,jrow)) enddo enddo write (stdout,*) ' "kmt" checksum = ', cksum !----------------------------------------------------------------------- ! zero integrated time average accumulators !----------------------------------------------------------------------- call ta_mom_tsi (0) !----------------------------------------------------------------------- ! initialize the time manager with specified initial conditions ! time, user reference time, model time, and how long to integrate. !----------------------------------------------------------------------- call tmngri (year0, month0, day0, hour0, min0, sec0 &, ryear, rmonth, rday, rhour, rmin, rsec &, irstdy, msrsdy, runlen, rununits, rundays, dtts) !----------------------------------------------------------------------- ! convert starting and ending longitudes for the stability tests ! to nearest model grid points. !----------------------------------------------------------------------- if (stabint .ge. c0) then iscfl = max(indp (cflons, xt, imt), 2) cflons = xt(iscfl) iecfl = min(indp (cflone, xt, imt), imt-1) cflone = xt(iecfl) jscfl = max(indp (cflats, yt, jmt), 2) cflats = yt(jscfl) jecfl = min(indp (cflate, yt, jmt), jmt-1) cflate = yt(jecfl) kscfl = indp (cfldps, zt, km) cfldps = zt(kscfl) kecfl = indp (cfldpe, zt, km) cfldpe = zt(kecfl) endif !----------------------------------------------------------------------- ! compute the surface area and volume of the ocean regions. index 0 ! refers to the sum of all regions. ! (values used in subroutine region are done in terms of meters, ! rather than centimetres) !----------------------------------------------------------------------- areag = c0 volgt = c0 do k=1,km volgk(k) = c0 enddo do n=0,numreg areat(n) = c0 areau(n) = c0 volt(n) = c0 volu(n) = c0 enddo do mask=1,nhreg areab(mask) = c0 volbt(mask) = c0 do k=1,km volbk(mask,k) = c0 enddo enddo do jrow=2,jmtm1 do i=2,imtm1 mask = mskhr(i,jrow) kz = kmt(i,jrow) if (kz .gt. 0 .and. mask .gt. 0) then boxat = cst(jrow) * dxt(i) * dyt(jrow) if (kmu(i,jrow) .ne. 0) then boxau = csu(jrow) * dxu(i) * dyu(jrow) else boxau = c0 endif areat(0) = areat(0) + boxat areau(0) = areau(0) + boxau areab(mask) = areab(mask) + boxat * 1.e-4 do k=1,kz volbk(mask,k) = volbk(mask,k) + boxat * dzt(k) * 1.e-6 dvolt = boxat * dzt(k) if (kmu(i,jrow) .ge. k) then dvolu = boxau * dzt(k) else dvolu = c0 endif n = nhreg*(mskvr(k)-1) + mask if (n .gt. 0) then volt(0) = volt(0) + dvolt volu(0) = volu(0) + dvolu volt(n) = volt(n) + dvolt volu(n) = volu(n) + dvolu do nv=1,nvreg ks = llvreg(nv,1) if (k .eq. ks) then areat(n) = areat(n) + boxat areau(n) = areau(n) + boxau endif enddo endif enddo endif enddo enddo do mask=0,numreg if (areat(mask) .eq. c0) then rareat(mask) = c0 else rareat(mask) = c1/areat(mask) endif if (areau(mask) .eq. c0) then rareau(mask) = c0 else rareau(mask) = c1/areau(mask) endif if (volt(mask) .eq. c0) then rvolt(mask) = c0 else rvolt(mask) = c1/volt(mask) endif if (volu(mask) .eq. c0) then rvolu(mask) = c0 else rvolu(mask) = c1/volu(mask) endif enddo do mask=1,nhreg do k=1,km volbt(mask) = volbt(mask) + volbk(mask,k) volgk(k) = volgk(k) + volbk(mask,k) enddo areag = areag + areab(mask) volgt = volgt + volbt(mask) enddo if (iotavg .ne. stdout .or. iotavg .lt. 0) then call getunit (iou, 'tracer_avg.dta' &, 'unformatted sequential append ieee') call reg1st (iou, .false., .true., .true., .false., .false.) call relunit (iou) endif if (iotavg .eq. stdout .or. iotavg .lt. 0) then call reg1st (stdout, .true., .true., .true., .false., .false.) endif ! compute and print statistics for regions sum = c0 do n=1,numreg sum = sum + volt(n) enddo sum = 100.0*sum/tcellv pctocn = 100.0*ocnp/float((imt-2)*(jmt-2)*km) diffa = 100.0 * (c1 - (tcella(1) - 10000.0*areag)/tcella(1)) write (stdout,9342) diffa, numreg, sum, pctocn 9342 format (' the horizontal regional masks cover',f8.3 &, '% of the total ocean surface area.'/ &, ' there are ', i6, ' regions over which tracer & ' &, 'momentum balances will be computed',/,' accounting for ' &, f6.2,'% of the total ocean volume.'/ &, 1x,f6.2,'% of the grid points lie within the ocean.'/) !----------------------------------------------------------------------- ! find all land mass perimeters for Poisson solvers !----------------------------------------------------------------------- auto_kmt_changes = .false. ! set mask for land mass perimeters on which to perform calculations ! imask(-n) = .false. [no equations ever on dry land mass n] ! imask(0) = .true. [equations at all mid ocean points] ! imask(n) = .true./.false [controls whether there will be ! equations on the ocean perimeter of ! land mass n] ! note: land mass 1 is the northwest-most land mass ! usually includes the "north pole", and at low resolutions, ! the "main continent" do isle=-mnisle,mnisle if (isle .ge. 0 .and. isle .le. nisle) then imask(isle) = .true. else imask(isle) = .false. endif enddo ! user-specified changes to island mask ! imask(1) = .true. ! imask(2) = .true. ! there are problems if imask is set .true. for a nonexistent ! island. ! print diagnostic information do isle=-mnisle,mnisle if (imask(isle)) then if (isle .eq. 0) then print '(a)','=> calculations enabled for mid ocean points' else print '(2a,i3)','=> calculations enabled for ocean ', & 'perimeter of land mass',isle endif endif enddo do isle=0,nisle if (.not. imask(isle)) then print '(2a,i3)','=> calculations disabled for ocean ', & 'perimeter of land mass',isle endif enddo ! imain is the land mass on which dpsi is normalized to 0 ! if imain is 0, then dpsi is not normalized. ! default value of imain is land mass with longest perimeter imain = min(2,nisle) do isle=1,nisle if (nippts(isle) .gt. nippts(imain)) then imain = isle endif enddo ! if any island perimeter is not calculated, imain must be one such do isle=1,nisle if (.not.imask(isle)) then imain = isle endif enddo if (imain .gt. 0 .and. imain .le. nisle) then print '(a,i4)', 'dpsi normalized to zero on land mass',imain elseif (imain .eq. 0) then print *, 'no normalization on dpsi' else print *, 'ERROR: illegal value for choice of normalization ', & 'land mass, imain =', imain endif print *,' (user may set "imain" to any valid land mass number)' !--------------------------------------------------------------------- ! compute checksum of density coefficients !--------------------------------------------------------------------- print *,' ' call print_checksum (c(1,1), km, 9 &, ' density coefficient checksum = ') !----------------------------------------------------------------------- ! checksum the starting stream function. !----------------------------------------------------------------------- call print_checksum (psi(1,1,1), imt, jmt &, ' checksum for psi(,,1) = ') call print_checksum (psi(1,1,2), imt, jmt &, ' checksum for psi(,,2) = ') !----------------------------------------------------------------------- ! compute an array to indicate "interior" stream function grid cells !----------------------------------------------------------------------- do jrow=1,jmt kmz(1,jrow) = 0 kmz(imt,jrow) = 0 enddo do i=1,imt kmz(i,1) = 0 kmz(i,jmt) = 0 enddo do jrow=2,jmtm1 do i=2,imt kmz(i,jrow) = min(kmu(i-1,jrow-1), kmu(i,jrow-1) &, kmu(i-1,jrow), kmu(i,jrow)) enddo enddo do jrow=1,jmt kmz(1,jrow) = kmz(imtm1,jrow) enddo !--------------------------------------------------------------------- ! find and print start & end indices for filtering !--------------------------------------------------------------------- write (stdout,9551) if (lsegf.gt.11) write (stdout,9552) write (stdout,9553) call findex (kmt, jmtfil, km, jft1, jft2, imt, istf, ietf) write (stdout,9554) call findex (kmu, jmtfil, km, jfu1, jfu2, imt, isuf, ieuf) write (stdout,9555) call findex (kmz, jmtfil, 1, jft1, jft2, imt, iszf, iezf) 9551 format (/' ==== start and end indices for Fourier filtering ====') 9552 format (' only 11 sets of indices fit across the page.', & ' others will not be printed.'/) 9553 format (/,' == filtering indices for t,s ==') 9554 format (/,' == filtering indices for u,v ==') 9555 format (/,' == filtering indices for stream function ==') !--------------------------------------------------------------------- ! print the timestep multipliers !--------------------------------------------------------------------- write (stdout,9601) (dtxcel(k),k=1,km) 9601 format(/,' "dtxcel(km)" tracer timestep multipliers:',/,10(f9.3)) !----------------------------------------------------------------------- ! initialize various things !----------------------------------------------------------------------- do jrow=1,jmt do i=1,imt ztd(i,jrow) = c0 zu(i,jrow,1) = c0 zu(i,jrow,2) = c0 enddo enddo ! Coriolis factors do jrow=1,jmt do i=1,imt cori(i,jrow,1) = c2*omega*sin(ulat(i,jrow)/radian) cori(i,jrow,2) = -cori(i,jrow,1) enddo enddo ! metric diffusion factors amix = am do jrow=1,jmt am3(jrow) = amix*(c1-tng(jrow)*tng(jrow))/(radius**2) am4(jrow,1) = -amix*c2*sine(jrow)/(radius*csu(jrow) & *csu(jrow)) am4(jrow,2) = -am4(jrow,1) enddo ! metric advection factors do jrow=1,jmt advmet(jrow,1) = tng(jrow)/radius advmet(jrow,2) = -advmet(jrow,1) enddo ! diffusive flux through bottom of cells do j=jsmw,jemw do k=0,km do i=1,imt diff_fb(i,k,j) = c0 adv_fb(i,k,j) = c0 enddo enddo enddo !----------------------------------------------------------------------- ! initialize diagnostics !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! initialize time mean "averaging" grid data !----------------------------------------------------------------------- call avgi !----------------------------------------------------------------------- ! initialize the isopycnal mixing !----------------------------------------------------------------------- call isopi (error, am, ah) !----------------------------------------------------------------------- ! initialize the model of ocean biogeochemistry and isotopes (MOBI) !----------------------------------------------------------------------- if (kpzd .ne. km) then print*, 'stop: kpzd is not equal to km' print*, ' currently inconsistent with use of mobi' endif do k=1,km if (dtxcel(k) .ne. 1.0) then print*, 'stop: dtxcel is not equal to one' print*, ' inconsistent with use of mobi' endif enddo call mobii !----------------------------------------------------------------------- ! do all consistency checks last !----------------------------------------------------------------------- call checks (error, vmixset, hmixset) return end subroutine depth_u (kmt, imt, jmt, zw, km, kmu, h, hr) !======================================================================= ! calculate depth arrays associated with "u" cells. ! input: ! kmt = number of oecan "t" cells from surface to bottom of ocean ! imt = longitudinal dimension of arrays ! jmt = latitudinal dimension of arrays ! zw = depth to bottom of "t" cells ! km = max number of depths ! output: ! kmu = number of ocean "u" cells from surface to bottom of ocean ! h = depth to ocean floor over "u" cells (cm) ! hr = reciprocal of "h" !======================================================================= implicit none integer i, imt, jmt, jrow, km integer kmt(imt,jmt), kmu(imt,jmt) real c0, c1 real h(imt,jmt), hr(imt,jmt), zw(km) !----------------------------------------------------------------------- ! set some constants !----------------------------------------------------------------------- c0 = 0.0 c1 = 1.0 !----------------------------------------------------------------------- ! compute number of vertical levels on the "u" grid !----------------------------------------------------------------------- do jrow=1,jmt kmu(imt,jrow) = 0 enddo do i=1,imt kmu(i,jmt) = 0 enddo do jrow=1,jmt-1 do i=1,imt-1 kmu(i,jrow) = min (kmt(i,jrow), kmt(i+1,jrow), kmt(i,jrow+1) &, kmt(i+1,jrow+1)) enddo enddo do jrow=1,jmt kmu(imt,jrow) = kmu(2,jrow) enddo !--------------------------------------------------------------------- ! compute depths and reciprocal depths over "u" cells !--------------------------------------------------------------------- do jrow=1,jmt do i=1,imt hr(i,jrow) = c0 h(i,jrow) = c0 if (kmu(i,jrow) .ne. 0) then hr(i,jrow) = c1/zw(kmu(i,jrow)) h (i,jrow) = zw(kmu(i,jrow)) endif enddo enddo return end subroutine rowi !----------------------------------------------------------------------- ! initialize prognostic quantities at "tau-1" and "tau" ! inputs: ! kmt = number of vertical levels on "t" cells ! yt = latitudes of "t" points ! zt = depths of "t" points !----------------------------------------------------------------------- implicit none character(120) :: fname, new_file_name, text integer i, itt, iou, ib(10), ic(10), j, jrow, n, k, kz, ln logical exists real dens, drodt, drods, p1, c100, c1e3, p001, p035, ucksum real vcksum, tcksum, scksum, theta0, tq, sq, s, d, tins, tinsit real checksum, C2K, jdi real geoheatflux include "size.h" include "param.h" include "npzd.h" ! rstd* include "pconst.h" include "stdunits.h" include "coord.h" include "csbc.h" include "iounit.h" include "levind.h" include "mw.h" include "diag.h" include "state.h" include "dens.h" real tmpik(imtm2,km), dmsk(imt,jmt) !----------------------------------------------------------------------- ! update pointers to tau-1, tau, & tau+1 data on disk. ! for latitude rows they point to latdisk(1) or latdisk(2) ! for 2D fields they point to records on kflds !----------------------------------------------------------------------- itt = 0 taum1disk = mod(itt+1,2) + 1 taudisk = mod(itt ,2) + 1 taup1disk = taum1disk p1 = 0.1 c100 = 100. c1e3 = 1000. p001 = 0.001 p035 = 0.035 C2K = 273.15 !----------------------------------------------------------------------- ! update pointers to tau-1, tau, & tau+1 data in the MW !----------------------------------------------------------------------- if (wide_open_mw) then ! rotate time levels instead of moving data taum1 = mod(itt+0,3) - 1 tau = mod(itt+1,3) - 1 taup1 = mod(itt+2,3) - 1 endif ucksum = 0.0 vcksum = 0.0 tcksum = 0.0 scksum = 0.0 !----------------------------------------------------------------------- ! initialize every latitude jrow either in the MW (when wide opened) ! or on disk (when jmw < jmt) !----------------------------------------------------------------------- do jrow=1,jmt if (wide_open_mw) then j = jrow else j = jmw endif jdi = min(max(jrow-1,1),jmtm2) !----------------------------------------------------------------------- ! zero out all variables. velocities are internal modes only !----------------------------------------------------------------------- do k=1,km do i=1,imt u(i,k,j,1,taup1) = c0 u(i,k,j,2,taup1) = c0 enddo enddo do n=1,nt do k=1,km do i=1,imt t(i,k,j,n,taup1) = c0 enddo enddo enddo !----------------------------------------------------------------------- ! set tracers !----------------------------------------------------------------------- do n=1,nt do i=1,imt do k=1,kmt(i,jrow) if (trim(mapt(n)) .eq. 'temp') then t(i,k,j,n,taup1) = theta0 (yt(jrow), zt(k)) elseif (trim(mapt(n)) .eq. 'salt') then t(i,k,j,n,taup1) = 0.03472 - 0.035 elseif (trim(mapt(n)) .eq. 'dic') then t(i,k,j,n,taup1) = 2.315 elseif (trim(mapt(n)) .eq. 'alk') then t(i,k,j,n,taup1) = 2.429 elseif (trim(mapt(n)) .eq. 'o2') then t(i,k,j,n,taup1) = 0.1692 elseif (trim(mapt(n)) .eq. 'po4') then if (k .eq. 1) then t(i,k,j,n,taup1) = 0.543 else t(i,k,j,n,taup1) = 2.165 endif elseif (trim(mapt(n)) .eq. 'dop') then if (k .le. 2) then t(i,k,j,n,taup1) = 0.156 elseif (k .gt. 2 .and. k .le. 12) then t(i,k,j,n,taup1) = 0.039 else t(i,k,j,n,taup1) = 0.0078 endif elseif (trim(mapt(n)) .eq. 'phyt') then t(i,k,j,n,taup1) = 0.14*exp((zt(1)-zt(k))/100.e2) elseif (trim(mapt(n)) .eq. 'zoop') then t(i,k,j,n,taup1) = 0.014*exp((zt(1)-zt(k))/100.e2) elseif (trim(mapt(n)) .eq. 'detr') then t(i,k,j,n,taup1) = 1.e-4 elseif (trim(mapt(n)) .eq. 'detrfe') then t(i,k,j,n,taup1) = 1.e-4*14.e-6*6.625 elseif (trim(mapt(n)) .eq. 'no3') then if (k .eq. 1) then t(i,k,j,n,taup1) = 5.30 else t(i,k,j,n,taup1) = 30.84 endif elseif (trim(mapt(n)) .eq. 'dfe') then if (k .eq. 1) then t(i,k,j,n,taup1) = 0.6e-3 else t(i,k,j,n,taup1) = 0.6e-3 endif elseif (trim(mapt(n)) .eq. 'don') then if (k .le. 2) then t(i,k,j,n,taup1) = 3.5 elseif (k .gt. 2 .and. k .le. 12) then t(i,k,j,n,taup1) = 1.5 else t(i,k,j,n,taup1) = 0.5 endif elseif (trim(mapt(n)) .eq. 'diaz') then t(i,k,j,n,taup1) = 0.014*exp((zt(1)-zt(k))/100.e2) elseif (trim(mapt(n)) .eq. 'c14') then t(i,k,j,n,taup1) = -150 ! permil (converted later) elseif (trim(mapt(n)) .eq. 'cfc11') then t(i,k,j,n,taup1) = 0.0 elseif (trim(mapt(n)) .eq. 'cfc12') then t(i,k,j,n,taup1) = 0.0 else if (k .eq. 1) then t(i,k,j,n,taup1) = c1 else t(i,k,j,n,taup1) = c0 endif endif enddo enddo ib(:) = 1 ib(2) = jdi ib(3) = 1 ic(:) = imtm2 ic(2) = 1 ic(3) = km if (trim(mapt(n)) .eq. 'temp') then ! expecting units of degrees K fname = new_file_name ("O_temp.nc") inquire (file=trim(fname), exist=exists) if (exists) then call openfile (trim(fname), iou) tmpik(1:imtm2,1:km) = t(2:imtm1,1:km,j,n,taup1) call getvara ('O_temp', iou, imtm2*km, ib, ic, tmpik &, c1, c0) text = "C" call getatttext (iou, 'O_temp', 'units', text) ! convert to model units (C) if (trim(text).eq."K") & where ((tmpik(1:imtm2,1:km)) .lt. 1.e30) & tmpik(1:imtm2,1:km) = tmpik(1:imtm2,1:km) - C2K t(2:imtm1,1:km,j,n,taup1) = tmpik(1:imtm2,1:km) call setbcx (t(1:imt,1:km,j,n,taup1),imt,km) else if (jrow .eq. 1) & print*, "Warning => Can not find ",trim(fname) endif elseif (trim(mapt(n)) .eq. 'salt') then ! expecting units of psu fname = new_file_name ("O_sal.nc") inquire (file=trim(fname), exist=exists) if (exists) then call openfile (trim(fname), iou) tmpik(1:imtm2,1:km) = t(2:imtm1,1:km,j,n,taup1) ! convert to model units ((psu-35)/1000) call getvara ('O_sal', iou, imtm2*km, ib, ic, tmpik &, p001, -p035) t(2:imtm1,1:km,j,n,taup1) = tmpik(1:imtm2,1:km) call setbcx (t(1:imt,1:km,j,n,taup1),imt,km) else if (jrow .eq. 1) & print*, "Warning => Can not find ",trim(fname) endif elseif (trim(mapt(n)) .eq. 'dic') then ! expecting units of mol m-3 (equals model units, umol cm-3) fname = new_file_name ("O_totcarb.nc") inquire (file=trim(fname), exist=exists) if (exists) then call openfile (trim(fname), iou) tmpik(1:imtm2,1:km) = t(2:imtm1,1:km,j,n,taup1) call getvara ('O_totcarb', iou, imtm2*km, ib, ic, tmpik &, c1, c0) t(2:imtm1,1:km,j,n,taup1) = tmpik(1:imtm2,1:km) call setbcx (t(1:imt,1:km,j,n,taup1),imt,km) else if (jrow .eq. 1) & print*, "Warning => Can not find ",trim(fname) endif elseif (trim(mapt(n)) .eq. 'alk') then ! expecting units of mol m-3 (equals model units, umol cm-3) fname = new_file_name ("O_alk.nc") inquire (file=trim(fname), exist=exists) if (exists) then call openfile (trim(fname), iou) tmpik(1:imtm2,1:km) = t(2:imtm1,1:km,j,n,taup1) call getvara ('O_alk', iou, imtm2*km, ib, ic, tmpik &, c1, c0) t(2:imtm1,1:km,j,n,taup1) = tmpik(1:imtm2,1:km) call setbcx (t(1:imt,1:km,j,n,taup1),imt,km) else if (jrow .eq. 1) & print*, "Warning => Can not find ",trim(fname) endif elseif (trim(mapt(n)) .eq. 'o2') then ! expecting units of mol m-3 (equals model units, umol cm-3) fname = new_file_name ("O_o2.nc") inquire (file=trim(fname), exist=exists) if (exists) then call openfile (trim(fname), iou) tmpik(1:imtm2,1:km) = t(2:imtm1,1:km,j,n,taup1) call getvara ('O_o2', iou, imtm2*km, ib, ic, tmpik &, c1, c0) t(2:imtm1,1:km,j,n,taup1) = tmpik(1:imtm2,1:km) call setbcx (t(1:imt,1:km,j,n,taup1),imt,km) else if (jrow .eq. 1) & print*, "Warning => Can not find ",trim(fname) endif elseif (trim(mapt(n)) .eq. 'po4') then ! expecting units of mol m-3 fname = new_file_name ("O_po4.nc") inquire (file=trim(fname), exist=exists) if (exists) then call openfile (trim(fname), iou) tmpik(1:imtm2,1:km) = t(2:imtm1,1:km,j,n,taup1) ! convert to model units (nmol cm-3) call getvara ('O_po4', iou, imtm2*km, ib, ic, tmpik &, c1e3, c0) t(2:imtm1,1:km,j,n,taup1) = tmpik(1:imtm2,1:km) call setbcx (t(1:imt,1:km,j,n,taup1),imt,km) else if (jrow .eq. 1) & print*, "Warning => Can not find ",trim(fname) endif elseif (trim(mapt(n)) .eq. 'dfe') then ! expecting units of mol m-3 fname = new_file_name ("O_dfe.nc") inquire (file=trim(fname), exist=exists) if (exists) then call openfile (trim(fname), iou) tmpik(1:imtm2,1:km) = t(2:imtm1,1:km,j,n,taup1) ! convert to model units (nmol cm-3) call getvara ('O_dfe', iou, imtm2*km, ib, ic, tmpik &, c1e3, c0) t(2:imtm1,1:km,j,n,taup1) = tmpik(1:imtm2,1:km) call setbcx (t(1:imt,1:km,j,n,taup1),imt,km) else if (jrow .eq. 1) & print*, "Warning => Can not find ",trim(fname) endif elseif (trim(mapt(n)) .eq. 'no3') then ! expecting units of mol m-3 fname = new_file_name ("O_no3.nc") inquire (file=trim(fname), exist=exists) if (exists) then call openfile (trim(fname), iou) tmpik(1:imtm2,1:km) = t(2:imtm1,1:km,j,n,taup1) ! convert to model units (nmol cm-3) call getvara ('O_no3', iou, imtm2*km, ib, ic, tmpik &, c1e3, c0) t(2:imtm1,1:km,j,n,taup1) = tmpik(1:imtm2,1:km) call setbcx (t(1:imt,1:km,j,n,taup1),imt,km) else if (jrow .eq. 1) & print*, "Warning => Can not find ",trim(fname) endif elseif (trim(mapt(n)) .eq. 'c14') then ! expecting units of permil fname = new_file_name ("O_dc14.nc") inquire (file=trim(fname), exist=exists) if (exists) then call openfile (trim(fname), iou) tmpik(1:imtm2,1:km) = t(2:imtm1,1:km,j,n,taup1) ! convert later call getvara ('O_dc14', iou, imtm2*km, ib, ic, tmpik &, c1, c0) t(2:imtm1,1:km,j,n,taup1) = tmpik(1:imtm2,1:km) call setbcx (t(1:imt,1:km,j,n,taup1),imt,km) else if (jrow .eq. 1) & print*, "Warning => Can not find ",trim(fname) endif endif do k=1,km do i=1,imt if (kmt(i,jrow) .lt. k) t(i,k,j,n,taup1) = c0 enddo enddo enddo do k=1,km do i=1,imt t(i,k,j,idin15,taup1) = 1.005*rn15std*t(i,k,j,ino3,taup1) & / (1 + 1.005*rn15std) t(i,k,j,idon15,taup1) = 1.0045*rn15std*t(i,k,j,idon,taup1) & / (1 + 1.0045*rn15std) t(i,k,j,idon15,taup1) = rn15std*t(i,k,j,idon,taup1) & / (1 + rn15std) t(i,k,j,iphytn15,taup1) = rn15std*t(i,k,j,iphyt,taup1) & / (1 + rn15std) t(i,k,j,izoopn15,taup1) = rn15std*t(i,k,j,izoop,taup1) & / (1 + rn15std) t(i,k,j,idetrn15,taup1) = rn15std*t(i,k,j,idetr,taup1) & / (1 + rn15std) t(i,k,j,idiazn15,taup1) = rn15std*t(i,k,j,idiaz,taup1) & / (1 + rn15std) enddo enddo do k=1,km do i=1,imt ! initialize at delta del13C = 0 permil t(i,k,j,idic13,taup1) = t(i,k,j,idic,taup1)*1.0004*rc13std & / (1. + 1.0004*rc13std) t(i,k,j,idoc13,taup1) = t(i,k,j,idon,taup1)*6.625e-3 & *rc13std/(1. + rc13std) t(i,k,j,iphytc13,taup1) = t(i,k,j,iphyt,taup1)*6.625e-3 & *rc13std/(1. + rc13std) t(i,k,j,izoopc13,taup1) = t(i,k,j,izoop,taup1)*6.625e-3 & *rc13std/(1. + rc13std) t(i,k,j,idetrc13,taup1) = t(i,k,j,idetr,taup1)*6.625e-3 & *rc13std/(1. + rc13std) t(i,k,j,idiazc13,taup1) = t(i,k,j,idiaz,taup1)*6.625e-3 & *rc13std/(1. + rc13std) enddo enddo ! convert c14 from permil to model units (umol cm-3) do k=1,km do i=1,imt t(i,k,j,ic14,taup1) = (t(i,k,j,ic14,taup1)*0.001 + 1) & *t(i,k,j,idic,taup1)*rc14std enddo enddo !----------------------------------------------------------------------- ! set sea level reference density field !----------------------------------------------------------------------- ! set reference density field to initial field do i=1,imt do k=1,km t_slh(i,jrow,k,1) = t(i,k,j,itemp,taup1) t_slh(i,jrow,k,2) = t(i,k,j,isalt,taup1) enddo enddo ! read reference density field if it exists fname = new_file_name ("O_slhref.nc") inquire (file=trim(fname), exist=exists) if (exists) then call openfile (fname, iou) ib(:) = 1 ib(2) = jdi ic(:) = 1 ic(1) = imtm2 ic(3) = km tmpik(1:imtm2,1:km) = t(2:imtm1,1:km,j,itemp,taup1) call getvara ('O_temp', iou, imtm2*km, ib, ic, tmpik &, c1, c0) text = "C" call getatttext (iou, 'O_temp', 'units', text) ! convert to model units (C) if (trim(text) .eq. "K") & where ((tmpik(1:imtm2,1:km)) .lt. 1.e30) & tmpik(1:imtm2,1:km) = tmpik(1:imtm2,1:km) - C2K do i=2,imtm1 do k=1,km t_slh(i,jrow,k,1) = tmpik(i-1,k) enddo enddo tmpik(1:imtm2,1:km) = t(2:imtm1,1:km,j,isalt,taup1) call getvara ('O_sal', iou, imtm2*km, ib, ic, tmpik &, p001, -p035) do i=2,imtm1 do k=1,km t_slh(i,jrow,k,2) = tmpik(i-1,k) enddo enddo endif !----------------------------------------------------------------------- ! calculate reference density field !----------------------------------------------------------------------- do i=2,imtm1 do k=1,kmt(i,jrow) s = (t_slh(i,jrow,k,2)*1000.) + 35. d = zt(k)*0.01 tins = tinsit(t_slh(i,jrow,k,1), s, d) d_slh(i,jrow,k) = dens(tins-to(k),t_slh(i,jrow,k,2)-so(k),k) enddo enddo ! zero for future accumulation t_slh(:,:,:,:) = 0. !----------------------------------------------------------------------- ! initialize surface boundary !----------------------------------------------------------------------- if (isst .ne. 0) sbc(1:imt,jrow,isst) = t(1:imt,1,j,itemp,taup1) if (isss .ne. 0) sbc(1:imt,jrow,isss) = t(1:imt,1,j,isalt,taup1) if (issdic .ne. 0) & sbc(1:imt,jrow,issdic) = t(1:imt,1,j,idic,taup1) if (issdic13 .ne. 0) & sbc(1:imt,jrow,issdic13) = t(1:imt,1,j,idic13,taup1) if (isso2 .ne. 0) sbc(1:imt,jrow,isso2) = t(1:imt,1,j,io2,taup1) if (issalk .ne. 0) & sbc(1:imt,jrow,issalk) = t(1:imt,1,j,ialk,taup1) if (isspo4 .ne. 0) & sbc(1:imt,jrow,isspo4) = t(1:imt,1,j,ipo4,taup1) if (issdop .ne. 0) & sbc(1:imt,jrow,issdop) = t(1:imt,1,j,idop,taup1) if (issphyt .ne. 0) & sbc(1:imt,jrow,issphyt) = t(1:imt,1,j,iphyt,taup1) if (isszoop .ne. 0) & sbc(1:imt,jrow,isszoop) = t(1:imt,1,j,izoop,taup1) if (issdetr .ne. 0) & sbc(1:imt,jrow,issdetr) = t(1:imt,1,j,idetr,taup1) if (issno3 .ne. 0) & sbc(1:imt,jrow,issno3) = t(1:imt,1,j,ino3,taup1) if (issdon .ne. 0) & sbc(1:imt,jrow,issdon) = t(1:imt,1,j,idon,taup1) if (issdiaz .ne. 0) & sbc(1:imt,jrow,issdiaz) = t(1:imt,1,j,idiaz,taup1) if (issdin15 .ne. 0) & sbc(1:imt,jrow,issdin15) = t(1:imt,1,j,idin15,taup1) if (issdon15 .ne. 0) & sbc(1:imt,jrow,issdon15) = t(1:imt,1,j,idon15,taup1) if (issphytn15 .ne. 0) & sbc(1:imt,jrow,issphytn15) = t(1:imt,1,j,iphytn15,taup1) if (isszoopn15 .ne. 0) & sbc(1:imt,jrow,isszoopn15) = t(1:imt,1,j,izoopn15,taup1) if (issdetrn15 .ne. 0) & sbc(1:imt,jrow,issdetrn15) = t(1:imt,1,j,idetrn15,taup1) if (issdiazn15 .ne. 0) & sbc(1:imt,jrow,issdiazn15) = t(1:imt,1,j,idiazn15,taup1) if (issdoc13 .ne. 0) & sbc(1:imt,jrow,issdoc13) = t(1:imt,1,j,idoc13,taup1) if (issphytc13 .ne. 0) & sbc(1:imt,jrow,issphytc13) = t(1:imt,1,j,iphytc13,taup1) if (isszoopc13 .ne. 0) & sbc(1:imt,jrow,isszoopc13) = t(1:imt,1,j,izoopc13,taup1) if (issdetrc13 .ne. 0) & sbc(1:imt,jrow,issdetrc13) = t(1:imt,1,j,idetrc13,taup1) if (issdiazc13 .ne. 0) & sbc(1:imt,jrow,issdiazc13) = t(1:imt,1,j,idiazc13,taup1) if (issc14 .ne. 0) & sbc(1:imt,jrow,issc14) = t(1:imt,1,j,ic14,taup1) if (isscfc11 .ne. 0) & sbc(1:imt,jrow,isscfc11) = t(1:imt,1,j,icfc11,taup1) if (isscfc12 .ne. 0) & sbc(1:imt,jrow,isscfc12) = t(1:imt,1,j,icfc12,taup1) if (issdfe .ne. 0) & sbc(1:imt,jrow,issdfe) = t(1:imt,1,j,idfe,taup1) if (issdetrfe .ne. 0) & sbc(1:imt,jrow,issdetrfe) = t(1:imt,1,j,idetrfe,taup1) ! bottom heat flux do i=1,imt bhf(i,jrow) = geoheatflux(xt(i),yt(jrow)) enddo !----------------------------------------------------------------------- ! zero out tracers in land points !----------------------------------------------------------------------- do i=1,imt kz = kmt(i,jrow) do k=1,km if (k .gt. kz) then do n=1,nt t(i,k,j,n,taup1) = c0 enddo endif enddo enddo !----------------------------------------------------------------------- ! checksum the initial conditions !----------------------------------------------------------------------- ucksum = ucksum + checksum (u(1,1,j,1,taup1), imt, km) vcksum = vcksum + checksum (u(1,1,j,2,taup1), imt, km) tcksum = tcksum + checksum (t(1,1,j,1,taup1), imt, km) scksum = scksum + checksum (t(1,1,j,2,taup1), imt, km) !----------------------------------------------------------------------- ! initialize every latitude jrow either on disk (when jmw < jmt) ! or in the MW (when the last jrow is complete and jmw = jmt) !----------------------------------------------------------------------- if (wide_open_mw) then if (jrow .eq. jmt) then call copy_all_rows (taup1, tau) call copy_all_rows (tau, taum1) endif else call putrow (latdisk(taudisk), nslab, jrow &, u(1,1,j,1,taup1), t(1,1,j,1,taup1)) call putrow (latdisk(taup1disk), nslab, jrow &, u(1,1,j,1,taup1), t(1,1,j,1,taup1)) endif enddo !----------------------------------------------------------------------- ! find inital average surface references !----------------------------------------------------------------------- print*, " " print*, "inital average surface references: " dmsk(:,:) = 1. where (kmt(:,:) .eq. 0) dmsk(:,:) = 0. gaost(:) = 0. if (isalt .ne. 0 .and. isss .ne. 0) then call areaavg (sbc(1,1,isss), dmsk, gaost(isalt)) gaost(isalt) = gaost(isalt) + 0.035 socn = gaost(isalt) print*, "global average sea surface salinity (psu) = " &, gaost(isalt)*1000. endif if (idic .ne. 0 .and. issdic .ne. 0) then call areaavg (sbc(1,1,issdic), dmsk, gaost(idic)) print*, "global average sea surface dic (mol m-3) = " &, gaost(idic) endif if (idic13 .ne. 0 .and. issdic13 .ne. 0) then call areaavg (sbc(1,1,issdic13), dmsk, gaost(idic13)) print*, "global average sea surface dic 13 (mol m-3) = " &, gaost(idic13) endif if (io2 .ne. 0 .and. isso2 .ne. 0) then call areaavg (sbc(1,1,isso2), dmsk, gaost(io2)) print*, "global average sea surface oxygen (mol m-3) = " &, gaost(io2) endif if (ialk .ne. 0 .and. issalk .ne. 0) then call areaavg (sbc(1,1,issalk), dmsk, gaost(ialk)) print*, "global average sea surface alkalinity (mol m-3) = " &, gaost(ialk) endif if (ipo4 .ne. 0 .and. isspo4 .ne. 0) then call areaavg (sbc(1,1,isspo4), dmsk, gaost(ipo4)) print*, "global average sea surface phosphate (mol m-3) = " &, gaost(ipo4)*0.001 endif if (idop .ne. 0 .and. issdop .ne. 0) then call areaavg (sbc(1,1,issdop), dmsk, gaost(idop)) print*, "global average sea surface DOP (mol m-3) = " &, gaost(idop)*0.001 endif if (ino3 .ne. 0 .and. issno3 .ne. 0) then call areaavg (sbc(1,1,issno3), dmsk, gaost(ino3)) print*, "global average sea surface nitrate (mol m-3) = " &, gaost(ino3)*0.001 endif if (idon .ne. 0 .and. issdon .ne. 0) then call areaavg (sbc(1,1,issdon), dmsk, gaost(idon)) print*, "global average sea surface DON (mol m-3) = " &, gaost(idon)*0.001 endif if (idin15 .ne. 0 .and. issdin15 .ne. 0) then call areaavg (sbc(1,1,issdin15), dmsk, gaost(idin15)) print*, "global average sea surface nitrate 15 (mol m-3) = " &, gaost(idin15)*0.001 endif if (idon15 .ne. 0 .and. issdon15 .ne. 0) then call areaavg (sbc(1,1,issdon15), dmsk, gaost(idon15)) print*, "global average sea surface DON15 (mol m-3) = " &, gaost(idon15)*0.001 endif if (idoc13 .ne. 0 .and. issdoc13 .ne. 0) then call areaavg (sbc(1,1,issdoc13), dmsk, gaost(idoc13)) print*, "global average sea surface DOC13" &, " (mol m-3) = ", gaost(idoc13)*0.001 endif if (iphytc13 .ne. 0 .and. issphytc13 .ne. 0) then call areaavg (sbc(1,1,issphytc13), dmsk, gaost(iphytc13)) print*, "global average sea surface phytoplankton C13" &, " (mol m-3) = ", gaost(iphytc13)*0.001 endif if (izoopc13 .ne. 0 .and. isszoopc13 .ne. 0) then call areaavg (sbc(1,1,isszoopc13), dmsk, gaost(izoopc13)) print*, "global average sea surface zooplankton C13" &, " (mol m-3) = ", gaost(izoopc13)*0.001 endif if (idetrc13 .ne. 0 .and. issdetrc13 .ne. 0) then call areaavg (sbc(1,1,issdetrc13), dmsk, gaost(idetrc13)) print*, "global average sea surface detritus c13" &, " (mol m-3) = ", gaost(idetrc13)*0.001 endif if (idiazc13 .ne. 0 .and. issdiazc13 .ne. 0) then call areaavg (sbc(1,1,issdiazc13), dmsk, gaost(idiazc13)) print*, "global average sea surface diazotrophs c13" &, " (mol m-3) = ", gaost(idiazc13)*0.001 endif if (idfe .ne. 0 .and. issdfe .ne. 0) then call areaavg (sbc(1,1,issdfe), dmsk, gaost(idfe)) print*, "global average sea surface iron (mol m-3) = " &, gaost(idfe)*0.001 endif if (ic14 .ne. 0 .and. issc14 .ne. 0) then call areaavg (sbc(1,1,issc14), dmsk, gaost(ic14)) print*, "global average sea surface carbon 14 (mol m-3) = " &, gaost(ic14) endif if (icfc11 .ne. 0 .and. isscfc11 .ne. 0) then call areaavg (sbc(1,1,isscfc11), dmsk, gaost(icfc11)) print*, "global average sea surface cfc 11 (mol m-3) = " &, gaost(icfc11) endif if (icfc12 .ne. 0 .and. isscfc12 .ne. 0) then call areaavg (sbc(1,1,isscfc12), dmsk, gaost(icfc12)) print*, "global average sea surface cfc 12 (mol m-3) = " &, gaost(icfc12) endif print*, " " write (stdout,*) ' I.C. checksum for t =',tcksum write (stdout,*) ' I.C. checksum for s =',scksum write (stdout,*) ' I.C. checksum for u =',ucksum write (stdout,*) ' I.C. checksum for v =',vcksum return end function theta0 (ydeg, depth) !======================================================================= ! this subroutine returns estimates of global mean potential ! temperature for model initialization as a function of depth. ! it is used to produce a reference thermal stratification for the ! upper 2000m of the MOM`s test case. below 2000m, the ! potential temperature returned is 2.0 degrees C. surface ! values are set slightly above 18.4 degrees C at the reference ! latitude "reflat". ! the estimates are produced from a 7th order polynomial fit to ! the annual mean world ocean potential temperature observations ! of Levitus (1982). ! input [units]: ! a latitude (ydeg): [degrees] ! a zt value (depth): [centimetres] ! output [units]: ! potential temperature estimate (est): [degrees centigrade] ! variables: ! coeft = coefficients for the polynomial fit of potential ! temperature vs. depth ! reflat = reference latitude at which observed surface ! temperatures approximately equal coeft(1) ! factor = the ratio of the cosine of the latitude requested ! ("ydeg") to the reference latitude ("reflat") ! used to scale the upper 2000 meters of the vertical ! temperature profile ! tmin,tmax = the minimum and maximum potential temperatures ! allowed at the time of model initialization ! reference: ! Levitus, S., Climatological atlas of the world ocean, NOAA ! Prof. Paper 13, US Gov`t printing Office, Washington, DC, 1982. !======================================================================= implicit none integer ndeg, nn parameter (ndeg=7) real c0, coslat, factor, pi, theta0, tmin, tmax, refcos, reflat real ydeg, z, depth, est, bb real coeft(ndeg+1) save coeft, factor, pi, tmin, tmax, reflat data coeft / 0.184231944E+02,-0.430306621E-01, 0.607121504E-04 & ,-0.523806281E-07, 0.272989082E-10,-0.833224666E-14 & , 0.136974583E-17,-0.935923382E-22/ data tmin, tmax, reflat /2.0, 25.0, 34.0/ c0 = 0.0 pi = atan(1.0) * 4.0 refcos = abs(cos(pi*reflat/180.)) coslat = abs(cos(pi*ydeg/180.)) factor = coslat/refcos z = depth * 0.01 if (z .gt. 2000.) then est = 2.0 else est = c0 bb = 1.0 do nn=1,ndeg+1 ! if (nn.gt.1) bb = z**(nn-1) est = est + coeft(nn)*bb bb = bb*z enddo est = est * factor endif if (est .gt. tmax) est = tmax if (est .lt. tmin) est = tmin theta0 = est return end