Changeset 1127


Ignore:
Timestamp:
Mar 19, 2009, 11:38:04 AM (15 years ago)
Author:
idelkadi
Message:

Corrections sur les wakes et la convection pour surmonter le probleme de l'eau negative

Location:
LMDZ4/branches/LMDZ4-dev/libf/phylmd
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • LMDZ4/branches/LMDZ4-dev/libf/phylmd/cv3_cine.F

    r879 r1127  
    3333      integer ifst(nloc),isublcl(nloc)
    3434      logical lswitch(nloc),lswitch1(nloc),lswitch2(nloc)
     35      logical exist_lfc(nloc)
     36      real plfc(nloc)
    3537      real dpmax
    3638      real deltap,dcin
    3739      real buoylcl(nloc),tvplcl(nloc),tvlcl(nloc)
    38       real plfc(nloc),p0(nloc)
     40      real p0(nloc)
    3941      real buoyz(nloc), buoy(nloc,nd)
    4042c
    … …  
    5052c      Recompute buoyancies
    5153c--------------------------------------------------------------
    52       DO k = 1,nl
     54      DO k = 1,nd
    5355        DO il = 1,ncum
     56!      print*,'tvp tv=',tvp(il,k),tv(il,k)
    5457          buoy(il,k) = tvp(il,k) - tv(il,k)
    5558        ENDDO
    5659      ENDDO
    57 c
    58 c---------------------------------------------------------------
    59 c premiere couche contenant un  niveau de flotabilite positive
    60 c et premiere couche contenant un  niveau de flotabilite negative
    61 c  au dessus du niveau de condensation
    62 c---------------------------------------------------------------
    63       do il = 1,ncum
    64         itop(il) =nl-1
    65         ineg(il) = nl-1
    66       enddo
    67       do 100 k=nl,1,-1
    68        do 110 il=1,ncum
    69         if (k .ge. icb(il)) then
    70          if (buoy(il,k) .gt. 0.) then
    71           itop(il)=k
    72          else
    73           ineg(il)=k
    74          endif
    75         endif
    76 110    continue
    77 100   continue
    78 c      print *,' itop, ineg, icb ',itop(1),ineg(1), icb(1)
    79 c
    8060c---------------------------------------------------------------
    8161c
    … …  
    10989c
    11090c---------------------------------------------------------------
     91c premiere couche contenant un  niveau de flotabilite positive
     92c et premiere couche contenant un  niveau de flotabilite negative
     93c  au dessus du niveau de condensation
     94c---------------------------------------------------------------
     95      do il = 1,ncum
     96        itop(il) =nl-1
     97        ineg(il) = nl-1
     98        exist_lfc(il) = .FALSE.
     99      enddo
     100      do 100 k=nl-1,1,-1
     101       do 110 il=1,ncum
     102        if (k .ge. ifst(il)) then
     103         if (buoy(il,k) .gt. 0.) then
     104          itop(il)=k
     105          exist_lfc(il) = .TRUE.
     106         else
     107          ineg(il)=k
     108         endif
     109        endif
     110110    continue
     111100   continue
     112c
     113c---------------------------------------------------------------
     114c When there is no positive buoyancy level, set Plfc, Cina and Cinb
     115c to arbitrary extreme values.
     116c---------------------------------------------------------------
     117      DO il = 1,ncum
     118       IF (.NOT.exist_lfc(il)) THEN
     119         Plfc(il) = 1.111
     120         Cinb(il) = -1111.
     121         Cina(il) = -1112.
     122       ENDIF
     123      ENDDO
     124c
     125c
     126c---------------------------------------------------------------
    111127c -- Two cases : BUOYlcl >= 0 and BUOYlcl < 0.
    112128c---------------------------------------------------------------
    … …  
    118134      DPMAX = 50.
    119135      DO il = 1,ncum
    120         lswitch1(il)=BUOYlcl(il) .GE. 0.
     136        lswitch1(il)=BUOYlcl(il) .GE. 0. .AND. exist_lfc(il)
    121137        lswitch(il) = lswitch1(il)
    122138      ENDDO
    … …  
    233249C
    234250      DO il = 1,ncum
    235         lswitch1(il)=BUOYlcl(il) .LT. 0.
     251        lswitch1(il)=BUOYlcl(il) .LT. 0. .AND. exist_lfc(il)
    236252        lswitch(il) = lswitch1(il)
    237253      ENDDO
    … …  
    239255c 2.0.1 Premiere  couche ou la flotabilite est negative au dessus du sol
    240256c ----------------------------------------------------
    241 c    au cas ou il existe  sinon ilow=1 (nk apres)
     257c    au cas ou elle existe  sinon ilow=1 (nk apres)
    242258c      on suppose que la parcelle part de la premiere couche
    243259c
    … …  
    248264      ENDDO
    249265c
    250       do 200 i=nl,1,-1
     266      do 200 k=nl,1,-1
    251267        DO il = 1,ncum
    252268        IF (lswitch(il) .AND. k .LE.icb(il)-1) THEN
    … …  
    292308        dcin = RD*(BUOYz(il)+BUOYlcl(il))*deltap/(P0(il)+Plcl(il))
    293309        CINB(il) = min(0.,dcin)
    294 cc        print *,'buoyz(il),buoylcl(il),deltap,p0(il),plcl(il),dcin ',
    295 cc     $           buoyz(il),buoylcl(il),deltap,p0(il),plcl(il),dcin
    296       ENDIF
    297       ENDDO
    298 c        print*, 'CINB ',CINB(1),'DCIN ',DCIN,I,BUOYz(1),BUOYlcl(1)
     310      ENDIF
     311      ENDDO
    299312c
    300313      DO il = 1,ncum
    … …  
    316329      ENDDO
    317330c
    318       IF (lswitch(1)) THEN
    319 c        print*,'ilow= ',ilow(1),'DCIN0 ',DCIN,P0(1),P(1,ilow(1))
    320 c        print*,'buoy',(BUOY(1,k),k=1,itop(1))
    321       ENDIF
    322331c
    323332C  Middle part of CINB : integral from P(ilow) to P(isublcl)
    … …  
    332341        ENDIF
    333342        ENDDO
    334 c      print*, 'CINB ', CINB(1), 'DCIN',DCIN,k,BUOY(1,k),BUOY(1,k+1)
    335343      ENDDO
    336344c
    … …  
    345353      ENDDO
    346354C
    347 c        print*, ' CINB ', CINB(1), 'Dcin ',dcin
    348355c
    349356cc      ENDIF
    … …  
    439446      ENDDO
    440447cc      ENDIF
    441 c       Print *,' Plcl,P(itop-1),P(itop),PLFC,BUOYlcl'
    442 c     $      ,Plcl(1),P(1,itop(1)-1),P(1,itop(1)),PLFC(1),BUOYlcl(1)
    443 C
    444 c        print*, 'CIN above', CINA(1), 'CIN below',CINB(1)
    445448c
    446449

  • LMDZ4/branches/LMDZ4-dev/libf/phylmd/cv3_routines.F

    r1044 r1127  
    11!
    2 ! $Header$
     2! $Header: /home/cvsroot/LMDZ4/libf/phylmd/cv3_routines.F,v 1.16 2008-11-06 16:29:35 lmdzadmin Exp $
    33!
    44c
    … …  
    120120
    121121c ori      do 110 k=1,nlp
    122       do 110 k=1,nl ! convect3
     122! abderr     do 110 k=1,nl ! convect3
     123       do 110 k=1,nlp
     124     
    123125        do 100 i=1,len
    124126cdebug          lv(i,k)= lv0-clmcpv*(t(i,k)-t0)
    … …  
    22562258      SUBROUTINE cv3_yield(nloc,ncum,nd,na,ntra
    22572259     :                    ,icb,inb,delt
    2258      :                    ,t,rr,t_wake,rr_wake,u,v,tra
     2260     :                    ,t,rr,t_wake,rr_wake,s_wake,u,v,tra
    22592261     :                    ,gz,p,ph,h,hp,lv,cpn,th,th_wake
    22602262     :                    ,ep,clw,m,tp,mp,rp,up,vp,trap
    … …  
    22832285      real t(nloc,nd), rr(nloc,nd), u(nloc,nd), v(nloc,nd)
    22842286      real t_wake(nloc,nd), rr_wake(nloc,nd)
     2287      real s_wake(nloc)
    22852288      real tra(nloc,nd,ntra), sig(nloc,nd)
    22862289      real gz(nloc,na), ph(nloc,nd+1), h(nloc,na), hp(nloc,na)
    … …  
    23272330      real esum(nloc), fsum(nloc), gsum(nloc), hsum(nloc)
    23282331      real th_wake(nloc,nd)
     2332      real alpha_qpos(nloc)
    23292333      real qcond(nloc,nd), nqcond(nloc,nd), wa(nloc,nd)  ! cld
    23302334      real siga(nloc,nd), sax(nloc,nd), mac(nloc,nd)      ! cld
    … …  
    29612965c   ***          integrated enthalpy and water tendencies         ***
    29622966c
     2967c Correction bug le 18-03-09
    29632968      do 503 il=1,ncum
    29642969      IF (iflag(il) .le. 1) THEN
    2965       ax=0.1*ment(il,inb(il),inb(il))*(hp(il,inb(il))-h(il,inb(il))
    2966      : +t(il,inb(il))*(cpv-cpd)
     2970       ax=0.1*ment(il,inb(il),inb(il))*(hp(il,inb(il))
     2971     : -h(il,inb(il))+t(il,inb(il))*(cpv-cpd)
    29672972     : *(rr(il,inb(il))-qent(il,inb(il),inb(il))))
    29682973     :  /(cpn(il,inb(il))*(ph(il,inb(il))-ph(il,inb(il)+1)))
    … …  
    30653070       enddo
    30663071      enddo
     3072
     3073c
     3074c   ***   Check that moisture stays positive. If not, scale tendencies
     3075c        in order to ensure moisture positivity
     3076      DO il = 1,ncum
     3077       IF (iflag(il) .le. 1) THEN
     3078        alpha_qpos(il) = max(1. , -delt*fr(il,1)/
     3079     :     (s_wake(il)*rr_wake(il,1)+(1.-s_wake(il))*rr(il,1)))
     3080       ENDIF
     3081      ENDDO
     3082      DO i = 2,nl
     3083       DO il = 1,ncum
     3084        IF (iflag(il) .le. 1) THEN
     3085        alpha_qpos(il) = max(alpha_qpos(il) , -delt*fr(il,i)/
     3086     :     (s_wake(il)*rr_wake(il,i)+(1.-s_wake(il))*rr(il,i)))
     3087        ENDIF
     3088       ENDDO
     3089      ENDDO
     3090      DO il = 1,ncum
     3091       IF (iflag(il) .le. 1 .and. alpha_qpos(il) .gt. 1.001) THEN
     3092        alpha_qpos(il) = alpha_qpos(il)*1.1
     3093       ENDIF
     3094      ENDDO
     3095      DO il = 1,ncum
     3096       IF (iflag(il) .le. 1) THEN
     3097        sigd(il) = sigd(il)/alpha_qpos(il)
     3098        precip(il) = precip(il)/alpha_qpos(il)
     3099       ENDIF
     3100      ENDDO
     3101      DO i = 1,nl
     3102       DO il = 1,ncum
     3103        IF (iflag(il) .le. 1) THEN
     3104         fr(il,i) = fr(il,i)/alpha_qpos(il)
     3105         ft(il,i) = ft(il,i)/alpha_qpos(il)
     3106         fqd(il,i) = fqd(il,i)/alpha_qpos(il)
     3107         ftd(il,i) = ftd(il,i)/alpha_qpos(il)
     3108         fu(il,i) = fu(il,i)/alpha_qpos(il)
     3109         fv(il,i) = fv(il,i)/alpha_qpos(il)
     3110         m(il,i) = m(il,i)/alpha_qpos(il)
     3111         mp(il,i) = mp(il,i)/alpha_qpos(il)
     3112         Vprecip(il,i) = Vprecip(il,i)/alpha_qpos(il)
     3113        ENDIF
     3114       ENDDO
     3115      ENDDO
     3116      DO i = 1,nl
     3117      DO j = 1,nl
     3118       DO il = 1,ncum
     3119        IF (iflag(il) .le. 1) THEN
     3120         ment(il,i,j) = ment(il,i,j)/alpha_qpos(il)
     3121        ENDIF
     3122       ENDDO
     3123      ENDDO
     3124      ENDDO
     3125      DO j = 1,ntra
     3126      DO i = 1,nl
     3127       DO il = 1,ncum
     3128        IF (iflag(il) .le. 1) THEN
     3129         ftra(il,i,j) = ftra(il,i,j)/alpha_qpos(il)
     3130        ENDIF
     3131       ENDDO
     3132      ENDDO
     3133      ENDDO
     3134
    30673135c
    30683136c   ***           reset counter and return           ***

  • LMDZ4/branches/LMDZ4-dev/libf/phylmd/cv3a_compress.F

    r972 r1127  
    33     :    ,plcl1,tnk1,qnk1,gznk1,hnk1,unk1,vnk1
    44     :    ,wghti1,pbase1,buoybase1
    5      :    ,t1,q1,qs1,t1_wake,q1_wake,qs1_wake,u1,v1,gz1,th1,th1_wake
     5     :    ,t1,q1,qs1,t1_wake,q1_wake,qs1_wake,s1_wake
     6     :    ,u1,v1,gz1,th1,th1_wake
    67     :    ,tra1
    78     :    ,h1     ,lv1     ,cpn1   ,p1,ph1,tv1    ,tp1,tvp1,clw1
    … …  
    1213     o    ,plcl,tnk,qnk,gznk,hnk,unk,vnk
    1314     o    ,wghti,pbase,buoybase
    14      o    ,t,q,qs,t_wake,q_wake,qs_wake,u,v,gz,th,th_wake
     15     o    ,t,q,qs,t_wake,q_wake,qs_wake,s_wake
     16     o    ,u,v,gz,th,th_wake
    1517     o    ,tra
    1618     o    ,h     ,lv     ,cpn    ,p,ph,tv    ,tp,tvp,clw
    … …  
    3941      real t1(len,nd),q1(len,nd),qs1(len,nd)
    4042      real t1_wake(len,nd),q1_wake(len,nd),qs1_wake(len,nd)
     43      real s1_wake(len)
    4144      real u1(len,nd),v1(len,nd)
    4245      real gz1(len,nd),th1(len,nd),th1_wake(len,nd)
    … …  
    5861      real t(len,nd),q(len,nd),qs(len,nd)
    5962      real t_wake(len,nd),q_wake(len,nd),qs_wake(len,nd)
     63      real s_wake(len)
    6064      real u(len,nd),v(len,nd)
    6165      real gz(len,nd),th(len,nd),th_wake(len,nd)
    … …  
    131135      if(iflag1(i).eq.0)then
    132136      nn=nn+1
     137      s_wake(nn)=s1_wake(i)
    133138      iflag(nn)=iflag1(i)
    134139      nk(nn)=nk1(i)

  • LMDZ4/branches/LMDZ4-dev/libf/phylmd/cva_driver.F

    r1062 r1127  
    22     &                   iflag_con,iflag_mix,
    33     &                   iflag_clos,delt,
    4      &                   t1,q1,qs1,t1_wake,q1_wake,qs1_wake,
     4     &                   t1,q1,qs1,t1_wake,q1_wake,qs1_wake,s1_wake,
    55     &                   u1,v1,tra1,
    66     &                   p1,ph1,
    … …  
    5050C      q1_wake       Real           Input        specific hum(unsat draught envt)
    5151C      qs1_wake      Real           Input        sat specific hum(unsat draughts envt)
     52C      s1_wake       Real           Input        fractionnal area covered by wakes
    5253C      u1            Real           Input        u-wind
    5354C      v1            Real           Input        v-wind
    … …  
    121122      real q1_wake(len,nd)
    122123      real qs1_wake(len,nd)
     124      real s1_wake(len)
    123125      real u1(len,nd)
    124126      real v1(len,nd)
    … …  
    198200!       Must be defined at same grid levels as T.
    199201!
     202!s_wake: Array of fractionnal area occupied by the wakes.
     203!
    200204!  u:   Array of zonal wind velocity (m/s) of dimension ND, witth first
    201205!       index corresponding with the lowest model level. Defined at
    … …  
    358362      real t(nloc,klev),q(nloc,klev),qs(nloc,klev)
    359363      real t_wake(nloc,klev),q_wake(nloc,klev),qs_wake(nloc,klev)
     364      real s_wake(nloc)
    360365      real u(nloc,klev),v(nloc,klev)
    361366      real gz(nloc,klev),h(nloc,klev)     ,hm(nloc,klev)
    … …  
    531536        print*,'Emanuel version 3 nouvelle'
    532537       endif
    533 
     538!       print*,'t1, q1 ',t1,q1
    534539       CALL cv3_prelim(len,nd,ndp1,t1,q1,p1,ph1      ! nd->na
    535540     o               ,lv1,cpn1,tv1,gz1,h1,hm1,th1)
    … …  
    668673
    669674      if (iflag_con.eq.3) then
    670      
     675!       print*,'ncum tv1 ',ncum,tv1
     676!       print*,'tvp1 ',tvp1
    671677       CALL cv3a_compress( len,nloc,ncum,nd,ntra
    672678     :    ,iflag1,nk1,icb1,icbs1
    673679     :    ,plcl1,tnk1,qnk1,gznk1,hnk1,unk1,vnk1
    674680     :    ,wghti1,pbase1,buoybase1
    675      :    ,t1,q1,qs1,t1_wake,q1_wake,qs1_wake,u1,v1,gz1,th1,th1_wake
     681     :    ,t1,q1,qs1,t1_wake,q1_wake,qs1_wake,s1_wake
     682     :    ,u1,v1,gz1,th1,th1_wake
    676683     :    ,tra1
    677684     :    ,h1     ,lv1     ,cpn1   ,p1,ph1,tv1    ,tp1,tvp1,clw1
    … …  
    682689     o    ,plcl,tnk,qnk,gznk,hnk,unk,vnk
    683690     o    ,wghti,pbase,buoybase
    684      o    ,t,q,qs,t_wake,q_wake,qs_wake,u,v,gz,th,th_wake
     691     o    ,t,q,qs,t_wake,q_wake,qs_wake,s_wake
     692     o    ,u,v,gz,th,th_wake
    685693     o    ,tra
    686694     o    ,h     ,lv     ,cpn    ,p,ph,tv    ,tp,tvp,clw
    … …  
    688696     o    ,sig,w0,ptop2
    689697     o    ,Ale,Alp  )
     698
     699!       print*,'tv ',tv
     700!       print*,'tvp ',tvp
    690701
    691702      endif
    … …  
    856867       CALL cv3_yield(nloc,ncum,nd,nd,ntra            ! na->nd
    857868     :                     ,icb,inb,delt
    858      :                     ,t,q,t_wake,q_wake,u,v,tra
     869     :                     ,t,q,t_wake,q_wake,s_wake,u,v,tra
    859870     :                     ,gz,p,ph,h,hp,lv,cpn,th,th_wake
    860871     :                     ,ep,clw,m,tp,mp,qp,up,vp,trap

  • LMDZ4/branches/LMDZ4-dev/libf/phylmd/wake.F

    r1059 r1127  
    1212     o                ,Cstar,d_deltat_gw
    1313     o                ,d_deltatw2,d_deltaqw2)
     14
    1415
    1516***************************************************************
    … …  
    157158      REAL alpk
    158159      REAL delta_t_min
    159       REAL Pupper
    160160      INTEGER nsub
    161161      REAL dtimesub
    162162      REAL sigmad, hwmin
     163      REAL :: sigmaw_max
    163164cIM 080208
    164165      LOGICAL, dimension(klon) :: gwake
    … …  
    183184      INTEGER, DIMENSION(klon) :: ktop, kupper
    184185
     186c Sub-timestep tendencies and related variables
     187       REAL d_deltatw(klon,klev),d_deltaqw(klon,klev)
     188       REAL d_te(klon,klev),d_qe(klon,klev)
     189       REAL d_sigmaw(klon),alpha(klon)
     190       REAL q0_min(klon),q1_min(klon)
     191       LOGICAL wk_adv(klon), OK_qx_qw(klon)
     192       REAL epsilon
     193       DATA epsilon/1.e-9/
     194
    185195c Autres variables internes
    186196      INTEGER isubstep, k, i
    … …  
    202212      REAL, DIMENSION(klon,klev) :: the, thu
    203213
    204       REAL, DIMENSION(klon,klev) :: d_deltatw, d_deltaqw
     214!      REAL, DIMENSION(klon,klev) :: d_deltatw, d_deltaqw
    205215
    206216      REAL, DIMENSION(klon,klev+1) :: omgbw
     217      REAL, DIMENSION(klon) :: pupper
    207218      REAL, DIMENSION(klon) :: omgtop
    208219      REAL, DIMENSION(klon,klev) :: dp_omgbw
    … …  
    279290        dqls(i,k) = 0.
    280291        d_deltat_gw(i,k)=0.
     292        d_te(i,k) = 0.
     293        d_qe(i,k) = 0.
     294        d_deltatw(i,k) = 0.
     295        d_deltaqw(i,k) = 0.
    281296!IM 060508 beg
    282297        d_deltatw2(i,k)=0.
    … …  
    294309      sigmaw(i) = amin1(sigmaw(i),0.99)
    295310      sigmaw0(i) = sigmaw(i)
     311      wape(i) = 0.
     312      wape2(i) = 0.
     313      d_sigmaw(i) = 0.
     314      ktopw(i) = 0
    296315      ENDDO
    297316C
    … …  
    406425c
    407426C       Pupper = 50000.  ! melting level
    408        Pupper = 60000.
     427c       Pupper = 60000.
    409428c       Pupper = 80000.  ! essais pour case_e
     429       DO i = 1,klon
     430ccc       Pupper(i) = 0.6*ph(i,1)
     431        Pupper(i) = 60000.
     432       ENDDO
     433
    410434C
    411435C    Determine Wake top pressure (Ptop) from buoyancy integral
    … …  
    481505      DO i=1,klon
    482506        IF (ph(i,k+1) .lt. ptop(i)) ktop(i)=k
    483         IF (ph(i,k+1) .lt. pupper) kupper(i)=k
     507        IF (ph(i,k+1) .lt. pupper(i)) kupper(i)=k
    484508      ENDDO
    485509      ENDDO
    … …  
    622646      ENDIF
    623647      ENDDO
    624 c
    625 C
     648
     649c
     650c Check qx and qw positivity
     651c --------------------------
     652      DO i = 1,klon
     653       q0_min(i)=min(  (qe(i,1)-sigmaw(i)*deltaqw(i,1)),
     654     $              (qe(i,1)+(1.-sigmaw(i))*deltaqw(i,1))  )
     655      ENDDO
     656      DO k = 2,klev
     657      DO i = 1,klon
     658        q1_min(i)=min(  (qe(i,k)-sigmaw(i)*deltaqw(i,k)),
     659     $              (qe(i,k)+(1.-sigmaw(i))*deltaqw(i,k))  )
     660        IF (q1_min(i).le.q0_min(i)) THEN
     661          q0_min(i)=q1_min(i)
     662        ENDIF
     663      ENDDO
     664      ENDDO
     665c
     666      DO i = 1,klon
     667       OK_qx_qw(i) = q0_min(i) .GE. 0.
     668       alpha(i) = 1.
     669      ENDDO
     670c
    626671CC -----------------------------------------------------------------
    627672C    Sub-time-stepping
    … …  
    634679      DO isubstep = 1,nsub
    635680c------------------------------------------------------------
    636       DO i=1,klon
     681c
     682c wk_adv is the logical flag enabling wake evolution in the time advance loop
     683      DO i = 1,klon
     684       wk_adv(i) = OK_qx_qw(i) .AND. alpha(i) .GE. 1.
     685      ENDDO
     686c
     687      DO i=1,klon
     688        IF (wk_adv(i)) THEN
    637689        gfl(i) = 2.*sqrt(3.14*wdens*sigmaw(i))
    638       ENDDO
    639       DO i=1,klon
    640         sigmaw(i) =sigmaw(i) + gfl(i)*Cstar(i)*dtimesub
    641         sigmaw(i) =amin1(sigmaw(i),0.99)     !!!!!!!!
     690        ENDIF
     691      ENDDO
     692      DO i=1,klon
     693        IF (wk_adv(i)) THEN
     694         d_sigmaw(i) = gfl(i)*Cstar(i)*dtimesub
     695c        sigmaw(i) =sigmaw(i) + gfl(i)*Cstar(i)*dtimesub
     696c        sigmaw(i) =min(sigmaw(i),0.99)     !!!!!!!!
    642697c        wdens = wdens0/(10.*sigmaw)
    643698c        sigmaw =max(sigmaw,sigd_con)
    644699c        sigmaw =max(sigmaw,sigmad)
     700        ENDIF
    645701      ENDDO
    646702C
    … …  
    650706cIM 060208 et omg sur les niveaux de 1 a klev+1, alors que avant l'on definit
    651707cIM 060208 au niveau k=1..?
    652       dp_deltomg(1:klon,1:klev)=0.
     708      DO k= 1,klev
     709      DO i = 1,klon
     710        dp_deltomg(i,k)=0.
     711      ENDDO
     712      ENDDO
    653713      DO k= 1,klev+1
    654714      DO i = 1,klon
    … …  
    658718c
    659719      DO i=1,klon
     720        IF (wk_adv(i)) THEN
    660721        z(i)= 0.
    661722        omg(i,1) = 0.
    662723        dp_deltomg(i,1) = -(gfl(i)*Cstar(i))/(sigmaw(i) * (1-sigmaw(i)))
     724        ENDIF
    663725      ENDDO
    664726c
    665727      DO k= 2,klev
    666728      DO i = 1,klon
    667        IF( k .LE. ktop(i)) THEN
     729       IF( wk_adv(i) .AND. k .LE. ktop(i)) THEN
    668730          dz(i) = -(ph(i,k)-ph(i,k-1))/(rho(i,k-1)*rg)
    669731          z(i) = z(i)+dz(i)
    … …  
    675737c
    676738      DO i = 1,klon
     739        IF (wk_adv(i)) THEN
    677740        dztop(i)=-(ptop(i)-ph(i,ktop(i)))/(rho(i,ktop(i))*rg)
    678741        ztop(i) = z(i)+dztop(i)
    679742        omgtop(i)=dp_deltomg(i,1)*ztop(i)
     743        ENDIF
    680744      ENDDO
    681745c
    … …  
    685749c
    686750       DO i=1,klon
     751        IF (wk_adv(i)) THEN
    687752        omgtop(i) = -rho(i,ktop(i))*rg*omgtop(i)
    688753        dp_deltomg(i,1) = omgtop(i)/(ptop(i)-ph(i,1))
     754        ENDIF
    689755       ENDDO
    690756c
    691757      DO k= 1,klev
    692758      DO i = 1,klon
    693        IF( k .LE. ktop(i)) THEN
     759       IF( wk_adv(i) .AND. k .LE. ktop(i)) THEN
    694760          omg(i,k) = - rho(i,k)*rg*omg(i,k)
    695761          dp_deltomg(i,k) = dp_deltomg(i,1)
    … …  
    701767
    702768      DO i=1,klon
    703       IF (kupper(i) .GT. ktop(i)) THEN
     769      IF ( wk_adv(i) .AND. kupper(i) .GT. ktop(i)) THEN
    704770        omg(i,kupper(i)+1) = - Rg*amdwn(i,kupper(i)+1)/sigmaw(i)
    705771     $                + Rg*amup(i,kupper(i)+1)/(1.-sigmaw(i))
    706772        dp_deltomg(i,kupper(i)) = (omgtop(i)-omg(i,kupper(i)+1))/
    707      $                     (ptop(i)-pupper)
     773     $                     (ptop(i)-pupper(i))
    708774      ENDIF
    709775      ENDDO
    … …  
    711777      DO k= 1,klev
    712778      DO i = 1,klon
    713        IF( k .GT. ktop(i) .AND. k .LE. kupper(i)) THEN
     779       IF( wk_adv(i) .AND. k .GT. ktop(i) .AND. k .LE. kupper(i)) THEN
    714780          dp_deltomg(i,k) = dp_deltomg(i,kupper(i))
    715781          omg(i,k) = omgtop(i)+(ph(i,k)-ptop(i))*dp_deltomg(i,kupper(i))
    … …  
    718784      ENDDO
    719785c
     786c
    720787c--    Compute wake average vertical velocity omgbw
    721788c
    … …  
    723790      DO k = 1,klev+1
    724791      DO i=1,klon
     792        IF ( wk_adv(i)) THEN
    725793        omgbw(i,k) = omgb(i,k)+(1.-sigmaw(i))*omg(i,k)
     794        ENDIF
    726795      ENDDO
    727796      ENDDO
    … …  
    730799      DO k = 1,klev
    731800      DO i=1,klon
     801        IF ( wk_adv(i)) THEN
    732802        dp_omgbw(i,k) = (omgbw(i,k+1)-omgbw(i,k))/(ph(i,k+1)-ph(i,k))
     803        ENDIF
    733804      ENDDO
    734805      ENDDO
    … …  
    739810      DO k = 1,klev
    740811      DO i=1,klon
    741        alpha_up(i,k) = 0.
    742        IF (omgb(i,k) .GT. 0.) alpha_up(i,k) = 1.
     812        IF ( wk_adv(i)) THEN
     813         alpha_up(i,k) = 0.
     814         IF (omgb(i,k) .GT. 0.) alpha_up(i,k) = 1.
     815        ENDIF
    743816      ENDDO
    744817      ENDDO
    … …  
    747820
    748821      DO i=1,klon
    749         RRe1(i) = 1.-sigmaw(i)
    750         RRe2(i) = sigmaw(i)
     822        IF ( wk_adv(i)) THEN
     823         RRe1(i) = 1.-sigmaw(i)
     824         RRe2(i) = sigmaw(i)
     825        ENDIF
    751826      ENDDO
    752827      RRd1 = -1.
    … …  
    757832      DO k= 1,klev
    758833      DO i = 1,klon
    759        IF(k .LE. kupper(i)+1) THEN
     834       IF( wk_adv(i) .AND. k .LE. kupper(i)+1) THEN
    760835        dth(i,k) = deltatw(i,k)/ppi(i,k)
    761836        Th1(i,k) = the(i,k) - sigmaw(i)     *dth(i,k)   ! undisturbed area
    … …  
    778853      DO k= 2,klev
    779854      DO i = 1,klon
    780        IF(k .LE. kupper(i)+1) THEN
     855       IF( wk_adv(i) .AND. k .LE. kupper(i)+1) THEN
    781856        D_Th1(i,k) = Th1(i,k-1)-Th1(i,k)
    782857        D_Th2(i,k) = Th2(i,k-1)-Th2(i,k)
    … …  
    790865
    791866      DO i=1,klon
    792         omgbdth(i,1) = 0.
    793         omgbdq(i,1) = 0.
     867        IF( wk_adv(i)) THEN
     868         omgbdth(i,1) = 0.
     869         omgbdq(i,1) = 0.
     870        ENDIF
    794871      ENDDO
    795872
    796873      DO k= 2,klev
    797874      DO i = 1,klon
    798        IF(k .LE. kupper(i)+1) THEN  !   loop on interfaces
     875       IF( wk_adv(i) .AND. k .LE. kupper(i)+1) THEN  !   loop on interfaces
    799876        omgbdth(i,k) = omgb(i,k)*(    dth(i,k-1) -     dth(i,k))
    800877        omgbdq(i,k)  = omgb(i,k)*(deltaqw(i,k-1) - deltaqw(i,k))
    … …  
    806883      DO k= 1,klev
    807884      DO i = 1,klon
    808        IF(k .LE. kupper(i)-1) THEN
     885       IF( wk_adv(i) .AND. k .LE. kupper(i)-1) THEN
    809886c-----------------------------------------------------------------
    810887c
    … …  
    829906c   and increment large scale tendencies
    830907c
    831          dtls(i,k) = dtls(i,k) +
    832      $               dtimesub*(
     908
     909c
     910C
     911CC -----------------------------------------------------------------
     912         d_te(i,k) =  dtimesub*(
    833913     $        ( RRe1(i)*omg(i,k  )*sigmaw(i)     *D_Th1(i,k)
    834914     $         -RRe2(i)*omg(i,k+1)*(1.-sigmaw(i))*D_Th2(i,k+1) )
    … …  
    836916     $         -sigmaw(i)*(1.-sigmaw(i))*dth(i,k)*dp_deltomg(i,k)
    837917     $                      )*ppi(i,k)
    838 c         print*,'dtls=',dtls(i,k)
    839 c
    840          dqls(i,k) = dqls(i,k) +
    841      $               dtimesub*(
     918c
     919         d_qe(i,k) =  dtimesub*(
    842920     $        ( RRe1(i)*omg(i,k  )*sigmaw(i)     *D_q1(i,k)
    843921     $         -RRe2(i)*omg(i,k+1)*(1.-sigmaw(i))*D_q2(i,k+1) )
    … …  
    845923     $         -sigmaw(i)*(1.-sigmaw(i))*deltaqw(i,k)*dp_deltomg(i,k)
    846924     $                      )
    847 c         print*,'dqls=',dqls(k)
    848        ENDIF
     925       ENDIF
     926
    849927c-------------------------------------------------------------------
    850928      ENDDO
    … …  
    856934      DO k= 1,klev
    857935      DO i = 1,klon
    858        IF(k .LE. kupper(i)-1) THEN
     936       IF( wk_adv(i) .AND. k .LE. kupper(i)-1) THEN
    859937c
    860938c Coefficient de répartition
    … …  
    912990
    913991        IF (dtimesub*Tgw(i,k).lt.1.e-10) THEN
    914           deltatw(i,k) = deltatw(i,k)+dtimesub*
    915      $          (ff(i)+dtKE(i,k)+dtPBL(i,k) 
     992          d_deltatw(i,k) = dtimesub*
     993     $          (ff(i)+dtKE(i,k)+dtPBL(i,k)
    916994     $          - spread(i,k)*deltatw(i,k)-Tgw(i,k)*deltatw(i,k))
    917995        ELSE
    918            deltatw(i,k) = deltatw(i,k)+1/Tgw(i,k)*(1-exp(-dtimesub*
     996           d_deltatw(i,k) = 1/Tgw(i,k)*(1-exp(-dtimesub*
    919997     $          Tgw(i,k)))*
    920998     $          (ff(i)+dtKE(i,k)+dtPBL(i,k)
    921999     $          - spread(i,k)*deltatw(i,k)-Tgw(i,k)*deltatw(i,k))
    9221000        ENDIF
    923    
     1001
    9241002        dth(i,k) = deltatw(i,k)/ppi(i,k)
    9251003
    9261004        gg(i)=d_deltaqw(i,k)/dtimesub
    9271005
    928        deltaqw(i,k) = deltaqw(i,k) +
    929      $         dtimesub*(gg(i)+ dqKE(i,k)+dqPBL(i,k) - spread(i,k)*
    930      $         deltaqw(i,k))
     1006       d_deltaqw(i,k) = dtimesub*(gg(i)+ dqKE(i,k)+dqPBL(i,k)
     1007     $                            - spread(i,k)*deltaqw(i,k))
    9311008
    9321009       d_deltatw2(i,k)=d_deltatw2(i,k)+d_deltatw(i,k)
    … …  
    9361013      ENDDO
    9371014
    938 C   And update large scale variables
     1015C
     1016C   Scale tendencies so that water vapour remains positive in w and x.
     1017C
     1018      call wake_vec_modulation(klon,klev,wk_adv,epsilon,qe,d_qe,deltaqw,
     1019     $                d_deltaqw,sigmaw,d_sigmaw,alpha)
     1020c
     1021      DO k = 1,klev
     1022      DO i = 1,klon
     1023       IF( wk_adv(i) .AND. k .LE. kupper(i)) THEN
     1024        d_te(i,k)=alpha(i)*d_te(i,k)
     1025        d_qe(i,k)=alpha(i)*d_qe(i,k)
     1026        d_deltatw(i,k)=alpha(i)*d_deltatw(i,k)
     1027        d_deltaqw(i,k)=alpha(i)*d_deltaqw(i,k)
     1028        d_deltat_gw(i,k)=alpha(i)*d_deltat_gw(i,k)
     1029       ENDIF
     1030      ENDDO
     1031      ENDDO
     1032      DO i = 1,klon
     1033       IF( wk_adv(i)) THEN
     1034        d_sigmaw(i)=alpha(i)*d_sigmaw(i)
     1035       ENDIF
     1036      ENDDO
     1037
     1038C   Update large scale variables and wake variables
    9391039cIM 060208 manque DO i + remplace DO k=1,kupper(i)
    9401040cIM 060208     DO k = 1,kupper(i)
    9411041      DO k= 1,klev
    9421042      DO i = 1,klon
    943        IF(k .LE. kupper(i)) THEN
     1043       IF( wk_adv(i) .AND. k .LE. kupper(i)) THEN
     1044        dtls(i,k)=dtls(i,k)+d_te(i,k)
     1045        dqls(i,k)=dqls(i,k)+d_qe(i,k)
     1046       ENDIF
     1047      ENDDO
     1048      ENDDO
     1049      DO k= 1,klev
     1050      DO i = 1,klon
     1051       IF( wk_adv(i) .AND. k .LE. kupper(i)) THEN
    9441052        te(i,k) = te0(i,k) + dtls(i,k)
    9451053        qe(i,k) = qe0(i,k) + dqls(i,k)
    9461054        the(i,k) = te(i,k)/ppi(i,k)
    947        ENDIF
    948       ENDDO
     1055        deltatw(i,k) = deltatw(i,k)+d_deltatw(i,k)
     1056        deltaqw(i,k) = deltaqw(i,k)+d_deltaqw(i,k)
     1057        dth(i,k) = deltatw(i,k)/ppi(i,k)
     1058       ENDIF
     1059      ENDDO
     1060      ENDDO
     1061      DO i = 1,klon
     1062       IF( wk_adv(i)) THEN
     1063        sigmaw(i) = sigmaw(i)+d_sigmaw(i)
     1064       ENDIF
    9491065      ENDDO
    9501066c
    … …  
    9561072c
    9571073      DO i=1,klon
    958       Ptop_provis(i)=ph(i,1)
     1074       IF ( wk_adv(i)) THEN
     1075        Ptop_provis(i)=ph(i,1)
     1076       ENDIF
    9591077      ENDDO
    9601078c
    9611079      DO k= 2,klev
    9621080      DO i=1,klon
    963         IF (Ptop_provis(i) .EQ. ph(i,1) .AND.
     1081        IF ( wk_adv(i) .AND.
     1082     $       Ptop_provis(i) .EQ. ph(i,1) .AND.
    9641083     $      dth(i,k) .GT. -delta_t_min .and.
    9651084     $      dth(i,k-1).LT. -delta_t_min) THEN
    … …  
    9811100      DO k = 1,klev
    9821101      DO i=1,klon
     1102       IF ( wk_adv(i)) THEN
    9831103        dz(i) = -(amax1(ph(i,k+1),Ptop_provis(i))-Ph(i,k))/(rho(i,k)*rg)
    9841104        IF (dz(i) .gt. 0) THEN
    … …  
    9871107         dthmin(i) = amin1(dthmin(i),dth(i,k))
    9881108        ENDIF
     1109       ENDIF
    9891110      ENDDO
    9901111      ENDDO
    … …  
    9931114
    9941115      DO i=1,klon
    995        hw(i) = 2.*sum_dth(i)/amin1(dthmin(i),-0.5)
    996        hw(i) = amax1(hwmin,hw(i))
     1116       IF ( wk_adv(i)) THEN
     1117         hw(i) = 2.*sum_dth(i)/amin1(dthmin(i),-0.5)
     1118         hw(i) = amax1(hwmin,hw(i))
     1119       ENDIF
    9971120      ENDDO
    9981121c
    … …  
    10061129      DO k = 1,klev
    10071130      DO i=1,klon
     1131       IF ( wk_adv(i)) THEN
    10081132        dz(i) = amin1(-(ph(i,k+1)-Ph(i,k))/(rho(i,k)*rg),hw(i)-z(i))
    10091133        IF (dz(i) .gt. 0) THEN
    … …  
    10121136         ktop(i) = k
    10131137        ENDIF
     1138       ENDIF
    10141139      ENDDO
    10151140      ENDDO
    … …  
    10181143c
    10191144      DO i=1,klon
     1145       IF ( wk_adv(i)) THEN
    10201146        Ptop_new(i)=ptop(i)
     1147       ENDIF
    10211148      ENDDO
    10221149c
    … …  
    10241151      DO i=1,klon
    10251152cIM v3JYG; IF (k .GE. ktop(i)
    1026        IF (k .LE. ktop(i) .AND.
     1153       IF ( wk_adv(i) .AND.
     1154     $      k .LE. ktop(i) .AND.
    10271155     $      ptop_new(i) .EQ. ptop(i) .AND.
    10281156     $      dth(i,k) .GT. -delta_t_min .and.
    … …  
    10371165c
    10381166      DO i=1,klon
    1039       ptop(i) = ptop_new(i)
     1167       IF ( wk_adv(i)) THEN
     1168        ptop(i) = ptop_new(i)
     1169       ENDIF
    10401170      ENDDO
    10411171
    … …  
    10501180      DO k = 1,klev
    10511181      DO i=1,klon
    1052         IF (k .GE. kupper(i)) THEN
     1182        IF ( wk_adv(i) .AND. k .GE. kupper(i)) THEN
    10531183         deltatw(i,k) = 0.
    10541184         deltaqw(i,k) = 0.
    … …  
    10581188c
    10591189C
     1190c-------------Cstar computation---------------------------------
     1191      DO i=1, klon
     1192      sum_thu(i) = 0.
     1193      sum_tu(i) = 0.
     1194      sum_qu(i) = 0.
     1195      sum_thvu(i) = 0.
     1196      sum_dth(i) = 0.
     1197      sum_dq(i) = 0.
     1198      sum_rho(i) = 0.
     1199      sum_dtdwn(i) = 0.
     1200      sum_dqdwn(i) = 0.
     1201
     1202      av_thu(i) = 0.
     1203      av_tu(i) =0.
     1204      av_qu(i) =0.
     1205      av_thvu(i) = 0.
     1206      av_dth(i) = 0.
     1207      av_dq(i) = 0.
     1208      av_rho(i) =0.
     1209      av_dtdwn(i) =0.
     1210      av_dqdwn(i) = 0.
     1211      ENDDO
     1212C
     1213C Integrals (and wake top level number)
     1214C --------------------------------------
     1215C
     1216C Initialize sum_thvu to 1st level virt. pot. temp.
     1217
     1218      DO i=1,klon
     1219      z(i) = 1.
     1220      dz(i) = 1.
     1221      sum_thvu(i) =  thu(i,1)*(1.+eps*qu(i,1))*dz(i)
     1222      sum_dth(i) = 0.
     1223      ENDDO
     1224
     1225      DO k = 1,klev
     1226      DO i=1,klon
     1227        dz(i) = -(max(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
     1228        IF (dz(i) .GT. 0) THEN
     1229         z(i) = z(i)+dz(i)
     1230         sum_thu(i) = sum_thu(i) + thu(i,k)*dz(i)
     1231         sum_tu(i) = sum_tu(i) + tu(i,k)*dz(i)
     1232         sum_qu(i) = sum_qu(i) + qu(i,k)*dz(i)
     1233         sum_thvu(i) = sum_thvu(i) + thu(i,k)*(1.+eps*qu(i,k))*dz(i)
     1234         sum_dth(i) = sum_dth(i) + dth(i,k)*dz(i)
     1235         sum_dq(i) = sum_dq(i) + deltaqw(i,k)*dz(i)
     1236         sum_rho(i) = sum_rho(i) + rhow(i,k)*dz(i)
     1237         sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i,k)*dz(i)
     1238         sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i,k)*dz(i)
     1239        ENDIF
     1240      ENDDO
     1241      ENDDO
     1242c
     1243      DO i=1,klon
     1244        hw0(i) = z(i)
     1245      ENDDO
     1246c
     1247C
     1248C - WAPE and mean forcing computation
     1249C ---------------------------------------
     1250C
     1251C ---------------------------------------
     1252C
     1253C Means
     1254
     1255      DO i=1,klon
     1256       av_thu(i) = sum_thu(i)/hw0(i)
     1257       av_tu(i) = sum_tu(i)/hw0(i)
     1258       av_qu(i) = sum_qu(i)/hw0(i)
     1259       av_thvu(i) = sum_thvu(i)/hw0(i)
     1260       av_dth(i) = sum_dth(i)/hw0(i)
     1261       av_dq(i) = sum_dq(i)/hw0(i)
     1262       av_rho(i) = sum_rho(i)/hw0(i)
     1263       av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
     1264       av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
     1265c
     1266       wape(i) = - rg*hw0(i)*(av_dth(i)
     1267     $     + eps*(av_thu(i)*av_dq(i)+av_dth(i)*av_qu(i)+av_dth(i)*
     1268     $     av_dq(i) ))/av_thvu(i)
     1269      ENDDO
     1270C
     1271C Filter out bad wakes
     1272
     1273      DO k = 1,klev
     1274       DO i=1,klon
     1275        IF ( wape(i) .LT. 0.) THEN
     1276          deltatw(i,k) = 0.
     1277          deltaqw(i,k) = 0.
     1278          dth(i,k) = 0.
     1279        ENDIF
     1280       ENDDO
     1281      ENDDO
     1282c
     1283      DO i=1,klon
     1284      IF ( wape(i) .LT. 0.) THEN
     1285        wape(i) = 0.
     1286        Cstar(i) = 0.
     1287        hw(i) = hwmin
     1288        sigmaw(i) = max(sigmad,sigd_con(i))
     1289        fip(i) = 0.
     1290        gwake(i) = .FALSE.
     1291      ELSE
     1292        Cstar(i) = stark*sqrt(2.*wape(i))
     1293        gwake(i) = .TRUE.
     1294      ENDIF
     1295      ENDDO
     1296
    10601297       ENDDO      ! end sub-timestep loop
    10611298C
    … …  
    10651302      DO k = 1,klev
    10661303      DO i=1,klon
    1067         IF (k .LE. kupper(i)-1) THEN
     1304        IF ( wk_adv(i) .AND. k .LE. kupper(i)-1) THEN
    10681305         dtls(i,k) = dtls(i,k)/dtime
    10691306         dqls(i,k) = dqls(i,k)/dtime
    … …  
    11111348      DO k =1,klev
    11121349      DO i=1,klon
     1350       IF ( wk_adv(i)) THEN
    11131351        rho(i,k) = p(i,k)/(rd*te(i,k))
    11141352        IF(k .eq. 1) THEN
    … …  
    11251363        rhow(i,k) = p(i,k)/(rd*(te(i,k)+deltatw(i,k)))
    11261364        dth(i,k) = deltatw(i,k)/ppi(i,k)
     1365       ENDIF
    11271366      ENDDO
    11281367      ENDDO
    … …  
    11341373
    11351374      DO i=1,klon
     1375       IF ( wk_adv(i)) THEN
    11361376        z(i) = 1.
    11371377        dz(i) = 1.
    11381378        sum_thvu(i) =  thu(i,1)*(1.+eps*qu(i,1))*dz(i)
    11391379        sum_dth(i) = 0.
     1380      ENDIF
    11401381      ENDDO
    11411382
    11421383      DO k = 1,klev
    11431384      DO i=1,klon
     1385       IF ( wk_adv(i)) THEN
    11441386        dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
    11451387        IF (dz(i) .GT. 0) THEN
    … …  
    11551397         sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i,k)*dz(i)
    11561398        ENDIF
    1157       ENDDO
    1158       ENDDO
    1159 c
    1160       DO i=1,klon
     1399       ENDIF
     1400      ENDDO
     1401      ENDDO
     1402c
     1403      DO i=1,klon
     1404       IF ( wk_adv(i)) THEN
    11611405        hw0(i) = z(i)
    1162       ENDDO
    1163 c
    1164 C 2.1 - WAPE and mean forcing computation
     1406       ENDIF
     1407      ENDDO
     1408c
     1409C - WAPE and mean forcing computation
    11651410C-------------------------------------------------------------
    11661411
    … …  
    11681413
    11691414      DO i=1, klon
     1415       IF ( wk_adv(i)) THEN
    11701416        av_thu(i) = sum_thu(i)/hw0(i)
    11711417        av_tu(i) = sum_tu(i)/hw0(i)
    … …  
    11811427     $     + eps*(av_thu(i)*av_dq(i)+av_dth(i)*av_qu(i)+
    11821428     $     av_dth(i)*av_dq(i) ))/av_thvu(i)
    1183       ENDDO
    1184 
    1185 C 2.2 Prognostic variable update
     1429       ENDIF
     1430      ENDDO
     1431
     1432C Prognostic variable update
    11861433C ------------------------------------------------------------
    11871434
    … …  
    11901437      DO k = 1,klev
    11911438      DO i=1,klon
    1192         IF ( wape2(i) .LT. 0.) THEN
     1439        IF ( wk_adv(i) .AND. wape2(i) .LT. 0.) THEN
    11931440          deltatw(i,k) = 0.
    11941441          deltaqw(i,k) = 0.
    … …  
    12001447
    12011448      DO i=1, klon
    1202       IF ( wape2(i) .LT. 0.) THEN
     1449       IF ( wk_adv(i)) THEN
     1450       IF ( wape2(i) .LT. 0.) THEN
    12031451        wape2(i) = 0.
    12041452        Cstar2(i) = 0.
    … …  
    12121460        gwake(i) = .TRUE.
    12131461      ENDIF
     1462      ENDIF
    12141463      ENDDO
    12151464c
    12161465      DO i=1, klon
     1466       IF ( wk_adv(i)) THEN
    12171467        ktopw(i) = ktop(i)
     1468       ENDIF
    12181469      ENDDO
    12191470c
    12201471      DO i=1, klon
    1221       IF (ktopw(i) .gt. 0 .and. gwake(i)) then
     1472       IF ( wk_adv(i)) THEN
     1473       IF (ktopw(i) .gt. 0 .and. gwake(i)) then
    12221474
    12231475Cjyg1     Utilisation d'un h_efficace constant ( ~ feeding layer)
    … …  
    12341486         FIP(i) = 0.
    12351487       ENDIF
     1488       ENDIF
    12361489      ENDDO
    12371490c
    … …  
    12411494C              alors il disparait en se mélangeant à la partie undisturbed
    12421495c
     1496      sigmaw_max = 0.9
    12431497      DO k = 1,klev
    12441498       DO i=1, klon
    1245          IF ((sigmaw(i).GT.0.9).or.
     1499c correction NICOLAS     $     ((wape(i).ge.wape2(i)).and.(wape2(i).le.1.0))) THEN
     1500!         print*,'wape wape2 ktopw OK_qx_qw =',
     1501!     $           wape(i),wape2(i),ktopw(i),OK_qx_qw(i)
     1502         IF ((sigmaw(i).GT.sigmaw_max).or.
    12461503     $     ((wape(i).ge.wape2(i)).and.(wape2(i).le.1.0)).or.
    1247      $      (ktopw(i).le.2)) THEN
     1504     $      (ktopw(i).le.2) .OR.
     1505     $     .not. OK_qx_qw(i)) THEN
    12481506cIM cf NR/JYG 251108  $     ((wape(i).ge.wape2(i)).and.(wape2(i).le.1.0))) THEN
    12491507ccc      IF (sigmaw(i).GT.0.9) THEN
    … …  
    12571515c
    12581516      DO i=1, klon
    1259          IF ((sigmaw(i).GT.0.9).or.
    1260      $     ((wape(i).ge.wape2(i)).and.(wape2(i).le.1.0))) THEN
     1517         IF ( (sigmaw(i).GT.sigmaw_max).or.
     1518     $      ((wape(i).ge.wape2(i)).and.(wape2(i).le.1.0)).or.
     1519     $      (ktopw(i).le.2) .OR.
     1520     $     .not. OK_qx_qw(i)) THEN
     1521! correction NICOLAS     $     ((wape(i).ge.wape2(i)).and.(wape2(i).le.1.0))) THEN
    12611522ccc      IF (sigmaw(i).GT.0.9) THEN
    12621523         wape(i) = 0.
    … …  
    12721533      RETURN
    12731534      END
     1535
     1536      SUBROUTINE wake_vec_modulation(nlon,nl,wk_adv,epsilon,qe,d_qe,
     1537     $           deltaqw,d_deltaqw,sigmaw,d_sigmaw,alpha)
     1538c------------------------------------------------------
     1539cDtermination du coefficient alpha tel que les tendances
     1540c corriges alpha*d_G, pour toutes les grandeurs G, correspondent
     1541c a une humidite positive dans la zone (x) et dans la zone (w).
     1542c------------------------------------------------------
     1543c
     1544 
     1545c  Input
     1546      REAL qe(nlon,nl),d_qe(nlon,nl)
     1547      REAL deltaqw(nlon,nl),d_deltaqw(nlon,nl)
     1548      REAL sigmaw(nlon),d_sigmaw(nlon)
     1549      LOGICAL wk_adv(nlon)
     1550      INTEGER nl,nlon
     1551c  Output
     1552      REAL alpha(nlon)
     1553c  Internal variables
     1554      REAL alpha1(nlon)
     1555      REAL x,a,b,c,discrim,zeta(nlon)
     1556      REAL epsilon
     1557!      DATA epsilon/1.e-9/
     1558c
     1559      DO k=1,nl
     1560      DO i = 1,nlon
     1561       IF (wk_adv(i)) THEN
     1562        IF ((deltaqw(i,k)+d_deltaqw(i,k)).ge.0.) then
     1563         zeta(i)=0.
     1564        ELSE
     1565         zeta(i)=1.
     1566        END IF
     1567       ENDIF
     1568      ENDDO
     1569      DO i = 1,nlon
     1570       IF (wk_adv(i)) THEN
     1571        x = qe(i,k)+(zeta(i)-sigmaw(i))*deltaqw(i,k)
     1572     $   +d_qe(i,k)+(zeta(i)-sigmaw(i))*d_deltaqw(i,k)
     1573     $   -d_sigmaw(i)*(deltaqw(i,k)+d_deltaqw(i,k))
     1574      a=-d_sigmaw(i)*d_deltaqw(i,k)
     1575      b=d_qe(i,k)+(zeta(i)-sigmaw(i))*d_deltaqw(i,k)
     1576     $           -deltaqw(i,k)*d_sigmaw(i)
     1577!      c=qe(i,k)+(zeta(i)-sigmaw(i))*deltaqw(i,k)-epsilon
     1578       c=qe(i,k)+(zeta(i)-sigmaw(i))*deltaqw(i,k)
     1579
     1580      discrim=b*b-4.*a*c
     1581!       print*,'ZETA *********************' 
     1582!       print*,'zeta sigmaw ',zeta(:)
     1583!       print*,'SIGMA *********************'
     1584!       print*,'sigmaw ',sigmaw(:)
     1585
     1586!       print*,' x ************************'
     1587!       print*,'x ',x
     1588!       print*,'  a+b ************************'
     1589!       print*,'a+b ',a+b
     1590
     1591!       print*,'a b c delta zeta ',a,b,c,discrim
     1592        IF (a+b .GE. 0.) THEN
     1593         alpha1(i)=1.
     1594        ELSE
     1595         IF (x .GE. 0.) THEN
     1596            alpha1(i)=1.
     1597         ELSE
     1598!              IF (a .GE. 0.) THEN
     1599              IF (a .GT. 0.) THEN
     1600!       print*,'a b c delta zeta ',a,b,c,discrim,zeta(i)
     1601!       print*,'-b+sqrt(discrim) ',-b+sqrt(discrim)
     1602                 alpha1(i)=0.9*min(   (2.*c)/(-b+sqrt(discrim)),
     1603     $                        (-b+sqrt(discrim))/(2.*a)   )
     1604              ELSE IF (a.eq.0.) THEN
     1605                 alpha1(i)=0.9*(-c/b)
     1606              ELSE
     1607!       print*,'a b c delta zeta ',a,b,c,discrim,zeta(i)
     1608!       print*,'-b+sqrt(discrim) ',-b+sqrt(discrim)
     1609                 alpha1(i)=0.9*max(   (2.*c)/(-b+sqrt(discrim)),
     1610     $                        (-b+sqrt(discrim))/(2.*a)   )
     1611              ENDIF
     1612         ENDIF
     1613        ENDIF
     1614       ENDIF
     1615      ENDDO
     1616      ENDDO
     1617c
     1618      DO i = 1,nlon
     1619       IF (wk_adv(i)) THEN
     1620        alpha(i) = min(alpha(i),alpha1(i))
     1621       ENDIF
     1622      ENDDO
     1623c
     1624      return
     1625      end
     1626
    12741627      Subroutine WAKE_scal (p,ph,ppi,dtime,sigd_con
    12751628     :                ,te0,qe0,omgb
    … …  
    23152668C              alors il disparait en se mélangeant à la partie undisturbed
    23162669
     2670! correction NICOLAS     .     ((wape.ge.wape2).and.(wape2.le.1.0))) THEN
    23172671      IF ((sigmaw.GT.0.9).or.
    23182672     .     ((wape.ge.wape2).and.(wape2.le.1.0)).or.(ktopw.le.2)) THEN
Note: See TracChangeset for help on using the changeset viewer.