CESM_Driver: shr_flux_mod.F90 Source File

00001 !===============================================================================
00002 ! SVN $Id: shr_flux_mod.F90 19407 2009-11-13 19:36:31Z tcraig $
00003 ! SVN $URL: https://svn-ccsm-models.cgd.ucar.edu/csm_share/branch_tags/cesm1_0_rel_tags/cesm1_0_rel01_share3_100616/shr/shr_flux_mod.F90 $
00004 !===============================================================================
00005 !BOP ===========================================================================
00006 !
00007 ! !MODULE: flux_mod -- CCSM shared flux calculations.
00008 !
00009 ! !DESCRIPTION:
00010 !
00011 !     CCSM shared flux calculations.
00012 !     
00013 ! !REVISION HISTORY:
00014 !     2006-Nov-07 - B. Kauffman - first version, code taken/migrated from cpl6
00015 !
00016 ! !INTERFACE: ------------------------------------------------------------------
00017 
00018 module shr_flux_mod
00019 
00020 ! !USES:
00021 
00022    use shr_kind_mod    ! shared kinds
00023    use shr_const_mod   ! shared constants
00024    use shr_sys_mod     ! shared system routines
00025    use shr_log_mod, only: s_loglev  => shr_log_Level
00026    use shr_log_mod, only: s_logunit => shr_log_Unit
00027 
00028    implicit none
00029 
00030    private ! default private
00031 
00032 ! !PUBLIC TYPES:
00033 
00034   ! none
00035 
00036 ! !PUBLIC MEMBER FUNCTIONS:
00037 
00038    public :: shr_flux_atmOcn      ! computes atm/ocn fluxes
00039    public :: shr_flux_atmIce      ! computes atm/ice fluxes
00040    public :: shr_flux_MOstability ! boundary layer stability scales/functions
00041 
00042 ! !PUBLIC DATA MEMBERS:
00043 
00044   integer(SHR_KIND_IN),parameter,public :: shr_flux_MOwScales   = 1 ! w scales  option
00045   integer(SHR_KIND_IN),parameter,public :: shr_flux_MOfunctions = 2 ! functions option
00046   real   (SHR_KIND_R8),parameter,public :: shr_flux_MOgammaM = 3.59_SHR_KIND_R8
00047   real   (SHR_KIND_R8),parameter,public :: shr_flux_MOgammaS = 7.86_SHR_KIND_R8
00048 
00049 !EOP
00050 
00051    !--- rename kinds for local readability only ---
00052    integer,parameter :: R8 = SHR_KIND_R8  ! 8 byte real
00053    integer,parameter :: IN = SHR_KIND_IN  ! native/default integer
00054 
00055    integer,parameter :: debug = 0 ! internal debug level
00056 
00057 !===============================================================================
00058 contains
00059 !===============================================================================
00060 !===============================================================================
00061 ! !BOP =========================================================================
00062 !
00063 ! !IROUTINE: shr_flux_atmOcn -- internal atm/ocn flux calculation
00064 !
00065 ! !DESCRIPTION:
00066 !
00067 !     Internal atm/ocn flux calculation
00068 !     
00069 ! !REVISION HISTORY:
00070 !     2002-Jun-10 - B. Kauffman - code migrated from cpl5 to cpl6
00071 !     2003-Apr-02 - B. Kauffman - taux & tauy now utilize ocn velocity
00072 !     2003-Apr-02 - B. Kauffman - tref,qref,duu10n mods as per Bill Large
00073 !     2006-Nov-07 - B. Kauffman - code migrated from cpl6 to share
00074 !
00075 ! !INTERFACE: ------------------------------------------------------------------
00076 
00077 SUBROUTINE shr_flux_atmOcn(nMax  ,zbot  ,ubot  ,vbot  ,thbot ,   & 
00078            &               qbot  ,rbot  ,tbot  ,us    ,vs    ,   &
00079            &               ts    ,mask  ,sen   ,lat   ,lwup  ,   &
00080            &               evap  ,taux  ,tauy  ,tref  ,qref  ,   &
00081            &               duu10n,  ustar_sv   ,re_sv ,ssq_sv    )
00082 
00083 ! !USES:
00084 
00085    implicit none
00086 
00087 ! !INPUT/OUTPUT PARAMETERS:
00088 
00089    !--- input arguments --------------------------------
00090    integer(IN),intent(in) ::       nMax  ! data vector length
00091    integer(IN),intent(in) :: mask (nMax) ! ocn domain mask       0 <=> out of domain
00092    real(R8)   ,intent(in) :: zbot (nMax) ! atm level height      (m)
00093    real(R8)   ,intent(in) :: ubot (nMax) ! atm u wind            (m/s)
00094    real(R8)   ,intent(in) :: vbot (nMax) ! atm v wind            (m/s)
00095    real(R8)   ,intent(in) :: thbot(nMax) ! atm potential T       (K)
00096    real(R8)   ,intent(in) :: qbot (nMax) ! atm specific humidity (kg/kg)
00097    real(R8)   ,intent(in) :: rbot (nMax) ! atm air density       (kg/m^3)
00098    real(R8)   ,intent(in) :: tbot (nMax) ! atm T                 (K) 
00099    real(R8)   ,intent(in) :: us   (nMax) ! ocn u-velocity        (m/s)
00100    real(R8)   ,intent(in) :: vs   (nMax) ! ocn v-velocity        (m/s)
00101    real(R8)   ,intent(in) :: ts   (nMax) ! ocn temperature       (K)
00102 
00103    !--- output arguments -------------------------------
00104    real(R8),intent(out)  ::  sen  (nMax) ! heat flux: sensible    (W/m^2)
00105    real(R8),intent(out)  ::  lat  (nMax) ! heat flux: latent      (W/m^2)
00106    real(R8),intent(out)  ::  lwup (nMax) ! heat flux: lw upward   (W/m^2)
00107    real(R8),intent(out)  ::  evap (nMax) ! water flux: evap  ((kg/s)/m^2)
00108    real(R8),intent(out)  ::  taux (nMax) ! surface stress, zonal      (N)
00109    real(R8),intent(out)  ::  tauy (nMax) ! surface stress, maridional (N)
00110    real(R8),intent(out)  ::  tref (nMax) ! diag:  2m ref height T     (K)
00111    real(R8),intent(out)  ::  qref (nMax) ! diag:  2m ref humidity (kg/kg)
00112    real(R8),intent(out)  :: duu10n(nMax) ! diag: 10m wind speed squared (m/s)^2
00113 
00114    real(R8),intent(out),optional :: ustar_sv(nMax) ! diag: ustar
00115    real(R8),intent(out),optional :: re_sv   (nMax) ! diag: sqrt of exchange coefficient (water)
00116    real(R8),intent(out),optional :: ssq_sv  (nMax) ! diag: sea surface humidity  (kg/kg)
00117  
00118 ! !EOP
00119 
00120    !--- local constants --------------------------------
00121    real(R8),parameter :: umin  =  0.5_R8 ! minimum wind speed       (m/s)
00122    real(R8),parameter :: zref  = 10.0_R8 ! reference height           (m)
00123    real(R8),parameter :: ztref =  2.0_R8 ! reference height for air T (m)
00124 
00125    !--- local variables --------------------------------
00126    integer(IN) :: n      ! vector loop index
00127    real(R8)    :: vmag   ! surface wind magnitude   (m/s)
00128    real(R8)    :: thvbot ! virtual temperature      (K)
00129    real(R8)    :: ssq    ! sea surface humidity     (kg/kg)
00130    real(R8)    :: delt   ! potential T difference   (K)
00131    real(R8)    :: delq   ! humidity difference      (kg/kg)
00132    real(R8)    :: stable ! stability factor
00133    real(R8)    :: rdn    ! sqrt of neutral exchange coeff (momentum) 
00134    real(R8)    :: rhn    ! sqrt of neutral exchange coeff (heat)     
00135    real(R8)    :: ren    ! sqrt of neutral exchange coeff (water)    
00136    real(R8)    :: rd     ! sqrt of exchange coefficient (momentum)         
00137    real(R8)    :: rh     ! sqrt of exchange coefficient (heat)             
00138    real(R8)    :: re     ! sqrt of exchange coefficient (water)            
00139    real(R8)    :: ustar  ! ustar             
00140    real(R8)    :: qstar  ! qstar             
00141    real(R8)    :: tstar  ! tstar             
00142    real(R8)    :: hol    ! H (at zbot) over L
00143    real(R8)    :: xsq    ! ?
00144    real(R8)    :: xqq    ! ?
00145    real(R8)    :: psimh  ! stability function at zbot (momentum)
00146    real(R8)    :: psixh  ! stability function at zbot (heat and water)
00147    real(R8)    :: psix2  ! stability function at ztref reference height
00148    real(R8)    :: alz    ! ln(zbot/zref)
00149    real(R8)    :: al2    ! ln(zref/ztref)
00150    real(R8)    :: u10n   ! 10m neutral wind 
00151    real(R8)    :: tau    ! stress at zbot
00152    real(R8)    :: cp     ! specific heat of moist air
00153    real(R8)    :: bn     ! exchange coef funct for interpolation
00154    real(R8)    :: bh     ! exchange coef funct for interpolation
00155    real(R8)    :: fac    ! vertical interpolation factor
00156 
00157    !--- local functions --------------------------------
00158    real(R8)    :: qsat   ! function: the saturation humididty of air (kg/m^3)
00159    real(R8)    :: cdn    ! function: neutral drag coeff at 10m
00160    real(R8)    :: psimhu ! function: unstable part of psimh
00161    real(R8)    :: psixhu ! function: unstable part of psimx
00162    real(R8)    :: Umps   ! dummy arg ~ wind velocity (m/s)
00163    real(R8)    :: Tk     ! dummy arg ~ temperature (K)
00164    real(R8)    :: xd     ! dummy arg ~ ?
00165  
00166    qsat(Tk)   = 640380.0_R8 / exp(5107.4_R8/Tk)
00167    cdn(Umps)  =   0.0027_R8 / Umps + 0.000142_R8 + 0.0000764_R8 * Umps
00168    psimhu(xd) = log((1.0_R8+xd*(2.0_R8+xd))*(1.0_R8+xd*xd)/8.0_R8) - 2.0_R8*atan(xd) + 1.571_R8
00169    psixhu(xd) = 2.0_R8 * log((1.0_R8 + xd*xd)/2.0_R8)
00170  
00171    !--- formats ----------------------------------------
00172    character(*),parameter :: subName = '(shr_flux_atmOcn) '
00173    character(*),parameter ::   F00 = "('(shr_flux_atmOcn) ',4a)"
00174 
00175 !-------------------------------------------------------------------------------
00176 ! PURPOSE:
00177 !   computes atm/ocn surface fluxes
00178 !
00179 ! NOTES: 
00180 !   o all fluxes are positive downward
00181 !   o net heat flux = net sw + lw up + lw down + sen + lat
00182 !   o here, tstar = <WT>/U*, and qstar = <WQ>/U*.
00183 !   o wind speeds should all be above a minimum speed (eg. 1.0 m/s)
00184 ! 
00185 ! ASSUMPTIONS:
00186 !   o Neutral 10m drag coeff: cdn = .0027/U10 + .000142 + .0000764 U10
00187 !   o Neutral 10m stanton number: ctn = .0327 sqrt(cdn), unstable
00188 !                                 ctn = .0180 sqrt(cdn), stable
00189 !   o Neutral 10m dalton number:  cen = .0346 sqrt(cdn)
00190 !   o The saturation humidity of air at T(K): qsat(T)  (kg/m^3)
00191 !-------------------------------------------------------------------------------
00192 
00193    if (debug > 0 .and. s_loglev > 0) write(s_logunit,F00) "enter"
00194  
00195    al2 = log(zref/ztref)
00196 
00197    DO n=1,nMax
00198      if (mask(n) /= 0) then
00199     
00200         !--- compute some needed quantities ---
00201         vmag   = max(umin, sqrt( (ubot(n)-us(n))**2 + (vbot(n)-vs(n))**2) )
00202         thvbot = thbot(n) * (1.0_R8 + shr_const_zvir * qbot(n)) ! virtual temp (K)
00203         ssq    = 0.98_R8 * qsat(ts(n)) / rbot(n)   ! sea surf hum (kg/kg)
00204         delt   = thbot(n) - ts(n)                  ! pot temp diff (K)
00205         delq   = qbot(n) - ssq                     ! spec hum dif (kg/kg)
00206         alz    = log(zbot(n)/zref) 
00207         cp     = shr_const_cpdair*(1.0_R8 + shr_const_cpvir*ssq) 
00208    
00209         !------------------------------------------------------------
00210         ! first estimate of Z/L and ustar, tstar and qstar
00211         !------------------------------------------------------------
00212    
00213         !--- neutral coefficients, z/L = 0.0 ---
00214         stable = 0.5_R8 + sign(0.5_R8 , delt)
00215         rdn    = sqrt(cdn(vmag))
00216         rhn    = (1.0_R8-stable) * 0.0327_R8 + stable * 0.018_R8 
00217         ren    = 0.0346_R8 
00218    
00219         !--- ustar, tstar, qstar ---
00220         ustar = rdn * vmag
00221         tstar = rhn * delt  
00222         qstar = ren * delq  
00223    
00224         !--- compute stability & evaluate all stability functions ---
00225         hol  = shr_const_karman*shr_const_g*zbot(n)*  &
00226                (tstar/thbot(n)+qstar/(1.0_R8/shr_const_zvir+qbot(n)))/ustar**2
00227         hol  = sign( min(abs(hol),10.0_R8), hol )
00228         stable = 0.5_R8 + sign(0.5_R8 , hol)
00229         xsq    = max(sqrt(abs(1.0_R8 - 16.0_R8*hol)) , 1.0_R8)
00230         xqq    = sqrt(xsq)
00231         psimh  = -5.0_R8*hol*stable + (1.0_R8-stable)*psimhu(xqq)
00232         psixh  = -5.0_R8*hol*stable + (1.0_R8-stable)*psixhu(xqq)
00233    
00234         !--- shift wind speed using old coefficient ---
00235         rd   = rdn / (1.0_R8 + rdn/shr_const_karman*(alz-psimh))
00236         u10n = vmag * rd / rdn 
00237    
00238         !--- update transfer coeffs at 10m and neutral stability ---
00239         rdn = sqrt(cdn(u10n))
00240         ren = 0.0346_R8
00241         rhn = (1.0_R8-stable)*0.0327_R8 + stable * 0.018_R8 
00242     
00243         !--- shift all coeffs to measurement height and stability ---
00244         rd = rdn / (1.0_R8 + rdn/shr_const_karman*(alz-psimh)) 
00245         rh = rhn / (1.0_R8 + rhn/shr_const_karman*(alz-psixh)) 
00246         re = ren / (1.0_R8 + ren/shr_const_karman*(alz-psixh)) 
00247    
00248         !--- update ustar, tstar, qstar using updated, shifted coeffs --
00249         ustar = rd * vmag 
00250         tstar = rh * delt 
00251         qstar = re * delq 
00252     
00253         !------------------------------------------------------------
00254         ! iterate to converge on Z/L, ustar, tstar and qstar
00255         !------------------------------------------------------------
00256     
00257         !--- compute stability & evaluate all stability functions ---
00258         hol  = shr_const_karman*shr_const_g*zbot(n)* &
00259                (tstar/thbot(n)+qstar/(1.0_R8/shr_const_zvir+qbot(n)))/ustar**2
00260         hol  = sign( min(abs(hol),10.0_R8), hol )
00261         stable = 0.5_R8 + sign(0.5_R8 , hol)
00262         xsq    = max(sqrt(abs(1.0_R8 - 16.0_R8*hol)) , 1.0_R8)
00263         xqq    = sqrt(xsq)
00264         psimh  = -5.0_R8*hol*stable + (1.0_R8-stable)*psimhu(xqq)
00265         psixh  = -5.0_R8*hol*stable + (1.0_R8-stable)*psixhu(xqq)
00266     
00267         !--- shift wind speed using old coeffs ---
00268         rd   = rdn / (1.0_R8 + rdn/shr_const_karman*(alz-psimh))
00269         u10n = vmag * rd/rdn 
00270     
00271         !--- update transfer coeffs at 10m and neutral stability ---
00272         rdn = sqrt(cdn(u10n))
00273         ren = 0.0346_R8
00274         rhn = (1.0_R8 - stable)*0.0327_R8 + stable * 0.018_R8 
00275    
00276         !--- shift all coeffs to measurement height and stability ---
00277         rd = rdn / (1.0_R8 + rdn/shr_const_karman*(alz-psimh)) 
00278         rh = rhn / (1.0_R8 + rhn/shr_const_karman*(alz-psixh)) 
00279         re = ren / (1.0_R8 + ren/shr_const_karman*(alz-psixh)) 
00280     
00281         !--- update ustar, tstar, qstar using updated, shifted coeffs ---
00282         ustar = rd * vmag 
00283         tstar = rh * delt 
00284         qstar = re * delq 
00285     
00286         !------------------------------------------------------------
00287         ! compute the fluxes
00288         !------------------------------------------------------------
00289     
00290         tau = rbot(n) * ustar * ustar 
00291        
00292         !--- momentum flux ---
00293         taux(n) = tau * (ubot(n)-us(n)) / vmag 
00294         tauy(n) = tau * (vbot(n)-vs(n)) / vmag 
00295         
00296         !--- heat flux ---
00297         sen (n) =                cp * tau * tstar / ustar 
00298         lat (n) =  shr_const_latvap * tau * qstar / ustar
00299         lwup(n) = -shr_const_stebol * ts(n)**4 
00300       
00301         !--- water flux ---
00302         evap(n) = lat(n)/shr_const_latvap 
00303     
00304         !------------------------------------------------------------
00305         ! compute diagnositcs: 2m ref T & Q, 10m wind speed squared
00306         !------------------------------------------------------------
00307         hol = hol*ztref/zbot(n)
00308         xsq = max( 1.0_R8, sqrt(abs(1.0_R8-16.0_R8*hol)) )
00309         xqq = sqrt(xsq)
00310         psix2   = -5.0_R8*hol*stable + (1.0_R8-stable)*psixhu(xqq)
00311         fac     = (rh/shr_const_karman) * (alz + al2 - psixh + psix2 )
00312         tref(n) = thbot(n) - delt*fac 
00313         tref(n) = tref(n) - 0.01_R8*ztref   ! pot temp to temp correction
00314         fac     = (re/shr_const_karman) * (alz + al2 - psixh + psix2 )
00315         qref(n) =  qbot(n) - delq*fac
00316     
00317         duu10n(n) = u10n*u10n ! 10m wind speed squared
00318 
00319         !------------------------------------------------------------
00320         ! optional diagnostics, needed for water tracer fluxes (dcn)
00321         !------------------------------------------------------------
00322         if (present(ustar_sv)) ustar_sv(n) = ustar
00323         if (present(re_sv   )) re_sv(n)    = re
00324         if (present(ssq_sv  )) ssq_sv(n)   = ssq
00325 
00326      else
00327         !------------------------------------------------------------
00328         ! no valid data here -- out of domain
00329         !------------------------------------------------------------
00330         sen   (n) = shr_const_spval  ! sensible         heat flux  (W/m^2)
00331         lat   (n) = shr_const_spval  ! latent           heat flux  (W/m^2)
00332         lwup  (n) = shr_const_spval  ! long-wave upward heat flux  (W/m^2)
00333         evap  (n) = shr_const_spval  ! evaporative water flux ((kg/s)/m^2)
00334         taux  (n) = shr_const_spval  ! x surface stress (N)
00335         tauy  (n) = shr_const_spval  ! y surface stress (N)
00336         tref  (n) = shr_const_spval  !  2m reference height temperature (K)
00337         qref  (n) = shr_const_spval  !  2m reference height humidity (kg/kg)
00338         duu10n(n) = shr_const_spval  ! 10m wind speed squared (m/s)^2
00339 
00340         if (present(ustar_sv)) ustar_sv(n) = shr_const_spval
00341         if (present(re_sv   )) re_sv   (n) = shr_const_spval
00342         if (present(ssq_sv  )) ssq_sv  (n) = shr_const_spval
00343      endif
00344    ENDDO 
00345 
00346 END subroutine shr_flux_atmOcn
00347 
00348 !BOP ===========================================================================
00349 !
00350 ! !IROUTINE: shr_flux_atmIce -- computes atm/ice fluxes
00351 !
00352 ! !DESCRIPTION:
00353 !    Computes atm/ice fluxes
00354 !
00355 ! !REVISION HISTORY:
00356 !    2006-Jun-12 - B. Kauffman, first version, adapted from dice6 code
00357 !
00358 ! !INTERFACE: ------------------------------------------------------------------
00359 
00360 subroutine shr_flux_atmIce(mask  ,zbot  ,ubot  ,vbot  ,thbot  &
00361                &          ,qbot  ,rbot  ,tbot  ,ts    ,sen    &
00362                &          ,lat   ,lwup  ,evap  ,taux  ,tauy   &
00363                &          ,tref  ,qref                        )
00364 
00365    implicit none
00366 
00367 ! !INPUT/OUTPUT PARAMETERS:
00368 
00369    !--- input arguments --------------------------------
00370    integer(IN),intent(in)  :: mask (:)    ! 0 <=> cell NOT in model domain
00371    real(R8)   ,intent(in)  :: zbot (:)    ! atm level height  (m)
00372    real(R8)   ,intent(in)  :: ubot (:)    ! atm u wind     (m/s)
00373    real(R8)   ,intent(in)  :: vbot (:)    ! atm v wind     (m/s)
00374    real(R8)   ,intent(in)  :: thbot(:)    ! atm potential T   (K)
00375    real(R8)   ,intent(in)  :: qbot (:)    ! atm specific humidity (kg/kg)
00376    real(R8)   ,intent(in)  :: rbot (:)    ! atm air density   (kg/m^3)
00377    real(R8)   ,intent(in)  :: tbot (:)    ! atm T       (K) 
00378    real(R8)   ,intent(in)  :: ts   (:)    ! surface temperature
00379 
00380    !--- output arguments -------------------------------
00381    real(R8)   ,intent(out) :: sen  (:)    ! sensible      heat flux  (W/m^2)
00382    real(R8)   ,intent(out) :: lat  (:)    ! latent        heat flux  (W/m^2)
00383    real(R8)   ,intent(out) :: lwup (:)    ! long-wave upward heat flux  (W/m^2)
00384    real(R8)   ,intent(out) :: evap (:)    ! evaporative water flux ((kg/s)/m^2)
00385    real(R8)   ,intent(out) :: taux (:)    ! x surface stress (N)
00386    real(R8)   ,intent(out) :: tauy (:)    ! y surface stress (N)
00387    real(R8)   ,intent(out) :: tref (:)    ! 2m reference height temperature
00388    real(R8)   ,intent(out) :: qref (:)    ! 2m reference height humidity
00389  
00390 !EOP
00391 
00392    !--- local constants --------------------------------
00393    real(R8),parameter :: umin   =  1.0_R8            ! minimum wind speed (m/s)
00394    real(R8),parameter :: zref   = 10.0_R8            ! ref height           ~ m
00395    real(R8),parameter :: ztref  =  2.0_R8            ! ref height for air T ~ m
00396    real(R8),parameter :: spval  = shr_const_spval    ! special value
00397    real(R8),parameter :: g      = shr_const_g        ! gravity
00398    real(R8),parameter :: cpdair = shr_const_cpdair   ! spec heat of dry air
00399    real(R8),parameter :: cpvir  = shr_const_cpvir    ! cpwv/cpdair - 1.0
00400    real(R8),parameter :: zvir   = shr_const_zvir     ! rh2o/rair   - 1.0
00401    real(R8),parameter :: latvap = shr_const_latvap   ! latent heat of evap
00402    real(R8),parameter :: latice = shr_const_latice   ! latent heat of fusion
00403    real(R8),parameter :: stebol = shr_const_stebol   ! Stefan-Boltzmann
00404    real(R8),parameter :: karman = shr_const_karman   ! Von Karman constant
00405    real(R8),parameter :: zzsice = 0.0005_R8          ! ice surface roughness
00406 
00407    !--- local variables --------------------------------
00408    integer(IN) :: lsize  ! array dimensions
00409    integer(IN) :: n      ! array indicies
00410    real(R8)    :: vmag   ! surface wind magnitude   (m/s)
00411    real(R8)    :: thvbot ! virtual temperature      (K)
00412    real(R8)    :: ssq    ! sea surface humidity     (kg/kg)
00413    real(R8)    :: dssqdt ! derivative of ssq wrt Ts (kg/kg/K)
00414    real(R8)    :: delt   ! potential T difference   (K)
00415    real(R8)    :: delq   ! humidity difference      (kg/kg)
00416    real(R8)    :: stable ! stability factor
00417    real(R8)    :: rdn    ! sqrt of neutral exchange coefficient (momentum)
00418    real(R8)    :: rhn    ! sqrt of neutral exchange coefficient (heat)
00419    real(R8)    :: ren    ! sqrt of neutral exchange coefficient (water)
00420    real(R8)    :: rd     ! sqrt of exchange coefficient (momentum)
00421    real(R8)    :: rh     ! sqrt of exchange coefficient (heat)
00422    real(R8)    :: re     ! sqrt of exchange coefficient (water)      
00423    real(R8)    :: ustar  ! ustar
00424    real(R8)    :: qstar  ! qstar
00425    real(R8)    :: tstar  ! tstar
00426    real(R8)    :: hol    ! H (at zbot) over L
00427    real(R8)    :: xsq    ! temporary variable
00428    real(R8)    :: xqq    ! temporary variable
00429    real(R8)    :: psimh  ! stability function at zbot (momentum)
00430    real(R8)    :: psixh  ! stability function at zbot (heat and water)
00431    real(R8)    :: alz    ! ln(zbot/z10)
00432    real(R8)    :: ltheat ! latent heat for surface
00433    real(R8)    :: tau    ! stress at zbot
00434    real(R8)    :: cp     ! specific heat of moist air
00435 
00436    real(R8)    :: bn     ! exchange coef funct for interpolation
00437    real(R8)    :: bh     ! exchange coef funct for interpolation
00438    real(R8)    :: fac    ! interpolation factor
00439    real(R8)    :: ln0    ! log factor for interpolation
00440    real(R8)    :: ln3    ! log factor for interpolation
00441 
00442    !--- local functions --------------------------------
00443    real(R8)   :: Tk      ! temperature (K)
00444    real(R8)   :: qsat    ! the saturation humididty of air (kg/m^3)
00445    real(R8)   :: dqsatdt ! derivivative of qsat wrt surface temperature
00446    real(R8)   :: xd      ! dummy argument  
00447    real(R8)   :: psimhu  ! unstable part of psimh
00448    real(R8)   :: psixhu  ! unstable part of psimx
00449 
00450    qsat(Tk)    = 627572.4_R8 / exp(5107.4_R8/Tk)
00451    dqsatdt(Tk) = (5107.4_R8 / Tk**2) * 627572.4_R8 / exp(5107.4_R8/Tk)
00452    psimhu(xd)  = log((1.0_R8+xd*(2.0_R8+xd))*(1.0_R8+xd*xd)/8.0_R8) - 2.0_R8*atan(xd) + 1.571_R8
00453    psixhu(xd)  =  2.0_R8 * log((1.0_R8 + xd*xd)/2.0_R8)
00454 
00455    !--- formats ----------------------------------------
00456    character(*),parameter :: subName =  "(shr_flux_atmIce) "
00457 
00458 !-------------------------------------------------------------------------------
00459 ! PURPOSE:
00460 !   using atm & ice state variables, compute atm/ice fluxes 
00461 !   and diagnostic 10m air temperature and humidity
00462 !
00463 ! NOTE: 
00464 !   o all fluxes are positive downward
00465 !   o net heat flux = net sw + lw up + lw down + sen + lat
00466 !   o here, tstar = <WT>/U*, and qstar = <WQ>/U*.
00467 !   o wind speeds should all be above a minimum speed (eg. 1.0 m/s)
00468 ! 
00469 ! ASSUME:
00470 !   o The saturation humidity of air at T(K): qsat(T)  (kg/m^3)
00471 !-------------------------------------------------------------------------------
00472 
00473    lsize = size(tbot)
00474 
00475    do n = 1,lsize
00476 
00477      if (mask(n) == 0) then
00478         sen  (n) = spval
00479         lat  (n) = spval
00480         lwup (n) = spval
00481         evap (n) = spval
00482         taux (n) = spval
00483         tauy (n) = spval
00484         tref (n) = spval
00485         qref (n) = spval
00486      else
00487         !--- define some needed variables ---
00488         vmag   = max(umin, sqrt(ubot(n)**2+vbot(n)**2))
00489         thvbot = thbot(n)*(1.0_R8 + zvir * qbot(n)) ! virtual pot temp (K)
00490          ssq   =  qsat  (ts(n)) / rbot(n)           ! sea surf hum (kg/kg)
00491         dssqdt = dqsatdt(ts(n)) / rbot(n)           ! deriv of ssq wrt Ts 
00492         delt   = thbot(n) - ts(n)                   ! pot temp diff (K)
00493         delq   = qbot(n) - ssq                        ! spec hum dif (kg/kg)
00494         alz    = log(zbot(n)/zref) 
00495         cp     = cpdair*(1.0_R8 + cpvir*ssq) 
00496         ltheat = latvap + latice
00497 
00498         !----------------------------------------------------------
00499         ! first estimate of Z/L and ustar, tstar and qstar
00500         !----------------------------------------------------------
00501 
00502         !--- neutral coefficients, z/L = 0.0 ---
00503         rdn = karman/log(zref/zzsice)
00504         rhn = rdn
00505         ren = rdn
00506 
00507         !--- ustar,tstar,qstar ----
00508         ustar = rdn * vmag
00509         tstar = rhn * delt  
00510         qstar = ren * delq  
00511 
00512         !--- compute stability & evaluate all stability functions ---
00513         hol    = karman * g * zbot(n) &
00514         &     * (tstar/thvbot+qstar/(1.0_R8/zvir+qbot(n))) / ustar**2
00515         hol    = sign( min(abs(hol),10.0_R8), hol )
00516         stable = 0.5_R8 + sign(0.5_R8 , hol)
00517         xsq    = max(sqrt(abs(1.0_R8 - 16.0_R8*hol)) , 1.0_R8)
00518         xqq    = sqrt(xsq)
00519         psimh  = -5.0_R8*hol*stable + (1.0_R8-stable)*psimhu(xqq)
00520         psixh  = -5.0_R8*hol*stable + (1.0_R8-stable)*psixhu(xqq)
00521 
00522         !--- shift all coeffs to measurement height and stability ---
00523         rd = rdn / (1.0_R8+rdn/karman*(alz-psimh))
00524         rh = rhn / (1.0_R8+rhn/karman*(alz-psixh))
00525         re = ren / (1.0_R8+ren/karman*(alz-psixh))
00526 
00527         !--- update ustar, tstar, qstar w/ updated, shifted coeffs --
00528         ustar = rd * vmag 
00529         tstar = rh * delt 
00530         qstar = re * delq 
00531 
00532         !----------------------------------------------------------
00533         ! iterate to converge on Z/L, ustar, tstar and qstar
00534         !----------------------------------------------------------
00535 
00536         !--- compute stability & evaluate all stability functions ---
00537         hol    = karman * g * zbot(n) &
00538         &      * (tstar/thvbot+qstar/(1.0_R8/zvir+qbot(n))) / ustar**2
00539         hol    = sign( min(abs(hol),10.0_R8), hol )
00540         stable = 0.5_R8 + sign(0.5_R8 , hol)
00541         xsq    = max(sqrt(abs(1.0_R8 - 16.0_R8*hol)) , 1.0_R8)
00542         xqq    = sqrt(xsq)
00543         psimh  = -5.0_R8*hol*stable + (1.0_R8-stable)*psimhu(xqq)
00544         psixh  = -5.0_R8*hol*stable + (1.0_R8-stable)*psixhu(xqq)
00545 
00546         !--- shift all coeffs to measurement height and stability ---
00547         rd = rdn / (1.0_R8+rdn/karman*(alz-psimh)) 
00548         rh = rhn / (1.0_R8+rhn/karman*(alz-psixh)) 
00549         re = ren / (1.0_R8+ren/karman*(alz-psixh)) 
00550 
00551         !--- update ustar, tstar, qstar w/ updated, shifted coeffs --
00552         ustar = rd * vmag 
00553         tstar = rh * delt 
00554         qstar = re * delq 
00555 
00556         !----------------------------------------------------------
00557         ! compute the fluxes
00558         !----------------------------------------------------------
00559 
00560         tau = rbot(n) * ustar * ustar 
00561     
00562         !--- momentum flux ---
00563         taux(n) = tau * ubot(n) / vmag 
00564         tauy(n) = tau * vbot(n) / vmag 
00565      
00566         !--- heat flux ---
00567         sen (n) =   cp * tau * tstar / ustar 
00568         lat (n) =  ltheat * tau * qstar / ustar
00569         lwup(n) = -stebol * ts(n)**4 
00570      
00571         !--- water flux ---
00572         evap(n) = lat(n)/ltheat 
00573 
00574         !----------------------------------------------------------
00575         ! compute diagnostic: 2m reference height temperature
00576         !----------------------------------------------------------
00577 
00578         !--- Compute function of exchange coefficients. Assume that 
00579         !--- cn = rdn*rdn, cm=rd*rd and ch=rh*rd, and therefore 
00580         !--- 1/sqrt(cn(n))=1/rdn and sqrt(cm(n))/ch(n)=1/rh 
00581         bn = karman/rdn
00582         bh = karman/rh
00583 
00584         !--- Interpolation factor for stable and unstable cases
00585         ln0 = log(1.0_R8 + (ztref/zbot(n))*(exp(bn) - 1.0_R8))
00586         ln3 = log(1.0_R8 + (ztref/zbot(n))*(exp(bn - bh) - 1.0_R8))
00587         fac = (ln0 - ztref/zbot(n)*(bn - bh))/bh * stable &
00588         &   + (ln0 - ln3)/bh * (1.0_R8-stable)
00589         fac = min(max(fac,0.0_R8),1.0_R8)
00590 
00591         !--- actual interpolation
00592         tref(n) = ts(n) + (tbot(n) - ts(n))*fac
00593         qref(n) = qbot(n) - delq*fac
00594 
00595      endif
00596    enddo 
00597 
00598 end subroutine shr_flux_atmIce
00599  
00600 !===============================================================================
00601 ! !BOP =========================================================================
00602 !
00603 ! !IROUTINE: shr_flux_MOstability -- Monin-Obukhov BL stability functions
00604 !
00605 ! !DESCRIPTION:
00606 !
00607 !    Monin-Obukhov boundary layer stability functions, two options:
00608 !    turbulent velocity scales or gradient and integral functions
00609 !    via option = shr_flux_MOwScales or shr_flux_MOfunctions
00610 !     
00611 ! !REVISION HISTORY:
00612 !    2007-Sep-19 - B. Kauffman, Bill Large - first version
00613 !
00614 ! !INTERFACE: ------------------------------------------------------------------
00615 
00616 subroutine shr_flux_MOstability(option,arg1,arg2,arg3,arg4,arg5)
00617 
00618 ! !USES:
00619 
00620    implicit none
00621 
00622 ! !INPUT/OUTPUT PARAMETERS:
00623 
00624    integer(IN),intent(in)           :: option ! shr_flux_MOwScales or MOfunctions 
00625    real(R8)   ,intent(in)           :: arg1   ! scales: uStar (in)  funct: zeta (in)
00626    real(R8)   ,intent(inout)        :: arg2   ! scales: zkB   (in)  funct: phim (out)
00627    real(R8)   ,intent(out)          :: arg3   ! scales: phim  (out) funct: phis (out)
00628    real(R8)   ,intent(out)          :: arg4   ! scales: phis  (out) funct: psim (out)
00629    real(R8)   ,intent(out),optional :: arg5   ! scales:    (unused) funct: psis (out)
00630  
00631 ! !EOP
00632 
00633    !----- local variables -----
00634    real(R8)           :: zeta  ! z/L
00635    real(R8)           :: uStar ! friction velocity
00636    real(R8)           :: zkB   ! (height)*(von Karman)*(surface bouyancy flux)
00637    real(R8)           :: phim  ! momentum    gradient function or scale
00638    real(R8)           :: phis  ! temperature gradient function or scale
00639    real(R8)           :: psim  ! momentum    integral function or scale
00640    real(R8)           :: psis  ! temperature integral function or scale
00641    real(R8)           :: temp  ! temporary-variable/partial calculation
00642 
00643    !----- local variables, stable case -----
00644    real(R8),parameter :: uStarMin = 0.001_R8 ! lower bound on uStar
00645    real(R8),parameter :: a = 1.000_R8  ! constant from Holtslag & de Bruin, equation 12
00646    real(R8),parameter :: b = 0.667_R8  ! constant from Holtslag & de Bruin, equation 12
00647    real(R8),parameter :: c = 5.000_R8  ! constant from Holtslag & de Bruin, equation 12
00648    real(R8),parameter :: d = 0.350_R8  ! constant from Holtslag & de Bruin, equation 12
00649 
00650    !----- local variables, unstable case -----
00651    real(R8),parameter :: a2 = 3.0_R8   ! constant from Wilson, equation 10 
00652 
00653    !----- formats -----
00654    character(*),parameter :: subName = '(shr_flux_MOstability) '
00655    character(*),parameter ::   F00 = "('(shr_flux_MOstability) ',4a)"
00656    character(*),parameter ::   F01 = "('(shr_flux_MOstability) ',a,i5)"
00657 
00658 !-------------------------------------------------------------------------------
00659 ! Notes::
00660 !   o this could be two routines, but are one to help keep them aligned
00661 !   o the stable calculation is taken from...
00662 !     A.A.M. HoltSlag and H.A.R. de Bruin, 1988:
00663 !     "Applied Modeling of the Nighttime Surface Energy Balance over Land",
00664 !     Journal of Applied Meteorology, Vol. 27, No. 6, June 1988, 659-704
00665 !   o the unstable calculation is taken from...
00666 !     D. Keith Wilson, 2001: "An Alternative Function for the Wind and 
00667 !     Temperature Gradients in Unstable Surface Layers", 
00668 !     Boundary-Layer Meteorology, 99 (2001), 151-158
00669 !-------------------------------------------------------------------------------
00670 
00671    !----- check for consistancy between option and arguments ------------------
00672    if (debug > 1 .and. s_loglev > 0) then
00673       if (debug > 2) write(s_logunit,F01) "enter, option = ",option
00674       if ( option == shr_flux_MOwScales .and. present(arg5) ) then
00675          write(s_logunit,F01) "ERROR: option1 must have four arguments"
00676          call shr_sys_abort(subName//"option inconsistant with arguments")
00677       else if ( option == shr_flux_MOfunctions .and. .not. present(arg5) ) then
00678          write(s_logunit,F01) "ERROR: option2 must have five arguments"
00679          call shr_sys_abort(subName//"option inconsistant with arguments")
00680       else 
00681          write(s_logunit,F01) "invalid option = ",option
00682          call shr_sys_abort(subName//"invalid option")
00683       end if
00684    end if
00685 
00686    !------ velocity scales option ----------------------------------------------
00687    if (option == shr_flux_MOwScales) then
00688 
00689       !--- input ---
00690       uStar = arg1
00691       zkB   = arg2
00692 
00693       if (zkB >= 0.0_R8) then ! ----- stable -----
00694          zeta = zkB/(max(uStar,uStarMin)**3)
00695          temp = exp(-d*zeta)
00696          phim = uStar/(1.0_R8 + zeta*(a + b*(1.0_R8 + c - d*zeta)*temp))
00697          phis = phim
00698       else                    ! ----- unstable -----
00699          temp = (zkB*zkB)**(1.0_R8/a2)   ! note: zkB < 0, zkB*zkB > 0
00700          phim = sqrt(uStar**2 + shr_flux_MOgammaM*temp)
00701          phis = sqrt(uStar**2 + shr_flux_MOgammaS*temp)
00702       end if
00703 
00704       !--- output ---
00705       arg3 = phim
00706       arg4 = phis
00707    !  arg5 = <unused>
00708 
00709    !------ stability function option -------------------------------------------
00710    else if (option == shr_flux_MOfunctions) then
00711 
00712       !--- input ---
00713       zeta  = arg1
00714 
00715       if (zeta >= 0.0_R8) then ! ----- stable -----
00716          temp = exp(-d*zeta)
00717          phim =        1.0_R8 + zeta*(a + b*(1.0_R8 + c - d*zeta)*temp)
00718          phis = phim
00719          psim = -a*zeta - b*(zeta - c/d)*temp - b*c/d
00720          psis = psim 
00721       else                    ! ----- unstable ----
00722          temp = (zeta*zeta)**(1.0_R8/a2)   ! note: zeta < 0, zeta*zeta > 0
00723          phim = 1.0_R8/sqrt(1.0_R8 + shr_flux_MOgammaM*temp)
00724          phis = 1.0_R8/sqrt(1.0_R8 + shr_flux_MOgammaS*temp)
00725          psim = a2*log(0.5_R8 + 0.5_R8/phim)
00726          psis = a2*log(0.5_R8 + 0.5_R8/phis)
00727       end if
00728 
00729       !--- output ---
00730       arg2 = phim
00731       arg3 = phis
00732       arg4 = psim
00733       arg5 = psis
00734    !----------------------------------------------------------------------------
00735    else 
00736       write(s_logunit,F01) "invalid option = ",option
00737       call shr_sys_abort(subName//"invalid option")
00738    endif
00739 
00740 end subroutine shr_flux_MOstability
00741 
00742 !===============================================================================
00743 !===============================================================================
00744 
00745 end module shr_flux_mod