! source file: /home/aschmitt/UVic/c13_lightbug_fixed/updates/gasbc.F subroutine gasbc (is, ie, js, je) !======================================================================= ! calculate boundary conditions for the atmospheric model ! based on code by: A. Fanning and M. Eby !======================================================================= c implicit none include "param.h" include "coord.h" include "csbc.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 "timeavgs.h" include "npzd.h" include "carbon.h" integer i, ie, iem1, is, isp1, j, je, jem1, js, jsp1, l real sss, sst, atmpres, xconv, ta_in real dic_in, ph, co2_in, co2star, dco2star, dco2surf, dpco2 real pco2surf, scco2, piston_vel, avgflxc, fa, calday, f, sco2 real o2sat, o2sato, o2surf, piston_o2, cfc11ccn, cfc12ccn, wt real sccfc, piston_cfc, sol_cfc, cfcsat, ao, tarea, tdc14ccn real h_r, d, f1, f2, f3, f4, f5 real ak, aaqg, adicg, dc13ccn, r13a, r13dic 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 ! here it is 100.*a*xconv (100 => m to cm, a=0.337, xconv=1/3.6e+05) xconv = 33.7/3.6e+05 cAAA modify air-sea gas exchange by constant factor xconv = xconv*.75 atmpres = 1.0 !atm co2flx = 0. ! fa is used in converting g carbon cm-2 => ppmv CO2 ! 4.138e-7 => 12e-6 g/umol carbon / 29 g/mol air fa = 1./(4.138e-7*rhoatm*shc) !----------------------------------------------------------------------- ! 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(:,:,idic13flx) = 0. sbc(:,:,ialkflx) = 0. sbc(:,:,io2flx) = 0. sbc(:,:,ipo4flx) = 0. sbc(:,:,ino3flx) = 0. !----------------------------------------------------------------------- ! set solar constant !----------------------------------------------------------------------- call solardata !----------------------------------------------------------------------- ! set c14 concentration !----------------------------------------------------------------------- call c14data tdc14ccn = 0. tarea = 0. !----------------------------------------------------------------------- ! 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) call zenith (i, c0, daylen, daylen, tlat, tlon, sindec, cosz) solins(:,:) = solarconst*eccf*cosz(:,:) !----------------------------------------------------------------------- ! update any atmospheric data !----------------------------------------------------------------------- call atmos !----------------------------------------------------------------------- ! calculate freezing point of sea water using UNESCO (1983) !----------------------------------------------------------------------- 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) ao = 1. - aice(i,j,1) !----------------------------------------------------------------------- ! calculate ocean carbon fluxes !----------------------------------------------------------------------- ta_in = sbc(i,j,issalk) dic_in = sbc(i,j,issdic) co2_in = co2ccn call co2calc_SWS (sst, sss, dic_in, ta_in, ph, co2_in & , atmpres, co2star, dco2star, pCO2surf, dpco2) ! 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 c fractionation factors from Zhang et al. 1995 Geochi. Cosmochi. Acta ak = 0.99915 ! kinetic fractionation aaqg = 0.998764 ! aquatic - gas fractionation adicg = 1.01051 - 1.05e-4*sst ! DIC-gas fractionation c ak = 1. c aaqg = 1. c adicg = 1. dc13ccn = -6.5 ! atmospheric d13C concentration c dc13ccn = 0. r13a = (dc13ccn*0.001 + 1.)*rc13std r13dic = sbc(i,j,issdic13)/sbc(i,j,issdic) c r13dic = r13std c sbc(i,j,idic13flx) = 0. sbc(i,j,idic13flx) = piston_vel*ak*aaqg & *((dco2star + co2star)*r13a - co2star*r13dic/adicg) c sbc(i,j,idic13flx) = sbc(i,j,idicflx)*rc13std !----------------------------------------------------------------------- ! calculate ocean c14 fluxes !----------------------------------------------------------------------- if (tlat(i,j) .gt. 20.) then dc14ccn = dc14ccnn elseif (tlat(i,j) .lt. -20.) then dc14ccn = dc14ccns else dc14ccn = dc14ccne endif tarea = tarea + dxt(i)*dyt(j)*cst(j) tdc14ccn = tdc14ccn + dc14ccn*dxt(i)*dyt(j)*cst(j) 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)/(273.15 + 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))*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)) & *rstd*0.1/12.e-6 endif enddo enddo !----------------------------------------------------------------------- ! set boundary conditions for carbon !----------------------------------------------------------------------- call setbcx (sbc(1,1,idicflx), imt, jmt) !----------------------------------------------------------------------- ! set boundary conditions for c14 !----------------------------------------------------------------------- call setbcx (sbc(1,1,ic14flx), imt, jmt) if (tarea .gt. 0) dc14ccn = tdc14ccn/tarea !----------------------------------------------------------------------- ! 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(1,1) = .true. newcoef(2,1) = .true. newcoef(1,2) = .true. newcoef(2,2) = .true. !----------------------------------------------------------------------- ! zero time averages if not in an averaging period !----------------------------------------------------------------------- if (.not. timavgperts) call ta_embm_snap (is, ie, js, je, 0) !----------------------------------------------------------------------- ! zero time step integrals if not in an averaging period !----------------------------------------------------------------------- if (.not. tsiperts) call ta_embm_tsi (0) return end