sfc_diff.f

 
 
      module sfc_diff
 
      use machine , only : kind_phys
 
      implicit none
 
      public :: sfc_diff_run
      public :: stability
 
      private
 
      real (kind=kind_phys), parameter :: ca=0.4_kind_phys  ! ca - von karman constant
 
      contains
 
      subroutine sfc_diff_run (im,rvrdm1,eps,epsm1,grav,                &  !intent(in)
     &                    ps,t1,q1,z1,garea,wind,                       &  !intent(in)
     &                    prsl1,prslki,prsik1,prslk1,                   &  !intent(in)
     &                    sigmaf,vegtype,shdmax,ivegsrc,                &  !intent(in)
     &                    z0pert,ztpert,                                &  ! mg, sfc-perts !intent(in)
     &                    flag_iter,redrag,                             &  !intent(in)
     &                    flag_lakefreeze,                              &  !intent(in)             
     &                    u10m,v10m,sfc_z0_type,                        &  !hafs,z0 type !intent(in)
     &                    u1,v1,usfco,vsfco,icplocn2atm,                &  
     &                    wet,dry,icy,                                  &  !intent(in)
     &                    thsfc_loc,                                    &  !intent(in)
     &                    tskin_wat, tskin_lnd, tskin_ice,              &  !intent(in)
     &                    tsurf_wat, tsurf_lnd, tsurf_ice,              &  !intent(in)
     &                     z0rl_wat,  z0rl_lnd,  z0rl_ice,              &  !intent(inout)
     &                     z0rl_wav,                                    &  !intent(inout)
     &                    ustar_wat, ustar_lnd, ustar_ice,              &  !intent(inout)
     &                       cm_wat,    cm_lnd,    cm_ice,              &  !intent(inout)
     &                       ch_wat,    ch_lnd,    ch_ice,              &  !intent(inout)
     &                       rb_wat,    rb_lnd,    rb_ice,              &  !intent(inout)
     &                   stress_wat,stress_lnd,stress_ice,              &  !intent(inout)
     &                       fm_wat,    fm_lnd,    fm_ice,              &  !intent(inout)
     &                       fh_wat,    fh_lnd,    fh_ice,              &  !intent(inout)
     &                     fm10_wat,  fm10_lnd,  fm10_ice,              &  !intent(inout)
     &                      fh2_wat,   fh2_lnd,   fh2_ice,              &  !intent(inout)
     &                    ztmax_wat, ztmax_lnd, ztmax_ice,              &  !intent(inout)
     &                    zvfun,                                        &  !intent(out)
     &                    errmsg, errflg)                                  !intent(out)
!
      implicit none
!
      integer, parameter  :: kp = kind_phys
      integer, intent(in) :: im, ivegsrc
      integer, intent(in) :: sfc_z0_type ! option for calculating surface roughness length over ocean
      integer, intent(in) :: icplocn2atm ! option for including ocean current in the computation of flux
 
      integer, dimension(:), intent(in) :: vegtype
 
      logical, intent(in) :: redrag ! reduced drag coeff. flag for high wind over sea (j.han)
      logical, dimension(:), intent(in) :: flag_iter, dry, icy
      logical, dimension(:), intent(in) :: flag_lakefreeze 
      logical, dimension(:), intent(inout) :: wet
 
      logical, intent(in) :: thsfc_loc ! Flag for reference pressure in theta calculation
 
      real(kind=kind_phys), dimension(:), intent(in)    :: u10m,v10m
      real(kind=kind_phys), dimension(:), intent(in)    :: u1,v1
      real(kind=kind_phys), dimension(:), intent(in)    :: usfco,vsfco
      real(kind=kind_phys), intent(in) :: rvrdm1, eps, epsm1, grav
      real(kind=kind_phys), dimension(:), intent(in)    ::              &
     &                    ps,t1,q1,z1,garea,prsl1,prslki,prsik1,prslk1, &
     &                    wind,sigmaf,shdmax,                           &
     &                    z0pert,ztpert ! mg, sfc-perts
      real(kind=kind_phys), dimension(:), intent(in)    ::              &
     &                    tskin_wat, tskin_lnd, tskin_ice,              &
     &                    tsurf_wat, tsurf_lnd, tsurf_ice
 
      real(kind=kind_phys), dimension(:), intent(in)    :: z0rl_wav
      real(kind=kind_phys), dimension(:), intent(inout) ::              &
     &                     z0rl_wat,  z0rl_lnd,  z0rl_ice,              &
     &                    ustar_wat, ustar_lnd, ustar_ice,              &
     &                       cm_wat,    cm_lnd,    cm_ice,              &
     &                       ch_wat,    ch_lnd,    ch_ice,              &
     &                       rb_wat,    rb_lnd,    rb_ice,              &
     &                   stress_wat,stress_lnd,stress_ice,              &
     &                       fm_wat,    fm_lnd,    fm_ice,              &
     &                       fh_wat,    fh_lnd,    fh_ice,              &
     &                     fm10_wat,  fm10_lnd,  fm10_ice,              &
     &                      fh2_wat,   fh2_lnd,   fh2_ice,              &
     &                    ztmax_wat, ztmax_lnd, ztmax_ice
      real(kind=kind_phys), dimension(:), intent(out)    :: zvfun
