Context Navigation

module_ForDiagnostics.f90 @ 1899

Last change on this file since 1899 was 1833, checked in by lfita, 8 years ago
Fixing calculation of afwa unstable indices
File size: 29.3 KB

Line
1	!! Fortran version of different diagnostics
2	! L. Fita. LMD May 2016
3	! gfortran module_generic.o module_ForDiagnosticsVars.o -c module_ForDiagnostics.F90
4	!
5	! f2py -m module_ForDiagnostics --f90exec=/usr/bin/gfortran-4.7 -c module_generic.F90 module_ForDiagnosticsVars.F90 module_ForDiagnostics.F90
6
7	MODULE module_ForDiagnostics
8
9	USE module_definitions
10	USE module_generic
11	USE module_ForDiagnosticsVars
12
13	CONTAINS
14
15	!!!!!!! Calculations
16	! compute_cape_afwa4D: Subroutine to use WRF phys/module_diag_afwa.F `buyoancy' subroutine to compute
17	! CAPE, CIN, ZLFC, PLFC, LI
18	! compute_cllmh4D3: Computation of low, medium and high cloudiness from a 4D CLDFRA and pressure being
19	! 3rd dimension the z-dim
20	! compute_cllmh3D3: Computation of low, medium and high cloudiness from a 3D CLDFRA and pressure being
21	! 3rd dimension the z-dim
22	! compute_cllmh: Computation of low, medium and high cloudiness
23	! compute_clt4D3: Computation of total cloudiness from a 4D CLDFRA being 3rd dimension the z-dim
24	! compute_clt3D3: Computation of total cloudiness from a 3D CLDFRA being 3rd dimension the z-dim
25	! compute_clt: Computation of total cloudiness
26	! compute_psl_ecmwf: Compute sea level pressure using ECMWF method following Mats Hamrud and Philippe Courtier [Pa]
27	! compute_massvertint1D: Subroutine to vertically integrate a 1D variable in eta vertical coordinates
28	! compute_vertint1D: Subroutine to vertically integrate a 1D variable in any vertical coordinates
29	! compute_zint4D: Subroutine to vertically integrate a 4D variable in any vertical coordinates
30	! compute_zmla_generic4D: Subroutine to compute pbl-height following a generic method
31	! compute_zwind4D: Subroutine to compute extrapolate the wind at a given height following the 'power law' methodology
32	! compute_zwind_log4D: Subroutine to compute extrapolate the wind at a given height following the 'logarithmic law' methodology
33	! compute_zwindMCO3D: Subroutine to compute extrapolate the wind at a given height following the 'power law' methodolog
34
35	!!!
36	! Calculations
37	!!!
38
39	SUBROUTINE compute_cllmh4D2(cldfra4D, pres4D, cllmh4D2, d1, d2, d3, d4)
40	! Subroutine to compute the low, medium and high cloudiness following 'newmicro.F90' from LMDZ from a 4D CLDFRA and pressure
41	! where zdim is the 2nd dimension (thus, cldfra4D(d1,d2,d3,d4) --> cllmh(3,d1,d3,d4) 1: low, 2: medium, 3: high
42	! It should be properly done via an 'INTERFACE', but...
43
44	IMPLICIT NONE
45
46	INTEGER, INTENT(in) :: d1, d2, d3, d4
47	REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in) :: cldfra4D, pres4D
48	REAL(r_k), DIMENSION(3,d1,d3,d4), INTENT(out) :: cllmh4D2
49
50	! Local
51	INTEGER :: i,j,k, zdim, Ndim
52
53	!!!!!!! Variables
54	! cldfra4D: 4D cloud fraction values [1]
55	! pres4D: 4D pressure values [Pa]
56	! Ndim: number of dimensions of the input data
57	! d[1-4]: dimensions of 'cldfra4D'
58	! zdim: number of the vertical-dimension within the matrix
59	! cltlmh4D2: low, medium, high cloudiness for the 4D cldfra and d2 being 'zdim'
60
61	fname = 'compute_cllmh4D2'
62	zdim = 2
63	Ndim = 4
64
65	DO i=1, d1
66	DO j=1, d3
67	DO k=1, d4
68	cllmh4D2(:,i,j,k) = var_cllmh(cldfra4D(i,:,j,k), pres4D(i,:,j,k), d2)
69	END DO
70	END DO
71	END DO
72
73	RETURN
74
75	END SUBROUTINE compute_cllmh4D2
76
77	SUBROUTINE compute_cllmh3D1(cldfra3D, pres3D, cllmh3D1, d1, d2, d3)
78	! Subroutine to compute the low, medium and high cloudiness following 'newmicro.F90' from LMDZ from a 3D CLDFRA and pressure
79	! where zdim is the 1st dimension (thus, cldfra3D(d1,d2,d3) --> cllmh(3,d2,d3) 1: low, 2: medium, 3: high
80	! It should be properly done via an 'INTERFACE', but...
81
82	IMPLICIT NONE
83
84	INTEGER, INTENT(in) :: d1, d2, d3
85	REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in) :: cldfra3D, pres3D
86	REAL(r_k), DIMENSION(3,d2,d3), INTENT(out) :: cllmh3D1
87
88	! Local
89	INTEGER :: i,j,k, zdim, Ndim
90
91	!!!!!!! Variables
92	! cldfra3D: 3D cloud fraction values [1]
93	! pres3D: 3D pressure values [Pa]
94	! Ndim: number of dimensions of the input data
95	! d[1-3]: dimensions of 'cldfra3D'
96	! zdim: number of the vertical-dimension within the matrix
97	! cltlmh3D1: low, medium, high cloudiness for the 3D cldfra and d1 being 'zdim'
98
99	fname = 'compute_cllmh3D1'
100	zdim = 1
101	Ndim = 3
102
103	DO i=1, d1
104	DO j=1, d2
105	cllmh3D1(:,i,j) = var_cllmh(cldfra3D(:,i,j), pres3D(:,i,j), d1)
106	END DO
107	END DO
108
109	RETURN
110
111	END SUBROUTINE compute_cllmh3D1
112
113	SUBROUTINE compute_cllmh(cldfra1D, cldfra2D, cldfra3D, cldfra4D, pres1D, pres2D, pres3D, pres4D, &
114	Ndim, zdim, cllmh1D, cllmh2D1, cllmh2D2, cllmh3D1, cllmh3D2, cllmh3D3, cllmh4D1, cllmh4D2, &
115	cllmh4D3, cllmh4D4, d1, d2, d3, d4)
116	! Subroutine to compute the low, medium and high cloudiness following 'newmicro.F90' from LMDZ
117
118	IMPLICIT NONE
119
120	INTEGER, INTENT(in) :: Ndim, d1, d2, d3, d4, zdim
121	REAL(r_k), DIMENSION(d1), OPTIONAL, INTENT(in) :: cldfra1D, pres1D
122	REAL(r_k), DIMENSION(d1,d2), OPTIONAL, INTENT(in) :: cldfra2D, pres2D
123	REAL(r_k), DIMENSION(d1,d2,d3), OPTIONAL, INTENT(in) :: cldfra3D, pres3D
124	REAL(r_k), DIMENSION(d1,d2,d3,d4), OPTIONAL, &
125	INTENT(in) :: cldfra4D, pres4D
126	REAL(r_k), DIMENSION(3), OPTIONAL, INTENT(out) :: cllmh1D
127	REAL(r_k), DIMENSION(d1,3), OPTIONAL, INTENT(out) :: cllmh2D1
128	REAL(r_k), DIMENSION(d2,3), OPTIONAL, INTENT(out) :: cllmh2D2
129	REAL(r_k), DIMENSION(d2,d3,3), OPTIONAL, INTENT(out) :: cllmh3D1
130	REAL(r_k), DIMENSION(d1,d3,3), OPTIONAL, INTENT(out) :: cllmh3D2
131	REAL(r_k), DIMENSION(d1,d2,3), OPTIONAL, INTENT(out) :: cllmh3D3
132	REAL(r_k), DIMENSION(d2,d3,d4,3), OPTIONAL, &
133	INTENT(out) :: cllmh4D1
134	REAL(r_k), DIMENSION(d1,d3,d4,3), OPTIONAL, &
135	INTENT(out) :: cllmh4D2
136	REAL(r_k), DIMENSION(d1,d2,d4,3), OPTIONAL, &
137	INTENT(out) :: cllmh4D3
138	REAL(r_k), DIMENSION(d1,d2,d3,3), OPTIONAL, &
139	INTENT(out) :: cllmh4D4
140
141	! Local
142	INTEGER :: i,j,k
143
144	!!!!!!! Variables
145	! cldfra[1-4]D: cloud fraction values [1]
146	! pres[1-4]D: pressure values [Pa]
147	! Ndim: number of dimensions of the input data
148	! d[1-4]: dimensions of 'cldfra'
149	! zdim: number of the vertical-dimension within the matrix
150	! cllmh1D: low, medium and high cloudiness for the 1D cldfra
151	! cllmh2D1: low, medium and high cloudiness for the 2D cldfra and d1 being 'zdim'
152	! cllmh2D2: low, medium and high cloudiness for the 2D cldfra and d2 being 'zdim'
153	! cllmh3D1: low, medium and high cloudiness for the 3D cldfra and d1 being 'zdim'
154	! cllmh3D2: low, medium and high cloudiness for the 3D cldfra and d2 being 'zdim'
155	! cllmh3D3: low, medium and high cloudiness for the 3D cldfra and d3 being 'zdim'
156	! cllmh4D1: low, medium and high cloudiness for the 4D cldfra and d1 being 'zdim'
157	! cllmh4D2: low, medium and high cloudiness for the 4D cldfra and d2 being 'zdim'
158	! cllmh4D3: low, medium and high cloudiness for the 4D cldfra and d3 being 'zdim'
159	! cllmh4D4: low, medium and high cloudiness for the 4D cldfra and d4 being 'zdim'
160
161	fname = 'compute_cllmh'
162
163	SELECT CASE (Ndim)
164	CASE (1)
165	cllmh1D = var_cllmh(cldfra1D, pres1D, d1)
166	CASE (2)
167	IF (zdim == 1) THEN
168	DO i=1, d2
169	cllmh2D1(i,:) = var_cllmh(cldfra2D(:,i), pres2D(:,i), d1)
170	END DO
171	ELSE IF (zdim == 2) THEN
172	DO i=1, d1
173	cllmh2D2(i,:) = var_cllmh(cldfra2D(:,i), pres2D(i,:), d2)
174	END DO
175	ELSE
176	PRINT *,TRIM(ErrWarnMsg('err'))
177	PRINT *,' ' // TRIM(fname) // ': wrong zdim:', zdim,' for Ndim=', Ndim, ' !!'
178	PRINT *,' accepted values: 1,2'
179	STOP
180	END IF
181	CASE (3)
182	IF (zdim == 1) THEN
183	DO i=1, d2
184	DO j=1, d3
185	cllmh3D1(i,j,:) = var_cllmh(cldfra3D(:,i,j), pres3D(:,i,j), d1)
186	END DO
187	END DO
188	ELSE IF (zdim == 2) THEN
189	DO i=1, d1
190	DO j=1, d3
191	cllmh3D2(i,j,:) = var_cllmh(cldfra3D(i,:,j), pres3D(i,:,j), d2)
192	END DO
193	END DO
194	ELSE IF (zdim == 3) THEN
195	DO i=1, d1
196	DO j=1, d2
197	cllmh3D3(i,j,:) = var_cllmh(cldfra3D(i,j,:), pres3D(i,j,:), d3)
198	END DO
199	END DO
200	ELSE
201	PRINT *,TRIM(ErrWarnMsg('err'))
202	PRINT *,' ' // TRIM(fname) // ': wrong zdim:', zdim,' for Ndim=', Ndim, ' !!'
203	PRINT *,' accepted values: 1,2,3'
204	STOP
205	END IF
206	CASE (4)
207	IF (zdim == 1) THEN
208	DO i=1, d2
209	DO j=1, d3
210	DO k=1, d4
211	cllmh4D1(i,j,k,:) = var_cllmh(cldfra4D(:,i,j,k), pres4D(:,i,j,k), d1)
212	END DO
213	END DO
214	END DO
215	ELSE IF (zdim == 2) THEN
216	DO i=1, d1
217	DO j=1, d3
218	DO k=1, d4
219	cllmh4D2(i,j,k,:) = var_cllmh(cldfra4D(i,:,j,k), pres4D(i,:,j,k), d2)
220	END DO
221	END DO
222	END DO
223	ELSE IF (zdim == 3) THEN
224	DO i=1, d2
225	DO j=1, d3
226	DO k=1, d4
227	cllmh4D3(i,j,k,:) = var_cllmh(cldfra4D(i,j,:,k), pres4D(i,j,:,k), d3)
228	END DO
229	END DO
230	END DO
231	ELSE IF (zdim == 4) THEN
232	DO i=1, d1
233	DO j=1, d2
234	DO k=1, d3
235	cllmh4D4(i,j,k,:) = var_cllmh(cldfra4D(i,j,k,:), pres4D(i,j,k,:), d4)
236	END DO
237	END DO
238	END DO
239	ELSE
240	PRINT *,TRIM(ErrWarnMsg('err'))
241	PRINT *,' ' // TRIM(fname) // ': wrong zdim:', zdim,' for Ndim=', Ndim, ' !!'
242	PRINT *,' accepted values: 1,2,3,4'
243	STOP
244	END IF
245	CASE DEFAULT
246	PRINT *,TRIM(ErrWarnMsg('err'))
247	PRINT *,' ' // TRIM(fname) // ': Ndim:', Ndim,' not ready !!'
248	STOP
249	END SELECT
250
251	RETURN
252
253	END SUBROUTINE compute_cllmh
254
255	SUBROUTINE compute_clt4D2(cldfra4D, clt4D2, d1, d2, d3, d4)
256	! Subroutine to compute the total cloudiness following 'newmicro.F90' from LMDZ from a 4D CLDFRA
257	! where zdim is the 2nd dimension (thus, cldfra4D(d1,d2,d3,d4) --> clt(d1,d3,d4)
258	! It should be properly done via an 'INTERFACE', but...
259
260	IMPLICIT NONE
261
262	INTEGER, INTENT(in) :: d1, d2, d3, d4
263	REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in) :: cldfra4D
264	REAL(r_k), DIMENSION(d1,d3,d4), INTENT(out) :: clt4D2
265
266	! Local
267	INTEGER :: i,j,k, zdim, Ndim
268
269	!!!!!!! Variables
270	! cldfra4D: 4D cloud fraction values [1]
271	! Ndim: number of dimensions of the input data
272	! d[1-4]: dimensions of 'cldfra4D'
273	! zdim: number of the vertical-dimension within the matrix
274	! clt4D2: total cloudiness for the 4D cldfra and d2 being 'zdim'
275
276	fname = 'compute_clt4D2'
277	zdim = 2
278	Ndim = 4
279
280	DO i=1, d1
281	DO j=1, d3
282	DO k=1, d4
283	clt4D2(i,j,k) = var_clt(cldfra4D(i,:,j,k), d2)
284	END DO
285	END DO
286	END DO
287
288	RETURN
289
290	END SUBROUTINE compute_clt4D2
291
292	SUBROUTINE compute_clt3D1(cldfra3D, clt3D1, d1, d2, d3)
293	! Subroutine to compute the total cloudiness following 'newmicro.F90' from LMDZ from a 3D CLDFRA
294	! where zdim is the 1st dimension (thus, cldfra4D(d1,d2,d3) --> clt(d2,d3)
295	! It should be properly done via an 'INTERFACE', but...
296
297	IMPLICIT NONE
298
299	INTEGER, INTENT(in) :: d1, d2, d3
300	REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in) :: cldfra3D
301	REAL(r_k), DIMENSION(d2,d3), INTENT(out) :: clt3D1
302
303	! Local
304	INTEGER :: i,j,k, zdim, Ndim
305
306	!!!!!!! Variables
307	! cldfra3D: 3D cloud fraction values [1]
308	! Ndim: number of dimensions of the input data
309	! d[1-3]: dimensions of 'cldfra3D'
310	! zdim: number of the vertical-dimension within the matrix
311	! clt3D1: total cloudiness for the 3D cldfra and d1 being 'zdim'
312
313	fname = 'compute_clt3D1'
314	zdim = 1
315	Ndim = 3
316
317	DO i=1, d2
318	DO j=1, d3
319	clt3D1(i,j) = var_clt(cldfra3D(:,i,j), d1)
320	END DO
321	END DO
322
323	RETURN
324
325	END SUBROUTINE compute_clt3D1
326
327	SUBROUTINE compute_clt(cldfra1D, cldfra2D, cldfra3D, cldfra4D, Ndim, zdim, clt1D, clt2D1, clt2D2, &
328	clt3D1, clt3D2, clt3D3, clt4D1, clt4D2, clt4D3, clt4D4, d1, d2, d3, d4)
329	! Subroutine to compute the total cloudiness following 'newmicro.F90' from LMDZ
330
331	IMPLICIT NONE
332
333	INTEGER, INTENT(in) :: Ndim, d1, d2, d3, d4, zdim
334	REAL(r_k), DIMENSION(d1), OPTIONAL, INTENT(in) :: cldfra1D
335	REAL(r_k), DIMENSION(d1,d2), OPTIONAL, INTENT(in) :: cldfra2D
336	REAL(r_k), DIMENSION(d1,d2,d3), OPTIONAL, INTENT(in) :: cldfra3D
337	REAL(r_k), DIMENSION(d1,d2,d3,d4), OPTIONAL, &
338	INTENT(in) :: cldfra4D
339	REAL(r_k), OPTIONAL, INTENT(out) :: clt1D
340	REAL(r_k), DIMENSION(d1), OPTIONAL, INTENT(out) :: clt2D1
341	REAL(r_k), DIMENSION(d2), OPTIONAL, INTENT(out) :: clt2D2
342	REAL(r_k), DIMENSION(d2,d3), OPTIONAL, INTENT(out) :: clt3D1
343	REAL(r_k), DIMENSION(d1,d3), OPTIONAL, INTENT(out) :: clt3D2
344	REAL(r_k), DIMENSION(d1,d2), OPTIONAL, INTENT(out) :: clt3D3
345	REAL(r_k), DIMENSION(d2,d3,d4), OPTIONAL,INTENT(out) :: clt4D1
346	REAL(r_k), DIMENSION(d1,d3,d4), OPTIONAL,INTENT(out) :: clt4D2
347	REAL(r_k), DIMENSION(d1,d2,d4), OPTIONAL,INTENT(out) :: clt4D3
348	REAL(r_k), DIMENSION(d1,d2,d3), OPTIONAL,INTENT(out) :: clt4D4
349
350	! Local
351	INTEGER :: i,j,k
352
353	!!!!!!! Variables
354	! cldfra[1-4]D: cloud fraction values [1]
355	! Ndim: number of dimensions of the input data
356	! d[1-4]: dimensions of 'cldfra'
357	! zdim: number of the vertical-dimension within the matrix
358	! clt1D: total cloudiness for the 1D cldfra
359	! clt2D1: total cloudiness for the 2D cldfra and d1 being 'zdim'
360	! clt2D2: total cloudiness for the 2D cldfra and d2 being 'zdim'
361	! clt3D1: total cloudiness for the 3D cldfra and d1 being 'zdim'
362	! clt3D2: total cloudiness for the 3D cldfra and d2 being 'zdim'
363	! clt3D3: total cloudiness for the 3D cldfra and d3 being 'zdim'
364	! clt4D1: total cloudiness for the 4D cldfra and d1 being 'zdim'
365	! clt4D2: total cloudiness for the 4D cldfra and d2 being 'zdim'
366	! clt4D3: total cloudiness for the 4D cldfra and d3 being 'zdim'
367	! clt4D4: total cloudiness for the 4D cldfra and d4 being 'zdim'
368
369	fname = 'compute_clt'
370
371	SELECT CASE (Ndim)
372	CASE (1)
373	clt1D = var_clt(cldfra1D, d1)
374	CASE (2)
375	IF (zdim == 1) THEN
376	DO i=1, d2
377	clt2D1(i) = var_clt(cldfra2D(:,i), d1)
378	END DO
379	ELSE IF (zdim == 2) THEN
380	DO i=1, d1
381	clt2D2(i) = var_clt(cldfra2D(:,i), d2)
382	END DO
383	ELSE
384	PRINT *,TRIM(ErrWarnMsg('err'))
385	PRINT *,' ' // TRIM(fname) // ': wrong zdim:', zdim,' for Ndim=', Ndim, ' !!'
386	PRINT *,' accepted values: 1,2'
387	STOP
388	END IF
389	CASE (3)
390	IF (zdim == 1) THEN
391	DO i=1, d2
392	DO j=1, d3
393	clt3D1(i,j) = var_clt(cldfra3D(:,i,j), d1)
394	END DO
395	END DO
396	ELSE IF (zdim == 2) THEN
397	DO i=1, d1
398	DO j=1, d3
399	clt3D2(i,j) = var_clt(cldfra3D(i,:,j), d2)
400	END DO
401	END DO
402	ELSE IF (zdim == 3) THEN
403	DO i=1, d1
404	DO j=1, d2
405	clt3D3(i,j) = var_clt(cldfra3D(i,j,:), d3)
406	END DO
407	END DO
408	ELSE
409	PRINT *,TRIM(ErrWarnMsg('err'))
410	PRINT *,' ' // TRIM(fname) // ': wrong zdim:', zdim,' for Ndim=', Ndim, ' !!'
411	PRINT *,' accepted values: 1,2,3'
412	STOP
413	END IF
414	CASE (4)
415	IF (zdim == 1) THEN
416	DO i=1, d2
417	DO j=1, d3
418	DO k=1, d4
419	clt4D1(i,j,k) = var_clt(cldfra4D(:,i,j,k), d1)
420	END DO
421	END DO
422	END DO
423	ELSE IF (zdim == 2) THEN
424	DO i=1, d1
425	DO j=1, d3
426	DO k=1, d4
427	clt4D2(i,j,k) = var_clt(cldfra4D(i,:,j,k), d2)
428	END DO
429	END DO
430	END DO
431	ELSE IF (zdim == 3) THEN
432	DO i=1, d2
433	DO j=1, d3
434	DO k=1, d4
435	clt4D3(i,j,k) = var_clt(cldfra4D(i,j,:,k), d3)
436	END DO
437	END DO
438	END DO
439	ELSE IF (zdim == 4) THEN
440	DO i=1, d1
441	DO j=1, d2
442	DO k=1, d3
443	clt4D4(i,j,k) = var_clt(cldfra4D(i,j,k,:), d4)
444	END DO
445	END DO
446	END DO
447	ELSE
448	PRINT *,TRIM(ErrWarnMsg('err'))
449	PRINT *,' ' // TRIM(fname) // ': wrong zdim:', zdim,' for Ndim=', Ndim, ' !!'
450	PRINT *,' accepted values: 1,2,3,4'
451	STOP
452	END IF
453	CASE DEFAULT
454	PRINT *,TRIM(ErrWarnMsg('err'))
455	PRINT *,' ' // TRIM(fname) // ': Ndim:', Ndim,' not ready !!'
456	STOP
457	END SELECT
458
459	RETURN
460
461	END SUBROUTINE compute_clt
462
463	SUBROUTINE compute_massvertint1D(var, mutot, dz, deta, integral)
464	! Subroutine to vertically integrate a 1D variable in eta vertical coordinates
465
466	IMPLICIT NONE
467
468	INTEGER, INTENT(in) :: dz
469	REAL(r_k), INTENT(in) :: mutot
470	REAL(r_k), DIMENSION(dz), INTENT(in) :: var, deta
471	REAL(r_k), INTENT(out) :: integral
472
473	! Local
474	INTEGER :: k
475
476	!!!!!!! Variables
477	! var: vertical variable to integrate (assuming kgkg-1)
478	! mutot: total dry-air mass in column
479	! dz: vertical dimension
480	! deta: eta-levels difference between full eta-layers
481
482	fname = 'compute_massvertint1D'
483
484	! integral=0.
485	! DO k=1,dz
486	! integral = integral + var(k)*deta(k)
487	! END DO
488	integral = SUM(var*deta)
489
490	integral=integral*mutot/g
491
492	RETURN
493
494	END SUBROUTINE compute_massvertint1D
495
496	SUBROUTINE compute_zint4D(var4D, dlev, zweight, d1, d2, d3, d4, int3D)
497	! Subroutine to vertically integrate a 4D variable in any vertical coordinates
498
499	IMPLICIT NONE
500
501	INTEGER, INTENT(in) :: d1,d2,d3,d4
502	REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in) :: var4D, dlev, zweight
503	REAL(r_k), DIMENSION(d1,d2,d4), INTENT(out) :: int3D
504
505	! Local
506	INTEGER :: i,j,l
507
508	!!!!!!! Variables
509	! var4D: vertical variable to integrate
510	! dlev: height of layers
511	! zweight: weight for each level to be applied (=1. for no effect)
512
513	fname = 'compute_zint4D'
514
515	DO i=1,d1
516	DO j=1,d2
517	DO l=1,d4
518	CALL compute_vertint1D(var4D(i,j,:,l),d3, dlev(i,j,:,l), zweight(i,j,:,l), &
519	int3D(i,j,l))
520	END DO
521	END DO
522	END DO
523
524	RETURN
525
526	END SUBROUTINE compute_zint4D
527
528	SUBROUTINE compute_vertint1D(var, dz, deta, zweight, integral)
529	! Subroutine to vertically integrate a 1D variable in any vertical coordinates
530
531	IMPLICIT NONE
532
533	INTEGER, INTENT(in) :: dz
534	REAL(r_k), DIMENSION(dz), INTENT(in) :: var, deta, zweight
535	REAL(r_k), INTENT(out) :: integral
536
537	! Local
538	INTEGER :: k
539
540	!!!!!!! Variables
541	! var: vertical variable to integrate
542	! dz: vertical dimension
543	! deta: eta-levels difference between layers
544	! zweight: weight for each level to be applied (=1. for no effect)
545
546	fname = 'compute_vertint1D'
547
548	! integral=0.
549	! DO k=1,dz
550	! integral = integral + var(k)*deta(k)
551	! END DO
552	integral = SUM(vardetazweight)
553
554	RETURN
555
556	END SUBROUTINE compute_vertint1D
557
558	SUBROUTINE compute_cape_afwa4D(ta, hur, press, zg, hgt, cape, cin, zlfc, plfc, li, parcelmethod, &
559	d1, d2, d3, d4)
560	! Subroutine to use WRF phys/module_diag_afwa.F `buyoancy' subroutine to compute CAPE, CIN, ZLFC, PLFC, LI
561
562	IMPLICIT NONE
563
564	INTEGER, INTENT(in) :: d1, d2, d3, d4, parcelmethod
565	REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in) :: ta, hur, press, zg
566	REAL(r_k), DIMENSION(d1,d2), INTENT(in) :: hgt
567	REAL(r_k), DIMENSION(d1,d2,d4), INTENT(out) :: cape, cin, zlfc, plfc, li
568
569	! Local
570	INTEGER :: i, j, it
571	INTEGER :: ofunc
572
573	!!!!!!! Variables
574	! ta: air temperature [K]
575	! hur: relative humidity [%]
576	! press: air pressure [Pa]
577	! zg: geopotential height [gpm]
578	! hgt: topographical height [m]
579	! cape: Convective available potential energy [Jkg-1]
580	! cin: Convective inhibition [Jkg-1]
581	! zlfc: height at the Level of free convection [m]
582	! plfc: pressure at the Level of free convection [Pa]
583	! li: lifted index [1]
584	! parcelmethod:
585	! Most Unstable = 1 (default)
586	! Mean layer = 2
587	! Surface based = 3
588
589	fname = 'compute_cape_afwa4D'
590
591	DO i=1, d1
592	DO j=1, d2
593	DO it=1, d4
594	ofunc = var_cape_afwa1D(d3, ta(i,j,:,it), hur(i,j,:,it), press(i,j,:,it), zg(i,j,:,it), &
595	1, cape(i,j,it), cin(i,j,it), zlfc(i,j,it), plfc(i,j,it), li(i,j,it), parcelmethod)
596	IF (zlfc(i,j,it) /= -1.) zlfc(i,j,it) = zlfc(i,j,it) - hgt(i,j)
597	END DO
598	END DO
599	END DO
600
601	RETURN
602
603	END SUBROUTINE compute_cape_afwa4D
604
605	SUBROUTINE compute_psl_ecmwf(ps, hgt, T, press, unpress, psl, d1, d2, d4)
606	! Subroutine to compute sea level pressure using ECMWF method following Mats Hamrud and Philippe Courtier [Pa]
607
608	IMPLICIT NONE
609
610	INTEGER, INTENT(in) :: d1, d2, d4
611	REAL(r_k), DIMENSION(d1,d2,d4), INTENT(in) :: ps, T, press, unpress
612	REAL(r_k), DIMENSION(d1,d2), INTENT(in) :: hgt
613	REAL(r_k), DIMENSION(d1,d2,d4), INTENT(out) :: psl
614
615	! Local
616	INTEGER :: i, j, it
617
618	!!!!!!! Variables
619	! ps: surface pressure [Pa]
620	! hgt: terrain height [m]
621	! T: temperature at first half-mass level [K]
622	! press: pressure at first full levels [Pa]
623	! unpress: pressure at first mass (half) levels [Pa]
624	! psl: sea-level pressure [Pa]
625
626	fname = 'compute_psl_ecmwf'
627
628	DO i=1, d1
629	DO j=1, d2
630	DO it=1, d4
631	CALL var_psl_ecmwf(ps(i,j,it), hgt(i,j), T(i,j,it), unpress(i,j,it), press(i,j,it), &
632	psl(i,j,it))
633	END DO
634	END DO
635	END DO
636
637	RETURN
638
639	END SUBROUTINE compute_psl_ecmwf
640
641	SUBROUTINE compute_zmla_generic4D(tpot, qratio, z, hgt, zmla3D, d1, d2, d3, d4)
642	! Subroutine to compute pbl-height following a generic method
643	! from Nielsen-Gammon et al., 2008 J. Appl. Meteor. Clim.
644	! applied also in Garcia-Diez et al., 2013, QJRMS
645	! where
646	! "The technique identifies the ML height as a threshold increase of potential temperature from
647	! its minimum value within the boundary layer."
648	! here applied similarly to Garcia-Diez et al. where
649	! zmla = "...first level where potential temperature exceeds the minimum potential temperature
650	! reached in the mixed layer by more than 1.5 K"
651
652	IMPLICIT NONE
653
654	INTEGER, INTENT(in) :: d1, d2, d3, d4
655	REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in) :: tpot, qratio, z
656	REAL(r_k), DIMENSION(d1,d2), INTENT(in) :: hgt
657	REAL(r_k), DIMENSION(d1,d2,d4), INTENT(out) :: zmla3D
658
659	! Local
660	INTEGER :: i, j, it
661
662	!!!!!!! Variables
663	! tpot: potential air temperature [K]
664	! qratio: water vapour mixing ratio [kgkg-1]
665	! z: height above sea level [m]
666	! hgt: terrain height [m]
667	! zmla3D: boundary layer height from surface [m]
668
669	fname = 'compute_zmla_generic4D'
670
671	DO i=1, d1
672	DO j=1, d2
673	DO it=1, d4
674	CALL var_zmla_generic(d3, qratio(i,j,:,it), tpot(i,j,:,it), z(i,j,:,it), hgt(i,j), &
675	zmla3D(i,j,it))
676	END DO
677	END DO
678	END DO
679
680	RETURN
681
682	END SUBROUTINE compute_zmla_generic4D
683
684	SUBROUTINE compute_zwind4D(ua, va, z, uas, vas, sina, cosa, zextrap, uaz, vaz, d1, d2, d3, d4)
685	! Subroutine to compute extrapolate the wind at a given height following the 'power law' methodology
686
687	IMPLICIT NONE
688
689	INTEGER, INTENT(in) :: d1, d2, d3, d4
690	REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in) :: ua, va, z
691	REAL(r_k), DIMENSION(d1,d2,d4), INTENT(in) :: uas, vas
692	REAL(r_k), DIMENSION(d1,d2), INTENT(in) :: sina, cosa
693	REAL(r_k), INTENT(in) :: zextrap
694	REAL(r_k), DIMENSION(d1,d2,d4), INTENT(out) :: uaz, vaz
695
696	! Local
697	INTEGER :: i, j, it
698
699	!!!!!!! Variables
700	! tpot: potential air temperature [K]
701	! qratio: water vapour mixing ratio [kgkg-1]
702	! z: height above surface [m]
703	! sina, cosa: local sine and cosine of map rotation [1.]
704	! zmla3D: boundary layer height from surface [m]
705
706	fname = 'compute_zwind4D'
707
708	DO i=1, d1
709	DO j=1, d2
710	DO it=1, d4
711	CALL var_zwind(d3, ua(i,j,:,it), va(i,j,:,it), z(i,j,:,it), uas(i,j,it), vas(i,j,it), &
712	sina(i,j), cosa(i,j), zextrap, uaz(i,j,it), vaz(i,j,it))
713	END DO
714	END DO
715	END DO
716
717	RETURN
718
719	END SUBROUTINE compute_zwind4D
720
721	SUBROUTINE compute_zwind_log4D(ua, va, z, uas, vas, sina, cosa, zextrap, uaz, vaz, d1, d2, d3, d4)
722	! Subroutine to compute extrapolate the wind at a given height following the 'logarithmic law' methodology
723
724	IMPLICIT NONE
725
726	INTEGER, INTENT(in) :: d1, d2, d3, d4
727	REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in) :: ua, va, z
728	REAL(r_k), DIMENSION(d1,d2,d4), INTENT(in) :: uas, vas
729	REAL(r_k), DIMENSION(d1,d2), INTENT(in) :: sina, cosa
730	REAL(r_k), INTENT(in) :: zextrap
731	REAL(r_k), DIMENSION(d1,d2,d4), INTENT(out) :: uaz, vaz
732
733	! Local
734	INTEGER :: i, j, it
735
736	!!!!!!! Variables
737	! tpot: potential air temperature [K]
738	! qratio: water vapour mixing ratio [kgkg-1]
739	! z: height above surface [m]
740	! sina, cosa: local sine and cosine of map rotation [1.]
741	! zmla3D: boundary layer height from surface [m]
742
743	fname = 'compute_zwind_log4D'
744
745	DO i=1, d1
746	DO j=1, d2
747	DO it=1, d4
748	CALL var_zwind_log(d3, ua(i,j,:,it), va(i,j,:,it), z(i,j,:,it), uas(i,j,it), vas(i,j,it), &
749	sina(i,j), cosa(i,j), zextrap, uaz(i,j,it), vaz(i,j,it))
750	END DO
751	END DO
752	END DO
753
754	RETURN
755
756	END SUBROUTINE compute_zwind_log4D
757
758	SUBROUTINE compute_zwindMO3D(d1, d2, d3, ust, znt, rmol, uas, vas, sina, cosa, newz, uznew, vznew)
759	! Subroutine to compute extrapolate the wind at a given height following the 'power law' methodology
760	! NOTE: only usefull for newz < 80. m
761
762	IMPLICIT NONE
763
764	INTEGER, INTENT(in) :: d1, d2, d3
765	REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in) :: ust, znt, rmol
766	REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in) :: uas, vas
767	REAL(r_k), DIMENSION(d1,d2), INTENT(in) :: sina, cosa
768	REAL(r_k), INTENT(in) :: newz
769	REAL(r_k), DIMENSION(d1,d2,d3), INTENT(out) :: uznew, vznew
770
771	! Local
772	INTEGER :: i, j, it
773
774	!!!!!!! Variables
775	! ust: u* in similarity theory [ms-1]
776	! znt: thermal time-varying roughness length [m]
777	! rmol: Inverse of the Obukhov length [m-1]
778	! uas: x-component 10-m wind speed [ms-1]
779	! vas: y-component 10-m wind speed [ms-1]
780	! sina, cosa: local sine and cosine of map rotation [1.]
781
782	fname = 'compute_zwindMO3D'
783
784	DO i=1, d1
785	DO j=1, d2
786	DO it=1, d3
787	CALL var_zwind_MOtheor(ust(i,j,it), znt(i,j,it), rmol(i,j,it), uas(i,j,it), vas(i,j,it), &
788	sina(i,j), cosa(i,j), newz, uznew(i,j,it), vznew(i,j,it))
789	END DO
790	END DO
791	END DO
792
793	RETURN
794
795	END SUBROUTINE compute_zwindMO3D
796
797	SUBROUTINE compute_potevap_orPM3D(d1, d2, d3, rho1, ust, uas, vas, tas, ps, qv1, potevap)
798	! Subroutine to compute potential evapotranspiration Penman-Monteith formulation implemented in
799	! ORCHIDEE in src_sechiba/enerbil.f90
800
801	IMPLICIT NONE
802
803	INTEGER, INTENT(in) :: d1, d2, d3
804	REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in) :: rho1, ust, uas, vas, tas, ps, qv1
805	REAL(r_k), DIMENSION(d1,d2,d3), INTENT(out) :: potevap
806
807	! Local
808	INTEGER :: i, j, it
809
810	!!!!!!! Variables
811	! rho1: atsmophere density at the first layer [kgm-3]
812	! ust: u* in similarity theory [ms-1]
813	! uas: x-component 10-m wind speed [ms-1]
814	! vas: y-component 10-m wind speed [ms-1]
815	! tas: 2-m atmosphere temperature [K]
816	! ps: surface pressure [Pa]
817	! qv1: 1st layer atmospheric mixing ratio [kgkg-1]
818	! potevap: potential evapo transpiration [kgm-2s-1]
819
820	fname = 'compute_potevap_orPM3D'
821
822	DO i=1, d1
823	DO j=1, d2
824	DO it=1, d3
825	CALL var_potevap_orPM(rho1(i,j,it), ust(i,j,it), uas(i,j,it), vas(i,j,it), tas(i,j,it), &
826	ps(i,j,it), qv1(i,j,it), potevap(i,j,it))
827	END DO
828	END DO
829	END DO
830
831	RETURN
832
833	END SUBROUTINE compute_potevap_orPM3D
834
835	END MODULE module_ForDiagnostics

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

Download in other formats:

Original Format