Changeset 974 for LMDZ4/trunk


Ignore:
Timestamp:
Jun 19, 2008, 12:26:15 PM (16 years ago)
Author:
lmdzadmin
Message:

Nouvelles versions vectorisees ; on garde versions scalaires; nom _scal
IM

Location:
LMDZ4/trunk/libf/phylmd
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • LMDZ4/trunk/libf/phylmd/calwake.F

    r953 r974  
    3030      IMPLICIT none
    3131c======================================================================
     32#include "dimensions.h"
     33!#include "dimphy.h"
     34#include "YOMCST.h"
     35
     36c Arguments
     37c----------
     38
     39      INTEGER  i,l,ktopw(klon)
     40      REAL   dtime
     41
     42      REAL paprs(klon,klev+1),pplay(klon,klev)
     43      REAL t(klon,klev), q(klon,klev), omgb(klon,klev)
     44      REAL dt_dwn(klon,klev), dq_dwn(klon,klev),M_dwn(klon,klev)
     45      REAL M_up(klon,klev)
     46      REAL dt_a(klon,klev), dq_a(klon,klev)
     47      REAL wdt_PBL(klon,klev), wdq_PBL(klon,klev)
     48      REAL udt_PBL(klon,klev), udq_PBL(klon,klev)
     49      REAL wake_deltat(klon,klev),wake_deltaq(klon,klev)
     50      REAL dt_wake(klon,klev),dq_wake(klon,klev)
     51      REAL wake_d_deltat_gw(klon,klev)
     52      REAL wake_h(klon),wake_s(klon)
     53      REAL wake_dth(klon,klev)
     54      REAL wake_pe(klon),wake_fip(klon),wake_gfl(klon)
     55      REAL undi_t(klon,klev),undi_q(klon,klev)
     56      REAL wake_omgbdth(klon,klev),wake_dp_omgb(klon,klev)
     57      REAL wake_dtKE(klon,klev),wake_dqKE(klon,klev)
     58      REAL wake_dtPBL(klon,klev),wake_dqPBL(klon,klev)
     59      REAL wake_omg(klon,klev+1),wake_dp_deltomg(klon,klev)
     60      REAL wake_spread(klon,klev),wake_Cstar(klon)
     61      REAL wake_ddeltat(klon,klev),wake_ddeltaq(klon,klev)
     62      REAL d_deltatw(klon,klev), d_deltaqw(klon,klev)
     63      INTEGER wake_k(klon)
     64      REAL sigd(klon)
     65      REAL wake_dens(klon)
     66
     67C  Variable internes
     68C  -----------------
     69
     70      REAL aire
     71      REAL p(klon,klev),ph(klon,klev+1),pi(klon,klev)
     72      REAL te(klon,klev),qe(klon,klev),omgbe(klon,klev)
     73      REAL dtdwn(klon,klev),dqdwn(klon,klev)
     74      REAL dta(klon,klev),dqa(klon,klev)
     75      REAL wdtPBL(klon,klev),wdqPBL(klon,klev)
     76      REAL udtPBL(klon,klev),udqPBL(klon,klev)
     77      REAL amdwn(klon,klev),amup(klon,klev)
     78      REAL dtw(klon,klev),dqw(klon,klev),dth(klon,klev)
     79      REAL d_deltat_gw(klon,klev)
     80      REAL dtls(klon,klev),dqls(klon,klev)
     81      REAL tu(klon,klev),qu(klon,klev)
     82      REAL hw(klon),sigmaw(klon),wape(klon),fip(klon),gfl(klon)
     83      REAL omgbdth(klon,klev),dp_omgb(klon,klev)
     84      REAL dtKE(klon,klev),dqKE(klon,klev)
     85      REAL dtPBL(klon,klev),dqPBL(klon,klev)
     86      REAL omg(klon,klev+1),dp_deltomg(klon,klev),spread(klon,klev)
     87      REAL Cstar(klon)
     88      REAL sigd0(klon),wdens(klon)
     89
     90      REAL RDCP
     91
     92c      print *, '-> calwake, wake_s ', wake_s(1)
     93
     94      RDCP=1./3.5
     95
     96c-----------------------------------------------------------
     97cIM 290108     DO 999 i=1,klon   ! a vectoriser
     98c----------------------------------------------------------
     99
     100
     101      DO l=1,klev
     102      DO i=1,klon
     103        p(i,l)= pplay(i,l)
     104        ph(i,l)= paprs(i,l)
     105        pi(i,l) = (pplay(i,l)/100000.)**RDCP
     106
     107        te(i,l) = t(i,l)
     108        qe(i,l) = q(i,l)
     109        omgbe(i,l) = omgb(i,l)
     110
     111        dtdwn(i,l)= dt_dwn(i,l)
     112        dqdwn(i,l)= dq_dwn(i,l)
     113        dta(i,l)= dt_a(i,l)
     114        dqa(i,l)= dq_a(i,l)
     115        wdtPBL(i,l)= wdt_PBL(i,l)
     116        wdqPBL(i,l)= wdq_PBL(i,l)
     117        udtPBL(i,l)= udt_PBL(i,l)
     118        udqPBL(i,l)= udq_PBL(i,l)
     119      ENDDO
     120      ENDDO
     121     
     122      DO i=1,klon
     123      sigd0(i)=sigd(i)
     124      ENDDO
     125c      print*, 'sigd0,sigd', sigd0, sigd(i)
     126      DO i=1,klon
     127      ph(i,klev+1)=0.
     128      ENDDO
     129
     130      DO i=1,klon
     131      ktopw(i) = wake_k(i)
     132      ENDDO
     133
     134      DO l=1,klev
     135      DO i=1,klon
     136        dtw(i,l) = wake_deltat(i,l)
     137        dqw(i,l) = wake_deltaq(i,l)
     138      ENDDO
     139      ENDDO
     140
     141      DO l=1,klev
     142      DO i=1,klon
     143        dtls(i,l)=dt_wake(i,l)
     144        dqls(i,l)=dq_wake(i,l)
     145      ENDDO
     146      ENDDO
     147
     148      DO i=1,klon
     149      hw(i) = wake_h(i)
     150      sigmaw(i)= wake_s(i)
     151      ENDDO
     152
     153cfkc les flux de masses sont evalues aux niveaux et valent 0 a la surface
     154cfkc  on veut le flux de masse au milieu des couches
     155
     156      DO l=1,klev-1
     157      DO i=1,klon
     158        amdwn(i,l)= 0.5*(M_dwn(i,l)+M_dwn(i,l+1))
     159        amdwn(i,l)= (M_dwn(i,l+1))
     160      ENDDO
     161      ENDDO
     162
     163c au sommet le flux de masse est nul
     164
     165      DO i=1,klon
     166      amdwn(i,klev)=0.5*M_dwn(i,klev)
     167      ENDDO
     168c
     169      DO l = 1,klev
     170      DO i=1,klon
     171        amup(i,l)=M_up(i,l)
     172      ENDDO
     173      ENDDO
     174
     175      call WAKE(p,ph,pi,dtime,sigd0
     176     $                ,te,qe,omgbe
     177     $                ,dtdwn,dqdwn,amdwn,amup,dta,dqa
     178     $                ,wdtPBL,wdqPBL,udtPBL,udqPBL
     179     $                ,dtw,dqw,dth,hw,sigmaw,wape,fip,gfl
     180     $                ,dtls,dqls,ktopw
     181     $                ,omgbdth,dp_omgb,wdens
     182     $                ,tu,qu
     183     $                ,dtKE,dqKE
     184     $                ,dtPBL,dqPBL
     185     $                ,omg,dp_deltomg,spread
     186     $                ,Cstar,d_deltat_gw
     187     $                ,d_deltatw,d_deltaqw)
     188
     189      DO i=1,klon
     190       IF (ktopw(i) .GT. 0) THEN
     191         DO l=1,klev
     192           wake_deltat(i,l)= dtw(i,l)
     193           wake_deltaq(i,l)= dqw(i,l)
     194           wake_d_deltat_gw(i,l)= d_deltat_gw(i,l)
     195           wake_omgbdth(i,l) = omgbdth(i,l)
     196           wake_dp_omgb(i,l) = dp_omgb(i,l)
     197           wake_dtKE(i,l) = dtKE(i,l)
     198           wake_dqKE(i,l) = dqKE(i,l)
     199           wake_dtPBL(i,l) = dtPBL(i,l)
     200           wake_dqPBL(i,l) = dqPBL(i,l)
     201           wake_omg(i,l) = omg(i,l)
     202           wake_dp_deltomg(i,l) = dp_deltomg(i,l)
     203           wake_spread(i,l) = spread(i,l)
     204           wake_dth(i,l) = dth(i,l)
     205           dt_wake(i,l) = dtls(i,l)
     206           dq_wake(i,l) = dqls(i,l)
     207           undi_t(i,l) = tu(i,l)
     208           undi_q(i,l) = qu(i,l)
     209           wake_ddeltat(i,l) = d_deltatw(i,l)
     210           wake_ddeltaq(i,l) = d_deltaqw(i,l)
     211         ENDDO
     212       ELSE
     213         DO l = 1,klev
     214           wake_deltat(i,l)= 0.
     215           wake_deltaq(i,l)= 0.
     216           wake_d_deltat_gw(i,l)= 0.
     217           wake_omgbdth(i,l) = 0.
     218           wake_dp_omgb(i,l) = 0.
     219           wake_dtKE(i,l) = 0.
     220           wake_dqKE(i,l) = 0.
     221           wake_omg(i,l) = 0.
     222           wake_dp_deltomg(i,l) = 0.
     223           wake_spread(i,l) = 0.
     224           wake_dth(i,l)=0.
     225           dt_wake(i,l)=0.
     226           dq_wake(i,l)=0.
     227           undi_t(i,l)=te(i,l)
     228           undi_q(i,l)=qe(i,l)
     229         ENDDO
     230       ENDIF
     231
     232       wake_h(i)= hw(i)
     233       wake_s(i)= sigmaw(i)
     234       wake_pe(i)= wape(i)
     235       wake_fip(i)= fip(i)
     236       wake_gfl(i) = gfl(i)
     237       wake_k(i) =ktopw(i)
     238       wake_Cstar(i) = Cstar(i)
     239       wake_dens(i) = wdens(i)
     240c
     241cIM 290108 999  CONTINUE
     242c
     243      ENDDO
     244      RETURN
     245      END
     246      SUBROUTINE CALWAKE_scal(paprs,pplay,dtime
     247     :             ,t,q,omgb
     248     :             ,dt_dwn,dq_dwn,M_dwn,M_up
     249     :             ,dt_a,dq_a,sigd
     250     :             ,wdt_PBL,wdq_PBL
     251     :             ,udt_PBL,udq_PBL
     252     o             ,wake_deltat,wake_deltaq,wake_dth
     253     o             ,wake_h,wake_s,wake_dens
     254     o             ,wake_pe,wake_fip,wake_gfl
     255     o             ,dt_wake,dq_wake
     256     o             ,wake_k
     257     o             ,undi_t,undi_q
     258     o             ,wake_omgbdth,wake_dp_omgb
     259     o             ,wake_dtKE,wake_dqKE
     260     o             ,wake_dtPBL,wake_dqPBL
     261     o             ,wake_omg,wake_dp_deltomg
     262     o             ,wake_spread,wake_Cstar,wake_d_deltat_gw
     263     o             ,wake_ddeltat,wake_ddeltaq)
     264***************************************************************
     265*                                                             *
     266* CALWAKE                                                     *
     267*           interface avec le schema de calcul de la poche    *
     268*           froide                                            *
     269*                                                             *
     270* written by   : CHERUY Frederique, 13/03/2000, 10.31.05      *
     271* modified by :  ROEHRIG Romain,    01/30/2007                *
     272***************************************************************
     273*
     274      USE dimphy
     275      IMPLICIT none
     276c======================================================================
    32277
    33278#include "dimensions.h"
    … …  
    151396      ENDDO
    152397
    153       call WAKE(p,ph,pi,dtime,sigd0
     398      call WAKE_scal(p,ph,pi,dtime,sigd0
    154399     $                ,te,qe,omgbe
    155400     $                ,dtdwn,dqdwn,amdwn,amup,dta,dqa

  • LMDZ4/trunk/libf/phylmd/wake.F

    r953 r974  
    2222***************************************************************
    2323c
    24       USE dimphy
     24      use dimphy
    2525      IMPLICIT none
    2626c============================================================================
    … …  
    112112
    113113#include "dimensions.h"
     114#include "YOMCST.h"
     115#include "cvthermo.h"
     116#include "iniprint.h"
     117
     118c Arguments en entree
     119c--------------------
     120
     121      REAL, dimension(klon,klev) :: p, ppi
     122      REAL, dimension(klon,klev+1) :: ph, omgb
     123      REAL dtime
     124      REAL, dimension(klon,klev) :: te0,qe0
     125      REAL, dimension(klon,klev) :: dtdwn, dqdwn
     126      REAL, dimension(klon,klev) :: wdtPBL,wdqPBL
     127      REAL, dimension(klon,klev) :: udtPBL,udqPBL
     128      REAL, dimension(klon,klev) :: amdwn, amup
     129      REAL, dimension(klon,klev) :: dta, dqa
     130      REAL, dimension(klon) :: sigd_con
     131
     132c Sorties
     133c--------
     134
     135      REAL, dimension(klon,klev) :: deltatw, deltaqw, dth
     136      REAL, dimension(klon,klev) :: tu, qu
     137      REAL, dimension(klon,klev) :: dtls, dqls
     138      REAL, dimension(klon,klev) :: dtKE, dqKE
     139      REAL, dimension(klon,klev) :: dtPBL, dqPBL
     140      REAL, dimension(klon,klev) :: spread
     141      REAL, dimension(klon,klev) :: d_deltatgw
     142      REAL, dimension(klon,klev) :: d_deltatw2, d_deltaqw2
     143      REAL, dimension(klon,klev+1) :: omgbdth, omg
     144      REAL, dimension(klon,klev) :: dp_omgb, dp_deltomg
     145      REAL, dimension(klon,klev) :: d_deltat_gw
     146      REAL, dimension(klon) :: hw, sigmaw, wape, fip, gfl, Cstar
     147      INTEGER, dimension(klon) :: ktopw
     148
     149c Variables internes
     150c-------------------
     151
     152c Variables à fixer
     153      REAL ALON
     154      REAL coefgw
     155      REAL :: wdens0, wdens
     156      REAL stark
     157      REAL alpk
     158      REAL delta_t_min
     159      REAL Pupper
     160      INTEGER nsub
     161      REAL dtimesub
     162      REAL sigmad, hwmin
     163cIM 080208
     164      LOGICAL, dimension(klon) :: gwake
     165
     166c Variables de sauvegarde
     167      REAL, dimension(klon,klev) :: deltatw0
     168      REAL, dimension(klon,klev) :: deltaqw0
     169      REAL, dimension(klon,klev) :: te, qe
     170      REAL, dimension(klon) :: sigmaw0, sigmaw1
     171
     172c Variables pour les GW
     173      REAL, DIMENSION(klon) :: LL
     174      REAL, dimension(klon,klev) :: N2
     175      REAL, dimension(klon,klev) :: Cgw
     176      REAL, dimension(klon,klev) :: Tgw
     177
     178c Variables liées au calcul de hw
     179      REAL, DIMENSION(klon) :: ptop_provis, ptop, ptop_new
     180      REAL, DIMENSION(klon) :: sum_dth
     181      REAL, DIMENSION(klon) :: dthmin
     182      REAL, DIMENSION(klon) :: z, dz, hw0
     183      INTEGER, DIMENSION(klon) :: ktop, kupper
     184
     185c Autres variables internes
     186      INTEGER isubstep, k, i
     187
     188      REAL, DIMENSION(klon) :: sum_thu, sum_tu, sum_qu,sum_thvu
     189      REAL, DIMENSION(klon) :: sum_dq, sum_rho
     190      REAL, DIMENSION(klon) :: sum_dtdwn, sum_dqdwn
     191      REAL, DIMENSION(klon) :: av_thu, av_tu, av_qu, av_thvu
     192      REAL, DIMENSION(klon) :: av_dth, av_dq, av_rho
     193      REAL, DIMENSION(klon) :: av_dtdwn, av_dqdwn
     194
     195      REAL, DIMENSION(klon,klev) :: rho, rhow
     196      REAL, DIMENSION(klon,klev+1) :: rhoh
     197      REAL, DIMENSION(klon,klev) :: rhow_moyen
     198      REAL, DIMENSION(klon,klev) :: zh
     199      REAL, DIMENSION(klon,klev+1) :: zhh
     200      REAL, DIMENSION(klon,klev) :: epaisseur1, epaisseur2
     201
     202      REAL, DIMENSION(klon,klev) :: the, thu
     203
     204      REAL, DIMENSION(klon,klev) :: d_deltatw, d_deltaqw
     205
     206      REAL, DIMENSION(klon,klev+1) :: omgbw
     207      REAL, DIMENSION(klon) :: omgtop
     208      REAL, DIMENSION(klon,klev) :: dp_omgbw
     209      REAL, DIMENSION(klon) :: ztop, dztop
     210      REAL, DIMENSION(klon,klev) :: alpha_up
     211     
     212      REAL, dimension(klon) :: RRe1, RRe2
     213      REAL :: RRd1, RRd2
     214      REAL, DIMENSION(klon,klev) :: Th1, Th2, q1, q2
     215      REAL, DIMENSION(klon,klev) :: D_Th1, D_Th2, D_dth
     216      REAL, DIMENSION(klon,klev) :: D_q1, D_q2, D_dq
     217      REAL, DIMENSION(klon,klev) :: omgbdq
     218
     219      REAL, dimension(klon) :: ff, gg
     220      REAL, dimension(klon) :: wape2, Cstar2, heff
     221
     222      REAL, DIMENSION(klon,klev) :: Crep
     223      REAL Crep_upper, Crep_sol
     224
     225C-------------------------------------------------------------------------
     226c         Initialisations
     227c-------------------------------------------------------------------------
     228
     229c      print*, 'wake initialisations'
     230
     231c   Essais d'initialisation avec sigmaw = 0.02 et hw = 10.
     232c-------------------------------------------------------------------------
     233
     234      DATA sigmad, hwmin /.02,10./
     235
     236C Longueur de maille (en m)
     237c-------------------------------------------------------------------------
     238
     239c      ALON = 3.e5
     240      ALON = 1.e6
     241
     242
     243C Configuration de coefgw,stark,wdens (22/02/06 by YU Jingmei)
     244c
     245c      coefgw : Coefficient pour les ondes de gravité
     246c       stark : Coefficient k dans Cstar=k*sqrt(2*WAPE)
     247c       wdens : Densité de poche froide par maille
     248c-------------------------------------------------------------------------
     249
     250      coefgw=10
     251c      coefgw=1
     252c      wdens0 = 1.0/(alon**2)   
     253      wdens = 1.0/(alon**2)       
     254      stark = 0.50
     255cCRtest
     256      alpk=0.1
     257c      alpk = 1.0
     258c      alpk = 0.5
     259c      alpk = 0.05
     260      Crep_upper=0.9
     261      Crep_sol=1.0
     262
     263
     264C Minimum value for |T_wake - T_undist|. Used for wake top definition
     265c-------------------------------------------------------------------------
     266
     267      delta_t_min = 0.2
     268
     269C 1. - Save initial values and initialize tendencies
     270C --------------------------------------------------
     271
     272      DO k=1,klev
     273      DO i=1, klon
     274        deltatw0(i,k) = deltatw(i,k)
     275        deltaqw0(i,k)= deltaqw(i,k)
     276        te(i,k) = te0(i,k)
     277        qe(i,k) = qe0(i,k)
     278        dtls(i,k) = 0.
     279        dqls(i,k) = 0.
     280        d_deltat_gw(i,k)=0.
     281!IM 060508 beg
     282        d_deltatw2(i,k)=0.
     283        d_deltaqw2(i,k)=0.
     284!IM 060508 end
     285      ENDDO
     286      ENDDO
     287c      sigmaw1=sigmaw
     288c      IF (sigd_con.GT.sigmaw1) THEN
     289c      print*, 'sigmaw,sigd_con', sigmaw, sigd_con
     290c      ENDIF
     291      DO i=1, klon
     292cc      sigmaw(i) = amax1(sigmaw(i),sigd_con(i))
     293      sigmaw(i) = amax1(sigmaw(i),sigmad)
     294      sigmaw(i) = amin1(sigmaw(i),0.99)
     295      sigmaw0(i) = sigmaw(i)
     296      ENDDO
     297C
     298C
     299C 2. - Prognostic part
     300C --------------------
     301C
     302C
     303C 2.1 - Undisturbed area and Wake integrals
     304C ---------------------------------------------------------
     305
     306      DO i=1, klon
     307      z(i) = 0.
     308      ktop(i)=0
     309      kupper(i) = 0
     310      sum_thu(i) = 0.
     311      sum_tu(i) = 0.
     312      sum_qu(i) = 0.
     313      sum_thvu(i) = 0.
     314      sum_dth(i) = 0.
     315      sum_dq(i) = 0.
     316      sum_rho(i) = 0.
     317      sum_dtdwn(i) = 0.
     318      sum_dqdwn(i) = 0.
     319
     320      av_thu(i) = 0.
     321      av_tu(i) =0.
     322      av_qu(i) =0.
     323      av_thvu(i) = 0.
     324      av_dth(i) = 0.
     325      av_dq(i) = 0.
     326      av_rho(i) =0.
     327      av_dtdwn(i) =0.
     328      av_dqdwn(i) = 0.
     329      ENDDO
     330c
     331c Distance between wakes
     332       DO i = 1,klon
     333        LL(i) = (1-sqrt(sigmaw(i)))/sqrt(wdens)
     334       ENDDO
     335C Potential temperatures and humidity
     336c----------------------------------------------------------
     337      DO k =1,klev
     338       DO i=1, klon
     339        rho(i,k) = p(i,k)/(rd*te(i,k))
     340        IF(k .eq. 1) THEN
     341          rhoh(i,k) = ph(i,k)/(rd*te(i,k))
     342          zhh(i,k)=0
     343        ELSE
     344          rhoh(i,k) = ph(i,k)*2./(rd*(te(i,k)+te(i,k-1)))
     345          zhh(i,k)=(ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*RG)+zhh(i,k-1)
     346        ENDIF
     347        the(i,k) = te(i,k)/ppi(i,k)
     348        thu(i,k) = (te(i,k) - deltatw(i,k)*sigmaw(i))/ppi(i,k)
     349        tu(i,k) = te(i,k) - deltatw(i,k)*sigmaw(i)
     350        qu(i,k)  =  qe(i,k) - deltaqw(i,k)*sigmaw(i)
     351        rhow(i,k) = p(i,k)/(rd*(te(i,k)+deltatw(i,k)))
     352        dth(i,k) = deltatw(i,k)/ppi(i,k)
     353       ENDDO
     354      ENDDO
     355       
     356      DO k = 1, klev-1
     357      DO i=1, klon
     358        IF(k.eq.1) THEN
     359          N2(i,k)=0
     360        ELSE
     361          N2(i,k)=amax1(0.,-RG**2/the(i,k)*rho(i,k)*(the(i,k+1)-
     362     $            the(i,k-1))/(p(i,k+1)-p(i,k-1)))
     363        ENDIF
     364        ZH(i,k)=(zhh(i,k)+zhh(i,k+1))/2
     365
     366        Cgw(i,k)=sqrt(N2(i,k))*ZH(i,k)
     367        Tgw(i,k)=coefgw*Cgw(i,k)/LL(i)
     368      ENDDO
     369      ENDDO
     370
     371      DO i=1, klon
     372      N2(i,klev)=0
     373      ZH(i,klev)=0
     374      Cgw(i,klev)=0
     375      Tgw(i,klev)=0
     376      ENDDO
     377
     378c  Calcul de la masse volumique moyenne de la colonne   (bdlmd)
     379c-----------------------------------------------------------------
     380
     381      DO k=1,klev
     382       DO i=1, klon
     383        epaisseur1(i,k)=0.
     384        epaisseur2(i,k)=0.
     385       ENDDO
     386      ENDDO
     387
     388      DO i=1, klon
     389      epaisseur1(i,1)= -(ph(i,2)-ph(i,1))/(rho(i,1)*rg)+1.
     390      epaisseur2(i,1)= -(ph(i,2)-ph(i,1))/(rho(i,1)*rg)+1.
     391      rhow_moyen(i,1) = rhow(i,1)
     392      ENDDO
     393
     394      DO k = 2, klev
     395      DO i=1, klon
     396        epaisseur1(i,k)= -(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg) +1.
     397        epaisseur2(i,k)=epaisseur2(i,k-1)+epaisseur1(i,k)
     398        rhow_moyen(i,k) = (rhow_moyen(i,k-1)*epaisseur2(i,k-1)+
     399     $                 rhow(i,k)*epaisseur1(i,k))/epaisseur2(i,k)
     400      ENDDO
     401      ENDDO
     402
     403C
     404C Choose an integration bound well above wake top
     405c-----------------------------------------------------------------
     406c
     407C       Pupper = 50000.  ! melting level
     408       Pupper = 60000.
     409c       Pupper = 80000.  ! essais pour case_e
     410C
     411C    Determine Wake top pressure (Ptop) from buoyancy integral
     412C    --------------------------------------------------------
     413c
     414c-1/ Pressure of the level where dth becomes less than delta_t_min.
     415
     416      DO i=1,klon
     417      ptop_provis(i)=ph(i,1)
     418      ENDDO
     419      DO k= 2,klev
     420      DO i=1,klon
     421c
     422cIM v3JYG; ptop_provis(i).LT. ph(i,1)
     423c
     424        IF (dth(i,k) .GT. -delta_t_min .and.
     425     $      dth(i,k-1).LT. -delta_t_min .and.
     426     $      ptop_provis(i).EQ. ph(i,1)) THEN
     427          ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)
     428     $          - (dth(i,k-1)+delta_t_min)*p(i,k)) /
     429     $          (dth(i,k) - dth(i,k-1))
     430        ENDIF
     431      ENDDO
     432      ENDDO
     433
     434c-2/ dth integral
     435
     436      DO i=1,klon
     437      sum_dth(i) = 0.
     438      dthmin(i) = -delta_t_min
     439      z(i) = 0.
     440      ENDDO
     441
     442      DO k = 1,klev
     443      DO i=1,klon
     444        dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-Ph(i,k))/(rho(i,k)*rg)
     445        IF (dz(i) .gt. 0) THEN
     446          z(i) = z(i)+dz(i)
     447          sum_dth(i) = sum_dth(i) + dth(i,k)*dz(i)
     448          dthmin(i) = amin1(dthmin(i),dth(i,k))
     449        ENDIF
     450      ENDDO
     451      ENDDO
     452
     453c-3/ height of triangle with area= sum_dth and base = dthmin
     454
     455      DO i=1,klon
     456      hw0(i) = 2.*sum_dth(i)/amin1(dthmin(i),-0.5)
     457      hw0(i) = amax1(hwmin,hw0(i))
     458      ENDDO
     459
     460c-4/ now, get Ptop
     461
     462      DO i=1,klon
     463      z(i) = 0.
     464      ptop(i) = ph(i,1)
     465      ENDDO
     466
     467      DO k = 1,klev
     468      DO i=1,klon
     469        dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg),hw0(i)-z(i))
     470        IF (dz(i) .gt. 0) THEN
     471         z(i) = z(i)+dz(i)
     472         ptop(i) = ph(i,k)-rho(i,k)*rg*dz(i)
     473        ENDIF
     474      ENDDO
     475      ENDDO
     476
     477
     478C-5/ Determination de ktop et kupper
     479
     480      DO k=klev,1,-1
     481      DO i=1,klon
     482        IF (ph(i,k+1) .lt. ptop(i)) ktop(i)=k
     483        IF (ph(i,k+1) .lt. pupper) kupper(i)=k
     484      ENDDO
     485      ENDDO
     486
     487c-6/ Correct ktop and ptop
     488
     489      DO i = 1,klon
     490        ptop_new(i)=ptop(i)
     491      ENDDO
     492      DO k= klev,2,-1
     493      DO i=1,klon
     494        IF (k .LE. ktop(i) .and.
     495     $      ptop_new(i) .EQ. ptop(i) .and.
     496     $      dth(i,k) .GT. -delta_t_min .and.
     497     $      dth(i,k-1).LT. -delta_t_min) THEN
     498          ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)
     499     $          - (dth(i,k-1)+delta_t_min)*p(i,k)) /
     500     $          (dth(i,k) - dth(i,k-1))
     501        ENDIF
     502      ENDDO
     503      ENDDO
     504
     505      DO i=1,klon
     506        ptop(i) = ptop_new(i)
     507      ENDDO
     508
     509      DO k=klev,1,-1
     510      DO i=1,klon
     511        IF (ph(i,k+1) .lt. ptop(i)) ktop(i)=k
     512      ENDDO
     513      ENDDO
     514c
     515c-5/ Set deltatw & deltaqw to 0 above kupper
     516c
     517      DO k = 1,klev
     518      DO i=1,klon
     519       IF (k.GE. kupper(i)) THEN
     520        deltatw(i,k) = 0.
     521        deltaqw(i,k) = 0.
     522       ENDIF
     523      ENDDO
     524      ENDDO
     525c
     526C
     527C Vertical gradient of LS omega
     528C
     529      DO k = 1,klev
     530      DO i=1,klon
     531       IF (k.LE. kupper(i)) THEN
     532        dp_omgb(i,k) = (omgb(i,k+1) - omgb(i,k))/(ph(i,k+1)-ph(i,k))
     533       ENDIF
     534      ENDDO
     535      ENDDO
     536C
     537C Integrals (and wake top level number)
     538C --------------------------------------
     539C
     540C Initialize sum_thvu to 1st level virt. pot. temp.
     541
     542      DO i=1,klon
     543      z(i) = 1.
     544      dz(i) = 1.
     545      sum_thvu(i) =  thu(i,1)*(1.+eps*qu(i,1))*dz(i)
     546      sum_dth(i) = 0.
     547      ENDDO
     548
     549      DO k = 1,klev
     550      DO i=1,klon
     551        dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
     552        IF (dz(i) .GT. 0) THEN
     553         z(i) = z(i)+dz(i)
     554         sum_thu(i) = sum_thu(i) + thu(i,k)*dz(i)
     555         sum_tu(i) = sum_tu(i) + tu(i,k)*dz(i)
     556         sum_qu(i) = sum_qu(i) + qu(i,k)*dz(i)
     557         sum_thvu(i) = sum_thvu(i) + thu(i,k)*(1.+eps*qu(i,k))*dz(i)
     558         sum_dth(i) = sum_dth(i) + dth(i,k)*dz(i)
     559         sum_dq(i) = sum_dq(i) + deltaqw(i,k)*dz(i)
     560         sum_rho(i) = sum_rho(i) + rhow(i,k)*dz(i)
     561         sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i,k)*dz(i)
     562         sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i,k)*dz(i)
     563        ENDIF
     564      ENDDO
     565      ENDDO
     566c
     567      DO i=1,klon
     568        hw0(i) = z(i)
     569      ENDDO
     570c
     571C
     572C 2.1 - WAPE and mean forcing computation
     573C ---------------------------------------
     574C
     575C ---------------------------------------
     576C
     577C Means
     578
     579      DO i=1,klon
     580      av_thu(i) = sum_thu(i)/hw0(i)
     581      av_tu(i) = sum_tu(i)/hw0(i)
     582      av_qu(i) = sum_qu(i)/hw0(i)
     583      av_thvu(i) = sum_thvu(i)/hw0(i)
     584c      av_thve = sum_thve/hw0
     585      av_dth(i) = sum_dth(i)/hw0(i)
     586      av_dq(i) = sum_dq(i)/hw0(i)
     587      av_rho(i) = sum_rho(i)/hw0(i)
     588      av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
     589      av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
     590
     591      wape(i) = - rg*hw0(i)*(av_dth(i)
     592     $     + eps*(av_thu(i)*av_dq(i)+av_dth(i)*av_qu(i)+av_dth(i)*
     593     $     av_dq(i) ))/av_thvu(i)
     594      ENDDO
     595C
     596C 2.2 Prognostic variable update
     597C ------------------------------
     598C
     599C Filter out bad wakes
     600
     601      DO k = 1,klev
     602       DO i=1,klon
     603        IF ( wape(i) .LT. 0.) THEN
     604          deltatw(i,k) = 0.
     605          deltaqw(i,k) = 0.
     606          dth(i,k) = 0.
     607        ENDIF
     608       ENDDO
     609      ENDDO
     610c
     611      DO i=1,klon
     612      IF ( wape(i) .LT. 0.) THEN
     613        wape(i) = 0.
     614        Cstar(i) = 0.
     615        hw(i) = hwmin
     616        sigmaw(i) = amax1(sigmad,sigd_con(i))
     617        fip(i) = 0.
     618        gwake(i) = .FALSE.
     619      ELSE
     620        Cstar(i) = stark*sqrt(2.*wape(i))
     621        gwake(i) = .TRUE.
     622      ENDIF
     623      ENDDO
     624c
     625C
     626CC -----------------------------------------------------------------
     627C    Sub-time-stepping
     628C    -----------------
     629C
     630      nsub=10
     631      dtimesub=dtime/nsub
     632c
     633c------------------------------------------------------------
     634      DO isubstep = 1,nsub
     635c------------------------------------------------------------
     636      DO i=1,klon
     637        gfl(i) = 2.*sqrt(3.14*wdens*sigmaw(i))
     638      ENDDO
     639      DO i=1,klon
     640        sigmaw(i) =sigmaw(i) + gfl(i)*Cstar(i)*dtimesub
     641        sigmaw(i) =amin1(sigmaw(i),0.99)     !!!!!!!!
     642c        wdens = wdens0/(10.*sigmaw)
     643c        sigmaw =max(sigmaw,sigd_con)
     644c        sigmaw =max(sigmaw,sigmad)
     645      ENDDO
     646C
     647C
     648c calcul de la difference de vitesse verticale poche - zone non perturbee
     649cIM 060208 differences par rapport au code initial; init. a 0 dp_deltomg
     650cIM 060208 et omg sur les niveaux de 1 a klev+1, alors que avant l'on definit
     651cIM 060208 au niveau k=1..?
     652      dp_deltomg(1:klon,1:klev)=0.
     653      DO k= 1,klev+1
     654      DO i = 1,klon
     655        omg(i,k)=0.
     656      ENDDO
     657      ENDDO
     658c
     659      DO i=1,klon
     660        z(i)= 0.
     661        omg(i,1) = 0.
     662        dp_deltomg(i,1) = -(gfl(i)*Cstar(i))/(sigmaw(i) * (1-sigmaw(i)))
     663      ENDDO
     664c
     665      DO k= 2,klev
     666      DO i = 1,klon
     667       IF( k .LE. ktop(i)) THEN
     668          dz(i) = -(ph(i,k)-ph(i,k-1))/(rho(i,k-1)*rg)
     669          z(i) = z(i)+dz(i)
     670          dp_deltomg(i,k)= dp_deltomg(i,1)
     671          omg(i,k)= dp_deltomg(i,1)*z(i)
     672       ENDIF
     673      ENDDO
     674      ENDDO
     675c
     676      DO i = 1,klon
     677        dztop(i)=-(ptop(i)-ph(i,ktop(i)))/(rho(i,ktop(i))*rg)
     678        ztop(i) = z(i)+dztop(i)
     679        omgtop(i)=dp_deltomg(i,1)*ztop(i)
     680      ENDDO
     681c
     682c        -----------------
     683c        From m/s to Pa/s
     684c        -----------------
     685c
     686       DO i=1,klon
     687        omgtop(i) = -rho(i,ktop(i))*rg*omgtop(i)
     688        dp_deltomg(i,1) = omgtop(i)/(ptop(i)-ph(i,1))
     689       ENDDO
     690c
     691      DO k= 1,klev
     692      DO i = 1,klon
     693       IF( k .LE. ktop(i)) THEN
     694          omg(i,k) = - rho(i,k)*rg*omg(i,k)
     695          dp_deltomg(i,k) = dp_deltomg(i,1)
     696       ENDIF
     697      ENDDO
     698      ENDDO
     699c
     700c   raccordement lineaire de omg de ptop a pupper
     701
     702      DO i=1,klon
     703      IF (kupper(i) .GT. ktop(i)) THEN
     704        omg(i,kupper(i)+1) = - Rg*amdwn(i,kupper(i)+1)/sigmaw(i)
     705     $                + Rg*amup(i,kupper(i)+1)/(1.-sigmaw(i))
     706        dp_deltomg(i,kupper(i)) = (omgtop(i)-omg(i,kupper(i)+1))/
     707     $                     (ptop(i)-pupper)
     708      ENDIF
     709      ENDDO
     710c
     711      DO k= 1,klev
     712      DO i = 1,klon
     713       IF( k .GT. ktop(i) .AND. k .LE. kupper(i)) THEN
     714          dp_deltomg(i,k) = dp_deltomg(i,kupper(i))
     715          omg(i,k) = omgtop(i)+(ph(i,k)-ptop(i))*dp_deltomg(i,kupper(i))
     716       ENDIF
     717      ENDDO
     718      ENDDO
     719c
     720c--    Compute wake average vertical velocity omgbw
     721c
     722c
     723      DO k = 1,klev+1
     724      DO i=1,klon
     725        omgbw(i,k) = omgb(i,k)+(1.-sigmaw(i))*omg(i,k)
     726      ENDDO
     727      ENDDO
     728c--    and its vertical gradient dp_omgbw
     729c
     730      DO k = 1,klev
     731      DO i=1,klon
     732        dp_omgbw(i,k) = (omgbw(i,k+1)-omgbw(i,k))/(ph(i,k+1)-ph(i,k))
     733      ENDDO
     734      ENDDO
     735C
     736c--    Upstream coefficients for omgb velocity
     737c--    (alpha_up(k) is the coefficient of the value at level k)
     738c--    (1-alpha_up(k) is the coefficient of the value at level k-1)
     739      DO k = 1,klev
     740      DO i=1,klon
     741       alpha_up(i,k) = 0.
     742       IF (omgb(i,k) .GT. 0.) alpha_up(i,k) = 1.
     743      ENDDO
     744      ENDDO
     745
     746c  Matrix expressing [The,deltatw] from  [Th1,Th2]
     747
     748      DO i=1,klon
     749        RRe1(i) = 1.-sigmaw(i)
     750        RRe2(i) = sigmaw(i)
     751      ENDDO
     752      RRd1 = -1.
     753      RRd2 = 1.
     754c
     755c--    Get [Th1,Th2], dth and [q1,q2]
     756c
     757      DO k= 1,klev
     758      DO i = 1,klon
     759       IF(k .LE. kupper(i)+1) THEN
     760        dth(i,k) = deltatw(i,k)/ppi(i,k)
     761        Th1(i,k) = the(i,k) - sigmaw(i)     *dth(i,k)   ! undisturbed area
     762        Th2(i,k) = the(i,k) + (1.-sigmaw(i))*dth(i,k)   ! wake
     763        q1(i,k) = qe(i,k) - sigmaw(i)     *deltaqw(i,k) ! undisturbed area
     764        q2(i,k) = qe(i,k) + (1.-sigmaw(i))*deltaqw(i,k) ! wake
     765       ENDIF
     766      ENDDO
     767      ENDDO
     768
     769      DO i=1,klon
     770       D_Th1(i,1) = 0.
     771       D_Th2(i,1) = 0.
     772       D_dth(i,1) = 0.
     773       D_q1(i,1) = 0.
     774       D_q2(i,1) = 0.
     775       D_dq(i,1) = 0.
     776      ENDDO
     777
     778      DO k= 2,klev
     779      DO i = 1,klon
     780       IF(k .LE. kupper(i)+1) THEN
     781        D_Th1(i,k) = Th1(i,k-1)-Th1(i,k)
     782        D_Th2(i,k) = Th2(i,k-1)-Th2(i,k)
     783        D_dth(i,k) = dth(i,k-1)-dth(i,k)
     784        D_q1(i,k) = q1(i,k-1)-q1(i,k)
     785        D_q2(i,k) = q2(i,k-1)-q2(i,k)
     786        D_dq(i,k) = deltaqw(i,k-1)-deltaqw(i,k)
     787       ENDIF
     788      ENDDO
     789      ENDDO
     790
     791      DO i=1,klon
     792        omgbdth(i,1) = 0.
     793        omgbdq(i,1) = 0.
     794      ENDDO
     795
     796      DO k= 2,klev
     797      DO i = 1,klon
     798       IF(k .LE. kupper(i)+1) THEN  !   loop on interfaces
     799        omgbdth(i,k) = omgb(i,k)*(    dth(i,k-1) -     dth(i,k))
     800        omgbdq(i,k)  = omgb(i,k)*(deltaqw(i,k-1) - deltaqw(i,k))
     801       ENDIF
     802      ENDDO
     803      ENDDO
     804c
     805c-----------------------------------------------------------------
     806      DO k= 1,klev
     807      DO i = 1,klon
     808       IF(k .LE. kupper(i)-1) THEN
     809c-----------------------------------------------------------------
     810c
     811c   Compute redistribution (advective) term
     812c
     813         d_deltatw(i,k) =
     814     $             dtimesub/(Ph(i,k)-Ph(i,k+1))*(
     815     $       RRd1*omg(i,k  )*sigmaw(i)     *D_Th1(i,k)
     816     $      -RRd2*omg(i,k+1)*(1.-sigmaw(i))*D_Th2(i,k+1)
     817     $      -(1.-alpha_up(i,k))*omgbdth(i,k) - alpha_up(i,k+1)*
     818     $      omgbdth(i,k+1))*ppi(i,k)
     819c         print*,'d_deltatw=',d_deltatw(i,k)
     820c
     821         d_deltaqw(i,k) =
     822     $             dtimesub/(Ph(i,k)-Ph(i,k+1))*(
     823     $       RRd1*omg(i,k  )*sigmaw(i)     *D_q1(i,k)
     824     $      -RRd2*omg(i,k+1)*(1.-sigmaw(i))*D_q2(i,k+1)
     825     $      -(1.-alpha_up(i,k))*omgbdq(i,k) - alpha_up(i,k+1)*
     826     $      omgbdq(i,k+1))
     827c         print*,'d_deltaqw=',d_deltaqw(i,k)
     828c
     829c   and increment large scale tendencies
     830c
     831         dtls(i,k) = dtls(i,k) +
     832     $               dtimesub*(
     833     $        ( RRe1(i)*omg(i,k  )*sigmaw(i)     *D_Th1(i,k)
     834     $         -RRe2(i)*omg(i,k+1)*(1.-sigmaw(i))*D_Th2(i,k+1) )
     835     $               /(Ph(i,k)-Ph(i,k+1))
     836     $         -sigmaw(i)*(1.-sigmaw(i))*dth(i,k)*dp_deltomg(i,k)
     837     $                      )*ppi(i,k)
     838c         print*,'dtls=',dtls(i,k)
     839c
     840         dqls(i,k) = dqls(i,k) +
     841     $               dtimesub*(
     842     $        ( RRe1(i)*omg(i,k  )*sigmaw(i)     *D_q1(i,k)
     843     $         -RRe2(i)*omg(i,k+1)*(1.-sigmaw(i))*D_q2(i,k+1) )
     844     $               /(Ph(i,k)-Ph(i,k+1))
     845     $         -sigmaw(i)*(1.-sigmaw(i))*deltaqw(i,k)*dp_deltomg(i,k)
     846     $                      )
     847c         print*,'dqls=',dqls(k)
     848       ENDIF
     849c-------------------------------------------------------------------
     850      ENDDO
     851      ENDDO
     852c------------------------------------------------------------------
     853C
     854C   Increment state variables
     855
     856      DO k= 1,klev
     857      DO i = 1,klon
     858       IF(k .LE. kupper(i)-1) THEN
     859c
     860c Coefficient de répartition
     861
     862        Crep(i,k)=Crep_sol*(ph(i,kupper(i))-ph(i,k))/(ph(i,kupper(i))
     863     $          -ph(i,1))
     864        Crep(i,k)=Crep(i,k)+Crep_upper*(ph(i,1)-ph(i,k))/(p(i,1)-
     865     $          ph(i,kupper(i)))
     866       
     867
     868c Reintroduce compensating subsidence term.
     869
     870c        dtKE(k)=(dtdwn(k)*Crep(k))/sigmaw
     871c        dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k))
     872c     .                   /(1-sigmaw)
     873c        dqKE(k)=(dqdwn(k)*Crep(k))/sigmaw
     874c        dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k))
     875c     .                   /(1-sigmaw)
     876c
     877c        dtKE(k)=(dtdwn(k)*Crep(k)+(1-Crep(k))*dta(k))/sigmaw
     878c        dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k)*Crep(k))
     879c     .                   /(1-sigmaw)
     880c        dqKE(k)=(dqdwn(k)*Crep(k)+(1-Crep(k))*dqa(k))/sigmaw
     881c        dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k)*Crep(k))
     882c     .                   /(1-sigmaw)
     883
     884        dtKE(i,k)=(dtdwn(i,k)/sigmaw(i) - dta(i,k)/(1.-sigmaw(i)))
     885        dqKE(i,k)=(dqdwn(i,k)/sigmaw(i) - dqa(i,k)/(1.-sigmaw(i)))
     886c        print*,'dtKE=',dtKE(k)
     887c        print*,'dqKE=',dqKE(k)
     888c
     889        dtPBL(i,k)=(wdtPBL(i,k)/sigmaw(i) - udtPBL(i,k)/(1.-sigmaw(i)))
     890        dqPBL(i,k)=(wdqPBL(i,k)/sigmaw(i) - udqPBL(i,k)/(1.-sigmaw(i)))
     891c
     892        spread(i,k) = (1.-sigmaw(i))*dp_deltomg(i,k)+gfl(i)*Cstar(i)/
     893     $  sigmaw(i)
     894
     895
     896c ajout d'un effet onde de gravité -Tgw(k)*deltatw(k) 03/02/06 YU Jingmei
     897
     898        d_deltat_gw(i,k)=d_deltat_gw(i,k)-Tgw(i,k)*deltatw(i,k)*
     899     $  dtimesub
     900        ff(i)=d_deltatw(i,k)/dtimesub
     901
     902c Sans GW
     903c
     904c        deltatw(k)=deltatw(k)+dtimesub*(ff+dtKE(k)-spread(k)*deltatw(k))
     905c
     906c GW formule 1
     907c
     908c        deltatw(k) = deltatw(k)+dtimesub*
     909c     $         (ff+dtKE(k) - spread(k)*deltatw(k)-Tgw(k)*deltatw(k))
     910c
     911c GW formule 2
     912
     913        IF (dtimesub*Tgw(i,k).lt.1.e-10) THEN
     914          deltatw(i,k) = deltatw(i,k)+dtimesub*
     915     $          (ff(i)+dtKE(i,k)+dtPBL(i,k)
     916     $          - spread(i,k)*deltatw(i,k)-Tgw(i,k)*deltatw(i,k))
     917        ELSE
     918           deltatw(i,k) = deltatw(i,k)+1/Tgw(i,k)*(1-exp(-dtimesub*
     919     $          Tgw(i,k)))*
     920     $          (ff(i)+dtKE(i,k)+dtPBL(i,k)
     921     $          - spread(i,k)*deltatw(i,k)-Tgw(i,k)*deltatw(i,k))
     922        ENDIF
     923   
     924        dth(i,k) = deltatw(i,k)/ppi(i,k)
     925
     926        gg(i)=d_deltaqw(i,k)/dtimesub
     927
     928       deltaqw(i,k) = deltaqw(i,k) +
     929     $         dtimesub*(gg(i)+ dqKE(i,k)+dqPBL(i,k) - spread(i,k)*
     930     $         deltaqw(i,k))
     931
     932       d_deltatw2(i,k)=d_deltatw2(i,k)+d_deltatw(i,k)
     933       d_deltaqw2(i,k)=d_deltaqw2(i,k)+d_deltaqw(i,k)
     934       ENDIF
     935      ENDDO
     936      ENDDO
     937
     938C   And update large scale variables
     939cIM 060208 manque DO i + remplace DO k=1,kupper(i)
     940cIM 060208     DO k = 1,kupper(i)
     941      DO k= 1,klev
     942      DO i = 1,klon
     943       IF(k .LE. kupper(i)) THEN
     944        te(i,k) = te0(i,k) + dtls(i,k)
     945        qe(i,k) = qe0(i,k) + dqls(i,k)
     946        the(i,k) = te(i,k)/ppi(i,k)
     947       ENDIF
     948      ENDDO
     949      ENDDO
     950c
     951C
     952c     Determine Ptop from buoyancy integral
     953c     ---------------------------------------
     954c
     955c-     1/ Pressure of the level where dth changes sign.
     956c
     957      DO i=1,klon
     958      Ptop_provis(i)=ph(i,1)
     959      ENDDO
     960c
     961      DO k= 2,klev
     962      DO i=1,klon
     963        IF (Ptop_provis(i) .EQ. ph(i,1) .AND.
     964     $      dth(i,k) .GT. -delta_t_min .and.
     965     $      dth(i,k-1).LT. -delta_t_min) THEN
     966          Ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)
     967     $          - (dth(i,k-1)+delta_t_min)*p(i,k)) /(dth(i,k)
     968     $          - dth(i,k-1))
     969        ENDIF
     970      ENDDO
     971      ENDDO
     972c
     973c-     2/ dth integral
     974c
     975      DO i=1,klon
     976       sum_dth(i) = 0.
     977       dthmin(i) = -delta_t_min
     978       z(i) = 0.
     979      ENDDO
     980
     981      DO k = 1,klev
     982      DO i=1,klon
     983        dz(i) = -(amax1(ph(i,k+1),Ptop_provis(i))-Ph(i,k))/(rho(i,k)*rg)
     984        IF (dz(i) .gt. 0) THEN
     985         z(i) = z(i)+dz(i)
     986         sum_dth(i) = sum_dth(i) + dth(i,k)*dz(i)
     987         dthmin(i) = amin1(dthmin(i),dth(i,k))
     988        ENDIF
     989      ENDDO
     990      ENDDO
     991c
     992c-     3/ height of triangle with area= sum_dth and base = dthmin
     993
     994      DO i=1,klon
     995       hw(i) = 2.*sum_dth(i)/amin1(dthmin(i),-0.5)
     996       hw(i) = amax1(hwmin,hw(i))
     997      ENDDO
     998c
     999c-     4/ now, get Ptop
     1000c
     1001      DO i=1,klon
     1002       ktop(i) = 0
     1003       z(i)=0.
     1004      ENDDO
     1005c
     1006      DO k = 1,klev
     1007      DO i=1,klon
     1008        dz(i) = amin1(-(ph(i,k+1)-Ph(i,k))/(rho(i,k)*rg),hw(i)-z(i))
     1009        IF (dz(i) .gt. 0) THEN
     1010         z(i) = z(i)+dz(i)
     1011         Ptop(i) = Ph(i,k)-rho(i,k)*rg*dz(i)
     1012         ktop(i) = k
     1013        ENDIF
     1014      ENDDO
     1015      ENDDO
     1016c
     1017c      4.5/Correct ktop and ptop
     1018c
     1019      DO i=1,klon
     1020        Ptop_new(i)=ptop(i)
     1021      ENDDO
     1022c
     1023      DO k= klev,2,-1
     1024      DO i=1,klon
     1025cIM v3JYG; IF (k .GE. ktop(i)
     1026       IF (k .LE. ktop(i) .AND.
     1027     $      ptop_new(i) .EQ. ptop(i) .AND.
     1028     $      dth(i,k) .GT. -delta_t_min .and.
     1029     $      dth(i,k-1).LT. -delta_t_min) THEN
     1030          Ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)
     1031     $          - (dth(i,k-1)+delta_t_min)*p(i,k)) /(dth(i,k)
     1032     $          - dth(i,k-1))
     1033        ENDIF
     1034      ENDDO
     1035      ENDDO
     1036c
     1037c
     1038      DO i=1,klon
     1039      ptop(i) = ptop_new(i)
     1040      ENDDO
     1041
     1042      DO k=klev,1,-1
     1043      DO i=1,klon
     1044        IF (ph(i,k+1) .LT. ptop(i)) ktop(i)=k
     1045      ENDDO
     1046      ENDDO
     1047c
     1048c      5/ Set deltatw & deltaqw to 0 above kupper
     1049c
     1050      DO k = 1,klev
     1051      DO i=1,klon
     1052        IF (k .GE. kupper(i)) THEN
     1053         deltatw(i,k) = 0.
     1054         deltaqw(i,k) = 0.
     1055        ENDIF
     1056      ENDDO
     1057      ENDDO
     1058c
     1059C
     1060       ENDDO      ! end sub-timestep loop
     1061C
     1062C -----------------------------------------------------------------
     1063c   Get back to tendencies per second
     1064c
     1065      DO k = 1,klev
     1066      DO i=1,klon
     1067        IF (k .LE. kupper(i)-1) THEN
     1068         dtls(i,k) = dtls(i,k)/dtime
     1069         dqls(i,k) = dqls(i,k)/dtime
     1070         d_deltatw2(i,k)=d_deltatw2(i,k)/dtime
     1071         d_deltaqw2(i,k)=d_deltaqw2(i,k)/dtime
     1072         d_deltat_gw(i,k) = d_deltat_gw(i,k)/dtime
     1073        ENDIF
     1074      ENDDO
     1075      ENDDO
     1076c
     1077c
     1078c----------------------------------------------------------
     1079c   Determine wake final state; recompute wape, cstar, ktop;
     1080c   filter out bad wakes.
     1081c----------------------------------------------------------
     1082c
     1083C 2.1 - Undisturbed area and Wake integrals
     1084C ---------------------------------------------------------
     1085
     1086      DO i=1,klon
     1087        z(i) = 0.
     1088        sum_thu(i) = 0.
     1089        sum_tu(i) = 0.
     1090        sum_qu(i) = 0.
     1091        sum_thvu(i) = 0.
     1092        sum_dth(i) = 0.
     1093        sum_dq(i) = 0.
     1094        sum_rho(i) = 0.
     1095        sum_dtdwn(i) = 0.
     1096        sum_dqdwn(i) = 0.
     1097
     1098        av_thu(i) = 0.
     1099        av_tu(i) =0.
     1100        av_qu(i) =0.
     1101        av_thvu(i) = 0.
     1102        av_dth(i) = 0.
     1103        av_dq(i) = 0.
     1104        av_rho(i) =0.
     1105        av_dtdwn(i) =0.
     1106        av_dqdwn(i) = 0.
     1107      ENDDO
     1108C Potential temperatures and humidity
     1109c----------------------------------------------------------
     1110
     1111      DO k =1,klev
     1112      DO i=1,klon
     1113        rho(i,k) = p(i,k)/(rd*te(i,k))
     1114        IF(k .eq. 1) THEN
     1115          rhoh(i,k) = ph(i,k)/(rd*te(i,k))
     1116          zhh(i,k)=0
     1117        ELSE
     1118          rhoh(i,k) = ph(i,k)*2./(rd*(te(i,k)+te(i,k-1)))
     1119          zhh(i,k)=(ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*RG)+zhh(i,k-1)
     1120        ENDIF
     1121        the(i,k) = te(i,k)/ppi(i,k)
     1122        thu(i,k) = (te(i,k) - deltatw(i,k)*sigmaw(i))/ppi(i,k)
     1123        tu(i,k) = te(i,k) - deltatw(i,k)*sigmaw(i)
     1124        qu(i,k)  =  qe(i,k) - deltaqw(i,k)*sigmaw(i)
     1125        rhow(i,k) = p(i,k)/(rd*(te(i,k)+deltatw(i,k)))
     1126        dth(i,k) = deltatw(i,k)/ppi(i,k)
     1127      ENDDO
     1128      ENDDO
     1129
     1130C Integrals (and wake top level number)
     1131C -----------------------------------------------------------
     1132
     1133C Initialize sum_thvu to 1st level virt. pot. temp.
     1134
     1135      DO i=1,klon
     1136        z(i) = 1.
     1137        dz(i) = 1.
     1138        sum_thvu(i) =  thu(i,1)*(1.+eps*qu(i,1))*dz(i)
     1139        sum_dth(i) = 0.
     1140      ENDDO
     1141
     1142      DO k = 1,klev
     1143      DO i=1,klon
     1144        dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
     1145        IF (dz(i) .GT. 0) THEN
     1146         z(i) = z(i)+dz(i)
     1147         sum_thu(i) = sum_thu(i) + thu(i,k)*dz(i)
     1148         sum_tu(i) = sum_tu(i) + tu(i,k)*dz(i)
     1149         sum_qu(i) = sum_qu(i) + qu(i,k)*dz(i)
     1150         sum_thvu(i) = sum_thvu(i) + thu(i,k)*(1.+eps*qu(i,k))*dz(i)
     1151         sum_dth(i) = sum_dth(i) + dth(i,k)*dz(i)
     1152         sum_dq(i) = sum_dq(i) + deltaqw(i,k)*dz(i)
     1153         sum_rho(i) = sum_rho(i) + rhow(i,k)*dz(i)
     1154         sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i,k)*dz(i)
     1155         sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i,k)*dz(i)
     1156        ENDIF
     1157      ENDDO
     1158      ENDDO
     1159c
     1160      DO i=1,klon
     1161        hw0(i) = z(i)
     1162      ENDDO
     1163c
     1164C 2.1 - WAPE and mean forcing computation
     1165C-------------------------------------------------------------
     1166
     1167C Means
     1168
     1169      DO i=1, klon
     1170        av_thu(i) = sum_thu(i)/hw0(i)
     1171        av_tu(i) = sum_tu(i)/hw0(i)
     1172        av_qu(i) = sum_qu(i)/hw0(i)
     1173        av_thvu(i) = sum_thvu(i)/hw0(i)
     1174        av_dth(i) = sum_dth(i)/hw0(i)
     1175        av_dq(i) = sum_dq(i)/hw0(i)
     1176        av_rho(i) = sum_rho(i)/hw0(i)
     1177        av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
     1178        av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
     1179
     1180        wape2(i) = - rg*hw0(i)*(av_dth(i)
     1181     $     + eps*(av_thu(i)*av_dq(i)+av_dth(i)*av_qu(i)+
     1182     $     av_dth(i)*av_dq(i) ))/av_thvu(i)
     1183      ENDDO
     1184
     1185C 2.2 Prognostic variable update
     1186C ------------------------------------------------------------
     1187
     1188C Filter out bad wakes
     1189c
     1190      DO k = 1,klev
     1191      DO i=1,klon
     1192        IF ( wape2(i) .LT. 0.) THEN
     1193          deltatw(i,k) = 0.
     1194          deltaqw(i,k) = 0.
     1195          dth(i,k) = 0.
     1196        ENDIF
     1197      ENDDO
     1198      ENDDO
     1199c
     1200
     1201      DO i=1, klon
     1202      IF ( wape2(i) .LT. 0.) THEN
     1203        wape2(i) = 0.
     1204        Cstar2(i) = 0.
     1205        hw(i) = hwmin
     1206        sigmaw(i) = amax1(sigmad,sigd_con(i))
     1207        fip(i) = 0.
     1208        gwake(i) = .FALSE.
     1209      ELSE
     1210        if(prt_level.ge.10) print*,'wape2>0'
     1211        Cstar2(i) = stark*sqrt(2.*wape2(i))
     1212        gwake(i) = .TRUE.
     1213      ENDIF
     1214      ENDDO
     1215c
     1216      DO i=1, klon
     1217        ktopw(i) = ktop(i)
     1218      ENDDO
     1219c
     1220      DO i=1, klon
     1221      IF (ktopw(i) .gt. 0 .and. gwake(i)) then
     1222
     1223Cjyg1     Utilisation d'un h_efficace constant ( ~ feeding layer)
     1224ccc       heff = 600.
     1225C      Utilisation de la hauteur hw
     1226cc       heff = 0.7*hw
     1227       heff(i) = hw(i)
     1228
     1229       FIP(i) = 0.5*rho(i,ktopw(i))*Cstar2(i)**3*heff(i)*2*
     1230     $      sqrt(sigmaw(i)*wdens*3.14)
     1231       FIP(i) = alpk * FIP(i)
     1232Cjyg2
     1233       ELSE
     1234         FIP(i) = 0.
     1235       ENDIF
     1236      ENDDO
     1237c
     1238C   Limitation de sigmaw
     1239c
     1240C   sécurité : si le wake occuppe plus de 90 % de la surface de la maille,
     1241C              alors il disparait en se mélangeant à la partie undisturbed
     1242c
     1243      DO k = 1,klev
     1244       DO i=1, klon
     1245         IF ((sigmaw(i).GT.0.9).or.
     1246     $     ((wape(i).ge.wape2(i)).and.(wape2(i).le.1.0))) THEN
     1247ccc      IF (sigmaw(i).GT.0.9) THEN
     1248          dtls(i,k) = 0.
     1249          dqls(i,k) = 0.
     1250          deltatw(i,k) = 0.
     1251          deltaqw(i,k) = 0.
     1252        ENDIF
     1253       ENDDO
     1254      ENDDO
     1255c
     1256      DO i=1, klon
     1257         IF ((sigmaw(i).GT.0.9).or.
     1258     $     ((wape(i).ge.wape2(i)).and.(wape2(i).le.1.0))) THEN
     1259ccc      IF (sigmaw(i).GT.0.9) THEN
     1260         wape(i) = 0.
     1261         hw(i) = hwmin
     1262         sigmaw(i) = sigmad
     1263         fip(i) = 0.
     1264        ELSE
     1265         wape(i) = wape2(i)
     1266        ENDIF
     1267      ENDDO
     1268c
     1269c
     1270      RETURN
     1271      END
     1272      Subroutine WAKE_scal (p,ph,ppi,dtime,sigd_con
     1273     :                ,te0,qe0,omgb
     1274     :                ,dtdwn,dqdwn,amdwn,amup,dta,dqa
     1275     :                ,wdtPBL,wdqPBL,udtPBL,udqPBL
     1276     o                ,deltatw,deltaqw,dth,hw,sigmaw,wape,fip,gfl
     1277     o                ,dtls,dqls
     1278     o                ,ktopw,omgbdth,dp_omgb,wdens
     1279     o                ,tu,qu
     1280     o                ,dtKE,dqKE
     1281     o                ,dtPBL,dqPBL
     1282     o                ,omg,dp_deltomg,spread
     1283     o                ,Cstar,d_deltat_gw
     1284     o                ,d_deltatw2,d_deltaqw2)
     1285
     1286***************************************************************
     1287*                                                             *
     1288* WAKE                                                        *
     1289*      retour a un Pupper fixe                                *
     1290*                                                             *
     1291* written by   :  GRANDPEIX Jean-Yves   09/03/2000            *
     1292* modified by :   ROEHRIG Romain        01/29/2007            *
     1293***************************************************************
     1294c
     1295      USE dimphy
     1296      IMPLICIT none
     1297c============================================================================
     1298C
     1299C
     1300C   But : Decrire le comportement des poches froides apparaissant dans les
     1301C        grands systemes convectifs, et fournir l'energie disponible pour
     1302C        le declenchement de nouvelles colonnes convectives.
     1303C
     1304C   Variables d'etat : deltatw    : ecart de temperature wake-undisturbed area
     1305C                      deltaqw    : ecart d'humidite wake-undisturbed area
     1306C                      sigmaw     : fraction d'aire occupee par la poche.
     1307C
     1308C   Variable de sortie :
     1309c
     1310c                        wape : WAke Potential Energy
     1311c                        fip  : Front Incident Power (W/m2) - ALP
     1312c                        gfl  : Gust Front Length per unit area (m-1)
     1313C                        dtls : large scale temperature tendency due to wake
     1314C                        dqls : large scale humidity tendency due to wake
     1315C                        hw   : hauteur de la poche
     1316C                     dp_omgb : vertical gradient of large scale omega
     1317C                      omgbdth: flux of Delta_Theta transported by LS omega
     1318C                      dtKE   : differential heating (wake - unpertubed)
     1319C                      dqKE   : differential moistening (wake - unpertubed)
     1320C                      omg    : Delta_omg =vertical velocity diff. wake-undist. (Pa/s)
     1321C                 dp_deltomg  : vertical gradient of omg (s-1)
     1322C                     spread  : spreading term in dt_wake and dq_wake
     1323C                 deltatw     : updated temperature difference (T_w-T_u).
     1324C                 deltaqw     : updated humidity difference (q_w-q_u).
     1325C                 sigmaw      : updated wake fractional area.
     1326C                 d_deltat_gw : delta T tendency due to GW
     1327c
     1328C   Variables d'entree :
     1329c
     1330c                        aire : aire de la maille
     1331c                        te0  : temperature dans l'environnement  (K)
     1332C                        qe0  : humidite dans l'environnement     (kg/kg)
     1333C                        omgb : vitesse verticale moyenne sur la maille (Pa/s)
     1334C                        dtdwn: source de chaleur due aux descentes (K/s)
     1335C                        dqdwn: source d'humidite due aux descentes (kg/kg/s)
     1336C                        dta  : source de chaleur due courants satures et detrain  (K/s)
     1337C                        dqa  : source d'humidite due aux courants satures et detra (kg/kg/s)
     1338C                        amdwn: flux de masse total des descentes, par unite de
     1339C                                surface de la maille (kg/m2/s)
     1340C                        amup : flux de masse total des ascendances, par unite de
     1341C                                surface de la maille (kg/m2/s)
     1342C                        p    : pressions aux milieux des couches (Pa)
     1343C                        ph   : pressions aux interfaces (Pa)
     1344C                        ppi  : (p/p_0)**kapa (adim)
     1345C                        dtime: increment temporel (s)
     1346c
     1347C   Variables internes :
     1348c
     1349c                        rhow : masse volumique de la poche froide
     1350C                        rho  : environment density at P levels
     1351C                        rhoh : environment density at Ph levels
     1352C                        te   : environment temperature | may change within
     1353C                        qe   : environment humidity    | sub-time-stepping
     1354C                        the  : environment potential temperature
     1355C                        thu  : potential temperature in undisturbed area
     1356C                        tu   :  temperature  in undisturbed area
     1357C                        qu   : humidity in undisturbed area
     1358C                      dp_omgb: vertical gradient og LS omega
     1359C                      omgbw  : wake average vertical omega
     1360C                     dp_omgbw: vertical gradient of omgbw
     1361C                      omgbdq : flux of Delta_q transported by LS omega
     1362C                        dth  : potential temperature diff. wake-undist.
     1363C                        th1  : first pot. temp. for vertical advection (=thu)
     1364C                        th2  : second pot. temp. for vertical advection (=thw)
     1365C                        q1   : first humidity for vertical advection
     1366C                        q2   : second humidity for vertical advection
     1367C                     d_deltatw   : terme de redistribution pour deltatw
     1368C                     d_deltaqw   : terme de redistribution pour deltaqw
     1369C                      deltatw0   : deltatw initial
     1370C                      deltaqw0   : deltaqw initial
     1371C                      hw0    : hw initial
     1372C                      sigmaw0: sigmaw initial
     1373C                      amflux : horizontal mass flux through wake boundary
     1374C                      wdens  : number of wakes per unit area (3D) or per
     1375C                               unit length (2D)
     1376C                      Tgw    : 1 sur la période de onde de gravité
     1377c                      Cgw    : vitesse de propagation de onde de gravité
     1378c                      LL     : distance entre 2 poches
     1379
     1380c-------------------------------------------------------------------------
     1381c          Déclaration de variables
     1382c-------------------------------------------------------------------------
     1383
     1384#include "dimensions.h"
    1141385cccc#include "dimphy.h"
    1151386#include "YOMCST.h"
Note: See TracChangeset for help on using the changeset viewer.