Changeset 2103 for trunk/LMDZ.GENERIC/libf
- Timestamp:
- Feb 18, 2019, 11:57:35 AM (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.GENERIC/libf/phystd/thermcell_flux.F90
r2101 r2103 56 56 ! ------ 57 57 58 INTEGER ig,l 58 INTEGER ig,l,k 59 59 INTEGER lout 60 60 … … 66 66 REAL ddd ! ddd0 - eee 67 67 REAL zzz ! temporary variable set to mass flux 68 68 REAL fact 69 69 REAL test 70 70 … … 75 75 INTEGER ncorecfm5 76 76 INTEGER ncorecfm6 77 INTEGER ncorecfm778 INTEGER ncorecfm879 77 INTEGER ncorecalpha 80 78 … … 95 93 ncorecfm5 = 0 96 94 ncorecfm6 = 0 97 ncorecfm7 = 098 ncorecfm8 = 099 95 100 96 ncorecalpha = 0 … … 178 174 IF (labort_physic) THEN 179 175 print *, 'ERROR: mass flux has negative value(s)!' 180 print *, 'ig,l,fm',ig,l,entr(igout,lout) 181 abort_message = 'negative incoming fm in thermcell_flux' 176 print *, 'ig,l,fm',igout, lout, fm(igout,lout) 177 print *, 'lmin,lmax', lmin(igout), lmax(igout) 178 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) 184 abort_message = 'Negative incoming fm in thermcell_flux' 185 CALL abort_physic(modname,abort_message,1) 186 ENDIF 187 188 !------------------------------------------------------------------------------ 189 ! Test entrainment positif 190 !------------------------------------------------------------------------------ 191 192 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 193 ! AB : According to the previous derivations, we MUST have positive entrainment 194 ! and detrainment everywhere! 195 ! Indeed, they are set to max between zero and a computed value. 196 ! Then tests are uncompromising. 197 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 198 199 labort_physic=.false. 200 201 DO ig=1,ngrid 202 IF (entr(ig,l).lt.0.) THEN 203 labort_physic = .true. 204 igout = ig 205 lout = l 206 ENDIF 207 ENDDO 208 209 IF (labort_physic) THEN 210 print *, 'ERROR: entrainment has negative value(s)!' 211 print *, 'ig,l,entr', igout, lout, entr(igout,lout) 212 print *, 'lmin,lmax', lmin(igout), lmax(igout) 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) 219 abort_message = 'Negative entrainment in thermcell_flux' 220 CALL abort_physic(modname,abort_message,1) 221 ENDIF 222 223 !------------------------------------------------------------------------------ 224 ! Test detrainment positif 225 !------------------------------------------------------------------------------ 226 227 DO ig=1,ngrid 228 IF (detr(ig,l).lt.0.) THEN 229 labort_physic = .true. 230 igout = ig 231 lout = l 232 ENDIF 233 ENDDO 234 235 IF (labort_physic) THEN 236 print *, 'ERROR: detrainment has negative value(s)!' 237 print *, 'ig,l,detr', igout, lout, detr(igout,lout) 238 print *, 'lmin,lmax', lmin(igout), lmax(igout) 239 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) 245 abort_message = 'Negative detrainment in thermcell_flux' 182 246 CALL abort_physic(modname,abort_message,1) 183 247 ENDIF … … 220 284 221 285 !------------------------------------------------------------------------------ 222 ! Test entrainment et detrainment positif223 !------------------------------------------------------------------------------224 225 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~226 ! AB : According to the previous derivations, we MUST have positive entrainment227 ! and detrainment everywhere!228 ! Indeed, they are set to max between zero and a computed value.229 ! Then tests are uncompromising.230 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~231 232 labort_physic=.false.233 234 DO ig=1,ngrid235 IF (entr(ig,l).lt.0.) THEN236 labort_physic = .true.237 igout = ig238 lout = l239 ENDIF240 ENDDO241 242 IF (labort_physic) THEN243 print *, 'ERROR: entrainment has negative value(s)!'244 print *, 'ig,l,entr', igout, lout, entr(igout,lout)245 abort_message = 'negative entrainment in thermcell_flux'246 CALL abort_physic(modname,abort_message,1)247 ENDIF248 249 DO ig=1,ngrid250 IF (detr(ig,l).lt.0.) THEN251 labort_physic = .true.252 igout = ig253 lout = l254 ENDIF255 ENDDO256 257 IF (labort_physic) THEN258 print *, 'ERROR: detrainment has negative value(s)!'259 print *, 'ig,l,detr', igout, lout, detr(igout,lout)260 abort_message = 'negative detrainment in thermcell_flux'261 CALL abort_physic(modname,abort_message,1)262 ENDIF263 264 !------------------------------------------------------------------------------265 286 ! Test detrainment inferieur au flux de masse 266 287 !------------------------------------------------------------------------------ … … 281 302 282 303 IF (prt_level.ge.10) THEN 283 print *, ' Detrainment is modified in thermcell_flux'304 print *, 'WARNING: Detrainment is modified in thermcell_flux!' 284 305 print *, 'ig,l,lmax', ig, l, lmax(ig) 285 306 ENDIF … … 299 320 ! 300 321 ! AB : Do we have to stop the plume here? 301 ! 302 ! AB : WARNING: if lmax is modified, be sure that fm, zw2, entr and detr are 322 ! AB : Attention, if lmax is modified, be sure that fm, zw2, entr and detr are 303 323 ! set to zero above this level! 304 324 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ … … 308 328 ! entr(ig,l)=0. 309 329 ! ELSE 310 ! ncorecfm 7=ncorecfm7+1330 ! ncorecfm2=ncorecfm2+1 311 331 ! ENDIF 312 332 ENDIF … … 332 352 ! print *, ' entrainment cannot compensate it!' 333 353 ! print *, 'ig,l,lmax(ig)', igout, l, lmax(igout) 334 ! print *, 'entr(ig,l)', entr(igout,l) 335 ! print *, 'fm(ig,l)', fm(igout,l) 354 ! print *, '--------------------' 355 ! print *, 'fm+', fm(igout,lout+1) 356 ! print *, 'entr,detr', entr(igout,lout), detr(igout,lout) 357 ! print *, 'fm-', fm(igout,lout) 336 358 ! abort_message = 'problem in thermcell_flux' 337 359 ! CALL abort_physic (modname,abort_message,1) … … 348 370 ! If it's not enough, we can increase entr in the layer above and decrease 349 371 ! the outgoing mass flux in the current layer. 350 ! If it's still unsufficient, code will abort. 372 ! If it's still unsufficient, code will abort (now commented). 373 ! Else we reset entr to its intial value and we renormalize entrainment, 374 ! detrainment and mass flux profiles such as we do not exceed the maximal 375 ! authorized entrained mass. 351 376 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 352 377 … … 363 388 ncorecfm3 = ncorecfm3 + 1 364 389 365 print *, 'CORRECTION: entr is to large!' 366 print *, 'ig,l,lmax', ig, l, lmax(ig) 367 368 IF (ddd.gt.0.) THEN 390 IF (prt_level.ge.10) THEN 391 print *, 'WARNING: Entrainment is modified in thermcell_flux!' 392 print *, 'ig,l,lmax', ig, l, lmax(ig) 393 ENDIF 394 395 IF (ddd.gt.0.) THEN ! on reduit le detr dans la couche 369 396 detr(ig,l) = ddd 397 ELSEIF (l.eq.lmax(ig)) THEN ! on ajuste le detr dans la couche 398 detr(ig,l) = fm(ig,l) + entr(ig,l) 399 ELSEIF (entr(ig,l+1).gt.ABS(ddd)) THEN ! on reduit le detr dans la couche et l'entr au-dessus 400 detr(ig,l) = 0. 401 fm(ig,l+1) = fm(ig,l) + entr(ig,l) 402 entr(ig,l+1) = entr(ig,l+1) + ddd 403 ! ELSE 404 ! entr(ig,l) = entr(ig,l) + eee 405 ! igout=ig 406 ! lout=l 407 ! labort_physic=.true. 370 408 ELSE 371 IF (l.eq.lmax(ig)) THEN 372 detr(ig,l) = fm(ig,l) + entr(ig,l) 373 ELSEIF (l.lt.lmax(ig)) THEN 374 entr(ig,l+1) = entr(ig,l+1) - ddd 375 detr(ig,l) = 0. 376 fm(ig,l+1) = fm(ig,l) + entr(ig,l) 377 ELSE 378 igout=ig 379 lout=l 380 labort_physic=.true. 409 fact = (eee0 - eee) / eee0 410 entr(ig,l) = eee0 411 detr(ig,:) = detr(ig,:) * fact 412 entr(ig,:) = entr(ig,:) * fact 413 fm(ig,:) = fm(ig,:) * fact 414 415 IF (prt_level.ge.1) THEN 416 print *, 'WARNING: Normalisation is modified in thermcell_flux!' 417 print *, 'old f, new f :', f(ig), fact * f(ig) 381 418 ENDIF 382 419 ENDIF … … 385 422 386 423 IF (labort_physic) THEN 387 print *, ' 388 print *, ' 389 print *, 'ig,l, lmax', igout, lout, lmax(igout)390 print *, ' entr :', eee0391 print *, ' detr :', ddd0424 print *, 'ERROR: Entrainment is greater than its maximal authorized value!' 425 print *, ' Nor detrainment neither entrainment can compensate it!' 426 print *, 'ig,l,entr', igout, lout, entr(igout,lout) 427 print *, 'lmin,lmax', lmin(igout), lmax(igout) 428 print *, '--------------------' 392 429 print *, 'e_max :', masse(igout,lout)*fomass_max/ptimestep 393 print *, ' --- fomass_max :', fomass_max 394 print *, ' --- masse :', masse(igout,lout) 395 print *, ' --- ptimestep :', ptimestep 396 print *, 'eee :', eee 397 print *, 'ddd :', ddd 398 print *, 'fm-', fm(ig,l) 399 print *, 'entr,detr', entr(ig,l), detr(ig,l) 400 print *, 'fm+', fm(ig,l+1) 430 print *, ' fomass_max :', fomass_max 431 print *, ' masse :', masse(igout,lout) 432 print *, ' ptimestep :', ptimestep 433 print *, 'norm :', f(igout) 434 print *, 'alim* :', alim_star(igout,lout) 435 print *, 'entr* :', entr_star(igout,lout) 436 print *, '--------------------' 437 print *, 'fm l+1 ', fm(igout,lout+2) 438 print *, 'entr <-- ', entr(igout,lout+1) 439 print *, 'detr -->', detr(igout,lout+1) 440 print *, 'fm l ', fm(igout,lout+1) 441 print *, 'entr <-- ', entr(igout,lout) 442 print *, 'detr -->', detr(igout,lout) 443 print *, 'fm l-1 ', fm(igout,lout) 444 call flush 401 445 abort_message = 'problem in thermcell_flux' 402 446 CALL abort_physic (modname,abort_message,1) … … 429 473 ! AB : The previous change doesn't observe the equation df/dz = e - d. 430 474 ! To avoid this issue, what is better to do? I choose to increase the 431 ! entrainment in the layer above. 432 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 433 entr(ig,l+1) = entr(ig,l+1) + f_old - fm(ig,l+1) 475 ! entrainment in the layer above when l<lmax. If l=lmax then fm(l+1)=0 and 476 ! there are never any problems. 477 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 478 IF (l.lt.lmax(ig)) THEN 479 entr(ig,l+1) = entr(ig,l+1) + f_old - fm(ig,l+1) 480 ENDIF 434 481 435 482 IF (prt_level.ge.10) THEN … … 437 484 print *, 'ig,l,lmax', ig, l, lmax(ig) 438 485 ENDIF 439 486 440 487 IF (prt_level.ge.100) THEN 441 488 print *, 'fm-', fm(ig,l) … … 444 491 ENDIF 445 492 446 ncorecalpha = ncorecalpha+1 447 448 ENDIF 449 ENDDO 450 451 !------------------------------------------------------------------------------ 452 ! Test sur fm sortant positif 453 !------------------------------------------------------------------------------ 454 455 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 456 ! AB : Is it usefull to do that? 457 ! We will check fm in the first test with the next value of l 458 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 459 ! DO ig=1,ngrid 460 ! IF (fm(ig,l+1).lt.0.) THEN 461 ! detr(ig,l) = detr(ig,l) + fm(ig,l+1) 462 ! fm(ig,l+1) = 0. 463 ! ncorecfm2 = ncorecfm2 + 1 464 ! ENDIF 465 ! ENDDO 466 ! 467 ! labort_physic=.false. 468 ! 469 ! DO ig=1,ngrid 470 ! IF (detr(ig,l).lt.0.) THEN 471 ! labort_physic = .true. 472 ! igout = ig 473 ! lout = l 474 ! ENDIF 475 ! ENDDO 476 ! 477 ! IF (labort_physic) THEN 478 ! print *, ' ERROR: mass flux has negative value(s) and detrainement' 479 ! print *, ' cannot compensate it!' 480 ! print*,'ig,l,entr', igout, lout, entr(igout,lout) 481 ! abort_message = 'negative outgoing fm in thermcell flux' 482 ! CALL abort_physic (modname,abort_message,1) 483 ! ENDIF 493 ncorecalpha = ncorecalpha + 1 494 495 ENDIF 496 ENDDO 484 497 485 498 ENDDO … … 491 504 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 492 505 ! AB : test added to avoid problem when lmax = 0, which is not a fortran index. 493 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 494 IF (lmax(ig).gt.0) THEN 495 DO ig=1,ngrid 506 ! Theoretically, we always get lmax >= lmin >= linf > 0 507 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 508 DO ig=1,ngrid 509 IF (lmax(ig).gt.0) THEN 496 510 detr(ig,lmax(ig)) = fm(ig,lmax(ig)) + entr(ig,lmax(ig)) 497 511 fm(ig,lmax(ig)+1) = 0. 498 512 entr(ig,lmax(ig)) = 0. 499 END DO500 END IF513 ENDIF 514 ENDDO 501 515 502 516 !------------------------------------------------------------------------------ … … 505 519 506 520 IF (prt_level.ge.1) THEN 507 508 ! AB : those outputs are commented for their uselessness.509 ! IF (ncorecfm2.ge.1) THEN510 ! print *, 'WARNING: Outgoing mass flux has negative value(s)!'511 ! print *, 'in', ncorecfm2, 'point(s)'512 ! ENDIF513 521 514 522 IF (ncorecfm4.ge.1) THEN … … 528 536 529 537 ! AB : those outputs are commented because corresponding test is also commented. 530 ! IF (ncorecfm 7.ge.1) THEN538 ! IF (ncorecfm2.ge.1) THEN 531 539 ! print *, 'WARNING: Detrainment is greater than mass flux!' 532 ! print *, 'in', ncorecfm 7, 'point(s)'540 ! print *, 'in', ncorecfm2, 'point(s)' 533 541 ! ENDIF 534 542 535 543 IF (ncorecfm3.ge.1) THEN 536 print *, 'WARNING: Entrained mass is greater than itsmaximal authorized value!'544 print *, 'WARNING: Entrained mass is greater than maximal authorized value!' 537 545 print *, 'in', ncorecfm3, 'point(s)' 538 546 ENDIF 539 547 540 548 IF (ncorecalpha.ge.1) THEN 541 print *, 'WARNING: Updraft fraction is greater than itsmaximal authorized value!'549 print *, 'WARNING: Updraft fraction is greater than maximal authorized value!' 542 550 print *, 'in', ncorecalpha, 'point(s)' 543 551 ENDIF … … 554 562 ! print *, 'WARNING: df/dz != entr - detr' 555 563 ! print *, 'ig,l', ig, l 556 ! print *, 'fm - :', fm(ig,l)564 ! print *, 'fm+ :', fm(ig,l+1) 557 565 ! print *, 'entr,detr', entr(ig,l), detr(ig,l) 558 ! print *, 'fm+ :', fm(ig,l+1) 566 ! print *, 'fm :', fm(ig,l) 567 ! print *, 'entr,detr', entr(ig,l-1), detr(ig,l-1) 568 ! print *, 'fm- :', fm(ig,l-1) 559 569 ! print *, 'err. :', test 560 570 ! ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
