subroutine gosbc (is, ie, js, je)

#if defined O_mom || defined O_embm
!=======================================================================
!     calculate the average fluxes for next ocean time step
!=======================================================================

      implicit none

      integer ie, is, je, js, i, j, nc
      integer iem1, isp1, jem1, jsp1, k, nn

      real f1, f1a, f1l, fh, fs, fwcflx, fwaflx, time
      real area, tarea, tsflx, rsocn, tmp, fg, fgs, sulphfac
# if defined O_oc_alk_Pg
      real area_aoa, tarea_aoa, aoa_add
# endif
# if defined O_oc_alk_as_mpi
      real area_aoa, tarea_aoa, aoa_add, tt, co2_target, ratio
      real month1,month2,month3,month4,month5,month6,month7,month8
      real month9,month10,month11,month12,ttt
# endif

      include "size.h"
      include "param.h"
      include "pconst.h"
      include "stdunits.h"
      include "calendar.h"
      include "csbc.h"
      include "coord.h"
      include "grdvar.h"
      include "tmngr.h"
      include "switch.h"
# if defined O_embm
      include "cembm.h"
      include "atm.h"
# endif
# 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
# if defined O_mtlm
      include "mtlm.h"
# endif
      include "levind.h"
# if defined O_fwa
      include "fwa.h"
      include "cregin.h"
      include "tmngr.h"
# endif
# if defined O_sed && !defined O_sed_uncoupled
      include "sed.h"
# endif
# if defined O_sulphate_data || defined O_sulphate_data_transient
      include "insolation.h"
      real tot, twt, cz(1), zero(1)
# endif

# if defined O_embm_read_sflx || defined O_embm_write_sflx
      integer ntrec
      save ntrec
      data ntrec /0/
# endif

      isp1 = is + 1
      iem1 = ie - 1
      jsp1 = js + 1
      jem1 = je - 1

# if defined O_mom && defined O_embm
      f1 = 1./atatm
      fh = 2.389e-8/atatm
      fs = -socn/atatm

!-----------------------------------------------------------------------
!     calculate average net fluxes. convert heat flux to cal/cm**2/s
!     from kW/m**2 and fresh water flux (cm/s) to an apparent salt
!     flux (g/cm**2/s) using global ocean average salinity, socn
!-----------------------------------------------------------------------

      do j=2,jmtm1
        do i=2,imtm1
          if (tmsk(i,j) .ge. 0.5) then
#  if defined O_plume
#   if !defined O_plume_brine
            subflux(i,j,1) = flux(i,j,isat)
            if (subflux(i,j,1) .gt. 0.) subflux(i,j,1) = 0.
            flux(i,j,isat) = flux(i,j,isat) - subflux(i,j,1)
            subflux(i,j,1) = fh*subflux(i,j,1)

            subflux(i,j,2) = flux(i,j,ishum)
#   endif
            if (subflux(i,j,2) .gt. 0.) subflux(i,j,2) = 0.
            flux(i,j,ishum) = flux(i,j,ishum) - subflux(i,j,2)
            subflux(i,j,2) = fs*subflux(i,j,2)
#  endif
#  if defined O_convect_brine
            cba0(i,j) = 0.
            do nc=0,ncat
              cba(i,j,nc) = f1*cba(i,j,nc)
              if (cba(i,j,nc) .gt. 0.) then
                cbf(i,j,nc) = fs*cbf(i,j,nc)/cba(i,j,nc)
                cba0(i,j) = cba0(i,j) + cba(i,j,nc)
              else
                cbf(i,j,nc) = 0.
                cba(i,j,nc) = 0.
              endif
            enddo
            if (cba0(i,j) .gt. 1.) then
              if (cba0(i,j) .gt. 1.000001) then
                print*, "==> Warning: ice area > 1: ", cba0(i,j)
                endif
              cba0(i,j) = 1.
            endif
            cba0(i,j) = 1. - cba0(i,j)
#  endif
            sbc(i,j,ihflx) = sbc(i,j,ihflx) + fh*flux(i,j,isat)
!           add virtual fluxes of salinity
            sbc(i,j,isflx) = sbc(i,j,isflx) + fs*flux(i,j,ishum)

          else
