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, c1e9
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_npzd_fe_limitation
real tmpijkm(1:100,1:100,1:3,1:12)
# end if
# 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
coa_grid(:,:)=0
! 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 ('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+6.0e-2*zt(k))/daylen/dzt(k) ! [s-1]
# if defined O_light_bug_in_Ref_Run
ztt(k) = -zt(k)+ dzt(k)/2.
# endif
rkwz(k) = 1./(kw*dzt(k))
enddo
# if defined O_light_bug_in_Ref_Run
!
# else
ztt(1)=0.0
do k=1,km
ztt(k+1)=(-1)*zw(k)
enddo
# endif
# if defined O_pipe_co2_with_fe
kpipe_fe(:,:) = 0.
# endif
# if defined O_npzd_fe_limitation
! read in the dissolved Fe conc from O_dissolved_fe.nc
fe_dissolved(:,:,:,:) = 0.
ib(:) = 1
ic(:) = imtm2
ic(2) = jmtm2
ic(3) = 3
ic(4) = 12
fname = new_file_name ("O_fe_dissolved.nc")
inquire (file=trim(fname), exist=exists)
if (exists) then
c1e9 = 1000000000
call openfile (trim(fname), iou)
call getvara ('O_dissolved_fe', iou, ic(1)*ic(2)*ic(3)*ic(4)
&, ib, ic, tmpijkm, c1e9, c0)
fe_dissolved(2:imtm1,2:jmtm1,:,:) = tmpijkm(1:imtm2
&, 1:jmtm2,:,:)
! print*,'fe_dissolved(50,50,1,1)=',fe_dissolved(50,50,1,1)
do m=1,12
do k=1,3
do j=1,jmt
fe_dissolved(1,j,k,m) = fe_dissolved(imtm1,j,k,m)
fe_dissolved(imt,j,k,m) = fe_dissolved(2,j,k,m)
enddo
do i=1,imt
fe_dissolved(i,1,k,m) = fe_dissolved(i,2,k,m)
fe_dissolved(i,jmt,k,m) = fe_dissolved(2,j,k,m)
enddo
enddo
enddo
else
print*,"Warning => Cannot find", trim(fname)
endif
# endif
# 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
&, 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
# 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