biogeophysics_lake.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

!

     ram    = 1./(ustar*ustar/um)
     rah    = 1./(temp1*ustar)
     raw    = 1./(temp2*ustar)

!
! Get derivative of fluxes with respect to ground temperature
!

     stftg3 = emg*sb*tgbef*tgbef*tgbef

     ax  = clm%sabg + emg*clm%forc_lwrad + 3.*stftg3*tgbef &
          + clm%forc_rho*cpair/rah*thm &
          - htvp*clm%forc_rho/raw*(qsatg-qsatgdT*tgbef - clm%forc_q) &
          + tksur*clm%t_lake(1)/dzsur

     bx  = 4.*stftg3 + clm%forc_rho*cpair/rah &
          + htvp*clm%forc_rho/raw*qsatgdT + tksur/dzsur

     clm%t_grnd = ax/bx

!
! Surface fluxes of momentum, sensible and latent heat
! using lake surface temperature from previous time step
!

     clm%eflx_sh_grnd = clm%forc_rho*cpair*(clm%t_grnd-thm)/rah
     clm%qflx_evap_soi = clm%forc_rho*(qsatg+qsatgdT*(clm%t_grnd-tgbef)-clm%forc_q)/raw

!
! Re-calculate saturated vapor pressure, specific humidity and their
! derivatives at lake surface
!

     call QSat(clm%t_grnd, clm%forc_pbot, eg, degdT, qsatg, &
               qsatgdT     )

     dth=thm-clm%t_grnd
     dqh=clm%forc_q-qsatg

     tstar = temp1*dth
     qstar = temp2*dqh

     dthv=dth*(1.+0.61*clm%forc_q)+0.61*clm%forc_th*dqh
     thvstar=tstar*(1.+0.61*clm%forc_q) + 0.61*clm%forc_th*qstar
     zeta=zldis*vkc * grav*thvstar/(ustar**2*thv)

     if (zeta >= 0.) then     !stable
        zeta = min(2._r8,max(zeta,0.01_r8))
        um = max(ur,0.1_r8)
     else                     !unstable
        zeta = max(-100._r8,min(zeta,-0.01_r8))
        wc = beta1*(-grav*ustar*thvstar*zii/thv)**0.333
        um = sqrt(ur*ur+wc*wc)
     endif
     obu = zldis/zeta

     if (obuold*obu < 0.) nmozsgn = nmozsgn+1
     if (nmozsgn >= 4) EXIT

  enddo

!
! If there is snow on the ground and t_grnd > tfrz: reset t_grnd = tfrz.
! Re-evaluate ground fluxes. Energy inbalance used to melt snow.  
! h2osno > 0.5 prevents spurious fluxes.
!

  if (clm%h2osno > 0.5 .AND. clm%t_grnd > tfrz) then
     clm%t_grnd = tfrz
     clm%eflx_sh_grnd = clm%forc_rho*cpair*(clm%t_grnd-thm)/rah
     clm%qflx_evap_soi = clm%forc_rho*(qsatg+qsatgdT*(clm%t_grnd-tgbef) & 
                       - clm%forc_q)/raw  !note: qsatg and qsatgdT should be f(tgbef)
  endif

!
! Net longwave from ground to atmosphere
!

  clm%eflx_lwrad_out = (1.-emg)*clm%forc_lwrad + stftg3*(-3.*tgbef+4.*clm%t_grnd)

!
! Radiative temperature
!

  clm%t_rad = (clm%eflx_lwrad_out/sb)**0.25

!
! Ground heat flux
!

  clm%eflx_soil_grnd = clm%sabg + clm%forc_lwrad - clm%eflx_lwrad_out - &
                       clm%eflx_sh_grnd - htvp*clm%qflx_evap_soi

  clm%taux   = -clm%forc_rho*clm%forc_u/ram
  clm%tauy   = -clm%forc_rho*clm%forc_v/ram

  clm%eflx_sh_tot   = clm%eflx_sh_grnd
  clm%qflx_evap_tot = clm%qflx_evap_soi
  clm%eflx_lh_tot   = htvp*clm%qflx_evap_soi
  clm%eflx_lh_grnd  = htvp*clm%qflx_evap_soi

!
! 2 m height air temperature
!

  clm%t_ref2m   = (clm%t_grnd + temp1*dth * 1./ &
                  vkc *log((2.+z0hg)/z0hg))

!
! Energy residual used for melting snow
!

  if (clm%h2osno > 0. .AND. clm%t_grnd >= tfrz) then
     hm = min( clm%h2osno*hfus/clm%dtime, max(clm%eflx_soil_grnd,0._r8) )
  else
     hm = 0.
  endif
  clm%qmelt = hm/hfus             ! snow melt (mm/s)

!
! [3] Lake layer temperature
!

!
! Lake density
!

  do j = 1, nlevlak
     rhow(j) = 1000.*( 1.0 - 1.9549e-05*(abs(clm%t_lake(j)-277.))**1.68 )
  enddo

