lake.f90

CLM集合卡曼滤波数据同化算法

! Evaluated stability-dependent variables using moz from prior iteration
         call obult(hu,ht,hq,0.,z0mg,z0hg,z0qg,obu,um,ustar,temp1,temp2,temp12m)

         obuold = obu
!
! Get derivative of fluxes with repect to ground temperature

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

         stftg3 = emg*stefnc*tgbef*tgbef*tgbef

         ax  = sabg + emg*frl + 3.*stftg3*tgbef &
             + rhoair*cp/rah*thm &
             - htvp*rhoair/raw*(qsatg-qsatgdT*tgbef - qm) &
             + tksur*tlak(1)/dzsur
 
         bx  = 4.*stftg3 + rhoair*cp/rah &
             + htvp*rhoair/raw*qsatgdT + tksur/dzsur
 
         tg = ax/bx

! surface fluxes of momentum, sensible and latent
! using ground temperatures from previous time step

         fseng = rhoair*cp*(tg-thm)/rah
         fevpg = rhoair*(qsatg+qsatgdT*(tg-tgbef)-qm)/raw
 
         call qsadv(tg,psrf,eg,degdT,qsatg,qsatgdT)
         dth=thm-tg
         dqh=qm-qsatg

         tstar = temp1*dth
         qstar = temp2*dqh

         thvstar=tstar+0.61*th*qstar
         zeta=zldis*vonkar*grv*thvstar/(ustar**2*thv)
         if(zeta >= 0.) then     !stable
           zeta = min(1.,max(zeta,1.e-6))
         else                    !unstable
           zeta = max(-100.,min(zeta,-1.e-6))
         endif
         obu = zldis/zeta

         if(zeta >= 0.)then
           um = max(ur,0.1)
         else
           wc = beta1*(-grv*ustar*thvstar*zii/thv)**0.333
           um = sqrt(ur*ur+wc*wc)
         endif

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

      enddo
! ----------------------------------------------------------------------

! if snow on ground and tg > tfrz: reset tg = tfrz. reevaluate ground fluxes.
! energy inbalance used to melt snow. scv > 0.5 prevents spurious fluxes

      if (scv > 0.5 .AND. tg > tfrz) then
         tg = tfrz
         fseng = rhoair*cp*(tg-thm)/rah
         fevpg = rhoair*(qsatg+qsatgdT*(tg-tgbef)-qm)/raw !*qsatg and qsatgdT
                                                          !*should be f(tgbef)
      end if

! net longwave from ground to atmosphere
      olrg = (1.-emg)*frl + stftg3*(-3.*tgbef+4.*tg)

! radiative temperature
      trad = (olrg/stefnc)**0.25

! ground heat flux
      fgrnd = sabg + frl - olrg - fseng - htvp*fevpg

      taux   = -rhoair*us/ram
      tauy   = -rhoair*vs/ram

      fsena  = fseng
      fevpa  = fevpg
      lfevpa = htvp*fevpg

! 2 m height air temperature
!     tref   = (tg + temp1*dth * 1./vonkar *alog((2.+z0hg)/z0hg))
      tref   = thm + temp1*dth * (1./temp12m - 1./temp1)

! energy residual for snow melting
      if (scv > 0. .AND. tg >= tfrz) then
         hm = min( scv*dlm/dtime, max(fgrnd,0.) )
      else
         hm = 0.
      end if
      qmelt = hm/dlm             ! snow melt (mm/s)

! ----------------------------------------------------------------------
!*[3] LAKE LAYER TEMPERATURE
! ----------------------------------------------------------------------

! lake density

      do j = 1, msl
         rhow(j) = 1000.*( 1.0 - 1.9549e-05*(abs(tlak(j)-277.))**1.68 )
      end do

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

      cwat = cl*rhowat
      km = tkwat/cwat

      fin = beta(idlak)*sabg + frl - (olrg+fsena+lfevpa+hm)
      u2m = max(1.0,ustar/vonkar*log(2./z0mg))

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

      do j = 1, msl-1
         drhodz = (rhow(j+1)-rhow(j)) / (zlak(j+1)-zlak(j))
         n2 = -grv / rhow(j) * drhodz
         num = 40. * n2 * (vonkar*zlak(j))**2
         den = max( (ws**2) * exp(-2.*ks*zlak(j)), 1.e-10 )
         ri = ( -1. + sqrt( max(1.+num/den, 0.) ) ) / 20.
         if (idlak == 1 .AND. tg > tfrz) then
            ke = vonkar*ws*zlak(j)/p0 * exp(-ks*zlak(j)) / (1.+37.*ri*ri)
         else
            ke = 0.
         end if
         kme(j) = km + ke 
      end do

      kme(msl) = kme(msl-1)

