sfc_diag.f

 
 
      module sfc_diag
      contains
 
      subroutine sfc_diag_run (im,xlat_d,xlon_d,                        &
     &                    lsm,lsm_ruc,grav,cp,eps,epsm1,con_rocp,       &
     &                    con_karman,                                   &
     &                    shflx,cdq,wind,                               &
     &                    usfco,vsfco,icplocn2atm,                      &
     &                    zf,ps,u1,v1,t1,q1,prslki,evap,fm,fh,fm10,fh2, &
     &                    ust,tskin,qsurf,thsfc_loc,diag_flux,diag_log, &
     &                    use_lake_model,iopt_lake,iopt_lake_clm,       &
     &                    lake_t2m,lake_q2m,use_lake2m,                 &
     &                    f10m,u10m,v10m,t2m,q2m,dpt2m,errmsg,errflg    &
     &                   )
!
      use machine , only : kind_phys, kind_dbl_prec
      use funcphys, only : fpvs
      use physcons, only : con_t0c
      implicit none
!
      integer, intent(in) :: im, lsm, lsm_ruc, iopt_lake, iopt_lake_clm
      logical, intent(in) :: use_lake2m
      integer, intent(in) :: icplocn2atm
      logical, intent(in) :: thsfc_loc  ! Flag for reference pot. temp.
      logical, intent(in) :: diag_flux  ! Flag for flux method in 2-m diagnostics
      logical, intent(in) :: diag_log   ! Flag for 2-m log diagnostics under stable conditions
      real(kind=kind_phys), intent(in) :: grav,cp,eps,epsm1,con_rocp
      real(kind=kind_phys), intent(in) :: con_karman
      real(kind=kind_phys), dimension(:), intent( in) ::                &
     &                      zf, ps, u1, v1, t1, q1, ust, tskin,         &
     &                      usfco, vsfco,                               &
     &                      qsurf, prslki, evap, fm, fh, fm10, fh2,     &
     &                      shflx, cdq, wind, xlat_d, xlon_d
      real(kind=kind_phys), dimension(:), intent(out) ::                &
     &                      f10m, u10m, v10m, t2m, q2m, dpt2m
      real(kind=kind_phys), dimension(:), intent(in), optional ::       &
     &                      lake_t2m, lake_q2m
      integer, dimension(:), intent(in) :: use_lake_model
      character(len=*), intent(out) :: errmsg
      integer,          intent(out) :: errflg
!
!     locals
!
      real (kind_phys), parameter :: zero     = 0._kind_dbl_prec
      real (kind_phys), parameter :: one      = 1._kind_dbl_prec
      real (kind_phys), parameter :: qmin=1.0e-8
      integer :: k,i
 
      logical :: debug_print
      real(kind=kind_phys) :: q1c, qv, tem, qv1, th2m, x2m, rho
      real(kind=kind_phys) :: dt, dq, qsfcmr, qsfcprox, ff, fac, dz1
      real(kind=kind_phys) :: t2_alt, q2_alt
      real(kind=kind_phys) :: thcon, cqs, chs, chs2, cqs2
      real(kind=kind_phys) :: testptlat, testptlon
      logical :: have_2m
!
      real(kind=kind_phys) :: fhi, qss, wrk
!     real(kind=kind_phys) sig2k, fhi, qss
!
!     real, parameter :: g=grav
!
      ! Initialize CCPP error handling variables
      errmsg = ''
      errflg = 0
 
      !--
      testptlat = 35.3_kind_phys 
      testptlon = 273.0_kind_phys
      !--
      debug_print = .false.
