Ignore:
Timestamp:
Mar 2, 2012, 12:04:08 PM (13 years ago)
Author:
acolaitis
Message:

Yamada4: update to latest terrestrial version, adapation to mars, bug fix and cosmetic changes

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/libf/phymars/yamada4.F

    r544 r554  
     1!************************************************************
     2!************************************************************
     3!
     4!      YAMADA4 EARTH =>>> MARS VERSION
     5!      Modifications by: A.C. 02-03-2012 (marked by 'MARS')
     6!      Original version given by F.H. 01-03-2012
     7!
     8!************************************************************
     9!************************************************************
    110      SUBROUTINE yamada4(ngrid,nlay,dt,g,rconst,plev,temp
    2      s   ,zlev,zlay,u,v,teta,cd,q2,km,kn,kq,ustar
     11     s   ,zlev,zlay,u,v,phc,pq,cd,q2,km,kn,kq,ustar
    312     s   ,iflag_pbl)
    413      IMPLICIT NONE
    5 c.......................................................................
     14!.......................................................................
     15! MARS
    616#include "dimensions.h"
    717#include "dimphys.h"
    818#include "tracer.h"
    919#include "callkeys.h"
    10 c.......................................................................
    11 c
    12 c
    13 c dt : pas de temps
    14 c g  : g
    15 c zlev : altitude a chaque niveau (interface inferieure de la couche
    16 c        de meme indice)
    17 c zlay : altitude au centre de chaque couche
    18 c u,v : vitesse au centre de chaque couche
    19 c       (en entree : la valeur au debut du pas de temps)
    20 c teta : temperature potentielle au centre de chaque couche
    21 c        (en entree : la valeur au debut du pas de temps)
    22 c cd : cdrag
    23 c      (en entree : la valeur au debut du pas de temps)
    24 c q2 : $q^2$ au bas de chaque couche
    25 c      (en entree : la valeur au debut du pas de temps)
    26 c      (en sortie : la valeur a la fin du pas de temps)
    27 c km : diffusivite turbulente de quantite de mouvement (au bas de chaque
    28 c      couche)
    29 c      (en sortie : la valeur a la fin du pas de temps)
    30 c kn : diffusivite turbulente des scalaires (au bas de chaque couche)
    31 c      (en sortie : la valeur a la fin du pas de temps)
    32 c
    33 c  iflag_pbl doit valoir entre 6 et 9
    34 c      l=6, on prend  systematiquement une longueur d equilibre
    35 c    iflag_pbl=6 : MY 2.0
    36 c    iflag_pbl=7 : MY 2.0.Fournier
    37 c    iflag_pbl=8 : MY 2.5
    38 c    iflag_pbl>=9 : MY 2.5 avec diffusion verticale
    39 
    40 c.......................................................................
    41       REAL dt,g,rconst
    42       real plev(ngrid,nlay+1),temp(ngrid,nlay)
    43       real ustar(ngrid)
    44       real kmin,qmin,pblhmin(ngrid),coriol(ngrid)
    45       REAL zlev(ngrid,nlay+1)
    46       REAL zlay(ngrid,nlay)
    47       REAL u(ngrid,nlay)
    48       REAL v(ngrid,nlay)
    49       REAL teta(ngrid,nlay)
    50       REAL cd(ngrid)
    51       REAL q2(ngrid,nlay+1),qpre
     20!.......................................................................
     21!
     22! dt : pas de temps
     23! g  : g
     24! zlev : altitude a chaque niveau (interface inferieure de la couche
     25!        de meme indice)
     26! zlay : altitude au centre de chaque couche
     27! u,v : vitesse au centre de chaque couche
     28!       (en entree : la valeur au debut du pas de temps)
     29! phc : temperature potentielle au centre de chaque couche
     30!        (en entree : la valeur au debut du pas de temps)
     31! cd : cdrag
     32!      (en entree : la valeur au debut du pas de temps)
     33! q2 : $q^2$ au bas de chaque couche
     34!      (en entree : la valeur au debut du pas de temps)
     35!      (en sortie : la valeur a la fin du pas de temps)
     36! km : diffusivite turbulente de quantite de mouvement (au bas de chaque
     37!      couche)
     38!      (en sortie : la valeur a la fin du pas de temps)
     39! kn : diffusivite turbulente des scalaires (au bas de chaque couche)
     40!      (en sortie : la valeur a la fin du pas de temps)
     41!
     42!  iflag_pbl doit valoir entre 6 et 9
     43!      l=6, on prend  systematiquement une longueur d equilibre
     44!    iflag_pbl=6 : MY 2.0
     45!    iflag_pbl=7 : MY 2.0.Fournier
     46!    iflag_pbl=8 : MY 2.5
     47!    iflag_pbl>=9 : MY 2.5 avec diffusion verticale
     48!
     49!.......................................................................
     50
     51      REAL,    INTENT(IN)    :: dt,g,rconst
     52      REAL,    INTENT(IN)    :: u(ngrid,nlay)
     53      REAL,    INTENT(IN)    :: v(ngrid,nlay)
     54      REAL,    INTENT(IN)    :: phc(ngrid,nlay)
     55      REAL,    INTENT(IN)    :: cd(ngrid)
     56      REAL,    INTENT(IN)    :: temp(ngrid,nlay)
     57      REAL,    INTENT(IN)    :: plev(ngrid,nlay+1)
     58      REAL,    INTENT(IN)    :: ustar(ngrid)
     59      REAL,    INTENT(IN)    :: zlev(ngrid,nlay+1)
     60      REAL,    INTENT(IN)    :: zlay(ngrid,nlay)
     61      INTEGER, INTENT(IN)    :: iflag_pbl,ngrid
     62      INTEGER, INTENT(IN)    :: nlay
     63      REAL,    INTENT(INOUT) :: q2(ngrid,nlay+1)
     64      REAL,    INTENT(OUT)   :: km(ngrid,nlay+1)
     65      REAL,    INTENT(OUT)   :: kn(ngrid,nlay+1)
     66      REAL,    INTENT(OUT)   :: kq(ngrid,nlay+1)
     67
     68      REAL kmin,qmin,pblhmin(ngrid),coriol(ngrid)
    5269      REAL unsdz(ngrid,nlay)
    5370      REAL unsdzdec(ngrid,nlay+1)
    54 
    55       REAL km(ngrid,nlay+1)
    5671      REAL kmpre(ngrid,nlay+1),tmp2
    5772      REAL mpre(ngrid,nlay+1)
    58       REAL kn(ngrid,nlay+1)
    59       REAL kq(ngrid,nlay+1)
    60       real ff(ngrid,nlay+1),delta(ngrid,nlay+1)
    61       real aa(ngrid,nlay+1),aa0,aa1
    62       integer iflag_pbl,ngrid,nlay,nlev
    63 
    64 
    65       logical first
    66       integer ipas
    67       save first,ipas
    68 cFH/IM     data first,ipas/.true.,0/
    69       data first,ipas/.false.,0/
    70 
    71       integer ig,k
    72 
    73 
    74       real ri,zrif,zalpha,zsm,zsn
    75       real rif(ngrid,nlay+1),sm(ngrid,nlay+1),alpha(ngrid,nlay)
    76 
    77       real m2(ngrid,nlay+1),dz(ngrid,nlay+1),zq,n2(ngrid,nlay+1)
    78       real dtetadz(ngrid,nlay+1)
    79       real m2cstat,mcstat,kmcstat
    80       real l(ngrid,nlay+1)
    81       real,allocatable,save :: l0(:)
    82       real sq(ngrid),sqz(ngrid),zz(ngrid,nlay+1)
    83       integer iter
    84 
    85       real ric,rifc,b1,kap
    86       save ric,rifc,b1,kap
    87       data ric,rifc,b1,kap/0.195,0.191,16.6,0.4/
    88       real frif,falpha,fsm
    89       real fl,zzz,zl0,zq2,zn2
    90 
    91 ! MARS
    92       REAL,SAVE :: q2min,q2max,knmin,kmmin
    93       DATA q2min,q2max,knmin,kmmin/1.E-10,1.E+2,1.E-5,1.E-5/
    94 
    95 ! mars q2min 0.001 mais trop fort !
    96 !      PARAMETER (gnmin=-10.E+0)
    97 !      PARAMETER (gnmax=0.0233E+0)
    98 !      PARAMETER (a1=0.92E+0)
    99 !      PARAMETER (a2=0.74E+0)
    100 !      PARAMETER (b2=10.1E+0)
    101 !      PARAMETER (c1=0.08E+0)
    102 
    103       real rino(ngrid,nlay+1),smyam(ngrid,nlay),styam(ngrid,nlay)
     73      REAL ff(ngrid,nlay+1),delta(ngrid,nlay+1)
     74      REAL aa(ngrid,nlay+1),aa0,aa1,qpre
     75
     76      LOGICAL first
     77      INTEGER ipas,nlev
     78      SAVE first,ipas
     79!FH/IM     DATA first,ipas/.true.,0/
     80      DATA first,ipas/.false.,0/
     81      INTEGER ig,k
     82
     83
     84      REAL ri,zrif,zalpha,zsm,zsn
     85      REAL rif(ngrid,nlay+1),sm(ngrid,nlay+1),alpha(ngrid,nlay)
     86
     87      REAL m2(ngrid,nlay+1),dz(ngrid,nlay+1),zq,n2(ngrid,nlay+1)
     88      REAL dtetadz(ngrid,nlay+1)
     89      REAL m2cstat,mcstat,kmcstat
     90      REAL l(ngrid,nlay+1)
     91      REAL,allocatable,SAVE :: l0(:)
     92      REAL sq(ngrid),sqz(ngrid),zz(ngrid,nlay+1)
     93      INTEGER iter
     94
     95      REAL ric,rifc,b1,kap
     96      SAVE ric,rifc,b1,kap
     97      DATA ric,rifc,b1,kap/0.195,0.191,16.6,0.4/
     98      REAL frif,falpha,fsm
     99      REAL fl,zzz,zl0,zq2,zn2
     100
     101      REAL rino(ngrid,nlay+1),smyam(ngrid,nlay),styam(ngrid,nlay)
    104102     s  ,lyam(ngrid,nlay),knyam(ngrid,nlay)
    105103     s  ,w2yam(ngrid,nlay),t2yam(ngrid,nlay)
    106       logical,save :: firstcall=.true.
     104      LOGICAL,SAVE :: firstcall=.true.
    107105      frif(ri)=0.6588*(ri+0.1776-sqrt(ri*ri-0.3221*ri+0.03156))
    108106      falpha(ri)=1.318*(0.2231-ri)/(0.2341-ri)
    … …  
    113111
    114112
     113! MARS
     114      REAL,SAVE :: q2min,q2max,knmin,kmmin
     115      DATA q2min,q2max,knmin,kmmin/1.E-10,1.E+2,1.E-5,1.E-5/
     116      INTEGER ico2,iq
     117      SAVE ico2
     118      REAL m_co2, m_noco2, A , B
     119      SAVE A, B
     120      REAL teta(ngrid,nlay)
     121      REAL pq(ngrid,nlay,nqmx)
     122
    115123      nlev=nlay+1
     124
     125c.......................................................................
     126c  Initialization
     127c.......................................................................
     128
     129      if(firstcall) then
     130        ico2=0
     131        if (tracer) then
     132!     Prepare Special treatment if one of the tracers is CO2 gas
     133           do iq=1,nqmx
     134             if (noms(iq).eq."co2") then
     135                ico2=iq
     136                m_co2 = 44.01E-3  ! CO2 molecular mass (kg/mol)
     137                m_noco2 = 33.37E-3  ! Non condensible mol mass (kg/mol)
     138!               Compute A and B coefficient use to compute
     139!               mean molecular mass Mair defined by
     140!               1/Mair = q(ico2)/m_co2 + (1-q(ico2))/m_noco2
     141!               1/Mair = A*q(ico2) + B
     142                A =(1/m_co2 - 1/m_noco2)
     143                B=1/m_noco2
     144             end if
     145           enddo
     146        endif
     147      allocate(l0(ngrid))
     148      firstcall=.false.
     149      endif !of if firstcall
     150
     151c.......................................................................
     152c  Special treatment for co2
     153c.......................................................................
     154
     155      if (ico2.ne.0) then
     156!     Special case if one of the tracers is CO2 gas
     157         DO k=1,nlay
     158           DO ig=1,ngrid
     159            teta(ig,k) = phc(ig,k)*(A*pq(ig,k,ico2)+B)
     160           ENDDO
     161         ENDDO
     162       else
     163          teta(:,:)=phc(:,:)
     164       end if
    116165     
    117       if (firstcall) then
    118         allocate(l0(ngrid))
    119         firstcall=.false.
    120       endif
    121 
    122 
    123166      if (.not.(iflag_pbl.ge.6.and.iflag_pbl.le.10)) then
    124167           stop'probleme de coherence dans appel a MY'
    … …  
    126169
    127170      ipas=ipas+1
    128       if (0.eq.1.and.first) then
    129       do ig=1,1000
    130          ri=(ig-800.)/500.
    131          if (ri.lt.ric) then
    132             zrif=frif(ri)
    133          else
    134             zrif=rifc
    135          endif
    136          if(zrif.lt.0.16) then
    137             zalpha=falpha(zrif)
    138             zsm=fsm(zrif)
    139          else
    140             zalpha=1.12
    141             zsm=0.085
    142          endif
    143 c     print*,ri,rif,zalpha,zsm
    144       enddo
    145       endif
    146 
    147 c.......................................................................
    148 c  les increments verticaux
    149 c.......................................................................
    150 c
    151 c!!!!! allerte !!!!!c
    152 c!!!!! zlev n'est pas declare a nlay !!!!!c
    153 c!!!!! ---->
    154                                                       DO ig=1,ngrid
    155             zlev(ig,nlev)=zlay(ig,nlay)
    156      &             +( zlay(ig,nlay) - zlev(ig,nlev-1) )
    157                                                       ENDDO
    158 c!!!!! <----
    159 c!!!!! allerte !!!!!c
    160 c
     171! MARS
     172!      if (0.eq.1.and.first) then
     173!      do ig=1,1000
     174!         ri=(ig-800.)/500.
     175!         if (ri.lt.ric) then
     176!            zrif=frif(ri)
     177!         else
     178!            zrif=rifc
     179!         endif
     180!         if(zrif.lt.0.16) then
     181!            zalpha=falpha(zrif)
     182!            zsm=fsm(zrif)
     183!         else
     184!            zalpha=1.12
     185!            zsm=0.085
     186!         endif
     187!     print*,ri,rif,zalpha,zsm
     188!      enddo
     189!      endif
     190
     191!.......................................................................
     192!  les increments verticaux
     193!.......................................................................
     194!
     195!!!!!! allerte !!!!!
     196!!!!!! zlev n'est pas declare a nlev !!!!!
     197!!!!!! ---->
     198! MARS
     199!
     200!                                                      DO ig=1,ngrid
     201!            zlev(ig,nlev)=zlay(ig,nlay)
     202!     &             +( zlay(ig,nlay) - zlev(ig,nlev-1) )
     203!                                                      ENDDO
     204!!!!! <----
     205!!!!! allerte !!!!!
     206
    161207      DO k=1,nlay
    162208                                                      DO ig=1,ngrid
    … …  
    175221      unsdzdec(ig,nlay+1)=1.E+0/(zlev(ig,nlay+1)-zlay(ig,nlay))
    176222                                                     ENDDO
    177 c
    178 c.......................................................................
     223!
     224!.......................................................................
    179225
    180226      do k=2,nlay
    … …  
    185231         dtetadz(ig,k)=(teta(ig,k)-teta(ig,k-1))/dz(ig,k)
    186232         n2(ig,k)=g*2.*dtetadz(ig,k)/(teta(ig,k-1)+teta(ig,k))
    187 c        n2(ig,k)=0.
     233!        n2(ig,k)=0.
    188234         ri=n2(ig,k)/max(m2(ig,k),1.e-10)
    189235         if (ri.lt.ric) then
    … …  
    200246         endif
    201247         zz(ig,k)=b1*m2(ig,k)*(1.-rif(ig,k))*sm(ig,k)
    202 c     print*,'RIF L=',k,rif(ig,k),ri*alpha(ig,k)
    203 
    204 
    205                                                           enddo
    206       enddo
    207 
    208 
    209 c====================================================================
    210 c   Au premier appel, on determine l et q2 de facon iterative.
    211 c iterration pour determiner la longueur de melange
     248!     print*,'RIF L=',k,rif(ig,k),ri*alpha(ig,k)
     249
     250
     251                                                          enddo
     252      enddo
     253
     254
     255!====================================================================
     256!   Au premier appel, on determine l et q2 de facon iterative.
     257! iterration pour determiner la longueur de melange
    212258
    213259
    214260      if (first.or.iflag_pbl.eq.6) then
    215261                                                          do ig=1,ngrid
    216 !      l0(ig)=10. !EARTH
     262! MARS
     263!      l0(ig)=10.
    217264      l0(ig)=160.
    218265                                                          enddo
    … …  
    239286                                                          do ig=1,ngrid
    240287         l0(ig)=0.2*sqz(ig)/sq(ig)
    241 c        l0(ig)=30.
    242                                                           enddo
    243 c      print*,'ITER=',iter,'  L0=',l0
    244 
    245       enddo
    246 
    247 c     print*,'Fin de l initialisation de q2 et l0'
     288!        l0(ig)=30.
     289                                                          enddo
     290!      print*,'ITER=',iter,'  L0=',l0
     291
     292      enddo
     293
     294!     print*,'Fin de l initialisation de q2 et l0'
    248295
    249296      endif ! first
    250297
    251 c====================================================================
    252 c  Calcul de la longueur de melange.
    253 c====================================================================
    254 
    255 c   Mise a jour de l0
     298!====================================================================
     299!  Calcul de la longueur de melange.
     300!====================================================================
     301
     302!   Mise a jour de l0
    256303                                                          do ig=1,ngrid
    257304      sq(ig)=1.e-10
    … …  
    267314                                                          do ig=1,ngrid
    268315      l0(ig)=0.2*sqz(ig)/sq(ig)
    269 c        l0(ig)=30.
    270                                                           enddo
    271 c      print*,'ITER=',iter,'  L0=',l0
    272 c   calcul de l(z)
     316!        l0(ig)=30.
     317                                                          enddo
     318!      print*,'ITER=',iter,'  L0=',l0
     319!   calcul de l(z)
    273320      do k=2,nlay
    274321                                                          do ig=1,ngrid
    … …  
    280327      enddo
    281328
    282 c====================================================================
    283 c   Yamada 2.0
    284 c====================================================================
     329!====================================================================
     330!   Yamada 2.0
     331!====================================================================
    285332      if (iflag_pbl.eq.6) then
    286333
    … …  
    293340
    294341      else if (iflag_pbl.eq.7) then
    295 c====================================================================
    296 c   Yamada 2.Fournier
    297 c====================================================================
    298 
    299 c  Calcul de l,  km, au pas precedent
     342!====================================================================
     343!   Yamada 2.Fournier
     344!====================================================================
     345
     346!  Calcul de l,  km, au pas precedent
    300347      do k=2,nlay
    301348                                                          do ig=1,ngrid
    … …  
    313360        mcstat=sqrt(m2cstat)
    314361
    315 c        print*,'M2 L=',k,mpre(ig,k),mcstat
    316 c
    317 c  -----{puis on ecrit la valeur de q qui annule l equation de m
    318 c        supposee en q3}
    319 c
     362!        print*,'M2 L=',k,mpre(ig,k),mcstat
     363!
     364!  -----{puis on ecrit la valeur de q qui annule l'equation de m
     365!        supposee en q3}
     366!
    320367        IF (k.eq.2) THEN
    321368          kmcstat=1.E+0 / mcstat
    … …  
    336383     &      /( unsdz(ig,k)+unsdz(ig,k-1) )
    337384        ENDIF
    338 c       print*,'T2 L=',k,tmp2
     385!       print*,'T2 L=',k,tmp2
    339386        tmp2=kmcstat
    340387     &      /( sm(ig,k)/q2(ig,k) )
    341388     &      /l(ig,k)
    342 !       q2(ig,k)=max(tmp2,1.E-10)**(2./3.)
    343 ! MARS
    344         q2(ig,k)=max(tmp2,q2min)**(2./3.)
    345 c       print*,'Q2 L=',k,q2(ig,k)
    346 c
     389
     390! MARS
     391!        q2(ig,k)=max(tmp2,1.e-12)**(2./3.)
     392        q2(ig,k)=max(q2min,max(tmp2,1.e-12)**(2./3.))
     393
     394!       print*,'Q2 L=',k,q2(ig,k)
     395!
    347396                                                          enddo
    348397      enddo
    349398
    350399      else if (iflag_pbl.ge.8) then
    351 c====================================================================
    352 c   Yamada 2.5 a la Didi
    353 c====================================================================
    354 
    355 
    356 c  Calcul de l,  km, au pas precedent
    357       do k=2,nlay
    358                                                           do ig=1,ngrid
    359 c        print*,'SMML=',sm(ig,k),l(ig,k)
     400!====================================================================
     401!   Yamada 2.5 a la Didi
     402!====================================================================
     403
     404
     405!  Calcul de l,  km, au pas precedent
     406      do k=2,nlay
     407                                                          do ig=1,ngrid
     408!        print*,'SMML=',sm(ig,k),l(ig,k)
    360409         delta(ig,k)=q2(ig,k)/(l(ig,k)**2*sm(ig,k))
    361410         if (delta(ig,k).lt.1.e-20) then
    362 c     print*,'ATTENTION   L=',k,'   Delta=',delta(ig,k)
     411!     print*,'ATTENTION   L=',k,'   Delta=',delta(ig,k)
    363412            delta(ig,k)=1.e-20
    364413         endif
    … …  
    368417         aa1=
    369418     s   (m2(ig,k)*(1.-rif(ig,k))-delta(ig,k)/b1)
    370 c abder      print*,'AA L=',k,aa0,aa1,aa1/max(m2(ig,k),1.e-20)
     419! abder      print*,'AA L=',k,aa0,aa1,aa1/max(m2(ig,k),1.e-20)
    371420         aa(ig,k)=aa1*dt/(delta(ig,k)*l(ig,k))
    372 c     print*,'0L=',k,l(ig,k),delta(ig,k),km(ig,k)
     421!     print*,'0L=',k,l(ig,k),delta(ig,k),km(ig,k)
    373422         qpre=sqrt(q2(ig,k))
    374423         if (iflag_pbl.eq.8 ) then
    … …  
    385434            endif
    386435         endif
     436
     437! MARS
     438         q2(ig,k)=min(max(q2(ig,k),q2min),q2max)
    387439!         q2(ig,k)=min(max(q2(ig,k),1.e-10),1.e4)
    388 ! MARS
    389          q2(ig,k)=min(max(q2(ig,k),q2min),q2max)
    390 c     print*,'Q2 L=',k,q2(ig,k),qpre*qpre
    391                                                           enddo
    392       enddo
    393 
    394 !MARS
    395 !      q2(:,nlay+1)=q2min
    396 ! OU ALORS MARS
     440
     441!     print*,'Q2 L=',k,q2(ig,k),qpre*qpre
     442                                                          enddo
     443      enddo
     444
     445! MARS
    397446      q2(:,nlay+1)=q2(:,nlay)
     447
    398448      endif ! Fin du cas 8
    399449
    400 c     print*,'OK8'
    401 
    402 c====================================================================
    403 c   Calcul des coefficients de m�ange
    404 c====================================================================
    405       do k=2,nlay
    406 c     print*,'k=',k
    407                                                           do ig=1,ngrid
    408 cabde      print*,'KML=',l(ig,k),q2(ig,k),sm(ig,k)
     450!     print*,'OK8'
     451
     452!====================================================================
     453!   Calcul des coefficients de melange
     454!====================================================================
     455      do k=2,nlay
     456!     print*,'k=',k
     457                                                          do ig=1,ngrid
     458!abde      print*,'KML=',l(ig,k),q2(ig,k),sm(ig,k)
    409459         zq=sqrt(q2(ig,k))
    410460         km(ig,k)=l(ig,k)*zq*sm(ig,k)
    411461         kn(ig,k)=km(ig,k)*alpha(ig,k)
    412 ! MARS
    413 !         km(ig,k)=max(l(ig,k)*zq*sm(ig,k),kmmin)
    414 !         kn(ig,k)=max(km(ig,k)*alpha(ig,k),knmin)
    415462         kq(ig,k)=l(ig,k)*zq*0.2
    416 c       print*,'KML=',km(ig,k),kn(ig,k)
    417                                                           enddo
    418       enddo
    419 ! MARS
    420 !      km(:,nlay+1)=kmmin
    421 !      kn(:,nlay+1)=knmin
    422 !      kq(:,nlay+1)=knmin
    423 ! OU ALORS MARS
     463!     print*,'KML=',km(ig,k),kn(ig,k)
     464                                                          enddo
     465      enddo
     466
     467! MARS
    424468      km(:,nlay+1)=km(:,nlay)
    425469      kn(:,nlay+1)=kn(:,nlay)
    … …  
    433477      endif
    434478
    435 c   Traitement des cas noctrunes avec l'introduction d'une longueur
    436 c   minimale.
    437 
    438 c====================================================================
    439 c   Traitement particulier pour les cas tres stables.
    440 c   D apres Holtslag Boville.
    441 
    442                                                           do ig=1,ngrid
    443       coriol(ig)=1.e-4
    444       pblhmin(ig)=0.07*ustar(ig)/max(abs(coriol(ig)),2.546e-5)
    445                                                           enddo
    446 
     479!   Traitement des cas noctrunes avec l'introduction d'une longueur
     480!   minilale.
     481!
     482!====================================================================
     483!   Traitement particulier pour les cas tres stables.
     484!   D'apres Holtslag Boville.
     485
     486! MARS : deactivating pblhmin following F.H. concerns
     487
     488!                                                          do ig=1,ngrid
     489!      coriol(ig)=1.e-4
     490!      pblhmin(ig)=0.07*ustar(ig)/max(abs(coriol(ig)),2.546e-5)
     491!                                                          enddo
     492!
    447493!      print*,'pblhmin ',pblhmin
    448 CTest a remettre 21 11 02
    449 c test abd 13 05 02      if(0.eq.1) then
    450       if(1.eq.1) then
    451       do k=2,nlay
    452          do ig=1,ngrid
    453             if (teta(ig,2).gt.teta(ig,1)) then
    454                qmin=ustar(ig)*(max(1.-zlev(ig,k)/pblhmin(ig),0.))**2
    455                kmin=kap*zlev(ig,k)*qmin
    456             else
    457                kmin=-1. ! kmin n'est utilise que pour les SL stables.
    458             endif
    459             if (kn(ig,k).lt.kmin.or.km(ig,k).lt.kmin) then
    460 c               print*,'Seuil min Km K=',k,kmin,km(ig,k),kn(ig,k)
    461 c     s           ,sqrt(q2(ig,k)),pblhmin(ig),qmin/sm(ig,k)
    462                kn(ig,k)=kmin
    463                km(ig,k)=kmin
    464                kq(ig,k)=kmin
    465 c   la longueur de melange est suposee etre l= kap z
    466 c   K=l q Sm d ou q2=(K/l Sm)**2
    467                q2(ig,k)=(qmin/sm(ig,k))**2
    468             endif
    469          enddo
    470       enddo
    471       endif
    472 
    473 c   Diagnostique pour stokage
     494!CTest a remettre 21 11 02
     495! test abd 13 05 02      if(0.eq.1) then
     496!      if(0.eq.1) then
     497!      do k=2,nlay
     498!         do ig=1,ngrid
     499!            if (teta(ig,2).gt.teta(ig,1)) then
     500!               qmin=ustar(ig)*(max(1.-zlev(ig,k)/pblhmin(ig),0.))**2
     501!               kmin=kap*zlev(ig,k)*qmin
     502!            else
     503!               kmin=-1. ! kmin n'est utilise que pour les SL stables.
     504!            endif
     505!            if (kn(ig,k).lt.kmin.or.km(ig,k).lt.kmin) then
     506!               print*,'Seuil min Km K=',k,kmin,km(ig,k),kn(ig,k)
     507!     s           ,sqrt(q2(ig,k)),pblhmin(ig),qmin/sm(ig,k)
     508!               kn(ig,k)=kmin
     509!               km(ig,k)=kmin
     510!               kq(ig,k)=kmin
     511!   la longueur de melange est suposee etre l= kap z
     512!   K=l q Sm d'ou q2=(K/l Sm)**2
     513!               q2(ig,k)=(qmin/sm(ig,k))**2
     514!            endif
     515!         enddo
     516!      enddo
     517!      endif
     518
     519!   Diagnostique pour stokage
    474520
    475521      if(1.eq.0)then
    … …  
    487533      knyam(1:ngrid,2:nlay)=kn(1:ngrid,2:nlay)
    488534
    489 c   Estimations de w'2 et T'2 d apres Abdela et McFarlane
     535!   Estimations de w'2 et T'2 d'apres Abdela et McFarlane
    490536
    491537      w2yam(1:ngrid,2:nlay)=q2(1:ngrid,2:nlay)*0.24
    … …  
    498544      endif
    499545
    500 c     print*,'OKFIN'
     546!     print*,'OKFIN'
    501547      first=.false.
    502548      return
    … …  
    505551     &  kmy,q2)
    506552      IMPLICIT NONE
    507 c.......................................................................
     553!.......................................................................
     554! MARS
    508555#include "dimensions.h"
    509556#include "dimphys.h"
    510557#include "tracer.h"
    511558#include "callkeys.h"
    512 c.......................................................................
    513 c
    514 c dt : pas de temps
    515 
    516       real plev(ngrid,nlay+1)
    517       real temp(ngrid,nlay)
    518       real timestep
    519       real gravity,rconst
    520       real kstar(ngrid,nlay+1),zz
    521       real kmy(ngrid,nlay+1)
    522       real q2(ngrid,nlay+1)
    523       real deltap(ngrid,nlay+1)
    524       real denom(ngrid,nlay+1),alpha(ngrid,nlay+1),beta(ngrid,nlay+1)
    525       integer ngrid,nlay
    526 
    527       integer i,k
     559!.......................................................................
     560!
     561! dt : pas de temps
     562!
     563      REAL plev(ngrid,nlay+1)
     564      REAL temp(ngrid,nlay)
     565      REAL timestep
     566      REAL gravity,rconst
     567      REAL kstar(ngrid,nlay+1),zz
     568      REAL kmy(ngrid,nlay+1)
     569      REAL q2(ngrid,nlay+1)
     570      REAL deltap(ngrid,nlay+1)
     571      REAL denom(ngrid,nlay+1),alpha(ngrid,nlay+1),beta(ngrid,nlay+1)
     572      INTEGER ngrid,nlay
     573
     574      INTEGER i,k
    528575
    529576!       print*,'RD=',rconst
    530577      do k=1,nlay
    531578         do i=1,ngrid
    532 c test
     579! test
    533580!       print*,'i,k',i,k
    534581!       print*,'temp(i,k)=',temp(i,k)
    … …  
    557604            denom(i,k)=deltap(i,k)+(1.-beta(i,k+1))*
    558605     s      kstar(i,k)+kstar(i,k-1)
    559 c   correction d un bug 10 01 2001
     606!   correction d'un bug 10 01 2001
    560607            alpha(i,k)=(q2(i,k)*deltap(i,k)
    561608     s      +kstar(i,k)*alpha(i,k+1))/denom(i,k)
    … …  
    564611      enddo
    565612
    566 c  Si on recalcule q2(1)
     613!  Si on recalcule q2(1)
    567614      if(1.eq.0) then
    568615      do i=1,ngrid
    … …  
    572619      enddo
    573620      endif
    574 c   sinon, on peut sauter cette boucle...
     621!   sinon, on peut sauter cette boucle...
    575622
    576623      do k=2,nlay+1
    … …  
    585632     &   plev,temp,kmy,q2)
    586633      IMPLICIT NONE
    587 c.......................................................................
     634!.......................................................................
     635! MARS
    588636#include "dimensions.h"
    589637#include "dimphys.h"
    590638#include "tracer.h"
    591639#include "callkeys.h"
    592 c.......................................................................
    593 c
    594 c dt : pas de temps
    595 
    596       real plev(ngrid,nlay+1)
    597       real temp(ngrid,nlay)
    598       real timestep
    599       real gravity,rconst
    600       real kstar(ngrid,nlay+1),zz
    601       real kmy(ngrid,nlay+1)
    602       real q2(ngrid,nlay+1)
    603       real deltap(ngrid,nlay+1)
    604       real denom(ngrid,nlay+1),alpha(ngrid,nlay+1),beta(ngrid,nlay+1)
    605       integer ngrid,nlay
    606 
    607       integer i,k
     640!.......................................................................
     641!
     642! dt : pas de temps
     643
     644      REAL plev(ngrid,nlay+1)
     645      REAL temp(ngrid,nlay)
     646      REAL timestep
     647      REAL gravity,rconst
     648      REAL kstar(ngrid,nlay+1),zz
     649      REAL kmy(ngrid,nlay+1)
     650      REAL q2(ngrid,nlay+1)
     651      REAL deltap(ngrid,nlay+1)
     652      REAL denom(ngrid,nlay+1),alpha(ngrid,nlay+1),beta(ngrid,nlay+1)
     653      INTEGER ngrid,nlay
     654
     655      INTEGER i,k
    608656
    609657      do k=1,nlay
Note: See TracChangeset for help on using the changeset viewer.