scale_land_phy_snow_diagnos.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
#include "scalelib.h"
module scale_land_phy_snow_diagnos
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_io
  use scale_prof
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: land_phy_snow_diags
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  subroutine land_phy_snow_diags( &
       LIA, LIS, LIE, LJA, LJS, LJE, &
       SNOW_frac,           &
       TMPA, PRSA,          &
       WA, UA, VA,          &
       RHOA, QVA,           &
       Z1, PBL,             &
       RHOS, PRSS, LST1,    &
       QVEF,                &
       Z0M, Z0H, Z0E,       &
       ZMFLX, XMFLX, YMFLX, &
       Ustar, Tstar, Qstar, &
       Wstar,               &
       RLmo,                &
       U10, V10,            &
       T2, Q2               )
    use scale_const, only: &
      eps   => const_eps,   &
      undef => const_undef, &
      pre00 => const_pre00, &
      rdry  => const_rdry,  &
      cpdry => const_cpdry, &
      rvap  => const_rvap,  &
      stb   => const_stb
    use scale_atmos_saturation, only: &
      qsat => atmos_saturation_dens2qsat_all
!      qsat => ATMOS_SATURATION_pres2qsat_all
    use scale_bulkflux, only: &
      bulkflux, &
      bulkflux_diagnose_surface
    implicit none
    integer, intent(in) :: lia, lis, lie
    integer, intent(in) :: lja, ljs, lje
 
    real(rp), intent(in) :: snow_frac(lia,lja) ! snow fraction
    real(rp), intent(in) :: tmpa(lia,lja) ! temperature at the lowest atmospheric layer [K]
    real(rp), intent(in) :: prsa(lia,lja) ! pressure at the lowest atmospheric layer [Pa]
    real(rp), intent(in) :: wa  (lia,lja) ! velocity w at the lowest atmospheric layer [m/s]
    real(rp), intent(in) :: ua  (lia,lja) ! velocity u at the lowest atmospheric layer [m/s]
    real(rp), intent(in) :: va  (lia,lja) ! velocity v at the lowest atmospheric layer [m/s]
    real(rp), intent(in) :: rhoa(lia,lja) ! density at the lowest atmospheric layer [kg/m3]
    real(rp), intent(in) :: qva (lia,lja) ! ratio of water vapor mass to total mass at the lowest atmospheric layer [kg/kg]
    real(rp), intent(in) :: z1  (lia,lja) ! cell center height at the lowest atmospheric layer [m]
    real(rp), intent(in) :: pbl (lia,lja) ! the top of atmospheric mixing layer [m]
    real(rp), intent(in) :: rhos(lia,lja) ! density  at the surface [kg/m3]
    real(rp), intent(in) :: prss(lia,lja) ! pressure at the surface [Pa]
 
    real(rp), intent(in) :: lst1  (lia,lja) ! land surface temperature [K]
    real(rp), intent(in) :: qvef  (lia,lja) ! efficiency of evaporation (0-1)
 
    real(rp), intent(in) :: z0m   (lia,lja) ! roughness length for momemtum [m]
    real(rp), intent(in) :: z0h   (lia,lja) ! roughness length for heat [m]
    real(rp), intent(in) :: z0e   (lia,lja) ! roughness length for vapor [m]
 
    real(rp), intent(out) :: zmflx(lia,lja) ! z-momentum flux at the surface [kg/m/s2]
    real(rp), intent(out) :: xmflx(lia,lja) ! x-momentum flux at the surface [kg/m/s2]
    real(rp), intent(out) :: ymflx(lia,lja) ! y-momentum flux at the surface [kg/m/s2]
    real(rp), intent(out) :: ustar(lia,lja) ! friction velocity         [m/s]
    real(rp), intent(out) :: tstar(lia,lja) ! temperature scale         [K]
    real(rp), intent(out) :: qstar(lia,lja) ! moisture scale            [kg/kg]
    real(rp), intent(out) :: wstar(lia,lja) ! convective velocity scale [m/s]
    real(rp), intent(out) :: rlmo (lia,lja) ! inversed Obukhov length   [1/m]
    real(rp), intent(out) :: u10  (lia,lja) ! velocity u at 10m [m/s]
    real(rp), intent(out) :: v10  (lia,lja) ! velocity v at 10m [m/s]
    real(rp), intent(out) :: t2   (lia,lja) ! temperature at 2m [K]
    real(rp), intent(out) :: q2   (lia,lja) ! water vapor at 2m [kg/kg]
 
 
    ! works
    real(rp) :: uabs ! modified absolute velocity [m/s]
    real(rp) :: ra   ! Aerodynamic resistance (=1/Ce) [1/s]
 
    real(rp) :: qvsat        ! saturation water vapor mixing ratio at surface [kg/kg]
    real(rp) :: qvs(lia,lja) ! water vapor mixing ratio at surface [kg/kg]
