subroutine gasbc (is, ie, js, je) !======================================================================= ! calculate boundary conditions for the atmospheric model !======================================================================= implicit none integer ie, is, je, js, i, iem1, isp1, j, jem1, jsp1, k, n real sss, sst, xconv, t_in, s_in, dic_in, ta_in, co2_in, pt_in real sit_in, atmpres, pHlo, pHhi, pH, co2star, dco2star, pCO2 real dpco2, CO3, Omega_c, Omega_a, scco2, piston_vel, avgflxc real calday, f, sco2, o2sat, o2sato, o2surf, piston_o2, cfc11ccn real cfc12ccn, wt, sccfc, piston_cfc, sol_cfc, cfcsat, ao, tarea real tdc14ccn, h_r, d, f1, f2, f3, f4, f5, area, C2K, tmp, zero real depth, hplus !======================= 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 $ ) parameter (nsrc=0 $ +1 $ +1 $ +1 $ +4 $ +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 & ) 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 "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 ! 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 ! { 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 ! { 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) !======================== include file "ice.h" ========================= ! arrays for the ice model ! hice = ice thickness (tau-1, tau, tau+1) ! aice = ice area (tau-1, tau, tau+1) ! tice = ice surface temperature ! hsno = effective thickness of snow ! psno = precipitation as snow (+ accumulation, - melt) ! frzpt = freezing temperature of sea water ! cpts ice model ("ice_cpts") ! A = area ! heff = effective thickness of ice ! Ts = ice surface temperature ! groice = thermo ice thickness change from diagnostic ! E = enthalpy of each ice layer ! strength = pressure (from Rothrock), must have ncat > 3 ! hseff = effective thickness of snow ! ice dynamics ("ice_evp") ! uice = u ice velocity ! vice = v ice velocity ! pice = pressure due to internal stress ! xint = x component of ice interaction ! yint = y component of ice interaction ! del = delta for ice dynamics models ! eI = divergence of ice velocity ! nseg = number of domains needed for ice calculations ! jsi = start of limited domain for ice calculations ! jei = end of limited domain for ice calculations ! brine convection ("convect_brine") ! cbf = flux from rejected brine over ice categories ! cba = area of brine rejection for each categories ! cba0 = area of the cell without brine rejection ! land ice ("landice_data") ! hicel = land ice thickness ! aicel = land ice area integer nseg, jsi, jei real hice, aice, tice, hsno, psno, frzpt real A, heff, Ts, groice, E, hseff, strength real uice, vice, pice, xint, yint, del, eI real cbf, cba, cba0, hicel, aicel common /ice_r/ hice(imt,jmt,2), aice(imt,jmt,2), tice(imt,jmt) common /ice_r/ hsno(imt,jmt,2), psno(imt,jmt), frzpt(imt,jmt) integer ncat parameter (ncat=1) common /ice_r/ uice(imt,jmt), vice(imt,jmt), pice(imt,jmt) common /ice_r/ xint(imt,jmt), yint(imt,jmt) common /ice_r/ del(imt,jmt), eI(imt,jmt) common /ice_i/ nseg, jsi(jmt), jei(jmt) !====================== include file "switch.h" ======================== ! all time dependent decisions are made by time manager "tmngr.F" ! and communicated elsewhere to the model via logical switches. ! inputs: (defaulted in "blkdta.F", optionally reset via namelist) ! runlen = integration period (see rununits). note "runlen" should ! be an integral number of density time steps. if not, ! then "runlen" is automatically adjusted to insure this. ! fractional days are supported but not fractional months ! or years. ! rununits= units of "runlen". may be "days", "months", or "years". ! tmngr will convert "runlen" which is in "rununits" ! to "rundays" in units of days. ! segtim = the integration time "runlen" is broken into a number of ! segments each of length "segtim" days. updated surface ! boundary conditions are applied to MOM every "segtim" ! days. this is useful when coupling to atmospheric models ! in which case both models exchange surface boundary ! conditions every "segtim" days where "segtim" ! is 1/(coupling frequency). without an atmospheric model, ! when getting surface boundary conditions from data, ! "segtim" is set to the time step (in days) by mom.F. in ! either case, "runlen" (in days) should be an integral ! number of "segtim". ! nmix = number of time steps between mixing timesteps. used ! to damp timestep splitting due to centred leapfrog. ! init = (true,false) indicates that this run is a ! (start from initial conditions, restart) ! initbgc = (true,false) indicates that this run is a ! (start from initial BGC conditions, restart) ! restrt = (true,false) = (do,don`t) write a restart at the end ! of the run ! eb = (true,false) configures for the use of a ! (euler backward,forward) type mixing timestep ! init_time = (true,false) sets restarts to initial time ! init_time_in = (true,false) sets input restart to initial time ! init_time_out = (true,false) sets output restart to initial time !----------------------------------------------------------------------- ! inputs to tmngr.F: diagnostic intervals !----------------------------------------------------------------------- ! note: switches are used to control the interval between doing ! diagnostics. units for all switches are in days. ! setting a switch < 0.0 disables whatever the switch is ! controlling. setting it = 0.0 causes the diagnostic to be ! done every time step, and setting it > 0.0 causes the ! diagnostic to be done repeatedly on the specified interval. ! cmixint = number of days between writing estimated mixing coeffs ! on faces of T cells and U cells ! crossint = number of days between writing diapycnal and isopycnal ! components of flow ! fctint = number of days between writing difference between ! FCT and leapfrog advection ! densityint = number of days between writing density ! exconvint = number of days between writing temperature rate of ! change due to explicit convection ! glenint = number of days between global energetics integrals. ! trmbint = number of days between momentum and tracer term ! balances (global and regional). ! itrmb = (true,false) = (do,don`t) write regional mask info for ! the term balance diagnostic. Typically set true ! at the beginning of a run; otherwise false since it is ! not necessary to keep writing a time independent field ! particularly when it may be a significant part of the ! time dependent part of the diagnostic. ! gyreint = number of days between calculation of tracer northward ! transport. ! igyre = (true,false) = (do,don`t) write regional mask info for ! the gyre diagnostic. Typically set true ! at the beginning of a run; otherwise false since it is ! not necessary to keep writing a time independent field ! particularly when it may be a significant part of the ! time dependent part of the diagnostic. ! vmsfint = number of days between calculation of vertical and ! meridional stream function. ! tyzint = number of days between calculation of zonally averaged ! tracer components. ! prxzint = number of days between printouts of x-z data. ! extint = number of days between printouts of external mode. ! dspint = number of days between surface pressure calculation. ! Note: only when "diagnostic_surface_height" is enabled. ! dspper = averaging period for "diagnostic_surface_height" ! tavgint = number of days between regional tracer averages (under ! horizontal regions). ! itavg = (true,false) = (do,don`t) write regional mask info for ! the tracer average diagnostic. Typically set true ! at the beginning of a run; otherwise false since it is ! not necessary to keep writing a time independent field ! particularly when it may be a significant part of the ! time dependent part of the diagnostic. ! tmbint = number of days over which tracer equation in averaged ! in depth and longitude to determine the meridional ! balance among storage, divergence, dissipation and ! forcing. ! tmbper = averaging period for "meridional_tracer_balance" ! itmb = (true,false) = (do,don`t) write "msktmb" for tracer ! the meridional balance diagnostic. Typically set true ! at the beginning of a run; otherwise false since it is ! not necessary to keep writing a time independent field ! particularly when it may be a significant part of the ! time dependent part of the diagnostic. ! tsiint = number of days between printing of time step integrals. ! tsiper = averaging period for "time_step_monitor" ! stabint = number of days between sampling for various stability ! criteria. ! timavgint= interval (days) for writing time mean data from ! the "averaging" grid (only when "time_averages" is ! enabled). if "timavgint" is not an integral number of ! density time steps,"timavgint" is automatically adjusted ! to insure this. if the number of days to integrate is ! not an integral number of "timavgint" then the last ! averaging period will be less than "timavgint" days.this ! may lead to one more averaging period than expected. ! see "iounit.h" for more details. ! timavgper= averaging period for "time_averages" ! xbtint = averaging period (days) for writing XBT data (only when ! "xbts" is enabled). if "xbtint" is not an integral ! number of density time steps, "xbtint" is automatically ! adjusted to insure this. if the number of days to ! integrate is not an integral number of "xbtint" then the ! last averaging period will be less than "xbtint" days. ! this may lead to one more averaging period than ! expected. see "iounit.h" for more details. ! xbtper = averaging period for "xbts" ! zmbcint = number of days between calculation of zonal mean ! surface boundary conditions (and related quantities) ! restint = number of days between saving restarts ! tbtint = averaging period (days) for writing term balances ! tbtper = averaging period for tracer term balances !----------------------------------------------------------------------- ! outputs from tmngr.F: logical switches !----------------------------------------------------------------------- ! rundays = integration time in days (from "runlen") ! the following are logical counterparts to the above switches are ! set within "tmngr" every time step. logical switches control all ! decisions about when to do things in MOM. ! cmixts = (false,true) = (don`t, do) do write estimated mixing ! coefficients on this time step. ! based on "cmixint". ! crossts = (false,true) = (don`t, do) write diapycnal and ! isopycnal components of flow on this time step. ! based on "crossint". ! fctts = (false,true) = (don`t, do) write difference between ! FCT and leapfrog advection on this time step. ! based on "fctint". ! densityts = (false,true) = (don`t, do) write density on this time ! step. based on "densityint". ! exconvts = (false,true) = (don`t, do) do write temperature change ! due to explicit convection on this time step. ! based on "exconvint". ! glents = (false,true) = (don`t, do) do calculation of global ! energy integrals on this time step. based on "glenint". ! trmbts = (false,true) = (don`t, do) do calculation of momentum & ! tracer term balance on this timestep. based on "trmbint" ! gyrets = (false,true) = (don`t, do) do calculation of tracer ! northward transport on this timestep. based on "gyreint" ! vmsfts = (false,true) = (don`t, do) do calculation of vertical ! and meridional stream function on this time step. ! based on "vmsfint" ! tyzts = (false,true) = (don`t, do) do calculation of zonally ! averaged tracer components on this time step. ! based on "tyzint" ! prxzts = (false,true) = (don`t, do) do printouts of x-z data ! on this time step. based on "prxzint" ! extts = (false,true) = (don`t, do) do printout of external mode ! on this time step. based on "extint" ! dspts = (false,true) = (don`t, do) do calculation of diagnostic ! surface pressure on this time step. based on "dspint" ! stabts = (false,true) = (don`t, do) test for stability on this ! time step. based on "stabint" ! tavgts = (false,true) = (don`t do) do tracer averages on this ! time step. based on "tavgint" ! tmbts = (false,true) = (don`t, do) write out tracer meridional . ! balance on this time step. based on "tmbint" ! tsits = (false,true) = (don`t, do) print time step integrals ! on this time step. based on "tsiint" ! zmbcts = (false,true) = (don`t, do) print zonal mean boundary ! conditions on this time step. based on "zmbcint" ! timavgts = (false,true) = (don`t, do) write time mean data ! on this time step. based on "timavgint" ! xbtts = (false,true) = (don`t, do) write averaged XBT data on ! this time step based on "xbtint" ! restts = (false,true) = (don`t, do) save a restart on this time ! step based on "restint" ! tbtts = (false,true) = (don`t, do) write averaged tracer term ! balance data on this time step based on "tbint" ! leapfrog= (false,true) on a (mixing, normal leapfrog) time step ! based on "nmix" ! euler1 = true on the 1st pass of an euler backward time step ! otherwise false. (applies when "eb" = true) ! euler2 = true on the 2nd pass of an euler backward time step ! otherwise false. (applies when "eb" = true) ! forward = true on a forward time step. otherwise false ! (applies when "eb" = false) ! the following logical switches are based on the model time step. ! first = (true,false) = when it`s (the first, not the first) ! time step of a run ! eots = end of a time step. always true except for first ! pass of an euler backward time step ! eorun = last time step of a run. always false except during the ! last time step of the run. ! eoday = true when within 1/2 time step of the end of a day ! else ... false ! eoweek = true when within 1/2 time step of the end of a 7 day ! week (referenced to the start of a year) else ...false ! eo2wks = true when within 1/2 time step of the end of two weeks ! (referenced to the start of a year) else ... false ! midmon = true when within 1/2 time step of the middle of a month ! else ... false ! eomon = true when within 1/2 time step of the end of a month ! else ... false ! eoyear = true when within 1/2 time step of the end of a year ! else ... false ! osegs = true on the 1st time step of an ocean segment in mom.F ! otherwise false. ! osege = true on the last time step of an ocean segment in mom.F ! otherwise false. character(8) :: rununits common /switc_c/ rununits integer nmix, ieoday,ieoweek,ieo2wks integer ieomon,imidmon,ieoyear,ieorun common /switc_i/ nmix, ieoday,ieoweek,ieo2wks common /switc_i/ ieomon,imidmon,ieoyear,ieorun logical eb, leapfrog, euler1, euler2, forward, eots logical init, first, restrt, itavg, itmb, itrmb, igyre logical eoday, eoweek, eo2wks, eomon, midmon, eoyear, eorun common /switc_l/ eb, leapfrog, euler1, euler2, forward, eots common /switc_l/ init, first, restrt common /switc_l/ itavg, itmb, itrmb, igyre common /switc_l/ eoday, eoweek, eo2wks common /switc_l/ eomon, midmon, eoyear, eorun logical init_time, init_time_in, init_time_out, initbgc common /switc_l/ init_time, init_time_in, init_time_out, initbgc real runlen, rundays common /switc_r/ runlen, rundays !----------------------------------------------------------------------- ! S W I T C H E S B A S E D O N A N I N T E R V A L ! each interval switch needs three variables in common. The ! following naming convention is used. ! 1) an interval (real) for diagnostic output (e.g,. glenint) ! 2) a switch (logical) for the interval (e.g., glents ) ! the third is an internal variable needed by the time manager ! to support calculation of the logical switch ! 3) an index (integer) (e.g., iglenint) ! the user must specify the interval [e.g., glenint] for diagnostic ! output in units of days. tmngr sets the corresponding logical ! switch [e.g., glents] every time step. It is set to true when ! within half a time step of the requested interval, otherwise it is ! false. All decisions relating to the interval [e.g., glenint] ! are based on the logical switch [e.g., glents]. ! internal time structures ! The switch index [e.g., iglenint] is used to subsrcipt into ! internal arrays maintained by tmngr.F. The switch index is ! allocated on the first call to function "alarm". ! The array entry [e.g., iinterval(iglenint)] is a time index to the ! internal representation of the interval [e.g., glenint]. ! The array entry [e.g., ialarm(iglenint)] is a time index to the ! next time the alarm will be true. !----------------------------------------------------------------------- integer itavgint, iglenint, itrmbint, iprxzint common /switc_i/ itavgint, iglenint, itrmbint, iprxzint logical tavgts, glents, trmbts, prxzts common /switc_l/ tavgts, glents, trmbts, prxzts real tavgint, glenint, trmbint, prxzint common /switc_r/ tavgint, glenint, trmbint, prxzint integer iextint, iexconvint, icmixint common /switc_i/ iextint, iexconvint, icmixint logical extts, exconvts, cmixts common /switc_l/ extts, exconvts, cmixts real extint, exconvint, cmixint common /switc_r/ extint, exconvint, cmixint integer ivmsfint, igyreint, ityzint, ifctint common /switc_i/ ivmsfint, igyreint, ityzint, ifctint logical vmsfts, gyrets, tyzts, fctts common /switc_l/ vmsfts, gyrets, tyzts, fctts real vmsfint, gyreint, tyzint, fctint common /switc_r/ vmsfint, gyreint, tyzint, fctint integer istabint, izmbcint, icrossint, idensityint common /switc_i/ istabint, izmbcint, icrossint, idensityint logical stabts, zmbcts, crossts, densityts common /switc_l/ stabts, zmbcts, crossts, densityts real stabint, zmbcint, crossint, densityint common /switc_r/ stabint, zmbcint, crossint, densityint integer iosegs, iosege common /switc_i/ iosegs, iosege logical osegs, osege common /switc_l/ osegs, osege real segtim common /switc_r/ segtim integer irestint common /switc_i/ irestint logical restts common /switcl/ restts real restint common /switcr/ restint !----------------------------------------------------------------------- ! S W I T C H E S B A S E D O N A N I N T E R V A L ! A N D A V E R A G I N G P E R I O D ! each averaging period switch needs five variables in common. The ! following naming convention is used. ! 1) an interval (real) for diagnostic output (e.g. xbtint ) ! 2) a switch (logical) for the interval (e.g. xbtts ) ! 3) an averaging period (real) (e.g. xbtper ) ! 4) a switch (logical) for accumulating (e.g. xbtperts) ! the third is an internal variable needed by the time manager ! to support calculation of the logical switches ! 5) an index (integer) (e.g. ixbtint ) ! The user must specify the interval [e.g., xbtint] for diagnostic ! output in units of days and the averaging period [e.g., xbtper] ! in units of days. The averaging period may be less than or equal ! to the interval. For example, if the interval is 30.0 days and the ! averaging period is 5.0 days, results will be averaged over all ! time steps within days 26, 27, 28, 29, and 30. An averaging period ! of 0.0 days averages over the last time step of the interval (as ! does xbtper = dt), and an averaging period less than zero turns ! the switches off for all time steps. ! The logical switch for writing output at the specified interval ! [e.g., xbtts] is set to true on the last time step of the ! averaging period. The logical switch for accumulating results ! [e.g., xbtperts] is true for all time steps within the averaging ! period, otherwise it is false. ! internal time structures ! The index [e.g., ixbtint] is allocated on the first call to ! function "avg_alarm". The array element iperiod(ixbtint) is an ! index to the time structure for the internal representation of ! "xbtper", and ilastsw(ixbtint) is the index of the switch that ! flags the last time step of the accumulation period. ! Depending on use, ilastsw(ixbtint) may either be the index ! of another "named" switch or the index of a new switch ! allocated on the first time step. ! In the latter case, iinterval(ilastsw(ixbtint)) is the index of ! the time structure where "xbtint" is stored in internal form, ! and ialarm(ilastsw(ixbtint)) is the index of the time when an ! accumulation period will next end. ! The variable nextts(ixbtint) is true whenever the next ! time step will begin the accumulation period. !----------------------------------------------------------------------- integer ixbtint, idspint, itmbint, itimavgint common /switc_i/ ixbtint, idspint, itmbint, itimavgint logical xbtts, dspts, tmbts, timavgts logical xbtperts, dspperts, tmbperts, timavgperts common /switc_l/ xbtts, dspts, tmbts, timavgts common /switc_l/ xbtperts, dspperts, tmbperts, timavgperts real xbtint, dspint, tmbint, timavgint real xbtper, dspper, tmbper, timavgper common /switc_r/ xbtint, dspint, tmbint, timavgint common /switc_r/ xbtper, dspper, tmbper, timavgper integer itbtint, itsiint common /switc_i/ itbtint, itsiint logical tbtts, tsits logical tbtperts, tsiperts common /switc_l/ tbtts, tsits common /switc_l/ tbtperts, tsiperts real tbtint, tsiint real tbtper, tsiper common /switc_r/ tbtint, tsiint common /switc_r/ tbtper, tsiper !----------------------------------------------------------------------- ! S W I T C H E S B A S E D O N ! C A L E N D A R O R P R E V I O U S S W I T C H ! A N D A V E R A G I N G P E R I O D ! the following logical switches are based on any calendar or ! interval switch and an averaging period (in days). The averaging ! period must be less than or equal to the interval. The last ! time step of the averaging period is at the end of the interval. ! If the averaging period is set to zero, the averaging period ! consists only of the last time period of the interval. If ! the averaging period is less than zero, these switches are always ! false. ! each averaging period switch needs four variables in common. For ! example, if the averaging period is before the end of each month ! then the calendar switch (eomon), and index (ieomon) are presumed ! to exist in common and need not be added. ! Additionally, four items are needed. ! 1) an averaging period (real) (e.g. testper ) ! 2) a switch (logical) for accumulating results (e.g. testperts) ! 3) a switch (logical) for the end of interval (e.g. testts ) ! the fourth is an internal variable needed by the time manager ! to support calculation of the logical switch ! 4) an index (integer) (e.g. itestper ) ! Suppose it is required to produce averages over all time steps ! during the last 5 days of each month. Then "testper" = 5.0 and ! the following will calculate the accumulating switch. ! testts = avg_alarm(itestper, ihalfstep, 0, testper, iref, ieomon) ! testperts = on(itestper) ! Note the use of "ieomon" to key off the months. The switch ! "testts" will be true whenever "eomon" is true. ! Also note that when an averaging switch is keyed off another ! switch, the switch interval argument is not used, but is ! retained for consistency with the form of other averaging ! switches. !----------------------------------------------------------------------- integer maxsw parameter (maxsw=100) integer itestper, nsw, ialarm, iinterval, iperiod, ilastsw common /switc_i/ itestper, nsw, ialarm(maxsw), iinterval(maxsw) common /switc_i/ iperiod(maxsw), ilastsw(maxsw) logical testts, testperts, cplts, tsicpl, tavgcpl, on, lastts logical nextts common /switc_l/ testts, testperts, cplts, tsicpl, tavgcpl common /switc_l/ on(maxsw), lastts(maxsw), nextts(maxsw) real testint, testper common /switc_r/ testint, testper !======================= include file "tmngr.h" ======================== ! time manager variables !----------------------------------------------------------------------- ! time manager inputs: !----------------------------------------------------------------------- ! how to choose a reference time: ! refrun = (true,false) to base calculation for diagnostic switches ! on (the start of each job, other reference time) ! example: ! suppose each job submission integrates ! for one month but the number of days per month changes. ! setting "refrun" = true and setting ! "timavgint" = (days in month)/3 will give 3 averaging ! periods per month of approximately 10 days each. the ! only restriction is that "timavgint"is an integral number ! of time steps (if not then "timavgint" is reset to insure ! this condition. other diagnostic switches do not have ! this restriction). ! refinit = (true, false) for basing calculation of logical switches ! on (initial conditions, other reference time) ! example: if term balances are desired every 20 days ! (trmbint=20.0) and refinit = true, then they ! will be done every 20 days starting from initial ! condition time. ! refuser = (true, false) to base calculations of logical switches ! on (user-chosen reference time, other reference time) ! if refuser = true, the user must also supply values for ! ryear, rmonth, rday, rhour, rmin, rsec (integer) ! example: if term balances are desired every 20 days ! (trmbint=20.0) and refuser = true, then they will be done ! every 20 days counting from reference time, ignoring the ! initial condition time. for comparing diagnostics from ! various experiments with different initial condition ! times, refuser = true will be more appropriate. setting ! refuser = true and choosing the reference time to be ! the initial condition time is the same as refinit = true. ! summary of how to choose the time for referencing calculations ! of logical switches ! refrun = T ==> referenced to the start of each run ! refinit = T ==> referenced to initial condition time given by: ! year0, month0, day0, hour0, min0, sec0 ! refuser = T ==> referenced to user specified reference time so ! must set: ryear, rmonth, rday, rhour, rmin, rsec !----------------------------------------------------------------------- ! time variable arrays ! arrays "iday" and "msday" contain the primary internal ! representation of all times within the time manager. they are ! referenced by using a subscript to indicate which time. ! iday = integer days (since Dec 31, 1899 when specifying a date) ! msday = non-negative integer milliseconds after midnight ! it is desirable to have time information expanded to include the ! following secondary time fields: ! year = ! month = ! day = ! hour = ! minute = ! second = ! tstamp = 32 character date and time stamp m/d/y h:m:s ! dayofyear = integer day of the year (1..yrlen) ! dayofweek = 1=sun - 7=sat ! daysinmon = days in the month ! daysinyear = days in the year ! those times for which primary and secondary information is ! maintained by the time manager are called "full times". those for ! which only primary information is kept are called "short times" ! indices to "full times" (including year, month ,day, etc). ! itime = simulation time corresponding to "itt" ! initial = time of the initial conditions ! irunstart = time of the start of the run ! iuser = user defined reference time ! iref = one of the three above selected by logicals ! (refinit, refrun, refuser) ! indices to "short times". ("iday", "msday" only) ! isunday = time of a sunday for week and two week switches ! ihalfstep = dt/2 beyond itime ! imodeltime = time since initial conditions ! iruntime = time since run start ! iusertime = time since user specified reference time ! idt = integer days and milliseconds of dt ! idtd2 = integer days and milliseconds of dt/2 ! ireftime = time used locally in alarm function ! for any time index (short or full) the internal representation ! may be converted to either real days or real seconds using ! the functions: ! realdays(index) ! realsecs(index) ! dayoyr = relative day number referenced to the beginning ! of the current year. (real) ! relyr = number of years (and fractional years) of model ! integration (for time tau+1 {itt}) relative to ! initial condition ! prelyr = relyr for previous time step ! stamp = 32 character date and time for current model timestep ! pstamp = 32 character date and time for previous model timestep ! calendar = calendar name ! itt = current time step counter (from initial cond.) ! itt0 = time step at start of current run ! variables used for initialization ! irstdy = integer number of days at start of run ! msrsdy = fractional day in millisec at start of run ! year0 = year of initial conditions ! month0 = month of initial conditions ! day0 = day of initial conditions ! hour0 = hour of initial conditions ! min0 = minute of initial conditions ! sec0 = second of initial conditions ! ryear = year of user specified reference time ! rmonth = month of user specified reference time ! rday = day of user specified reference time ! rhour = hour of user specified reference time ! rmin = minute of user specified reference time ! rsec = second of user specified reference time !----------------------------------------------------------------------- integer ntimes, nfulltimes parameter (ntimes = 100, nfulltimes = 30) character(32) :: tstamp, stamp, pstamp character(120) :: calendar common /tmngrc/ tstamp(nfulltimes), stamp, pstamp, calendar integer nextfulltime, nexttime integer initial, iref, irunstart, itime, iuser integer iruntime, imodeltime, ireftime, iusertime integer ihalfstep, isunday integer itemptime,itemptime2,itmptime,itmptime2,itmptime3 integer idt, idtd2 integer iday, msday, year, month, day, hour, minute, second integer dayofyear, dayofweek, daysinmon, daysinyear integer itt0, itt, irstdy, msrsdy integer year0, month0, day0, hour0, min0, sec0 integer ryear, rmonth, rday, rhour, rmin, rsec common /tmngri_i/ nextfulltime, nexttime common /tmngri_i/ initial, iref, irunstart, itime, iuser common /tmngri_i/ iruntime, imodeltime, ireftime, iusertime common /tmngri_i/ ihalfstep, isunday common /tmngri_i/ itemptime, itemptime2, itmptime, itmptime2 common /tmngri_i/ itmptime3 common /tmngri_i/ idt, idtd2 common /tmngri_i/ iday(ntimes), msday(ntimes) common /tmngri_i/ year(nfulltimes), month(nfulltimes) common /tmngri_i/ day(nfulltimes), hour(nfulltimes) common /tmngri_i/ minute(nfulltimes), second(nfulltimes) common /tmngri_i/ dayofyear(nfulltimes), dayofweek(nfulltimes) common /tmngri_i/ daysinmon(nfulltimes), daysinyear(nfulltimes) common /tmngri_i/ itt0, itt, irstdy, msrsdy common /tmngri_i/ year0, month0, day0, hour0, min0, sec0 common /tmngri_i/ ryear, rmonth, rday, rhour, rmin, rsec logical refrun, refinit, refuser common /tmngr_l/ refrun, refinit, refuser real dayoyr, relyr, prelyr common /tmngr_r/ dayoyr, relyr, prelyr !----------------------------------------------------------------------- ! acceleration parameters for long term forcing !----------------------------------------------------------------------- ! an acceleration factor can be used to accelerate long term ! forcing. forcing and output at time intervals of less than one ! year are not affected (eg. seasonal cycle of insolation). ! averaging intervals are also not changed. only longer term ! forcings are accelerated (eg. CO2 or orbital changes). ! runlen, runstep and most output intervals will will be divided by ! this factor (accel). output intervals are divided by accel if the ! interval divided by accel is larger or equal to a year. both ! runlen and runstep divided accel must still be evenly divisible ! by the coupling time (segtim). output intervals divided by accel ! should still be evenly divisible by the coupling time and the ! averaging period. ! times in various output files may appear to be inconsistent. ! for example the year may switch part way through a seasonal cycle. ! be cautious when looking at output averaged over less than a year. ! it is probably best to use an integer value for accel that will ! leave the runlen and runstep as integer multiples of a year. the ! initial year can be set to offset any existing time in a restart. ! choosing an inappropriate acceleration factor may mess up output ! intervals and averaging periods. think before you use. ! accel = factor for accelerating forcing transient data ! accel_yr0 = relyr to start acceleration real accel, accel_yr0 common /tmngr_r/ accel, accel_yr0 !======================= include file "cembm.h" ======================== ! parameters for use in the energy balance model (also see atm.h) ! namix = time steps between mixing (set in atmos.in) ! lf = time step flag (1=>leapfrog, 2=>forward) ! niats = number of ice advection sub time steps ! nivts = time steps between recalculating ice velocities ! nivc = time step counter for nivts ! ns = number of subcycles for advection and diffusion ! pyear = default paleo calendar year (-/+ = BC/AD) ! dts = time step (2*dtatm=>leapfrog, dtatm=>forward) ! co2ccn = atmospheric CO2 concentration (ppmv) ! co2emit = atmospheric CO2 emissions flux (g cm-2 s-1) ! co2emit_fuel = emissions flux from fossil fuels (g cm-2 s-1) ! co2emit_land = emissions flux from land use change (g cm-2 s-1) ! anthro = radiative forcing by atmospheric CO2 ! co2for = atmospheric CO2 forcing term (in units of heat flux) ! c14ccn = atmospheric C14 concentration (ppmv) ! dc14ccn = atmospheric dC14 concentration (permil) ! dc14ccnn = northern hemisphere atmospheric dC14 concentration ! dc14ccne = equatorial atmospheric dC14 concentration ! dc14ccns = southern hemisphere atmospheric dC14 concentration ! c14prod = atmospheric C14 production (umol cm-2 s-1) ! cfc11ccnn = northern hemisphere atmospheric CFC11 concentration ! cfc11ccns = southern hemisphere atmospheric CFC11 concentration ! cfc12ccnn = northern hemisphere atmospheric CFC12 concentration ! cfc12ccns = southern hemisphere atmospheric CFC12 concentration ! scatter = proportion of solar scattered by the atmosphere ! solarconst = solar constant (g s-3) ! cssh = constant used in calculation of ssh (g g-1) ! cdatm = drag coefficient (dimensionless) ! cpatm = atmospheric heat capacity (cm**2 s-2 K-1) ! sht = scale height for temperature ! shq = scale height for specific humidity ! shc = scale height for carbon ! rhoatm = density of air at sea surface (g cm-3) ! esatm = atmosphere emissivity times Stefan's constant ! rhoocn = representative sea surface density ! esocn = ocean emissivity times Stefan's constant ! vlocn = latent heat of vapourization of water ! cdice = drag coefficient (dimensionless) ! dampice = time scale for freezing first layer under ice (days) ! rhoice = ice density (g cm-3) ! rhosno = snow density (g cm-3) ! esice = ice emissivity times Stefan's constant ! slice = latent heat of sublimation of ice ! flice = latent heat of fusion of ice (cm2 s-2) ! condice = ice conductivity (g*cm s-3 K-1) ! tsno = air temperature for accumulating snow ! hsno_max = maximum snow depth ! totaltime = total time for long term averages ! rlapse = lapse rate (K cm-1) ! soilmax = soil water field capacity (cm) ! eslnd = land emissivity time Stefan's constant ! pass = atmospheric transmission coefficient (%/100) ! ice_calb = ice coalbedo (%/100) ! sno_calb = snow coalbedo (%/100) ! pcfactor = precip - cloud correlation factor ( %/100) ! rf1 = factor used in calculating lapse rate reduction ! rf2 = factor used in calculating lapse rate reduction ! co2_yr = co2 forcing year (-/+ = BC/AD) ! agric_yr = agricutural forcing year (-/+ = BC/AD) ! landice_yr = ice sheet forcing year (-/+ = BC/AD) ! solar_yr = solar forcing year (-/+ = BC/AD) ! orbit_yr = orbital forcing year (-/+ = BC/AD) ! volcano_yr = volcanic forcing year (-/+ = BC/AD) ! sulph_yr = sulphate forcing year (-/+ = BC/AD) ! aggfor_yr = additional greenhouse gas forcing year (-/+ = BC/AD) ! cfcs_yr = cfc forcing year (-/+ = BC/AD) ! c14_yr = c14 forcing year (-/+ = BC/AD) ! ice_yr = redundant (same as landice_yr, kept for compatibility) ! dalt_v = dalton number over vegetation (no vegetation model) ! dalt_o = dalton number over ocean ! dalt_i = dalton number over ice ! rhmax = maximum relative humidity ! volcfor = anomalous volcanic forcing (g/s**3) ! aggfor = anomalous additional greenhouse gas forcing (g s-3) ! aggfor_os = additional greenhouse gas forcing offset (g s-3) ! atmsa = atmospheric surface area (cm 2) ! ocnsa = exposed ocean surface area (cm 2) ! sealev = anomalous sea level change (cm) ! dsealev = anomalous sea level change in segtim (cm) ! sealev_yr = sea level forcing year (-/+ = BC/AD) ! itrack_co2 = index for averaging co2 array ! ntrack_co2 = number of points for average co2 ! itrack_sat = index for averaging sat array ! ntrack_sat = number of points for average sat ! vcsref = climate sensitivity reference temperature (C) ! vcsfac = climate sensitivity factor (mW/m2/C) ! vcsyri = year to start changing the climate sensitivity (year) ! gtoppm = conversion from g cm-2 to ppmv ! carbemit = accumulated co2 emissions (Pg) ! adiff = anomolous diffusion factor (percent/100/C) ! dtbar = global average sat anomoly (C) ! cropf = scaler for crop area ! pastf = scaler for pasture area integer namix, lf, niats, nivts, nivc, ns, itrack_co2 integer ntrack_co2, itrack_sat, ntrack_sat common /cembm_i/ namix, lf, niats, nivts, nivc, ns, itrack_co2 common /cembm_i/ ntrack_co2, itrack_sat, ntrack_sat real pyear, dts, co2ccn, co2emit, co2emit_fuel, co2emit_land real anthro, co2for, c14ccn, dc14ccn, dc14ccnn, dc14ccne, dc14ccns real c14prod, cfc11ccnn, cfc11ccns, cfc12ccnn, cfc12ccns, scatter real solarconst, cssh, cdatm, cpatm, sht, shq, shc, rhoatm, esatm real rhoocn, esocn, vlocn, cdice, dampice, rhoice, rhosno, esice real slice, flice, condice, tsno, hsno_max, totaltime, rlapse real soilmax, eslnd, pass, ice_calb, sno_calb, pcfactor, rf1, rf2 real co2_yr, agric_yr, landice_yr, solar_yr, orbit_yr, volcano_yr real sulph_yr, aggfor_yr, cfcs_yr, c14_yr, ice_yr, dalt_v, dalt_o real dalt_i, rhmax, volcfor, aggfor, aggfor_os, atmsa, ocnsa real sealev, dsealev, sealev_yr, vcsref, vcsfac, vcsyri, gtoppm real carbemit, adiff, dtbar, cropf, pastf, anthro_srm common /cembm_r/ pyear, dts, co2ccn, co2emit, co2emit_fuel common /cembm_r/ co2emit_land, anthro, co2for, c14ccn, dc14ccn common /cembm_r/ dc14ccnn, dc14ccne, dc14ccns, c14prod, cfc11ccnn common /cembm_r/ cfc11ccns, cfc12ccnn, cfc12ccns, cssh, cdatm common /cembm_r/ cpatm, sht, shq, shc, rhoatm, esatm, rhoocn common /cembm_r/ scatter, solarconst, esocn, dampice, rhoice common /cembm_r/ rhosno, esice, slice, flice, condice, vlocn common /cembm_r/ cdice, tsno, hsno_max, totaltime, rlapse, soilmax common /cembm_r/ eslnd, pass, ice_calb, sno_calb, pcfactor, rf1 common /cembm_r/ rf2, co2_yr, agric_yr, landice_yr, solar_yr common /cembm_r/ orbit_yr, volcano_yr, sulph_yr, aggfor_yr common /cembm_r/ cfcs_yr, c14_yr, ice_yr, dalt_v, dalt_o, dalt_i common /cembm_r/ rhmax, volcfor, aggfor, aggfor_os, atmsa, ocnsa common /cembm_r/ sealev, dsealev, sealev_yr, vcsref, vcsfac common /cembm_r/ vcsyri, gtoppm, carbemit, adiff, dtbar, cropf common /cembm_r/ pastf, anthro_srm ! ntatsa = time step counter for time averaging ! ntatia = number of time averaged time step integrals ! tai_sat = average integrated elevated surface air temperature ! tai_shum = average integrated surface specific humidity ! tai_precip = average integrated precipitation ! tai_evap = average integrated evaporation ! tai_ohice = total integrated sea ice volume ! tai_oaice = total integrated sea ice area ! tai_hsno = total integrated snow volume ! tai_lhice = total integrated land ice volume ! tai_laice = total integrated land ice area ! tai_co2ccn = average integrated CO2 concentration ! tai_co2emit = average integrated CO2 emissions ! tai_dc14ccn = average integrated dC14 concentration ! tai_cfc11ccn = average integrated CFC11 concentration ! tai_cfc12ccn = average integrated CFC12 concentration ! tai_maxit = average maximum iterations for atmospheric solver ! tai_nsat = average northern hemisphere surface air temperature ! tai_ssat = average southern hemisphere surface air temperature ! tai_nshum = average northern hemisphere surface specific humidity ! tai_sshum = average southern hemisphere surface specific humidity ! tai_nprecip = average northern hemisphere precipitation ! tai_sprecip = average southern hemisphere precipitation ! tai_nevap = average northern hemisphere evaporation ! tai_sevap = average southern hemisphere evaporation ! tai_nohice = total northern hemisphere sea ice volume ! tai_sohice = total southern hemisphere sea ice volume ! tai_noaice = total northern hemisphere sea ice area ! tai_soaice = total southern hemisphere sea ice area ! tai_nhsno = total northern hemisphere snow volume ! tai_shsno = total southern hemisphere snow volume ! tai_nlhice = total northern hemisphere land ice volume ! tai_slhice = total southern hemisphere land ice volume ! tai_nlaice = total northern hemisphere land ice area ! tai_slaice = total southern hemisphere land ice area ! tai_lsat = average sat over land ! tai_osat = average sat over ocean ! tai_lprecip = average precip over land ! tai_oprecip = average precip over ocean ! tai_levap = average evap over land ! tai_oevap = average evap over ocean ! tai_solins = average incoming solar ! tai_upsens = average surface sensible heat ! tai_uplwr = average surface upward longwave radiation ! tai_outlwr = average outgoing longwave radiation ! tai_dnswr = average downward (absorbed) shortwave ! tai_absswr = net absorbed shortwave radiation ! tai_netrad = net radiation at the top of the atmosphere ! tai_palb = average planetary albedo ! tai_aalb = average atmospheric albedo ! tai_salb = average surface albedo ! tai_lsalb = average surface albedo over land ! tai_osalb = average surface albedo over ocean ! tai_sst = average sea surface temperature ! tai_sss = average sea surface salinity ! tai_ssdic = average sea surface dic ! tai_ssc14 = average sea surface c14 ! tai_ssalk = average sea surface alkalinity ! tai_sso2 = average sea surface o2 ! tai_sspo4 = average sea surface po4 ! tai_ssno3 = average sea surface no3 !{ wyao 20170703 ! wyao 20170703} ! tai_ssdfe = average sea surface dfe ! tai_sscfc11 = average sea surface cfc11 ! tai_sscfc12 = average sea surface cfc12 ! tai_sulph = average tropospheric sulphate forcing ! tai_volc = average volcanic forcing ! tai_agg = average additional greenhouse gas forcing ! tai_catm = average total carbon in atmosphere ! tai_carbemit = average taccumulated co2 emissions integer ntatsa, ntatia common /cembm_i/ ntatsa, ntatia real tai_sat, tai_shum, tai_precip, tai_evap, tai_ohice, tai_oaice real tai_hsno, tai_lhice, tai_laice, tai_co2ccn, tai_co2emit real tai_dc14ccn, tai_cfc11ccn, tai_cfc12ccn, tai_maxit, tai_nsat real tai_ssat, tai_nshum, tai_sshum, tai_nprecip, tai_sprecip real tai_nevap, tai_sevap, tai_nohice, tai_sohice, tai_noaice real tai_soaice, tai_nhsno, tai_shsno, tai_nlhice, tai_slhice real tai_nlaice, tai_slaice, tai_lsat, tai_osat, tai_lprecip real tai_oprecip, tai_levap, tai_oevap, tai_solins, tai_upsens real tai_uplwr, tai_outlwr, tai_dnswr, tai_absswr, tai_netrad real tai_palb, tai_aalb, tai_salb, tai_lsalb, tai_osalb, tai_sst real tai_sss, tai_ssdic, tai_ssc14, tai_ssalk, tai_sso2, tai_sspo4 real tai_ssno3, tai_sscfc11, tai_sscfc12, tai_sulph, tai_volc real tai_agg, tai_catm, tai_ssdfe, tai_carbemit !{ wyao 20170703 ! wyao 20170703} common /cembm_r/ tai_sat, tai_shum, tai_precip, tai_evap common /cembm_r/ tai_ohice, tai_oaice, tai_hsno, tai_lhice common /cembm_r/ tai_laice, tai_co2ccn, tai_co2emit, tai_dc14ccn common /cembm_r/ tai_cfc11ccn, tai_cfc12ccn, tai_maxit, tai_nsat common /cembm_r/ tai_ssat, tai_nshum, tai_sshum, tai_nprecip common /cembm_r/ tai_sprecip, tai_nevap, tai_sevap, tai_nohice common /cembm_r/ tai_sohice, tai_noaice, tai_soaice, tai_nhsno common /cembm_r/ tai_shsno, tai_nlhice, tai_slhice, tai_nlaice common /cembm_r/ tai_slaice, tai_lsat, tai_osat, tai_lprecip common /cembm_r/ tai_oprecip, tai_levap, tai_oevap, tai_solins common /cembm_r/ tai_upsens, tai_uplwr, tai_outlwr, tai_dnswr common /cembm_r/ tai_absswr, tai_netrad, tai_palb, tai_aalb common /cembm_r/ tai_salb, tai_lsalb, tai_osalb, tai_sst common /cembm_r/ tai_sss, tai_ssdic, tai_ssc14, tai_ssalk common /cembm_r/ tai_sso2, tai_sspo4, tai_ssno3, tai_sscfc11 common /cembm_r/ tai_sscfc12, tai_sulph, tai_volc, tai_agg common /cembm_r/ tai_catm, tai_ssdfe, tai_carbemit !{ wyao 20170703 ! wyao 20170703} !======================== include file "atm.h" ========================= ! arrays for the energy-moisture balance model ! note: units for heat flux are: g s-3 (or mW m-2) ! units for fresh water flux are in g cm-2 s-1 ! downward is into the top surface (ocean, ice or land) ! upward is into the bottom of the atmosphere ! outward is out of the top of the atmosphere ! inward is into the top of the atmosphere ! at = tracers (previous and most recent) ! elev = elevations (cm) ! flux = flux accumulator (1=heat, 2=water, 3=carbon, ...) ! rh = relative humidity (%/100) ! tmsk = tracer grid ocean mask (0 = land, 1 = ocean) ! umsk = velocity grid ocean mask (0 = land, 1 = ocean) ! dn = northward atmospheric tracer diffusivity ! de = eastward atmospheric tracer diffusivity ! fcor = Coriolis factor ! heat fluxes ! solins = incoming short wave radiation (g s-3) ! dnswr = downward surface shortwave flux absorbed (g s-3) ! uplwr = upward surface longwave flux (g s-3) ! upsens = upward surface sensible heat flux (g s-3) ! upltnt = upward surface latent heat flux (g s-3) ! outlwr = outgoing atmosphere longwave flux (g s-3) ! fresh water fluxes ! precip = precipitation (g cm-2 s-1) ! evap = evaporation (g cm-2 s-1) ! disch = discharge from continents (g cm-2 s-1) ! vflux = normalized virtual flux ! land model ! soilm = soil moisture, as depth in bucket (g cm -2) ! runoff = water runoff from continents (g cm-2 s-1) ! surf = land surface temperature (C) ! annual average temperature ("embm_awind" and "embm_adiff") ! rtbar = running average atmospheric temperature (C) ! atbar = average temperature accumulator (C) ! tbar = climatological average atmospheric temperature (C) ! wind feedback ("embm_awind") ! awx = anomalous x component of wind (cm s-1) ! awy = anomalous y component of wind (cm s-1) ! apress = anomalous pressure (g cm-2 s-2) ! sulphates ("sulphate_data_transient") ! sulph = sulphate forcing ! distributed CO2 emissions ("co2emit_data" and "carbon_co2_2d") ! co2dist = distributed co2 forcing ! sea level ("sealev" or "sealev_data") ! elev_sealev = elevation anomaly due to sea level changes ! flux adjustment ("save_flxadj") ! flxadj = flux adjustment ! co2 emissions from tracking co2 ("co2emit_track_co2") ! track_co2 = array of global average co2 values ! co2 emissions from tracking sat("co2emit_track_sat" or "embm_vcs") ! track_sat = array of global average sat values ! indicies for atmospheric tracer array ! isat = index for surface air temperature ! ishum = index for surface specific humidity ! ico2 = index for atmospheric co2 ! mapat = map for atmospheric tracer names ! mtrack = maximum size for tracking arrays character(10) :: mapat(nat) common /embm_c/ mapat integer isat, ishum, ico2, mtrack parameter (mtrack=3650) common /embm_i/ isat, ishum, ico2 real at, elev, flux, rh, tmsk, umsk, dn, de, fcor, solins, dnswr real uplwr, upsens, upltnt, outlwr, precip, evap, disch, vflux real soilm, runoff, surf, rtbar, atbar, tbar, awx, awy, apress real sulph, co2dist, elev_sealev, flxadj, track_co2, track_sat common /embm_r/ at(imt,jmt,2,nat), elev(imt,jmt) common /embm_r/ flux(imt,jmt,nat+2), rh(imt,jmt) common /embm_r/ tmsk(imt,jmt), umsk(imt,jmt), dn(imt,jmt,nat) common /embm_r/ de(imt,jmt,nat), fcor(imt,jmt), solins(imt,jmt) common /embm_r/ dnswr(imt,jmt), uplwr(imt,jmt), upsens(imt,jmt) common /embm_r/ upltnt(imt,jmt), outlwr(imt,jmt), precip(imt,jmt) common /embm_r/ evap(imt,jmt), disch(imt,jmt), vflux(imt,jmt) common /embm_r/ soilm(imt,jmt,2), runoff(imt,jmt) common /embm_r/ surf(imt,jmt) common /embm_r/ rtbar(imt,jmt), atbar(imt,jmt), tbar(imt,jmt) !===================== include file "insolation.h" ===================== ! parameters used in calculating solar insolation ! eccen = orbital eccentricity ! obliq = obliquity (degrees) ! mvelp = moving vernal equinox longitude of perihelion (degrees) ! lambm0 = Mean long of perihelion at vernal equinox (radians) ! sindec = sine of solar declination angle in rad ! eccf = Earth-sun distance factor (ie. (1/r)**2) real eccen, obliq , mvelp, lambm0, sindec, eccf common /insolation_r/ eccen, obliq , mvelp, lambm0, sindec, eccf !====================== include file "calendar.h" ====================== ! calendar specification arrays !----------------------------------------------------------------------- ! eqyear = true to select a calendar in which each year ! has the same number of days (i.e., no leap years) ! false selects a Julian calendar ! eqmon = true to force all months to have the same number of days ! false => the usual 31, 28, 31, 30, ..., days per month. ! only used when eqyear = true ! dayname = character names of days ! monname = character names of months ! monlen = the length of each month (in days) when eqmon is true ! yrlen = the length of a typical (non-leap) year in days ! daypm = array of month lengths in days (non-leap) ! msum = array of cumulative days preceding each month ! (again, non-leap) ! daylen = the length of a day in seconds !----------------------------------------------------------------------- character(10) :: dayname character(12) :: monname common /calen_c/ dayname(7), monname(12) logical eqyear, eqmon common /calen_l/ eqyear, eqmon integer daypm, msum, yrlen, monlen common /calen_i/ daypm(12), msum(12), yrlen, monlen real daylen common /calen_r/ daylen !====================== 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 "solve.h" ======================== ! variables needed for solving atmospheric advection and diffusion ! newcoef = logical flag for calculating new coefficients logical newcoef(2,nat) common /solver_l/ newcoef integer iimtm2, jjmtm2, nord, nelm parameter (iimtm2 = imtm2) parameter (jjmtm2 = jmtm2) parameter (nord = iimtm2*jjmtm2) ! itin = requested maximum iterations ! itout = actual iterations ! newcoef = logical flag for calculating new coefficients ! bv = right hand side vector (b) ! xv = left hand side vector (x) ! epsin = requested maximum error ! epsout = actual error integer itin(nat), itout(nat) common /solver_i/ itin, itout ! for mgrid routine storage is by compass coefficient ! ap, an, as, ae, aw = centre, north, south, east, and west coef ! ie. ap*xp = an*xn + as*xs + ae*xe + aw*xw + bp integer levelin, levelout common /solver_i/ levelin, levelout real(kind=8) bv, xv, epsin, epsout, an, as, ae, aw, ap common /solver_r/ bv(nord), xv(nord), epsin(nat), epsout(nat) common /solver_r/ an(nord,2,nat), as(nord,2,nat), ae(nord,2,nat) common /solver_r/ aw(nord,2,nat), ap(nord,2,nat) ! grid terms for the atmospheric solver real dwgrd, degrd, azgrd, dsgrd, dngrd, asgrd, angrd common /solve_r/ dwgrd(2:imtm1), degrd(2:imtm1), azgrd(2:imtm1) common /solve_r/ dsgrd(2:jmtm1), dngrd(2:jmtm1), asgrd(2:jmtm1) common /solve_r/ angrd(2:jmtm1) !====================== include file "diaga.h" ========================= ! variables used for computing diagnostics: ! totalk = total number of levels involved in convection ! vdepth = ventilation depth (cm) ! pe = potential energy lost due to explicit convection (g/s2) real totalk, vdepth, pe common /cdiaga_r/ totalk(imt,jsmw:jemw), vdepth(imt,jsmw:jemw) common /cdiaga_r/ pe(imt,jsmw:jemw) ! sspH = sea surface pH ! ssCO3 = sea surface CO3 concentration ! ssOc = sea surface Omega_c ! ssOa = sea surface Omega_a ! sspCO2 = sea surface pCO2 real sspH, ssCO3, ssOc, ssOa, sspCO2 common /cdiaga_r/ sspH(imt,jmt), ssCO3(imt,jmt), ssOc(imt,jmt) common /cdiaga_r/ ssOa(imt,jmt), sspCO2(imt,jmt) real dmsk(is:ie,js:je) isp1 = is iem1 = ie jsp1 = js jem1 = je ! ! xconv is constant to convert piston_vel from cm/hr -> cm/s ! ! The following updates the gas transfer velocity paramter (a) to the ! most recent estimates of Wanninkhof 2014 LnO Methods, which is a similar ! value to the one used in recent publications by A. Schmittner ! here it is 100.*a*xconv (100 => m to cm, a=0.251, xconv=1/3.6e+05) xconv = 25.1/3.6e+05 ! The old param. is here as 100.*a*xconv (100 => m to cm, a=0.337, xconv=1/3.6e+05) ! uncomment to use old param. xconv = 33.7/3.6e+05 C2K = 273.15 pHlo = 6. pHhi = 10. sit_in = 7.6875e-03 !mol/m^3 pt_in = 0.5125e-3 !mol/m^3 atmpres = 1.0 !atm zero = 0. !----------------------------------------------------------------------- ! zero totals for new accumulation !----------------------------------------------------------------------- atatm = 0. flux(:,:,:) = 0. sbc(:,:,ihflx) = 0. sbc(:,:,isflx) = 0. sbc(:,:,iro) = 0. sbc(:,:,idicflx) = 0. sbc(:,:,ialkflx) = 0. sbc(:,:,io2flx) = 0. sbc(:,:,ipo4flx) = 0. !{ wyao 20170704 ! wyao 20170704} sbc(:,:,ino3flx) = 0. !{ wyao 20170704 ! wyao 20170704} !----------------------------------------------------------------------- ! update any atmospheric data !----------------------------------------------------------------------- call data (is, ie, js, je) !----------------------------------------------------------------------- ! calculate freezing point of sea water using UNESCO (1983) !----------------------------------------------------------------------- sspH(:,:) = 0. ssCO3(:,:) = 0. ssOc(:,:) = 0. ssOa(:,:) = 0. sspCO2(:,:) = 0. do j=jsp1,jem1 do i=isp1,iem1 if (tmsk(i,j) .ge. 0.5) then sss = 1000.0*sbc(i,j,isss) + 35.0 frzpt(i,j) = -.0575*sss + 1.71e-3*sss**1.5 - 2.155e-4*sss**2 sst = sbc(i,j,isst) ! put reasonable limits on sst and sss for chemistry flux calculations sst = min(35.,max(sst,-2.)) sss = min(45.,max(sss,0.)) ao = 1. - aice(i,j,2) !----------------------------------------------------------------------- ! calculate ocean carbon fluxes !----------------------------------------------------------------------- t_in = sst s_in = sss dic_in = sbc(i,j,issdic) ta_in = sbc(i,j,issalk) co2_in = co2ccn call co2calc_SWS (t_in, s_in, dic_in, ta_in, co2_in, pt_in &, sit_in, atmpres, zero, pHlo, pHhi, pH &, co2star, dco2star, pCO2, hplus, CO3 &, Omega_c, Omega_a) ! sspH is really H+ until the final averaging sspH(i,j) = hplus ssCO3(i,j) = CO3 ssOc(i,j) = Omega_c ssOa(i,j) = Omega_a sspCO2(i,j) = pCO2 ! Schmidt number for CO2 ! The Schmidt numbers have been update to follow Wanninkhof 2014 ! The old parameterization is commented out below for reference scco2 = 2116.8 - 136.25*sst + 4.7353*sst**2 & - 0.092307*sst**3 + 0.0007555*sst**4 ! scco2 = 2073.1 - 125.62*sst + 3.6276*sst**2 ! & - 0.043219*sst**3 piston_vel = ao*xconv*((sbc(i,j,iws)*0.01)**2) & *((scco2/660.)**(-0.5)) ! dic in umol cm-3 or (mol m-3) => flux in umol cm-2 s-1 sbc(i,j,idicflx) = piston_vel*dco2star !----------------------------------------------------------------------- ! calculate ocean oxygen fluxes !----------------------------------------------------------------------- ! Schmidt number for O2 ! The Schmidt numbers have been update to follow Wanninkhof 2014 ! The old parameterization is commented out below for reference sco2 = 1920.4 - 135.6*sst + 5.2122*sst**2 - 0.10939*sst**3 & + 0.00093777*sst**4 ! sco2 = 1638.0 - 81.83*sst + 1.483*sst**2 - 0.008004*sst**3 ! piston velocity for O2 piston_o2 = ao*xconv*((sbc(i,j,iws)*0.01)**2) & *(sco2/660.0)**(-0.5) ! oxygen saturation concentration [mol/m^3] f1 = alog((298.15 - sst)/(C2K + sst)) f2 = f1*f1 f3 = f2*f1 f4 = f3*f1 f5 = f4*f1 o2sat = exp (2.00907 + 3.22014*f1 + 4.05010*f2 & + 4.94457*f3 - 2.56847E-1*f4 + 3.88767*f5 & + sss*(-6.24523e-3 - 7.37614e-3*f1 - 1.03410e-2*f2 & - 8.17083E-3*f3) - 4.88682E-7*sss*sss) ! Convert from ml/l to mol/m^3 o2sat = o2sat/22391.6*1000.0 sbc(i,j,io2flx) = piston_o2*(o2sat - sbc(i,j,isso2)) endif enddo enddo !----------------------------------------------------------------------- ! set boundary conditions for carbon !----------------------------------------------------------------------- call setbcx (sbc(1,1,idicflx), imt, jmt) dmsk(:,:) = 1. call areaavg (sbc(1,1,idicflx), dmsk, avgflxc) carbemit = carbemit + co2emit*atmsa*segtim*daylen*1e-15 !----------------------------------------------------------------------- ! set boundary conditions for oxygen !----------------------------------------------------------------------- call setbcx (sbc(1,1,io2flx), imt, jmt) !----------------------------------------------------------------------- ! calculate CO2 forcing !----------------------------------------------------------------------- call co2forc !----------------------------------------------------------------------- ! set flags to calculate new coefficients !----------------------------------------------------------------------- newcoef(:,:) = .true. return end