Changeset 878


Ignore:
Timestamp:
Jan 14, 2008, 1:03:39 PM (16 years ago)
Author:
Laurent Fairhead
Message:

Bascule de la physique du LMD vers la physique avec thermiques
LF

Location:
LMDZ4/trunk/libf/phylmd
Files:
15 added
17 edited

Legend:

Unmodified
Added
Removed

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

    r766 r878  
    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/phylmd/clouds_gno.F

    r559 r878  
    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/phylmd/coef_diff_turb_mod.F90

    r793 r878  
    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/phylmd/conf_phys.F90

    r840 r878  
    66
    77  subroutine conf_phys(ocean, ok_veget, ok_journe, ok_mensuel, ok_instan, ok_hf, &
     8 &                     seuil_inversion, &
    89 &                     fact_cldcon, facttemps,ok_newmicro,iflag_cldcon, &
    9 !IM&                   ratqsbas,ratqshaut,ip_ebil_phy, &
    10  &                     ratqsbas,ratqshaut, &
     10 &                     iflag_ratqs,ratqsbas,ratqshaut, &
    1111 &                     ok_ade, ok_aie, &
    1212 &                     bl95_b0, bl95_b1,&
    … …  
    5151  real                 :: fact_cldcon, facttemps,ratqsbas,ratqshaut
    5252  integer              :: iflag_cldcon
     53  integer              :: iflag_ratqs
    5354
    5455  character (len = 6),SAVE  :: ocean_omp
    … …  
    6162  real,SAVE           :: ratqshaut_omp
    6263  integer,SAVE        :: iflag_cldcon_omp, ip_ebil_phy_omp
     64  integer,SAVE        :: iflag_ratqs_omp
    6365
    6466! Local
    6567  integer              :: numout = 6
    6668  real                 :: zzz
     69
     70  real :: seuil_inversion
     71  real,save :: seuil_inversion_omp
    6772
    6873  integer :: iflag_thermals,nsplit_thermals
    … …  
    8287  REAL,SAVE :: cdmmax_omp,cdhmax_omp,ksta_omp,ksta_ter_omp
    8388  LOGICAL,SAVE :: ok_kzmin_omp
    84   REAL, SAVE ::  fmagic_omp
     89  REAL, SAVE :: fmagic_omp
    8590  INTEGER,SAVE :: iflag_pbl_omp,lev_histhf_omp,lev_histday_omp,lev_histmth_omp
    8691  CHARACTER*4, SAVE :: type_run_omp
    … …  
    207212!Config Help =
    208213!               
    209 !
    210214  ip_ebil_phy_omp = 0
    211215  call getin('ip_ebil_phy', ip_ebil_phy_omp)
     216!
     217!Config Key  = seuil_inversion
     218!Config Desc = Seuil ur dTh pour le choix entre les schemas de CL
     219!Config Def  = -0.1
     220!Config Help =
     221!               
     222  seuil_inversion_omp = -0.1
     223  call getin('seuil_inversion', seuil_inversion_omp)
     224
    212225!!
    213226!! Constante solaire & Parametres orbitaux & taux gaz effet de serre BEG
    … …  
    438451  reevap_ice_omp = .false.
    439452  call getin('reevap_ice',reevap_ice_omp)
     453
     454!Config Key  = iflag_ratqs
     455!Config Desc =
     456!Config Def  = 1
     457!Config Help =
     458!
     459  iflag_ratqs_omp = 1
     460  call getin('iflag_ratqs',iflag_ratqs_omp)
     461
    440462!
    441463!Config Key  = iflag_cldcon
    … …  
    835857    ratqshaut = ratqshaut_omp
    836858    iflag_cldcon = iflag_cldcon_omp
     859    iflag_ratqs = iflag_ratqs_omp
    837860    ip_ebil_phy = ip_ebil_phy_omp
    838861    iflag_thermals = iflag_thermals_omp
    … …  
    840863    type_run = type_run_omp
    841864    ok_isccp = ok_isccp_omp
     865    seuil_inversion=seuil_inversion_omp
    842866    lonmin_ins = lonmin_ins_omp
    843867    lonmax_ins = lonmax_ins_omp
    … …  
    891915  write(numout,*)' iflag_pdf = ', iflag_pdf
    892916  write(numout,*)' iflag_cldcon = ', iflag_cldcon
     917  write(numout,*)' iflag_ratqs = ', iflag_ratqs
     918  write(numout,*)' seuil_inversion = ', seuil_inversion
    893919  write(numout,*)' fact_cldcon = ', fact_cldcon
    894920  write(numout,*)' facttemps = ', facttemps

  • LMDZ4/trunk/libf/phylmd/dimphy.h

    r524 r878  
    22! $Header$
    33!
    4 c-----------------------------------------------------------------------
     4!-----------------------------------------------------------------------
    55      INTEGER KIDIA, KFDIA, KLON, KLEV
    6       PARAMETER (KIDIA=1,KFDIA=iim*(jjm-1)+2-1/jjm,
    7      .           KLON=KFDIA-KIDIA+1,KLEV=llm)
    8 c-----------------------------------------------------------------------
     6      PARAMETER (KIDIA=1,KFDIA=iim*(jjm-1)+2-1/jjm)
     7      PARAMETER (KLON=KFDIA-KIDIA+1,KLEV=llm)
     8!-----------------------------------------------------------------------
    99      INTEGER nbtr ! nombre de vrais traceurs
    1010      PARAMETER (nbtr=nqmx-2+1/(nqmx-1))
    11 c-----------------------------------------------------------------------
     11!-----------------------------------------------------------------------
    1212      REAL zmasq(KLON)
    1313      COMMON/terreoce/zmasq

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

    r860 r878  
    5353c Champs 2D:
    5454c
     55         CALL histdef(nid_day, "weakinv", "Weak inversion", "",
     56     s           iim,jjmp1,nhori, 1,1,1, -99, 32,
     57     s           "ave(X)", zstophy,zout)
     58
     59         CALL histdef(nid_day, "dthmin", "dTheta mini", "K/m",
     60     s           iim,jjmp1,nhori, 1,1,1, -99, 32,
     61     s           "ave(X)", zstophy,zout)
     62
     63
    5564         CALL histdef(nid_day, "tsol", "Surface Temperature", "K",
    5665     .                iim,jj_nb,nhori, 1,1,1, -99, 32,
    … …  
    538547     $         "ave(X)", zstophy,zout)
    539548C
     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
    540565         END DO
    541566C           

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

    r864 r878  
    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",
    … …  
    879910         if (nqmax.GE.3) THEN
    880911           DO iq=3,nqmax
    881            iiq=niadv(iq)
     912             iiq=niadv(iq)
    882913             CALL histdef(nid_mth, tnom(iq), ttext(iiq), "-",
    883914     .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
    884915     .                "ave(X)", zstophy,zout)
    885            ENDDO
    886          ENDIF
     916             ENDDO
     917           ENDIF
    887918#endif
    888919c
    … …  
    960991c
    961992          CALL histdef(nid_mth, "dtcon","dtcon","K/s",
     993     .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
     994     .                "ave(X)", zstophy,zout)
     995c
     996         CALL histdef(nid_mth, "dtthe", "Dry adjust. dT", "K/s",
     997     .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
     998     .                "ave(X)", zstophy,zout)
     999
     1000         CALL histdef(nid_mth,"dqthe","Dry adjust. dQ","(kg/kg)/s",
    9621001     .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
    9631002     .                "ave(X)", zstophy,zout)
    … …  
    16681707      IF(lev_histmth.GE.4) THEN
    16691708c
     1709!FH Sorties pour la couche limite
     1710      if (iflag_pbl>1) then
     1711         CALL histdef(nid_mth, "tke","TKE","m2/s2",
     1712     .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
     1713     .                "ave(X)", zstophy,zout)
     1714c
     1715         CALL histdef(nid_mth, "tke_max","TKE max","m2/s2",
     1716     .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
     1717     .                "t_max(X)", zstophy,zout)
     1718      endif
     1719c
     1720         CALL histdef(nid_mth, "kz","Kz melange","m2/s",
     1721     .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
     1722     .                "ave(X)", zstophy,zout)
     1723c
     1724         CALL histdef(nid_mth, "kz_max","Kz melange max","m2/s",
     1725     .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
     1726     .                "t_max(X)", zstophy,zout)
     1727
     1728
     1729
     1730
    16701731         CALL histdef(nid_mth, "clwcon",
    16711732     .                "Convective Cloud Liquid water content"
    … …  
    18311892#ifndef INCA
    18321893         if (nqmax.GE.3) THEN
    1833            DO iq=3,nqmax
    1834              iiq=niadv(iq)
    1835              CALL histdef(nid_mth, tnom(iq), ttext(iiq), "-",
    1836      .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
    1837      .                "ave(X)", zstophy,zout)
    1838              ENDDO
    1839            ENDIF
     1894         DO iq=3,nqmax
     1895         iiq=niadv(iq)
     1896         CALL histdef(nid_mth, tnom(iq), ttext(iiq), "-",
     1897     .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
     1898     .                "ave(X)", zstophy,zout)
     1899         ENDDO
     1900         ENDIF
    18401901#endif
    18411902

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

    r803 r878  
    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/phylmd/phyetat0.F

    r782 r878  
    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/phylmd/phyredem.F

    r782 r878  
    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/phylmd/physiq.F

    r864 r878  
    177177      REAL fm_therm(klon,klev+1)
    178178      REAL entr_therm(klon,klev)
     179      real,allocatable,save :: clwcon0th(:,:),rnebcon0th(:,:)
     180c$OMP THREADPRIVATE(clwcon0th,rnebcon0th)
     181
     182      real weak_inversion(klon),dthmin(klon)
     183      real seuil_inversion
     184      save seuil_inversion
     185c$OMP THREADPRIVATE(seuil_inversion)
     186      integer iflag_ratqs
     187      save iflag_ratqs
     188c$OMP THREADPRIVATE(iflag_ratqs)
     189
     190      integer lmax_th(klon)
     191      integer limbas(klon)
     192      real ratqscth(klon,klev)
     193      real ratqsdiff(klon,klev)
     194      real zqsatth(klon,klev)
     195
    179196c======================================================================
    180197c
    … …  
    929946      REAL fluxu(klon,klev, nbsrf)   ! flux turbulent de vitesse u
    930947      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
    931952c
    932953      REAL zxfluxt(klon, klev)
    … …  
    10501071c$OMP THREADPRIVATE(d_u_con,d_v_con)
    10511072      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)
    10521074      REAL d_t_ajs(klon,klev), d_q_ajs(klon,klev)
    10531075      REAL d_u_ajs(klon,klev), d_v_ajs(klon,klev)
    … …  
    14711493#endif
    14721494      allocate( clwcon0(klon,klev),rnebcon0(klon,klev))
     1495      allocate( clwcon0th(klon,klev),rnebcon0th(klon,klev))
    14731496      allocate( tau_ae(klon,klev,2), piz_ae(klon,klev,2))
    14741497      allocate( cg_ae(klon,klev,2))
    … …  
    15071530         CALL suphec ! initialiser constantes et parametres phys.
    15081531      ENDIF
     1532
     1533      print*,'CONVERGENCE PHYSIQUE THERM 1 '
    15091534
    15101535
    … …  
    15671592c
    15681593         call conf_phys(ocean, ok_veget, ok_journe, ok_mensuel,
    1569      .                  ok_instan, ok_hf,
     1594     .                  ok_instan, ok_hf, seuil_inversion,
    15701595     .                  fact_cldcon, facttemps,ok_newmicro,
    1571 cIM  .                  iflag_cldcon,ratqsbas,ratqshaut, if_ebil,
    1572      .                  iflag_cldcon,ratqsbas,ratqshaut,
     1596     .                  iflag_cldcon,iflag_ratqs,ratqsbas,ratqshaut,
    15731597     .                  ok_ade, ok_aie,
    15741598     .                  bl95_b0, bl95_b1,
    … …  
    15811605         itap    = 0
    15821606         itaprad = 0
     1607
     1608!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     1609!! Un petit travail à faire ici.
     1610!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     1611
     1612         if (iflag_pbl>1) then
     1613             PRINT*, "Using method MELLOR&YAMADA"
     1614         endif
     1615          ! NB! pbl_tke could/should be read and written from (re)startphy.nc
     1616          ALLOCATE(pbl_tke(klon,klev+1,nbsrf)) 
     1617          pbl_tke(:,:,:) = 1.e-8
     1618
    15831619
    15841620         CALL phyetat0 ("startphy.nc",dtime,co2_ppm_etat0,solaire_etat0,
    … …  
    15891625     .       radsol,clesphy0,
    15901626     .       zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,tabcntr0,
    1591      .       t_ancien, q_ancien, ancien_ok, rnebcon, ratqs,clwcon)
     1627     .       t_ancien, q_ancien, ancien_ok, rnebcon, ratqs,clwcon,
     1628     .       pbl_tke)
     1629
     1630
     1631
     1632
     1633!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     1634
    15921635
    15931636       DO i=1,klon
    … …  
    16201663     .                        pdtphys
    16211664            abort_message='Pas physique n est pas correct '
    1622             call abort_gcm(modname,abort_message,1)
     1665!           call abort_gcm(modname,abort_message,1)
     1666            dtime=pdtphys
    16231667         ENDIF
    16241668         IF (nlon .NE. klon) THEN
    … …  
    16701714         ENDIF
    16711715
     1716           rugoro=0.
    16721717c34EK
    16731718         IF (ok_orodr) THEN
    1674            DO i=1,klon
    1675              rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
    1676            ENDDO
     1719
     1720           rugoro=0.
     1721
     1722!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     1723! FH sans doute a enlever de finitivement ou, si on le garde, l'activer
     1724! justement quand ok_orodr = false.
     1725! ce rugoro est utilise par la couche limite et fait double emploi
     1726! avec les paramétrisations spécifiques de Francois Lott.
     1727!           DO i=1,klon
     1728!             rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
     1729!           ENDDO
     1730!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     1731
    16771732           CALL SUGWD(klon,klev,paprs,pplay)
    16781733           DO i=1,klon
    … …  
    20882143     d     wfbils,    wfbilo,    fluxt,   fluxu,  fluxv,
    20892144     -     dsens,     devap,     zxsnow,
    2090      -     zxfluxt,   zxfluxq,   q2m,     fluxq )
     2145     -     zxfluxt,   zxfluxq,   q2m,     fluxq, pbl_tke )
    20912146c
    20922147c
    … …  
    22602315          DO i = 1, klon
    22612316            ema_pct(i)  = paprs(i,itop_con(i))
     2317            if (itop_con(i).gt.klev-3) then
     2318               print*,'La convection monte trop haut '
     2319               print*,'itop_con(,',i,',)=',itop_con(i)
     2320            endif
    22622321          ENDDO
    22632322          DO i = 1, klon
    … …  
    23432402c===================================================================
    23442403c
     2404       call stratocu_if(klon,klev,pctsrf,paprs, pplay,t_seri
     2405     s ,seuil_inversion,weak_inversion,dthmin)
     2406
     2407
     2408
     2409      d_t_ajsb(:,:)=0.
     2410      d_q_ajsb(:,:)=0.
    23452411      d_t_ajs(:,:)=0.
    23462412      d_u_ajs(:,:)=0.
    23472413      d_v_ajs(:,:)=0.
    23482414      d_q_ajs(:,:)=0.
     2415      clwcon0th(:,:)=0.
    23492416      fm_therm(:,:)=0.
    23502417      entr_therm(:,:)=0.
    … …  
    23572424c  ====
    23582425         IF(prt_level>9)WRITE(lunout,*)'pas de convection'
    2359       else if(iflag_thermals.eq.0) then
    2360 
    2361 c  Ajustement sec
    2362 c  ==============
    2363          IF(prt_level>9)WRITE(lunout,*)'ajsec'
    2364          CALL ajsec(paprs, pplay, t_seri,q_seri, d_t_ajs, d_q_ajs)
    2365          t_seri(:,:) = t_seri(:,:) + d_t_ajs(:,:)
    2366          q_seri(:,:) = q_seri(:,:) + d_q_ajs(:,:)
     2426
     2427
    23672428      else
     2429
    23682430c  Thermiques
    23692431c  ==========
    23702432         IF(prt_level>9)WRITE(lunout,*)'JUSTE AVANT , iflag_thermals='
    23712433     s   ,iflag_thermals,'   nsplit_thermals=',nsplit_thermals
     2434         print*,'JUSTE AVANT , iflag_thermals='
     2435     s   ,iflag_thermals,'   nsplit_thermals=',nsplit_thermals
     2436
     2437
     2438         if (iflag_thermals.gt.1) then
    23722439         call calltherm(pdtphys
    2373      s      ,pplay,paprs,pphi
    2374      s      ,u_seri,v_seri,t_seri,q_seri
     2440     s      ,pplay,paprs,pphi,weak_inversion
     2441     s      ,u_seri,v_seri,t_seri,q_seri,zqsat,debut
    23752442     s      ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs
    2376      s      ,fm_therm,entr_therm)
     2443     s      ,fm_therm,entr_therm,zqasc,clwcon0th,lmax_th,ratqscth,
     2444     s       ratqsdiff,zqsatth)
     2445         endif
     2446
     2447!        call calltherm(pdtphys
     2448!    s      ,pplay,paprs,pphi
     2449!    s      ,u_seri,v_seri,t_seri,q_seri
     2450!    s      ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs
     2451!    s      ,fm_therm,entr_therm)
     2452
     2453c  Ajustement sec
     2454c  ==============
     2455
     2456! Dans le cas où on active les thermiques, on fait partir l'ajustement
     2457! a partir du sommet des thermiques.
     2458! Dans le cas contraire, on demarre au niveau 1.
     2459
     2460         if (iflag_thermals.ge.13.or.iflag_thermals.eq.0) then
     2461
     2462         if(iflag_thermals.eq.0) then
     2463            IF(prt_level>9)WRITE(lunout,*)'ajsec'
     2464            limbas(:)=1
     2465         else
     2466            limbas(:)=lmax_th(:)
     2467         endif
     2468 
     2469! Attention : le call ajsec_convV2 n'est maintenu que momentanneement
     2470! pour des test de convergence numerique.
     2471! Le nouveau ajsec est a priori mieux, meme pour le cas
     2472! iflag_thermals = 0 (l'ancienne version peut faire des tendances
     2473! non nulles numeriquement pour des mailles non concernees.
     2474
     2475         if (iflag_thermals.eq.0) then
     2476            CALL ajsec_convV2(paprs, pplay, t_seri,q_seri
     2477     s      , d_t_ajsb, d_q_ajsb)
     2478         else
     2479            CALL ajsec(paprs, pplay, t_seri,q_seri,limbas
     2480     s      , d_t_ajsb, d_q_ajsb)
     2481         endif
     2482
     2483         t_seri(:,:) = t_seri(:,:) + d_t_ajsb(:,:)
     2484         q_seri(:,:) = q_seri(:,:) + d_q_ajsb(:,:)
     2485         d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
     2486         d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
     2487
     2488         endif
     2489
    23772490      endif
    23782491c
    … …  
    24062519         enddo
    24072520         enddo
     2521      else if (iflag_cldcon.eq.4) then
     2522         ptconvth(:,:)=.false.
     2523         ratqsc(:,:)=0.
     2524         print*,'avant clouds_gno thermique'
     2525         call clouds_gno
     2526     s   (klon,klev,q_seri,zqsat,clwcon0th,ptconvth,ratqsc,rnebcon0th)
     2527         print*,' CLOUDS_GNO OK'
    24082528      endif
    24092529
    24102530c   ratqs stables
    24112531c   -------------
    2412       do k=1,klev
    2413 cIM RAJOUT boucle do=i
    2414        do i=1, klon
    2415 cIM         ratqss(:,k)=ratqsbas+(ratqshaut-ratqsbas)*
    2416 cIM     s   min((paprs(:,1)-pplay(:,k))/(paprs(:,1)-30000.),1.)
    2417          ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
    2418      s   min((paprs(i,1)-pplay(i,k))/(paprs(i,1)-30000.),1.)
    2419 cIM      print*,' IMratqs STABLE i, k',i,k,ratqss(i,k)
    2420        enddo
    2421       enddo
     2532
     2533      if (iflag_ratqs.eq.0) then
     2534
     2535! Le cas iflag_ratqs=0 correspond a la version IPCC 2005 du modele.
     2536         do k=1,klev
     2537            do i=1, klon
     2538               ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
     2539     s         min((paprs(i,1)-pplay(i,k))/(paprs(i,1)-30000.),1.)
     2540            enddo
     2541         enddo
     2542
     2543! Pour iflag_ratqs=1 ou 2, le ratqs est constant au dessus de
     2544! 300 hPa (ratqshaut), varie lineariement en fonction de la pression
     2545! entre 600 et 300 hPa et est soit constant (ratqsbas) pour iflag_ratqs=1
     2546! soit lineaire (entre 0 a la surface et ratqsbas) pour iflag_ratqs=2
     2547! Il s'agit de differents tests dans la phase de reglage du modele
     2548! avec thermiques.
     2549
     2550      else if (iflag_ratqs.eq.1) then
     2551
     2552         do k=1,klev
     2553            do i=1, klon
     2554               if (pplay(i,k).ge.60000.) then
     2555                  ratqss(i,k)=ratqsbas
     2556               else if ((pplay(i,k).ge.30000.).and.
     2557     s            (pplay(i,k).lt.60000.)) then
     2558                  ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
     2559     s            (60000.-pplay(i,k))/(60000.-30000.)
     2560               else
     2561                  ratqss(i,k)=ratqshaut
     2562               endif
     2563            enddo
     2564         enddo
     2565
     2566      else
     2567
     2568         do k=1,klev
     2569            do i=1, klon
     2570               if (pplay(i,k).ge.60000.) then
     2571                  ratqss(i,k)=ratqsbas
     2572     s            *(paprs(i,1)-pplay(i,k))/(paprs(i,1)-60000.)
     2573               else if ((pplay(i,k).ge.30000.).and.
     2574     s             (pplay(i,k).lt.60000.)) then
     2575                    ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
     2576     s              (60000.-pplay(i,k))/(60000.-30000.)
     2577               else
     2578                    ratqss(i,k)=ratqshaut
     2579               endif
     2580            enddo
     2581         enddo
     2582      endif
     2583
     2584
    24222585
    24232586
    24242587c  ratqs final
    24252588c  -----------
    2426       if (iflag_cldcon.eq.1 .or.iflag_cldcon.eq.2) then
    2427 c   les ratqs sont une conbinaison de ratqss et ratqsc
    2428 c   ratqs final
    2429 c   1e4 (en gros 3 heures), en dur pour le moment, est le temps de
    2430 c   relaxation des ratqs
    2431 c         facttemps=exp(-pdtphys/1.e4)
    2432          facteur=exp(-pdtphys*facttemps)
    2433          ratqs(:,:)=max(ratqs(:,:)*facteur,ratqss(:,:))
    2434          ratqs(:,:)=max(ratqs(:,:),ratqsc(:,:))
    2435 c         print*,'calcul des ratqs fini'
     2589
     2590      if (iflag_cldcon.eq.1 .or.iflag_cldcon.eq.2
     2591     s    .or.iflag_cldcon.eq.4) then
     2592
     2593! On ajoute une constante au ratqsc*2 pour tenir compte de
     2594! fluctuations turbulentes de petite echelle
     2595
     2596         do k=1,klev
     2597            do i=1,klon
     2598               if ((fm_therm(i,k).gt.1.e-10)) then
     2599                  ratqsc(i,k)=sqrt(ratqsc(i,k)**2+0.05**2)
     2600               endif
     2601            enddo
     2602         enddo
     2603
     2604!   les ratqs sont une conbinaison de ratqss et ratqsc
     2605!   1800s, en dur pour le moment, est le temps de
     2606!   relaxation des ratqs
     2607
     2608         facteur=exp(-pdtphys/1800.)
     2609
     2610         print*,'WARNING ratqs a revoir '
     2611         ratqs(:,:)=ratqsc(:,:)*(1.-facteur)+ratqs(:,:)*facteur
     2612         ratqs(:,:)=sqrt(ratqs(:,:)**2+ratqss(:,:)**2)
     2613
    24362614      else
    2437 c   on ne prend que le ratqs stable pour fisrtilp
     2615!   on ne prend que le ratqs stable pour fisrtilp
    24382616         ratqs(:,:)=ratqss(:,:)
    24392617      endif
    … …  
    25052683cIM cf FH
    25062684c     IF (iflag_cldcon.eq.-1) THEN ! seulement pour Tiedtke
    2507        IF (iflag_cldcon.le.-1) THEN ! seulement pour Tiedtke
     2685      IF (iflag_cldcon.le.-1) THEN ! seulement pour Tiedtke
    25082686       snow_tiedtke=0.
    25092687c     print*,'avant calcul de la pseudo precip '
    … …  
    25422720      ENDDO
    25432721
    2544       ELSE IF (iflag_cldcon.eq.3) THEN
     2722      ELSE IF (iflag_cldcon.ge.3) THEN
    25452723c  On prend pour les nuages convectifs le max du calcul de la
    25462724c  convection et du calcul du pas de temps precedent diminue d'un facteur
    25472725c  facttemps
    2548 c      facttemps=pdtphys/1.e4
    25492726      facteur = pdtphys *facttemps
    25502727      do k=1,klev
    … …  
    33143491     .      radsol,
    33153492     .      zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,
    3316      .      t_ancien, q_ancien, rnebcon, ratqs, clwcon)
     3493     .      t_ancien, q_ancien, rnebcon, ratqs, clwcon,
     3494     .      pbl_tke)
    33173495      ENDIF
    33183496     

  • LMDZ4/trunk/libf/phylmd/thermcell.F

    r776 r878  
    1 !
    2 ! $Header$
    3 !
    4 
    5       SUBROUTINE thermcell(ngrid,nlay,ptimestep
    6      s                  ,pplay,pplev,pphi
     1      SUBROUTINE calcul_sec(ngrid,nlay,ptimestep
     2     s                  ,pplay,pplev,pphi,zlev
    73     s                  ,pu,pv,pt,po
    8      s                  ,pduadj,pdvadj,pdtadj,pdoadj
    9      s                  ,fm0,entr0
     4     s                  ,zmax,wmax,zw2,lmix
    105c    s                  ,pu_therm,pv_therm
    116     s                  ,r_aspect,l_mix,w2di,tho)
    12       USE dimphy
     7
    138      IMPLICIT NONE
    149
    … …  
    1813c   de "thermiques" explicitement representes
    1914c
    20 c   Reecriture a partir d'un listing papier à Habas, le 14/02/00
    21 c
    22 c   le thermique est suppose homogene et dissipe par melange avec
    23 c   son environnement. la longueur l_mix controle l'efficacite du
    24 c   melange
    25 c
    26 c   Le calcul du transport des differentes especes se fait en prenant
     15c   Réécriture à partir d'un listing papier à Habas, le 14/02/00
     16c
     17c   le thermique est supposé homogène et dissipé par mélange avec
     18c   son environnement. la longueur l_mix contrôle l'efficacité du
     19c   mélange
     20c
     21c   Le calcul du transport des différentes espèces se fait en prenant
    2722c   en compte:
    2823c     1. un flux de masse montant
    … …  
    3732c   -------------
    3833
    39 cym#include "dimensions.h"
    40 cym#include "dimphy.h"
     34#include "dimensions.h"
     35#include "dimphy.h"
    4136#include "YOMCST.h"
    4237
    … …  
    5651      save idetr
    5752      data idetr/3/
    58 c$OMP THREADPRIVATE(idetr)
     53
    5954c   local:
    6055c   ------
    … …  
    8984      real fracc(klon,klev+1)
    9085      real zf,zf2
    91       real,allocatable,save :: thetath2(:,:),wth2(:,:)
    92 c$OMP THREADPRIVATE(thetath2,wth2)
    93 cym      common/comtherm/thetath2,wth2
     86      real thetath2(klon,klev),wth2(klon,klev)
     87      common/comtherm/thetath2,wth2
    9488
    9589      real count_time
    … …  
    9892      data isplit/0/
    9993      save isplit
    100 c$OMP THREADPRIVATE(isplit)
     94
    10195      logical sorties
    10296      real rho(klon,klev),rhobarz(klon,klev+1),masse(klon,klev)
    … …  
    116110      real fm(klon,klev+1),entr(klon,klev)
    117111      real fmc(klon,klev+1)
    118 
     112 
    119113cCR:nouvelles variables
    120114      real f_star(klon,klev+1),entr_star(klon,klev)
    121115      real entr_star_tot(klon),entr_star2(klon)
     116      real zalim(klon)
     117      integer lalim(klon)
     118      real norme(klon)
    122119      real f(klon), f0(klon)
    123120      real zlevinter(klon)
     121      logical therm
    124122      logical first
    125123      data first /.false./
    126124      save first
    127 c$OMP THREADPRIVATE(first)
    128125cRC
    129126
    … …  
    138135      save ncorrec
    139136      data ncorrec/0/
    140 c$OMP THREADPRIVATE(ncorrec)
    141       logical,save :: firstCall=.true.
    142 c$OMP THREADPRIVATE(firstCall)
     137     
    143138c
    144139c-----------------------------------------------------------------------
    … …  
    146141c   ---------------
    147142c
    148       if (firstcall) then
    149         allocate(thetath2(klon,klev),wth2(klon,klev))
    150         thetath2(:,:)=0.
    151         wth2(:,:)=0.
    152         firstcall=.false.
    153       endif
    154      
    155143       sorties=.true.
    156144      IF(ngrid.NE.klon) THEN
    … …  
    165153c   ---------------------------------------------------
    166154
    167        print*,'0 OK convect8'
     155c       print*,'0 OK convect8'
    168156
    169157      DO 1010 l=1,nlay
    … …  
    1781661010  CONTINUE
    179167
    180        print*,'1 OK convect8'
     168c       print*,'1 OK convect8'
    181169c                       --------------------
    182170c
    … …  
    291279
    292280con ne considere que les premieres couches instables
     281      therm=.false.
    293282      do k=nlay-2,1,-1
    294283         do ig=1,ngrid
    295284            if (ztv(ig,k).gt.ztv(ig,k+1).and.
    296285     s          ztv(ig,k+1).le.ztv(ig,k+2)) then
    297                lentr(ig)=k
     286               lentr(ig)=k+1
     287               therm=.true.
    298288            endif
    299289          enddo
    300290      enddo
    301    
     291climitation de la valeur du lentr
     292c      do ig=1,ngrid
     293c         lentr(ig)=min(5,lentr(ig))
     294c      enddo
    302295c determination du lmin: couche d ou provient le thermique
    303296      do ig=1,ngrid
    … …  
    309302               lmin(ig)=l-1
    310303            endif
     304         enddo
     305      enddo
     306cinitialisations
     307      do ig=1,ngrid
     308         zalim(ig)=0.
     309         norme(ig)=0.
     310         lalim(ig)=1
     311      enddo
     312      do k=1,klev-1
     313         do ig=1,ngrid
     314       zalim(ig)=zalim(ig)+zlev(ig,k)*MAX(0.,(ztv(ig,k)-ztv(ig,k+1))
     315     s          /(zlev(ig,k+1)-zlev(ig,k)))
     316c     s         *(zlev(ig,k+1)-zlev(ig,k))
     317       norme(ig)=norme(ig)+MAX(0.,(ztv(ig,k)-ztv(ig,k+1))
     318     s          /(zlev(ig,k+1)-zlev(ig,k)))
     319c    s          *(zlev(ig,k+1)-zlev(ig,k))
     320         enddo
     321       enddo
     322       do ig=1,ngrid
     323          if (norme(ig).gt.1.e-10) then
     324             zalim(ig)=max(10.*zalim(ig)/norme(ig),zlev(ig,2))
     325c             zalim(ig)=min(zalim(ig),zlev(ig,lentr(ig)))
     326          endif
     327       enddo
     328cdétermination du lalim correspondant
     329      do k=1,klev-1
     330         do ig=1,ngrid
     331      if ((zalim(ig).gt.zlev(ig,k)).and.(zalim(ig).le.zlev(ig,k+1)))
     332     s   then
     333         lalim(ig)=k
     334      endif
    311335         enddo
    312336      enddo
    … …  
    316340         do ig=1,ngrid
    317341            if (ztv(ig,l).gt.ztv(ig,l+1).and.
    318      s          l.ge.lmin(ig).and.l.le.lentr(ig)) then
    319                  entr_star(ig,l)=(ztv(ig,l)-ztv(ig,l+1))*
    320      s                           (zlev(ig,l+1)-zlev(ig,l))
    321             endif
    322          enddo
    323       enddo
    324 c pas de thermique si couches 1->5 stables
     342     s          l.ge.lmin(ig).and.l.lt.lentr(ig)) then
     343                 entr_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.)
     344c     s                           *(zlev(ig,l+1)-zlev(ig,l))
     345     s                           *sqrt(zlev(ig,l+1))
     346cautre def
     347c                entr_star(ig,l)=zlev(ig,l+1)*(1.-(zlev(ig,l+1)
     348c     s                         /zlev(ig,lentr(ig)+2)))**(3./2.)
     349            endif
     350         enddo
     351      enddo
     352cnouveau test
     353c      if (therm) then
     354      do l=1,klev-1
     355         do ig=1,ngrid
     356            if (ztv(ig,l).gt.ztv(ig,l+1).and.
     357     s          l.ge.lmin(ig).and.l.le.lalim(ig)
     358     s          .and.zalim(ig).gt.1.e-10) then
     359c            if (l.le.lentr(ig)) then
     360c               entr_star(ig,l)=zlev(ig,l+1)*(1.-(zlev(ig,l+1)
     361c     s                         /zalim(ig)))**(3./2.)
     362c               write(10,*)zlev(ig,l),entr_star(ig,l)
     363            endif
     364         enddo
     365      enddo
     366c      endif
     367c pas de thermique si couche 1 stable
    325368      do ig=1,ngrid
    326369         if (lmin(ig).gt.5) then
    … …  
    339382         enddo
    340383      enddo
    341 c
    342       print*,'fin calcul entr_star'
     384c Calcul entrainement normalise
     385      do ig=1,ngrid
     386         if (entr_star_tot(ig).gt.1.e-10) then
     387c         do l=1,lentr(ig)
     388          do l=1,klev
     389cdef possibles pour entr_star: zdthetadz, dthetadz, zdtheta
     390            entr_star(ig,l)=entr_star(ig,l)/entr_star_tot(ig)
     391         enddo
     392         endif
     393      enddo
     394c
     395c      print*,'fin calcul entr_star'
    343396      do k=1,klev
    344397         do ig=1,ngrid
    … …  
    394447ctest
    395448               if (abs(zw2(ig,l+1)-zw2(ig,l)).lt.1e-10) then
    396                   print*,'pb linter'
     449c                  print*,'pb linter'
    397450               endif
    398451               linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))
    … …  
    402455            else
    403456               if (zw2(ig,l+1).lt.0.) then
    404                   print*,'pb1 zw2<0'
     457c                  print*,'pb1 zw2<0'
    405458               endif
    406459               wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
    … …  
    413466         enddo
    414467      enddo
    415       print*,'fin calcul zw2'
     468c      print*,'fin calcul zw2'
    416469c
    417470c Calcul de la couche correspondant a la hauteur du thermique
    418471      do ig=1,ngrid
    419472         lmax(ig)=lentr(ig)
     473c         lmax(ig)=lalim(ig)
    420474      enddo
    421475      do ig=1,ngrid
    422476         do l=nlay,lentr(ig)+1,-1
     477c         do l=nlay,lalim(ig)+1,-1
    423478            if (zw2(ig,l).le.1.e-10) then
    424479               lmax(ig)=l-1
    … …  
    426481         enddo
    427482      enddo
    428 c pas de thermique si couches 1->5 stables
     483c pas de thermique si couche 1 stable
    429484      do ig=1,ngrid
    430485         if (lmin(ig).gt.5) then
    431486            lmax(ig)=1
    432487            lmin(ig)=1
     488            lentr(ig)=1
     489            lalim(ig)=1
    433490         endif
    434491      enddo
    … …  
    443500            if (l.le.lmax(ig)) then
    444501                if (zw2(ig,l).lt.0.)then
    445                   print*,'pb2 zw2<0'
     502c                  print*,'pb2 zw2<0'
    446503                endif
    447504                zw2(ig,l)=sqrt(zw2(ig,l))
    … …  
    465522       zmax(ig)=max(zmax(ig),zlevinter(ig)-zlev(ig,lmin(ig)))
    466523      enddo
    467 
    468       print*,'avant fermeture'
     524      do ig=1,ngrid
     525c      write(8,*)zmax(ig),lmax(ig),lentr(ig),lmin(ig)
     526      enddo
     527con stoppe après les calculs de zmax et wmax
     528      RETURN
     529
     530c      print*,'avant fermeture'
    469531c Fermeture,determination de f
     532cAttention! entrainement normalisé ou pas?
    470533      do ig=1,ngrid
    471534         entr_star2(ig)=0.
    … …  
    476539         else
    477540             do k=lmin(ig),lentr(ig)
     541c             do k=lmin(ig),lalim(ig)
    478542                entr_star2(ig)=entr_star2(ig)+entr_star(ig,k)**2
    479543     s                    /(rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k)))
    … …  
    481545c Nouvelle fermeture
    482546             f(ig)=wmax(ig)/(max(500.,zmax(ig))*r_aspect
    483      s             *entr_star2(ig))*entr_star_tot(ig)
     547     s             *entr_star2(ig))
     548c     s            *entr_star_tot(ig)
    484549ctest
    485550c             if (first) then
    486 c             f(ig)=f(ig)+(f0(ig)-f(ig))*exp(-ptimestep/zmax(ig)
    487 c     s             *wmax(ig))
     551             f(ig)=f(ig)+(f0(ig)-f(ig))*exp(-ptimestep/zmax(ig)
     552     s             *wmax(ig))
    488553c             endif
    489554         endif
    490 c         f0(ig)=f(ig)
     555         f0(ig)=f(ig)
    491556c         first=.true.
    492557      enddo
    493       print*,'apres fermeture'
    494 
     558c      print*,'apres fermeture'
     559con stoppe après la fermeture
     560      RETURN
    495561c Calcul de l'entrainement
    496562       do k=1,klev
    … …  
    499565         enddo
    500566      enddo
     567con stoppe après le calcul de entr
     568c      RETURN
     569cCR:test pour entrainer moins que la masse
     570c       do ig=1,ngrid
     571c          do l=1,lentr(ig)
     572c             if ((entr(ig,l)*ptimestep).gt.(0.9*masse(ig,l))) then
     573c                entr(ig,l+1)=entr(ig,l+1)+entr(ig,l)
     574c     s                       -0.9*masse(ig,l)/ptimestep
     575c                entr(ig,l)=0.9*masse(ig,l)/ptimestep
     576c             endif
     577c          enddo
     578c       enddo
     579cCR: fin test
    501580c Calcul des flux
    502581      do ig=1,ngrid
    … …  
    516595c   calcul de la largeur de chaque ascendance dans le cas conservatif.
    517596c   dans ce cas simple, on suppose que la largeur de l'ascendance provenant
    518 c   d'une couche est egale a la hauteur de la couche alimentante.
     597c   d'une couche est égale à la hauteur de la couche alimentante.
    519598c   La vitesse maximale dans l'ascendance est aussi prise comme estimation
    520599c   de la vitesse d'entrainement horizontal dans la couche alimentante.
    … …  
    536615c  cette option est finalement en dur.
    537616                  if ((l_mix*zlev(ig,l)).lt.0.)then
    538                    print*,'pb l_mix*zlev<0'
     617c                   print*,'pb l_mix*zlev<0'
    539618                  endif
     619cCR: test: nouvelle def de lambda
     620c                  larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
     621                  if (zw2(ig,l).gt.1.e-10) then
     622                  larg_detr(ig,l)=sqrt((l_mix/zw2(ig,l))*zlev(ig,l))
     623                  else
    540624                  larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
     625                  endif
     626cRC
    541627c              else if (idetr.eq.1) then
    542628c                 larg_detr(ig,l)=larg_cons(ig,l)
    … …  
    555641c      print*,'10 OK convect8'
    556642c     print*,'WA2 ',wa_moy
    557 c   calcul de la fraction de la maille concerne par l'ascendance en tenant
     643c   calcul de la fraction de la maille concernée par l'ascendance en tenant
    558644c   compte de l'epluchage du thermique.
    559645c
    … …  
    578664            else
    579665            zmix(ig)=zlev(ig,lmix(ig))
    580             print*,'pb zmix'
     666c            print*,'pb zmix'
    581667            endif
    582668         else
    … …  
    657743      enddo
    658744     
    659       print*,'fin calcul fraca'
     745c      print*,'fin calcul fraca'
    660746c      print*,'11 OK convect8'
    661747c     print*,'Ea3 ',wa_moy
    … …  
    700786      enddo
    701787
    702       print*,'12 OK convect8'
     788c      print*,'12 OK convect8'
    703789c     print*,'WA4 ',wa_moy
    704790cc------------------------------------------------------------------
    … …  
    760846            detr(ig,l)=fm(ig,l)+entr(ig,l)-fm(ig,l+1)
    761847            if (detr(ig,l).lt.0.) then
    762                 entr(ig,l)=entr(ig,l)-detr(ig,l)
     848c                entr(ig,l)=entr(ig,l)-detr(ig,l)
     849                fm(ig,l+1)=fm(ig,l)+entr(ig,l)
    763850                detr(ig,l)=0.
    764851c     print*,'WARNING !!! detrainement negatif ',ig,l
    … …  
    831918c      enddo
    832919
    833       print*,'14 OK convect8'
     920c      print*,'14 OK convect8'
    834921c------------------------------------------------------------------
    835922c   Calculs pour les sorties
    … …  
    901988#endif
    902989123   continue
    903 c#define troisD
     990#define troisD
    904991#ifdef troisD
    905992c       if (sorties) then
    … …  
    9881075c     if(wa_moy(1,4).gt.1.e-10) stop
    9891076
    990       print*,'19 OK convect8'
     1077c      print*,'19 OK convect8'
    9911078      return
    9921079      end
    9931080
    994       subroutine dqthermcell(ngrid,nlay,ptimestep,fm,entr,masse
    995      .    ,q,dq,qa)
    996       use dimphy
    997       implicit none
    998 
    999 c=======================================================================
    1000 c
    1001 c   Calcul du transport verticale dans la couche limite en presence
    1002 c   de "thermiques" explicitement representes
    1003 c   calcul du dq/dt une fois qu'on connait les ascendances
    1004 c
    1005 c=======================================================================
    1006 
    1007 cym#include "dimensions.h"
    1008 cym#include "dimphy.h"
    1009 
    1010       integer ngrid,nlay
    1011 
    1012       real ptimestep
    1013       real masse(ngrid,nlay),fm(ngrid,nlay+1)
    1014       real entr(ngrid,nlay)
    1015       real q(ngrid,nlay)
    1016       real dq(ngrid,nlay)
    1017 
    1018       real qa(klon,klev),detr(klon,klev),wqd(klon,klev+1)
    1019 
    1020       integer ig,k
    1021 
    1022 c   calcul du detrainement
    1023 
    1024       do k=1,nlay
    1025          do ig=1,ngrid
    1026             detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
    1027          enddo
    1028       enddo
    1029 
    1030 c   calcul de la valeur dans les ascendances
    1031       do ig=1,ngrid
    1032          qa(ig,1)=q(ig,1)
    1033       enddo
    1034 
    1035       do k=2,nlay
    1036          do ig=1,ngrid
    1037             if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.
    1038      s         1.e-5*masse(ig,k)) then
    1039                qa(ig,k)=(fm(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k))
    1040      s         /(fm(ig,k+1)+detr(ig,k))
     1081      SUBROUTINE fermeture_seche(ngrid,nlay
     1082     s                ,pplay,pplev,pphi,zlev,rhobarz,f0,zpspsk
     1083     s                ,alim_star,zh,zo,lentr,lmin,nu_min,nu_max,r_aspect
     1084     s                ,zmax,wmax)
     1085
     1086      IMPLICIT NONE
     1087
     1088#include "dimensions.h"
     1089#include "dimphy.h"
     1090#include "YOMCST.h"
     1091
     1092      INTEGER ngrid,nlay
     1093      REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
     1094      real pphi(ngrid,nlay)
     1095      real zlev(klon,klev+1)
     1096      real alim_star(klon,klev)
     1097      real f0(klon)
     1098      integer lentr(klon)
     1099      integer lmin(klon)
     1100      real zmax(klon)
     1101      real wmax(klon)
     1102      real nu_min
     1103      real nu_max
     1104      real r_aspect
     1105      real rhobarz(klon,klev+1)
     1106      REAL zh(klon,klev)
     1107      real zo(klon,klev)
     1108      real zpspsk(klon,klev)
     1109
     1110      integer ig,l
     1111
     1112      real f_star(klon,klev+1)
     1113      real detr_star(klon,klev)
     1114      real entr_star(klon,klev)
     1115      real zw2(klon,klev+1)
     1116      real linter(klon)
     1117      integer lmix(klon)
     1118      integer lmax(klon)
     1119      real zlevinter(klon)
     1120      real wa_moy(klon,klev+1)
     1121      real wmaxa(klon)
     1122      REAL ztv(klon,klev)
     1123      REAL ztva(klon,klev)
     1124      real nu(klon,klev)
     1125      real zmax0_sec(klon)
     1126      save zmax0_sec
     1127
     1128      do l=1,nlay
     1129         do ig=1,ngrid
     1130      ztv(ig,l)=zh(ig,l)/zpspsk(ig,l)
     1131      ztv(ig,l)=ztv(ig,l)*(1.+RETV*zo(ig,l))
     1132         enddo
     1133      enddo
     1134      do l=1,nlay-2
     1135         do ig=1,ngrid
     1136            if (ztv(ig,l).gt.ztv(ig,l+1)
     1137     s         .and.alim_star(ig,l).gt.1.e-10
     1138     s         .and.zw2(ig,l).lt.1e-10) then
     1139               f_star(ig,l+1)=alim_star(ig,l)
     1140ctest:calcul de dteta
     1141               zw2(ig,l+1)=2.*RG*(ztv(ig,l)-ztv(ig,l+1))/ztv(ig,l+1)
     1142     s                     *(zlev(ig,l+1)-zlev(ig,l))
     1143     s                     *0.4*pphi(ig,l)/(pphi(ig,l+1)-pphi(ig,l))
     1144            else if ((zw2(ig,l).ge.1e-10).and.
     1145     s         (f_star(ig,l)+alim_star(ig,l)).gt.1.e-10) then
     1146cestimation du detrainement a partir de la geometrie du pas precedent
     1147ctests sur la definition du detr
     1148             nu(ig,l)=(nu_min+nu_max)/2.
     1149     s         *(1.-(nu_max-nu_min)/(nu_max+nu_min)
     1150     s  *tanh((((ztva(ig,l-1)-ztv(ig,l))/ztv(ig,l))/0.0005)))
     1151         
     1152             detr_star(ig,l)=rhobarz(ig,l)
     1153     s                      *sqrt(zw2(ig,l))
     1154     s                       /(r_aspect*zmax0_sec(ig))*
     1155c     s                       /(r_aspect*zmax0(ig))*
     1156     s                      (sqrt(nu(ig,l)*zlev(ig,l+1)
     1157     s                /sqrt(zw2(ig,l)))
     1158     s                     -sqrt(nu(ig,l)*zlev(ig,l)
     1159     s                /sqrt(zw2(ig,l))))
     1160         detr_star(ig,l)=detr_star(ig,l)/f0(ig)
     1161         if ((detr_star(ig,l)).gt.f_star(ig,l)) then
     1162              detr_star(ig,l)=f_star(ig,l)
     1163         endif
     1164         entr_star(ig,l)=0.9*detr_star(ig,l)
     1165             if ((l.lt.lentr(ig))) then
     1166                 entr_star(ig,l)=0.
     1167c                 detr_star(ig,l)=0.
     1168             endif
     1169            print*,'ok detr_star'
     1170cprise en compte du detrainement dans le calcul du flux
     1171             f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)
     1172     s                      -detr_star(ig,l)
     1173ctest sur le signe de f_star
     1174       if ((f_star(ig,l+1)+detr_star(ig,l)).gt.1.e-10) then
     1175cAM on melange Tl et qt du thermique
     1176          ztva(ig,l)=(f_star(ig,l)*ztva(ig,l-1)+(entr_star(ig,l)
     1177     s                    +alim_star(ig,l))
     1178     s                    *ztv(ig,l))/(f_star(ig,l+1)+detr_star(ig,l))
     1179               zw2(ig,l+1)=zw2(ig,l)*(f_star(ig,l)
     1180     s                     /(f_star(ig,l+1)+detr_star(ig,l)))**2+
     1181     s                     2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
     1182     s                     *(zlev(ig,l+1)-zlev(ig,l))
     1183            endif
     1184        endif
     1185c
     1186            if (zw2(ig,l+1).lt.0.) then
     1187               linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))
     1188     s           -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
     1189               zw2(ig,l+1)=0.
     1190               print*,'linter=',linter(ig)
    10411191            else
    1042                qa(ig,k)=q(ig,k)
    1043             endif
    1044          enddo
    1045       enddo
    1046 
    1047       do k=2,nlay
    1048          do ig=1,ngrid
    1049 c             wqd(ig,k)=fm(ig,k)*0.5*(q(ig,k-1)+q(ig,k))
    1050             wqd(ig,k)=fm(ig,k)*q(ig,k)
    1051          enddo
    1052       enddo
    1053       do ig=1,ngrid
    1054          wqd(ig,1)=0.
    1055          wqd(ig,nlay+1)=0.
    1056       enddo
    1057 
    1058       do k=1,nlay
    1059          do ig=1,ngrid
    1060             dq(ig,k)=(detr(ig,k)*qa(ig,k)-entr(ig,k)*q(ig,k)
    1061      s               -wqd(ig,k)+wqd(ig,k+1))
    1062      s               /masse(ig,k)
    1063          enddo
    1064       enddo
    1065 
    1066       return
    1067       end
    1068       subroutine dvthermcell(ngrid,nlay,ptimestep,fm,entr,masse
    1069      .    ,fraca,larga
    1070      .    ,u,v,du,dv,ua,va)
    1071        use dimphy
    1072       implicit none
    1073 
    1074 c=======================================================================
    1075 c
    1076 c   Calcul du transport verticale dans la couche limite en presence
    1077 c   de "thermiques" explicitement representes
    1078 c   calcul du dq/dt une fois qu'on connait les ascendances
    1079 c
    1080 c=======================================================================
    1081 
    1082 cym#include "dimensions.h"
    1083 cym#include "dimphy.h"
    1084 
    1085       integer ngrid,nlay
    1086 
    1087       real ptimestep
    1088       real masse(ngrid,nlay),fm(ngrid,nlay+1)
    1089       real fraca(ngrid,nlay+1)
    1090       real larga(ngrid)
    1091       real entr(ngrid,nlay)
    1092       real u(ngrid,nlay)
    1093       real ua(ngrid,nlay)
    1094       real du(ngrid,nlay)
    1095       real v(ngrid,nlay)
    1096       real va(ngrid,nlay)
    1097       real dv(ngrid,nlay)
    1098 
    1099       real qa(klon,klev),detr(klon,klev)
    1100       real wvd(klon,klev+1),wud(klon,klev+1)
    1101       real gamma0,gamma(klon,klev+1)
    1102       real dua,dva
    1103       integer iter
    1104 
    1105       integer ig,k
    1106 
    1107 c   calcul du detrainement
    1108 
    1109       do k=1,nlay
    1110          do ig=1,ngrid
    1111             detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
    1112          enddo
    1113       enddo
    1114 
    1115 c   calcul de la valeur dans les ascendances
    1116       do ig=1,ngrid
    1117          ua(ig,1)=u(ig,1)
    1118          va(ig,1)=v(ig,1)
    1119       enddo
    1120 
    1121       do k=2,nlay
    1122          do ig=1,ngrid
    1123             if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.
    1124      s         1.e-5*masse(ig,k)) then
    1125 c   On itere sur la valeur du coeff de freinage.
    1126 c              gamma0=rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k))
    1127                gamma0=masse(ig,k)
    1128      s         *sqrt( 0.5*(fraca(ig,k+1)+fraca(ig,k)) )
    1129      s         *0.5/larga(ig)
    1130 c              gamma0=0.
    1131 c   la premiere fois on multiplie le coefficient de freinage
    1132 c   par le module du vent dans la couche en dessous.
    1133                dua=ua(ig,k-1)-u(ig,k-1)
    1134                dva=va(ig,k-1)-v(ig,k-1)
    1135                do iter=1,5
    1136                   gamma(ig,k)=gamma0*sqrt(dua**2+dva**2)
    1137                   ua(ig,k)=(fm(ig,k)*ua(ig,k-1)
    1138      s               +(entr(ig,k)+gamma(ig,k))*u(ig,k))
    1139      s               /(fm(ig,k+1)+detr(ig,k)+gamma(ig,k))
    1140                   va(ig,k)=(fm(ig,k)*va(ig,k-1)
    1141      s               +(entr(ig,k)+gamma(ig,k))*v(ig,k))
    1142      s               /(fm(ig,k+1)+detr(ig,k)+gamma(ig,k))
    1143 c                 print*,k,ua(ig,k),va(ig,k),u(ig,k),v(ig,k),dua,dva
    1144                   dua=ua(ig,k)-u(ig,k)
    1145                   dva=va(ig,k)-v(ig,k)
    1146                enddo
     1192               wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
     1193            endif
     1194            if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
     1195c   lmix est le niveau de la couche ou w (wa_moy) est maximum
     1196               lmix(ig)=l+1
     1197               wmaxa(ig)=wa_moy(ig,l+1)
     1198            endif
     1199         enddo
     1200      enddo
     1201      print*,'fin calcul zw2'
     1202c
     1203c Calcul de la couche correspondant a la hauteur du thermique
     1204      do ig=1,ngrid
     1205         lmax(ig)=lentr(ig)
     1206      enddo
     1207      do ig=1,ngrid
     1208         do l=nlay,lentr(ig)+1,-1
     1209            if (zw2(ig,l).le.1.e-10) then
     1210               lmax(ig)=l-1
     1211            endif
     1212         enddo
     1213      enddo
     1214c pas de thermique si couche 1 stable
     1215      do ig=1,ngrid
     1216         if (lmin(ig).gt.1) then
     1217            lmax(ig)=1
     1218             lmin(ig)=1
     1219             lentr(ig)=1
     1220         endif
     1221      enddo
     1222c   
     1223c Determination de zw2 max
     1224      do ig=1,ngrid
     1225         wmax(ig)=0.
     1226      enddo
     1227
     1228      do l=1,nlay
     1229         do ig=1,ngrid
     1230            if (l.le.lmax(ig)) then
     1231                if (zw2(ig,l).lt.0.)then
     1232                  print*,'pb2 zw2<0'
     1233                endif
     1234                zw2(ig,l)=sqrt(zw2(ig,l))
     1235                wmax(ig)=max(wmax(ig),zw2(ig,l))
    11471236            else
    1148                ua(ig,k)=u(ig,k)
    1149                va(ig,k)=v(ig,k)
    1150                gamma(ig,k)=0.
    1151             endif
    1152          enddo
    1153       enddo
    1154 
    1155       do k=2,nlay
    1156          do ig=1,ngrid
    1157             wud(ig,k)=fm(ig,k)*u(ig,k)
    1158             wvd(ig,k)=fm(ig,k)*v(ig,k)
    1159          enddo
    1160       enddo
    1161       do ig=1,ngrid
    1162          wud(ig,1)=0.
    1163          wud(ig,nlay+1)=0.
    1164          wvd(ig,1)=0.
    1165          wvd(ig,nlay+1)=0.
    1166       enddo
    1167 
    1168       do k=1,nlay
    1169          do ig=1,ngrid
    1170             du(ig,k)=((detr(ig,k)+gamma(ig,k))*ua(ig,k)
    1171      s               -(entr(ig,k)+gamma(ig,k))*u(ig,k)
    1172      s               -wud(ig,k)+wud(ig,k+1))
    1173      s               /masse(ig,k)
    1174             dv(ig,k)=((detr(ig,k)+gamma(ig,k))*va(ig,k)
    1175      s               -(entr(ig,k)+gamma(ig,k))*v(ig,k)
    1176      s               -wvd(ig,k)+wvd(ig,k+1))
    1177      s               /masse(ig,k)
    1178          enddo
    1179       enddo
    1180 
    1181       return
    1182       end
    1183       subroutine dqthermcell2(ngrid,nlay,ptimestep,fm,entr,masse,frac
    1184      .    ,q,dq,qa)
    1185       use dimphy
    1186       implicit none
    1187 
    1188 c=======================================================================
    1189 c
    1190 c   Calcul du transport verticale dans la couche limite en presence
    1191 c   de "thermiques" explicitement representes
    1192 c   calcul du dq/dt une fois qu'on connait les ascendances
    1193 c
    1194 c=======================================================================
    1195 
    1196 cym#include "dimensions.h"
    1197 cym#include "dimphy.h"
    1198 
    1199       integer ngrid,nlay
    1200 
    1201       real ptimestep
    1202       real masse(ngrid,nlay),fm(ngrid,nlay+1)
    1203       real entr(ngrid,nlay),frac(ngrid,nlay)
    1204       real q(ngrid,nlay)
    1205       real dq(ngrid,nlay)
    1206 
    1207       real qa(klon,klev),detr(klon,klev),wqd(klon,klev+1)
    1208       real qe(klon,klev),zf,zf2
    1209 
    1210       integer ig,k
    1211 
    1212 c   calcul du detrainement
    1213 
    1214       do k=1,nlay
    1215          do ig=1,ngrid
    1216             detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
    1217          enddo
    1218       enddo
    1219 
    1220 c   calcul de la valeur dans les ascendances
    1221       do ig=1,ngrid
    1222          qa(ig,1)=q(ig,1)
    1223          qe(ig,1)=q(ig,1)
    1224       enddo
    1225 
    1226       do k=2,nlay
    1227          do ig=1,ngrid
    1228             if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.
    1229      s         1.e-5*masse(ig,k)) then
    1230                zf=0.5*(frac(ig,k)+frac(ig,k+1))
    1231                zf2=1./(1.-zf)
    1232                qa(ig,k)=(fm(ig,k)*qa(ig,k-1)+zf2*entr(ig,k)*q(ig,k))
    1233      s         /(fm(ig,k+1)+detr(ig,k)+entr(ig,k)*zf*zf2)
    1234                qe(ig,k)=(q(ig,k)-zf*qa(ig,k))*zf2
    1235             else
    1236                qa(ig,k)=q(ig,k)
    1237                qe(ig,k)=q(ig,k)
    1238             endif
    1239          enddo
    1240       enddo
    1241 
    1242       do k=2,nlay
    1243          do ig=1,ngrid
    1244 c             wqd(ig,k)=fm(ig,k)*0.5*(q(ig,k-1)+q(ig,k))
    1245             wqd(ig,k)=fm(ig,k)*qe(ig,k)
    1246          enddo
    1247       enddo
    1248       do ig=1,ngrid
    1249          wqd(ig,1)=0.
    1250          wqd(ig,nlay+1)=0.
    1251       enddo
    1252 
    1253       do k=1,nlay
    1254          do ig=1,ngrid
    1255             dq(ig,k)=(detr(ig,k)*qa(ig,k)-entr(ig,k)*qe(ig,k)
    1256      s               -wqd(ig,k)+wqd(ig,k+1))
    1257      s               /masse(ig,k)
    1258          enddo
    1259       enddo
    1260 
    1261       return
    1262       end
    1263       subroutine dvthermcell2(ngrid,nlay,ptimestep,fm,entr,masse
    1264      .    ,fraca,larga
    1265      .    ,u,v,du,dv,ua,va)
    1266       use dimphy
    1267       implicit none
    1268 
    1269 c=======================================================================
    1270 c
    1271 c   Calcul du transport verticale dans la couche limite en presence
    1272 c   de "thermiques" explicitement representes
    1273 c   calcul du dq/dt une fois qu'on connait les ascendances
    1274 c
    1275 c=======================================================================
    1276 
    1277 cym#include "dimensions.h"
    1278 cym#include "dimphy.h"
    1279 
    1280       integer ngrid,nlay
    1281 
    1282       real ptimestep
    1283       real masse(ngrid,nlay),fm(ngrid,nlay+1)
    1284       real fraca(ngrid,nlay+1)
    1285       real larga(ngrid)
    1286       real entr(ngrid,nlay)
    1287       real u(ngrid,nlay)
    1288       real ua(ngrid,nlay)
    1289       real du(ngrid,nlay)
    1290       real v(ngrid,nlay)
    1291       real va(ngrid,nlay)
    1292       real dv(ngrid,nlay)
    1293 
    1294       real qa(klon,klev),detr(klon,klev),zf,zf2
    1295       real wvd(klon,klev+1),wud(klon,klev+1)
    1296       real gamma0,gamma(klon,klev+1)
    1297       real ue(klon,klev),ve(klon,klev)
    1298       real dua,dva
    1299       integer iter
    1300 
    1301       integer ig,k
    1302 
    1303 c   calcul du detrainement
    1304 
    1305       do k=1,nlay
    1306          do ig=1,ngrid
    1307             detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
    1308          enddo
    1309       enddo
    1310 
    1311 c   calcul de la valeur dans les ascendances
    1312       do ig=1,ngrid
    1313          ua(ig,1)=u(ig,1)
    1314          va(ig,1)=v(ig,1)
    1315          ue(ig,1)=u(ig,1)
    1316          ve(ig,1)=v(ig,1)
    1317       enddo
    1318 
    1319       do k=2,nlay
    1320          do ig=1,ngrid
    1321             if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.
    1322      s         1.e-5*masse(ig,k)) then
    1323 c   On itere sur la valeur du coeff de freinage.
    1324 c              gamma0=rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k))
    1325                gamma0=masse(ig,k)
    1326      s         *sqrt( 0.5*(fraca(ig,k+1)+fraca(ig,k)) )
    1327      s         *0.5/larga(ig)
    1328      s         *1.
    1329 c    s         *0.5
    1330 c              gamma0=0.
    1331                zf=0.5*(fraca(ig,k)+fraca(ig,k+1))
    1332                zf=0.
    1333                zf2=1./(1.-zf)
    1334 c   la premiere fois on multiplie le coefficient de freinage
    1335 c   par le module du vent dans la couche en dessous.
    1336                dua=ua(ig,k-1)-u(ig,k-1)
    1337                dva=va(ig,k-1)-v(ig,k-1)
    1338                do iter=1,5
    1339 c   On choisit une relaxation lineaire.
    1340                   gamma(ig,k)=gamma0
    1341 c   On choisit une relaxation quadratique.
    1342                   gamma(ig,k)=gamma0*sqrt(dua**2+dva**2)
    1343                   ua(ig,k)=(fm(ig,k)*ua(ig,k-1)
    1344      s               +(zf2*entr(ig,k)+gamma(ig,k))*u(ig,k))
    1345      s               /(fm(ig,k+1)+detr(ig,k)+entr(ig,k)*zf*zf2
    1346      s                 +gamma(ig,k))
    1347                   va(ig,k)=(fm(ig,k)*va(ig,k-1)
    1348      s               +(zf2*entr(ig,k)+gamma(ig,k))*v(ig,k))
    1349      s               /(fm(ig,k+1)+detr(ig,k)+entr(ig,k)*zf*zf2
    1350      s                 +gamma(ig,k))
    1351 c                 print*,k,ua(ig,k),va(ig,k),u(ig,k),v(ig,k),dua,dva
    1352                   dua=ua(ig,k)-u(ig,k)
    1353                   dva=va(ig,k)-v(ig,k)
    1354                   ue(ig,k)=(u(ig,k)-zf*ua(ig,k))*zf2
    1355                   ve(ig,k)=(v(ig,k)-zf*va(ig,k))*zf2
    1356                enddo
    1357             else
    1358                ua(ig,k)=u(ig,k)
    1359                va(ig,k)=v(ig,k)
    1360                ue(ig,k)=u(ig,k)
    1361                ve(ig,k)=v(ig,k)
    1362                gamma(ig,k)=0.
    1363             endif
    1364          enddo
    1365       enddo
    1366 
    1367       do k=2,nlay
    1368          do ig=1,ngrid
    1369             wud(ig,k)=fm(ig,k)*ue(ig,k)
    1370             wvd(ig,k)=fm(ig,k)*ve(ig,k)
    1371          enddo
    1372       enddo
    1373       do ig=1,ngrid
    1374          wud(ig,1)=0.
    1375          wud(ig,nlay+1)=0.
    1376          wvd(ig,1)=0.
    1377          wvd(ig,nlay+1)=0.
    1378       enddo
    1379 
    1380       do k=1,nlay
    1381          do ig=1,ngrid
    1382             du(ig,k)=((detr(ig,k)+gamma(ig,k))*ua(ig,k)
    1383      s               -(entr(ig,k)+gamma(ig,k))*ue(ig,k)
    1384      s               -wud(ig,k)+wud(ig,k+1))
    1385      s               /masse(ig,k)
    1386             dv(ig,k)=((detr(ig,k)+gamma(ig,k))*va(ig,k)
    1387      s               -(entr(ig,k)+gamma(ig,k))*ve(ig,k)
    1388      s               -wvd(ig,k)+wvd(ig,k+1))
    1389      s               /masse(ig,k)
    1390          enddo
    1391       enddo
    1392 
    1393       return
    1394       end
     1237                 zw2(ig,l)=0.
     1238            endif
     1239          enddo
     1240      enddo
     1241
     1242c   Longueur caracteristique correspondant a la hauteur des thermiques.
     1243      do  ig=1,ngrid
     1244         zmax(ig)=0.
     1245         zlevinter(ig)=zlev(ig,1)
     1246      enddo
     1247      do  ig=1,ngrid
     1248c calcul de zlevinter
     1249          zlevinter(ig)=(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*
     1250     s    linter(ig)+zlev(ig,lmax(ig))-lmax(ig)*(zlev(ig,lmax(ig)+1)
     1251     s    -zlev(ig,lmax(ig)))
     1252cpour le cas ou on prend tjs lmin=1
     1253c       zmax(ig)=max(zmax(ig),zlevinter(ig)-zlev(ig,lmin(ig)))
     1254       zmax(ig)=max(zmax(ig),zlevinter(ig)-zlev(ig,1))
     1255       zmax0_sec(ig)=zmax(ig)
     1256      enddo
     1257       
     1258      RETURN
     1259      END

  • LMDZ4/trunk/libf/phylmd/thermcell.h

    r776 r878  
    1 !
    2 ! $Header$
    3 !
    4 
    51      integer iflag_thermals,nsplit_thermals
    62      real r_aspect_thermals,l_mix_thermals,tho_thermals
    73      integer w2di_thermals,isplit
    84
    9       common/ctherm/iflag_thermals,nsplit_thermals
    10      s              ,r_aspect_thermals,l_mix_thermals,tho_thermals
    11      s              ,w2di_thermals
     5      common/ctherm1/iflag_thermals,nsplit_thermals
     6      common/ctherm2/r_aspect_thermals,l_mix_thermals,tho_thermals
     7      common/ctherm3/w2di_thermals
    128
    13 c$OMP THREADPRIVATE(/ctherm/)

  • LMDZ4/trunk/libf/phylmd/write_histISCCP.h

    r827 r878  
    2020cIM: champ 3d : (lon,lat,pres) pour un tau fixe
    2121c
    22        CALL histwrite_phy(nid_isccp,"cldISCCP_"//taulev(k)//verticaxe(n),
     22      CALL histwrite_phy(nid_isccp,"cldISCCP_"//taulev(k)//verticaxe(n),
    2323     .                  itau_w,zx_tmp_fi3d)
    2424        ENDDO !k

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

    r782 r878  
    3030     &                   pctsrf_new(:,is_ter))
    3131c
     32
     33      CALL histwrite_phy(nid_day,"dthmin",itau_w,dthmin)
     34      CALL histwrite_phy(nid_day,"weakinv",itau_w,weak_inversion)
     35
     36
    3237cym      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zxtsol,zx_tmp_2d)
    3338      CALL histwrite_phy(nid_day,"tsol",itau_w,zxtsol)
    … …  
    486491     $                     zx_tmp_fi2d)
    487492C
     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
    488503      END DO 
    489504c=================================================================

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

    r864 r878  
    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
    … …  
    951978      CALL histwrite_phy(nid_mth,"dtcon",itau_w,zx_tmp_fi3d)
    952979c
    953 c temperature tendency due to dry convective processes
    954        DO l=1, klev
    955        DO i=1, klon
    956         zx_tmp_fi3d(i,l)=d_t_ajs(i,l)/pdtphys
    957        ENDDO !i
    958        ENDDO !l
    959 c
    960 cym      CALL gr_fi_ecrit(klev, klon,iim,jjmp1, zx_tmp_fi3d,zx_tmp_3d)
    961       CALL histwrite_phy(nid_mth,"dtajs",itau_w,zx_tmp_fi3d)
     980
     981
     982
     983
    962984c
    963985c  temperature tendency due to large scale precipitation
    … …  
    16901712      IF(lev_histmth.GE.4) THEN
    16911713c
     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
    16921732cym      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, clwcon0, zx_tmp_3d)
    16931733      CALL histwrite_phy(nid_mth,"clwcon",itau_w,clwcon0)
    … …  
    17561796cym      CALL gr_fi_ecrit(klev,klon,iim,jjmp1, ratqs, zx_tmp_3d)
    17571797      CALL histwrite_phy(nid_mth,"ratqs",itau_w,ratqs)
     1798cIM: 101003 : K/30min ==> K/s
     1799      zx_tmp_fi3d(1:klon,1:klev)=d_t_ajsb(1:klon,1:klev)/pdtphys
     1800cym      CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
     1801      CALL histwrite_phy(nid_mth,"dtajs",itau_w,zx_tmp_fi3d)
     1802c
     1803      zx_tmp_fi3d(1:klon,1:klev)=d_q_ajsb(1:klon,1:klev)/pdtphys
     1804cym      CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
     1805      CALL histwrite_phy(nid_mth,"dqajs",itau_w,zx_tmp_fi3d)
     1806c
     1807cIM: 101003 : K/30min ==> K/s
     1808      zx_tmp_fi3d(1:klon,1:klev)=d_t_ajs(1:klon,1:klev)/pdtphys
     1809cym      CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
     1810      CALL histwrite_phy(nid_mth,"dtthe",itau_w,zx_tmp_fi3d)
    17581811c
    17591812      zx_tmp_fi3d(1:klon,1:klev)=d_q_ajs(1:klon,1:klev)/pdtphys
    17601813cym      CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
    1761       CALL histwrite_phy(nid_mth,"dqajs",itau_w,zx_tmp_fi3d)
     1814      CALL histwrite_phy(nid_mth,"dqthe",itau_w,zx_tmp_fi3d)
     1815c
    17621816c
    17631817cIM: 101003 : K/day ==> K/s

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

    r789 r878  
    429429                                                          enddo
    430430
    431        print*,'pblhmin ',pblhmin
     431!      print*,'pblhmin ',pblhmin
    432432CTest a remettre 21 11 02
    433433c test abd 13 05 02      if(0.eq.1) then
    … …  
    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.