Ignore:
Timestamp:
Dec 1, 2016, 3:34:58 PM (8 years ago)
Author:
fhourdin
Message:

Modification de yamada4 pour corriger quelques erreurs et
autoriser de jouer avec les fonctions de stabilité.
Par defaut
new_yamada4=n donne les mêmes resultats numeriques qu'avant.
new_yamada4=y permet d'activer une nombre de Richardson critique
a partir du quel on ne fait plus décroitre les fonction de
stabilité Sm et Alpha de Mellor et Yamada.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/phylmd/yamada4.F90

    r2574 r2721  
    22
    33SUBROUTINE yamada4(ni, nsrf, ngrid, dt, g, rconst, plev, temp, zlev, zlay, u, v, teta, &
    4     cd, q2, km, kn, kq, ustar, iflag_pbl)
     4    cd, tke, km, kn, kq, ustar, iflag_pbl)
    55
    66  USE dimphy
    … …  
    5656  !             iflag_pbl=11 -> the model starts with NP from start files created by ce0l
    5757  !                          -> the model can run with longer time-steps.
     58  !             2016/11/30 (EV etienne.vignon@univ-grenoble-alpes.fr)
     59  !               On met tke (=q2/2) en entrée plutôt que q2
     60  !               On corrige l'update de la tke
    5861  !
    5962  ! Inpout/Output :
    6063  !==============
    61   ! q2 : $q^2$ au bas de chaque couche
     64  ! tke : tke au bas de chaque couche
    6265  ! (en entree : la valeur au debut du pas de temps)
    6366  ! (en sortie : la valeur a la fin du pas de temps)
    … …  
    9093  REAL teta(klon, klev)
    9194  REAL cd(klon)
    92   REAL q2(klon, klev+1)
     95  REAL tke(klon, klev+1)
    9396  REAL unsdz(klon, klev)
    9497  REAL unsdzdec(klon, klev+1)
    … …  
    104107  INCLUDE "clesphys.h"
    105108
    106 
     109  REAL q2(klon, klev+1)
    107110  REAL kmpre(klon, klev+1), tmp2, qpre
    108111  REAL mpre(klon, klev+1)
    … …  
    127130  REAL zz(klon, klev+1)
    128131  INTEGER iter
    129   REAL ric, rifc, b1, kap
    130   SAVE ric, rifc, b1, kap
    131   DATA ric, rifc, b1, kap/0.195, 0.191, 16.6, 0.4/
    132   !$OMP THREADPRIVATE(ric,rifc,b1,kap)
    133   REAL seuilsm, seuilalpha
     132  REAL,SAVE :: ric0,ric,rifc, b1, kap
     133  !$OMP THREADPRIVATE(ric0,ric,rifc,b1,kap)
     134  DATA b1, kap/16.6, 0.4/
     135  REAL,SAVE :: seuilsm, seuilalpha
     136  !$OMP THREADPRIVATE(seuilsm, seuilalpha)
    134137  REAL,SAVE :: lmixmin
    135138  !$OMP THREADPRIVATE(lmixmin)
     139  LOGICAL, SAVE :: new_yamada4
     140  !$OMP THREADPRIVATE(new_yamada4)
     141  REAL, SAVE :: yun,ydeux
     142  !$OMP THREADPRIVATE(yun,ydeux)
    136143  REAL frif, falpha, fsm
    137144  REAL rino(klon, klev+1), smyam(klon, klev), styam(klon, klev), &
    … …  
    141148
    142149
     150
    143151  ! Fonctions utilis??es
    144152  !--------------------
    … …  
    150158
    151159  IF (firstcall) THEN
     160! Seuil dans le code de turbulence
     161    new_yamada4=.false.
     162    CALL getin_p('new_yamada4',new_yamada4)
     163! Pour garantir la continuite dans la mise au point de CMIP6.
     164! Eliminer l'option new_yamada4 des le printemps 2016
     165    IF (new_yamada4) THEN
     166       ric=0.143 ! qui donne des valeurs proches des seuils proposes
     167                 ! dans YAMADA 1983 : sm=0.0845 (0.085 dans Y83)
     168                 !                    sm=1.1213 (1.12  dans Y83)
     169       CALL getin_p('yamada4_ric',ric)
     170       ric0=0.19489      ! ric=0.195 originalement, mais produisait sm<0
     171       ric=min(ric,ric0) ! Au delà de ric0, sm devient négatif
     172       rifc=frif(ric)
     173       seuilsm=fsm(frif(ric))
     174       seuilalpha=falpha(frif(ric))
     175       yun=1.
     176       ydeux=2.
     177!      yun=2.
     178!      ydeux=1.
     179    ELSE
     180  !!  DATA ric, rifc, b1, kap/0.195, 0.191, 16.6, 0.4/
     181       ric=0.195
     182       rifc=0.191
     183       seuilalpha=1.12
     184       seuilsm=0.085
     185       yun=2.
     186       ydeux=1.
     187    ENDIF
     188    PRINT*,'YAMADA4 RIc, RIfc, Sm_min, Alpha_min',ric,rifc,seuilsm,seuilalpha
    152189    firstcall = .FALSE.
    153190    lmixmin=1.
    … …  
    169206  END IF
    170207
    171 ! Seuil dans le code de turbulence
    172  seuilalpha=1.12
    173  seuilsm=0.085
    174208
    175209  nlay = klev
    … …  
    231265        rif(ig, k) = rifc
    232266      END IF
    233 
     267if (new_yamada4) then
     268        alpha(ig, k) = max(falpha(rif(ig,k)),seuilalpha)
     269        sm(ig, k) = max(fsm(rif(ig,k)),seuilsm)
     270else
    234271      IF (rif(ig,k)<0.16) THEN
    235272        alpha(ig, k) = falpha(rif(ig,k))
    … …  
    239276        sm(ig, k) = seuilsm
    240277      END IF
     278
     279end if
    241280      zz(ig, k) = b1*m2(ig, k)*(1.-rif(ig,k))*sm(ig, k)
    242281    END DO
    … …  
    244283
    245284
     285
     286
     287
     288  !=======================================================================
     289  !     DIFFERENT TYPES  DE SCHEMA de  YAMADA
     290  !=======================================================================
     291
     292  ! On commence par calculé q2 à partir de la tke
     293
     294  IF (new_yamada4) THEN
     295      DO k=1,klev+1
     296         q2(1:ngrid,k)=tke(1:ngrid,k)*ydeux
     297      ENDDO
     298  ELSE
     299      DO k=1,klev+1
     300         q2(1:ngrid,k)=tke(1:ngrid,k)
     301      ENDDO
     302  ENDIF
    246303
    247304! ====================================================================
    … …  
    252309  CALL mixinglength(ni,nsrf,ngrid,iflag_pbl,pbl_lmixmin_alpha,lmixmin,zlay,zlev,u,v,q2,n2, l)
    253310
    254 
    255 
    256   !=======================================================================
    257   !     DIFFERENT TYPES  DE SCHEMA de  YAMADA
    258   !=======================================================================
    259311
    260312  !--------------
    … …  
    350402    DO k = 2, klev - 1
    351403      km(1:ngrid, k) = l(1:ngrid, k)*sqrt(q2(1:ngrid,k))*sm(1:ngrid, k)
    352       q2(1:ngrid, k) = q2(1:ngrid, k) + dt*km(1:ngrid, k)*m2(1:ngrid, k)*(1.-rif(1:ngrid,k))
     404      q2(1:ngrid, k) = q2(1:ngrid, k) + ydeux*dt*km(1:ngrid, k)*m2(1:ngrid, k)*(1.-rif(1:ngrid,k))
     405!     q2(1:ngrid, k) = q2(1:ngrid, k) + dt*km(1:ngrid, k)*m2(1:ngrid, k)*(1.-rif(1:ngrid,k))
    353406      q2(1:ngrid, k) = min(max(q2(1:ngrid,k),1.E-10), 1.E4)
    354       q2(1:ngrid, k) = 1./(1./sqrt(q2(1:ngrid,k))+dt/(2*l(1:ngrid,k)*b1))
     407       q2(1:ngrid, k) = 1./(1./sqrt(q2(1:ngrid,k))+dt/(yun*l(1:ngrid,k)*b1))
     408!     q2(1:ngrid, k) = 1./(1./sqrt(q2(1:ngrid,k))+dt/(2*l(1:ngrid,k)*b1))
    355409      q2(1:ngrid, k) = q2(1:ngrid, k)*q2(1:ngrid, k)
    356410    END DO
    … …  
    517571      lyam(1:ngrid, 2:klev)
    518572  END IF
     573
     574
     575
     576!============================================================================
     577! Mise à jour de la tke
     578!============================================================================
     579
     580  IF (new_yamada4) THEN
     581     DO k=1,klev+1
     582        tke(1:ngrid,k)=q2(1:ngrid,k)/ydeux
     583     ENDDO
     584  ELSE
     585     DO k=1,klev+1
     586        tke(1:ngrid,k)=q2(1:ngrid,k)
     587     ENDDO
     588  ENDIF
     589
    519590
    520591!============================================================================
Note: See TracChangeset for help on using the changeset viewer.