Changeset 852 for LMDZ4/trunk/libf


Ignore:
Timestamp:
Oct 11, 2007, 3:43:42 PM (17 years ago)
Author:
Laurent Fairhead
Message:

Mise a jour de la physique avec thermiques avec la version de FH d'aout 2007
LF

Location:
LMDZ4/trunk/libf/phytherm
Files:
17 edited

Legend:

Unmodified
Added
Removed

  • LMDZ4/trunk/libf/phytherm/ajsec.F

    r814 r852  
    22! $Header$
    33!
    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
     7c======================================================================
     8c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
     9c Objet: ajustement sec (adaptation du GCM du LMD)
     10c======================================================================
     11c Arguments:
     12c t-------input-R- Temperature
     13c
     14c d_t-----output-R-Incrementation de la temperature
     15c======================================================================
     16cym#include "dimensions.h"
     17cym#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)
     22c
     23      INTEGER limbas(klon), limhau ! les couches a ajuster
     24c
     25      LOGICAL mixq
     26ccc      PARAMETER (mixq=.TRUE.)
     27      PARAMETER (mixq=.FALSE.)
     28c
     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
     37c
     38c Initialisation:
     39c
     40cym
     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
     49c------------------------------------- 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
     58c
     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
     64c
     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
     75c------------------------------------- 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
     107c------------ 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
     118c
     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
     125c
     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
     137c
     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
     146c
     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
     154c
     155      RETURN
     156      END
     157
     158      SUBROUTINE ajsec_convV2(paprs, pplay, t,q, d_t,d_q)
    5159      USE dimphy
    6160      IMPLICIT none

  • LMDZ4/trunk/libf/phytherm/clouds_gno.F

    r814 r852  
    219219c -- test numerical convergence:
    220220
    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
    222234          if ( ABS(dist(i)/beta(i)) .LT. epsilon ) then
    223235c           print*,'v-u **2',(v(i)-u(i))**2

  • LMDZ4/trunk/libf/phytherm/coef_diff_turb_mod.F90

    r814 r852  
    1818  SUBROUTINE coef_diff_turb(dtime, nsrf, knon, ni, &
    1919       ypaprs, ypplay, yu, yv, yq, yt, yts, yrugos, yqsurf, &
    20        ycoefm, ycoefh)
     20       ycoefm, ycoefh ,yq2)
    2121!
    2222! Calculate coefficients(ycoefm, ycoefh) for turbulent diffusion in the
    … …  
    3535    REAL, DIMENSION(klon,klev), INTENT(IN)     :: yq, yt
    3636    REAL, DIMENSION(klon), INTENT(IN)          :: yts, yrugos, yqsurf
     37    REAL, DIMENSION(klon,klev+1)               :: yq2
    3738 
    3839! Output arguments
    … …  
    4142    REAL, DIMENSION(klon,klev), INTENT(OUT)    :: ycoefm
    4243
    43 ! Local variables with attribute SAVE
    44 !****************************************************************************************
    45     REAL, DIMENSION(:,:,:), ALLOCATABLE, SAVE  :: q2save
    46     !$OMP THREADPRIVATE(q2save)
    47     LOGICAL, SAVE                              :: first_call=.TRUE.
    48     !$OMP THREADPRIVATE(first_call)   
    49 
    5044! Other local variables
    5145!****************************************************************************************
    5246    INTEGER                                    :: k, i, j
    5347    REAL, DIMENSION(klon,klev)                 :: ycoefm0, ycoefh0, yzlay, yteta
    54     REAL, DIMENSION(klon,klev+1)               :: yzlev, yq2, q2diag, ykmm, ykmn, ykmq
     48    REAL, DIMENSION(klon,klev+1)               :: yzlev, q2diag, ykmm, ykmn, ykmq
    5549    REAL, DIMENSION(klon)                      :: y_cd_h, y_cd_m
    5650    REAL, DIMENSION(klon)                      :: yustar
    … …  
    9488            ycoefm0, ycoefh0)
    9589
    96        DO k = 1, klev
     90       DO k = 2, klev
    9791          DO i = 1, knon
    9892             ycoefm(i,k) = MAX(ycoefm(i,k),ycoefm0(i,k))
    … …  
    121115            ycoefm0,ycoefh0)
    122116       
    123        DO k = 1, klev
     117       DO k = 2, klev
    124118          DO i = 1, knon
    125119             ycoefm(i,k) = MAX(ycoefm(i,k),ycoefm0(i,k))
    … …  
    137131
    138132    IF (iflag_pbl.GE.3) THEN
    139 
    140        IF (first_call) THEN
    141           PRINT*, "Using method MELLOR&YAMADA"
    142           ! NB! q2save could/should be read and written from (re)startphy.nc
    143           ALLOCATE(q2save(klon,klev+1,nbsrf))
    144           q2save(:,:,:) = 1.e-8
    145           first_call=.FALSE.
    146        END IF
    147133
    148134       yzlay(1:knon,1)= &
    … …  
    173159       END DO
    174160
    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)
    192170       
    193171       CALL ustarhb(knon,yu,yv,y_cd_m, yustar)
    … …  
    215193       ycoefh(1:knon,2:klev)=ykmn(1:knon,2:klev)
    216194               
    217 ! update q2save with yq2
    218        DO j=1,knon
    219           DO k=1,klev+1
    220              i=ni(j)
    221              q2save(i,k,nsrf) = yq2(j,k)
    222           END DO
    223        END DO
    224        
    225195    ENDIF !(iflag_pbl.ge.3)
    226196

  • LMDZ4/trunk/libf/phytherm/ini_histday.h

    r814 r852  
    5555         CALL histdef(nid_day, "weakinv", "Weak inversion", "",
    5656     s           iim,jjmp1,nhori, 1,1,1, -99, 32,
    57      s           "ave(X)", zsto,zout)
     57     s           "ave(X)", zstophy,zout)
    5858
    5959         CALL histdef(nid_day, "dthmin", "dTheta mini", "K/m",
    6060     s           iim,jjmp1,nhori, 1,1,1, -99, 32,
    61      s           "ave(X)", zsto,zout)
     61     s           "ave(X)", zstophy,zout)
    6262
    6363
    … …  
    547547     $         "ave(X)", zstophy,zout)
    548548C
     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
     564C
    549565         END DO
    550566C           

  • LMDZ4/trunk/libf/phytherm/ini_histmth.h

    r841 r852  
    285285c
    286286         CALL histdef(nid_mth, "sens", "Sensible heat flux", "W/m2",
    287      .                iim,jjmp1,nhori, 1,1,1, -99, 32,
     287     .                iim,jj_nb,nhori, 1,1,1, -99, 32,
    288288     .                "ave(X)", zstophy,zout)
    289289c
    … …  
    310310         CALL histdef(nid_mth, "fqfonte","Land ice melt",
    311311     .                "kg/m2/s",iim,jj_nb,nhori, 1,1,1, -99, 32,
    312      .                "ave(X)", zsto,zout)
     312     .                "ave(X)", zstophy,zout)
    313313
    314314         DO nsrf = 1, nbsrf
    … …  
    700700      IF(lev_histmth.GE.4) THEN
    701701c
     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)
     707c
     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
     712c
     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)
     716c
     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
    702722         CALL histdef(nid_mth, "clwcon",
    703723     .                "Convective Cloud Liquid water content"
    … …  
    730750     .                "ave(X)", zstophy,zout)
    731751c
     752
     753
     754c
    732755c        CALL histdef(nid_mth, "ducon", "Convection du", "m/s2",
    733756c    .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
    … …  
    774797     .                "ave(X)", zstophy,zout)
    775798
     799c
     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)
    776807c
    777808         CALL histdef(nid_mth, "dtajs", "Dry adjust. dT", "K/s",
    … …  
    914945c
    915946          CALL histdef(nid_mth, "dtcon","dtcon","K/s",
     947     .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
     948     .                "ave(X)", zstophy,zout)
     949c
     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",
    916955     .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
    917956     .                "ave(X)", zstophy,zout)
    … …  
    16221661      IF(lev_histmth.GE.4) THEN
    16231662c
     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)
     1668c
     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
     1673c
     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)
     1677c
     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
    16241685         CALL histdef(nid_mth, "clwcon",
    16251686     .                "Convective Cloud Liquid water content"

  • LMDZ4/trunk/libf/phytherm/pbl_surface_mod.F90

    r815 r852  
    197197       wfbils,    wfbilo,    flux_t,   flux_u, flux_v,&
    198198       dflux_t,   dflux_q,   zxsnow,                  &
    199        zxfluxt,   zxfluxq,   q2m,      flux_q )
     199       zxfluxt,   zxfluxq,   q2m,      flux_q, tke    )
    200200!****************************************************************************************
    201201! Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
    … …  
    240240! flux_t---output-R- flux de chaleur sensible (CpT) J/m**2/s (W/m**2)
    241241!                    (orientation positive vers le bas)
     242! tke---input/output-R- tke (kg/m**2/s)
    242243! flux_q---output-R- flux de vapeur d'eau (kg/m**2/s)
    243244! flux_u---output-R- tension du vent X: (kg m/s)/(m**2 s) ou Pascal
    … …  
    348349    REAL, DIMENSION(klon, nbsrf),INTENT(OUT)        :: q2m
    349350    REAL, DIMENSION(klon, klev, nbsrf), INTENT(OUT) :: flux_q
     351
     352! Input/output
     353    REAL, DIMENSION(klon, klev+1, nbsrf), INTENT(INOUT) :: tke
    350354
    351355
    … …  
    421425    REAL, DIMENSION(klon,klev)         :: delp
    422426    REAL, DIMENSION(klon,klev+1)       :: ypaprs
     427    REAL, DIMENSION(klon,klev+1)       :: ytke
    423428    REAL, DIMENSION(klon,nsoilmx)      :: ytsoil
    424429    REAL, DIMENSION(klon,nbsrf)        :: pctsrf_pot
    … …  
    451456
    452457
     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   
    453463!****************************************************************************************
    454464! End of declarations
    … …  
    514524    flux_u = 0.0  ; flux_v = 0.0     ; d_t = 0.0       ; d_q = 0.0     
    515525    d_u = 0.0     ; d_v = 0.0        ; zcoefh = 0.0    ; yqsol = 0.0   
    516     ytherm = 0.0
     526    ytherm = 0.0  ; ytke=0.
    517527     
    518528    ytsoil = 999999.
    … …  
    657667             ypplay(j,k) = pplay(i,k)
    658668             ydelp(j,k) = delp(i,k)
     669             ytke(j,k)=tke(i,k,nsrf)
    659670             yu(j,k) = u(i,k)
    660671             yv(j,k) = v(i,k)
    … …  
    687698       CALL coef_diff_turb(dtime, nsrf, knon, ni,  &
    688699            ypaprs, ypplay, yu, yv, yq, yt, yts, yrugos, yqsurf,  &
    689             ycoefm, ycoefh)
     700            ycoefm, ycoefh,ytke)
    690701       
    691702!****************************************************************************************
    … …  
    839850!****************************************************************************************
    840851
     852       tke(:,:,nsrf)=0.
    841853       DO k = 1, klev
    842854          DO j = 1, knon
    … …  
    853865             flux_u(i,k,nsrf) = y_flux_u(j,k)
    854866             flux_v(i,k,nsrf) = y_flux_v(j,k)
     867
     868             tke(i,k,nsrf)=ytke(j,k)
     869
    855870          ENDDO
    856871       ENDDO

  • LMDZ4/trunk/libf/phytherm/phyetat0.F

    r814 r852  
    1313     .           rugsrel_p,tabcntr0,
    1414     .           t_ancien_p,q_ancien_p,ancien_ok_p, rnebcon_p, ratqs_p,
    15      .           clwcon_p)
     15     .           clwcon_p,pbl_tke_p)
    1616
    1717      USE dimphy
    … …  
    3737#include "clesphys.h"
    3838#include "temps.h"
     39#include "thermcell.h"
     40#include "compbl.h"
    3941c======================================================================
    4042      CHARACTER*(*) fichnom
    … …  
    4547      REAL solaire_etat0
    4648      REAL tsol_p(klon,nbsrf)
     49      REAL pbl_tke_p(klon,klev,nbsrf)
    4750      REAL tsoil_p(klon,nsoilmx,nbsrf)
    4851      REAL tslab_p(klon), seaice_p(klon)
    … …  
    8386      REAL rlat(klon_glo), rlon(klon_glo)
    8487      REAL tsol(klon_glo,nbsrf)
     88      REAL pbl_tke(klon_glo,klev,nbsrf)
    8589      REAL tsoil(klon_glo,nsoilmx,nbsrf)
    8690cIM "slab" ocean
    … …  
    127131c
    128132      INTEGER nid, nvarid
    129       INTEGER ierr, i, nsrf, isoil
     133      INTEGER ierr, i, nsrf, isoil ,k
    130134      INTEGER length
    131135      PARAMETER (length=100)
    … …  
    14921496      xmax = MAXval(run_off_lic_0)
    14931497      PRINT*,'(ecart-type) run_off_lic_0:', xmin, xmax
     1498
     1499
     1500c Lecture de l'energie cinetique turbulente
     1501c
     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
    14941540c
    14951541c Fermer le fichier:
    … …  
    15161562      call Scatter( rlon,rlon_p)
    15171563      call Scatter( tsol,tsol_p)
     1564      IF (iflag_pbl>1) then
     1565         call Scatter( pbl_tke,pbl_tke_p)
     1566      endif
    15181567      call Scatter( tsoil,tsoil_p)
    15191568      call Scatter( tslab,tslab_p)

  • LMDZ4/trunk/libf/phytherm/phyredem.F

    r814 r852  
    99     .           radsol_p,zmea_p,zstd_p,zsig_p,
    1010     .           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)
    1213
    1314      USE dimphy
    … …  
    2930#include "control.h"
    3031#include "temps.h"
     32#include "thermcell.h"
     33#include "compbl.h"
    3134c======================================================================
    3235      CHARACTER*(*) fichnom
    … …  
    3538      REAL rlat_p(klon), rlon_p(klon)
    3639      REAL tsol_p(klon,nbsrf)
     40      REAL pbl_tke_p(klon,klev,nbsrf)
    3741      REAL tsoil_p(klon,nsoilmx,nbsrf)
    3842      CHARACTER*6 ocean
    … …  
    7074      REAL rlat(klon_glo), rlon(klon_glo)
    7175      REAL tsol(klon_glo,nbsrf)
     76      REAL pbl_tke(klon_glo,klev,nbsrf)
    7277      REAL tsoil(klon_glo,nsoilmx,nbsrf)
    7378      REAL tslab(klon_glo), seaice(klon_glo)
    … …  
    135140      call Gather( rlon_p,rlon)
    136141      call Gather( tsol_p,tsol)
     142      call Gather( pbl_tke_p,pbl_tke)
    137143      call Gather( tsoil_p,tsoil)
    138144      call Gather( tslab_p,tslab)
    … …  
    899905c
    900906c
     907!!!!!!!!!!!!!!!!!!!! DEB TKE PBL !!!!!!!!!!!!!!!!!!!!!!!!!
     908c
     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 !!!!!!!!!!!!!!!!!!!!!!!!!
     937c
    901938      ierr = NF_CLOSE(nid)
    902939c

  • LMDZ4/trunk/libf/phytherm/physiq.F

    r839 r852  
    189189
    190190      integer lmax_th(klon)
     191      integer limbas(klon)
    191192      real ratqscth(klon,klev)
    192193      real ratqsdiff(klon,klev)
    … …  
    945946      REAL fluxu(klon,klev, nbsrf)   ! flux turbulent de vitesse u
    946947      REAL fluxv(klon,klev, nbsrf)   ! flux turbulent de vitesse v
     948
     949      REAL, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: pbl_tke ! turb kinetic energy
     950c   !$OMP THREADPRIVATE(pbl_tke)
     951
    947952c
    948953      REAL zxfluxt(klon, klev)
    … …  
    10661071c$OMP THREADPRIVATE(d_u_con,d_v_con)
    10671072      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)
    10681074      REAL d_t_ajs(klon,klev), d_q_ajs(klon,klev)
    10691075      REAL d_u_ajs(klon,klev), d_v_ajs(klon,klev)
    … …  
    15541560         itaprad = 0
    15551561
     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
    15561574         CALL phyetat0 ("startphy.nc",dtime,co2_ppm_etat0,solaire_etat0,
    15571575     .       rlat,rlon,pctsrf, ftsol,
    … …  
    15611579     .       radsol,clesphy0,
    15621580     .       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
    15641589
    15651590       DO i=1,klon
    … …  
    16431668         ENDIF
    16441669
     1670           rugoro=0.
    16451671c34EK
    16461672         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
    16501686           CALL SUGWD(klon,klev,paprs,pplay)
    16511687           DO i=1,klon
    … …  
    20612097     d     wfbils,    wfbilo,    fluxt,   fluxu,  fluxv,
    20622098     -     dsens,     devap,     zxsnow,
    2063      -     zxfluxt,   zxfluxq,   q2m,     fluxq )
     2099     -     zxfluxt,   zxfluxq,   q2m,     fluxq, pbl_tke )
    20642100c
    20652101c
    … …  
    22332269          DO i = 1, klon
    22342270            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
    22352275          ENDDO
    22362276          DO i = 1, klon
    … …  
    23212361
    23222362
     2363      d_t_ajsb(:,:)=0.
     2364      d_q_ajsb(:,:)=0.
    23232365      d_t_ajs(:,:)=0.
    23242366      d_u_ajs(:,:)=0.
    … …  
    23362378c  ====
    23372379         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
    23462382      else
     2383
    23472384c  Thermiques
    23482385c  ==========
    … …  
    23512388         print*,'JUSTE AVANT , iflag_thermals='
    23522389     s   ,iflag_thermals,'   nsplit_thermals=',nsplit_thermals
     2390
     2391
     2392         if (iflag_thermals.gt.1) then
    23532393         call calltherm(pdtphys
    23542394     s      ,pplay,paprs,pphi,weak_inversion
    … …  
    23572397     s      ,fm_therm,entr_therm,zqasc,clwcon0th,lmax_th,ratqscth,
    23582398     s       ratqsdiff,zqsatth)
     2399         endif
    23592400
    23602401!        call calltherm(pdtphys
    … …  
    23632404!    s      ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs
    23642405!    s      ,fm_therm,entr_therm)
     2406
     2407c  Ajustement sec
     2408c  ==============
     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
    23652444      endif
    23662445c
    … …  
    23972476         ptconvth(:,:)=.false.
    23982477         ratqsc(:,:)=0.
    2399          print*,'avant clouds_gno'
     2478         print*,'avant clouds_gno thermique'
    24002479         call clouds_gno
    24012480     s   (klon,klev,q_seri,zqsat,clwcon0th,ptconvth,ratqsc,rnebcon0th)
    … …  
    33213400     .      radsol,
    33223401     .      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)
    33243404      ENDIF
    33253405     

  • LMDZ4/trunk/libf/phytherm/thermcell_closure.F90

    r814 r852  
    1818      REAL zmax(ngrid),zmax_sec(ngrid)
    1919      REAL wmax(ngrid),wmax_sec(ngrid)
     20      real zdenom
    2021
    2122      REAL alim_star2(ngrid)
    … …  
    3536     &                    /(rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k)))
    3637             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
    3750             if ((zmax_sec(ig).gt.1.e-10).and.(1.eq.1)) then
    3851             f(ig)=wmax_sec(ig)/(max(500.,zmax_sec(ig))*r_aspect  &
    … …  
    4154     &                     zmax_sec(ig))*wmax_sec(ig))
    4255             else
    43              f(ig)=wmax(ig)/(max(500.,zmax(ig))*r_aspect*alim_star2(ig))
     56             f(ig)=wmax(ig)/zdenom
    4457            f(ig)=f(ig)+(f0(ig)-f(ig))*exp((-ptimestep/  &
    4558     &                     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      &                            lentr,lmin,zmax,wmax,lev_out)
     1       SUBROUTINE thermcell_dry(ngrid,nlay,zlev,pphi,ztv,alim_star,  &
     2     &                            lalim,lmin,zmax,wmax,lev_out)
    33
    44!--------------------------------------------------------------------------
    … …  
    1313       REAL pphi(ngrid,nlay)
    1414       REAl ztv(ngrid,nlay)
    15        REAL entr_star(ngrid,nlay)
    16        INTEGER lentr(ngrid)
     15       REAL alim_star(ngrid,nlay)
     16       INTEGER lalim(ngrid)
    1717      integer lev_out                           ! niveau pour les print
    1818
    … …  
    3333          do l=1,nlay+1
    3434             zw2(ig,l)=0.
    35              f_star(ig,l)=0.
    3635             wa_moy(ig,l)=0.
    3736          enddo
    … …  
    4746       enddo
    4847!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
    4977       do l=1,nlay-2
    5078         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!------------------------------------------------------------------------
    5685               zw2(ig,l+1)=2.*RG*(ztv(ig,l)-ztv(ig,l+1))/ztv(ig,l+1)  &
    5786     &                     *(zlev(ig,l+1)-zlev(ig,l))  &
    5887     &                     *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)  &
    63112     &                    *ztv(ig,l))/f_star(ig,l+1)
    64113               zw2(ig,l+1)=zw2(ig,l)*(f_star(ig,l)/f_star(ig,l+1))**2+  &
    … …  
    67116            endif
    68117! 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
    69128            if (zw2(ig,l+1).lt.0.) then
    70129               linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))  &
    71130     &           -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
    72131               zw2(ig,l+1)=0.
    73             else
     132               lmax(ig)=l
     133            endif
     134
    74135               wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
    75             endif
    76136
    77137            if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
    … …  
    84144       if (lev_out.ge.1) print*,'fin calcul zw2'
    85145!
     146!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     147! A eliminer :
     148! Ce calcul de lmax est fait en meme temps que celui de linter, plus haut
    86149! 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
    97162!   
    98163! Determination de zw2 max
    … …  
    119184      do  ig=1,ngrid
    120185! 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)))
    124197           zmax(ig)=max(zmax(ig),zlevinter(ig)-zlev(ig,lmin(ig)))
    125198      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
    127207     
    128208      RETURN

  • LMDZ4/trunk/libf/phytherm/thermcell_flux.F90

    r814 r852  
    2828
    2929      integer ncorecfm1,ncorecfm2,ncorecfm3,ncorecalpha
    30       integer ncorecfm4,ncorecfm5,ncorecfm6
     30      integer ncorecfm4,ncorecfm5,ncorecfm6,ncorecfm7,ncorecfm8
    3131     
    3232
    3333      REAL entr(ngrid,klev),detr(ngrid,klev)
    3434      REAL fm(ngrid,klev+1)
     35      REAL zfm
    3536
    3637      integer igout
    … …  
    4041      REAL f_old,ddd0,eee0,ddd,eee,zzz
    4142
    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
    4648
    4749      ncorecfm1=0
    … …  
    5153      ncorecfm5=0
    5254      ncorecfm6=0
     55      ncorecfm7=0
     56      ncorecfm8=0
    5357      ncorecalpha=0
    5458
    … …  
    5660      fm(:,:)=0.
    5761     
     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
    5873
    5974!-------------------------------------------------------------------------
    … …  
    6681            if (l.le.lmax(ig)) then
    6782               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)
    6992                    print*,'entr_star(ig,l)',entr_star(ig,l)
    7093                    print*,'alim_star(ig,l)',alim_star(ig,l)
    … …  
    7295                    stop
    7396               endif
    74             else
    75                if (abs(entr_star(ig,l))+abs(alim_star(ig,l))+abs(detr_star(ig,l)).gt.0) then
    76                     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                     stop
    81                endif
    8297            endif
    8398         enddo
    … …  
    92107         detr(:,l)=f(:)*detr_star(:,l)
    93108      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
    94120
    95121      fm(:,1)=0.
    … …  
    107133      enddo
    108134
    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  ')
    117153
    118154!-------------------------------------------------------------------------
    … …  
    130166         enddo
    131167      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  ')
    132172
    133173!-------------------------------------------------------------------------
    … …  
    151191      enddo
    152192
     193      if (lev_out.ge.10) &
     194    &    call printflux(ngrid,klev,lunout,igout,f,lmax,lalim, &
     195    &    ptimestep,masse,entr,detr,fm,'4  ')
     196
    153197
    154198!-------------------------------------------------------------------------
    … …  
    167211      enddo
    168212
     213      if (lev_out.ge.10) &
     214    &    call printflux(ngrid,klev,lunout,igout,f,lmax,lalim, &
     215    &    ptimestep,masse,entr,detr,fm,'5  ')
     216
    169217!-------------------------------------------------------------------------
    170218!detr ne peut pas etre superieur a fm
    … …  
    172220
    173221      if(1.eq.1) then
     222
    174223      do l=1,klev
    175224         do ig=1,ngrid
    … …  
    181230               ncorecfm6=ncorecfm6+1
    182231               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
    185251            if (entr(ig,l).lt.0.) then
    186252               print*,'ig,l,lmax(ig)',ig,l,lmax(ig)
    … …  
    194260
    195261
     262      if (lev_out.ge.10) &
     263    &    call printflux(ngrid,klev,lunout,igout,f,lmax,lalim, &
     264    &    ptimestep,masse,entr,detr,fm,'6  ')
     265
    196266!-------------------------------------------------------------------------
    197267!fm ne peut pas etre negatif
    … …  
    207277            endif
    208278            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)
    210280               print*,'detr(ig,l)',detr(ig,l)
    211281               print*,'fm(ig,l)',fm(ig,l)
    … …  
    215285     enddo
    216286
     287      if (lev_out.ge.10) &
     288    &    call printflux(ngrid,klev,lunout,igout,f,lmax,lalim, &
     289    &    ptimestep,masse,entr,detr,fm,'7  ')
     290
    217291!-----------------------------------------------------------------------
    218292!la fraction couverte ne peut pas etre superieure a 1           
    219293!-----------------------------------------------------------------------
     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
    220312     do l=1,klev
    221313        do ig=1,ngrid
    222314           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.)) then
     315           zfm=rhobarz(ig,l+1)*zw2(ig,l+1)*alphamax
     316           if ( fm(ig,l+1) .gt. zfm) then
    225317              f_old=fm(ig,l+1)
    226               fm(ig,l+1)=rhobarz(ig,l+1)*zw2(ig,l+1)
    227               zw2(ig,l+1)=0.
    228               zqla(ig,l+1)=0.
     318              fm(ig,l+1)=zfm
     319!             zw2(ig,l+1)=0.
     320!             zqla(ig,l+1)=0.
    229321              detr(ig,l)=detr(ig,l)+f_old-fm(ig,l+1)
    230               lmax(ig)=l+1
    231               zmax(ig)=zlev(ig,lmax(ig))
     322!             lmax(ig)=l+1
     323!             zmax(ig)=zlev(ig,lmax(ig))
    232324!             print*,'alpha>1',l+1,lmax(ig)
    233325              ncorecalpha=ncorecalpha+1
    … …  
    237329     enddo
    238330!
     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
    239336
    240337!-----------------------------------------------------------------------
    … …  
    248345            eee0=entr(ig,l)
    249346            ddd0=detr(ig,l)
    250             eee=entr(ig,l)-masse(ig,l)*fracmax/ptimestep
     347            eee=entr(ig,l)-masse(ig,l)*fomass_max/ptimestep
    251348            ddd=detr(ig,l)-eee
    252349            if (eee.gt.0.) then
    … …  
    272369                         print*,'detr',detr(ig,l)
    273370                         print*,'masse',masse(ig,l)
    274                          print*,'fracmax',fracmax
    275                          print*,'masse(ig,l)*fracmax/ptimestep',masse(ig,l)*fracmax/ptimestep
     371                         print*,'fomass_max',fomass_max
     372                         print*,'masse(ig,l)*fomass_max/ptimestep',masse(ig,l)*fomass_max/ptimestep
    276373                         print*,'ptimestep',ptimestep
    277374                         print*,'lmax(ig)',lmax(ig)
    … …  
    309406    &     ncorecfm4,'x fm4',ncorecfm5,'x fm5 et', &
    310407    &     ncorecfm6,'x fm6', &
     408    &     ncorecfm7,'x fm7', &
     409    &     ncorecfm8,'x fm8', &
    311410    &     ncorecalpha,'x alpha'
    312411      endif
    313412
     413      if (lev_out.ge.10) &
     414    &    call printflux(ngrid,klev,lunout,igout,f,lmax,lalim, &
     415    &    ptimestep,masse,entr,detr,fm,'fin')
     416
    314417
    315418      return
    316419      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!
    14      SUBROUTINE thermcell_main(ngrid,nlay,ptimestep  &
    25     &                  ,pplay,pplev,pphi,debut  &
    … …  
    5558!   ------
    5659
    57       integer,save :: igout=871
     60!      integer,save :: igout=4259
     61      integer,save :: igout=2928
    5862      integer,save :: lunout=6
    59       integer,save :: lev_out=0
     63      integer,save :: lev_out=10
    6064
    6165      INTEGER ig,k,l,ll
    … …  
    129133      real zlevinter(klon)
    130134      logical debut
     135       real seuil
    131136
    132137!
    … …  
    142147!   ---------------
    143148!
     149
     150      seuil=0.25
     151
    144152      if (lev_out.ge.1) print*,'thermcell_main V4'
    145153
    … …  
    228236
    229237!------------------------------------------------------------------
    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
    256289!------------------------------------------------------------------
    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
    259293!------------------------------------------------------------------
    260294!
    … …  
    262296      CALL thermcell_init(ngrid,nlay,ztv,zlev,  &
    263297     &                    lalim,lmin,alim_star,alim_star_tot,lev_out)
     298
     299call test_ltherm(ngrid,nlay,pplev,pplay,lmin,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_init lmin  ')
     300call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_init lalim ')
     301
    264302
    265303      if (lev_out.ge.1) print*,'thermcell_main apres thermcell_init'
    … …  
    268306         write(lunout,*) 'lmin ',lmin(igout)
    269307         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)
    270311      endif
    271312
    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!------------------------------------------------------------------------------
    275319!
    276320      CALL thermcell_dry(ngrid,nlay,zlev,pphi,ztv,alim_star,  &
    277321     &                      lalim,lmin,zmax_sec,wmax_sec,lev_out)
    278322
     323call test_ltherm(ngrid,nlay,pplev,pplay,lmin,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry  lmin  ')
     324call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry  lalim ')
     325
    279326      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
    280335
    281336
    … …  
    289344     &           lalim,zmax_sec,f0,detr_star,entr_star,f_star,ztva,  &
    290345     &           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  ')
    291348
    292349      if (lev_out.ge.1) print*,'thermcell_main apres thermcell_plume'
    … …  
    297354         write(lunout,*) ' ig l alim_star entr_star detr_star f_star '
    298355         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)
    300357      endif
    301358
    … …  
    307364     &           zlev,lmax,zmax,zmax0,zmix,wmax,lev_out)
    308365
     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
    309372      if (lev_out.ge.1) print*,'thermcell_main apres thermcell_height'
    310373
    … …  
    322385!-------------------------------------------------------------------------------
    323386
    324       CALL thermcell_flux(ngrid,nlay,ptimestep,masse, &
     387      CALL thermcell_flux2(ngrid,nlay,ptimestep,masse, &
    325388     &       lalim,lmax,alim_star,  &
    326389     &       entr_star,detr_star,f,rhobarz,zlev,zw2,fm,entr,  &
    … …  
    328391
    329392      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  ')
    330395
    331396!c------------------------------------------------------------------
    … …  
    391456      enddo
    392457     
     458      print*,'14a OK convect8'
    393459! calcul du niveau de condensation
    394460! initialisation
    … …  
    410476         enddo
    411477      enddo
     478      print*,'14b OK convect8'
    412479      do k=nlay,1,-1
    413480         do ig=1,ngrid
    … …  
    418485         enddo
    419486      enddo
     487      print*,'14c OK convect8'
    420488!calcul des moments
    421489!initialisation
    … …  
    429497         enddo
    430498      enddo     
     499      print*,'14d OK convect8'
    431500      do l=1,nlay
    432501         do ig=1,ngrid
    … …  
    440509            q2(ig,l)=zf2*(zqta(ig,l)*1000.-po(ig,l)*1000.)**2
    441510!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.))
    443512         enddo
    444513      enddo
    445514!calcul du ratqscdiff
     515      print*,'14e OK convect8'
    446516      var=0.
    447517      vardiff=0.
    … …  
    452522         enddo
    453523      enddo
     524      print*,'14f OK convect8'
    454525      do ig=1,ngrid
    455526          do l=1,lalim(ig)
    … …  
    462533          enddo
    463534      enddo
     535      print*,'14g OK convect8'
    464536      do l=1,nlay
    465537         do ig=1,ngrid
    … …  
    496568
    497569!-----------------------------------------------------------------------------
     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  
    11      SUBROUTINE thermcell_plume(ngrid,klev,ztv,zthl,po,zl,rhobarz,  &
    22     &           zlev,pplev,pphi,zpspsk,l_mix,r_aspect,alim_star,  &
    3      &           lentr,zmax_sec,f0,detr_star,entr_star,f_star,ztva,  &
     3     &           lalim,zmax_sec,f0,detr_star,entr_star,f_star,ztva,  &
    44     &           ztla,zqla,zqta,zha,zw2,zqsatth,lmix,linter,lev_out)
    55
    … …  
    2929      REAL l_mix
    3030      REAL r_aspect
    31       INTEGER lentr(ngrid)
     31      INTEGER lalim(ngrid)
    3232      integer lev_out                           ! niveau pour les print
    3333
    … …  
    133133!tests sur la definition du detr
    134134!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
    140148      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
    141156!---------------------------------------------------------------------------
    142157!calcul de l entrainement et du detrainement lateral
    … …  
    217232!calcul de entr_star
    218233!
    219         if (detr_star(ig,l).gt.f_star(ig,l)) then
    220            detr_star(ig,l)=f_star(ig,l)
    221 !a decommenter ou pas?
    222 !           entr_star(ig,l)=0.
    223         endif
    224 
    225234! Deplacement du calcul de entr_star pour eviter d'avoir aussi
    226235! entr_star > fstar.
     236! Redeplacer suite a la transformation du cas detr>f
    227237! 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
    230274      else
    231275         detr_star(ig,l)=0.
    232276         entr_star(ig,l)=0.
    233277      endif
     278
     279!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     280! FH inutile si on conserve comme on l'a fait plus haut entr=detr
     281! dans la couche d'alimentation
    234282!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!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    238287!
    239288!prise en compte du detrainement et de l entrainement dans le calcul du flux
    … …  
    308357!
    309358!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
    310368        if (zw2(ig,l+1).lt.0.) then
    311369           linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))  &
    312370     &               -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
    313371           zw2(ig,l+1)=0.
    314         else
     372        endif
     373
    315374           wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
    316         endif
     375
    317376        if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
    318377!   lmix est le niveau de la couche ou w (wa_moy) est maximum

  • LMDZ4/trunk/libf/phytherm/write_histday.h

    r814 r852  
    491491     $                     zx_tmp_fi2d)
    492492C
     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
    493503      END DO 
    494504c=================================================================

  • LMDZ4/trunk/libf/phytherm/write_histmth.h

    r814 r852  
    681681      IF(lev_histmth.GE.4) THEN
    682682c
     683c  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
    683701cym      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, clwcon0, zx_tmp_3d)
    684702      CALL histwrite_phy(nid_mth,"clwcon",itau_w,clwcon0)
    … …  
    757775c
    758776cIM: 101003 : K/30min ==> K/s
     777      zx_tmp_fi3d(1:klon,1:klev)=d_t_ajsb(1:klon,1:klev)/pdtphys
     778cym      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)
     780c
     781      zx_tmp_fi3d(1:klon,1:klev)=d_q_ajsb(1:klon,1:klev)/pdtphys
     782cym      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)
     784c
     785cIM: 101003 : K/30min ==> K/s
    759786      zx_tmp_fi3d(1:klon,1:klev)=d_t_ajs(1:klon,1:klev)/pdtphys
    760787cym      CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
    761       CALL histwrite_phy(nid_mth,"dtajs",itau_w,zx_tmp_fi3d)
     788      CALL histwrite_phy(nid_mth,"dtthe",itau_w,zx_tmp_fi3d)
    762789c
    763790      zx_tmp_fi3d(1:klon,1:klev)=d_q_ajs(1:klon,1:klev)/pdtphys
    764791cym      CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
    765       CALL histwrite_phy(nid_mth,"dqajs",itau_w,zx_tmp_fi3d)
     792      CALL histwrite_phy(nid_mth,"dqthe",itau_w,zx_tmp_fi3d)
    766793c
    767794cIM: 101003 : K/day ==> K/s
    … …  
    933960      CALL histwrite_phy(nid_mth,"dtcon",itau_w,zx_tmp_fi3d)
    934961c
    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
     965cIM: 101003 : K/30min ==> K/s
     966      zx_tmp_fi3d(1:klon,1:klev)=d_t_ajsb(1:klon,1:klev)/pdtphys
     967cym      CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
    943968      CALL histwrite_phy(nid_mth,"dtajs",itau_w,zx_tmp_fi3d)
     969c
     970      zx_tmp_fi3d(1:klon,1:klev)=d_q_ajsb(1:klon,1:klev)/pdtphys
     971cym      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)
     973c
     974cIM: 101003 : K/30min ==> K/s
     975      zx_tmp_fi3d(1:klon,1:klev)=d_t_ajs(1:klon,1:klev)/pdtphys
     976cym      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)
     978c
     979      zx_tmp_fi3d(1:klon,1:klev)=d_q_ajs(1:klon,1:klev)/pdtphys
     980cym      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)
     982c
     983
    944984c
    945985c  temperature tendency due to large scale precipitation
    … …  
    16721712      IF(lev_histmth.GE.4) THEN
    16731713c
     1714c  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
    16741732cym      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, clwcon0, zx_tmp_3d)
    16751733      CALL histwrite_phy(nid_mth,"clwcon",itau_w,clwcon0)
    … …  
    17391797      CALL histwrite_phy(nid_mth,"ratqs",itau_w,ratqs)
    17401798c
    1741       zx_tmp_fi3d(1:klon,1:klev)=d_q_ajs(1:klon,1:klev)/pdtphys
    1742 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 c
    17451799cIM: 101003 : K/day ==> K/s
    17461800      zx_tmp_fi3d(1:klon,1:klev)=heat(1:klon,1:klev)/RDAY

  • LMDZ4/trunk/libf/phytherm/yamada4.F

    r814 r852  
    458458c   Diagnostique pour stokage
    459459
     460      if(1.eq.0)then
    460461      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)
    465473
    466474c   Estimations de w'2 et T'2 d'apres Abdela et McFarlane
    467475
    468         if(1.eq.0)then
    469       w2yam=q2(:,1:klev)*0.24
    470      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         endif
     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
    476484
    477485c     print*,'OKFIN'
Note: See TracChangeset for help on using the changeset viewer.