#  if defined O_plume
            subflux(i,j,1) = 0.
            subflux(i,j,2) = 0.
#  endif
#  if defined O_convect_brine
            cbf(i,j,:) = 0.
            cba(i,j,:) = 0.
            cba0(i,j) = 1.
#  endif
            sbc(i,j,ihflx) = 0.
            sbc(i,j,isflx) = 0.
          endif
          if (umsk(i,j) .ge. 0.5) then
#  if defined O_ice_evp || defined O_embm_awind
            sbc(i,j,itaux) = f1*flux(i,j,nat+1)
            sbc(i,j,itauy) = f1*flux(i,j,nat+2)
#  endif
          else
            sbc(i,j,itaux) = 0.
            sbc(i,j,itauy) = 0.
          endif
        enddo
      enddo

      call setbcx (sbc(1,1,ihflx), imt, jmt)
      call setbcx (sbc(1,1,isflx), imt, jmt)
      call setbcx (sbc(1,1,itaux), imt, jmt)
      call setbcx (sbc(1,1,itauy), imt, jmt)
#  if defined O_convect_brine
      do nc=0,ncat
        call setbcx (cbf(1,1,nc), imt, jmt)
        call setbcx (cba(1,1,nc), imt, jmt)
      enddo
      call setbcx (cba0, imt, jmt)
#  endif
# endif

!-----------------------------------------------------------------------
!     update boundary conditions from the land model
!     do this now instead of in gasbc so fields can be written out
!-----------------------------------------------------------------------

      f1l = 0.
      f1a = 0.
# if defined O_embm
      if (atatm .ne. 0.) f1a = 1.0/atatm
# endif
# if defined O_mtlm
      if (atlnd .ne. 0.) f1l = 1.0/atlnd
      do j=2,jmtm1
        do i=2,imtm1
          if (land_map(i,j) .ne. 0) then
            sbc(i,j,iro) = sbc(i,j,iro)*f1l
            sbc(i,j,isca) = sbc(i,j,isca)*f1l
            sbc(i,j,ievap) = sbc(i,j,ievap)*f1l
            sbc(i,j,ilwr) = sbc(i,j,ilwr)*f1l
            sbc(i,j,isens) = sbc(i,j,isens)*f1l
#  if defined O_carbon
            sbc(i,j,inpp) =  sbc(i,j,inpp)*f1l
            sbc(i,j,isr) =  sbc(i,j,isr)*f1l
#   if defined O_reproduce_updates_01
#    if defined O_mtlm_segday
            sbc(i,j,iburn) =  sbc(i,j,iburn)*f1l/segtim
#    else
            sbc(i,j,iburn) =  sbc(i,j,iburn)*f1l
#    endif
#   endif
#  endif
          else
            sbc(i,j,iro) = sbc(i,j,iro)*f1a
            sbc(i,j,ievap) = 0.
            sbc(i,j,ilwr) = 0.
            sbc(i,j,isens) = 0.
#  if defined O_mtlm && defined O_carbon
            sbc(i,j,inpp) = 0.
            sbc(i,j,isr) = 0.
#   if defined O_reproduce_updates_01
            sbc(i,j,iburn) = 0.
#   endif
#  endif
          endif
        enddo
      enddo
      call setbcx (sbc(1,1,isca), imt, jmt)
      call setbcx (sbc(1,1,ievap), imt, jmt)
      call setbcx (sbc(1,1,ilwr), imt, jmt)
      call setbcx (sbc(1,1,isens), imt, jmt)
#  if defined O_carbon
      call setbcx (sbc(1,1,inpp), imt, jmt)
      call setbcx (sbc(1,1,isr), imt, jmt)
#  endif
# else
      sbc(:,:,iro) = sbc(:,:,iro)*f1a
# endif
      call setbcx (sbc(1,1,iro), imt, jmt)

# if defined O_embm
!-----------------------------------------------------------------------
!     zero diagnostic for river discharge and call river model
!-----------------------------------------------------------------------
      disch(:,:) = 0.
      call rivmodel

