My Project
 All Classes Files Functions Variables Macros
wake.F
Go to the documentation of this file.
1 !
2 ! $Id: wake.F 1649 2012-08-22 15:32:25Z emillour $
3 !
4  Subroutine wake (p,ph,pi,dtime,sigd_con
5  : ,te0,qe0,omgb
6  : ,dtdwn,dqdwn,amdwn,amup,dta,dqa
7  : ,wdtpbl,wdqpbl,udtpbl,udqpbl
8  o ,deltatw,deltaqw,dth,hw,sigmaw,wape,fip,gfl
9  o ,dtls,dqls
10  o ,ktopw,omgbdth,dp_omgb,wdens
11  o ,tu,qu
12  o ,dtke,dqke
13  o ,dtpbl,dqpbl
14  o ,omg,dp_deltomg,spread
15  o ,cstar,d_deltat_gw
16  o ,d_deltatw2,d_deltaqw2)
17 
18 
19 ***************************************************************
20 * *
21 * WAKE *
22 * retour a un Pupper fixe *
23 * *
24 * written by : GRANDPEIX Jean-Yves 09/03/2000 *
25 * modified by : ROEHRIG Romain 01/29/2007 *
26 ***************************************************************
27 c
28  use dimphy
29  IMPLICIT none
30 c============================================================================
31 C
32 C
33 C But : Decrire le comportement des poches froides apparaissant dans les
34 C grands systemes convectifs, et fournir l'energie disponible pour
35 C le declenchement de nouvelles colonnes convectives.
36 C
37 C Variables d'etat : deltatw : ecart de temperature wake-undisturbed area
38 C deltaqw : ecart d'humidite wake-undisturbed area
39 C sigmaw : fraction d'aire occupee par la poche.
40 C
41 C Variable de sortie :
42 c
43 c wape : WAke Potential Energy
44 c fip : Front Incident Power (W/m2) - ALP
45 c gfl : Gust Front Length per unit area (m-1)
46 C dtls : large scale temperature tendency due to wake
47 C dqls : large scale humidity tendency due to wake
48 C hw : hauteur de la poche
49 C dp_omgb : vertical gradient of large scale omega
50 C wdens : densite de poches
51 C omgbdth: flux of Delta_Theta transported by LS omega
52 C dtKE : differential heating (wake - unpertubed)
53 C dqKE : differential moistening (wake - unpertubed)
54 C omg : Delta_omg =vertical velocity diff. wake-undist. (Pa/s)
55 C dp_deltomg : vertical gradient of omg (s-1)
56 C spread : spreading term in dt_wake and dq_wake
57 C deltatw : updated temperature difference (T_w-T_u).
58 C deltaqw : updated humidity difference (q_w-q_u).
59 C sigmaw : updated wake fractional area.
60 C d_deltat_gw : delta T tendency due to GW
61 c
62 C Variables d'entree :
63 c
64 c aire : aire de la maille
65 c te0 : temperature dans l'environnement (K)
66 C qe0 : humidite dans l'environnement (kg/kg)
67 C omgb : vitesse verticale moyenne sur la maille (Pa/s)
68 C dtdwn: source de chaleur due aux descentes (K/s)
69 C dqdwn: source d'humidite due aux descentes (kg/kg/s)
70 C dta : source de chaleur due courants satures et detrain (K/s)
71 C dqa : source d'humidite due aux courants satures et detra (kg/kg/s)
72 C amdwn: flux de masse total des descentes, par unite de
73 C surface de la maille (kg/m2/s)
74 C amup : flux de masse total des ascendances, par unite de
75 C surface de la maille (kg/m2/s)
76 C p : pressions aux milieux des couches (Pa)
77 C ph : pressions aux interfaces (Pa)
78 C pi : (p/p_0)**kapa (adim)
79 C dtime: increment temporel (s)
80 c
81 C Variables internes :
82 c
83 c rhow : masse volumique de la poche froide
84 C rho : environment density at P levels
85 C rhoh : environment density at Ph levels
86 C te : environment temperature | may change within
87 C qe : environment humidity | sub-time-stepping
88 C the : environment potential temperature
89 C thu : potential temperature in undisturbed area
90 C tu : temperature in undisturbed area
91 C qu : humidity in undisturbed area
92 C dp_omgb: vertical gradient og LS omega
93 C omgbw : wake average vertical omega
94 C dp_omgbw: vertical gradient of omgbw
95 C omgbdq : flux of Delta_q transported by LS omega
96 C dth : potential temperature diff. wake-undist.
97 C th1 : first pot. temp. for vertical advection (=thu)
98 C th2 : second pot. temp. for vertical advection (=thw)
99 C q1 : first humidity for vertical advection
100 C q2 : second humidity for vertical advection
101 C d_deltatw : terme de redistribution pour deltatw
102 C d_deltaqw : terme de redistribution pour deltaqw
103 C deltatw0 : deltatw initial
104 C deltaqw0 : deltaqw initial
105 C hw0 : hw initial
106 C sigmaw0: sigmaw initial
107 C amflux : horizontal mass flux through wake boundary
108 C wdens_ref: initial number of wakes per unit area (3D) or per
109 C unit length (2D), at the beginning of each time step
110 C Tgw : 1 sur la période de onde de gravité
111 c Cgw : vitesse de propagation de onde de gravité
112 c LL : distance entre 2 poches
113 
114 c-------------------------------------------------------------------------
115 c Déclaration de variables
116 c-------------------------------------------------------------------------
117 
118 #include "dimensions.h"
119 #include "YOMCST.h"
120 #include "cvthermo.h"
121 #include "iniprint.h"
122 
123 c Arguments en entree
124 c--------------------
125 
126  REAL, dimension(klon,klev) :: p, pi
127  REAL, dimension(klon,klev+1) :: ph, omgb
128  REAL dtime
129  REAL, dimension(klon,klev) :: te0,qe0
130  REAL, dimension(klon,klev) :: dtdwn, dqdwn
131  REAL, dimension(klon,klev) :: wdtpbl,wdqpbl
132  REAL, dimension(klon,klev) :: udtpbl,udqpbl
133  REAL, dimension(klon,klev) :: amdwn, amup
134  REAL, dimension(klon,klev) :: dta, dqa
135  REAL, dimension(klon) :: sigd_con
136 
137 c Sorties
138 c--------
139 
140  REAL, dimension(klon,klev) :: deltatw, deltaqw, dth
141  REAL, dimension(klon,klev) :: tu, qu
142  REAL, dimension(klon,klev) :: dtls, dqls
143  REAL, dimension(klon,klev) :: dtke, dqke
144  REAL, dimension(klon,klev) :: dtpbl, dqpbl
145  REAL, dimension(klon,klev) :: spread
146  REAL, dimension(klon,klev) :: d_deltatgw
147  REAL, dimension(klon,klev) :: d_deltatw2, d_deltaqw2
148  REAL, dimension(klon,klev+1) :: omgbdth, omg
149  REAL, dimension(klon,klev) :: dp_omgb, dp_deltomg
150  REAL, dimension(klon,klev) :: d_deltat_gw
151  REAL, dimension(klon) :: hw, sigmaw, wape, fip, gfl, cstar
152  REAL, dimension(klon) :: wdens
153  INTEGER, dimension(klon) :: ktopw
154 
155 c Variables internes
156 c-------------------
157 
158 c Variables à fixer
159  REAL alon
160  REAL coefgw
161  REAL :: wdens_ref
162  REAL stark
163  REAL alpk
164  REAL delta_t_min
165  INTEGER nsub
166  REAL dtimesub
167  REAL sigmad, hwmin,wapecut
168  REAL :: sigmaw_max
169  REAL :: dens_rate
170  REAL wdens0
171 cIM 080208
172  LOGICAL, dimension(klon) :: gwake
173 
174 c Variables de sauvegarde
175  REAL, dimension(klon,klev) :: deltatw0
176  REAL, dimension(klon,klev) :: deltaqw0
177  REAL, dimension(klon,klev) :: te, qe
178  REAL, dimension(klon) :: sigmaw0, sigmaw1
179 
180 c Variables pour les GW
181  REAL, DIMENSION(klon) :: ll
182  REAL, dimension(klon,klev) :: n2
183  REAL, dimension(klon,klev) :: cgw
184  REAL, dimension(klon,klev) :: tgw
185 
186 c Variables liées au calcul de hw
187  REAL, DIMENSION(klon) :: ptop_provis, ptop, ptop_new
188  REAL, DIMENSION(klon) :: sum_dth
189  REAL, DIMENSION(klon) :: dthmin
190  REAL, DIMENSION(klon) :: z, dz, hw0
191  INTEGER, DIMENSION(klon) :: ktop, kupper
192 
193 c Sub-timestep tendencies and related variables
194  REAL d_deltatw(klon,klev),d_deltaqw(klon,klev)
195  REAL d_te(klon,klev),d_qe(klon,klev)
196  REAL d_sigmaw(klon),alpha(klon)
197  REAL q0_min(klon),q1_min(klon)
198  LOGICAL wk_adv(klon), ok_qx_qw(klon)
199  REAL epsilon
200  DATA epsilon/1.e-15/
201 
202 c Autres variables internes
203  INTEGER isubstep, k, i
204 
205  REAL, DIMENSION(klon) :: sum_thu, sum_tu, sum_qu,sum_thvu
206  REAL, DIMENSION(klon) :: sum_dq, sum_rho
207  REAL, DIMENSION(klon) :: sum_dtdwn, sum_dqdwn
208  REAL, DIMENSION(klon) :: av_thu, av_tu, av_qu, av_thvu
209  REAL, DIMENSION(klon) :: av_dth, av_dq, av_rho
210  REAL, DIMENSION(klon) :: av_dtdwn, av_dqdwn
211 
212  REAL, DIMENSION(klon,klev) :: rho, rhow
213  REAL, DIMENSION(klon,klev+1) :: rhoh
214  REAL, DIMENSION(klon,klev) :: rhow_moyen
215  REAL, DIMENSION(klon,klev) :: zh
216  REAL, DIMENSION(klon,klev+1) :: zhh
217  REAL, DIMENSION(klon,klev) :: epaisseur1, epaisseur2
218 
219  REAL, DIMENSION(klon,klev) :: the, thu
220 
221 ! REAL, DIMENSION(klon,klev) :: d_deltatw, d_deltaqw
222 
223  REAL, DIMENSION(klon,klev+1) :: omgbw
224  REAL, DIMENSION(klon) :: pupper
225  REAL, DIMENSION(klon) :: omgtop
226  REAL, DIMENSION(klon,klev) :: dp_omgbw
227  REAL, DIMENSION(klon) :: ztop, dztop
228  REAL, DIMENSION(klon,klev) :: alpha_up
229 
230  REAL, dimension(klon) :: rre1, rre2
231  REAL :: rrd1, rrd2
232  REAL, DIMENSION(klon,klev) :: th1, th2, q1, q2
233  REAL, DIMENSION(klon,klev) :: d_th1, d_th2, d_dth
234  REAL, DIMENSION(klon,klev) :: d_q1, d_q2, d_dq
235  REAL, DIMENSION(klon,klev) :: omgbdq
236 
237  REAL, dimension(klon) :: ff, gg
238  REAL, dimension(klon) :: wape2, cstar2, heff
239 
240  REAL, DIMENSION(klon,klev) :: crep
241  REAL crep_upper, crep_sol
242 
243  REAL, DIMENSION(klon,klev) :: ppi
244 
245 ccc nrlmd
246  real, dimension(klon) :: death_rate,nat_rate
247  real, dimension(klon,klev) :: entr
248  real, dimension(klon,klev) :: detr
249 
250 C-------------------------------------------------------------------------
251 c Initialisations
252 c-------------------------------------------------------------------------
253 
254 c print*, 'wake initialisations'
255 
256 c Essais d'initialisation avec sigmaw = 0.02 et hw = 10.
257 c-------------------------------------------------------------------------
258 
259  DATA wapecut,sigmad, hwmin /5.,.02,10./
260 ccc nrlmd
261  DATA sigmaw_max /0.4/
262  DATA dens_rate /0.1/
263 ccc
264 C Longueur de maille (en m)
265 c-------------------------------------------------------------------------
266 
267 c ALON = 3.e5
268  alon = 1.e6
269 
270 
271 C Configuration de coefgw,stark,wdens (22/02/06 by YU Jingmei)
272 c
273 c coefgw : Coefficient pour les ondes de gravité
274 c stark : Coefficient k dans Cstar=k*sqrt(2*WAPE)
275 c wdens : Densité de poche froide par maille
276 c-------------------------------------------------------------------------
277 
278 ccc nrlmd coefgw=10
279 c coefgw=1
280 c wdens0 = 1.0/(alon**2)
281 ccc nrlmd wdens = 1.0/(alon**2)
282 ccc nrlmd stark = 0.50
283 cCRtest
284 ccc nrlmd alpk=0.1
285 c alpk = 1.0
286 c alpk = 0.5
287 c alpk = 0.05
288 c
289  stark = 0.33
290  alpk = 0.25
291  wdens_ref = 8.e-12
292  coefgw = 4.
293  crep_upper=0.9
294  crep_sol=1.0
295 
296 ccc nrlmd Lecture du fichier wake_param.data
297  OPEN(99,file='wake_param.data',status='old',
298  $ form='formatted',err=9999)
299  READ(99,*,end=9998) stark
300  READ(99,*,end=9998) alpk
301  READ(99,*,end=9998) wdens_ref
302  READ(99,*,end=9998) coefgw
303 9998 Continue
304  CLOSE(99)
305 9999 Continue
306 c
307 c Initialisation de toutes des densites a wdens_ref.
308 c Les densites peuvent evoluer si les poches debordent
309 c (voir au tout debut de la boucle sur les substeps)
310  wdens = wdens_ref
311 c
312 c print*,'stark',stark
313 c print*,'alpk',alpk
314 c print*,'wdens',wdens
315 c print*,'coefgw',coefgw
316 ccc
317 C Minimum value for |T_wake - T_undist|. Used for wake top definition
318 c-------------------------------------------------------------------------
319 
320  delta_t_min = 0.2
321 
322 C 1. - Save initial values and initialize tendencies
323 C --------------------------------------------------
324 
325  DO k=1,klev
326  DO i=1, klon
327  ppi(i,k)=pi(i,k)
328  deltatw0(i,k) = deltatw(i,k)
329  deltaqw0(i,k)= deltaqw(i,k)
330  te(i,k) = te0(i,k)
331  qe(i,k) = qe0(i,k)
332  dtls(i,k) = 0.
333  dqls(i,k) = 0.
334  d_deltat_gw(i,k)=0.
335  d_te(i,k) = 0.
336  d_qe(i,k) = 0.
337  d_deltatw(i,k) = 0.
338  d_deltaqw(i,k) = 0.
339 !IM 060508 beg
340  d_deltatw2(i,k)=0.
341  d_deltaqw2(i,k)=0.
342 !IM 060508 end
343  ENDDO
344  ENDDO
345 c sigmaw1=sigmaw
346 c IF (sigd_con.GT.sigmaw1) THEN
347 c print*, 'sigmaw,sigd_con', sigmaw, sigd_con
348 c ENDIF
349  DO i=1, klon
350 cc sigmaw(i) = amax1(sigmaw(i),sigd_con(i))
351  sigmaw(i) = amax1(sigmaw(i),sigmad)
352  sigmaw(i) = amin1(sigmaw(i),0.99)
353  sigmaw0(i) = sigmaw(i)
354  wape(i) = 0.
355  wape2(i) = 0.
356  d_sigmaw(i) = 0.
357  ktopw(i) = 0
358  ENDDO
359 C
360 C
361 C 2. - Prognostic part
362 C --------------------
363 C
364 C
365 C 2.1 - Undisturbed area and Wake integrals
366 C ---------------------------------------------------------
367 
368  DO i=1, klon
369  z(i) = 0.
370  ktop(i)=0
371  kupper(i) = 0
372  sum_thu(i) = 0.
373  sum_tu(i) = 0.
374  sum_qu(i) = 0.
375  sum_thvu(i) = 0.
376  sum_dth(i) = 0.
377  sum_dq(i) = 0.
378  sum_rho(i) = 0.
379  sum_dtdwn(i) = 0.
380  sum_dqdwn(i) = 0.
381 
382  av_thu(i) = 0.
383  av_tu(i) =0.
384  av_qu(i) =0.
385  av_thvu(i) = 0.
386  av_dth(i) = 0.
387  av_dq(i) = 0.
388  av_rho(i) =0.
389  av_dtdwn(i) =0.
390  av_dqdwn(i) = 0.
391  ENDDO
392 c
393 c Distance between wakes
394  DO i = 1,klon
395  ll(i) = (1-sqrt(sigmaw(i)))/sqrt(wdens(i))
396  ENDDO
397 C Potential temperatures and humidity
398 c----------------------------------------------------------
399  DO k =1,klev
400  DO i=1, klon
401 ! write(*,*)'wake 1',i,k,rd,te(i,k)
402  rho(i,k) = p(i,k)/(rd*te(i,k))
403 ! write(*,*)'wake 2',rho(i,k)
404  IF(k .eq. 1) THEN
405 ! write(*,*)'wake 3',i,k,rd,te(i,k)
406  rhoh(i,k) = ph(i,k)/(rd*te(i,k))
407 ! write(*,*)'wake 4',i,k,rd,te(i,k)
408  zhh(i,k)=0
409  ELSE
410 ! write(*,*)'wake 5',rd,(te(i,k)+te(i,k-1))
411  rhoh(i,k) = ph(i,k)*2./(rd*(te(i,k)+te(i,k-1)))
412 ! write(*,*)'wake 6',(-rhoh(i,k)*RG)+zhh(i,k-1)
413  zhh(i,k)=(ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*rg)+zhh(i,k-1)
414  ENDIF
415 ! write(*,*)'wake 7',ppi(i,k)
416  the(i,k) = te(i,k)/ppi(i,k)
417  thu(i,k) = (te(i,k) - deltatw(i,k)*sigmaw(i))/ppi(i,k)
418  tu(i,k) = te(i,k) - deltatw(i,k)*sigmaw(i)
419  qu(i,k) = qe(i,k) - deltaqw(i,k)*sigmaw(i)
420 ! write(*,*)'wake 8',(rd*(te(i,k)+deltatw(i,k)))
421  rhow(i,k) = p(i,k)/(rd*(te(i,k)+deltatw(i,k)))
422  dth(i,k) = deltatw(i,k)/ppi(i,k)
423  ENDDO
424  ENDDO
425 
426  DO k = 1, klev-1
427  DO i=1, klon
428  IF(k.eq.1) THEN
429  n2(i,k)=0
430  ELSE
431  n2(i,k)=amax1(0.,-rg**2/the(i,k)*rho(i,k)*(the(i,k+1)-
432  $ the(i,k-1))/(p(i,k+1)-p(i,k-1)))
433  ENDIF
434  zh(i,k)=(zhh(i,k)+zhh(i,k+1))/2
435 
436  cgw(i,k)=sqrt(n2(i,k))*zh(i,k)
437  tgw(i,k)=coefgw*cgw(i,k)/ll(i)
438  ENDDO
439  ENDDO
440 
441  DO i=1, klon
442  n2(i,klev)=0
443  zh(i,klev)=0
444  cgw(i,klev)=0
445  tgw(i,klev)=0
446  ENDDO
447 
448 c Calcul de la masse volumique moyenne de la colonne (bdlmd)
449 c-----------------------------------------------------------------
450 
451  DO k=1,klev
452  DO i=1, klon
453  epaisseur1(i,k)=0.
454  epaisseur2(i,k)=0.
455  ENDDO
456  ENDDO
457 
458  DO i=1, klon
459  epaisseur1(i,1)= -(ph(i,2)-ph(i,1))/(rho(i,1)*rg)+1.
460  epaisseur2(i,1)= -(ph(i,2)-ph(i,1))/(rho(i,1)*rg)+1.
461  rhow_moyen(i,1) = rhow(i,1)
462  ENDDO
463 
464  DO k = 2, klev
465  DO i=1, klon
466  epaisseur1(i,k)= -(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg) +1.
467  epaisseur2(i,k)=epaisseur2(i,k-1)+epaisseur1(i,k)
468  rhow_moyen(i,k) = (rhow_moyen(i,k-1)*epaisseur2(i,k-1)+
469  $ rhow(i,k)*epaisseur1(i,k))/epaisseur2(i,k)
470  ENDDO
471  ENDDO
472 
473 C
474 C Choose an integration bound well above wake top
475 c-----------------------------------------------------------------
476 c
477 C Pupper = 50000. ! melting level
478 c Pupper = 60000.
479 c Pupper = 80000. ! essais pour case_e
480  DO i = 1,klon
481  pupper(i) = 0.6*ph(i,1)
482  pupper(i) = max(pupper(i), 45000.)
483 ccc Pupper(i) = 60000.
484  ENDDO
485 
486 C
487 C Determine Wake top pressure (Ptop) from buoyancy integral
488 C --------------------------------------------------------
489 c
490 c-1/ Pressure of the level where dth becomes less than delta_t_min.
491 
492  DO i=1,klon
493  ptop_provis(i)=ph(i,1)
494  ENDDO
495  DO k= 2,klev
496  DO i=1,klon
497 c
498 cIM v3JYG; ptop_provis(i).LT. ph(i,1)
499 c
500  IF (dth(i,k) .GT. -delta_t_min .and.
501  $ dth(i,k-1).LT. -delta_t_min .and.
502  $ ptop_provis(i).EQ. ph(i,1)) THEN
503  ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)
504  $ - (dth(i,k-1)+delta_t_min)*p(i,k)) /
505  $ (dth(i,k) - dth(i,k-1))
506  ENDIF
507  ENDDO
508  ENDDO
509 
510 c-2/ dth integral
511 
512  DO i=1,klon
513  sum_dth(i) = 0.
514  dthmin(i) = -delta_t_min
515  z(i) = 0.
516  ENDDO
517 
518  DO k = 1,klev
519  DO i=1,klon
520  dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*rg)
521  IF (dz(i) .gt. 0) THEN
522  z(i) = z(i)+dz(i)
523  sum_dth(i) = sum_dth(i) + dth(i,k)*dz(i)
524  dthmin(i) = amin1(dthmin(i),dth(i,k))
525  ENDIF
526  ENDDO
527  ENDDO
528 
529 c-3/ height of triangle with area= sum_dth and base = dthmin
530 
531  DO i=1,klon
532  hw0(i) = 2.*sum_dth(i)/amin1(dthmin(i),-0.5)
533  hw0(i) = amax1(hwmin,hw0(i))
534  ENDDO
535 
536 c-4/ now, get Ptop
537 
538  DO i=1,klon
539  z(i) = 0.
540  ptop(i) = ph(i,1)
541  ENDDO
542 
543  DO k = 1,klev
544  DO i=1,klon
545  dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg),hw0(i)-z(i))
546  IF (dz(i) .gt. 0) THEN
547  z(i) = z(i)+dz(i)
548  ptop(i) = ph(i,k)-rho(i,k)*rg*dz(i)
549  ENDIF
550  ENDDO
551  ENDDO
552 
553 
554 C-5/ Determination de ktop et kupper
555 
556  DO k=klev,1,-1
557  DO i=1,klon
558  IF (ph(i,k+1) .lt. ptop(i)) ktop(i)=k
559  IF (ph(i,k+1) .lt. pupper(i)) kupper(i)=k
560  ENDDO
561  ENDDO
562 
563 c On evite kupper = 1 et kupper = klev
564  DO i=1,klon
565  kupper(i) = max(kupper(i),2)
566  kupper(i) = min(kupper(i),klev-1)
567  ENDDO
568 
569 
570 c-6/ Correct ktop and ptop
571 
572  DO i = 1,klon
573  ptop_new(i)=ptop(i)
574  ENDDO
575  DO k= klev,2,-1
576  DO i=1,klon
577  IF (k .LE. ktop(i) .and.
578  $ ptop_new(i) .EQ. ptop(i) .and.
579  $ dth(i,k) .GT. -delta_t_min .and.
580  $ dth(i,k-1).LT. -delta_t_min) THEN
581  ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)
582  $ - (dth(i,k-1)+delta_t_min)*p(i,k)) /
583  $ (dth(i,k) - dth(i,k-1))
584  ENDIF
585  ENDDO
586  ENDDO
587 
588  DO i=1,klon
589  ptop(i) = ptop_new(i)
590  ENDDO
591 
592  DO k=klev,1,-1
593  DO i=1,klon
594  IF (ph(i,k+1) .lt. ptop(i)) ktop(i)=k
595  ENDDO
596  ENDDO
597 c
598 c-5/ Set deltatw & deltaqw to 0 above kupper
599 c
600  DO k = 1,klev
601  DO i=1,klon
602  IF (k.GE. kupper(i)) THEN
603  deltatw(i,k) = 0.
604  deltaqw(i,k) = 0.
605  ENDIF
606  ENDDO
607  ENDDO
608 c
609 C
610 C Vertical gradient of LS omega
611 C
612  DO k = 1,klev
613  DO i=1,klon
614  IF (k.LE. kupper(i)) THEN
615  dp_omgb(i,k) = (omgb(i,k+1) - omgb(i,k))/(ph(i,k+1)-ph(i,k))
616  ENDIF
617  ENDDO
618  ENDDO
619 C
620 C Integrals (and wake top level number)
621 C --------------------------------------
622 C
623 C Initialize sum_thvu to 1st level virt. pot. temp.
624 
625  DO i=1,klon
626  z(i) = 1.
627  dz(i) = 1.
628  sum_thvu(i) = thu(i,1)*(1.+eps*qu(i,1))*dz(i)
629  sum_dth(i) = 0.
630  ENDDO
631 
632  DO k = 1,klev
633  DO i=1,klon
634  dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
635  IF (dz(i) .GT. 0) THEN
636  z(i) = z(i)+dz(i)
637  sum_thu(i) = sum_thu(i) + thu(i,k)*dz(i)
638  sum_tu(i) = sum_tu(i) + tu(i,k)*dz(i)
639  sum_qu(i) = sum_qu(i) + qu(i,k)*dz(i)
640  sum_thvu(i) = sum_thvu(i) + thu(i,k)*(1.+eps*qu(i,k))*dz(i)
641  sum_dth(i) = sum_dth(i) + dth(i,k)*dz(i)
642  sum_dq(i) = sum_dq(i) + deltaqw(i,k)*dz(i)
643  sum_rho(i) = sum_rho(i) + rhow(i,k)*dz(i)
644  sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i,k)*dz(i)
645  sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i,k)*dz(i)
646  ENDIF
647  ENDDO
648  ENDDO
649 c
650  DO i=1,klon
651  hw0(i) = z(i)
652  ENDDO
653 c
654 C
655 C 2.1 - WAPE and mean forcing computation
656 C ---------------------------------------
657 C
658 C ---------------------------------------
659 C
660 C Means
661 
662  DO i=1,klon
663  av_thu(i) = sum_thu(i)/hw0(i)
664  av_tu(i) = sum_tu(i)/hw0(i)
665  av_qu(i) = sum_qu(i)/hw0(i)
666  av_thvu(i) = sum_thvu(i)/hw0(i)
667 c av_thve = sum_thve/hw0
668  av_dth(i) = sum_dth(i)/hw0(i)
669  av_dq(i) = sum_dq(i)/hw0(i)
670  av_rho(i) = sum_rho(i)/hw0(i)
671  av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
672  av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
673 
674  wape(i) = - rg*hw0(i)*(av_dth(i)
675  $ + eps*(av_thu(i)*av_dq(i)+av_dth(i)*av_qu(i)+av_dth(i)*
676  $ av_dq(i) ))/av_thvu(i)
677  ENDDO
678 C
679 C 2.2 Prognostic variable update
680 C ------------------------------
681 C
682 C Filter out bad wakes
683 
684  DO k = 1,klev
685  DO i=1,klon
686  IF ( wape(i) .LT. 0.) THEN
687  deltatw(i,k) = 0.
688  deltaqw(i,k) = 0.
689  dth(i,k) = 0.
690  ENDIF
691  ENDDO
692  ENDDO
693 c
694  DO i=1,klon
695  IF ( wape(i) .LT. 0.) THEN
696  wape(i) = 0.
697  cstar(i) = 0.
698  hw(i) = hwmin
699  sigmaw(i) = amax1(sigmad,sigd_con(i))
700  fip(i) = 0.
701  gwake(i) = .false.
702  ELSE
703  cstar(i) = stark*sqrt(2.*wape(i))
704  gwake(i) = .true.
705  ENDIF
706  ENDDO
707 
708 c
709 c Check qx and qw positivity
710 c --------------------------
711  DO i = 1,klon
712  q0_min(i)=min( (qe(i,1)-sigmaw(i)*deltaqw(i,1)),
713  $ (qe(i,1)+(1.-sigmaw(i))*deltaqw(i,1)) )
714  ENDDO
715  DO k = 2,klev
716  DO i = 1,klon
717  q1_min(i)=min( (qe(i,k)-sigmaw(i)*deltaqw(i,k)),
718  $ (qe(i,k)+(1.-sigmaw(i))*deltaqw(i,k)) )
719  IF (q1_min(i).le.q0_min(i)) THEN
720  q0_min(i)=q1_min(i)
721  ENDIF
722  ENDDO
723  ENDDO
724 c
725  DO i = 1,klon
726  ok_qx_qw(i) = q0_min(i) .GE. 0.
727  alpha(i) = 1.
728  ENDDO
729 c
730 CC -----------------------------------------------------------------
731 C Sub-time-stepping
732 C -----------------
733 C
734  nsub=10
735  dtimesub=dtime/nsub
736 c
737 c------------------------------------------------------------
738  DO isubstep = 1,nsub
739 c------------------------------------------------------------
740 c
741 c wk_adv is the logical flag enabling wake evolution in the time advance loop
742  DO i = 1,klon
743  wk_adv(i) = ok_qx_qw(i) .AND. alpha(i) .GE. 1.
744  ENDDO
745 c
746 ccc nrlmd Ajout d'un recalcul de wdens dans le cas d'un entrainement négatif de ktop à kupper --------
747 ccc On calcule pour cela une densité wdens0 pour laquelle on aurait un entrainement nul ---
748  DO i=1,klon
749 cc print *,' isubstep,wk_adv(i),cstar(i),wape(i) ',
750 cc $ isubstep,wk_adv(i),cstar(i),wape(i)
751  IF (wk_adv(i) .AND. cstar(i).GT.0.01) THEN
752  omg(i,kupper(i)+1) = - rg*amdwn(i,kupper(i)+1)/sigmaw(i)
753  $ + rg*amup(i,kupper(i)+1)/(1.-sigmaw(i))
754  wdens0 = ( sigmaw(i) / (4.*3.14) ) *
755  $ ( (1.-sigmaw(i)) * omg(i,kupper(i)+1) /
756  $ ( (ph(i,1)-pupper(i)) * cstar(i) ) ) **(2)
757  IF ( wdens(i) .LE. wdens0*1.1 ) THEN
758  wdens(i) = wdens0
759  ENDIF
760 cc print*,'omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i)
761 cc $ ,ph(i,1)-pupper(i)',
762 cc $ omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i)
763 cc $ ,ph(i,1)-pupper(i)
764  ENDIF
765  ENDDO
766 
767 ccc nrlmd
768 
769  DO i=1,klon
770  IF (wk_adv(i)) THEN
771  gfl(i) = 2.*sqrt(3.14*wdens(i)*sigmaw(i))
772  sigmaw(i)=amin1(sigmaw(i),sigmaw_max)
773  ENDIF
774  ENDDO
775  DO i=1,klon
776  IF (wk_adv(i)) THEN
777 ccc nrlmd Introduction du taux de mortalité des poches et test sur sigmaw_max=0.4
778 ccc d_sigmaw(i) = gfl(i)*Cstar(i)*dtimesub
779  IF (sigmaw(i).ge.sigmaw_max) THEN
780  death_rate(i)=gfl(i)*cstar(i)/sigmaw(i)
781  ELSE
782  death_rate(i)=0.
783  END IF
784  d_sigmaw(i) = gfl(i)*cstar(i)*dtimesub
785  $ - death_rate(i)*sigmaw(i)*dtimesub
786 c $ - nat_rate(i)*sigmaw(i)*dtimesub
787 cc print*, 'd_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i),
788 cc $ death_rate(i),ktop(i),kupper(i)',
789 cc $ d_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i),
790 cc $ death_rate(i),ktop(i),kupper(i)
791 
792 c sigmaw(i) =sigmaw(i) + gfl(i)*Cstar(i)*dtimesub
793 c sigmaw(i) =min(sigmaw(i),0.99) !!!!!!!!
794 c wdens = wdens0/(10.*sigmaw)
795 c sigmaw =max(sigmaw,sigd_con)
796 c sigmaw =max(sigmaw,sigmad)
797  ENDIF
798  ENDDO
799 C
800 C
801 c calcul de la difference de vitesse verticale poche - zone non perturbee
802 cIM 060208 differences par rapport au code initial; init. a 0 dp_deltomg
803 cIM 060208 et omg sur les niveaux de 1 a klev+1, alors que avant l'on definit
804 cIM 060208 au niveau k=1..?
805  DO k= 1,klev
806  DO i = 1,klon
807  if (wk_adv(i)) THEN !!! nrlmd
808  dp_deltomg(i,k)=0.
809  end if
810  ENDDO
811  ENDDO
812  DO k= 1,klev+1
813  DO i = 1,klon
814  if (wk_adv(i)) THEN !!! nrlmd
815  omg(i,k)=0.
816  end if
817  ENDDO
818  ENDDO
819 c
820  DO i=1,klon
821  IF (wk_adv(i)) THEN
822  z(i)= 0.
823  omg(i,1) = 0.
824  dp_deltomg(i,1) = -(gfl(i)*cstar(i))/(sigmaw(i) * (1-sigmaw(i)))
825  ENDIF
826  ENDDO
827 c
828  DO k= 2,klev
829  DO i = 1,klon
830  IF( wk_adv(i) .AND. k .LE. ktop(i)) THEN
831  dz(i) = -(ph(i,k)-ph(i,k-1))/(rho(i,k-1)*rg)
832  z(i) = z(i)+dz(i)
833  dp_deltomg(i,k)= dp_deltomg(i,1)
834  omg(i,k)= dp_deltomg(i,1)*z(i)
835  ENDIF
836  ENDDO
837  ENDDO
838 c
839  DO i = 1,klon
840  IF (wk_adv(i)) THEN
841  dztop(i)=-(ptop(i)-ph(i,ktop(i)))/(rho(i,ktop(i))*rg)
842  ztop(i) = z(i)+dztop(i)
843  omgtop(i)=dp_deltomg(i,1)*ztop(i)
844  ENDIF
845  ENDDO
846 c
847 c -----------------
848 c From m/s to Pa/s
849 c -----------------
850 c
851  DO i=1,klon
852  IF (wk_adv(i)) THEN
853  omgtop(i) = -rho(i,ktop(i))*rg*omgtop(i)
854  dp_deltomg(i,1) = omgtop(i)/(ptop(i)-ph(i,1))
855  ENDIF
856  ENDDO
857 c
858  DO k= 1,klev
859  DO i = 1,klon
860  IF( wk_adv(i) .AND. k .LE. ktop(i)) THEN
861  omg(i,k) = - rho(i,k)*rg*omg(i,k)
862  dp_deltomg(i,k) = dp_deltomg(i,1)
863  ENDIF
864  ENDDO
865  ENDDO
866 c
867 c raccordement lineaire de omg de ptop a pupper
868 
869  DO i=1,klon
870  IF ( wk_adv(i) .AND. kupper(i) .GT. ktop(i)) THEN
871  omg(i,kupper(i)+1) = - rg*amdwn(i,kupper(i)+1)/sigmaw(i)
872  $ + rg*amup(i,kupper(i)+1)/(1.-sigmaw(i))
873  dp_deltomg(i,kupper(i)) = (omgtop(i)-omg(i,kupper(i)+1))/
874  $ (ptop(i)-pupper(i))
875  ENDIF
876  ENDDO
877 c
878 cc DO i=1,klon
879 cc print*,'Pente entre 0 et kupper (référence)'
880 cc $ ,omg(i,kupper(i)+1)/(pupper(i)-ph(i,1))
881 cc print*,'Pente entre ktop et kupper'
882 cc $ ,(omg(i,kupper(i)+1)-omgtop(i))/(pupper(i)-ptop(i))
883 cc ENDDO
884 cc
885  DO k= 1,klev
886  DO i = 1,klon
887  IF( wk_adv(i) .AND. k .GT. ktop(i) .AND. k .LE. kupper(i)) THEN
888  dp_deltomg(i,k) = dp_deltomg(i,kupper(i))
889  omg(i,k) = omgtop(i)+(ph(i,k)-ptop(i))*dp_deltomg(i,kupper(i))
890  ENDIF
891  ENDDO
892  ENDDO
893 ccc nrlmd
894 cc DO i=1,klon
895 cc print*,'deltaw_ktop,deltaw_conv',omgtop(i),omg(i,kupper(i)+1)
896 cc END DO
897 ccc
898 c
899 c
900 c-- Compute wake average vertical velocity omgbw
901 c
902 c
903  DO k = 1,klev+1
904  DO i=1,klon
905  IF ( wk_adv(i)) THEN
906  omgbw(i,k) = omgb(i,k)+(1.-sigmaw(i))*omg(i,k)
907  ENDIF
908  ENDDO
909  ENDDO
910 c-- and its vertical gradient dp_omgbw
911 c
912  DO k = 1,klev
913  DO i=1,klon
914  IF ( wk_adv(i)) THEN
915  dp_omgbw(i,k) = (omgbw(i,k+1)-omgbw(i,k))/(ph(i,k+1)-ph(i,k))
916  ENDIF
917  ENDDO
918  ENDDO
919 C
920 c-- Upstream coefficients for omgb velocity
921 c-- (alpha_up(k) is the coefficient of the value at level k)
922 c-- (1-alpha_up(k) is the coefficient of the value at level k-1)
923  DO k = 1,klev
924  DO i=1,klon
925  IF ( wk_adv(i)) THEN
926  alpha_up(i,k) = 0.
927  IF (omgb(i,k) .GT. 0.) alpha_up(i,k) = 1.
928  ENDIF
929  ENDDO
930  ENDDO
931 
932 c Matrix expressing [The,deltatw] from [Th1,Th2]
933 
934  DO i=1,klon
935  IF ( wk_adv(i)) THEN
936  rre1(i) = 1.-sigmaw(i)
937  rre2(i) = sigmaw(i)
938  ENDIF
939  ENDDO
940  rrd1 = -1.
941  rrd2 = 1.
942 c
943 c-- Get [Th1,Th2], dth and [q1,q2]
944 c
945  DO k= 1,klev
946  DO i = 1,klon
947  IF( wk_adv(i) .AND. k .LE. kupper(i)+1) THEN
948  dth(i,k) = deltatw(i,k)/ppi(i,k)
949  th1(i,k) = the(i,k) - sigmaw(i) *dth(i,k) ! undisturbed area
950  th2(i,k) = the(i,k) + (1.-sigmaw(i))*dth(i,k) ! wake
951  q1(i,k) = qe(i,k) - sigmaw(i) *deltaqw(i,k) ! undisturbed area
952  q2(i,k) = qe(i,k) + (1.-sigmaw(i))*deltaqw(i,k) ! wake
953  ENDIF
954  ENDDO
955  ENDDO
956 
957  DO i=1,klon
958  if (wk_adv(i)) then !!! nrlmd
959  d_th1(i,1) = 0.
960  d_th2(i,1) = 0.
961  d_dth(i,1) = 0.
962  d_q1(i,1) = 0.
963  d_q2(i,1) = 0.
964  d_dq(i,1) = 0.
965  end if
966  ENDDO
967 
968  DO k= 2,klev
969  DO i = 1,klon
970  IF( wk_adv(i) .AND. k .LE. kupper(i)+1) THEN
971  d_th1(i,k) = th1(i,k-1)-th1(i,k)
972  d_th2(i,k) = th2(i,k-1)-th2(i,k)
973  d_dth(i,k) = dth(i,k-1)-dth(i,k)
974  d_q1(i,k) = q1(i,k-1)-q1(i,k)
975  d_q2(i,k) = q2(i,k-1)-q2(i,k)
976  d_dq(i,k) = deltaqw(i,k-1)-deltaqw(i,k)
977  ENDIF
978  ENDDO
979  ENDDO
980 
981  DO i=1,klon
982  IF( wk_adv(i)) THEN
983  omgbdth(i,1) = 0.
984  omgbdq(i,1) = 0.
985  ENDIF
986  ENDDO
987 
988  DO k= 2,klev
989  DO i = 1,klon
990  IF( wk_adv(i) .AND. k .LE. kupper(i)+1) THEN ! loop on interfaces
991  omgbdth(i,k) = omgb(i,k)*( dth(i,k-1) - dth(i,k))
992  omgbdq(i,k) = omgb(i,k)*(deltaqw(i,k-1) - deltaqw(i,k))
993  ENDIF
994  ENDDO
995  ENDDO
996 c
997 c-----------------------------------------------------------------
998  DO k= 1,klev
999  DO i = 1,klon
1000  IF( wk_adv(i) .AND. k .LE. kupper(i)-1) THEN
1001 c-----------------------------------------------------------------
1002 c
1003 c Compute redistribution (advective) term
1004 c
1005  d_deltatw(i,k) =
1006  $ dtimesub/(ph(i,k)-ph(i,k+1))*(
1007  $ rrd1*omg(i,k )*sigmaw(i) *d_th1(i,k)
1008  $ -rrd2*omg(i,k+1)*(1.-sigmaw(i))*d_th2(i,k+1)
1009  $ -(1.-alpha_up(i,k))*omgbdth(i,k) - alpha_up(i,k+1)*
1010  $ omgbdth(i,k+1))*ppi(i,k)
1011 c print*,'d_deltatw=',d_deltatw(i,k)
1012 c
1013  d_deltaqw(i,k) =
1014  $ dtimesub/(ph(i,k)-ph(i,k+1))*(
1015  $ rrd1*omg(i,k )*sigmaw(i) *d_q1(i,k)
1016  $ -rrd2*omg(i,k+1)*(1.-sigmaw(i))*d_q2(i,k+1)
1017  $ -(1.-alpha_up(i,k))*omgbdq(i,k) - alpha_up(i,k+1)*
1018  $ omgbdq(i,k+1))
1019 c print*,'d_deltaqw=',d_deltaqw(i,k)
1020 c
1021 c and increment large scale tendencies
1022 c
1023 
1024 c
1025 C
1026 CC -----------------------------------------------------------------
1027  d_te(i,k) = dtimesub*(
1028  $ ( rre1(i)*omg(i,k )*sigmaw(i) *d_th1(i,k)
1029  $ -rre2(i)*omg(i,k+1)*(1.-sigmaw(i))*d_th2(i,k+1) )
1030  $ /(ph(i,k)-ph(i,k+1))
1031 ccc nrlmd $ -sigmaw(i)*(1.-sigmaw(i))*dth(i,k)*dp_deltomg(i,k)
1032  $ -sigmaw(i)*(1.-sigmaw(i))*dth(i,k)
1033  $ *(omg(i,k)-omg(i,k+1))/(ph(i,k)-ph(i,k+1))
1034 ccc
1035  $ )*ppi(i,k)
1036 c
1037  d_qe(i,k) = dtimesub*(
1038  $ ( rre1(i)*omg(i,k )*sigmaw(i) *d_q1(i,k)
1039  $ -rre2(i)*omg(i,k+1)*(1.-sigmaw(i))*d_q2(i,k+1) )
1040  $ /(ph(i,k)-ph(i,k+1))
1041 ccc nrlmd $ -sigmaw(i)*(1.-sigmaw(i))*deltaqw(i,k)*dp_deltomg(i,k)
1042  $ -sigmaw(i)*(1.-sigmaw(i))*deltaqw(i,k)
1043  $ *(omg(i,k)-omg(i,k+1))/(ph(i,k)-ph(i,k+1))
1044 ccc
1045  $ )
1046 ccc nrlmd
1047  ELSE IF(wk_adv(i) .AND. k .EQ. kupper(i)) THEN
1048  d_te(i,k) = dtimesub*(
1049  $ ( rre1(i)*omg(i,k )*sigmaw(i) *d_th1(i,k)
1050  $ /(ph(i,k)-ph(i,k+1)))
1051  $ )*ppi(i,k)
1052 
1053  d_qe(i,k) = dtimesub*(
1054  $ ( rre1(i)*omg(i,k )*sigmaw(i) *d_q1(i,k)
1055  $ /(ph(i,k)-ph(i,k+1)))
1056  $ )
1057 
1058  ENDIF
1059 ccc
1060  ENDDO
1061  ENDDO
1062 c------------------------------------------------------------------
1063 C
1064 C Increment state variables
1065 
1066  DO k= 1,klev
1067  DO i = 1,klon
1068 ccc nrlmd IF( wk_adv(i) .AND. k .LE. kupper(i)-1) THEN
1069  IF( wk_adv(i) .AND. k .LE. kupper(i)) THEN
1070 ccc
1071 
1072 
1073 c
1074 c Coefficient de répartition
1075 
1076  crep(i,k)=crep_sol*(ph(i,kupper(i))-ph(i,k))/(ph(i,kupper(i))
1077  $ -ph(i,1))
1078  crep(i,k)=crep(i,k)+crep_upper*(ph(i,1)-ph(i,k))/(p(i,1)-
1079  $ ph(i,kupper(i)))
1080 
1081 
1082 c Reintroduce compensating subsidence term.
1083 
1084 c dtKE(k)=(dtdwn(k)*Crep(k))/sigmaw
1085 c dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k))
1086 c . /(1-sigmaw)
1087 c dqKE(k)=(dqdwn(k)*Crep(k))/sigmaw
1088 c dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k))
1089 c . /(1-sigmaw)
1090 c
1091 c dtKE(k)=(dtdwn(k)*Crep(k)+(1-Crep(k))*dta(k))/sigmaw
1092 c dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k)*Crep(k))
1093 c . /(1-sigmaw)
1094 c dqKE(k)=(dqdwn(k)*Crep(k)+(1-Crep(k))*dqa(k))/sigmaw
1095 c dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k)*Crep(k))
1096 c . /(1-sigmaw)
1097 
1098  dtke(i,k)=(dtdwn(i,k)/sigmaw(i) - dta(i,k)/(1.-sigmaw(i)))
1099  dqke(i,k)=(dqdwn(i,k)/sigmaw(i) - dqa(i,k)/(1.-sigmaw(i)))
1100 c print*,'dtKE= ',dtKE(i,k),' dqKE= ',dqKE(i,k)
1101 c
1102  dtpbl(i,k)=(wdtpbl(i,k)/sigmaw(i) - udtpbl(i,k)/(1.-sigmaw(i)))
1103  dqpbl(i,k)=(wdqpbl(i,k)/sigmaw(i) - udqpbl(i,k)/(1.-sigmaw(i)))
1104 c print*,'dtPBL= ',dtPBL(i,k),' dqPBL= ',dqPBL(i,k)
1105 c
1106 ccc nrlmd Prise en compte du taux de mortalité
1107 ccc Définitions de entr, detr
1108  detr(i,k)=0.
1109 
1110  entr(i,k)=detr(i,k)+gfl(i)*cstar(i)+
1111  $ sigmaw(i)*(1.-sigmaw(i))*dp_deltomg(i,k)
1112 
1113  spread(i,k) = (entr(i,k)-detr(i,k))/sigmaw(i)
1114 ccc spread(i,k) = (1.-sigmaw(i))*dp_deltomg(i,k)+gfl(i)*Cstar(i)/
1115 ccc $ sigmaw(i)
1116 
1117 
1118 c ajout d'un effet onde de gravité -Tgw(k)*deltatw(k) 03/02/06 YU Jingmei
1119 
1120 ! write(lunout,*)'wake.F ',i,k, dtimesub,d_deltat_gw(i,k),
1121 ! & Tgw(i,k),deltatw(i,k)
1122  d_deltat_gw(i,k)=d_deltat_gw(i,k)-tgw(i,k)*deltatw(i,k)*
1123  $ dtimesub
1124 ! write(lunout,*)'wake.F ',i,k, dtimesub,d_deltatw(i,k)
1125  ff(i)=d_deltatw(i,k)/dtimesub
1126 
1127 c Sans GW
1128 c
1129 c deltatw(k)=deltatw(k)+dtimesub*(ff+dtKE(k)-spread(k)*deltatw(k))
1130 c
1131 c GW formule 1
1132 c
1133 c deltatw(k) = deltatw(k)+dtimesub*
1134 c $ (ff+dtKE(k) - spread(k)*deltatw(k)-Tgw(k)*deltatw(k))
1135 c
1136 c GW formule 2
1137 
1138  IF (dtimesub*tgw(i,k).lt.1.e-10) THEN
1139  d_deltatw(i,k) = dtimesub*
1140  $ (ff(i)+dtke(i,k)+dtpbl(i,k)
1141 ccc $ -spread(i,k)*deltatw(i,k)
1142  $ - entr(i,k)*deltatw(i,k)/sigmaw(i)
1143  $ - (death_rate(i)*sigmaw(i)+detr(i,k))*deltatw(i,k)
1144  $ / (1.-sigmaw(i))
1145 ccc
1146  $ -tgw(i,k)*deltatw(i,k))
1147  ELSE
1148  d_deltatw(i,k) = 1/tgw(i,k)*(1-exp(-dtimesub*
1149  $ tgw(i,k)))*
1150  $ (ff(i)+dtke(i,k)+dtpbl(i,k)
1151 ccc $ -spread(i,k)*deltatw(i,k)
1152  $ - entr(i,k)*deltatw(i,k)/sigmaw(i)
1153  $ - (death_rate(i)*sigmaw(i)+detr(i,k))*deltatw(i,k)
1154  $ / (1.-sigmaw(i))
1155 ccc
1156  $ -tgw(i,k)*deltatw(i,k))
1157  ENDIF
1158 
1159  dth(i,k) = deltatw(i,k)/ppi(i,k)
1160 
1161  gg(i)=d_deltaqw(i,k)/dtimesub
1162 
1163  d_deltaqw(i,k) = dtimesub*(gg(i)+ dqke(i,k)+dqpbl(i,k)
1164 ccc $ -spread(i,k)*deltaqw(i,k))
1165  $ - entr(i,k)*deltaqw(i,k)/sigmaw(i)
1166  $ - (death_rate(i)*sigmaw(i)+detr(i,k))*deltaqw(i,k)
1167  $ /(1.-sigmaw(i)))
1168 ccc
1169 
1170 ccc nrlmd
1171 ccc d_deltatw2(i,k)=d_deltatw2(i,k)+d_deltatw(i,k)
1172 ccc d_deltaqw2(i,k)=d_deltaqw2(i,k)+d_deltaqw(i,k)
1173 ccc
1174  ENDIF
1175  ENDDO
1176  ENDDO
1177 
1178 C
1179 C Scale tendencies so that water vapour remains positive in w and x.
1180 C
1181  call wake_vec_modulation(klon,klev,wk_adv,epsilon,qe,d_qe,deltaqw,
1182  $ d_deltaqw,sigmaw,d_sigmaw,alpha)
1183 c
1184 ccc nrlmd
1185 cc print*,'alpha'
1186 cc do i=1,klon
1187 cc print*,alpha(i)
1188 cc end do
1189 ccc
1190  DO k = 1,klev
1191  DO i = 1,klon
1192  IF( wk_adv(i) .AND. k .LE. kupper(i)) THEN
1193  d_te(i,k)=alpha(i)*d_te(i,k)
1194  d_qe(i,k)=alpha(i)*d_qe(i,k)
1195  d_deltatw(i,k)=alpha(i)*d_deltatw(i,k)
1196  d_deltaqw(i,k)=alpha(i)*d_deltaqw(i,k)
1197  d_deltat_gw(i,k)=alpha(i)*d_deltat_gw(i,k)
1198  ENDIF
1199  ENDDO
1200  ENDDO
1201  DO i = 1,klon
1202  IF( wk_adv(i)) THEN
1203  d_sigmaw(i)=alpha(i)*d_sigmaw(i)
1204  ENDIF
1205  ENDDO
1206 
1207 C Update large scale variables and wake variables
1208 cIM 060208 manque DO i + remplace DO k=1,kupper(i)
1209 cIM 060208 DO k = 1,kupper(i)
1210  DO k= 1,klev
1211  DO i = 1,klon
1212  IF( wk_adv(i) .AND. k .LE. kupper(i)) THEN
1213  dtls(i,k)=dtls(i,k)+d_te(i,k)
1214  dqls(i,k)=dqls(i,k)+d_qe(i,k)
1215 ccc nrlmd
1216  d_deltatw2(i,k)=d_deltatw2(i,k)+d_deltatw(i,k)
1217  d_deltaqw2(i,k)=d_deltaqw2(i,k)+d_deltaqw(i,k)
1218 ccc
1219  ENDIF
1220  ENDDO
1221  ENDDO
1222  DO k= 1,klev
1223  DO i = 1,klon
1224  IF( wk_adv(i) .AND. k .LE. kupper(i)) THEN
1225  te(i,k) = te0(i,k) + dtls(i,k)
1226  qe(i,k) = qe0(i,k) + dqls(i,k)
1227  the(i,k) = te(i,k)/ppi(i,k)
1228  deltatw(i,k) = deltatw(i,k)+d_deltatw(i,k)
1229  deltaqw(i,k) = deltaqw(i,k)+d_deltaqw(i,k)
1230  dth(i,k) = deltatw(i,k)/ppi(i,k)
1231 cc print*,'k,qx,qw',k,qe(i,k)-sigmaw(i)*deltaqw(i,k)
1232 cc $ ,qe(i,k)+(1-sigmaw(i))*deltaqw(i,k)
1233  ENDIF
1234  ENDDO
1235  ENDDO
1236  DO i = 1,klon
1237  IF( wk_adv(i)) THEN
1238  sigmaw(i) = sigmaw(i)+d_sigmaw(i)
1239  ENDIF
1240  ENDDO
1241 c
1242 C
1243 c Determine Ptop from buoyancy integral
1244 c ---------------------------------------
1245 c
1246 c- 1/ Pressure of the level where dth changes sign.
1247 c
1248  DO i=1,klon
1249  IF ( wk_adv(i)) THEN
1250  ptop_provis(i)=ph(i,1)
1251  ENDIF
1252  ENDDO
1253 c
1254  DO k= 2,klev
1255  DO i=1,klon
1256  IF ( wk_adv(i) .AND.
1257  $ ptop_provis(i) .EQ. ph(i,1) .AND.
1258  $ dth(i,k) .GT. -delta_t_min .and.
1259  $ dth(i,k-1).LT. -delta_t_min) THEN
1260  ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)
1261  $ - (dth(i,k-1)+delta_t_min)*p(i,k)) /(dth(i,k)
1262  $ - dth(i,k-1))
1263  ENDIF
1264  ENDDO
1265  ENDDO
1266 c
1267 c- 2/ dth integral
1268 c
1269  DO i=1,klon
1270  if (wk_adv(i)) then !!! nrlmd
1271  sum_dth(i) = 0.
1272  dthmin(i) = -delta_t_min
1273  z(i) = 0.
1274  end if
1275  ENDDO
1276 
1277  DO k = 1,klev
1278  DO i=1,klon
1279  IF ( wk_adv(i)) THEN
1280  dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*rg)
1281  IF (dz(i) .gt. 0) THEN
1282  z(i) = z(i)+dz(i)
1283  sum_dth(i) = sum_dth(i) + dth(i,k)*dz(i)
1284  dthmin(i) = amin1(dthmin(i),dth(i,k))
1285  ENDIF
1286  ENDIF
1287  ENDDO
1288  ENDDO
1289 c
1290 c- 3/ height of triangle with area= sum_dth and base = dthmin
1291 
1292  DO i=1,klon
1293  IF ( wk_adv(i)) THEN
1294  hw(i) = 2.*sum_dth(i)/amin1(dthmin(i),-0.5)
1295  hw(i) = amax1(hwmin,hw(i))
1296  ENDIF
1297  ENDDO
1298 c
1299 c- 4/ now, get Ptop
1300 c
1301  DO i=1,klon
1302  if (wk_adv(i)) then !!! nrlmd
1303  ktop(i) = 0
1304  z(i)=0.
1305  end if
1306  ENDDO
1307 c
1308  DO k = 1,klev
1309  DO i=1,klon
1310  IF ( wk_adv(i)) THEN
1311  dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg),hw(i)-z(i))
1312  IF (dz(i) .gt. 0) THEN
1313  z(i) = z(i)+dz(i)
1314  ptop(i) = ph(i,k)-rho(i,k)*rg*dz(i)
1315  ktop(i) = k
1316  ENDIF
1317  ENDIF
1318  ENDDO
1319  ENDDO
1320 c
1321 c 4.5/Correct ktop and ptop
1322 c
1323  DO i=1,klon
1324  IF ( wk_adv(i)) THEN
1325  ptop_new(i)=ptop(i)
1326  ENDIF
1327  ENDDO
1328 c
1329  DO k= klev,2,-1
1330  DO i=1,klon
1331 cIM v3JYG; IF (k .GE. ktop(i)
1332  IF ( wk_adv(i) .AND.
1333  $ k .LE. ktop(i) .AND.
1334  $ ptop_new(i) .EQ. ptop(i) .AND.
1335  $ dth(i,k) .GT. -delta_t_min .and.
1336  $ dth(i,k-1).LT. -delta_t_min) THEN
1337  ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)
1338  $ - (dth(i,k-1)+delta_t_min)*p(i,k)) /(dth(i,k)
1339  $ - dth(i,k-1))
1340  ENDIF
1341  ENDDO
1342  ENDDO
1343 c
1344 c
1345  DO i=1,klon
1346  IF ( wk_adv(i)) THEN
1347  ptop(i) = ptop_new(i)
1348  ENDIF
1349  ENDDO
1350 
1351  DO k=klev,1,-1
1352  DO i=1,klon
1353  if (wk_adv(i)) then !!! nrlmd
1354  IF (ph(i,k+1) .LT. ptop(i)) ktop(i)=k
1355  end if
1356  ENDDO
1357  ENDDO
1358 c
1359 c 5/ Set deltatw & deltaqw to 0 above kupper
1360 c
1361  DO k = 1,klev
1362  DO i=1,klon
1363  IF ( wk_adv(i) .AND. k .GE. kupper(i)) THEN
1364  deltatw(i,k) = 0.
1365  deltaqw(i,k) = 0.
1366  ENDIF
1367  ENDDO
1368  ENDDO
1369 c
1370 C
1371 c-------------Cstar computation---------------------------------
1372  DO i=1, klon
1373  if (wk_adv(i)) then !!! nrlmd
1374  sum_thu(i) = 0.
1375  sum_tu(i) = 0.
1376  sum_qu(i) = 0.
1377  sum_thvu(i) = 0.
1378  sum_dth(i) = 0.
1379  sum_dq(i) = 0.
1380  sum_rho(i) = 0.
1381  sum_dtdwn(i) = 0.
1382  sum_dqdwn(i) = 0.
1383 
1384  av_thu(i) = 0.
1385  av_tu(i) =0.
1386  av_qu(i) =0.
1387  av_thvu(i) = 0.
1388  av_dth(i) = 0.
1389  av_dq(i) = 0.
1390  av_rho(i) =0.
1391  av_dtdwn(i) =0.
1392  av_dqdwn(i) = 0.
1393  end if
1394  ENDDO
1395 C
1396 C Integrals (and wake top level number)
1397 C --------------------------------------
1398 C
1399 C Initialize sum_thvu to 1st level virt. pot. temp.
1400 
1401  DO i=1,klon
1402  if (wk_adv(i)) then !!! nrlmd
1403  z(i) = 1.
1404  dz(i) = 1.
1405  sum_thvu(i) = thu(i,1)*(1.+eps*qu(i,1))*dz(i)
1406  sum_dth(i) = 0.
1407  end if
1408  ENDDO
1409 
1410  DO k = 1,klev
1411  DO i=1,klon
1412  if (wk_adv(i)) then !!! nrlmd
1413  dz(i) = -(max(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
1414  IF (dz(i) .GT. 0) THEN
1415  z(i) = z(i)+dz(i)
1416  sum_thu(i) = sum_thu(i) + thu(i,k)*dz(i)
1417  sum_tu(i) = sum_tu(i) + tu(i,k)*dz(i)
1418  sum_qu(i) = sum_qu(i) + qu(i,k)*dz(i)
1419  sum_thvu(i) = sum_thvu(i) + thu(i,k)*(1.+eps*qu(i,k))*dz(i)
1420  sum_dth(i) = sum_dth(i) + dth(i,k)*dz(i)
1421  sum_dq(i) = sum_dq(i) + deltaqw(i,k)*dz(i)
1422  sum_rho(i) = sum_rho(i) + rhow(i,k)*dz(i)
1423  sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i,k)*dz(i)
1424  sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i,k)*dz(i)
1425  ENDIF
1426  end if
1427  ENDDO
1428  ENDDO
1429 c
1430  DO i=1,klon
1431  if (wk_adv(i)) then !!! nrlmd
1432  hw0(i) = z(i)
1433  end if
1434  ENDDO
1435 c
1436 C
1437 C - WAPE and mean forcing computation
1438 C ---------------------------------------
1439 C
1440 C ---------------------------------------
1441 C
1442 C Means
1443 
1444  DO i=1,klon
1445  if (wk_adv(i)) then !!! nrlmd
1446  av_thu(i) = sum_thu(i)/hw0(i)
1447  av_tu(i) = sum_tu(i)/hw0(i)
1448  av_qu(i) = sum_qu(i)/hw0(i)
1449  av_thvu(i) = sum_thvu(i)/hw0(i)
1450  av_dth(i) = sum_dth(i)/hw0(i)
1451  av_dq(i) = sum_dq(i)/hw0(i)
1452  av_rho(i) = sum_rho(i)/hw0(i)
1453  av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
1454  av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
1455 c
1456  wape(i) = - rg*hw0(i)*(av_dth(i)
1457  $ + eps*(av_thu(i)*av_dq(i)+av_dth(i)*av_qu(i)+av_dth(i)*
1458  $ av_dq(i) ))/av_thvu(i)
1459  end if
1460  ENDDO
1461 C
1462 C Filter out bad wakes
1463 
1464  DO k = 1,klev
1465  DO i=1,klon
1466  if (wk_adv(i)) then !!! nrlmd
1467  IF ( wape(i) .LT. 0.) THEN
1468  deltatw(i,k) = 0.
1469  deltaqw(i,k) = 0.
1470  dth(i,k) = 0.
1471  ENDIF
1472  end if
1473  ENDDO
1474  ENDDO
1475 c
1476  DO i=1,klon
1477  if (wk_adv(i)) then !!! nrlmd
1478  IF ( wape(i) .LT. 0.) THEN
1479  wape(i) = 0.
1480  cstar(i) = 0.
1481  hw(i) = hwmin
1482  sigmaw(i) = max(sigmad,sigd_con(i))
1483  fip(i) = 0.
1484  gwake(i) = .false.
1485  ELSE
1486  cstar(i) = stark*sqrt(2.*wape(i))
1487  gwake(i) = .true.
1488  ENDIF
1489  end if
1490  ENDDO
1491 
1492  ENDDO ! end sub-timestep loop
1493 C
1494 C -----------------------------------------------------------------
1495 c Get back to tendencies per second
1496 c
1497  DO k = 1,klev
1498  DO i=1,klon
1499 
1500 ccc nrlmd IF ( wk_adv(i) .AND. k .LE. kupper(i)) THEN
1501  IF ( ok_qx_qw(i) .AND. k .LE. kupper(i)) THEN
1502 ccc
1503  dtls(i,k) = dtls(i,k)/dtime
1504  dqls(i,k) = dqls(i,k)/dtime
1505  d_deltatw2(i,k)=d_deltatw2(i,k)/dtime
1506  d_deltaqw2(i,k)=d_deltaqw2(i,k)/dtime
1507  d_deltat_gw(i,k) = d_deltat_gw(i,k)/dtime
1508 cc print*,'k,dqls,omg,entr,detr',k,dqls(i,k),omg(i,k),entr(i,k)
1509 cc $ ,death_rate(i)*sigmaw(i)
1510  ENDIF
1511  ENDDO
1512  ENDDO
1513 
1514 c
1515 c----------------------------------------------------------
1516 c Determine wake final state; recompute wape, cstar, ktop;
1517 c filter out bad wakes.
1518 c----------------------------------------------------------
1519 c
1520 C 2.1 - Undisturbed area and Wake integrals
1521 C ---------------------------------------------------------
1522 
1523  DO i=1,klon
1524 ccc nrlmd if (wk_adv(i)) then !!! nrlmd
1525  if (ok_qx_qw(i)) then
1526 ccc
1527  z(i) = 0.
1528  sum_thu(i) = 0.
1529  sum_tu(i) = 0.
1530  sum_qu(i) = 0.
1531  sum_thvu(i) = 0.
1532  sum_dth(i) = 0.
1533  sum_dq(i) = 0.
1534  sum_rho(i) = 0.
1535  sum_dtdwn(i) = 0.
1536  sum_dqdwn(i) = 0.
1537 
1538  av_thu(i) = 0.
1539  av_tu(i) =0.
1540  av_qu(i) =0.
1541  av_thvu(i) = 0.
1542  av_dth(i) = 0.
1543  av_dq(i) = 0.
1544  av_rho(i) =0.
1545  av_dtdwn(i) =0.
1546  av_dqdwn(i) = 0.
1547  end if
1548  ENDDO
1549 C Potential temperatures and humidity
1550 c----------------------------------------------------------
1551 
1552  DO k =1,klev
1553  DO i=1,klon
1554 ccc nrlmd IF ( wk_adv(i)) THEN
1555  if (ok_qx_qw(i)) then
1556 ccc
1557  rho(i,k) = p(i,k)/(rd*te(i,k))
1558  IF(k .eq. 1) THEN
1559  rhoh(i,k) = ph(i,k)/(rd*te(i,k))
1560  zhh(i,k)=0
1561  ELSE
1562  rhoh(i,k) = ph(i,k)*2./(rd*(te(i,k)+te(i,k-1)))
1563  zhh(i,k)=(ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*rg)+zhh(i,k-1)
1564  ENDIF
1565  the(i,k) = te(i,k)/ppi(i,k)
1566  thu(i,k) = (te(i,k) - deltatw(i,k)*sigmaw(i))/ppi(i,k)
1567  tu(i,k) = te(i,k) - deltatw(i,k)*sigmaw(i)
1568  qu(i,k) = qe(i,k) - deltaqw(i,k)*sigmaw(i)
1569  rhow(i,k) = p(i,k)/(rd*(te(i,k)+deltatw(i,k)))
1570  dth(i,k) = deltatw(i,k)/ppi(i,k)
1571  ENDIF
1572  ENDDO
1573  ENDDO
1574 
1575 C Integrals (and wake top level number)
1576 C -----------------------------------------------------------
1577 
1578 C Initialize sum_thvu to 1st level virt. pot. temp.
1579 
1580  DO i=1,klon
1581 ccc nrlmd IF ( wk_adv(i)) THEN
1582  if (ok_qx_qw(i)) then
1583 ccc
1584  z(i) = 1.
1585  dz(i) = 1.
1586  sum_thvu(i) = thu(i,1)*(1.+eps*qu(i,1))*dz(i)
1587  sum_dth(i) = 0.
1588  ENDIF
1589  ENDDO
1590 
1591  DO k = 1,klev
1592  DO i=1,klon
1593 ccc nrlmd IF ( wk_adv(i)) THEN
1594  if (ok_qx_qw(i)) then
1595 ccc
1596  dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg)
1597  IF (dz(i) .GT. 0) THEN
1598  z(i) = z(i)+dz(i)
1599  sum_thu(i) = sum_thu(i) + thu(i,k)*dz(i)
1600  sum_tu(i) = sum_tu(i) + tu(i,k)*dz(i)
1601  sum_qu(i) = sum_qu(i) + qu(i,k)*dz(i)
1602  sum_thvu(i) = sum_thvu(i) + thu(i,k)*(1.+eps*qu(i,k))*dz(i)
1603  sum_dth(i) = sum_dth(i) + dth(i,k)*dz(i)
1604  sum_dq(i) = sum_dq(i) + deltaqw(i,k)*dz(i)
1605  sum_rho(i) = sum_rho(i) + rhow(i,k)*dz(i)
1606  sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i,k)*dz(i)
1607  sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i,k)*dz(i)
1608  ENDIF
1609  ENDIF
1610  ENDDO
1611  ENDDO
1612 c
1613  DO i=1,klon
1614 ccc nrlmd IF ( wk_adv(i)) THEN
1615  if (ok_qx_qw(i)) then
1616 ccc
1617  hw0(i) = z(i)
1618  ENDIF
1619  ENDDO
1620 c
1621 C - WAPE and mean forcing computation
1622 C-------------------------------------------------------------
1623 
1624 C Means
1625 
1626  DO i=1, klon
1627 ccc nrlmd IF ( wk_adv(i)) THEN
1628  if (ok_qx_qw(i)) then
1629 ccc
1630  av_thu(i) = sum_thu(i)/hw0(i)
1631  av_tu(i) = sum_tu(i)/hw0(i)
1632  av_qu(i) = sum_qu(i)/hw0(i)
1633  av_thvu(i) = sum_thvu(i)/hw0(i)
1634  av_dth(i) = sum_dth(i)/hw0(i)
1635  av_dq(i) = sum_dq(i)/hw0(i)
1636  av_rho(i) = sum_rho(i)/hw0(i)
1637  av_dtdwn(i) = sum_dtdwn(i)/hw0(i)
1638  av_dqdwn(i) = sum_dqdwn(i)/hw0(i)
1639 
1640  wape2(i) = - rg*hw0(i)*(av_dth(i)
1641  $ + eps*(av_thu(i)*av_dq(i)+av_dth(i)*av_qu(i)
1642  $ + av_dth(i)*av_dq(i) ))/av_thvu(i)
1643  ENDIF
1644  ENDDO
1645 
1646 C Prognostic variable update
1647 C ------------------------------------------------------------
1648 
1649 C Filter out bad wakes
1650 c
1651  DO k = 1,klev
1652  DO i=1,klon
1653 ccc nrlmd IF ( wk_adv(i) .AND. wape2(i) .LT. 0.) THEN
1654  if (ok_qx_qw(i) .AND. wape2(i) .LT. 0.) then
1655 ccc
1656  deltatw(i,k) = 0.
1657  deltaqw(i,k) = 0.
1658  dth(i,k) = 0.
1659  ENDIF
1660  ENDDO
1661  ENDDO
1662 c
1663 
1664  DO i=1, klon
1665 ccc nrlmd IF ( wk_adv(i)) THEN
1666  if (ok_qx_qw(i)) then
1667 ccc
1668  IF ( wape2(i) .LT. 0.) THEN
1669  wape2(i) = 0.
1670  cstar2(i) = 0.
1671  hw(i) = hwmin
1672  sigmaw(i) = amax1(sigmad,sigd_con(i))
1673  fip(i) = 0.
1674  gwake(i) = .false.
1675  ELSE
1676  if(prt_level.ge.10) print*,'wape2>0'
1677  cstar2(i) = stark*sqrt(2.*wape2(i))
1678  gwake(i) = .true.
1679  ENDIF
1680  ENDIF
1681  ENDDO
1682 c
1683  DO i=1, klon
1684 ccc nrlmd IF ( wk_adv(i)) THEN
1685  if (ok_qx_qw(i)) then
1686 ccc
1687  ktopw(i) = ktop(i)
1688  ENDIF
1689  ENDDO
1690 c
1691  DO i=1, klon
1692 ccc nrlmd IF ( wk_adv(i)) THEN
1693  if (ok_qx_qw(i)) then
1694 ccc
1695  IF (ktopw(i) .gt. 0 .and. gwake(i)) then
1696 
1697 Cjyg1 Utilisation d'un h_efficace constant ( ~ feeding layer)
1698 ccc heff = 600.
1699 C Utilisation de la hauteur hw
1700 cc heff = 0.7*hw
1701  heff(i) = hw(i)
1702 
1703  fip(i) = 0.5*rho(i,ktopw(i))*cstar2(i)**3*heff(i)*2*
1704  $ sqrt(sigmaw(i)*wdens(i)*3.14)
1705  fip(i) = alpk * fip(i)
1706 Cjyg2
1707  ELSE
1708  fip(i) = 0.
1709  ENDIF
1710  ENDIF
1711  ENDDO
1712 c
1713 C Limitation de sigmaw
1714 
1715 ccc nrlmd
1716 c DO i=1,klon
1717 c IF (OK_qx_qw(i)) THEN
1718 c IF (sigmaw(i).GE.sigmaw_max) sigmaw(i)=sigmaw_max
1719 c ENDIF
1720 c ENDDO
1721 ccc
1722  DO k = 1,klev
1723  DO i=1, klon
1724 
1725 ccc nrlmd On maintient désormais constant sigmaw en régime permanent
1726 ccc IF ((sigmaw(i).GT.sigmaw_max).or.
1727  IF ( ((wape(i).ge.wape2(i)).and.(wape2(i).le.1.0)).or.
1728  $ (ktopw(i).le.2) .OR.
1729  $ .not. ok_qx_qw(i) ) THEN
1730 ccc
1731  dtls(i,k) = 0.
1732  dqls(i,k) = 0.
1733  deltatw(i,k) = 0.
1734  deltaqw(i,k) = 0.
1735  ENDIF
1736  ENDDO
1737  ENDDO
1738 c
1739 ccc nrlmd On maintient désormais constant sigmaw en régime permanent
1740  DO i=1, klon
1741  IF ( ((wape(i).ge.wape2(i)).and.(wape2(i).le.1.0)).or.
1742  $ (ktopw(i).le.2) .OR.
1743  $ .not. ok_qx_qw(i) ) THEN
1744  wape(i) = 0.
1745  cstar(i)=0.
1746  hw(i) = hwmin
1747  sigmaw(i) = sigmad
1748  fip(i) = 0.
1749  ELSE
1750  wape(i) = wape2(i)
1751  cstar(i)=cstar2(i)
1752  ENDIF
1753 cc print*,'wape wape2 ktopw OK_qx_qw =',
1754 cc $ wape(i),wape2(i),ktopw(i),OK_qx_qw(i)
1755  ENDDO
1756 c
1757 c
1758  RETURN
1759  END
1760 
1761  SUBROUTINE wake_vec_modulation(nlon,nl,wk_adv,epsilon,qe,d_qe,
1762  $ deltaqw,d_deltaqw,sigmaw,d_sigmaw,alpha)
1763 c------------------------------------------------------
1764 cDtermination du coefficient alpha tel que les tendances
1765 c corriges alpha*d_G, pour toutes les grandeurs G, correspondent
1766 c a une humidite positive dans la zone (x) et dans la zone (w).
1767 c------------------------------------------------------
1768 c
1769 
1770 c Input
1771  REAL qe(nlon,nl),d_qe(nlon,nl)
1772  REAL deltaqw(nlon,nl),d_deltaqw(nlon,nl)
1773  REAL sigmaw(nlon),d_sigmaw(nlon)
1774  LOGICAL wk_adv(nlon)
1775  INTEGER nl,nlon
1776 c Output
1777  REAL alpha(nlon)
1778 c Internal variables
1779  REAL zeta(nlon,nl)
1780  REAL alpha1(nlon)
1781  REAL x,a,b,c,discrim
1782  REAL epsilon
1783 ! DATA epsilon/1.e-15/
1784 c
1785  DO k=1,nl
1786  DO i = 1,nlon
1787  IF (wk_adv(i)) THEN
1788  IF ((deltaqw(i,k)+d_deltaqw(i,k)).ge.0.) then
1789  zeta(i,k)=0.
1790  ELSE
1791  zeta(i,k)=1.
1792  END IF
1793  ENDIF
1794  ENDDO
1795  DO i = 1,nlon
1796  IF (wk_adv(i)) THEN
1797  x = qe(i,k)+(zeta(i,k)-sigmaw(i))*deltaqw(i,k)
1798  $ + d_qe(i,k)+(zeta(i,k)-sigmaw(i))*d_deltaqw(i,k)
1799  $ - d_sigmaw(i)*(deltaqw(i,k)+d_deltaqw(i,k))
1800  a = -d_sigmaw(i)*d_deltaqw(i,k)
1801  b = d_qe(i,k)+(zeta(i,k)-sigmaw(i))*d_deltaqw(i,k)
1802  $ - deltaqw(i,k)*d_sigmaw(i)
1803  c = qe(i,k)+(zeta(i,k)-sigmaw(i))*deltaqw(i,k)+epsilon
1804  discrim = b*b-4.*a*c
1805 c print*, 'x, a, b, c, discrim', x, a, b, c, discrim
1806  IF (a+b .GE. 0.) THEN !! Condition suffisante pour la positivité de ovap
1807  alpha1(i)=1.
1808  ELSE
1809  IF (x .GE. 0.) THEN
1810  alpha1(i)=1.
1811  ELSE
1812  IF (a .GT. 0.) THEN
1813  alpha1(i)=0.9*min( (2.*c)/(-b+sqrt(discrim)),
1814  $ (-b+sqrt(discrim))/(2.*a) )
1815  ELSE IF (a .eq. 0.) then
1816  alpha1(i)=0.9*(-c/b)
1817  ELSE
1818 c print*,'a,b,c discrim',a,b,c discrim
1819  alpha1(i)=0.9*max( (2.*c)/(-b+sqrt(discrim)),
1820  $ (-b+sqrt(discrim))/(2.*a) )
1821  ENDIF
1822  ENDIF
1823  ENDIF
1824  alpha(i) = min(alpha(i),alpha1(i))
1825  ENDIF
1826  ENDDO
1827  ENDDO
1828 !
1829  return
1830  end
1831 
1832  Subroutine wake_scal (p,ph,ppi,dtime,sigd_con
1833  : ,te0,qe0,omgb
1834  : ,dtdwn,dqdwn,amdwn,amup,dta,dqa
1835  : ,wdtpbl,wdqpbl,udtpbl,udqpbl
1836  o ,deltatw,deltaqw,dth,hw,sigmaw,wape,fip,gfl
1837  o ,dtls,dqls
1838  o ,ktopw,omgbdth,dp_omgb,wdens
1839  o ,tu,qu
1840  o ,dtke,dqke
1841  o ,dtpbl,dqpbl
1842  o ,omg,dp_deltomg,spread
1843  o ,cstar,d_deltat_gw
1844  o ,d_deltatw2,d_deltaqw2)
1845 
1846 ***************************************************************
1847 * *
1848 * WAKE *
1849 * retour a un Pupper fixe *
1850 * *
1851 * written by : GRANDPEIX Jean-Yves 09/03/2000 *
1852 * modified by : ROEHRIG Romain 01/29/2007 *
1853 ***************************************************************
1854 c
1855  USE dimphy
1856  IMPLICIT none
1857 c============================================================================
1858 C
1859 C
1860 C But : Decrire le comportement des poches froides apparaissant dans les
1861 C grands systemes convectifs, et fournir l'energie disponible pour
1862 C le declenchement de nouvelles colonnes convectives.
1863 C
1864 C Variables d'etat : deltatw : ecart de temperature wake-undisturbed area
1865 C deltaqw : ecart d'humidite wake-undisturbed area
1866 C sigmaw : fraction d'aire occupee par la poche.
1867 C
1868 C Variable de sortie :
1869 c
1870 c wape : WAke Potential Energy
1871 c fip : Front Incident Power (W/m2) - ALP
1872 c gfl : Gust Front Length per unit area (m-1)
1873 C dtls : large scale temperature tendency due to wake
1874 C dqls : large scale humidity tendency due to wake
1875 C hw : hauteur de la poche
1876 C dp_omgb : vertical gradient of large scale omega
1877 C omgbdth: flux of Delta_Theta transported by LS omega
1878 C dtKE : differential heating (wake - unpertubed)
1879 C dqKE : differential moistening (wake - unpertubed)
1880 C omg : Delta_omg =vertical velocity diff. wake-undist. (Pa/s)
1881 C dp_deltomg : vertical gradient of omg (s-1)
1882 C spread : spreading term in dt_wake and dq_wake
1883 C deltatw : updated temperature difference (T_w-T_u).
1884 C deltaqw : updated humidity difference (q_w-q_u).
1885 C sigmaw : updated wake fractional area.
1886 C d_deltat_gw : delta T tendency due to GW
1887 c
1888 C Variables d'entree :
1889 c
1890 c aire : aire de la maille
1891 c te0 : temperature dans l'environnement (K)
1892 C qe0 : humidite dans l'environnement (kg/kg)
1893 C omgb : vitesse verticale moyenne sur la maille (Pa/s)
1894 C dtdwn: source de chaleur due aux descentes (K/s)
1895 C dqdwn: source d'humidite due aux descentes (kg/kg/s)
1896 C dta : source de chaleur due courants satures et detrain (K/s)
1897 C dqa : source d'humidite due aux courants satures et detra (kg/kg/s)
1898 C amdwn: flux de masse total des descentes, par unite de
1899 C surface de la maille (kg/m2/s)
1900 C amup : flux de masse total des ascendances, par unite de
1901 C surface de la maille (kg/m2/s)
1902 C p : pressions aux milieux des couches (Pa)
1903 C ph : pressions aux interfaces (Pa)
1904 C ppi : (p/p_0)**kapa (adim)
1905 C dtime: increment temporel (s)
1906 c
1907 C Variables internes :
1908 c
1909 c rhow : masse volumique de la poche froide
1910 C rho : environment density at P levels
1911 C rhoh : environment density at Ph levels
1912 C te : environment temperature | may change within
1913 C qe : environment humidity | sub-time-stepping
1914 C the : environment potential temperature
1915 C thu : potential temperature in undisturbed area
1916 C tu : temperature in undisturbed area
1917 C qu : humidity in undisturbed area
1918 C dp_omgb: vertical gradient og LS omega
1919 C omgbw : wake average vertical omega
1920 C dp_omgbw: vertical gradient of omgbw
1921 C omgbdq : flux of Delta_q transported by LS omega
1922 C dth : potential temperature diff. wake-undist.
1923 C th1 : first pot. temp. for vertical advection (=thu)
1924 C th2 : second pot. temp. for vertical advection (=thw)
1925 C q1 : first humidity for vertical advection
1926 C q2 : second humidity for vertical advection
1927 C d_deltatw : terme de redistribution pour deltatw
1928 C d_deltaqw : terme de redistribution pour deltaqw
1929 C deltatw0 : deltatw initial
1930 C deltaqw0 : deltaqw initial
1931 C hw0 : hw initial
1932 C sigmaw0: sigmaw initial
1933 C amflux : horizontal mass flux through wake boundary
1934 C wdens : number of wakes per unit area (3D) or per
1935 C unit length (2D)
1936 C Tgw : 1 sur la période de onde de gravité
1937 c Cgw : vitesse de propagation de onde de gravité
1938 c LL : distance entre 2 poches
1939 
1940 c-------------------------------------------------------------------------
1941 c Déclaration de variables
1942 c-------------------------------------------------------------------------
1943 
1944 #include "dimensions.h"
1945 cccc#include "dimphy.h"
1946 #include "YOMCST.h"
1947 #include "cvthermo.h"
1948 #include "iniprint.h"
1949 
1950 c Arguments en entree
1951 c--------------------
1952 
1953  REAL p(klev),ph(klev+1),ppi(klev)
1954  REAL dtime
1955  REAL te0(klev),qe0(klev)
1956  REAL omgb(klev+1)
1957  REAL dtdwn(klev), dqdwn(klev)
1958  REAL wdtpbl(klev),wdqpbl(klev)
1959  REAL udtpbl(klev),udqpbl(klev)
1960  REAL amdwn(klev), amup(klev)
1961  REAL dta(klev), dqa(klev)
1962  REAL sigd_con
1963 
1964 c Sorties
1965 c--------
1966 
1967  REAL deltatw(klev), deltaqw(klev), dth(klev)
1968  REAL tu(klev), qu(klev)
1969  REAL dtls(klev), dqls(klev)
1970  REAL dtke(klev), dqke(klev)
1971  REAL dtpbl(klev), dqpbl(klev)
1972  REAL spread(klev)
1973  REAL d_deltatgw(klev)
1974  REAL d_deltatw2(klev), d_deltaqw2(klev)
1975  REAL omgbdth(klev+1), omg(klev+1)
1976  REAL dp_omgb(klev), dp_deltomg(klev)
1977  REAL d_deltat_gw(klev)
1978  REAL hw, sigmaw, wape, fip, gfl, cstar
1979  INTEGER ktopw
1980 
1981 c Variables internes
1982 c-------------------
1983 
1984 c Variables à fixer
1985  REAL alon
1986  REAL coefgw
1987  REAL wdens0, wdens
1988  REAL stark
1989  REAL alpk
1990  REAL delta_t_min
1991  REAL pupper
1992  INTEGER nsub
1993  REAL dtimesub
1994  REAL sigmad, hwmin
1995 
1996 c Variables de sauvegarde
1997  REAL deltatw0(klev)
1998  REAL deltaqw0(klev)
1999  REAL te(klev), qe(klev)
2000  REAL sigmaw0, sigmaw1
2001 
2002 c Variables pour les GW
2003  REAL ll
2004  REAL n2(klev)
2005  REAL cgw(klev)
2006  REAL tgw(klev)
2007 
2008 c Variables liées au calcul de hw
2009  REAL ptop_provis, ptop, ptop_new
2010  REAL sum_dth
2011  REAL dthmin
2012  REAL z, dz, hw0
2013  INTEGER ktop, kupper
2014 
2015 c Autres variables internes
2016  INTEGER isubstep, k
2017 
2018  REAL sum_thu, sum_tu, sum_qu,sum_thvu
2019  REAL sum_dq, sum_rho
2020  REAL sum_dtdwn, sum_dqdwn
2021  REAL av_thu, av_tu, av_qu, av_thvu
2022  REAL av_dth, av_dq, av_rho
2023  REAL av_dtdwn, av_dqdwn
2024 
2025  REAL rho(klev), rhoh(klev+1), rhow(klev)
2026  REAL rhow_moyen(klev)
2027  REAL zh(klev), zhh(klev+1)
2028  REAL epaisseur1(klev), epaisseur2(klev)
2029 
2030  REAL the(klev), thu(klev)
2031 
2032  REAL d_deltatw(klev), d_deltaqw(klev)
2033 
2034  REAL omgbw(klev+1), omgtop
2035  REAL dp_omgbw(klev)
2036  REAL ztop, dztop
2037  REAL alpha_up(klev)
2038 
2039  REAL rre1, rre2, rrd1, rrd2
2040  REAL th1(klev), th2(klev), q1(klev), q2(klev)
2041  REAL d_th1(klev), d_th2(klev), d_dth(klev)
2042  REAL d_q1(klev), d_q2(klev), d_dq(klev)
2043  REAL omgbdq(klev)
2044 
2045  REAL ff, gg
2046  REAL wape2, cstar2, heff
2047 
2048  REAL crep(klev)
2049  REAL crep_upper, crep_sol
2050 
2051 C-------------------------------------------------------------------------
2052 c Initialisations
2053 c-------------------------------------------------------------------------
2054 
2055 c print*, 'wake initialisations'
2056 
2057 c Essais d'initialisation avec sigmaw = 0.02 et hw = 10.
2058 c-------------------------------------------------------------------------
2059 
2060  DATA sigmad, hwmin /.02,10./
2061 
2062 C Longueur de maille (en m)
2063 c-------------------------------------------------------------------------
2064 
2065 c ALON = 3.e5
2066  alon = 1.e6
2067 
2068 
2069 C Configuration de coefgw,stark,wdens (22/02/06 by YU Jingmei)
2070 c
2071 c coefgw : Coefficient pour les ondes de gravité
2072 c stark : Coefficient k dans Cstar=k*sqrt(2*WAPE)
2073 c wdens : Densité de poche froide par maille
2074 c-------------------------------------------------------------------------
2075 
2076  coefgw=10
2077 c coefgw=1
2078 c wdens0 = 1.0/(alon**2)
2079  wdens = 1.0/(alon**2)
2080  stark = 0.50
2081 cCRtest
2082  alpk=0.1
2083 c alpk = 1.0
2084 c alpk = 0.5
2085 c alpk = 0.05
2086  crep_upper=0.9
2087  crep_sol=1.0
2088 
2089 
2090 C Minimum value for |T_wake - T_undist|. Used for wake top definition
2091 c-------------------------------------------------------------------------
2092 
2093  delta_t_min = 0.2
2094 
2095 
2096 C 1. - Save initial values and initialize tendencies
2097 C --------------------------------------------------
2098 
2099  DO k=1,klev
2100  deltatw0(k) = deltatw(k)
2101  deltaqw0(k)= deltaqw(k)
2102  te(k) = te0(k)
2103  qe(k) = qe0(k)
2104  dtls(k) = 0.
2105  dqls(k) = 0.
2106  d_deltat_gw(k)=0.
2107  d_deltatw2(k)=0.
2108  d_deltaqw2(k)=0.
2109  ENDDO
2110 c sigmaw1=sigmaw
2111 c IF (sigd_con.GT.sigmaw1) THEN
2112 c print*, 'sigmaw,sigd_con', sigmaw, sigd_con
2113 c ENDIF
2114  sigmaw = max(sigmaw,sigd_con)
2115  sigmaw = max(sigmaw,sigmad)
2116  sigmaw = min(sigmaw,0.99)
2117  sigmaw0 = sigmaw
2118 c wdens=wdens0/(10.*sigmaw)
2119 c IF (sigd_con.GT.sigmaw1) THEN
2120 c print*, 'sigmaw1,sigd1', sigmaw, sigd_con
2121 c ENDIF
2122 
2123 C 2. - Prognostic part
2124 C =========================================================
2125 
2126 c print *, 'prognostic wake computation'
2127 
2128 
2129 C 2.1 - Undisturbed area and Wake integrals
2130 C ---------------------------------------------------------
2131 
2132  z = 0.
2133  ktop=0
2134  kupper = 0
2135  sum_thu = 0.
2136  sum_tu = 0.
2137  sum_qu = 0.
2138  sum_thvu = 0.
2139  sum_dth = 0.
2140  sum_dq = 0.
2141  sum_rho = 0.
2142  sum_dtdwn = 0.
2143  sum_dqdwn = 0.
2144 
2145  av_thu = 0.
2146  av_tu =0.
2147  av_qu =0.
2148  av_thvu = 0.
2149  av_dth = 0.
2150  av_dq = 0.
2151  av_rho =0.
2152  av_dtdwn =0.
2153  av_dqdwn = 0.
2154 
2155 C Potential temperatures and humidity
2156 c----------------------------------------------------------
2157 
2158  DO k =1,klev
2159  rho(k) = p(k)/(rd*te(k))
2160  IF(k .eq. 1) THEN
2161  rhoh(k) = ph(k)/(rd*te(k))
2162  zhh(k)=0
2163  ELSE
2164  rhoh(k) = ph(k)*2./(rd*(te(k)+te(k-1)))
2165  zhh(k)=(ph(k)-ph(k-1))/(-rhoh(k)*rg)+zhh(k-1)
2166  ENDIF
2167  the(k) = te(k)/ppi(k)
2168  thu(k) = (te(k) - deltatw(k)*sigmaw)/ppi(k)
2169  tu(k) = te(k) - deltatw(k)*sigmaw
2170  qu(k) = qe(k) - deltaqw(k)*sigmaw
2171  rhow(k) = p(k)/(rd*(te(k)+deltatw(k)))
2172  dth(k) = deltatw(k)/ppi(k)
2173  ll = (1-sqrt(sigmaw))/sqrt(wdens)
2174  ENDDO
2175 
2176  DO k = 1, klev-1
2177  IF(k.eq.1) THEN
2178  n2(k)=0
2179  ELSE
2180  n2(k)=max(0.,-rg**2/the(k)*rho(k)*(the(k+1)-the(k-1))
2181  $ /(p(k+1)-p(k-1)))
2182  ENDIF
2183  zh(k)=(zhh(k)+zhh(k+1))/2
2184 
2185  cgw(k)=sqrt(n2(k))*zh(k)
2186  tgw(k)=coefgw*cgw(k)/ll
2187  ENDDO
2188 
2189  n2(klev)=0
2190  zh(klev)=0
2191  cgw(klev)=0
2192  tgw(klev)=0
2193 
2194 c Calcul de la masse volumique moyenne de la colonne
2195 c-----------------------------------------------------------------
2196 
2197  DO k=1,klev
2198  epaisseur1(k)=0.
2199  epaisseur2(k)=0.
2200  ENDDO
2201 
2202  epaisseur1(1)= -(ph(2)-ph(1))/(rho(1)*rg)+1.
2203  epaisseur2(1)= -(ph(2)-ph(1))/(rho(1)*rg)+1.
2204  rhow_moyen(1) = rhow(1)
2205 
2206  DO k = 2, klev
2207  epaisseur1(k)= -(ph(k+1)-ph(k))/(rho(k)*rg) +1.
2208  epaisseur2(k)=epaisseur2(k-1)+epaisseur1(k)
2209  rhow_moyen(k) = (rhow_moyen(k-1)*epaisseur2(k-1)+
2210  $ rhow(k)*epaisseur1(k))/epaisseur2(k)
2211  ENDDO
2212 
2213 
2214 C Choose an integration bound well above wake top
2215 c-----------------------------------------------------------------
2216 
2217 c Pupper = 50000. ! melting level
2218  pupper = 60000.
2219 c Pupper = 70000.
2220 
2221 
2222 C Determine Wake top pressure (Ptop) from buoyancy integral
2223 C-----------------------------------------------------------------
2224 
2225 c-1/ Pressure of the level where dth becomes less than delta_t_min.
2226 
2227  ptop_provis=ph(1)
2228  DO k= 2,klev
2229  IF (dth(k) .GT. -delta_t_min .and.
2230  $ dth(k-1).LT. -delta_t_min) THEN
2231  ptop_provis = ((dth(k)+delta_t_min)*p(k-1)
2232  $ - (dth(k-1)+delta_t_min)*p(k)) /(dth(k) - dth(k-1))
2233  go to 25
2234  ENDIF
2235  ENDDO
2236 25 CONTINUE
2237 
2238 c-2/ dth integral
2239 
2240  sum_dth = 0.
2241  dthmin = -delta_t_min
2242  z = 0.
2243 
2244  DO k = 1,klev
2245  dz = -(max(ph(k+1),ptop_provis)-ph(k))/(rho(k)*rg)
2246  IF (dz .le. 0) go to 40
2247  z = z+dz
2248  sum_dth = sum_dth + dth(k)*dz
2249  dthmin = min(dthmin,dth(k))
2250  ENDDO
2251 40 CONTINUE
2252 
2253 c-3/ height of triangle with area= sum_dth and base = dthmin
2254 
2255  hw0 = 2.*sum_dth/min(dthmin,-0.5)
2256  hw0 = max(hwmin,hw0)
2257 
2258 c-4/ now, get Ptop
2259 
2260  z = 0.
2261  ptop = ph(1)
2262 
2263  DO k = 1,klev
2264  dz = min(-(ph(k+1)-ph(k))/(rho(k)*rg),hw0-z)
2265  IF (dz .le. 0) go to 45
2266  z = z+dz
2267  ptop = ph(k)-rho(k)*rg*dz
2268  ENDDO
2269 45 CONTINUE
2270 
2271 
2272 C-5/ Determination de ktop et kupper
2273 
2274  DO k=klev,1,-1
2275  IF (ph(k+1) .lt. ptop) ktop=k
2276  IF (ph(k+1) .lt. pupper) kupper=k
2277  ENDDO
2278 
2279 c-6/ Correct ktop and ptop
2280 
2281  ptop_new=ptop
2282  DO k= ktop,2,-1
2283  IF (dth(k) .GT. -delta_t_min .and.
2284  $ dth(k-1).LT. -delta_t_min) THEN
2285  ptop_new = ((dth(k)+delta_t_min)*p(k-1)
2286  $ - (dth(k-1)+delta_t_min)*p(k)) /(dth(k) - dth(k-1))
2287  go to 225
2288  ENDIF
2289  ENDDO
2290 225 CONTINUE
2291 
2292  ptop = ptop_new
2293 
2294  DO k=klev,1,-1
2295  IF (ph(k+1) .lt. ptop) ktop=k
2296  ENDDO
2297 
2298 c Set deltatw & deltaqw to 0 above kupper
2299 c-----------------------------------------------------------
2300 
2301  DO k = kupper,klev
2302  deltatw(k) = 0.
2303  deltaqw(k) = 0.
2304  ENDDO
2305 
2306 
2307 C Vertical gradient of LS omega
2308 C------------------------------------------------------------
2309 
2310  DO k = 1,kupper
2311  dp_omgb(k) = (omgb(k+1) - omgb(k))/(ph(k+1)-ph(k))
2312  ENDDO
2313 
2314 
2315 C Integrals (and wake top level number)
2316 C -----------------------------------------------------------
2317 
2318 C Initialize sum_thvu to 1st level virt. pot. temp.
2319 
2320  z = 1.
2321  dz = 1.
2322  sum_thvu = thu(1)*(1.+eps*qu(1))*dz
2323  sum_dth = 0.
2324 
2325  DO k = 1,klev
2326  dz = -(max(ph(k+1),ptop)-ph(k))/(rho(k)*rg)
2327  IF (dz .LE. 0) go to 50
2328  z = z+dz
2329  sum_thu = sum_thu + thu(k)*dz
2330  sum_tu = sum_tu + tu(k)*dz
2331  sum_qu = sum_qu + qu(k)*dz
2332  sum_thvu = sum_thvu + thu(k)*(1.+eps*qu(k))*dz
2333  sum_dth = sum_dth + dth(k)*dz
2334  sum_dq = sum_dq + deltaqw(k)*dz
2335  sum_rho = sum_rho + rhow(k)*dz
2336  sum_dtdwn = sum_dtdwn + dtdwn(k)*dz
2337  sum_dqdwn = sum_dqdwn + dqdwn(k)*dz
2338  ENDDO
2339 50 CONTINUE
2340 
2341  hw0 = z
2342 
2343 C 2.1 - WAPE and mean forcing computation
2344 C-------------------------------------------------------------
2345 
2346 C Means
2347 
2348  av_thu = sum_thu/hw0
2349  av_tu = sum_tu/hw0
2350  av_qu = sum_qu/hw0
2351  av_thvu = sum_thvu/hw0
2352 c av_thve = sum_thve/hw0
2353  av_dth = sum_dth/hw0
2354  av_dq = sum_dq/hw0
2355  av_rho = sum_rho/hw0
2356  av_dtdwn = sum_dtdwn/hw0
2357  av_dqdwn = sum_dqdwn/hw0
2358 
2359  wape = - rg*hw0*(av_dth
2360  $ + eps*(av_thu*av_dq+av_dth*av_qu+av_dth*av_dq ))/av_thvu
2361 
2362 C 2.2 Prognostic variable update
2363 C ------------------------------------------------------------
2364 
2365 C Filter out bad wakes
2366 
2367  IF ( wape .LT. 0.) THEN
2368  if(prt_level.ge.10) print*,'wape<0'
2369  wape = 0.
2370  hw = hwmin
2371  sigmaw = max(sigmad,sigd_con)
2372  fip = 0.
2373  DO k = 1,klev
2374  deltatw(k) = 0.
2375  deltaqw(k) = 0.
2376  dth(k) = 0.
2377  ENDDO
2378  ELSE
2379  if(prt_level.ge.10) print*,'wape>0'
2380  cstar = stark*sqrt(2.*wape)
2381  ENDIF
2382 
2383 C------------------------------------------------------------------
2384 C Sub-time-stepping
2385 C------------------------------------------------------------------
2386 
2387 c nsub=36
2388  nsub=10
2389  dtimesub=dtime/nsub
2390 
2391 c------------------------------------------------------------
2392  DO isubstep = 1,nsub
2393 c------------------------------------------------------------
2394 
2395 c print*,'---------------','substep=',isubstep,'-------------'
2396 
2397 c Evolution of sigmaw
2398 
2399 
2400  gfl = 2.*sqrt(3.14*wdens*sigmaw)
2401 
2402  sigmaw =sigmaw + gfl*cstar*dtimesub
2403  sigmaw =min(sigmaw,0.99) !!!!!!!!
2404 c wdens = wdens0/(10.*sigmaw)
2405 c sigmaw =max(sigmaw,sigd_con)
2406 c sigmaw =max(sigmaw,sigmad)
2407 
2408 c calcul de la difference de vitesse verticale poche - zone non perturbee
2409 
2410  z= 0.
2411  dp_deltomg(1:klev)=0.
2412  omg(1:klev+1)=0.
2413 
2414  omg(1) = 0.
2415  dp_deltomg(1) = -(gfl*cstar)/(sigmaw * (1-sigmaw))
2416 
2417  DO k=2,ktop
2418  dz = -(ph(k)-ph(k-1))/(rho(k-1)*rg)
2419  z = z+dz
2420  dp_deltomg(k)= dp_deltomg(1)
2421  omg(k)= dp_deltomg(1)*z
2422  ENDDO
2423 
2424  dztop=-(ptop-ph(ktop))/(rho(ktop)*rg)
2425  ztop = z+dztop
2426  omgtop=dp_deltomg(1)*ztop
2427 
2428 
2429 c Conversion de la vitesse verticale de m/s a Pa/s
2430 
2431  omgtop = -rho(ktop)*rg*omgtop
2432  dp_deltomg(1) = omgtop/(ptop-ph(1))
2433 
2434  DO k = 1,ktop
2435  omg(k) = - rho(k)*rg*omg(k)
2436  dp_deltomg(k) = dp_deltomg(1)
2437  ENDDO
2438 
2439 c raccordement lineaire de omg de ptop a pupper
2440 
2441  IF (kupper .GT. ktop) THEN
2442  omg(kupper+1) = - rg*amdwn(kupper+1)/sigmaw
2443  $ + rg*amup(kupper+1)/(1.-sigmaw)
2444  dp_deltomg(kupper) = (omgtop-omg(kupper+1))/(ptop-pupper)
2445  DO k=ktop+1,kupper
2446  dp_deltomg(k) = dp_deltomg(kupper)
2447  omg(k) = omgtop+(ph(k)-ptop)*dp_deltomg(kupper)
2448  ENDDO
2449  ENDIF
2450 
2451 c Compute wake average vertical velocity omgbw
2452 
2453  DO k = 1,klev+1
2454  omgbw(k) = omgb(k)+(1.-sigmaw)*omg(k)
2455  ENDDO
2456 
2457 c and its vertical gradient dp_omgbw
2458 
2459  DO k = 1,klev
2460  dp_omgbw(k) = (omgbw(k+1)-omgbw(k))/(ph(k+1)-ph(k))
2461  ENDDO
2462 
2463 
2464 c Upstream coefficients for omgb velocity
2465 c-- (alpha_up(k) is the coefficient of the value at level k)
2466 c-- (1-alpha_up(k) is the coefficient of the value at level k-1)
2467 
2468  DO k = 1,klev
2469  alpha_up(k) = 0.
2470  IF (omgb(k) .GT. 0.) alpha_up(k) = 1.
2471  ENDDO
2472 
2473 c Matrix expressing [The,deltatw] from [Th1,Th2]
2474 
2475  rre1 = 1.-sigmaw
2476  rre2 = sigmaw
2477  rrd1 = -1.
2478  rrd2 = 1.
2479 
2480 c Get [Th1,Th2], dth and [q1,q2]
2481 
2482  DO k = 1,kupper+1
2483  dth(k) = deltatw(k)/ppi(k)
2484  th1(k) = the(k) - sigmaw *dth(k) ! undisturbed area
2485  th2(k) = the(k) + (1.-sigmaw)*dth(k) ! wake
2486  q1(k) = qe(k) - sigmaw *deltaqw(k) ! undisturbed area
2487  q2(k) = qe(k) + (1.-sigmaw)*deltaqw(k) ! wake
2488  ENDDO
2489 
2490  d_th1(1) = 0.
2491  d_th2(1) = 0.
2492  d_dth(1) = 0.
2493  d_q1(1) = 0.
2494  d_q2(1) = 0.
2495  d_dq(1) = 0.
2496 
2497  DO k = 2,kupper+1 ! loop on interfaces
2498  d_th1(k) = th1(k-1)-th1(k)
2499  d_th2(k) = th2(k-1)-th2(k)
2500  d_dth(k) = dth(k-1)-dth(k)
2501  d_q1(k) = q1(k-1)-q1(k)
2502  d_q2(k) = q2(k-1)-q2(k)
2503  d_dq(k) = deltaqw(k-1)-deltaqw(k)
2504  ENDDO
2505 
2506  omgbdth(1) = 0.
2507  omgbdq(1) = 0.
2508 
2509  DO k = 2,kupper+1 ! loop on interfaces
2510  omgbdth(k) = omgb(k)*( dth(k-1) - dth(k))
2511  omgbdq(k) = omgb(k)*(deltaqw(k-1) - deltaqw(k))
2512  ENDDO
2513 
2514 
2515 c-----------------------------------------------------------------
2516  DO k=1,kupper-1
2517 c-----------------------------------------------------------------
2518 c
2519 c Compute redistribution (advective) term
2520 c
2521  d_deltatw(k) =
2522  $ dtimesub/(ph(k)-ph(k+1))*(
2523  $ rrd1*omg(k )*sigmaw *d_th1(k)
2524  $ -rrd2*omg(k+1)*(1.-sigmaw)*d_th2(k+1)
2525  $ -(1.-alpha_up(k))*omgbdth(k) - alpha_up(k+1)*omgbdth(k+1)
2526  $ )*ppi(k)
2527 c print*,'d_deltatw=',d_deltatw(k)
2528 c
2529  d_deltaqw(k) =
2530  $ dtimesub/(ph(k)-ph(k+1))*(
2531  $ rrd1*omg(k )*sigmaw *d_q1(k)
2532  $ -rrd2*omg(k+1)*(1.-sigmaw)*d_q2(k+1)
2533  $ -(1.-alpha_up(k))*omgbdq(k) - alpha_up(k+1)*omgbdq(k+1)
2534  $ )
2535 c print*,'d_deltaqw=',d_deltaqw(k)
2536 c
2537 c and increment large scale tendencies
2538 c
2539  dtls(k) = dtls(k) +
2540  $ dtimesub*(
2541  $ ( rre1*omg(k )*sigmaw *d_th1(k)
2542  $ -rre2*omg(k+1)*(1.-sigmaw)*d_th2(k+1) )
2543  $ /(ph(k)-ph(k+1))
2544  $ -sigmaw*(1.-sigmaw)*dth(k)*dp_deltomg(k)
2545  $ )*ppi(k)
2546 c print*,'dtls=',dtls(k)
2547 c
2548  dqls(k) = dqls(k) +
2549  $ dtimesub*(
2550  $ ( rre1*omg(k )*sigmaw *d_q1(k)
2551  $ -rre2*omg(k+1)*(1.-sigmaw)*d_q2(k+1) )
2552  $ /(ph(k)-ph(k+1))
2553  $ -sigmaw*(1.-sigmaw)*deltaqw(k)*dp_deltomg(k)
2554  $ )
2555 c print*,'dqls=',dqls(k)
2556 
2557 c-------------------------------------------------------------------
2558  ENDDO
2559 c------------------------------------------------------------------
2560 
2561 C Increment state variables
2562 
2563  DO k = 1,kupper-1
2564 
2565 c Coefficient de répartition
2566 
2567  crep(k)=crep_sol*(ph(kupper)-ph(k))/(ph(kupper)-ph(1))
2568  crep(k)=crep(k)+crep_upper*(ph(1)-ph(k))/(p(1)-ph(kupper))
2569 
2570 
2571 c Reintroduce compensating subsidence term.
2572 
2573 c dtKE(k)=(dtdwn(k)*Crep(k))/sigmaw
2574 c dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k))
2575 c . /(1-sigmaw)
2576 c dqKE(k)=(dqdwn(k)*Crep(k))/sigmaw
2577 c dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k))
2578 c . /(1-sigmaw)
2579 c
2580 c dtKE(k)=(dtdwn(k)*Crep(k)+(1-Crep(k))*dta(k))/sigmaw
2581 c dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k)*Crep(k))
2582 c . /(1-sigmaw)
2583 c dqKE(k)=(dqdwn(k)*Crep(k)+(1-Crep(k))*dqa(k))/sigmaw
2584 c dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k)*Crep(k))
2585 c . /(1-sigmaw)
2586 
2587  dtke(k)=(dtdwn(k)/sigmaw - dta(k)/(1.-sigmaw))
2588  dqke(k)=(dqdwn(k)/sigmaw - dqa(k)/(1.-sigmaw))
2589 c print*,'dtKE=',dtKE(k)
2590 c print*,'dqKE=',dqKE(k)
2591 c
2592  dtpbl(k)=(wdtpbl(k)/sigmaw - udtpbl(k)/(1.-sigmaw))
2593  dqpbl(k)=(wdqpbl(k)/sigmaw - udqpbl(k)/(1.-sigmaw))
2594 c
2595  spread(k) = (1.-sigmaw)*dp_deltomg(k)+gfl*cstar/sigmaw
2596 c print*,'spread=',spread(k)
2597 
2598 
2599 c ajout d'un effet onde de gravité -Tgw(k)*deltatw(k) 03/02/06 YU Jingmei
2600 
2601  d_deltat_gw(k)=d_deltat_gw(k)-tgw(k)*deltatw(k)* dtimesub
2602 c print*,'d_delta_gw=',d_deltat_gw(k)
2603  ff=d_deltatw(k)/dtimesub
2604 
2605 c Sans GW
2606 c
2607 c deltatw(k)=deltatw(k)+dtimesub*(ff+dtKE(k)-spread(k)*deltatw(k))
2608 c
2609 c GW formule 1
2610 c
2611 c deltatw(k) = deltatw(k)+dtimesub*
2612 c $ (ff+dtKE(k) - spread(k)*deltatw(k)-Tgw(k)*deltatw(k))
2613 c
2614 c GW formule 2
2615 
2616  IF (dtimesub*tgw(k).lt.1.e-10) THEN
2617  deltatw(k) = deltatw(k)+dtimesub*
2618  $ (ff+dtke(k)+dtpbl(k)
2619  $ - spread(k)*deltatw(k)-tgw(k)*deltatw(k))
2620  ELSE
2621  deltatw(k) = deltatw(k)+1/tgw(k)*(1-exp(-dtimesub*tgw(k)))*
2622  $ (ff+dtke(k)+dtpbl(k)
2623  $ - spread(k)*deltatw(k)-tgw(k)*deltatw(k))
2624  ENDIF
2625 
2626  dth(k) = deltatw(k)/ppi(k)
2627 
2628  gg=d_deltaqw(k)/dtimesub
2629 
2630  deltaqw(k) = deltaqw(k) +
2631  $ dtimesub*(gg+ dqke(k)+dqpbl(k) - spread(k)*deltaqw(k))
2632 
2633  d_deltatw2(k)=d_deltatw2(k)+d_deltatw(k)
2634  d_deltaqw2(k)=d_deltaqw2(k)+d_deltaqw(k)
2635  ENDDO
2636 
2637 C And update large scale variables
2638 
2639  DO k = 1,kupper
2640  te(k) = te0(k) + dtls(k)
2641  qe(k) = qe0(k) + dqls(k)
2642  the(k) = te(k)/ppi(k)
2643  ENDDO
2644 
2645 c Determine Ptop from buoyancy integral
2646 c----------------------------------------------------------------------
2647 
2648 c-1/ Pressure of the level where dth changes sign.
2649 
2650  ptop_provis=ph(1)
2651 
2652  DO k= 2,klev
2653  IF (dth(k) .GT. -delta_t_min .and.
2654  $ dth(k-1).LT. -delta_t_min) THEN
2655  ptop_provis = ((dth(k)+delta_t_min)*p(k-1)
2656  $ - (dth(k-1)+delta_t_min)*p(k)) /(dth(k) - dth(k-1))
2657  go to 65
2658  ENDIF
2659  ENDDO
2660 65 CONTINUE
2661 
2662 c-2/ dth integral
2663 
2664  sum_dth = 0.
2665  dthmin = -delta_t_min
2666  z = 0.
2667 
2668  DO k = 1,klev
2669  dz = -(max(ph(k+1),ptop_provis)-ph(k))/(rho(k)*rg)
2670  IF (dz .le. 0) go to 70
2671  z = z+dz
2672  sum_dth = sum_dth + dth(k)*dz
2673  dthmin = min(dthmin,dth(k))
2674  ENDDO
2675 70 CONTINUE
2676 
2677 c-3/ height of triangle with area= sum_dth and base = dthmin
2678 
2679  hw = 2.*sum_dth/min(dthmin,-0.5)
2680  hw = max(hwmin,hw)
2681 
2682 c-4/ now, get Ptop
2683 
2684  ktop = 0
2685  z=0.
2686 
2687  DO k = 1,klev
2688  dz = min(-(ph(k+1)-ph(k))/(rho(k)*rg),hw-z)
2689  IF (dz .le. 0) go to 75
2690  z = z+dz
2691  ptop = ph(k)-rho(k)*rg*dz
2692  ktop = k
2693  ENDDO
2694 75 CONTINUE
2695 
2696 c-5/Correct ktop and ptop
2697 
2698  ptop_new=ptop
2699 
2700  DO k= ktop,2,-1
2701  IF (dth(k) .GT. -delta_t_min .and.
2702  $ dth(k-1).LT. -delta_t_min) THEN
2703  ptop_new = ((dth(k)+delta_t_min)*p(k-1)
2704  $ - (dth(k-1)+delta_t_min)*p(k)) /(dth(k) - dth(k-1))
2705  go to 275
2706  ENDIF
2707  ENDDO
2708 275 CONTINUE
2709 
2710  ptop = ptop_new
2711 
2712  DO k=klev,1,-1
2713  IF (ph(k+1) .LT. ptop) ktop=k
2714  ENDDO
2715 
2716 c-6/ Set deltatw & deltaqw to 0 above kupper
2717 
2718  DO k = kupper,klev
2719  deltatw(k) = 0.
2720  deltaqw(k) = 0.
2721  ENDDO
2722 
2723 c------------------------------------------------------------------
2724  ENDDO ! end sub-timestep loop
2725 C -----------------------------------------------------------------
2726 
2727 c Get back to tendencies per second
2728 
2729  DO k = 1,kupper-1
2730  dtls(k) = dtls(k)/dtime
2731  dqls(k) = dqls(k)/dtime
2732  d_deltatw2(k)=d_deltatw2(k)/dtime
2733  d_deltaqw2(k)=d_deltaqw2(k)/dtime
2734  d_deltat_gw(k) = d_deltat_gw(k)/dtime
2735  ENDDO
2736 
2737 C 2.1 - Undisturbed area and Wake integrals
2738 C ---------------------------------------------------------
2739 
2740  z = 0.
2741  sum_thu = 0.
2742  sum_tu = 0.
2743  sum_qu = 0.
2744  sum_thvu = 0.
2745  sum_dth = 0.
2746  sum_dq = 0.
2747  sum_rho = 0.
2748  sum_dtdwn = 0.
2749  sum_dqdwn = 0.
2750 
2751  av_thu = 0.
2752  av_tu =0.
2753  av_qu =0.
2754  av_thvu = 0.
2755  av_dth = 0.
2756  av_dq = 0.
2757  av_rho =0.
2758  av_dtdwn =0.
2759  av_dqdwn = 0.
2760 
2761 C Potential temperatures and humidity
2762 c----------------------------------------------------------
2763 
2764  DO k =1,klev
2765  rho(k) = p(k)/(rd*te(k))
2766  IF(k .eq. 1) THEN
2767  rhoh(k) = ph(k)/(rd*te(k))
2768  zhh(k)=0
2769  ELSE
2770  rhoh(k) = ph(k)*2./(rd*(te(k)+te(k-1)))
2771  zhh(k)=(ph(k)-ph(k-1))/(-rhoh(k)*rg)+zhh(k-1)
2772  ENDIF
2773  the(k) = te(k)/ppi(k)
2774  thu(k) = (te(k) - deltatw(k)*sigmaw)/ppi(k)
2775  tu(k) = te(k) - deltatw(k)*sigmaw
2776  qu(k) = qe(k) - deltaqw(k)*sigmaw
2777  rhow(k) = p(k)/(rd*(te(k)+deltatw(k)))
2778  dth(k) = deltatw(k)/ppi(k)
2779 
2780  ENDDO
2781 
2782 C Integrals (and wake top level number)
2783 C -----------------------------------------------------------
2784 
2785 C Initialize sum_thvu to 1st level virt. pot. temp.
2786 
2787  z = 1.
2788  dz = 1.
2789  sum_thvu = thu(1)*(1.+eps*qu(1))*dz
2790  sum_dth = 0.
2791 
2792  DO k = 1,klev
2793  dz = -(max(ph(k+1),ptop)-ph(k))/(rho(k)*rg)
2794 
2795  IF (dz .LE. 0) go to 51
2796  z = z+dz
2797  sum_thu = sum_thu + thu(k)*dz
2798  sum_tu = sum_tu + tu(k)*dz
2799  sum_qu = sum_qu + qu(k)*dz
2800  sum_thvu = sum_thvu + thu(k)*(1.+eps*qu(k))*dz
2801  sum_dth = sum_dth + dth(k)*dz
2802  sum_dq = sum_dq + deltaqw(k)*dz
2803  sum_rho = sum_rho + rhow(k)*dz
2804  sum_dtdwn = sum_dtdwn + dtdwn(k)*dz
2805  sum_dqdwn = sum_dqdwn + dqdwn(k)*dz
2806  ENDDO
2807  51 CONTINUE
2808 
2809  hw0 = z
2810 
2811 C 2.1 - WAPE and mean forcing computation
2812 C-------------------------------------------------------------
2813 
2814 C Means
2815 
2816  av_thu = sum_thu/hw0
2817  av_tu = sum_tu/hw0
2818  av_qu = sum_qu/hw0
2819  av_thvu = sum_thvu/hw0
2820  av_dth = sum_dth/hw0
2821  av_dq = sum_dq/hw0
2822  av_rho = sum_rho/hw0
2823  av_dtdwn = sum_dtdwn/hw0
2824  av_dqdwn = sum_dqdwn/hw0
2825 
2826  wape2 = - rg*hw0*(av_dth
2827  $ + eps*(av_thu*av_dq+av_dth*av_qu+av_dth*av_dq ))/av_thvu
2828 
2829 
2830 C 2.2 Prognostic variable update
2831 C ------------------------------------------------------------
2832 
2833 C Filter out bad wakes
2834 
2835  IF ( wape2 .LT. 0.) THEN
2836  if(prt_level.ge.10) print*,'wape2<0'
2837  wape2 = 0.
2838  hw = hwmin
2839  sigmaw = max(sigmad,sigd_con)
2840  fip = 0.
2841  DO k = 1,klev
2842  deltatw(k) = 0.
2843  deltaqw(k) = 0.
2844  dth(k) = 0.
2845  ENDDO
2846  ELSE
2847  if(prt_level.ge.10) print*,'wape2>0'
2848  cstar2 = stark*sqrt(2.*wape2)
2849 
2850  ENDIF
2851 
2852  ktopw = ktop
2853 
2854  IF (ktopw .gt. 0) then
2855 
2856 Cjyg1 Utilisation d'un h_efficace constant ( ~ feeding layer)
2857 ccc heff = 600.
2858 C Utilisation de la hauteur hw
2859 cc heff = 0.7*hw
2860  heff = hw
2861 
2862  fip = 0.5*rho(ktopw)*cstar2**3*heff*2*sqrt(sigmaw*wdens*3.14)
2863  fip = alpk * fip
2864 Cjyg2
2865  ELSE
2866  fip = 0.
2867  ENDIF
2868 
2869 
2870 C Limitation de sigmaw
2871 c
2872 C sécurité : si le wake occuppe plus de 90 % de la surface de la maille,
2873 C alors il disparait en se mélangeant à la partie undisturbed
2874 
2875 ! correction NICOLAS . ((wape.ge.wape2).and.(wape2.le.1.0))) THEN
2876  IF ((sigmaw.GT.0.9).or.
2877  . ((wape.ge.wape2).and.(wape2.le.1.0)).or.(ktopw.le.2)) THEN
2878 cIM cf NR/JYG 251108 . ((wape.ge.wape2).and.(wape2.le.1.0))) THEN
2879 c IF (sigmaw.GT.0.9) THEN
2880  DO k = 1,klev
2881  dtls(k) = 0.
2882  dqls(k) = 0.
2883  deltatw(k) = 0.
2884  deltaqw(k) = 0.
2885  ENDDO
2886  wape = 0.
2887  hw = hwmin
2888  sigmaw = sigmad
2889  fip = 0.
2890  ENDIF
2891 
2892  RETURN
2893  END
2894 
2895 
2896