Changeset 1322 for LMDZ4/branches


Ignore:
Timestamp:
Mar 12, 2010, 11:54:11 AM (14 years ago)
Author:
Laurent Fairhead
Message:

Improvements concerning wake parametrisation (from JYG, NR, IT, with more to come).
Alp_offset is read in form physiq.def file

Améliorations à la paramétrisation des poches froides (de JYG, NR, IT, d'autres
sont à venir)
Alp_offset est rajouté à la liste des paramètres lus dans physiq.def

Location:
LMDZ4/branches/LMDZ4V5.0-dev/libf
Files:
14 edited

Legend:

Unmodified
Added
Removed

  • LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3d/etat0_netcdf.F

    r1299 r1322  
    144144      REAL      :: solarlong0
    145145      real :: seuil_inversion
     146      real :: alp_offset
    146147
    147148      integer  read_climoz ! read ozone climatology
    … …  
    184185     &                 iflag_thermals,nsplit_thermals,tau_thermals,     &
    185186     &                 iflag_thermals_ed,iflag_thermals_optflux,        &
    186      &                 iflag_coupl,iflag_clos,iflag_wake, read_climoz )
     187     &                 iflag_coupl,iflag_clos,iflag_wake, read_climoz,  &
     188     &                 alp_offset)
    187189
    188190! co2_ppm0 : initial value of atmospheric CO2 from .def file (co2_ppm value)

  • LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/etat0_netcdf.F

    r1299 r1322  
    144144      REAL      :: solarlong0
    145145      real :: seuil_inversion
     146      REAL :: alp_offset
    146147
    147148      integer  read_climoz ! read ozone climatology
    … …  
    184185     &                 iflag_thermals,nsplit_thermals,tau_thermals,     &
    185186     &                 iflag_thermals_ed,iflag_thermals_optflux,        &
    186      &                 iflag_coupl,iflag_clos,iflag_wake, read_climoz )
     187     &                 iflag_coupl,iflag_clos,iflag_wake, read_climoz,  &
     188     &                 alp_offset)
    187189
    188190! co2_ppm0 : initial value of atmospheric CO2 from .def file (co2_ppm value)

  • LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/leapfrog_p.F

    r1299 r1322  
    677677         call suspend_timer(timer_caldyn)
    678678
     679        if (prt_level >= 10) then
    679680         write(lunout,*)
    680681     &   'leapfrog_p: Entree dans la physique : Iteration No ',true_itau
     682        endif
    681683c$OMP END MASTER
    682684

  • LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/calwake.F

    r990 r1322  
     1!
     2! $Id$
     3!
    14      SUBROUTINE CALWAKE(paprs,pplay,dtime
    25     :             ,t,q,omgb
    … …  
    189192     $                ,Cstar,d_deltat_gw
    190193     $                ,d_deltatw,d_deltaqw)
    191 
    192       DO i=1,klon
    193        IF (ktopw(i) .GT. 0) THEN
    194          DO l=1,klev
     194c
     195      DO l=1,klev
     196       DO i=1,klon
     197        IF (ktopw(i) .GT. 0) THEN
    195198           wake_deltat(i,l)= dtw(i,l)
    196199           wake_deltaq(i,l)= dqw(i,l)
    … …  
    212215           wake_ddeltat(i,l) = d_deltatw(i,l)
    213216           wake_ddeltaq(i,l) = d_deltaqw(i,l)
    214          ENDDO
    215        ELSE
    216          DO l = 1,klev
     217        ELSE
    217218           wake_deltat(i,l)= 0.
    218219           wake_deltaq(i,l)= 0.
    … …  
    222223           wake_dtKE(i,l) = 0.
    223224           wake_dqKE(i,l) = 0.
     225           wake_dtPBL(i,l) = 0.
     226           wake_dqPBL(i,l) = 0.
    224227           wake_omg(i,l) = 0.
    225228           wake_dp_deltomg(i,l) = 0.
    … …  
    230233           undi_t(i,l)=te(i,l)
    231234           undi_q(i,l)=qe(i,l)
    232          ENDDO
    233        ENDIF
    234 
     235           wake_ddeltat(i,l) = 0.
     236           wake_ddeltaq(i,l) = 0.
     237        ENDIF
     238       ENDDO
     239      ENDDO
     240c
     241      DO i=1,klon
    235242       wake_h(i)= hw(i)
    236243       wake_s(i)= sigmaw(i)
    … …  
    241248       wake_Cstar(i) = Cstar(i)
    242249       wake_dens(i) = wdens(i)
    243 c
    244 cIM 290108 999  CONTINUE
    245 c
    246       ENDDO
     250      ENDDO
     251c
    247252      RETURN
    248253      END
     254
    249255      SUBROUTINE CALWAKE_scal(paprs,pplay,dtime
    250256     :             ,t,q,omgb

  • LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/conf_phys.F90

    r1299 r1322  
    2222                       iflag_thermals,nsplit_thermals,tau_thermals, &
    2323                       iflag_thermals_ed,iflag_thermals_optflux, &
    24                        iflag_coupl,iflag_clos,iflag_wake, read_climoz)
     24                       iflag_coupl,iflag_clos,iflag_wake, read_climoz, &
     25                       alp_offset)
    2526
    2627   use IOIPSL
    … …  
    110111  integer :: iflag_clos
    111112  integer :: iflag_wake
     113  real :: alp_offset
     114  REAL, SAVE :: alp_offset_omp
    112115  integer,SAVE :: iflag_coupl_omp,iflag_clos_omp,iflag_wake_omp
    113116  integer,SAVE :: iflag_cvl_sigd_omp
    … …  
    10321035  iflag_wake_omp = 0
    10331036  call getin('iflag_wake',iflag_wake_omp)
     1037
     1038!Config Key  = alp_offset
     1039!Config Desc = 
     1040!Config Def  = 0
     1041!Config Help =
     1042!
     1043  alp_offset_omp = 0.01
     1044  call getin('alp_offset',alp_offset_omp)
    10341045
    10351046!
    … …  
    14491460    iflag_clos = iflag_clos_omp
    14501461    iflag_wake = iflag_wake_omp
     1462    alp_offset = alp_offset_omp
    14511463    iflag_cvl_sigd = iflag_cvl_sigd_omp
    14521464    type_run = type_run_omp
    … …  
    16161628  write(numout,*)' Fmax = ', Fmax
    16171629  write(numout,*)' alphas = ', alphas
     1630  write(numout,*)' iflag_wake = ', iflag_wake
     1631  write(numout,*)' alp_offset = ', alp_offset
    16181632
    16191633  write(numout,*)' lonmin lonmax latmin latmax bilKP_ins =',&

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

    r1146 r1322  
     1!
     2! $Id$
     3!
    14        SUBROUTINE cv3_cine(nloc,ncum,nd,icb,inb
    25     :                      ,pbase,plcl,p,ph,tv,tvp
    3      :                      ,cina,cinb)
     6     :                      ,cina,cinb,plfc)
    47
    58***************************************************************
    … …  
    2629c
    2730c output
    28       real cina(nloc),cinb(nloc)
     31      real cina(nloc),cinb(nloc),plfc(nloc)
    2932c
    3033c local variables
    … …  
    3437      logical lswitch(nloc),lswitch1(nloc),lswitch2(nloc)
    3538      logical exist_lfc(nloc)
    36       real plfc(nloc)
    3739      real dpmax
    3840      real deltap,dcin

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

    r1299 r1322  
    7474c      dtcrit = -5.0
    7575c      tau    = 3000.
    76 cc      tau = 1800.
    77 c     tau= 2800.
    78       tau=8000.
     76      tau = 1800.
     77cc      tau=8000.
    7978      beta   = 1.0 - delt/tau
    8079      alpha1 = 1.5e-3
    … …  
    20942093cc---end jyg---
    20952094c
    2096 c--------retour à la formulation originale d''Emanuel.
     2095c--------retour à la formulation originale d'Emanuel.
    20972096      b6=bfac*50.*sigd(il)*(ph(il,i)-ph(il,i+1))*sigt*afac
    20982097      c6=water(il,i+1)+bfac*wdtrain(il)
    … …  
    21002099      if(c6.gt.0.0)then
    21012100       revap=0.5*(-b6+sqrt(b6*b6+4.*c6))
    2102        water(il,i)=revap*revap      !equation de conservation
     2101cjyg    Dans sa formulation originale, Emanuel calcule l'evaporation par:
     2102cc             evap(il,i)=sigt*afac*revap
     2103c     ce qui n'est pas correct. Dans cv_routines, la formulation a été modifiee.
     2104c     Ici,l'evaporation evap est simplement calculee par l'equation de
     2105c     conservation.
     2106       water(il,i)=revap*revap
    21032107      else
     2108cjyg----   Correction : si c6 <= 0, water(il,i)=0.
    21042109       water(il,i) = 0.
    21052110      endif
    … …  
    23452350      real esum(nloc), fsum(nloc), gsum(nloc), hsum(nloc)
    23462351      real th_wake(nloc,nd)
    2347       real alpha_qpos(nloc)
     2352      real alpha_qpos(nloc),alpha_qpos1(nloc)
    23482353      real qcond(nloc,nd), nqcond(nloc,nd), wa(nloc,nd)  ! cld
    23492354      real siga(nloc,nd), sax(nloc,nd), mac(nloc,nd)      ! cld
    … …  
    30503055       do il=1,ncum
    30513056        IF (iflag(il) .le. 1) THEN
     3057        IF (cvflag_grav) then
     3058        ex=0.01*grav*ment(il,inb(il),inb(il))
     3059     :      *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
     3060     :      /(ph(i l,inb(il))-ph(il,inb(il)+1))
     3061        ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
     3062        ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
     3063     :       +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
     3064     :          /(ph(il,inb(il)-1)-ph(il,inb(il)))
     3065        else
    30523066        ex=0.1*ment(il,inb(il),inb(il))
    30533067     :      *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
    … …  
    30573071     :       +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
    30583072     :          /(ph(il,inb(il)-1)-ph(il,inb(il)))
     3073        ENDIF   !cvflag grav
    30593074        ENDIF    !iflag
    30603075       enddo
    … …  
    31223137c        in order to ensure moisture positivity
    31233138      DO il = 1,ncum
     3139      alpha_qpos(il)=1.
    31243140       IF (iflag(il) .le. 1) THEN
    3125         alpha_qpos(il) = max(1. , -delt*fr(il,1)/
     3141        if (fr(il,1) .le. 0.) then
     3142            alpha_qpos(il) = max(alpha_qpos(il) ,
     3143     :     (-delt*fr(il,1))/
    31263144     :     (s_wake(il)*rr_wake(il,1)+(1.-s_wake(il))*rr(il,1)))
     3145        end if
    31273146       ENDIF
    31283147      ENDDO
    … …  
    31303149       DO il = 1,ncum
    31313150        IF (iflag(il) .le. 1) THEN
    3132         alpha_qpos(il) = max(alpha_qpos(il) , -delt*fr(il,i)/
     3151          IF (fr(il,i) .le. 0.) THEN
     3152           alpha_qpos1(il)=max(1. , (-delt*fr(il,i))/
    31333153     :     (s_wake(il)*rr_wake(il,i)+(1.-s_wake(il))*rr(il,i)))
     3154             IF (alpha_qpos1(il) .ge. alpha_qpos(il))
     3155     :           alpha_qpos(il)=alpha_qpos1(il)
     3156          ENDIF
    31343157        ENDIF
    31353158       ENDDO

  • LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/cv3p1_closure.F

    r1299 r1322  
    6363      integer nsupmax(nloc)
    6464      real supcrit,temp(nloc,nd)
    65       real P1(nloc),Pmin(nloc)
     65      real P1(nloc),Pmin(nloc),plfc(nloc)
    6666      real asupmax0(nloc)
    6767      logical ok(nloc)
    … …  
    385385      CALL cv3_cine (nloc,ncum,nd,icb,inb
    386386     :                      ,pbase,plcl,p,ph,tv,tvp
    387      :                      ,cina,cinb)
     387     :                      ,cina,cinb,plfc)
    388388c
    389389      DO il = 1,ncum
    … …  
    495495      do k= 1,nl
    496496       do il = 1,ncum
    497 !IM       IF (k .ge. icb(il) .and. k .le. inb(il)) THEN
    498         IF (k .ge. icb(il)+1 .and. k .le. inb(il)) THEN
     497!old       IF (k .ge. icb(il) .and. k .le. inb(il)) THEN
     498!IM        IF (k .ge. icb(il)+1 .and. k .le. inb(il)) THEN
     499       IF (k .ge. icb(il) .and. k .le. inb(il)         !cor jyg
     500     $     .and. icb(il)+1 .le. inb(il)) THEN          !cor jyg
    499501         cbmflim(il) = cbmflim(il)+MLIM(il,k)
    500502        ENDIF
    … …  
    548550        do il = 1,ncum
    549551         IF ( k .ge. icb(il)+1 .AND. k .le. inb(il)) THEN
    550          sig(il,k) = beta*sig(il,k)+(1.-beta)*coef(il)*siglim(il,k)
    551 cc         sig(il,k) = beta*sig(il,k)+siglim(il,k)
    552          w0(il,k) = beta*w0(il,k)  +(1.-beta)*wlim(il,k)
    553          AMU=SIG(il,k)*W0(il,k)
     552         amu=beta*sig(il,k)*w0(il,k)+
     553     :   (1.-beta)*coef(il)*siglim(il,k)*wlim(il,k)
     554         w0(il,k) = wlim(il,k)
     555         w0(il,k) =max(w0(il,k),1.e-10)
     556         sig(il,k)=amu/w0(il,k)
     557         sig(il,k)=min(sig(il,k),1.)
    554558cc         amu = 0.5*(SIG(il,k)+sigold(il,k))*W0(il,k)
    555559         M(il,k)=AMU*0.007*P(il,k)*(PH(il,k)-PH(il,k+1))/TV(il,k)

  • LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/ini_wake.F

    r970 r1322  
     1!
     2! $Id$
     3!
    14      SUBROUTINE INI_WAKE(wape,fip,it_wape_prescr,
    25     :     wape_prescr, fip_prescr, alp_bl_prescr, ale_bl_prescr)
    … …  
    2326c   wape_prescr    : valeur prescrite de la WAPE.
    2427c   fip_prescr     : valeur prescrite de la FIP.
     28c   ale_bl_prescr  : valeur prescrite de la Ale de PBL.
     29c   alp_bl_prescr  : valeur prescrite de la Alp de PBL.
    2530c
    2631c Variables internes
    … …  
    2934c   w  = WAPE lue
    3035c   f  = FIP lue
     36c   alebl  = Ale de PBL lue
     37c   alpbl  = Alp de PBL lue
    3138c
     39      include 'iniprint.h'
    3240cdeclarations
    3341      real ale_bl_prescr
    3442      real alp_bl_prescr
    3543      real it
    36 cCR: on rajoute ale et alp de la PBL precrits
    37 c     open (99,file='wake.data',form='formatted')
    38 c     read (99,*) it
    39 c     read (99,*) w
    40 c     read (99,*) f
    41 c     read (99,*) u
    42 c     read (99,*) p
    43 c     close (99)
    4444
    4545! FH A mettre si besoin dans physiq.def
    … …  
    4848      w=4.
    4949      f=0.1
    50       u=0.1
    51       p=4.
     50      alebl=4.
     51      alpbl=0.1
    5252c
    53       print *,' it,w ',it,w
     53cCR: on rajoute ale et alp de la PBL precrits
     54      open (99,file='wake.data',form='formatted',status='old',err=902)
     55      read (99,*) it
     56      read (99,*) w
     57      read (99,*) f
     58      read (99,*,end=901) alebl
     59      read (99,*,end=901) alpbl
     60901   close (99)
     61902   continue
     62c
     63      write(lunout,*)' it,wape ',it,wape
    5464      it_wape_prescr = it
    5565      if (w .lt. 0) then
    … …  
    6171      endif
    6272c
    63       print *,' u,p ',u,p
    64       alp_bl_prescr=u
    65       ale_bl_prescr=p
     73      write(lunout,*)' alebl, alpbl ',alebl,alpbl
     74      ale_bl_prescr=alebl
     75      alp_bl_prescr=alpbl
    6676      print *,'Initialisation de la poche : WAPE, FIP imposees ='
    6777     $               ,wape_prescr, fip_prescr

  • LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/phyetat0.F

    r1311 r1322  
    11!
    2 ! $Header$
     2! $Id$
    33!
    44c
    … …  
    228228     $          'coherente ', i, zmasq(i), pctsrf(i, is_ter)
    229229     $          ,pctsrf(i, is_lic)
     230            WRITE(*,*) 'Je force la coherence zmasq=fractint'
    230231            zmasq(i) = fractint(i)
    231232        ENDIF
    … …  
    238239     $          'coherente ', i, zmasq(i) , pctsrf(i, is_oce)
    239240     $          ,pctsrf(i, is_sic)
     241            WRITE(*,*) 'Je force la coherence zmasq=fractint'
    240242            zmasq(i) = fractint(i)
    241243        ENDIF
    … …  
    987989      PRINT*,'(ecart-type) wake_cstar:', xmin, xmax
    988990c
     991c wake_pe
     992c
     993      CALL get_field("WAKE_PE",wake_pe,found)
     994      IF (.NOT. found) THEN
     995         PRINT*, "phyetat0: Le champ <WAKE_PE> est absent"
     996         PRINT*, "Depart legerement fausse. Mais je continue"
     997         wake_pe=0.
     998      ENDIF
     999      xmin = 1.0E+20
     1000      xmax = -1.0E+20
     1001      xmin = MINval(wake_pe)
     1002      xmax = MAXval(wake_pe)
     1003      PRINT*,'(ecart-type) wake_pe:', xmin, xmax
     1004c
    9891005c wake_fip
    9901006c

  • LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/phyredem.F

    r1302 r1322  
    299299!!!!!!!!!!!!!!!!!!!! FIN TKE PBL !!!!!!!!!!!!!!!!!!!!!!!!!
    300300cIM ajout zmax0, f0, ema_work1, ema_work2
    301 cIM wake_deltat, wake_deltaq, wake_s, wake_cstar, wake_fip
     301cIM wake_deltat, wake_deltaq, wake_s, wake_cstar, wake_pe, wake_fip
    302302     
    303303      CALL put_field("ZMAX0","ZMAX0",zmax0)
    … …  
    318318      CALL put_field("WAKE_CSTAR","WAKE_CSTAR",wake_cstar)
    319319     
     320      CALL put_field("WAKE_PE","WAKE_PE",wake_pe)
     321
    320322      CALL put_field("WAKE_FIP","WAKE_FIP",wake_fip)
    321323

  • LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/phys_state_var_mod.F90

    r1299 r1322  
     1!
     2! $Id$
     3!
    14      MODULE phys_state_var_mod
    25! Variables sauvegardees pour le startphy.nc
    … …  
    169172! wake_Cstar  : vitesse d'etalement de la poche
    170173! wake_s      : fraction surfacique occupee par la poche froide
     174! wake_pe     : wake potential energy - WAPE
    171175! wake_fip    : Gust Front Impinging power - ALP
    172176! dt_wake, dq_wake: LS tendencies due to wake
    … …  
    179183      REAL,ALLOCATABLE,SAVE :: wake_s(:)
    180184!$OMP THREADPRIVATE(wake_s)
     185      REAL,ALLOCATABLE,SAVE :: wake_pe(:)
     186!$OMP THREADPRIVATE(wake_pe)
    181187      REAL,ALLOCATABLE,SAVE :: wake_fip(:)
    182188!$OMP THREADPRIVATE(wake_fip)
    … …  
    370376      ALLOCATE(wght_th(klon,klev))
    371377      ALLOCATE(wake_deltat(klon,klev), wake_deltaq(klon,klev))
    372       ALLOCATE(wake_Cstar(klon), wake_s(klon), wake_fip(klon))
     378      ALLOCATE(wake_Cstar(klon), wake_s(klon))
     379      ALLOCATE(wake_pe(klon), wake_fip(klon))
    373380      ALLOCATE(dt_wake(klon,klev), dq_wake(klon,klev))
    374381      ALLOCATE(pfrac_impa(klon,klev), pfrac_nucl(klon,klev))
    … …  
    464471      deallocate(lalim_conv, wght_th)
    465472      deallocate(wake_deltat, wake_deltaq)
    466       deallocate(wake_Cstar, wake_s, wake_fip)
     473      deallocate(wake_Cstar, wake_s, wake_pe, wake_fip)
    467474      deallocate(dt_wake, dq_wake)
    468475      deallocate(pfrac_impa, pfrac_nucl)

  • LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/physiq.F

    r1311 r1322  
    625625      REAL q_undi(klon,klev)               ! humidite moyenne dans la zone non perturbee
    626626c
    627       REAL wake_pe(klon)              ! Wake potential energy - WAPE
     627cjyg
     628ccc      REAL wake_pe(klon)              ! Wake potential energy - WAPE
    628629
    629630      REAL wake_gfl(klon)             ! Gust Front Length
    … …  
    641642      REAL dt_a(klon,klev)
    642643      REAL dq_a(klon,klev)
     644      REAL, SAVE :: alp_offset
     645c$OMP THREADPRIVATE(alp_offset)
     646
    643647c
    644648cRR:fin declarations poches froides
    … …  
    11881192     .     iflag_thermals_ed,iflag_thermals_optflux,
    11891193c     nv flags pour la convection et les poches froides
    1190      .     iflag_coupl,iflag_clos,iflag_wake, read_climoz)
     1194     .     iflag_coupl,iflag_clos,iflag_wake, read_climoz,
     1195     &     alp_offset)
    11911196      call phys_state_var_init(read_climoz)
    11921197      print*, '================================================='
    … …  
    20212026c
    20222027ccalcul de ale_wake et alp_wake
    2023        do i = 1,klon
    2024           if (iflag_wake.eq.1) then
    2025           ale_wake(i) = 0.5*wake_cstar(i)**2
    2026           alp_wake(i) = wake_fip(i)
    2027           else
    2028           ale_wake(i) = 0.
    2029           alp_wake(i) = 0.
    2030           endif
    2031        enddo
     2028       if (iflag_wake.eq.1) then
     2029         if (itap .le. it_wape_prescr) then
     2030          do i = 1,klon
     2031           ale_wake(i) = wape_prescr
     2032           alp_wake(i) = fip_prescr
     2033          enddo
     2034         else
     2035          do i = 1,klon
     2036cjyg  ALE=WAPE au lieu de ALE = 1/2 Cstar**2
     2037ccc           ale_wake(i) = 0.5*wake_cstar(i)**2
     2038           ale_wake(i) = wake_pe(i)
     2039           alp_wake(i) = wake_fip(i)
     2040          enddo
     2041         endif
     2042       else
     2043         do i = 1,klon
     2044           ale_wake(i) = 0.
     2045           alp_wake(i) = 0.
     2046         enddo
     2047       endif
    20322048ccombinaison avec ale et alp de couche limite: constantes si pas de couplage, valeurs calculees
    20332049cdans le thermique sinon
    20342050       if (iflag_coupl.eq.0) then
    2035           if (debut) print*,'ALE et ALP imposes'
     2051          if (debut.and.prt_level.gt.9)
     2052     $                     WRITE(lunout,*)'ALE et ALP imposes'
    20362053          do i = 1,klon
    20372054con ne couple que ale
    … …  
    20462063          do i = 1,klon
    20472064              ALE(i) = max(ale_wake(i),Ale_bl(i))
    2048               ALP(i) = alp_wake(i) + Alp_bl(i)
     2065c avant        ALP(i) = alp_wake(i) + Alp_bl(i)
     2066              ALP(i) = alp_wake(i) + Alp_bl(i) + alp_offset ! modif sb
    20492067c         write(20,*)'ALE',ALE(i),Ale_bl(i),ale_wake(i)
    20502068c         write(21,*)'ALP',ALP(i),Alp_bl(i),alp_wake(i)
    … …  
    22892307
    22902308      endif
     2309c
     2310c===================================================================
     2311cJYG
     2312      IF (ip_ebil_phy.ge.2) THEN
     2313        ztit='after wake'
     2314        CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
     2315     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
     2316     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
     2317        call diagphy(airephy,ztit,ip_ebil_phy
     2318     e      , zero_v, zero_v, zero_v, zero_v, zero_v
     2319     e      , zero_v, zero_v, zero_v, ztsol
     2320     e      , d_h_vcol, d_qt, d_ec
     2321     s      , fs_bound, fq_bound )
     2322      END IF
     2323
    22912324c      print*,'apres callwake iflag_cldcon=', iflag_cldcon
    22922325c
    … …  
    24022435     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
    24032436     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
     2437        call diagphy(airephy,ztit,ip_ebil_phy
     2438     e      , zero_v, zero_v, zero_v, zero_v, zero_v
     2439     e      , zero_v, zero_v, zero_v, ztsol
     2440     e      , d_h_vcol, d_qt, d_ec
     2441     s      , fs_bound, fq_bound )
    24042442      END IF
    24052443
    … …  
    24502488         enddo
    24512489         tau_overturning_th(:)=zmax_th(:)/max(0.5*wmax_th(:),0.1)
    2452          print*,'TAU TH OK ',tau_overturning_th(1),detr_therm(1,3)
     2490         if(prt_level.ge.9)
     2491     &     write(lunout,*)'TAU TH OK ',
     2492     &     tau_overturning_th(1),detr_therm(1,3)
    24532493
    24542494c On impose que l'air autour de la fraction couverte par le thermique
    … …  
    25612601
    25622602!   les ratqs sont une combinaison de ratqss et ratqsc
    2563 !       write(lunout,*)'PHYLMD NOUVEAU TAU_RATQS ',tau_ratqs
     2603       if(prt_level.ge.9)
     2604     $       write(lunout,*)'PHYLMD NOUVEAU TAU_RATQS ',tau_ratqs
    25642605
    25652606         if (tau_ratqs>1.e-10) then
    … …  
    27852826     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
    27862827     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
     2828        call diagphy(airephy,ztit,ip_ebil_phy
     2829     e      , zero_v, zero_v, zero_v, zero_v, zero_v
     2830     e      , zero_v, zero_v, zero_v, ztsol
     2831     e      , d_h_vcol, d_qt, d_ec
     2832     s      , fs_bound, fq_bound )
    27872833      END IF
    27882834c
    … …  
    29923038      itaprad = itaprad + 1
    29933039
    2994       if (iflag_radia.eq.0) then
     3040      if (iflag_radia.eq.0 .and. prt_level.ge.9) then
    29953041      print *,'--------------------------------------------------'
    29963042      print *,'>>>> ATTENTION rayonnement desactive pour ce cas'
    … …  
    31953241     e      , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
    31963242     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
     3243         call diagphy(airephy,ztit,ip_ebil_phy
     3244     e      , zero_v, zero_v, zero_v, zero_v, zero_v
     3245     e      , zero_v, zero_v, zero_v, ztsol
     3246     e      , d_h_vcol, d_qt, d_ec
     3247     s      , fs_bound, fq_bound )
    31973248      END IF
    31983249c
    … …  
    32623313      IF (offline) THEN
    32633314
    3264          print*,'Attention on met a 0 les thermiques pour phystoke'
    3265          call phystokenc (
     3315       IF (prt_level.ge.9)
     3316     $    print*,'Attention on met a 0 les thermiques pour phystoke'
     3317         call phystokenc (
    32663318     I                   nlon,klev,pdtphys,rlon,rlat,
    32673319     I                   t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,

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

    r1277 r1322  
    1       Subroutine WAKE (p,ph,ppi,dtime,sigd_con
     1!
     2! $Id$
     3!
     4      Subroutine WAKE (p,ph,pi,dtime,sigd_con
    25     :                ,te0,qe0,omgb
    36     :                ,dtdwn,dqdwn,amdwn,amup,dta,dqa
    … …  
    2427c
    2528      use dimphy
     29      IMPLICIT none
     30c============================================================================
     31C
     32C
     33C   But : Decrire le comportement des poches froides apparaissant dans les
     34C        grands systemes convectifs, et fournir l'energie disponible pour
     35C        le declenchement de nouvelles colonnes convectives.
     36C
     37C   Variables d'etat : deltatw    : ecart de temperature wake-undisturbed area
     38C                      deltaqw    : ecart d'humidite wake-undisturbed area
     39C                      sigmaw     : fraction d'aire occupee par la poche.
     40C
     41C   Variable de sortie :
     42c
     43c                        wape : WAke Potential Energy
     44c                        fip  : Front Incident Power (W/m2) - ALP
     45c                        gfl  : Gust Front Length per unit area (m-1)
     46C                        dtls : large scale temperature tendency due to wake
     47C                        dqls : large scale humidity tendency due to wake
     48C                        hw   : hauteur de la poche
     49C                     dp_omgb : vertical gradient of large scale omega
     50C                     wdens   : densite de poches
     51C                      omgbdth: flux of Delta_Theta transported by LS omega
     52C                      dtKE   : differential heating (wake - unpertubed)
     53C                      dqKE   : differential moistening (wake - unpertubed)
     54C                      omg    : Delta_omg =vertical velocity diff. wake-undist. (Pa/s)
     55C                 dp_deltomg  : vertical gradient of omg (s-1)
     56C                     spread  : spreading term in dt_wake and dq_wake
     57C                 deltatw     : updated temperature difference (T_w-T_u).
     58C                 deltaqw     : updated humidity difference (q_w-q_u).
     59C                 sigmaw      : updated wake fractional area.
     60C                 d_deltat_gw : delta T tendency due to GW
     61c
     62C   Variables d'entree :
     63c
     64c                        aire : aire de la maille
     65c                        te0  : temperature dans l'environnement  (K)
     66C                        qe0  : humidite dans l'environnement     (kg/kg)
     67C                        omgb : vitesse verticale moyenne sur la maille (Pa/s)
     68C                        dtdwn: source de chaleur due aux descentes (K/s)
     69C                        dqdwn: source d'humidite due aux descentes (kg/kg/s)
     70C                        dta  : source de chaleur due courants satures et detrain  (K/s)
     71C                        dqa  : source d'humidite due aux courants satures et detra (kg/kg/s)
     72C                        amdwn: flux de masse total des descentes, par unite de
     73C                                surface de la maille (kg/m2/s)
     74C                        amup : flux de masse total des ascendances, par unite de
     75C                                surface de la maille (kg/m2/s)
     76C                        p    : pressions aux milieux des couches (Pa)
     77C                        ph   : pressions aux interfaces (Pa)
     78C                        pi  : (p/p_0)**kapa (adim)
     79C                        dtime: increment temporel (s)
     80c
     81C   Variables internes :
     82c
     83c                        rhow : masse volumique de la poche froide
     84C                        rho  : environment density at P levels
     85C                        rhoh : environment density at Ph levels
     86C                        te   : environment temperature | may change within
     87C                        qe   : environment humidity    | sub-time-stepping
     88C                        the  : environment potential temperature
     89C                        thu  : potential temperature in undisturbed area
     90C                        tu   :  temperature  in undisturbed area
     91C                        qu   : humidity in undisturbed area
     92C                      dp_omgb: vertical gradient og LS omega
     93C                      omgbw  : wake average vertical omega
     94C                     dp_omgbw: vertical gradient of omgbw
     95C                      omgbdq : flux of Delta_q transported by LS omega
     96C                        dth  : potential temperature diff. wake-undist.
     97C                        th1  : first pot. temp. for vertical advection (=thu)
     98C                        th2  : second pot. temp. for vertical advection (=thw)
     99C                        q1   : first humidity for vertical advection
     100C                        q2   : second humidity for vertical advection
     101C                     d_deltatw   : terme de redistribution pour deltatw
     102C                     d_deltaqw   : terme de redistribution pour deltaqw
     103C                      deltatw0   : deltatw initial
     104C                      deltaqw0   : deltaqw initial
     105C                      hw0    : hw initial
     106C                      sigmaw0: sigmaw initial
     107C                      amflux : horizontal mass flux through wake boundary
     108C                      wdens_ref: initial number of wakes per unit area (3D) or per
     109C                               unit length (2D), at the beginning of each time step
     110C                      Tgw    : 1 sur la période de onde de gravité
     111c                      Cgw    : vitesse de propagation de onde de gravité
     112c                      LL     : distance entre 2 poches
     113
     114c-------------------------------------------------------------------------
     115c          Déclaration de variables
     116c-------------------------------------------------------------------------
     117
     118#include "dimensions.h"
     119#include "YOMCST.h"
     120#include "cvthermo.h"
     121#include "iniprint.h"
     122
     123c Arguments en entree
     124c--------------------
     125
     126      REAL, dimension(klon,klev) :: p, pi
     127      REAL, dimension(klon,klev+1) :: ph, omgb
     128      REAL dtime
     129      REAL, dimension(klon,klev) :: te0,qe0
     130      REAL, dimension(klon,klev) :: dtdwn, dqdwn
     131      REAL, dimension(klon,klev) :: wdtPBL,wdqPBL
     132      REAL, dimension(klon,klev) :: udtPBL,udqPBL
     133      REAL, dimension(klon,klev) :: amdwn, amup
     134      REAL, dimension(klon,klev) :: dta, dqa
     135      REAL, dimension(klon) :: sigd_con
     136
     137c Sorties
     138c--------
     139
     140      REAL, dimension(klon,klev) :: deltatw, deltaqw, dth
     141      REAL, dimension(klon,klev) :: tu, qu
     142      REAL, dimension(klon,klev) :: dtls, dqls
     143      REAL, dimension(klon,klev) :: dtKE, dqKE
     144      REAL, dimension(klon,klev) :: dtPBL, dqPBL
     145      REAL, dimension(klon,klev) :: spread
     146      REAL, dimension(klon,klev) :: d_deltatgw
     147      REAL, dimension(klon,klev) :: d_deltatw2, d_deltaqw2
     148      REAL, dimension(klon,klev+1) :: omgbdth, omg
     149      REAL, dimension(klon,klev) :: dp_omgb, dp_deltomg
     150      REAL, dimension(klon,klev) :: d_deltat_gw
     151      REAL, dimension(klon) :: hw, sigmaw, wape, fip, gfl, Cstar
     152      REAL, dimension(klon) :: wdens
     153      INTEGER, dimension(klon) :: ktopw
     154
     155c Variables internes
     156c-------------------
     157
     158c Variables à fixer
     159      REAL ALON
     160      REAL coefgw
     161      REAL :: wdens_ref
     162      REAL stark
     163      REAL alpk
     164      REAL delta_t_min
     165      INTEGER nsub
     166      REAL dtimesub
     167      REAL sigmad, hwmin,wapecut
     168      REAL :: sigmaw_max
     169      REAL :: dens_rate
     170      REAL wdens0
     171cIM 080208
     172      LOGICAL, dimension(klon) :: gwake
     173
     174c Variables de sauvegarde
     175      REAL, dimension(klon,klev) :: deltatw0
     176      REAL, dimension(klon,klev) :: deltaqw0
     177      REAL, dimension(klon,klev) :: te, qe
     178      REAL, dimension(klon) :: sigmaw0, sigmaw1
     179
     180c Variables pour les GW
     181      REAL, DIMENSION(klon) :: LL
     182      REAL, dimension(klon,klev) :: N2
     183      REAL, dimension(klon,klev) :: Cgw
     184      REAL, dimension(klon,klev) :: Tgw
     185
     186c Variables liées au calcul de hw
     187      REAL, DIMENSION(klon) :: ptop_provis, ptop, ptop_new
     188      REAL, DIMENSION(klon) :: sum_dth
     189      REAL, DIMENSION(klon) :: dthmin
     190      REAL, DIMENSION(klon) :: z, dz, hw0
     191      INTEGER, DIMENSION(klon) :: ktop, kupper
     192
     193c Sub-timestep tendencies and related variables
     194       REAL d_deltatw(klon,klev),d_deltaqw(klon,klev)
     195       REAL d_te(klon,klev),d_qe(klon,klev)
     196       REAL d_sigmaw(klon),alpha(klon)
     197       REAL q0_min(klon),q1_min(klon)
     198       LOGICAL wk_adv(klon), OK_qx_qw(klon)
     199       REAL epsilon
     200       DATA epsilon/1.e-15/
     201
     202c Autres variables internes
     203      INTEGER isubstep, k, i
     204
     205      REAL, DIMENSION(klon) :: sum_thu, sum_tu, sum_qu,sum_thvu
     206      REAL, DIMENSION(klon) :: sum_dq, sum_rho
     207      REAL, DIMENSION(klon) :: sum_dtdwn, sum_dqdwn
     208      REAL, DIMENSION(klon) :: av_thu, av_tu, av_qu, av_thvu
     209      REAL, DIMENSION(klon) :: av_dth, av_dq, av_rho
     210      REAL, DIMENSION(klon) :: av_dtdwn, av_dqdwn
     211
     212      REAL, DIMENSION(klon,klev) :: rho, rhow
     213      REAL, DIMENSION(klon,klev+1) :: rhoh
     214      REAL, DIMENSION(klon,klev) :: rhow_moyen
     215      REAL, DIMENSION(klon,klev) :: zh
     216      REAL, DIMENSION(klon,klev+1) :: zhh
     217      REAL, DIMENSION(klon,klev) :: epaisseur1, epaisseur2
     218
     219      REAL, DIMENSION(klon,klev) :: the, thu
     220
     221!      REAL, DIMENSION(klon,klev) :: d_deltatw, d_deltaqw
     222
     223      REAL, DIMENSION(klon,klev+1) :: omgbw
     224      REAL, DIMENSION(klon) :: pupper
     225      REAL, DIMENSION(klon) :: omgtop
     226      REAL, DIMENSION(klon,klev) :: dp_omgbw
     227      REAL, DIMENSION(klon) :: ztop, dztop
     228      REAL, DIMENSION(klon,klev) :: alpha_up
     229     
     230      REAL, dimension(klon) :: RRe1, RRe2
     231      REAL :: RRd1, RRd2
     232      REAL, DIMENSION(klon,klev) :: Th1, Th2, q1, q2
     233      REAL, DIMENSION(klon,klev) :: D_Th1, D_Th2, D_dth
     234      REAL, DIMENSION(klon,klev) :: D_q1, D_q2, D_dq
     235      REAL, DIMENSION(klon,klev) :: omgbdq
     236
     237      REAL, dimension(klon) :: ff, gg
     238      REAL, dimension(klon) :: wape2, Cstar2, heff
     239
     240      REAL, DIMENSION(klon,klev) :: Crep
     241      REAL Crep_upper, Crep_sol
     242
     243      REAL, DIMENSION(klon,klev) :: ppi
     244
     245ccc nrlmd
     246      real, dimension(klon) :: death_rate,nat_rate
     247      real, dimension(klon,klev) :: entr
     248      real, dimension(klon,klev) :: detr
     249
     250C-------------------------------------------------------------------------
     251c         Initialisations
     252c-------------------------------------------------------------------------
     253
     254c      print*, 'wake initialisations'
     255
     256c   Essais d'initialisation avec sigmaw = 0.02 et hw = 10.
     257c-------------------------------------------------------------------------
     258
     259      DATA wapecut,sigmad, hwmin /5.,.02,10./
     260ccc nrlmd
     261      DATA sigmaw_max /0.4/
     262      DATA dens_rate /0.1/
     263ccc
     264C Longueur de maille (en m)
     265c-------------------------------------------------------------------------
     266
     267c      ALON = 3.e5
     268      ALON = 1.e6
     269
     270
     271C Configuration de coefgw,stark,wdens (22/02/06 by YU Jingmei)
     272c
     273c      coefgw : Coefficient pour les ondes de gravité
     274c       stark : Coefficient k dans Cstar=k*sqrt(2*WAPE)
     275c       wdens : Densité de poche froide par maille
     276c-------------------------------------------------------------------------
     277
     278ccc nrlmd      coefgw=10
     279c      coefgw=1
     280c      wdens0 = 1.0/(alon**2)
     281ccc nrlmd      wdens = 1.0/(alon**2)
     282ccc nrlmd      stark = 0.50
     283cCRtest
     284ccc nrlmd      alpk=0.1
     285c      alpk = 1.0
     286c      alpk = 0.5
     287c      alpk = 0.05
     288c
     289       stark  = 0.33
     290       Alpk   = 0.25
     291       wdens_ref  = 8.e-12
     292       coefgw = 4.
     293      Crep_upper=0.9
     294      Crep_sol=1.0
     295
     296ccc nrlmd Lecture du fichier wake_param.data
     297c$$$      OPEN(99,file='wake_param.data',status='old',
     298c$$$     $          form='formatted',err=9999)
     299c$$$      READ(99,*,end=9998) stark
     300c$$$      READ(99,*,end=9998) Alpk
     301c$$$      READ(99,*,end=9998) wdens_ref
     302c$$$      READ(99,*,end=9998) coefgw
     303c$$$9998  Continue
     304c$$$      CLOSE(99)
     305c$$$9999  Continue
     306c
     307c   Initialisation de toutes des densites a wdens_ref.
     308c   Les densites peuvent evoluer si les poches debordent
     309c   (voir au tout debut de la boucle sur les substeps)
     310      wdens = wdens_ref
     311c
     312c      print*,'stark',stark
     313c      print*,'alpk',alpk
     314c      print*,'wdens',wdens
     315c      print*,'coefgw',coefgw
     316ccc
     317C Minimum value for |T_wake - T_undist|. Used for wake top definition
     318c-------------------------------------------------------------------------
     319
     320      delta_t_min = 0.2
     321
     322C 1. - Save initial values and initialize tendencies
     323C --------------------------------------------------
     324
     325      DO k=1,klev
     326      DO i=1, klon
     327        ppi(i,k)=pi(i,k)
     328        deltatw0(i,k) = deltatw(i,k)
     329        deltaqw0(i,k)= deltaqw(i,k)
     330        te(i,k) = te0(i,k)
     331        qe(i,k) = qe0(i,k)
     332        dtls(i,k) = 0.
     333        dqls(i,k) = 0.
     334        d_deltat_gw(i,k)=0.
     335        d_te(i,k) = 0.
     336        d_qe(i,k) = 0.
     337        d_deltatw(i,k) = 0.
     338        d_deltaqw(i,k) = 0.
     339!IM 060508 beg
     340        d_deltatw2(i,k)=0.
     341        d_deltaqw2(i,k)=0.
     342!IM 060508 end
     343      ENDDO
     344      ENDDO
     345c      sigmaw1=sigmaw
     346c      IF (sigd_con.GT.sigmaw1) THEN
     347c      print*, 'sigmaw,sigd_con', sigmaw, sigd_con
     348c      ENDIF
     349      DO i=1, klon
     350cc      sigmaw(i) = amax1(sigmaw(i),sigd_con(i))
     351      sigmaw(i) = amax1(sigmaw(i),sigmad)
     352      sigmaw(i) = amin1(sigmaw(i),0.99)
     353      sigmaw0(i) = sigmaw(i)
     354      wape(i) = 0.
     355      wape2(i) = 0.
     356      d_sigmaw(i) = 0.
     357      ktopw(i) = 0
     358      ENDDO
     359C
     360C
     361C 2. - Prognostic part
     362C --------------------
     363C
     364C
     365C 2.1 - Undisturbed area and Wake integrals
     366C ---------------------------------------------------------
     367
     368      DO i=1, klon
     369      z(i) = 0.
     370      ktop(i)=0
     371      kupper(i) = 0
     372      sum_thu(i) = 0.
     373      sum_tu(i) = 0.
     374      sum_qu(i) = 0.
     375      sum_thvu(i) = 0.
     376      sum_dth(i) = 0.
     377      sum_dq(i) = 0.
     378      sum_rho(i) = 0.
     379      sum_dtdwn(i) = 0.
     380      sum_dqdwn(i) = 0.
     381
     382      av_thu(i) = 0.
     383      av_tu(i) =0.
     384      av_qu(i) =0.
     385      av_thvu(i) = 0.
     386      av_dth(i) = 0.
     387      av_dq(i) = 0.
     388      av_rho(i) =0.
     389      av_dtdwn(i) =0.
     390      av_dqdwn(i) = 0.
     391      ENDDO
     392c
     393c Distance between wakes
     394       DO i = 1,klon
     395        LL(i) = (1-sqrt(sigmaw(i)))/sqrt(wdens(i))
     396       ENDDO
     397C Potential temperatures and humidity
     398c----------------------------------------------------------
     399      DO k =1,klev
     400       DO i=1, klon
     401!        write(*,*)'wake 1',i,k,rd,te(i,k)
     402        rho(i,k) = p(i,k)/(rd*te(i,k))
     403!        write(*,*)'wake 2',rho(i,k)
     404        IF(k .eq. 1) THEN
     405!        write(*,*)'wake 3',i,k,rd,te(i,k)
     406          rhoh(i,k) = ph(i,k)/(rd*te(i,k))
     407!        write(*,*)'wake 4',i,k,rd,te(i,k)
     408          zhh(i,k)=0
     409        ELSE
     410!          write(*,*)'wake 5',rd,(te(i,k)+te(i,k-1))
     411          rhoh(i,k) = ph(i,k)*2./(rd*(te(i,k)+te(i,k-1)))
     412!          write(*,*)'wake 6',(-rhoh(i,k)*RG)+zhh(i,k-1)
     413          zhh(i,k)=(ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*RG)+zhh(i,k-1)
     414        ENDIF
     415!          write(*,*)'wake 7',ppi(i,k)
     416        the(i,k) = te(i,k)/ppi(i,k)
     417        thu(i,k) = (te(i,k) - deltatw(i,k)*sigmaw(i))/ppi(i,k)
     418        tu(i,k) = te(i,k) - deltatw(i,k)*sigmaw(i)
     419        qu(i,k)  =  qe(i,k) - deltaqw(i,k)*sigmaw(i)
     420!          write(*,*)'wake 8',(rd*(te(i,k)+deltatw(i,k)))
     421        rhow(i,k) = p(i,k)/(rd*(te(i,k)+deltatw(i,k)))
     422        dth(i,k) = deltatw(i,k)/ppi(i,k)
     423       ENDDO
     424      ENDDO
     425       
     426      DO k = 1, klev-1
     427      DO i=1, klon
     428        IF(k.eq.1) THEN
     429          N2(i,k)=0
     430        ELSE
     431          N2(i,k)=amax1(0.,-RG**2/the(i,k)*rho(i,k)*(the(i,k+1)-
     432     $            the(i,k-1))/(p(i,k+1)-p(i,k-1)))
     433        ENDIF
     434        ZH(i,k)=(zhh(i,k)+zhh(i,k+1))/2
     435
     436        Cgw(i,k)=sqrt(N2(i,k))*ZH(i,k)
     437        Tgw(i,k)=coefgw*Cgw(i,k)/LL(i)
     438      ENDDO
     439      ENDDO
     440
     441      DO i=1, klon
     442      N2(i,klev)=0
     443      ZH(i,klev)=0
     444      Cgw(i,klev)=0
     445      Tgw(i,klev)=0
     446      ENDDO
     447
     448c  Calcul de la masse volumique moyenne de la colonne   (bdlmd)
     449c-----------------------------------------------------------------
     450
     451      DO k=1,klev
     452       DO i=1, klon
     453        epaisseur1(i,k)=0.
     454        epaisseur2(i,k)=0.
     455       ENDDO
     456      ENDDO
     457
     458      DO i=1, klon
     459      epaisseur1(i,1)= -(ph(i,2)-ph(i,1))/(rho(i,1)*rg)+1.
     460      epaisseur2(i,1)= -(ph(i,2)-ph(i,1))/(rho(i,1)*rg)+1.
     461      rhow_moyen(i,1) = rhow(i,1)
     462      ENDDO
     463
     464      DO k = 2, klev
     465      DO i=1, klon
     466        epaisseur1(i,k)= -(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg) +1.
     467        epaisseur2(i,k)=epaisseur2(i,k-1)+epaisseur1(i,k)
     468        rhow_moyen(i,k) = (rhow_moyen(i,k-1)*epaisseur2(i,k-1)+
     469     $                 rhow(i,k)*epaisseur1(i,k))/epaisseur2(i,k)
     470      ENDDO
     471      ENDDO
     472
     473C
     474C Choose an integration bound well above wake top
     475c-----------------------------------------------------------------
     476c
     477C       Pupper = 50000.  ! melting level
     478c       Pupper = 60000.
     479c       Pupper = 80000.  ! essais pour case_e
     480       DO i = 1,klon
     481        Pupper(i) = 0.6*ph(i,1)
     482        Pupper(i) = max(Pupper(i), 45000.)
     483ccc        Pupper(i) = 60000.
     484       ENDDO
     485
     486C
     487C    Determine Wake top pressure (Ptop) from buoyancy integral
     488C    --------------------------------------------------------
     489c
     490c-1/ Pressure of the level where dth becomes less than delta_t_min.
     491
     492      DO i=1,klon
     493      ptop_provis(i)=ph(i,1)
     494      ENDDO
     495      DO k= 2,klev
     496      DO i=1,klon
     497c
     498cIM v3JYG; ptop_provis(i).LT. ph(i,1)
     499c
     500        IF (dth(i,k) .GT. -delta_t_min .and.
     501     $      dth(i,k-1).LT. -delta_t_min .and.
     502     $      ptop_provis(i).EQ. ph(i,1)) THEN
     503          ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)
     504     $          - (dth(i,k-1)+delta_t_min)*p(i,k)) /
     505     $          (dth(i,k) - dth(i,k-1))
     506        ENDIF
     507      ENDDO
     508      ENDDO
     509
     510c-2/ dth integral
     511
     512      DO i=1,klon
     513      sum_dth(i) = 0.
     514      dthmin(i) = -delta_t_min
     515      z(i) = 0.
     516      ENDDO
     517
     518      DO k = 1,klev
     519      DO i=1,klon
     520        dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-Ph(i,k))/(rho(i,k)*rg)
     521        IF (dz(i) .gt. 0) THEN
     522          z(i) = z(i)+dz(i)
     523          sum_dth(i) = sum_dth(i) + dth(i,k)*dz(i)
     524          dthmin(i) = amin1(dthmin(i),dth(i,k))
     525        ENDIF
     526      ENDDO
     527      ENDDO
     528
     529c-3/ height of triangle with area= sum_dth and base = dthmin
     530
     531      DO i=1,klon
     532      hw0(i) = 2.*sum_dth(i)/amin1(dthmin(i),-0.5)
     533      hw0(i) = amax1(hwmin,hw0(i))
     534      ENDDO
     535
     536c-4/ now, get Ptop
     537
     538      DO i=1,klon
     539      z(i) = 0.
     540      ptop(i) = ph(i,1)
     541      ENDDO
     542
     543      DO k = 1,klev
     544      DO i=1,klon
     545        dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg),hw0(i)-z(i))
     546        IF (dz(i) .gt. 0) THEN
     547         z(i) = z(i)+dz(i)
     548         ptop(i) = ph(i,k)-rho(i,k)*rg*dz(i)
     549        ENDIF
     550      ENDDO
     551      ENDDO
     552
     553
     554C-5/ Determination de ktop et kupper
     555
     556      DO k=klev,1,-1
     557      DO i=1,klon
     558        IF (ph(i,k+1) .lt. ptop(i)) ktop(i)=k
     559        IF (ph(i,k+1) .lt. pupper(i)) kupper(i)=k
     560      ENDDO
     561      ENDDO
     562
     563c      On evite kupper = 1
     564      DO i=1,klon
     565        kupper(i) = max(kupper(i),2)
     566      ENDDO
     567
     568
     569c-6/ Correct ktop and ptop
     570
     571      DO i = 1,klon
     572        ptop_new(i)=ptop(i)
     573      ENDDO
     574      DO k= klev,2,-1
     575      DO i=1,klon
     576        IF (k .LE. ktop(i) .and.
     577     $      ptop_new(i) .EQ. ptop(i) .and.
     578     $      dth(i,k) .GT. -delta_t_min .and.
     579     $      dth(i,k-1).LT. -delta_t_min) THEN
     580          ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)
     581     $          - (dth(i,k-1)+delta_t_min)*p(i,k)) /
     582     $          (dth(i,k) - dth(i,k-1))
     583        ENDIF
     584      ENDDO
     585      ENDDO
     586
     587      DO i=1,klon
     588        ptop(i) = ptop_new(i)
     589      ENDDO
     590
     591      DO k=klev,1,-1
     592      DO i=1,klon
     593        IF (ph(i,k+1) .lt. ptop(i)) ktop(i)=k
     594      ENDDO
     595      ENDDO
     596c
     597c-5/ Set deltatw & deltaqw to 0 above kupper
     598c
     599      DO k = 1,klev
     600      DO i=1,klon
     601       IF (k.GE. kupper(i)) THEN
     602        deltatw(i,k) = 0.
     603        deltaqw(i,k) = 0.
     604       ENDIF
     605      ENDDO
     606      ENDDO
     607c
     608C
     609C Vertical gradient of LS omega
     610C
     611      DO k = 1,klev
     612      DO i=1,klon
     613       IF (k.LE. kupper(i)) THEN
     614        dp_omgb(i,k) = (omgb(i,k+1) - omgb(i,k))/(ph(i,k+1)-ph(i,k))
     615       ENDIF
     616      ENDDO
     617      ENDDO
     618C
     619C Integrals (and wake top level number)
     620C --------------------------------------
     621C
     622C Initialize sum_thvu to 1st level virt. pot. temp.
     623
     624      DO i=1,klon
     625      z(i) = 1.
     626      dz(i) = 1.
     627      sum_thvu(i) =  thu(i,1)*(1.+eps*qu(i,1))*dz(i)
     628      sum_dth(i) = 0.
     629      ENDDO
     630
     631      DO k = 1,klev
     632      DO i=1,klon
     633        dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
     634        IF (dz(i) .GT. 0) THEN
     635         z(i) = z(i)+dz(i)
     636         sum_thu(i) = sum_thu(i) + thu(i,k)*dz(i)
     637         sum_tu(i) = sum_tu(i) + tu(i,k)*dz(i)
     638         sum_qu(i) = sum_qu(i) + qu(i,k)*dz(i)
     639         sum_thvu(i) = sum_thvu(i) + thu(i,k)*(1.+eps*qu(i,k))*dz(i)
     640         sum_dth(i) = sum_dth(i) + dth(i,k)*dz(i)
     641         sum_dq(i) = sum_dq(i) + deltaqw(i,k)*dz(i)
     642         sum_rho(i) = sum_rho(i) + rhow(i,k)*dz(i)
     643         sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i,k)*dz(i)
     644         sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i,k)*dz(i)
     645        ENDIF
     646      ENDDO
     647      ENDDO
     648c
     649      DO i=1,klon
     650        hw0(i) = z(i)
     651      ENDDO
     652c
     653C
     654C 2.1 - WAPE and mean forcing computation
     655C ---------------------------------------
     656C
     657C ---------------------------------------
     658C
     659C Means
     660
     661      DO i=1,klon
     662      av_thu(i) = sum_thu(i)/hw0(i)
     663      av_tu(i) = sum_tu(i)/hw0(i)
     664      av_qu(i) = sum_qu(i)/hw0(i)
     665      av_thvu(i) = sum_thvu(i)/hw0(i)
     666c      av_thve = sum_thve/hw0
     667      av_dth(i) = sum_dth(i)/hw0(i)
     668      av_dq(i) = sum_dq(i)/hw0(i)
     669      av_rho(i) = sum_rho(i)/hw0(i)
     670      av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
     671      av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
     672
     673      wape(i) = - rg*hw0(i)*(av_dth(i)
     674     $     + eps*(av_thu(i)*av_dq(i)+av_dth(i)*av_qu(i)+av_dth(i)*
     675     $     av_dq(i) ))/av_thvu(i)
     676      ENDDO
     677C
     678C 2.2 Prognostic variable update
     679C ------------------------------
     680C
     681C Filter out bad wakes
     682
     683      DO k = 1,klev
     684       DO i=1,klon
     685        IF ( wape(i) .LT. 0.) THEN
     686          deltatw(i,k) = 0.
     687          deltaqw(i,k) = 0.
     688          dth(i,k) = 0.
     689        ENDIF
     690       ENDDO
     691      ENDDO
     692c
     693      DO i=1,klon
     694      IF ( wape(i) .LT. 0.) THEN
     695        wape(i) = 0.
     696        Cstar(i) = 0.
     697        hw(i) = hwmin
     698        sigmaw(i) = amax1(sigmad,sigd_con(i))
     699        fip(i) = 0.
     700        gwake(i) = .FALSE.
     701      ELSE
     702        Cstar(i) = stark*sqrt(2.*wape(i))
     703        gwake(i) = .TRUE.
     704      ENDIF
     705      ENDDO
     706
     707c
     708c Check qx and qw positivity
     709c --------------------------
     710      DO i = 1,klon
     711       q0_min(i)=min(  (qe(i,1)-sigmaw(i)*deltaqw(i,1)),
     712     $              (qe(i,1)+(1.-sigmaw(i))*deltaqw(i,1))  )
     713      ENDDO
     714      DO k = 2,klev
     715      DO i = 1,klon
     716        q1_min(i)=min(  (qe(i,k)-sigmaw(i)*deltaqw(i,k)),
     717     $              (qe(i,k)+(1.-sigmaw(i))*deltaqw(i,k))  )
     718        IF (q1_min(i).le.q0_min(i)) THEN
     719          q0_min(i)=q1_min(i)
     720        ENDIF
     721      ENDDO
     722      ENDDO
     723c
     724      DO i = 1,klon
     725       OK_qx_qw(i) = q0_min(i) .GE. 0.
     726       alpha(i) = 1.
     727      ENDDO
     728c
     729CC -----------------------------------------------------------------
     730C    Sub-time-stepping
     731C    -----------------
     732C
     733      nsub=10
     734      dtimesub=dtime/nsub
     735c
     736c------------------------------------------------------------
     737      DO isubstep = 1,nsub
     738c------------------------------------------------------------
     739c
     740c wk_adv is the logical flag enabling wake evolution in the time advance loop
     741      DO i = 1,klon
     742       wk_adv(i) = OK_qx_qw(i) .AND. alpha(i) .GE. 1.
     743      ENDDO
     744c
     745ccc nrlmd   Ajout d'un recalcul de wdens dans le cas d'un entrainement négatif de ktop à kupper --------
     746ccc           On calcule pour cela une densité wdens0 pour laquelle on aurait un entrainement nul ---
     747      DO i=1,klon
     748cc       print *,' isubstep,wk_adv(i),cstar(i),wape(i) ',
     749cc     $           isubstep,wk_adv(i),cstar(i),wape(i)
     750        IF (wk_adv(i) .AND. cstar(i).GT.0.01) THEN
     751           omg(i,kupper(i)+1) = - Rg*amdwn(i,kupper(i)+1)/sigmaw(i)
     752     $                + Rg*amup(i,kupper(i)+1)/(1.-sigmaw(i))
     753           wdens0 = ( sigmaw(i) / (4.*3.14) ) *
     754     $     ( (1.-sigmaw(i)) * omg(i,kupper(i)+1) /
     755     $     ( (ph(i,1)-pupper(i)) * cstar(i) )  ) **(2)
     756         IF ( wdens(i) .LE. wdens0*1.1 ) THEN
     757            wdens(i) = wdens0
     758         ENDIF
     759cc         print*,'omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i)
     760cc     $     ,ph(i,1)-pupper(i)',
     761cc     $             omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i)
     762cc     $     ,ph(i,1)-pupper(i)
     763        ENDIF
     764      ENDDO
     765
     766ccc nrlmd
     767
     768      DO i=1,klon
     769       IF (wk_adv(i)) THEN
     770        gfl(i) = 2.*sqrt(3.14*wdens(i)*sigmaw(i))
     771        sigmaw(i)=amin1(sigmaw(i),sigmaw_max)
     772       ENDIF
     773      ENDDO
     774      DO i=1,klon
     775        IF (wk_adv(i)) THEN
     776ccc nrlmd          Introduction du taux de mortalité des poches et test sur sigmaw_max=0.4
     777ccc         d_sigmaw(i) = gfl(i)*Cstar(i)*dtimesub
     778           IF (sigmaw(i).ge.sigmaw_max) THEN
     779           death_rate(i)=gfl(i)*Cstar(i)/sigmaw(i)
     780           ELSE
     781             death_rate(i)=0.
     782           END IF
     783        d_sigmaw(i) = gfl(i)*Cstar(i)*dtimesub
     784     $               - death_rate(i)*sigmaw(i)*dtimesub
     785c     $              - nat_rate(i)*sigmaw(i)*dtimesub
     786cc        print*, 'd_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i),
     787cc     $  death_rate(i),ktop(i),kupper(i)',
     788cc     $                 d_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i),
     789cc     $  death_rate(i),ktop(i),kupper(i)
     790
     791c        sigmaw(i) =sigmaw(i) + gfl(i)*Cstar(i)*dtimesub
     792c        sigmaw(i) =min(sigmaw(i),0.99)     !!!!!!!!
     793c        wdens = wdens0/(10.*sigmaw)
     794c        sigmaw =max(sigmaw,sigd_con)
     795c        sigmaw =max(sigmaw,sigmad)
     796        ENDIF
     797      ENDDO
     798C
     799C
     800c calcul de la difference de vitesse verticale poche - zone non perturbee
     801cIM 060208 differences par rapport au code initial; init. a 0 dp_deltomg
     802cIM 060208 et omg sur les niveaux de 1 a klev+1, alors que avant l'on definit
     803cIM 060208 au niveau k=1..?
     804      DO k= 1,klev
     805      DO i = 1,klon
     806      if (wk_adv(i)) THEN !!! nrlmd
     807        dp_deltomg(i,k)=0.
     808      end if
     809      ENDDO
     810      ENDDO
     811      DO k= 1,klev+1
     812      DO i = 1,klon
     813      if (wk_adv(i)) THEN !!! nrlmd
     814        omg(i,k)=0.
     815      end if
     816      ENDDO
     817      ENDDO
     818c
     819      DO i=1,klon
     820        IF (wk_adv(i)) THEN
     821        z(i)= 0.
     822        omg(i,1) = 0.
     823        dp_deltomg(i,1) = -(gfl(i)*Cstar(i))/(sigmaw(i) * (1-sigmaw(i)))
     824        ENDIF
     825      ENDDO
     826c
     827      DO k= 2,klev
     828      DO i = 1,klon
     829       IF( wk_adv(i) .AND. k .LE. ktop(i)) THEN
     830          dz(i) = -(ph(i,k)-ph(i,k-1))/(rho(i,k-1)*rg)
     831          z(i) = z(i)+dz(i)
     832          dp_deltomg(i,k)= dp_deltomg(i,1)
     833          omg(i,k)= dp_deltomg(i,1)*z(i)
     834       ENDIF
     835      ENDDO
     836      ENDDO
     837c
     838      DO i = 1,klon
     839        IF (wk_adv(i)) THEN
     840        dztop(i)=-(ptop(i)-ph(i,ktop(i)))/(rho(i,ktop(i))*rg)
     841        ztop(i) = z(i)+dztop(i)
     842        omgtop(i)=dp_deltomg(i,1)*ztop(i)
     843        ENDIF
     844      ENDDO
     845c
     846c        -----------------
     847c        From m/s to Pa/s
     848c        -----------------
     849c
     850       DO i=1,klon
     851        IF (wk_adv(i)) THEN
     852        omgtop(i) = -rho(i,ktop(i))*rg*omgtop(i)
     853        dp_deltomg(i,1) = omgtop(i)/(ptop(i)-ph(i,1))
     854        ENDIF
     855       ENDDO
     856c
     857      DO k= 1,klev
     858      DO i = 1,klon
     859       IF( wk_adv(i) .AND. k .LE. ktop(i)) THEN
     860          omg(i,k) = - rho(i,k)*rg*omg(i,k)
     861          dp_deltomg(i,k) = dp_deltomg(i,1)
     862       ENDIF
     863      ENDDO
     864      ENDDO
     865c
     866c   raccordement lineaire de omg de ptop a pupper
     867
     868      DO i=1,klon
     869      IF ( wk_adv(i) .AND. kupper(i) .GT. ktop(i)) THEN
     870        omg(i,kupper(i)+1) = - Rg*amdwn(i,kupper(i)+1)/sigmaw(i)
     871     $                + Rg*amup(i,kupper(i)+1)/(1.-sigmaw(i))
     872        dp_deltomg(i,kupper(i)) = (omgtop(i)-omg(i,kupper(i)+1))/
     873     $                     (ptop(i)-pupper(i))
     874      ENDIF
     875      ENDDO
     876c
     877cc      DO i=1,klon
     878cc        print*,'Pente entre 0 et kupper (référence)'
     879cc     $        ,omg(i,kupper(i)+1)/(pupper(i)-ph(i,1))
     880cc        print*,'Pente entre ktop et kupper'
     881cc     $        ,(omg(i,kupper(i)+1)-omgtop(i))/(pupper(i)-ptop(i))
     882cc      ENDDO
     883cc
     884      DO k= 1,klev
     885      DO i = 1,klon
     886       IF( wk_adv(i) .AND. k .GT. ktop(i) .AND. k .LE. kupper(i)) THEN
     887          dp_deltomg(i,k) = dp_deltomg(i,kupper(i))
     888          omg(i,k) = omgtop(i)+(ph(i,k)-ptop(i))*dp_deltomg(i,kupper(i))
     889       ENDIF
     890      ENDDO
     891      ENDDO
     892ccc nrlmd
     893cc      DO i=1,klon
     894cc      print*,'deltaw_ktop,deltaw_conv',omgtop(i),omg(i,kupper(i)+1)
     895cc      END DO
     896ccc
     897c
     898c
     899c--    Compute wake average vertical velocity omgbw
     900c
     901c
     902      DO k = 1,klev+1
     903      DO i=1,klon
     904        IF ( wk_adv(i)) THEN
     905        omgbw(i,k) = omgb(i,k)+(1.-sigmaw(i))*omg(i,k)
     906        ENDIF
     907      ENDDO
     908      ENDDO
     909c--    and its vertical gradient dp_omgbw
     910c
     911      DO k = 1,klev
     912      DO i=1,klon
     913        IF ( wk_adv(i)) THEN
     914        dp_omgbw(i,k) = (omgbw(i,k+1)-omgbw(i,k))/(ph(i,k+1)-ph(i,k))
     915        ENDIF
     916      ENDDO
     917      ENDDO
     918C
     919c--    Upstream coefficients for omgb velocity
     920c--    (alpha_up(k) is the coefficient of the value at level k)
     921c--    (1-alpha_up(k) is the coefficient of the value at level k-1)
     922      DO k = 1,klev
     923      DO i=1,klon
     924        IF ( wk_adv(i)) THEN
     925         alpha_up(i,k) = 0.
     926         IF (omgb(i,k) .GT. 0.) alpha_up(i,k) = 1.
     927        ENDIF
     928      ENDDO
     929      ENDDO
     930
     931c  Matrix expressing [The,deltatw] from  [Th1,Th2]
     932
     933      DO i=1,klon
     934        IF ( wk_adv(i)) THEN
     935         RRe1(i) = 1.-sigmaw(i)
     936         RRe2(i) = sigmaw(i)
     937        ENDIF
     938      ENDDO
     939      RRd1 = -1.
     940      RRd2 = 1.
     941c
     942c--    Get [Th1,Th2], dth and [q1,q2]
     943c
     944      DO k= 1,klev
     945      DO i = 1,klon
     946       IF( wk_adv(i) .AND. k .LE. kupper(i)+1) THEN
     947        dth(i,k) = deltatw(i,k)/ppi(i,k)
     948        Th1(i,k) = the(i,k) - sigmaw(i)     *dth(i,k)   ! undisturbed area
     949        Th2(i,k) = the(i,k) + (1.-sigmaw(i))*dth(i,k)   ! wake
     950        q1(i,k) = qe(i,k) - sigmaw(i)     *deltaqw(i,k) ! undisturbed area
     951        q2(i,k) = qe(i,k) + (1.-sigmaw(i))*deltaqw(i,k) ! wake
     952       ENDIF
     953      ENDDO
     954      ENDDO
     955
     956      DO i=1,klon
     957       if (wk_adv(i)) then !!! nrlmd
     958       D_Th1(i,1) = 0.
     959       D_Th2(i,1) = 0.
     960       D_dth(i,1) = 0.
     961       D_q1(i,1) = 0.
     962       D_q2(i,1) = 0.
     963       D_dq(i,1) = 0.
     964       end if
     965      ENDDO
     966
     967      DO k= 2,klev
     968      DO i = 1,klon
     969       IF( wk_adv(i) .AND. k .LE. kupper(i)+1) THEN
     970        D_Th1(i,k) = Th1(i,k-1)-Th1(i,k)
     971        D_Th2(i,k) = Th2(i,k-1)-Th2(i,k)
     972        D_dth(i,k) = dth(i,k-1)-dth(i,k)
     973        D_q1(i,k) = q1(i,k-1)-q1(i,k)
     974        D_q2(i,k) = q2(i,k-1)-q2(i,k)
     975        D_dq(i,k) = deltaqw(i,k-1)-deltaqw(i,k)
     976       ENDIF
     977      ENDDO
     978      ENDDO
     979
     980      DO i=1,klon
     981        IF( wk_adv(i)) THEN
     982         omgbdth(i,1) = 0.
     983         omgbdq(i,1) = 0.
     984        ENDIF
     985      ENDDO
     986
     987      DO k= 2,klev
     988      DO i = 1,klon
     989       IF( wk_adv(i) .AND. k .LE. kupper(i)+1) THEN  !   loop on interfaces
     990        omgbdth(i,k) = omgb(i,k)*(    dth(i,k-1) -     dth(i,k))
     991        omgbdq(i,k)  = omgb(i,k)*(deltaqw(i,k-1) - deltaqw(i,k))
     992       ENDIF
     993      ENDDO
     994      ENDDO
     995c
     996c-----------------------------------------------------------------
     997      DO k= 1,klev
     998      DO i = 1,klon
     999       IF( wk_adv(i) .AND. k .LE. kupper(i)-1) THEN
     1000c-----------------------------------------------------------------
     1001c
     1002c   Compute redistribution (advective) term
     1003c
     1004         d_deltatw(i,k) =
     1005     $             dtimesub/(Ph(i,k)-Ph(i,k+1))*(
     1006     $       RRd1*omg(i,k  )*sigmaw(i)     *D_Th1(i,k)
     1007     $      -RRd2*omg(i,k+1)*(1.-sigmaw(i))*D_Th2(i,k+1)
     1008     $      -(1.-alpha_up(i,k))*omgbdth(i,k) - alpha_up(i,k+1)*
     1009     $      omgbdth(i,k+1))*ppi(i,k)
     1010c         print*,'d_deltatw=',d_deltatw(i,k)
     1011c
     1012         d_deltaqw(i,k) =
     1013     $             dtimesub/(Ph(i,k)-Ph(i,k+1))*(
     1014     $       RRd1*omg(i,k  )*sigmaw(i)     *D_q1(i,k)
     1015     $      -RRd2*omg(i,k+1)*(1.-sigmaw(i))*D_q2(i,k+1)
     1016     $      -(1.-alpha_up(i,k))*omgbdq(i,k) - alpha_up(i,k+1)*
     1017     $      omgbdq(i,k+1))
     1018c         print*,'d_deltaqw=',d_deltaqw(i,k)
     1019c
     1020c   and increment large scale tendencies
     1021c
     1022
     1023c
     1024C
     1025CC -----------------------------------------------------------------
     1026         d_te(i,k) =  dtimesub*(
     1027     $        ( RRe1(i)*omg(i,k  )*sigmaw(i)     *D_Th1(i,k)
     1028     $         -RRe2(i)*omg(i,k+1)*(1.-sigmaw(i))*D_Th2(i,k+1) )
     1029     $               /(Ph(i,k)-Ph(i,k+1))
     1030ccc nrlmd     $         -sigmaw(i)*(1.-sigmaw(i))*dth(i,k)*dp_deltomg(i,k)
     1031     $         -sigmaw(i)*(1.-sigmaw(i))*dth(i,k)
     1032     $         *(omg(i,k)-omg(i,k+1))/(Ph(i,k)-Ph(i,k+1))
     1033ccc
     1034     $                      )*ppi(i,k)
     1035c
     1036         d_qe(i,k) =  dtimesub*(
     1037     $        ( RRe1(i)*omg(i,k  )*sigmaw(i)     *D_q1(i,k)
     1038     $         -RRe2(i)*omg(i,k+1)*(1.-sigmaw(i))*D_q2(i,k+1) )
     1039     $               /(Ph(i,k)-Ph(i,k+1))
     1040ccc nrlmd     $         -sigmaw(i)*(1.-sigmaw(i))*deltaqw(i,k)*dp_deltomg(i,k)
     1041     $         -sigmaw(i)*(1.-sigmaw(i))*deltaqw(i,k)
     1042     $         *(omg(i,k)-omg(i,k+1))/(Ph(i,k)-Ph(i,k+1))
     1043ccc
     1044     $                      )
     1045ccc nrlmd
     1046       ELSE IF(wk_adv(i) .AND. k .EQ. kupper(i)) THEN
     1047         d_te(i,k) =  dtimesub*(
     1048     $       ( RRe1(i)*omg(i,k  )*sigmaw(i)     *D_Th1(i,k)
     1049     $        /(Ph(i,k)-Ph(i,k+1)))
     1050     $                       )*ppi(i,k)
     1051
     1052         d_qe(i,k) =  dtimesub*(
     1053     $       ( RRe1(i)*omg(i,k  )*sigmaw(i)     *D_q1(i,k)       
     1054     $        /(Ph(i,k)-Ph(i,k+1)))
     1055     $                       )
     1056
     1057       ENDIF
     1058ccc
     1059      ENDDO
     1060      ENDDO
     1061c------------------------------------------------------------------
     1062C
     1063C   Increment state variables
     1064
     1065      DO k= 1,klev
     1066      DO i = 1,klon
     1067ccc nrlmd       IF( wk_adv(i) .AND. k .LE. kupper(i)-1) THEN
     1068        IF( wk_adv(i) .AND. k .LE. kupper(i)) THEN
     1069ccc
     1070
     1071
     1072c
     1073c Coefficient de répartition
     1074
     1075        Crep(i,k)=Crep_sol*(ph(i,kupper(i))-ph(i,k))/(ph(i,kupper(i))
     1076     $          -ph(i,1))
     1077        Crep(i,k)=Crep(i,k)+Crep_upper*(ph(i,1)-ph(i,k))/(p(i,1)-
     1078     $          ph(i,kupper(i)))
     1079       
     1080
     1081c Reintroduce compensating subsidence term.
     1082
     1083c        dtKE(k)=(dtdwn(k)*Crep(k))/sigmaw
     1084c        dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k))
     1085c     .                   /(1-sigmaw)
     1086c        dqKE(k)=(dqdwn(k)*Crep(k))/sigmaw
     1087c        dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k))
     1088c     .                   /(1-sigmaw)
     1089c
     1090c        dtKE(k)=(dtdwn(k)*Crep(k)+(1-Crep(k))*dta(k))/sigmaw
     1091c        dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k)*Crep(k))
     1092c     .                   /(1-sigmaw)
     1093c        dqKE(k)=(dqdwn(k)*Crep(k)+(1-Crep(k))*dqa(k))/sigmaw
     1094c        dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k)*Crep(k))
     1095c     .                   /(1-sigmaw)
     1096
     1097        dtKE(i,k)=(dtdwn(i,k)/sigmaw(i) - dta(i,k)/(1.-sigmaw(i)))
     1098        dqKE(i,k)=(dqdwn(i,k)/sigmaw(i) - dqa(i,k)/(1.-sigmaw(i)))
     1099c        print*,'dtKE= ',dtKE(i,k),' dqKE= ',dqKE(i,k)
     1100c
     1101        dtPBL(i,k)=(wdtPBL(i,k)/sigmaw(i) - udtPBL(i,k)/(1.-sigmaw(i)))
     1102        dqPBL(i,k)=(wdqPBL(i,k)/sigmaw(i) - udqPBL(i,k)/(1.-sigmaw(i)))
     1103c        print*,'dtPBL= ',dtPBL(i,k),' dqPBL= ',dqPBL(i,k)
     1104c
     1105ccc nrlmd          Prise en compte du taux de mortalité
     1106ccc               Définitions de entr, detr
     1107        detr(i,k)=0.
     1108
     1109        entr(i,k)=detr(i,k)+gfl(i)*cstar(i)+
     1110     $          sigmaw(i)*(1.-sigmaw(i))*dp_deltomg(i,k)
     1111
     1112        spread(i,k) = (entr(i,k)-detr(i,k))/sigmaw(i)
     1113ccc        spread(i,k) = (1.-sigmaw(i))*dp_deltomg(i,k)+gfl(i)*Cstar(i)/
     1114ccc     $  sigmaw(i)
     1115
     1116
     1117c ajout d'un effet onde de gravité -Tgw(k)*deltatw(k) 03/02/06 YU Jingmei
     1118
     1119!      write(lunout,*)'wake.F ',i,k, dtimesub,d_deltat_gw(i,k),
     1120!     &  Tgw(i,k),deltatw(i,k)
     1121        d_deltat_gw(i,k)=d_deltat_gw(i,k)-Tgw(i,k)*deltatw(i,k)*
     1122     $  dtimesub
     1123!      write(lunout,*)'wake.F ',i,k, dtimesub,d_deltatw(i,k)
     1124        ff(i)=d_deltatw(i,k)/dtimesub
     1125
     1126c Sans GW
     1127c
     1128c        deltatw(k)=deltatw(k)+dtimesub*(ff+dtKE(k)-spread(k)*deltatw(k))
     1129c
     1130c GW formule 1
     1131c
     1132c        deltatw(k) = deltatw(k)+dtimesub*
     1133c     $         (ff+dtKE(k) - spread(k)*deltatw(k)-Tgw(k)*deltatw(k))
     1134c
     1135c GW formule 2
     1136
     1137        IF (dtimesub*Tgw(i,k).lt.1.e-10) THEN
     1138          d_deltatw(i,k) = dtimesub*
     1139     $       (ff(i)+dtKE(i,k)+dtPBL(i,k)
     1140ccc     $       -spread(i,k)*deltatw(i,k)
     1141     $       - entr(i,k)*deltatw(i,k)/sigmaw(i)
     1142     $       - (death_rate(i)*sigmaw(i)+detr(i,k))*deltatw(i,k)
     1143     $       / (1.-sigmaw(i))
     1144ccc
     1145     $       -Tgw(i,k)*deltatw(i,k))
     1146        ELSE
     1147           d_deltatw(i,k) = 1/Tgw(i,k)*(1-exp(-dtimesub*
     1148     $       Tgw(i,k)))*
     1149     $       (ff(i)+dtKE(i,k)+dtPBL(i,k)
     1150ccc     $       -spread(i,k)*deltatw(i,k)
     1151     $       - entr(i,k)*deltatw(i,k)/sigmaw(i)
     1152     $       - (death_rate(i)*sigmaw(i)+detr(i,k))*deltatw(i,k)
     1153     $       / (1.-sigmaw(i))
     1154ccc
     1155     $       -Tgw(i,k)*deltatw(i,k))
     1156        ENDIF
     1157
     1158        dth(i,k) = deltatw(i,k)/ppi(i,k)
     1159
     1160        gg(i)=d_deltaqw(i,k)/dtimesub
     1161
     1162       d_deltaqw(i,k) = dtimesub*(gg(i)+ dqKE(i,k)+dqPBL(i,k)
     1163ccc     $     -spread(i,k)*deltaqw(i,k))
     1164     $        - entr(i,k)*deltaqw(i,k)/sigmaw(i)
     1165     $        - (death_rate(i)*sigmaw(i)+detr(i,k))*deltaqw(i,k)
     1166     $        /(1.-sigmaw(i)))
     1167ccc
     1168
     1169ccc nrlmd
     1170ccc       d_deltatw2(i,k)=d_deltatw2(i,k)+d_deltatw(i,k)
     1171ccc       d_deltaqw2(i,k)=d_deltaqw2(i,k)+d_deltaqw(i,k)
     1172ccc
     1173       ENDIF
     1174      ENDDO
     1175      ENDDO
     1176
     1177C
     1178C   Scale tendencies so that water vapour remains positive in w and x.
     1179C
     1180      call wake_vec_modulation(klon,klev,wk_adv,epsilon,qe,d_qe,deltaqw,
     1181     $                d_deltaqw,sigmaw,d_sigmaw,alpha)
     1182c
     1183ccc nrlmd
     1184cc      print*,'alpha'
     1185cc      do i=1,klon
     1186cc         print*,alpha(i)
     1187cc      end do
     1188ccc
     1189      DO k = 1,klev
     1190      DO i = 1,klon
     1191       IF( wk_adv(i) .AND. k .LE. kupper(i)) THEN
     1192        d_te(i,k)=alpha(i)*d_te(i,k)
     1193        d_qe(i,k)=alpha(i)*d_qe(i,k)
     1194        d_deltatw(i,k)=alpha(i)*d_deltatw(i,k)
     1195        d_deltaqw(i,k)=alpha(i)*d_deltaqw(i,k)
     1196        d_deltat_gw(i,k)=alpha(i)*d_deltat_gw(i,k)
     1197       ENDIF
     1198      ENDDO
     1199      ENDDO
     1200      DO i = 1,klon
     1201       IF( wk_adv(i)) THEN
     1202        d_sigmaw(i)=alpha(i)*d_sigmaw(i)
     1203       ENDIF
     1204      ENDDO
     1205
     1206C   Update large scale variables and wake variables
     1207cIM 060208 manque DO i + remplace DO k=1,kupper(i)
     1208cIM 060208     DO k = 1,kupper(i)
     1209      DO k= 1,klev
     1210      DO i = 1,klon
     1211       IF( wk_adv(i) .AND. k .LE. kupper(i)) THEN
     1212        dtls(i,k)=dtls(i,k)+d_te(i,k)
     1213        dqls(i,k)=dqls(i,k)+d_qe(i,k)
     1214ccc nrlmd
     1215        d_deltatw2(i,k)=d_deltatw2(i,k)+d_deltatw(i,k)
     1216        d_deltaqw2(i,k)=d_deltaqw2(i,k)+d_deltaqw(i,k)
     1217ccc
     1218       ENDIF
     1219      ENDDO
     1220      ENDDO
     1221      DO k= 1,klev
     1222      DO i = 1,klon
     1223       IF( wk_adv(i) .AND. k .LE. kupper(i)) THEN
     1224        te(i,k) = te0(i,k) + dtls(i,k)
     1225        qe(i,k) = qe0(i,k) + dqls(i,k)
     1226        the(i,k) = te(i,k)/ppi(i,k)
     1227        deltatw(i,k) = deltatw(i,k)+d_deltatw(i,k)
     1228        deltaqw(i,k) = deltaqw(i,k)+d_deltaqw(i,k)
     1229        dth(i,k) = deltatw(i,k)/ppi(i,k)
     1230cc      print*,'k,qx,qw',k,qe(i,k)-sigmaw(i)*deltaqw(i,k)
     1231cc     $        ,qe(i,k)+(1-sigmaw(i))*deltaqw(i,k)
     1232       ENDIF
     1233      ENDDO
     1234      ENDDO
     1235      DO i = 1,klon
     1236       IF( wk_adv(i)) THEN
     1237        sigmaw(i) = sigmaw(i)+d_sigmaw(i)
     1238       ENDIF
     1239      ENDDO
     1240c
     1241C
     1242c     Determine Ptop from buoyancy integral
     1243c     ---------------------------------------
     1244c
     1245c-     1/ Pressure of the level where dth changes sign.
     1246c
     1247      DO i=1,klon
     1248       IF ( wk_adv(i)) THEN
     1249        Ptop_provis(i)=ph(i,1)
     1250       ENDIF
     1251      ENDDO
     1252c
     1253      DO k= 2,klev
     1254      DO i=1,klon
     1255        IF ( wk_adv(i) .AND.
     1256     $       Ptop_provis(i) .EQ. ph(i,1) .AND.
     1257     $      dth(i,k) .GT. -delta_t_min .and.
     1258     $      dth(i,k-1).LT. -delta_t_min) THEN
     1259          Ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)
     1260     $          - (dth(i,k-1)+delta_t_min)*p(i,k)) /(dth(i,k)
     1261     $          - dth(i,k-1))
     1262        ENDIF
     1263      ENDDO
     1264      ENDDO
     1265c
     1266c-     2/ dth integral
     1267c
     1268      DO i=1,klon
     1269       if (wk_adv(i)) then !!! nrlmd
     1270       sum_dth(i) = 0.
     1271       dthmin(i) = -delta_t_min
     1272       z(i) = 0.
     1273       end if
     1274      ENDDO
     1275
     1276      DO k = 1,klev
     1277      DO i=1,klon
     1278       IF ( wk_adv(i)) THEN
     1279        dz(i) = -(amax1(ph(i,k+1),Ptop_provis(i))-Ph(i,k))/(rho(i,k)*rg)
     1280        IF (dz(i) .gt. 0) THEN
     1281         z(i) = z(i)+dz(i)
     1282         sum_dth(i) = sum_dth(i) + dth(i,k)*dz(i)
     1283         dthmin(i) = amin1(dthmin(i),dth(i,k))
     1284        ENDIF
     1285       ENDIF
     1286      ENDDO
     1287      ENDDO
     1288c
     1289c-     3/ height of triangle with area= sum_dth and base = dthmin
     1290
     1291      DO i=1,klon
     1292       IF ( wk_adv(i)) THEN
     1293         hw(i) = 2.*sum_dth(i)/amin1(dthmin(i),-0.5)
     1294         hw(i) = amax1(hwmin,hw(i))
     1295       ENDIF
     1296      ENDDO
     1297c
     1298c-     4/ now, get Ptop
     1299c
     1300      DO i=1,klon
     1301       if (wk_adv(i)) then !!! nrlmd
     1302       ktop(i) = 0
     1303       z(i)=0.
     1304       end if
     1305      ENDDO
     1306c
     1307      DO k = 1,klev
     1308      DO i=1,klon
     1309       IF ( wk_adv(i)) THEN
     1310        dz(i) = amin1(-(ph(i,k+1)-Ph(i,k))/(rho(i,k)*rg),hw(i)-z(i))
     1311        IF (dz(i) .gt. 0) THEN
     1312         z(i) = z(i)+dz(i)
     1313         Ptop(i) = Ph(i,k)-rho(i,k)*rg*dz(i)
     1314         ktop(i) = k
     1315        ENDIF
     1316       ENDIF
     1317      ENDDO
     1318      ENDDO
     1319c
     1320c      4.5/Correct ktop and ptop
     1321c
     1322      DO i=1,klon
     1323       IF ( wk_adv(i)) THEN
     1324        Ptop_new(i)=ptop(i)
     1325       ENDIF
     1326      ENDDO
     1327c
     1328      DO k= klev,2,-1
     1329      DO i=1,klon
     1330cIM v3JYG; IF (k .GE. ktop(i)
     1331       IF ( wk_adv(i) .AND.
     1332     $      k .LE. ktop(i) .AND.
     1333     $      ptop_new(i) .EQ. ptop(i) .AND.
     1334     $      dth(i,k) .GT. -delta_t_min .and.
     1335     $      dth(i,k-1).LT. -delta_t_min) THEN
     1336          Ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)
     1337     $          - (dth(i,k-1)+delta_t_min)*p(i,k)) /(dth(i,k)
     1338     $          - dth(i,k-1))
     1339        ENDIF
     1340      ENDDO
     1341      ENDDO
     1342c
     1343c
     1344      DO i=1,klon
     1345       IF ( wk_adv(i)) THEN
     1346        ptop(i) = ptop_new(i)
     1347       ENDIF
     1348      ENDDO
     1349
     1350      DO k=klev,1,-1
     1351      DO i=1,klon
     1352      if (wk_adv(i)) then !!! nrlmd
     1353        IF (ph(i,k+1) .LT. ptop(i)) ktop(i)=k
     1354      end if
     1355      ENDDO
     1356      ENDDO
     1357c
     1358c      5/ Set deltatw & deltaqw to 0 above kupper
     1359c
     1360      DO k = 1,klev
     1361      DO i=1,klon
     1362        IF ( wk_adv(i) .AND. k .GE. kupper(i)) THEN
     1363         deltatw(i,k) = 0.
     1364         deltaqw(i,k) = 0.
     1365        ENDIF
     1366      ENDDO
     1367      ENDDO
     1368c
     1369C
     1370c-------------Cstar computation---------------------------------
     1371      DO i=1, klon
     1372       if (wk_adv(i)) then !!! nrlmd
     1373      sum_thu(i) = 0.
     1374      sum_tu(i) = 0.
     1375      sum_qu(i) = 0.
     1376      sum_thvu(i) = 0.
     1377      sum_dth(i) = 0.
     1378      sum_dq(i) = 0.
     1379      sum_rho(i) = 0.
     1380      sum_dtdwn(i) = 0.
     1381      sum_dqdwn(i) = 0.
     1382
     1383      av_thu(i) = 0.
     1384      av_tu(i) =0.
     1385      av_qu(i) =0.
     1386      av_thvu(i) = 0.
     1387      av_dth(i) = 0.
     1388      av_dq(i) = 0.
     1389      av_rho(i) =0.
     1390      av_dtdwn(i) =0.
     1391      av_dqdwn(i) = 0.
     1392       end if
     1393      ENDDO
     1394C
     1395C Integrals (and wake top level number)
     1396C --------------------------------------
     1397C
     1398C Initialize sum_thvu to 1st level virt. pot. temp.
     1399
     1400      DO i=1,klon
     1401       if (wk_adv(i)) then !!! nrlmd
     1402      z(i) = 1.
     1403      dz(i) = 1.
     1404      sum_thvu(i) =  thu(i,1)*(1.+eps*qu(i,1))*dz(i)
     1405      sum_dth(i) = 0.
     1406       end if
     1407      ENDDO
     1408
     1409      DO k = 1,klev
     1410      DO i=1,klon
     1411       if (wk_adv(i)) then !!! nrlmd
     1412        dz(i) = -(max(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
     1413        IF (dz(i) .GT. 0) THEN
     1414         z(i) = z(i)+dz(i)
     1415         sum_thu(i) = sum_thu(i) + thu(i,k)*dz(i)
     1416         sum_tu(i) = sum_tu(i) + tu(i,k)*dz(i)
     1417         sum_qu(i) = sum_qu(i) + qu(i,k)*dz(i)
     1418         sum_thvu(i) = sum_thvu(i) + thu(i,k)*(1.+eps*qu(i,k))*dz(i)
     1419         sum_dth(i) = sum_dth(i) + dth(i,k)*dz(i)
     1420         sum_dq(i) = sum_dq(i) + deltaqw(i,k)*dz(i)
     1421         sum_rho(i) = sum_rho(i) + rhow(i,k)*dz(i)
     1422         sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i,k)*dz(i)
     1423         sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i,k)*dz(i)
     1424        ENDIF
     1425       end if
     1426      ENDDO
     1427      ENDDO
     1428c
     1429      DO i=1,klon
     1430       if (wk_adv(i)) then !!! nrlmd
     1431        hw0(i) = z(i)
     1432       end if
     1433      ENDDO
     1434c
     1435C
     1436C - WAPE and mean forcing computation
     1437C ---------------------------------------
     1438C
     1439C ---------------------------------------
     1440C
     1441C Means
     1442
     1443      DO i=1,klon
     1444       if (wk_adv(i)) then !!! nrlmd
     1445       av_thu(i) = sum_thu(i)/hw0(i)
     1446       av_tu(i) = sum_tu(i)/hw0(i)
     1447       av_qu(i) = sum_qu(i)/hw0(i)
     1448       av_thvu(i) = sum_thvu(i)/hw0(i)
     1449       av_dth(i) = sum_dth(i)/hw0(i)
     1450       av_dq(i) = sum_dq(i)/hw0(i)
     1451       av_rho(i) = sum_rho(i)/hw0(i)
     1452       av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
     1453       av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
     1454c
     1455       wape(i) = - rg*hw0(i)*(av_dth(i)
     1456     $     + eps*(av_thu(i)*av_dq(i)+av_dth(i)*av_qu(i)+av_dth(i)*
     1457     $     av_dq(i) ))/av_thvu(i)
     1458       end if
     1459      ENDDO
     1460C
     1461C Filter out bad wakes
     1462
     1463      DO k = 1,klev
     1464       DO i=1,klon
     1465        if (wk_adv(i)) then !!! nrlmd
     1466        IF ( wape(i) .LT. 0.) THEN
     1467          deltatw(i,k) = 0.
     1468          deltaqw(i,k) = 0.
     1469          dth(i,k) = 0.
     1470        ENDIF
     1471        end if
     1472       ENDDO
     1473      ENDDO
     1474c
     1475      DO i=1,klon
     1476      if (wk_adv(i)) then !!! nrlmd
     1477      IF ( wape(i) .LT. 0.) THEN
     1478        wape(i) = 0.
     1479        Cstar(i) = 0.
     1480        hw(i) = hwmin
     1481        sigmaw(i) = max(sigmad,sigd_con(i))
     1482        fip(i) = 0.
     1483        gwake(i) = .FALSE.
     1484      ELSE
     1485        Cstar(i) = stark*sqrt(2.*wape(i))
     1486        gwake(i) = .TRUE.
     1487      ENDIF
     1488      end if
     1489      ENDDO
     1490
     1491       ENDDO      ! end sub-timestep loop
     1492C
     1493C -----------------------------------------------------------------
     1494c   Get back to tendencies per second
     1495c
     1496      DO k = 1,klev
     1497      DO i=1,klon
     1498
     1499ccc nrlmd        IF ( wk_adv(i) .AND. k .LE. kupper(i)) THEN
     1500        IF ( OK_qx_qw(i) .AND. k .LE. kupper(i)) THEN
     1501ccc
     1502         dtls(i,k) = dtls(i,k)/dtime
     1503         dqls(i,k) = dqls(i,k)/dtime
     1504         d_deltatw2(i,k)=d_deltatw2(i,k)/dtime
     1505         d_deltaqw2(i,k)=d_deltaqw2(i,k)/dtime
     1506         d_deltat_gw(i,k) = d_deltat_gw(i,k)/dtime
     1507cc      print*,'k,dqls,omg,entr,detr',k,dqls(i,k),omg(i,k),entr(i,k)
     1508cc     $         ,death_rate(i)*sigmaw(i)
     1509        ENDIF
     1510      ENDDO
     1511      ENDDO
     1512
     1513c
     1514c----------------------------------------------------------
     1515c   Determine wake final state; recompute wape, cstar, ktop;
     1516c   filter out bad wakes.
     1517c----------------------------------------------------------
     1518c
     1519C 2.1 - Undisturbed area and Wake integrals
     1520C ---------------------------------------------------------
     1521
     1522      DO i=1,klon
     1523ccc nrlmd       if (wk_adv(i)) then !!! nrlmd
     1524      if (OK_qx_qw(i)) then
     1525ccc
     1526        z(i) = 0.
     1527        sum_thu(i) = 0.
     1528        sum_tu(i) = 0.
     1529        sum_qu(i) = 0.
     1530        sum_thvu(i) = 0.
     1531        sum_dth(i) = 0.
     1532        sum_dq(i) = 0.
     1533        sum_rho(i) = 0.
     1534        sum_dtdwn(i) = 0.
     1535        sum_dqdwn(i) = 0.
     1536
     1537        av_thu(i) = 0.
     1538        av_tu(i) =0.
     1539        av_qu(i) =0.
     1540        av_thvu(i) = 0.
     1541        av_dth(i) = 0.
     1542        av_dq(i) = 0.
     1543        av_rho(i) =0.
     1544        av_dtdwn(i) =0.
     1545        av_dqdwn(i) = 0.
     1546       end if   
     1547      ENDDO
     1548C Potential temperatures and humidity
     1549c----------------------------------------------------------
     1550
     1551      DO k =1,klev
     1552      DO i=1,klon
     1553ccc nrlmd       IF ( wk_adv(i)) THEN
     1554       if (OK_qx_qw(i)) then
     1555ccc
     1556        rho(i,k) = p(i,k)/(rd*te(i,k))
     1557        IF(k .eq. 1) THEN
     1558          rhoh(i,k) = ph(i,k)/(rd*te(i,k))
     1559          zhh(i,k)=0
     1560        ELSE
     1561          rhoh(i,k) = ph(i,k)*2./(rd*(te(i,k)+te(i,k-1)))
     1562          zhh(i,k)=(ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*RG)+zhh(i,k-1)
     1563        ENDIF
     1564        the(i,k) = te(i,k)/ppi(i,k)
     1565        thu(i,k) = (te(i,k) - deltatw(i,k)*sigmaw(i))/ppi(i,k)
     1566        tu(i,k) = te(i,k) - deltatw(i,k)*sigmaw(i)
     1567        qu(i,k)  =  qe(i,k) - deltaqw(i,k)*sigmaw(i)
     1568        rhow(i,k) = p(i,k)/(rd*(te(i,k)+deltatw(i,k)))
     1569        dth(i,k) = deltatw(i,k)/ppi(i,k)
     1570       ENDIF
     1571      ENDDO
     1572      ENDDO
     1573
     1574C Integrals (and wake top level number)
     1575C -----------------------------------------------------------
     1576
     1577C Initialize sum_thvu to 1st level virt. pot. temp.
     1578
     1579      DO i=1,klon
     1580ccc nrlmd       IF ( wk_adv(i)) THEN
     1581      if (OK_qx_qw(i)) then
     1582ccc
     1583        z(i) = 1.
     1584        dz(i) = 1.
     1585        sum_thvu(i) =  thu(i,1)*(1.+eps*qu(i,1))*dz(i)
     1586        sum_dth(i) = 0.
     1587      ENDIF
     1588      ENDDO
     1589
     1590      DO k = 1,klev
     1591      DO i=1,klon
     1592ccc nrlmd       IF ( wk_adv(i)) THEN
     1593       if (OK_qx_qw(i)) then
     1594ccc
     1595        dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
     1596        IF (dz(i) .GT. 0) THEN
     1597         z(i) = z(i)+dz(i)
     1598         sum_thu(i) = sum_thu(i) + thu(i,k)*dz(i)
     1599         sum_tu(i) = sum_tu(i) + tu(i,k)*dz(i)
     1600         sum_qu(i) = sum_qu(i) + qu(i,k)*dz(i)
     1601         sum_thvu(i) = sum_thvu(i) + thu(i,k)*(1.+eps*qu(i,k))*dz(i)
     1602         sum_dth(i) = sum_dth(i) + dth(i,k)*dz(i)
     1603         sum_dq(i) = sum_dq(i) + deltaqw(i,k)*dz(i)
     1604         sum_rho(i) = sum_rho(i) + rhow(i,k)*dz(i)
     1605         sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i,k)*dz(i)
     1606         sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i,k)*dz(i)
     1607        ENDIF
     1608       ENDIF
     1609      ENDDO
     1610      ENDDO
     1611c
     1612      DO i=1,klon
     1613ccc nrlmd       IF ( wk_adv(i)) THEN
     1614       if (OK_qx_qw(i)) then
     1615ccc
     1616        hw0(i) = z(i)
     1617       ENDIF
     1618      ENDDO
     1619c
     1620C - WAPE and mean forcing computation
     1621C-------------------------------------------------------------
     1622
     1623C Means
     1624
     1625      DO i=1, klon
     1626ccc nrlmd       IF ( wk_adv(i)) THEN
     1627      if (OK_qx_qw(i)) then
     1628ccc
     1629        av_thu(i) = sum_thu(i)/hw0(i)
     1630        av_tu(i) = sum_tu(i)/hw0(i)
     1631        av_qu(i) = sum_qu(i)/hw0(i)
     1632        av_thvu(i) = sum_thvu(i)/hw0(i)
     1633        av_dth(i) = sum_dth(i)/hw0(i)
     1634        av_dq(i) = sum_dq(i)/hw0(i)
     1635        av_rho(i) = sum_rho(i)/hw0(i)
     1636        av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
     1637        av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
     1638
     1639        wape2(i) = - rg*hw0(i)*(av_dth(i)
     1640     $     + eps*(av_thu(i)*av_dq(i)+av_dth(i)*av_qu(i)
     1641     $     + av_dth(i)*av_dq(i) ))/av_thvu(i)
     1642       ENDIF
     1643      ENDDO
     1644
     1645C Prognostic variable update
     1646C ------------------------------------------------------------
     1647
     1648C Filter out bad wakes
     1649c
     1650      DO k = 1,klev
     1651      DO i=1,klon
     1652ccc nrlmd        IF ( wk_adv(i) .AND. wape2(i) .LT. 0.) THEN
     1653      if (OK_qx_qw(i) .AND. wape2(i) .LT. 0.) then
     1654ccc
     1655          deltatw(i,k) = 0.
     1656          deltaqw(i,k) = 0.
     1657          dth(i,k) = 0.
     1658        ENDIF
     1659      ENDDO
     1660      ENDDO
     1661c
     1662
     1663      DO i=1, klon
     1664ccc nrlmd       IF ( wk_adv(i)) THEN
     1665      if (OK_qx_qw(i)) then
     1666ccc
     1667       IF ( wape2(i) .LT. 0.) THEN
     1668        wape2(i) = 0.
     1669        Cstar2(i) = 0.
     1670        hw(i) = hwmin
     1671        sigmaw(i) = amax1(sigmad,sigd_con(i))
     1672        fip(i) = 0.
     1673        gwake(i) = .FALSE.
     1674      ELSE
     1675        if(prt_level.ge.10) print*,'wape2>0'
     1676        Cstar2(i) = stark*sqrt(2.*wape2(i))
     1677        gwake(i) = .TRUE.
     1678      ENDIF
     1679      ENDIF
     1680      ENDDO
     1681c
     1682      DO i=1, klon
     1683ccc nrlmd       IF ( wk_adv(i)) THEN
     1684       if (OK_qx_qw(i)) then
     1685ccc
     1686        ktopw(i) = ktop(i)
     1687       ENDIF
     1688      ENDDO
     1689c
     1690      DO i=1, klon
     1691ccc nrlmd       IF ( wk_adv(i)) THEN
     1692       if (OK_qx_qw(i)) then
     1693ccc
     1694       IF (ktopw(i) .gt. 0 .and. gwake(i)) then
     1695
     1696Cjyg1     Utilisation d'un h_efficace constant ( ~ feeding layer)
     1697ccc       heff = 600.
     1698C      Utilisation de la hauteur hw
     1699cc       heff = 0.7*hw
     1700       heff(i) = hw(i)
     1701
     1702       FIP(i) = 0.5*rho(i,ktopw(i))*Cstar2(i)**3*heff(i)*2*
     1703     $      sqrt(sigmaw(i)*wdens(i)*3.14)
     1704       FIP(i) = alpk * FIP(i)
     1705Cjyg2
     1706       ELSE
     1707         FIP(i) = 0.
     1708       ENDIF
     1709       ENDIF
     1710      ENDDO
     1711c
     1712C   Limitation de sigmaw
     1713
     1714ccc nrlmd
     1715c       DO i=1,klon
     1716c         IF (OK_qx_qw(i)) THEN
     1717c          IF (sigmaw(i).GE.sigmaw_max) sigmaw(i)=sigmaw_max
     1718c        ENDIF
     1719c       ENDDO
     1720ccc
     1721      DO k = 1,klev
     1722       DO i=1, klon
     1723
     1724ccc nrlmd      On maintient désormais constant sigmaw en régime permanent
     1725ccc      IF ((sigmaw(i).GT.sigmaw_max).or.
     1726        IF     ( ((wape(i).ge.wape2(i)).and.(wape2(i).le.1.0)).or.
     1727     $         (ktopw(i).le.2) .OR.
     1728     $         .not. OK_qx_qw(i) ) THEN
     1729ccc
     1730          dtls(i,k) = 0.
     1731          dqls(i,k) = 0.
     1732          deltatw(i,k) = 0.
     1733          deltaqw(i,k) = 0.
     1734        ENDIF
     1735       ENDDO
     1736      ENDDO
     1737c
     1738ccc nrlmd      On maintient désormais constant sigmaw en régime permanent
     1739      DO i=1, klon
     1740        IF  ( ((wape(i).ge.wape2(i)).and.(wape2(i).le.1.0)).or.
     1741     $        (ktopw(i).le.2) .OR.
     1742     $        .not. OK_qx_qw(i)   ) THEN
     1743         wape(i) = 0.
     1744         cstar(i)=0.
     1745         hw(i) = hwmin
     1746         sigmaw(i) = sigmad
     1747         fip(i) = 0.
     1748        ELSE
     1749         wape(i) = wape2(i)
     1750         cstar(i)=cstar2(i)
     1751        ENDIF
     1752cc        print*,'wape wape2 ktopw OK_qx_qw =',
     1753cc     $          wape(i),wape2(i),ktopw(i),OK_qx_qw(i)
     1754      ENDDO
     1755c
     1756c
     1757      RETURN
     1758      END
     1759
     1760      SUBROUTINE wake_vec_modulation(nlon,nl,wk_adv,epsilon,qe,d_qe,
     1761     $           deltaqw,d_deltaqw,sigmaw,d_sigmaw,alpha)
     1762c------------------------------------------------------
     1763cDtermination du coefficient alpha tel que les tendances
     1764c corriges alpha*d_G, pour toutes les grandeurs G, correspondent
     1765c a une humidite positive dans la zone (x) et dans la zone (w).
     1766c------------------------------------------------------
     1767c
     1768 
     1769c  Input
     1770      REAL qe(nlon,nl),d_qe(nlon,nl)
     1771      REAL deltaqw(nlon,nl),d_deltaqw(nlon,nl)
     1772      REAL sigmaw(nlon),d_sigmaw(nlon)
     1773      LOGICAL wk_adv(nlon)
     1774      INTEGER nl,nlon
     1775c  Output
     1776      REAL alpha(nlon)
     1777c  Internal variables
     1778      REAL zeta(nlon,nl)
     1779      REAL alpha1(nlon)
     1780      REAL x,a,b,c,discrim
     1781      REAL epsilon
     1782!      DATA epsilon/1.e-15/
     1783c
     1784      DO k=1,nl
     1785      DO i = 1,nlon
     1786       IF (wk_adv(i)) THEN
     1787        IF ((deltaqw(i,k)+d_deltaqw(i,k)).ge.0.) then
     1788         zeta(i,k)=0.
     1789        ELSE
     1790         zeta(i,k)=1.
     1791        END IF
     1792       ENDIF
     1793      ENDDO
     1794      DO i = 1,nlon
     1795       IF (wk_adv(i)) THEN
     1796        x = qe(i,k)+(zeta(i,k)-sigmaw(i))*deltaqw(i,k)
     1797     $    + d_qe(i,k)+(zeta(i,k)-sigmaw(i))*d_deltaqw(i,k)
     1798     $    - d_sigmaw(i)*(deltaqw(i,k)+d_deltaqw(i,k))
     1799        a = -d_sigmaw(i)*d_deltaqw(i,k)
     1800        b = d_qe(i,k)+(zeta(i,k)-sigmaw(i))*d_deltaqw(i,k)
     1801     $    - deltaqw(i,k)*d_sigmaw(i)
     1802        c = qe(i,k)+(zeta(i,k)-sigmaw(i))*deltaqw(i,k)+epsilon
     1803        discrim = b*b-4.*a*c
     1804c      print*, 'x, a, b, c, discrim', x, a, b, c, discrim
     1805        IF (a+b .GE. 0.) THEN !! Condition suffisante pour la positivité de ovap
     1806         alpha1(i)=1.
     1807        ELSE
     1808         IF (x .GE. 0.) THEN
     1809            alpha1(i)=1.
     1810         ELSE
     1811              IF (a .GT. 0.) THEN
     1812                 alpha1(i)=0.9*min(   (2.*c)/(-b+sqrt(discrim)),
     1813     $                        (-b+sqrt(discrim))/(2.*a)   )
     1814              ELSE IF (a .eq. 0.) then
     1815                 alpha1(i)=0.9*(-c/b)
     1816              ELSE
     1817c         print*,'a,b,c discrim',a,b,c discrim
     1818                 alpha1(i)=0.9*max(   (2.*c)/(-b+sqrt(discrim)),
     1819     $                        (-b+sqrt(discrim))/(2.*a)   )
     1820              ENDIF
     1821         ENDIF
     1822        ENDIF
     1823       alpha(i) = min(alpha(i),alpha1(i))
     1824       ENDIF
     1825      ENDDO
     1826      ENDDO
     1827!
     1828      return
     1829      end
     1830
     1831      Subroutine WAKE_scal (p,ph,ppi,dtime,sigd_con
     1832     :                ,te0,qe0,omgb
     1833     :                ,dtdwn,dqdwn,amdwn,amup,dta,dqa
     1834     :                ,wdtPBL,wdqPBL,udtPBL,udqPBL
     1835     o                ,deltatw,deltaqw,dth,hw,sigmaw,wape,fip,gfl
     1836     o                ,dtls,dqls
     1837     o                ,ktopw,omgbdth,dp_omgb,wdens
     1838     o                ,tu,qu
     1839     o                ,dtKE,dqKE
     1840     o                ,dtPBL,dqPBL
     1841     o                ,omg,dp_deltomg,spread
     1842     o                ,Cstar,d_deltat_gw
     1843     o                ,d_deltatw2,d_deltaqw2)
     1844
     1845***************************************************************
     1846*                                                             *
     1847* WAKE                                                        *
     1848*      retour a un Pupper fixe                                *
     1849*                                                             *
     1850* written by   :  GRANDPEIX Jean-Yves   09/03/2000            *
     1851* modified by :   ROEHRIG Romain        01/29/2007            *
     1852***************************************************************
     1853c
     1854      USE dimphy
    261855      IMPLICIT none
    271856c============================================================================
    … …  
    1131942
    1141943#include "dimensions.h"
    115 #include "YOMCST.h"
    116 #include "cvthermo.h"
    117 #include "iniprint.h"
    118 
    119 c Arguments en entree
    120 c--------------------
    121 
    122       REAL, dimension(klon,klev) :: p, ppi
    123       REAL, dimension(klon,klev+1) :: ph, omgb
    124       REAL dtime
    125       REAL, dimension(klon,klev) :: te0,qe0
    126       REAL, dimension(klon,klev) :: dtdwn, dqdwn
    127       REAL, dimension(klon,klev) :: wdtPBL,wdqPBL
    128       REAL, dimension(klon,klev) :: udtPBL,udqPBL
    129       REAL, dimension(klon,klev) :: amdwn, amup
    130       REAL, dimension(klon,klev) :: dta, dqa
    131       REAL, dimension(klon) :: sigd_con
    132 
    133 c Sorties
    134 c--------
    135 
    136       REAL, dimension(klon,klev) :: deltatw, deltaqw, dth
    137       REAL, dimension(klon,klev) :: tu, qu
    138       REAL, dimension(klon,klev) :: dtls, dqls
    139       REAL, dimension(klon,klev) :: dtKE, dqKE
    140       REAL, dimension(klon,klev) :: dtPBL, dqPBL
    141       REAL, dimension(klon,klev) :: spread
    142       REAL, dimension(klon,klev) :: d_deltatgw
    143       REAL, dimension(klon,klev) :: d_deltatw2, d_deltaqw2
    144       REAL, dimension(klon,klev+1) :: omgbdth, omg
    145       REAL, dimension(klon,klev) :: dp_omgb, dp_deltomg
    146       REAL, dimension(klon,klev) :: d_deltat_gw
    147       REAL, dimension(klon) :: hw, sigmaw, wape, fip, gfl, Cstar
    148       INTEGER, dimension(klon) :: ktopw
    149 
    150 c Variables internes
    151 c-------------------
    152 
    153 c Variables à fixer
    154       REAL ALON
    155       REAL coefgw
    156       REAL :: wdens0, wdens
    157       REAL stark
    158       REAL alpk
    159       REAL delta_t_min
    160       INTEGER nsub
    161       REAL dtimesub
    162       REAL sigmad, hwmin
    163       REAL :: sigmaw_max
    164 cIM 080208
    165       LOGICAL, dimension(klon) :: gwake
    166 
    167 c Variables de sauvegarde
    168       REAL, dimension(klon,klev) :: deltatw0
    169       REAL, dimension(klon,klev) :: deltaqw0
    170       REAL, dimension(klon,klev) :: te, qe
    171       REAL, dimension(klon) :: sigmaw0, sigmaw1
    172 
    173 c Variables pour les GW
    174       REAL, DIMENSION(klon) :: LL
    175       REAL, dimension(klon,klev) :: N2
    176       REAL, dimension(klon,klev) :: Cgw
    177       REAL, dimension(klon,klev) :: Tgw
    178 
    179 c Variables liées au calcul de hw
    180       REAL, DIMENSION(klon) :: ptop_provis, ptop, ptop_new
    181       REAL, DIMENSION(klon) :: sum_dth
    182       REAL, DIMENSION(klon) :: dthmin
    183       REAL, DIMENSION(klon) :: z, dz, hw0
    184       INTEGER, DIMENSION(klon) :: ktop, kupper
    185 
    186 c 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 
    193 c Autres variables internes
    194       INTEGER isubstep, k, i
    195 
    196       REAL, DIMENSION(klon) :: sum_thu, sum_tu, sum_qu,sum_thvu
    197       REAL, DIMENSION(klon) :: sum_dq, sum_rho
    198       REAL, DIMENSION(klon) :: sum_dtdwn, sum_dqdwn
    199       REAL, DIMENSION(klon) :: av_thu, av_tu, av_qu, av_thvu
    200       REAL, DIMENSION(klon) :: av_dth, av_dq, av_rho
    201       REAL, DIMENSION(klon) :: av_dtdwn, av_dqdwn
    202 
    203       REAL, DIMENSION(klon,klev) :: rho, rhow
    204       REAL, DIMENSION(klon,klev+1) :: rhoh
    205       REAL, DIMENSION(klon,klev) :: rhow_moyen
    206       REAL, DIMENSION(klon,klev) :: zh
    207       REAL, DIMENSION(klon,klev+1) :: zhh
    208       REAL, DIMENSION(klon,klev) :: epaisseur1, epaisseur2
    209 
    210       REAL, DIMENSION(klon,klev) :: the, thu
    211 
    212 !      REAL, DIMENSION(klon,klev) :: d_deltatw, d_deltaqw
    213 
    214       REAL, DIMENSION(klon,klev+1) :: omgbw
    215       REAL, DIMENSION(klon) :: pupper
    216       REAL, DIMENSION(klon) :: omgtop
    217       REAL, DIMENSION(klon,klev) :: dp_omgbw
    218       REAL, DIMENSION(klon) :: ztop, dztop
    219       REAL, DIMENSION(klon,klev) :: alpha_up
    220      
    221       REAL, dimension(klon) :: RRe1, RRe2
    222       REAL :: RRd1, RRd2
    223       REAL, DIMENSION(klon,klev) :: Th1, Th2, q1, q2, T1
    224       REAL, DIMENSION(klon,klev) :: D_Th1, D_Th2, D_dth
    225       REAL, DIMENSION(klon,klev) :: D_q1, D_q2, D_dq
    226       REAL, DIMENSION(klon,klev) :: omgbdq
    227 
    228       REAL, dimension(klon) :: ff, gg
    229       REAL, dimension(klon) :: wape2, Cstar2, heff
    230 
    231       REAL, DIMENSION(klon,klev) :: Crep
    232       REAL Crep_upper, Crep_sol
    233 
    234 C-------------------------------------------------------------------------
    235 c         Initialisations
    236 c-------------------------------------------------------------------------
    237 
    238 c      print*, 'wake initialisations'
    239 
    240 c   Essais d'initialisation avec sigmaw = 0.02 et hw = 10.
    241 c-------------------------------------------------------------------------
    242 
    243       DATA sigmad, hwmin /.02,10./
    244 
    245 C Longueur de maille (en m)
    246 c-------------------------------------------------------------------------
    247 
    248 c      ALON = 3.e5
    249       ALON = 1.e6
    250 
    251 
    252 C Configuration de coefgw,stark,wdens (22/02/06 by YU Jingmei)
    253 c
    254 c      coefgw : Coefficient pour les ondes de gravité
    255 c       stark : Coefficient k dans Cstar=k*sqrt(2*WAPE)
    256 c       wdens : Densité de poche froide par maille
    257 c-------------------------------------------------------------------------
    258 
    259       coefgw=10
    260 c      coefgw=1
    261 c      wdens0 = 1.0/(alon**2)   
    262       wdens = 1.0/(alon**2)       
    263       stark = 0.50
    264 cCRtest
    265       alpk=0.1
    266 c      alpk = 1.0
    267 c      alpk = 0.5
    268 c      alpk = 0.05
    269       Crep_upper=0.9
    270       Crep_sol=1.0
    271 
    272 C Minimum value for |T_wake - T_undist|. Used for wake top definition
    273 c-------------------------------------------------------------------------
    274 
    275       delta_t_min = 0.2
    276 
    277 C 1. - Save initial values and initialize tendencies
    278 C --------------------------------------------------
    279 
    280       DO k=1,klev
    281       DO i=1, klon
    282         deltatw0(i,k) = deltatw(i,k)
    283         deltaqw0(i,k)= deltaqw(i,k)
    284         te(i,k) = te0(i,k)
    285         qe(i,k) = qe0(i,k)
    286         dtls(i,k) = 0.
    287         dqls(i,k) = 0.
    288         d_deltat_gw(i,k)=0.
    289         d_te(i,k) = 0.
    290         d_qe(i,k) = 0.
    291         d_deltatw(i,k) = 0.
    292         d_deltaqw(i,k) = 0.
    293 !IM 060508 beg
    294         d_deltatw2(i,k)=0.
    295         d_deltaqw2(i,k)=0.
    296 !IM 060508 end
    297       ENDDO
    298       ENDDO
    299 c      sigmaw1=sigmaw
    300 c      IF (sigd_con.GT.sigmaw1) THEN
    301 c      print*, 'sigmaw,sigd_con', sigmaw, sigd_con
    302 c      ENDIF
    303       DO i=1, klon
    304 cc      sigmaw(i) = amax1(sigmaw(i),sigd_con(i))
    305       sigmaw(i) = amax1(sigmaw(i),sigmad)
    306       sigmaw(i) = amin1(sigmaw(i),0.99)
    307       sigmaw0(i) = sigmaw(i)
    308       wape(i) = 0.
    309       wape2(i) = 0.
    310       d_sigmaw(i) = 0.
    311       ktopw(i) = 0
    312       ENDDO
    313 C
    314 C
    315 C 2. - Prognostic part
    316 C --------------------
    317 C
    318 C
    319 C 2.1 - Undisturbed area and Wake integrals
    320 C ---------------------------------------------------------
    321 
    322       DO i=1, klon
    323       z(i) = 0.
    324       ktop(i)=0
    325       kupper(i) = 0
    326       sum_thu(i) = 0.
    327       sum_tu(i) = 0.
    328       sum_qu(i) = 0.
    329       sum_thvu(i) = 0.
    330       sum_dth(i) = 0.
    331       sum_dq(i) = 0.
    332       sum_rho(i) = 0.
    333       sum_dtdwn(i) = 0.
    334       sum_dqdwn(i) = 0.
    335 
    336       av_thu(i) = 0.
    337       av_tu(i) =0.
    338       av_qu(i) =0.
    339       av_thvu(i) = 0.
    340       av_dth(i) = 0.
    341       av_dq(i) = 0.
    342       av_rho(i) =0.
    343       av_dtdwn(i) =0.
    344       av_dqdwn(i) = 0.
    345       ENDDO
    346 c
    347 c Distance between wakes
    348        DO i = 1,klon
    349         LL(i) = (1-sqrt(sigmaw(i)))/sqrt(wdens)
    350        ENDDO
    351 C Potential temperatures and humidity
    352 c----------------------------------------------------------
    353       DO k =1,klev
    354        DO i=1, klon
    355         rho(i,k) = p(i,k)/(rd*te(i,k))
    356         IF(k .eq. 1) THEN
    357           rhoh(i,k) = ph(i,k)/(rd*te(i,k))
    358           zhh(i,k)=0
    359         ELSE
    360           rhoh(i,k) = ph(i,k)*2./(rd*(te(i,k)+te(i,k-1)))
    361           zhh(i,k)=(ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*RG)+zhh(i,k-1)
    362         ENDIF
    363         the(i,k) = te(i,k)/ppi(i,k)
    364         thu(i,k) = (te(i,k) - deltatw(i,k)*sigmaw(i))/ppi(i,k)
    365         tu(i,k) = te(i,k) - deltatw(i,k)*sigmaw(i)
    366         qu(i,k)  =  qe(i,k) - deltaqw(i,k)*sigmaw(i)
    367         rhow(i,k) = p(i,k)/(rd*(tu(i,k)+deltatw(i,k)))
    368         dth(i,k) = deltatw(i,k)/ppi(i,k)
    369        ENDDO
    370       ENDDO
    371        
    372       DO k = 1, klev-1
    373       DO i=1, klon
    374         IF(k.eq.1) THEN
    375           N2(i,k)=0
    376         ELSE
    377           N2(i,k)=amax1(0.,-RG**2/the(i,k)*rho(i,k)*(the(i,k+1)-
    378      $            the(i,k-1))/(p(i,k+1)-p(i,k-1)))
    379         ENDIF
    380         ZH(i,k)=(zhh(i,k)+zhh(i,k+1))/2
    381 
    382         Cgw(i,k)=sqrt(N2(i,k))*ZH(i,k)
    383         Tgw(i,k)=coefgw*Cgw(i,k)/LL(i)
    384       ENDDO
    385       ENDDO
    386 
    387       DO i=1, klon
    388       N2(i,klev)=0
    389       ZH(i,klev)=0
    390       Cgw(i,klev)=0
    391       Tgw(i,klev)=0
    392       ENDDO
    393 
    394 c  Calcul de la masse volumique moyenne de la colonne   (bdlmd)
    395 c-----------------------------------------------------------------
    396 
    397       DO k=1,klev
    398        DO i=1, klon
    399         epaisseur1(i,k)=0.
    400         epaisseur2(i,k)=0.
    401        ENDDO
    402       ENDDO
    403 
    404       DO i=1, klon
    405       epaisseur1(i,1)= -(ph(i,2)-ph(i,1))/(rho(i,1)*rg)+1.
    406       epaisseur2(i,1)= -(ph(i,2)-ph(i,1))/(rho(i,1)*rg)+1.
    407       rhow_moyen(i,1) = rhow(i,1)
    408       ENDDO
    409 
    410       DO k = 2, klev
    411       DO i=1, klon
    412         epaisseur1(i,k)= -(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg) +1.
    413         epaisseur2(i,k)=epaisseur2(i,k-1)+epaisseur1(i,k)
    414         rhow_moyen(i,k) = (rhow_moyen(i,k-1)*epaisseur2(i,k-1)+
    415      $                 rhow(i,k)*epaisseur1(i,k))/epaisseur2(i,k)
    416       ENDDO
    417       ENDDO
    418 
    419 C
    420 C Choose an integration bound well above wake top
    421 c-----------------------------------------------------------------
    422 c
    423 C       Pupper = 50000.  ! melting level
    424 c       Pupper = 60000.
    425 c       Pupper = 80000.  ! essais pour case_e
    426        DO i = 1,klon
    427 ccc       Pupper(i) = 0.6*ph(i,1)
    428         Pupper(i) = 60000.
    429        ENDDO
    430 
    431 C
    432 C    Determine Wake top pressure (Ptop) from buoyancy integral
    433 C    --------------------------------------------------------
    434 c
    435 c-1/ Pressure of the level where dth becomes less than delta_t_min.
    436 
    437       DO i=1,klon
    438       ptop_provis(i)=ph(i,1)
    439       ENDDO
    440       DO k= 2,klev
    441       DO i=1,klon
    442 c
    443 cIM v3JYG; ptop_provis(i).LT. ph(i,1)
    444 c
    445         IF (dth(i,k) .GT. -delta_t_min .and.
    446      $      dth(i,k-1).LT. -delta_t_min .and.
    447      $      ptop_provis(i).EQ. ph(i,1)) THEN
    448           ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)
    449      $          - (dth(i,k-1)+delta_t_min)*p(i,k)) /
    450      $          (dth(i,k) - dth(i,k-1))
    451         ENDIF
    452       ENDDO
    453       ENDDO
    454 
    455 c-2/ dth integral
    456 
    457       DO i=1,klon
    458       sum_dth(i) = 0.
    459       dthmin(i) = -delta_t_min
    460       z(i) = 0.
    461       ENDDO
    462 
    463       DO k = 1,klev
    464       DO i=1,klon
    465         dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-Ph(i,k))/(rho(i,k)*rg)
    466         IF (dz(i) .gt. 0) THEN
    467           z(i) = z(i)+dz(i)
    468           sum_dth(i) = sum_dth(i) + dth(i,k)*dz(i)
    469           dthmin(i) = amin1(dthmin(i),dth(i,k))
    470         ENDIF
    471       ENDDO
    472       ENDDO
    473 
    474 c-3/ height of triangle with area= sum_dth and base = dthmin
    475 
    476       DO i=1,klon
    477       hw0(i) = 2.*sum_dth(i)/amin1(dthmin(i),-0.5)
    478       hw0(i) = amax1(hwmin,hw0(i))
    479       ENDDO
    480 
    481 c-4/ now, get Ptop
    482 
    483       DO i=1,klon
    484       z(i) = 0.
    485       ptop(i) = ph(i,1)
    486       ENDDO
    487 
    488       DO k = 1,klev
    489       DO i=1,klon
    490         dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg),hw0(i)-z(i))
    491         IF (dz(i) .gt. 0) THEN
    492          z(i) = z(i)+dz(i)
    493          ptop(i) = ph(i,k)-rho(i,k)*rg*dz(i)
    494         ENDIF
    495       ENDDO
    496       ENDDO
    497 
    498 
    499 C-5/ Determination de ktop et kupper
    500 
    501       DO k=klev,1,-1
    502       DO i=1,klon
    503         IF (ph(i,k+1) .lt. ptop(i)) ktop(i)=k
    504         IF (ph(i,k+1) .lt. pupper(i)) kupper(i)=k
    505       ENDDO
    506       ENDDO
    507 
    508 c-6/ Correct ktop and ptop
    509 
    510       DO i = 1,klon
    511         ptop_new(i)=ptop(i)
    512       ENDDO
    513       DO k= klev,2,-1
    514       DO i=1,klon
    515         IF (k .LE. ktop(i) .and.
    516      $      ptop_new(i) .EQ. ptop(i) .and.
    517      $      dth(i,k) .GT. -delta_t_min .and.
    518      $      dth(i,k-1).LT. -delta_t_min) THEN
    519           ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)
    520      $          - (dth(i,k-1)+delta_t_min)*p(i,k)) /
    521      $          (dth(i,k) - dth(i,k-1))
    522         ENDIF
    523       ENDDO
    524       ENDDO
    525 
    526       DO i=1,klon
    527         ptop(i) = ptop_new(i)
    528       ENDDO
    529 
    530       DO k=klev,1,-1
    531       DO i=1,klon
    532         IF (ph(i,k+1) .lt. ptop(i)) ktop(i)=k
    533       ENDDO
    534       ENDDO
    535 c
    536 c-5/ Set deltatw & deltaqw to 0 above kupper
    537 c
    538       DO k = 1,klev
    539       DO i=1,klon
    540        IF (k.GE. kupper(i)) THEN
    541         deltatw(i,k) = 0.
    542         deltaqw(i,k) = 0.
    543        ENDIF
    544       ENDDO
    545       ENDDO
    546 c
    547 C
    548 C Vertical gradient of LS omega
    549 C
    550       DO k = 1,klev
    551       DO i=1,klon
    552        IF (k.LE. kupper(i)) THEN
    553         dp_omgb(i,k) = (omgb(i,k+1) - omgb(i,k))/(ph(i,k+1)-ph(i,k))
    554        ENDIF
    555       ENDDO
    556       ENDDO
    557 C
    558 C Integrals (and wake top level number)
    559 C --------------------------------------
    560 C
    561 C Initialize sum_thvu to 1st level virt. pot. temp.
    562 
    563       DO i=1,klon
    564       z(i) = 1.
    565       dz(i) = 1.
    566       sum_thvu(i) =  thu(i,1)*(1.+eps*qu(i,1))*dz(i)
    567       sum_dth(i) = 0.
    568       ENDDO
    569 
    570       DO k = 1,klev
    571       DO i=1,klon
    572         dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
    573         IF (dz(i) .GT. 0) THEN
    574          z(i) = z(i)+dz(i)
    575          sum_thu(i) = sum_thu(i) + thu(i,k)*dz(i)
    576          sum_tu(i) = sum_tu(i) + tu(i,k)*dz(i)
    577          sum_qu(i) = sum_qu(i) + qu(i,k)*dz(i)
    578          sum_thvu(i) = sum_thvu(i) + thu(i,k)*(1.+eps*qu(i,k))*dz(i)
    579          sum_dth(i) = sum_dth(i) + dth(i,k)*dz(i)
    580          sum_dq(i) = sum_dq(i) + deltaqw(i,k)*dz(i)
    581          sum_rho(i) = sum_rho(i) + rhow(i,k)*dz(i)
    582          sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i,k)*dz(i)
    583          sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i,k)*dz(i)
    584         ENDIF
    585       ENDDO
    586       ENDDO
    587 c
    588       DO i=1,klon
    589         hw0(i) = z(i)
    590       ENDDO
    591 c
    592 C
    593 C 2.1 - WAPE and mean forcing computation
    594 C ---------------------------------------
    595 C
    596 C ---------------------------------------
    597 C
    598 C Means
    599 
    600       DO i=1,klon
    601       av_thu(i) = sum_thu(i)/hw0(i)
    602       av_tu(i) = sum_tu(i)/hw0(i)
    603       av_qu(i) = sum_qu(i)/hw0(i)
    604       av_thvu(i) = sum_thvu(i)/hw0(i)
    605 c      av_thve = sum_thve/hw0
    606       av_dth(i) = sum_dth(i)/hw0(i)
    607       av_dq(i) = sum_dq(i)/hw0(i)
    608       av_rho(i) = sum_rho(i)/hw0(i)
    609       av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
    610       av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
    611 
    612       wape(i) = - rg*hw0(i)*(av_dth(i)
    613      $     + eps*(av_thu(i)*av_dq(i)+av_dth(i)*av_qu(i)+av_dth(i)*
    614      $     av_dq(i) ))/av_thvu(i)
    615       ENDDO
    616 C
    617 C 2.2 Prognostic variable update
    618 C ------------------------------
    619 C
    620 C Filter out bad wakes
    621 
    622       DO k = 1,klev
    623        DO i=1,klon
    624         IF ( wape(i) .LT. 0.) THEN
    625           deltatw(i,k) = 0.
    626           deltaqw(i,k) = 0.
    627           dth(i,k) = 0.
    628         ENDIF
    629        ENDDO
    630       ENDDO
    631 c
    632       DO i=1,klon
    633       IF ( wape(i) .LT. 0.) THEN
    634         wape(i) = 0.
    635         Cstar(i) = 0.
    636         hw(i) = hwmin
    637         sigmaw(i) = amax1(sigmad,sigd_con(i))
    638         fip(i) = 0.
    639         gwake(i) = .FALSE.
    640       ELSE
    641         Cstar(i) = stark*sqrt(2.*wape(i))
    642         gwake(i) = .TRUE.
    643       ENDIF
    644       ENDDO
    645 
    646 c
    647 c Check qx and qw positivity
    648 c --------------------------
    649       DO i = 1,klon
    650        q0_min(i)=min(  (qe(i,1)-sigmaw(i)*deltaqw(i,1)),
    651      $              (qe(i,1)+(1.-sigmaw(i))*deltaqw(i,1))  )
    652       ENDDO
    653       DO k = 2,klev
    654       DO i = 1,klon
    655         q1_min(i)=min(  (qe(i,k)-sigmaw(i)*deltaqw(i,k)),
    656      $              (qe(i,k)+(1.-sigmaw(i))*deltaqw(i,k))  )
    657         IF (q1_min(i).le.q0_min(i)) THEN
    658           q0_min(i)=q1_min(i)
    659         ENDIF
    660       ENDDO
    661       ENDDO
    662 c
    663       DO i = 1,klon
    664        OK_qx_qw(i) = q0_min(i) .GE. 0.
    665        alpha(i) = 1.
    666       ENDDO
    667 c
    668 CC -----------------------------------------------------------------
    669 C    Sub-time-stepping
    670 C    -----------------
    671 C
    672       nsub=10
    673       dtimesub=dtime/nsub
    674 c
    675 c------------------------------------------------------------
    676       DO isubstep = 1,nsub
    677 c------------------------------------------------------------
    678 c
    679 c wk_adv is the logical flag enabling wake evolution in the time advance loop
    680       DO i = 1,klon
    681        wk_adv(i) = OK_qx_qw(i) .AND. alpha(i) .GE. 1.
    682       ENDDO
    683 c
    684       DO i=1,klon
    685         IF (wk_adv(i)) THEN
    686         gfl(i) = 2.*sqrt(3.14*wdens*sigmaw(i))
    687         ENDIF
    688       ENDDO
    689       DO i=1,klon
    690         IF (wk_adv(i)) THEN
    691          d_sigmaw(i) = gfl(i)*Cstar(i)*dtimesub
    692 c        sigmaw(i) =sigmaw(i) + gfl(i)*Cstar(i)*dtimesub
    693 c        sigmaw(i) =min(sigmaw(i),0.99)     !!!!!!!!
    694 c        wdens = wdens0/(10.*sigmaw)
    695 c        sigmaw =max(sigmaw,sigd_con)
    696 c        sigmaw =max(sigmaw,sigmad)
    697         ENDIF
    698       ENDDO
    699 C
    700 C
    701 c calcul de la difference de vitesse verticale poche - zone non perturbee
    702 cIM 060208 differences par rapport au code initial; init. a 0 dp_deltomg
    703 cIM 060208 et omg sur les niveaux de 1 a klev+1, alors que avant l'on definit
    704 cIM 060208 au niveau k=1..?
    705       DO k= 1,klev
    706       DO i = 1,klon
    707         dp_deltomg(i,k)=0.
    708       ENDDO
    709       ENDDO
    710       DO k= 1,klev+1
    711       DO i = 1,klon
    712         omg(i,k)=0.
    713       ENDDO
    714       ENDDO
    715 c
    716       DO i=1,klon
    717         IF (wk_adv(i)) THEN
    718         z(i)= 0.
    719         omg(i,1) = 0.
    720         dp_deltomg(i,1) = -(gfl(i)*Cstar(i))/(sigmaw(i) * (1-sigmaw(i)))
    721         ENDIF
    722       ENDDO
    723 c
    724       DO k= 2,klev
    725       DO i = 1,klon
    726        IF( wk_adv(i) .AND. k .LE. ktop(i)) THEN
    727           dz(i) = -(ph(i,k)-ph(i,k-1))/(rho(i,k-1)*rg)
    728           z(i) = z(i)+dz(i)
    729           dp_deltomg(i,k)= dp_deltomg(i,1)
    730           omg(i,k)= dp_deltomg(i,1)*z(i)
    731        ENDIF
    732       ENDDO
    733       ENDDO
    734 c
    735       DO i = 1,klon
    736         IF (wk_adv(i)) THEN
    737         dztop(i)=-(ptop(i)-ph(i,ktop(i)))/(rho(i,ktop(i))*rg)
    738         ztop(i) = z(i)+dztop(i)
    739         omgtop(i)=dp_deltomg(i,1)*ztop(i)
    740         ENDIF
    741       ENDDO
    742 c
    743 c        -----------------
    744 c        From m/s to Pa/s
    745 c        -----------------
    746 c
    747        DO i=1,klon
    748         IF (wk_adv(i)) THEN
    749         omgtop(i) = -rho(i,ktop(i))*rg*omgtop(i)
    750         dp_deltomg(i,1) = omgtop(i)/(ptop(i)-ph(i,1))
    751         ENDIF
    752        ENDDO
    753 c
    754       DO k= 1,klev
    755       DO i = 1,klon
    756        IF( wk_adv(i) .AND. k .LE. ktop(i)) THEN
    757           omg(i,k) = - rho(i,k)*rg*omg(i,k)
    758           dp_deltomg(i,k) = dp_deltomg(i,1)
    759        ENDIF
    760       ENDDO
    761       ENDDO
    762 c
    763 c   raccordement lineaire de omg de ptop a pupper
    764 
    765       DO i=1,klon
    766       IF ( wk_adv(i) .AND. kupper(i) .GT. ktop(i)) THEN
    767         omg(i,kupper(i)+1) = - Rg*amdwn(i,kupper(i)+1)/sigmaw(i)
    768      $                + Rg*amup(i,kupper(i)+1)/(1.-sigmaw(i))
    769         dp_deltomg(i,kupper(i)) = (omgtop(i)-omg(i,kupper(i)+1))/
    770      $                     (ptop(i)-pupper(i))
    771       ENDIF
    772       ENDDO
    773 c
    774       DO k= 1,klev
    775       DO i = 1,klon
    776        IF( wk_adv(i) .AND. k .GT. ktop(i) .AND. k .LE. kupper(i)) THEN
    777           dp_deltomg(i,k) = dp_deltomg(i,kupper(i))
    778           omg(i,k) = omgtop(i)+(ph(i,k)-ptop(i))*dp_deltomg(i,kupper(i))
    779        ENDIF
    780       ENDDO
    781       ENDDO
    782 c
    783 c
    784 c--    Compute wake average vertical velocity omgbw
    785 c
    786 c
    787       DO k = 1,klev+1
    788       DO i=1,klon
    789         IF ( wk_adv(i)) THEN
    790         omgbw(i,k) = omgb(i,k)+(1.-sigmaw(i))*omg(i,k)
    791         ENDIF
    792       ENDDO
    793       ENDDO
    794 c--    and its vertical gradient dp_omgbw
    795 c
    796       DO k = 1,klev
    797       DO i=1,klon
    798         IF ( wk_adv(i)) THEN
    799         dp_omgbw(i,k) = (omgbw(i,k+1)-omgbw(i,k))/(ph(i,k+1)-ph(i,k))
    800         ENDIF
    801       ENDDO
    802       ENDDO
    803 C
    804 c--    Upstream coefficients for omgb velocity
    805 c--    (alpha_up(k) is the coefficient of the value at level k)
    806 c--    (1-alpha_up(k) is the coefficient of the value at level k-1)
    807       DO k = 1,klev
    808       DO i=1,klon
    809         IF ( wk_adv(i)) THEN
    810          alpha_up(i,k) = 0.
    811          IF (omgb(i,k) .GT. 0.) alpha_up(i,k) = 1.
    812         ENDIF
    813       ENDDO
    814       ENDDO
    815 
    816 c  Matrix expressing [The,deltatw] from  [Th1,Th2]
    817 
    818       DO i=1,klon
    819         IF ( wk_adv(i)) THEN
    820          RRe1(i) = 1.-sigmaw(i)
    821          RRe2(i) = sigmaw(i)
    822         ENDIF
    823       ENDDO
    824       RRd1 = -1.
    825       RRd2 = 1.
    826 c
    827 c--    Get [Th1,Th2], dth and [q1,q2]
    828 c
    829       DO k= 1,klev
    830       DO i = 1,klon
    831        IF( wk_adv(i) .AND. k .LE. kupper(i)+1) THEN
    832         dth(i,k) = deltatw(i,k)/ppi(i,k)
    833         Th1(i,k) = the(i,k) - sigmaw(i)     *dth(i,k)   ! undisturbed area
    834         Th2(i,k) = the(i,k) + (1.-sigmaw(i))*dth(i,k)   ! wake
    835         q1(i,k) = qe(i,k) - sigmaw(i)     *deltaqw(i,k) ! undisturbed area
    836         q2(i,k) = qe(i,k) + (1.-sigmaw(i))*deltaqw(i,k) ! wake
    837         T1(i,k) = te(i,k) - sigmaw(i)*deltatw(i,k)! undisturb itlmd
    838        ENDIF
    839       ENDDO
    840       ENDDO
    841 
    842       DO i=1,klon
    843        D_Th1(i,1) = 0.   !!!itlmd : ne pas mettre if wk_adv cf nrlmd?
    844        D_Th2(i,1) = 0.
    845        D_dth(i,1) = 0.
    846        D_q1(i,1) = 0.
    847        D_q2(i,1) = 0.
    848        D_dq(i,1) = 0.
    849       ENDDO
    850 
    851       DO k= 2,klev
    852       DO i = 1,klon
    853        IF( wk_adv(i) .AND. k .LE. kupper(i)+1) THEN
    854         D_Th1(i,k) = Th1(i,k-1)-Th1(i,k)
    855         D_Th2(i,k) = Th2(i,k-1)-Th2(i,k)
    856         D_dth(i,k) = dth(i,k-1)-dth(i,k)
    857         D_q1(i,k) = q1(i,k-1)-q1(i,k)
    858         D_q2(i,k) = q2(i,k-1)-q2(i,k)
    859         D_dq(i,k) = deltaqw(i,k-1)-deltaqw(i,k)
    860        ENDIF
    861       ENDDO
    862       ENDDO
    863 
    864       DO i=1,klon
    865         IF( wk_adv(i)) THEN
    866          omgbdth(i,1) = 0.
    867          omgbdq(i,1) = 0.
    868         ENDIF
    869       ENDDO
    870 
    871       DO k= 2,klev
    872       DO i = 1,klon
    873        IF( wk_adv(i) .AND. k .LE. kupper(i)+1) THEN  !   loop on interfaces
    874         omgbdth(i,k) = omgb(i,k)*(    dth(i,k-1) -     dth(i,k))
    875         omgbdq(i,k)  = omgb(i,k)*(deltaqw(i,k-1) - deltaqw(i,k))
    876        ENDIF
    877       ENDDO
    878       ENDDO
    879 c
    880 c-----------------------------------------------------------------
    881       DO k= 1,klev
    882       DO i = 1,klon
    883        IF( wk_adv(i) .AND. k .LE. kupper(i)-1) THEN
    884 c-----------------------------------------------------------------
    885 c
    886 c   Compute redistribution (advective) term
    887 c
    888          d_deltatw(i,k) =
    889      $             dtimesub/(Ph(i,k)-Ph(i,k+1))*(
    890      $       RRd1*omg(i,k  )*sigmaw(i)     *D_Th1(i,k)
    891      $      -RRd2*omg(i,k+1)*(1.-sigmaw(i))*D_Th2(i,k+1)
    892      $      -(1.-alpha_up(i,k))*omgbdth(i,k) - alpha_up(i,k+1)*
    893      $      omgbdth(i,k+1))*ppi(i,k)
    894 c         print*,'d_deltatw=',d_deltatw(i,k)
    895 c
    896          d_deltaqw(i,k) =
    897      $             dtimesub/(Ph(i,k)-Ph(i,k+1))*(
    898      $       RRd1*omg(i,k  )*sigmaw(i)     *D_q1(i,k)
    899      $      -RRd2*omg(i,k+1)*(1.-sigmaw(i))*D_q2(i,k+1)
    900      $      -(1.-alpha_up(i,k))*omgbdq(i,k) - alpha_up(i,k+1)*
    901      $      omgbdq(i,k+1))
    902 c         print*,'d_deltaqw=',d_deltaqw(i,k)
    903 c
    904 c   and increment large scale tendencies
    905 c
    906 
    907 c
    908 C
    909 CC -----------------------------------------------------------------
    910          d_te(i,k) =  dtimesub*(
    911      $        ( RRe1(i)*omg(i,k  )*sigmaw(i)     *D_Th1(i,k)
    912      $         -RRe2(i)*omg(i,k+1)*(1.-sigmaw(i))*D_Th2(i,k+1) )
    913      $               /(Ph(i,k)-Ph(i,k+1))
    914      $         -sigmaw(i)*(1.-sigmaw(i))*dth(i,k)*
    915      $            (omg(i,k)-omg(i,k+1))/(Ph(i,k)-Ph(i,k+1)) !instead of dp_deltomg(i,k)
    916      $                      )*ppi(i,k)
    917 c
    918          d_qe(i,k) =  dtimesub*(
    919      $        ( RRe1(i)*omg(i,k  )*sigmaw(i)     *D_q1(i,k)
    920      $         -RRe2(i)*omg(i,k+1)*(1.-sigmaw(i))*D_q2(i,k+1) )
    921      $               /(Ph(i,k)-Ph(i,k+1))
    922      $         -sigmaw(i)*(1.-sigmaw(i))*deltaqw(i,k)*
    923      $           (omg(i,k)-omg(i,k+1))/(Ph(i,k)-Ph(i,k+1))!instead of dp_deltomg(i,k)
    924      $                      )
    925         ELSE IF(wk_adv(i) .AND. k .EQ. kupper(i)) THEN ! corr pour conserver l'eau
    926 
    927          d_te(i,k) =  dtimesub*(
    928      $        ( RRe1(i)*omg(i,k  )*sigmaw(i)     *D_Th1(i,k))
    929      $               /(Ph(i,k)-Ph(i,k+1))
    930      $                      )*ppi(i,k)
    931 
    932          d_qe(i,k) =  dtimesub*(
    933      $        ( RRe1(i)*omg(i,k  )*sigmaw(i)     *D_q1(i,k))
    934      $               /(Ph(i,k)-Ph(i,k+1))
    935      $                      )
    936        ENDIF
    937 
    938 c-------------------------------------------------------------------
    939       ENDDO
    940       ENDDO
    941 c------------------------------------------------------------------
    942 C
    943 C   Increment state variables
    944 
    945       DO k= 1,klev
    946       DO i = 1,klon
    947        IF( wk_adv(i) .AND. k .LE. kupper(i)-1) THEN
    948 c
    949 c Coefficient de répartition
    950 
    951         Crep(i,k)=Crep_sol*(ph(i,kupper(i))-ph(i,k))/(ph(i,kupper(i))
    952      $          -ph(i,1))
    953         Crep(i,k)=Crep(i,k)+Crep_upper*(ph(i,1)-ph(i,k))/(p(i,1)-
    954      $          ph(i,kupper(i)))
    955        
    956 
    957 c Reintroduce compensating subsidence term.
    958 
    959 c        dtKE(k)=(dtdwn(k)*Crep(k))/sigmaw
    960 c        dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k))
    961 c     .                   /(1-sigmaw)
    962 c        dqKE(k)=(dqdwn(k)*Crep(k))/sigmaw
    963 c        dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k))
    964 c     .                   /(1-sigmaw)
    965 c
    966 c        dtKE(k)=(dtdwn(k)*Crep(k)+(1-Crep(k))*dta(k))/sigmaw
    967 c        dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k)*Crep(k))
    968 c     .                   /(1-sigmaw)
    969 c        dqKE(k)=(dqdwn(k)*Crep(k)+(1-Crep(k))*dqa(k))/sigmaw
    970 c        dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k)*Crep(k))
    971 c     .                   /(1-sigmaw)
    972 
    973         dtKE(i,k)=(dtdwn(i,k)/sigmaw(i) - dta(i,k)/(1.-sigmaw(i)))
    974         dqKE(i,k)=(dqdwn(i,k)/sigmaw(i) - dqa(i,k)/(1.-sigmaw(i)))
    975 c        print*,'dtKE=',dtKE(k)
    976 c        print*,'dqKE=',dqKE(k)
    977 c
    978         dtPBL(i,k)=(wdtPBL(i,k)/sigmaw(i) - udtPBL(i,k)/(1.-sigmaw(i)))
    979         dqPBL(i,k)=(wdqPBL(i,k)/sigmaw(i) - udqPBL(i,k)/(1.-sigmaw(i)))
    980 c
    981         spread(i,k) = (1.-sigmaw(i))*dp_deltomg(i,k)+gfl(i)*Cstar(i)/
    982      $  sigmaw(i)
    983 
    984 
    985 c ajout d'un effet onde de gravité -Tgw(k)*deltatw(k) 03/02/06 YU Jingmei
    986 
    987         d_deltat_gw(i,k)=d_deltat_gw(i,k)-Tgw(i,k)*deltatw(i,k)*
    988      $  dtimesub
    989         ff(i)=d_deltatw(i,k)/dtimesub
    990 
    991 c Sans GW
    992 c
    993 c        deltatw(k)=deltatw(k)+dtimesub*(ff+dtKE(k)-spread(k)*deltatw(k))
    994 c
    995 c GW formule 1
    996 c
    997 c        deltatw(k) = deltatw(k)+dtimesub*
    998 c     $         (ff+dtKE(k) - spread(k)*deltatw(k)-Tgw(k)*deltatw(k))
    999 c
    1000 c GW formule 2
    1001 
    1002         IF (dtimesub*Tgw(i,k).lt.1.e-10) THEN
    1003           d_deltatw(i,k) = dtimesub*
    1004      $          (ff(i)+dtKE(i,k)+dtPBL(i,k)
    1005      $          - spread(i,k)*deltatw(i,k)-Tgw(i,k)*deltatw(i,k))
    1006         ELSE
    1007            d_deltatw(i,k) = 1/Tgw(i,k)*(1-exp(-dtimesub*
    1008      $          Tgw(i,k)))*
    1009      $          (ff(i)+dtKE(i,k)+dtPBL(i,k)
    1010      $          - spread(i,k)*deltatw(i,k)-Tgw(i,k)*deltatw(i,k))
    1011         ENDIF
    1012 
    1013         dth(i,k) = deltatw(i,k)/ppi(i,k)
    1014 
    1015         gg(i)=d_deltaqw(i,k)/dtimesub
    1016 
    1017        d_deltaqw(i,k) = dtimesub*(gg(i)+ dqKE(i,k)+dqPBL(i,k)
    1018      $                            - spread(i,k)*deltaqw(i,k))
    1019 
    1020        d_deltatw2(i,k)=d_deltatw2(i,k)+d_deltatw(i,k)
    1021        d_deltaqw2(i,k)=d_deltaqw2(i,k)+d_deltaqw(i,k)
    1022        ENDIF
    1023       ENDDO
    1024       ENDDO
    1025 
    1026 C
    1027 C   Scale tendencies so that water vapour remains positive in w and x.
    1028 C
    1029       call wake_vec_modulation(klon,klev,wk_adv,qe,d_qe,deltaqw,
    1030      $                d_deltaqw,sigmaw,d_sigmaw,alpha)
    1031 c
    1032       DO k = 1,klev
    1033       DO i = 1,klon
    1034        IF( wk_adv(i) .AND. k .LE. kupper(i)) THEN
    1035         d_te(i,k)=alpha(i)*d_te(i,k)
    1036         d_qe(i,k)=alpha(i)*d_qe(i,k)
    1037         d_deltatw(i,k)=alpha(i)*d_deltatw(i,k)
    1038         d_deltaqw(i,k)=alpha(i)*d_deltaqw(i,k)
    1039         d_deltat_gw(i,k)=alpha(i)*d_deltat_gw(i,k)
    1040        ENDIF
    1041       ENDDO
    1042       ENDDO
    1043       DO i = 1,klon
    1044        IF( wk_adv(i)) THEN
    1045         d_sigmaw(i)=alpha(i)*d_sigmaw(i)
    1046        ENDIF
    1047       ENDDO
    1048 
    1049 C   Update large scale variables and wake variables
    1050 cIM 060208 manque DO i + remplace DO k=1,kupper(i)
    1051 cIM 060208     DO k = 1,kupper(i)
    1052       DO k= 1,klev
    1053       DO i = 1,klon
    1054        IF( wk_adv(i) .AND. k .LE. kupper(i)) THEN
    1055         dtls(i,k)=dtls(i,k)+d_te(i,k)
    1056         dqls(i,k)=dqls(i,k)+d_qe(i,k)
    1057        ENDIF
    1058       ENDDO
    1059       ENDDO
    1060       DO k= 1,klev
    1061       DO i = 1,klon
    1062        IF( wk_adv(i) .AND. k .LE. kupper(i)) THEN
    1063         te(i,k) = te0(i,k) + dtls(i,k)
    1064         qe(i,k) = qe0(i,k) + dqls(i,k)
    1065         the(i,k) = te(i,k)/ppi(i,k)
    1066         deltatw(i,k) = deltatw(i,k)+d_deltatw(i,k)
    1067         deltaqw(i,k) = deltaqw(i,k)+d_deltaqw(i,k)
    1068         dth(i,k) = deltatw(i,k)/ppi(i,k)
    1069        ENDIF
    1070       ENDDO
    1071       ENDDO
    1072       DO i = 1,klon
    1073        IF( wk_adv(i)) THEN
    1074         sigmaw(i) = sigmaw(i)+d_sigmaw(i)
    1075        ENDIF
    1076       ENDDO
    1077 c
    1078 C
    1079 c     Determine Ptop from buoyancy integral
    1080 c     ---------------------------------------
    1081 c
    1082 c-     1/ Pressure of the level where dth changes sign.
    1083 c
    1084       DO i=1,klon
    1085        IF ( wk_adv(i)) THEN
    1086         Ptop_provis(i)=ph(i,1)
    1087        ENDIF
    1088       ENDDO
    1089 c
    1090       DO k= 2,klev
    1091       DO i=1,klon
    1092         IF ( wk_adv(i) .AND.
    1093      $       Ptop_provis(i) .EQ. ph(i,1) .AND.
    1094      $      dth(i,k) .GT. -delta_t_min .and.
    1095      $      dth(i,k-1).LT. -delta_t_min) THEN
    1096           Ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)
    1097      $          - (dth(i,k-1)+delta_t_min)*p(i,k)) /(dth(i,k)
    1098      $          - dth(i,k-1))
    1099         ENDIF
    1100       ENDDO
    1101       ENDDO
    1102 c
    1103 c-     2/ dth integral
    1104 c
    1105       DO i=1,klon
    1106        sum_dth(i) = 0.
    1107        dthmin(i) = -delta_t_min
    1108        z(i) = 0.
    1109       ENDDO
    1110 
    1111       DO k = 1,klev
    1112       DO i=1,klon
    1113        IF ( wk_adv(i)) THEN
    1114         dz(i) = -(amax1(ph(i,k+1),Ptop_provis(i))-Ph(i,k))/(rho(i,k)*rg)
    1115         IF (dz(i) .gt. 0) THEN
    1116          z(i) = z(i)+dz(i)
    1117          sum_dth(i) = sum_dth(i) + dth(i,k)*dz(i)
    1118          dthmin(i) = amin1(dthmin(i),dth(i,k))
    1119         ENDIF
    1120        ENDIF
    1121       ENDDO
    1122       ENDDO
    1123 c
    1124 c-     3/ height of triangle with area= sum_dth and base = dthmin
    1125 
    1126       DO i=1,klon
    1127        IF ( wk_adv(i)) THEN
    1128          hw(i) = 2.*sum_dth(i)/amin1(dthmin(i),-0.5)
    1129          hw(i) = amax1(hwmin,hw(i))
    1130        ENDIF
    1131       ENDDO
    1132 c
    1133 c-     4/ now, get Ptop
    1134 c
    1135       DO i=1,klon
    1136        ktop(i) = 0
    1137        z(i)=0.
    1138       ENDDO
    1139 c
    1140       DO k = 1,klev
    1141       DO i=1,klon
    1142        IF ( wk_adv(i)) THEN
    1143         dz(i) = amin1(-(ph(i,k+1)-Ph(i,k))/(rho(i,k)*rg),hw(i)-z(i))
    1144         IF (dz(i) .gt. 0) THEN
    1145          z(i) = z(i)+dz(i)
    1146          Ptop(i) = Ph(i,k)-rho(i,k)*rg*dz(i)
    1147          ktop(i) = k
    1148         ENDIF
    1149        ENDIF
    1150       ENDDO
    1151       ENDDO
    1152 c
    1153 c      4.5/Correct ktop and ptop
    1154 c
    1155       DO i=1,klon
    1156        IF ( wk_adv(i)) THEN
    1157         Ptop_new(i)=ptop(i)
    1158        ENDIF
    1159       ENDDO
    1160 c
    1161       DO k= klev,2,-1
    1162       DO i=1,klon
    1163 cIM v3JYG; IF (k .GE. ktop(i)
    1164        IF ( wk_adv(i) .AND.
    1165      $      k .LE. ktop(i) .AND.
    1166      $      ptop_new(i) .EQ. ptop(i) .AND.
    1167      $      dth(i,k) .GT. -delta_t_min .and.
    1168      $      dth(i,k-1).LT. -delta_t_min) THEN
    1169           Ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)
    1170      $          - (dth(i,k-1)+delta_t_min)*p(i,k)) /(dth(i,k)
    1171      $          - dth(i,k-1))
    1172         ENDIF
    1173       ENDDO
    1174       ENDDO
    1175 c
    1176 c
    1177       DO i=1,klon
    1178        IF ( wk_adv(i)) THEN
    1179         ptop(i) = ptop_new(i)
    1180        ENDIF
    1181       ENDDO
    1182 
    1183       DO k=klev,1,-1
    1184       DO i=1,klon
    1185         IF (ph(i,k+1) .LT. ptop(i)) ktop(i)=k
    1186       ENDDO
    1187       ENDDO
    1188 c
    1189 c      5/ Set deltatw & deltaqw to 0 above kupper
    1190 c
    1191       DO k = 1,klev
    1192       DO i=1,klon
    1193         IF ( wk_adv(i) .AND. k .GE. kupper(i)) THEN
    1194          deltatw(i,k) = 0.
    1195          deltaqw(i,k) = 0.
    1196         ENDIF
    1197       ENDDO
    1198       ENDDO
    1199 c
    1200 C
    1201 c-------------Cstar computation---------------------------------
    1202       DO i=1, klon
    1203       sum_thu(i) = 0.
    1204       sum_tu(i) = 0.
    1205       sum_qu(i) = 0.
    1206       sum_thvu(i) = 0.
    1207       sum_dth(i) = 0.
    1208       sum_dq(i) = 0.
    1209       sum_rho(i) = 0.
    1210       sum_dtdwn(i) = 0.
    1211       sum_dqdwn(i) = 0.
    1212 
    1213       av_thu(i) = 0.
    1214       av_tu(i) =0.
    1215       av_qu(i) =0.
    1216       av_thvu(i) = 0.
    1217       av_dth(i) = 0.
    1218       av_dq(i) = 0.
    1219       av_rho(i) =0.
    1220       av_dtdwn(i) =0.
    1221       av_dqdwn(i) = 0.
    1222       ENDDO
    1223 C
    1224 C Integrals (and wake top level number)
    1225 C --------------------------------------
    1226 C
    1227 C Initialize sum_thvu to 1st level virt. pot. temp.
    1228 
    1229       DO i=1,klon
    1230       z(i) = 1.
    1231       dz(i) = 1.
    1232       sum_thvu(i) =  thu(i,1)*(1.+eps*qu(i,1))*dz(i)
    1233       sum_dth(i) = 0.
    1234       ENDDO
    1235 
    1236       DO k = 1,klev
    1237       DO i=1,klon
    1238         dz(i) = -(max(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
    1239         IF (dz(i) .GT. 0) THEN
    1240          z(i) = z(i)+dz(i)
    1241          sum_thu(i) = sum_thu(i) + th1(i,k)*dz(i)
    1242          sum_tu(i) = sum_tu(i) + t1(i,k)*dz(i)
    1243          sum_qu(i) = sum_qu(i) + q1(i,k)*dz(i)
    1244          sum_thvu(i) = sum_thvu(i) + th1(i,k)*(1.+eps*q1(i,k))*dz(i)!itlmd
    1245 
    1246          sum_dth(i) = sum_dth(i) + dth(i,k)*dz(i)
    1247          sum_dq(i) = sum_dq(i) + deltaqw(i,k)*dz(i)
    1248          sum_rho(i) = sum_rho(i) + rhow(i,k)*dz(i)
    1249          sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i,k)*dz(i)
    1250          sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i,k)*dz(i)
    1251         ENDIF
    1252       ENDDO
    1253       ENDDO
    1254 c
    1255       DO i=1,klon
    1256         hw0(i) = z(i)
    1257       ENDDO
    1258 c
    1259 C
    1260 C - WAPE and mean forcing computation
    1261 C ---------------------------------------
    1262 C
    1263 C ---------------------------------------
    1264 C
    1265 C Means
    1266 
    1267       DO i=1,klon
    1268        av_thu(i) = sum_thu(i)/hw0(i)
    1269        av_tu(i) = sum_tu(i)/hw0(i)
    1270        av_qu(i) = sum_qu(i)/hw0(i)
    1271        av_thvu(i) = sum_thvu(i)/hw0(i)
    1272        av_dth(i) = sum_dth(i)/hw0(i)
    1273        av_dq(i) = sum_dq(i)/hw0(i)
    1274        av_rho(i) = sum_rho(i)/hw0(i)
    1275        av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
    1276        av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
    1277 c
    1278        wape(i) = - rg*hw0(i)*(av_dth(i)
    1279      $     + eps*(av_thu(i)*av_dq(i)+av_dth(i)*av_qu(i)+av_dth(i)*
    1280      $     av_dq(i) ))/av_thvu(i)
    1281       ENDDO
    1282 C
    1283 C Filter out bad wakes
    1284 
    1285       DO k = 1,klev
    1286        DO i=1,klon
    1287         IF ( wape(i) .LT. 0.) THEN
    1288           deltatw(i,k) = 0.
    1289           deltaqw(i,k) = 0.
    1290           dth(i,k) = 0.
    1291         ENDIF
    1292        ENDDO
    1293       ENDDO
    1294 c
    1295       DO i=1,klon
    1296       IF ( wape(i) .LT. 0.) THEN
    1297         wape(i) = 0.
    1298         Cstar(i) = 0.
    1299         hw(i) = hwmin
    1300         sigmaw(i) = max(sigmad,sigd_con(i))
    1301         fip(i) = 0.
    1302         gwake(i) = .FALSE.
    1303       ELSE
    1304         Cstar(i) = stark*sqrt(2.*wape(i))
    1305         gwake(i) = .TRUE.
    1306       ENDIF
    1307       ENDDO
    1308 
    1309        ENDDO      ! end sub-timestep loop
    1310 C
    1311 C -----------------------------------------------------------------
    1312 c   Get back to tendencies per second
    1313 c
    1314       DO k = 1,klev
    1315       DO i=1,klon
    1316        IF ( wk_adv(i) .AND. k .LE. kupper(i)) THEN  !! corr conservation eau
    1317          dtls(i,k) = dtls(i,k)/dtime
    1318          dqls(i,k) = dqls(i,k)/dtime
    1319          d_deltatw2(i,k)=d_deltatw2(i,k)/dtime
    1320          d_deltaqw2(i,k)=d_deltaqw2(i,k)/dtime
    1321          d_deltat_gw(i,k) = d_deltat_gw(i,k)/dtime
    1322         ENDIF
    1323       ENDDO
    1324       ENDDO
    1325 c
    1326 c
    1327 c----------------------------------------------------------
    1328 c   Determine wake final state; recompute wape, cstar, ktop;
    1329 c   filter out bad wakes.
    1330 c----------------------------------------------------------
    1331 c
    1332 C 2.1 - Undisturbed area and Wake integrals
    1333 C ---------------------------------------------------------
    1334 
    1335       DO i=1,klon
    1336         z(i) = 0.
    1337         sum_thu(i) = 0.
    1338         sum_tu(i) = 0.
    1339         sum_qu(i) = 0.
    1340         sum_thvu(i) = 0.
    1341         sum_dth(i) = 0.
    1342         sum_dq(i) = 0.
    1343         sum_rho(i) = 0.
    1344         sum_dtdwn(i) = 0.
    1345         sum_dqdwn(i) = 0.
    1346 
    1347         av_thu(i) = 0.
    1348         av_tu(i) =0.
    1349         av_qu(i) =0.
    1350         av_thvu(i) = 0.
    1351         av_dth(i) = 0.
    1352         av_dq(i) = 0.
    1353         av_rho(i) =0.
    1354         av_dtdwn(i) =0.
    1355         av_dqdwn(i) = 0.
    1356       ENDDO
    1357 C Potential temperatures and humidity
    1358 c----------------------------------------------------------
    1359 
    1360       DO k =1,klev
    1361       DO i=1,klon
    1362        IF ( wk_adv(i)) THEN
    1363         rho(i,k) = p(i,k)/(rd*te(i,k))
    1364         IF(k .eq. 1) THEN
    1365           rhoh(i,k) = ph(i,k)/(rd*te(i,k))
    1366           zhh(i,k)=0
    1367         ELSE
    1368           rhoh(i,k) = ph(i,k)*2./(rd*(te(i,k)+te(i,k-1)))
    1369           zhh(i,k)=(ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*RG)+zhh(i,k-1)
    1370         ENDIF
    1371         the(i,k) = te(i,k)/ppi(i,k)
    1372         thu(i,k) = (te(i,k) - deltatw(i,k)*sigmaw(i))/ppi(i,k)
    1373         tu(i,k) = te(i,k) - deltatw(i,k)*sigmaw(i)
    1374         qu(i,k)  =  qe(i,k) - deltaqw(i,k)*sigmaw(i)
    1375         rhow(i,k) = p(i,k)/(rd*(tu(i,k)+deltatw(i,k)))
    1376         dth(i,k) = deltatw(i,k)/ppi(i,k)
    1377        ENDIF
    1378       ENDDO
    1379       ENDDO
    1380 
    1381 C Integrals (and wake top level number)
    1382 C -----------------------------------------------------------
    1383 
    1384 C Initialize sum_thvu to 1st level virt. pot. temp.
    1385 
    1386       DO i=1,klon
    1387        IF ( wk_adv(i)) THEN
    1388         z(i) = 1.
    1389         dz(i) = 1.
    1390         sum_thvu(i) =  thu(i,1)*(1.+eps*qu(i,1))*dz(i)
    1391         sum_dth(i) = 0.
    1392       ENDIF
    1393       ENDDO
    1394 
    1395       DO k = 1,klev
    1396       DO i=1,klon
    1397        IF ( wk_adv(i)) THEN
    1398         dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
    1399         IF (dz(i) .GT. 0) THEN
    1400          z(i) = z(i)+dz(i)
    1401          sum_thu(i) = sum_thu(i) + thu(i,k)*dz(i)
    1402          sum_tu(i) = sum_tu(i) + tu(i,k)*dz(i)
    1403          sum_qu(i) = sum_qu(i) + qu(i,k)*dz(i)
    1404          sum_thvu(i) = sum_thvu(i) + thu(i,k)*(1.+eps*qu(i,k))*dz(i)
    1405          sum_dth(i) = sum_dth(i) + dth(i,k)*dz(i)
    1406          sum_dq(i) = sum_dq(i) + deltaqw(i,k)*dz(i)
    1407          sum_rho(i) = sum_rho(i) + rhow(i,k)*dz(i)
    1408          sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i,k)*dz(i)
    1409          sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i,k)*dz(i)
    1410         ENDIF
    1411        ENDIF
    1412       ENDDO
    1413       ENDDO
    1414 c
    1415       DO i=1,klon
    1416        IF ( wk_adv(i)) THEN
    1417         hw0(i) = z(i)
    1418        ENDIF
    1419       ENDDO
    1420 c
    1421 C - WAPE and mean forcing computation
    1422 C-------------------------------------------------------------
    1423 
    1424 C Means
    1425 
    1426       DO i=1, klon
    1427        IF ( wk_adv(i)) THEN
    1428         av_thu(i) = sum_thu(i)/hw0(i)
    1429         av_tu(i) = sum_tu(i)/hw0(i)
    1430         av_qu(i) = sum_qu(i)/hw0(i)
    1431         av_thvu(i) = sum_thvu(i)/hw0(i)
    1432         av_dth(i) = sum_dth(i)/hw0(i)
    1433         av_dq(i) = sum_dq(i)/hw0(i)
    1434         av_rho(i) = sum_rho(i)/hw0(i)
    1435         av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
    1436         av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
    1437 
    1438         wape2(i) = - rg*hw0(i)*(av_dth(i)
    1439      $     + eps*(av_thu(i)*av_dq(i)+av_dth(i)*av_qu(i)+
    1440      $     av_dth(i)*av_dq(i) ))/av_thvu(i)
    1441        ENDIF
    1442       ENDDO
    1443 
    1444 C Prognostic variable update
    1445 C ------------------------------------------------------------
    1446 
    1447 C Filter out bad wakes
    1448 c
    1449       DO k = 1,klev
    1450       DO i=1,klon
    1451         IF ( wk_adv(i) .AND. wape2(i) .LT. 0.) THEN
    1452           deltatw(i,k) = 0.
    1453           deltaqw(i,k) = 0.
    1454           dth(i,k) = 0.
    1455         ENDIF
    1456       ENDDO
    1457       ENDDO
    1458 c
    1459 
    1460       DO i=1, klon
    1461        IF ( wk_adv(i)) THEN
    1462        IF ( wape2(i) .LT. 0.) THEN
    1463         wape2(i) = 0.
    1464         Cstar2(i) = 0.
    1465         hw(i) = hwmin
    1466         sigmaw(i) = amax1(sigmad,sigd_con(i))
    1467         fip(i) = 0.
    1468         gwake(i) = .FALSE.
    1469       ELSE
    1470         if(prt_level.ge.10) print*,'wape2>0'
    1471         Cstar2(i) = stark*sqrt(2.*wape2(i))
    1472         gwake(i) = .TRUE.
    1473       ENDIF
    1474       ENDIF
    1475       ENDDO
    1476 c
    1477       DO i=1, klon
    1478        IF ( wk_adv(i)) THEN
    1479         ktopw(i) = ktop(i)
    1480        ENDIF
    1481       ENDDO
    1482 c
    1483       DO i=1, klon
    1484        IF ( wk_adv(i)) THEN
    1485        IF (ktopw(i) .gt. 0 .and. gwake(i)) then
    1486 
    1487 Cjyg1     Utilisation d'un h_efficace constant ( ~ feeding layer)
    1488 ccc       heff = 600.
    1489 C      Utilisation de la hauteur hw
    1490 cc       heff = 0.7*hw
    1491        heff(i) = hw(i)
    1492 
    1493        FIP(i) = 0.5*rho(i,ktopw(i))*Cstar2(i)**3*heff(i)*2*
    1494      $      sqrt(sigmaw(i)*wdens*3.14)
    1495        FIP(i) = alpk * FIP(i)
    1496 Cjyg2
    1497        ELSE
    1498          FIP(i) = 0.
    1499        ENDIF
    1500        ENDIF
    1501       ENDDO
    1502 c
    1503 C   Limitation de sigmaw
    1504 c
    1505 C   sécurité : si le wake occuppe plus de 90 % de la surface de la maille,
    1506 C              alors il disparait en se mélangeant à la partie undisturbed
    1507 c
    1508       sigmaw_max = 0.9
    1509       DO k = 1,klev
    1510        DO i=1, klon
    1511 c correction NICOLAS     $     ((wape(i).ge.wape2(i)).and.(wape2(i).le.1.0))) THEN
    1512 !         print*,'wape wape2 ktopw OK_qx_qw =',
    1513 !     $           wape(i),wape2(i),ktopw(i),OK_qx_qw(i)
    1514          IF ((sigmaw(i).GT.sigmaw_max).or.
    1515      $     ((wape(i).ge.wape2(i)).and.(wape2(i).le.1.0)).or.
    1516      $      (ktopw(i).le.2) .OR.
    1517      $     .not. OK_qx_qw(i)) THEN
    1518 cIM cf NR/JYG 251108  $     ((wape(i).ge.wape2(i)).and.(wape2(i).le.1.0))) THEN
    1519 ccc      IF (sigmaw(i).GT.0.9) THEN
    1520           dtls(i,k) = 0.
    1521           dqls(i,k) = 0.
    1522           deltatw(i,k) = 0.
    1523           deltaqw(i,k) = 0.
    1524         ENDIF
    1525        ENDDO
    1526       ENDDO
    1527 c
    1528       DO i=1, klon
    1529          IF ( (sigmaw(i).GT.sigmaw_max).or.
    1530      $      ((wape(i).ge.wape2(i)).and.(wape2(i).le.1.0)).or.
    1531      $      (ktopw(i).le.2) .OR.
    1532      $     .not. OK_qx_qw(i)) THEN
    1533 ! correction NICOLAS     $     ((wape(i).ge.wape2(i)).and.(wape2(i).le.1.0))) THEN
    1534 ccc      IF (sigmaw(i).GT.0.9) THEN
    1535          wape(i) = 0.
    1536          cstar(i)= 0.  !!corr itlmd
    1537          hw(i) = hwmin
    1538          sigmaw(i) = sigmad
    1539          fip(i) = 0.
    1540         ELSE
    1541          wape(i) = wape2(i)
    1542          cstar(i)= cstar2(i) !!corr itlmd
    1543         ENDIF
    1544       ENDDO
    1545 c
    1546 c
    1547       RETURN
    1548       END
    1549 
    1550       SUBROUTINE wake_vec_modulation(nlon,nl,wk_adv,qe,d_qe,
    1551      $           deltaqw,d_deltaqw,sigmaw,d_sigmaw,alpha)
    1552 c------------------------------------------------------
    1553 cDtermination du coefficient alpha tel que les tendances
    1554 c corriges alpha*d_G, pour toutes les grandeurs G, correspondent
    1555 c a une humidite positive dans la zone (x) et dans la zone (w).
    1556 c------------------------------------------------------
    1557 c
    1558  
    1559 c  Input
    1560       REAL qe(nlon,nl),d_qe(nlon,nl)
    1561       REAL deltaqw(nlon,nl),d_deltaqw(nlon,nl)
    1562       REAL sigmaw(nlon),d_sigmaw(nlon)
    1563       LOGICAL wk_adv(nlon)
    1564       INTEGER nl,nlon
    1565 c  Output
    1566       REAL alpha(nlon)
    1567 c  Internal variables
    1568       REAL alpha1(nlon)
    1569       REAL x,a,b,c,discrim,zeta(nlon)
    1570       REAL epsilon
    1571       DATA epsilon/1.e-15/
    1572 c
    1573       DO k=1,nl
    1574       DO i = 1,nlon
    1575        IF (wk_adv(i)) THEN
    1576         IF ((deltaqw(i,k)+d_deltaqw(i,k)).ge.0.) then
    1577          zeta(i)=0.
    1578         ELSE
    1579          zeta(i)=1.
    1580         END IF
    1581        ENDIF
    1582       ENDDO
    1583       DO i = 1,nlon
    1584        IF (wk_adv(i)) THEN
    1585         x = qe(i,k)+(zeta(i)-sigmaw(i))*deltaqw(i,k)
    1586      $   +d_qe(i,k)+(zeta(i)-sigmaw(i))*d_deltaqw(i,k)
    1587      $   -d_sigmaw(i)*(deltaqw(i,k)+d_deltaqw(i,k))
    1588       a=-d_sigmaw(i)*d_deltaqw(i,k)
    1589       b=d_qe(i,k)+(zeta(i)-sigmaw(i))*d_deltaqw(i,k)
    1590      $           -deltaqw(i,k)*d_sigmaw(i)
    1591       c=qe(i,k)+(zeta(i)-sigmaw(i))*deltaqw(i,k)-epsilon
    1592 !       c=qe(i,k)+(zeta(i)-sigmaw(i))*deltaqw(i,k)
    1593 
    1594       discrim=b*b-4.*a*c
    1595 !       print*,'ZETA *********************' 
    1596 !       print*,'zeta sigmaw ',zeta(:)
    1597 !       print*,'SIGMA *********************'
    1598 !       print*,'sigmaw ',sigmaw(:)
    1599 
    1600 !       print*,' x ************************'
    1601 !       print*,'x ',x
    1602 !       print*,'  a+b ************************'
    1603 !       print*,'a+b ',a+b
    1604 
    1605 !       print*,'a b c delta zeta ',a,b,c,discrim
    1606         IF (a+b .GE. 0.) THEN
    1607          alpha1(i)=1.
    1608         ELSE
    1609          IF (x .GE. 0.) THEN
    1610             alpha1(i)=1.
    1611          ELSE
    1612 !              IF (a .GE. 0.) THEN
    1613               IF (a .GT. 0.) THEN
    1614 !       print*,'a b c delta zeta ',a,b,c,discrim,zeta(i)
    1615 !       print*,'-b+sqrt(discrim) ',-b+sqrt(discrim)
    1616                  alpha1(i)=0.9*min(   (2.*c)/(-b+sqrt(discrim)),
    1617      $                        (-b+sqrt(discrim))/(2.*a)   )
    1618               ELSE IF (a.eq.0.) THEN
    1619                  alpha1(i)=0.9*(-c/b)
    1620               ELSE
    1621 !       print*,'a b c delta zeta ',a,b,c,discrim,zeta(i)
    1622 !       print*,'-b+sqrt(discrim) ',-b+sqrt(discrim)
    1623                  alpha1(i)=0.9*max(   (2.*c)/(-b+sqrt(discrim)),
    1624      $                        (-b+sqrt(discrim))/(2.*a)   )
    1625               ENDIF
    1626          ENDIF
    1627         ENDIF
    1628        ENDIF
    1629       ENDDO
    1630       ENDDO
    1631 c
    1632       DO i = 1,nlon
    1633        IF (wk_adv(i)) THEN
    1634         alpha(i) = min(alpha(i),alpha1(i))
    1635        ENDIF
    1636       ENDDO
    1637 c
    1638       return
    1639       end
    1640 
    1641       Subroutine WAKE_scal (p,ph,ppi,dtime,sigd_con
    1642      :                ,te0,qe0,omgb
    1643      :                ,dtdwn,dqdwn,amdwn,amup,dta,dqa
    1644      :                ,wdtPBL,wdqPBL,udtPBL,udqPBL
    1645      o                ,deltatw,deltaqw,dth,hw,sigmaw,wape,fip,gfl
    1646      o                ,dtls,dqls
    1647      o                ,ktopw,omgbdth,dp_omgb,wdens
    1648      o                ,tu,qu
    1649      o                ,dtKE,dqKE
    1650      o                ,dtPBL,dqPBL
    1651      o                ,omg,dp_deltomg,spread
    1652      o                ,Cstar,d_deltat_gw
    1653      o                ,d_deltatw2,d_deltaqw2)
    1654 
    1655 ***************************************************************
    1656 *                                                             *
    1657 * WAKE                                                        *
    1658 *      retour a un Pupper fixe                                *
    1659 *                                                             *
    1660 * written by   :  GRANDPEIX Jean-Yves   09/03/2000            *
    1661 * modified by :   ROEHRIG Romain        01/29/2007            *
    1662 ***************************************************************
    1663 c
    1664       USE dimphy
    1665       IMPLICIT none
    1666 c============================================================================
    1667 C
    1668 C
    1669 C   But : Decrire le comportement des poches froides apparaissant dans les
    1670 C        grands systemes convectifs, et fournir l'energie disponible pour
    1671 C        le declenchement de nouvelles colonnes convectives.
    1672 C
    1673 C   Variables d'etat : deltatw    : ecart de temperature wake-undisturbed area
    1674 C                      deltaqw    : ecart d'humidite wake-undisturbed area
    1675 C                      sigmaw     : fraction d'aire occupee par la poche.
    1676 C
    1677 C   Variable de sortie :
    1678 c
    1679 c                        wape : WAke Potential Energy
    1680 c                        fip  : Front Incident Power (W/m2) - ALP
    1681 c                        gfl  : Gust Front Length per unit area (m-1)
    1682 C                        dtls : large scale temperature tendency due to wake
    1683 C                        dqls : large scale humidity tendency due to wake
    1684 C                        hw   : hauteur de la poche
    1685 C                     dp_omgb : vertical gradient of large scale omega
    1686 C                      omgbdth: flux of Delta_Theta transported by LS omega
    1687 C                      dtKE   : differential heating (wake - unpertubed)
    1688 C                      dqKE   : differential moistening (wake - unpertubed)
    1689 C                      omg    : Delta_omg =vertical velocity diff. wake-undist. (Pa/s)
    1690 C                 dp_deltomg  : vertical gradient of omg (s-1)
    1691 C                     spread  : spreading term in dt_wake and dq_wake
    1692 C                 deltatw     : updated temperature difference (T_w-T_u).
    1693 C                 deltaqw     : updated humidity difference (q_w-q_u).
    1694 C                 sigmaw      : updated wake fractional area.
    1695 C                 d_deltat_gw : delta T tendency due to GW
    1696 c
    1697 C   Variables d'entree :
    1698 c
    1699 c                        aire : aire de la maille
    1700 c                        te0  : temperature dans l'environnement  (K)
    1701 C                        qe0  : humidite dans l'environnement     (kg/kg)
    1702 C                        omgb : vitesse verticale moyenne sur la maille (Pa/s)
    1703 C                        dtdwn: source de chaleur due aux descentes (K/s)
    1704 C                        dqdwn: source d'humidite due aux descentes (kg/kg/s)
    1705 C                        dta  : source de chaleur due courants satures et detrain  (K/s)
    1706 C                        dqa  : source d'humidite due aux courants satures et detra (kg/kg/s)
    1707 C                        amdwn: flux de masse total des descentes, par unite de
    1708 C                                surface de la maille (kg/m2/s)
    1709 C                        amup : flux de masse total des ascendances, par unite de
    1710 C                                surface de la maille (kg/m2/s)
    1711 C                        p    : pressions aux milieux des couches (Pa)
    1712 C                        ph   : pressions aux interfaces (Pa)
    1713 C                        ppi  : (p/p_0)**kapa (adim)
    1714 C                        dtime: increment temporel (s)
    1715 c
    1716 C   Variables internes :
    1717 c
    1718 c                        rhow : masse volumique de la poche froide
    1719 C                        rho  : environment density at P levels
    1720 C                        rhoh : environment density at Ph levels
    1721 C                        te   : environment temperature | may change within
    1722 C                        qe   : environment humidity    | sub-time-stepping
    1723 C                        the  : environment potential temperature
    1724 C                        thu  : potential temperature in undisturbed area
    1725 C                        tu   :  temperature  in undisturbed area
    1726 C                        qu   : humidity in undisturbed area
    1727 C                      dp_omgb: vertical gradient og LS omega
    1728 C                      omgbw  : wake average vertical omega
    1729 C                     dp_omgbw: vertical gradient of omgbw
    1730 C                      omgbdq : flux of Delta_q transported by LS omega
    1731 C                        dth  : potential temperature diff. wake-undist.
    1732 C                        th1  : first pot. temp. for vertical advection (=thu)
    1733 C                        th2  : second pot. temp. for vertical advection (=thw)
    1734 C                        q1   : first humidity for vertical advection
    1735 C                        q2   : second humidity for vertical advection
    1736 C                     d_deltatw   : terme de redistribution pour deltatw
    1737 C                     d_deltaqw   : terme de redistribution pour deltaqw
    1738 C                      deltatw0   : deltatw initial
    1739 C                      deltaqw0   : deltaqw initial
    1740 C                      hw0    : hw initial
    1741 C                      sigmaw0: sigmaw initial
    1742 C                      amflux : horizontal mass flux through wake boundary
    1743 C                      wdens  : number of wakes per unit area (3D) or per
    1744 C                               unit length (2D)
    1745 C                      Tgw    : 1 sur la période de onde de gravité
    1746 c                      Cgw    : vitesse de propagation de onde de gravité
    1747 c                      LL     : distance entre 2 poches
    1748 
    1749 c-------------------------------------------------------------------------
    1750 c          Déclaration de variables
    1751 c-------------------------------------------------------------------------
    1752 
    1753 #include "dimensions.h"
    17541944cccc#include "dimphy.h"
    17551945#include "YOMCST.h"
Note: See TracChangeset for help on using the changeset viewer.