Changeset 2127 for trunk/LMDZ.GENERIC/libf/phystd/thermcell_flux.F90
- Timestamp:
- Apr 29, 2019, 10:07:47 AM (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.GENERIC/libf/phystd/thermcell_flux.F90
r2103 r2127 3 3 ! 4 4 SUBROUTINE thermcell_flux(ngrid,nlay,ptimestep,masse, & 5 lalim,lmin,lmax,alim_star,entr_star,detr_star, & 6 f,rhobarz,zlev,zw2,fm,entr, & 7 detr,zqla,lev_out,lunout1,igout) 8 9 10 !============================================================================== 11 ! thermcell_flux: deduction des flux 12 !============================================================================== 5 lalim,lmin,lmax,entr_star,detr_star, & 6 f,rhobarz,zlev,zw2,fm,entr,detr,zqla) 7 8 9 !=============================================================================== 10 ! Purpose: deduction des flux 11 ! 12 ! Modif 2019/04 (AB alexandre.boissinot@lmd.jussieu.fr) 13 ! 14 !=============================================================================== 13 15 14 16 USE print_control_mod, ONLY: prt_level … … 18 20 19 21 20 !============================================================================== 22 !=============================================================================== 21 23 ! Declaration 22 !============================================================================== 23 24 ! inputs: 25 ! ------- 26 27 INTEGER ngrid 28 INTEGER nlay 29 INTEGER igout 30 INTEGER lev_out 31 INTEGER lunout1 24 !=============================================================================== 25 26 ! Inputs: 27 ! ------- 28 29 INTEGER ngrid, nlay 32 30 INTEGER lmin(ngrid) 33 31 INTEGER lmax(ngrid) 34 32 INTEGER lalim(ngrid) 35 33 36 REAL alim_star(ngrid,nlay) 34 ! AB : I remove alimentation, which is in fact included in entr_star. EN COURS 35 ! REAL alim_star(ngrid,nlay) 37 36 REAL entr_star(ngrid,nlay) 38 37 REAL detr_star(ngrid,nlay) … … 46 45 REAL zmax(ngrid) 47 46 48 ! outputs:49 ! 47 ! Outputs: 48 ! -------- 50 49 51 50 REAL entr(ngrid,nlay) … … 53 52 REAL fm(ngrid,nlay+1) 54 53 55 ! local:56 ! 57 58 INTEGER ig, l,k59 INTEGER lout60 61 REAL zfm ! mass flux such as updraft fraction is equal to its maximal authorized value alphamax62 REAL f_old ! save initial value of mass flux63 REAL eee0 ! save initial value of entrainment64 REAL ddd0 ! save initial value of detrainment54 ! Local: 55 ! ------ 56 57 INTEGER ig, l, k 58 INTEGER igout, lout ! Error grid point number and level 59 60 REAL zfm ! Mass flux such as updraft fraction is equal to its maximal authorized value alphamax 61 REAL f_old ! Save initial value of mass flux 62 REAL eee0 ! Save initial value of entrainment 63 REAL ddd0 ! Save initial value of detrainment 65 64 REAL eee ! eee0 - layer mass * maximal authorized mass fraction / dt 66 65 REAL ddd ! ddd0 - eee 67 REAL zzz ! temporary variable set to mass flux66 REAL zzz ! Temporary variable set to mass flux 68 67 REAL fact 69 68 REAL test … … 83 82 CHARACTER (len=80) :: abort_message 84 83 85 !============================================================================== 84 !=============================================================================== 86 85 ! Initialization 87 !============================================================================== 86 !=============================================================================== 88 87 89 88 ncorecfm1 = 0 … … 100 99 fm(:,:) = 0. 101 100 102 IF (prt_level.ge.10) THEN 103 write(lunout1,*) 'Dans thermcell_flux 0' 104 write(lunout1,*) 'flux base ', f(igout) 105 write(lunout1,*) 'lmax ', lmax(igout) 106 write(lunout1,*) 'lalim ', lalim(igout) 107 write(lunout1,*) 'ig= ', igout 108 write(lunout1,*) ' l E* A* D* ' 109 write(lunout1,'(i4,3e15.5)') (l,entr_star(igout,l), & 110 & alim_star(igout,l),detr_star(igout,l), & 111 & l=1,lmax(igout)) 112 ENDIF 113 114 !============================================================================== 101 !=============================================================================== 115 102 ! Calcul de l'entrainement, du detrainement et du flux de masse 116 !============================================================================== 117 118 !------------------------------------------------------------------------------ 103 !=============================================================================== 104 105 !------------------------------------------------------------------------------- 119 106 ! Multiplication par la norme issue de la relation de fermeture 120 !------------------------------------------------------------------------------ 107 !------------------------------------------------------------------------------- 121 108 122 109 DO l=1,nlay 123 entr(:,l) = f(:) * (entr_star(:,l) + alim_star(:,l)) 110 ! AB : I remove alimentation, which is in fact included in entr_star. EN COURS 111 ! entr(:,l) = f(:) * (entr_star(:,l) + alim_star(:,l)) 112 entr(:,l) = f(:) * entr_star(:,l) 124 113 detr(:,l) = f(:) * detr_star(:,l) 125 114 ENDDO 126 115 127 !------------------------------------------------------------------------------ 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 !------------------------------------------------------------------------------- 128 129 ! Calcul du flux de masse 129 !------------------------------------------------------------------------------ 130 !------------------------------------------------------------------------------- 130 131 131 132 DO l=1,nlay … … 144 145 ENDDO 145 146 146 !============================================================================== 147 !=============================================================================== 147 148 ! Tests de validite 148 !============================================================================== 149 !=============================================================================== 149 150 150 151 DO l=1,nlay 151 152 152 !------------------------------------------------------------------------------ 153 !------------------------------------------------------------------------------- 153 154 ! Verification de la positivite du flux de masse entrant 154 !------------------------------------------------------------------------------ 155 156 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 155 !------------------------------------------------------------------------------- 156 157 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 157 158 ! AB : I remove the correction and replace it by an uncompromising test. 158 159 ! According to the previous derivations, we MUST have positive mass flux … … 160 161 ! Then the only value which can be negative is the mass flux at the top of 161 162 ! the plume. However, this value was set to zero a few lines above. 162 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 163 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 163 164 164 165 labort_physic=.false. … … 177 178 print *, 'lmin,lmax', lmin(igout), lmax(igout) 178 179 print *, '-------------------------------' 179 print *, 'fm+', fm(igout,lout+1) 180 print *, 'entr,detr', entr(igout,lout), detr(igout,lout) 181 print *, 'fm ', fm(igout,lout) 182 print *, 'entr,detr', entr(igout,lout-1), detr(igout,lout-1) 183 print *, 'fm-', fm(igout,lout-1) 180 DO k=nlay,1,-1 181 print *, 'fm ', fm(igout,k+1) 182 print *, 'entr,detr', entr(igout,k), detr(igout,k) 183 ENDDO 184 print *, 'fm-', fm(igout,1) 185 print *, '-------------------------------' 184 186 abort_message = 'Negative incoming fm in thermcell_flux' 185 187 CALL abort_physic(modname,abort_message,1) 186 188 ENDIF 187 189 188 !------------------------------------------------------------------------------ 190 !------------------------------------------------------------------------------- 189 191 ! Test entrainment positif 190 !------------------------------------------------------------------------------ 191 192 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 192 !------------------------------------------------------------------------------- 193 194 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 193 195 ! AB : According to the previous derivations, we MUST have positive entrainment 194 196 ! and detrainment everywhere! 195 197 ! Indeed, they are set to max between zero and a computed value. 196 198 ! Then tests are uncompromising. 197 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 198 199 labort_physic=.false. 199 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 200 200 201 201 DO ig=1,ngrid … … 212 212 print *, 'lmin,lmax', lmin(igout), lmax(igout) 213 213 print *, '-------------------------------' 214 print *, 'fm+', fm(igout,lout+1) 215 print *, 'entr,detr', entr(igout,lout), detr(igout,lout) 216 print *, 'fm ', fm(igout,lout) 217 print *, 'entr,detr', entr(igout,lout-1), detr(igout,lout-1) 218 print *, 'fm-', fm(igout,lout-1) 214 DO k=nlay,1,-1 215 print *, 'fm ', fm(igout,k+1) 216 print *, 'entr,detr', entr(igout,k), detr(igout,k) 217 ENDDO 218 print *, 'fm-', fm(igout,1) 219 print *, '-------------------------------' 219 220 abort_message = 'Negative entrainment in thermcell_flux' 220 221 CALL abort_physic(modname,abort_message,1) 221 222 ENDIF 222 223 223 !------------------------------------------------------------------------------ 224 !------------------------------------------------------------------------------- 224 225 ! Test detrainment positif 225 !------------------------------------------------------------------------------ 226 !------------------------------------------------------------------------------- 226 227 227 228 DO ig=1,ngrid … … 238 239 print *, 'lmin,lmax', lmin(igout), lmax(igout) 239 240 print *, '-------------------------------' 240 print *, 'fm+', fm(igout,lout+1) 241 print *, 'entr,detr', entr(igout,lout), detr(igout,lout) 242 print *, 'fm ', fm(igout,lout) 243 print *, 'entr,detr', entr(igout,lout-1), detr(igout,lout-1) 244 print *, 'fm-', fm(igout,lout-1) 241 DO k=nlay,1,-1 242 print *, 'fm ', fm(igout,k+1) 243 print *, 'entr,detr', entr(igout,k), detr(igout,k) 244 ENDDO 245 print *, 'fm-', fm(igout,1) 246 print *, '-------------------------------' 245 247 abort_message = 'Negative detrainment in thermcell_flux' 246 248 CALL abort_physic(modname,abort_message,1) 247 249 ENDIF 248 250 249 !------------------------------------------------------------------------------ 251 !------------------------------------------------------------------------------- 250 252 ! Test sur fraca constante ou croissante au-dessus de lalim 251 !------------------------------------------------------------------------------ 253 !------------------------------------------------------------------------------- 252 254 253 255 ! AB : Do we have to decree that? If so, set iflag_thermals_optflux to 0 … … 267 269 ENDIF 268 270 269 !------------------------------------------------------------------------------ 271 !------------------------------------------------------------------------------- 270 272 ! Test sur flux de masse constant ou decroissant au-dessus de lalim 271 !------------------------------------------------------------------------------ 273 !------------------------------------------------------------------------------- 272 274 273 275 ! AB : Do we have to decree that? If so, set iflag_thermals_optflux to 0 … … 283 285 ENDIF 284 286 285 !------------------------------------------------------------------------------ 287 !------------------------------------------------------------------------------- 286 288 ! Test detrainment inferieur au flux de masse 287 !------------------------------------------------------------------------------ 288 289 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 289 !------------------------------------------------------------------------------- 290 291 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 290 292 ! AB : Even if fm has no negative value, it can be lesser than detr. 291 293 ! That's not suitable because when we will mix the plume with the … … 294 296 ! otherwise we get a negative outgoing mass flux (cf. above). 295 297 ! That is why we correct the detrainment as follows. 296 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 298 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 297 299 298 300 DO ig=1,ngrid … … 314 316 ncorecfm6 = ncorecfm6 + 1 315 317 316 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 318 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 317 319 ! Dans le cas ou on est au dessus de la couche d'alimentation et que le 318 320 ! detrainement est plus fort que le flux de masse, on stope le thermique. … … 320 322 ! 321 323 ! AB : Do we have to stop the plume here? 322 ! AB : Attention, if lmax is modified, be sure that fm, zw2, entr and detr are323 ! setto zero above this level!324 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 324 ! AB : Attention, if lmax is modified, be sure that fm, zw2, entr, detr are set 325 ! to zero above this level! 326 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 325 327 ! IF (l.gt.lalim(ig)) THEN 326 328 ! lmax(ig)=l … … 331 333 ! ENDIF 332 334 ENDIF 333 !335 334 336 ! IF (l.gt.lmax(ig)) THEN 335 337 ! detr(ig,l) = 0. … … 360 362 ! ENDIF 361 363 362 !------------------------------------------------------------------------------ 364 !------------------------------------------------------------------------------- 363 365 ! Test fraction masse entrainee inferieure a fomass_max 364 !------------------------------------------------------------------------------ 365 366 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 366 !------------------------------------------------------------------------------- 367 368 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 367 369 ! AB : Entrainment is bigger than the maximal authorized value. 368 370 ! If we consider that the excess entrainement is in fact plume air which … … 374 376 ! detrainment and mass flux profiles such as we do not exceed the maximal 375 377 ! authorized entrained mass. 376 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 378 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 377 379 378 380 labort_physic=.false. … … 432 434 print *, ' ptimestep :', ptimestep 433 435 print *, 'norm :', f(igout) 434 print *, 'alim* :', alim_star(igout,lout)436 ! print *, 'alim* :', alim_star(igout,lout) 435 437 print *, 'entr* :', entr_star(igout,lout) 436 438 print *, '--------------------' … … 447 449 ENDIF 448 450 449 !------------------------------------------------------------------------------ 451 !------------------------------------------------------------------------------- 450 452 ! Test fraction couverte inferieure a alphamax 451 !------------------------------------------------------------------------------ 452 453 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 453 !------------------------------------------------------------------------------- 454 455 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 454 456 ! FH : Partie a revisiter. 455 457 ! Il semble qu'etaient codees ici deux optiques dans le cas F/(rho*w) > 1 … … 461 463 ! semble de toutes facons deraisonable d'avoir des fractions de 1... 462 464 ! Ci-dessous, et dans l'etat actuel, le premier des deux if est sans doute inutile. 463 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 465 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 464 466 465 467 DO ig=1,ngrid 466 468 zfm = rhobarz(ig,l+1) * zw2(ig,l+1) * alphamax 469 467 470 IF (fm(ig,l+1) .gt. zfm) THEN 468 471 f_old = fm(ig,l+1) … … 470 473 detr(ig,l) = detr(ig,l) + f_old - fm(ig,l+1) 471 474 472 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 475 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 473 476 ! AB : The previous change doesn't observe the equation df/dz = e - d. 474 477 ! To avoid this issue, what is better to do? I choose to increase the 475 478 ! entrainment in the layer above when l<lmax. If l=lmax then fm(l+1)=0 and 476 479 ! there are never any problems. 477 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 480 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 478 481 IF (l.lt.lmax(ig)) THEN 479 482 entr(ig,l+1) = entr(ig,l+1) + f_old - fm(ig,l+1) … … 498 501 ENDDO 499 502 500 !------------------------------------------------------------------------------ 503 !------------------------------------------------------------------------------- 501 504 ! We check if fm, entr and detr vanish correctly in layer lmax 502 !------------------------------------------------------------------------------ 503 504 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 505 !------------------------------------------------------------------------------- 506 507 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 505 508 ! AB : test added to avoid problem when lmax = 0, which is not a fortran index. 506 509 ! Theoretically, we always get lmax >= lmin >= linf > 0 507 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 510 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 508 511 DO ig=1,ngrid 509 512 IF (lmax(ig).gt.0) THEN … … 514 517 ENDDO 515 518 516 !------------------------------------------------------------------------------ 519 !------------------------------------------------------------------------------- 517 520 ! Impression des bidouilles utilisees pour corriger les problemes 518 !------------------------------------------------------------------------------ 521 !------------------------------------------------------------------------------- 519 522 520 523 IF (prt_level.ge.1) THEN … … 553 556 ENDIF 554 557 555 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 558 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 556 559 ! AB : temporary test added to check the validity of eq. df/dz = e - d 557 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 560 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 558 561 ! DO l=1,nlay 559 562 ! DO ig=1,ngrid
Note: See TracChangeset
for help on using the changeset viewer.
