subroutine gasbc (is, ie, js, je)

#if defined O_embm
!=======================================================================
!     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

      include "size.h"
      include "param.h"
      include "pconst.h"
      include "stdunits.h"
      include "coord.h"
      include "csbc.h"
# if defined O_mom
      include "mw.h"
# endif
# if defined O_ice
#  if defined O_ice_cpts
      include "cpts.h"
#  endif
      include "ice.h"
# endif
      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"
# if defined O_mtlm
      include "mtlm.h"
# endif
# if defined O_save_carbon_carbonate_chem
      include "diaga.h"
# endif

# if !defined O_embm_annual
      real cosz(is:ie,js:je)
# endif
      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
# if defined O_carbon
      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.
# endif
# if defined O_mom
#  if !defined O_constant_flux_reference || defined O_cfcs_data_transient

      dmsk(:,:) = 1.
      where (kmt(:,:) .eq. 0) dmsk(:,:) = 0.
#  endif
#  if !defined O_constant_flux_reference

!-----------------------------------------------------------------------
!     calculate new global average sea surface flux references
!-----------------------------------------------------------------------
      if (issdic .ne. 0)
     &  call areaavg (sbc(1,1,issdic), dmsk, gaost(idic))
      if (issalk .ne. 0)
     &  call areaavg (sbc(1,1,issalk), dmsk, gaost(ialk))
      if (isso2 .ne. 0)
     &  call areaavg (sbc(1,1,isso2), dmsk, gaost(io2))
      if (isspo4 .ne. 0)
     &  call areaavg (sbc(1,1,isspo4), dmsk, gaost(ipo4))
#   if !defined O_npzd_no_vflux
      if (issphyt .ne. 0)
     &  call areaavg (sbc(1,1,issphyt), dmsk, gaost(iphyt))
      if (isszoop .ne. 0)
     &  call areaavg (sbc(1,1,isszoop), dmsk, gaost(izoop))
      if (issdetr .ne. 0)
     &  call areaavg (sbc(1,1,issdetr), dmsk, gaost(idetr))
#   endif
      if (issno3 .ne. 0)
     &  call areaavg (sbc(1,1,issno3), dmsk, gaost(ino3))
#   if !defined O_npzd_no_vflux
      if (issdiaz .ne. 0)
     &  call areaavg (sbc(1,1,issdiaz), dmsk, gaost(idiaz))
#   endif
      if (issc14 .ne. 0)
     &  call areaavg (sbc(1,1,issc14), dmsk, gaost(ic14))
#  endif
#  if defined O_cfcs_data_transient

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

!-----------------------------------------------------------------------
!     zero totals for new accumulation
!-----------------------------------------------------------------------
      atatm = 0.
      flux(:,:,:) = 0.
# if defined O_plume
      subflux(:,:,:) = 0.
# endif
# if defined O_convect_brine
      cbf(:,:,:) = 0.
      cba(:,:,:) = 0.
# endif
      sbc(:,:,ihflx) = 0.
      sbc(:,:,isflx) = 0.
      sbc(:,:,iro) = 0.
# if defined O_mtlm
      sbc(:,:,iat) = 0.
      sbc(:,:,irh) = 0.
      sbc(:,:,ipr) = 0.
      sbc(:,:,ips) = 0.
      sbc(:,:,iaws) = 0.
      sbc(:,:,iswr) = 0.
# endif
# if defined O_carbon
      sbc(:,:,idicflx) = 0.
#  if defined O_carbon_14
      sbc(:,:,ic14flx) = 0.
#  endif
# endif
# if defined O_npzd_alk
      sbc(:,:,ialkflx) = 0.
# endif
# if defined O_npzd_o2
      sbc(:,:,io2flx) = 0.
# endif
# if defined O_npzd
      sbc(:,:,ipo4flx) = 0.
#  if !defined O_npzd_no_vflux
      sbc(:,:,iphytflx) = 0.
      sbc(:,:,izoopflx) = 0.
      sbc(:,:,idetrflx) = 0.
#  endif
#  if defined O_npzd_nitrogen
      sbc(:,:,ino3flx) = 0.
#   if !defined O_npzd_no_vflux
      sbc(:,:,idiazflx) = 0.
#   endif
#  endif
# endif
# if defined O_cfcs_data_transient
      sbc(:,:,icfc11flx) = 0.
      sbc(:,:,icfc12flx) = 0.
# endif
# if defined O_solar_data || defined O_solar_data_transient

!-----------------------------------------------------------------------
!     set solar constant
!-----------------------------------------------------------------------
      call solardata
