Context Navigation

module_ForDiagnostics.f90 @ 2734

Last change on this file since 2734 was 2675, checked in by lfita, 6 years ago

Adding:

'frontogenesis': Calculation of the frontogenesis diagnostics

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

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

Download in other formats:

Original Format