Context Navigation

module_ForDiagnostics.f90 @ 2381

Last change on this file since 2381 was 2345, checked in by lfita, 6 years ago
Adding right filling value in variable `range' after homogenization
File size: 47.1 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_scientific
12	USE module_ForDiagnosticsVars
13
14	CONTAINS
15
16	!!!!!!! Calculations
17	! compute_cape_afwa4D: Subroutine to use WRF phys/module_diag_afwa.F `buyoancy' subroutine to compute
18	! CAPE, CIN, ZLFC, PLFC, LI
19	! compute_cellbnds: Subroutine to compute cellboundaries using wind-staggered lon, lats as
20	! intersection of their related parallels and meridians
21	! compute_cellbndsreg: Subroutine to compute cellboundaries using lon, lat from a reglar lon/lat
22	! projection as intersection of their related parallels and meridians
23	! compute_cllmh4D3: Computation of low, medium and high cloudiness from a 4D CLDFRA and pressure being
24	! 3rd dimension the z-dim
25	! compute_cllmh3D3: Computation of low, medium and high cloudiness from a 3D CLDFRA and pressure being
26	! 3rd dimension the z-dim
27	! compute_cllmh: Computation of low, medium and high cloudiness
28	! compute_clt4D3: Computation of total cloudiness from a 4D CLDFRA being 3rd dimension the z-dim
29	! compute_clt3D3: Computation of total cloudiness from a 3D CLDFRA being 3rd dimension the z-dim
30	! compute_clt: Computation of total cloudiness
31	! compute_fog_K84: Computation of fog and visibility following Kunkel, (1984)
32	! compute_fog_RUC: Computation of fog and visibility following RUC method Smirnova, (2000)
33	! compute_fog_FRAML50: fog and visibility following Gultepe and Milbrandt, (2010)
34	! compute_massvertint1D: Subroutine to vertically integrate a 1D variable in eta vertical coordinates
35	! compute_psl_ecmwf: Compute sea level pressure using ECMWF method following Mats Hamrud and Philippe Courtier [Pa]
36	! compute_range_faces: Subroutine to compute faces [uphill, valleys, downhill] of a mountain range along a given face
37	! compute_vertint1D: Subroutine to vertically integrate a 1D variable in any vertical coordinates
38	! compute_zint4D: Subroutine to vertically integrate a 4D variable in any vertical coordinates
39	! compute_zmla_generic4D: Subroutine to compute pbl-height following a generic method
40	! compute_zwind4D: Subroutine to compute extrapolate the wind at a given height following the 'power law' methodology
41	! compute_zwind_log4D: Subroutine to compute extrapolate the wind at a given height following the 'logarithmic law' methodology
42	! compute_zwindMCO3D: Subroutine to compute extrapolate the wind at a given height following the 'power law' methodolog
43
44	!!!
45	! Calculations
46	!!!
47
48	SUBROUTINE compute_cllmh4D2(cldfra4D, pres4D, cllmh4D2, d1, d2, d3, d4)
49	! Subroutine to compute the low, medium and high cloudiness following 'newmicro.F90' from LMDZ from a 4D CLDFRA and pressure
50	! where zdim is the 2nd dimension (thus, cldfra4D(d1,d2,d3,d4) --> cllmh(3,d1,d3,d4) 1: low, 2: medium, 3: high
51	! It should be properly done via an 'INTERFACE', but...
52
53	IMPLICIT NONE
54
55	INTEGER, INTENT(in) :: d1, d2, d3, d4
56	REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in) :: cldfra4D, pres4D
57	REAL(r_k), DIMENSION(3,d1,d3,d4), INTENT(out) :: cllmh4D2
58
59	! Local
60	INTEGER :: i,j,k, zdim, Ndim
61
62	!!!!!!! Variables
63	! cldfra4D: 4D cloud fraction values [1]
64	! pres4D: 4D pressure values [Pa]
65	! Ndim: number of dimensions of the input data
66	! d[1-4]: dimensions of 'cldfra4D'
67	! zdim: number of the vertical-dimension within the matrix
68	! cltlmh4D2: low, medium, high cloudiness for the 4D cldfra and d2 being 'zdim'
69
70	fname = 'compute_cllmh4D2'
71	zdim = 2
72	Ndim = 4
73
74	DO i=1, d1
75	DO j=1, d3
76	DO k=1, d4
77	cllmh4D2(:,i,j,k) = var_cllmh(cldfra4D(i,:,j,k), pres4D(i,:,j,k), d2)
78	END DO
79	END DO
80	END DO
81
82	RETURN
83
84	END SUBROUTINE compute_cllmh4D2
85
86	SUBROUTINE compute_cllmh3D1(cldfra3D, pres3D, cllmh3D1, d1, d2, d3)
87	! Subroutine to compute the low, medium and high cloudiness following 'newmicro.F90' from LMDZ from a 3D CLDFRA and pressure
88	! where zdim is the 1st dimension (thus, cldfra3D(d1,d2,d3) --> cllmh(3,d2,d3) 1: low, 2: medium, 3: high
89	! It should be properly done via an 'INTERFACE', but...
90
91	IMPLICIT NONE
92
93	INTEGER, INTENT(in) :: d1, d2, d3
94	REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in) :: cldfra3D, pres3D
95	REAL(r_k), DIMENSION(3,d2,d3), INTENT(out) :: cllmh3D1
96
97	! Local
98	INTEGER :: i,j,k, zdim, Ndim
99
100	!!!!!!! Variables
101	! cldfra3D: 3D cloud fraction values [1]
102	! pres3D: 3D pressure values [Pa]
103	! Ndim: number of dimensions of the input data
104	! d[1-3]: dimensions of 'cldfra3D'
105	! zdim: number of the vertical-dimension within the matrix
106	! cltlmh3D1: low, medium, high cloudiness for the 3D cldfra and d1 being 'zdim'
107
108	fname = 'compute_cllmh3D1'
109	zdim = 1
110	Ndim = 3
111
112	DO i=1, d1
113	DO j=1, d2
114	cllmh3D1(:,i,j) = var_cllmh(cldfra3D(:,i,j), pres3D(:,i,j), d1)
115	END DO
116	END DO
117
118	RETURN
119
120	END SUBROUTINE compute_cllmh3D1
121
122	SUBROUTINE compute_cllmh(cldfra1D, cldfra2D, cldfra3D, cldfra4D, pres1D, pres2D, pres3D, pres4D, &
123	Ndim, zdim, cllmh1D, cllmh2D1, cllmh2D2, cllmh3D1, cllmh3D2, cllmh3D3, cllmh4D1, cllmh4D2, &
124	cllmh4D3, cllmh4D4, d1, d2, d3, d4)
125	! Subroutine to compute the low, medium and high cloudiness following 'newmicro.F90' from LMDZ
126
127	IMPLICIT NONE
128
129	INTEGER, INTENT(in) :: Ndim, d1, d2, d3, d4, zdim
130	REAL(r_k), DIMENSION(d1), OPTIONAL, INTENT(in) :: cldfra1D, pres1D
131	REAL(r_k), DIMENSION(d1,d2), OPTIONAL, INTENT(in) :: cldfra2D, pres2D
132	REAL(r_k), DIMENSION(d1,d2,d3), OPTIONAL, INTENT(in) :: cldfra3D, pres3D
133	REAL(r_k), DIMENSION(d1,d2,d3,d4), OPTIONAL, &
134	INTENT(in) :: cldfra4D, pres4D
135	REAL(r_k), DIMENSION(3), OPTIONAL, INTENT(out) :: cllmh1D
136	REAL(r_k), DIMENSION(d1,3), OPTIONAL, INTENT(out) :: cllmh2D1
137	REAL(r_k), DIMENSION(d2,3), OPTIONAL, INTENT(out) :: cllmh2D2
138	REAL(r_k), DIMENSION(d2,d3,3), OPTIONAL, INTENT(out) :: cllmh3D1
139	REAL(r_k), DIMENSION(d1,d3,3), OPTIONAL, INTENT(out) :: cllmh3D2
140	REAL(r_k), DIMENSION(d1,d2,3), OPTIONAL, INTENT(out) :: cllmh3D3
141	REAL(r_k), DIMENSION(d2,d3,d4,3), OPTIONAL, &
142	INTENT(out) :: cllmh4D1
143	REAL(r_k), DIMENSION(d1,d3,d4,3), OPTIONAL, &
144	INTENT(out) :: cllmh4D2
145	REAL(r_k), DIMENSION(d1,d2,d4,3), OPTIONAL, &
146	INTENT(out) :: cllmh4D3
147	REAL(r_k), DIMENSION(d1,d2,d3,3), OPTIONAL, &
148	INTENT(out) :: cllmh4D4
149
150	! Local
151	INTEGER :: i,j,k
152
153	!!!!!!! Variables
154	! cldfra[1-4]D: cloud fraction values [1]
155	! pres[1-4]D: pressure values [Pa]
156	! Ndim: number of dimensions of the input data
157	! d[1-4]: dimensions of 'cldfra'
158	! zdim: number of the vertical-dimension within the matrix
159	! cllmh1D: low, medium and high cloudiness for the 1D cldfra
160	! cllmh2D1: low, medium and high cloudiness for the 2D cldfra and d1 being 'zdim'
161	! cllmh2D2: low, medium and high cloudiness for the 2D cldfra and d2 being 'zdim'
162	! cllmh3D1: low, medium and high cloudiness for the 3D cldfra and d1 being 'zdim'
163	! cllmh3D2: low, medium and high cloudiness for the 3D cldfra and d2 being 'zdim'
164	! cllmh3D3: low, medium and high cloudiness for the 3D cldfra and d3 being 'zdim'
165	! cllmh4D1: low, medium and high cloudiness for the 4D cldfra and d1 being 'zdim'
166	! cllmh4D2: low, medium and high cloudiness for the 4D cldfra and d2 being 'zdim'
167	! cllmh4D3: low, medium and high cloudiness for the 4D cldfra and d3 being 'zdim'
168	! cllmh4D4: low, medium and high cloudiness for the 4D cldfra and d4 being 'zdim'
169
170	fname = 'compute_cllmh'
171
172	SELECT CASE (Ndim)
173	CASE (1)
174	cllmh1D = var_cllmh(cldfra1D, pres1D, d1)
175	CASE (2)
176	IF (zdim == 1) THEN
177	DO i=1, d2
178	cllmh2D1(i,:) = var_cllmh(cldfra2D(:,i), pres2D(:,i), d1)
179	END DO
180	ELSE IF (zdim == 2) THEN
181	DO i=1, d1
182	cllmh2D2(i,:) = var_cllmh(cldfra2D(:,i), pres2D(i,:), d2)
183	END DO
184	ELSE
185	PRINT *,TRIM(ErrWarnMsg('err'))
186	PRINT *,' ' // TRIM(fname) // ': wrong zdim:', zdim,' for Ndim=', Ndim, ' !!'
187	PRINT *,' accepted values: 1,2'
188	STOP
189	END IF
190	CASE (3)
191	IF (zdim == 1) THEN
192	DO i=1, d2
193	DO j=1, d3
194	cllmh3D1(i,j,:) = var_cllmh(cldfra3D(:,i,j), pres3D(:,i,j), d1)
195	END DO
196	END DO
197	ELSE IF (zdim == 2) THEN
198	DO i=1, d1
199	DO j=1, d3
200	cllmh3D2(i,j,:) = var_cllmh(cldfra3D(i,:,j), pres3D(i,:,j), d2)
201	END DO
202	END DO
203	ELSE IF (zdim == 3) THEN
204	DO i=1, d1
205	DO j=1, d2
206	cllmh3D3(i,j,:) = var_cllmh(cldfra3D(i,j,:), pres3D(i,j,:), d3)
207	END DO
208	END DO
209	ELSE
210	PRINT *,TRIM(ErrWarnMsg('err'))
211	PRINT *,' ' // TRIM(fname) // ': wrong zdim:', zdim,' for Ndim=', Ndim, ' !!'
212	PRINT *,' accepted values: 1,2,3'
213	STOP
214	END IF
215	CASE (4)
216	IF (zdim == 1) THEN
217	DO i=1, d2
218	DO j=1, d3
219	DO k=1, d4
220	cllmh4D1(i,j,k,:) = var_cllmh(cldfra4D(:,i,j,k), pres4D(:,i,j,k), d1)
221	END DO
222	END DO
223	END DO
224	ELSE IF (zdim == 2) THEN
225	DO i=1, d1
226	DO j=1, d3
227	DO k=1, d4
228	cllmh4D2(i,j,k,:) = var_cllmh(cldfra4D(i,:,j,k), pres4D(i,:,j,k), d2)
229	END DO
230	END DO
231	END DO
232	ELSE IF (zdim == 3) THEN
233	DO i=1, d2
234	DO j=1, d3
235	DO k=1, d4
236	cllmh4D3(i,j,k,:) = var_cllmh(cldfra4D(i,j,:,k), pres4D(i,j,:,k), d3)
237	END DO
238	END DO
239	END DO
240	ELSE IF (zdim == 4) THEN
241	DO i=1, d1
242	DO j=1, d2
243	DO k=1, d3
244	cllmh4D4(i,j,k,:) = var_cllmh(cldfra4D(i,j,k,:), pres4D(i,j,k,:), d4)
245	END DO
246	END DO
247	END DO
248	ELSE
249	PRINT *,TRIM(ErrWarnMsg('err'))
250	PRINT *,' ' // TRIM(fname) // ': wrong zdim:', zdim,' for Ndim=', Ndim, ' !!'
251	PRINT *,' accepted values: 1,2,3,4'
252	STOP
253	END IF
254	CASE DEFAULT
255	PRINT *,TRIM(ErrWarnMsg('err'))
256	PRINT *,' ' // TRIM(fname) // ': Ndim:', Ndim,' not ready !!'
257	STOP
258	END SELECT
259
260	RETURN
261
262	END SUBROUTINE compute_cllmh
263
264	SUBROUTINE compute_clt4D2(cldfra4D, clt4D2, d1, d2, d3, d4)
265	! Subroutine to compute the total cloudiness following 'newmicro.F90' from LMDZ from a 4D CLDFRA
266	! where zdim is the 2nd dimension (thus, cldfra4D(d1,d2,d3,d4) --> clt(d1,d3,d4)
267	! It should be properly done via an 'INTERFACE', but...
268
269	IMPLICIT NONE
270
271	INTEGER, INTENT(in) :: d1, d2, d3, d4
272	REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in) :: cldfra4D
273	REAL(r_k), DIMENSION(d1,d3,d4), INTENT(out) :: clt4D2
274
275	! Local
276	INTEGER :: i,j,k, zdim, Ndim
277
278	!!!!!!! Variables
279	! cldfra4D: 4D cloud fraction values [1]
280	! Ndim: number of dimensions of the input data
281	! d[1-4]: dimensions of 'cldfra4D'
282	! zdim: number of the vertical-dimension within the matrix
283	! clt4D2: total cloudiness for the 4D cldfra and d2 being 'zdim'
284
285	fname = 'compute_clt4D2'
286	zdim = 2
287	Ndim = 4
288
289	DO i=1, d1
290	DO j=1, d3
291	DO k=1, d4
292	clt4D2(i,j,k) = var_clt(cldfra4D(i,:,j,k), d2)
293	END DO
294	END DO
295	END DO
296
297	RETURN
298
299	END SUBROUTINE compute_clt4D2
300
301	SUBROUTINE compute_clt3D1(cldfra3D, clt3D1, d1, d2, d3)
302	! Subroutine to compute the total cloudiness following 'newmicro.F90' from LMDZ from a 3D CLDFRA
303	! where zdim is the 1st dimension (thus, cldfra4D(d1,d2,d3) --> clt(d2,d3)
304	! It should be properly done via an 'INTERFACE', but...
305
306	IMPLICIT NONE
307
308	INTEGER, INTENT(in) :: d1, d2, d3
309	REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in) :: cldfra3D
310	REAL(r_k), DIMENSION(d2,d3), INTENT(out) :: clt3D1
311
312	! Local
313	INTEGER :: i,j,k, zdim, Ndim
314
315	!!!!!!! Variables
316	! cldfra3D: 3D cloud fraction values [1]
317	! Ndim: number of dimensions of the input data
318	! d[1-3]: dimensions of 'cldfra3D'
319	! zdim: number of the vertical-dimension within the matrix
320	! clt3D1: total cloudiness for the 3D cldfra and d1 being 'zdim'
321
322	fname = 'compute_clt3D1'
323	zdim = 1
324	Ndim = 3
325
326	DO i=1, d2
327	DO j=1, d3
328	clt3D1(i,j) = var_clt(cldfra3D(:,i,j), d1)
329	END DO
330	END DO
331
332	RETURN
333
334	END SUBROUTINE compute_clt3D1
335
336	SUBROUTINE compute_clt(cldfra1D, cldfra2D, cldfra3D, cldfra4D, Ndim, zdim, clt1D, clt2D1, clt2D2, &
337	clt3D1, clt3D2, clt3D3, clt4D1, clt4D2, clt4D3, clt4D4, d1, d2, d3, d4)
338	! Subroutine to compute the total cloudiness following 'newmicro.F90' from LMDZ
339
340	IMPLICIT NONE
341
342	INTEGER, INTENT(in) :: Ndim, d1, d2, d3, d4, zdim
343	REAL(r_k), DIMENSION(d1), OPTIONAL, INTENT(in) :: cldfra1D
344	REAL(r_k), DIMENSION(d1,d2), OPTIONAL, INTENT(in) :: cldfra2D
345	REAL(r_k), DIMENSION(d1,d2,d3), OPTIONAL, INTENT(in) :: cldfra3D
346	REAL(r_k), DIMENSION(d1,d2,d3,d4), OPTIONAL, &
347	INTENT(in) :: cldfra4D
348	REAL(r_k), OPTIONAL, INTENT(out) :: clt1D
349	REAL(r_k), DIMENSION(d1), OPTIONAL, INTENT(out) :: clt2D1
350	REAL(r_k), DIMENSION(d2), OPTIONAL, INTENT(out) :: clt2D2
351	REAL(r_k), DIMENSION(d2,d3), OPTIONAL, INTENT(out) :: clt3D1
352	REAL(r_k), DIMENSION(d1,d3), OPTIONAL, INTENT(out) :: clt3D2
353	REAL(r_k), DIMENSION(d1,d2), OPTIONAL, INTENT(out) :: clt3D3
354	REAL(r_k), DIMENSION(d2,d3,d4), OPTIONAL,INTENT(out) :: clt4D1
355	REAL(r_k), DIMENSION(d1,d3,d4), OPTIONAL,INTENT(out) :: clt4D2
356	REAL(r_k), DIMENSION(d1,d2,d4), OPTIONAL,INTENT(out) :: clt4D3
357	REAL(r_k), DIMENSION(d1,d2,d3), OPTIONAL,INTENT(out) :: clt4D4
358
359	! Local
360	INTEGER :: i,j,k
361
362	!!!!!!! Variables
363	! cldfra[1-4]D: cloud fraction values [1]
364	! Ndim: number of dimensions of the input data
365	! d[1-4]: dimensions of 'cldfra'
366	! zdim: number of the vertical-dimension within the matrix
367	! clt1D: total cloudiness for the 1D cldfra
368	! clt2D1: total cloudiness for the 2D cldfra and d1 being 'zdim'
369	! clt2D2: total cloudiness for the 2D cldfra and d2 being 'zdim'
370	! clt3D1: total cloudiness for the 3D cldfra and d1 being 'zdim'
371	! clt3D2: total cloudiness for the 3D cldfra and d2 being 'zdim'
372	! clt3D3: total cloudiness for the 3D cldfra and d3 being 'zdim'
373	! clt4D1: total cloudiness for the 4D cldfra and d1 being 'zdim'
374	! clt4D2: total cloudiness for the 4D cldfra and d2 being 'zdim'
375	! clt4D3: total cloudiness for the 4D cldfra and d3 being 'zdim'
376	! clt4D4: total cloudiness for the 4D cldfra and d4 being 'zdim'
377
378	fname = 'compute_clt'
379
380	SELECT CASE (Ndim)
381	CASE (1)
382	clt1D = var_clt(cldfra1D, d1)
383	CASE (2)
384	IF (zdim == 1) THEN
385	DO i=1, d2
386	clt2D1(i) = var_clt(cldfra2D(:,i), d1)
387	END DO
388	ELSE IF (zdim == 2) THEN
389	DO i=1, d1
390	clt2D2(i) = var_clt(cldfra2D(:,i), d2)
391	END DO
392	ELSE
393	PRINT *,TRIM(ErrWarnMsg('err'))
394	PRINT *,' ' // TRIM(fname) // ': wrong zdim:', zdim,' for Ndim=', Ndim, ' !!'
395	PRINT *,' accepted values: 1,2'
396	STOP
397	END IF
398	CASE (3)
399	IF (zdim == 1) THEN
400	DO i=1, d2
401	DO j=1, d3
402	clt3D1(i,j) = var_clt(cldfra3D(:,i,j), d1)
403	END DO
404	END DO
405	ELSE IF (zdim == 2) THEN
406	DO i=1, d1
407	DO j=1, d3
408	clt3D2(i,j) = var_clt(cldfra3D(i,:,j), d2)
409	END DO
410	END DO
411	ELSE IF (zdim == 3) THEN
412	DO i=1, d1
413	DO j=1, d2
414	clt3D3(i,j) = var_clt(cldfra3D(i,j,:), d3)
415	END DO
416	END DO
417	ELSE
418	PRINT *,TRIM(ErrWarnMsg('err'))
419	PRINT *,' ' // TRIM(fname) // ': wrong zdim:', zdim,' for Ndim=', Ndim, ' !!'
420	PRINT *,' accepted values: 1,2,3'
421	STOP
422	END IF
423	CASE (4)
424	IF (zdim == 1) THEN
425	DO i=1, d2
426	DO j=1, d3
427	DO k=1, d4
428	clt4D1(i,j,k) = var_clt(cldfra4D(:,i,j,k), d1)
429	END DO
430	END DO
431	END DO
432	ELSE IF (zdim == 2) THEN
433	DO i=1, d1
434	DO j=1, d3
435	DO k=1, d4
436	clt4D2(i,j,k) = var_clt(cldfra4D(i,:,j,k), d2)
437	END DO
438	END DO
439	END DO
440	ELSE IF (zdim == 3) THEN
441	DO i=1, d2
442	DO j=1, d3
443	DO k=1, d4
444	clt4D3(i,j,k) = var_clt(cldfra4D(i,j,:,k), d3)
445	END DO
446	END DO
447	END DO
448	ELSE IF (zdim == 4) THEN
449	DO i=1, d1
450	DO j=1, d2
451	DO k=1, d3
452	clt4D4(i,j,k) = var_clt(cldfra4D(i,j,k,:), d4)
453	END DO
454	END DO
455	END DO
456	ELSE
457	PRINT *,TRIM(ErrWarnMsg('err'))
458	PRINT *,' ' // TRIM(fname) // ': wrong zdim:', zdim,' for Ndim=', Ndim, ' !!'
459	PRINT *,' accepted values: 1,2,3,4'
460	STOP
461	END IF
462	CASE DEFAULT
463	PRINT *,TRIM(ErrWarnMsg('err'))
464	PRINT *,' ' // TRIM(fname) // ': Ndim:', Ndim,' not ready !!'
465	STOP
466	END SELECT
467
468	RETURN
469
470	END SUBROUTINE compute_clt
471
472	SUBROUTINE compute_massvertint1D(var, mutot, dz, deta, integral)
473	! Subroutine to vertically integrate a 1D variable in eta vertical coordinates
474
475	IMPLICIT NONE
476
477	INTEGER, INTENT(in) :: dz
478	REAL(r_k), INTENT(in) :: mutot
479	REAL(r_k), DIMENSION(dz), INTENT(in) :: var, deta
480	REAL(r_k), INTENT(out) :: integral
481
482	! Local
483	INTEGER :: k
484
485	!!!!!!! Variables
486	! var: vertical variable to integrate (assuming kgkg-1)
487	! mutot: total dry-air mass in column
488	! dz: vertical dimension
489	! deta: eta-levels difference between full eta-layers
490
491	fname = 'compute_massvertint1D'
492
493	! integral=0.
494	! DO k=1,dz
495	! integral = integral + var(k)*deta(k)
496	! END DO
497	integral = SUM(var*deta)
498
499	integral=integral*mutot/g
500
501	RETURN
502
503	END SUBROUTINE compute_massvertint1D
504
505	SUBROUTINE compute_zint4D(var4D, dlev, zweight, d1, d2, d3, d4, int3D)
506	! Subroutine to vertically integrate a 4D variable in any vertical coordinates
507
508	IMPLICIT NONE
509
510	INTEGER, INTENT(in) :: d1,d2,d3,d4
511	REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in) :: var4D, dlev, zweight
512	REAL(r_k), DIMENSION(d1,d2,d4), INTENT(out) :: int3D
513
514	! Local
515	INTEGER :: i,j,l
516
517	!!!!!!! Variables
518	! var4D: vertical variable to integrate
519	! dlev: height of layers
520	! zweight: weight for each level to be applied (=1. for no effect)
521
522	fname = 'compute_zint4D'
523
524	DO i=1,d1
525	DO j=1,d2
526	DO l=1,d4
527	CALL compute_vertint1D(var4D(i,j,:,l),d3, dlev(i,j,:,l), zweight(i,j,:,l), &
528	int3D(i,j,l))
529	END DO
530	END DO
531	END DO
532
533	RETURN
534
535	END SUBROUTINE compute_zint4D
536
537	SUBROUTINE compute_vertint1D(var, dz, deta, zweight, integral)
538	! Subroutine to vertically integrate a 1D variable in any vertical coordinates
539
540	IMPLICIT NONE
541
542	INTEGER, INTENT(in) :: dz
543	REAL(r_k), DIMENSION(dz), INTENT(in) :: var, deta, zweight
544	REAL(r_k), INTENT(out) :: integral
545
546	! Local
547	INTEGER :: k
548
549	!!!!!!! Variables
550	! var: vertical variable to integrate
551	! dz: vertical dimension
552	! deta: eta-levels difference between layers
553	! zweight: weight for each level to be applied (=1. for no effect)
554
555	fname = 'compute_vertint1D'
556
557	! integral=0.
558	! DO k=1,dz
559	! integral = integral + var(k)*deta(k)
560	! END DO
561	integral = SUM(vardetazweight)
562
563	RETURN
564
565	END SUBROUTINE compute_vertint1D
566
567	SUBROUTINE compute_cape_afwa4D(ta, hur, press, zg, hgt, cape, cin, zlfc, plfc, li, parcelmethod, &
568	d1, d2, d3, d4)
569	! Subroutine to use WRF phys/module_diag_afwa.F `buyoancy' subroutine to compute CAPE, CIN, ZLFC, PLFC, LI
570
571	IMPLICIT NONE
572
573	INTEGER, INTENT(in) :: d1, d2, d3, d4, parcelmethod
574	REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in) :: ta, hur, press, zg
575	REAL(r_k), DIMENSION(d1,d2), INTENT(in) :: hgt
576	REAL(r_k), DIMENSION(d1,d2,d4), INTENT(out) :: cape, cin, zlfc, plfc, li
577
578	! Local
579	INTEGER :: i, j, it
580	INTEGER :: ofunc
581
582	!!!!!!! Variables
583	! ta: air temperature [K]
584	! hur: relative humidity [%]
585	! press: air pressure [Pa]
586	! zg: geopotential height [gpm]
587	! hgt: topographical height [m]
588	! cape: Convective available potential energy [Jkg-1]
589	! cin: Convective inhibition [Jkg-1]
590	! zlfc: height at the Level of free convection [m]
591	! plfc: pressure at the Level of free convection [Pa]
592	! li: lifted index [1]
593	! parcelmethod:
594	! Most Unstable = 1 (default)
595	! Mean layer = 2
596	! Surface based = 3
597
598	fname = 'compute_cape_afwa4D'
599
600	DO i=1, d1
601	DO j=1, d2
602	DO it=1, d4
603	ofunc = var_cape_afwa1D(d3, ta(i,j,:,it), hur(i,j,:,it), press(i,j,:,it), zg(i,j,:,it), &
604	1, cape(i,j,it), cin(i,j,it), zlfc(i,j,it), plfc(i,j,it), li(i,j,it), parcelmethod)
605	IF (zlfc(i,j,it) /= -1.) zlfc(i,j,it) = zlfc(i,j,it) - hgt(i,j)
606	END DO
607	END DO
608	END DO
609
610	RETURN
611
612	END SUBROUTINE compute_cape_afwa4D
613
614	SUBROUTINE compute_psl_ecmwf(ps, hgt, T, press, unpress, psl, d1, d2, d4)
615	! Subroutine to compute sea level pressure using ECMWF method following Mats Hamrud and Philippe Courtier [Pa]
616
617	IMPLICIT NONE
618
619	INTEGER, INTENT(in) :: d1, d2, d4
620	REAL(r_k), DIMENSION(d1,d2,d4), INTENT(in) :: ps, T, press, unpress
621	REAL(r_k), DIMENSION(d1,d2), INTENT(in) :: hgt
622	REAL(r_k), DIMENSION(d1,d2,d4), INTENT(out) :: psl
623
624	! Local
625	INTEGER :: i, j, it
626
627	!!!!!!! Variables
628	! ps: surface pressure [Pa]
629	! hgt: terrain height [m]
630	! T: temperature at first half-mass level [K]
631	! press: pressure at first full levels [Pa]
632	! unpress: pressure at first mass (half) levels [Pa]
633	! psl: sea-level pressure [Pa]
634
635	fname = 'compute_psl_ecmwf'
636
637	DO i=1, d1
638	DO j=1, d2
639	DO it=1, d4
640	CALL var_psl_ecmwf(ps(i,j,it), hgt(i,j), T(i,j,it), unpress(i,j,it), press(i,j,it), &
641	psl(i,j,it))
642	END DO
643	END DO
644	END DO
645
646	RETURN
647
648	END SUBROUTINE compute_psl_ecmwf
649
650	SUBROUTINE compute_zmla_generic4D(tpot, qratio, z, hgt, zmla3D, d1, d2, d3, d4)
651	! Subroutine to compute pbl-height following a generic method
652	! from Nielsen-Gammon et al., 2008 J. Appl. Meteor. Clim.
653	! applied also in Garcia-Diez et al., 2013, QJRMS
654	! where
655	! "The technique identifies the ML height as a threshold increase of potential temperature from
656	! its minimum value within the boundary layer."
657	! here applied similarly to Garcia-Diez et al. where
658	! zmla = "...first level where potential temperature exceeds the minimum potential temperature
659	! reached in the mixed layer by more than 1.5 K"
660
661	IMPLICIT NONE
662
663	INTEGER, INTENT(in) :: d1, d2, d3, d4
664	REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in) :: tpot, qratio, z
665	REAL(r_k), DIMENSION(d1,d2), INTENT(in) :: hgt
666	REAL(r_k), DIMENSION(d1,d2,d4), INTENT(out) :: zmla3D
667
668	! Local
669	INTEGER :: i, j, it
670
671	!!!!!!! Variables
672	! tpot: potential air temperature [K]
673	! qratio: water vapour mixing ratio [kgkg-1]
674	! z: height above sea level [m]
675	! hgt: terrain height [m]
676	! zmla3D: boundary layer height from surface [m]
677
678	fname = 'compute_zmla_generic4D'
679
680	DO i=1, d1
681	DO j=1, d2
682	DO it=1, d4
683	CALL var_zmla_generic(d3, qratio(i,j,:,it), tpot(i,j,:,it), z(i,j,:,it), hgt(i,j), &
684	zmla3D(i,j,it))
685	END DO
686	END DO
687	END DO
688
689	RETURN
690
691	END SUBROUTINE compute_zmla_generic4D
692
693	SUBROUTINE compute_zwind4D(ua, va, z, uas, vas, sina, cosa, zextrap, uaz, vaz, d1, d2, d3, d4)
694	! Subroutine to compute extrapolate the wind at a given height following the 'power law' methodology
695
696	IMPLICIT NONE
697
698	INTEGER, INTENT(in) :: d1, d2, d3, d4
699	REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in) :: ua, va, z
700	REAL(r_k), DIMENSION(d1,d2,d4), INTENT(in) :: uas, vas
701	REAL(r_k), DIMENSION(d1,d2), INTENT(in) :: sina, cosa
702	REAL(r_k), INTENT(in) :: zextrap
703	REAL(r_k), DIMENSION(d1,d2,d4), INTENT(out) :: uaz, vaz
704
705	! Local
706	INTEGER :: i, j, it
707
708	!!!!!!! Variables
709	! tpot: potential air temperature [K]
710	! qratio: water vapour mixing ratio [kgkg-1]
711	! z: height above surface [m]
712	! sina, cosa: local sine and cosine of map rotation [1.]
713	! zmla3D: boundary layer height from surface [m]
714
715	fname = 'compute_zwind4D'
716
717	DO i=1, d1
718	DO j=1, d2
719	DO it=1, d4
720	CALL var_zwind(d3, ua(i,j,:,it), va(i,j,:,it), z(i,j,:,it), uas(i,j,it), vas(i,j,it), &
721	sina(i,j), cosa(i,j), zextrap, uaz(i,j,it), vaz(i,j,it))
722	END DO
723	END DO
724	END DO
725
726	RETURN
727
728	END SUBROUTINE compute_zwind4D
729
730	SUBROUTINE compute_zwind_log4D(ua, va, z, uas, vas, sina, cosa, zextrap, uaz, vaz, d1, d2, d3, d4)
731	! Subroutine to compute extrapolate the wind at a given height following the 'logarithmic law' methodology
732
733	IMPLICIT NONE
734
735	INTEGER, INTENT(in) :: d1, d2, d3, d4
736	REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in) :: ua, va, z
737	REAL(r_k), DIMENSION(d1,d2,d4), INTENT(in) :: uas, vas
738	REAL(r_k), DIMENSION(d1,d2), INTENT(in) :: sina, cosa
739	REAL(r_k), INTENT(in) :: zextrap
740	REAL(r_k), DIMENSION(d1,d2,d4), INTENT(out) :: uaz, vaz
741
742	! Local
743	INTEGER :: i, j, it
744
745	!!!!!!! Variables
746	! tpot: potential air temperature [K]
747	! qratio: water vapour mixing ratio [kgkg-1]
748	! z: height above surface [m]
749	! sina, cosa: local sine and cosine of map rotation [1.]
750	! zmla3D: boundary layer height from surface [m]
751
752	fname = 'compute_zwind_log4D'
753
754	DO i=1, d1
755	DO j=1, d2
756	DO it=1, d4
757	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), &
758	sina(i,j), cosa(i,j), zextrap, uaz(i,j,it), vaz(i,j,it))
759	END DO
760	END DO
761	END DO
762
763	RETURN
764
765	END SUBROUTINE compute_zwind_log4D
766
767	SUBROUTINE compute_zwindMO3D(d1, d2, d3, ust, znt, rmol, uas, vas, sina, cosa, newz, uznew, vznew)
768	! Subroutine to compute extrapolate the wind at a given height following the 'power law' methodology
769	! NOTE: only usefull for newz < 80. m
770
771	IMPLICIT NONE
772
773	INTEGER, INTENT(in) :: d1, d2, d3
774	REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in) :: ust, znt, rmol
775	REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in) :: uas, vas
776	REAL(r_k), DIMENSION(d1,d2), INTENT(in) :: sina, cosa
777	REAL(r_k), INTENT(in) :: newz
778	REAL(r_k), DIMENSION(d1,d2,d3), INTENT(out) :: uznew, vznew
779
780	! Local
781	INTEGER :: i, j, it
782
783	!!!!!!! Variables
784	! ust: u* in similarity theory [ms-1]
785	! znt: thermal time-varying roughness length [m]
786	! rmol: Inverse of the Obukhov length [m-1]
787	! uas: x-component 10-m wind speed [ms-1]
788	! vas: y-component 10-m wind speed [ms-1]
789	! sina, cosa: local sine and cosine of map rotation [1.]
790
791	fname = 'compute_zwindMO3D'
792
793	DO i=1, d1
794	DO j=1, d2
795	DO it=1, d3
796	CALL var_zwind_MOtheor(ust(i,j,it), znt(i,j,it), rmol(i,j,it), uas(i,j,it), vas(i,j,it), &
797	sina(i,j), cosa(i,j), newz, uznew(i,j,it), vznew(i,j,it))
798	END DO
799	END DO
800	END DO
801
802	RETURN
803
804	END SUBROUTINE compute_zwindMO3D
805
806	SUBROUTINE compute_potevap_orPM3D(d1, d2, d3, rho1, ust, uas, vas, tas, ps, qv1, potevap)
807	! Subroutine to compute potential evapotranspiration Penman-Monteith formulation implemented in
808	! ORCHIDEE in src_sechiba/enerbil.f90
809
810	IMPLICIT NONE
811
812	INTEGER, INTENT(in) :: d1, d2, d3
813	REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in) :: rho1, ust, uas, vas, tas, ps, qv1
814	REAL(r_k), DIMENSION(d1,d2,d3), INTENT(out) :: potevap
815
816	! Local
817	INTEGER :: i, j, it
818
819	!!!!!!! Variables
820	! rho1: atsmophere density at the first layer [kgm-3]
821	! ust: u* in similarity theory [ms-1]
822	! uas: x-component 10-m wind speed [ms-1]
823	! vas: y-component 10-m wind speed [ms-1]
824	! tas: 2-m atmosphere temperature [K]
825	! ps: surface pressure [Pa]
826	! qv1: 1st layer atmospheric mixing ratio [kgkg-1]
827	! potevap: potential evapo transpiration [kgm-2s-1]
828
829	fname = 'compute_potevap_orPM3D'
830
831	DO i=1, d1
832	DO j=1, d2
833	DO it=1, d3
834	CALL var_potevap_orPM(rho1(i,j,it), ust(i,j,it), uas(i,j,it), vas(i,j,it), tas(i,j,it), &
835	ps(i,j,it), qv1(i,j,it), potevap(i,j,it))
836	END DO
837	END DO
838	END DO
839
840	RETURN
841
842	END SUBROUTINE compute_potevap_orPM3D
843
844	SUBROUTINE compute_fog_K84(d1, d2, d3, qc, qi, fog, vis)
845	! Subroutine to compute fog: qcloud + qice /= 0.
846	! And visibility following Kunkel, B. A., (1984): Parameterization of droplet terminal velocity and
847	! extinction coefficient in fog models. J. Climate Appl. Meteor., 23, 34â41.
848
849	IMPLICIT NONE
850
851	INTEGER, INTENT(in) :: d1, d2, d3
852	REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in) :: qc, qi
853	INTEGER, DIMENSION(d1,d2,d3), INTENT(out) :: fog
854	REAL(r_k), DIMENSION(d1,d2,d3), INTENT(out) :: vis
855
856	! Local
857	INTEGER :: i, j, it
858
859	!!!!!!! Variables
860	! qc: cloud mixing ratio [kgkg-1]
861	! qi, ice mixing ratio [kgkg-1]
862	! fog: presence of fog (1: yes, 0: no)
863	! vis: visibility within fog [km]
864
865	fname = 'compute_fog_K84'
866
867	DO i=1, d1
868	DO j=1, d2
869	DO it=1, d3
870	CALL var_fog_K84(qc(i,j,it), qi(i,j,it), fog(i,j,it), vis(i,j,it))
871	END DO
872	END DO
873	END DO
874
875	RETURN
876
877	END SUBROUTINE compute_fog_K84
878
879	SUBROUTINE compute_fog_RUC(d1, d2, d3, qv, ta, pres, fog, vis)
880	! Subroutine to compute fog: qcloud + qice /= 0.
881	! And visibility following RUC method Smirnova, T. G., S. G. Benjamin, and J. M. Brown, 2000: Case
882	! study verification of RUC/MAPS fog and visibility forecasts. Preprints, 9 th Conference on
883	! Aviation, Range, and Aerospace Meteorlogy, AMS, Orlando, FL, Sep. 2000. Paper#2.3, 6 pp.
884
885	IMPLICIT NONE
886
887	INTEGER, INTENT(in) :: d1, d2, d3
888	REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in) :: qv, ta, pres
889	INTEGER, DIMENSION(d1,d2,d3), INTENT(out) :: fog
890	REAL(r_k), DIMENSION(d1,d2,d3), INTENT(out) :: vis
891
892	! Local
893	INTEGER :: i, j, it
894
895	!!!!!!! Variables
896	! qv: water vapor mixing ratio [kgkg-1]
897	! ta: temperature [K]
898	! pres: pressure [Pa]
899	! fog: presence of fog (1: yes, 0: no)
900	! vis: visibility within fog [km]
901
902	fname = 'compute_fog_RUC'
903
904	DO i=1, d1
905	DO j=1, d2
906	DO it=1, d3
907	CALL var_fog_RUC(qv(i,j,it), ta(i,j,it), pres(i,j,it), fog(i,j,it), vis(i,j,it))
908	END DO
909	END DO
910	END DO
911
912	RETURN
913
914	END SUBROUTINE compute_fog_RUC
915
916	SUBROUTINE compute_fog_FRAML50(d1, d2, d3, qv, ta, pres, fog, vis)
917	! Subroutine to compute fog (vis < 1 km) and visibility following
918	! Gultepe, I. and J.A. Milbrandt, 2010: Probabilistic Parameterizations of Visibility Using
919	! Observations of Rain Precipitation Rate, Relative Humidity, and Visibility. J. Appl. Meteor.
920	! Climatol., 49, 36-46, https://doi.org/10.1175/2009JAMC1927.1
921	! Interest is focused on a 'general' fog/visibilty approach, thus the fit at 50 % of probability is
922	! chosen
923	! Effects from precipitation are not considered
924
925	IMPLICIT NONE
926
927	INTEGER, INTENT(in) :: d1, d2, d3
928	REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in) :: qv, ta, pres
929	INTEGER, DIMENSION(d1,d2,d3), INTENT(out) :: fog
930	REAL(r_k), DIMENSION(d1,d2,d3), INTENT(out) :: vis
931
932	! Local
933	INTEGER :: i, j, it
934
935	!!!!!!! Variables
936	! qv: mixing ratio in [kgkg-1]
937	! ta: temperature [K]
938	! pres: pressure field [Pa]
939	! fog: presence of fog (1: yes, 0: no)
940	! vis: visibility within fog [km]
941
942	fname = 'compute_fog_FRAML50'
943
944	DO i=1, d1
945	DO j=1, d2
946	DO it=1, d3
947	CALL var_fog_FRAML50(qv(i,j,it), ta(i,j,it), pres(i,j,it), fog(i,j,it), vis(i,j,it))
948	END DO
949	END DO
950	END DO
951
952	RETURN
953
954	END SUBROUTINE compute_fog_FRAML50
955
956	SUBROUTINE compute_range_faces(d1, d2, lon, lat, hgt, xdist, ydist, dist, face, dsfilt, dsnewrange, &
957	hvalrng, hgtmax, pthgtmax, derivhgt, peaks, valleys, origfaces, filtfaces, ranges, rangeshgtmax, &
958	ptrangeshgtmax)
959	! Subroutine to compute faces [uphill, valleys, downhill] of a mountain range along a given face
960
961	IMPLICIT NONE
962
963	INTEGER, INTENT(in) :: d1, d2
964	REAL(r_k), INTENT(in) :: dsfilt, dsnewrange, hvalrng
965	REAL(r_k), DIMENSION(d1,d2), INTENT(in) :: lon, lat, hgt, xdist, ydist, dist
966	CHARACTER(len=*) :: face
967	REAL(r_k), DIMENSION(d1,d2), INTENT(out) :: derivhgt, hgtmax, rangeshgtmax
968	INTEGER, DIMENSION(d1,d2), INTENT(out) :: pthgtmax, origfaces, filtfaces, peaks, &
969	valleys, ranges, ptrangeshgtmax
970	! Local
971	INTEGER :: i, j
972	INTEGER :: pthgtmax1, Npeaks, Nvalleys, Nranges
973	REAL(r_k) :: hgtmax1
974	INTEGER, DIMENSION(d1) :: ipeaks1, ivalleys1, irangeshgtmax1
975	INTEGER, DIMENSION(d2) :: ipeaks2, ivalleys2, irangeshgtmax2
976	REAL(r_k), DIMENSION(d1) :: rangeshgtmax1
977	REAL(r_k), DIMENSION(d2) :: rangeshgtmax2
978	INTEGER, DIMENSION(2,d1) :: ranges1
979	INTEGER, DIMENSION(2,d2) :: ranges2
980	INTEGER, DIMENSION(d1,d2) :: iranges
981	LOGICAL, DIMENSION(d1,d2) :: Lranges
982
983	!!!!!!! Variables
984	! lon: longitude [degrees east]
985	! lat: latitude [degrees north]
986	! hgt: topograpical height [m]
987	! face: which face (axis along which produce slices) to use to compute the faces: WE, SN
988	! dsfilt: distance to filter orography smaller scale of it [m]
989	! dsnewrange: distance to start a new mountain range [m]
990	! hvalrng: maximum height of a valley to mark change of range [m]
991	! hgtmax: maximum height of the face [m]
992	! pthgtmax: grid point of the maximum height [1]
993	! derivhgt: topograpic derivative along axis [m deg-1]
994	! peaks: peak point
995	! valleys: valley point
996	! origfaces: original faces (-1, downhill; 0: valley; 1: uphill)
997	! filtfaces: filtered faces (-1, downhill; 0: valley; 1: uphill)
998	! ranges: number of range
999	! rangeshgtmax: maximum height for each individual range [m]
1000	! ptrangeshgtmax: grid point maximum height for each individual range [1]
1001
1002	fname = 'compute_range_faces'
1003
1004	peaks = 0
1005	valleys = 0
1006	pthgtmax = 0
1007	rangeshgtmax = fillVal64
1008	IF (TRIM(face) == 'WE') THEN
1009	DO j=1, d2
1010	!PRINT ,'Lluis:', j-1, '**'
1011	CALL var_range_faces(d1, lon(:,j), lat(:,j), hgt(:,j), xdist(:,j), dsfilt, &
1012	dsnewrange, hvalrng, hgtmax1, pthgtmax1, derivhgt(:,j), Npeaks, ipeaks1, &
1013	Nvalleys, ivalleys1, origfaces(:,j), filtfaces(:,j), Nranges, ranges1, &
1014	rangeshgtmax1, irangeshgtmax1)
1015	hgtmax(:,j) = hgtmax1
1016	pthgtmax(pthgtmax1,j) = 1
1017	DO i=1, Npeaks
1018	peaks(ipeaks1(i),j) = 1
1019	END DO
1020	DO i=1, Nvalleys
1021	valleys(ivalleys1(i),j) = 1
1022	END DO
1023	DO i=1, Nranges
1024	iranges(ranges1(1,i):ranges1(2,i),j) = i
1025	rangeshgtmax(ranges1(1,i):ranges1(2,i),j) = rangeshgtmax1(i)
1026	ptrangeshgtmax(irangeshgtmax1(i),j) = 1
1027	END DO
1028	END DO
1029	ELSE IF (TRIM(face) == 'SN') THEN
1030	DO i=1, d1
1031	CALL var_range_faces(d2, lon(i,:), lat(i,:), hgt(i,:), ydist(i,:), dsfilt, &
1032	dsnewrange, hvalrng, hgtmax1, pthgtmax1, derivhgt(i,:), Npeaks, ipeaks2, &
1033	Nvalleys, ivalleys2, origfaces(i,:), filtfaces(i,:), Nranges, ranges2, &
1034	rangeshgtmax2, irangeshgtmax2)
1035	hgtmax(i,:) = hgtmax1
1036	pthgtmax(i,pthgtmax1) = 1
1037	DO j=1, Npeaks
1038	peaks(i,ipeaks2(j)) = 1
1039	END DO
1040	DO j=1, Nvalleys
1041	valleys(i,ivalleys2(j)) = 1
1042	END DO
1043	DO j=1, Nranges
1044	iranges(i,ranges2(1,j):ranges2(2,j)) = j
1045	rangeshgtmax(i,ranges2(1,j):ranges2(2,j)) = rangeshgtmax2(j)
1046	ptrangeshgtmax(i,irangeshgtmax2(j)) = 1
1047	END DO
1048	END DO
1049	ELSE
1050	PRINT *,TRIM(ErrWarnMsg('err'))
1051	PRINT *,' ' // TRIM(fname) // ": wrong face: '" // TRIM(face) // "' !!"
1052	PRINT *,' accepted ones: WE, SN'
1053	STOP
1054	END IF
1055
1056	! Homogenizing indices of the ranges
1057	CALL continguos_homogene_zones(d1, d2, iranges, Nranges, ranges)
1058	WHERE (ranges == -1)
1059	ranges = fillValI
1060	END WHERE
1061
1062	RETURN
1063
1064	END SUBROUTINE compute_range_faces
1065
1066	SUBROUTINE compute_cellbnds(dx, dy, sdx, sdy, ulon, ulat, vlon, vlat, xbnds, ybnds)
1067	! Subroutine to compute cellboundaries using wind-staggered lon, lats as intersection of their related
1068	! parallels and meridians
1069
1070	IMPLICIT NONE
1071
1072	INTEGER, INTENT(in) :: dx, dy, sdx, sdy
1073	REAL(r_k), DIMENSION(sdx, dy), INTENT(in) :: ulon, ulat
1074	REAL(r_k), DIMENSION(dx, sdy), INTENT(in) :: vlon, vlat
1075	REAL(r_k), DIMENSION(dx, dy, 4), INTENT(out) :: xbnds, ybnds
1076
1077	! Local
1078	INTEGER :: i,j,iv
1079	INTEGER :: ix,ex,iy,ey
1080	REAL(r_k) :: tmpval1, tmpval2
1081	CHARACTER(len=2), DIMENSION(4) :: Svertex
1082	INTEGER, DIMENSION(4,2,2,2) :: indices
1083	REAL(r_k), DIMENSION(2) :: ptintsct
1084	REAL(r_k), DIMENSION(2,2) :: merid, paral
1085	LOGICAL :: intsct
1086
1087	!!!!!!! Variables
1088	! dx, dy: un-staggered dimensions
1089	! sdx, sdy: staggered dimensions
1090	! ulon, ulat: x-wind staggered longitudes and latitudes
1091	! vlon, vlat: y-wind staggered longitudes and latitudes
1092	! xbnds, ybnds: x and y cell boundaries
1093
1094	fname = 'compute_cellbnds'
1095
1096	! Indices to use indices[SW/NW/NE/SE, m/p, x/y, i/e]
1097	Svertex = (/ 'SW', 'NW', 'NE', 'SE' /)
1098
1099	! SW
1100	indices(1,1,1,1) = 0
1101	indices(1,1,1,2) = 0
1102	indices(1,1,2,1) = -1
1103	indices(1,1,2,2) = 0
1104	indices(1,2,1,1) = -1
1105	indices(1,2,1,2) = 0
1106	indices(1,2,2,1) = -1
1107	indices(1,2,2,2) = -1
1108	! NW
1109	indices(2,1,1,1) = 0
1110	indices(2,1,1,2) = 0
1111	indices(2,1,2,1) = 0
1112	indices(2,1,2,2) = 1
1113	indices(2,2,1,1) = -1
1114	indices(2,2,1,2) = 0
1115	indices(2,2,2,1) = 1
1116	indices(2,2,2,2) = 1
1117	! NE
1118	indices(3,1,1,1) = 1
1119	indices(3,1,1,2) = 1
1120	indices(3,1,2,1) = 0
1121	indices(3,1,2,2) = 1
1122	indices(3,2,1,1) = 0
1123	indices(3,2,1,2) = 1
1124	indices(3,2,2,1) = 1
1125	indices(3,2,2,2) = 1
1126	! SE
1127	indices(4,1,1,1) = 1
1128	indices(4,1,1,2) = 1
1129	indices(4,1,2,1) = -1
1130	indices(4,1,2,2) = 0
1131	indices(4,2,1,1) = 0
1132	indices(4,2,1,2) = 1
1133	indices(4,2,2,1) = -1
1134	indices(4,2,2,2) = -1
1135
1136	DO i=2,dx-1
1137	DO j=1,dy
1138	DO iv=1,4
1139
1140	ix = MAX(i+indices(iv,1,1,1),1)
1141	!ex = MIN(i+indices(iv,1,1,2),dx)
1142	ex = MAX(i+indices(iv,1,1,2),1)
1143	iy = MAX(j+indices(iv,1,2,1),1)
1144	ey = MIN(j+indices(iv,1,2,2),dy)
1145
1146	merid(1,1) = ulon(ix,iy)
1147	merid(1,2) = ulat(ix,iy)
1148	merid(2,1) = ulon(ex,ey)
1149	merid(2,2) = ulat(ex,ey)
1150
1151	ix = MAX(i+indices(iv,2,1,1),1)
1152	ex = MIN(i+indices(iv,2,1,2),dx)
1153	iy = MAX(j+indices(iv,2,2,1),1)
1154	!ey = MIN(i+indices(iv,2,2,2),dy)
1155	ey = MAX(j+indices(iv,2,2,2),1)
1156	paral(1,1) = vlon(ix,iy)
1157	paral(1,2) = vlat(ix,iy)
1158	paral(2,1) = vlon(ex,ey)
1159	paral(2,2) = vlat(ex,ey)
1160
1161	CALL intersection_2Dlines(merid, paral, intsct, ptintsct)
1162	IF (.NOT.intsct) THEN
1163	msg = 'not intersection found for ' // Svertex(iv) // ' vertex'
1164	CALL ErrMsg(msg, fname, -1)
1165	END IF
1166	xbnds(i,j,iv) = ptintsct(1)
1167	ybnds(i,j,iv) = ptintsct(2)
1168
1169	END DO
1170	END DO
1171	END DO
1172
1173	! Dealing with the boundary values
1174	i = 1
1175	DO j=1,dy
1176	DO iv=1,4
1177
1178	ix = MAX(i+indices(iv,1,1,1),1)
1179	!ex = MIN(i+indices(iv,1,1,2),dx)
1180	ex = MAX(i+indices(iv,1,1,2),1)
1181	iy = MAX(j+indices(iv,1,2,1),1)
1182	ey = MIN(j+indices(iv,1,2,2),dy)
1183	merid(1,1) = ulon(ix,iy)
1184	merid(1,2) = ulat(ix,iy)
1185	merid(2,1) = ulon(ex,ey)
1186	merid(2,2) = ulat(ex,ey)
1187
1188	ix = MAX(i+indices(iv,2,1,1),1)
1189	ex = MIN(i+indices(iv,2,1,2),dx)
1190	iy = MAX(j+indices(iv,2,2,1),1)
1191	!ey = MIN(i+indices(iv,2,2,2),dy)
1192	ey = MAX(j+indices(iv,2,2,2),1)
1193	IF (iv == 1 .OR. iv == 2) THEN
1194	! Projecting values using dx from next grid point
1195	tmpval1 = vlon(2,iy)
1196	paral(2,1) = vlon(ex,ey)
1197	tmpval2 = tmpval1 - paral(2,1)
1198	paral(1,1) = paral(2,1) - tmpval2
1199	tmpval1 = vlat(2,iy)
1200	paral(2,2) = vlat(ex,ey)
1201	tmpval2 = tmpval1 - paral(2,2)
1202	paral(1,2) = paral(2,2) - tmpval2
1203	ELSE
1204	paral(1,1) = vlon(ix,iy)
1205	paral(1,2) = vlat(ix,iy)
1206	paral(2,1) = vlon(ex,ey)
1207	paral(2,2) = vlat(ex,ey)
1208	END IF
1209
1210	CALL intersection_2Dlines(merid, paral, intsct, ptintsct)
1211	IF (.NOT.intsct) THEN
1212	msg = 'not intersection found for ' // Svertex(iv) // ' vertex'
1213	CALL ErrMsg(msg, fname, -1)
1214	END IF
1215	xbnds(i,j,iv) = ptintsct(1)
1216	ybnds(i,j,iv) = ptintsct(2)
1217
1218	END DO
1219	END DO
1220
1221	i = dx
1222	DO j=1,dy
1223	DO iv=1,4
1224
1225	ix = MAX(i+indices(iv,1,1,1),1)
1226	!ex = MIN(i+indices(iv,1,1,2),dx)
1227	ex = MAX(i+indices(iv,1,1,2),1)
1228	iy = MAX(j+indices(iv,1,2,1),1)
1229	ey = MIN(j+indices(iv,1,2,2),dy)
1230	merid(1,1) = ulon(ix,iy)
1231	merid(1,2) = ulat(ix,iy)
1232	merid(2,1) = ulon(ex,ey)
1233	merid(2,2) = ulat(ex,ey)
1234
1235	ix = MAX(i+indices(iv,2,1,1),1)
1236	ex = MIN(i+indices(iv,2,1,2),dx)
1237	iy = MAX(j+indices(iv,2,2,1),1)
1238	!ey = MIN(i+indices(iv,2,2,2),dy)
1239	ey = MAX(j+indices(iv,2,2,2),1)
1240	IF (iv == 3 .OR. iv == 4) THEN
1241	! Projecting values using dx from previous grid point
1242	tmpval1 = vlon(dx-1,iy)
1243	paral(2,1) = vlon(ex,ey)
1244	tmpval2 = tmpval1 - paral(2,1)
1245	paral(1,1) = paral(2,1) - tmpval2
1246	tmpval1 = vlat(dx-1,iy)
1247	paral(2,2) = vlat(ex,ey)
1248	tmpval2 = tmpval1 - paral(2,2)
1249	paral(1,2) = paral(2,2) - tmpval2
1250	ELSE
1251	paral(1,1) = vlon(ix,iy)
1252	paral(1,2) = vlat(ix,iy)
1253	paral(2,1) = vlon(ex,ey)
1254	paral(2,2) = vlat(ex,ey)
1255	END IF
1256
1257	CALL intersection_2Dlines(merid, paral, intsct, ptintsct)
1258	IF (.NOT.intsct) THEN
1259	msg = 'not intersection found for ' // Svertex(iv) // ' vertex'
1260	CALL ErrMsg(msg, fname, -1)
1261	END IF
1262	xbnds(i,j,iv) = ptintsct(1)
1263	ybnds(i,j,iv) = ptintsct(2)
1264
1265	END DO
1266	END DO
1267
1268	END SUBROUTINE compute_cellbnds
1269
1270	SUBROUTINE compute_cellbndsreg(dx, dy, lon, lat, xbnds, ybnds)
1271	! Subroutine to compute cellboundaries using lon, lat from a reglar lon/lat projection as intersection
1272	! of their related parallels and meridians
1273
1274	IMPLICIT NONE
1275
1276	INTEGER, INTENT(in) :: dx, dy
1277	REAL(r_k), DIMENSION(dx, dy), INTENT(in) :: lon, lat
1278	REAL(r_k), DIMENSION(dx, dy, 4), INTENT(out) :: xbnds, ybnds
1279
1280	! Local
1281	INTEGER :: i,j,iv
1282	INTEGER :: ix,ex,iy,ey
1283	CHARACTER(len=2), DIMENSION(4) :: Svertex
1284	INTEGER, DIMENSION(4,2,2) :: indices
1285
1286	!!!!!!! Variables
1287	! dx, dy: un-staggered dimensions
1288	! lon, lat: longitudes and latitudes
1289	! xbnds, ybnds: x and y cell boundaries
1290
1291	fname = 'compute_cellbndsreg'
1292
1293	! Indices to use indices[SW/NW/NE/SE, m/p, x/y, i/e]
1294	Svertex = (/ 'SW', 'NW', 'NE', 'SE' /)
1295
1296	! SW
1297	indices(1,1,1) = -1
1298	indices(1,1,2) = 0
1299	indices(1,2,1) = -1
1300	indices(1,2,2) = 0
1301	! NW
1302	indices(2,1,1) = -1
1303	indices(2,1,2) = 0
1304	indices(2,2,1) = 0
1305	indices(2,2,2) = 1
1306	! NE
1307	indices(3,1,1) = 0
1308	indices(3,1,2) = 1
1309	indices(3,2,1) = 0
1310	indices(3,2,2) = 1
1311	! SE
1312	indices(4,1,1) = 0
1313	indices(4,1,2) = 1
1314	indices(4,2,1) = -1
1315	indices(4,2,2) = 0
1316
1317	DO i=1,dx
1318	DO j=1,dy
1319	DO iv=1,4
1320
1321	ix = MAX(i+indices(iv,1,1),1)
1322	ix = MIN(ix,dx)
1323	ex = MAX(i+indices(iv,1,2),1)
1324	ex = MIN(ex,dx)
1325	iy = MAX(j+indices(iv,2,1),1)
1326	iy = MIN(iy,dy)
1327	ey = MAX(j+indices(iv,2,2),1)
1328	ey = MIN(ey,dy)
1329
1330	xbnds(i,j,iv) = 0.5*(lon(ix,iy) + lon(ex,ey))
1331	ybnds(i,j,iv) = 0.5*(lat(ix,iy) + lat(ex,ey))
1332
1333	END DO
1334	END DO
1335	END DO
1336
1337	END SUBROUTINE
1338
1339	SUBROUTINE compute_Koeppen_Geiger_climates(dx, dy, dt, pr, tas, climatesI, climatesS, climlegend)
1340	! Subroutine to compute the Koeppen-Geiger climates after:
1341	! Kottek et al., Meteorologische Zeitschrift, Vol. 15, No. 3, 259-263
1342
1343	IMPLICIT NONE
1344
1345	INTEGER, INTENT(in) :: dx, dy, dt
1346	REAL(r_k), DIMENSION(dx, dy, dt), INTENT(in) :: pr, tas
1347	INTEGER, DIMENSION(dy, dt), INTENT(out) :: climatesI
1348	CHARACTER(LEN=3), DIMENSION(dy, dt), INTENT(out) :: climatesS
1349	CHARACTER(LEN=1000), INTENT(out) :: climlegend
1350
1351	! Local
1352	INTEGER :: i,j
1353
1354
1355
1356
1357	END SUBROUTINE compute_Koeppen_Geiger_climates
1358
1359	END MODULE module_ForDiagnostics

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

Download in other formats:

Original Format