Changeset 3851 for LMDZ6/branches/LMDZ-tracers/libf/phylmd/yamada4.F90
- Timestamp:
- Feb 22, 2021, 12:44:07 PM (4 years ago)
- Location:
- LMDZ6/branches/LMDZ-tracers
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/branches/LMDZ-tracers
- Property svn:mergeinfo changed
-
LMDZ6/branches/LMDZ-tracers/libf/phylmd/yamada4.F90
r3531 r3851 6 6 USE dimphy 7 7 USE ioipsl_getin_p_mod, ONLY : getin_p 8 8 USE phys_local_var_mod, only: tke_dissip 9 9 10 IMPLICIT NONE 10 11 include "iniprint.h" … … 56 57 ! iflag_pbl=11 -> the model starts with NP from start files created by ce0l 57 58 ! -> the model can run with longer time-steps. 58 ! 2016/11/30 (EV etienne.vignon@ univ-grenoble-alpes.fr)59 ! 2016/11/30 (EV etienne.vignon@lmd.ipsl.fr) 59 60 ! On met tke (=q2/2) en entr??e plut??t que q2 60 61 ! On corrige l'update de la tke 61 ! 62 ! 2020/10/01 (EV) 63 ! On ajoute la dissipation de la TKE en diagnostique de sortie 64 ! 62 65 ! Inpout/Output : 63 66 !============== … … 121 124 REAL,SAVE :: viscom,viscoh 122 125 !$OMP THREADPRIVATE( hboville,viscom,viscoh) 123 INTEGER ig, k126 INTEGER ig, jg, k 124 127 REAL ri, zrif, zalpha, zsm, zsn 125 128 REAL rif(klon, klev+1), sm(klon, klev+1), alpha(klon, klev) … … 186 189 viscom=1.46E-5 187 190 viscoh=2.06E-5 191 !lmixmin=1.0E-3 188 192 lmixmin=0. 189 193 yamada4_num=5 … … 416 420 ELSE IF (iflag_pbl>=10) THEN 417 421 422 shear(:,:)=0. 423 buoy(:,:)=0. 424 dissip(:,:)=0. 425 km(:,:)=0. 426 418 427 IF (yamada4_num>=1) THEN 419 428 … … 424 433 shear(ig,k)=km(ig, k)*m2(ig, k) 425 434 buoy(ig,k)=km(ig, k)*m2(ig, k)*(-1.*rif(ig,k)) 426 dissip(ig,k)=((sqrt(q2(ig,k)))**3)/(b1*l(ig,k)) 435 ! dissip(ig,k)=min(max(((sqrt(q2(ig,k)))**3)/(b1*l(ig,k)),1.E-12),1.E4) 436 dissip(ig,k)=((sqrt(q2(ig,k)))**3)/(b1*l(ig,k)) 427 437 ENDDO 428 438 ENDDO 429 439 430 440 IF (yamada4_num==1) THEN ! Schema du MAR tel quel 431 441 DO k = 2, klev - 1 … … 478 488 ENDDO 479 489 ENDDO 490 491 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 492 !! Attention, yamada4_num=5 est inexacte car néglige les termes de flottabilité 493 !! en conditions instables 494 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 480 495 ELSE IF (yamada4_num==5) THEN ! version modifiee avec integration exacte pour la dissipation 481 496 DO k = 2, klev - 1 … … 507 522 DO k = 2, klev - 1 508 523 DO ig=1,ngrid 524 !tkeprov=q2(ig,k)/ydeux 525 !tkeprov=tkeprov+max(buoy(ig,k)+shear(ig,k),0.)*dt 526 !disseff=dissip(ig,k)-min(0.,buoy(ig,k)+shear(ig,k)) 527 !tkeexp=exp(-dt*disseff/tkeprov) 528 !tkeprov= tkeprov*tkeexp 529 !q2(ig,k)=tkeprov*ydeux 530 ! En cas stable, on traite la flotabilite comme la 531 ! dissipation, en supposant que dissipeff/TKE est constant. 532 ! Puis on prend la solution exacte 533 ! 534 ! With drag and dissipation from high vegetation (EV & FC, 05/10/2020) 535 winds(ig,k)=sqrt(u(ig,k)**2+v(ig,k)**2) 509 536 tkeprov=q2(ig,k)/ydeux 510 tkeprov=tkeprov+max(buoy(ig,k)+shear(ig,k) ,0.)*dt511 disseff=dissip(ig,k)-min(0.,buoy(ig,k)+shear(ig,k) )537 tkeprov=tkeprov+max(buoy(ig,k)+shear(ig,k)+drgpro(ig,k)*(winds(ig,k))**3,0.)*dt 538 disseff=dissip(ig,k)-min(0.,buoy(ig,k)+shear(ig,k)+drgpro(ig,k)*(winds(ig,k))**3) + drgpro(ig,k)*tkeprov 512 539 tkeexp=exp(-dt*disseff/tkeprov) 513 540 tkeprov= tkeprov*tkeexp 514 541 q2(ig,k)=tkeprov*ydeux 515 ! En cas stable, on traite la flotabilite comme la 516 ! dissipation, en supposant que buoy/q2^3 est constant. 517 ! Puis on prend la solution exacte 542 518 543 ENDDO 519 544 ENDDO … … 725 750 726 751 !============================================================================ 752 ! Diagnostique de la dissipation 753 !============================================================================ 754 755 ! Diagnostics 756 tke_dissip(1:ngrid,:,nsrf)=0. 757 ! DO k=2,klev 758 ! DO ig=1,ngrid 759 ! jg=ni(ig) 760 ! tke_dissip(jg,k,nsrf)=dissip(ig,k) 761 ! ENDDO 762 ! ENDDO 763 764 !============================================================================= 727 765 728 766 RETURN … … 1017 1055 !===================================================================== 1018 1056 1019 1057 l1(1:ngrid,:)=0. 1020 1058 IF (iflag_pbl==8 .OR. iflag_pbl==10) THEN 1021 1059 … … 1135 1173 1136 1174 1137 DO k= 2,klev1175 DO k=1,klev+1 1138 1176 DO ig=1,ngrid 1139 1177 lmix(ig,k)=MAX(MAX(l1(ig,k), l2(ig,k)),lmixmin)
Note: See TracChangeset
for help on using the changeset viewer.