mr98.F90

Go to the documentation of this file.

!SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
!SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
!SFX_LIC for details. version 1.
!     #########
      SUBROUTINE mr98      (PZ0SEA,                                         &
                              pta, pexna, prhoa, psst, pexns, pqa,            &
                              ptt,                                            &
                              pvmod, pzref, puref,                            &
                              pps, pqsat,                                     &
                              psfth, psftq, pustar,                           &
                              pcd, pcdn, pch, pri, presa, pz0hsea             )  
!     #######################################################################
!
!
!!****  *MR98*  
!!
!!    PURPOSE
!!    -------
!      Calculate the surface fluxes of heat, moisture, and momentum over
!      water surfaces.  
!     
!!**  METHOD
!!    ------
!
!!    EXTERNAL
!!    --------
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------ 
!!
!!      
!!    REFERENCE
!!    ---------
!!      
!!    AUTHOR
!!    ------
!!      V. Masson           * Meteo-France *
!!      from part of water_flux.f90 written by M. tomasini
!!
!!    MODIFICATIONS
!!    -------------
!!      Original      11/2004
!!      (M.Tomasini)  25/07/00  Add an iterative method for the fluxes calculation
!!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
!
USE modd_csts,     ONLY : xpi, xcpd, xg, xkarman
USE modd_surf_par, ONLY : xundef
!
USE mode_thermos
!
!
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
!
REAL, DIMENSION(:), INTENT(IN)       :: pta   ! air temperature at atm. level
REAL, DIMENSION(:), INTENT(IN)       :: pqa   ! air humidity at atm. level (kg/kg)
REAL, DIMENSION(:), INTENT(IN)       :: pexna ! Exner function at atm. level
REAL, DIMENSION(:), INTENT(IN)       :: prhoa ! air density at atm. level
REAL, DIMENSION(:), INTENT(IN)       :: pvmod ! module of wind at atm. wind level
REAL, DIMENSION(:), INTENT(IN)       :: pzref ! atm. level for temp. and humidity
REAL, DIMENSION(:), INTENT(IN)       :: puref ! atm. level for wind
REAL, DIMENSION(:), INTENT(IN)       :: psst  ! Sea Surface Temperature
REAL, DIMENSION(:), INTENT(IN)       :: pexns ! Exner function at sea surface
REAL, DIMENSION(:), INTENT(IN)       :: pps   ! air pressure at sea surface
REAL,               INTENT(IN)       :: ptt   ! temperature of freezing point
!
REAL, DIMENSION(:), INTENT(INOUT)    :: pz0sea! roughness length over the ocean
!                                         
!                                         
!  surface fluxes : latent heat, sensible heat, friction fluxes
REAL, DIMENSION(:), INTENT(OUT)      :: psfth ! heat flux  (W/m2)
REAL, DIMENSION(:), INTENT(OUT)      :: psftq ! water flux (kg/m2/s)
REAL, DIMENSION(:), INTENT(OUT)      :: pustar! friction velocity (m/s)
!
! diagnostics
REAL, DIMENSION(:), INTENT(OUT)      :: pqsat ! humidity at saturation
REAL, DIMENSION(:), INTENT(OUT)      :: pcd   ! heat drag coefficient
REAL, DIMENSION(:), INTENT(OUT)      :: pcdn  ! momentum drag coefficient
REAL, DIMENSION(:), INTENT(OUT)      :: pch   ! neutral momentum drag coefficient
REAL, DIMENSION(:), INTENT(OUT)      :: pri   ! Richardson number
REAL, DIMENSION(:), INTENT(OUT)      :: presa     ! aerodynamical resistance
REAL, DIMENSION(:), INTENT(OUT)      :: pz0hsea   ! roughness length for heat
!                                         
!
!
!*      0.2    declarations of local variables
!
!
REAL, DIMENSION(SIZE(PTA)) :: zz0seah,        &!roughness length heat
                                zz0seaq,        &!roughness length moisture
                                zthvi  
!
REAL, DIMENSION(SIZE(PTA)) :: zqsa, ztha, zdu, zdt, zdq, znu, zwg
!                             ZQSA = specific humidity at the
!                                     sea surface
!                             ZTHA = potential temperature at the
!                                     the first level
!                             ZDU, ZDT, ZDQ = wind, potential
!                                     temperature and specific 
!                                     humidity difference between
!                                     the first level and the sea
!                             ZNU = air kinematic visosity
!                             ZWG = wind subgrid "gustiness" effects
REAL, DIMENSION(SIZE(PTA)) :: ztstar, zqstar, ztvstar,          &
!                             T*, Q*, Tv* usefull in fluxes computation 
                                zbuflx  
