Ignore:
Timestamp:
Apr 2, 2013, 11:48:33 AM (12 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

File:
1 edited

Legend:

Unmodified
Added
Removed

  • 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.