!
      character(len=*), intent(out) :: errmsg
      integer,          intent(out) :: errflg
!
!     locals
!
      integer   i
      real(kind=kind_phys)  :: windrel
!
      real(kind=kind_phys) :: rat, tv1, thv1, restar, wind10m,
     &                        czilc, tem1, tem2, virtfac
!
 
      real(kind=kind_phys) :: tvs, z0, z0max, ztmax, gdx
!
      real(kind=kind_phys), parameter :: z0lo=0.1, z0up=1.0
!
      real(kind=kind_phys), parameter ::
     &        one=1.0_kp, zero=0.0_kp, half=0.5_kp, qmin=1.0e-8_kp
     &,       charnock=.018_kp, z0s_max=.317e-2_kp                      &! a limiting value at high winds over sea
     &,       zmin=1.0e-6_kp                                            &
     &,       vis=1.4e-5_kp, rnu=1.51e-5_kp, visi=one/vis               &
     &,       log01=log(0.01_kp), log05=log(0.05_kp), log07=log(0.07_kp)
 
!     parameter (charnock=.014,ca=.4)!c ca is the von karman constant
!     parameter (alpha=5.,a0=-3.975,a1=12.32,b1=-7.755,b2=6.041)
!     parameter (a0p=-7.941,a1p=24.75,b1p=-8.705,b2p=7.899,vis=1.4e-5)
 
!     real(kind=kind_phys) aa1,bb1,bb2,cc,cc1,cc2,arnu
!     parameter (aa1=-1.076,bb1=.7045,cc1=-.05808)
!     parameter (bb2=-.1954,cc2=.009999)
!     parameter (arnu=.135*rnu)
!
!    z0s_max=.196e-2 for u10_crit=25 m/s
!    z0s_max=.317e-2 for u10_crit=30 m/s
!    z0s_max=.479e-2 for u10_crit=35 m/s
!
! mbek -- toga-coare flux algorithm
!     parameter (rnu=1.51e-5,arnu=0.11*rnu)
 
! Initialize CCPP error handling variables
      errmsg = ''
      errflg = 0
 
!  initialize variables. all units are supposedly m.k.s. unless specified
!  ps is in pascals, wind is wind speed,
!  surface roughness length is converted to m from cm
!
!       write(0,*)'in sfc_diff, sfc_z0_type=',sfc_z0_type
 
      do i=1,im
        if(flag_iter(i) .or. flag_lakefreeze(i)) then
 
          ! Need to initialize ztmax arrays
          ztmax_lnd(i) = 1. ! log(1) = 0
          ztmax_ice(i) = 1. ! log(1) = 0
          ztmax_wat(i) = 1. ! log(1) = 0
 
          virtfac = one + rvrdm1 * max(q1(i),qmin)
 
          tv1 = t1(i) * virtfac ! Virtual temperature in middle of lowest layer
          if(thsfc_loc) then ! Use local potential temperature
            thv1  = t1(i) * prslki(i) * virtfac
          else ! Use potential temperature reference to 1000 hPa
            thv1    = t1(i) / prslk1(i) * virtfac
          endif
 
          zvfun(i) = zero
          gdx = sqrt(garea(i))
 
!  compute stability dependent exchange coefficients
!  this portion of the code is presently suppressed
!
          if (dry(i)) then ! Some land
 
            if(thsfc_loc) then ! Use local potential temperature
              tvs   = half * (tsurf_lnd(i)+tskin_lnd(i)) * virtfac
            else ! Use potential temperature referenced to 1000 hPa
              tvs   = half * (tsurf_lnd(i)+tskin_lnd(i))/prsik1(i)
     &                     * virtfac
            endif
 
            z0max = max(zmin, min(0.01_kp * z0rl_lnd(i), z1(i)))
!** xubin's new z0  over land
            tem1  = one - shdmax(i)
            tem2  = tem1 * tem1
            tem1  = one  - tem2
 
            if( ivegsrc == 1 ) then
 
              if (vegtype(i) == 10) then
                z0max = exp( tem2*log01 + tem1*log07 )
              elseif (vegtype(i) == 6) then
                z0max = exp( tem2*log01 + tem1*log05 )
              elseif (vegtype(i) == 7) then
!               z0max = exp( tem2*log01 + tem1*log01 )
                z0max = 0.01_kp
              elseif (vegtype(i) == 16) then
