Changeset 2132 for trunk/LMDZ.GENERIC/libf/phystd/thermcell_flux.F90
- Timestamp:
- Apr 29, 2019, 3:13:04 PM (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.GENERIC/libf/phystd/thermcell_flux.F90
r2127 r2132 3 3 ! 4 4 SUBROUTINE thermcell_flux(ngrid,nlay,ptimestep,masse, & 5 l alim,lmin,lmax,entr_star,detr_star,&5 lmin,lmax,entr_star,detr_star, & 6 6 f,rhobarz,zlev,zw2,fm,entr,detr,zqla) 7 7 8 8 9 9 !=============================================================================== 10 ! Purpose: deduction des flux10 ! Purpose: fluxes deduction 11 11 ! 12 12 ! Modif 2019/04 (AB alexandre.boissinot@lmd.jussieu.fr) … … 15 15 16 16 USE print_control_mod, ONLY: prt_level 17 USE thermcell_mod 17 USE thermcell_mod, ONLY: fomass_max, alpha_max 18 18 19 19 IMPLICIT NONE … … 30 30 INTEGER lmin(ngrid) 31 31 INTEGER lmax(ngrid) 32 INTEGER lalim(ngrid) 33 34 ! AB : I remove alimentation, which is in fact included in entr_star. EN COURS 35 ! REAL alim_star(ngrid,nlay) 32 36 33 REAL entr_star(ngrid,nlay) 37 34 REAL detr_star(ngrid,nlay) … … 58 55 INTEGER igout, lout ! Error grid point number and level 59 56 60 REAL zfm ! Mass flux such as updraft fraction is equal to its maximal authorized value alpha max57 REAL zfm ! Mass flux such as updraft fraction is equal to its maximal authorized value alpha_max 61 58 REAL f_old ! Save initial value of mass flux 62 59 REAL eee0 ! Save initial value of entrainment … … 65 62 REAL ddd ! ddd0 - eee 66 63 REAL zzz ! Temporary variable set to mass flux 67 REAL fact 68 REAL test 69 70 INTEGER ncorecfm1 71 INTEGER ncorecfm2 72 INTEGER ncorecfm3 73 INTEGER ncorecfm4 74 INTEGER ncorecfm5 75 INTEGER ncorecfm6 76 INTEGER ncorecalpha 77 78 LOGICAL check_debug 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 69 79 70 LOGICAL labort_physic 80 71 … … 86 77 !=============================================================================== 87 78 88 ncorecfm1 = 0 89 ncorecfm2 = 0 90 ncorecfm3 = 0 91 ncorecfm4 = 0 92 ncorecfm5 = 0 93 ncorecfm6 = 0 94 79 ncorecdetr = 0 80 ncorecentr = 0 95 81 ncorecalpha = 0 96 82 … … 100 86 101 87 !=============================================================================== 102 ! Calcul de l'entrainement, du detrainement et du flux de masse88 ! Mass flux, entrainment and detrainment 103 89 !=============================================================================== 104 90 … … 108 94 109 95 DO l=1,nlay 110 ! AB : I remove alimentation, which is in fact included in entr_star. EN COURS111 ! entr(:,l) = f(:) * (entr_star(:,l) + alim_star(:,l))112 96 entr(:,l) = f(:) * entr_star(:,l) 113 97 detr(:,l) = f(:) * detr_star(:,l) 114 98 ENDDO 115 99 116 ! AB : temporary test 117 DO l=1,nlay 118 DO ig=1,ngrid 119 IF (detr(ig,l)<0.) THEN 120 print *, 'ERROR: detr is negative in thermcell_flux!' 121 print *, 'ig,l', ig, l 122 print *, 'detr_star', detr_star(ig,l) 123 print *, 'detr', detr(ig,l) 124 ENDIF 125 ENDDO 126 ENDDO 127 128 !------------------------------------------------------------------------------- 129 ! Calcul du flux de masse 100 !------------------------------------------------------------------------------- 101 ! Mass flux and boundary conditions 130 102 !------------------------------------------------------------------------------- 131 103 … … 146 118 147 119 !=============================================================================== 148 ! Tests de validite120 ! Validity tests and corrections 149 121 !=============================================================================== 150 122 … … 152 124 153 125 !------------------------------------------------------------------------------- 154 ! Verification de la positivite du flux de masse entrant126 ! Is incoming mass flux positive ? 155 127 !------------------------------------------------------------------------------- 156 128 … … 189 161 190 162 !------------------------------------------------------------------------------- 191 ! Test entrainment positif163 ! Is entrainment positive ? 192 164 !------------------------------------------------------------------------------- 193 165 … … 216 188 print *, 'entr,detr', entr(igout,k), detr(igout,k) 217 189 ENDDO 218 print *, 'fm -', fm(igout,1)190 print *, 'fm ', fm(igout,1) 219 191 print *, '-------------------------------' 220 192 abort_message = 'Negative entrainment in thermcell_flux' … … 223 195 224 196 !------------------------------------------------------------------------------- 225 ! Test detrainment positif197 ! Is detrainment positive ? 226 198 !------------------------------------------------------------------------------- 227 199 … … 243 215 print *, 'entr,detr', entr(igout,k), detr(igout,k) 244 216 ENDDO 245 print *, 'fm -', fm(igout,1)217 print *, 'fm ', fm(igout,1) 246 218 print *, '-------------------------------' 247 219 abort_message = 'Negative detrainment in thermcell_flux' … … 250 222 251 223 !------------------------------------------------------------------------------- 252 ! Test sur fraca constante ou croissante au-dessus de lalim 253 !------------------------------------------------------------------------------- 254 255 ! AB : Do we have to decree that? If so, set iflag_thermals_optflux to 0 256 IF (iflag_thermals_optflux==0) THEN 257 DO ig=1,ngrid 258 IF (l.ge.lalim(ig).and.l.le.lmax(ig).and.(zw2(ig,l+1).gt.1.e-10).and.(zw2(ig,l).gt.1.e-10) ) THEN 259 zzz = fm(ig,l) * rhobarz(ig,l+1) * zw2(ig,l+1) & 260 & / (rhobarz(ig,l) * zw2(ig,l)) 261 262 IF (fm(ig,l+1).gt.zzz) THEN 263 detr(ig,l) = detr(ig,l) + fm(ig,l+1) - zzz 264 fm(ig,l+1) = zzz 265 ncorecfm4 = ncorecfm4 + 1 266 ENDIF 267 ENDIF 268 ENDDO 269 ENDIF 270 271 !------------------------------------------------------------------------------- 272 ! Test sur flux de masse constant ou decroissant au-dessus de lalim 273 !------------------------------------------------------------------------------- 274 275 ! AB : Do we have to decree that? If so, set iflag_thermals_optflux to 0 276 IF (iflag_thermals_optflux==0) THEN 277 DO ig=1,ngrid 278 IF ((fm(ig,l+1).gt.fm(ig,l)).and.(l.gt.lalim(ig))) then 279 f_old = fm(ig,l+1) 280 fm(ig,l+1) = fm(ig,l) 281 detr(ig,l) = detr(ig,l) + f_old - fm(ig,l+1) 282 ncorecfm5 = ncorecfm5 + 1 283 ENDIF 284 ENDDO 285 ENDIF 286 287 !------------------------------------------------------------------------------- 288 ! Test detrainment inferieur au flux de masse 224 ! Is detrainment lesser than incoming mass flux ? 289 225 !------------------------------------------------------------------------------- 290 226 … … 309 245 310 246 IF (prt_level.ge.100) THEN 247 print *, 'fm+', fm(ig,l+1) 248 print *, 'entr,detr', entr(ig,l), detr(ig,l) 311 249 print *, 'fm-', fm(ig,l) 312 print *, 'entr,detr', entr(ig,l), detr(ig,l) 313 print *, 'fm+', fm(ig,l+1) 314 ENDIF 315 316 ncorecfm6 = ncorecfm6 + 1 317 318 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 319 ! Dans le cas ou on est au dessus de la couche d'alimentation et que le 320 ! detrainement est plus fort que le flux de masse, on stope le thermique. 321 ! test : on commente 322 ! 323 ! AB : Do we have to stop the plume here? 324 ! AB : Attention, if lmax is modified, be sure that fm, zw2, entr, detr are set 325 ! to zero above this level! 326 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 327 ! IF (l.gt.lalim(ig)) THEN 328 ! lmax(ig)=l 329 ! fm(ig,l+1)=0. 330 ! entr(ig,l)=0. 331 ! ELSE 332 ! ncorecfm2=ncorecfm2+1 333 ! ENDIF 334 ENDIF 335 336 ! IF (l.gt.lmax(ig)) THEN 337 ! detr(ig,l) = 0. 338 ! fm(ig,l+1) = 0. 339 ! entr(ig,l) = 0. 340 ! ENDIF 341 ENDDO 342 343 ! labort_physic=.false. 344 ! 345 ! DO ig=1,ngrid 346 ! IF (entr(ig,l).lt.0.) THEN 347 ! labort_physic=.true. 348 ! igout=ig 349 ! ENDIF 350 ! ENDDO 351 ! 352 ! IF (labort_physic) THEN 353 ! print *, 'ERROR: detrainment is greater than mass flux and' 354 ! print *, ' entrainment cannot compensate it!' 355 ! print *, 'ig,l,lmax(ig)', igout, l, lmax(igout) 356 ! print *, '--------------------' 357 ! print *, 'fm+', fm(igout,lout+1) 358 ! print *, 'entr,detr', entr(igout,lout), detr(igout,lout) 359 ! print *, 'fm-', fm(igout,lout) 360 ! abort_message = 'problem in thermcell_flux' 361 ! CALL abort_physic (modname,abort_message,1) 362 ! ENDIF 363 364 !------------------------------------------------------------------------------- 365 ! Test fraction masse entrainee inferieure a fomass_max 250 ENDIF 251 252 ncorecdetr = ncorecdetr + 1 253 254 ENDIF 255 ENDDO 256 257 !------------------------------------------------------------------------------- 258 ! Is entrainment mass fraction lesser than fomass_max ? 366 259 !------------------------------------------------------------------------------- 367 260 … … 388 281 IF (eee.gt.0.) THEN 389 282 entr(ig,l) = entr(ig,l) - eee 390 ncorec fm3 = ncorecfm3+ 1283 ncorecentr = ncorecentr + 1 391 284 392 285 IF (prt_level.ge.10) THEN … … 395 288 ENDIF 396 289 397 IF (ddd.gt.0.) THEN ! on reduit le detr dans la couche290 IF (ddd.gt.0.) THEN ! detr in the current layer is reduced 398 291 detr(ig,l) = ddd 399 ELSEIF (l.eq.lmax(ig)) THEN ! on ajuste le detr dans la couche292 ELSEIF (l.eq.lmax(ig)) THEN ! detr in the last layer is adjusted 400 293 detr(ig,l) = fm(ig,l) + entr(ig,l) 401 ELSEIF (entr(ig,l+1).gt.ABS(ddd)) THEN ! on reduit le detr dans la couche et l'entr au-dessus294 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 reduced 402 295 detr(ig,l) = 0. 403 296 fm(ig,l+1) = fm(ig,l) + entr(ig,l) … … 423 316 ENDDO 424 317 425 IF (labort_physic) THEN 426 print *, 'ERROR: Entrainment is greater than its maximal authorized value!' 427 print *, ' Nor detrainment neither entrainment can compensate it!' 428 print *, 'ig,l,entr', igout, lout, entr(igout,lout) 429 print *, 'lmin,lmax', lmin(igout), lmax(igout) 430 print *, '--------------------' 431 print *, 'e_max :', masse(igout,lout)*fomass_max/ptimestep 432 print *, ' fomass_max :', fomass_max 433 print *, ' masse :', masse(igout,lout) 434 print *, ' ptimestep :', ptimestep 435 print *, 'norm :', f(igout) 436 ! print *, 'alim* :', alim_star(igout,lout) 437 print *, 'entr* :', entr_star(igout,lout) 438 print *, '--------------------' 439 print *, 'fm l+1 ', fm(igout,lout+2) 440 print *, 'entr <-- ', entr(igout,lout+1) 441 print *, 'detr -->', detr(igout,lout+1) 442 print *, 'fm l ', fm(igout,lout+1) 443 print *, 'entr <-- ', entr(igout,lout) 444 print *, 'detr -->', detr(igout,lout) 445 print *, 'fm l-1 ', fm(igout,lout) 446 call flush 447 abort_message = 'problem in thermcell_flux' 448 CALL abort_physic (modname,abort_message,1) 449 ENDIF 450 451 !------------------------------------------------------------------------------- 452 ! Test fraction couverte inferieure a alphamax 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 ? 453 343 !------------------------------------------------------------------------------- 454 344 … … 466 356 467 357 DO ig=1,ngrid 468 zfm = rhobarz(ig,l+1) * zw2(ig,l+1) * alpha max358 zfm = rhobarz(ig,l+1) * zw2(ig,l+1) * alpha_max 469 359 470 360 IF (fm(ig,l+1) .gt. zfm) THEN … … 495 385 496 386 ncorecalpha = ncorecalpha + 1 497 498 387 ENDIF 499 388 ENDDO … … 502 391 503 392 !------------------------------------------------------------------------------- 504 ! We check if fm, entr and detr vanish correctly in layer lmax393 ! Boundary conditions 505 394 !------------------------------------------------------------------------------- 506 395 … … 518 407 519 408 !------------------------------------------------------------------------------- 520 ! Impression des bidouilles utilisees pour corriger les problemes 521 !------------------------------------------------------------------------------- 522 523 IF (prt_level.ge.1) THEN 524 525 IF (ncorecfm4.ge.1) THEN 526 print *, 'WARNING: Deacreasing updraft fraction above lalim!' 527 print *, 'in', ncorecfm4, 'point(s)' 528 ENDIF 529 530 IF (ncorecfm5.ge.1) THEN 531 print *, 'WARNING: Increasing mass flux above lalim!' 532 print *, 'in', ncorecfm5, 'point(s)' 533 ENDIF 534 535 IF (ncorecfm6.ge.1) THEN 409 ! Corrections display 410 !------------------------------------------------------------------------------- 411 412 IF (prt_level.GE.1) THEN 413 414 IF (ncorecdetr.GE.1) THEN 536 415 print *, 'WARNING: Detrainment is greater than mass flux!' 537 print *, 'in', ncorecfm6, 'point(s)' 538 ENDIF 539 540 ! AB : those outputs are commented because corresponding test is also commented. 541 ! IF (ncorecfm2.ge.1) THEN 542 ! print *, 'WARNING: Detrainment is greater than mass flux!' 543 ! print *, 'in', ncorecfm2, 'point(s)' 544 ! ENDIF 545 546 IF (ncorecfm3.ge.1) THEN 416 print *, ' In', ncorecdetr, 'point(s)' 417 ENDIF 418 419 IF (ncorecentr.GE.1) THEN 547 420 print *, 'WARNING: Entrained mass is greater than maximal authorized value!' 548 print *, ' in', ncorecfm3, 'point(s)'549 ENDIF 550 551 IF (ncorecalpha. ge.1) THEN421 print *, ' In', ncorecentr, 'point(s)' 422 ENDIF 423 424 IF (ncorecalpha.GE.1) THEN 552 425 print *, 'WARNING: Updraft fraction is greater than maximal authorized value!' 553 print *, ' in', ncorecalpha, 'point(s)'426 print *, ' In', ncorecalpha, 'point(s)' 554 427 ENDIF 555 428 556 429 ENDIF 557 430 558 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~559 431 ! AB : temporary test added to check the validity of eq. df/dz = e - d 560 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~561 432 ! DO l=1,nlay 562 433 ! DO ig=1,ngrid
Note: See TracChangeset
for help on using the changeset viewer.
