Changeset 2143
- Timestamp:
- Jun 20, 2019, 4:11:54 PM (6 years ago)
- Location:
- trunk/LMDZ.GENERIC
- Files:
-
- 1 added
- 9 edited
-
README (modified) (1 diff)
-
libf/phystd/thermcell_closure.F90 (modified) (9 diffs)
-
libf/phystd/thermcell_dq.F90 (modified) (4 diffs)
-
libf/phystd/thermcell_dv2.F90 (added)
-
libf/phystd/thermcell_env.F90 (modified) (4 diffs)
-
libf/phystd/thermcell_flux.F90 (modified) (17 diffs)
-
libf/phystd/thermcell_height.F90 (modified) (15 diffs)
-
libf/phystd/thermcell_main.F90 (modified) (12 diffs)
-
libf/phystd/thermcell_mod.F90 (modified) (5 diffs)
-
libf/phystd/thermcell_plume.F90 (modified) (13 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.GENERIC/README
r2138 r2143 1483 1483 but would require to be able to have writediag in 5D) 1484 1484 EDIT (22/05/19) : To have correct calculations, output is now exp(-tau), you need to postproc it to have tau. 1485 1486 == 20/06/2019 == AB 1487 - update the thermal plume model (check formulae consistency between thermcell_plume and thercell_closure, 1488 compute correctly thermal plume height, fix alimentation computation in the first unstable layer) 1489 -
trunk/LMDZ.GENERIC/libf/phystd/thermcell_closure.F90
r2127 r2143 3 3 ! 4 4 SUBROUTINE thermcell_closure(ngrid,nlay,ptimestep,rho,zlev, & 5 lmax,alim_star,zmax,wmax,f) 5 lmax,alim_star,zmin,zmax,wmax,f) 6 6 7 7 8 !=============================================================================== … … 13 14 ! 2. On ne garde qu'une version des couples wmax,zmax et wmax_sec,zmax_sec 14 15 ! l'idee etant que le choix se fasse a l'appel de thermcell_closure 15 ! 3. Vectorisation en mettant les boucles en l l'exterieur avec des if16 ! 3. Vectorisation en mettant les boucles en l a l'exterieur avec des if 16 17 !=============================================================================== 17 18 … … 35 36 REAL zlev(ngrid,nlay) 36 37 REAL alim_star(ngrid,nlay) 38 REAL zmin(ngrid) 37 39 REAL zmax(ngrid) 38 40 REAL wmax(ngrid) … … 51 53 REAL alim_star_tot(ngrid) 52 54 REAL alim_star2(ngrid) 53 REAL zdenom(ngrid)55 REAL plume_height(ngrid) 54 56 55 57 !=============================================================================== … … 64 66 llmax = 1 65 67 68 DO ig=1,ngrid 69 plume_height(ig) = zmax(ig) - zmin(ig) 70 ENDDO 71 66 72 !=============================================================================== 67 73 ! Closure … … 73 79 74 80 DO ig=1,ngrid 75 IF (lmax(ig) .GT.llmax) THEN81 IF (lmax(ig) > llmax) THEN 76 82 llmax = lmax(ig) 77 83 ENDIF … … 86 92 IF (l < lmax(ig)) THEN 87 93 alim_star2(ig) = alim_star2(ig) + alim_star(ig,l)**2 & 88 & / (rho(ig,l) * (zlev(ig,l+1) - zlev(ig,l))) 94 & / (rho(ig,l) * (zlev(ig,l+1) - zlev(ig,l))) ! => intergration is ok because alim_star = a* dz 89 95 alim_star_tot(ig) = alim_star_tot(ig) + alim_star(ig,l) 90 96 ENDIF … … 93 99 94 100 DO ig=1,ngrid 95 IF ( alim_star2(ig).GT.0.) THEN96 f(ig) = wmax(ig) * alim_star_tot(ig) & 97 & / ( max(1.,zmax(ig)) * r_aspect_thermals * alim_star2(ig))101 IF ((alim_star2(ig) > 0.).and.(plume_height(ig) > 0.)) THEN 102 f(ig) = wmax(ig) * alim_star_tot(ig) & ! => normalization is ok 103 & / (plume_height(ig) * r_aspect_thermals * alim_star2(ig)) 98 104 ELSE 99 105 f(ig) = 0. … … 101 107 ENDDO 102 108 109 print *, 'A*2 ', alim_star2(1) 110 print *, 'A* ', alim_star_tot(1) 111 print *, 'H ', plume_height(1) 112 print *, 'wmax', wmax(1) 103 113 104 114 RETURN -
trunk/LMDZ.GENERIC/libf/phystd/thermcell_dq.F90
r2127 r2143 53 53 ! ------ 54 54 55 INTEGER ig, l 55 INTEGER ig, l, k 56 56 INTEGER niter, iter 57 57 … … 59 59 REAL qold(ngrid,nlay) 60 60 REAL fqa(ngrid,nlay+1) 61 REAL wqd(ngrid,nlay+1)62 61 REAL zzm 63 62 … … 86 85 print *, 'ERROR: entrainment is greater than the layer mass!' 87 86 print *, 'ig,l,entr', ig, l, entr(ig,l) 87 print *, 'lmin', lmin(ig) 88 88 print *, '-------------------------------' 89 89 print *, 'entr*dt,mass', entr(ig,l)*ptimestep, masse(ig,l) 90 90 print *, '-------------------------------' 91 print *, 'fm+', fm(ig,l+1) 92 print *, 'entr,detr', entr(ig,l), detr(ig,l) 93 print *, 'fm ', fm(ig,l) 94 print *, 'entr,detr', entr(ig,l-1), detr(ig,l-1) 95 print *, 'fm-', fm(ig,l-1) 91 DO k=nlay,1,-1 92 print *, 'fm ', fm(ig,k+1) 93 print *, 'entr,detr', entr(ig,k), detr(ig,k) 94 ENDDO 95 print *, 'fm ', fm(ig,1) 96 print *, '-------------------------------' 96 97 CALL abort 97 98 ENDIF … … 117 118 qa(ig,l) = q(ig,l) 118 119 ENDIF 119 120 ! IF (qa(ig,l).lt.0.) THEN121 ! print *, 'WARNING: qa is negative!'122 ! print *, 'qa', qa(ig,l)123 ! ENDIF124 125 ! IF (q(ig,l).lt.0.) THEN126 ! print *, 'WARNING: q is negative!'127 ! print *, 'q', q(ig,l)128 ! ENDIF129 120 ENDDO 130 121 ENDDO -
trunk/LMDZ.GENERIC/libf/phystd/thermcell_env.F90
r2130 r2143 51 51 ! -------- 52 52 53 REAL zqt(ngrid,nlay) ! q 54 REAL zql(ngrid,nlay) ! q 53 REAL zqt(ngrid,nlay) ! qt environment 54 REAL zql(ngrid,nlay) ! ql environment 55 55 REAL zt(ngrid,nlay) ! T environment 56 56 REAL ztv(ngrid,nlay) ! TRPV environment … … 76 76 zhl(:,:) = pt(:,:) / zpopsk(:,:) 77 77 78 zqt(:,:) = pq(:,:,igcm_h2o_vap) 79 zql(:,:) = 0. 80 78 81 !=============================================================================== 79 82 ! Condensation and latent heat release … … 81 84 82 85 IF (water) THEN 83 84 zqt(:,:) = pq(:,:,igcm_h2o_vap)85 86 86 87 DO l=1,nlay … … 101 102 ELSE 102 103 103 zqt(:,:) = 0.104 zql(:,:) = 0.105 104 zt(:,:) = pt(:,:) 106 ztv(:,:) = zhl(:,:)105 ztv(:,:) = pt(:,:) / zpopsk(:,:) 107 106 108 107 ENDIF -
trunk/LMDZ.GENERIC/libf/phystd/thermcell_flux.F90
r2132 r2143 4 4 SUBROUTINE thermcell_flux(ngrid,nlay,ptimestep,masse, & 5 5 lmin,lmax,entr_star,detr_star, & 6 f,rhobarz,zlev,zw2,fm,entr,detr ,zqla)7 8 9 !=============================================================================== 10 ! Purpose: fluxes deduction6 f,rhobarz,zlev,zw2,fm,entr,detr) 7 8 9 !=============================================================================== 10 ! Purpose: deduction des flux 11 11 ! 12 12 ! Modif 2019/04 (AB alexandre.boissinot@lmd.jussieu.fr) … … 39 39 REAL rhobarz(ngrid,nlay) 40 40 REAL f(ngrid) 41 REAL zqla(ngrid,nlay)42 REAL zmax(ngrid)43 41 44 42 ! Outputs: … … 53 51 54 52 INTEGER ig, l, k 55 INTEGER igout, lout ! Error grid point numberand level53 INTEGER igout, lout ! Error grid point and level 56 54 57 55 REAL zfm ! Mass flux such as updraft fraction is equal to its maximal authorized value alpha_max … … 62 60 REAL ddd ! ddd0 - eee 63 61 REAL zzz ! Temporary variable set to mass flux 64 REAL fact ! Factor between old norm and new norm 65 66 INTEGER ncorecdetr ! detr > fm- counter 67 INTEGER ncorecentr ! entr > e_max counter 68 INTEGER ncorecalpha ! fm > zfm counter 62 REAL fact 63 REAL test 64 65 INTEGER ncorecentr 66 INTEGER ncorecdetr 67 INTEGER nerrorequa 68 INTEGER ncorecfact 69 INTEGER ncorecalpha 69 70 70 71 LOGICAL labort_physic 71 72 72 CHARACTER (len=20) :: modname='thermcell_flux'73 CHARACTER (len=80) :: abort_message74 75 73 !=============================================================================== 76 74 ! Initialization 77 75 !=============================================================================== 78 76 79 ncorecdetr = 0 80 ncorecentr = 0 77 nerrorequa = 0 78 ncorecentr = 0 79 ncorecdetr = 0 80 ncorecfact = 0 81 81 ncorecalpha = 0 82 82 … … 85 85 fm(:,:) = 0. 86 86 87 !=============================================================================== 88 ! Mass flux, entrainment and detrainment 87 labort_physic = .false. 88 89 fact = 0. 90 91 !=============================================================================== 92 ! Calcul de l'entrainement, du detrainement et du flux de masse 89 93 !=============================================================================== 90 94 … … 99 103 100 104 !------------------------------------------------------------------------------- 101 ! Mass flux and boundary conditions105 ! Calcul du flux de masse 102 106 !------------------------------------------------------------------------------- 103 107 104 108 DO l=1,nlay 105 109 DO ig=1,ngrid 106 IF (l .lt.lmax(ig) .and. l.ge.lmin(ig)) THEN110 IF (l < lmax(ig) .and. l >= lmin(ig)) THEN 107 111 fm(ig,l+1) = fm(ig,l) + entr(ig,l) - detr(ig,l) 108 ELSEIF (l .eq.lmax(ig)) THEN112 ELSEIF (l == lmax(ig)) THEN 109 113 fm(ig,l+1) = 0. 110 114 detr(ig,l) = fm(ig,l) + entr(ig,l) … … 118 122 119 123 !=============================================================================== 120 ! Validity tests and corrections124 ! Checking 121 125 !=============================================================================== 122 126 … … 127 131 !------------------------------------------------------------------------------- 128 132 129 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 130 ! AB : I remove the correction and replace it by an uncompromising test. 131 ! According to the previous derivations, we MUST have positive mass flux 132 ! everywhere! Indeed, as soon as fm becomes negative, the plume stops. 133 ! Then the only value which can be negative is the mass flux at the top of 134 ! the plume. However, this value was set to zero a few lines above. 135 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 136 137 labort_physic=.false. 138 139 DO ig=1,ngrid 140 IF (fm(ig,l).lt.0.) THEN 133 DO ig=1,ngrid 134 IF (fm(ig,l) < 0.) THEN 141 135 labort_physic = .true. 142 136 igout = ig … … 145 139 ENDDO 146 140 147 IF (labort_physic) THEN148 print *, 'ERROR: mass flux has negative value(s)!'149 print *, 'ig,l,fm',igout, lout, fm(igout,lout)150 print *, 'lmin,lmax', lmin(igout), lmax(igout)151 print *, '-------------------------------'152 DO k=nlay,1,-1153 print *, 'fm ', fm(igout,k+1)154 print *, 'entr,detr', entr(igout,k), detr(igout,k)155 ENDDO156 print *, 'fm-', fm(igout,1)157 print *, '-------------------------------'158 abort_message = 'Negative incoming fm in thermcell_flux'159 CALL abort_physic(modname,abort_message,1)160 ENDIF161 162 141 !------------------------------------------------------------------------------- 163 142 ! Is entrainment positive ? 164 143 !------------------------------------------------------------------------------- 165 144 166 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 167 ! AB : According to the previous derivations, we MUST have positive entrainment 168 ! and detrainment everywhere! 169 ! Indeed, they are set to max between zero and a computed value. 170 ! Then tests are uncompromising. 171 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 172 173 DO ig=1,ngrid 174 IF (entr(ig,l).lt.0.) THEN 145 DO ig=1,ngrid 146 IF (entr(ig,l) < 0.) THEN 175 147 labort_physic = .true. 176 148 igout = ig … … 179 151 ENDDO 180 152 153 !------------------------------------------------------------------------------- 154 ! Is detrainment positive ? 155 !------------------------------------------------------------------------------- 156 157 DO ig=1,ngrid 158 IF (detr(ig,l) < 0.) THEN 159 labort_physic = .true. 160 igout = ig 161 lout = l 162 ENDIF 163 ENDDO 164 165 !------------------------------------------------------------------------------- 166 ! Abort 167 !------------------------------------------------------------------------------- 168 181 169 IF (labort_physic) THEN 182 print *, 'ERROR: entrainment has negative value(s)!' 170 print *, '---------------------------------------------------------' 171 print *, 'ERROR: mass flux has negative value(s)!' 183 172 print *, 'ig,l,entr', igout, lout, entr(igout,lout) 184 173 print *, 'lmin,lmax', lmin(igout), lmax(igout) 185 print *, '- ------------------------------'174 print *, '- - - - - - - - - - - - - - - - - - - - - - - - - - - - -' 186 175 DO k=nlay,1,-1 187 176 print *, 'fm ', fm(igout,k+1) … … 189 178 ENDDO 190 179 print *, 'fm ', fm(igout,1) 191 print *, '-------------------------------' 192 abort_message = 'Negative entrainment in thermcell_flux' 193 CALL abort_physic(modname,abort_message,1) 194 ENDIF 195 196 !------------------------------------------------------------------------------- 197 ! Is detrainment positive ? 198 !------------------------------------------------------------------------------- 199 200 DO ig=1,ngrid 201 IF (detr(ig,l).lt.0.) THEN 202 labort_physic = .true. 203 igout = ig 204 lout = l 205 ENDIF 206 ENDDO 207 208 IF (labort_physic) THEN 209 print *, 'ERROR: detrainment has negative value(s)!' 210 print *, 'ig,l,detr', igout, lout, detr(igout,lout) 211 print *, 'lmin,lmax', lmin(igout), lmax(igout) 212 print *, '-------------------------------' 213 DO k=nlay,1,-1 214 print *, 'fm ', fm(igout,k+1) 215 print *, 'entr,detr', entr(igout,k), detr(igout,k) 216 ENDDO 217 print *, 'fm ', fm(igout,1) 218 print *, '-------------------------------' 219 abort_message = 'Negative detrainment in thermcell_flux' 220 CALL abort_physic(modname,abort_message,1) 221 ENDIF 222 223 !------------------------------------------------------------------------------- 224 ! Is detrainment lesser than incoming mass flux ? 180 print *, '---------------------------------------------------------' 181 CALL abort 182 ENDIF 183 184 !------------------------------------------------------------------------------- 185 ! Is detrainment lessser than incoming mass flux ? 225 186 !------------------------------------------------------------------------------- 226 187 … … 238 199 detr(ig,l) = fm(ig,l) 239 200 entr(ig,l) = fm(ig,l+1) 240 241 IF (prt_level.ge.10) THEN242 print *, 'WARNING: Detrainment is modified in thermcell_flux!'243 print *, 'ig,l,lmax', ig, l, lmax(ig)244 ENDIF245 246 IF (prt_level.ge.100) THEN247 print *, 'fm+', fm(ig,l+1)248 print *, 'entr,detr', entr(ig,l), detr(ig,l)249 print *, 'fm-', fm(ig,l)250 ENDIF251 252 201 ncorecdetr = ncorecdetr + 1 253 254 ENDIF 255 ENDDO 256 257 !------------------------------------------------------------------------------- 258 ! Is entrainment mass fraction lesser than fomass_max ? 202 ENDIF 203 ENDDO 204 205 !------------------------------------------------------------------------------- 206 ! Is entrained mass lesser than fomass_max ? 259 207 !------------------------------------------------------------------------------- 260 208 … … 271 219 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 272 220 273 labort_physic=.false.274 275 221 DO ig=1,ngrid 276 222 eee0 = entr(ig,l) … … 278 224 eee = entr(ig,l) - masse(ig,l) * fomass_max / ptimestep 279 225 ddd = detr(ig,l) - eee 280 281 IF (eee.gt.0.) THEN 226 IF (eee > 0.) THEN 282 227 entr(ig,l) = entr(ig,l) - eee 283 228 ncorecentr = ncorecentr + 1 284 285 IF (prt_level.ge.10) THEN 286 print *, 'WARNING: Entrainment is modified in thermcell_flux!' 287 print *, 'ig,l,lmax', ig, l, lmax(ig) 288 ENDIF 289 290 IF (ddd.gt.0.) THEN ! detr in the current layer is reduced 229 IF (ddd > 0.) THEN 291 230 detr(ig,l) = ddd 292 ELSEIF (l .eq.lmax(ig)) THEN ! detr in the last layer is adjusted231 ELSEIF (l == lmax(ig)) THEN 293 232 detr(ig,l) = fm(ig,l) + entr(ig,l) 294 ELSEIF (entr(ig,l+1) .gt.ABS(ddd)) THEN ! detr in the current layer is set to 0 and entr in the layer above is reduced233 ELSEIF (entr(ig,l+1) > ABS(ddd)) THEN 295 234 detr(ig,l) = 0. 296 235 fm(ig,l+1) = fm(ig,l) + entr(ig,l) … … 298 237 ! ELSE 299 238 ! entr(ig,l) = entr(ig,l) + eee 300 ! igout =ig301 ! lout =l302 ! labort_physic =.true.239 ! igout = ig 240 ! lout = l 241 ! labort_physic = .true. 303 242 ELSE 304 fact = (eee0 - eee) / eee0243 fact = max(fact, eee0 / (eee0 - eee)) 305 244 entr(ig,l) = eee0 306 detr(ig,:) = detr(ig,:) * fact 307 entr(ig,:) = entr(ig,:) * fact 308 fm(ig,:) = fm(ig,:) * fact 309 310 IF (prt_level.ge.1) THEN 311 print *, 'WARNING: Normalisation is modified in thermcell_flux!' 312 print *, 'old f, new f :', f(ig), fact * f(ig) 313 ENDIF 245 ncorecfact = ncorecfact + 1 314 246 ENDIF 315 247 ENDIF 316 248 ENDDO 317 249 318 ! IF (labort_physic) THEN 319 ! print *, 'ERROR: Entrainment is greater than its maximal authorized value!' 320 ! print *, ' Nor detrainment neither entrainment can compensate it!' 321 ! print *, 'ig,l,entr', igout, lout, entr(igout,lout) 322 ! print *, 'lmin,lmax', lmin(igout), lmax(igout) 323 ! print *, '--------------------' 324 ! print *, 'e_max :', masse(igout,lout)*fomass_max/ptimestep 325 ! print *, ' fomass_max :', fomass_max 326 ! print *, ' masse :', masse(igout,lout) 327 ! print *, ' ptimestep :', ptimestep 328 ! print *, 'norm :', f(igout) 329 ! print *, 'entr* :', entr_star(igout,lout) 330 ! print *, '--------------------' 331 ! DO k=nlay,1,-1 332 ! print *, 'fm ', fm(igout,k+1) 333 ! print *, 'entr,detr', entr(igout,k), detr(igout,k) 334 ! ENDDO 335 ! print *, 'fm ', fm(igout,1) 336 ! print *, '-------------------------------' 337 ! abort_message = 'Entrainment is too large in thermcell_flux' 338 ! CALL abort_physic (modname,abort_message,1) 339 ! ENDIF 340 341 !------------------------------------------------------------------------------- 342 ! Is updraft fraction lesser than alpha_max ? 343 !------------------------------------------------------------------------------- 344 345 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 346 ! FH : Partie a revisiter. 347 ! Il semble qu'etaient codees ici deux optiques dans le cas F/(rho*w) > 1 348 ! - soit limiter la hauteur du thermique en considerant que c'est 349 ! la derniere chouche 350 ! - soit limiter F a rho w. 351 ! Dans le second cas, il faut en fait limiter a un peu moins que ca parce 352 ! qu'on a des 1/(1-alpha) un peu plus loin dans thermcell_main et qu'il 353 ! semble de toutes facons deraisonable d'avoir des fractions de 1... 354 ! Ci-dessous, et dans l'etat actuel, le premier des deux if est sans doute inutile. 355 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 250 IF (labort_physic) THEN 251 print *, '---------------------------------------------------------' 252 print *, 'ERROR: Entrainment is greater than maximal authorized value!' 253 print *, ' Nor detrainment neither entrainment can compensate it!' 254 print *, 'ig,l,entr', igout, lout, entr(igout,lout) 255 print *, 'lmin,lmax', lmin(igout), lmax(igout) 256 print *, '- - - - - - - - - - - - - - - - - - - - - - - - - - - - -' 257 print *, 'e_max :', masse(igout,lout)*fomass_max/ptimestep 258 print *, ' fomass_max :', fomass_max 259 print *, ' masse :', masse(igout,lout) 260 print *, ' ptimestep :', ptimestep 261 print *, 'norm :', f(igout) 262 print *, 'entr* :', entr_star(igout,lout) 263 print *, '- - - - - - - - - - - - - - - - - - - - - - - - - - - - -' 264 DO k=nlay,1,-1 265 print *, 'fm ', fm(igout,k+1) 266 print *, 'entr,detr', entr(igout,k), detr(igout,k) 267 ENDDO 268 print *, 'fm ', fm(igout,1) 269 print *, '---------------------------------------------------------' 270 CALL abort 271 ENDIF 272 273 !------------------------------------------------------------------------------- 274 ! Is Updraft fraction lesser than alpha_max ? 275 !------------------------------------------------------------------------------- 356 276 357 277 DO ig=1,ngrid 358 278 zfm = rhobarz(ig,l+1) * zw2(ig,l+1) * alpha_max 359 279 360 IF (fm(ig,l+1) .gt.zfm) THEN280 IF (fm(ig,l+1) > zfm) THEN 361 281 f_old = fm(ig,l+1) 362 282 fm(ig,l+1) = zfm 363 283 detr(ig,l) = detr(ig,l) + f_old - fm(ig,l+1) 364 284 ncorecalpha = ncorecalpha + 1 365 285 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 366 286 ! AB : The previous change doesn't observe the equation df/dz = e - d. 367 ! To avoid this issue, what is better to do? I choose to increase the 368 ! entrainment in the layer above when l<lmax. If l=lmax then fm(l+1)=0 and 369 ! there are never any problems. 370 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 371 IF (l.lt.lmax(ig)) THEN 372 entr(ig,l+1) = entr(ig,l+1) + f_old - fm(ig,l+1) 373 ENDIF 374 375 IF (prt_level.ge.10) THEN 376 print *, 'Mass flux is modified in thermcell_flux' 377 print *, 'ig,l,lmax', ig, l, lmax(ig) 378 ENDIF 379 380 IF (prt_level.ge.100) THEN 381 print *, 'fm-', fm(ig,l) 382 print *, 'entr,detr', entr(ig,l), detr(ig,l) 383 print *, 'fm+', fm(ig,l+1) 384 ENDIF 385 386 ncorecalpha = ncorecalpha + 1 287 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 288 entr(ig,l+1) = entr(ig,l+1) + f_old - fm(ig,l+1) 387 289 ENDIF 388 290 ENDDO … … 390 292 ENDDO 391 293 392 !------------------------------------------------------------------------------- 393 ! Boundary conditions 394 !------------------------------------------------------------------------------- 395 396 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 397 ! AB : test added to avoid problem when lmax = 0, which is not a fortran index. 398 ! Theoretically, we always get lmax >= lmin >= linf > 0 399 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 294 IF (fact > 0.) THEN 295 entr(:,:) = entr(:,:) / fact 296 detr(:,:) = entr(:,:) / fact 297 fm(:,:) = entr(:,:) / fact 298 ENDIF 299 300 !------------------------------------------------------------------------------- 301 ! Is equation df/dz = e - d still verified ? 302 !------------------------------------------------------------------------------- 303 304 ! DO l=1,nlay 305 ! DO ig=1,ngrid 306 ! test = abs(fm(ig,l) + entr(ig,l) - detr(ig,l) - fm(ig,l+1)) 307 ! IF (test > 1.e-10) THEN 308 ! nerrorequa = nerrorequa + 1 309 ! ENDIF 310 ! ENDDO 311 ! ENDDO 312 313 !------------------------------------------------------------------------------- 314 ! Reset top boundary conditions 315 !------------------------------------------------------------------------------- 316 400 317 DO ig=1,ngrid 401 IF (lmax(ig). gt.0) THEN318 IF (lmax(ig).GT.0) THEN 402 319 detr(ig,lmax(ig)) = fm(ig,lmax(ig)) + entr(ig,lmax(ig)) 403 320 fm(ig,lmax(ig)+1) = 0. … … 406 323 ENDDO 407 324 408 ! -------------------------------------------------------------------------------409 ! Corrections display410 ! -------------------------------------------------------------------------------411 412 IF (prt_level .GE.1) THEN413 414 IF (ncorecdetr. GE.1) THEN325 !=============================================================================== 326 ! Outputs 327 !=============================================================================== 328 329 IF (prt_level > 0) THEN 330 331 IF (ncorecdetr.ge.1) THEN 415 332 print *, 'WARNING: Detrainment is greater than mass flux!' 416 print *, ' In', ncorecdetr, 'point(s)'417 ENDIF 418 419 IF (ncorecentr. GE.1) THEN333 print *, 'In', ncorecdetr, 'grid point(s) over', nlay, 'x', ngrid 334 ENDIF 335 336 IF (ncorecentr.ge.1) THEN 420 337 print *, 'WARNING: Entrained mass is greater than maximal authorized value!' 421 print *, ' In', ncorecentr, 'point(s)' 422 ENDIF 423 424 IF (ncorecalpha.GE.1) THEN 338 print *, 'in', ncorecentr, 'grid point(s) over', nlay, 'x', ngrid 339 ENDIF 340 341 IF (ncorecfact.ge.1) THEN 342 print *, 'WARNING: Entrained mass needs renormalization to be ok!' 343 print *, 'in', ncorecfact, 'grid point(s) over', nlay, 'x', ngrid 344 ENDIF 345 346 ! IF (nerrorequa.ge.1) THEN 347 ! print *, 'WARNING: !' 348 ! print *, 'in', nerrorequa, 'grid point(s) over', nlay, 'x', ngrid 349 ! ENDIF 350 351 IF (ncorecalpha.ge.1) THEN 425 352 print *, 'WARNING: Updraft fraction is greater than maximal authorized value!' 426 print *, ' In', ncorecalpha, 'point(s)'353 print *, 'in', ncorecalpha, 'grid point(s) over', nlay, 'x', ngrid 427 354 ENDIF 428 355 429 356 ENDIF 430 431 ! AB : temporary test added to check the validity of eq. df/dz = e - d432 ! DO l=1,nlay433 ! DO ig=1,ngrid434 ! test = abs(fm(ig,l) + entr(ig,l) - detr(ig,l) - fm(ig,l+1))435 ! IF (test.gt.1.e-10) THEN436 ! print *, 'WARNING: df/dz != entr - detr'437 ! print *, 'ig,l', ig, l438 ! print *, 'fm+ :', fm(ig,l+1)439 ! print *, 'entr,detr', entr(ig,l), detr(ig,l)440 ! print *, 'fm :', fm(ig,l)441 ! print *, 'entr,detr', entr(ig,l-1), detr(ig,l-1)442 ! print *, 'fm- :', fm(ig,l-1)443 ! print *, 'err. :', test444 ! ENDIF445 ! ENDDO446 ! ENDDO447 357 448 358 -
trunk/LMDZ.GENERIC/libf/phystd/thermcell_height.F90
r2127 r2143 2 2 ! 3 3 ! 4 SUBROUTINE thermcell_height(ngrid,nlay,lmin,linter,lmix, zw2,&5 zlev, lmax,zmax,zmix,wmax,f_star)4 SUBROUTINE thermcell_height(ngrid,nlay,lmin,linter,lmix,lmax,zw2, & 5 zlev,zmin,zmix,zmax,wmax,f_star) 6 6 7 7 … … 9 9 ! Purpose: calcul des caracteristiques du thermique: zmax,wmax,zmix 10 10 !=============================================================================== 11 12 USE thermcell_mod13 USE print_control_mod, ONLY: prt_level14 11 15 12 IMPLICIT NONE … … 36 33 37 34 REAL linter(ngrid) 38 REAL zmax(ngrid) ! maximal vertical speed 39 REAL zmix(ngrid) ! altitude of maximal vertical speed 40 REAL wmax(ngrid) ! maximal vertical speed 41 REAL zw2(ngrid,nlay+1) ! square vertical speed (becomes its squere root) 35 REAL zmin(ngrid) ! Altitude of the plume bottom 36 REAL zmax(ngrid) ! Altitude of the plume top 37 REAL zmix(ngrid) ! Altitude of maximal vertical speed 38 REAL wmax(ngrid) ! Maximal vertical speed 39 REAL zw2(ngrid,nlay+1) ! Square vertical speed (becomes its square root) 42 40 43 41 ! Local: … … 46 44 INTEGER ig, l 47 45 48 REAL num(ngrid)49 REAL denom(ngrid)50 46 REAL zlevinter(ngrid) 51 REAL zlev2(ngrid,nlay+1)52 47 53 48 !=============================================================================== … … 69 64 ENDDO 70 65 71 DO ig=1,ngrid 72 zlev2(ig,:) = zlev(ig,:) - zlev(ig,lmin(ig)) 73 ENDDO 74 75 !=============================================================================== 76 ! Calcul de la hauteur max du thermique 77 !=============================================================================== 66 !=============================================================================== 67 ! Thermal plume height 68 !=============================================================================== 69 70 !------------------------------------------------------------------------------- 71 ! Calcul de zmin 72 !------------------------------------------------------------------------------- 73 74 DO l=1,nlay 75 DO ig=1,ngrid 76 zmin(ig) = zlev(ig,lmin(ig)) 77 ENDDO 78 ENDDO 78 79 79 80 !------------------------------------------------------------------------------- … … 83 84 DO ig=1,ngrid 84 85 DO l=nlay,lmin(ig)+1,-1 85 IF ( zw2(ig,l).le.0..or.f_star(ig,l).le.0.) THEN86 IF ((zw2(ig,l) <= 0.).or.(f_star(ig,l) <= 0.)) THEN 86 87 lmax(ig) = l - 1 87 88 ENDIF … … 90 91 91 92 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 92 ! On traite le cas particulier qu'il faudrait eviter ou le thermique 93 ! atteind le haut du modele ... 94 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 95 DO ig=1,ngrid 96 IF (zw2(ig,nlay).gt.0.) THEN 93 ! AB: Problematic case where thermal plume reaches the top of the model... 94 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 95 DO ig=1,ngrid 96 IF (zw2(ig,nlay) > 0.) THEN 97 97 print *, 'WARNING: a thermal plume reaches the highest layer!' 98 98 print *, 'ig,l', ig, nlay … … 110 110 ! Otherwise we have lmax>lmin. 111 111 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 112 113 112 DO ig=1,ngrid 114 113 l = lmax(ig) 115 116 IF (l==lmin(ig)) THEN 114 IF (l == lmin(ig)) THEN 117 115 linter(ig) = l 118 116 ELSE … … 121 119 ENDDO 122 120 123 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 124 ! AB : zmax is now equal to zlevinter-zlev(lmin) because lmin can be different 125 ! from 1. 126 ! We have to substract this level because zmax must be the plume height (to 127 ! derive the good closure), not the plume highest altitude. 128 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 129 130 DO ig=1,ngrid 131 zlevinter(ig) = (zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*linter(ig) & 132 & + zlev(ig,lmax(ig))-lmax(ig)*(zlev(ig,lmax(ig)+1) & 133 & - zlev(ig,lmax(ig))) 134 zmax(ig) = max(zmax(ig), zlevinter(ig) - zlev(ig,lmin(ig))) 135 ENDDO 136 137 !=============================================================================== 138 ! Calcul de wmax et du niveau d'inversion. 121 DO ig=1,ngrid 122 zlevinter(ig) = zlev(ig,lmax(ig)) + (linter(ig) - lmax(ig)) & 123 & * (zlev(ig,lmax(ig)+1) - zlev(ig,lmax(ig))) 124 zmax(ig) = max(zmax(ig), zlevinter(ig)) 125 ENDDO 126 127 !=============================================================================== 128 ! Thermal plume maximal speed and inversion layer 139 129 !=============================================================================== 140 130 … … 145 135 DO l=1,nlay 146 136 DO ig=1,ngrid 147 IF( l.le.lmax(ig) .and. l.gt.lmin(ig)) THEN137 IF((l <= lmax(ig)).and.(l > lmin(ig))) THEN 148 138 IF (zw2(ig,l).lt.0.) THEN 149 139 print *, 'ERROR: zw2 has negative value(s)!' … … 158 148 zw2(ig,l) = sqrt(zw2(ig,l)) 159 149 160 IF (zw2(ig,l) .gt.wmax(ig)) THEN150 IF (zw2(ig,l) > wmax(ig)) THEN 161 151 wmax(ig) = zw2(ig,l) 162 152 lmix(ig) = l … … 172 162 !------------------------------------------------------------------------------- 173 163 174 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 175 ! AB : We substract zlev(lmin) in zmix formula because we have to 176 ! compare it with zmax which is zlev(linter)-zlev(lmin). 177 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 178 179 DO ig=1,ngrid 180 IF (lmix(ig).gt.lmin(ig)) THEN 164 DO ig=1,ngrid 165 IF (lmix(ig) > lmin(ig)) THEN 181 166 IF (((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig))) & 182 167 & *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1))) & 183 168 & -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1)) & 184 & *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))) .gt.1e-10)&169 & *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))) > 1e-10) & 185 170 & THEN 186 171 zmix(ig) = ((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig))) & … … 191 176 & *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1))) & 192 177 & -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1)) & 193 & *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig)))))) & 194 & - zlev(ig,lmin(ig)) 178 & *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig)))))) 195 179 ELSE 196 zmix(ig) = zlev(ig,lmix(ig)) - zlev(ig,lmin(ig))180 zmix(ig) = zlev(ig,lmix(ig)) 197 181 print *, 'WARNING: problematic zmix value!' 198 182 ENDIF … … 203 187 204 188 !=============================================================================== 205 ! Verification zmax > zmix 206 !=============================================================================== 207 208 DO ig=1,ngrid 209 IF ((zmax(ig)-zmix(ig)).lt.0.) THEN 189 ! Check consistence between zmax and zmix 190 !=============================================================================== 191 192 !------------------------------------------------------------------------------- 193 ! 194 !------------------------------------------------------------------------------- 195 196 DO ig=1,ngrid 197 IF ((zmax(ig)-zmix(ig)) < 0.) THEN 198 print *, 'WARNING: we have zmix > zmax!' 199 print *, 'zmax,zmix_old,zmix_new', zmax(ig), zmix(ig), 0.9 * zmax(ig) 210 200 zmix(ig) = 0.9 * zmax(ig) 211 print *, 'WARNING: we have zmix > zmax!'212 print *, 'zmax,zmix', zmax(ig), zmix(ig)213 201 ENDIF 214 202 ENDDO … … 220 208 DO ig=1,ngrid 221 209 DO l=1,nlay 222 IF ( zmix(ig).ge.zlev2(ig,l).and.zmix(ig).lt.zlev2(ig,l+1)) THEN210 IF ((zmix(ig) >= zlev(ig,l)).and.(zmix(ig) < zlev(ig,l+1))) THEN 223 211 lmix(ig) = l 224 212 ENDIF -
trunk/LMDZ.GENERIC/libf/phystd/thermcell_main.F90
r2132 r2143 49 49 50 50 USE thermcell_mod 51 USE tracer_h, ONLY: igcm_h2o_vap 52 USE print_control_mod, ONLY: lunout, prt_level 51 USE print_control_mod, ONLY: prt_level 53 52 54 53 IMPLICIT NONE … … 96 95 INTEGER lmix(ngrid) ! Layer in which plume vertical speed is maximal 97 96 INTEGER lmin(ngrid) ! First unstable layer 97 98 REAL zmix(ngrid) ! Altitude of maximal vertical speed 99 REAL zmax(ngrid) ! Maximal altitudes where plumes are active 100 REAL zmin(ngrid) ! Minimal altitudes where plumes are active 98 101 99 102 REAL zlay(ngrid,nlay) ! Layers altitudes … … 124 127 REAL zqa(ngrid,nlay,nq) ! q plume (ql=0, qv=qt) 125 128 129 REAL f_star(ngrid,nlay+1) ! Normalized mass flux 130 REAL entr_star(ngrid,nlay) ! Normalized entrainment (E* = e* dz) 131 REAL detr_star(ngrid,nlay) ! Normalized detrainment (D* = d* dz) 132 133 REAL fm(ngrid,nlay+1) ! Mass flux 134 REAL entr(ngrid,nlay) ! Entrainment (E = e dz) 135 REAL detr(ngrid,nlay) ! Detrainment (D = d dz) 136 137 REAL f(ngrid) ! Mass flux norm 138 REAL lambda ! Time relaxation coefficent 139 REAL fraca(ngrid,nlay+1) ! Updraft fraction 126 140 REAL linter(ngrid) ! Level (continuous) of maximal vertical speed 127 REAL zmix(ngrid) ! Altitude of maximal vertical speed128 REAL zmax(ngrid) ! Plume height129 141 REAL wmax(ngrid) ! Maximal vertical speed 130 142 REAL zw2(ngrid,nlay+1) ! Plume vertical speed 131 132 REAL fraca(ngrid,nlay+1) ! Updraft fraction133 134 REAL f_star(ngrid,nlay+1) ! Normalized mass flux135 REAL entr_star(ngrid,nlay) ! Normalized entrainment136 REAL detr_star(ngrid,nlay) ! Normalized detrainment137 138 REAL f(ngrid) ! Mass flux norm139 REAL fm(ngrid,nlay+1) ! Mass flux140 REAL entr(ngrid,nlay) ! Entrainment141 REAL detr(ngrid,nlay) ! Detrainment142 143 REAL lambda ! Time relaxation coefficent144 145 143 REAL zdthladj(ngrid,nlay) ! Potential temperature variations 146 144 REAL dummy(ngrid,nlay) ! Dummy argument for thermcell_dq() 147 148 CHARACTER (len=20) :: modname='thermcell_main'149 CHARACTER (len=80) :: abort_message150 145 151 146 !=============================================================================== … … 178 173 pdqadj(:,:,:) = 0.0 179 174 180 ! AB: Careful, hard-coded value from Earth tuned version of the thermal plume model!181 175 DO ig=1,ngrid 182 f0(ig) = max(f0(ig), 1.e-2) 183 ENDDO 184 185 IF (prt_level.ge.20) then 186 DO ig=1,ngrid 187 print *, 'ig,f0', ig, f0(ig) 188 ENDDO 189 ENDIF 176 ! AB: Careful: Hard-coded value from Earth version! 177 ! f0(ig) = max(f0(ig), 1.e-2) 178 ! AB: No pescribed minimal value for f0 179 f0(ig) = max(f0(ig), 0.) 180 ENDDO 190 181 191 182 !=============================================================================== … … 212 203 213 204 DO l=1,nlay 214 zlay(:,l) = pphi(:,l) /RG205 zlay(:,l) = pphi(:,l) / RG 215 206 ENDDO 216 207 … … 222 213 223 214 IF (prt_level.ge.10) THEN 224 write(lunout,*) 'WARNING: thermcell_main rhobarz(:,1)=rho(:,1)' 215 print *, 'WARNING: density in the first layer is equal to density at the first level!' 216 print *, 'rhobarz(:,1)', rhobarz(:,1) 217 print *, 'rho(:,1) ', rho(:,1) 225 218 ENDIF 226 219 … … 319 312 & zhl,zqt,zql,rhobarz,zlev,pplev,pphi,zpopsk, & 320 313 & detr_star,entr_star,f_star, & 321 & ztva,zhla,zqla,zqta,zta, 322 & zw2, zqsa,lmix,lmin)323 314 & ztva,zhla,zqla,zqta,zta,zqsa, & 315 & zw2,lmix,lmin) 316 324 317 !------------------------------------------------------------------------------- 325 318 ! Thermal plumes characteristics: zmax, zmix, wmax … … 327 320 328 321 ! AB: Careful, zw2 became its square root in thermcell_height! 329 CALL thermcell_height(ngrid,nlay,lmin,linter,lmix, zw2,&330 & zlev, lmax,zmax,zmix,wmax,f_star)322 CALL thermcell_height(ngrid,nlay,lmin,linter,lmix,lmax,zw2, & 323 & zlev,zmin,zmix,zmax,wmax,f_star) 331 324 332 325 !=============================================================================== … … 339 332 340 333 CALL thermcell_closure(ngrid,nlay,ptimestep,rho,zlev, & 341 & lmax,entr_star,zm ax,wmax,f)334 & lmax,entr_star,zmin,zmax,wmax,f) 342 335 343 336 IF (tau_thermals>1.) THEN … … 352 345 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 353 346 IF (.not. (f0(1).ge.0.) ) THEN 354 abort_message = '.not. (f0(1).ge.0.)'347 print *, 'ERROR: mass flux norm is not positive!' 355 348 print *, 'f0 =', f0(1) 356 CALL abort _physic(modname,abort_message,1)349 CALL abort 357 350 ENDIF 358 351 … … 393 386 DO l=2,nlay 394 387 DO ig=1,ngrid 395 IF (zw2(ig,l) .gt.0.) THEN388 IF (zw2(ig,l) > 0.) THEN 396 389 fraca(ig,l) = fm(ig,l) / (rhobarz(ig,l) * zw2(ig,l)) 397 390 ELSE … … 433 426 CALL thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,detr0,masse, & 434 427 & zu,pduadj,zua,lmin) 435 428 436 429 CALL thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,detr0,masse, & 437 430 & zv,pdvadj,zva,lmin) 438 431 432 ! CALL thermcell_dv2(ngrid,nlay,ptimestep,fm0,entr0,detr0,masse,fraca, & 433 ! & zmax,zmin,pu,pv,pduadj,pdvadj,zua,zva) 439 434 440 435 RETURN -
trunk/LMDZ.GENERIC/libf/phystd/thermcell_mod.F90
r2132 r2143 5 5 6 6 ! Flags for computations 7 8 INTEGER, PARAMETER :: dqimpl = 1! 1 flag for thermcell_dq version (1 : implicit scheme || 0 : explicit scheme)7 ! default 8 INTEGER,SAVE :: dqimpl ! 1 flag for thermcell_dq version (1 : implicit scheme || 0 : explicit scheme) 9 9 10 10 ! Physical parameters 11 11 12 REAL, PARAMETER :: r_aspect_thermals = 1.0 !Aspect ratio of the thermals (width / height)13 REAL, PARAMETER :: tau_thermals = 0.! 0. Relaxation time14 REAL, PARAMETER :: betalpha = 1.0 ! 0.9 - factorbetween 0 (e=d) and 1 (rho*fraca=cst)15 REAL, PARAMETER :: afact = 1. ! 2./3. Buoyancy contribution - factorbetween 0 and 116 REAL, PARAMETER :: fact_epsilon = 1.e-4 ! 2.e-3 Friction at plume borders - exponential decrease17 REAL, PARAMETER :: nu = 0.000 !Geometrical contributions to entrainment and detrainment18 REAL, PARAMETER :: alpha_max = 0.7! 0.7 Maximal permitted updraft fraction19 REAL, PARAMETER :: fomass_max = 0.5! 0.5 Maximal permitted outgoing layer mass fraction20 REAL, PARAMETER :: pres_limit = 2.e5 ! 1.e5 Minimal permitted pressure to trigger a thermal plume12 REAL,SAVE :: r_aspect_thermals ! 1.0 Aspect ratio of the thermals (width / height) 13 REAL,SAVE :: tau_thermals ! 0. Relaxation time 14 REAL,SAVE :: betalpha ! 0.9 - between 0 (e=d) and 1 (rho*fraca=cst) 15 REAL,SAVE :: afact ! 2./3. - buoyancy acceleration efficiency, between 0 and 1 16 REAL,SAVE :: fact_epsilon ! 2.e-3 - turbulence and friction deceleration 17 REAL,SAVE :: nu ! 0. Geometrical contributions to entrainment and detrainment 18 REAL,SAVE :: alpha_max ! 0.7 Maximal permitted updraft fraction 19 REAL,SAVE :: fomass_max ! 0.5 Maximal permitted outgoing layer mass fraction 20 REAL,SAVE :: pres_limit ! 1.e5 - 21 21 22 22 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ … … 29 29 ! ngrid. 30 30 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 31 INTEGER, PARAMETER :: linf = 1! 131 INTEGER,SAVE :: linf ! 1 32 32 33 33 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ … … 35 35 ! layer linf which can be used to force convection to begin in it. 36 36 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 37 REAL, PARAMETER :: d_temp = 0.! 0.37 REAL,SAVE :: d_temp ! 0. 38 38 39 39 ! Physical constants … … 52 52 SUBROUTINE init_thermcell_mod(g, rcp, r, pi, T_h2o_ice_liq, RV) 53 53 54 USE ioipsl_getin_p_mod, ONLY: getin_p 55 54 56 IMPLICIT NONE 55 57 … … 67 69 RD = r 68 70 71 print *, 'Plume aspect ratio ?' 72 r_aspect_thermals = 1. 73 call getin_p('r_aspect_thermals', r_aspect_thermals) 74 print *, 'r_aspect_thermals = ', r_aspect_thermals 75 76 print *, 'Relaxation time ?' 77 tau_thermals = 0. 78 call getin_p('tau_thermals', tau_thermals) 79 print *, 'tau_thermals = ', tau_thermals 80 81 print *, 'betalpha ?' 82 betalpha = 0.9 83 call getin_p('betalpha', betalpha) 84 print *, 'betalpha = ', betalpha 85 86 print *, 'Buoyancy acceleration efficiency ?' 87 afact = 2./3. 88 call getin_p('afact', afact) 89 print *, 'afact = ', afact 90 91 print *, 'Turbulence and friction factor ?' 92 fact_epsilon = 2.e-3 93 call getin_p('fact_epsilon', fact_epsilon) 94 print *, 'fact_epsilon = ', fact_epsilon 95 96 print *, 'Constant entrainment and detrainment term ?' 97 nu = 0. 98 call getin_p('nu', nu) 99 print *, 'nu = ', nu 100 101 print *, 'Maximal authorized updraft fraction ?' 102 alpha_max = 0.7 103 call getin_p('alpha_max', alpha_max) 104 print *, 'alpha_max = ', alpha_max 105 106 print *, 'Maximal authorized entrained mass flux ?' 107 fomass_max = 0.5 108 call getin_p('fomass_max', fomass_max) 109 print *, 'fomass_max = ', fomass_max 110 111 print *, 'Minimal pressure below which plume can no longer trigger ?' 112 pres_limit = 1.e5 113 call getin_p('pres_limit', pres_limit) 114 print *, 'pres_limit = ', pres_limit 115 116 print *, 'Deepest layer from which a plume can rise ?' 117 linf = 1 118 call getin_p('linf', linf) 119 print *, 'linf = ', linf 120 121 print *, 'd_temp ?' 122 d_temp = 0. 123 call getin_p('d_temp', d_temp) 124 print *, 'd_temp = ', d_temp 125 69 126 RETURN 70 127 END SUBROUTINE init_thermcell_mod -
trunk/LMDZ.GENERIC/libf/phystd/thermcell_plume.F90
r2130 r2143 5 5 zhl,zqt,zql,rhobarz,zlev,pplev,pphi,zpopsk, & 6 6 detr_star,entr_star,f_star, & 7 ztva,zhla,zqla,zqta,zta, 8 zw2, zqsa,lmix,lmin)7 ztva,zhla,zqla,zqta,zta,zqsa, & 8 zw2,lmix,lmin) 9 9 10 10 … … 62 62 63 63 REAL ztva(ngrid,nlay) ! TRPV plume (after mixing) 64 REAL zhla(ngrid,nlay) ! TP plume (after mixing)64 REAL zhla(ngrid,nlay) ! TP plume ? 65 65 REAL zqla(ngrid,nlay) ! ql plume (after mixing) 66 REAL zqta(ngrid,nlay) ! qt plume (after mixing)66 REAL zqta(ngrid,nlay) ! qt plume ? 67 67 REAL zqsa(ngrid,nlay) ! qsat plume (after mixing) 68 REAL zw2(ngrid,nlay+1) ! w 2plume (after mixing)68 REAL zw2(ngrid,nlay+1) ! w plume (after mixing) 69 69 70 70 ! Local: … … 77 77 REAL zta_est(ngrid,nlay) ! TR plume (before mixing) 78 78 REAL zqsa_est(ngrid) ! qsat plume (before mixing) 79 REAL zw2_est(ngrid,nlay+1) ! w 2plume (before mixing)79 REAL zw2_est(ngrid,nlay+1) ! w plume (before mixing) 80 80 81 81 REAL zta(ngrid,nlay) ! TR plume (after mixing) … … 128 128 ztv2(:,linf) = ztv(:,linf) + d_temp 129 129 130 active(:) = .false. 131 130 132 !=============================================================================== 131 133 ! First layer computation … … 133 135 134 136 DO ig=1,ngrid 135 active(ig) = .false.136 137 l = linf 137 DO WHILE (.not.active(ig).and.(pplev(ig,l+1).GT.pres_limit).and.(l.LT.nlay)) 138 zdz = (zlev(ig,l+1) - zlev(ig,l)) 139 zw2(ig,l+1) = 2. * afact * RG * zdz & ! Do we have to divide by 1+betalpha (consider entrainment) ? 140 & * (ztv2(ig,l) - ztv2(ig,l+1)) & 141 & / (ztv2(ig,l+1) * (1. + betalpha)) 142 zw2m = zw2(ig,l+1) / 2. 143 entr_star(ig,l) = zdz * zbetalpha * (afact * RG * (ztv2(ig,l) & 144 & - ztv2(ig,l+1)) / ztv2(ig,l+1) / zw2m - fact_epsilon) 145 IF ((ztv2(ig,l).gt.ztv2(ig,l+1)).and.(entr_star(ig,l).GT.0.)) THEN 146 active(ig) = .true. 138 DO WHILE (.not.active(ig).and.(pplev(ig,l+1) > pres_limit).and.(l < nlay)) 139 zbuoy(ig,l) = RG * (ztv2(ig,l) - ztv2(ig,l+1)) / ztv2(ig,l+1) 140 zdz = zlev(ig,l+1) - zlev(ig,l) 141 zw2m = afact * zbuoy(ig,l) * zdz / (1. + betalpha) 142 ! gamma = afact * zbuoy(ig,l) - fact_epsilon * zw2m 143 ! test = gamma / zw2m - nu 144 test = zbuoy(ig,l) 145 IF (test > 0.) THEN 147 146 lmin(ig) = l 148 f_star(ig,l+1) = entr_star(ig,l) 149 zw2_est(ig,l+1) = zw2(ig,l+1) 150 ELSE 151 zw2(ig,l+1) = 0. 152 entr_star(ig,l) = 0. 147 ! entr_star(ig,l) = zdz * f_star(ig,l) * zbetalpha * gamma / zw2m - nu ! Problem because f*(ig,l) = 0 148 ! detr_star(ig,l) = f_star(ig,l) * nu ! Problem because f*(ig,l) = 0 149 ! f_star(ig,l+1) = entr_star(ig,l) - detr_star(ig,l) 150 entr_star(ig,l) = 1. 151 f_star(ig,l+1) = 1. 152 zw2_est(ig,l+1) = zw2m * 2. 153 zw2(ig,l+1) = zw2_est(ig,l+1) 153 154 ENDIF 154 155 l = l + 1 … … 163 164 164 165 !------------------------------------------------------------------------------- 165 ! Check if thermal plume is (still) active 166 !------------------------------------------------------------------------------- 167 168 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 169 ! AB: we decide here if the plume is still active or not. When the plume's 170 ! first level is reached, we set active to "true". Otherwise, it is given 171 ! by zw2, f_star and entr_star. 172 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 173 DO ig=1,ngrid 174 IF (l==lmin(ig)+1) THEN 175 active(ig) = .true. 176 ENDIF 177 178 active(ig) = active(ig) & 179 & .and. (zw2(ig,l).GT.1.e-10) & 180 & .and. (f_star(ig,l) + entr_star(ig,l)).GT.1.e-10 181 ENDDO 166 ! Is thermal plume (still) active ? 167 !------------------------------------------------------------------------------- 168 169 DO ig=1,ngrid 170 active(ig) = (active(ig).or.(l == lmin(ig)+1)) & 171 & .and.(zw2(ig,l) > 1.e-10) & 172 & .and.(f_star(ig,l) > 1.e-10) 173 ENDDO 174 175 !------------------------------------------------------------------------------- 176 ! Latent heat release (before mixing) 177 !------------------------------------------------------------------------------- 182 178 183 179 ztemp(:) = zpopsk(:,l) * zhla(:,l-1) … … 188 184 189 185 !------------------------------------------------------------------------------- 190 ! Vertical speed (before mixing between plume and environment)186 ! Vertical speed (before mixing) 191 187 !------------------------------------------------------------------------------- 192 188 193 189 DO ig=1,ngrid 194 190 IF (active(ig)) THEN 195 zqla_est(ig,l) = max(0.,zqta(ig,l-1)-zqsa_est(ig))! zqla_est set to ql plume191 zqla_est(ig,l) = MAX(0.,zqta(ig,l-1) - zqsa_est(ig)) ! zqla_est set to ql plume 196 192 zta_est(ig,l) = zhla(ig,l-1) * zpopsk(ig,l) & ! zta_est set to TR plume 197 193 & + RLvCp * zqla_est(ig,l) … … 229 225 230 226 zdz = zlev(ig,l+1) - zlev(ig,l) 231 zw2m = (zw2_est(ig,l+1) + zw2_est(ig,l)) / 2. ! AB: est-ce la bonne vitesse a utiliser ?227 zw2m = (zw2_est(ig,l+1) + zw2_est(ig,l)) / 2. 232 228 gamma = afact * zbuoy(ig,l) - fact_epsilon * zw2m 233 229 234 IF (zw2_est(ig,l) .GT.0.) THEN230 IF (zw2_est(ig,l) > 0.) THEN 235 231 test = gamma / zw2m - nu 236 232 ELSE … … 241 237 ENDIF 242 238 243 IF (test .gt.0.) THEN244 detr_star(ig,l) = zdz * nu245 entr_star(ig,l) = zdz * (zbetalpha * gamma / zw2m + nu)239 IF (test > 0.) THEN 240 detr_star(ig,l) = zdz * f_star(ig,l) * nu 241 entr_star(ig,l) = zdz * f_star(ig,l) * (zbetalpha * gamma / zw2m + nu) 246 242 ELSE 247 detr_star(ig,l) = zdz * (nu - betalpha * gamma / zw2m)248 entr_star(ig,l) = zdz * nu243 detr_star(ig,l) = zdz * f_star(ig,l) * (nu - betalpha * gamma / zw2m) 244 entr_star(ig,l) = zdz * f_star(ig,l) * nu 249 245 ENDIF 250 246 251 ! IF (detr_star(ig,l).lt.0.) THEN252 ! print *, 'WARNING: detrainment is negative!'253 ! print *, 'l,detr', l, detr_star(ig,l)254 ! ENDIF255 256 ! IF (entr_star(ig,l).lt.0.) THEN257 ! print *, 'WARNING: entrainment is negative!'258 ! print *, 'l,entr', l, entr_star(ig,l)259 ! ENDIF260 261 247 f_star(ig,l+1) = f_star(ig,l) + entr_star(ig,l) - detr_star(ig,l) 262 248 263 ! IF (f_star(ig,l+1).le.0.) THEN264 ! print *, 'WARNING: mass flux is negative!'265 ! print *, 'l,f_star', l+1, f_star(ig,l+1)266 ! ENDIF267 268 249 ENDIF 269 250 ENDDO … … 273 254 !------------------------------------------------------------------------------- 274 255 275 activetmp(:) = active(:) .and. (f_star(:,l+1).GT.1.e-10)256 activetmp(:) = active(:).and.(f_star(:,l+1) > 1.e-10) 276 257 277 258 DO ig=1,ngrid … … 286 267 ENDDO 287 268 269 !------------------------------------------------------------------------------- 270 ! Latent heat release (after mixing) 271 !------------------------------------------------------------------------------- 272 288 273 ztemp(:) = zpopsk(:,l) * zhla(:,l) 289 274 … … 295 280 296 281 !------------------------------------------------------------------------------- 297 ! Vertical speed (after mixing between plume and environment)282 ! Vertical speed (after mixing) 298 283 !------------------------------------------------------------------------------- 299 284 300 285 DO ig=1,ngrid 301 286 IF (activetmp(ig)) THEN 302 zqla(ig,l) = max(0.,zqta(ig,l)-zqsa(ig,l))! zqla is set to ql plume (mixed)303 304 zta(ig,l) = zhla(ig,l) * zpopsk(ig,l) + RLvCp * zqla(ig,l) ! ztva is set to TR plume (mixed)287 zqla(ig,l) = MAX(0.,zqta(ig,l) - zqsa(ig,l)) ! zqla is set to ql plume (mixed) 288 zta(ig,l) = zhla(ig,l) * zpopsk(ig,l) & ! ztva is set to TR plume (mixed) 289 & + RLvCp * zqla(ig,l) 305 290 ztva(ig,l) = zta(ig,l) / zpopsk(ig,l) & ! ztva is set to TRPV plume (mixed) 306 291 & * (1. + RETV*(zqta(ig,l)-zqla(ig,l)) - zqla(ig,l)) … … 329 314 330 315 331 RETURN 316 RETURN 332 317 END
Note: See TracChangeset
for help on using the changeset viewer.