!               z0max = exp( tem2*log01 + tem1*log01 )
                z0max = 0.01_kp
              else
                z0max = exp( tem2*log01 + tem1*log(z0max) )
              endif
 
            elseif (ivegsrc == 2 ) then
 
              if (vegtype(i) == 7) then
                z0max = exp( tem2*log01 + tem1*log07 )
              elseif (vegtype(i) == 8) then
                z0max = exp( tem2*log01 + tem1*log05 )
              elseif (vegtype(i) == 9) then
!               z0max = exp( tem2*log01 + tem1*log01 )
                z0max = 0.01_kp
              elseif (vegtype(i) == 11) then
!               z0max = exp( tem2*log01 + tem1*log01 )
                z0max = 0.01_kp
              else
                z0max = exp( tem2*log01 + tem1*log(z0max) )
              endif
 
            endif
! mg, sfc-perts: add surface perturbations to z0max over land
            if (z0pert(i) /= zero ) then
              z0max = z0max * (10.0_kp**z0pert(i))
            endif
 
            z0max = max(z0max, zmin)
 
!!          czilc = 10.0 ** (- (0.40/0.07) * z0) ! fei's canopy height dependance of czil
!           czilc = 0.8_kp
!
!           tem1  = 1.0_kp - sigmaf(i)
!           ztmax_lnd(i) = z0max*exp( - tem1*tem1
!    &              * czilc*ca*sqrt(ustar_lnd(i)*(0.01/1.5e-05)))
!
            czilc = 10.0_kp ** (- 4.0_kp * z0max) ! Trier et al. (2011,WAF)
            czilc = max(min(czilc, 0.8_kp), 0.08_kp)
            tem1 = 1.0_kp - sigmaf(i)
            czilc = czilc * tem1 * tem1
            ztmax_lnd(i) = z0max * exp( - czilc * ca
     &            * 258.2_kp * sqrt(ustar_lnd(i)*z0max) )
!
! mg, sfc-perts: add surface perturbations to ztmax/z0max ratio over land
            if (ztpert(i) /= zero) then
              ztmax_lnd(i) = ztmax_lnd(i) * (10.0_kp**ztpert(i))
            endif
            ztmax_lnd(i) = max(ztmax_lnd(i), zmin)
!
! compute a function of surface roughness & green vegetation fraction (zvfun)       
!
            tem1 = (z0max - z0lo) / (z0up - z0lo)
            tem1 = min(max(tem1, zero), 1.0_kp)
            tem2 = max(sigmaf(i), 0.1_kp)
            zvfun(i) = sqrt(tem1 * tem2)
!
            call stability
!  ---  inputs:
     &       (z1(i), zvfun(i), gdx, tv1, thv1, wind(i),
     &        z0max, ztmax_lnd(i), tvs, grav, thsfc_loc,
!  ---  outputs:
     &        rb_lnd(i), fm_lnd(i), fh_lnd(i), fm10_lnd(i), fh2_lnd(i),
     &        cm_lnd(i), ch_lnd(i), stress_lnd(i), ustar_lnd(i))
          endif ! Dry points
 
          if (icy(i)) then ! Some ice
 
            zvfun(i) = zero
 
            if(thsfc_loc) then ! Use local potential temperature
              tvs   = half * (tsurf_ice(i)+tskin_ice(i)) * virtfac
            else ! Use potential temperature referenced to 1000 hPa
              tvs   = half * (tsurf_ice(i)+tskin_ice(i))/prsik1(i)
     &                     * virtfac 
            endif
 
            z0max = max(zmin, min(0.01_kp * z0rl_ice(i), z1(i)))
!** xubin's new z0  over land and sea ice
            tem1  = one - shdmax(i)
            tem2  = tem1 * tem1
            tem1  = one  - tem2
 
! Removed the following lines by W. Zheng, for effective z0m (z0max) is applied only
! for land.
!wz         if( ivegsrc == 1 ) then
!wz
!wz           z0max = exp( tem2*log01 + tem1*log(z0max) )
!wz         elseif (ivegsrc == 2 ) then
!wz           z0max = exp( tem2*log01 + tem1*log(z0max) )
!wz         endif
 
            z0max = max(z0max, zmin)
 
!!          czilc = 10.0 ** (- (0.40/0.07) * z0) ! fei's canopy height
!           dependance of czil
!           czilc = 0.8_kp
!
!           tem1  = 1.0_kp - sigmaf(i)
!           ztmax_ice(i) = z0max*exp( - tem1*tem1
!    &              * czilc*ca*sqrt(ustar_ice(i)*(0.01/1.5e-05)))
!
            czilc = 10.0_kp ** (- 4.0_kp * z0max)
            czilc = max(min(czilc, 0.8_kp), 0.08_kp)
            tem1 = 1.0_kp - sigmaf(i)
            czilc = czilc * tem1 * tem1
            ztmax_ice(i) = z0max * exp( - czilc * ca
     &            * 258.2_kp * sqrt(ustar_ice(i)*z0max) )
!
            ztmax_ice(i) = max(ztmax_ice(i), 1.0e-6)
!
            call stability
