Changeset 1738


Ignore:
Timestamp:
Apr 2, 2013, 11:48:33 AM (11 years ago)
Author:
idelkadi
Message:

Modifications for numerical stability of the boundary layer parameterizations.
Concerning yamada4.F :

  • option with asymptotic mixing length l0 imposed (iflag_pbl=8 and 9)
  • option with a new temporal scheme (iflag_pbl=10 and 11)
  • Correction for very stable PBLs (iflag_pbl=10 and 11) iflag_pbl=8 converges numerically with NPv3.1 iflag_pbl=11 -> now starts with NP from start files created by ce0l

-> the model can run with longer time-steps.

Concerning thermals :

Introduction of an implicit computation of vertical advection in
the environment of thermal plumes in thermcell_dq
impl = 0 : explicit, 1 : implicit, -1 : old version
controled by iflag_thermals =

15, 16 run with impl=-1 : numerical convergence with NPv3
17, 18 run with impl=1 : more stable

15 and 17 correspond to the activation of the stratocumulus "bidouille"

Modified routines (phylmd/):
calltherm.F90 : for managing the various options of thermcell_dq
coef_diff_turb_mod.F90 : yamada4 called for iflag_pbl<= 18 instead of 11
physiq.F : desactivation of the vertical diffusion of TKE for iflag_pbl=10
thermcellV0_main.F90 : calling thermcell_dq with implicit scheme
thermcell_dq.F90 : thermcell_dq with optional explicit/implicit scheme
thermcell_main.F90 : various option for vertical transport by thermals
yamada4.F : Yamada scheme with new options

Location:
LMDZ5/trunk/libf/phylmd
Files:
7 edited

Legend:

