! source file: /sfs/fs6/home-geomar/smomw258/UVOK_1.1/Kiel_Feb_2019/source/embm/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)
return
end