# endif
# if !defined O_embm_annual

!-----------------------------------------------------------------------
!     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)
# endif
# if defined O_volcano_data_transient

!-----------------------------------------------------------------------
!     set anomalous volcanic forcing (modify solins)
!-----------------------------------------------------------------------
      call volcdata
      dmsk(:,:) = 1.
      call areaavg (solins, dmsk, tmp)
      if (tmp .gt. 0.) solins(:,:) = solins(:,:)*(tmp - volcfor)/tmp
# endif
# if defined O_co2emit_data_transient

!-----------------------------------------------------------------------
!     set co2 emissions
!-----------------------------------------------------------------------
      call co2emitdata
#  if defined O_carbon_co2_2d

!-----------------------------------------------------------------------
!     set co2 emissions distribution
!-----------------------------------------------------------------------
      call co2distdata
#  endif
# endif
# if defined O_co2ccn_data || defined O_co2ccn_data_transient || defined O_co2emit_track_co2

!-----------------------------------------------------------------------
!     set co2 concentration or emissions by tracking average co2
!-----------------------------------------------------------------------
      call co2ccndata
# endif
# if defined O_co2emit_track_sat || defined O_embm_vcs

!-----------------------------------------------------------------------
!     set co2 emissions by tracking average surface air temperature
!-----------------------------------------------------------------------
      call satdata
# endif
# if defined O_carbon_14
#  if defined O_c14ccn_data || defined O_c14ccn_data_transient

!-----------------------------------------------------------------------
!     set c14 concentration
!-----------------------------------------------------------------------
      call c14data
#   if defined O_c14ccn_data
      tdc14ccn = 0.
      tarea = 0.
#   endif
#  endif
# endif
# if defined O_cfcs_data_transient

!-----------------------------------------------------------------------
!     set CFC concentration
!-----------------------------------------------------------------------
      call cfcdata
# endif
# if defined O_aggfor_data_transient

!-----------------------------------------------------------------------
!     set additional greenhouse gas forcing
!-----------------------------------------------------------------------
      call aggdata
# endif

!-----------------------------------------------------------------------
!     update any atmospheric data
!-----------------------------------------------------------------------
      call data (is, ie, js, je)
# if defined O_embm_awind

!-----------------------------------------------------------------------
!     calculate winds with new feedback
!-----------------------------------------------------------------------
      call add_awind (is, ie, js, je)
# endif
# if defined O_embm && defined O_sealev_data_transient

!-----------------------------------------------------------------------
!     set anomalous sea level
!-----------------------------------------------------------------------
#  if defined O_sealev_data_transient &&  defined O_sealev_salinity
      dsealev = sealev
      call sealevdata
      dsealev = sealev - dsealev
#  else
      call sealevdata
#  endif

      area = ocnsa/atmsa
      do j=2,jmtm1
        do i=2,imtm1
          elev_sealev(i,j) = sealev*(1. - area)*tmsk(i,j)
     &                     - sealev*area*(1. - tmsk(i,j))
        enddo
      enddo
# endif

!-----------------------------------------------------------------------
!     calculate freezing point of sea water using UNESCO (1983)
!-----------------------------------------------------------------------

# if defined O_save_carbon_carbonate_chem
      sspH(:,:) = 0.
      ssCO3(:,:) = 0.
      ssOc(:,:) = 0.
      ssOa(:,:) = 0.
      sspCO2(:,:) = 0.
# endif
      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
# if defined O_mom
#  if defined O_carbon || defined O_npzd_o2 || defined O_cfcs_data_transient
            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.))
#   if defined O_ice
#    if defined O_ice_cpts
            ao = 1.
            do n=1,ncat
              ao =  ao - A(i,j,2,n)
            enddo
#    else
            ao = 1. - aice(i,j,2)
#    endif
#   else
            ao = 1.
#   endif
#  endif
#  if defined O_carbon

!-----------------------------------------------------------------------
!           calculate ocean carbon fluxes
!-----------------------------------------------------------------------
            t_in = sst
            s_in = sss
            dic_in = sbc(i,j,issdic)
#   if defined O_npzd_alk
            ta_in = sbc(i,j,issalk)
#   else
            ta_in = 2.36775*sss/(socn*1000.)
#   endif
#   if defined O_carbon_co2_2d
            co2_in = at(i,j,2,ico2)
#   else
            co2_in = co2ccn
#   endif
            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, dpco2, CO3
     &,                       Omega_c, Omega_a)
