!     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   = spvalend subroutine Biogeophysics_Lake