!  ---  inputs:
     &     (z1(i), zvfun(i), gdx, tv1, thv1, wind(i),
     &      z0max, ztmax_ice(i), tvs, grav, thsfc_loc,
!  ---  outputs:
     &      rb_ice(i), fm_ice(i), fh_ice(i), fm10_ice(i), fh2_ice(i),
     &      cm_ice(i), ch_ice(i), stress_ice(i), ustar_ice(i))
      endif ! Icy points
 
! BWG: Everything from here to end of subroutine was after
!      the stuff now put into "stability"
 
          if (wet(i)) then ! Some open ocean
  
            zvfun(i) = zero
 
            if(thsfc_loc) then ! Use local potential temperature
              tvs        = half * (tsurf_wat(i)+tskin_wat(i)) * virtfac
            else
              tvs        = half * (tsurf_wat(i)+tskin_wat(i))/prsik1(i)
     &                          * virtfac
            endif
 
            if (icplocn2atm == 0) then 
              wind10m=sqrt(u10m(i)*u10m(i)+v10m(i)*v10m(i))
              windrel=wind(i)
            else if (icplocn2atm ==1) then
              wind10m=sqrt((u10m(i)-usfco(i))**2+(v10m(i)-vsfco(i))**2)
              windrel=sqrt((u1(i)-usfco(i))**2+(v1(i)-vsfco(i))**2) 
            endif
 
            if (sfc_z0_type == -1) then    ! using wave model derived momentum roughness
              tem1 = 0.11 * vis / ustar_wat(i)
              z0 = tem1 +  0.01_kp * z0rl_wav(i)
 
              if (redrag) then
                z0max = max(min(z0, z0s_max),1.0e-7_kp)
              else
                z0max = max(min(z0,0.1_kp), 1.0e-7_kp)
              endif
              z0rl_wat(i) = 100.0_kp * z0max   ! cm
            else
              z0    = 0.01_kp * z0rl_wat(i)
              z0max = max(zmin, min(z0,z1(i)))
            endif
!
!**  test xubin's new z0
 
!           ztmax  = z0max
 
            restar = max(ustar_wat(i)*z0max*visi, 0.000001_kp)
 
!           restar = log(restar)
!           restar = min(restar,5.)
!           restar = max(restar,-5.)
!           rat    = aa1 + (bb1 + cc1*restar) * restar
!           rat    = rat    / (1. + (bb2 + cc2*restar) * restar))
!  rat taken from zeng, zhao and dickinson 1997
 
            rat   = min(7.0_kp, 2.67_kp * sqrt(sqrt(restar)) - 2.57_kp)
            ztmax_wat(i) = max(z0max * exp(-rat), zmin)
!
            if (sfc_z0_type == 6) then
              call znot_t_v6(wind10m, ztmax_wat(i))   ! 10-m wind,m/s, ztmax(m)
            else if (sfc_z0_type == 7) then
              call znot_t_v7(wind10m, ztmax_wat(i))   ! 10-m wind,m/s, ztmax(m)
            else if (sfc_z0_type > 0) then
              write(0,*)'no option for sfc_z0_type=',sfc_z0_type
              errflg = 1
              errmsg = 'ERROR(sfc_diff_run): no option for sfc_z0_type'
              return
            endif
!
            call stability
!  ---  inputs:
     &       (z1(i), zvfun(i), gdx, tv1, thv1, windrel,
     &        z0max, ztmax_wat(i), tvs, grav, thsfc_loc,
!  ---  outputs:
     &        rb_wat(i), fm_wat(i), fh_wat(i), fm10_wat(i), fh2_wat(i),
     &        cm_wat(i), ch_wat(i), stress_wat(i), ustar_wat(i))
!
!  update z0 over ocean
!
            if (sfc_z0_type >= 0) then
              if (sfc_z0_type == 0) then
!               z0 = (charnock / grav) * ustar_wat(i) * ustar_wat(i)
                tem1 = 0.11 * vis / ustar_wat(i)
                z0 = tem1 + (charnock/grav)*ustar_wat(i)*ustar_wat(i)
 
 
! mbek -- toga-coare flux algorithm
!               z0 = (charnock / grav) * ustar(i)*ustar(i) +  arnu/ustar(i)
!  new implementation of z0
!               cc = ustar(i) * z0 / rnu
!               pp = cc / (1. + cc)
!               ff = grav * arnu / (charnock * ustar(i) ** 3)
!               z0 = arnu / (ustar(i) * ff ** pp)
 
                if (redrag) then
                  z0rl_wat(i) = 100.0_kp * max(min(z0, z0s_max),        &
     &                                                 1.0e-7_kp)
                else
                  z0rl_wat(i) = 100.0_kp * max(min(z0,0.1_kp), 1.e-7_kp)
                endif
 
              elseif (sfc_z0_type == 6) then   ! wang
                 call znot_m_v6(wind10m, z0)   ! wind, m/s, z0, m
                 z0rl_wat(i) = 100.0_kp * z0   ! cm
              elseif (sfc_z0_type == 7) then   ! wang
                 call znot_m_v7(wind10m, z0)   ! wind, m/s, z0, m
                 z0rl_wat(i) = 100.0_kp * z0   ! cm
              else
                 z0rl_wat(i) = 1.0e-4_kp
              endif
 
           elseif (z0rl_wav(i) <= 1.0e-7_kp .or.
     &             z0rl_wav(i) > 1.0_kp) then