#   if defined O_save_carbon_carbonate_chem
            sspH(i,j) = pH
            ssCO3(i,j) = CO3
            ssOc(i,j) = Omega_c
            ssOa(i,j) = Omega_a
            sspCO2(i,j) = pCO2
#   endif
!           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
#   if defined O_carbon_co2_2d
!           convert from umol cm-2 s-1 => g cm-2 s-1
            flux(i,j,ico2) = sbc(i,j,idicflx)*12.e-6 + flux(i,j,ico2)
#   endif
#   if defined O_carbon_14

!-----------------------------------------------------------------------
!           calculate ocean c14 fluxes
!-----------------------------------------------------------------------
#    if defined O_c14ccn_data
            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)
#    endif

            sbc(i,j,ic14flx) = piston_vel*((dco2star + co2star)
     &                         *(1 + dc14ccn*0.001)*rstd
     &                       - co2star*sbc(i,j,issc14)/sbc(i,j,issdic))
#   endif
#  endif
#  if defined O_npzd_o2

!-----------------------------------------------------------------------
!           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))
#  endif
#  if defined O_cfcs_data_transient

!-----------------------------------------------------------------------
!           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
# endif
# if defined O_carbon && defined O_mtlm

          else
!-----------------------------------------------------------------------
!           calculate land carbon fluxes
!-----------------------------------------------------------------------
#  if defined O_carbon_co2_2d
!           convert from kg m-2 s-1 => g cm-2 s-1
            flux(i,j,ico2) = (sbc(i,j,inpp) - sbc(i,j,isr)
     &                     - sbc(i,j,iburn))*0.1 + flux(i,j,ico2)
#  else
!           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
#  endif
#  if defined O_carbon_14
!           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
# endif
          endif
        enddo
      enddo

# if defined O_carbon

!-----------------------------------------------------------------------
!     set boundary conditions for carbon
!-----------------------------------------------------------------------
      call setbcx (sbc(1,1,idicflx), imt, jmt)
#  if defined O_carbon_co2_2d
#   if !defined O_carbon_uncoupled && !defined O_co2ccn_data && !defined O_co2ccn_data_transient
      flux(:,:,ico2) = flux(:,:,ico2) - co2emit*co2dist(:,:,2)
#   endif
      call setbcx (flux(1,1,ico2), imt, jmt)
#  else
      dmsk(:,:) = 1.
      call areaavg (sbc(1,1,idicflx), dmsk, avgflxc)
#   if !defined O_carbon_uncoupled && !defined O_co2ccn_data && !defined O_co2ccn_data_transient
      co2ccn = co2ccn + (co2emit - avgflxc*12.e-6)*segtim*daylen*gtoppm
#   endif
#  endif
#  if defined O_global_sums
!     convert from g cm-2 s-1 to umol s-1 for conversion later
      dtoic = dtoic - co2emit*atmsa*segtim*daylen/12e-6
#  endif
#  if defined O_carbon_14
      dmsk(:,:) = 1.
      call areaavg (sbc(1,1,ic14flx), dmsk, avgflxc)
#   if defined O_carbon_14_coupled
      c14ccn = c14ccn + (c14prod - avgflxc)*12.e-6*segtim*daylen*gtoppm
!     calculate dc14ccn from c14ccn and co2ccn
      dc14ccn = 1000.*(c14ccn/co2ccn/rstd - 1.)
#   endif

!-----------------------------------------------------------------------
!     set boundary conditions for c14
!-----------------------------------------------------------------------
      call setbcx (sbc(1,1,ic14flx), imt, jmt)
#   if defined O_c14ccn_data
      if (tarea .gt. 0) dc14ccn = tdc14ccn/tarea
#   endif
#  endif
# endif
# if defined O_npzd_o2

!-----------------------------------------------------------------------
!     set boundary conditions for oxygen
!-----------------------------------------------------------------------
      call setbcx (sbc(1,1,io2flx), imt, jmt)
# endif
# if defined O_cfcs_data_transient

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

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

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

!-----------------------------------------------------------------------
!     update boundary conditions over vegetation
!-----------------------------------------------------------------------
      call gvsbc
# endif
# if defined O_time_averages

!-----------------------------------------------------------------------
!     zero time averages if not in an averaging period
!-----------------------------------------------------------------------
      if (.not. timavgperts) call ta_embm_tavg (is, ie, js, je, 0)
# endif
# if defined O_time_step_monitor

!-----------------------------------------------------------------------
!     zero time step integrals if not in an averaging period
!-----------------------------------------------------------------------
      if (.not. tsiperts) call ta_embm_tsi (is, ie, js, je, 0)
# endif
#endif

      return
      end