!                             ZBUFLX = surface buoyancy flux
!
REAL, DIMENSION(SIZE(PTA)) :: zlmon, zzeta, zpsim, zpsih, zpsiq,  &
                                zkhi, zkhic, zpsic, zf,             &
                                zpsihstand,zpsimstand,zzetast  
!                             ZLMON = Monin Obukhov length
!                             ZZETA = Monin Obukhov Stability
!                                     parameter : Zeta=Zref/Lmon 
!                             ZPSIM, ZPSIH, ZPSIQ = Stability
!                                     functions for Momentum,
!                                     Heat and Moisture
!                             ZKHI = variable for Stability functions
!                             ZPSIC = Stability function correction for
!                                     very unstable conditions
!                             ZKHIC = variable for ZPSIC
!                             ZF = weight parameter in psi function
!                                     computation
!
!  constants for wind gustiness effect (ZWG) calculation
!                                            
REAL                      ::  zbeta
!                             ZBETA = coef.  for the gustiness
!                                      effect on the wind 

!                                     Wg=(ZBETA).W*
REAL                      ::  zzbl
!                             ZZBL = atmospheric boundary 
!                                    layer depth (m)
!
REAL                      ::  zstand        
!
!  loop parameters for iterative case
!
INTEGER                   ::  jiter
!                             JITER = loop parameter for
!                                     iterations
!
INTEGER, PARAMETER        ::  iitermax = 10
REAL(KIND=JPRB) :: zhook_handle
!                             IITERMAX = number of iterations
!
!-------------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('MR98',0,zhook_handle)
pri(:) = xundef
pch(:) = xundef
pcd(:) = xundef
pcdn(:) = xundef
presa(:) = xundef
!
psfth(:) =xundef
psftq(:)=xundef
pustar(:)=xundef
!
!                                    Although it is not necessary, fluxes over
!                                    water surfaces are calculated at every
!                                    grid-points, even for those over the
!                                    continent.  When the "WHERE" function of
!                                    the CRAY F90 will be vectorized, it can
!                                    be used to discriminate points over the
!                                    ocean and the continent.  But for now, 
!                                    the current approach is believed the
!                                    simplest, without being too penalizing.
!                                    It mainly consists in averaging the
!                                    surface fluxes over water and land using
!                                    the cover type variables.
!
!
!       1.     Saturated specific humidity near the water surface
!       ---------------------------------------------------------
!
pqsat(:) = qsat(psst(:),pps(:))
!
!-------------------------------------------------------------------------------
!
!         Begining of the fluxes calculation
!            by an iterative method coming from the one
!                                  proposed by Fairall et al (1996) for
!                                  the TOGA-COARE Bulk Flux Algorithm
!
!
!
!       3.1            initialisations
!       ----------------------------------------------------
!
zzbl=650.       
zstand=15.         
zbeta = 0.6
!
  zqsa(:)  = 0.98 * pqsat(:)
  ztha(:)  = pta(:) / pexna(:)
  zthvi(:) = ( 1 + 0.61 * pqa(:) ) * ztha(:)
!
  zdu(:) = pvmod(:)
  zdt(:) = ztha(:) - psst(:)  /pexns(:)
  zdq(:) = pqa(:)  - zqsa(:)
!
  pustar(:) = max( 0.04 * zdu(:) , 5e-3)
  ztstar(:)  = 0.04 * zdt(:)
  zqstar(:)  = 0.04 * zdq(:)
  ztvstar(:) = ztstar(:)*(1+0.61*pqa(:)) + 0.61*zqstar(:)*ztha(:)
!
  znu(:) =  1.318e-5 + 9.282e-8 * (pta(:) - ptt)
!
!
!
!           3.2    begining of iterations 
!           -------------------------------------------------
!
  DO jiter = 1 , iitermax
!
!           3.2.1   compute stability functions
!           -------------------------------------------------
!
    WHERE (abs(ztvstar(:)) < 1e-6)
      zzeta(:)   = 0.0
      zzetast(:) = 0.0
    ELSEWHERE
      zlmon(:)   =  zthvi(:)*pustar(:)*pustar(:)/(ztvstar(:)*xg*xkarman)
      zzeta(:)   =  max( pzref(:) / zlmon(:) , -20000.)
      zzetast(:) =  max( zstand / zlmon(:) , -20000.)
    END WHERE
!
!
    WHERE(zzeta(:) >= 0.0)
      zpsim(:) = -4.7*zzeta(:)
      zpsih(:) = zpsim(:)
      zpsiq(:) = zpsih(:)
      zpsihstand(:) = -4.7*zzetast(:)
      zpsimstand(:) = zpsihstand(:)
    ELSEWHERE
      zkhi(:) = (1 - 16 * zzeta(:))**0.25
      zpsim(:) =  2*log((1+zkhi(:))/2)                            &
                    + log((1+zkhi(:)*zkhi(:))/2)                    &
                    - 2*atan(zkhi(:))                               &
                    + xpi/2  
      zpsih(:) = 2 * log((1+zkhi(:)*zkhi(:))/2) 