!            z0 = (charnock / grav) * ustar_wat(i) * ustar_wat(i)
             tem1 = 0.11 * vis / ustar_wat(i)
             z0 = tem1 + (charnock/grav)*ustar_wat(i)*ustar_wat(i)
 
             if (redrag) then
               z0rl_wat(i) = 100.0_kp * max(min(z0, z0s_max),1.0e-7_kp)
             else
               z0rl_wat(i) = 100.0_kp * max(min(z0,0.1_kp), 1.0e-7_kp)
             endif
 
            endif
 
          endif              ! end of if(open ocean)
!
        endif                ! end of if(flagiter) loop
      enddo
 
      return
      end subroutine sfc_diff_run
 
!----------------------------------------
      subroutine stability                                              &
!  ---  inputs:
     &     ( z1, zvfun, gdx, tv1, thv1, wind, z0max, ztmax, tvs, grav,  &
     &       thsfc_loc,                                                 &
!  ---  outputs:
     &       rb, fm, fh, fm10, fh2, cm, ch, stress, ustar)
!-----
 
      integer, parameter :: kp = kind_phys
!  ---  inputs:
      real(kind=kind_phys), intent(in) ::                               &
     &       z1, zvfun, gdx, tv1, thv1, wind, z0max, ztmax, tvs, grav
      logical,              intent(in) :: thsfc_loc
 
!  ---  outputs:
      real(kind=kind_phys), intent(out) ::                              &
     &       rb, fm, fh, fm10, fh2, cm, ch, stress, ustar
 
!  ---  locals:
      real(kind=kind_phys), parameter :: alpha=5.0_kp, a0=-3.975_kp     &
     &,         a1=12.32_kp,   alpha4=4.0_kp*alpha                      &
     &,         b1=-7.755_kp,  b2=6.041_kp                              &
     &,         xkrefsqr=0.3_kp,  xkmin=0.05_kp                         & 
     &,         xkgdx=3000.0_kp                                         &
     &,         a0p=-7.941_kp, a1p=24.75_kp, b1p=-8.705_kp, b2p=7.899_kp&
     &,         zolmin=-10.0_kp, zero=0.0_kp, one=1.0_kp
 
      real(kind=kind_phys) aa,     aa0,    bb,     bb0, dtv,   adtv,
     &                     hl1,    hl12,   pm,     ph,  pm10,  ph2,
     &                     z1i,
     &                     fms,    fhs,    hl0,    hl0inf, hlinf,
     &                     hl110,  hlt,    hltinf, olinf,
     &                     tem1,   tem2,   zolmax
 
      real(kind=kind_phys) xkzo
 
          z1i = one / z1
 
!
!  set background diffusivities with one for gdx >= xkgdx and 
!   as a function of horizontal grid size for gdx < xkgdx 
!   (i.e., gdx/xkgdx for gdx < xkgdx)
!
          if(gdx >= xkgdx) then
            xkzo = one
          else
            xkzo = gdx / xkgdx
          endif
 
          tem1 = tv1 - tvs
          if(tem1 > zero) then
            tem2 = xkzo * zvfun
            xkzo = min(max(tem2, xkmin), xkzo)
          endif
 
          zolmax = xkrefsqr / sqrt(xkzo)
 
!  compute stability indices (rb and hlinf)
 
          dtv     = thv1 - tvs
          adtv    = max(abs(dtv),0.001_kp)
          dtv     = sign(1.0_kp,dtv) * adtv
 
          if(thsfc_loc) then ! Use local potential temperature
            rb      = max(-5000.0_kp, (grav+grav) * dtv * z1
     &              / ((thv1 + tvs) * wind * wind))
          else ! Use potential temperature referenced to 1000 hPa
            rb      = max(-5000.0_kp, grav * dtv * z1
     &              / (tv1 * wind * wind))
          endif
 
          tem1    = one / z0max
          tem2    = one / ztmax
          fm      = log((z0max+z1)  * tem1)
          fh      = log((ztmax+z1)  * tem2)
          fm10    = log((z0max+10.0_kp) * tem1)
          fh2     = log((ztmax+2.0_kp)  * tem2)
          hlinf   = rb * fm * fm / fh
          hlinf   = min(max(hlinf,zolmin),zolmax)
