! source file: /Users/nmengis/Documents/UVic_ESCM/2.10.backup/updates/02_CE_permafrost_merge_CMIP6forcing/source/common/gasbc.F 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, T_ph_srm, T_pc_srm, T_period, k_adg, time real delta include "size.h" include "param.h" include "pconst.h" include "stdunits.h" include "coord.h" include "csbc.h" include "mw.h" include "ice.h" include "switch.h" include "tmngr.h" include "cembm.h" include "atm.h" include "insolation.h" include "calendar.h" include "grdvar.h" include "levind.h" include "solve.h" include "mtlm.h" include "diaga.h" real cosz(is:ie,js:je) real dmsk(is:ie,js:je) isp1 = is + 1 iem1 = ie - 1 jsp1 = js + 1 jem1 = je - 1 ! 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 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(:,:,iat) = 0. sbc(:,:,irh) = 0. sbc(:,:,ipr) = 0. sbc(:,:,ips) = 0. sbc(:,:,iaws) = 0. sbc(:,:,iswr) = 0. sbc(:,:,idicflx) = 0. sbc(:,:,ic14flx) = 0. sbc(:,:,ialkflx) = 0. sbc(:,:,io2flx) = 0. sbc(:,:,ipo4flx) = 0. sbc(:,:,ino3flx) = 0. !----------------------------------------------------------------------- ! set solar constant !----------------------------------------------------------------------- call solardata !----------------------------------------------------------------------- ! update insolation for the current day !----------------------------------------------------------------------- ! subroutine decl is expecting a 365.25 day year calday = dayoyr*365.25/yrlen call decl (calday, eccen, obliq, mvelp, lambm0, sindec, eccf) i = (ie-is+1)*(je-js+1) ! get average zenith angle call zenith (i, c0, daylen, daylen, tlat, tlon, sindec, cosz) solins(is:ie,js:je) = solarconst*eccf*cosz(is:ie,js:je) !----------------------------------------------------------------------- ! set co2 concentration or emissions by tracking average co2 !----------------------------------------------------------------------- call co2ccndata !----------------------------------------------------------------------- ! set additional greenhouse gas forcing !----------------------------------------------------------------------- call aggdata !----------------------------------------------------------------------- ! update any atmospheric data !----------------------------------------------------------------------- call data (is, ie, js, je) !----------------------------------------------------------------------- ! calculate winds with new feedback !----------------------------------------------------------------------- call add_awind (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 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 c14 fluxes !----------------------------------------------------------------------- sbc(i,j,ic14flx) = piston_vel*((dco2star + co2star) & *(1 + dc14ccn*0.001)*rstd & - co2star*sbc(i,j,issc14)/sbc(i,j,issdic)) !----------------------------------------------------------------------- ! calculate ocean oxygen fluxes !----------------------------------------------------------------------- ! Schmidt number for O2 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)) else !----------------------------------------------------------------------- ! calculate land carbon fluxes !----------------------------------------------------------------------- ! convert from kg m-2 s-1 => umol cm-2 s-1 sbc(i,j,idicflx) = (sbc(i,j,inpp) - sbc(i,j,isr) & - sbc(i,j,iburn))*0.1/12.e-6 ! use the carbon flux scaled by rstd for c14 sbc(i,j,ic14flx) = (sbc(i,j,inpp) - sbc(i,j,isr) & - sbc(i,j,iburn))*rstd*0.1/12.e-6 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 dmsk(:,:) = 1. call areaavg (sbc(1,1,ic14flx), dmsk, avgflxc) !----------------------------------------------------------------------- ! set boundary conditions for c14 !----------------------------------------------------------------------- call setbcx (sbc(1,1,ic14flx), imt, jmt) !----------------------------------------------------------------------- ! 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. !----------------------------------------------------------------------- ! zero time averages if not in an averaging period !----------------------------------------------------------------------- if (.not. timavgperts) call ta_embm_tavg (is, ie, js, je, 0) !----------------------------------------------------------------------- ! zero time step integrals if not in an averaging period !----------------------------------------------------------------------- if (.not. tsiperts) call ta_embm_tsi (is, ie, js, je, 0) return end