source: LMDZ4/trunk/libf/phylmd/wake.F @ 942

Last change on this file since 942 was 940, checked in by Laurent Fairhead, 17 years ago

On remplace le fichier include dimphy.h par le module dimphy.F90i pour etre
coherent avec le partout
LF
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 30.1 KB
Line 
1      Subroutine WAKE (p,ph,ppi,dtime,sigd_con
2     :                ,te0,qe0,omgb,ibas
3     :                ,dtdwn,dqdwn,amdwn,amup,dta,dqa
4     :                ,wdtPBL,wdqPBL,udtPBL,udqPBL
5     o                ,deltatw,deltaqw,dth,hw,sigmaw,wape,fip,gfl
6     o                ,dtls,dqls
7     o                ,ktopw,omgbdth,dp_omgb,wdens
8     o                ,tu,qu
9     o                ,dtKE,dqKE
10     o                ,dtPBL,dqPBL
11     o                ,omg,dp_deltomg,spread
12     o                ,Cstar,d_deltat_gw
13     o                ,d_deltatw2,d_deltaqw2)
14
15***************************************************************
16*                                                             *
17* WAKE                                                        *
18*      retour a un Pupper fixe                                *
19*                                                             *
20* written by   :  GRANDPEIX Jean-Yves   09/03/2000            *
21* modified by :   ROEHRIG Romain        01/29/2007            *
22***************************************************************
23c
24      USE dimphy
25      IMPLICIT none
26c============================================================================
27C
28C
29C   But : Decrire le comportement des poches froides apparaissant dans les
30C        grands systemes convectifs, et fournir l'energie disponible pour
31C        le declenchement de nouvelles colonnes convectives.
32C
33C   Variables d'etat : deltatw    : ecart de temperature wake-undisturbed area
34C                      deltaqw    : ecart d'humidite wake-undisturbed area
35C                      sigmaw     : fraction d'aire occupee par la poche.
36C
37C   Variable de sortie :
38c
39c                        wape : WAke Potential Energy
40c                        fip  : Front Incident Power (W/m2) - ALP
41c                        gfl  : Gust Front Length per unit area (m-1)
42C                        dtls : large scale temperature tendency due to wake
43C                        dqls : large scale humidity tendency due to wake
44C                        hw   : hauteur de la poche
45C                     dp_omgb : vertical gradient of large scale omega
46C                      omgbdth: flux of Delta_Theta transported by LS omega
47C                      dtKE   : differential heating (wake - unpertubed)
48C                      dqKE   : differential moistening (wake - unpertubed)
49C                      omg    : Delta_omg =vertical velocity diff. wake-undist. (Pa/s)
50C                 dp_deltomg  : vertical gradient of omg (s-1)
51C                     spread  : spreading term in dt_wake and dq_wake
52C                 deltatw     : updated temperature difference (T_w-T_u).
53C                 deltaqw     : updated humidity difference (q_w-q_u).
54C                 sigmaw      : updated wake fractional area.
55C                 d_deltat_gw : delta T tendency due to GW
56c
57C   Variables d'entree :
58c
59c                        aire : aire de la maille
60c                        te0  : temperature dans l'environnement  (K)
61C                        qe0  : humidite dans l'environnement     (kg/kg)
62C                        omgb : vitesse verticale moyenne sur la maille (Pa/s)
63C                        ibas : cloud base level number
64C                        dtdwn: source de chaleur due aux descentes (K/s)
65C                        dqdwn: source d'humidite due aux descentes (kg/kg/s)
66C                        dta  : source de chaleur due courants satures et detrain  (K/s)
67C                        dqa  : source d'humidite due aux courants satures et detra (kg/kg/s)
68C                        amdwn: flux de masse total des descentes, par unite de
69C                                surface de la maille (kg/m2/s)
70C                        amup : flux de masse total des ascendances, par unite de
71C                                surface de la maille (kg/m2/s)
72C                        p    : pressions aux milieux des couches (Pa)
73C                        ph   : pressions aux interfaces (Pa)
74C                        ppi  : (p/p_0)**kapa (adim)
75C                        dtime: increment temporel (s)
76c
77C   Variables internes :
78c
79c                        rhow : masse volumique de la poche froide
80C                        rho  : environment density at P levels
81C                        rhoh : environment density at Ph levels
82C                        te   : environment temperature | may change within
83C                        qe   : environment humidity    | sub-time-stepping
84C                        the  : environment potential temperature
85C                        thu  : potential temperature in undisturbed area
86C                        tu   :  temperature  in undisturbed area
87C                        qu   : humidity in undisturbed area
88C                      dp_omgb: vertical gradient og LS omega
89C                      omgbw  : wake average vertical omega
90C                     dp_omgbw: vertical gradient of omgbw
91C                      omgbdq : flux of Delta_q transported by LS omega
92C                        dth  : potential temperature diff. wake-undist.
93C                        th1  : first pot. temp. for vertical advection (=thu)
94C                        th2  : second pot. temp. for vertical advection (=thw)
95C                        q1   : first humidity for vertical advection
96C                        q2   : second humidity for vertical advection
97C                     d_deltatw   : terme de redistribution pour deltatw
98C                     d_deltaqw   : terme de redistribution pour deltaqw
99C                      deltatw0   : deltatw initial
100C                      deltaqw0   : deltaqw initial
101C                      hw0    : hw initial
102C                      sigmaw0: sigmaw initial
103C                      amflux : horizontal mass flux through wake boundary
104C                      wdens  : number of wakes per unit area (3D) or per
105C                               unit length (2D)
106C                      Tgw    : 1 sur la période de onde de gravité
107c                      Cgw    : vitesse de propagation de onde de gravité
108c                      LL     : distance entre 2 poches
109
110c-------------------------------------------------------------------------
111c          Déclaration de variables
112c-------------------------------------------------------------------------
113
114#include "dimensions.h"
115cccc#include "dimphy.h"
116#include "YOMCST.h"
117#include "cvthermo.h"
118
119c Arguments en entree
120c--------------------
121
122      REAL p(klev),ph(klev+1),ppi(klev)
123      REAL dtime
124      REAL te0(klev),qe0(klev)
125      REAL omgb(klev+1)
126      INTEGER ibas
127      REAL dtdwn(klev), dqdwn(klev)
128      REAL wdtPBL(klev),wdqPBL(klev)
129      REAL udtPBL(klev),udqPBL(klev)
130      REAL amdwn(klev), amup(klev)
131      REAL dta(klev), dqa(klev)
132      REAL sigd_con
133
134c Sorties
135c--------
136
137      REAL deltatw(klev), deltaqw(klev), dth(klev)
138      REAL tu(klev), qu(klev)
139      REAL dtls(klev), dqls(klev)
140      REAL dtKE(klev), dqKE(klev)
141      REAL dtPBL(klev), dqPBL(klev)
142      REAL spread(klev)
143      REAL d_deltatgw(klev)
144      REAL d_deltatw2(klev), d_deltaqw2(klev)
145      REAL omgbdth(klev+1), omg(klev+1)
146      REAL dp_omgb(klev), dp_deltomg(klev)
147      REAL d_deltat_gw(klev)
148      REAL hw, sigmaw, wape, fip, gfl, Cstar
149      INTEGER ktopw
150
151c Variables internes
152c-------------------
153
154c Variables à fixer
155      REAL ALON
156      REAL coefgw
157      REAL wdens0, wdens
158      REAL stark
159      REAL alpk
160      REAL delta_t_min
161      REAL Pupper
162      INTEGER nsub
163      REAL dtimesub
164      REAL sigmad, hwmin
165
166c Variables de sauvegarde
167      REAL deltatw0(klev)
168      REAL deltaqw0(klev)
169      REAL te(klev), qe(klev)
170      REAL sigmaw0, sigmaw1
171
172c Variables pour les GW
173      REAL LL
174      REAL N2(klev)
175      REAL Cgw(klev)
176      REAL Tgw(klev)
177
178c Variables liées au calcul de hw
179      REAL ptop_provis, ptop, ptop_new
180      REAL sum_dth
181      REAL dthmin
182      REAL z, dz, hw0
183      INTEGER ktop, kupper
184
185c Autres variables internes
186      INTEGER isubstep, k
187
188      REAL sum_thu, sum_tu, sum_qu,sum_thvu
189      REAL sum_dq, sum_rho
190      REAL sum_dtdwn, sum_dqdwn
191      REAL av_thu, av_tu, av_qu, av_thvu
192      REAL av_dth, av_dq, av_rho
193      REAL av_dtdwn, av_dqdwn
194
195      REAL rho(klev), rhoh(klev+1), rhow(klev)
196      REAL rhow_moyen(klev)
197      REAL zh(klev), zhh(klev+1)
198      REAL epaisseur1(klev), epaisseur2(klev)
199
200      REAL pbase
201     
202      REAL the(klev), thu(klev)
203
204      REAL d_deltatw(klev), d_deltaqw(klev)
205
206      REAL omgbw(klev+1), omgtop
207      REAL dp_omgbw(klev)
208      REAL ztop, dztop
209      REAL alpha_up(klev)
210     
211      REAL RRe1, RRe2, RRd1, RRd2
212      REAL Th1(klev), Th2(klev), q1(klev), q2(klev)
213      REAL D_Th1(klev), D_Th2(klev), D_dth(klev)
214      REAL D_q1(klev), D_q2(klev), D_dq(klev)
215      REAL omgbdq(klev)
216
217      REAL ff, gg
218      REAL wape2, Cstar2, heff
219
220      REAL Crep(klev)
221      REAL Crep_upper, Crep_sol
222
223C-------------------------------------------------------------------------
224c         Initialisations
225c-------------------------------------------------------------------------
226
227c      print*, 'wake initialisations'
228
229c   Essais d'initialisation avec sigmaw = 0.02 et hw = 10.
230c-------------------------------------------------------------------------
231
232      DATA sigmad, hwmin /.02,10./
233
234C Longueur de maille (en m)
235c-------------------------------------------------------------------------
236
237c      ALON = 3.e5
238      ALON = 1.e6
239
240
241C Configuration de coefgw,stark,wdens (22/02/06 by YU Jingmei)
242c
243c      coefgw : Coefficient pour les ondes de gravité
244c       stark : Coefficient k dans Cstar=k*sqrt(2*WAPE)
245c       wdens : Densité de poche froide par maille
246c-------------------------------------------------------------------------
247
248      coefgw=10
249c      coefgw=1
250c      wdens0 = 1.0/(alon**2)   
251      wdens = 1.0/(alon**2)       
252      stark = 0.50
253cCRtest
254      alpk=0.1
255c      alpk = 1.0
256c      alpk = 0.5
257c      alpk = 0.05
258      Crep_upper=0.9
259      Crep_sol=1.0
260
261
262C Minimum value for |T_wake - T_undist|. Used for wake top definition
263c-------------------------------------------------------------------------
264
265      delta_t_min = 0.2
266
267
268C Cloud base
269c-------------------------------------------------------------------------
270
271      Pbase = P(ibas)
272
273
274C 1. - Save initial values and initialize tendencies
275C --------------------------------------------------
276
277      DO k=1,klev
278        deltatw0(k) = deltatw(k)
279        deltaqw0(k)= deltaqw(k)
280        te(k) = te0(k)
281        qe(k) = qe0(k)
282        dtls(k) = 0.
283        dqls(k) = 0.
284        d_deltat_gw(k)=0.
285        d_deltatw2(k)=0.
286        d_deltaqw2(k)=0.
287      ENDDO
288c      sigmaw1=sigmaw
289c      IF (sigd_con.GT.sigmaw1) THEN
290c      print*, 'sigmaw,sigd_con', sigmaw, sigd_con
291c      ENDIF
292      sigmaw = max(sigmaw,sigd_con)
293      sigmaw = max(sigmaw,sigmad)
294      sigmaw = min(sigmaw,0.99)
295      sigmaw0 = sigmaw
296c      wdens=wdens0/(10.*sigmaw)
297c      IF (sigd_con.GT.sigmaw1) THEN
298c      print*, 'sigmaw1,sigd1', sigmaw, sigd_con
299c      ENDIF
300
301C 2. - Prognostic part
302C =========================================================
303
304c      print *, 'prognostic wake computation'
305
306
307C 2.1 - Undisturbed area and Wake integrals
308C ---------------------------------------------------------
309
310      z = 0.
311      ktop=0
312      kupper = 0
313      sum_thu = 0.
314      sum_tu = 0.
315      sum_qu = 0.
316      sum_thvu = 0.
317      sum_dth = 0.
318      sum_dq = 0.
319      sum_rho = 0.
320      sum_dtdwn = 0.
321      sum_dqdwn = 0.
322
323      av_thu = 0.
324      av_tu =0.
325      av_qu =0.
326      av_thvu = 0.
327      av_dth = 0.
328      av_dq = 0.
329      av_rho =0.
330      av_dtdwn =0.
331      av_dqdwn = 0.
332
333C Potential temperatures and humidity
334c----------------------------------------------------------
335
336      DO k =1,klev
337        rho(k) = p(k)/(rd*te(k))
338        IF(k .eq. 1) THEN
339          rhoh(k) = ph(k)/(rd*te(k))
340          zhh(k)=0
341        ELSE
342          rhoh(k) = ph(k)*2./(rd*(te(k)+te(k-1)))
343          zhh(k)=(ph(k)-ph(k-1))/(-rhoh(k)*RG)+zhh(k-1)
344        ENDIF
345        the(k) = te(k)/ppi(k)
346        thu(k) = (te(k) - deltatw(k)*sigmaw)/ppi(k)
347        tu(k) = te(k) - deltatw(k)*sigmaw
348        qu(k)  =  qe(k) - deltaqw(k)*sigmaw
349        rhow(k) = p(k)/(rd*(te(k)+deltatw(k)))
350        dth(k) = deltatw(k)/ppi(k)
351        LL = (1-sqrt(sigmaw))/sqrt(wdens)       
352      ENDDO
353       
354      DO k = 1, klev-1
355        IF(k.eq.1) THEN
356          N2(k)=0
357        ELSE
358          N2(k)=max(0.,-RG**2/the(k)*rho(k)*(the(k+1)-the(k-1))
359     $           /(p(k+1)-p(k-1)))
360        ENDIF
361        ZH(k)=(zhh(k)+zhh(k+1))/2
362
363        Cgw(k)=sqrt(N2(k))*ZH(k)
364        Tgw(k)=coefgw*Cgw(k)/LL
365      ENDDO
366         
367      N2(klev)=0
368      ZH(klev)=0
369      Cgw(klev)=0
370      Tgw(klev)=0
371
372c  Calcul de la masse volumique moyenne de la colonne
373c-----------------------------------------------------------------
374
375      DO k=1,klev
376        epaisseur1(k)=0.
377        epaisseur2(k)=0.
378      ENDDO
379
380      epaisseur1(1)= -(Ph(2)-Ph(1))/(rho(1)*rg)+1.
381      epaisseur2(1)= -(Ph(2)-Ph(1))/(rho(1)*rg)+1.
382      rhow_moyen(1) = rhow(1)
383
384      DO k = 2, klev
385        epaisseur1(k)= -(Ph(k+1)-Ph(k))/(rho(k)*rg) +1.
386        epaisseur2(k)=epaisseur2(k-1)+epaisseur1(k)
387        rhow_moyen(k) = (rhow_moyen(k-1)*epaisseur2(k-1)+
388     $                 rhow(k)*epaisseur1(k))/epaisseur2(k)
389      ENDDO
390
391
392C Choose an integration bound well above wake top
393c-----------------------------------------------------------------
394
395c       Pupper = 50000.  ! melting level
396       Pupper = 60000.
397c       Pupper = 70000.
398
399
400C    Determine Wake top pressure (Ptop) from buoyancy integral
401C-----------------------------------------------------------------
402
403c-1/ Pressure of the level where dth becomes less than delta_t_min.
404
405      Ptop_provis=ph(1)
406      DO k= 2,klev
407        IF (dth(k) .GT. -delta_t_min .and.
408     $      dth(k-1).LT. -delta_t_min) THEN
409          Ptop_provis = ((dth(k)+delta_t_min)*p(k-1)
410     $          - (dth(k-1)+delta_t_min)*p(k)) /(dth(k) - dth(k-1))
411          GO TO 25
412        ENDIF
413      ENDDO
41425    CONTINUE
415
416c-2/ dth integral
417
418      sum_dth = 0.
419      dthmin = -delta_t_min
420      z = 0.
421
422      DO k = 1,klev
423        dz = -(max(ph(k+1),Ptop_provis)-Ph(k))/(rho(k)*rg)
424        IF (dz .le. 0) GO TO 40
425        z = z+dz
426        sum_dth = sum_dth + dth(k)*dz
427        dthmin = min(dthmin,dth(k))
428      ENDDO
42940    CONTINUE
430
431c-3/ height of triangle with area= sum_dth and base = dthmin
432
433      hw0 = 2.*sum_dth/min(dthmin,-0.5)
434      hw0 = max(hwmin,hw0)
435
436c-4/ now, get Ptop
437
438      z = 0.
439      ptop = ph(1)
440
441      DO k = 1,klev
442        dz = min(-(ph(k+1)-Ph(k))/(rho(k)*rg),hw0-z)
443        IF (dz .le. 0) GO TO 45
444        z = z+dz
445        Ptop = Ph(k)-rho(k)*rg*dz
446      ENDDO
44745    CONTINUE
448
449
450C-5/ Determination de ktop et kupper
451
452      DO k=klev,1,-1
453        IF (ph(k+1) .lt. ptop) ktop=k
454        IF (ph(k+1) .lt. pupper) kupper=k
455      ENDDO
456
457c-6/ Correct ktop and ptop
458
459      Ptop_new=ptop
460      DO k= ktop,2,-1
461        IF (dth(k) .GT. -delta_t_min .and.
462     $      dth(k-1).LT. -delta_t_min) THEN
463          Ptop_new = ((dth(k)+delta_t_min)*p(k-1)
464     $          - (dth(k-1)+delta_t_min)*p(k)) /(dth(k) - dth(k-1))
465          GO TO 225
466        ENDIF
467      ENDDO
468225   CONTINUE
469
470      ptop = ptop_new
471
472      DO k=klev,1,-1
473        IF (ph(k+1) .lt. ptop) ktop=k
474      ENDDO
475
476c Set deltatw & deltaqw to 0 above kupper
477c-----------------------------------------------------------
478
479      DO k = kupper,klev
480        deltatw(k) = 0.
481        deltaqw(k) = 0.
482      ENDDO
483
484
485C Vertical gradient of LS omega
486C------------------------------------------------------------
487
488      DO k = 1,kupper
489        dp_omgb(k) = (omgb(k+1) - omgb(k))/(ph(k+1)-ph(k))
490      ENDDO
491
492
493C Integrals (and wake top level number)
494C -----------------------------------------------------------
495
496C Initialize sum_thvu to 1st level virt. pot. temp.
497
498      z = 1.
499      dz = 1.
500      sum_thvu =  thu(1)*(1.+eps*qu(1))*dz
501      sum_dth = 0.
502
503      DO k = 1,klev
504        dz = -(max(ph(k+1),Ptop)-Ph(k))/(rho(k)*rg)
505        IF (dz .LE. 0) GO TO 50
506        z = z+dz
507        sum_thu = sum_thu + thu(k)*dz
508        sum_tu = sum_tu + tu(k)*dz
509        sum_qu = sum_qu + qu(k)*dz
510        sum_thvu = sum_thvu + thu(k)*(1.+eps*qu(k))*dz
511        sum_dth = sum_dth + dth(k)*dz
512        sum_dq = sum_dq + deltaqw(k)*dz
513        sum_rho = sum_rho + rhow(k)*dz
514        sum_dtdwn = sum_dtdwn + dtdwn(k)*dz
515        sum_dqdwn = sum_dqdwn + dqdwn(k)*dz
516      ENDDO
51750    CONTINUE
518
519      hw0 = z
520
521C 2.1 - WAPE and mean forcing computation
522C-------------------------------------------------------------
523
524C Means
525
526      av_thu = sum_thu/hw0
527      av_tu = sum_tu/hw0
528      av_qu = sum_qu/hw0
529      av_thvu = sum_thvu/hw0
530c      av_thve = sum_thve/hw0
531      av_dth = sum_dth/hw0
532      av_dq = sum_dq/hw0
533      av_rho = sum_rho/hw0
534      av_dtdwn = sum_dtdwn/hw0
535      av_dqdwn = sum_dqdwn/hw0
536
537      wape = - rg*hw0*(av_dth
538     $     + eps*(av_thu*av_dq+av_dth*av_qu+av_dth*av_dq ))/av_thvu
539
540C 2.2 Prognostic variable update
541C ------------------------------------------------------------
542
543C Filter out bad wakes
544
545      IF ( wape .LT. 0.) THEN
546        print*,'wape<0'
547        wape = 0.
548        hw = hwmin
549        sigmaw = max(sigmad,sigd_con)
550        fip = 0.
551        DO k = 1,klev
552          deltatw(k) = 0.
553          deltaqw(k) = 0.
554          dth(k) = 0.
555        ENDDO
556      ELSE
557        print*,'wape>0'
558        Cstar = stark*sqrt(2.*wape)
559      ENDIF
560
561C------------------------------------------------------------------
562C    Sub-time-stepping
563C------------------------------------------------------------------
564
565c      nsub=36
566      nsub=10
567      dtimesub=dtime/nsub
568
569c------------------------------------------------------------
570      DO isubstep = 1,nsub
571c------------------------------------------------------------
572
573c        print*,'---------------','substep=',isubstep,'-------------'
574
575c  Evolution of sigmaw
576
577
578        gfl = 2.*sqrt(3.14*wdens*sigmaw)           
579
580        sigmaw =sigmaw + gfl*Cstar*dtimesub
581        sigmaw =min(sigmaw,0.99)     !!!!!!!!
582c        wdens = wdens0/(10.*sigmaw)
583c        sigmaw =max(sigmaw,sigd_con)
584c        sigmaw =max(sigmaw,sigmad)
585
586c calcul de la difference de vitesse verticale poche - zone non perturbee
587
588        z= 0.
589        dp_deltomg(1:klev)=0.
590        omg(1:klev+1)=0.
591
592        omg(1) = 0.
593        dp_deltomg(1) = -(gfl*Cstar)/(sigmaw * (1-sigmaw))
594
595        DO k=2,ktop
596          dz = -(Ph(k)-Ph(k-1))/(rho(k-1)*rg)
597          z = z+dz
598          dp_deltomg(k)= dp_deltomg(1)
599          omg(k)= dp_deltomg(1)*z
600        ENDDO
601
602        dztop=-(Ptop-Ph(ktop))/(rho(ktop)*rg)
603        ztop = z+dztop
604        omgtop=dp_deltomg(1)*ztop
605
606
607c Conversion de la vitesse verticale de m/s a Pa/s
608
609        omgtop = -rho(ktop)*rg*omgtop
610        dp_deltomg(1) = omgtop/(ptop-ph(1))
611
612        DO k = 1,ktop
613          omg(k) = - rho(k)*rg*omg(k)
614          dp_deltomg(k) = dp_deltomg(1)
615        ENDDO
616
617c   raccordement lineaire de omg de ptop a pupper
618
619      IF (kupper .GT. ktop) THEN
620        omg(kupper+1) = - Rg*amdwn(kupper+1)/sigmaw
621     $                + Rg*amup(kupper+1)/(1.-sigmaw)
622        dp_deltomg(kupper) = (omgtop-omg(kupper+1))/(Ptop-Pupper)
623        DO k=ktop+1,kupper
624          dp_deltomg(k) = dp_deltomg(kupper)
625          omg(k) = omgtop+(ph(k)-Ptop)*dp_deltomg(kupper)
626        ENDDO
627      ENDIF
628
629c   Compute wake average vertical velocity omgbw
630
631      DO k = 1,klev+1
632        omgbw(k) = omgb(k)+(1.-sigmaw)*omg(k)
633      ENDDO
634
635c  and its vertical gradient dp_omgbw
636
637      DO k = 1,klev
638        dp_omgbw(k) = (omgbw(k+1)-omgbw(k))/(ph(k+1)-ph(k))
639      ENDDO
640
641
642c  Upstream coefficients for omgb velocity
643c--    (alpha_up(k) is the coefficient of the value at level k)
644c--    (1-alpha_up(k) is the coefficient of the value at level k-1)
645
646      DO k = 1,klev
647       alpha_up(k) = 0.
648       IF (omgb(k) .GT. 0.) alpha_up(k) = 1.
649      ENDDO
650
651c  Matrix expressing [The,deltatw] from  [Th1,Th2]
652
653      RRe1 = 1.-sigmaw
654      RRe2 = sigmaw
655      RRd1 = -1.
656      RRd2 = 1.
657
658c Get [Th1,Th2], dth and [q1,q2]
659
660      DO k = 1,kupper+1
661        dth(k) = deltatw(k)/ppi(k)
662        Th1(k) = the(k) - sigmaw     *dth(k)   ! undisturbed area
663        Th2(k) = the(k) + (1.-sigmaw)*dth(k)   ! wake
664        q1(k) = qe(k) - sigmaw     *deltaqw(k) ! undisturbed area
665        q2(k) = qe(k) + (1.-sigmaw)*deltaqw(k) ! wake
666      ENDDO
667
668      D_Th1(1) = 0.
669      D_Th2(1) = 0.
670      D_dth(1) = 0.
671      D_q1(1) = 0.
672      D_q2(1) = 0.
673      D_dq(1) = 0.
674
675      DO k = 2,kupper+1 !   loop on interfaces
676        D_Th1(k) = Th1(k-1)-Th1(k)
677        D_Th2(k) = Th2(k-1)-Th2(k)
678        D_dth(k) = dth(k-1)-dth(k)
679        D_q1(k) = q1(k-1)-q1(k)
680        D_q2(k) = q2(k-1)-q2(k)
681        D_dq(k) = deltaqw(k-1)-deltaqw(k)
682      ENDDO
683
684      omgbdth(1) = 0.
685      omgbdq(1) = 0.
686
687      DO k = 2,kupper+1  !   loop on interfaces
688        omgbdth(k) = omgb(k)*(    dth(k-1) -     dth(k))
689        omgbdq(k)  = omgb(k)*(deltaqw(k-1) - deltaqw(k))
690      ENDDO
691
692
693c-----------------------------------------------------------------
694      DO k=1,kupper-1
695c-----------------------------------------------------------------
696c
697c   Compute redistribution (advective) term
698c
699         d_deltatw(k) =
700     $             dtimesub/(Ph(k)-Ph(k+1))*(
701     $       RRd1*omg(k  )*sigmaw     *D_Th1(k)
702     $      -RRd2*omg(k+1)*(1.-sigmaw)*D_Th2(k+1)
703     $      -(1.-alpha_up(k))*omgbdth(k) - alpha_up(k+1)*omgbdth(k+1)
704     $                      )*ppi(k)
705c         print*,'d_deltatw=',d_deltatw(k)
706c
707         d_deltaqw(k) =
708     $             dtimesub/(Ph(k)-Ph(k+1))*(
709     $       RRd1*omg(k  )*sigmaw     *D_q1(k)
710     $      -RRd2*omg(k+1)*(1.-sigmaw)*D_q2(k+1)
711     $      -(1.-alpha_up(k))*omgbdq(k) - alpha_up(k+1)*omgbdq(k+1)
712     $                      )
713c         print*,'d_deltaqw=',d_deltaqw(k)
714c
715c   and increment large scale tendencies
716c
717         dtls(k) = dtls(k) +
718     $               dtimesub*(
719     $        ( RRe1*omg(k  )*sigmaw     *D_Th1(k)
720     $         -RRe2*omg(k+1)*(1.-sigmaw)*D_Th2(k+1) )
721     $               /(Ph(k)-Ph(k+1))
722     $         -sigmaw*(1.-sigmaw)*dth(k)*dp_deltomg(k)
723     $                      )*ppi(k)
724c         print*,'dtls=',dtls(k)
725c
726         dqls(k) = dqls(k) +
727     $               dtimesub*(
728     $        ( RRe1*omg(k  )*sigmaw     *D_q1(k)
729     $         -RRe2*omg(k+1)*(1.-sigmaw)*D_q2(k+1) )
730     $               /(Ph(k)-Ph(k+1))
731     $         -sigmaw*(1.-sigmaw)*deltaqw(k)*dp_deltomg(k)
732     $                      )
733c         print*,'dqls=',dqls(k)
734
735c-------------------------------------------------------------------
736      ENDDO
737c------------------------------------------------------------------
738
739C   Increment state variables
740
741      DO k = 1,kupper-1
742
743c Coefficient de répartition
744
745        Crep(k)=Crep_sol*(ph(kupper)-ph(k))/(ph(kupper)-ph(1))
746        Crep(k)=Crep(k)+Crep_upper*(ph(1)-ph(k))/(p(1)-ph(kupper))
747       
748
749c Reintroduce compensating subsidence term.
750
751c        dtKE(k)=(dtdwn(k)*Crep(k))/sigmaw
752c        dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k))
753c     .                   /(1-sigmaw)
754c        dqKE(k)=(dqdwn(k)*Crep(k))/sigmaw
755c        dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k))
756c     .                   /(1-sigmaw)
757c
758c        dtKE(k)=(dtdwn(k)*Crep(k)+(1-Crep(k))*dta(k))/sigmaw
759c        dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k)*Crep(k))
760c     .                   /(1-sigmaw)
761c        dqKE(k)=(dqdwn(k)*Crep(k)+(1-Crep(k))*dqa(k))/sigmaw
762c        dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k)*Crep(k))
763c     .                   /(1-sigmaw)
764
765        dtKE(k)=(dtdwn(k)/sigmaw - dta(k)/(1.-sigmaw))
766        dqKE(k)=(dqdwn(k)/sigmaw - dqa(k)/(1.-sigmaw))
767c        print*,'dtKE=',dtKE(k)
768c        print*,'dqKE=',dqKE(k)
769c
770        dtPBL(k)=(wdtPBL(k)/sigmaw - udtPBL(k)/(1.-sigmaw))
771        dqPBL(k)=(wdqPBL(k)/sigmaw - udqPBL(k)/(1.-sigmaw))
772c
773        spread(k) = (1.-sigmaw)*dp_deltomg(k)+gfl*Cstar/sigmaw
774c        print*,'spread=',spread(k)
775
776
777c ajout d'un effet onde de gravité -Tgw(k)*deltatw(k) 03/02/06 YU Jingmei
778
779        d_deltat_gw(k)=d_deltat_gw(k)-Tgw(k)*deltatw(k)* dtimesub
780c        print*,'d_delta_gw=',d_deltat_gw(k)
781        ff=d_deltatw(k)/dtimesub
782
783c Sans GW
784c
785c        deltatw(k)=deltatw(k)+dtimesub*(ff+dtKE(k)-spread(k)*deltatw(k))
786c
787c GW formule 1
788c
789c        deltatw(k) = deltatw(k)+dtimesub*
790c     $         (ff+dtKE(k) - spread(k)*deltatw(k)-Tgw(k)*deltatw(k))
791c
792c GW formule 2
793
794        IF (dtimesub*Tgw(k).lt.1.e-10) THEN
795          deltatw(k) = deltatw(k)+dtimesub*
796     $          (ff+dtKE(k)+dtPBL(k)
797     $          - spread(k)*deltatw(k)-Tgw(k)*deltatw(k))
798        ELSE
799           deltatw(k) = deltatw(k)+1/Tgw(k)*(1-exp(-dtimesub*Tgw(k)))*
800     $          (ff+dtKE(k)+dtPBL(k)
801     $          - spread(k)*deltatw(k)-Tgw(k)*deltatw(k))
802        ENDIF
803   
804        dth(k) = deltatw(k)/ppi(k)
805
806        gg=d_deltaqw(k)/dtimesub
807
808       deltaqw(k) = deltaqw(k) +
809     $         dtimesub*(gg+ dqKE(k)+dqPBL(k) - spread(k)*deltaqw(k))
810
811       d_deltatw2(k)=d_deltatw2(k)+d_deltatw(k)
812       d_deltaqw2(k)=d_deltaqw2(k)+d_deltaqw(k)
813      ENDDO
814
815C   And update large scale variables
816
817      DO k = 1,kupper
818        te(k) = te0(k) + dtls(k)
819        qe(k) = qe0(k) + dqls(k)
820        the(k) = te(k)/ppi(k)
821      ENDDO
822
823c     Determine Ptop from buoyancy integral
824c----------------------------------------------------------------------
825
826c-1/ Pressure of the level where dth changes sign.
827
828      Ptop_provis=ph(1)
829
830      DO k= 2,klev
831        IF (dth(k) .GT. -delta_t_min .and.
832     $      dth(k-1).LT. -delta_t_min) THEN
833          Ptop_provis = ((dth(k)+delta_t_min)*p(k-1)
834     $          - (dth(k-1)+delta_t_min)*p(k)) /(dth(k) - dth(k-1))
835        GO TO 65
836        ENDIF
837      ENDDO
83865    CONTINUE
839
840c-2/ dth integral
841
842      sum_dth = 0.
843      dthmin = -delta_t_min
844      z = 0.
845
846      DO k = 1,klev
847        dz = -(max(ph(k+1),Ptop_provis)-Ph(k))/(rho(k)*rg)
848        IF (dz .le. 0) GO TO 70
849        z = z+dz
850        sum_dth = sum_dth + dth(k)*dz
851        dthmin = min(dthmin,dth(k))
852      ENDDO
85370    CONTINUE
854
855c-3/ height of triangle with area= sum_dth and base = dthmin
856
857      hw = 2.*sum_dth/min(dthmin,-0.5)
858      hw = max(hwmin,hw)
859
860c-4/ now, get Ptop
861
862      ktop = 0
863      z=0.
864
865      DO k = 1,klev
866        dz = min(-(ph(k+1)-Ph(k))/(rho(k)*rg),hw-z)
867        IF (dz .le. 0) GO TO 75
868        z = z+dz
869        Ptop = Ph(k)-rho(k)*rg*dz
870        ktop = k
871      ENDDO
87275    CONTINUE
873
874c-5/Correct ktop and ptop
875
876      Ptop_new=ptop
877
878      DO k= ktop,2,-1
879        IF (dth(k) .GT. -delta_t_min .and.
880     $      dth(k-1).LT. -delta_t_min) THEN
881          Ptop_new = ((dth(k)+delta_t_min)*p(k-1)
882     $          - (dth(k-1)+delta_t_min)*p(k)) /(dth(k) - dth(k-1))
883          GO TO 275
884        ENDIF
885      ENDDO
886275   CONTINUE
887
888      ptop = ptop_new
889
890      DO k=klev,1,-1
891        IF (ph(k+1) .LT. ptop) ktop=k
892      ENDDO
893
894c-6/ Set deltatw & deltaqw to 0 above kupper
895
896      DO k = kupper,klev
897        deltatw(k) = 0.
898        deltaqw(k) = 0.
899      ENDDO
900
901c------------------------------------------------------------------
902      ENDDO      ! end sub-timestep loop
903C -----------------------------------------------------------------
904
905c   Get back to tendencies per second
906
907      DO k = 1,kupper-1
908        dtls(k) = dtls(k)/dtime
909        dqls(k) = dqls(k)/dtime
910        d_deltatw2(k)=d_deltatw2(k)/dtime
911        d_deltaqw2(k)=d_deltaqw2(k)/dtime
912        d_deltat_gw(k) = d_deltat_gw(k)/dtime
913      ENDDO
914
915C 2.1 - Undisturbed area and Wake integrals
916C ---------------------------------------------------------
917
918      z = 0.
919      sum_thu = 0.
920      sum_tu = 0.
921      sum_qu = 0.
922      sum_thvu = 0.
923      sum_dth = 0.
924      sum_dq = 0.
925      sum_rho = 0.
926      sum_dtdwn = 0.
927      sum_dqdwn = 0.
928
929      av_thu = 0.
930      av_tu =0.
931      av_qu =0.
932      av_thvu = 0.
933      av_dth = 0.
934      av_dq = 0.
935      av_rho =0.
936      av_dtdwn =0.
937      av_dqdwn = 0.
938
939C Potential temperatures and humidity
940c----------------------------------------------------------
941
942      DO k =1,klev
943        rho(k) = p(k)/(rd*te(k))
944        IF(k .eq. 1) THEN
945          rhoh(k) = ph(k)/(rd*te(k))
946          zhh(k)=0
947        ELSE
948          rhoh(k) = ph(k)*2./(rd*(te(k)+te(k-1)))
949          zhh(k)=(ph(k)-ph(k-1))/(-rhoh(k)*RG)+zhh(k-1)
950        ENDIF
951        the(k) = te(k)/ppi(k)
952        thu(k) = (te(k) - deltatw(k)*sigmaw)/ppi(k)
953        tu(k) = te(k) - deltatw(k)*sigmaw
954        qu(k)  =  qe(k) - deltaqw(k)*sigmaw
955        rhow(k) = p(k)/(rd*(te(k)+deltatw(k)))
956        dth(k) = deltatw(k)/ppi(k)
957       
958      ENDDO
959
960C Integrals (and wake top level number)
961C -----------------------------------------------------------
962
963C Initialize sum_thvu to 1st level virt. pot. temp.
964
965      z = 1.
966      dz = 1.
967      sum_thvu =  thu(1)*(1.+eps*qu(1))*dz
968      sum_dth = 0.
969
970      DO k = 1,klev
971        dz = -(max(ph(k+1),Ptop)-Ph(k))/(rho(k)*rg)
972
973        IF (dz .LE. 0) GO TO 51
974        z = z+dz
975        sum_thu = sum_thu + thu(k)*dz
976        sum_tu = sum_tu + tu(k)*dz
977        sum_qu = sum_qu + qu(k)*dz
978        sum_thvu = sum_thvu + thu(k)*(1.+eps*qu(k))*dz
979        sum_dth = sum_dth + dth(k)*dz
980        sum_dq = sum_dq + deltaqw(k)*dz
981        sum_rho = sum_rho + rhow(k)*dz
982        sum_dtdwn = sum_dtdwn + dtdwn(k)*dz
983        sum_dqdwn = sum_dqdwn + dqdwn(k)*dz
984      ENDDO
985 51   CONTINUE
986
987      hw0 = z
988
989C 2.1 - WAPE and mean forcing computation
990C-------------------------------------------------------------
991
992C Means
993
994      av_thu = sum_thu/hw0
995      av_tu = sum_tu/hw0
996      av_qu = sum_qu/hw0
997      av_thvu = sum_thvu/hw0
998      av_dth = sum_dth/hw0
999      av_dq = sum_dq/hw0
1000      av_rho = sum_rho/hw0
1001      av_dtdwn = sum_dtdwn/hw0
1002      av_dqdwn = sum_dqdwn/hw0
1003
1004      wape2 = - rg*hw0*(av_dth
1005     $     + eps*(av_thu*av_dq+av_dth*av_qu+av_dth*av_dq ))/av_thvu
1006
1007
1008C 2.2 Prognostic variable update
1009C ------------------------------------------------------------
1010
1011C Filter out bad wakes
1012
1013      IF ( wape2 .LT. 0.) THEN
1014        print*,'wape2<0'
1015        wape2 = 0.
1016        hw = hwmin
1017        sigmaw = max(sigmad,sigd_con)
1018        fip = 0.
1019        DO k = 1,klev
1020          deltatw(k) = 0.
1021          deltaqw(k) = 0.
1022          dth(k) = 0.
1023        ENDDO
1024      ELSE
1025        print*,'wape2>0'
1026        Cstar2 = stark*sqrt(2.*wape2)
1027
1028      ENDIF
1029
1030      ktopw = ktop
1031
1032      IF (ktopw .gt. 0) then
1033
1034Cjyg1     Utilisation d'un h_efficace constant ( ~ feeding layer)
1035ccc       heff = 600.
1036C      Utilisation de la hauteur hw
1037cc       heff = 0.7*hw
1038       heff = hw
1039
1040       FIP = 0.5*rho(ktopw)*Cstar2**3*heff*2*sqrt(sigmaw*wdens*3.14)
1041       FIP = alpk * FIP
1042Cjyg2
1043       ELSE
1044         FIP = 0.
1045       ENDIF
1046
1047
1048C   Limitation de sigmaw
1049c
1050C   sécurité : si le wake occuppe plus de 90 % de la surface de la maille,
1051C              alors il disparait en se mélangeant à la partie undisturbed
1052
1053      IF ((sigmaw.GT.0.9).or.
1054     .     ((wape.ge.wape2).and.(wape2.le.1.0))) THEN
1055c      IF (sigmaw.GT.0.9) THEN
1056        DO k = 1,klev
1057          dtls(k) = 0.
1058          dqls(k) = 0.
1059          deltatw(k) = 0.
1060          deltaqw(k) = 0.
1061        ENDDO
1062        wape = 0.
1063        hw = hwmin
1064        sigmaw = sigmad
1065        fip = 0.
1066      ENDIF
1067
1068      RETURN
1069      END
1070
1071
1072
Note: See TracBrowser for help on using the repository browser.