Context Navigation

module_ForDiagnostics.f90 @ 2424

Last change on this file since 2424 was 2387, checked in by lfita, 6 years ago

Adding:

`tws': Wet Bulb temperature after Stull, 2011

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

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

Download in other formats:

Original Format