Context Navigation

wake.F @ 922

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format