subroutine setvbc (joff, js, je, is, ie) !======================================================================= ! set momentum and tracer vertical boundary conditions ! input: ! joff = offset relating "j" in the MW to latitude "jrow" ! js = starting row in the MW ! je = ending row in the MW ! is = starting longitude index in the MW ! ie = ending longitude index in the MW !======================================================================= implicit none integer js, je, istrt, is, iend, ie, n, j, i, jrow, joff, kz real uvmag !======================= include file "size.h" ========================= !----------------------------------------------------------------------- ! USER INPUT: !----------------------------------------------------------------------- ! imt = number of grid points in the longitudinal direction ! (calculated points are from 2 through imt-1. End points ! are boundaries) ! jmt = number of grid points (latitude rows) in the latitudinal ! direction (calculated points are from 2 through jmt-1. ! End points are boundaries) ! km = number of grid points in the vertical direction ! (calculated points are from 1 through km) ! nt = number of tracers (temperature, salinity, ...) ! nsrc = number of tracer with sources ! kpzd = depth for limited npzd model ! jmz = size for "unrotated" zonal averages ! jmzm1 = jmz minus one ! mnisle = maximum number of islands (unconnected land masses) ! maxipp = maximum number of all island perimeter points !----------------------------------------------------------------------- integer imt, jmt, km, nt, nsrc, kpzd, nat, jmz, jmzm1, mnisle integer maxipp, jmw, jsmw, jemw parameter (imt= 1, jmt= 1, km= 19) parameter (nt=2 $ +1 $ +1 $ +1 $ +4 $ +2 $ +2 $ ) parameter (nsrc=0 $ +1 $ +1 $ +1 $ +4 $ +2 $ +2 $ ) parameter (kpzd=km) parameter (nat=2 $, jmz=jmt, jmzm1=jmz-1) parameter (mnisle=50, maxipp=5000) parameter (jmw=jmt) !----------------------------------------------------------------------- ! set first and last calculated row within the MW. other rows ! are used as buffers !----------------------------------------------------------------------- ! jsmw = 1st calculated row within the MW ! jemw = last calculated row within the MW parameter (jsmw=1, jemw=1) ! Moses-Triffid land model ! POINTS = Maximum number of points in grid. ! STEPSM = Maximum number of timesteps in a day. ! klmax = maximum ocean depth levels over which the land model can exist integer POINTS, STEPSM, klmax parameter (POINTS=14300, STEPSM=24, klmax=0) ! NNVG = Number of non-vegetation surface types. ! NPFT = Number of plant functional types. ! NTYPE = Number of surface types. ! SOIL = Index of the surface type 'Soil' ! Land surface types : ! 1 - Broadleaf Tree ! 2 - Needleleaf Tree ! 3 - C3 Grass ! 4 - C4 Grass ! 5 - Shrub ! 6 - Soil integer NNVG, NPFT, NTYPE, SOIL parameter (NNVG=4, NPFT=5, NTYPE=6, SOIL=6) !======================= include file "param.h" ======================== ! main parameter file which sets ocean characteristics: ! nvar = number of prognostic variables ! lseg = maximum number of longitudinal stream function segments ! nlatpr = maximum number of latitudes for matrix printouts ! on diagnostic time steps ! nhreg = number of regions in the horizontal used for averaging ! tracers. ! nvreg = number of regions in the vertical used for term balance ! calculations. note "nvreg" is not used for tracer ! averages ! numreg = total number of regions ( = product of nhreg & nvreg) ! used for term balance calculations ! nvarbh = number of prognostic variables using biharmonic mixing ! ncrows = number of calculated rows within the MW. ! (the remaining rows are buffer rows). ! rstd = standard c14/c12 ratio integer lseg, nlatpr, nhreg, nvreg, numreg, nvar, nvarbh integer imtm1, kmm1, imtp1, imtm2, jmtp1, jmtm1, jmtm2, jscan integer kmp1, kmp2, imtkm, nwds, nkflds, nslab, ntmin2, ncrows real rstd parameter (lseg=5, nlatpr=10) parameter (nhreg=3, nvreg=1, numreg=nhreg*nvreg) parameter (nvar=nt+2) parameter (imtm1=imt-1, kmm1=km-1) parameter (imtp1=imt+1, imtm2=imt-2 &, jmtp1=jmt+1, jmtm1=jmt-1, jmtm2=jmt-2 &, jscan=jmtm2 &, kmp1=km+1, kmp2=km+2 &, imtkm=imt*km, nwds=imt*jmt, nkflds=2 &, nslab=imt*nvar*km, ntmin2=nt+1/nt) parameter (ncrows = jmw - 3 + jmw/jmt) !====================== include file "pconst.h" ======================== ! rules for parameter constants ! use prefix of "c" for whole real numbers (eg: c57 for 57.0) ! use "m" after prefix to designate negative values (minus sign) ! (eg: cm7 for -7.0) ! use prefix of "p" for non repeating fractions (eg: p5 for 0.5) ! use prefix of "r" for reciprocals (eg: r3 for 1/3.0) ! combine use of prefix above and "e" for scientific notation, with ! (eg: c5e4 for 5.0e4, c1em10 for 1.0e-10) real c0, c1, c2, c3, c4, c5, c7, c8, c14, c16, c360, p125, p25 real p5, p75, epsln, c24, c60, c1440, r24, r60, r1440, secday parameter (c0=0.0, c1=1.0, c2=2.0, c3=3.0, c4=4.0, c5=5.0, c7=7.0) parameter (c8=8.0) parameter (c14=14.0, c16=16.0, c360=360.0) parameter (p125=0.125, p25=0.25, p5=0.5, p75=0.75) parameter (epsln=1.0e-20) parameter (c24=24.0, c60=60.0, c1440=1440.0) parameter (r24=c1/c24, r60=c1/c60, r1440=c1/c1440) parameter (secday=c1/(c60*c1440)) !====================== include file "stdunits.h" ====================== ! stdin = unit number for standard input. ! stdout = unit number for standard output. ! stderr = unit number for standard error. integer stdin, stdout, stderr parameter (stdin = 5, stdout = 6, stderr = 6) !======================= include file "coord.h" ======================== ! model grid point coordinates ! grid definition: ! the model uses a staggered arakawa "b" grid which is setup and ! generated by the "grids.F" module. ! xt(i) = longitude of the ith "t" point in degrees. i=1..imt ! xu(i) = longitude of the ith "u,v" point in degrees. i=1..imt ! yt(j) = latitude of the jth "t" point in degrees. j=1..jmt ! yu(j) = latitude of the jth "u,v" point in degrees. j=1..jmt ! zt(k) = distance from surface down to centre of level k (in cm) ! (for depth of "t" and "u,v" grid points: k=1,km) ! zw(k) = distance from surface down to bottom of level k (in cm) ! (for depth of "t" and "u,v" grid points: k=1,km) ! dxtdeg = widths for "t" grid cells (degrees) ! dytdeg = heights for "t" grid cells (degrees) ! dxudeg = widths for "u" grid cells (degrees) ! dyudeg = heights for "u" grid cells (degrees) ! dzt(k) = vertical resolution of "t" and "u" grid cells (in cm) ! dzw(k) = vertical resolution of "w" grid cells (in cm) ! "i" increases in an eastward direction, "j" increases in a ! northward direction, and "k" increases downward. real xt, yt, xu, yu, zw, zt real dxtdeg, dytdeg, dzt real dxudeg, dyudeg, dzw common /coord/ xt(imt), yt(jmt), xu(imt), yu(jmt), zw(km), zt(km) common /coord/ dxtdeg(imt), dytdeg(jmt), dzt(km) common /coord/ dxudeg(imt), dyudeg(jmt), dzw(0:km) !======================= include file "csbc.h" ========================= ! surface boundary conditions (S.B.C.) ! numsbc = total number of surface boundary conditions ! indices and names set in UVic_ESCM.F: ! itaux = x component of wind stress (dynes cm-2) ! itauy = y component of wind stress (dynes cm-2) ! iws = surface wind speed (cm s-1) ! iaca = atmospheric coalbedo ! isca = surface coalbedo (%/100) ! ihflx = heat flux ! isflx = salt flux ! isst = ocean model SST (C) ! isss = ocean model SSS (PSU*0.001-0.035) ! iwa = surface wind angle (degrees) ! iro = surface runoff (g cm-2 s-1) ! iwxq = x component of wind for moisture advection (cm s-1) ! iwyq = y component of wind for moisture advection (cm s-1) ! iwxt = x component of wind for temperature advection (cm s-1) ! iwyt = y component of wind for temperature advection (cm s-1) ! iwxc = x component of wind for carbon advection (cm s-1) ! iwyc = y component of wind for carbon advection (cm s-1) ! ipsw = penetrating shortwave (cal cm-2 s-1) ! isu = x component of ocean first layer velocity (cm s-1) ! isv = y component of ocean first layer velocity (cm s-1) ! igu = x component of ocean second layer velocity (cm s-1) ! igv = y component of ocean second layer velocity (cm s-1) ! issdic = sea surface concentration of dic (umol cm-3) ! idicflx = sea surface flux of carbon (umol cm-2 s-1) ! issalk = sea surface concentration of alkalinity (umol cm-3) ! ialkflx = sea surface flux of alkalinity (umol cm-2 s-1) ! isso2 = sea surface concentration of oxygen (umol cm-3) ! io2flx = sea surface flux of oxygen (umol cm-2 s-1) ! isspo4 = sea surface concentration of phosphate (nmol cm-3) ! ipo4flx = sea surface flux of phosphate (nmol cm-2 s-1) !{ wyao 20170704 ! wyao 20170704} ! issphyt = sea surface concentration of phytoplankton (nmol P cm-3) ! iphytflx = sea surface flux of phytoplankton (nmol P cm-2 s-1) ! isszoop = sea surface concentration of zooplankton (nmol P cm-3) ! izoopflx = sea surface flux of zooplankton (nmol P cm-2 s-1) ! issdetr = sea surface concentration of detritus (nmol P cm-3) ! idetrflx = sea surface flux of detritus (nmol P cm-2 s-1) ! issdetrfe = sea surface concentration of particualte Fe (nmol Fe cm-3) ! idetrfeflx= sea surface flux of particulate Fe (nmol Fe cm-2 s-1) ! issno3 = sea surface concentration of nitrate (nmol P cm-3) ! ino3flx = sea surface flux of nitrate (nmol P cm-2 s-1) !{ wyao 20170704 ! wyao 20170704} ! issdfe = sea surface concentration of iron (nmol Fe cm-3) ! idfeflx = sea surface flux of iron (nmol Fe cm-2 s-1) ! idfeadep = sea surface flux of iron (nmol Fe cm-2 s-1) ! issdiaz = sea surface concentration of diazotraphs (nmol P cm-3) ! idiazflx = sea surface flux of diazotraphs (nmol P cm-2 s-1) ! issc14 = sea surface concentration of carbon 14 (umol cm-3) ! ic14flx = sea surface flux of carbon 14 (umol cm-2 s-1) ! isscfc11 = sea surface concentration of cfc11 (umol cm-3) ! icfc11flx = sea surface flux of cfc11 (umol cm-2 s-1) ! isscfc12 = sea surface concentration of cfc12 (umol cm-3) ! icfc12flx = sea surface flux of cfc12 (umol cm-2 s-1) ! iat = surface air temperature (C) ! irh = surface relative humidity (%/100) ! ipr = precipitation as rain (g cm-2 s-1) ! ips = precipitation as snow (g cm-2 s-1) ! iaws = averaged surface wind speed (cm s-1) ! iswr = surface shortwave radiation (g s-3) ! ilwr = surface longwave radiation (g s-3) ! isens = surface sensible heat flux (g s-3) ! ievap = surface evaporation (g cm-2 s-1) ! idtr = diurnal temperature range (C) ! inpp = carbon flux from net primary production (kg C m-2 s-1) ! isr = carbon flux from soil respiration (kg m-2 s-1) ! iburn = carbon flux from burning vegetation (kg m-2 s-1) ! ibtemp = ocean bottom temperature (C) ! ibsalt = ocean bottom temperature (PSU*0.001-0.035) ! ibdic = ocean bottom dissolved organic carbon (umol cm-3) ! ibdicfx = ocean bottom flux of carbon (umol cm-2 s-1) ! ibalk = ocean bottom alkalinity (umol cm-3) ! ibalkfx = ocean bottom flux of alkalinity (umol cm-2 s-1) ! ibo2 = ocean bottom oxygen (umol cm-3) ! ircal = rain rate of calcite (umol cm-2 s-1) ! irorg = rain rate of organic carbon (umol cm-2 s-1) ! mapsbc = surface boundary conditions names ! sbc = surface boundary condition data. ! gaost = global average ocean surface tracer ! socn = average ocean sea surface salinity ! ntspas = the number of ocean time steps per atmosphere segment ! ntspls = the number of ocean time steps per land segment ! ntspos = the number of ocean time steps per ocean segment ! nats = number of atm time steps since mixing time step ! nots = number of ocn time steps since mixing time step ! addflxa = logical flag for adding only even mode fluxes from atm ! addflxo = logical flag for adding only even mode fluxes from ocn ! dtatm = atmosphere time step (s) ! dtism = icesheet time step (s) ! dtismyr = icesheet time step (years) ! nismacc = icesheet acceleration time (years) ! dtlnd = land time step (s) ! dtocn = ocean time step (s) ! dtsed = sediment time step (s) ! dtsedyr = sediment time step (years) ! nsedacc = sediment acceleration time (years) ! atatm = atmosphere boundary condition averaging time (s) ! atlnd = land boundary condition averaging time (s) ! atocn = ocean boundary condition averaging time (s) ! atsed = sediment boundary condition averaging time (s) ! atism = icesheet boundary condition averaging time (s) ! dampts = time scale for damping surface tracers (days) to data ! dampdz = ocean level thickness for converting Newtonian damping ! to a surface flux ! land_map = map with indices for coupling to land arrays ! dtoih = total system heat lost minus heat gained ! dtoic = total system carbon lost minus carbon gained ! avgpertavg = averaging period for time averages ! avgtimtavg = averaging time for time averages ! avgpertsi = averaging period for time integrals ! avgtimtsi = averaging time for time integrals integer numsbc parameter (numsbc = 14 & + 4 & + 2 & + 2 & + 2 & + 2 & + 2 & + 3 & ) integer itaux, itauy, iws, iaca, isca, ihflx, isflx, isst, isss integer iwa, iro, iwxq, iwyq, iwxt, iwyt, iwxc, iwyc, ipsw, isu integer isv, igu, igv, issdic, idicflx, issalk, ialkflx, isso2 integer io2flx, isspo4, ipo4flx, issphyt, iphytflx, isszoop integer izoopflx, issdetr, idetrflx, issno3, ino3flx, issdiaz integer idiazflx, issc14, ic14flx, isscfc11, icfc11flx, isscfc12 integer icfc12flx, iat, irh, ipr, ips, iaws, iswr, ilwr, isens integer ievap, idtr, inpp, isr, iburn, ibtemp, ibsalt, ibdic integer ibdicfx, ibalk, ibalkfx, ibo2, ircal, irorg !{ wyao 20170704 integer issdfe ! wyao 20170704} integer idfeflx, issdetrfe, idetrfeflx, idfeadep common /csbc_i/ itaux, itauy, iws, iaca, isca, ihflx, isflx, isst common /csbc_i/ isss, iwa, iro, iwxq, iwyq, iwxt, iwyt, iwxc, iwyc common /csbc_i/ ipsw, isu, isv, igu, igv, issdic, idicflx, issalk common /csbc_i/ ialkflx, isso2, io2flx, isspo4, ipo4flx, issphyt common /csbc_i/ iphytflx, isszoop, izoopflx, issdetr, idetrflx common /csbc_i/ issno3, ino3flx, issdiaz, idiazflx, issc14 common /csbc_i/ ic14flx, isscfc11, icfc11flx, isscfc12, icfc12flx common /csbc_i/ iat, irh, ipr, ips, iaws, iswr, ilwr, isens, ievap common /csbc_i/ idtr, inpp, isr, iburn, ibtemp, ibsalt, ibdic common /csbc_i/ ibdicfx, ibalk, ibalkfx, ibo2, ircal, irorg !{ wyao 20170704 ! wyao 20170704} common /csbc_i/ issdfe, idfeflx, issdetrfe, idetrfeflx common /csbc_i/ idfeadep character(20) :: mapsbc common /csbc_c/ mapsbc(numsbc) real sbc common /csbc_r/ sbc(imt,jmt,numsbc) real gaost, socn common /csbc_r/ gaost(nt), socn integer ntspas, ntspls, ntspos, nats, nots, nismacc, nsedacc common /csbc_i/ ntspas, ntspls, ntspos, nats, nots, nismacc common /csbc_i/ nsedacc logical addflxa, addflxo common /csbc_l/ addflxa, addflxo real dtatm, dtism, dtismyr, dtlnd, dtocn, dtsed, dtsedyr, atatm real atism, atlnd, atocn, atsed common /csbc_r/ dtatm, dtism, dtismyr, dtlnd, dtocn, dtsed common /csbc_r/ dtsedyr, atatm, atism, atlnd, atocn, atsed real dampts, dampdz, subflux, subz common /csbc_r/ dampts(nt), dampdz(nt) integer ntlbc common /csbc_i/ ntlbc real dtoih, dtoic common /csbc_r/ dtoih, dtoic real avgpertavg, avgtimtavg, avgpertsi, avgtimtsi common /csbc_r/ avgpertavg, avgtimtavg, avgpertsi, avgtimtsi !====================== include file "grdvar.h" ======================== ! variables which are functions of the grid defined by "coord.h" ! dxt = longitudinal width of "t" grid box at the ! equator (in cm) ! dxtr = reciprocal of "dxt" ! dxt2r = reciprocal of "2*dxt" ! dxt4r = reciprocal of "4*dxt" ! dxu = longitudinal width of "u,v" grid box at the ! equator (in cm) ! dxur = reciprocal of "dxu" ! dxu2r = reciprocal of "2*dxu" ! dxu4r = reciprocal of "4*dxu" ! dxmetr = reciprocal of "(dxt(i)+dxt(i+1))" ! duw = xu(i) - xt(i) ! due = xt(i+1) - xu(i) ! dus = yu(jrow) - yt(jrow) ! dun = yt(jrow+1) - yu(jrow) ! dyt = latitudinal height of "t" grid box (in cm) ! dytr = reciprocal of "dyt" ! dyt2r = reciprocal of "2*dyt" ! dyt4r = reciprocal of "4*dyt" ! dyu = latitudinal height of "u,v" grid box (in cm) ! dyur = reciprocal of "dyu" ! dyu2r = reciprocal of "2*dyu" ! dyu4r = reciprocal of "4*dyu" ! csu = cosine of "u,v" grid point latitude ! csur = reciprocal of "csu" ! cst = cosine of "t" grid point latitude ! cstr = reciprocal of "cst" ! phi = latitude of "u,v" grid point in radians ! phit = latitude of "t" grid point in radians ! sine = sine of "u,v" grid point latitude ! tng = tan of "u,v" grid point latitude ! c2dzt(k)= "2*dzt" ! dztr(k) = reciprocal of dzt ("t" cell vertical resolution) ! dzt2r(k)= reciprocal of "2*dzt" ! dzwr(k) = reciprocal of dzw ("w" cell vertical resolution) ! dzw2r(k)= reciprocal of "2*dzw" ! dztur(k)= upper diffusion grid factor = 1.0/(dzw(k-1)*dzt(k)) ! dztlr(k)= lower diffusion grid factor = 1.0/(dzw(k)*dzt(k)) ! tlat = tracer grid geographic latitude (degrees) ! ulat = velocity grid geographic latitude (degrees) ! tlon = tracer grid geographic longitude (degrees) ! ulon = velocity grid geographic longitude (degrees) ! tgarea = tracer grid geographic area (cm2) ! ugarea = velocity grid geographic area (cm2) real dxt, dxtr, dxt2r, dxu, dxur, dxu2r, dxu4r, dxt4r real dyt, dytr, dyt2r, dyu, dyur, dyu2r, dyu4r, dyt4r real csu, csur, cst, cstr, cstdytr, cstdyt2r real csudyur, csudyu2r, cst_dytr, csu_dyur real phi, phit, sine, tng, c2dzt, dztr, dzt2r real dzwr, dzw2r, dxmetr, duw, due, dun, dus, dztur, dztlr real quick_x, curv_xp, curv_xn, quick_y, curv_yp, curv_yn real quick_z, curv_zp, curv_zn real tlat, tlon, ulat, ulon, tgarea, ugarea common /grdvar/ dxt(imt), dxtr(imt), dxt2r(imt), dxu(imt) common /grdvar/ dxur(imt), dxu2r(imt), dxu4r(imt), dxt4r(imt) common /grdvar/ dyt(jmt), dytr(jmt), dyt2r(jmt), dyu(jmt) common /grdvar/ dyur(jmt), dyu2r(jmt), dyu4r(jmt), dyt4r(jmt) common /grdvar/ csu(jmt), csur(jmt), cst(jmt), cstr(jmt) common /grdvar/ cstdytr(jmt), cstdyt2r(jmt) common /grdvar/ csudyur(jmt), csudyu2r(jmt) common /grdvar/ cst_dytr(jmt), csu_dyur(jmt) common /grdvar/ phi(jmt), phit(jmt), sine(jmt), tng(jmt) common /grdvar/ c2dzt(km), dztr(km), dzt2r(km) common /grdvar/ dzwr(0:km), dzw2r(0:km) common /grdvar/ dxmetr(imt), duw(imt), due(imt) common /grdvar/ dun(jmt), dus(jmt) common /grdvar/ dztur(km), dztlr(km) common /grdvar/ tlat(imt,jmt), tlon(imt,jmt), tgarea(imt,jmt) common /grdvar/ ulat(imt,jmt), ulon(imt,jmt), ugarea(imt,jmt) ! tcella = "t" cell surface area cm**2 (entire ocean) ! ucella = "u" cell surface area cm**2 (entire ocean) ! tcellv = "t" cell volume cm**3 (entire ocean) ! ucellv = "u" cell volume cm**3 (entire ocean) real tcellv, ucellv, tcella, ucella common /grdvar/ tcellv, ucellv, tcella(km), ucella(km) !====================== include file "levind.h" ======================== ! vertical level indicators which define model geometry & bottom ! topography: ! kmt = number of vertical boxes over "t" points ! kmu = number of vertical boxes over "u,v" points integer kmt, kmu real sg_bathy common /levind/ kmt(imt,jmt), kmu(imt,jmt) common /levind/ sg_bathy(imt,jmt,km) !====================== include file "scalar.h" ======================== ! various scalar quantities: ! dtts = time step for density & tracers (in seconds) ! dtuv = time step for baroclinic velocity (in seconds) ! dtsf = time step for barotropic velocity (in seconds) ! c2dtts = 2*dtts ! c2dtuv = 2*dtuv ! c2dtsf = 2*dtsf ! acor = (>0, 0) = (implicit, explicit) treatment of Coriolis ! term for internal and external modes. ! rho0 = mean density for Bousinessq approximation ! rho0r = 1/rho0 ! omega = earth`s rotation rate (radians/sec) ! radius = earth`s radius (cm) ! grav = earth`s gravitational acceleration (cm/sec**2) ! cdbot = bottom drag coefficient ! ncon = number of passes through convective code in tracer ! gcor = time centring for Coriolis term ! taux0 = constant zonal windstress (dynes/cm**2) for idealized ! equatorial studies ! tauy0 = constant meridional windstress (dynes/cm**2) for ! idealized equatorial studies integer ncon common /scalar_i/ ncon real dtts, dtuv, dtsf, c2dtts, c2dtuv, c2dtsf, acor, rho0 real rho0r, omega, radius, grav, cdbot, gcor, taux0, tauy0 common /scalar_r/ dtts, dtuv, dtsf, c2dtts, c2dtuv, c2dtsf, acor common /scalar_r/ rho0, rho0r, omega, radius, grav, cdbot, gcor common /scalar_r/ taux0, tauy0 ! various non dimensional quantities: ! radian = degrees per radian ! pi = something good to eat real radian, pi common /ndcon/ radian, pi !======================= include file "mw.h" =========================== ! M E M O R Y W I N D O W ! Data layout on disk ! On disk, a latitude consists of a row of prognostic quantities. ! Latitudes are ordered from south to north on two logical disk ! units: one for time level "tau" and one for time level "tau-1". ! Newly computed "tau+1" quantities are written to the "tau-1" disk ! which serves double duty as a "tau+1" disk. ! Data layout in memory ! A memory window "MW" is constructed to hold a group of "jmw" ! adjacent latitude rows of prognostic quantities from disks "tau-1" ! and "tau". Parameter "jmw" controls the size of the MW and can be ! set anywhere from a minimum of three latitude rows to all "jmt" ! latitudes. In addition, the MW holds diagnostic quantities ! (density, hydrostatic pressure gradients, and advective ! velocities) along with work arrays for constructing intermediate ! fluxes used in solving the tracer and momentum equations. A ! latitude row in the MW is referred to by index "j" and corresponds ! to the physical latitude row "jrow" on disk. ! Data flow between disk and memory ! The MW is loaded with prognostic data from the first "jww" ! latitude rows on disks "tau" and "tau-1". As the tracer and ! momentum equations are solved for rows j=2 through jmw-1 in the ! MW, the solutions are written to disk "tau+1". When finished, the ! MW is moved northward by moving quantities from the northernmost ! two rows into the southernmost two rows within the MW. The ! remaining MW rows are then loaded with more latitude rows from ! disk. The process continues until all latitude rows 2 through ! jmt-1 on disk "tau+1" have been updated. !======================================================================= ! taum1 = tau-1 time level for variables in MW ! tau = tau time level for variables in MW ! taup1 = tau+1 time level for variables in MW integer taum1, tau, taup1 common /mw_i/ taum1, tau, taup1 !----------------------------------------------------------------------- ! MW arrays for prognostic equations: !----------------------------------------------------------------------- ! u(i,k,j,n,tau) = total velocity where: ! i = index for longitude ! k = index for depth ! j = index for latitude row within MW ! n = component (1 = zonal, 2 = meridional) ! tau = time level (tau-1, tau, tau+1) ! (only internal modes are on disk and at tau+1 in the MW) ! t(i,k,j,n,tau) = tracer where: ! i = index for longitude ! k = index for depth ! j = index for latitude row within MW ! n = component (1 = temperature, 2 = salinity) ! if nt > 2 then other tracers are allowed. ! tau = time level (tau-1, tau, tau+1) ! note: temperature is potential temperature in degrees Celsius and ! salinity is in "model units", the deviation from 0.035 grams ! of salt/cm**3 of water, or, assuming a water density of ! 1 gram/cm**3, the deviation from 0.035 g of salt/g of water. ! one can convert model units to the more common units of parts ! per thousand (ppt) by adding 0.035 grams/cm**3 to the model ! units and then multiplying by 1000. real u,t common /mw_r/ u(imt,km,jmw,2,-1:1), t(imt,km,jmw,nt,-1:1) ! indicies for ocean tracer array ! itemp = index for temperature ! isalt = index for salinity ! idic = index for dissolved inorganic carbon ! ic14 = index for carbon 14 ! icfc11 = index for cfc11 ! icfc12 = index for cfc12 ! ialk = index for alkalinity ! ipo4 = index for phosphate ! iphyt = index for phytoplankton ! izoop = index for zooplankton ! idetr = index for detritus ! idfe = index for iron ! idetrfe= index for particulate iron ! io2 = index for oxygen ! ino3 = index for nitrate ! idiaz = index for diazotrophs ! mapt = map for ocean tracer names ! { wyao 20170627 ! wyao 20170627} character(10) :: mapt common /mw_c/ mapt(nt) integer itemp, isalt, idic, ic14, icfc11, icfc12, ialk, ipo4 integer iphyt, izoop, idetr, io2, ino3, idiaz common /mw_i/ itemp, isalt, idic, ic14, icfc11, icfc12, ialk common /mw_i/ ipo4, iphyt, izoop, idetr, io2, ino3, idiaz integer idfe, idetrfe common /mw_i/ idfe, idetrfe ! { wyao 20170627 ! wyao 20170627} ! indicies for ocean tracer source array ! istemp = index for temperature ! issalt = index for salinity ! isdic = index for carbon ! isc14 = index for carbon 14 ! isalk = index for alkalinity ! ispo4 = index for phosphate ! isphyt = index for phytoplankton ! iszoop = index for zooplankton ! isdetr = index for detritus ! isdfe = index for iron ! isdetrfe= index for particulate iron ! { wyao 20170627 ! wyao 20170627} ! itrc = index of tracer sources for all tracers (0 = no source) ! mapst = map for ocean tracer source names character(10) :: mapst common /mw_c/ mapst(nt) integer itrc common /mw_i/ itrc(nt) integer istemp, issalt, isdic, isc14, isalk, ispo4, isphyt integer iszoop, isdetr, iso2, isno3, isdiaz common /mw_i/ istemp, issalt, isdic, isc14,isalk, ispo4, isphyt common /mw_i/ iszoop, isdetr, iso2, isno3, isdiaz integer isdfe, isdetrfe common /mw_i/ isdfe, isdetrfe ! { wyao 20170627 ! wyao 20170627} !----------------------------------------------------------------------- ! MW arrays for diagnostic equations and workspace: !----------------------------------------------------------------------- ! diagnostic advective velocities are in units of cm/sec ! adv_vet = advective velocity on the eastern face of a "T" cell ! adv_vnt = advective velocity on the northern face of a "T" cell ! adv_veu = advective velocity on the eastern face of a "u" cell ! adv_vnu = advective velocity on the northern face of a "u" cell ! adv_vbt = advective velocity on the bottom face of a "T" cell ! adv_vbu = advective velocity on the bottom face of a "u" cell ! rho = density at centre of a "T" cell in units of gm/cm**3 ! note: there is an arbitrary constant which is only a ! function of depth in "rho". It is related to ! subtracting a reference level density for purposes of ! accuracy. ! grad_p = hydrostatic pressure gradient for "u" cell. There are ! two components: (1,2) is for (dp/dx, dp/dy) real adv_vet, adv_vnt, adv_veu, adv_vnu, adv_vbt, adv_vbu, rho real rhotaum1, rhotaup1, rhotilde, grad_p common /mw_r/ adv_vet(imt,km,jsmw:jmw), adv_vnt(imt,km,1:jmw) common /mw_r/ adv_veu(imt,km,jsmw:jemw) common /mw_r/ adv_vnu(imt,km,1:jemw) common /mw_r/ adv_vbt(imt,0:km,jsmw:jmw) common /mw_r/ adv_vbu(imt,0:km,jsmw:jemw) common /mw_r/ rho(imt,km,jsmw:jmw) common /mw_r/ grad_p(imt,km,jsmw:jemw,2) ! tmask = tracer cell land/sea mask = (0.0, 1.0) on (land, sea) ! umask = velocity cell land/sea mask = (0.0, 1.0) on (land, sea) real tmask, umask common /mw_r/ tmask(imt,km,1:jmw), umask(imt,km,1:jmw) ! these workspace arrays are recalculated for each component of the ! equations so do not have to be moved as the MW moves northward. ! adv_fe = advective flux across the eastern face of a cell ! adv_fn = advective flux across the northern face of a cell ! (removed in most cases and put directly into the ! statement functions for speed optimization.) ! adv_fb = advective flux across the bottom face of a cell ! diff_fe = diffusive flux across the eastern face of a cell ! diff_fn = diffusive flux across the northern face of a cell ! diff_fb = diffusive flux across the bottom face of a cell ! source = source term real adv_fe common /mw_r/ adv_fe(imt,km,jsmw:jemw) real adv_fn common /mw_r/ adv_fn(imt,km,1:jemw) real adv_fb, diff_fe common /mw_r/ adv_fb(imt,0:km,jsmw:jemw) common /mw_r/ diff_fe(imt,km,jsmw:jemw) real diff_fn common /mw_r/ diff_fn(imt,km,1:jemw) real diff_fb common /mw_r/ diff_fb(imt,0:km,jsmw:jemw) real diff_fbiso common /mw_r/ diff_fbiso(imt,0:km,jsmw:jemw) real source common /mw_r/ source(imt,km,jsmw:jemw) real zzi common /mw_r/ zzi(imt,km,jsmw:jemw) ! these grid factors are for optimizations and are recalculated as ! the MW moves northward so they do not have to be moved. real cstdxtr, cstdxur, cstdxt2r, ah_cstdxur, csudxur, csudxtr real csudxu2r, am_csudxtr common /mw_r/ cstdxtr(imt,jsmw:jmw), cstdxur(imt,jsmw:jmw) common /mw_r/ cstdxt2r(imt,jsmw:jmw), ah_cstdxur(imt,jsmw:jmw) common /mw_r/ csudxur(imt,jsmw:jmw) common /mw_r/ csudxtr(imt,jsmw:jmw) common /mw_r/ csudxu2r(imt,jsmw:jmw) common /mw_r/ am_csudxtr(imt,km,jsmw:jmw) ! these variables are either constant or globally dimensioned by ! "jmt", so they do not need to be moved as the MW moves northward ! advmet = coeff for metric advection. ! cori = Coriolis parameter for velocity component "n" real advmet, cori, cf common /advec_r/ advmet(jmt,2) common /coriol_r/ cori(imt,jmt,2) common /extwrk_r/ cf(imt,jmt,-1:1,-1:1) ! smf = surface momentum flux ! 1 => zonal wind stress (dynes/cm**2) ! 2 => meridional wind stress (dynes/cm**2) ! bmf = bottom momentum flux ! 1 => zonal bottom drag (dynes/cm**2) ! 2 => meridional bottom drag (dynes/cm**2) ! stf = surface tracer flux ! 1 => surface heat flux (cal/cm**2/sec = cm*degC/sec = ly/sec) ! 2 => surface salt flux (grams of salt/cm**2/sec) ! btf = bottom tracer flux (for consistency but normally zero!) ! 1 => bottom heat flux (cal/cm**2/sec = cm*degC/sec = ly/sec) ! 2 => bottom salt flux (grams of salt/cm**2/sec) real smf, bmf, stf, btf common /mw_r/ smf(imt,1:jmw,2), bmf(imt,1:jmw,2) common /mw_r/ stf(imt,1:jmw,nt), btf(imt,1:jmw,nt) ! "antidiffusive" flux as in Zalesak, 1989( see FCTstm.h) ! same for R+, R- ! anti_fe = antidiffusive flux across the eastern face of a T cell ! anti_fn = antidiffusive flux across the northern face of a T cell ! anti_fb = antidiffusive flux across the bottom face of a T cell ! R_plusY = ratio of maximal feasible to maximal possible change ! of tracer T in subroutine tracer.F, N-S dimension delimiter ! R_minusY = ratio of minimal feasible to minimal possible change ! of tracer T in subroutine tracer.F, N-S dimension delimiter real anti_fe, anti_fn, anti_fb, R_plusY, R_minusY common /mw_r/ anti_fe(imt,km,jsmw:jmw,nt) common /mw_r/ anti_fn(imt,km,1:jmw-1+jmw/jmt,nt) common /mw_r/ anti_fb(imt,0:km,jsmw:jmw,nt) common /mw_r/ R_plusY(imt,km,1:jmw-1+jmw/jmt,nt) common /mw_r/ R_minusY(imt,km,1:jmw-1+jmw/jmt,nt) !----------------------------------------------------------------------- ! bail out if starting row exceeds ending row !----------------------------------------------------------------------- if (js .gt. je) return istrt = is iend = ie !---------------------------------------------------------------------- ! set no flux condition for all tracers at surface & bottom. !---------------------------------------------------------------------- do n=1,nt do j=js,je do i=istrt,iend stf(i,j,n) = c0 btf(i,j,n) = c0 enddo enddo enddo !---------------------------------------------------------------------- ! apply surface tracer and momentum fluxes from the atmosphere ! code is for 2 tracer and 2 momentum fluxes. !---------------------------------------------------------------------- do j=js,je jrow = j + joff do i=istrt,iend stf(i,j,itemp) = sbc(i,jrow,ihflx)*tmask(i,1,j) stf(i,j,isalt) = sbc(i,jrow,isflx)*tmask(i,1,j) stf(i,j,idic) = sbc(i,jrow,idicflx)*tmask(i,1,j) stf(i,j,ialk) = sbc(i,jrow,ialkflx)*tmask(i,1,j) stf(i,j,io2) = sbc(i,jrow,io2flx)*tmask(i,1,j) stf(i,j,ipo4) = sbc(i,jrow,ipo4flx)*tmask(i,1,j) !{ wyao 20170714 ! wyao 20170714} stf(i,j,ino3) = sbc(i,jrow,ino3flx)*tmask(i,1,j) !{ wyao 20170714 ! wyao 20170714} stf(i,j,idfe) = sbc(i,jrow,idfeflx)*tmask(i,1,j) smf(i,j,1) = sbc(i,jrow,itaux)*umask(i,1,j) smf(i,j,2) = sbc(i,jrow,itauy)*umask(i,1,j) enddo enddo return end