! source file: /raid23/csomes/UVic/2.9/c_n_isotopes/last_glacial_experiments/dlgp15/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"

      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.3829382e2
      b10 = -3.47968e1
      b20 = 1.02790e1
      b01 = 2.60065
      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

!-----------------------------------------------------------------------
!         set the incoming short wave
!-----------------------------------------------------------------------
          dnswr(i,j) = solins(i,j)*sbc(i,j,iaca)*pass*sbc(i,j,isca)

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

          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)
!-----------------------------------------------------------------------
      anthro = co2for*alog(co2ccn/280.0)-0.6

      return
      end