#  if defined O_sed && !defined O_sed_uncoupled
!-----------------------------------------------------------------------
!     apply CaCO3 weathering flux proportional to discharge
!-----------------------------------------------------------------------
#   if defined O_sed_constrain_weath
      if (weathflx .lt. 0.) weathflx = 0.
#   endif
#   if defined O_save_carbon_totals
      dicwflx =  weathflx
#   endif
      tmp = 0.
      do j=2,jmtm1
        do i=2,imtm1
          if (kmt(i,j) .ne. 0.) then
            tmp = tmp + disch(i,j)*dxt(i)*dyt(j)*cst(j)
          endif
        enddo
      enddo
      tmp = weathflx/tmp
      do j=2,jmtm1
#   if defined O_global_sums || defined O_save_carbon_totals
        fgs = dyt(j)*cst(j)*segtim/secday
#   endif
        do i=2,imtm1
          if (kmt(i,j) .ne. 0.) then
#   if defined O_carbon
            sbc(i,j,idicflx) = sbc(i,j,idicflx) + disch(i,j)*tmp
#   if defined O_npzd_alk
            sbc(i,j,ialkflx) = sbc(i,j,ialkflx) + disch(i,j)*2.*tmp
#   endif
#    if defined O_save_carbon_totals
            carblith = carblith - disch(i,j)*dxt(i)*tmp*fgs
     &              - sbc(i,j,ibdicfx)*dxt(i)*fgs*dzt(kmt(i,j))
#    endif
#    if defined O_global_sums
            dtoic = dtoic - disch(i,j)*dxt(i)*tmp*fgs
#    endif
#   endif
          endif
        enddo
      enddo

#  endif
# endif
# if defined O_sealev_data_transient && defined O_sealev_salinity
!-----------------------------------------------------------------------
!     add in flux from sea level change
!-----------------------------------------------------------------------
      do j=2,jmtm1
        do i=2,imtm1
          if (kmt(i,j) .gt. 0)
     &      sbc(i,j,isflx) = sbc(i,j,isflx) + fs*dsealev
        enddo
      enddo

# endif
# if defined O_mom && defined O_embm
#  if defined O_fwa
!-----------------------------------------------------------------------
!     add additional fresh water flux anomaly
!-----------------------------------------------------------------------

      time = year0 + accel_yr0 + (relyr - accel_yr0)*accel
      if (time .ge. fwayri .and. time .le. fwayrf) then

        if (areafwa .gt. 0) then
!         fwaflxi is in Sv (1e6 m3 s-1) and fwarate is in Sv/yr
          fwaflx = fwaflxi + (time - fwayri)*fwarate
!         convert to flux in g salt cm-2 s-1
          fwaflx = -socn*fwaflx*1.e12/areafwa
#   if defined O_fwa_precip
          call areaavg (precip, fwawt, tmp)
          if (tmp .gt. 0) fwaflx = fwaflx/tmp
#   endif
          do j=2,jmtm1
            do i=2,imtm1
              sbc(i,j,isflx) = sbc(i,j,isflx) + fwaflx*fwawt(i,j)
#   if defined O_fwa_precip
     &                         *precip(i,j)
#   endif
            enddo
          enddo
        endif

        if (compensate .and. areafwc .gt. 0) then
!         fwaflxi is in Sv (1e6 m3 s-1) and fwarate is in Sv/yr
          fwcflx = fwaflxi + (time - fwayri)*fwarate
!         convert to opposite flux in g salt cm-2 s-1
          fwcflx = socn*fwcflx*1.e12/areafwc
#   if defined O_fwa_compevap
          call areaavg (evap, fwcwt, tmp)
          if (tmp .gt. 0) fwcflx = fwcflx/tmp
#   endif
          do j=2,jmtm1
            do i=2,imtm1
              sbc(i,j,isflx) = sbc(i,j,isflx) + fwcflx*fwcwt(i,j)
#   if defined O_fwa_compevap
     &                         *evap(i,j)
#   endif
            enddo
          enddo
        endif

      endif

