Changeset 2143 for trunk/LMDZ.GENERIC/libf/phystd/thermcell_flux.F90
- Timestamp:
- Jun 20, 2019, 4:11:54 PM (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
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
Note: See TracChangeset
for help on using the changeset viewer.