!
!  stable case
!
          if (dtv >= zero) then
            hl1 = hlinf
            if(hlinf > 0.25_kp) then
              tem1   = hlinf * z1i
              hl0inf = z0max * tem1
              hltinf = ztmax * tem1
              aa     = sqrt(one + alpha4 * hlinf)
              aa0    = sqrt(one + alpha4 * hl0inf)
              bb     = aa
              bb0    = sqrt(one + alpha4 * hltinf)
              pm     = aa0 - aa + log( (aa + one)/(aa0 + one) )
              ph     = bb0 - bb + log( (bb + one)/(bb0 + one) )
              fms    = fm - pm
              fhs    = fh - ph
              hl1    = fms * fms * rb / fhs
              hl1    = min(hl1, zolmax)
            endif
!
!  second iteration
!
            tem1  = hl1 * z1i
            hl0   = z0max * tem1
            hlt   = ztmax * tem1
            aa    = sqrt(one + alpha4 * hl1)
            aa0   = sqrt(one + alpha4 * hl0)
            bb    = aa
            bb0   = sqrt(one + alpha4 * hlt)
            pm    = aa0 - aa + log( (one+aa)/(one+aa0) )
            ph    = bb0 - bb + log( (one+bb)/(one+bb0) )
            hl110 = hl1 * 10.0_kp * z1i
            aa    = sqrt(one + alpha4 * hl110)
            pm10  = aa0 - aa + log( (one+aa)/(one+aa0) )
            hl12  = (hl1+hl1) * z1i
!           aa    = sqrt(one + alpha4 * hl12)
            bb    = sqrt(one + alpha4 * hl12)
            ph2   = bb0 - bb + log( (one+bb)/(one+bb0) )
!
!  unstable case - check for unphysical obukhov length
!
          else                          ! dtv < 0 case
            olinf = z1 / hlinf
            tem1  = 50.0_kp * z0max
            if(abs(olinf) <= tem1) then
              hlinf = -z1 / tem1
              hlinf = max(hlinf, zolmin)
            endif
!
!  get pm and ph
!
            if (hlinf >= -0.5_kp) then
              hl1   = hlinf
              pm    = (a0  + a1*hl1)  * hl1   / (one+ (b1+b2*hl1)  *hl1)
              ph    = (a0p + a1p*hl1) * hl1   / (one+ (b1p+b2p*hl1)*hl1)
              hl110 = hl1 * 10.0_kp * z1i
              pm10  = (a0 + a1*hl110) * hl110/(one+(b1+b2*hl110)*hl110)
              hl12  = (hl1+hl1) * z1i
              ph2   = (a0p + a1p*hl12) * hl12/(one+(b1p+b2p*hl12)*hl12)
            else                       ! hlinf < 0.05
              hl1   = -hlinf
              tem1  = one / sqrt(hl1)
              pm    = log(hl1) + 2.0_kp * sqrt(tem1) - .8776_kp
              ph    = log(hl1) + 0.5_kp * tem1 + 1.386_kp
!             pm    = log(hl1) + 2.0 * hl1 ** (-.25) - .8776
!             ph    = log(hl1) + 0.5 * hl1 ** (-.5) + 1.386
              hl110 = hl1 * 10.0_kp * z1i
              pm10  = log(hl110) + 2.0_kp/sqrt(sqrt(hl110)) - 0.8776_kp
!             pm10  = log(hl110) + 2. * hl110 ** (-.25) - .8776
              hl12  = (hl1+hl1) * z1i
              ph2   = log(hl12) + 0.5_kp / sqrt(hl12) + 1.386_kp
!             ph2   = log(hl12) + .5 * hl12 ** (-.5) + 1.386
            endif
 
          endif          ! end of if (dtv >= 0 ) then loop
!
!  finish the exchange coefficient computation to provide fm and fh
!
          fm        = fm - pm
          fh        = fh - ph
          fm10      = fm10 - pm10
          fh2       = fh2 - ph2
          cm        = ca * ca / (fm * fm)
          ch        = ca * ca / (fm * fh)
          tem1      = 0.00001_kp/z1
          cm        = max(cm, tem1)
          ch        = max(ch, tem1)
          stress    = cm * wind * wind
          ustar     = sqrt(stress)
 
      return
!.................................
      end subroutine stability
!---------------------------------
 
 
      SUBROUTINE znot_m_v6(uref, znotm)
      use machine , only : kind_phys
      IMPLICIT NONE
