scale_cpl_phy_sfc_fixed_temp.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
#include "scalelib.h"
module scale_cpl_phy_sfc_fixed_temp
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_io
  use scale_prof
  use scale_cpl_sfc_index
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: cpl_phy_sfc_fixed_temp_setup
  public :: cpl_phy_sfc_fixed_temp

  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  logical :: initialized = .false.

  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  subroutine cpl_phy_sfc_fixed_temp_setup
    implicit none
    !---------------------------------------------------------------------------

    if ( initialized ) return

    log_newline
    log_info("CPL_PHY_SFC_FIXED_TEMP_setup",*) 'Setup'

    initialized = .true.

    return
  end subroutine cpl_phy_sfc_fixed_temp_setup

  !-----------------------------------------------------------------------------
  subroutine cpl_phy_sfc_fixed_temp( &
       IA, IS, IE,          &
       JA, JS, JE,          &
       TMPA, PRSA,          &
       WA, UA, VA,          &
       RHOA, QVA, LH,       &
       Z1, PBL,             &
       RHOS, PRSS,          &
       RFLXD,               &
       TMPS, QVEF,          &
       ALBEDO,              &
       Rb, Z0M, Z0H, Z0E,   &
       calc_flag, dt,       &
       ZMFLX, XMFLX, YMFLX, &
       SHFLX, QVFLX, GFLX,  &
       U10, V10, T2, Q2     )
    use scale_const, only: &
       pre00 => const_pre00, &
       rdry  => const_rdry,  &
       cpdry => const_cpdry, &
       rvap  => const_rvap,  &
       stb   => const_stb
    use scale_atmos_saturation, only: &
       qsat => atmos_saturation_pres2qsat_all
    use scale_bulkflux, only: &
       bulkflux
    implicit none

    integer,  intent(in)  :: ia, is, ie
    integer,  intent(in)  :: JA, JS, JE
    real(RP), intent(in)  :: TMPA     (ia,ja)                     ! temperature at the lowest atmospheric layer [K]
    real(RP), intent(in)  :: PRSA     (ia,ja)                     ! pressure    at the lowest atmospheric layer [Pa]
    real(RP), intent(in)  :: WA       (ia,ja)                     ! velocity w  at the lowest atmospheric layer [m/s]
    real(RP), intent(in)  :: UA       (ia,ja)                     ! velocity u  at the lowest atmospheric layer [m/s]
    real(RP), intent(in)  :: VA       (ia,ja)                     ! velocity v  at the lowest atmospheric layer [m/s]
    real(RP), intent(in)  :: RHOA     (ia,ja)                     ! density     at the lowest atmospheric layer [kg/m3]
    real(RP), intent(in)  :: QVA      (ia,ja)                     ! ratio of water vapor mass to total mass at the lowest atmospheric layer [kg/kg]
    real(RP), intent(in)  :: LH       (ia,ja)                     ! latent heat at the lowest atmospheric layer [J/kg]
    real(RP), intent(in)  :: Z1       (ia,ja)                     ! cell center height at the lowest atmospheric layer [m]
    real(RP), intent(in)  :: PBL      (ia,ja)                     ! the top of atmospheric mixing layer [m]
    real(RP), intent(in)  :: RHOS     (ia,ja)                     ! density  at the surface [kg/m3]
    real(RP), intent(in)  :: PRSS     (ia,ja)                     ! pressure at the surface [Pa]
    real(RP), intent(in)  :: RFLXD    (ia,ja,n_rad_dir,n_rad_rgn) ! downward radiation flux at the surface (direct/diffuse,IR/near-IR/VIS) [J/m2/s]
    real(RP), intent(in)  :: TMPS     (ia,ja)                     ! surface temperature [K]
    real(RP), intent(in)  :: QVEF     (ia,ja)                     ! efficiency of evaporation (0-1)
    real(RP), intent(in)  :: ALBEDO   (ia,ja,n_rad_dir,n_rad_rgn) ! surface albedo (direct/diffuse,IR/near-IR/VIS) (0-1)
    real(RP), intent(in)  :: Rb       (ia,ja)                     ! stomata resistance [1/s]
    real(RP), intent(in)  :: Z0M      (ia,ja)                     ! roughness length for momemtum [m]
    real(RP), intent(in)  :: Z0H      (ia,ja)                     ! roughness length for heat     [m]
    real(RP), intent(in)  :: Z0E      (ia,ja)                     ! roughness length for vapor    [m]
    logical,  intent(in)  :: calc_flag(ia,ja)                     ! to decide calculate or not
    real(DP), intent(in)  :: dt                                   ! delta time
    real(RP), intent(out) :: ZMFLX    (ia,ja)                     ! z-momentum      flux at the surface [kg/m/s2]
    real(RP), intent(out) :: XMFLX    (ia,ja)                     ! x-momentum      flux at the surface [kg/m/s2]
    real(RP), intent(out) :: YMFLX    (ia,ja)                     ! y-momentum      flux at the surface [kg/m/s2]
    real(RP), intent(out) :: SHFLX    (ia,ja)                     ! sensible heat   flux at the surface [J/m2/s]
    real(RP), intent(out) :: QVFLX    (ia,ja)                     ! water vapor     flux at the surface [kg/m2/s]
    real(RP), intent(out) :: GFLX     (ia,ja)                     ! subsurface heat flux at the surface [J/m2/s]
    real(RP), intent(out) :: U10      (ia,ja)                     ! velocity u  at 10m [m/s]
    real(RP), intent(out) :: V10      (ia,ja)                     ! velocity v  at 10m [m/s]
    real(RP), intent(out) :: T2       (ia,ja)                     ! temperature at 2m  [K]
    real(RP), intent(out) :: Q2       (ia,ja)                     ! water vapor at 2m  [kg/kg]

    real(RP) :: emis    ! surface longwave emission                 [J/m2/s]
    real(RP) :: LWD     ! surface downward longwave  radiation flux [J/m2/s]
    real(RP) :: LWU     ! surface upward   longwave  radiation flux [J/m2/s]
    real(RP) :: SWD     ! surface downward shortwave radiation flux [J/m2/s]
    real(RP) :: SWU     ! surface upward   shortwave radiation flux [J/m2/s]
    real(RP) :: res     ! residual

    real(RP) :: Ustar   ! friction velocity               [m]
    real(RP) :: Tstar   ! friction potential temperature  [K]
    real(RP) :: Qstar   ! friction water vapor mass ratio [kg/kg]
    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     ! 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 ! calculation parameter for U10 [1]
    real(RP) :: FracT2  ! calculation parameter for T2  [1]
    real(RP) :: FracQ2  ! calculation parameter for Q2  [1]

    integer  :: i, j
    !---------------------------------------------------------------------------

    log_progress(*) 'coupler / physics / surface / FIXED-TEMP'

    ! calculate surface flux
