Changeset 852
- Timestamp:
- Oct 11, 2007, 3:43:42 PM (17 years ago)
- Location:
- LMDZ4/trunk/libf/phytherm
- Files:
-
- 17 edited
-
ajsec.F (modified) (1 diff)
-
clouds_gno.F (modified) (1 diff)
-
coef_diff_turb_mod.F90 (modified) (8 diffs)
-
ini_histday.h (modified) (2 diffs)
-
ini_histmth.h (modified) (7 diffs)
-
pbl_surface_mod.F90 (modified) (10 diffs)
-
phyetat0.F (modified) (7 diffs)
-
phyredem.F (modified) (6 diffs)
-
physiq.F (modified) (15 diffs)
-
thermcell_closure.F90 (modified) (3 diffs)
-
thermcell_dry.F90 (modified) (7 diffs)
-
thermcell_flux.F90 (modified) (20 diffs)
-
thermcell_main.F90 (modified) (20 diffs)
-
thermcell_plume.F90 (modified) (5 diffs)
-
write_histday.h (modified) (1 diff)
-
write_histmth.h (modified) (5 diffs)
-
yamada4.F (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ4/trunk/libf/phytherm/ajsec.F
r814 r852 2 2 ! $Header$ 3 3 ! 4 SUBROUTINE ajsec(paprs, pplay, t,q, d_t,d_q) 4 SUBROUTINE ajsec(paprs, pplay, t,q,limbas,d_t,d_q) 5 USE dimphy 6 IMPLICIT none 7 c====================================================================== 8 c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818 9 c Objet: ajustement sec (adaptation du GCM du LMD) 10 c====================================================================== 11 c Arguments: 12 c t-------input-R- Temperature 13 c 14 c d_t-----output-R-Incrementation de la temperature 15 c====================================================================== 16 cym#include "dimensions.h" 17 cym#include "dimphy.h" 18 #include "YOMCST.h" 19 REAL paprs(klon,klev+1), pplay(klon,klev) 20 REAL t(klon,klev), q(klon,klev) 21 REAL d_t(klon,klev), d_q(klon,klev) 22 c 23 INTEGER limbas(klon), limhau ! les couches a ajuster 24 c 25 LOGICAL mixq 26 ccc PARAMETER (mixq=.TRUE.) 27 PARAMETER (mixq=.FALSE.) 28 c 29 REAL zh(klon,klev) 30 REAL zho(klon,klev) 31 REAL zq(klon,klev) 32 REAL zpk(klon,klev) 33 REAL zpkdp(klon,klev) 34 REAL hm, sm, qm 35 LOGICAL modif(klon), down 36 INTEGER i, k, k1, k2 37 c 38 c Initialisation: 39 c 40 cym 41 limhau=klev 42 43 DO k = 1, klev 44 DO i = 1, klon 45 d_t(i,k) = 0.0 46 d_q(i,k) = 0.0 47 ENDDO 48 ENDDO 49 c------------------------------------- detection des profils a modifier 50 DO k = 1, limhau 51 DO i = 1, klon 52 zpk(i,k) = pplay(i,k)**RKAPPA 53 zh(i,k) = RCPD * t(i,k)/ zpk(i,k) 54 zho(i,k) = zh(i,k) 55 zq(i,k) = q(i,k) 56 ENDDO 57 ENDDO 58 c 59 DO k = 1, limhau 60 DO i = 1, klon 61 zpkdp(i,k) = zpk(i,k) * (paprs(i,k)-paprs(i,k+1)) 62 ENDDO 63 ENDDO 64 c 65 DO i = 1, klon 66 modif(i) = .FALSE. 67 ENDDO 68 DO k = 2, limhau 69 DO i = 1, klon 70 IF (.NOT.modif(i).and.k-1>limbas(i)) THEN 71 IF ( zh(i,k).LT.zh(i,k-1) ) modif(i) = .TRUE. 72 ENDIF 73 ENDDO 74 ENDDO 75 c------------------------------------- correction des profils instables 76 DO 1080 i = 1, klon 77 IF (modif(i)) THEN 78 k2 = limbas(i) 79 8000 CONTINUE 80 k2 = k2 + 1 81 IF (k2 .GT. limhau) goto 8001 82 IF (zh(i,k2) .LT. zh(i,k2-1)) THEN 83 k1 = k2 - 1 84 k = k1 85 sm = zpkdp(i,k2) 86 hm = zh(i,k2) 87 qm = zq(i,k2) 88 8020 CONTINUE 89 sm = sm +zpkdp(i,k) 90 hm = hm +zpkdp(i,k) * (zh(i,k)-hm) / sm 91 qm = qm +zpkdp(i,k) * (zq(i,k)-qm) / sm 92 down = .FALSE. 93 IF (k1 .ne. limbas(i)) THEN 94 IF (hm .LT. zh(i,k1-1)) down = .TRUE. 95 ENDIF 96 IF (down) THEN 97 k1 = k1 - 1 98 k = k1 99 ELSE 100 IF ((k2 .EQ. limhau)) GOTO 8021 101 IF ((zh(i,k2+1).GE.hm)) GOTO 8021 102 k2 = k2 + 1 103 k = k2 104 ENDIF 105 GOTO 8020 106 8021 CONTINUE 107 c------------ nouveau profil : constant (valeur moyenne) 108 DO k = k1, k2 109 zh(i,k) = hm 110 zq(i,k) = qm 111 ENDDO 112 k2 = k2 + 1 113 ENDIF 114 GOTO 8000 115 8001 CONTINUE 116 ENDIF 117 1080 CONTINUE 118 c 119 DO k = 1, limhau 120 DO i = 1, klon 121 d_t(i,k) = (zh(i,k)-zho(i,k))*zpk(i,k)/RCPD 122 d_q(i,k) = zq(i,k) - q(i,k) 123 ENDDO 124 ENDDO 125 c 126 ! FH : les d_q et d_t sont maintenant calcules de facon a valoir 127 ! effectivement 0. si on ne fait rien. 128 ! 129 ! IF (limbas.GT.1) THEN 130 ! DO k = 1, limbas-1 131 ! DO i = 1, klon 132 ! d_t(i,k) = 0.0 133 ! d_q(i,k) = 0.0 134 ! ENDDO 135 ! ENDDO 136 ! ENDIF 137 c 138 ! IF (limhau.LT.klev) THEN 139 ! DO k = limhau+1, klev 140 ! DO i = 1, klon 141 ! d_t(i,k) = 0.0 142 ! d_q(i,k) = 0.0 143 ! ENDDO 144 ! ENDDO 145 ! ENDIF 146 c 147 IF (.NOT.mixq) THEN 148 DO k = 1, klev 149 DO i = 1, klon 150 d_q(i,k) = 0.0 151 ENDDO 152 ENDDO 153 ENDIF 154 c 155 RETURN 156 END 157 158 SUBROUTINE ajsec_convV2(paprs, pplay, t,q, d_t,d_q) 5 159 USE dimphy 6 160 IMPLICIT none -
LMDZ4/trunk/libf/phytherm/clouds_gno.F
r814 r852 219 219 c -- test numerical convergence: 220 220 221 c print*,'avant test ',i,k,lconv(i),u(i),v(i) 221 if (beta(i).lt.1.e-10) then 222 print*,'avant test ',i,k,lconv(i),u(i),v(i),beta(i) 223 stop 224 endif 225 if (abs(fprime(i)).lt.1.e-11) then 226 print*,'avant test fprime<.e-11 ' 227 s ,i,k,lconv(i),u(i),v(i),beta(i),fprime(i) 228 print*,'klon,ND,R,RS,QSUB', 229 s klon,ND,R(i,k),rs(i,k),qsub(i,k) 230 fprime(i)=sign(1.e-11,fprime(i)) 231 endif 232 233 222 234 if ( ABS(dist(i)/beta(i)) .LT. epsilon ) then 223 235 c print*,'v-u **2',(v(i)-u(i))**2 -
LMDZ4/trunk/libf/phytherm/coef_diff_turb_mod.F90
r814 r852 18 18 SUBROUTINE coef_diff_turb(dtime, nsrf, knon, ni, & 19 19 ypaprs, ypplay, yu, yv, yq, yt, yts, yrugos, yqsurf, & 20 ycoefm, ycoefh )20 ycoefm, ycoefh ,yq2) 21 21 ! 22 22 ! Calculate coefficients(ycoefm, ycoefh) for turbulent diffusion in the … … 35 35 REAL, DIMENSION(klon,klev), INTENT(IN) :: yq, yt 36 36 REAL, DIMENSION(klon), INTENT(IN) :: yts, yrugos, yqsurf 37 REAL, DIMENSION(klon,klev+1) :: yq2 37 38 38 39 ! Output arguments … … 41 42 REAL, DIMENSION(klon,klev), INTENT(OUT) :: ycoefm 42 43 43 ! Local variables with attribute SAVE44 !****************************************************************************************45 REAL, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: q2save46 !$OMP THREADPRIVATE(q2save)47 LOGICAL, SAVE :: first_call=.TRUE.48 !$OMP THREADPRIVATE(first_call)49 50 44 ! Other local variables 51 45 !**************************************************************************************** 52 46 INTEGER :: k, i, j 53 47 REAL, DIMENSION(klon,klev) :: ycoefm0, ycoefh0, yzlay, yteta 54 REAL, DIMENSION(klon,klev+1) :: yzlev, yq2,q2diag, ykmm, ykmn, ykmq48 REAL, DIMENSION(klon,klev+1) :: yzlev, q2diag, ykmm, ykmn, ykmq 55 49 REAL, DIMENSION(klon) :: y_cd_h, y_cd_m 56 50 REAL, DIMENSION(klon) :: yustar … … 94 88 ycoefm0, ycoefh0) 95 89 96 DO k = 1, klev90 DO k = 2, klev 97 91 DO i = 1, knon 98 92 ycoefm(i,k) = MAX(ycoefm(i,k),ycoefm0(i,k)) … … 121 115 ycoefm0,ycoefh0) 122 116 123 DO k = 1, klev117 DO k = 2, klev 124 118 DO i = 1, knon 125 119 ycoefm(i,k) = MAX(ycoefm(i,k),ycoefm0(i,k)) … … 137 131 138 132 IF (iflag_pbl.GE.3) THEN 139 140 IF (first_call) THEN141 PRINT*, "Using method MELLOR&YAMADA"142 ! NB! q2save could/should be read and written from (re)startphy.nc143 ALLOCATE(q2save(klon,klev+1,nbsrf))144 q2save(:,:,:) = 1.e-8145 first_call=.FALSE.146 END IF147 133 148 134 yzlay(1:knon,1)= & … … 173 159 END DO 174 160 175 ! Compresse q2save 176 DO k = 1, klev+1 177 DO j = 1, knon 178 i = ni(j) 179 yq2(j,k)=q2save(i,k,nsrf) 180 ENDDO 181 ENDDO 182 183 ! Bug introduit volontairement pour converger avec les resultats 184 ! du papier sur les thermiques. 185 IF (1.EQ.1) THEN 186 y_cd_m(1:knon) = ycoefm(1:knon,1) 187 y_cd_h(1:knon) = ycoefh(1:knon,1) 188 ELSE 189 y_cd_h(1:knon) = ycoefm(1:knon,1) 190 y_cd_m(1:knon) = ycoefh(1:knon,1) 191 ENDIF 161 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 162 ! Pour memoire, le papier Hourdin et al. 2002 a ete obtenur avec un 163 ! bug sur les coefficients de surface : 164 ! y_cd_h(1:knon) = ycoefm(1:knon,1) 165 ! y_cd_m(1:knon) = ycoefh(1:knon,1) 166 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 167 168 y_cd_m(1:knon) = ycoefm(1:knon,1) 169 y_cd_h(1:knon) = ycoefh(1:knon,1) 192 170 193 171 CALL ustarhb(knon,yu,yv,y_cd_m, yustar) … … 215 193 ycoefh(1:knon,2:klev)=ykmn(1:knon,2:klev) 216 194 217 ! update q2save with yq2218 DO j=1,knon219 DO k=1,klev+1220 i=ni(j)221 q2save(i,k,nsrf) = yq2(j,k)222 END DO223 END DO224 225 195 ENDIF !(iflag_pbl.ge.3) 226 196 -
LMDZ4/trunk/libf/phytherm/ini_histday.h
r814 r852 55 55 CALL histdef(nid_day, "weakinv", "Weak inversion", "", 56 56 s iim,jjmp1,nhori, 1,1,1, -99, 32, 57 s "ave(X)", zsto ,zout)57 s "ave(X)", zstophy,zout) 58 58 59 59 CALL histdef(nid_day, "dthmin", "dTheta mini", "K/m", 60 60 s iim,jjmp1,nhori, 1,1,1, -99, 32, 61 s "ave(X)", zsto ,zout)61 s "ave(X)", zstophy,zout) 62 62 63 63 … … 547 547 $ "ave(X)", zstophy,zout) 548 548 C 549 550 551 ! FH Sorties specifiques pour Mellor et Yamada 552 if (iflag_pbl>1 .and. lev_histday.gt.10 ) then 553 call histdef(nid_day, "tke_"//clnsurf(nsrf), 554 $ "Max Turb. Kinetic Energy "//clnsurf(nsrf), "-", 555 $ iim,jj_nb,nhori, klev,1,klev,nvert, 32, 556 $ "ave(X)", zstophy,zout) 557 558 call histdef(nid_day, "tke_max_"//clnsurf(nsrf), 559 $ "Max Turb. Kinetic Energy "//clnsurf(nsrf), "-", 560 $ iim,jj_nb,nhori, klev,1,klev,nvert, 32, 561 $ "t_max(X)", zstophy,zout) 562 endif 563 564 C 549 565 END DO 550 566 C -
LMDZ4/trunk/libf/phytherm/ini_histmth.h
r841 r852 285 285 c 286 286 CALL histdef(nid_mth, "sens", "Sensible heat flux", "W/m2", 287 . iim,jj mp1,nhori, 1,1,1, -99, 32,287 . iim,jj_nb,nhori, 1,1,1, -99, 32, 288 288 . "ave(X)", zstophy,zout) 289 289 c … … 310 310 CALL histdef(nid_mth, "fqfonte","Land ice melt", 311 311 . "kg/m2/s",iim,jj_nb,nhori, 1,1,1, -99, 32, 312 . "ave(X)", zsto ,zout)312 . "ave(X)", zstophy,zout) 313 313 314 314 DO nsrf = 1, nbsrf … … 700 700 IF(lev_histmth.GE.4) THEN 701 701 c 702 !FH Sorties pour la couche limite 703 if (iflag_pbl>1) then 704 CALL histdef(nid_mth, "tke","TKE","m2/s2", 705 . iim,jj_nb,nhori, klev,1,klev,nvert, 32, 706 . "ave(X)", zstophy,zout) 707 c 708 CALL histdef(nid_mth, "tke_max","TKE max","m2/s2", 709 . iim,jj_nb,nhori, klev,1,klev,nvert, 32, 710 . "t_max(X)", zstophy,zout) 711 endif 712 c 713 CALL histdef(nid_mth, "kz","Kz melange","m2/s", 714 . iim,jj_nb,nhori, klev,1,klev,nvert, 32, 715 . "ave(X)", zstophy,zout) 716 c 717 CALL histdef(nid_mth, "kz_max","Kz melange max","m2/s", 718 . iim,jj_nb,nhori, klev,1,klev,nvert, 32, 719 . "t_max(X)", zstophy,zout) 720 721 702 722 CALL histdef(nid_mth, "clwcon", 703 723 . "Convective Cloud Liquid water content" … … 730 750 . "ave(X)", zstophy,zout) 731 751 c 752 753 754 c 732 755 c CALL histdef(nid_mth, "ducon", "Convection du", "m/s2", 733 756 c . iim,jj_nb,nhori, klev,1,klev,nvert, 32, … … 774 797 . "ave(X)", zstophy,zout) 775 798 799 c 800 CALL histdef(nid_mth, "dtthe", "Dry adjust. dT", "K/s", 801 . iim,jj_nb,nhori, klev,1,klev,nvert, 32, 802 . "ave(X)", zstophy,zout) 803 804 CALL histdef(nid_mth,"dqthe","Dry adjust. dQ","(kg/kg)/s", 805 . iim,jj_nb,nhori, klev,1,klev,nvert, 32, 806 . "ave(X)", zstophy,zout) 776 807 c 777 808 CALL histdef(nid_mth, "dtajs", "Dry adjust. dT", "K/s", … … 914 945 c 915 946 CALL histdef(nid_mth, "dtcon","dtcon","K/s", 947 . iim,jj_nb,nhori, klev,1,klev,nvert, 32, 948 . "ave(X)", zstophy,zout) 949 c 950 CALL histdef(nid_mth, "dtthe", "Dry adjust. dT", "K/s", 951 . iim,jj_nb,nhori, klev,1,klev,nvert, 32, 952 . "ave(X)", zstophy,zout) 953 954 CALL histdef(nid_mth,"dqthe","Dry adjust. dQ","(kg/kg)/s", 916 955 . iim,jj_nb,nhori, klev,1,klev,nvert, 32, 917 956 . "ave(X)", zstophy,zout) … … 1622 1661 IF(lev_histmth.GE.4) THEN 1623 1662 c 1663 !FH Sorties pour la couche limite 1664 if (iflag_pbl>1) then 1665 CALL histdef(nid_mth, "tke","TKE","m2/s2", 1666 . iim,jj_nb,nhori, klev,1,klev,nvert, 32, 1667 . "ave(X)", zstophy,zout) 1668 c 1669 CALL histdef(nid_mth, "tke_max","TKE max","m2/s2", 1670 . iim,jj_nb,nhori, klev,1,klev,nvert, 32, 1671 . "t_max(X)", zstophy,zout) 1672 endif 1673 c 1674 CALL histdef(nid_mth, "kz","Kz melange","m2/s", 1675 . iim,jj_nb,nhori, klev,1,klev,nvert, 32, 1676 . "ave(X)", zstophy,zout) 1677 c 1678 CALL histdef(nid_mth, "kz_max","Kz melange max","m2/s", 1679 . iim,jj_nb,nhori, klev,1,klev,nvert, 32, 1680 . "t_max(X)", zstophy,zout) 1681 1682 1683 1684 1624 1685 CALL histdef(nid_mth, "clwcon", 1625 1686 . "Convective Cloud Liquid water content" -
LMDZ4/trunk/libf/phytherm/pbl_surface_mod.F90
r815 r852 197 197 wfbils, wfbilo, flux_t, flux_u, flux_v,& 198 198 dflux_t, dflux_q, zxsnow, & 199 zxfluxt, zxfluxq, q2m, flux_q )199 zxfluxt, zxfluxq, q2m, flux_q, tke ) 200 200 !**************************************************************************************** 201 201 ! Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818 … … 240 240 ! flux_t---output-R- flux de chaleur sensible (CpT) J/m**2/s (W/m**2) 241 241 ! (orientation positive vers le bas) 242 ! tke---input/output-R- tke (kg/m**2/s) 242 243 ! flux_q---output-R- flux de vapeur d'eau (kg/m**2/s) 243 244 ! flux_u---output-R- tension du vent X: (kg m/s)/(m**2 s) ou Pascal … … 348 349 REAL, DIMENSION(klon, nbsrf),INTENT(OUT) :: q2m 349 350 REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_q 351 352 ! Input/output 353 REAL, DIMENSION(klon, klev+1, nbsrf), INTENT(INOUT) :: tke 350 354 351 355 … … 421 425 REAL, DIMENSION(klon,klev) :: delp 422 426 REAL, DIMENSION(klon,klev+1) :: ypaprs 427 REAL, DIMENSION(klon,klev+1) :: ytke 423 428 REAL, DIMENSION(klon,nsoilmx) :: ytsoil 424 429 REAL, DIMENSION(klon,nbsrf) :: pctsrf_pot … … 451 456 452 457 458 !jg+ temporary test 459 REAL, DIMENSION(klon,klev) :: y_flux_u_old, y_flux_v_old 460 REAL, DIMENSION(klon,klev) :: y_d_u_old, y_d_v_old 461 !jg- 462 453 463 !**************************************************************************************** 454 464 ! End of declarations … … 514 524 flux_u = 0.0 ; flux_v = 0.0 ; d_t = 0.0 ; d_q = 0.0 515 525 d_u = 0.0 ; d_v = 0.0 ; zcoefh = 0.0 ; yqsol = 0.0 516 ytherm = 0.0 526 ytherm = 0.0 ; ytke=0. 517 527 518 528 ytsoil = 999999. … … 657 667 ypplay(j,k) = pplay(i,k) 658 668 ydelp(j,k) = delp(i,k) 669 ytke(j,k)=tke(i,k,nsrf) 659 670 yu(j,k) = u(i,k) 660 671 yv(j,k) = v(i,k) … … 687 698 CALL coef_diff_turb(dtime, nsrf, knon, ni, & 688 699 ypaprs, ypplay, yu, yv, yq, yt, yts, yrugos, yqsurf, & 689 ycoefm, ycoefh )700 ycoefm, ycoefh,ytke) 690 701 691 702 !**************************************************************************************** … … 839 850 !**************************************************************************************** 840 851 852 tke(:,:,nsrf)=0. 841 853 DO k = 1, klev 842 854 DO j = 1, knon … … 853 865 flux_u(i,k,nsrf) = y_flux_u(j,k) 854 866 flux_v(i,k,nsrf) = y_flux_v(j,k) 867 868 tke(i,k,nsrf)=ytke(j,k) 869 855 870 ENDDO 856 871 ENDDO -
LMDZ4/trunk/libf/phytherm/phyetat0.F
r814 r852 13 13 . rugsrel_p,tabcntr0, 14 14 . t_ancien_p,q_ancien_p,ancien_ok_p, rnebcon_p, ratqs_p, 15 . clwcon_p )15 . clwcon_p,pbl_tke_p) 16 16 17 17 USE dimphy … … 37 37 #include "clesphys.h" 38 38 #include "temps.h" 39 #include "thermcell.h" 40 #include "compbl.h" 39 41 c====================================================================== 40 42 CHARACTER*(*) fichnom … … 45 47 REAL solaire_etat0 46 48 REAL tsol_p(klon,nbsrf) 49 REAL pbl_tke_p(klon,klev,nbsrf) 47 50 REAL tsoil_p(klon,nsoilmx,nbsrf) 48 51 REAL tslab_p(klon), seaice_p(klon) … … 83 86 REAL rlat(klon_glo), rlon(klon_glo) 84 87 REAL tsol(klon_glo,nbsrf) 88 REAL pbl_tke(klon_glo,klev,nbsrf) 85 89 REAL tsoil(klon_glo,nsoilmx,nbsrf) 86 90 cIM "slab" ocean … … 127 131 c 128 132 INTEGER nid, nvarid 129 INTEGER ierr, i, nsrf, isoil 133 INTEGER ierr, i, nsrf, isoil ,k 130 134 INTEGER length 131 135 PARAMETER (length=100) … … 1492 1496 xmax = MAXval(run_off_lic_0) 1493 1497 PRINT*,'(ecart-type) run_off_lic_0:', xmin, xmax 1498 1499 1500 c Lecture de l'energie cinetique turbulente 1501 c 1502 1503 IF (iflag_pbl>1) then 1504 PRINT*, 'phyetat0: Le champ <TKE> est absent' 1505 PRINT*, ' Mais je vais essayer de lire TKE**' 1506 DO nsrf = 1, nbsrf 1507 IF (nsrf.GT.99) THEN 1508 PRINT*, "Trop de sous-mailles" 1509 CALL abort 1510 ENDIF 1511 WRITE(str2,'(i2.2)') nsrf 1512 ierr = NF_INQ_VARID (nid, "TKE"//str2, nvarid) 1513 IF (ierr.NE.NF_NOERR) THEN 1514 PRINT*, "WARNING phyetat0: <TKE"//str2//"> est absent" 1515 pbl_tke(:,:,nsrf)=1.e-8 1516 ELSE 1517 #ifdef NC_DOUBLE 1518 ierr = NF_GET_VAR_DOUBLE(nid, nvarid, pbl_tke(1,1,nsrf)) 1519 #else 1520 ierr = NF_GET_VAR_REAL(nid, nvarid, pbl_tke(1,1,nsrf)) 1521 #endif 1522 IF (ierr.NE.NF_NOERR) THEN 1523 PRINT*, "WARNING phyetat0: echec lect <TKE"//str2//">" 1524 CALL abort 1525 ENDIF 1526 ENDIF 1527 1528 xmin = 1.0E+20 1529 xmax = -1.0E+20 1530 DO k = 1, klev 1531 DO i = 1, klon_glo 1532 xmin = MIN(pbl_tke(i,k,nsrf),xmin) 1533 xmax = MAX(pbl_tke(i,k,nsrf),xmax) 1534 ENDDO 1535 ENDDO 1536 PRINT*,'Temperature du sol TKE**:', nsrf, xmin, xmax 1537 ENDDO 1538 ENDIF 1539 1494 1540 c 1495 1541 c Fermer le fichier: … … 1516 1562 call Scatter( rlon,rlon_p) 1517 1563 call Scatter( tsol,tsol_p) 1564 IF (iflag_pbl>1) then 1565 call Scatter( pbl_tke,pbl_tke_p) 1566 endif 1518 1567 call Scatter( tsoil,tsoil_p) 1519 1568 call Scatter( tslab,tslab_p) -
LMDZ4/trunk/libf/phytherm/phyredem.F
r814 r852 9 9 . radsol_p,zmea_p,zstd_p,zsig_p, 10 10 . zgam_p,zthe_p,zpic_p,zval_p,rugsrel_p, 11 . t_ancien_p, q_ancien_p, rnebcon_p, ratqs_p, clwcon_p) 11 . t_ancien_p, q_ancien_p, rnebcon_p, ratqs_p, clwcon_p, 12 . pbl_tke_p) 12 13 13 14 USE dimphy … … 29 30 #include "control.h" 30 31 #include "temps.h" 32 #include "thermcell.h" 33 #include "compbl.h" 31 34 c====================================================================== 32 35 CHARACTER*(*) fichnom … … 35 38 REAL rlat_p(klon), rlon_p(klon) 36 39 REAL tsol_p(klon,nbsrf) 40 REAL pbl_tke_p(klon,klev,nbsrf) 37 41 REAL tsoil_p(klon,nsoilmx,nbsrf) 38 42 CHARACTER*6 ocean … … 70 74 REAL rlat(klon_glo), rlon(klon_glo) 71 75 REAL tsol(klon_glo,nbsrf) 76 REAL pbl_tke(klon_glo,klev,nbsrf) 72 77 REAL tsoil(klon_glo,nsoilmx,nbsrf) 73 78 REAL tslab(klon_glo), seaice(klon_glo) … … 135 140 call Gather( rlon_p,rlon) 136 141 call Gather( tsol_p,tsol) 142 call Gather( pbl_tke_p,pbl_tke) 137 143 call Gather( tsoil_p,tsoil) 138 144 call Gather( tslab_p,tslab) … … 899 905 c 900 906 c 907 !!!!!!!!!!!!!!!!!!!! DEB TKE PBL !!!!!!!!!!!!!!!!!!!!!!!!! 908 c 909 IF (iflag_pbl>1) then 910 DO nsrf = 1, nbsrf 911 IF (nsrf.LE.99) THEN 912 WRITE(str2,'(i2.2)') nsrf 913 ierr = NF_REDEF (nid) 914 #ifdef NC_DOUBLE 915 ierr = NF_DEF_VAR (nid,"TKE"//str2,NF_DOUBLE,1,idim3 916 $ ,nvarid) 917 #else 918 ierr = NF_DEF_VAR (nid,"TKE"//str2,NF_FLOAT,1,idim3 919 $ ,nvarid) 920 #endif 921 ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 15, 922 . "Energ. Cineti. Turb."//str2) 923 ierr = NF_ENDDEF(nid) 924 ELSE 925 PRINT*, "Trop de sous-mailles" 926 CALL abort 927 ENDIF 928 #ifdef NC_DOUBLE 929 ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,pbl_tke(:,:,nsrf)) 930 #else 931 ierr = NF_PUT_VAR_REAL (nid,nvarid,pbl_tke(:,:,nsrf)) 932 #endif 933 ENDDO 934 ENDIF 935 936 !!!!!!!!!!!!!!!!!!!! FIN TKE PBL !!!!!!!!!!!!!!!!!!!!!!!!! 937 c 901 938 ierr = NF_CLOSE(nid) 902 939 c -
LMDZ4/trunk/libf/phytherm/physiq.F
r839 r852 189 189 190 190 integer lmax_th(klon) 191 integer limbas(klon) 191 192 real ratqscth(klon,klev) 192 193 real ratqsdiff(klon,klev) … … 945 946 REAL fluxu(klon,klev, nbsrf) ! flux turbulent de vitesse u 946 947 REAL fluxv(klon,klev, nbsrf) ! flux turbulent de vitesse v 948 949 REAL, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: pbl_tke ! turb kinetic energy 950 c !$OMP THREADPRIVATE(pbl_tke) 951 947 952 c 948 953 REAL zxfluxt(klon, klev) … … 1066 1071 c$OMP THREADPRIVATE(d_u_con,d_v_con) 1067 1072 REAL d_t_lsc(klon,klev),d_q_lsc(klon,klev),d_ql_lsc(klon,klev) 1073 REAL d_t_ajsb(klon,klev), d_q_ajsb(klon,klev) 1068 1074 REAL d_t_ajs(klon,klev), d_q_ajs(klon,klev) 1069 1075 REAL d_u_ajs(klon,klev), d_v_ajs(klon,klev) … … 1554 1560 itaprad = 0 1555 1561 1562 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1563 !! Un petit travail à faire ici. 1564 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1565 1566 if (iflag_pbl>1) then 1567 PRINT*, "Using method MELLOR&YAMADA" 1568 endif 1569 ! NB! pbl_tke could/should be read and written from (re)startphy.nc 1570 ALLOCATE(pbl_tke(klon,klev+1,nbsrf)) 1571 pbl_tke(:,:,:) = 1.e-8 1572 1573 1556 1574 CALL phyetat0 ("startphy.nc",dtime,co2_ppm_etat0,solaire_etat0, 1557 1575 . rlat,rlon,pctsrf, ftsol, … … 1561 1579 . radsol,clesphy0, 1562 1580 . zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,tabcntr0, 1563 . t_ancien, q_ancien, ancien_ok, rnebcon, ratqs,clwcon) 1581 . t_ancien, q_ancien, ancien_ok, rnebcon, ratqs,clwcon, 1582 . pbl_tke) 1583 1584 1585 1586 1587 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1588 1564 1589 1565 1590 DO i=1,klon … … 1643 1668 ENDIF 1644 1669 1670 rugoro=0. 1645 1671 c34EK 1646 1672 IF (ok_orodr) THEN 1647 DO i=1,klon 1648 rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0) 1649 ENDDO 1673 1674 rugoro=0. 1675 1676 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1677 ! FH sans doute a enlever de finitivement ou, si on le garde, l'activer 1678 ! justement quand ok_orodr = false. 1679 ! ce rugoro est utilise par la couche limite et fait double emploi 1680 ! avec les paramétrisations spécifiques de Francois Lott. 1681 ! DO i=1,klon 1682 ! rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0) 1683 ! ENDDO 1684 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1685 1650 1686 CALL SUGWD(klon,klev,paprs,pplay) 1651 1687 DO i=1,klon … … 2061 2097 d wfbils, wfbilo, fluxt, fluxu, fluxv, 2062 2098 - dsens, devap, zxsnow, 2063 - zxfluxt, zxfluxq, q2m, fluxq )2099 - zxfluxt, zxfluxq, q2m, fluxq, pbl_tke ) 2064 2100 c 2065 2101 c … … 2233 2269 DO i = 1, klon 2234 2270 ema_pct(i) = paprs(i,itop_con(i)) 2271 if (itop_con(i).gt.klev-3) then 2272 print*,'La convection monte trop haut ' 2273 print*,'itop_con(,',i,',)=',itop_con(i) 2274 endif 2235 2275 ENDDO 2236 2276 DO i = 1, klon … … 2321 2361 2322 2362 2363 d_t_ajsb(:,:)=0. 2364 d_q_ajsb(:,:)=0. 2323 2365 d_t_ajs(:,:)=0. 2324 2366 d_u_ajs(:,:)=0. … … 2336 2378 c ==== 2337 2379 IF(prt_level>9)WRITE(lunout,*)'pas de convection' 2338 else if(iflag_thermals.eq.0) then 2339 2340 c Ajustement sec 2341 c ============== 2342 IF(prt_level>9)WRITE(lunout,*)'ajsec' 2343 CALL ajsec(paprs, pplay, t_seri,q_seri, d_t_ajs, d_q_ajs) 2344 t_seri(:,:) = t_seri(:,:) + d_t_ajs(:,:) 2345 q_seri(:,:) = q_seri(:,:) + d_q_ajs(:,:) 2380 2381 2346 2382 else 2383 2347 2384 c Thermiques 2348 2385 c ========== … … 2351 2388 print*,'JUSTE AVANT , iflag_thermals=' 2352 2389 s ,iflag_thermals,' nsplit_thermals=',nsplit_thermals 2390 2391 2392 if (iflag_thermals.gt.1) then 2353 2393 call calltherm(pdtphys 2354 2394 s ,pplay,paprs,pphi,weak_inversion … … 2357 2397 s ,fm_therm,entr_therm,zqasc,clwcon0th,lmax_th,ratqscth, 2358 2398 s ratqsdiff,zqsatth) 2399 endif 2359 2400 2360 2401 ! call calltherm(pdtphys … … 2363 2404 ! s ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs 2364 2405 ! s ,fm_therm,entr_therm) 2406 2407 c Ajustement sec 2408 c ============== 2409 2410 ! Dans le cas où on active les thermiques, on fait partir l'ajustement 2411 ! a partir du sommet des thermiques. 2412 ! Dans le cas contraire, on demarre au niveau 1. 2413 2414 if (iflag_thermals.ge.13.or.iflag_thermals.eq.0) then 2415 2416 if(iflag_thermals.eq.0) then 2417 IF(prt_level>9)WRITE(lunout,*)'ajsec' 2418 limbas(:)=1 2419 else 2420 limbas(:)=lmax_th(:) 2421 endif 2422 2423 ! Attention : le call ajsec_convV2 n'est maintenu que momentanneement 2424 ! pour des test de convergence numerique. 2425 ! Le nouveau ajsec est a priori mieux, meme pour le cas 2426 ! iflag_thermals = 0 (l'ancienne version peut faire des tendances 2427 ! non nulles numeriquement pour des mailles non concernees. 2428 2429 if (iflag_thermals.eq.0) then 2430 CALL ajsec_convV2(paprs, pplay, t_seri,q_seri 2431 s , d_t_ajsb, d_q_ajsb) 2432 else 2433 CALL ajsec(paprs, pplay, t_seri,q_seri,limbas 2434 s , d_t_ajsb, d_q_ajsb) 2435 endif 2436 2437 t_seri(:,:) = t_seri(:,:) + d_t_ajsb(:,:) 2438 q_seri(:,:) = q_seri(:,:) + d_q_ajsb(:,:) 2439 d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:) 2440 d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:) 2441 2442 endif 2443 2365 2444 endif 2366 2445 c … … 2397 2476 ptconvth(:,:)=.false. 2398 2477 ratqsc(:,:)=0. 2399 print*,'avant clouds_gno '2478 print*,'avant clouds_gno thermique' 2400 2479 call clouds_gno 2401 2480 s (klon,klev,q_seri,zqsat,clwcon0th,ptconvth,ratqsc,rnebcon0th) … … 3321 3400 . radsol, 3322 3401 . zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro, 3323 . t_ancien, q_ancien, rnebcon, ratqs, clwcon) 3402 . t_ancien, q_ancien, rnebcon, ratqs, clwcon, 3403 . pbl_tke) 3324 3404 ENDIF 3325 3405 -
LMDZ4/trunk/libf/phytherm/thermcell_closure.F90
r814 r852 18 18 REAL zmax(ngrid),zmax_sec(ngrid) 19 19 REAL wmax(ngrid),wmax_sec(ngrid) 20 real zdenom 20 21 21 22 REAL alim_star2(ngrid) … … 35 36 & /(rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k))) 36 37 enddo 38 zdenom=max(500.,zmax(ig))*r_aspect*alim_star2(ig) 39 if (zdenom<1.e-14) then 40 print*,'ig=',ig 41 print*,'alim_star2',alim_star2(ig) 42 print*,'zmax',zmax(ig) 43 print*,'r_aspect',r_aspect 44 print*,'zdenom',zdenom 45 print*,'alim_star',alim_star(ig,:) 46 print*,'zmax_sec',zmax_sec(ig) 47 print*,'wmax_sec',wmax_sec(ig) 48 stop 49 endif 37 50 if ((zmax_sec(ig).gt.1.e-10).and.(1.eq.1)) then 38 51 f(ig)=wmax_sec(ig)/(max(500.,zmax_sec(ig))*r_aspect & … … 41 54 & zmax_sec(ig))*wmax_sec(ig)) 42 55 else 43 f(ig)=wmax(ig)/ (max(500.,zmax(ig))*r_aspect*alim_star2(ig))56 f(ig)=wmax(ig)/zdenom 44 57 f(ig)=f(ig)+(f0(ig)-f(ig))*exp((-ptimestep/ & 45 58 & zmax(ig))*wmax(ig)) -
LMDZ4/trunk/libf/phytherm/thermcell_dry.F90
r814 r852 1 SUBROUTINE thermcell_dry(ngrid,nlay,zlev,pphi,ztv, entr_star, &2 & l entr,lmin,zmax,wmax,lev_out)1 SUBROUTINE thermcell_dry(ngrid,nlay,zlev,pphi,ztv,alim_star, & 2 & lalim,lmin,zmax,wmax,lev_out) 3 3 4 4 !-------------------------------------------------------------------------- … … 13 13 REAL pphi(ngrid,nlay) 14 14 REAl ztv(ngrid,nlay) 15 REAL entr_star(ngrid,nlay)16 INTEGER l entr(ngrid)15 REAL alim_star(ngrid,nlay) 16 INTEGER lalim(ngrid) 17 17 integer lev_out ! niveau pour les print 18 18 … … 33 33 do l=1,nlay+1 34 34 zw2(ig,l)=0. 35 f_star(ig,l)=0.36 35 wa_moy(ig,l)=0. 37 36 enddo … … 47 46 enddo 48 47 !calcul de la vitesse a partir de la CAPE en melangeant thetav 48 49 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 50 ! A eliminer 51 ! Ce if complique etait fait pour reperer la premiere couche instable 52 ! Ici, c'est lmin. 53 ! 54 ! do l=1,nlay-2 55 ! do ig=1,ngrid 56 ! if (ztv(ig,l).gt.ztv(ig,l+1) & 57 ! & .and.alim_star(ig,l).gt.1.e-10 & 58 ! & .and.zw2(ig,l).lt.1e-10) then 59 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 60 61 62 ! Calcul des F^*, integrale verticale de E^* 63 f_star(:,1)=0. 64 do l=1,nlay 65 f_star(:,l+1)=f_star(:,l)+alim_star(:,l) 66 enddo 67 68 ! niveau (reel) auquel zw2 s'annule FH :n'etait pas initialise 69 linter(:)=0. 70 71 ! couche la plus haute concernee par le thermique. 72 lmax(:)=1 73 74 ! Le niveau linter est une variable continue qui se trouve dans la couche 75 ! lmax 76 49 77 do l=1,nlay-2 50 78 do ig=1,ngrid 51 if (ztv(ig,l).gt.ztv(ig,l+1) & 52 & .and.entr_star(ig,l).gt.1.e-10 & 53 & .and.zw2(ig,l).lt.1e-10) then 54 f_star(ig,l+1)=entr_star(ig,l) 55 ! 79 if (l.eq.lmin(ig).and.lalim(ig).gt.1) then 80 81 !------------------------------------------------------------------------ 82 ! Calcul de la vitesse en haut de la premiere couche instable. 83 ! Premiere couche du panache thermique 84 !------------------------------------------------------------------------ 56 85 zw2(ig,l+1)=2.*RG*(ztv(ig,l)-ztv(ig,l+1))/ztv(ig,l+1) & 57 86 & *(zlev(ig,l+1)-zlev(ig,l)) & 58 87 & *0.4*pphi(ig,l)/(pphi(ig,l+1)-pphi(ig,l)) 59 else if ((zw2(ig,l).ge.1e-10).and. & 60 & (f_star(ig,l)+entr_star(ig,l).gt.1.e-10)) then 61 f_star(ig,l+1)=f_star(ig,l)+entr_star(ig,l) 62 ztva(ig,l)=(f_star(ig,l)*ztva(ig,l-1)+entr_star(ig,l) & 88 89 !------------------------------------------------------------------------ 90 ! Tant que la vitesse en bas de la couche et la somme du flux de masse 91 ! et de l'entrainement (c'est a dire le flux de masse en haut) sont 92 ! positifs, on calcul 93 ! 1. le flux de masse en haut f_star(ig,l+1) 94 ! 2. la temperature potentielle virtuelle dans la couche ztva(ig,l) 95 ! 3. la vitesse au carré en haut zw2(ig,l+1) 96 !------------------------------------------------------------------------ 97 98 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 99 ! A eliminer : dans cette version, si zw2 est > 0 on a un therique. 100 ! et donc, au dessus, f_star(ig,l+1) est forcement suffisamment 101 ! grand puisque on n'a pas de detrainement. 102 ! f_star est une fonction croissante. 103 ! c'est donc vraiment sur zw2 uniquement qu'il faut faire le test. 104 ! else if ((zw2(ig,l).ge.1e-10).and. & 105 ! & (f_star(ig,l)+alim_star(ig,l).gt.1.e-10)) then 106 ! f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l) 107 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 108 109 else if (zw2(ig,l).ge.1e-10) then 110 111 ztva(ig,l)=(f_star(ig,l)*ztva(ig,l-1)+alim_star(ig,l) & 63 112 & *ztv(ig,l))/f_star(ig,l+1) 64 113 zw2(ig,l+1)=zw2(ig,l)*(f_star(ig,l)/f_star(ig,l+1))**2+ & … … 67 116 endif 68 117 ! determination de zmax continu par interpolation lineaire 118 !------------------------------------------------------------------------ 119 120 if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then 121 ! stop'On tombe sur le cas particulier de thermcell_dry' 122 print*,'On tombe sur le cas particulier de thermcell_dry' 123 zw2(ig,l+1)=0. 124 linter(ig)=l+1 125 lmax(ig)=l 126 endif 127 69 128 if (zw2(ig,l+1).lt.0.) then 70 129 linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l)) & 71 130 & -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l)) 72 131 zw2(ig,l+1)=0. 73 else 132 lmax(ig)=l 133 endif 134 74 135 wa_moy(ig,l+1)=sqrt(zw2(ig,l+1)) 75 endif76 136 77 137 if (wa_moy(ig,l+1).gt.wmaxa(ig)) then … … 84 144 if (lev_out.ge.1) print*,'fin calcul zw2' 85 145 ! 146 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 147 ! A eliminer : 148 ! Ce calcul de lmax est fait en meme temps que celui de linter, plus haut 86 149 ! Calcul de la couche correspondant a la hauteur du thermique 87 do ig=1,ngrid 88 lmax(ig)=lentr(ig) 89 enddo 90 do ig=1,ngrid 91 do l=nlay,lentr(ig)+1,-1 92 if (zw2(ig,l).le.1.e-10) then 93 lmax(ig)=l-1 94 endif 95 enddo 96 enddo 150 ! do ig=1,ngrid 151 ! lmax(ig)=lalim(ig) 152 ! enddo 153 ! do ig=1,ngrid 154 ! do l=nlay,lalim(ig)+1,-1 155 ! if (zw2(ig,l).le.1.e-10) then 156 ! lmax(ig)=l-1 157 ! endif 158 ! enddo 159 ! enddo 160 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 161 97 162 ! 98 163 ! Determination de zw2 max … … 119 184 do ig=1,ngrid 120 185 ! calcul de zlevinter 121 zlevinter(ig)=(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))* & 122 & linter(ig)+zlev(ig,lmax(ig))-lmax(ig)*(zlev(ig,lmax(ig)+1) & 123 & -zlev(ig,lmax(ig))) 186 187 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 188 ! FH A eliminer 189 ! Simplification 190 ! zlevinter(ig)=(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))* & 191 ! & linter(ig)+zlev(ig,lmax(ig))-lmax(ig)*(zlev(ig,lmax(ig)+1) & 192 ! & -zlev(ig,lmax(ig))) 193 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 194 195 zlevinter(ig)=zlev(ig,lmax(ig)) + & 196 & (linter(ig)-lmax(ig))*(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig))) 124 197 zmax(ig)=max(zmax(ig),zlevinter(ig)-zlev(ig,lmin(ig))) 125 198 enddo 126 !on stoppe après les calculs de zmax et wmax 199 200 ! Verification que lalim<=lmax 201 do ig=1,ngrid 202 if(lalim(ig)>lmax(ig)) then 203 print*,'WARNING thermcell_dry ig=',ig,' lalim=',lalim(ig),' lmax(ig)=',lmax(ig) 204 lmax(ig)=lalim(ig) 205 endif 206 enddo 127 207 128 208 RETURN -
LMDZ4/trunk/libf/phytherm/thermcell_flux.F90
r814 r852 28 28 29 29 integer ncorecfm1,ncorecfm2,ncorecfm3,ncorecalpha 30 integer ncorecfm4,ncorecfm5,ncorecfm6 30 integer ncorecfm4,ncorecfm5,ncorecfm6,ncorecfm7,ncorecfm8 31 31 32 32 33 33 REAL entr(ngrid,klev),detr(ngrid,klev) 34 34 REAL fm(ngrid,klev+1) 35 REAL zfm 35 36 36 37 integer igout … … 40 41 REAL f_old,ddd0,eee0,ddd,eee,zzz 41 42 42 REAL fracmax 43 save fracmax 44 45 fracmax=0.5 43 REAL fomass_max,alphamax 44 save fomass_max,alphamax 45 46 fomass_max=0.5 47 alphamax=0.7 46 48 47 49 ncorecfm1=0 … … 51 53 ncorecfm5=0 52 54 ncorecfm6=0 55 ncorecfm7=0 56 ncorecfm8=0 53 57 ncorecalpha=0 54 58 … … 56 60 fm(:,:)=0. 57 61 62 if (lev_out.ge.10) then 63 write(lunout,*) 'Dans thermcell_flux 0' 64 write(lunout,*) 'flux base ',f(igout) 65 write(lunout,*) 'lmax ',lmax(igout) 66 write(lunout,*) 'lalim ',lalim(igout) 67 write(lunout,*) 'ig= ',igout 68 write(lunout,*) ' l E* A* D* ' 69 write(lunout,'(i4,3e15.5)') (l,entr_star(igout,l),alim_star(igout,l),detr_star(igout,l) & 70 & ,l=1,lmax(igout)) 71 endif 72 58 73 59 74 !------------------------------------------------------------------------- … … 66 81 if (l.le.lmax(ig)) then 67 82 if (entr_star(ig,l).gt.1.) then 68 print*,'ig,l,lmax(ig)',ig,l,lmax(ig) 83 print*,'WARNING thermcell_flux 1 ig,l,lmax(ig)',ig,l,lmax(ig) 84 print*,'entr_star(ig,l)',entr_star(ig,l) 85 print*,'alim_star(ig,l)',alim_star(ig,l) 86 print*,'detr_star(ig,l)',detr_star(ig,l) 87 ! stop 88 endif 89 else 90 if (abs(entr_star(ig,l))+abs(alim_star(ig,l))+abs(detr_star(ig,l)).gt.0.) then 91 print*,'cas 1 : ig,l,lmax(ig)',ig,l,lmax(ig) 69 92 print*,'entr_star(ig,l)',entr_star(ig,l) 70 93 print*,'alim_star(ig,l)',alim_star(ig,l) … … 72 95 stop 73 96 endif 74 else75 if (abs(entr_star(ig,l))+abs(alim_star(ig,l))+abs(detr_star(ig,l)).gt.0) then76 print*,'ig,l,lmax(ig)',ig,l,lmax(ig)77 print*,'entr_star(ig,l)',entr_star(ig,l)78 print*,'alim_star(ig,l)',alim_star(ig,l)79 print*,'detr_star(ig,l)',detr_star(ig,l)80 stop81 endif82 97 endif 83 98 enddo … … 92 107 detr(:,l)=f(:)*detr_star(:,l) 93 108 enddo 109 110 if (lev_out.ge.10) then 111 write(lunout,*) 'Dans thermcell_flux 1' 112 write(lunout,*) 'flux base ',f(igout) 113 write(lunout,*) 'lmax ',lmax(igout) 114 write(lunout,*) 'lalim ',lalim(igout) 115 write(lunout,*) 'ig= ',igout 116 write(lunout,*) ' l E D W2' 117 write(lunout,'(i4,3e15.5)') (l,entr(igout,l),detr(igout,l) & 118 & ,zw2(igout,l+1),l=1,lmax(igout)) 119 endif 94 120 95 121 fm(:,1)=0. … … 107 133 enddo 108 134 109 if (lev_out.ge.10) then 110 write(lunout,*) 'Dans thermcell_flux 1' 111 write(lunout,*) 'flux base ',f(igout) 112 write(lunout,*) 'lmax ',lmax(igout) 113 write(lunout,*) 'lalim ',lalim(igout) 114 write(lunout,'(i6,i4,3e15.5)') (ig,l,entr(igout,l),detr(igout,l) & 115 & ,fm(igout,l+1),l=1,lmax(igout)) 116 endif 135 136 137 ! Test provisoire : pour comprendre pourquoi on corrige plein de fois 138 ! le cas fm6, on commence par regarder une premiere fois avant les 139 ! autres corrections. 140 141 do l=1,klev 142 do ig=1,ngrid 143 if (detr(ig,l).gt.fm(ig,l)) then 144 ncorecfm8=ncorecfm8+1 145 ! igout=ig 146 endif 147 enddo 148 enddo 149 150 if (lev_out.ge.10) & 151 & call printflux(ngrid,klev,lunout,igout,f,lmax,lalim, & 152 & ptimestep,masse,entr,detr,fm,'2 ') 117 153 118 154 !------------------------------------------------------------------------- … … 130 166 enddo 131 167 enddo 168 169 if (lev_out.ge.10) & 170 & call printflux(ngrid,klev,lunout,igout,f,lmax,lalim, & 171 & ptimestep,masse,entr,detr,fm,'3 ') 132 172 133 173 !------------------------------------------------------------------------- … … 151 191 enddo 152 192 193 if (lev_out.ge.10) & 194 & call printflux(ngrid,klev,lunout,igout,f,lmax,lalim, & 195 & ptimestep,masse,entr,detr,fm,'4 ') 196 153 197 154 198 !------------------------------------------------------------------------- … … 167 211 enddo 168 212 213 if (lev_out.ge.10) & 214 & call printflux(ngrid,klev,lunout,igout,f,lmax,lalim, & 215 & ptimestep,masse,entr,detr,fm,'5 ') 216 169 217 !------------------------------------------------------------------------- 170 218 !detr ne peut pas etre superieur a fm … … 172 220 173 221 if(1.eq.1) then 222 174 223 do l=1,klev 175 224 do ig=1,ngrid … … 181 230 ncorecfm6=ncorecfm6+1 182 231 detr(ig,l)=fm(ig,l) 183 entr(ig,l)=fm(ig,l+1) 184 endif 232 ! entr(ig,l)=fm(ig,l+1) 233 234 ! Dans le cas ou on est au dessus de la couche d'alimentation et que le 235 ! detrainement est plus fort que le flux de masse, on stope le thermique. 236 if (l.gt.lalim(ig)) then 237 lmax(ig)=l 238 fm(ig,l+1)=0. 239 entr(ig,l)=0. 240 else 241 ncorecfm7=ncorecfm7+1 242 endif 243 endif 244 245 if(l.gt.lmax(ig)) then 246 detr(ig,l)=0. 247 fm(ig,l+1)=0. 248 entr(ig,l)=0. 249 endif 250 185 251 if (entr(ig,l).lt.0.) then 186 252 print*,'ig,l,lmax(ig)',ig,l,lmax(ig) … … 194 260 195 261 262 if (lev_out.ge.10) & 263 & call printflux(ngrid,klev,lunout,igout,f,lmax,lalim, & 264 & ptimestep,masse,entr,detr,fm,'6 ') 265 196 266 !------------------------------------------------------------------------- 197 267 !fm ne peut pas etre negatif … … 207 277 endif 208 278 if (detr(ig,l).lt.0.) then 209 print*,' ig,l,lmax(ig)',ig,l,lmax(ig)279 print*,'cas 2 : ig,l,lmax(ig)',ig,l,lmax(ig) 210 280 print*,'detr(ig,l)',detr(ig,l) 211 281 print*,'fm(ig,l)',fm(ig,l) … … 215 285 enddo 216 286 287 if (lev_out.ge.10) & 288 & call printflux(ngrid,klev,lunout,igout,f,lmax,lalim, & 289 & ptimestep,masse,entr,detr,fm,'7 ') 290 217 291 !----------------------------------------------------------------------- 218 292 !la fraction couverte ne peut pas etre superieure a 1 219 293 !----------------------------------------------------------------------- 294 295 296 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 297 ! FH Partie a revisiter. 298 ! Il semble qu'etaient codees ici deux optiques dans le cas 299 ! F/ (rho *w) > 1 300 ! soit limiter la hauteur du thermique en considerant que c'est 301 ! la derniere chouche, soit limiter F a rho w. 302 ! Dans le second cas, il faut en fait limiter a un peu moins 303 ! que ca parce qu'on a des 1 / ( 1 -alpha) un peu plus loin 304 ! dans thermcell_main et qu'il semble de toutes facons deraisonable 305 ! d'avoir des fractions de 1.. 306 ! Ci dessous, et dans l'etat actuel, le premier des deux if est 307 ! sans doute inutile. 308 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 309 310 if(1.eq.0) then 311 220 312 do l=1,klev 221 313 do ig=1,ngrid 222 314 if (zw2(ig,l+1).gt.1.e-10) then 223 if ((((fm(ig,l+1))/(rhobarz(ig,l+1)*zw2(ig,l+1))).gt. &224 & 1.)) then315 zfm=rhobarz(ig,l+1)*zw2(ig,l+1)*alphamax 316 if ( fm(ig,l+1) .gt. zfm) then 225 317 f_old=fm(ig,l+1) 226 fm(ig,l+1)= rhobarz(ig,l+1)*zw2(ig,l+1)227 228 318 fm(ig,l+1)=zfm 319 ! zw2(ig,l+1)=0. 320 ! zqla(ig,l+1)=0. 229 321 detr(ig,l)=detr(ig,l)+f_old-fm(ig,l+1) 230 231 322 ! lmax(ig)=l+1 323 ! zmax(ig)=zlev(ig,lmax(ig)) 232 324 ! print*,'alpha>1',l+1,lmax(ig) 233 325 ncorecalpha=ncorecalpha+1 … … 237 329 enddo 238 330 ! 331 if (lev_out.ge.10) & 332 & call printflux(ngrid,klev,lunout,igout,f,lmax,lalim, & 333 & ptimestep,masse,entr,detr,fm,'8 ') 334 335 endif 239 336 240 337 !----------------------------------------------------------------------- … … 248 345 eee0=entr(ig,l) 249 346 ddd0=detr(ig,l) 250 eee=entr(ig,l)-masse(ig,l)*f racmax/ptimestep347 eee=entr(ig,l)-masse(ig,l)*fomass_max/ptimestep 251 348 ddd=detr(ig,l)-eee 252 349 if (eee.gt.0.) then … … 272 369 print*,'detr',detr(ig,l) 273 370 print*,'masse',masse(ig,l) 274 print*,'f racmax',fracmax275 print*,'masse(ig,l)*f racmax/ptimestep',masse(ig,l)*fracmax/ptimestep371 print*,'fomass_max',fomass_max 372 print*,'masse(ig,l)*fomass_max/ptimestep',masse(ig,l)*fomass_max/ptimestep 276 373 print*,'ptimestep',ptimestep 277 374 print*,'lmax(ig)',lmax(ig) … … 309 406 & ncorecfm4,'x fm4',ncorecfm5,'x fm5 et', & 310 407 & ncorecfm6,'x fm6', & 408 & ncorecfm7,'x fm7', & 409 & ncorecfm8,'x fm8', & 311 410 & ncorecalpha,'x alpha' 312 411 endif 313 412 413 if (lev_out.ge.10) & 414 & call printflux(ngrid,klev,lunout,igout,f,lmax,lalim, & 415 & ptimestep,masse,entr,detr,fm,'fin') 416 314 417 315 418 return 316 419 end 420 421 subroutine printflux(ngrid,klev,lunout,igout,f,lmax,lalim, & 422 & ptimestep,masse,entr,detr,fm,descr) 423 424 implicit none 425 426 integer ngrid,klev,lunout,igout,l,lm 427 428 integer lmax(klev),lalim(klev) 429 real ptimestep,masse(ngrid,klev),entr(ngrid,klev),detr(ngrid,klev) 430 real fm(ngrid,klev+1),f(ngrid) 431 432 character*3 descr 433 434 lm=lmax(igout)+5 435 if(lm.gt.klev) lm=klev 436 437 print*,'Impression jusque lm=',lm 438 439 write(lunout,*) 'Dans thermcell_flux '//descr 440 write(lunout,*) 'flux base ',f(igout) 441 write(lunout,*) 'lmax ',lmax(igout) 442 write(lunout,*) 'lalim ',lalim(igout) 443 write(lunout,*) 'ig= ',igout 444 write(lunout,'(a3,4a14)') 'l','M','E','D','F' 445 write(lunout,'(i4,4e14.4)') (l,masse(igout,l)/ptimestep, & 446 & entr(igout,l),detr(igout,l) & 447 & ,fm(igout,l+1),l=1,lm) 448 449 450 do l=lmax(igout)+1,klev 451 if (abs(entr(igout,l))+abs(detr(igout,l))+abs(fm(igout,l)).gt.0.) then 452 print*,'cas 1 : igout,l,lmax(igout)',igout,l,lmax(igout) 453 print*,'entr(igout,l)',entr(igout,l) 454 print*,'detr(igout,l)',detr(igout,l) 455 print*,'fm(igout,l)',fm(igout,l) 456 stop 457 endif 458 enddo 459 460 return 461 end 462 -
LMDZ4/trunk/libf/phytherm/thermcell_main.F90
r814 r852 1 ! 2 ! $Header$ 3 ! 1 4 SUBROUTINE thermcell_main(ngrid,nlay,ptimestep & 2 5 & ,pplay,pplev,pphi,debut & … … 55 58 ! ------ 56 59 57 integer,save :: igout=871 60 ! integer,save :: igout=4259 61 integer,save :: igout=2928 58 62 integer,save :: lunout=6 59 integer,save :: lev_out= 063 integer,save :: lev_out=10 60 64 61 65 INTEGER ig,k,l,ll … … 129 133 real zlevinter(klon) 130 134 logical debut 135 real seuil 131 136 132 137 ! … … 142 147 ! --------------- 143 148 ! 149 150 seuil=0.25 151 144 152 if (lev_out.ge.1) print*,'thermcell_main V4' 145 153 … … 228 236 229 237 !------------------------------------------------------------------ 230 ! Calcul de w2, carre de w a partir de la cape 231 ! 232 ! Indicages: 233 ! l'ascendance provenant du niveau k traverse l'interface l avec 234 ! une vitesse wa(k,l). 235 ! 236 ! -------------------- 237 ! 238 ! + + + + + + + + + + 239 ! 240 ! wa(k,l) ---- -------------------- l 241 ! /\ 242 ! /||\ + + + + + + + + + + 243 ! || 244 ! || -------------------- 245 ! || 246 ! || + + + + + + + + + + 247 ! || 248 ! || -------------------- 249 ! ||__ 250 ! |___ + + + + + + + + + + k 251 ! 252 ! -------------------- 253 ! 254 ! 255 ! 238 ! 239 ! /|\ 240 ! -------- | F_k+1 ------- 241 ! ----> D_k 242 ! /|\ <---- E_k , A_k 243 ! -------- | F_k --------- 244 ! ----> D_k-1 245 ! <---- E_k-1 , A_k-1 246 ! 247 ! 248 ! 249 ! 250 ! 251 ! --------------------------- 252 ! 253 ! ----- F_lmax+1=0 ---------- \ 254 ! lmax (zmax) | 255 ! --------------------------- | 256 ! | 257 ! --------------------------- | 258 ! | 259 ! --------------------------- | 260 ! | 261 ! --------------------------- | 262 ! | 263 ! --------------------------- | 264 ! | E 265 ! --------------------------- | D 266 ! | 267 ! --------------------------- | 268 ! | 269 ! --------------------------- \ | 270 ! lalim | | 271 ! --------------------------- | | 272 ! | | 273 ! --------------------------- | | 274 ! | A | 275 ! --------------------------- | | 276 ! | | 277 ! --------------------------- | | 278 ! lmin (=1 pour le moment) | | 279 ! ----- F_lmin=0 ------------ / / 280 ! 281 ! --------------------------- 282 ! ////////////////////////// 283 ! 284 ! 285 !============================================================================= 286 ! Calculs initiaux ne faisant pas intervenir les changements de phase 287 !============================================================================= 288 256 289 !------------------------------------------------------------------ 257 !definition du profil d alimentation a partir de la flottabilite: 258 !alim_star, alim_star_tot, lalim et lmin 290 ! 1. alim_star est le profil vertical de l'alimentation à la base du 291 ! panache thermique, calculé à partir de la flotabilité de l'air sec 292 ! 2. lmin et lalim sont les indices inferieurs et superieurs de alim_star 259 293 !------------------------------------------------------------------ 260 294 ! … … 262 296 CALL thermcell_init(ngrid,nlay,ztv,zlev, & 263 297 & lalim,lmin,alim_star,alim_star_tot,lev_out) 298 299 call test_ltherm(ngrid,nlay,pplev,pplay,lmin,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_init lmin ') 300 call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_init lalim ') 301 264 302 265 303 if (lev_out.ge.1) print*,'thermcell_main apres thermcell_init' … … 268 306 write(lunout,*) 'lmin ',lmin(igout) 269 307 write(lunout,*) 'lalim ',lalim(igout) 308 write(lunout,*) ' ig l alim_star thetav' 309 write(lunout,'(i6,i4,2e15.5)') (igout,l,alim_star(igout,l) & 310 & ,ztv(igout,l),l=1,lalim(igout)+4) 270 311 endif 271 312 272 !----------------------------------------------------------------------------------- 273 !calcul des caracteristiques du thermique sec pour le calcul de detr et la fermeture 274 !----------------------------------------------------------------------------------- 313 !----------------------------------------------------------------------------- 314 ! 3. wmax_sec et zmax_sec sont les vitesses et altitudes maximum d'un 315 ! panache sec conservatif (e=d=0) alimente selon alim_star 316 ! Il s'agit d'un calcul de type CAPE 317 ! zmax_sec est utilisé pour déterminer la géométrie du thermique. 318 !------------------------------------------------------------------------------ 275 319 ! 276 320 CALL thermcell_dry(ngrid,nlay,zlev,pphi,ztv,alim_star, & 277 321 & lalim,lmin,zmax_sec,wmax_sec,lev_out) 278 322 323 call test_ltherm(ngrid,nlay,pplev,pplay,lmin,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry lmin ') 324 call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry lalim ') 325 279 326 if (lev_out.ge.1) print*,'thermcell_main apres thermcell_dry' 327 if (lev_out.ge.10) then 328 write(lunout,*) 'Dans thermcell_main 1b' 329 write(lunout,*) 'lmin ',lmin(igout) 330 write(lunout,*) 'lalim ',lalim(igout) 331 write(lunout,*) ' ig l alim_star entr_star detr_star f_star ' 332 write(lunout,'(i6,i4,e15.5)') (igout,l,alim_star(igout,l) & 333 & ,l=1,lalim(igout)+4) 334 endif 280 335 281 336 … … 289 344 & lalim,zmax_sec,f0,detr_star,entr_star,f_star,ztva, & 290 345 & ztla,zqla,zqta,zha,zw2,zqsatth,lmix,linter,lev_out) 346 call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lalim ') 347 call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lmix ') 291 348 292 349 if (lev_out.ge.1) print*,'thermcell_main apres thermcell_plume' … … 297 354 write(lunout,*) ' ig l alim_star entr_star detr_star f_star ' 298 355 write(lunout,'(i6,i4,4e15.5)') (igout,l,alim_star(igout,l),entr_star(igout,l),detr_star(igout,l) & 299 & ,f_star(igout,l+1),l=1, lmax(igout))356 & ,f_star(igout,l+1),l=1,nint(linter(igout))+5) 300 357 endif 301 358 … … 307 364 & zlev,lmax,zmax,zmax0,zmix,wmax,lev_out) 308 365 366 367 call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lalim ') 368 call test_ltherm(ngrid,nlay,pplev,pplay,lmin ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmin ') 369 call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmix ') 370 call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmax ') 371 309 372 if (lev_out.ge.1) print*,'thermcell_main apres thermcell_height' 310 373 … … 322 385 !------------------------------------------------------------------------------- 323 386 324 CALL thermcell_flux (ngrid,nlay,ptimestep,masse, &387 CALL thermcell_flux2(ngrid,nlay,ptimestep,masse, & 325 388 & lalim,lmax,alim_star, & 326 389 & entr_star,detr_star,f,rhobarz,zlev,zw2,fm,entr, & … … 328 391 329 392 if (lev_out.ge.1) print*,'thermcell_main apres thermcell_flux' 393 call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lalim ') 394 call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lmax ') 330 395 331 396 !c------------------------------------------------------------------ … … 391 456 enddo 392 457 458 print*,'14a OK convect8' 393 459 ! calcul du niveau de condensation 394 460 ! initialisation … … 410 476 enddo 411 477 enddo 478 print*,'14b OK convect8' 412 479 do k=nlay,1,-1 413 480 do ig=1,ngrid … … 418 485 enddo 419 486 enddo 487 print*,'14c OK convect8' 420 488 !calcul des moments 421 489 !initialisation … … 429 497 enddo 430 498 enddo 499 print*,'14d OK convect8' 431 500 do l=1,nlay 432 501 do ig=1,ngrid … … 440 509 q2(ig,l)=zf2*(zqta(ig,l)*1000.-po(ig,l)*1000.)**2 441 510 !test: on calcul q2/po=ratqsc 442 ratqscth(ig,l)=sqrt( q2(ig,l))/(po(ig,l)*1000.)511 ratqscth(ig,l)=sqrt(max(q2(ig,l),1.e-6)/(po(ig,l)*1000.)) 443 512 enddo 444 513 enddo 445 514 !calcul du ratqscdiff 515 print*,'14e OK convect8' 446 516 var=0. 447 517 vardiff=0. … … 452 522 enddo 453 523 enddo 524 print*,'14f OK convect8' 454 525 do ig=1,ngrid 455 526 do l=1,lalim(ig) … … 462 533 enddo 463 534 enddo 535 print*,'14g OK convect8' 464 536 do l=1,nlay 465 537 do ig=1,ngrid … … 496 568 497 569 !----------------------------------------------------------------------------- 570 571 subroutine test_ltherm(klon,klev,pplev,pplay,long,seuil,ztv,po,ztva,zqla,f_star,zw2,comment) 572 573 integer klon,klev 574 real pplev(klon,klev+1),pplay(klon,klev) 575 real ztv(klon,klev) 576 real po(klon,klev) 577 real ztva(klon,klev) 578 real zqla(klon,klev) 579 real f_star(klon,klev) 580 real zw2(klon,klev) 581 integer long(klon) 582 real seuil 583 character*21 comment 584 585 print*,'TEST ',comment 586 ! test sur la hauteur des thermiques ... 587 do i=1,klon 588 if (pplay(i,long(i)).lt.seuil*pplev(i,1)) then 589 print*,'WARNING ',comment,' au point ',i,' K= ',long(i) 590 print*,' K P(MB) THV(K) Qenv(g/kg)THVA QLA(g/kg) F* W2' 591 do k=1,klev 592 write(6,'(i3,7f10.3)') k,pplay(i,k),ztv(i,k),1000*po(i,k),ztva(i,k),1000*zqla(i,k),f_star(i,k),zw2(i,k) 593 enddo 594 ! stop 595 endif 596 enddo 597 598 599 return 600 end 601 -
LMDZ4/trunk/libf/phytherm/thermcell_plume.F90
r814 r852 1 1 SUBROUTINE thermcell_plume(ngrid,klev,ztv,zthl,po,zl,rhobarz, & 2 2 & zlev,pplev,pphi,zpspsk,l_mix,r_aspect,alim_star, & 3 & l entr,zmax_sec,f0,detr_star,entr_star,f_star,ztva, &3 & lalim,zmax_sec,f0,detr_star,entr_star,f_star,ztva, & 4 4 & ztla,zqla,zqta,zha,zw2,zqsatth,lmix,linter,lev_out) 5 5 … … 29 29 REAL l_mix 30 30 REAL r_aspect 31 INTEGER l entr(ngrid)31 INTEGER lalim(ngrid) 32 32 integer lev_out ! niveau pour les print 33 33 … … 133 133 !tests sur la definition du detr 134 134 !calcul de detr_star et entr_star 135 w_est(ig,3)=zw2(ig,2)* & 136 & ((f_star(ig,2))**2) & 137 & /(f_star(ig,2)+alim_star(ig,2))**2+ & 138 & 2.*RG*(ztva(ig,1)-ztv(ig,2))/ztv(ig,2) & 139 & *(zlev(ig,3)-zlev(ig,2)) 135 136 137 138 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 139 ! FH le test miraculeux de Catherine ? Le bout du tunel ? 140 ! w_est(ig,3)=zw2(ig,2)* & 141 ! & ((f_star(ig,2))**2) & 142 ! & /(f_star(ig,2)+alim_star(ig,2))**2+ & 143 ! & 2.*RG*(ztva(ig,1)-ztv(ig,2))/ztv(ig,2) & 144 ! & *(zlev(ig,3)-zlev(ig,2)) 145 ! if (l.gt.2) then 146 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 147 140 148 if (l.gt.2) then 149 w_est(ig,3)=zw2(ig,2)* & 150 & ((f_star(ig,2))**2) & 151 & /(f_star(ig,2)+alim_star(ig,2))**2+ & 152 & 2.*RG*(ztva(ig,2)-ztv(ig,2))/ztv(ig,2) & 153 & *(zlev(ig,3)-zlev(ig,2)) 154 155 141 156 !--------------------------------------------------------------------------- 142 157 !calcul de l entrainement et du detrainement lateral … … 217 232 !calcul de entr_star 218 233 ! 219 if (detr_star(ig,l).gt.f_star(ig,l)) then220 detr_star(ig,l)=f_star(ig,l)221 !a decommenter ou pas?222 ! entr_star(ig,l)=0.223 endif224 225 234 ! Deplacement du calcul de entr_star pour eviter d'avoir aussi 226 235 ! entr_star > fstar. 236 ! Redeplacer suite a la transformation du cas detr>f 227 237 ! FH 228 entr_star(ig,l)=0.4*detr_star(ig,l) 229 ! 238 if(l.gt.lalim(ig)) then 239 entr_star(ig,l)=0.4*detr_star(ig,l) 240 else 241 242 ! FH : 243 ! Cette ligne doit permettre de garantir qu'on a toujours un flux = 1 244 ! en haut de la couche d'alimentation. 245 ! A remettre en questoin a la premiere occasion mais ca peut aider a 246 ! ecrire un code robuste. 247 ! Que ce soit avec ca ou avec l'ancienne facon de faire (e* = 0 mais 248 ! d* non nul) on a une discontinuité de e* ou d* en haut de la couche 249 ! d'alimentation, ce qui n'est pas forcement heureux. 250 entr_star(ig,l)=max(detr_star(ig,l)-alim_star(ig,l),0.) 251 detr_star(ig,l)=entr_star(ig,l) 252 endif 253 254 ! 255 if (detr_star(ig,l).gt.f_star(ig,l)) then 256 257 ! Ce test est là entre autres parce qu'on passe par des valeurs 258 ! delirantes de detr_star. 259 ! ca vaut sans doute le coup de verifier pourquoi. 260 261 detr_star(ig,l)=f_star(ig,l) 262 if (l.gt.lalim(ig)+1) then 263 entr_star(ig,l)=0. 264 alim_star(ig,l)=0. 265 ! FH ajout pour forcer a stoper le thermique juste sous le sommet 266 ! de la couche (voir calcul de finter) 267 zw2(ig,l+1)=-1.e-10 268 linter(ig)=l+1 269 else 270 entr_star(ig,l)=detr_star(ig,l) 271 endif 272 endif 273 230 274 else 231 275 detr_star(ig,l)=0. 232 276 entr_star(ig,l)=0. 233 277 endif 278 279 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 280 ! FH inutile si on conserve comme on l'a fait plus haut entr=detr 281 ! dans la couche d'alimentation 234 282 !pas d entrainement dans la couche alim 235 if ((l.lt.lentr(ig))) then 236 entr_star(ig,l)=0. 237 endif 283 ! if ((l.le.lalim(ig))) then 284 ! entr_star(ig,l)=0. 285 ! endif 286 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 238 287 ! 239 288 !prise en compte du detrainement et de l entrainement dans le calcul du flux … … 308 357 ! 309 358 !initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max 359 360 if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then 361 ! stop'On tombe sur le cas particulier de thermcell_dry' 362 print*,'On tombe sur le cas particulier de thermcell_plume' 363 zw2(ig,l+1)=0. 364 linter(ig)=l+1 365 endif 366 367 310 368 if (zw2(ig,l+1).lt.0.) then 311 369 linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l)) & 312 370 & -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l)) 313 371 zw2(ig,l+1)=0. 314 else 372 endif 373 315 374 wa_moy(ig,l+1)=sqrt(zw2(ig,l+1)) 316 endif 375 317 376 if (wa_moy(ig,l+1).gt.wmaxa(ig)) then 318 377 ! lmix est le niveau de la couche ou w (wa_moy) est maximum -
LMDZ4/trunk/libf/phytherm/write_histday.h
r814 r852 491 491 $ zx_tmp_fi2d) 492 492 C 493 ! FH Sorties specifiques pour Mellor et Yamada 494 if (iflag_pbl>1 .and. lev_histday.gt.10 ) then 495 496 CALL histwrite_phy(nid_day,"tke_"//clnsurf(nsrf),itau_w, 497 $ pbl_tke(:,1:klev,nsrf)) 498 499 CALL histwrite_phy(nid_day,"tke_max_"//clnsurf(nsrf),itau_w, 500 $ pbl_tke(:,1:klev,nsrf)) 501 endif 502 493 503 END DO 494 504 c================================================================= -
LMDZ4/trunk/libf/phytherm/write_histmth.h
r814 r852 681 681 IF(lev_histmth.GE.4) THEN 682 682 c 683 c FH Sorties pour la couche limite 684 CALL histwrite_phy(nid_mth,"kz",itau_w,ycoefh) 685 CALL histwrite_phy(nid_mth,"kz_max",itau_w,ycoefh) 686 687 if(iflag_pbl>1) then 688 zx_tmp_fi3d=0. 689 do nsrf=1,nbsrf 690 do k=1,klev 691 zx_tmp_fi3d(:,k)=zx_tmp_fi3d(:,k) 692 , +pctsrf(:,nsrf)*pbl_tke(:,k,nsrf) 693 enddo 694 enddo 695 CALL histwrite_phy(nid_mth,"tke",itau_w,zx_tmp_fi3d) 696 CALL histwrite_phy(nid_mth,"tke_max",itau_w,zx_tmp_fi3d) 697 endif 698 699 700 683 701 cym CALL gr_fi_ecrit(klev,klon,iim,jjmp1, clwcon0, zx_tmp_3d) 684 702 CALL histwrite_phy(nid_mth,"clwcon",itau_w,clwcon0) … … 757 775 c 758 776 cIM: 101003 : K/30min ==> K/s 777 zx_tmp_fi3d(1:klon,1:klev)=d_t_ajsb(1:klon,1:klev)/pdtphys 778 cym CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d) 779 CALL histwrite_phy(nid_mth,"dtajs",itau_w,zx_tmp_fi3d) 780 c 781 zx_tmp_fi3d(1:klon,1:klev)=d_q_ajsb(1:klon,1:klev)/pdtphys 782 cym CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d) 783 CALL histwrite_phy(nid_mth,"dqajs",itau_w,zx_tmp_fi3d) 784 c 785 cIM: 101003 : K/30min ==> K/s 759 786 zx_tmp_fi3d(1:klon,1:klev)=d_t_ajs(1:klon,1:klev)/pdtphys 760 787 cym CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d) 761 CALL histwrite_phy(nid_mth,"dt ajs",itau_w,zx_tmp_fi3d)788 CALL histwrite_phy(nid_mth,"dtthe",itau_w,zx_tmp_fi3d) 762 789 c 763 790 zx_tmp_fi3d(1:klon,1:klev)=d_q_ajs(1:klon,1:klev)/pdtphys 764 791 cym CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d) 765 CALL histwrite_phy(nid_mth,"dq ajs",itau_w,zx_tmp_fi3d)792 CALL histwrite_phy(nid_mth,"dqthe",itau_w,zx_tmp_fi3d) 766 793 c 767 794 cIM: 101003 : K/day ==> K/s … … 933 960 CALL histwrite_phy(nid_mth,"dtcon",itau_w,zx_tmp_fi3d) 934 961 c 935 c temperature tendency due to dry convective processes 936 DO l=1, klev 937 DO i=1, klon 938 zx_tmp_fi3d(i,l)=d_t_ajs(i,l)/pdtphys 939 ENDDO !i 940 ENDDO !l 941 c 942 cym CALL gr_fi_ecrit(klev, klon,iim,jjmp1, zx_tmp_fi3d,zx_tmp_3d) 962 963 964 965 cIM: 101003 : K/30min ==> K/s 966 zx_tmp_fi3d(1:klon,1:klev)=d_t_ajsb(1:klon,1:klev)/pdtphys 967 cym CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d) 943 968 CALL histwrite_phy(nid_mth,"dtajs",itau_w,zx_tmp_fi3d) 969 c 970 zx_tmp_fi3d(1:klon,1:klev)=d_q_ajsb(1:klon,1:klev)/pdtphys 971 cym CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d) 972 CALL histwrite_phy(nid_mth,"dqajs",itau_w,zx_tmp_fi3d) 973 c 974 cIM: 101003 : K/30min ==> K/s 975 zx_tmp_fi3d(1:klon,1:klev)=d_t_ajs(1:klon,1:klev)/pdtphys 976 cym CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d) 977 CALL histwrite_phy(nid_mth,"dtthe",itau_w,zx_tmp_fi3d) 978 c 979 zx_tmp_fi3d(1:klon,1:klev)=d_q_ajs(1:klon,1:klev)/pdtphys 980 cym CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d) 981 CALL histwrite_phy(nid_mth,"dqthe",itau_w,zx_tmp_fi3d) 982 c 983 944 984 c 945 985 c temperature tendency due to large scale precipitation … … 1672 1712 IF(lev_histmth.GE.4) THEN 1673 1713 c 1714 c FH Sorties pour la couche limite 1715 CALL histwrite_phy(nid_mth,"kz",itau_w,ycoefh) 1716 CALL histwrite_phy(nid_mth,"kz_max",itau_w,ycoefh) 1717 1718 if(iflag_pbl>1) then 1719 zx_tmp_fi3d=0. 1720 do nsrf=1,nbsrf 1721 do k=1,klev 1722 zx_tmp_fi3d(:,k)=zx_tmp_fi3d(:,k) 1723 , +pctsrf(:,nsrf)*pbl_tke(:,k,nsrf) 1724 enddo 1725 enddo 1726 CALL histwrite_phy(nid_mth,"tke",itau_w,zx_tmp_fi3d) 1727 CALL histwrite_phy(nid_mth,"tke_max",itau_w,zx_tmp_fi3d) 1728 endif 1729 1730 1731 1674 1732 cym CALL gr_fi_ecrit(klev,klon,iim,jjmp1, clwcon0, zx_tmp_3d) 1675 1733 CALL histwrite_phy(nid_mth,"clwcon",itau_w,clwcon0) … … 1739 1797 CALL histwrite_phy(nid_mth,"ratqs",itau_w,ratqs) 1740 1798 c 1741 zx_tmp_fi3d(1:klon,1:klev)=d_q_ajs(1:klon,1:klev)/pdtphys1742 cym CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)1743 CALL histwrite_phy(nid_mth,"dqajs",itau_w,zx_tmp_fi3d)1744 c1745 1799 cIM: 101003 : K/day ==> K/s 1746 1800 zx_tmp_fi3d(1:klon,1:klev)=heat(1:klon,1:klev)/RDAY -
LMDZ4/trunk/libf/phytherm/yamada4.F
r814 r852 458 458 c Diagnostique pour stokage 459 459 460 if(1.eq.0)then 460 461 rino=rif 461 smyam(:,1:klev)=sm(:,1:klev) 462 styam=sm(:,1:klev)*alpha(:,1:klev) 463 lyam(1:klon,1:klev)=l(:,1:klev) 464 knyam(1:klon,1:klev)=kn(:,1:klev) 462 smyam(1:ngrid,1)=0. 463 styam(1:ngrid,1)=0. 464 lyam(1:ngrid,1)=0. 465 knyam(1:ngrid,1)=0. 466 w2yam(1:ngrid,1)=0. 467 t2yam(1:ngrid,1)=0. 468 469 smyam(1:ngrid,2:klev)=sm(1:ngrid,2:klev) 470 styam(1:ngrid,2:klev)=sm(1:ngrid,2:klev)*alpha(1:ngrid,2:klev) 471 lyam(1:ngrid,2:klev)=l(1:ngrid,2:klev) 472 knyam(1:ngrid,2:klev)=kn(1:ngrid,2:klev) 465 473 466 474 c Estimations de w'2 et T'2 d'apres Abdela et McFarlane 467 475 468 if(1.eq.0)then469 w2yam=q2(:,1:klev)*0.24470 s +lyam(:,1:klev)*5.17*kn(:,1:klev)*n2(:,1:klev)471 s /sqrt(q2(:,1:klev)) 472 473 t2yam=9.1*kn(:,1:klev)*dtetadz(:,1:klev)**2/sqrt(q2(:,1:klev))474 s *lyam(:,1:klev)475 476 w2yam(1:ngrid,2:klev)=q2(1:ngrid,2:klev)*0.24 477 s +lyam(1:ngrid,2:klev)*5.17*kn(1:ngrid,2:klev) 478 s *n2(1:ngrid,2:klev)/sqrt(q2(1:ngrid,2:klev)) 479 480 t2yam(1:ngrid,2:klev)=9.1*kn(1:ngrid,2:klev) 481 s *dtetadz(1:ngrid,2:klev)**2 482 s /sqrt(q2(1:ngrid,2:klev))*lyam(1:ngrid,2:klev) 483 endif 476 484 477 485 c print*,'OKFIN'
Note: See TracChangeset
for help on using the changeset viewer.
