Context Navigation

module_ForDiagnostics.f90 @ 2329

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

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

Download in other formats:

Original Format