! source file: /sfs/fs6/home-geomar/smomw258/UVOK_1.1/Kiel_Feb_2019/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

      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"

      include "npzd.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
!     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(:,:,iat) = 0.
      sbc(:,:,irh) = 0.
      sbc(:,:,ipr) = 0.
      sbc(:,:,ips) = 0.
      sbc(:,:,iaws) = 0.
      sbc(:,:,iswr) = 0.

      sbc(:,:,idicflx) = 0.

      sbc(:,:,ic14flx) = 0.

      sbc(:,:,impflx) = 0.
      sbc(:,:,impaflx) = 0.
      sbc(:,:,imppflx) = 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

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

          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)

! plastic enters the ocean
! Geyer et al. (2017) plastic production growth rate 8.4%/y
! 1950 = 2MT, 2015 = 380 MT ever produced, 3% of this in the ocean
! Koelmans et al. (2017) and possibly 1% as MP???
      call getplastic
      do j=jsp1,jem1
        do i=isp1,iem1
          if (tmsk(i,j) .ge. 0.5) then
            sbc(i,j,impflx) = plastic(i,j,2)*pinp
          endif
        enddo
      enddo

      call setbcx (sbc(1,1,impflx), 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.

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