#include <misc.h>#include <preproc.h>subroutine Biogeophysics_Lake (clm) !-----------------------------------------------------------------------!!  CLMCLMCLMCLMCLMCLMCLMCLMCLMCL  A community developed and sponsored, freely!  L                           M  available land surface process model.!  M --COMMUNITY LAND MODEL--  C!  C                           L!  LMCLMCLMCLMCLMCLMCLMCLMCLMCLM!!-----------------------------------------------------------------------! Purpose:! Calculates lake temperatures and surface fluxes.!! Method:! Lake temperatures are determined from a one-dimensional thermal! stratification model based on eddy diffusion concepts to ! represent vertical mixing of heat.!! d ts    d            d ts     1 ds! ---- = -- [(km + ke) ----] + -- --!  dt    dz             dz     cw dz   !! where: ts = temperature (kelvin)!         t = time (s)!         z = depth (m)!        km = molecular diffusion coefficient (m**2/s)!        ke = eddy diffusion coefficient (m**2/s)!        cw = heat capacity (j/m**3/kelvin)!         s = heat source term (w/m**2)!! There are two types of lakes: !   Deep lakes are 50 m. !   Shallow lakes are 10 m deep.!!   For unfrozen deep lakes:    ke > 0 and    convective mixing!   For unfrozen shallow lakes: ke = 0 and no convective mixing!! Use the Crank-Nicholson method to set up tridiagonal system of equations to! solve for ts at time n+1, where the temperature equation for layer i is! r_i = a_i [ts_i-1] n+1 + b_i [ts_i] n+1 + c_i [ts_i+1] n+1!! The solution conserves energy as:!! cw*([ts(      1)] n+1 - [ts(      1)] n)*dz(      1)/dt + ... +! cw*([ts(nlevlak)] n+1 - [ts(nlevlak)] n)*dz(nlevlak)/dt = fin!! where:! [ts] n   = old temperature (kelvin)! [ts] n+1 = new temperature (kelvin)! fin      = heat flux into lake (w/m**2)!          = beta*sabg + forc_lwrad - eflx_lwrad_out - eflx_sh_tot - eflx_lh_tot!            - hm + phi(1) + ... + phi(nlevlak) !! Author:! Gordon Bonan! 15 September 1999: Yongjiu Dai; Initial code! 15 December 1999:  Paul Houser and Jon Radakovich; F90 Revision ! April 2002: Vertenstein/Oleson/Levis; Final form!!-----------------------------------------------------------------------! $Id: Biogeophysics_Lake.F90,v 1.6.6.2 2002/04/27 15:38:36 erik Exp $!-----------------------------------------------------------------------  use precision  use clmtype  use clm_varpar, only : nlevlak  use clm_varcon, only : hvap, hsub, hfus, cpair, cpliq, tkwat, tkice, &                        sb, vkc, grav, denh2o, tfrz, spval  implicit none!----Arguments----------------------------------------------------------  type (clm1d), intent(inout) :: clm      ! CLM 1-D Module!----Local Variables----------------------------------------------------  integer i,j          ! do loop or array index  integer :: idlak = 1 ! index of lake, 1 = deep lake, 2 = shallow lake  integer niters       ! maximum number of iterations for surface temperature  integer iter         ! iteration index  integer nmozsgn      ! number of times moz changes sign  real(r8) ax      ! temporary variable  real(r8) bx      ! temporary variable  real(r8) beta1   ! coefficient of convective velocity [-]  real(r8) degdT   ! d(eg)/dT  real(r8) dqh     ! diff of humidity between ref. height and surface  real(r8) dth     ! diff of virtual temp. between ref. height and surface  real(r8) dthv    ! diff of vir. poten. temp. between ref. height and surface  real(r8) dzsur   ! thickness of snow cover plus top lake layer [m]  real(r8) eg      ! water vapor pressure at temperature T [pa]  real(r8) emg     ! ground emissivity  real(r8) hm      ! energy residual [W/m2]  real(r8) htvp    ! latent heat of vaporization (or sublimation) [J/kg]  real(r8) obu     ! monin-obukhov length [m]  real(r8) obuold  ! monin-obukhov length of previous iteration [m]  real(r8) qsatg   ! saturated humidity [kg/kg]  real(r8) qsatgdT ! d(qsatg)/dT  real(r8) qstar   ! moisture scaling parameter  real(r8) ram     ! aerodynamic resistance [s/m]  real(r8) rah     ! thermal resistance [s/m]  real(r8) raw     ! moisture resistance [s/m]  real(r8) stftg3  ! temporary variable  real(r8) temp1   ! relation for potential temperature profile  real(r8) temp2   ! relation for specific humidity profile  real(r8) tgbef   ! lake surface temperature from previous iteration [K]  real(r8) thm     ! intermediate variable (forc_t+0.0098*forc_hgt_t) [K]  real(r8) thv     ! virtual potential temperature [K]  real(r8) thvstar ! virtual potential temperature scaling parameter [K]  real(r8) tksur   ! thermal conductivity of snow/soil [W/m/K]  real(r8) tstar   ! temperature scaling parameter [K]  real(r8) um      ! wind speed including the stability effect [m/s]  real(r8) ur      ! wind speed at reference height [m/s]  real(r8) ustar   ! friction velocity [m/s]  real(r8) wc      ! convective velocity [m/s]  real(r8) zeta    ! dimensionless height used in Monin-Obukhov theory  real(r8) zii     ! convective boundary height [m]  real(r8) zldis   ! reference height "minus" zero displacement height [m]  real(r8) displa  ! displacement height [m]  real(r8) z0mg    ! roughness length over ground, momentum [m]  real(r8) z0hg    ! roughness length over ground, sensible heat [m]  real(r8) z0qg    ! roughness length over ground, latent heat [m]  real(r8) beta(2) ! fraction solar rad absorbed at surface:                   !  depends on lake type  real(r8) za(2)   ! base of surface absorption layer (m): depends on lake type  real(r8) eta(2)  ! light extinction coefficient (m-1): depends on lake type  real(r8) p0      ! neutral value of turbulent prandtl number  real(r8) a(nlevlak)    ! "a" vector for tridiagonal matrix  real(r8) b(nlevlak)    ! "b" vector for tridiagonal matrix  real(r8) c(nlevlak)    ! "c" vector for tridiagonal matrix  real(r8) r(nlevlak)    ! "r" vector for tridiagonal solution  real(r8) rhow(nlevlak) ! density of water [kg/m3]  real(r8) phi(nlevlak)  ! solar radiation absorbed by layer [W/m2]  real(r8) kme(nlevlak)  ! molecular + eddy diffusion coefficient [m2/s]  real(r8) cwat    ! specific heat capacity of water [J/m3/K]  real(r8) ws      ! surface friction velocity [m/s]  real(r8) ks      ! coefficient  real(r8) in      ! relative flux of solar radiation into layer  real(r8) out     ! relative flux of solar radiation out of layer  real(r8) ri      ! richardson number  real(r8) fin     ! heat flux into lake - flux out of lake [W/m2]  real(r8) ocvts   ! (cwat*(t_lake[n  ])*dz  real(r8) ncvts   ! (cwat*(t_lake[n+1])*dz  real(r8) m1      ! intermediate variable for calculating r, a, b, c  real(r8) m2      ! intermediate variable for calculating r, a, b, c  real(r8) m3      ! intermediate variable for calculating r, a, b, c  real(r8) ke      ! eddy diffusion coefficient [m2/s]  real(r8) km      ! molecular diffusion coefficient [m2/s]  real(r8) zin     ! depth at top of layer [m]  real(r8) zout    ! depth at bottom of layer [m]  real(r8) drhodz  ! d [rhow] /dz [kg/m4]  real(r8) n2      ! brunt-vaisala frequency [s-2]  real(r8) num     ! used in calculating ri  real(r8) den     ! used in calculating ri  real(r8) tav     ! used in aver temp for convectively mixed layers  real(r8) nav     ! used in aver temp for convectively mixed layers  real(r8) phidum  ! temporary value of phi  real(r8) u2m     ! 2 m wind speed [m/s]  real(r8) cf      ! s m2/umol -> s/m!----End Variable List--------------------------------------------------!! Surface Radiation!  call SurfaceRadiation (clm)!! Determine beginning water balance!  clm%begwb = clm%h2osno  !! [1] Constants and model parameters!!! Constants for lake temperature model!  beta = (/0.4, 0.4/)                            ! (deep lake, shallow lake)  za   = (/0.6, 0.5/)      eta  = (/0.1, 0.5/)    p0   = 1.  !! Roughness lengths!  if (clm%t_grnd >= tfrz)then                    ! for unfrozen lake     z0mg = 0.01  else                                           ! for frozen lake     z0mg = 0.04  endif  z0hg = z0mg  z0qg = z0mg!! Latent heat !  if (clm%forc_t > tfrz) then     htvp = hvap  else			     htvp = hsub  endif!! Switch between vaporization and sublimation causes rapid solution! separation in perturbation growth test!#if (defined PERGRO)  htvp = hvap#endif!! Emissivity!  emg = 0.97!! [2] Surface temperature and fluxes!  dzsur = clm%dz(1) + clm%snowdp!! Saturated vapor pressure, specific humidity and their derivatives! at lake surface!  call QSat(clm%t_grnd, clm%forc_pbot, eg, degdT, qsatg, &            qsatgdT     )!! Potential, virtual potential temperature, and wind speed at the! reference height!  beta1=1.       ! -  (in computing W_*)  zii = 1000.    ! m  (pbl height)  thm = clm%forc_t + 0.0098*clm%forc_hgt_t         ! intermediate variable   thv = clm%forc_th*(1.+0.61*clm%forc_q)           ! virtual potential T  ur = max(1.0_r8,sqrt(clm%forc_u*clm%forc_u+clm%forc_v*clm%forc_v))!! Initialize stability variables!  nmozsgn = 0  obuold = 0.  dth   = thm-clm%t_grnd  dqh   = clm%forc_q-qsatg  dthv  = dth*(1.+0.61*clm%forc_q)+0.61*clm%forc_th*dqh  zldis = clm%forc_hgt_u-0.!! Initialize Monin-Obukhov length and wind speed including stability effect!  call MoninObukIni(ur, thv, dthv, zldis, z0mg, &                    um, obu  )!! Begin stability iteration!  niters = 3  do iter = 1, niters     tgbef = clm%t_grnd     if (clm%t_grnd > tfrz) then        tksur = tkwat     else        tksur = tkice     endif!! Determine friction velocity, and potential temperature and humidity! profiles of the surface boundary layer!     displa = 0.0_r8     call FrictionVelocity (displa, z0mg,  z0hg,  z0qg,  obu, &                            iter, ur, um, ustar, temp1, temp2, clm)     obuold = obu!! Determine aerodynamic resistances