SURFEX v7.3
General documentation of Surfex
|
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