SURFEX v7.3
General documentation of Surfex
 All Classes Files Functions Variables Typedefs
/home/dasprezs/EXPORT_v7_3/src/SURFEX/mr98.F90
Go to the documentation of this file.
00001 !     #########
00002       SUBROUTINE MR98      (PZ0SEA,                                         &
00003                               PTA, PEXNA, PRHOA, PSST, PEXNS, PQA,            &
00004                               PTT,                                            &
00005                               PVMOD, PZREF, PUREF,                            &
00006                               PPS, PQSAT,                                     &
00007                               PSFTH, PSFTQ, PUSTAR,                           &
00008                               PCD, PCDN, PCH, PRI, PRESA, PZ0HSEA             )  
00009 !     #######################################################################
00010 !
00011 !
00012 !!****  *MR98*  
00013 !!
00014 !!    PURPOSE
00015 !!    -------
00016 !      Calculate the surface fluxes of heat, moisture, and momentum over
00017 !      water surfaces.  
00018 !     
00019 !!**  METHOD
00020 !!    ------
00021 !
00022 !!    EXTERNAL
00023 !!    --------
00024 !!
00025 !!    IMPLICIT ARGUMENTS
00026 !!    ------------------ 
00027 !!
00028 !!      
00029 !!    REFERENCE
00030 !!    ---------
00031 !!      
00032 !!    AUTHOR
00033 !!    ------
00034 !!      V. Masson           * Meteo-France *
00035 !!      from part of water_flux.f90 written by M. tomasini
00036 !!
00037 !!    MODIFICATIONS
00038 !!    -------------
00039 !!      Original      11/2004
00040 !!      (M.Tomasini)  25/07/00  Add an iterative method for the fluxes calculation
00041 !!-------------------------------------------------------------------------------
00042 !
00043 !*       0.     DECLARATIONS
00044 !               ------------
00045 !
00046 USE MODD_CSTS,     ONLY : XPI, XCPD, XG, XKARMAN
00047 USE MODD_SURF_PAR, ONLY : XUNDEF
00048 !
00049 USE MODE_THERMOS
00050 !
00051 !
00052 !
00053 USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
00054 USE PARKIND1  ,ONLY : JPRB
00055 !
00056 IMPLICIT NONE
00057 !
00058 !*      0.1    declarations of arguments
00059 !
00060 !
00061 REAL, DIMENSION(:), INTENT(IN)       :: PTA   ! air temperature at atm. level
00062 REAL, DIMENSION(:), INTENT(IN)       :: PQA   ! air humidity at atm. level (kg/kg)
00063 REAL, DIMENSION(:), INTENT(IN)       :: PEXNA ! Exner function at atm. level
00064 REAL, DIMENSION(:), INTENT(IN)       :: PRHOA ! air density at atm. level
00065 REAL, DIMENSION(:), INTENT(IN)       :: PVMOD ! module of wind at atm. wind level
00066 REAL, DIMENSION(:), INTENT(IN)       :: PZREF ! atm. level for temp. and humidity
00067 REAL, DIMENSION(:), INTENT(IN)       :: PUREF ! atm. level for wind
00068 REAL, DIMENSION(:), INTENT(IN)       :: PSST  ! Sea Surface Temperature
00069 REAL, DIMENSION(:), INTENT(IN)       :: PEXNS ! Exner function at sea surface
00070 REAL, DIMENSION(:), INTENT(IN)       :: PPS   ! air pressure at sea surface
00071 REAL,               INTENT(IN)       :: PTT   ! temperature of freezing point
00072 !
00073 REAL, DIMENSION(:), INTENT(INOUT)    :: PZ0SEA! roughness length over the ocean
00074 !                                         
00075 !                                         
00076 !  surface fluxes : latent heat, sensible heat, friction fluxes
00077 REAL, DIMENSION(:), INTENT(OUT)      :: PSFTH ! heat flux  (W/m2)
00078 REAL, DIMENSION(:), INTENT(OUT)      :: PSFTQ ! water flux (kg/m2/s)
00079 REAL, DIMENSION(:), INTENT(OUT)      :: PUSTAR! friction velocity (m/s)
00080 !
00081 ! diagnostics
00082 REAL, DIMENSION(:), INTENT(OUT)      :: PQSAT ! humidity at saturation
00083 REAL, DIMENSION(:), INTENT(OUT)      :: PCD   ! heat drag coefficient
00084 REAL, DIMENSION(:), INTENT(OUT)      :: PCDN  ! momentum drag coefficient
00085 REAL, DIMENSION(:), INTENT(OUT)      :: PCH   ! neutral momentum drag coefficient
00086 REAL, DIMENSION(:), INTENT(OUT)      :: PRI   ! Richardson number
00087 REAL, DIMENSION(:), INTENT(OUT)      :: PRESA     ! aerodynamical resistance
00088 REAL, DIMENSION(:), INTENT(OUT)      :: PZ0HSEA   ! roughness length for heat
00089 !                                         
00090 !
00091 !
00092 !*      0.2    declarations of local variables
00093 !
00094 !
00095 REAL, DIMENSION(SIZE(PTA)) :: ZZ0SEAH,        !roughness length heat
00096                                 ZZ0SEAQ,        !roughness length moisture
00097                                 ZTHVI  
00098 !
00099 REAL, DIMENSION(SIZE(PTA)) :: ZQSA, ZTHA, ZDU, ZDT, ZDQ, ZNU, ZWG
00100 !                             ZQSA = specific humidity at the
00101 !                                     sea surface
00102 !                             ZTHA = potential temperature at the
00103 !                                     the first level
00104 !                             ZDU, ZDT, ZDQ = wind, potential
00105 !                                     temperature and specific 
00106 !                                     humidity difference between
00107 !                                     the first level and the sea
00108 !                             ZNU = air kinematic visosity
00109 !                             ZWG = wind subgrid "gustiness" effects
00110 REAL, DIMENSION(SIZE(PTA)) :: ZTSTAR, ZQSTAR, ZTVSTAR,          
00111 !                             T*, Q*, Tv* usefull in fluxes computation 
00112                                 ZBUFLX  
00113 !                             ZBUFLX = surface buoyancy flux
00114 !
00115 REAL, DIMENSION(SIZE(PTA)) :: ZLMON, ZZETA, ZPSIM, ZPSIH, ZPSIQ,  
00116                                 ZKHI, ZKHIC, ZPSIC, ZF,             
00117                                 ZPSIHSTAND,ZPSIMSTAND,ZZETAST  
00118 !                             ZLMON = Monin Obukhov length
00119 !                             ZZETA = Monin Obukhov Stability
00120 !                                     parameter : Zeta=Zref/Lmon 
00121 !                             ZPSIM, ZPSIH, ZPSIQ = Stability
00122 !                                     functions for Momentum,
00123 !                                     Heat and Moisture
00124 !                             ZKHI = variable for Stability functions
00125 !                             ZPSIC = Stability function correction for
00126 !                                     very unstable conditions
00127 !                             ZKHIC = variable for ZPSIC
00128 !                             ZF = weight parameter in psi function
00129 !                                     computation
00130 !
00131 !  constants for wind gustiness effect (ZWG) calculation
00132 !                                            
00133 REAL                      ::  ZBETA
00134 !                             ZBETA = coef.  for the gustiness
00135 !                                      effect on the wind 
00136 
00137 !                                     Wg=(ZBETA).W*
00138 REAL                      ::  ZZBL
00139 !                             ZZBL = atmospheric boundary 
00140 !                                    layer depth (m)
00141 !
00142 REAL                      ::  ZSTAND        
00143 !
00144 !  loop parameters for iterative case
00145 !
00146 INTEGER                   ::  JITER
00147 !                             JITER = loop parameter for
00148 !                                     iterations
00149 !
00150 INTEGER, PARAMETER        ::  IITERMAX = 10
00151 REAL(KIND=JPRB) :: ZHOOK_HANDLE
00152 !                             IITERMAX = number of iterations
00153 !
00154 !-------------------------------------------------------------------------------
00155 !
00156 IF (LHOOK) CALL DR_HOOK('MR98',0,ZHOOK_HANDLE)
00157 PRI(:) = XUNDEF
00158 PCH(:) = XUNDEF
00159 PCD(:) = XUNDEF
00160 PCDN(:) = XUNDEF
00161 PRESA(:) = XUNDEF
00162 !
00163 PSFTH(:) =XUNDEF
00164 PSFTQ(:)=XUNDEF
00165 PUSTAR(:)=XUNDEF
00166 !
00167 !                                    Although it is not necessary, fluxes over
00168 !                                    water surfaces are calculated at every
00169 !                                    grid-points, even for those over the
00170 !                                    continent.  When the "WHERE" function of
00171 !                                    the CRAY F90 will be vectorized, it can
00172 !                                    be used to discriminate points over the
00173 !                                    ocean and the continent.  But for now, 
00174 !                                    the current approach is believed the
00175 !                                    simplest, without being too penalizing.
00176 !                                    It mainly consists in averaging the
00177 !                                    surface fluxes over water and land using
00178 !                                    the cover type variables.
00179 !
00180 !
00181 !       1.     Saturated specific humidity near the water surface
00182 !       ---------------------------------------------------------
00183 !
00184 PQSAT(:) = QSAT(PSST(:),PPS(:))
00185 !
00186 !-------------------------------------------------------------------------------
00187 !
00188 !         Begining of the fluxes calculation
00189 !            by an iterative method coming from the one
00190 !                                  proposed by Fairall et al (1996) for
00191 !                                  the TOGA-COARE Bulk Flux Algorithm
00192 !
00193 !
00194 !
00195 !       3.1            initialisations
00196 !       ----------------------------------------------------
00197 !
00198 ZZBL=650.       
00199 ZSTAND=15.         
00200 ZBETA = 0.6
00201 !
00202   ZQSA(:)  = 0.98 * PQSAT(:)
00203   ZTHA(:)  = PTA(:) / PEXNA(:)
00204   ZTHVI(:) = ( 1 + 0.61 * PQA(:) ) * ZTHA(:)
00205 !
00206   ZDU(:) = PVMOD(:)
00207   ZDT(:) = ZTHA(:) - PSST(:)  /PEXNS(:)
00208   ZDQ(:) = PQA(:)  - ZQSA(:)
00209 !
00210   PUSTAR(:) = MAX( 0.04 * ZDU(:) , 5E-3)
00211   ZTSTAR(:)  = 0.04 * ZDT(:)
00212   ZQSTAR(:)  = 0.04 * ZDQ(:)
00213   ZTVSTAR(:) = ZTSTAR(:)*(1+0.61*PQA(:)) + 0.61*ZQSTAR(:)*ZTHA(:)
00214 !
00215   ZNU(:) =  1.318E-5 + 9.282E-8 * (PTA(:) - PTT)
00216 !
00217 !
00218 !
00219 !           3.2    begining of iterations 
00220 !           -------------------------------------------------
00221 !
00222   DO JITER = 1 , IITERMAX
00223 !
00224 !           3.2.1   compute stability functions
00225 !           -------------------------------------------------
00226 !
00227     WHERE (ABS(ZTVSTAR(:)) < 1E-6)
00228       ZZETA(:)   = 0.0
00229       ZZETAST(:) = 0.0
00230     ELSEWHERE
00231       ZLMON(:)   =  ZTHVI(:)*PUSTAR(:)*PUSTAR(:)/(ZTVSTAR(:)*XG*XKARMAN)
00232       ZZETA(:)   =  MAX( PZREF(:) / ZLMON(:) , -20000.)
00233       ZZETAST(:) =  MAX( ZSTAND / ZLMON(:) , -20000.)
00234     END WHERE
00235 !
00236 !
00237     WHERE(ZZETA(:) >= 0.0)
00238       ZPSIM(:) = -4.7*ZZETA(:)
00239       ZPSIH(:) = ZPSIM(:)
00240       ZPSIQ(:) = ZPSIH(:)
00241       ZPSIHSTAND(:) = -4.7*ZZETAST(:)
00242       ZPSIMSTAND(:) = ZPSIHSTAND(:)
00243     ELSEWHERE
00244       ZKHI(:) = (1 - 16 * ZZETA(:))**0.25
00245       ZPSIM(:) =  2*LOG((1+ZKHI(:))/2)                            &
00246                     + LOG((1+ZKHI(:)*ZKHI(:))/2)                    &
00247                     - 2*ATAN(ZKHI(:))                               &
00248                     + XPI/2  
00249       ZPSIH(:) = 2 * LOG((1+ZKHI(:)*ZKHI(:))/2) 
00250 !
00251       ZKHI(:)  = (1 - 16 * ZZETAST(:))**0.25
00252       ZPSIHSTAND(:) =  2 * LOG((1+ZKHI(:)*ZKHI(:))/2) 
00253       ZPSIMSTAND(:) =   2*LOG((1+ZKHI(:))/2)                                   &
00254                         + LOG((1+ZKHI(:)*ZKHI(:))/2)                             &
00255                         - 2*ATAN(ZKHI(:))                                        &
00256                         + XPI/2  
00257 !
00258 ! to match very unstable conditions
00259 !
00260       ZKHIC(:) = (1 - 12.87 * ZZETA(:))**0.33
00261       ZPSIC(:) =  1.5 * LOG ((1+ZKHIC(:)+ZKHIC(:)*ZKHIC(:))/3)                 &
00262                     - (3**0.5)*ATAN((2*ZKHIC(:)+1)/(3**0.5))                     &
00263                     + XPI/(3**0.5)  
00264 !
00265       ZF(:) = 1 / (1+ ZZETA(:)*ZZETA(:)) 
00266       ZPSIM(:) =  ZPSIM(:)*ZF(:) + ZPSIC(:)*(1-ZF(:))
00267       ZPSIH(:) =  ZPSIH(:)*ZF(:) + ZPSIC(:)*(1-ZF(:))
00268       ZPSIQ(:) =  ZPSIH(:) 
00269 !
00270       ZKHIC(:) = (1 - 12.87 * ZZETAST(:))**0.33
00271       ZPSIC(:) =  1.5 * LOG ((1+ZKHIC(:)+ZKHIC(:)*ZKHIC(:))/3)                 &
00272                     - (3**0.5)*ATAN((2*ZKHIC(:)+1)/(3**0.5))                     &
00273                     + XPI/(3**0.5)  
00274 !
00275       ZF(:) = 1 / (1+ ZZETAST(:)*ZZETAST(:)) 
00276       ZPSIMSTAND(:) =  ZPSIMSTAND(:)*ZF(:) + ZPSIC(:)*(1-ZF(:))
00277       ZPSIHSTAND(:) =  ZPSIHSTAND(:)*ZF(:) + ZPSIC(:)*(1-ZF(:))
00278 !
00279     END WHERE
00280 !
00281 !           3.2.2 compute roughness length Zo, Zoh and Zoq
00282 !           ----------------------------------------------
00283 !
00284     PZ0SEA(:) = 0.011*PUSTAR(:)*PUSTAR(:)/XG                                 &
00285                    +0.11*ZNU(:)/(PUSTAR(:))  
00286 !
00287 !
00288     WHERE(PUSTAR(:) > 0.23)
00289       ZZ0SEAH(:) = 0.14*ZNU(:)/(PUSTAR(:)-0.2) + 7E-6
00290       ZZ0SEAQ(:) = 0.20*ZNU(:)/(PUSTAR(:)-0.2) + 9E-6
00291     ELSEWHERE
00292       ZZ0SEAH(:) = 0.015*PUSTAR(:)*PUSTAR(:)/XG                             &
00293                       + 0.18*ZNU(:)/(PUSTAR(:))  
00294       ZZ0SEAQ(:) = 0.0205*PUSTAR(:)*PUSTAR(:)/XG                            &
00295                       + 0.294*ZNU(:)/(PUSTAR(:))  
00296     END WHERE
00297 !
00298 !
00299 !          3.2.3  compute friction velocity u*, T*, Q* and Tv*
00300 !          ---------------------------------------------------
00301 !
00302     PUSTAR(:) = MAX(                                                         &
00303                     ZDU(:)*XKARMAN / (LOG(PZREF(:)/PZ0SEA(:)) - ZPSIM(:))       &
00304                     , 5E-3 )  
00305 !
00306     ZTSTAR(:) = ZDT(:)*XKARMAN / (LOG(PZREF(:)/ZZ0SEAH(:)) - ZPSIH(:))
00307 !
00308     ZQSTAR(:) = ZDQ(:)*XKARMAN / (LOG(PZREF(:)/ZZ0SEAQ(:)) - ZPSIQ(:))
00309 !
00310     ZTVSTAR(:) = ZTSTAR(:)*(1+0.61*PQA(:))+0.61*ZQSTAR(:)*ZTHA(:)
00311 !
00312 !
00313 !         3.2.4  compute subgrid correction for wind due to convective motions
00314 !         --------------------------------------------------------------------
00315 !
00316     ZBUFLX(:) = -XG*ZTVSTAR(:)*PUSTAR(:)/ZTHVI(:)
00317 !
00318     WHERE(ZBUFLX(:) > 0.)
00319       ZWG(:) = ZBETA*(ZBUFLX(:)*ZZBL)**0.333
00320     ELSEWHERE
00321       ZWG(:) = 0.
00322     END WHERE
00323 !
00324     ZDU(:) = (PVMOD(:)*PVMOD(:)+ZWG(:)*ZWG(:))**0.5
00325 !
00326   END DO
00327 !
00328 PZ0HSEA = ZZ0SEAH
00329 !-------------------------------------------------------------------------------
00330 !
00331 !        3.3  here are the fluxes of iterative case
00332 !        ------------------------------------------
00333 !
00334   PSFTH(:)  = - PRHOA(:)  * XCPD  * PUSTAR(:) * ZTSTAR(:)
00335   PSFTQ(:)  = - PRHOA(:)  *         PUSTAR(:) * ZQSTAR(:)
00336 IF (LHOOK) CALL DR_HOOK('MR98',1,ZHOOK_HANDLE)
00337 !
00338 !-------------------------------------------------------------------------------
00339 !
00340 END SUBROUTINE MR98