meltf.f90

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

  SUBROUTINE meltf ( lb, iub, dtime, &
                     fact, brr, hs, dhsdT, &
                     tssbef, tss, wliq, wice, imelt, &
                     scv, snowdp, sm, xmf )

!=======================================================================
!      Source file: meltf.f90
! Original version: Yongjiu Dai, September 15, 1999
!
! calculation of the phase change within snow and soil layers:
! 
! (1) check the conditions which the phase change may take place,
!     i.e., the layer temperature is great than the freezing point
!     and the ice mass is not equal to zero (i.e., melting),
!     or layer temperature is less than the freezing point
!     and the liquid water mass is not equal to zero (i.e., freezing);
! (2) assess the rate of phase change from the energy excess (or deficit)
!     after setting the layer temperature to freezing point;
! (3) re-adjust the ice and liquid mass, and the layer temperature
!
!=======================================================================

  USE PHYCON_MODULE ! physical constant
  IMPLICIT NONE

! Dummy argument

  integer, INTENT(in) :: & 
        iub,             &! upper bound of array (i.e., soil layers)
        lb                ! lower bound of array (i.e., snl +1)

  real, INTENT(in) :: &
        dtime             ! time step [second]

  real, INTENT(in) :: &
        tssbef(lb : iub),&! temperature at previous time step [K]
        brr   (lb : iub),&! 
        fact  (lb : iub),&! temporary variables
        hs,              &! net ground heat flux into the surface
        dhsdT             ! temperature derivative of "hs"

  real, INTENT(inout) :: &
        tss   (lb : iub),&! temperature at current time step [K]
        wice  (lb : iub),&! ice lens [kg/m2]
        wliq  (lb : iub),&! liquid water [kg/m2]
        scv,             &! snow mass [kg/m2]
        snowdp            ! snow depth [m]

  real, INTENT(out) :: &
        sm,              &! rate of snowmelt [mm/s, kg/(m2 s)]
        xmf               ! total latent heat of phase change

  integer, INTENT(out) :: & 
        imelt(lb : iub)   ! flag for melting or freezing [-]

!-----------------------------------------------------------------------
! Local 
  integer j
  real  hm(lb : iub),    &! energy residual [W/m2]
        xm(lb : iub),    &! metling or freezing within a time step [kg/m2]
        heatr,           &! energy residual or loss after melting or freezing
        temp1,           &! temporary variables [kg/m2]
        temp2             ! temporary variables [kg/m2]

  real, dimension(lb : iub) :: wmass0, wice0, wliq0
  real  propor,tinc  
!
!=======================================================================
!
      sm = 0.
      xmf = 0.
      do j = lb, iub
         imelt(j) = 0
         hm(j) = 0.
         xm(j) = 0.
         wice0(j) = wice(j)
         wliq0(j) = wliq(j)
         wmass0(j) = wice(j) + wliq(j)
      enddo

! Melting identification
! if ice exists above melt point, melt some to liquid.
!-----------------------------------------------------------------------
      do j = lb, iub
         if(wice(j) > 0. .AND. tss(j) > tfrz)then
            imelt(j) = 1
            tss(j) = tfrz
         endif

! Freezing identification
! if liquid exists below melt point, freeze some to ice.
!-----------------------------------------------------------------------
         if(wliq(j) > 0. .AND. tss(j) < tfrz) then
            imelt(j) = 2
            tss(j) = tfrz
         endif
      enddo

! If snow exists, but its thickness less than the critical value (0.01 m)
!-----------------------------------------------------------------------
      if(lb == 1 .AND. scv > 0.)then
         if(tss(1) > tfrz)then
            imelt(1) = 1
            tss(1) = tfrz
         endif
      endif

!-----------------------------------------------------------------------
! Calculate the energy surplus and loss for melting and freezing
      do j = lb, iub
         if(imelt(j) > 0)then
            tinc = tss(j)-tssbef(j)
            if(j > lb)then
               hm(j) = brr(j) - tinc/fact(j) 
            else
               hm(j) = hs + dhsdT*tinc + brr(j) - tinc/fact(j) 
            endif
         endif
      enddo

      do j = lb, iub
         if(imelt(j) == 1 .AND. hm(j) < 0.) then
           hm(j) = 0.
           imelt(j) = 0
         endif

! this error was checked carefully, it results from the the computed error
! of "Tridiagonal-Matrix" in subroutine "thermal".

         if(imelt(j) == 2 .AND. hm(j) > 0.) then
           hm(j) = 0.
           imelt(j) = 0
         endif
      enddo

!-----------------------------------------------------------------------
! The rate of melting and freezing
      do j = lb, iub
      if(imelt(j) > 0 .AND. abs(hm(j)) > .0) then
         xm(j) = hm(j)*dtime/dlm                        ! kg/m2

! if snow exists, but its thickness less than the critical value (1 cm)
! Note: more work is need on how to tune the snow depth at this case
!-----------------------------------------------------------------------
         if(j == 1)then
         if ((lb == 1) .AND. (scv > 0.) .AND. (xm(j) > 0.))then
            temp1 = scv                                 ! kg/m2
            scv = max(0.,temp1-xm(j))
            propor = scv/temp1
            snowdp = propor * snowdp
            heatr = hm(j) - dlm*(temp1-scv)/dtime       ! W/m2
            if(heatr > 0.) then
               xm(j) = heatr*dtime/dlm                  ! kg/m2
               hm(j) = heatr                            ! W/m2
            else
               xm(j) = 0.
               hm(j) = 0.
            endif
            sm = max(0.,(temp1-scv))/dtime              ! kg/(m2 s)
            xmf = dlm*sm
         endif
         endif
!
!-----------------------------------------------------------------------
!
         heatr = 0.
         if(xm(j) > 0.) then
            wice(j) = max(0., wice0(j)-xm(j))
            heatr = hm(j) - dlm*(wice0(j)-wice(j))/dtime
         else if(xm(j) < 0.) then
            wice(j) = min(wmass0(j), wice0(j)-xm(j))
            heatr = hm(j) - dlm*(wice0(j)-wice(j))/dtime  
         endif

         wliq(j) = max(0.,wmass0(j)-wice(j))

         if(abs(heatr) > 0.)then
            if(j > lb)then
               tss(j) = tss(j) + fact(j)*heatr
            else
               tss(j) = tss(j) + fact(j)*heatr/(1.-fact(j)*dhsdT)
            endif
            if(wliq(j)*wice(j) > 0.) tss(j) = tfrz
         endif

         xmf = xmf + dlm * (wice0(j)-wice(j))/dtime

         if(imelt(j) == 1 .AND. j < 1) &
         sm = sm + max(0.,(wice0(j)-wice(j)))/dtime  

      endif
      enddo

  END SUBROUTINE meltf