! Calculate areodynamical roughness over water with input 10-m wind
! For low-to-moderate winds, try to match the Cd-U10 relationship from COARE V3.5 (Edson et al. 2013)
! For high winds, try to fit available observational data
!
! Bin Liu, NOAA/NCEP/EMC 2017
! 
! uref(m/s)   :   wind speed at 10-m height
! znotm(meter):   areodynamical roughness scale over water
!
 
      REAL(kind=kind_phys), INTENT(IN) :: uref
      REAL(kind=kind_phys), INTENT(OUT):: znotm
      real(kind=kind_phys), parameter  :: p13 = -1.296521881682694e-02,
     &      p12 =  2.855780863283819e-01, p11 = -1.597898515251717e+00,
     &      p10 = -8.396975715683501e+00,
 
     &      p25 =  3.790846746036765e-10, p24 =  3.281964357650687e-09,
     &      p23 =  1.962282433562894e-07, p22 = -1.240239171056262e-06,
     &      p21 =  1.739759082358234e-07, p20 =  2.147264020369413e-05,
 
     &      p35 =  1.840430200185075e-07, p34 = -2.793849676757154e-05,
     &      p33 =  1.735308193700643e-03, p32 = -6.139315534216305e-02,
     &      p31 =  1.255457892775006e+00, p30 = -1.663993561652530e+01,
 
     &      p40 =  4.579369142033410e-04
  
 
       if (uref >= 0.0 .and.  uref <= 6.5 ) then
        znotm = exp(p10 + uref * (p11 + uref * (p12 + uref*p13))) 
       elseif (uref > 6.5 .and. uref <= 15.7) then
        znotm = p20 + uref * (p21 + uref * (p22 + uref * (p23
     &              + uref * (p24 + uref * p25))))
       elseif (uref > 15.7 .and. uref <= 53.0) then
        znotm = exp( p30 + uref * (p31 + uref * (p32 + uref * (p33
     &                   + uref * (p34 + uref * p35)))))
       elseif ( uref > 53.0) then
         znotm = p40
       else
        print*, 'Wrong input uref value:',uref
       endif
 
      END SUBROUTINE znot_m_v6
 
!
!! uref(m/s)   :   wind speed at 10-m height
!! znott(meter):   scalar roughness scale over water
      SUBROUTINE znot_t_v6(uref, znott)
      use machine , only : kind_phys
      IMPLICIT NONE
 
 
 
      REAL(kind=kind_phys), INTENT(IN) :: uref
      REAL(kind=kind_phys), INTENT(OUT):: znott
      real(kind=kind_phys), parameter  :: p00 =  1.100000000000000e-04,
     &      p15 = -9.144581627678278e-10, p14 =  7.020346616456421e-08,
     &      p13 = -2.155602086883837e-06, p12 =  3.333848806567684e-05,
     &      p11 = -2.628501274963990e-04, p10 =  8.634221567969181e-04,
 
     &      p25 = -8.654513012535990e-12, p24 =  1.232380050058077e-09,
     &      p23 = -6.837922749505057e-08, p22 =  1.871407733439947e-06,
     &      p21 = -2.552246987137160e-05, p20 =  1.428968311457630e-04,
 
     &      p35 =  3.207515102100162e-12, p34 = -2.945761895342535e-10,
     &      p33 =  8.788972147364181e-09, p32 = -3.814457439412957e-08,
     &      p31 = -2.448983648874671e-06, p30 =  3.436721779020359e-05,
 
     &      p45 = -3.530687797132211e-11, p44 =  3.939867958963747e-09,
     &      p43 = -1.227668406985956e-08, p42 = -1.367469811838390e-05,
     &      p41 =  5.988240863928883e-04, p40 = -7.746288511324971e-03,
 
     &      p56 = -1.187982453329086e-13, p55 =  4.801984186231693e-11,
     &      p54 = -8.049200462388188e-09, p53 =  7.169872601310186e-07,
     &      p52 = -3.581694433758150e-05, p51 =  9.503919224192534e-04,
     &      p50 = -1.036679430885215e-02,
 
     &      p60 =  4.751256171799112e-05
 
      if (uref >= 0.0 .and. uref < 5.9 ) then
         znott = p00
      elseif (uref >= 5.9 .and. uref <= 15.4) then
         znott = p10 + uref * (p11 + uref * (p12 + uref * (p13
     &               + uref * (p14 + uref * p15))))
      elseif (uref > 15.4 .and. uref <= 21.6) then
         znott = p20 + uref * (p21 + uref * (p22 + uref * (p23 
     &               + uref * (p24 + uref * p25))))
      elseif (uref > 21.6 .and. uref <= 42.2) then
         znott = p30 + uref * (p31 + uref * (p32 + uref * (p33 
     &               + uref * (p34 + uref * p35))))
      elseif ( uref > 42.2 .and. uref <= 53.3) then
         znott = p40 + uref * (p41 + uref * (p42 + uref * (p43 
     &               + uref * (p44 + uref * p45))))
      elseif ( uref > 53.3 .and. uref <= 80.0) then
         znott = p50 + uref * (p51 + uref * (p52 + uref * (p53 
     &               + uref * (p54 + uref * (p55 + uref * p56)))))
      elseif ( uref > 80.0) then
         znott = p60
      else
         print*, 'Wrong input uref value:',uref
      endif
 
      END SUBROUTINE znot_t_v6
 
 
