! source file: /sfs/fs1/work-geomar6/smomw113/UVic_ESCM/2.9_dk_0814/updates/07_aouassim/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

      real sst_ase_anomaly

      real current_time

      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 "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
!     here it is 100.*a*xconv (100 => m to cm, a=0.337, xconv=1/3.6e+05)
      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.

      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(:,:,idicflx) = 0.

      sbc(:,:,ialkflx) = 0.

      sbc(:,:,io2flx) = 0.

      sbc(:,:,io2abiotflx) = 0.

      sbc(:,:,ipo4flx) = 0.

      sbc(:,:,ino3flx) = 0.

      sbc(:,:,icfc11flx) = 0.
      sbc(:,:,icfc12flx) = 0.

      sst_ase_anomaly = 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 emissions
!-----------------------------------------------------------------------
      call co2emitdata

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

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

            current_time = year0 + accel_yr0
     &         + (relyr - accel_yr0)*accel

! wkoeve@geomar.de, 22.06.2016

! Note, wk, 22.12.17
! Why do I call i-1, j-1?
! In sstdata, the dimenstion is datasst(1:imtm2,1:jmtm2,1:ln),
! i.e. without the edge cells, ergo i, j run 1:100, each.
! Here I loop over j=jsp1,jem1 and i=isp1,iem1
! hence from 2:101, for both.
            call sstdata(i-1,j-1
     &,       sst_ase_anomaly)

            t_in = sst - sst_ase_anomaly

            if ((i.eq.50) .and. (j.eq.50)) then
               write(*,*) current_time,'  ',sst_ase_anomaly
            endif

! Note, wk, 22.12.2017: The dimensions of
! freeout1(imt,jmt) are 1:102, earch.
! Since only 2:101 will be saved to file, I actually
! have to right into i,j since their min, max values are
! 2 and 101, respectively
            ase_sst(i,j) = t_in

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

!           use the t_in read in for CO2-fluxes
            f1 = alog((298.15 - t_in)/(C2K + t_in))

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

            sbc(i,j,io2abiotflx) = piston_o2 * (o2sat
     &         - sbc(i,j,isso2abiot))

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

          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

!-----------------------------------------------------------------------
!     set boundary conditions for oxygen
!-----------------------------------------------------------------------
      call setbcx (sbc(1,1,io2flx), imt, jmt)

!-----------------------------------------------------------------------
!     set zonal boundary conditions for abiotic oxygen (in common/util.F)
!-----------------------------------------------------------------------
      call setbcx (sbc(1,1,io2abiotflx), 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.

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