! source file: /home/dullman/UVic/PI_LGMalb_tests/CCSMb/updates/fluxes.F

      subroutine fluxes (is, ie, js, je)

!=======================================================================
!     calculate energy and moisture fluxes

!     Note: evaporation and precipitation are in g cm-2 s-1
!           and humidities are in g g-1

!     for Thompson and Warren outgoing radiation (see: Thompson S.J.,
!     and S.G. Warren 'parameterization of outgoing ...'J. Atmos. Sci.,
!     39, 2667-2680, 1982.
!=======================================================================

      implicit none

      integer i, ie, iem1, imax, is, isp1, iter
      integer j, je, jem1, jmax, js, jsp1, maxit, n

      logical track

      real b00, b10, b20, b01, b11, b21, b02, b12, b22, b03, b13, b23
      real delta, df, dt, dultnt, dulwr, dusens, emax, f, fb, ff, fg
      real fh, fl, fm, qair, qlnd, rhrh, sr, scrit, ssh, tair, teff
      real telev, tlnd, tlold, tol, tol2, ultnt, ulwr, usens, wspd
      real vcs, avg_sat, C2K

      include "size.h"
      include "param.h"
      include "pconst.h"
      include "stdunits.h"
      include "cembm.h"
      include "atm.h"
      include "csbc.h"
      include "ice.h"
      include "veg.h"

      real sat(imt,jmt)
      real totelev(imt,jmt)
      real sat_anom(imt,jmt)
      real acalb_adj(imt,jmt)
!      real acalb_trim(imt,jmt)
      
      isp1 = is + 1
      iem1 = ie - 1
      jsp1 = js + 1
      jem1 = je - 1

!-----------------------------------------------------------------------
!     set appropriate constants
!-----------------------------------------------------------------------
      fb = 0.94*rhoatm*cpatm
      maxit = 10
      tol = 0.01
      emax = 0.0
      imax = 0
      jmax = 0
      scrit = 0.75*soilmax
      ff = rhoatm*vlocn
      C2K = 273.15

!     Thomson and Warren constants

      b00 = 2.37944e2 !adjust for clim sens
      b10 = -3.47968e1
      b20 = 1.02790e1
      b01 = 2.627 !adjust for clim sens
      b11 = -1.62064
      b21 = 6.34856e-1
      b02 = 4.40272e-3
      b12 = -2.26092e-2
      b22 = 1.12265e-2
      b03 = -2.05237e-5
      b13 = -9.67e-5
      b23 = 5.62925e-5

      tol2 = tol*2.0

      do j=jsp1,jem1
        do i=isp1,iem1
# if defined O_embm_albslope || defined O_embm_LWslope
          totelev(i,j) = elev(i,j)
     &          + hicel(i,j,2)
          sat(i,j) = at(i,j,2,isat) - totelev(i,j)*rlapse
          sat_anom(i,j) = sat(i,j) - sbc(i,j,isatbase)
!old way of calculating sat_anom, WRONG! NO LAPSE!          
!          sat_anom(i,j) = at(i,j,2,isat) - sbc(i,j,isatbase)
# endif
!-----------------------------------------------------------------------
!         set the incoming short wave
!-----------------------------------------------------------------------
# if defined O_embm_albslope
          acalb_adj(i,j) = (1. - (1. - sbc(i,j,iaca))
     &                     *(sbc(i,j,iaslope)*sat_anom(i,j)+1.))
          acalb_trim(i,j) = max(acalb_adj(i,j),acalb_min)
          acalb_trim(i,j) = min(acalb_trim(i,j),acalb_max)

          dnswr_noalb(i,j) = solins(i,j)*sbc(i,j,iaca)*pass
     &                       *sbc(i,j,isca)
#   if defined O_embm_dust
          dnswr(i,j) = max(0.,(solins(i,j)*acalb_trim(i,j)
     &                 +sbc(i,j,idsw))*pass*sbc(i,j,isca))
          dnswr_cloud(i,j) = dnswr(i,j) - dnswr_noalb(i,j)
#   else
          dnswr(i,j) = solins(i,j)*acalb_trim(i,j)*pass*sbc(i,j,isca)
          dnswr_cloud(i,j) = dnswr(i,j) - dnswr_noalb(i,j)
#   endif

# else 

#   if defined O_embm_dust
          dnswr(i,j) = max(0.,(solins(i,j)*sbc(i,j,iaca)
     &                 +sbc(i,j,idsw))*pass*sbc(i,j,isca))
#   else
          dnswr(i,j) = solins(i,j)*sbc(i,j,iaca)*pass*sbc(i,j,isca)
#   endif

# endif

!-----------------------------------------------------------------------
!         set wind speed and effective elevated air temperature
!-----------------------------------------------------------------------
          wspd = sbc(i,j,iws)
          telev = elev(i,j)
     &          + hicel(i,j,2)
          teff = at(i,j,2,isat)
     &         - telev*rlapse*rf1*exp(max(-1.,-telev/rf2))

!-----------------------------------------------------------------------
!         calculate outgoing longwave radiation
!-----------------------------------------------------------------------
          rhrh = rh(i,j)*rh(i,j)
# if defined O_embm_LWslope
          outlwr(i,j) = 1.0e3*(b00 + b10*rh(i,j) + b20*rhrh
     &                + (b01 + b11*rh(i,j) + b21*rhrh)*teff
     &                + (b02 + b12*rh(i,j) + b22*rhrh)*teff**2
     &                + (b03 + b13*rh(i,j) + b23*rhrh)*teff**3
     &                + (sbc(i,j,ilwslope)*sat_anom(i,j)))
     &                - anthro
#   if defined O_embm_dust
     &                - sbc(i,j,idlw)
#   endif
          outlwr_cloud(i,j)=1.0e3*(sbc(i,j,ilwslope)*sat_anom(i,j))
# else
          outlwr(i,j) = 1.0e3*(b00 + b10*rh(i,j) + b20*rhrh
     &                + (b01 + b11*rh(i,j) + b21*rhrh)*teff
     &                + (b02 + b12*rh(i,j) + b22*rhrh)*teff**2
     &                + (b03 + b13*rh(i,j) + b23*rhrh)*teff**3)
     &                - anthro
#   if defined O_embm_dust
     &                - sbc(i,j,idlw)            
#   endif

# endif
          tair = at(i,j,2,isat) - telev*rlapse

          if (tmsk(i,j) .ge. 0.5) then

!-----------------------------------------------------------------------
!           calculations only for ocean points
!-----------------------------------------------------------------------
            dt = sbc(i,j,isst) - tair
            fg = dalt_o*wspd

!-----------------------------------------------------------------------
!           calculate evaporation or sublimation (ensure it is positive)
!-----------------------------------------------------------------------
            ssh = cssh*exp(17.67*sbc(i,j,isst)/(sbc(i,j,isst) + 243.5))
            evap(i,j) = max(c0, rhoatm*fg*(ssh - at(i,j,2,ishum)))
            upltnt(i,j) = vlocn*evap(i,j)

!-----------------------------------------------------------------------
!           calculate upward sensible heat flux
!-----------------------------------------------------------------------
            upsens(i,j) = fb*fg*(dt)

!-----------------------------------------------------------------------
!           calculate upward longwave re-radiation
!-----------------------------------------------------------------------
            uplwr(i,j) = esocn*(sbc(i,j,isst) + C2K)**4
     &                 - esatm*(tair + C2K)**4

          elseif (land_map(i,j) .ne. 0) then

!----------------------------------------------------------------------
!           set fluxes over land from the land model
!---------------------------------------------------------------------
            upltnt(i,j) = 0.0
            evap(i,j) = sbc(i,j,ievap)
            upsens(i,j) = sbc(i,j,isens)
            uplwr(i,j) = sbc(i,j,ilwr)

          else

!-----------------------------------------------------------------------
!            calculations only for land points

!           find land temperature by balancing the surface heat budget
!             dwsr = ultnt + usens + ulwr
!           using Newton's method:
!             t(i+1) = t(i) - f(t(i))/df(t(i))
!           where:
!             f(t(i)) = dwsr - ultnt - usens - ulwr
!             -df(t(i)) = dultnt - dusens - dulwr
!-----------------------------------------------------------------------
            tlnd = surf(i,j)
            tlold = tlnd
            fm = esatm*(tair + C2K)**4
            fg = rhoatm
!           calculate stomatal resistance
            sr = (1.-agric(i,j,2))*veg_rs(iveg(i,j))
     &         + agric(i,j,2)*veg_rs(iagric)
            dalt_v = veg_dalt(i,j)
!           add in aerodynamic resistance
            sr = sr + 1.0/(dalt_v*wspd + epsln)
!           set beta parameter for calculating actual evaporation
            fh = min(max(c0+epsln, (soilm(i,j,lf)/soilmax)**(0.25)),c1)
!           set coefficients for latent heat (fl) and evaporation (fg)
            fl = fh*ff/(sr)
            fg = fh*fg/(sr)
            dusens = fb*dalt_v*wspd

!-----------------------------------------------------------------------
!           start loop for all land grid points
!-----------------------------------------------------------------------
            qair = rh(i,j)*cssh*exp(17.67*tair/(tair + 243.5))
            iter = 0
            delta = tol2
            do while (abs(delta) .gt. tol .and. iter .le. maxit)
              iter = iter + 1
              qlnd = cssh*exp(17.67*tlnd/(tlnd + 243.5))
              if (qlnd .gt. qair) then
                ultnt = fl*(qlnd - qair)
                dultnt = fl*qlnd*17.67*243.5/(tlnd + 243.5)**2
              else
                ultnt = 0.0
                dultnt = 0.0
              endif
              usens = dusens*(tlnd - tair)
              ulwr = eslnd*(tlnd + C2K)**4 - fm
              dulwr = 4.0*eslnd*(tlnd + C2K)**3
              f = dnswr(i,j) - ultnt - usens - ulwr
              df = dultnt + dusens + dulwr
              delta = f/df
              tlnd = tlnd + delta
            enddo
            if (iter .gt. maxit) then
!             if not converged, set to last converged temperature
              if (abs(delta) .gt. emax) then
                emax = abs(delta)
                imax = i
                jmax = j
                tlnd = tlold
              endif
            endif

!-----------------------------------------------------------------------
!           calculate fluxes on land
!-----------------------------------------------------------------------
            surf(i,j) = tlnd
            qlnd = cssh*exp(17.67*tlnd/(tlnd + 243.5))
            evap(i,j) = max(c0, fg*(qlnd - qair))
            evap(i,j) = max(c0, min(soilm(i,j,lf)/dts, evap(i,j)))
            flux(i,j,ishum) = evap(i,j)*(1.0 - aice(i,j,2))
            upltnt(i,j) = vlocn*evap(i,j)
            upsens(i,j) = dusens*(tlnd - tair)
            uplwr(i,j) = eslnd*(tlnd + C2K)**4 - fm

!           ensure fluxes are balanced since land can't absorb error
            upsens(i,j) = dnswr(i,j) - upltnt(i,j) - uplwr(i,j)

          endif

        enddo
      enddo

      if (emax .gt. 0.0) write (stdout,*)
     &  '==> Warning: land surface temperature not converging: '
     &, 'emax, i, j, soilm:', emax, imax, jmax, soilm(imax,jmax,2)

      return
      end

      subroutine precipitate (is, ie, js, je)

!=======================================================================
!     calculate precipitation explicitly and update humidity
!=======================================================================

      implicit none

      integer i, ie, iem1, is, isp1, j, je, jem1, js, jsp1, k, n, negq

      real fb, fc, qmax, rate, tair, teff, telev, soiltemp, pson, psot
      real ssh, tmp

      include "size.h"
      include "param.h"
      include "pconst.h"
      include "stdunits.h"
      include "scalar.h"
      include "cembm.h"
      include "atm.h"
      include "switch.h"
      include "csbc.h"
      include "ice.h"
      include "mtlm.h"
      real hs(imt,jmt)

      data negq /0/
      save negq

      if (eoyear) negq = 0

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

!-----------------------------------------------------------------------
!     set appropriate constants
!-----------------------------------------------------------------------
      fb = rhoatm*shq/dts
      fc = dts/rhosno
!     maximum relative humidity after rain
      call unloadland (POINTS, LYING_SNOW, imt, jmt, land_map, hs)
!     convert from kg/m2 to cm
      hs(:,:) = hs(:,:)*0.1/rhosno

      precip(:,:) = 0.
      rh(:,:) = 0.
      k = 0
      do j=jsp1,jem1
        do i=isp1,iem1

!-----------------------------------------------------------------------
!         check if specific humidity is greater than rhmax of saturation
!-----------------------------------------------------------------------
          telev = elev(i,j)
     &          + hicel(i,j,2)
          teff = at(i,j,2,isat)
     &         - telev*rlapse*rf1*exp(max(-1.,-telev/rf2))
          ssh = cssh*exp(17.67*teff/(teff + 243.5))
          qmax = rhmax*ssh
          if (at(i,j,2,ishum) .gt. qmax) then
            tmp = fb*(at(i,j,2,ishum) - qmax)
            precip(i,j) = precip(i,j) + tmp
            at(i,j,2,ishum) = at(i,j,2,ishum) - tmp/fb
            rh(i,j) = rhmax
          endif
          rh(i,j) = at(i,j,2,ishum)/(ssh + epsln)
          rh(i,j) = max(c0, min(c1, rh(i,j)))

!-----------------------------------------------------------------------
!         calculate snowfall (hsno at tau was set in the ice model)
!-----------------------------------------------------------------------
!         tair may be adjusted by a snowfall offset temperature tsno
          tair = at(i,j,2,isat) - tsno - telev*rlapse
          psno(i,j) = 0.0
          if (land_map(i,j) .eq. 0) hs(i,j) = hsno(i,j,2)
          if (tair .le. c0 .and. hs(i,j) .lt. hsno_max)
     &      psno(i,j) = min((hsno_max - hs(i,j))/fc, precip(i,j))
          if (tmsk(i,j) .ge. 0.5) then
!           only allow snow where there is sea ice
            psno(i,j) = psno(i,j)*aice(i,j,2)
            hsno(i,j,2) = hsno(i,j,2) + fc*psno(i,j)
            if (addflxa) flux(i,j,ishum) = flux(i,j,ishum)
     &                                   - dts*psno(i,j)
          elseif (land_map(i,j) .eq. 0) then
            hsno(i,j,2) = hsno(i,j,2) + fc*psno(i,j)
          endif

!-----------------------------------------------------------------------
!         update soilm and allocate surplus soil moisture to runoff
!-----------------------------------------------------------------------
          if (tmsk(i,j) .lt. 0.5 .and. land_map(i,j) .eq. 0) then
            flux(i,j,ishum) = flux(i,j,ishum) - precip(i,j) + psno(i,j)
            soiltemp = soilm(i,j,2)
            soilm(i,j,2) = soilm(i,j,lf) - dts*flux(i,j,ishum)
            soilm(i,j,2) = max(c0, soilm(i,j,2))
            soilm(i,j,1) = soiltemp
            if (soilm(i,j,2) .gt. soilmax) then
              runoff(i,j) = (soilm(i,j,2) - soilmax)/dts
              soilm(i,j,2) = soilmax
            else
              runoff(i,j) = c0
            endif
          endif
        enddo
      enddo

      call embmbc (psno)
      call embmbc (hsno(1,1,2))

      return
      end

      subroutine co2forc
!=======================================================================
!     calculate global average CO2 forcing
!=======================================================================

      implicit none

      real yr

      include "cembm.h"

!-----------------------------------------------------------------------
!     relative forcing from 280 ppmv (anthro)
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
!     relative forcing (anthro) is added to current forcing
!-----------------------------------------------------------------------
cAAA CH4 and N2O LGM forcing added (-0.577 W/m^2), calculated from equations in IPCC report, using PMIP/CMIP concentrations, included for LGM by turning on the O_GHG_LGM flag in the mk.in file
      anthro = co2for*alog(co2ccn/280.0)
# if O_GHG_LGM
     &         - 0.577e3   
# endif

      return
      end