Unmodified
Added
Removed

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

    r1638 r1738  
    234234     &      ,tau_thermals,Ale,Alp,lalim_conv,wght_th &
    235235     &      ,zmax0,f0,zw2,fraca)
    236           else if (iflag_thermals==15.or.iflag_thermals==16) then
     236          else if (iflag_thermals>=15.and.iflag_thermals<=18) then
    237237
    238238!            print*,'THERM iflag_thermas_ed=',iflag_thermals_ed
    … …  
    271271! fait bien ce qu'on croit.
    272272
    273        flag_bidouille_stratocu=iflag_thermals<=12.or.iflag_thermals==14.or.iflag_thermals==16
     273       flag_bidouille_stratocu=iflag_thermals<=12.or.iflag_thermals==14.or.iflag_thermals==16.or.iflag_thermals==18
    274274
    275275      if (iflag_thermals<=12) then

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

    r1604 r1738  
    158158         
    159159!   iflag_pbl peut etre utilise comme longuer de melange
    160        IF (iflag_pbl.GE.11) THEN
     160       IF (iflag_pbl.GE.18) THEN
    161161          CALL vdif_kcay(knon,dtime,RG,RD,ypaprs,yt, &
    162162               yzlev,yzlay,yu,yv,yteta, &

  • LMDZ5/trunk/libf/phylmd/physiq.F

    r1737 r1738  
    26752675! Transport de la TKE par les panaches thermiques.
    26762676! FH : 2010/02/01
    2677       if (iflag_pbl.eq.10) then
    2678       call thermcell_dtke(klon,klev,nbsrf,pdtphys,fm_therm,entr_therm,
    2679      s           rg,paprs,pbl_tke)
    2680       endif
     2677!     if (iflag_pbl.eq.10) then
     2678!     call thermcell_dtke(klon,klev,nbsrf,pdtphys,fm_therm,entr_therm,
     2679!    s           rg,paprs,pbl_tke)
     2680!     endif
    26812681! ----------------------------------------------------------------------
    26822682!IM/FH: 2011/02/23

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

    r1403 r1738  
    519519!------------------------------------------------------------------
    520520
    521       call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse,  &
     521      call thermcell_dq(ngrid,nlay,1,ptimestep,fm0,entr0,masse,  &
    522522     &                    zthl,zdthladj,zta,lev_out)
    523       call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse,  &
     523      call thermcell_dq(ngrid,nlay,1,ptimestep,fm0,entr0,masse,  &
    524524     &                   po,pdoadj,zoa,lev_out)
    525525
    … …  
    561561
    562562! calcul purement conservatif pour le transport de V
    563          call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse  &
     563         call thermcell_dq(ngrid,nlay,1,ptimestep,fm0,entr0,masse  &
    564564     &    ,zu,pduadj,zua,lev_out)
    565          call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse  &
     565         call thermcell_dq(ngrid,nlay,1,ptimestep,fm0,entr0,masse  &
    566566     &    ,zv,pdvadj,zva,lev_out)
    567567      endif

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

    r1403 r1738  
    1       subroutine thermcell_dq(ngrid,nlay,ptimestep,fm,entr,  &
     1      subroutine thermcell_dq(ngrid,nlay,impl,ptimestep,fm,entr,  &
    22     &           masse,q,dq,qa,lev_out)
    33      implicit none
    … …  
    1010!   calcul du dq/dt une fois qu'on connait les ascendances
    1111!
     12! Modif 2013/01/04 (FH hourdin@lmd.jussieu.fr)
     13!  Introduction of an implicit computation of vertical advection in
     14!  the environment of thermal plumes in thermcell_dq
     15!  impl =     0 : explicit, 1 : implicit, -1 : old version
     16!
    1217!=======================================================================
    1318
    14       integer ngrid,nlay
     19      integer ngrid,nlay,impl
    1520
    1621      real ptimestep
    … …  
    2833      real cfl
    2934
    30       real qold(ngrid,nlay)
    31       real ztimestep
     35      real qold(ngrid,nlay),fqa(ngrid,nlay+1)
    3236      integer niter,iter
    3337      CHARACTER (LEN=20) :: modname='thermcell_dq'
    … …  
    3539
    3640
     41! Old explicite scheme
     42      if (impl==-1) then
     43         call thermcell_dq_o(ngrid,nlay,ptimestep,fm,entr,  &
     44     &           masse,q,dq,qa,lev_out)
     45         return
     46      endif
    3747
    3848! Calcul du critere CFL pour l'advection dans la subsidence
    … …  
    5060      enddo
    5161
    52 !IM 090508     print*,'CFL CFL CFL CFL ',cfl
    53 
    54 #undef CFL
    55 #ifdef CFL
    56 ! On subdivise le calcul en niter pas de temps.
    57       niter=int(cfl)+1
    58 #else
    59       niter=1
    60 #endif
    61 
    62       ztimestep=ptimestep/niter
    6362      qold=q
    6463
    6564
    66 do iter=1,niter
    6765      if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
    6866
    … …  
    8886      enddo
    8987
     88! Computation of tracer concentrations in the ascending plume
     89      do ig=1,ngrid
     90         qa(ig,1)=q(ig,1)
     91      enddo
     92
     93      do k=2,nlay
     94         do ig=1,ngrid
     95            if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.  &
     96     &         1.e-5*masse(ig,k)) then
     97         qa(ig,k)=(fm(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k))  &
     98     &         /(fm(ig,k+1)+detr(ig,k))
     99            else
     100               qa(ig,k)=q(ig,k)
     101            endif
     102            if (qa(ig,k).lt.0.) then
     103!               print*,'qa<0!!!'
     104            endif
     105            if (q(ig,k).lt.0.) then
     106!               print*,'q<0!!!'
     107            endif
     108         enddo
     109      enddo
     110
     111! Plume vertical flux
     112      do k=2,nlay-1
     113         fqa(:,k)=fm(:,k)*qa(:,k-1)
     114      enddo
     115      fqa(:,1)=0. ; fqa(:,nlay)=0.
     116
     117
     118! Trace species evolution
     119   if (impl==0) then
     120      do k=1,nlay-1
     121         q(:,k)=q(:,k)+(fqa(:,k)-fqa(:,k+1)-fm(:,k)*q(:,k)+fm(:,k+1)*q(:,k+1)) &
     122     &               *ptimestep/masse(:,k)
     123      enddo
     124   else
     125      do k=nlay-1,1,-1
     126         q(:,k)=(masse(:,k)*q(:,k)/ptimestep+fqa(:,k)-fqa(:,k+1)+fm(:,k+1)*q(:,k+1)) &
     127     &               /(fm(:,k)+masse(:,k)/ptimestep)
     128      enddo
     129   endif
     130
     131! Tendencies
     132      do k=1,nlay
     133         do ig=1,ngrid
     134            dq(ig,k)=(q(ig,k)-qold(ig,k))/ptimestep
     135            q(ig,k)=qold(ig,k)
     136         enddo
     137      enddo
     138
     139return
     140end
     141
     142!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     143! Obsolete version kept for convergence with Cmip5 NPv3.1 simulations
     144!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     145
     146      subroutine thermcell_dq_o(ngrid,nlay,ptimestep,fm,entr,  &
     147     &           masse,q,dq,qa,lev_out)
     148      implicit none
     149
     150#include "iniprint.h"
     151!=======================================================================
     152!
     153!   Calcul du transport verticale dans la couche limite en presence
     154!   de "thermiques" explicitement representes
     155!   calcul du dq/dt une fois qu'on connait les ascendances
     156!
     157!=======================================================================
     158
     159      integer ngrid,nlay
     160
     161      real ptimestep
     162      real masse(ngrid,nlay),fm(ngrid,nlay+1)
     163      real entr(ngrid,nlay)
     164      real q(ngrid,nlay)
     165      real dq(ngrid,nlay)
     166      integer lev_out                           ! niveau pour les print
     167
     168      real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
     169
     170      real zzm
     171
     172      integer ig,k
     173      real cfl
     174
     175      real qold(ngrid,nlay)
     176      real ztimestep
     177      integer niter,iter
     178      CHARACTER (LEN=20) :: modname='thermcell_dq'
     179      CHARACTER (LEN=80) :: abort_message
     180
     181
     182
     183! Calcul du critere CFL pour l'advection dans la subsidence
     184      cfl = 0.
     185      do k=1,nlay
     186         do ig=1,ngrid
     187            zzm=masse(ig,k)/ptimestep
     188            cfl=max(cfl,fm(ig,k)/zzm)
     189            if (entr(ig,k).gt.zzm) then
     190               print*,'entr dt > m ',entr(ig,k)*ptimestep,masse(ig,k)
     191               abort_message = ''
     192               CALL abort_gcm (modname,abort_message,1)
     193            endif
     194         enddo
     195      enddo
     196
     197!IM 090508     print*,'CFL CFL CFL CFL ',cfl
     198
     199#undef CFL
     200#ifdef CFL
     201! On subdivise le calcul en niter pas de temps.
     202      niter=int(cfl)+1
     203#else
     204      niter=1
     205#endif
     206
     207      ztimestep=ptimestep/niter
     208      qold=q
     209
     210
     211do iter=1,niter
     212      if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
     213
     214!   calcul du detrainement
     215      do k=1,nlay
     216         do ig=1,ngrid
     217            detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
     218!           print*,'Q2 DQ ',detr(ig,k),fm(ig,k),entr(ig,k)
     219!test
     220            if (detr(ig,k).lt.0.) then
     221               entr(ig,k)=entr(ig,k)-detr(ig,k)
     222               detr(ig,k)=0.
     223!               print*,'detr2<0!!!','ig=',ig,'k=',k,'f=',fm(ig,k),
     224!     s         'f+1=',fm(ig,k+1),'e=',entr(ig,k),'d=',detr(ig,k)
     225            endif
     226            if (fm(ig,k+1).lt.0.) then
     227!               print*,'fm2<0!!!'
     228            endif
     229            if (entr(ig,k).lt.0.) then
     230!               print*,'entr2<0!!!'
     231            endif
     232         enddo
     233      enddo
     234
    90235!   calcul de la valeur dans les ascendances
    91236      do ig=1,ngrid

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

    r1638 r1738  
    2222
    2323      USE dimphy
     24      USE ioipsl
    2425      USE comgeomphy , ONLY:rlond,rlatd
    2526      IMPLICIT NONE
    … …  
    4445!     4. un detrainement
    4546!
     47! Modif 2013/01/04 (FH hourdin@lmd.jussieu.fr)
     48!    Introduction of an implicit computation of vertical advection in
     49!    the environment of thermal plumes in thermcell_dq
     50!    impl =     0 : explicit, 1 : implicit, -1 : old version
     51!    controled by iflag_thermals =
     52!       15, 16 run with impl=-1 : numerical convergence with NPv3
     53!       17, 18 run with impl=1  : more stable
     54!    15 and 17 correspond to the activation of the stratocumulus "bidouille"
     55!
    4656!=======================================================================
     57
    4758
    4859!-----------------------------------------------------------------------
    … …  
    7990
    8091      integer icount
     92
     93      integer, save :: dvdq=1,dqimpl=-1
     94!$OMP THREADPRIVATE(dvdq,dqimpl)
    8195      data icount/0/
    8296      save icount
    … …  
    247261
    248262      if (debut)  then
     263!        call getin('dvdq',dvdq)
     264!        call getin('dqimpl',dqimpl)
     265
     266         if (iflag_thermals==15.or.iflag_thermals==16) then
     267            dvdq=0
     268            dqimpl=-1
     269         else
     270            dvdq=1
     271            dqimpl=1
     272         endif
     273
    249274         fm0=0.
    250275         entr0=0.
    … …  
    593618!------------------------------------------------------------------
    594619
    595       call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse,  &
     620      call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse,  &
    596621     &                    zthl,zdthladj,zta,lev_out)
    597       call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse,  &
     622      call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse,  &
    598623     &                   po,pdoadj,zoa,lev_out)
    599624
    … …  
    620645
    621646!IM 090508 
    622       if (1.eq.1) then
    623 !IM 070508 vers. _dq       
    624 !     if (1.eq.0) then
    625 
     647      if (dvdq == 0 ) then
    626648
    627649! Calcul du transport de V tenant compte d'echange par gradient
    … …  
    629651
    630652         call thermcell_dv2(ngrid,nlay,ptimestep,fm0,entr0,masse  &
    631      &    ,fraca,zmax  &
     653!    &    ,fraca*dvdq,zmax &
     654     &    ,fraca,zmax &
    632655     &    ,zu,zv,pduadj,pdvadj,zua,zva,lev_out)
    633656
    … …  
    635658
    636659! calcul purement conservatif pour le transport de V
    637          call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse  &
     660         call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse  &
    638661     &    ,zu,pduadj,zua,lev_out)
    639          call thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,masse  &
     662         call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse  &
    640663     &    ,zv,pdvadj,zva,lev_out)
     664
    641665      endif
    642666

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

    r1403 r1738  
    3737c    iflag_pbl=6 : MY 2.0
    3838c    iflag_pbl=7 : MY 2.0.Fournier
    39 c    iflag_pbl=8 : MY 2.5
    40 c    iflag_pbl>=9 : MY 2.5 avec diffusion verticale
    41 
    42 c.......................................................................
     39c    iflag_pbl=8/9 : MY 2.5
     40c       iflag_pbl=8 with special obsolete treatments for convergence
     41c       with Cmpi5 NPv3.1 simulations
     42c    iflag_pbl=10/11 :  New scheme M2 and N2 explicit and dissiptation exact
     43c    iflag_pbl=12 = 11 with vertical diffusion off q2
     44c
     45c  2013/04/01 (FH hourdin@lmd.jussieu.fr)
     46c     Correction for very stable PBLs (iflag_pbl=10 and 11)
     47c     iflag_pbl=8 converges numerically with NPv3.1
     48c     iflag_pbl=11 -> the model starts with NP from start files created by ce0l
     49c                  -> the model can run with longer time-steps.
     50c.......................................................................
     51
    4352      REAL dt,g,rconst
    4453      real plev(klon,klev+1),temp(klon,klev)
    … …  
    6372      real aa(klon,klev+1),aa0,aa1
    6473      integer iflag_pbl,ngrid
    65 
    66 
    6774      integer nlay,nlev
    6875
    … …  
    118125
    119126
    120       if (.not.(iflag_pbl.ge.6.and.iflag_pbl.le.10)) then
     127      if (.not.(iflag_pbl.ge.6.and.iflag_pbl.le.12)) then
    121128           stop'probleme de coherence dans appel a MY'
    122129      endif
    123130
    124131      ipas=ipas+1
    125       if (0.eq.1.and.first) then
    126       do ig=1,1000
    127          ri=(ig-800.)/500.
    128          if (ri.lt.ric) then
    129             zrif=frif(ri)
    130          else
    131             zrif=rifc
    132          endif
    133          if(zrif.lt.0.16) then
    134             zalpha=falpha(zrif)
    135             zsm=fsm(zrif)
    136          else
    137             zalpha=1.12
    138             zsm=0.085
    139          endif
    140 c     print*,ri,rif,zalpha,zsm
    141       enddo
    142       endif
     132
    143133
    144134c.......................................................................
    … …  
    173163                                                     ENDDO
    174164c
    175 c.......................................................................
     165!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     166! Computing M^2, N^2, Richardson numbers, stability functions
     167!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    176168
    177169      do k=2,klev
    … …  
    197189         endif
    198190         zz(ig,k)=b1*m2(ig,k)*(1.-rif(ig,k))*sm(ig,k)
    199 c     print*,'RIF L=',k,rif(ig,k),ri*alpha(ig,k)
    200 
    201 
    202                                                           enddo
    203       enddo
    204 
    205 
    206 c====================================================================
    207 c   Au premier appel, on determine l et q2 de facon iterative.
    208 c iterration pour determiner la longueur de melange
    209 
    210 
    211       if (first.or.iflag_pbl.eq.6) then
    212                                                           do ig=1,ngrid
    213       l0(ig)=10.
    214                                                           enddo
    215       do k=2,klev-1
    216                                                           do ig=1,ngrid
    217         l(ig,k)=l0(ig)*kap*zlev(ig,k)/(kap*zlev(ig,k)+l0(ig))
    218                                                           enddo
    219       enddo
    220 
    221       do iter=1,10
    222                                                           do ig=1,ngrid
    223          sq(ig)=1.e-10
    224          sqz(ig)=1.e-10
    225                                                           enddo
    226          do k=2,klev-1
    227                                                           do ig=1,ngrid
    228            q2(ig,k)=l(ig,k)**2*zz(ig,k)
    229            l(ig,k)=fl(zlev(ig,k),l0(ig),q2(ig,k),n2(ig,k))
    230            zq=sqrt(q2(ig,k))
    231            sqz(ig)=sqz(ig)+zq*zlev(ig,k)*(zlay(ig,k)-zlay(ig,k-1))
    232            sq(ig)=sq(ig)+zq*(zlay(ig,k)-zlay(ig,k-1))
    233                                                           enddo
    234          enddo
    235                                                           do ig=1,ngrid
    236          l0(ig)=0.2*sqz(ig)/sq(ig)
    237 c        l0(ig)=30.
    238                                                           enddo
    239 c      print*,'ITER=',iter,'  L0=',l0
    240 
    241       enddo
    242 
    243 c     print*,'Fin de l initialisation de q2 et l0'
    244 
    245       endif ! first
    246 
    247 c====================================================================
    248 c  Calcul de la longueur de melange.
     191                                                          enddo
     192      enddo
     193
     194
     195c====================================================================
     196c  Computing the mixing length
    249197c====================================================================
    250198
    251199c   Mise a jour de l0
     200      if (iflag_pbl==8.or.iflag_pbl==10) then
     201
     202!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     203! Iterative computation of l0
     204! This version is kept for iflag_pbl only for convergence
     205! with NPv3.1 Cmip5 simulations
     206!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     207
    252208                                                          do ig=1,ngrid
    253209      sq(ig)=1.e-10
    … …  
    263219                                                          do ig=1,ngrid
    264220      l0(ig)=0.2*sqz(ig)/sq(ig)
    265 c        l0(ig)=30.
    266                                                           enddo
    267 c      print*,'ITER=',iter,'  L0=',l0
    268 c   calcul de l(z)
     221                                                          enddo
    269222      do k=2,klev
    270223                                                          do ig=1,ngrid
    271224         l(ig,k)=fl(zlev(ig,k),l0(ig),q2(ig,k),n2(ig,k))
    272          if(first) then
    273            q2(ig,k)=l(ig,k)**2*zz(ig,k)
    274          endif
    275                                                           enddo
    276       enddo
     225                                                          enddo
     226      enddo
     227!     print*,'L0 cas 8 ou 10 ',l0
     228
     229      else
     230
     231!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     232! In all other case, the assymptotic mixing length l0 is imposed (100m)
     233!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     234
     235          l0(:)=150.
     236          do k=2,klev
     237                                                          do ig=1,ngrid
     238             l(ig,k)=fl(zlev(ig,k),l0(ig),q2(ig,k),n2(ig,k))
     239                                                          enddo
     240          enddo
     241!     print*,'L0 cas autres ',l0
     242
     243      endif
     244
    277245
    278246c====================================================================
    … …  
    282250
    283251      do k=2,klev
    284                                                           do ig=1,ngrid
    285          q2(ig,k)=l(ig,k)**2*zz(ig,k)
    286                                                           enddo
     252         q2(:,k)=l(:,k)**2*zz(:,k)
    287253      enddo
    288254
    … …  
    342308      enddo
    343309
    344       else if (iflag_pbl.ge.8) then
     310      else if (iflag_pbl==8.or.iflag_pbl==9) then
    345311c====================================================================
    346312c   Yamada 2.5 a la Didi
    … …  
    366332c     print*,'0L=',k,l(ig,k),delta(ig,k),km(ig,k)
    367333         qpre=sqrt(q2(ig,k))
    368          if (iflag_pbl.eq.8 ) then
     334!        if (iflag_pbl.eq.8 ) then
    369335            if (aa(ig,k).gt.0.) then
    370336               q2(ig,k)=(qpre+aa(ig,k)*qpre*qpre)**2
    … …  
    372338               q2(ig,k)=(qpre/(1.-aa(ig,k)*qpre))**2
    373339            endif
    374          else ! iflag_pbl=9
    375             if (aa(ig,k)*qpre.gt.0.9) then
    376                q2(ig,k)=(qpre*10.)**2
    377             else
    378                q2(ig,k)=(qpre/(1.-aa(ig,k)*qpre))**2
    379             endif
    380          endif
     340!        else ! iflag_pbl=9
     341!           if (aa(ig,k)*qpre.gt.0.9) then
     342!              q2(ig,k)=(qpre*10.)**2
     343!           else
     344!              q2(ig,k)=(qpre/(1.-aa(ig,k)*qpre))**2
     345!           endif
     346!        endif
    381347         q2(ig,k)=min(max(q2(ig,k),1.e-10),1.e4)
    382348c     print*,'Q2 L=',k,q2(ig,k),qpre*qpre
    383349                                                          enddo
    384350      enddo
     351
     352      else if (iflag_pbl>=10) then
     353
     354!        print*,'Schema mixte D'
     355!        print*,'Longueur ',l(:,:)
     356         do k=2,klev-1
     357            l(:,k)=max(l(:,k),1.)
     358            km(:,k)=l(:,k)*sqrt(q2(:,k))*sm(:,k)
     359            q2(:,k)=q2(:,k)+dt*km(:,k)*m2(:,k)*(1.-rif(:,k))
     360            q2(:,k)=min(max(q2(:,k),1.e-10),1.e4)
     361            q2(:,k)=1./(1./sqrt(q2(:,k))+dt/(2*l(:,k)*b1))
     362            q2(:,k)=q2(:,k)*q2(:,k)
     363         enddo
     364
     365
     366      else
     367         stop'Cas nom prevu dans yamada4'
    385368
    386369      endif ! Fin du cas 8
    … …  
    404387
    405388! Transport diffusif vertical de la TKE.
    406       if (iflag_pbl.ge.9) then
     389      if (iflag_pbl.ge.12) then
    407390!       print*,'YAMADA VDIF'
    408391        q2(:,1)=q2(:,2)
    … …  
    425408                                                          enddo
    426409
    427 !      print*,'pblhmin ',pblhmin
     410!     print*,'pblhmin ',pblhmin
    428411CTest a remettre 21 11 02
    429412c test abd 13 05 02      if(0.eq.1) then
    430       if(1.eq.1) then
     413      if(1==1) then
     414      if(iflag_pbl==8.or.iflag_pbl==10) then
     415
    431416      do k=2,klev
    432417         do ig=1,ngrid
    … …  
    449434         enddo
    450435      enddo
     436
     437      else
     438
     439      do k=2,klev
     440         do ig=1,ngrid
     441            if (teta(ig,2).gt.teta(ig,1)) then
     442               qmin=ustar(ig)*(max(1.-zlev(ig,k)/pblhmin(ig),0.))**2
     443               kmin=kap*zlev(ig,k)*qmin
     444            else
     445               kmin=-1. ! kmin n'est utilise que pour les SL stables.
     446            endif
     447            if (kn(ig,k).lt.kmin.or.km(ig,k).lt.kmin) then
     448c               print*,'Seuil min Km K=',k,kmin,km(ig,k),kn(ig,k)
     449c     s           ,sqrt(q2(ig,k)),pblhmin(ig),qmin/sm(ig,k)
     450               kn(ig,k)=kmin
     451               km(ig,k)=kmin
     452               kq(ig,k)=kmin
     453c   la longueur de melange est suposee etre l= kap z
     454c   K=l q Sm d'ou q2=(K/l Sm)**2
     455               sm(ig,k)=1.
     456               alpha(ig,k)=1.
     457               q2(ig,k)=min((qmin/sm(ig,k))**2,10.)
     458               zq=sqrt(q2(ig,k))
     459               km(ig,k)=l(ig,k)*zq*sm(ig,k)
     460               kn(ig,k)=km(ig,k)*alpha(ig,k)
     461               kq(ig,k)=l(ig,k)*zq*0.2
     462            endif
     463         enddo
     464      enddo
     465      endif
     466
    451467      endif
    452468
Note: See TracChangeset for help on using the changeset viewer.