!
!     estimate sigma ** k at 2 m
!
!     sig2k = 1. - 4. * g * 2. / (cp * 280.)
!
!  initialize variables. all units are supposedly m.k.s. unless specified
!  ps is in pascals
!
!!
 
      do i = 1, im
        f10m(i) = fm10(i) / fm(i)
        if (icplocn2atm ==0) then
          u10m(i) = f10m(i) * u1(i)
          v10m(i) = f10m(i) * v1(i)
        else if (icplocn2atm ==1) then
          u10m(i) = usfco(i)+f10m(i) * (u1(i)-usfco(i))
          v10m(i) = vsfco(i)+f10m(i) * (v1(i)-vsfco(i))
        endif
        have_2m = use_lake_model(i)>0 .and. use_lake2m .and.            &
     &                iopt_lake==iopt_lake_clm
        if(have_2m) then
           t2m(i) = lake_t2m(i)
           q2m(i) = lake_q2m(i)
        endif
        fhi     = fh2(i) / fh(i)
        wrk     = 1.0 - fhi
 
        if(lsm /= lsm_ruc) then
        !-- original method
          if(have_2m) then
             ! already have 2m T & Q from lake
          else
           if(thsfc_loc) then ! Use local potential temperature
            t2m(i)=tskin(i)*wrk + t1(i)*prslki(i)*fhi - (grav+grav)/cp
           else ! Use potential temperature referenced to 1000 hPa
            t2m(i) = tskin(i)*wrk + t1(i)*fhi - (grav+grav)/cp
           endif
           if(evap(i) >= 0.) then !  for evaporation>0, use inferred qsurf to deduce q2m
            q2m(i) = qsurf(i)*wrk + max(qmin,q1(i))*fhi
           else                   !  for dew formation, use saturated q at tskin
            qss    = fpvs(tskin(i))
            qss    = eps * qss/(ps(i) + epsm1 * qss)
            q2m(i) = qss*wrk + max(qmin,q1(i))*fhi
           endif
          endif
          qss    = fpvs(t2m(i))
          qss    = eps * qss / (ps(i) + epsm1 * qss)
          q2m(i) = min(q2m(i),qss)
        else
        !RUC lsm
          thcon    = (1.e5_kind_phys/ps(i))**con_rocp
          !-- make sure 1st level q is not higher than saturated value
          qss    = fpvs(t1(i))
          qss    = eps * qss / (ps(i) + epsm1 * qss)
          q1c    = min(q1(i),qss) ! lev 1 spec. humidity
 
          qv1    = q1c / (one - q1c) ! lev 1 mixing ratio
          qsfcmr   = qsurf(i)/(one - qsurf(i)) ! surface mixing ratio
          chs = cdq(i) * wind(i)
          cqs = chs
          chs2 = ust(i)*con_karman/fh2(i)
          cqs2 = chs2
          qsfcprox = max(qmin,qv1 + evap(i)/cqs) ! surface mix. ratio computed from the flux
 
         ruc_have_2m: if(.not.have_2m) then
           if(diag_flux) then
           !-- flux method
            th2m = tskin(i)*thcon - shflx(i)/chs2
            t2m(i) = th2m/thcon
            x2m = max(qmin,qsfcprox - evap(i)/cqs2) ! mix. ratio
            q2m(i) = x2m/(one + x2m) ! spec. humidity
           else
            t2m(i) = tskin(i)*wrk + t1(i)*fhi - (grav+grav)/cp
            q2m(i) = qsurf(i)*wrk + max(qmin,q1c)*fhi
           endif ! flux method
 
           if(diag_log) then
           !-- Alternative logarithmic diagnostics:
            dt = t1(i) - tskin(i)
            dq = qv1 - qsfcmr
            dz1= zf(i) ! level of atm. forcing
            IF (dt > zero) THEN
              ff  = min(max(one-dt/10._kind_phys,0.01_kind_phys), one)
              !for now, set zt = 0.05
              fac = log((2._kind_phys +.05_kind_phys)/(0.05_kind_phys + &
     &          ff))/log((dz1 + .05_kind_phys)/(0.05_kind_phys + ff))
              t2_alt = tskin(i) + fac * dt
            ELSE
            !no alternatives (yet) for unstable conditions
              t2_alt = t2m(i)
            ENDIF
 
            IF (dq > zero) THEN
              ff  = min(max(one-dq/0.003_kind_phys,0.01_kind_phys), one)
              !-- for now, set zt = 0.05
              fac = log((2._kind_phys +.05_kind_phys)/(0.05_kind_phys + &
     &            ff))/log((dz1 + .05_kind_phys)/(0.05_kind_phys + ff))
              q2_alt = qsfcmr + fac * dq ! mix. ratio
              q2_alt = q2_alt/(one + q2_alt) ! spec. humidity
            ELSE
            !no alternatives (yet) for unstable conditions
              q2_alt = q2m(i)
            ENDIF
            !-- Note: use of alternative diagnostics will make 
            !   it cooler and drier with stable stratification
            t2m(i) = t2_alt  
            q2m(i) = q2_alt
           endif ! log method for stable regime
 
           !-- check that T2m values lie in the range between tskin and t1
           x2m     = max(min(tskin(i),t1(i)) , t2m(i))
           t2m(i)  = min(max(tskin(i),t1(i)) , x2m)
           !-- check that Q2m values lie in the range between qsurf and q1
           x2m    = max(min(qsurf(i),q1c) , q2m(i))
           q2m(i) = min(max(qsurf(i),q1c) , x2m)
 
           !-- make sure q2m is not oversaturated
           qss    = fpvs(t2m(i))
           qss    = eps * qss/(ps(i) + epsm1 * qss)
           q2m(i) = min(q2m(i),qss)
 
           if(diag_flux) then
         !-- from WRF, applied in HRRR - Jimy Dudhia
         !   Limit Q2m diagnostic to no more than 5 percent higher than lowest level value
         !   This prevents unrealistic values when QFX is not mostly surface
         !   flux because calculation is based on surface flux only.
         !   Problems occurred in transition periods and weak winds and vegetation source
            q2m(i) = min(q2m(i),1.05_kind_dbl_prec*q1c) ! works if qsurf > q1c, evaporation
           endif
         endif ruc_have_2m
 
 
         !-- Compute dew point, using vapor pressure
         qv  = max(qmin,(q2m(i)/(1.-q2m(i))))
         tem = max(ps(i) * qv/( eps - epsm1 *qv), qmin)
         dpt2m(i) = 243.5_kind_dbl_prec/( ( 17.67_kind_dbl_prec /       &
     &             log(tem/611.2_kind_dbl_prec) ) - one) + con_t0c
         dpt2m(i) = min(dpt2m(i),t2m(i))
       
 
         if (debug_print) then
         !-- diagnostics for a test point with known lat/lon
           if (abs(xlat_d(i)-testptlat).lt.0.2 .and.                    &
     &      abs(xlon_d(i)-testptlon).lt.0.2)then
            print 100,'(ruc_lsm_diag)  i=',i,                           &
     &      '  lat,lon=',xlat_d(i),xlon_d(i),'zf ',zf(i),               &
     &      'tskin ',tskin(i),'t2m ',t2m(i),'t1',t1(i),'shflx',shflx(i),&
     &      'qsurf ',qsurf(i),'qsfcprox ',qsfcprox,'q2m ',q2m(i),       &
     &      'q1 ',q1(i),'evap ',evap(i),'dpt2m ',dpt2m(i),              &
     &      'chs2 ',chs2,'cqs2 ',cqs2,'cqs ',cqs,'cdq',cdq(i)
           endif
         endif
 100    format (";;; ",a,i4,a,2f14.7/(4(a10,'='es11.4)))
      endif ! RUC LSM
      enddo
 
      return
      end subroutine sfc_diag_run
 
      end module sfc_diag