Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

PHY_Atm_CP_RUN.f90 @ 2100

Last change on this file since 2100 was 2089, checked in by Laurent Fairhead, 10 years ago

Inclusion de la physique de MAR

Integration of MAR physics

File size: 18.6 KB

Line
1	subroutine PHY_Atm_CP_RUN(mzc,kcolc)
2
3	!------------------------------------------------------------------------------+
4	! Mon 17-Jun-2013 MAR \|
5	! subroutine PHY_Atm_CP_RUN interfaces \|
6	! Bechtold et al. (2001) Convective Parameterization with MAR \|
7	! \|
8	! version 3.p.4.1 created by H. Gallee, Mon 8-Apr-2013 \|
9	! Last Modification by H. Gallee, Mon 17-Jun-2013 \|
10	! \|
11	!------------------------------------------------------------------------------+
12	! \|
13	! INPUT \|
14	! ^^^^^ it_EXP : Experiment Iteration Counter \|
15	! dt__CP : Mass Flux Scheme: Time Step \|
16	! dxHOST : grid spacing of HOST MODEL [m] \|
17	! Ta__DY(kcolp,mzpp) : air temperature [K] \|
18	! qs__CM(kcolp,mzpp) : air snow Particl. concentr. [kg/kg] \|
19	! qr__CM(kcolp,mzp ) : air rain drops concentr. [kg/kg] \|
20	! \|
21	! INPUT/OUTPUT \|
22	! ^^^^^^^^^^^^ qv__DY(kcolp,mzpp) : Specific Humidity [kg/kg] \|
23	! qw__CM(kcolp,mzp ) : air cloud droplets concentr. [kg/kg] \|
24	! qi__CM(kcolp,mzp ) : air cloud crystals concentr. [kg/kg] \|
25	! snowCP(kcolp ) : Snow (convective) [m] \|
26	! snowCM(kcolp ) : Snow (convective + stratiform) [m] \|
27	! rainCP(kcolp ) : Rain (convective) [m] \|
28	! rainCM(kcolp ) : Rain (convective + stratiform) [m] \|
29	! \|
30	! OUTPUT dpktCP(kcolp,mzp ) : Reduc. Pot.Temperat.Tendency [K/X/s] \|
31	! ^^^^^^ dqv_CP(kcolp,mzp ) : Specific Humidity Tendency [kg/kg/s] \|
32	! dqw_CP(kcolp,mzp ) : cloud dropl.Concent.Tendency [kg/kg/s] \|
33	! dqi_CP(kcolp,mzp ) : cloud cryst.Concent.Tendency [kg/kg/s] \|
34	! \|
35	! dss_CP(kcolp ) : Snow (convective) Tendency [m/s] \|
36	! drr_CP(kcolp ) : Rain (convective) Tendency [m/s] \|
37	! \|
38	! pkt_DY(kcolp,mzp ) : Reduced Potential Temperature [K/X] \|
39	! \|
40	! CAPECP(kcolp ) : Convective Avail.Potent.Energy \|
41	! \|
42	! REFER. : MesoNH MASS FLUX CONVECTION Routine \|
43	! ^^^^^^^^ (Bechtold et al., 2001, QJRMS 127, pp 869-886) \|
44	! \|
45	! # OPTION : #EW Energy and Water Conservation \|
46	! # ^^^^^^^^ \|
47	! \|
48	! MODIF. HGallee: 18-11-2004: Adaptation to CVAmnh.f90.laurent \|
49	! ^^^^^^ (Argument kensbl of CONVECTION removed) \|
50	! \|
51	!------------------------------------------------------------------------------+
52
53	use Mod_Real
54	use Mod_PHY____dat
55	use Mod_PHY____grd
56	use Mod_PHY_CP_ctr
57	use Mod_PHY_CP_dat
58	use Mod_PHY_CP_grd
59	use Mod_PHY_CP_kkl
60	use Mod_PHY_DY_kkl
61	use Mod_PHY_AT_kkl
62	use Mod_PHY_CM_kkl
63
64	use Mod_Atm_CP_RUN
65
66	IMPLICIT NONE
67
68
69
70	! Global Variables
71	! ================
72
73	integer :: mzc ! Nb of levels
74	integer :: kcolc ! Nb of columns
75
76
77
78
79	! Local Variables
80	! ================
81
82	real(kind=real8) :: bANA !
83	real(kind=real8) :: zANA !
84
85	integer :: i ! x-Axis Index
86	integer :: j ! y-Axis Index
87	integer :: k ! Level Index (from top to bottom)
88	integer :: klc ! Level Index (from bottom to top )
89	integer :: ikl ! Column Index
90
91	REAL :: Pdxdy0(kcolc)
92	REAL :: P_pa_0(kcolc,mzc)
93	REAL :: P_za_0(kcolc,mzc)
94	REAL :: P_Ta_0(kcolc,mzc)
95	REAL :: P_Qa_0(kcolc,mzc)
96	REAL :: P_Qw_0(kcolc,mzc)
97	REAL :: P_Qi_0(kcolc,mzc)
98	REAL :: P_Ua_0(kcolc,mzc)
99	REAL :: P_Va_0(kcolc,mzc)
100	REAL :: P_Wa_0(kcolc,mzc)
101
102	integer :: Kstep1(kcolc) ! convective counter
103	integer :: K_CbT1(kcolc) ! cloud top level
104	integer :: K_CbB1(kcolc) ! cloud base level
105
106	integer, parameter :: KTCCH0=1 !
107	REAL :: P_CH_0(kcolc,mzc,KTCCH0)
108	REAL :: PdCH_1(kcolc,mzc,KTCCH0)
109
110	REAL :: PdTa_1(kcolc,mzc)
111	REAL :: PdQa_1(kcolc,mzc)
112	REAL :: PdQw_1(kcolc,mzc)
113	REAL :: PdQi_1(kcolc,mzc)
114	REAL :: Pdrr_1(kcolc)
115	REAL :: Pdss_1(kcolc)
116	REAL :: PuMF_1(kcolc,mzc) ! Upward Mass Flux
117	REAL :: PdMF_1(kcolc,mzc) ! Downward Mass Flux
118	REAL :: Pfrr_1(kcolc,mzc) ! Liquid Precipitation Flux
119	REAL :: Pfss_1(kcolc,mzc) ! Solid Precipitation Flux
120	REAL :: Pcape1(kcolc) ! CAPE [J/kg]
121
122
123	! Diagnostic Variables
124	! --------------------
125
126	! #EW integer :: irmx ,jrmx ,iter_0 !
127	! #EW real(kind=real8) :: rr_max,temp_r,energ0 !
128	! #EW real(kind=real8) :: water0,waterb !
129
130
131
132	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
133	! !
134	! ALLOCATION
135	! ----------
136
137	IF (it_RUN.EQ.1 .OR. FlagDALLOC) THEN !
138	allocate ( wa_ANA(kcolc,mzc) ) ! ANAbatic Wind Speed [m/s]
139	END IF
140	! !
141	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
142
143
144
145
146	! Update Convective Mass Flux
147	! ===========================
148
149	IF ( mod(iitCV0,jjtCV0).EQ.0 ) THEN
150
151	! Martin control
152	PRINT*,'Dans PHY_Atm_CP_RUN'
153	PRINT*,'size(ii__AP)',size(ii__AP)
154	PRINT*,'size(jj__AP)',size(jj__AP)
155	PRINT*,'size(wa_ANA)=',size(wa_ANA)
156	PRINT*,'ii__AP(1)=',ii__AP(1)
157	PRINT*,'jj__AP(1)=',jj__AP(1)
158	PRINT*,'ii__AP(kcolp)=',ii__AP(kcolp)
159	PRINT*,'jj__AP(kcolp)=',jj__AP(kcolp)
160	PRINT*,'mzp=',mzp
161	! Martin control
162
163	DO ikl = 1,kcolp
164
165	! PRINT*,'ikl=',ikl
166
167	i = ii__AP(ikl)
168	j = jj__AP(ikl)
169
170	! PRINT*,'ii__AP(',ikl,')=',ii__AP(ikl)
171	! PRINT*,'jj__AP(',ikl,')=',jj__AP(ikl)
172
173
174	! Contribution from Subgrid Mountain Breeze
175	! -----------------------------------------
176
177	DO k=1,mzp
178	wa_ANA(ikl,k ) = 0.0
179	END DO
180	IF(Lo_ANA) THEN
181	DO k=1,mzp
182	bANA = min(Z___DY(ikl,k ), zi__AT(ikl))
183	zANA = hANA(ikl) + bANA * 2.0
184	IF (Z___DY(ikl,k ) .LE. zANA .AND. &
185	& Ta__DY(ikl,mzpp) .GT. Ta__DY(ikl,mzp)) THEN
186	wa_ANA(ikl,k ) = rANA * bANA * 0.5 ! Half Integrated Horizontal Divergence
187
188	END IF
189	END DO
190	END IF
191
192
193
194	! Mass Flux convective Scheme: Set Up Vertical Profiles
195	! -----------------------------------------------------
196
197	Pdxdy0(ikl) = dxHOST * dxHOST ! grid area [m2]
198	DO klc= 1,mzp
199	k = mzpp-klc
200	P_pa_0(ikl,klc) = (psa_DY( ikl)sigma(k) + pt__DY) 1.e3 ! pressure in layer [Pa]
201	P_za_0(ikl,klc) = Z___DY(ikl,k) ! height of model layer [m]
202	P_Ta_0(ikl,klc) = Ta__DY(ikl,k) ! grid scale T at time t [K]
203	P_Qa_0(ikl,klc) = qv__DY(ikl,k) ! grid scale water vapor at time t [kg/kg]
204	P_Qw_0(ikl,klc) = qw__CM(ikl,k) / (1.0-qw__CM(ikl,k)) ! grid scale Cloud drops at time t [kg/kg]
205	P_Qi_0(ikl,klc) = qi__CM(ikl,k) / (1.0-qi__CM(ikl,k)) ! grid scale Cloud ice at time t [kg/kg]
206	P_Ua_0(ikl,klc) = ua__DY(ikl,k) ! grid scale hor. wind u at time t [m/s]
207	P_Va_0(ikl,klc) = va__DY(ikl,k) ! grid scale hor. wind v at time t [m/s]
208	P_Wa_0(ikl,klc) = wa__DY(ikl,k) + wa_ANA(ikl,k) ! grid scale vertic.wind at time t [m/s]
209	END DO
210
211	END DO ! ikl = 1,kcolp
212
213
214	! Mass Flux convective Scheme: Bechtold et al. (2001) Convective Parameterization
215	! -------------------------------------------------------------------------------
216
217	! ***************
218	call CONVECTION( &
219	& kcolc , mzc , kidia0, kfdia0, kbdia0, ktdia0, &
220	& pdtCV , Odeep , Oshal , Orset0, Odown0, kIce_0, &
221	& OsetA0, PTdcv , PTscv , &
222	& kensbl, &
223	& P_pa_0, P_za_0, Pdxdy0, &
224	& P_Ta_0, P_Qa_0, P_Qw_0, P_Qi_0, P_Ua_0, P_Va_0, P_Wa_0, &
225	& Kstep1, PdTa_1, PdQa_1, PdQw_1, PdQi_1, &
226	& Pdrr_1, Pdss_1, &
227	& PuMF_1, PdMF_1, Pfrr_1, Pfss_1, Pcape1, K_CbT1, K_CbB1, &
228	& OCvTC0, KTCCH0, P_CH_0, PdCH_1)
229	! ***************
230
231
232
233	! Mass Flux convective Scheme: products
234	! -------------------------------------
235
236	DO ikl = 1,kcolp
237	CAPECP(ikl) = Pcape1(ikl)
238	timeCP(ikl) = Kstep1(ikl) *dt__CP
239	drr_CP(ikl) = Pdrr_1(ikl)
240	dss_CP(ikl) = Pdss_1(ikl)
241
242	DO klc= 1,mzp
243	k = mzpp - klc
244	dpktCP(ikl,k) = PdTa_1(ikl,klc) /ExnrDY(ikl,k)
245	dqv_CP(ikl,k) = PdQa_1(ikl,klc)
246	dqw_CP(ikl,k) = PdQw_1(ikl,klc)
247	dqi_CP(ikl,k) = PdQi_1(ikl,klc)
248
249	END DO
250
251	END DO
252
253
254	END IF ! mod(iitCV0,jjtCV0).EQ.0 (UPDATE) THEN
255
256
257
258
259	! Vertical Integrated Energy and Water Content
260	! ============================================
261
262	! #EW irmx = i_x0
263	! #EW jrmx = j_y0
264	! #EW rr_max = 0.0
265
266	DO ikl = 1,kcolp
267	i = ii__AP(ikl)
268	j = jj__AP(ikl)
269
270	! #EW enr0EW(ikl) = 0.0
271	! #EW wat0EW(ikl) = 0.0
272
273	! #EW DO k=1,mzp
274	! #EW temp_r = pkt_DY(ikl,k)*ExnrDY(ikl,k)
275	! #EW enr0EW(ikl) = enr0EW(ikl) &
276	! #EW& +(temp_r &
277	! #EW& -(qw__CM(ikl,k)+qr__CM(ikl,k)) * Lv_CPd &
278	! #EW& -(qi__CM(ikl,k)+qs__CM(ikl,k)) * Ls_CPd)*dsigmi(k)
279	! #EW wat0EW(ikl) = wat0EW(ikl) &
280	! #EW& +(qv__DY(ikl,k) &
281	! #EW& + qw__CM(ikl,k)+qr__CM(ikl,k) &
282	! #EW& + qi__CM(ikl,k)+qs__CM(ikl,k) )*dsigmi(k)
283	! #EW END DO
284
285	! #EW enr0EW(ikl) = enr0EW(ikl) * psa_DY( ikl) * Grav_I
286	! #EW wat0EW(ikl) = wat0EW(ikl) * psa_DY( ikl) * Grav_I
287	! .. wat0EW [m] contains implicit factor 1.d3 [kPa-->Pa] /ro_Wat
288
289	! #EW energ0 = energ0 - enr0EW(ikl)
290	! #EW water0 = water0 - wat0EW(ikl)
291
292
293
294
295	! Update of Mass Flux convective Tendencies
296	! =========================================
297
298	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
299
300	IF (timeCP(ikl) .GT. 0.) THEN
301
302
303	! Temporal tendencies on pkt_DY, qv__DY and rainCM
304	! ------------------------------------------------
305
306	DO k=1,mzp
307	dqv_CP(ikl,k) = min(dqv_CP(ikl,k),(qv__DY(ikl,k)-epsq)/dt__CP)
308	dqw_CP(ikl,k) = min(dqw_CP(ikl,k),(qw__CM(ikl,k)-epsq)/dt__CP)
309	dqi_CP(ikl,k) = min(dqi_CP(ikl,k),(qi__CM(ikl,k)-epsq)/dt__CP)
310	ENDDO
311
312	rainCM(ikl) = rainCM(ikl) + drr_CP(ikl) *dt__CP
313	rainCP(ikl) = rainCP(ikl) + drr_CP(ikl) *dt__CP
314
315	snowCM(ikl) = snowCM(ikl) + dss_CP(ikl) *dt__CP
316	snowCP(ikl) = snowCP(ikl) + dss_CP(ikl) *dt__CP
317
318	ELSE
319
320	DO k=1,mzp
321	dpktCP(ikl,k) = 0.00
322	dqv_CP(ikl,k) = 0.00
323	dqw_CP(ikl,k) = 0.00
324	dqi_CP(ikl,k) = 0.00
325	ENDDO
326
327
328	ENDIF ! {timeCP(ikl) .gt. 0}
329
330	timeCP(ikl) = max(timeCP(ikl) - dt__CP,zer0)
331	! .. ^^^^ remaining time before the end of convection
332	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
333
334
335
336
337	! Vertical Integrated Energy and Water Content
338	! ============================================
339
340	! #EW enr1EW(ikl) = 0.0
341	! #EW wat1EW(ikl) = 0.0
342	! #EW watfEW(ikl) =-drr_CP(ikl)
343
344	! #EW DO k=1,mz
345	! #EW temp_r = pkt_DY(ikl,k)*ExnrDY(ikl,k)
346	! #EW enr1EW(ikl) = enr1EW(ikl) &
347	! #EW& +(temp_r &
348	! #EW& -(qw__CM(ikl,k)+qr__CM(ikl,k)) *Lv_CPd &
349	! #EW& -(qi__CM(ikl,k)+qs__CM(ikl,k)) Ls_CPd)dsigmi(k)
350	! #EW wat1EW(ikl) = wat1EW(ikl) &
351	! #EW& +(qv__DY(ikl,k) &
352	! #EW& + qw__CM(ikl,k)+qr__CM(ikl,k) &
353	! #EW& + qi__CM(ikl,k)+qs__CM(ikl,k) )*dsigmi(k)
354	! #EW END DO
355
356	! #EW enr1EW(ikl) = enr1EW(ikl) * psa_DY( ikl) *Grav_I &
357	! #EW& - drr_CP(ikl) *Lv_CPd
358	! #EW wat1EW(ikl) = wat1EW(ikl) * psa_DY( ikl) *Grav_I
359	! .. wat1EW [m] contains implicit factor 1.d3 [kPa-->Pa] /rhoWat
360
361	! #EW energ0 = energ0 + enr1EW(ikl)
362	! #EW water0 = water0 + wat1EW(ikl)
363	! #EW iter_0 = iter_0 + 1
364
365
366
367
368	! Vertical Integrated Energy and Water Content: OUTPUT
369	! ====================================================
370
371	! #EW IF (drr_CP(ikl).gt.rr_max) THEN
372	! #EW rr_max = drr_CP(ikl)
373	! #EW irmx = i
374	! #EW jrmx = j
375	! #EW END IF
376
377	END DO
378
379	! #EW waterb = wat1EW(irmx,jrmx)-wat0EW(irmx,jrmx)-watfEW(irmx,jrmx)
380	! #EW write(6,606) it_EXP,enr0EW(irmx,jrmx),1.d3*wat0EW(irmx,jrmx), &
381	! #EW& irmx,jrmx,enr1EW(irmx,jrmx),1.d3*wat1EW(irmx,jrmx), &
382	! #EW& 1.d3*watfEW(irmx,jrmx), &
383	! #EW& 1.d3*waterb , &
384	! #EW& energ0/iter_0 , water0/iter_0
385	606 format(i9,' Before CVAj: E0 =',f12.6,' W0 = ',f9.6, &
386	& /,i5,i4,' After CVAj: E1 =',f12.6,' W1 = ',f9.6, &
387	& ' W Flux =',f9.6, &
388	& ' Div(W) =',e9.3, &
389	& /,9x,' Mean dE =',f12.9,' dW = ',e9.3)
390
391
392
393
394	! Incrementation Step Nb
395	! ======================
396
397	iitCV0 = iitCV0 + 1
398
399
400	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
401	! !
402	! DE-ALLOCATION
403	! =============
404
405	IF (FlagDALLOC) THEN !
406	deallocate ( wa_ANA ) ! ANAbatic Wind Speed [m/s]
407	END IF
408
409	! Bug corrigé par Martin:
410	! IF (it_RUN.EQ.1 .OR. FlagDALLOC) THEN !
411	! deallocate ( wa_ANA ) ! ANAbatic Wind Speed [m/s]
412	! END IF
413	! !
414	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
415
416
417
418	return
419	end subroutine PHY_Atm_CP_RUN

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

Download in other formats:

Original Format