!
! Eddy diffusion + molecular diffusion coefficient:
! eddy diffusion coefficient used for unfrozen deep lakes only
!

  cwat = cpliq*denh2o
  km = tkwat/cwat

  fin = beta(idlak) * clm%sabg + clm%forc_lwrad - (clm%eflx_lwrad_out + &
        clm%eflx_sh_tot + clm%eflx_lh_tot + hm)
  u2m = max(1.0_r8,ustar/vkc*log(2./z0mg))

  ws = 1.2e-03 * u2m
  ks = 6.6*sqrt(abs(sin(clm%lat)))*(u2m**(-1.84))

  do j = 1, nlevlak-1
     drhodz = (rhow(j+1)-rhow(j)) / (clm%z(j+1)-clm%z(j))
     n2 = -grav / rhow(j) * drhodz
     num = 40. * n2 * (vkc*clm%z(j))**2
     den = max( (ws**2) * exp(-2.*ks*clm%z(j)), 1.e-10_r8 )
     ri = ( -1. + sqrt( max(1.+num/den, 0._r8) ) ) / 20.
     if (idlak == 1 .AND. clm%t_grnd > tfrz) then
        ke = vkc*ws*clm%z(j)/p0 * exp(-ks*clm%z(j)) / (1.+37.*ri*ri)
     else
        ke = 0.
     endif
     kme(j) = km + ke 
  enddo

  kme(nlevlak) = kme(nlevlak-1)

!
! Heat source term: unfrozen lakes only
!

  do j = 1, nlevlak
     zin  = clm%z(j) - 0.5*clm%dz(j)
     zout = clm%z(j) + 0.5*clm%dz(j)
     in  = exp( -eta(idlak)*max(  zin-za(idlak),0._r8 ) )
     out = exp( -eta(idlak)*max( zout-za(idlak),0._r8 ) )

!
! Assume solar absorption is only in the considered depth
!

     if (j == nlevlak) out = 0.  
     if (clm%t_grnd > tfrz) then
        phidum = (in-out) * clm%sabg * (1.-beta(idlak))
     else if (j == 1) then
        phidum= clm%sabg * (1.-beta(idlak))
     else
        phidum = 0.
     endif
     phi(j) = phidum
  enddo

!
! Sum cwat*t_lake*dz for energy check
!

  ocvts = 0.
  do j = 1, nlevlak
     ocvts = ocvts + cwat*clm%t_lake(j)*clm%dz(j) 
  enddo

!
! Set up vector r and vectors a, b, c that define tridiagonal matrix
!

  j = 1
  m2 = clm%dz(j)/kme(j) + clm%dz(j+1)/kme(j+1)
  m3 = clm%dtime/clm%dz(j)
  r(j) = clm%t_lake(j) + (fin+phi(j))*m3/cwat - (clm%t_lake(j)-clm%t_lake(j+1))*m3/m2
  a(j) = 0.
  b(j) = 1. + m3/m2
  c(j) = -m3/m2

  j = nlevlak
  m1 = clm%dz(j-1)/kme(j-1) + clm%dz(j)/kme(j)
  m3 = clm%dtime/clm%dz(j)
  r(j) = clm%t_lake(j) + phi(j)*m3/cwat + (clm%t_lake(j-1)-clm%t_lake(j))*m3/m1
  a(j) = -m3/m1
  b(j) = 1. + m3/m1
  c(j) = 0.

  do j = 2, nlevlak-1
     m1 = clm%dz(j-1)/kme(j-1) + clm%dz(j  )/kme(j  )
     m2 = clm%dz(j  )/kme(j  ) + clm%dz(j+1)/kme(j+1)
     m3 = clm%dtime/clm%dz(j)
     r(j) = clm%t_lake(j) + phi(j)*m3/cwat + &
            (clm%t_lake(j-1) - clm%t_lake(j))*m3/m1 - &
            (clm%t_lake(j)-clm%t_lake(j+1))*m3/m2

     a(j) = -m3/m1
     b(j) = 1. + m3/m1 + m3/m2
     c(j) = -m3/m2
  enddo

!
! Solve for t_lake: a, b, c, r, u 
!

  call Tridiagonal (nlevlak, a, b, c, r, clm%t_lake(1:nlevlak)) 

!
! Convective mixing: make sure cwat*dz*ts is conserved.
!

  if (idlak == 1 .AND. clm%t_grnd > tfrz) then
     do j = 1, nlevlak-1
        if (rhow(j) > rhow(j+1)) then
           tav = 0.
           nav = 0.
           do i = 1, j+1
              tav = tav + clm%t_lake(i)*clm%dz(i)
              nav = nav + clm%dz(i)
           enddo
           tav = tav/nav
           do i = 1, j+1
              clm%t_lake(i) = tav
              rhow(i) = 1000.*( 1.0 - 1.9549e-05*(abs(clm%t_lake(i)-277.))**1.68 )
           enddo
        endif
     enddo
  endif
 
!
! Sum cwat*t_lake*dz and total energy into lake for energy check
!

  ncvts = 0.
  do j = 1, nlevlak
     ncvts = ncvts + cwat*clm%t_lake(j)*clm%dz(j) 
     fin = fin + phi(j)
  enddo

  clm%errsoi = (ncvts-ocvts) / clm%dtime - fin

!
! [4] Set other clm values for lake points
!

!
! The following are needed for global average on history tape.
! Note: time invariant variables set in initialization phase:
! z, dz, snl, h2osoi_liq, and h2osoi_ice
!

  clm%t_veg = clm%forc_t  ! to be consistent with treatment of t_veg for bare soil points

  clm%eflx_sh_veg     = 0.
  clm%eflx_lh_vegt    = 0.
  clm%eflx_lh_vege    = 0.
  clm%eflx_lwrad_net  = clm%eflx_lwrad_out -  clm%forc_lwrad
  
! Components that are not displayed over lake on history tape and 
! therefore need to be set to spval here

  clm%rssun    = spval
  clm%rssha    = spval
  clm%t_snow   = spval

end subroutine Biogeophysics_Lake