Context Navigation

module_ForDiagnostics.f90 @ 2805

Last change on this file since 2805 was 2765, checked in by lfita, 6 years ago

Adding:

`1st order 2D horizontal gradient'

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

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

Download in other formats:

Original Format