obult.f90

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

  SUBROUTINE obult ( hu,ht,hq,displa,z0m,z0h,z0q,obu,um,ustar,temp1,temp2,temp12m )

! ======================================================================
!      Source file: obult.f90
! Original version: Yongjiu Dai, September 15, 1999
!
! calculation of friction velocity, relation for potential temperatur
! and humidity profiles of surface boundary layer.
! the scheme is based on the work of Zeng et al. (1998):
! Intercomparison of bulk aerodynamic algorithms for the computation
! of sea surface fluxes using TOGA CORE and TAO data. J. Climate, Vol. 11: 2628-2644
! ======================================================================

  USE PHYCON_MODULE    ! physical constant
  IMPLICIT NONE
!
! ---------------------- dummy argument --------------------------------
!
  real, INTENT(in) :: &!
        hu,           &! observational height of wind [m]
        ht,           &! observational height of temperature [m]
        hq,           &! observational height of humidity [m]
        displa,       &! displacement height [m]
        z0m,          &! roughness length, momentum [m]
        z0h,          &! roughness length, sensible heat [m]
        z0q,          &! roughness length, latent heat [m]
        obu,          &! monin-obukhov length (m)
        um             ! wind speed including the stablity effect [m/s]

  real, INTENT(out) ::&!
        ustar,        &! friction velocity [m/s]
        temp1,        &! relation for potential temperature profile
        temp12m,      &!
        temp2          ! relation for specific humidity profile
!
!------------------------ local variables ------------------------------
!
  real  zldis,        &! reference height "minus" zero displacement heght [m]
        psi,          &! stability function for unstable case
        zetam,        &! transition point of flux-gradient relation (wind profile)
        zetat,        &! transition point of flux-gradient relation (temp. profile)
        zeta           ! dimensionless height used in Monin-Obukhov theory
!
!-----------------------------------------------------------------------
! adjustment factors for unstable (moz < 0) or stable (moz > 0) conditions.

! wind profile
        zldis=hu-displa
        zeta=zldis/obu
        zetam=1.574
        if(zeta < -zetam)then           ! zeta < -1
          ustar = vonkar*um / ( alog(-zetam*obu/z0m) - psi(1,-zetam) &
                + psi(1,z0m/obu) + 1.14*((-zeta)**0.333-(zetam)**0.333) )
        else if(zeta < 0.)then          ! -1 <= zeta < 0
          ustar = vonkar*um / ( alog(zldis/z0m) - psi(1,zeta) + psi(1,z0m/obu) )
        else if(zeta <= 1.)then         !  0 <= ztea <= 1
          ustar = vonkar*um / ( alog(zldis/z0m) + 5.*zeta - 5.*z0m/obu )
        else                            !  1 < zeta, phi=5+zeta
          ustar = vonkar*um / ( alog(obu/z0m) + 5. - 5.*z0m/obu + (5.*alog(zeta)+zeta-1.) )
        endif

! temperature profile
        zldis=ht-displa
        zeta=zldis/obu
        zetat=0.465
        if(zeta < -zetat)then           ! zeta < -1
          temp1 = vonkar / ( alog(-zetat*obu/z0h)-psi(2,-zetat) &
                + psi(2,z0h/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) )
        else if(zeta < 0.)then          ! -1 <= zeta < 0
          temp1 = vonkar / ( alog(zldis/z0h) - psi(2,zeta) + psi(2,z0h/obu) )
        else if(zeta <= 1.)then         !  0 <= ztea <= 1
          temp1 = vonkar / ( alog(zldis/z0h) + 5.*zeta - 5.*z0h/obu )
        else                            !  1 < zeta, phi=5+zeta
          temp1 = vonkar / ( alog(obu/z0h) + 5. - 5.*z0h/obu + (5.*alog(zeta)+zeta-1.) )
        endif

        ! for 2 meter screen temperature
        zldis=2.+z0h  ! ht-displa
        zeta=zldis/obu
        zetat=0.465
        if(zeta < -zetat)then           ! zeta < -1
          temp12m = vonkar / ( alog(-zetat*obu/z0h)-psi(2,-zetat) &
                + psi(2,z0h/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) )
        else if(zeta < 0.)then          ! -1 <= zeta < 0
          temp12m = vonkar / ( alog(zldis/z0h) - psi(2,zeta) + psi(2,z0h/obu) )
        else if(zeta <= 1.)then         !  0 <= ztea <= 1
          temp12m = vonkar / ( alog(zldis/z0h) + 5.*zeta - 5.*z0h/obu )
        else                            !  1 < zeta, phi=5+zeta
          temp12m = vonkar / ( alog(obu/z0h) + 5. - 5.*z0h/obu + (5.*alog(zeta)+zeta-1.) )
        endif