!
      zkhi(:)  = (1 - 16 * zzetast(:))**0.25
      zpsihstand(:) =  2 * log((1+zkhi(:)*zkhi(:))/2) 
      zpsimstand(:) =   2*log((1+zkhi(:))/2)                                   &
                        + log((1+zkhi(:)*zkhi(:))/2)                             &
                        - 2*atan(zkhi(:))                                        &
                        + xpi/2  
!
! to match very unstable conditions
!
      zkhic(:) = (1 - 12.87 * zzeta(:))**0.33
      zpsic(:) =  1.5 * log((1+zkhic(:)+zkhic(:)*zkhic(:))/3)                 &
                    - (3**0.5)*atan((2*zkhic(:)+1)/(3**0.5))                     &
                    + xpi/(3**0.5)  
!
      zf(:) = 1 / (1+ zzeta(:)*zzeta(:)) 
      zpsim(:) =  zpsim(:)*zf(:) + zpsic(:)*(1-zf(:))
      zpsih(:) =  zpsih(:)*zf(:) + zpsic(:)*(1-zf(:))
      zpsiq(:) =  zpsih(:) 
!
      zkhic(:) = (1 - 12.87 * zzetast(:))**0.33
      zpsic(:) =  1.5 * log((1+zkhic(:)+zkhic(:)*zkhic(:))/3)                 &
                    - (3**0.5)*atan((2*zkhic(:)+1)/(3**0.5))                     &
                    + xpi/(3**0.5)  
!
      zf(:) = 1 / (1+ zzetast(:)*zzetast(:)) 
      zpsimstand(:) =  zpsimstand(:)*zf(:) + zpsic(:)*(1-zf(:))
      zpsihstand(:) =  zpsihstand(:)*zf(:) + zpsic(:)*(1-zf(:))
!
    END WHERE
!
!           3.2.2 compute roughness length Zo, Zoh and Zoq
!           ----------------------------------------------
!
    pz0sea(:) = 0.011*pustar(:)*pustar(:)/xg                                 &
                   +0.11*znu(:)/(pustar(:))  
!
!
    WHERE(pustar(:) > 0.23)
      zz0seah(:) = 0.14*znu(:)/(pustar(:)-0.2) + 7e-6
      zz0seaq(:) = 0.20*znu(:)/(pustar(:)-0.2) + 9e-6
    ELSEWHERE
      zz0seah(:) = 0.015*pustar(:)*pustar(:)/xg                             &
                      + 0.18*znu(:)/(pustar(:))  
      zz0seaq(:) = 0.0205*pustar(:)*pustar(:)/xg                            &
                      + 0.294*znu(:)/(pustar(:))  
    END WHERE
!
!
!          3.2.3  compute friction velocity u*, T*, Q* and Tv*
!          ---------------------------------------------------
!
    pustar(:) = max(                                                         &
                    zdu(:)*xkarman / (log(pzref(:)/pz0sea(:)) - zpsim(:))       &
                    , 5e-3 )  
!
    ztstar(:) = zdt(:)*xkarman / (log(pzref(:)/zz0seah(:)) - zpsih(:))
!
    zqstar(:) = zdq(:)*xkarman / (log(pzref(:)/zz0seaq(:)) - zpsiq(:))
!
    ztvstar(:) = ztstar(:)*(1+0.61*pqa(:))+0.61*zqstar(:)*ztha(:)
!
!
!         3.2.4  compute subgrid correction for wind due to convective motions
!         --------------------------------------------------------------------
!
    zbuflx(:) = -xg*ztvstar(:)*pustar(:)/zthvi(:)
!
    WHERE(zbuflx(:) > 0.)
      zwg(:) = zbeta*(zbuflx(:)*zzbl)**0.333
    ELSEWHERE
      zwg(:) = 0.
    END WHERE
!
    zdu(:) = (pvmod(:)*pvmod(:)+zwg(:)*zwg(:))**0.5
!
  END DO
!
pz0hsea = zz0seah
!-------------------------------------------------------------------------------
!
!        3.3  here are the fluxes of iterative case
!        ------------------------------------------
!
  psfth(:)  = - prhoa(:)  * xcpd  * pustar(:) * ztstar(:)
  psftq(:)  = - prhoa(:)  *         pustar(:) * zqstar(:)
IF (lhook) CALL dr_hook('MR98',1,zhook_handle)
!
!-------------------------------------------------------------------------------
!
END SUBROUTINE mr98