! source file: /Users/oschlies/UVIC/master/testcase/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
!=======================================================================

      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"

      integer i, ie, iem1, is, isp1, j, je, jem1, js, jsp1, l

      real sss, sst, phlo, phhi, sit_in, pt_in, 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 n2osat, a1, a2, a3, a4, b1, b2, b3, beta, atn2o
      real sccfc, piston_cfc, sol_cfc, cfcsat, ao, tarea, tdc14ccn
      real h_r, d, f1, f2, f3, f4, f5
      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
      phlo = 6.
      phhi = 10.
      sit_in = 7.6875e-03 !mol/m^3
      pt_in = 0.5125e-3   !mol/m^3
      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(:,:,ialkflx) = 0.
      sbc(:,:,io2flx) = 0.
      sbc(:,:,in2oflx) = 0.
      sbc(:,:,iabioto2flx) = 0.
      sbc(:,:,ipo4flx) = 0.
      sbc(:,:,ino3flx) = 0.

!-----------------------------------------------------------------------
!     set solar constant
!-----------------------------------------------------------------------
      call solardata

!-----------------------------------------------------------------------
!     set co2 concentration
!-----------------------------------------------------------------------
      call co2data

!-----------------------------------------------------------------------
!     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, pt_in, sit_in
     &,                       phlo, phhi, 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

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

!-----------------------------------------------------------------------
!           calculate ocean N2O fluxes
!           Andreas Oschlies, 16.7.08
!-----------------------------------------------------------------------
!           Schmidt number for N2O assumed to be identical to that for O2
            sco2 = 1638.0 - 81.83*sst + 1.483*sst**2 - 0.008004*sst**3
!           piston velocity for N2O assumed to be identical to that for O2
            piston_o2 = ao*xconv*((sbc(i,j,iws)*0.01)**2)
     &                  *(sco2/660.0)**(-0.5)
! solubility coefficient for N2O in seawater in mol/l/atm (Weiss & Preiss, Mar.Chem., 8, 347-359, 1980)
! after Matlab script by A. Freing
            a1 = -165.8806
            a2 = 222.8743
            a3 = 92.0792
            a4 = -1.48425
            b1 = -0.056235
            b2 = 0.031619
            b3 = -0.0048472
            beta = exp(a1 + a2*(100./(sst+273.15))
     &          + a3*alog((sst+273.15)/100.) + a4*((sst+273.15)/100)**2
     &          +sss*(b1 + b2*((sst+273.15)/100.)
     &               + b3*((sst+273.15)/100.)**2))
! atmospheric N2O is assumed constant as 316 ppb (IPCC, 2001)
            atn2o = 3.16e-7
            n2osat = atn2o*beta
!           Convert from mol/l to mol/m^3
            n2osat = n2osat*1000.0
            sbc(i,j,in2oflx) = piston_o2*(n2osat - sbc(i,j,issn2o))

!-----------------------------------------------------------------------
!           calculate ocean abiotic O2 fluxes
!           Andreas Oschlies, 16.8.11
!-----------------------------------------------------------------------
!           Schmidt number for abiotic O2 identical to that for O2
            sco2 = 1638.0 - 81.83*sst + 1.483*sst**2 - 0.008004*sst**3
!           piston velocity for abiotic O2 identical to that 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,iabioto2flx) = piston_o2*(o2sat - sbc(i,j,issabioto2))

          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)
      dmsk(:,:) = 1.
      call areaavg (sbc(1,1,idicflx), dmsk, avgflxc)
      co2ccn = co2ccn - (co2flx + avgflxc)*segtim*12.e-6*daylen*fa

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

!-----------------------------------------------------------------------
!     set boundary conditions for N2O
!-----------------------------------------------------------------------
      call setbcx (sbc(1,1,in2oflx), imt, jmt)

!-----------------------------------------------------------------------
!     set boundary conditions for abiotic O2
!-----------------------------------------------------------------------
      call setbcx (sbc(1,1,iabioto2flx), 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