! humidity profile
        zldis=hq-displa
        zeta=zldis/obu
        zetat=0.465
        if(zeta < -zetat)then           ! zeta < -1
          temp2 = vonkar / ( alog(-zetat*obu/z0q) - psi(2,-zetat) &
                + psi(2,z0q/obu) + 0.8*((zetat)**(-0.333)-(-zeta)**(-0.333)) )
        else if(zeta < 0.)then          ! -1 <= zeta < 0
          temp2 = vonkar / ( alog(zldis/z0q) - psi(2,zeta) + psi(2,z0q/obu) )
        else if(zeta <= 1.)then         !  0 <= ztea <= 1
          temp2 = vonkar / ( alog(zldis/z0q) + 5.*zeta - 5.*z0q/obu )
        else                            !  1 < zeta, phi=5+zeta
          temp2 = vonkar / ( alog(obu/z0q) + 5. - 5.*z0q/obu + (5.*alog(zeta)+zeta-1.) )
        endif

   END SUBROUTINE obult




   SUBROUTINE obuini(ur,th,thm,thv,dth,dqh,dthv,zldis,z0m,um,obu)

! ======================================================================
!      Source file: obuini.f90
! Original version: Yongjiu Dai, September 15, 1999
!
! initialzation of Monin-Obukhov length,
! the scheme is based on the work of Zeng et al. (1998):
! Intercomparison of bulk aerodynamic algorithms for the computation
! of sea surface fluxes using TOGA CORE and TAO data. J. Climate, Vol. 11: 2628-2644
! ======================================================================

   USE PHYCON_MODULE ! physical constant
   IMPLICIT NONE

! Dummy argument
  real, INTENT(in) :: &!
        ur,           &! wind speed at reference height [m/s]
        thm,          &! intermediate variable (tm+0.0098*ht)
        th,           &! potential temperature [kelvin]
        thv,          &! virtual potential temperature (kelvin)
        dth,          &! diff of virtual temp. between ref. height and surface
        dthv,         &! diff of vir. poten. temp. between ref. height and surface
        dqh,          &! diff of humidity between ref. height and surface
        zldis,        &! reference height "minus" zero displacement heght [m]
        z0m            ! roughness length, momentum [m]

  real, INTENT(out) ::&!
        um,           &! wind speed including the stablity effect [m/s]
        obu            ! monin-obukhov length (m)

! Local       
  real  wc,           &! convective velocity [m/s]
        rib,          &! bulk Richardson number
        zeta,         &! dimensionless height used in Monin-Obukhov theory
        ustar          ! friction velocity [m/s]

!-----------------------------------------------------------------------
! Initial values of u* and convective velocity
!
        ustar=0.06
        wc=0.5
        if(dthv >= 0.)then
          um=max(ur,0.1)
        else
          um=sqrt(ur*ur+wc*wc)
        endif

        rib=grv*zldis*dthv/(thv*um*um)

        if(rib >= 0.)then      ! neutral or stable
          zeta = rib*alog(zldis/z0m)/(1.-5.*min(rib,0.19))
          zeta = min(2.,max(zeta,1.e-6))
        else                   ! unstable
          zeta = rib*alog(zldis/z0m)
          zeta = max(-100.,min(zeta,-1.e-6))
        endif
        obu=zldis/zeta

   END SUBROUTINE obuini



   FUNCTION psi(k,zeta)

!=======================================================================
! stability function for rib < 0

      IMPLICIT NONE
      
      integer k    !
      real zeta,  &! dimensionless height used in Monin-Obukhov theory
           psi,   &! stability function for unstable case
           chik    ! 

      chik = (1.-16.*zeta)**0.25
      if(k == 1)then
        psi = 2.*alog((1.+chik)*0.5)+alog((1.+chik*chik)*0.5)-2.*atan(chik)+2.*atan(1.)
      else
        psi = 2.*alog((1.+chik*chik)*0.5)
      endif

   END FUNCTION psi