subroutine setmom (is, ie, js, je) #if defined O_mom !======================================================================= ! set up everything which must be done only once per run !======================================================================= implicit none character (120) :: fname, new_file_name, logfile 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, inqvardef 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 # if defined O_neptune include "cnep.h" # endif include "coord.h" # if defined O_fourfil || defined O_firfil include "cpolar.h" # endif include "cprnts.h" include "cregin.h" include "csbc.h" # if defined O_shortwave include "cshort.h" # endif 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" # if defined O_isopycmix include "isopyc.h" # endif include "levind.h" include "mw.h" include "scalar.h" include "stab.h" include "state.h" include "switch.h" include "tmngr.h" include "vmixc.h" # if defined O_mom_tbt include "tbt.h" # endif # if defined O_fourfil || defined O_firfil integer kmz(imt,jmt) # endif # if defined O_tidal_kv include "tidal_kv.h" # endif # if defined O_fwa include "fwa.h" # endif # if defined O_npzd include "npzd.h" # endif real tmpij(imtm2,jmtm2) # if defined O_tidal_kv 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) # endif ! 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 # if defined O_rigid_lid_surface_pressure alph = 1.0 gam = 0.0 theta = 1.0 # endif # if defined O_implicit_free_surface alph = 0.3333333 gam = 0.3333333 theta = 0.5 # endif ! 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)) # if defined O_tidal_kv !----------------------------------------------------------------------- ! 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 !----------------------------------------------------------------------- edr(:,:) = 0. ib(:) = 1 ic(:) = 1 ic(1) = imtm2 ic(2) = jmtm2 fname = new_file_name ("O_tidenrg.nc") inquire (file=trim(fname), exist=exists) if (exists) then c1e3 = 1.e3 call openfile (trim(fname), iou) call getvara ('O_tidenrg', iou, imtm2*jmtm2, ib, ic, tmpij &, c1e3, c0) edr(2:imtm1,2:jmtm1) = tmpij(1:imtm2,1:jmtm2) else print*, "Warning => Can not find ",trim(fname) endif # 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 defined O_linearized_density ! eliminate the pressure effect when using linearized density do k=1,km ro0(k) = c0 to(k) = c0 so(k) = c0 enddo # endif # if defined O_linearized_advection ! set the initial idealized stratification as a function of ! temperature only. Salinity is fixed at 35 ppt (tbar(k,2)=0) but ! passive tracers (if nt>2) are 1.0 at k=1 and 0.0 for k>1 tzero = 7.5 tone = 10.0 z0 = 30.0e2 bigl = 80.0e2 write (stdout,'(//a/)') & 'Initial Tracer profile as a function of depth' do k=1,km tbarz(k,1) = tzero*(1.0-tanh((zt(k)-bigl)/z0)) + & tone*(1.0-zt(k)/zt(km)) tbarz(k,2) = c0 if (nt .gt. 2) then do n=3,nt if (k .eq. 1) then tbarz(k,n) = c1 else tbarz(k,n) = c0 endif enddo endif write (stdout,'(1x,"k=",i3,", zt(k)=",f8.1 &, "cm, tbarz(k,n),n=1,nt =",10e14.7)') & k, zt(k),(tbarz(k,n),n=1,nt) enddo # endif !----------------------------------------------------------------------- ! 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 !----------------------------------------------------------------------- # if defined O_neptune ! calculate "neptune" velocities call neptune # endif ! initialize two dimensional fields on disk # if defined O_rigid_lid_surface_pressure || defined O_implicit_free_surface ! initialize surface pressure fields in memory do jrow=1,jmt do i=1,imt ps(i,jrow,1) = c0 ps(i,jrow,2) = c0 pguess(i,jrow) = c0 ubar(i,jrow,1) = c0 ubar(i,jrow,2) = c0 ubarm1(i,jrow,1) = c0 ubarm1(i,jrow,2) = c0 enddo enddo # endif # if defined O_stream_function ! initialize stream function fields in memory do n=1,2 do jrow=1,jmt do i=1,imt # if defined O_neptune ! initialize to "neptune" streamfunction psi(i,jrow,n) = pnep(i,jrow) ptd(i,jrow) = c0 # else psi(i,jrow,n) = c0 ptd(i,jrow) = c0 # endif 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 # endif ! 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) # if defined O_restart_2 fname = new_file_name ("restart_2_mom.nc") inquire (file=trim(fname), exist=exists) if (exists) call mom_rest_in (fname, is, ie, js, je) # endif 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 # if defined O_time_averages !----------------------------------------------------------------------- ! zero integrated time average accumulators !----------------------------------------------------------------------- call ta_mom_tsi (0) # endif !----------------------------------------------------------------------- ! 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) # if defined O_stability_tests !----------------------------------------------------------------------- ! 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 # endif # if defined O_restorst && !defined O_replacst !----------------------------------------------------------------------- ! damp surface tracers back to data. schematically, the restoring ! term will be = (dampdz/(dampts*daylen)) * (data - t) ! where dampdz is the thickness (cm) and dampts is the time ! scale (days) !----------------------------------------------------------------------- do n=1,nt write (stdout,'(/,1x,a,i2,a,a,/,1x,1pe14.7,a,1pe14.7,a,/)') & ' Surface tracer #',n,' will be damped back to data using a' &, ' Newtonian restoring time scale of ' &, dampts(n),' days. and a level thickness =', dampdz(n),' cm.' enddo # endif # if defined O_shortwave !----------------------------------------------------------------------- ! Solar Shortwave energy penetrates below the ocean surface. Clear ! water assumes energy partitions between two exponentials as ! follows: ! 58% of the energy decays with a 35 cm e-folding scale ! 42% of the energy decays with a 23 m e-folding scale ! if the thickness of the first ocean level "dzt(1)" is 50 meters, ! then shortwave penetration wouldn't matter. however, for ! dzt(1) = 10 meters, the effect can be significant. this can be ! particularly noticeable in the summer hemisphere. ! Paulson and Simpson (1977 Irradiance measurements in the upper ! ocean JPO 7, 952-956) ! Also see ... Jerlov (1968 Optical oceanography. Elsevier) ! A General Circulation Model for Upper Ocean ! Simulation (Rosati and Miyakoda JPO vol 18,Nov 1988) !----------------------------------------------------------------------- write (stdout,'(/a/)') &' => Shortwave penetration is a double exponential as follows:' rpart = 0.58 efold1 = 35.0e0 efold2 = 23.0e2 rscl1 = 1.0/efold1 rscl2 = 1.0/efold2 ! note that pen(0) is set 0.0 instead of 1.0 to compensate for the ! shortwave part of the total surface flux in "stf(i,1)" pen(0) = c0 do k=1,km swarg1 = -min(zw(k)*rscl1,70.0) swarg2 = -min(zw(k)*rscl2,70.0) pen(k) = rpart*exp(swarg1) + (c1-rpart)*exp(swarg2) divpen(k) = (pen(k-1) - pen(k))/dzt(k) write (stdout,9200) k, zw(k)*1.e-2, pen(k), zt(k)*1.e-2 &, divpen(k) enddo write (stdout,*) ' ' 9200 format (1x,'k=',i3,' zw(k)=',f8.2,'(m) pen(k)=',e14.7 &, ' zt(k)=',f8.2,'(m) divpen(k)=',e14.7) # 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 defined O_tracer_averages 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 # 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 # if defined O_symmetry ! do not perform island integrals for land masses whose perimeters ! touch the equator in models symmetric about the equator do isle=1,nisle do n=1,nippts(isle) j = jperm(iofs(isle)+n) if (j .eq. jmt-1) then imask(isle) = .false. endif enddo enddo # endif ! 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 defined O_symmetry ! do not normalize dpsi when using symmetry about the equator imain = 0 # endif 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 = ') # if defined O_stream_function !----------------------------------------------------------------------- ! 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) = ') # endif # if defined O_fourfil || defined O_firfil !----------------------------------------------------------------------- ! 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 # if defined O_cyclic do jrow=1,jmt kmz(1,jrow) = kmz(imtm1,jrow) enddo # endif !--------------------------------------------------------------------- ! 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) # if defined O_rigid_lid_surface_pressure || defined O_implicit_free_surface call findex (kmu, jmtfil, 1, jfu1, jfu2, imt, iszf, iezf) # endif # if defined O_stream_function call findex (kmz, jmtfil, 1, jft1, jft2, imt, iszf, iezf) # endif 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 ==') # endif !--------------------------------------------------------------------- ! 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 # if defined O_rigid_lid_surface_pressure || defined O_implicit_free_surface uhat(i,jrow,1) = c0 uhat(i,jrow,2) = c0 divf(i,jrow) = c0 # endif # if defined O_stream_function ztd(i,jrow) = c0 # endif 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 # if defined O_consthmix # if defined O_biharmonic amix = sqrt(abs(ambi)) # else amix = am # endif 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 # endif ! 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 !----------------------------------------------------------------------- # if defined O_meridional_tracer_budget !----------------------------------------------------------------------- ! set basins and initialize arrays !----------------------------------------------------------------------- if (tmbint .ge. c0) then ! set all points to signify basin # 1 do jrow=1,jmt do i=2,imtm1 msktmb(i,jrow) = 1 enddo msktmb(1,jrow) = 0 msktmb(imt,jrow) = 0 enddo ! divide "msktmb" into other basins (2 .. ntmbb) here if desired. ! set all land points to basin # 0 do jrow=1,jmt do i=1,imt if (kmt(i,jrow) .eq. 0) msktmb(i,jrow) = 0 enddo enddo if (ntmbb .gt. 1) then write (stdout,*) & ' Basin arrangement for Meridional Tracer Budget diagnostic' call iplot (msktmb, imt, imt, jmt) endif ! write out the meridional tracer budget basin mask if ((iotmb .ne. stdout .or. iotmb .lt. 0) .and. itmb) then call getunit (iu, 'tracer_bud.dta' &, 'unformatted sequential append ieee') iotext = & ' read (iotmb) imt, jmt, ((msktmb(i,j),i=1,imt),j=1,jmt)' write (iu) stamp, iotext, expnam write (iu) imt, jmt, ((msktmb(i,j),i=1,imt),j=1,jmt) call relunit (iu) endif endif numtmb = 0 do mask=0,ntmbb do jrow=1,jmt smdvol(jrow,mask) = c0 enddo enddo do mask=0,ntmbb do n=1,nt do jrow=1,jmt tstor(jrow,n,mask) = c0 tdiv(jrow,n,mask) = c0 tflux(jrow,n,mask) = c0 tdif(jrow,n,mask) = c0 tsorc(jrow,n,mask) = c0 enddo enddo enddo # endif # if defined O_bryan_lewis_vertical || defined O_bryan_lewis_horizontal !----------------------------------------------------------------------- ! initialize Bryan_Lewis tracer diffusion coefficients !----------------------------------------------------------------------- call blmixi # endif # if defined O_time_averages !----------------------------------------------------------------------- ! initialize time mean "averaging" grid data !----------------------------------------------------------------------- call avgi # endif # if defined O_xbts !----------------------------------------------------------------------- ! initialize XBT locations and averaging arrays !----------------------------------------------------------------------- call xbti # endif # if defined O_mom_tbt !----------------------------------------------------------------------- ! initialize tbt averaging arrays !----------------------------------------------------------------------- do j=1,jmt do i=1,imt do n=1,nt do k=1,kmt(i,j) do m=1,11 tbt(i,j,k,n,m) = 0. enddo enddo tbtsf(i,j,n) = 0. enddo enddo enddo ntbtts = 0 # endif # if defined O_ppvmix !----------------------------------------------------------------------- ! initialize pacanowski-philander vertical mixing scheme !----------------------------------------------------------------------- call ppmixi (vmixset) # endif # if defined O_smagnlmix && !defined O_consthmix !----------------------------------------------------------------------- ! initialize Smagorinsky nonlinear horizontal mixing scheme !----------------------------------------------------------------------- call smagnli (hmixset) # endif # if defined O_isopycmix !----------------------------------------------------------------------- ! initialize the isopycnal mixing !----------------------------------------------------------------------- call isopi (error, am, ah) # endif # if defined O_npzd ! convert units of NPZD parameters to MOM units redctn = redctn*1.e-3 redotn = redotn*1.e-3 redotp = redotn/redptn redctp = redctn/redptn redntp = 1./redptn k1p = k1n*redptn kw = kw*1.e-2 kc = kc*1.e-2 ki = ki*1.e-2 wd0 = wd0*1.e2 alpha = alpha/daylen abio = abio/daylen nup = nup/daylen nupt0 = nupt0/daylen gbio = gbio/daylen epsbio = epsbio/daylen nuz = nuz/daylen gamma2 = gamma2/daylen nud0 = nud0/daylen ! calculate sinking speed of detritus divided by grid width do k=1,km ! linear increase wd0-200m with depth wd(k) = (wd0+4.0e-2*zt(k))/daylen/dzt(k) ! [s-1] ztt(k) = -zt(k)+ dzt(k)/2. rkwz(k) = 1./(kw*dzt(k)) enddo # if defined O_carbon || defined O_npzd_alk !--------------------------------------------------------------------- ! calculate variables used in calcite remineralization !--------------------------------------------------------------------- rcak(1) = -(exp(-zw(1)/dcaco3)-1.0)/dzt(1) rcab(1) = -1./dzt(1) do k=2,km rcak(k) = -(exp(-zw(k)/dcaco3))/dzt(k) & + (exp(-zw(k-1)/dcaco3))/dzt(k) rcab(k) = (exp(-zw(k-1)/dcaco3))/dzt(k) enddo # endif # endif # if defined O_fwa !----------------------------------------------------------------------- ! calculate areas and masks for fresh water anomalies !----------------------------------------------------------------------- fname = new_file_name ("O_fwawt.nc") inquire (file=trim(fname), exist=exists) if (.not. exists) then print*, "Warning => ", trim(fname), " does not exist." fwawt(:,:) = 0. if (mrfwa .gt. 0 .and. mrfwa .le. nhreg) then where (mskhr(:,:) .eq. mrfwa) fwawt(:,:) = 1. print*, " using regional mask to define area" else print*, " using box corners to define area" fwawt(isfwa1:iefwa1,jsfwa:jefwa) = 1. fwawt(isfwa2:iefwa2,jsfwa:jefwa) = 1. endif else call openfile (fname, iou) ib(:) = 1 ic(:) = 1 ic(1) = imtm2 ic(2) = jmtm2 call getvara ('O_fwawt', iou, imtm2*jmtm2, ib, ic, tmpij &, c1, c0) fwawt(2:imtm1,2:jmtm1) = tmpij(1:imtm2,1:jmtm2) call embmbc (fwawt) endif areafwc = 0. areafwa = 0. do j=2,jmtm1 do i=2,imtm1 if (kmt(i,j) .gt. 0) then areafwa = areafwa + dxt(i)*dyt(j)*cst(j)*fwawt(i,j) fwcwt(i,j) = 1. - fwawt(i,j) areafwc = areafwc + dxt(i)*dyt(j)*cst(j)*fwcwt(i,j) else fwawt(i,j) = 0. fwcwt(i,j) = 0. endif enddo enddo # endif !----------------------------------------------------------------------- ! 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 # if defined O_cyclic do jrow=1,jmt kmu(imt,jrow) = kmu(2,jrow) enddo # endif # if defined O_symmetry do i=1,imt kmu(i,jmt) = kmu(i,jmt-2) enddo # endif !--------------------------------------------------------------------- ! 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 include "size.h" include "param.h" include "pconst.h" include "stdunits.h" include "coord.h" include "csbc.h" include "iounit.h" include "levind.h" include "mw.h" # if defined O_tai_slh include "diag.h" include "state.h" include "dens.h" # endif 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. '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. '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. '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 # if !defined O_idealized_ic 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 # if defined O_reproduce_updates_01 &, c1, c0) # else &, c1e3, c0) # endif 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 # if defined O_reproduce_updates_01 &, c1, c0) # else &, c1e3, c0) # endif 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 # endif enddo # if defined O_carbon && defined O_carbon_14 ! 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)*rstd & *t(i,k,j,idic,taup1) enddo enddo # endif # if defined O_linearized_advection ! initialize density profile with tbarz do n=1,nt do i=1,imt do k=1,km t(i,k,j,n,taup1) = tbarz(k,n) enddo enddo enddo if (jrow .eq. jmt) then write (stdout,'(a,a)') & '=> Note: initialized T & S with idealized profile T(z), S=35ppt' endif # endif # if defined O_tai_slh !----------------------------------------------------------------------- ! 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. # endif # if defined O_sed !----------------------------------------------------------------------- ! initialize bottom boundary for missing ocean tracers from data !----------------------------------------------------------------------- ib(:) = 1 ib(2) = jdi ib(3) = 1 ic(:) = imtm2 ic(2) = 1 ic(3) = km # if !defined O_carbon fname = new_file_name ("O_totcarb.nc") inquire (file=trim(fname), exist=exists) if (exists) then call openfile (trim(fname), iou) call getvara ('O_totcarb', iou, imtm2*km, ib, ic, tmpik &, c1, c0) do i=2,imtm1 k = kmt(i,jrow) sbc(i,jrow,ibdic) = tmpik(i-1,k) enddo call setbcx (sbc(:,jrow,ibdic), imt, 1) else print*, "Error => Can not find ",trim(fname) stop '=> setmom sed dic' endif # endif # if !defined O_npzd_alk fname = new_file_name ("O_alk.nc") inquire (file=trim(fname), exist=exists) if (exists) then call openfile (trim(fname), iou) call getvara ('O_alk', iou, imtm2*km, ib, ic, tmpik, c1, c0) do i=2,imtm1 k = kmt(i,jrow) sbc(i,jrow,ibalk) = tmpik(i,k) enddo call setbcx (sbc(:,jrow,ibalk), imt, 1) else print*, "Error => Can not find ",trim(fname) stop '=> setmom sed alk' endif # endif # if !defined O_npzd_o2 fname = new_file_name ("O_o2.nc") inquire (file=trim(fname), exist=exists) if (exists) then call openfile (trim(fname), iou) call getvara ('O_o2', iou, imtm2*km, ib, ic, tmpik, c1, c0) do i=2,imtm1 k = kmt(i,jrow) sbc(i,jrow,ibo2) = tmpik(i,k) enddo call setbcx (sbc(:,jrow,ibo2), imt, 1) else print*, "Error => Can not find ",trim(fname) stop '=> setmom sed o2' endif # endif # endif !----------------------------------------------------------------------- ! 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 (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 (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 (issdiaz .ne. 0) & sbc(1:imt,jrow,issdiaz) = t(1:imt,1,j,idiaz,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) !----------------------------------------------------------------------- ! 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 (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 !defined O_npzd_no_vflux if (iphyt .ne. 0 .and. issphyt .ne. 0) then call areaavg (sbc(1,1,issphyt), dmsk, gaost(iphyt)) print*, "global average sea surface phytoplankton (mol m-3) = " &, gaost(iphyt)*0.001 endif if (izoop .ne. 0 .and. isszoop .ne. 0) then call areaavg (sbc(1,1,isszoop), dmsk, gaost(izoop)) print*, "global average sea surface zooplankton (mol m-3) = " &, gaost(izoop)*0.001 endif if (idetr .ne. 0 .and. issdetr .ne. 0) then call areaavg (sbc(1,1,issdetr), dmsk, gaost(idetr)) print*, "global average sea surface detritus (mol m-3) = " &, gaost(idetr)*0.001 endif # 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 !defined O_npzd_no_vflux if (idiaz .ne. 0 .and. issdiaz .ne. 0) then call areaavg (sbc(1,1,issdiaz), dmsk, gaost(idiaz)) print*, "global average sea surface diazotrophs (mol m-3) = " &, gaost(idiaz)*0.001 endif # 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 #endif return end