#  endif
#  if defined O_carbon || defined O_npzd_alk || defined O_npzd_o2 || defined O_npzd || defined O_cfcs_data || defined O_cfcs_data_transient
!-----------------------------------------------------------------------
!     add normalized virtual fluxes to other tracers
!-----------------------------------------------------------------------
! Note that the virtual flux is not written out in the diagnostic file unless
! the O_save_fluxes_with_virtual option is turned on.  The virtual flux is still
! calculated in the model, but the flux is subtracted from the F_XXX diagnostic  
! output so that without this option on there appears to be no flux.
!-----------------------------------------------------------------------
      tarea = 0.
      tsflx = 0.
      rsocn = 1./socn
      do j=2,jmtm1
        do i=2,imtm1
          if (tmsk(i,j) .ge. 0.5) then
            area = dxt(i)*dyt(j)*cst(j)
            tarea = tarea + area
            tsflx = tsflx + sbc(i,j,isflx)*area
          endif
        enddo
      enddo
      tsflx = tsflx/tarea
      do j=2,jmtm1
        do i=2,imtm1
          if (tmsk(i,j) .ge. 0.5) then
            tmp = (sbc(i,j,isflx) - tsflx)*rsocn
            vflux(i,j) = tmp
#   if defined O_carbon
            sbc(i,j,idicflx) = sbc(i,j,idicflx) + gaost(idic)*tmp
#    if defined O_carbon_14
            sbc(i,j,ic14flx) = sbc(i,j,ic14flx) + gaost(ic14)*tmp
#    endif
#   endif
#   if defined O_npzd_alk
            sbc(i,j,ialkflx) = sbc(i,j,ialkflx) + gaost(ialk)*tmp
#   endif
#   if defined O_npzd_o2
            sbc(i,j,io2flx) = sbc(i,j,io2flx) + gaost(io2)*tmp
#   endif
#   if defined O_npzd
            sbc(i,j,ipo4flx) = sbc(i,j,ipo4flx) + gaost(ipo4)*tmp
#    if !defined O_npzd_no_vflux
            sbc(i,j,iphytflx) = sbc(i,j,iphytflx) + gaost(iphyt)*tmp
            sbc(i,j,izoopflx) = sbc(i,j,izoopflx) + gaost(izoop)*tmp
            sbc(i,j,idetrflx) = sbc(i,j,idetrflx) + gaost(idetr)*tmp
#    endif
#    if defined O_npzd_nitrogen
            sbc(i,j,ino3flx) = sbc(i,j,ino3flx) + gaost(ino3)*tmp
#     if !defined O_npzd_no_vflux
            sbc(i,j,idiazflx) = sbc(i,j,idiazflx) + gaost(idiaz)*tmp
#     endif
#    endif
#   endif
#   if defined O_cfcs_data || defined O_cfcs_data_transient
            sbc(i,j,icfc11flx) = sbc(i,j,icfc11flx) + gaost(icfc11)*tmp
            sbc(i,j,icfc12flx) = sbc(i,j,icfc12flx) + gaost(icfc12)*tmp
#   endif
#   if defined O_oc_alk_as_a_function_of_emissions
!------------------------------------------------------------------------------------
! Artifically ocean alkalinization (aoa) adds alkalinity to the surface ocean. 
! This option increase alkalinity as a function of CO2 emissions.
! nn is the intensity for aoa with co2emit
! output, the unit is umol/s/cm2, default nn is "1", change the number as you need to.
!--------------------------------------------------------------------------------------
            nn=1  
            sbc(i,j,ialkflx) = sbc(i,j,ialkflx) + nn*co2emit/(12.e-6)
!
#   endif
          endif
        enddo
      enddo
      
