! 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.

      dmsk(:,:) = 1.
      where (kmt(:,:) .eq. 0) dmsk(:,:) = 0.

!     always update cfc flux references since initial is usually zero
      if (isscfc11 .ne. 0)
     &  call areaavg (sbc(1,1,isscfc11), dmsk, gaost(icfc11))
      if (isscfc12 .ne. 0)
     &  call areaavg (sbc(1,1,isscfc12), dmsk, gaost(icfc12))

!-----------------------------------------------------------------------
!     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.

      sbc(:,:,icfc11flx) = 0.
      sbc(:,:,icfc12flx) = 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 anomalous volcanic forcing (modify solins)
!-----------------------------------------------------------------------
      call volcdata
      dmsk(:,:) = 1.
      call areaavg (solins, dmsk, tmp)
      if (tmp .gt. 0.) solins(:,:) = solins(:,:)*(tmp - volcfor)/tmp

!-----------------------------------------------------------------------
!     set co2 emissions
!-----------------------------------------------------------------------
      call co2emitdata

!-----------------------------------------------------------------------
!     set CFC concentration
!-----------------------------------------------------------------------
      call cfcdata

!-----------------------------------------------------------------------
!     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))

!-----------------------------------------------------------------------
!           calculate ocean CFC11 fluxes
!-----------------------------------------------------------------------
            if (tlat(i,j) .gt. 10.) then
              cfc11ccn = cfc11ccnn
            elseif (tlat(i,j) .lt. -10.) then
              cfc11ccn = cfc11ccns
            else
              wt = (tlat(i,j) + 10.)/20.
              cfc11ccn = cfc11ccnn*wt + cfc11ccns*(1. - wt)
            endif
!           Schmidt number for CFC11
            sccfc = 3501.8 -210.31*sst + 6.1851*sst**2 -0.07513*sst**3
!           piston velocity for CFC
            piston_cfc = ao*xconv*((sbc(i,j,iws)*0.01)**2)
     &                 *(sccfc/660.0)**(-0.5)
!           cfc saturation concentration [mol/m^3]
            f1 = (sst + 273.16)*0.01
            d = (0.091459 - 0.0157274*f1)*f1 - 0.142382
            sol_cfc = exp(-229.9261 + 319.6552/f1 + 119.4471*alog(f1)
     $              - 1.39165*f1*f1 + sss*d )
!           conversion from mol/(l * atm) to mol/(m3 * pptv)

            cfcsat =  1.0e-12 *1000.*sol_cfc*cfc11ccn
            sbc(i,j,icfc11flx) = piston_cfc*(cfcsat - sbc(i,j,isscfc11))

!-----------------------------------------------------------------------
!           calculate ocean CFC12 fluxes
!-----------------------------------------------------------------------
            if (tlat(i,j) .gt. 10.) then
              cfc12ccn = cfc12ccnn
            elseif (tlat(i,j) .lt. -10.) then
              cfc12ccn = cfc12ccns
            else
              wt = (tlat(i,j) + 10.)/20.
              cfc12ccn = cfc12ccnn*wt + cfc12ccns*(1. - wt)
            endif
!           Schmidt number for CFC12
            sccfc = 3845.4 -228.95*sst + 6.1908*sst**2 -0.067430*sst**3
 !          piston velocity for CFC12
            piston_cfc = ao*xconv*((sbc(i,j,iws)*0.01)**2)
     &                 *(sccfc/660.0)**(-0.5)
!           cfc saturation concentration [mol/m^3]
            f1 = (sst + 273.16)*0.01
            d = (0.091015 - 0.0153924*f1)*f1 - 0.143566
            sol_cfc = exp(-218.0971 + 298.9702/f1 + 113.8049*alog(f1)
     $              - 1.39165*f1*f1 + sss*d )
!           conversion from mol/(l * atm) to mol/(m3 * pptv)
            cfcsat =  1.0e-12 *1000.*sol_cfc*cfc12ccn
            sbc(i,j,icfc12flx) = piston_cfc*(cfcsat - sbc(i,j,isscfc12))

          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)

      co2ccn = co2ccn + (co2emit - avgflxc*12.e-6)*segtim*daylen*gtoppm

      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)

!-----------------------------------------------------------------------
!     set boundary conditions for CFC11
!-----------------------------------------------------------------------
      call setbcx (sbc(1,1,icfc11flx), imt, jmt)

!-----------------------------------------------------------------------
!     set boundary conditions for CFC12
!-----------------------------------------------------------------------
      call setbcx (sbc(1,1,icfc12flx), imt, jmt)

!-----------------------------------------------------------------------
!     calculate CO2 forcing
!-----------------------------------------------------------------------
      call co2forc

!-----------------------------------------------------------------------
!     set flags to calculate new coefficients
!-----------------------------------------------------------------------
      newcoef(:,:) = .true.

!-----------------------------------------------------------------------
!     update boundary conditions over vegetation
!-----------------------------------------------------------------------
      call gvsbc

!-----------------------------------------------------------------------
!     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