#ifndef __GFORTRAN__
    !$omp parallel do default(none) &
    !$omp private(qdry,Rtot,QVsat,QVS,Ustar,Tstar,Qstar,Uabs,Ra,FracU10,FracT2,FracQ2,res,emis,LWD,LWU,SWD,SWU) &
    !$omp shared(IS,IE,JS,JE,Rdry,CPdry,bulkflux, &
    !$omp        calc_flag,TMPA,QVA,LH,UA,VA,WA,Z1,PBL,PRSA,TMPS,PRSS,RHOS,QVEF,Z0M,Z0H,Z0E,ALBEDO,RFLXD,Rb, &
    !$omp        SHFLX,QVFLX,GFLX,ZMFLX,XMFLX,YMFLX,U10,V10,T2,Q2)
#else
    !$omp parallel do default(shared) &
    !$omp private(qdry,Rtot,QVsat,QVS,Ustar,Tstar,Qstar,Uabs,Ra,FracU10,FracT2,FracQ2,res,emis,LWD,LWU,SWD,SWU)
#endif
    do j = js, je
    do i = is, ie
       if ( calc_flag(i,j) ) then

          qdry = 1.0_rp - qva(i,j)
          rtot = qdry * rdry + qva(i,j) * rvap

          call qsat( tmps(i,j), prss(i,j), qdry, qvsat )

          qvs = ( 1.0_rp-qvef(i,j) ) * qva(i,j) &
              + (        qvef(i,j) ) * qvsat

          call bulkflux( ustar,     & ! [OUT]
                         tstar,     & ! [OUT]
                         qstar,     & ! [OUT]
                         uabs,      & ! [OUT]
                         ra,        & ! [OUT]
                         fracu10,   & ! [OUT]
                         fract2,    & ! [OUT]
                         fracq2,    & ! [OUT]
                         tmpa(i,j), & ! [IN]
                         tmps(i,j), & ! [IN]
                         prsa(i,j), & ! [IN]
                         prss(i,j), & ! [IN]
                         qva(i,j), & ! [IN]
                         qvs,       & ! [IN]
                         ua(i,j), & ! [IN]
                         va(i,j), & ! [IN]
                         z1(i,j), & ! [IN]
                         pbl(i,j), & ! [IN]
                         z0m(i,j), & ! [IN]
                         z0h(i,j), & ! [IN]
                         z0e(i,j)  ) ! [IN]

          zmflx(i,j) = -rhos(i,j) * ustar * ustar / uabs * wa(i,j)
          xmflx(i,j) = -rhos(i,j) * ustar * ustar / uabs * ua(i,j)
          ymflx(i,j) = -rhos(i,j) * ustar * ustar / uabs * va(i,j)
          shflx(i,j) = -rhos(i,j) * ustar * tstar * cpdry
          qvflx(i,j) = -rhos(i,j) * ustar * qstar * ra / ( ra+rb(i,j) )

          emis = ( 1.0_rp-albedo(i,j,i_r_diffuse,i_r_ir) ) * stb * tmps(i,j)**4

          lwd  = rflxd(i,j,i_r_diffuse,i_r_ir)
          lwu  = rflxd(i,j,i_r_diffuse,i_r_ir)  * albedo(i,j,i_r_diffuse,i_r_ir) + emis
          swd  = rflxd(i,j,i_r_direct ,i_r_nir) &
               + rflxd(i,j,i_r_diffuse,i_r_nir) &
               + rflxd(i,j,i_r_direct ,i_r_vis) &
               + rflxd(i,j,i_r_diffuse,i_r_vis)
          swu  = rflxd(i,j,i_r_direct ,i_r_nir) * albedo(i,j,i_r_direct ,i_r_nir) &
               + rflxd(i,j,i_r_diffuse,i_r_nir) * albedo(i,j,i_r_diffuse,i_r_nir) &
               + rflxd(i,j,i_r_direct ,i_r_vis) * albedo(i,j,i_r_direct ,i_r_vis) &
               + rflxd(i,j,i_r_diffuse,i_r_vis) * albedo(i,j,i_r_diffuse,i_r_vis)

          ! calculation for residual
          res = swd - swu + lwd - lwu - shflx(i,j) - qvflx(i,j) * lh(i,j)

          ! put residual in ground heat flux (positive for upward)
          gflx(i,j) = -res

          ! diagnostic variables considering unstable/stable state
          !U10(i,j) = FracU10 * UA(i,j)
          !V10(i,j) = FracU10 * VA(i,j)
          !T2 (i,j) = ( 1.0_RP - FracT2 ) * TMPS(i,j) + FracT2 * TMPA(i,j)
          !Q2 (i,j) = ( 1.0_RP - FracQ2 ) * QVS       + FracQ2 * QVA (i,j)

          ! diagnostic variables for neutral state
          u10(i,j) = ua(i,j) * log( 10.0_rp / z0m(i,j) ) / log( z1(i,j) / z0m(i,j) )
          v10(i,j) = va(i,j) * log( 10.0_rp / z0m(i,j) ) / log( z1(i,j) / z0m(i,j) )
          t2(i,j) = tmps(i,j) + ( tmpa(i,j) - tmps(i,j) ) * log( 2.0_rp / z0h(i,j) ) / log( z1(i,j) / z0h(i,j) )
          q2(i,j) = qvs       + (  qva(i,j) - qvs       ) * log( 2.0_rp / z0e(i,j) ) / log( z1(i,j) / z0e(i,j) )

       else ! not calculate surface flux
          zmflx(i,j) = 0.0_rp
          xmflx(i,j) = 0.0_rp
          ymflx(i,j) = 0.0_rp
          shflx(i,j) = 0.0_rp
          qvflx(i,j) = 0.0_rp
          gflx(i,j) = 0.0_rp
          u10(i,j) = 0.0_rp
          v10(i,j) = 0.0_rp
          t2(i,j) = 0.0_rp
          q2(i,j) = 0.0_rp
       endif
    enddo
    enddo

    return
  end subroutine cpl_phy_sfc_fixed_temp

end module scale_cpl_phy_sfc_fixed_temp