#   if defined O_npzd_alk
#   if defined O_oc_alk_Pg
!-------------------------------------------------------------------------------------
! This method should add XX Pg yr-1 of Ca(OH)2 to the non-arctic surface ocean (j=18,90).  To do  
! this convert 10 Pg to mol and divide by the area of the ocean. Then add it. Remember that:
! 1 g Ca(OH)2 = 74.1 mol.  10 Pg yr-1 equals ~ 4279324154000 umol s-1 cm-2 Ca(OH)2 and that for 
! each mole of Ca(OH)2 added alkalinity increases by 2 (this is why the # is multiplied by 2 below).
! The 10 Pg option is the one use in the Keller et. al. 2014 Nat. Comm. paper.
!-------------------------------------------------------------------------------------
      tarea_aoa = 0.
      do j=18,90
        do i=2,imtm1
          if (tmsk(i,j) .ge. 0.5) then
             area_aoa = dxt(i)*dyt(j)*cst(j)
             tarea_aoa = tarea_aoa + area_aoa
          endif
        enddo
      enddo
# if defined O_oc_alk_inc_w_time
! This option will start off adding XX Pg and then increase that amount
! as time progresses.  This was not used in Keller et al., 2014
      time = (year0 + accel_yr0 + (relyr - accel_yr0)*accel)-1565
!      aoa_add = (4279324154000)/tarea_oa  ! 5 Pg addition
      aoa_add = (2*(4279324154000+33924*time))/tarea_aoa  ! 10 Pg addition
!      aoa_add = (4*4279324154000)/tarea_oa   ! 20 Pg addition
# else
# if defined O_CDR-MIP_C4
! To add Pmol of alkalinity per year convert the value to umol s-1 cm-2
      aoa_add = (4439370877727)/tarea_aoa  ! 0.14 Pmol TA yr-1 addition
# else
!      aoa_add = (4279324154000)/tarea_aoa  ! 5 Pg addition
      aoa_add = (2*4279324154000)/tarea_aoa  ! 10 Pg addition
!      aoa_add = (4*4279324154000)/tarea_aoa   ! 20 Pg addition
# endif
# endif
      do j=18,90
        do i=2,imtm1
          if (tmsk(i,j) .ge. 0.5) then
!             area_aoa = dxt(i)*dyt(j)*cst(j)
             sbc(i,j,ialkflx) = sbc(i,j,ialkflx) + aoa_add
          endif
        enddo
      enddo
#   endif

#   if defined O_oc_alk_as_mpi
! This code will add alkalinity in an attempt to pull the CO2 emissions down from
! RCP 8.5 to 4.5 starting in the year 2006.  It should mimic what is done for the
! ComparCE project with the MPI-ESM.
      tt=year0 + accel_yr0 + (relyr - accel_yr0)*accel-2006.
      ttt=tt-int(tt)

      if(tt.lt.70)then
      co2_target=1.597e-7*tt**5-2.8627e-5*tt**4+0.0013*tt**3-
     & 0.0142*tt**2+3.2492*tt+377.5062
      else
      co2_target=5.6466e-4*tt**3-0.1237*tt**2+9.0482*tt+341.0481
      endif

      ratio=co2ccn/co2_target
!      write(*,*) "co2_target",co2_target

! Jan
      if((ttt .lt. 0.028) .and. (ttt .gt. 0.027) .and. (ratio.gt.1.01))
     & then
        do j=2,jmtm1
         do i=2,imtm1
          if (tmsk(i,j) .ge. 0.5) then
            sbc(i,j,ialkflx) = sbc(i,j,ialkflx) + 5*co2emit/(12.e-7)
!             sbc(i,j,ialkflx) = sbc(i,j,ialkflx)
          endif
         enddo
        enddo
      endif

! Feb
      if((ttt .lt. 0.11) .and. (ttt .gt. 0.109) .and. (ratio.gt.1.01))
     & then
        do j=2,jmtm1
         do i=2,imtm1
          if (tmsk(i,j) .ge. 0.5) then
            sbc(i,j,ialkflx) = sbc(i,j,ialkflx) + 5*co2emit/(12.e-7)
!             sbc(i,j,ialkflx) = sbc(i,j,ialkflx)
          endif
         enddo
        enddo
      endif

! Mar
      if((ttt .lt. 0.192) .and. (ttt .gt. 0.191) .and. (ratio.gt.1.01))
     & then
        do j=2,jmtm1
         do i=2,imtm1
          if (tmsk(i,j) .ge. 0.5) then
            sbc(i,j,ialkflx) = sbc(i,j,ialkflx) + 5*co2emit/(12.e-7)
!             sbc(i,j,ialkflx) = sbc(i,j,ialkflx)
          endif
         enddo
        enddo
      endif

! Apr
      if((ttt .lt. 0.274) .and. (ttt .gt. 0.273) .and. (ratio.gt.1.01))
     & then
        do j=2,jmtm1
         do i=2,imtm1
          if (tmsk(i,j) .ge. 0.5) then
            sbc(i,j,ialkflx) = sbc(i,j,ialkflx) + 5*co2emit/(12.e-7)
!             sbc(i,j,ialkflx) = sbc(i,j,ialkflx)
          endif
         enddo
        enddo
      endif

! May
      if((ttt .lt. 0.357) .and. (ttt .gt. 0.356) .and. (ratio.gt.1.01))
     & then
        do j=2,jmtm1
         do i=2,imtm1
          if (tmsk(i,j) .ge. 0.5) then
            sbc(i,j,ialkflx) = sbc(i,j,ialkflx) + 5*co2emit/(12.e-7)
!             sbc(i,j,ialkflx) = sbc(i,j,ialkflx)
          endif
         enddo
        enddo
      endif

! Jun
      if((ttt .lt. 0.439) .and. (ttt .gt. 0.438) .and. (ratio.gt.1.01))
     & then
        do j=2,jmtm1
         do i=2,imtm1
          if (tmsk(i,j) .ge. 0.5) then
            sbc(i,j,ialkflx) = sbc(i,j,ialkflx) + 5*co2emit/(12.e-7)
!             sbc(i,j,ialkflx) = sbc(i,j,ialkflx)
          endif
         enddo
        enddo
      endif

! July
      if((ttt .lt. 0.521) .and. (ttt .gt. 0.520) .and. (ratio.gt.1.01))
     & then
        do j=2,jmtm1
         do i=2,imtm1
          if (tmsk(i,j) .ge. 0.5) then
            sbc(i,j,ialkflx) = sbc(i,j,ialkflx) + 5*co2emit/(12.e-7)
!             sbc(i,j,ialkflx) = sbc(i,j,ialkflx)
          endif
         enddo
        enddo
      endif

! Aug
      if((ttt .lt. 0.603) .and. (ttt .gt. 0.602) .and. (ratio.gt.1.01))
     & then
        do j=2,jmtm1
         do i=2,imtm1
          if (tmsk(i,j) .ge. 0.5) then
            sbc(i,j,ialkflx) = sbc(i,j,ialkflx) + 5*co2emit/(12.e-7)
!             sbc(i,j,ialkflx) = sbc(i,j,ialkflx)
          endif
         enddo
        enddo
      endif

! Sep
      if((ttt .lt. 0.699) .and. (ttt .gt. 0.698) .and. (ratio.gt.1.01))
     & then
        do j=2,jmtm1
         do i=2,imtm1
          if (tmsk(i,j) .ge. 0.5) then
            sbc(i,j,ialkflx) = sbc(i,j,ialkflx) + 5*co2emit/(12.e-7)
!             sbc(i,j,ialkflx) = sbc(i,j,ialkflx)
          endif
         enddo
        enddo
      endif

! Oct
      if((ttt .lt. 0.781) .and. (ttt .gt. 0.780) .and. (ratio.gt.1.01))
     & then
        do j=2,jmtm1
         do i=2,imtm1
          if (tmsk(i,j) .ge. 0.5) then
            sbc(i,j,ialkflx) = sbc(i,j,ialkflx) + 5*co2emit/(12.e-7)
!             sbc(i,j,ialkflx) = sbc(i,j,ialkflx)
          endif
         enddo
        enddo
      endif

! Nov
      if((ttt .lt. 0.864) .and. (ttt .gt. 0.863) .and. (ratio.gt.1.01))
     & then
        do j=2,jmtm1
         do i=2,imtm1
          if (tmsk(i,j) .ge. 0.5) then
            sbc(i,j,ialkflx) = sbc(i,j,ialkflx) + 5*co2emit/(12.e-7)
!             sbc(i,j,ialkflx) = sbc(i,j,ialkflx)
          endif
         enddo
        enddo
      endif

! Dec
      if((ttt .lt. 0.946) .and. (ttt .gt. 0.945) .and. (ratio.gt.1.01))
     & then
        do j=2,jmtm1
         do i=2,imtm1
          if (tmsk(i,j) .ge. 0.5) then
            sbc(i,j,ialkflx) = sbc(i,j,ialkflx) + 5*co2emit/(12.e-7)
!             sbc(i,j,ialkflx) = sbc(i,j,ialkflx)
          endif
         enddo
        enddo
      endif
#   endif
#   endif
#   if defined O_carbon
      call setbcx (sbc(1,1,idicflx), imt, jmt)
#    if defined O_carbon_14
      call setbcx (sbc(1,1,ic14flx), imt, jmt)
#    endif
#   endif
#   if defined O_npzd_alk
      call setbcx (sbc(1,1,ialkflx), imt, jmt)
#   endif
#   if defined O_npzd_o2
      call setbcx (sbc(1,1,io2flx), imt, jmt)
#   endif
#   if defined O_npzd
      call setbcx (sbc(1,1,ipo4flx), imt, jmt)
#    if !defined O_npzd_no_vflux
      call setbcx (sbc(1,1,iphytflx), imt, jmt)
      call setbcx (sbc(1,1,izoopflx), imt, jmt)
      call setbcx (sbc(1,1,idetrflx), imt, jmt)
#    endif
#    if defined O_npzd_nitrogen
      call setbcx (sbc(1,1,ino3flx), imt, jmt)
#     if !defined O_npzd_no_vflux
      call setbcx (sbc(1,1,idiazflx), imt, jmt)
#     endif
#    endif
#   endif
#   if defined O_cfcs_data || defined O_cfcs_data_transient
      call setbcx (sbc(1,1,icfc11flx), imt, jmt)
      call setbcx (sbc(1,1,icfc12flx), imt, jmt)
#   endif
#  endif
# endif
# if defined O_embm_read_sflx
      call read_var ("F_salt.nc", "F_salt", sbc(1,1,isflx), ntrec)
# elif defined O_embm_write_sflx
      call write_var ("F_salt.nc", "F_salt", sbc(1,1,isflx), ntrec)
# endif
#endif
# if defined O_embm && defined O_sulphate_data || defined O_sulphate_data_transient

!-----------------------------------------------------------------------
!     set anthropogenic sulphate forcing (increase in surface albedo)
!     for the next atmospheric step from updated surface albedo
!-----------------------------------------------------------------------
      call sulphdata

      zero(1) = 0.
      cz(1) = 0.
!     full direct and indirect sulphate forcing
      sulphfac = 0.29*aero_scal
      
!      print*, "aero_scal", aero_scal
!      print*, "sulphfac", sulphfac

      
#  if defined O_sulphate_data_direct && !defined O_sulphate_data_indirect
!     reduce sulphate forcing to just the direct effect (0.4/1.1)
!     updated using the mean values for the 2004-2014 from the Forester 
!     dataset, 05 June 2019 (0.343/1.241) by N. Mengis
      sulphfac = sulphfac*.343/1.241
#  endif
#  if defined O_sulphate_data_indirect && !defined O_sulphate_data_direct
!     reduce sulphate forcing to just the indirect effect (0.7/1.1)
!     updated using the mean values for the 2004-2014 from the Forester 
!     dataset, 05 June 2019 (0.898/1.241) by N. Mengis
      sulphfac = sulphfac*.898/1.241
#  endif
      do j=jsp1,jem1
        do i=isp1,iem1
          if (sulph(i,j,2) .gt. 0. .and. solins(i,j) .gt. 0.) then
            tot = 0.
            twt = 0.
!           calculate sunlight weighted daily average zenith angle
            do k=1,24
                call zenith (1, (k-1)*3600., 3600., daylen, tlat(i,j)
     &,                      zero(1), sindec, cz(1))
              tot = tot + cz(1)*cz(1)
              twt = twt + cz(1)
            enddo
            cz(1) = 0.
            if (twt .gt. 0) cz(1) = tot/twt
!           convert from optical depth to albedo
            if (cz(1) .gt. 0.05) then
!             reonstruct surface coalbedo
              tmp = solins(i,j) - solins(i,j)*(1. - sbc(i,j,iaca))
     &            - solins(i,j)*sbc(i,j,iaca)*scatter
              if (tmp .gt. 0.) then
                tmp = dnswr(i,j)/tmp
                sulph(i,j,2) = sulphfac*sulph(i,j,2)*tmp**2/cz(1)
              else
                sulph(i,j,2) = 0.
              endif
            else
              sulph(i,j,2) = 0.
            endif
          endif
        enddo
      enddo
# endif

      return
      end

#if defined O_embm_read_sflx
      subroutine read_var (fname, vname, var, ntrec)

      implicit none

      character (*) :: fname, vname

      integer ib(10), ic(10), iou, ntrec

      include "size.h"
      include "param.h"
      include "pconst.h"
      include "stdunits.h"
      include "tmngr.h"

      real var(imt,jmt), tmpij(imtm2,jmtm2)

      ntrec = ntrec + 1
      call openfile (trim(fname), iou)
      ib(:) = 1
      ib(3) = ntrec
      ic(:) = 1
      ic(1) = imtm2
      ic(2) = jmtm2
      call getvara (trim(vname), iou, imtm2*jmtm2, ib, ic, tmpij
     &, 1., 0.)
      var(2:imtm1,2:jmtm1) = tmpij(1:imtm2:1:jmtm2)
      embmbc (var)

      return
      end
#endif

#if defined O_embm_write_sflx
      subroutine write_var (fname, vname, var, ntrec)

      implicit none

      character (*) :: fname, vname

      integer ib(10), ic(10), it(10), id_time, id_xt, id_yt, iou, ntrec
      integer nyear, nmonth, nday, nhour, nmin, nsec

      real tmp

      include "size.h"
      include "param.h"
      include "pconst.h"
      include "stdunits.h"
      include "coord.h"
      include "tmngr.h"

      real var(imt,jmt), tmpij(imtm2,jmtm2), tmpi(imtm2), tmpj(jmtm2)

      ntrec = ntrec + 1
      call openfile (trim(fname), iou)
      if (ntrec .eq. 1) then
        call redef (iou)

        call defdim ('time', iou, 0, id_time)
        call defdim ('longitude', iou, imtm2, id_xt)
        call defdim ('latitude', iou, jmtm2, id_yt)

        it(:) = id_time
        call defvar ('time', iou, 1, it, c0, c0, 'T', 'D'
     &,   'time', 'time', 'years since 0-1-1')
        call putatttext (iou, 'time', 'calendar', calendar)

        it(1) = id_xt
        call defvar ('longitude', iou, 1, it, 0., 0., 'X', 'D'
     &,   'longitude', 'longitude', 'degrees_east')
        it(1) = id_yt
        call defvar ('latitude', iou, 1, it, 0., 0., 'Y', 'D'
     &,   'latitude', 'latitude', 'degrees_north')

        it(:) = id_xt
        it(2) = id_yt
        it(3) = id_time
        call defvar (trim(vname), iou, 3, it, -1.e20, 1.e20, ' ', 'D'
     &,   ' ',' ', ' ')

        call enddef (iou)

!-----------------------------------------------------------------------
!       write 1d data (t)
!-----------------------------------------------------------------------
        call rdstmp (stamp, nyear, nmonth, nday, nhour, nmin, nsec)

        ib(:) = 1
        ic(:) = imtm2
        tmpi(1:imtm2) = xt(2:imtm1)
        call putvara ('longitude', iou, imtm2, ib, ic, tmpi, c1, c0)
        ic(:) = jmtm2
        tmpj(1:jmtm2) = xt(2:jmtm1)
        call putvara ('latitude', iou, jmtm2, ib, ic, tmpj, c1, c0)

      endif
      call putvars ('time', iou, ntrec, relyr, c1, c0)
      ib(:) = 1
      ib(3) = ntrec
      ic(:) = 1
      ic(1) = imtm2
      ic(2) = jmtm2
      tmpij(1:imtm2:1:jmtm2) = var(2:imtm1,2:jmtm1)
      call putvara (trim(vname), iou, imtm2*jmtm2, ib, ic, tmpij
     &, 1., 0.)

      return
      end
#endif