!    real(RP) :: Rtot         ! total gas constant
!    real(RP) :: qdry         ! dry air mass ratio [kg/kg]
 
    real(rp) :: fracu10(lia,lja) ! calculation parameter for U10 [-]
    real(rp) :: fract2 (lia,lja) ! calculation parameter for T2 [-]
    real(rp) :: fracq2 (lia,lja) ! calculation parameter for Q2 [-]
 
    real(rp) :: mflux
 
    logical :: calc_flag(lia,lja)
 
    integer  :: i, j
    !---------------------------------------------------------------------------
 
    log_info("LAND_PHY_SNOW_DIAGS",*) 'Snow surface diagnostic'
 
    ! calculate surface flux
    !$omp parallel do &
    !$omp private(QVsat,Uabs,MFLUX)
    do j = ljs, lje
    do i = lis, lie
 
      if( snow_frac(i,j) > 0.0_rp ) then
 
         calc_flag(i,j) = .true.
 
!        qdry = 1.0_RP - QVA(i,j)
!        Rtot = qdry * Rdry + QVA(i,j) * Rvap
 
!        call qsat( LST1(i,j),  PRSS(i,j), qdry, & ! [IN]
!                   QVsat                        ) ! [OUT]
        call qsat( lst1(i,j),  rhos(i,j), & ! [IN]
                   qvsat                  ) ! [OUT]
 
        qvs(i,j)  = ( 1.0_rp - qvef(i,j) ) * qva(i,j) + qvef(i,j) * qvsat
 
        uabs = sqrt( wa(i,j)**2 + ua(i,j)**2 + va(i,j)**2 )
 
        call bulkflux( tmpa(i,j), lst1(i,j),                  & ! [IN]
                       prsa(i,j), prss(i,j),                  & ! [IN]
                       qva(i,j), qvs(i,j),                  & ! [IN]
                       uabs, z1(i,j), pbl(i,j),               & ! [IN]
                       z0m(i,j), z0h(i,j), z0e(i,j),          & ! [IN]
                       ustar(i,j), tstar(i,j), qstar(i,j),    & ! [OUT]
                       wstar(i,j), rlmo(i,j), ra,             & ! [OUT]
                       fracu10(i,j), fract2(i,j), fracq2(i,j) ) ! [OUT]
 
        if ( uabs < eps ) then
           zmflx(i,j) = 0.0_rp
           xmflx(i,j) = 0.0_rp
           ymflx(i,j) = 0.0_rp
        else
           mflux = - rhos(i,j) * ustar(i,j)**2
           zmflx(i,j) = mflux * wa(i,j) / uabs
           xmflx(i,j) = mflux * ua(i,j) / uabs
           ymflx(i,j) = mflux * va(i,j) / uabs
        end if
 
      else
 
         ! not calculate surface flux
 
         calc_flag(i,j) = .false.
 
         zmflx(i,j) = undef
         xmflx(i,j) = undef
         ymflx(i,j) = undef
         ustar(i,j) = undef
         tstar(i,j) = undef
         qstar(i,j) = undef
         wstar(i,j) = undef
         rlmo(i,j) = undef
         u10(i,j) = undef
         v10(i,j) = undef
         t2(i,j) = undef
         q2(i,j) = undef
 
      end if
 
    end do
    end do
 
    call bulkflux_diagnose_surface( lia, lis, lie, lja, ljs, lje, &
                                    ua(:,:), va(:,:),                      & ! (in)
                                    tmpa(:,:), qva(:,:),                   & ! (in)
                                    lst1(:,:), qvs(:,:),                   & ! (in)
                                    z1(:,:), z0m(:,:), z0h(:,:), z0e(:,:), & ! (in)
                                    u10(:,:), v10(:,:), t2(:,:), q2(:,:),  & ! (out)
                                    mask = calc_flag(:,:),                 & ! (in)
                                    fracu10 = fracu10(:,:),                & ! (in)
                                    fract2 = fract2(:,:),                  & ! (in)
                                    fracq2 = fracq2(:,:)                   ) ! (in)
 
    return
  end subroutine land_phy_snow_diags
 
end module scale_land_phy_snow_diagnos