!
!! Bin Liu, NOAA/NCEP/EMC 2018
!
!! uref(m/s)   :   wind speed at 10-m height
!! znotm(meter):   areodynamical roughness scale over water
      SUBROUTINE znot_m_v7(uref, znotm)
      use machine , only : kind_phys
      IMPLICIT NONE
 
 
 
      REAL(kind=kind_phys), INTENT(IN) :: uref
      REAL(kind=kind_phys), INTENT(OUT):: znotm
 
      real(kind=kind_phys), parameter  :: p13 = -1.296521881682694e-02,
     &      p12 =  2.855780863283819e-01, p11 = -1.597898515251717e+00,
     &      p10 = -8.396975715683501e+00,
 
     &      p25 =  3.790846746036765e-10, p24 =  3.281964357650687e-09,
     &      p23 =  1.962282433562894e-07, p22 = -1.240239171056262e-06,
     &      p21 =  1.739759082358234e-07, p20 =  2.147264020369413e-05,
 
     &      p35 =  1.897534489606422e-07, p34 = -3.019495980684978e-05,
     &      p33 =  1.931392924987349e-03, p32 = -6.797293095862357e-02,
     &      p31 =  1.346757797103756e+00, p30 = -1.707846930193362e+01,
 
     &      p40 =  3.371427455376717e-04
 
      if (uref >= 0.0 .and.  uref <= 6.5 ) then
        znotm = exp( p10 + uref * (p11 + uref * (p12 + uref * p13)))
      elseif (uref > 6.5 .and. uref <= 15.7) then
        znotm = p20 + uref * (p21 + uref * (p22 + uref * (p23
     &              + uref * (p24 + uref * p25))))
      elseif (uref > 15.7 .and. uref <= 53.0) then
        znotm = exp( p30 + uref * (p31 + uref * (p32 + uref * (p33
     &                   + uref * (p34 + uref * p35)))))
      elseif ( uref > 53.0) then
        znotm = p40
      else
        print*, 'Wrong input uref value:',uref
      endif
 
      END SUBROUTINE znot_m_v7
 
      SUBROUTINE znot_t_v7(uref, znott)
      use machine , only : kind_phys
      IMPLICIT NONE
 
!
 
      REAL(kind=kind_phys), INTENT(IN) :: uref
      REAL(kind=kind_phys), INTENT(OUT):: znott
 
      real(kind=kind_phys), parameter  :: p00 =  1.100000000000000e-04,
 
     &      p15 = -9.193764479895316e-10, p14 =  7.052217518653943e-08,
     &      p13 = -2.163419217747114e-06, p12 =  3.342963077911962e-05,
     &      p11 = -2.633566691328004e-04, p10 =  8.644979973037803e-04,
 
     &      p25 = -9.402722450219142e-12, p24 =  1.325396583616614e-09,
     &      p23 = -7.299148051141852e-08, p22 =  1.982901461144764e-06,
     &      p21 = -2.680293455916390e-05, p20 =  1.484341646128200e-04,
 
     &      p35 =  7.921446674311864e-12, p34 = -1.019028029546602e-09,
     &      p33 =  5.251986927351103e-08, p32 = -1.337841892062716e-06,
     &      p31 =  1.659454106237737e-05, p30 = -7.558911792344770e-05,
 
     &      p45 = -2.694370426850801e-10, p44 =  5.817362913967911e-08,
     &      p43 = -5.000813324746342e-06, p42 =  2.143803523428029e-04,
     &      p41 = -4.588070983722060e-03, p40 =  3.924356617245624e-02,
 
     &      p56 = -1.663918773476178e-13, p55 =  6.724854483077447e-11,
     &      p54 = -1.127030176632823e-08, p53 =  1.003683177025925e-06,
     &      p52 = -5.012618091180904e-05, p51 =  1.329762020689302e-03,
     &      p50 = -1.450062148367566e-02, p60 =  6.840803042788488e-05
 
        if (uref >= 0.0 .and. uref < 5.9 ) then
           znott = p00
        elseif (uref >= 5.9 .and. uref <= 15.4) then
           znott = p10 + uref * (p11 + uref * (p12 + uref * (p13
     &                 + uref * (p14 + uref * p15))))
        elseif (uref > 15.4 .and. uref <= 21.6) then
           znott = p20 + uref * (p21 + uref * (p22 + uref * (p23
     &                 + uref * (p24 + uref * p25))))
        elseif (uref > 21.6 .and. uref <= 42.6) then
           znott = p30 + uref * (p31 + uref * (p32 + uref * (p33
     &                 + uref * (p34 + uref * p35))))
        elseif ( uref > 42.6 .and. uref <= 53.0) then
           znott = p40 + uref * (p41 + uref * (p42 + uref * (p43
     &                 + uref * (p44 + uref * p45))))
        elseif ( uref > 53.0 .and. uref <= 80.0) then
           znott = p50 + uref * (p51 + uref * (p52 + uref * (p53
     &                 + uref * (p54 + uref * (p55 + uref * p56)))))
        elseif ( uref > 80.0) then
           znott = p60
        else
           print*, 'Wrong input uref value:',uref
        endif
 
        END SUBROUTINE znot_t_v7
      end module sfc_diff