! heat source term: unfrozen lakes only

      do j = 1, msl
         zin  = zlak(j) - 0.5*dzlak(j)
         zout = zlak(j) + 0.5*dzlak(j)
         in  = exp( -eta(idlak)*max(  zin-za(idlak),0. ) )
         out = exp( -eta(idlak)*max( zout-za(idlak),0. ) )
         !assumed solar absorption is only in the considered depth
         if(j == msl) out = 0.  
         if (tg > tfrz) then
            phidum = (in-out) * sabg * (1.-beta(idlak))
         else if (j == 1) then
            phidum= sabg * (1.-beta(idlak))
         else
            phidum = 0.
         end if
         phi(j) = phidum
      end do

! sum cwat*tlak*dzlak for energy check

      ocvts = 0.
      do j = 1, msl
         ocvts = ocvts + cwat*tlak(j)*dzlak(j) 
      end do

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

      j = 1
      m2 = dzlak(j)/kme(j) + dzlak(j+1)/kme(j+1)
      m3 = dtime/dzlak(j)
      r(j) = tlak(j) + (fin+phi(j))*m3/cwat - (tlak(j)-tlak(j+1))*m3/m2
      a(j) = 0.
      b(j) = 1. + m3/m2
      c(j) = -m3/m2

      j = msl
      m1 = dzlak(j-1)/kme(j-1) + dzlak(j)/kme(j)
      m3 = dtime/dzlak(j)
      r(j) = tlak(j) + phi(j)*m3/cwat + (tlak(j-1)-tlak(j))*m3/m1
      a(j) = -m3/m1
      b(j) = 1. + m3/m1
      c(j) = 0.

      do j = 2, msl-1
         m1 = dzlak(j-1)/kme(j-1) + dzlak(j  )/kme(j  )
         m2 = dzlak(j  )/kme(j  ) + dzlak(j+1)/kme(j+1)
         m3 = dtime/dzlak(j)
         r(j) = tlak(j) + phi(j)*m3/cwat &
              +(tlak(j-1)-tlak(j))*m3/m1 - (tlak(j)-tlak(j+1))*m3/m2

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

! solve for tlak: a, b, c, r, u go from 1 to nsoi. tlak = 1 to npt

      call tridia (msl ,a ,b ,c ,r ,tlak) 

! convective mixing: make sure cwat*dzlak*ts is conserved. mixing
! is only allowed for unfrozen deep lakes. mix every 3 time steps

      if (mod(istep,3) == 0) then
         if(idlak == 1 .AND. tg > tfrz) then
            do j = 1, msl-1
               if(rhow(j) > rhow(j+1)) then

                  tav = 0.
                  nav = 0.
                  do i = 1, j+1
                     tav = tav + tlak(i)*dzlak(i)
                     nav = nav + dzlak(i)
                  end do
                  tav = tav/nav
   
                  do i = 1, j+1
                     tlak(i) = tav
                     rhow(i) = 1000.*( 1.0 &
                             - 1.9549e-05*(abs(tlak(i)-277.))**1.68 )
                  end do

               end if
            end do
         end if
      end if

! sum cwat*tlak*dzlak and total energy into lake for energy check

      ncvts = 0.
      do j = 1, msl
         ncvts = ncvts + cwat*tlak(j)*dzlak(j) 
         fin = fin + phi(j)
      end do

      errore = (ncvts-ocvts) / dtime - fin

! ----------------------------------------------------------------------
! [4] snow on the lake ice 
! ----------------------------------------------------------------------

      qseva = 0.
      qsubl = 0.
      qfros = 0.
      qsdew = 0.

      if(fevpg >= 0.)then
! sublimation. do not allow for more sublimation than there is snow
! after melt. remaining surface evaporation used for infiltration
         qsubl = min( fevpg, scv/dtime-qmelt )
         qseva = fevpg - qsubl
      else
         if(tg < tfrz-0.1)then
            qfros = abs(fevpg)
         else
            qsdew = abs(fevpg)
         endif
      endif

! update snow pack
      scv = scv + (snowrate-qmelt-qsubl+qfros)*dtime
      scv = max( scv, 0. )

! no snow if lake unfrozen
      if (tg > tfrz) scv = 0.

! snow height and fractional coverage
      snowdp = scv/250.       !assumed a constant snow bulk density = 250.

! water mass
      do j = 1, msl
         wice(j) = 0.
         wliq(j) = 1000.*dzlak(j)
      enddo

      END SUBROUTINE lake