Context Navigation

module_ForDiagnostics.f90 @ 1766

Last change on this file since 1766 was 1762, checked in by lfita, 8 years ago

Adding:

`mrso': total soil moisture
`slw': total liquid water

File size: 50.5 KB

Line
1	!! Fortran version of different diagnostics
2	! L. Fita. LMD May 2016
3	! gfortran module_generic.o module_ForDiagnosticsVars.o -c module_ForDiagnostics.F90
4	!
5	! f2py -m module_ForDiagnostics --f90exec=/usr/bin/gfortran-4.7 -c module_generic.F90 module_ForDiagnosticsVars.F90 module_ForDiagnostics.F90
6
7	MODULE module_ForDiagnostics
8
9	USE module_definitions
10	USE module_generic
11	USE module_ForDiagnosticsVars
12
13	CONTAINS
14
15	!!!!!!! Calculations
16	! compute_cllmh4D3: Computation of low, medium and high cloudiness from a 4D CLDFRA and pressure being
17	! 3rd dimension the z-dim
18	! compute_cllmh3D3: Computation of low, medium and high cloudiness from a 3D CLDFRA and pressure being
19	! 3rd dimension the z-dim
20	! compute_cllmh: Computation of low, medium and high cloudiness
21	! compute_clt4D3: Computation of total cloudiness from a 4D CLDFRA being 3rd dimension the z-dim
22	! compute_clt3D3: Computation of total cloudiness from a 3D CLDFRA being 3rd dimension the z-dim
23	! compute_clt: Computation of total cloudiness
24	! compute_psl_ptarget4d2: Compute sea level pressure using a target pressure. Similar to the Benjamin
25	! and Miller (1990). Method found in p_interp.F90
26	! compute_tv4d: 4D calculation of virtual temperaure
27	! VirtualTemp1D: Function for 1D calculation of virtual temperaure
28
29	!!!
30	! Calculations
31	!!!
32
33	SUBROUTINE compute_cllmh4D2(cldfra4D, pres4D, cllmh4D2, d1, d2, d3, d4)
34	! Subroutine to compute the low, medium and high cloudiness following 'newmicro.F90' from LMDZ from a 4D CLDFRA and pressure
35	! where zdim is the 2nd dimension (thus, cldfra4D(d1,d2,d3,d4) --> cllmh(3,d1,d3,d4) 1: low, 2: medium, 3: high
36	! It should be properly done via an 'INTERFACE', but...
37
38	IMPLICIT NONE
39
40	INTEGER, INTENT(in) :: d1, d2, d3, d4
41	REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in) :: cldfra4D, pres4D
42	REAL(r_k), DIMENSION(3,d1,d3,d4), INTENT(out) :: cllmh4D2
43
44	! Local
45	INTEGER :: i,j,k, zdim, Ndim
46
47	!!!!!!! Variables
48	! cldfra4D: 4D cloud fraction values [1]
49	! pres4D: 4D pressure values [Pa]
50	! Ndim: number of dimensions of the input data
51	! d[1-4]: dimensions of 'cldfra4D'
52	! zdim: number of the vertical-dimension within the matrix
53	! cltlmh4D2: low, medium, high cloudiness for the 4D cldfra and d2 being 'zdim'
54
55	fname = 'compute_cllmh4D2'
56	zdim = 2
57	Ndim = 4
58
59	DO i=1, d1
60	DO j=1, d3
61	DO k=1, d4
62	cllmh4D2(:,i,j,k) = var_cllmh(cldfra4D(i,:,j,k), pres4D(i,:,j,k), d2)
63	END DO
64	END DO
65	END DO
66
67	RETURN
68
69	END SUBROUTINE compute_cllmh4D2
70
71	SUBROUTINE compute_cllmh3D1(cldfra3D, pres3D, cllmh3D1, d1, d2, d3)
72	! Subroutine to compute the low, medium and high cloudiness following 'newmicro.F90' from LMDZ from a 3D CLDFRA and pressure
73	! where zdim is the 1st dimension (thus, cldfra3D(d1,d2,d3) --> cllmh(3,d2,d3) 1: low, 2: medium, 3: high
74	! It should be properly done via an 'INTERFACE', but...
75
76	IMPLICIT NONE
77
78	INTEGER, INTENT(in) :: d1, d2, d3
79	REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in) :: cldfra3D, pres3D
80	REAL(r_k), DIMENSION(3,d2,d3), INTENT(out) :: cllmh3D1
81
82	! Local
83	INTEGER :: i,j,k, zdim, Ndim
84
85	!!!!!!! Variables
86	! cldfra3D: 3D cloud fraction values [1]
87	! pres3D: 3D pressure values [Pa]
88	! Ndim: number of dimensions of the input data
89	! d[1-3]: dimensions of 'cldfra3D'
90	! zdim: number of the vertical-dimension within the matrix
91	! cltlmh3D1: low, medium, high cloudiness for the 3D cldfra and d1 being 'zdim'
92
93	fname = 'compute_cllmh3D1'
94	zdim = 1
95	Ndim = 3
96
97	DO i=1, d1
98	DO j=1, d2
99	cllmh3D1(:,i,j) = var_cllmh(cldfra3D(:,i,j), pres3D(:,i,j), d1)
100	END DO
101	END DO
102
103	RETURN
104
105	END SUBROUTINE compute_cllmh3D1
106
107	SUBROUTINE compute_cllmh(cldfra1D, cldfra2D, cldfra3D, cldfra4D, pres1D, pres2D, pres3D, pres4D, &
108	Ndim, zdim, cllmh1D, cllmh2D1, cllmh2D2, cllmh3D1, cllmh3D2, cllmh3D3, cllmh4D1, cllmh4D2, &
109	cllmh4D3, cllmh4D4, d1, d2, d3, d4)
110	! Subroutine to compute the low, medium and high cloudiness following 'newmicro.F90' from LMDZ
111
112	IMPLICIT NONE
113
114	INTEGER, INTENT(in) :: Ndim, d1, d2, d3, d4, zdim
115	REAL(r_k), DIMENSION(d1), OPTIONAL, INTENT(in) :: cldfra1D, pres1D
116	REAL(r_k), DIMENSION(d1,d2), OPTIONAL, INTENT(in) :: cldfra2D, pres2D
117	REAL(r_k), DIMENSION(d1,d2,d3), OPTIONAL, INTENT(in) :: cldfra3D, pres3D
118	REAL(r_k), DIMENSION(d1,d2,d3,d4), OPTIONAL, &
119	INTENT(in) :: cldfra4D, pres4D
120	REAL(r_k), DIMENSION(3), OPTIONAL, INTENT(out) :: cllmh1D
121	REAL(r_k), DIMENSION(d1,3), OPTIONAL, INTENT(out) :: cllmh2D1
122	REAL(r_k), DIMENSION(d2,3), OPTIONAL, INTENT(out) :: cllmh2D2
123	REAL(r_k), DIMENSION(d2,d3,3), OPTIONAL, INTENT(out) :: cllmh3D1
124	REAL(r_k), DIMENSION(d1,d3,3), OPTIONAL, INTENT(out) :: cllmh3D2
125	REAL(r_k), DIMENSION(d1,d2,3), OPTIONAL, INTENT(out) :: cllmh3D3
126	REAL(r_k), DIMENSION(d2,d3,d4,3), OPTIONAL, &
127	INTENT(out) :: cllmh4D1
128	REAL(r_k), DIMENSION(d1,d3,d4,3), OPTIONAL, &
129	INTENT(out) :: cllmh4D2
130	REAL(r_k), DIMENSION(d1,d2,d4,3), OPTIONAL, &
131	INTENT(out) :: cllmh4D3
132	REAL(r_k), DIMENSION(d1,d2,d3,3), OPTIONAL, &
133	INTENT(out) :: cllmh4D4
134
135	! Local
136	INTEGER :: i,j,k
137
138	!!!!!!! Variables
139	! cldfra[1-4]D: cloud fraction values [1]
140	! pres[1-4]D: pressure values [Pa]
141	! Ndim: number of dimensions of the input data
142	! d[1-4]: dimensions of 'cldfra'
143	! zdim: number of the vertical-dimension within the matrix
144	! cllmh1D: low, medium and high cloudiness for the 1D cldfra
145	! cllmh2D1: low, medium and high cloudiness for the 2D cldfra and d1 being 'zdim'
146	! cllmh2D2: low, medium and high cloudiness for the 2D cldfra and d2 being 'zdim'
147	! cllmh3D1: low, medium and high cloudiness for the 3D cldfra and d1 being 'zdim'
148	! cllmh3D2: low, medium and high cloudiness for the 3D cldfra and d2 being 'zdim'
149	! cllmh3D3: low, medium and high cloudiness for the 3D cldfra and d3 being 'zdim'
150	! cllmh4D1: low, medium and high cloudiness for the 4D cldfra and d1 being 'zdim'
151	! cllmh4D2: low, medium and high cloudiness for the 4D cldfra and d2 being 'zdim'
152	! cllmh4D3: low, medium and high cloudiness for the 4D cldfra and d3 being 'zdim'
153	! cllmh4D4: low, medium and high cloudiness for the 4D cldfra and d4 being 'zdim'
154
155	fname = 'compute_cllmh'
156
157	SELECT CASE (Ndim)
158	CASE (1)
159	cllmh1D = var_cllmh(cldfra1D, pres1D, d1)
160	CASE (2)
161	IF (zdim == 1) THEN
162	DO i=1, d2
163	cllmh2D1(i,:) = var_cllmh(cldfra2D(:,i), pres2D(:,i), d1)
164	END DO
165	ELSE IF (zdim == 2) THEN
166	DO i=1, d1
167	cllmh2D2(i,:) = var_cllmh(cldfra2D(:,i), pres2D(i,:), d2)
168	END DO
169	ELSE
170	PRINT *,TRIM(ErrWarnMsg('err'))
171	PRINT *,' ' // TRIM(fname) // ': wrong zdim:', zdim,' for Ndim=', Ndim, ' !!'
172	PRINT *,' accepted values: 1,2'
173	STOP
174	END IF
175	CASE (3)
176	IF (zdim == 1) THEN
177	DO i=1, d2
178	DO j=1, d3
179	cllmh3D1(i,j,:) = var_cllmh(cldfra3D(:,i,j), pres3D(:,i,j), d1)
180	END DO
181	END DO
182	ELSE IF (zdim == 2) THEN
183	DO i=1, d1
184	DO j=1, d3
185	cllmh3D2(i,j,:) = var_cllmh(cldfra3D(i,:,j), pres3D(i,:,j), d2)
186	END DO
187	END DO
188	ELSE IF (zdim == 3) THEN
189	DO i=1, d1
190	DO j=1, d2
191	cllmh3D3(i,j,:) = var_cllmh(cldfra3D(i,j,:), pres3D(i,j,:), d3)
192	END DO
193	END DO
194	ELSE
195	PRINT *,TRIM(ErrWarnMsg('err'))
196	PRINT *,' ' // TRIM(fname) // ': wrong zdim:', zdim,' for Ndim=', Ndim, ' !!'
197	PRINT *,' accepted values: 1,2,3'
198	STOP
199	END IF
200	CASE (4)
201	IF (zdim == 1) THEN
202	DO i=1, d2
203	DO j=1, d3
204	DO k=1, d4
205	cllmh4D1(i,j,k,:) = var_cllmh(cldfra4D(:,i,j,k), pres4D(:,i,j,k), d1)
206	END DO
207	END DO
208	END DO
209	ELSE IF (zdim == 2) THEN
210	DO i=1, d1
211	DO j=1, d3
212	DO k=1, d4
213	cllmh4D2(i,j,k,:) = var_cllmh(cldfra4D(i,:,j,k), pres4D(i,:,j,k), d2)
214	END DO
215	END DO
216	END DO
217	ELSE IF (zdim == 3) THEN
218	DO i=1, d2
219	DO j=1, d3
220	DO k=1, d4
221	cllmh4D3(i,j,k,:) = var_cllmh(cldfra4D(i,j,:,k), pres4D(i,j,:,k), d3)
222	END DO
223	END DO
224	END DO
225	ELSE IF (zdim == 4) THEN
226	DO i=1, d1
227	DO j=1, d2
228	DO k=1, d3
229	cllmh4D4(i,j,k,:) = var_cllmh(cldfra4D(i,j,k,:), pres4D(i,j,k,:), d4)
230	END DO
231	END DO
232	END DO
233	ELSE
234	PRINT *,TRIM(ErrWarnMsg('err'))
235	PRINT *,' ' // TRIM(fname) // ': wrong zdim:', zdim,' for Ndim=', Ndim, ' !!'
236	PRINT *,' accepted values: 1,2,3,4'
237	STOP
238	END IF
239	CASE DEFAULT
240	PRINT *,TRIM(ErrWarnMsg('err'))
241	PRINT *,' ' // TRIM(fname) // ': Ndim:', Ndim,' not ready !!'
242	STOP
243	END SELECT
244
245	RETURN
246
247	END SUBROUTINE compute_cllmh
248
249	SUBROUTINE compute_clt4D2(cldfra4D, clt4D2, d1, d2, d3, d4)
250	! Subroutine to compute the total cloudiness following 'newmicro.F90' from LMDZ from a 4D CLDFRA
251	! where zdim is the 2nd dimension (thus, cldfra4D(d1,d2,d3,d4) --> clt(d1,d3,d4)
252	! It should be properly done via an 'INTERFACE', but...
253
254	IMPLICIT NONE
255
256	INTEGER, INTENT(in) :: d1, d2, d3, d4
257	REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in) :: cldfra4D
258	REAL(r_k), DIMENSION(d1,d3,d4), INTENT(out) :: clt4D2
259
260	! Local
261	INTEGER :: i,j,k, zdim, Ndim
262
263	!!!!!!! Variables
264	! cldfra4D: 4D cloud fraction values [1]
265	! Ndim: number of dimensions of the input data
266	! d[1-4]: dimensions of 'cldfra4D'
267	! zdim: number of the vertical-dimension within the matrix
268	! clt4D2: total cloudiness for the 4D cldfra and d2 being 'zdim'
269
270	fname = 'compute_clt4D2'
271	zdim = 2
272	Ndim = 4
273
274	DO i=1, d1
275	DO j=1, d3
276	DO k=1, d4
277	clt4D2(i,j,k) = var_clt(cldfra4D(i,:,j,k), d2)
278	END DO
279	END DO
280	END DO
281
282	RETURN
283
284	END SUBROUTINE compute_clt4D2
285
286	SUBROUTINE compute_clt3D1(cldfra3D, clt3D1, d1, d2, d3)
287	! Subroutine to compute the total cloudiness following 'newmicro.F90' from LMDZ from a 3D CLDFRA
288	! where zdim is the 1st dimension (thus, cldfra4D(d1,d2,d3) --> clt(d2,d3)
289	! It should be properly done via an 'INTERFACE', but...
290
291	IMPLICIT NONE
292
293	INTEGER, INTENT(in) :: d1, d2, d3
294	REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in) :: cldfra3D
295	REAL(r_k), DIMENSION(d2,d3), INTENT(out) :: clt3D1
296
297	! Local
298	INTEGER :: i,j,k, zdim, Ndim
299
300	!!!!!!! Variables
301	! cldfra3D: 3D cloud fraction values [1]
302	! Ndim: number of dimensions of the input data
303	! d[1-3]: dimensions of 'cldfra3D'
304	! zdim: number of the vertical-dimension within the matrix
305	! clt3D1: total cloudiness for the 3D cldfra and d1 being 'zdim'
306
307	fname = 'compute_clt3D1'
308	zdim = 1
309	Ndim = 3
310
311	DO i=1, d2
312	DO j=1, d3
313	clt3D1(i,j) = var_clt(cldfra3D(:,i,j), d1)
314	END DO
315	END DO
316
317	RETURN
318
319	END SUBROUTINE compute_clt3D1
320
321	SUBROUTINE compute_clt(cldfra1D, cldfra2D, cldfra3D, cldfra4D, Ndim, zdim, clt1D, clt2D1, clt2D2, &
322	clt3D1, clt3D2, clt3D3, clt4D1, clt4D2, clt4D3, clt4D4, d1, d2, d3, d4)
323	! Subroutine to compute the total cloudiness following 'newmicro.F90' from LMDZ
324
325	IMPLICIT NONE
326
327	INTEGER, INTENT(in) :: Ndim, d1, d2, d3, d4, zdim
328	REAL(r_k), DIMENSION(d1), OPTIONAL, INTENT(in) :: cldfra1D
329	REAL(r_k), DIMENSION(d1,d2), OPTIONAL, INTENT(in) :: cldfra2D
330	REAL(r_k), DIMENSION(d1,d2,d3), OPTIONAL, INTENT(in) :: cldfra3D
331	REAL(r_k), DIMENSION(d1,d2,d3,d4), OPTIONAL, &
332	INTENT(in) :: cldfra4D
333	REAL(r_k), OPTIONAL, INTENT(out) :: clt1D
334	REAL(r_k), DIMENSION(d1), OPTIONAL, INTENT(out) :: clt2D1
335	REAL(r_k), DIMENSION(d2), OPTIONAL, INTENT(out) :: clt2D2
336	REAL(r_k), DIMENSION(d2,d3), OPTIONAL, INTENT(out) :: clt3D1
337	REAL(r_k), DIMENSION(d1,d3), OPTIONAL, INTENT(out) :: clt3D2
338	REAL(r_k), DIMENSION(d1,d2), OPTIONAL, INTENT(out) :: clt3D3
339	REAL(r_k), DIMENSION(d2,d3,d4), OPTIONAL,INTENT(out) :: clt4D1
340	REAL(r_k), DIMENSION(d1,d3,d4), OPTIONAL,INTENT(out) :: clt4D2
341	REAL(r_k), DIMENSION(d1,d2,d4), OPTIONAL,INTENT(out) :: clt4D3
342	REAL(r_k), DIMENSION(d1,d2,d3), OPTIONAL,INTENT(out) :: clt4D4
343
344	! Local
345	INTEGER :: i,j,k
346
347	!!!!!!! Variables
348	! cldfra[1-4]D: cloud fraction values [1]
349	! Ndim: number of dimensions of the input data
350	! d[1-4]: dimensions of 'cldfra'
351	! zdim: number of the vertical-dimension within the matrix
352	! clt1D: total cloudiness for the 1D cldfra
353	! clt2D1: total cloudiness for the 2D cldfra and d1 being 'zdim'
354	! clt2D2: total cloudiness for the 2D cldfra and d2 being 'zdim'
355	! clt3D1: total cloudiness for the 3D cldfra and d1 being 'zdim'
356	! clt3D2: total cloudiness for the 3D cldfra and d2 being 'zdim'
357	! clt3D3: total cloudiness for the 3D cldfra and d3 being 'zdim'
358	! clt4D1: total cloudiness for the 4D cldfra and d1 being 'zdim'
359	! clt4D2: total cloudiness for the 4D cldfra and d2 being 'zdim'
360	! clt4D3: total cloudiness for the 4D cldfra and d3 being 'zdim'
361	! clt4D4: total cloudiness for the 4D cldfra and d4 being 'zdim'
362
363	fname = 'compute_clt'
364
365	SELECT CASE (Ndim)
366	CASE (1)
367	clt1D = var_clt(cldfra1D, d1)
368	CASE (2)
369	IF (zdim == 1) THEN
370	DO i=1, d2
371	clt2D1(i) = var_clt(cldfra2D(:,i), d1)
372	END DO
373	ELSE IF (zdim == 2) THEN
374	DO i=1, d1
375	clt2D2(i) = var_clt(cldfra2D(:,i), d2)
376	END DO
377	ELSE
378	PRINT *,TRIM(ErrWarnMsg('err'))
379	PRINT *,' ' // TRIM(fname) // ': wrong zdim:', zdim,' for Ndim=', Ndim, ' !!'
380	PRINT *,' accepted values: 1,2'
381	STOP
382	END IF
383	CASE (3)
384	IF (zdim == 1) THEN
385	DO i=1, d2
386	DO j=1, d3
387	clt3D1(i,j) = var_clt(cldfra3D(:,i,j), d1)
388	END DO
389	END DO
390	ELSE IF (zdim == 2) THEN
391	DO i=1, d1
392	DO j=1, d3
393	clt3D2(i,j) = var_clt(cldfra3D(i,:,j), d2)
394	END DO
395	END DO
396	ELSE IF (zdim == 3) THEN
397	DO i=1, d1
398	DO j=1, d2
399	clt3D3(i,j) = var_clt(cldfra3D(i,j,:), d3)
400	END DO
401	END DO
402	ELSE
403	PRINT *,TRIM(ErrWarnMsg('err'))
404	PRINT *,' ' // TRIM(fname) // ': wrong zdim:', zdim,' for Ndim=', Ndim, ' !!'
405	PRINT *,' accepted values: 1,2,3'
406	STOP
407	END IF
408	CASE (4)
409	IF (zdim == 1) THEN
410	DO i=1, d2
411	DO j=1, d3
412	DO k=1, d4
413	clt4D1(i,j,k) = var_clt(cldfra4D(:,i,j,k), d1)
414	END DO
415	END DO
416	END DO
417	ELSE IF (zdim == 2) THEN
418	DO i=1, d1
419	DO j=1, d3
420	DO k=1, d4
421	clt4D2(i,j,k) = var_clt(cldfra4D(i,:,j,k), d2)
422	END DO
423	END DO
424	END DO
425	ELSE IF (zdim == 3) THEN
426	DO i=1, d2
427	DO j=1, d3
428	DO k=1, d4
429	clt4D3(i,j,k) = var_clt(cldfra4D(i,j,:,k), d3)
430	END DO
431	END DO
432	END DO
433	ELSE IF (zdim == 4) THEN
434	DO i=1, d1
435	DO j=1, d2
436	DO k=1, d3
437	clt4D4(i,j,k) = var_clt(cldfra4D(i,j,k,:), d4)
438	END DO
439	END DO
440	END DO
441	ELSE
442	PRINT *,TRIM(ErrWarnMsg('err'))
443	PRINT *,' ' // TRIM(fname) // ': wrong zdim:', zdim,' for Ndim=', Ndim, ' !!'
444	PRINT *,' accepted values: 1,2,3,4'
445	STOP
446	END IF
447	CASE DEFAULT
448	PRINT *,TRIM(ErrWarnMsg('err'))
449	PRINT *,' ' // TRIM(fname) // ': Ndim:', Ndim,' not ready !!'
450	STOP
451	END SELECT
452
453	RETURN
454
455	END SUBROUTINE compute_clt
456
457	SUBROUTINE compute_psl_ptarget4d2(press, ps, hgt, ta, qv, ptarget, psl, d1, d2, d3, d4)
458	! Subroutine to compute sea level pressure using a target pressure. Similar to the Benjamin
459	! and Miller (1990). Method found in p_interp.F90
460
461	IMPLICIT NONE
462
463	INTEGER, INTENT(in) :: d1, d2, d3, d4
464	REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in) :: press, ta, qv
465	REAL(r_k), DIMENSION(d1,d2,d4), INTENT(in) :: ps
466	REAL(r_k), DIMENSION(d1,d2), INTENT(in) :: hgt
467	REAL(r_k), INTENT(in) :: ptarget
468	REAL(r_k), DIMENSION(d1,d2,d4), INTENT(out) :: psl
469
470	! Local
471	INTEGER :: i, j, it
472	INTEGER :: kin
473	INTEGER :: kupper
474	REAL(r_k) :: dpmin, dp, tbotextrap, &
475	tvbotextrap, virtual
476	! Exponential related to standard atmosphere lapse rate r_d*gammav/g
477	REAL(r_k), PARAMETER :: expon=r_d*gammav/grav
478
479	!!!!!!! Variables
480	! press: Atmospheric pressure [Pa]
481	! ps: surface pressure [Pa]
482	! hgt: surface height
483	! ta: temperature [K]
484	! qv: water vapor mixing ratio
485	! dz: number of vertical levels
486	! psl: sea-level pressure
487
488	fname = 'compute_psl_ptarget4d2'
489
490	! Minimal distance between pressures [Pa]
491	dpmin=1.e4
492	psl=0.
493
494	DO i=1,d1
495	DO j=1,d2
496	IF (hgt(i,j) /= 0.) THEN
497	DO it=1,d4
498
499	! target pressure to be used for the extrapolation [Pa] (defined in namelist.input)
500	! ptarget = 70000. default value
501
502	! We are below both the ground and the lowest data level.
503
504	! First, find the model level that is closest to a "target" pressure
505	! level, where the "target" pressure is delta-p less that the local
506	! value of a horizontally smoothed surface pressure field. We use
507	! delta-p = 150 hPa here. A standard lapse rate temperature profile
508	! passing through the temperature at this model level will be used
509	! to define the temperature profile below ground. This is similar
510	! to the Benjamin and Miller (1990) method, using
511	! 700 hPa everywhere for the "target" pressure.
512
513	kupper = 0
514	loop_kIN: DO kin=d3,1,-1
515	kupper = kin
516	dp=abs( press(i,j,kin,it) - ptarget )
517	IF (dp .GT. dpmin) EXIT loop_kIN
518	dpmin=min(dpmin,dp)
519	ENDDO loop_kIN
520
521	tbotextrap=ta(i,j,kupper,it)(ps(i,j,it)/ptarget)*expon
522	tvbotextrap=virtualTemp1D(tbotextrap,qv(i,j,kupper,it))
523
524	psl(i,j,it) = ps(i,j,it)((tvbotextrap+gammavhgt(i,j))/tvbotextrap)**(1/expon)
525	END DO
526	ELSE
527	psl(i,j,:) = ps(i,j,:)
528	END IF
529	END DO
530	END DO
531
532	RETURN
533
534	END SUBROUTINE compute_psl_ptarget4d2
535
536	SUBROUTINE compute_massvertint1D(var, mutot, dz, deta, integral)
537	! Subroutine to vertically integrate a 1D variable in eta vertical coordinates
538
539	IMPLICIT NONE
540
541	INTEGER, INTENT(in) :: dz
542	REAL(r_k), INTENT(in) :: mutot
543	REAL(r_k), DIMENSION(dz), INTENT(in) :: var, deta
544	REAL(r_k), INTENT(out) :: integral
545
546	! Local
547	INTEGER :: k
548
549	!!!!!!! Variables
550	! var: vertical variable to integrate (assuming kgkg-1)
551	! mutot: total dry-air mass in column
552	! dz: vertical dimension
553	! deta: eta-levels difference between full eta-layers
554
555	fname = 'compute_massvertint1D'
556
557	! integral=0.
558	! DO k=1,dz
559	! integral = integral + var(k)*deta(k)
560	! END DO
561	integral = SUM(var*deta)
562
563	integral=integral*mutot/g
564
565	RETURN
566
567	END SUBROUTINE compute_massvertint1D
568
569	SUBROUTINE compute_zint4D(var4D, dlev, zweight, d1, d2, d3, d4, int3D)
570	! Subroutine to vertically integrate a 4D variable in any vertical coordinates
571
572	IMPLICIT NONE
573
574	INTEGER, INTENT(in) :: d1,d2,d3,d4
575	REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in) :: var4D, dlev, zweight
576	REAL(r_k), DIMENSION(d1,d2,d4), INTENT(out) :: int3D
577
578	! Local
579	INTEGER :: i,j,l
580
581	!!!!!!! Variables
582	! var4D: vertical variable to integrate
583	! dlev: height of layers
584	! zweight: weight for each level to be applied (=1. for no effect)
585
586	fname = 'compute_zint4D'
587
588	DO i=1,d1
589	DO j=1,d2
590	DO l=1,d4
591	CALL compute_vertint1D(var4D(i,j,:,l),d3, dlev(i,j,:,l), zweight(i,j,:,l), &
592	int3D(i,j,l))
593	END DO
594	END DO
595	END DO
596
597	RETURN
598
599	END SUBROUTINE compute_zint4D
600
601	SUBROUTINE compute_vertint1D(var, dz, deta, zweight, integral)
602	! Subroutine to vertically integrate a 1D variable in any vertical coordinates
603
604	IMPLICIT NONE
605
606	INTEGER, INTENT(in) :: dz
607	REAL(r_k), DIMENSION(dz), INTENT(in) :: var, deta, zweight
608	REAL(r_k), INTENT(out) :: integral
609
610	! Local
611	INTEGER :: k
612
613	!!!!!!! Variables
614	! var: vertical variable to integrate
615	! dz: vertical dimension
616	! deta: eta-levels difference between layers
617	! zweight: weight for each level to be applied (=1. for no effect)
618
619	fname = 'compute_vertint1D'
620
621	! integral=0.
622	! DO k=1,dz
623	! integral = integral + var(k)*deta(k)
624	! END DO
625	integral = SUM(vardetazweight)
626
627	RETURN
628
629	END SUBROUTINE compute_vertint1D
630
631	SUBROUTINE compute_tv4d(ta,qv,tv,d1,d2,d3,d4)
632	! 4D calculation of virtual temperaure
633
634	IMPLICIT NONE
635
636	INTEGER, INTENT(in) :: d1, d2, d3, d4
637	REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in) :: ta, qv
638	REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(out) :: tv
639
640	! Variables
641	! ta: temperature [K]
642	! qv: mixing ratio [kgkg-1]
643	! tv: virtual temperature
644
645	tv = ta*(oneRK+(qv/epsilonv))/(oneRK+qv)
646
647	END SUBROUTINE compute_tv4d
648
649	FUNCTION VirtualTemp1D (ta,qv) result (tv)
650	! 1D calculation of virtual temperaure
651
652	IMPLICIT NONE
653
654	REAL(r_k), INTENT(in) :: ta, qv
655	REAL(r_k) :: tv
656
657	! Variables
658	! ta: temperature [K]
659	! qv: mixing ratio [kgkg-1]
660
661	tv = ta*(oneRK+(qv/epsilonv))/(oneRK+qv)
662
663	END FUNCTION VirtualTemp1D
664
665	! ---- BEGIN modified from module_diag_afwa.F ---- !
666
667	!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
668	!~
669	!~ Name:
670	!~ Theta
671	!~
672	!~ Description:
673	!~ This function calculates potential temperature as defined by
674	!~ Poisson's equation, given temperature and pressure ( hPa ).
675	!~
676	!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
677	FUNCTION Theta ( t, p )
678
679	IMPLICIT NONE
680
681	!~ Variable declaration
682	! --------------------
683	REAL(r_k), INTENT ( IN ) :: t
684	REAL(r_k), INTENT ( IN ) :: p
685	REAL(r_k) :: theta
686
687	! Using WRF values
688	!REAL :: Rd ! Dry gas constant
689	!REAL :: Cp ! Specific heat of dry air at constant pressure
690	!REAL :: p00 ! Standard pressure ( 1000 hPa )
691	REAL(r_k) :: Rd, p00
692
693	!Rd = 287.04
694	!Cp = 1004.67
695	!p00 = 1000.00
696
697	Rd = r_d
698	p00 = p1000mb/100.
699
700	!~ Poisson's equation
701	! ------------------
702	theta = t * ( (p00/p)**(Rd/Cp) )
703
704	END FUNCTION Theta
705
706	!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
707	!~
708	!~ Name:
709	!~ Thetae
710	!~
711	!~ Description:
712	!~ This function returns equivalent potential temperature using the
713	!~ method described in Bolton 1980, Monthly Weather Review, equation 43.
714	!~
715	!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
716	FUNCTION Thetae ( tK, p, rh, mixr )
717
718	IMPLICIT NONE
719
720	!~ Variable Declarations
721	! ---------------------
722	REAL(r_k) :: tK ! Temperature ( K )
723	REAL(r_k) :: p ! Pressure ( hPa )
724	REAL(r_k) :: rh ! Relative humidity
725	REAL(r_k) :: mixr ! Mixing Ratio ( kg kg^-1)
726	REAL(r_k) :: te ! Equivalent temperature ( K )
727	REAL(r_k) :: thetae ! Equivalent potential temperature
728
729	! Using WRF values
730	!REAL, PARAMETER :: R = 287.04 ! Universal gas constant (J/deg kg)
731	!REAL, PARAMETER :: P0 = 1000.0 ! Standard pressure at surface (hPa)
732	REAL(r_k) :: R, p00, Lv
733	!REAL, PARAMETER :: lv = 2.54(10*6) ! Latent heat of vaporization
734	! (J kg^-1)
735	!REAL, PARAMETER :: cp = 1004.67 ! Specific heat of dry air constant
736	! at pressure (J/deg kg)
737	REAL(r_k) :: tlc ! LCL temperature
738
739	R = r_d
740	p00 = p1000mb/100.
741	lv = XLV
742
743	!~ Calculate the temperature of the LCL
744	! ------------------------------------
745	tlc = TLCL ( tK, rh )
746
747	!~ Calculate theta-e
748	! -----------------
749	thetae = (tK * (p00/p)*( (R/Cp)(1.- ( (.28E-3)mixr1000.) ) ) )* &
750	exp( (((3.376/tlc)-.00254))*&
751	(mixr1000.(1.+(.81E-3)mixr1000.)) )
752
753	END FUNCTION Thetae
754
755	!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
756	!~
757	!~ Name:
758	!~ The2T.f90
759	!~
760	!~ Description:
761	!~ This function returns the temperature at any pressure level along a
762	!~ saturation adiabat by iteratively solving for it from the parcel
763	!~ thetae.
764	!~
765	!~ Dependencies:
766	!~ function thetae.f90
767	!~
768	!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
769	FUNCTION The2T ( thetaeK, pres, flag ) result ( tparcel )
770
771	IMPLICIT NONE
772
773	!~ Variable Declaration
774	! --------------------
775	REAL(r_k), INTENT ( IN ) :: thetaeK
776	REAL(r_k), INTENT ( IN ) :: pres
777	LOGICAL, INTENT ( INOUT ) :: flag
778	REAL(r_k) :: tparcel
779
780	REAL(r_k) :: thetaK
781	REAL(r_k) :: tovtheta
782	REAL(r_k) :: tcheck
783	REAL(r_k) :: svpr, svpr2
784	REAL(r_k) :: smixr, smixr2
785	REAL(r_k) :: thetae_check, thetae_check2
786	REAL(r_k) :: tguess_2, correction
787
788	LOGICAL :: found
789	INTEGER :: iter
790
791	! Using WRF values
792	!REAL :: R ! Dry gas constant
793	!REAL :: Cp ! Specific heat for dry air
794	!REAL :: kappa ! Rd / Cp
795	!REAL :: Lv ! Latent heat of vaporization at 0 deg. C
796	REAL(r_k) :: R, kappa, Lv
797
798	R = r_d
799	Lv = XLV
800	!R = 287.04
801	!Cp = 1004.67
802	Kappa = R/Cp
803	!Lv = 2.500E+6
804
805	!~ Make initial guess for temperature of the parcel
806	! ------------------------------------------------
807	tovtheta = (pres/100000.0)**(r/cp)
808	tparcel = thetaeK/exp(lv.012/(cp295.))*tovtheta
809
810	iter = 1
811	found = .false.
812	flag = .false.
813
814	DO
815	IF ( iter > 105 ) EXIT
816
817	tguess_2 = tparcel + REAL ( 1 )
818
819	svpr = 6.122 * exp ( (17.67*(tparcel-273.15)) / (tparcel-29.66) )
820	smixr = ( 0.622*svpr ) / ( (pres/100.0)-svpr )
821	svpr2 = 6.122 * exp ( (17.67*(tguess_2-273.15)) / (tguess_2-29.66) )
822	smixr2 = ( 0.622*svpr2 ) / ( (pres/100.0)-svpr2 )
823
824	! ------------------------------------------------------------------ ~!
825	!~ When this function was orinially written, the final parcel ~!
826	!~ temperature check was based off of the parcel temperature and ~!
827	!~ not the theta-e it produced. As there are multiple temperature- ~!
828	!~ mixing ratio combinations that can produce a single theta-e value, ~!
829	!~ we change the check to be based off of the resultant theta-e ~!
830	!~ value. This seems to be the most accurate way of backing out ~!
831	!~ temperature from theta-e. ~!
832	!~ ~!
833	!~ Rentschler, April 2010 ~!
834	! ------------------------------------------------------------------ !
835
836	!~ Old way...
837	!thetaK = thetaeK / EXP (lv * smixr /(cp*tparcel) )
838	!tcheck = thetaK * tovtheta
839
840	!~ New way
841	thetae_check = Thetae ( tparcel, pres/100., 100., smixr )
842	thetae_check2 = Thetae ( tguess_2, pres/100., 100., smixr2 )
843
844	!~ Whew doggies - that there is some accuracy...
845	!IF ( ABS (tparcel-tcheck) < .05) THEN
846	IF ( ABS (thetaeK-thetae_check) < .001) THEN
847	found = .true.
848	flag = .true.
849	EXIT
850	END IF
851
852	!~ Old
853	!tparcel = tparcel + (tcheck - tparcel)*.3
854
855	!~ New
856	correction = ( thetaeK-thetae_check ) / ( thetae_check2-thetae_check )
857	tparcel = tparcel + correction
858
859	iter = iter + 1
860	END DO
861
862	!IF ( .not. found ) THEN
863	! print*, "Warning! Thetae to temperature calculation did not converge!"
864	! print*, "Thetae ", thetaeK, "Pressure ", pres
865	!END IF
866
867	END FUNCTION The2T
868
869	!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
870	!~
871	!~ Name:
872	!~ VirtualTemperature
873	!~
874	!~ Description:
875	!~ This function returns virtual temperature given temperature ( K )
876	!~ and mixing ratio.
877	!~
878	!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
879	FUNCTION VirtualTemperature ( tK, w ) result ( Tv )
880
881	IMPLICIT NONE
882
883	!~ Variable declaration
884	real(r_k), intent ( in ) :: tK !~ Temperature
885	real(r_k), intent ( in ) :: w !~ Mixing ratio ( kg kg^-1 )
886	real(r_k) :: Tv !~ Virtual temperature
887
888	Tv = tK * ( 1.0 + (w/0.622) ) / ( 1.0 + w )
889
890	END FUNCTION VirtualTemperature
891
892	!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
893	!~
894	!~ Name:
895	!~ SaturationMixingRatio
896	!~
897	!~ Description:
898	!~ This function calculates saturation mixing ratio given the
899	!~ temperature ( K ) and the ambient pressure ( Pa ). Uses
900	!~ approximation of saturation vapor pressure.
901	!~
902	!~ References:
903	!~ Bolton (1980), Monthly Weather Review, pg. 1047, Eq. 10
904	!~
905	!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
906	FUNCTION SaturationMixingRatio ( tK, p ) result ( ws )
907
908	IMPLICIT NONE
909
910	REAL(r_k), INTENT ( IN ) :: tK
911	REAL(r_k), INTENT ( IN ) :: p
912	REAL(r_k) :: ws
913
914	REAL(r_k) :: es
915
916	es = 6.122 * exp ( (17.67*(tK-273.15))/ (tK-29.66) )
917	ws = ( 0.622*es ) / ( (p/100.0)-es )
918
919	END FUNCTION SaturationMixingRatio
920
921	!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
922	!~
923	!~ Name:
924	!~ tlcl
925	!~
926	!~ Description:
927	!~ This function calculates the temperature of a parcel of air would have
928	!~ if lifed dry adiabatically to it's lifting condensation level (lcl).
929	!~
930	!~ References:
931	!~ Bolton (1980), Monthly Weather Review, pg. 1048, Eq. 22
932	!~
933	!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
934	FUNCTION TLCL ( tk, rh )
935
936	IMPLICIT NONE
937
938	REAL(r_k), INTENT ( IN ) :: tK !~ Temperature ( K )
939	REAL(r_k), INTENT ( IN ) :: rh !~ Relative Humidity ( % )
940	REAL(r_k) :: tlcl
941
942	REAL(r_k) :: denom, term1, term2
943
944	term1 = 1.0 / ( tK - 55.0 )
945	!! Lluis
946	! IF ( rh > REAL (0) ) THEN
947	IF ( rh > zeroRK ) THEN
948	term2 = ( LOG (rh/100.0) / 2840.0 )
949	ELSE
950	term2 = ( LOG (0.001/oneRK) / 2840.0 )
951	END IF
952	denom = term1 - term2
953	!! Lluis
954	! tlcl = ( 1.0 / denom ) + REAL ( 55 )
955	tlcl = ( oneRK / denom ) + 55*oneRK
956
957	END FUNCTION TLCL
958
959	FUNCTION var_cape_afwa1D(nz, tk, rhv, p, hgt, sfc, cape, cin, zlfc, plfc, lidx, parcel) RESULT (ostat)
960	! Function to compute cape on a 1D column following implementation in phys/module_diag_afwa.F
961
962	IMPLICIT NONE
963
964	INTEGER, INTENT(in) :: nz, sfc
965	REAL(r_k), DIMENSION(nz), INTENT(in) :: tk, rhv, p, hgt
966	REAL(r_k), INTENT(out) :: cape, cin, zlfc, plfc, lidx
967	INTEGER :: ostat
968	INTEGER, INTENT(in) :: parcel
969
970	! Local
971	!~ Derived profile variables
972	! -------------------------
973	REAL(r_k), DIMENSION(nz) :: rh, ws, w, dTvK, buoy
974	REAL(r_k) :: tlclK, plcl, nbuoy, pbuoy
975
976	!~ Source parcel information
977	! -------------------------
978	REAL(r_k) :: srctK, srcrh, srcws, srcw, srcp, &
979	srctheta, srcthetaeK
980	INTEGER :: srclev
981	REAL(r_k) :: spdiff
982
983	!~ Parcel variables
984	! ----------------
985	REAL(r_k) :: ptK, ptvK, tvK, pw
986
987	!~ Other utility variables
988	! -----------------------
989	INTEGER :: i, j, k
990	INTEGER :: lfclev
991	INTEGER :: prcl
992	INTEGER :: mlev
993	INTEGER :: lyrcnt
994	LOGICAL :: flag
995	LOGICAL :: wflag
996	REAL(r_k) :: freeze
997	REAL(r_k) :: pdiff
998	REAL(r_k) :: pm, pu, pd
999	REAL(r_k) :: lidxu
1000	REAL(r_k) :: lidxd
1001
1002	REAL(r_k), PARAMETER :: Rd = r_d
1003	REAL(r_k), PARAMETER :: RUNDEF = -9.999E30
1004
1005	!!!!!!! Variables
1006	! nz: Number of vertical levels
1007	! sfc: Surface level in the profile
1008	! tk: Temperature profile [K]
1009	! rhv: Relative Humidity profile [1]
1010	! rh: Relative Humidity profile [%]
1011	! p: Pressure profile [Pa]
1012	! hgt: Geopotential height profile [gpm]
1013	! cape: CAPE [Jkg-1]
1014	! cin: CIN [Jkg-1]
1015	! zlfc: LFC Height [gpm]
1016	! plfc: LFC Pressure [Pa]
1017	! lidx: Lifted index
1018	! FROM: https://en.wikipedia.org/wiki/Lifted_index
1019	! lidx >= 6: Very Stable Conditions
1020	! 6 > lidx > 1: Stable Conditions, Thunderstorms Not Likely
1021	! 0 > lidx > -2: Slightly Unstable, Thunderstorms Possible, With Lifting Mechanism (i.e., cold front, daytime heating, ...)
1022	! -2 > lidx > -6: Unstable, Thunderstorms Likely, Some Severe With Lifting Mechanism
1023	! -6 > lidx: Very Unstable, Severe Thunderstorms Likely With Lifting Mechanism
1024	! ostat: Function return status (Nonzero is bad)
1025	! parcel:
1026	! Most Unstable = 1 (default)
1027	! Mean layer = 2
1028	! Surface based = 3
1029	!~ Derived profile variables
1030	! -------------------------
1031	! ws: Saturation mixing ratio
1032	! w: Mixing ratio
1033	! dTvK: Parcel / ambient Tv difference
1034	! buoy: Buoyancy
1035	! tlclK: LCL temperature [K]
1036	! plcl: LCL pressure [Pa]
1037	! nbuoy: Negative buoyancy
1038	! pbuoy: Positive buoyancy
1039
1040	!~ Source parcel information
1041	! -------------------------
1042	! srctK: Source parcel temperature [K]
1043	! srcrh: Source parcel rh [%]
1044	! srcws: Source parcel sat. mixing ratio
1045	! srcw: Source parcel mixing ratio
1046	! srcp: Source parcel pressure [Pa]
1047	! srctheta: Source parcel theta [K]
1048	! srcthetaeK: Source parcel theta-e [K]
1049	! srclev: Level of the source parcel
1050	! spdiff: Pressure difference
1051
1052	!~ Parcel variables
1053	! ----------------
1054	! ptK: Parcel temperature [K]
1055	! ptvK: Parcel virtual temperature [K]
1056	! tvK: Ambient virtual temperature [K]
1057	! pw: Parcel mixing ratio
1058
1059	!~ Other utility variables
1060	! -----------------------
1061	! lfclev: Level of LFC
1062	! prcl: Internal parcel type indicator
1063	! mlev: Level for ML calculation
1064	! lyrcnt: Number of layers in mean layer
1065	! flag: Dummy flag
1066	! wflag: Saturation flag
1067	! freeze: Water loading multiplier
1068	! pdiff: Pressure difference between levs
1069	! pm, pu, pd: Middle, upper, lower pressures
1070	! lidxu: Lifted index at upper level
1071	! lidxd: Lifted index at lower level
1072
1073	fname = 'var_cape_afwa'
1074
1075	!~ Initialize variables
1076	! --------------------
1077	rh = rhv*100.
1078	ostat = 0
1079	CAPE = zeroRK
1080	CIN = zeroRK
1081	ZLFC = RUNDEF
1082	PLFC = RUNDEF
1083
1084	!~ Look for submitted parcel definition
1085	!~ 1 = Most unstable
1086	!~ 2 = Mean layer
1087	!~ 3 = Surface based
1088	! -------------------------------------
1089	IF ( parcel > 3 .or. parcel < 1 ) THEN
1090	prcl = 1
1091	ELSE
1092	prcl = parcel
1093	END IF
1094
1095	!~ Initalize our parcel to be (sort of) surface based. Because of
1096	!~ issues we've been observing in the WRF model, specifically with
1097	!~ excessive surface moisture values at the surface, using a true
1098	!~ surface based parcel is resulting a more unstable environment
1099	!~ than is actually occuring. To address this, our surface parcel
1100	!~ is now going to be defined as the parcel between 25-50 hPa
1101	!~ above the surface. UPDATE - now that this routine is in WRF,
1102	!~ going to trust surface info. GAC 20140415
1103	! ----------------------------------------------------------------
1104
1105	!~ Compute mixing ratio values for the layer
1106	! -----------------------------------------
1107	DO k = sfc, nz
1108	ws ( k ) = SaturationMixingRatio ( tK(k), p(k) )
1109	w ( k ) = ( rh(k)/100.0 ) * ws ( k )
1110	END DO
1111
1112	srclev = sfc
1113	srctK = tK ( sfc )
1114	srcrh = rh ( sfc )
1115	srcp = p ( sfc )
1116	srcws = ws ( sfc )
1117	srcw = w ( sfc )
1118	srctheta = Theta ( tK(sfc), p(sfc)/100.0 )
1119
1120	!~ Compute the profile mixing ratio. If the parcel is the MU parcel,
1121	!~ define our parcel to be the most unstable parcel within the lowest
1122	!~ 180 mb.
1123	! -------------------------------------------------------------------
1124	mlev = sfc + 1
1125	DO k = sfc + 1, nz
1126
1127	!~ Identify the last layer within 100 hPa of the surface
1128	! -----------------------------------------------------
1129	pdiff = ( p (sfc) - p (k) ) / REAL ( 100 )
1130	IF ( pdiff <= REAL (100) ) mlev = k
1131
1132	!~ If we've made it past the lowest 180 hPa, exit the loop
1133	! -------------------------------------------------------
1134	IF ( pdiff >= REAL (180) ) EXIT
1135
1136	IF ( prcl == 1 ) THEN
1137	!IF ( (p(k) > 70000.0) .and. (w(k) > srcw) ) THEN
1138	IF ( (w(k) > srcw) ) THEN
1139	srctheta = Theta ( tK(k), p(k)/100.0 )
1140	srcw = w ( k )
1141	srclev = k
1142	srctK = tK ( k )
1143	srcrh = rh ( k )
1144	srcp = p ( k )
1145	END IF
1146	END IF
1147
1148	END DO
1149
1150	!~ If we want the mean layer parcel, compute the mean values in the
1151	!~ lowest 100 hPa.
1152	! ----------------------------------------------------------------
1153	lyrcnt = mlev - sfc + 1
1154	IF ( prcl == 2 ) THEN
1155
1156	srclev = sfc
1157	srctK = SUM ( tK (sfc:mlev) ) / REAL ( lyrcnt )
1158	srcw = SUM ( w (sfc:mlev) ) / REAL ( lyrcnt )
1159	srcrh = SUM ( rh (sfc:mlev) ) / REAL ( lyrcnt )
1160	srcp = SUM ( p (sfc:mlev) ) / REAL ( lyrcnt )
1161	srctheta = Theta ( srctK, srcp/100. )
1162
1163	END IF
1164
1165	srcthetaeK = Thetae ( srctK, srcp/100.0, srcrh, srcw )
1166
1167	!~ Calculate temperature and pressure of the LCL
1168	! ---------------------------------------------
1169	tlclK = TLCL ( tK(srclev), rh(srclev) )
1170	plcl = p(srclev) * ( (tlclK/tK(srclev))**(Cp/Rd) )
1171
1172	!~ Now lift the parcel
1173	! -------------------
1174
1175	buoy = REAL ( 0 )
1176	pw = srcw
1177	wflag = .false.
1178	DO k = srclev, nz
1179	IF ( p (k) <= plcl ) THEN
1180
1181	!~ The first level after we pass the LCL, we're still going to
1182	!~ lift the parcel dry adiabatically, as we haven't added the
1183	!~ the required code to switch between the dry adiabatic and moist
1184	!~ adiabatic cooling. Since the dry version results in a greater
1185	!~ temperature loss, doing that for the first step so we don't over
1186	!~ guesstimate the instability.
1187	! ----------------------------------------------------------------
1188
1189	IF ( wflag ) THEN
1190	flag = .false.
1191
1192	!~ Above the LCL, our parcel is now undergoing moist adiabatic
1193	!~ cooling. Because of the latent heating being undergone as
1194	!~ the parcel rises above the LFC, must iterative solve for the
1195	!~ parcel temperature using equivalant potential temperature,
1196	!~ which is conserved during both dry adiabatic and
1197	!~ pseudoadiabatic displacements.
1198	! --------------------------------------------------------------
1199	ptK = The2T ( srcthetaeK, p(k), flag )
1200
1201	!~ Calculate the parcel mixing ratio, which is now changing
1202	!~ as we condense moisture out of the parcel, and is equivalent
1203	!~ to the saturation mixing ratio, since we are, in theory, at
1204	!~ saturation.
1205	! ------------------------------------------------------------
1206	pw = SaturationMixingRatio ( ptK, p(k) )
1207
1208	!~ Now we can calculate the virtual temperature of the parcel
1209	!~ and the surrounding environment to assess the buoyancy.
1210	! ----------------------------------------------------------
1211	ptvK = VirtualTemperature ( ptK, pw )
1212	tvK = VirtualTemperature ( tK (k), w (k) )
1213
1214	!~ Modification to account for water loading
1215	! -----------------------------------------
1216	freeze = 0.033 * ( 263.15 - pTvK )
1217	IF ( freeze > 1.0 ) freeze = 1.0
1218	IF ( freeze < 0.0 ) freeze = 0.0
1219
1220	!~ Approximate how much of the water vapor has condensed out
1221	!~ of the parcel at this level
1222	! ---------------------------------------------------------
1223	freeze = freeze * 333700.0 * ( srcw - pw ) / 1005.7
1224
1225	pTvK = pTvK - pTvK * ( srcw - pw ) + freeze
1226	dTvK ( k ) = ptvK - tvK
1227	buoy ( k ) = g * ( dTvK ( k ) / tvK )
1228
1229	ELSE
1230
1231	!~ Since the theta remains constant whilst undergoing dry
1232	!~ adiabatic processes, can back out the parcel temperature
1233	!~ from potential temperature below the LCL
1234	! --------------------------------------------------------
1235	ptK = srctheta / ( 100000.0/p(k) )**(Rd/Cp)
1236
1237	!~ Grab the parcel virtual temperture, can use the source
1238	!~ mixing ratio since we are undergoing dry adiabatic cooling
1239	! ----------------------------------------------------------
1240	ptvK = VirtualTemperature ( ptK, srcw )
1241
1242	!~ Virtual temperature of the environment
1243	! --------------------------------------
1244	tvK = VirtualTemperature ( tK (k), w (k) )
1245
1246	!~ Buoyancy at this level
1247	! ----------------------
1248	dTvK ( k ) = ptvK - tvK
1249	buoy ( k ) = g * ( dtvK ( k ) / tvK )
1250
1251	wflag = .true.
1252
1253	END IF
1254
1255	ELSE
1256
1257	!~ Since the theta remains constant whilst undergoing dry
1258	!~ adiabatic processes, can back out the parcel temperature
1259	!~ from potential temperature below the LCL
1260	! --------------------------------------------------------
1261	ptK = srctheta / ( 100000.0/p(k) )**(Rd/Cp)
1262
1263	!~ Grab the parcel virtual temperture, can use the source
1264	!~ mixing ratio since we are undergoing dry adiabatic cooling
1265	! ----------------------------------------------------------
1266	ptvK = VirtualTemperature ( ptK, srcw )
1267
1268	!~ Virtual temperature of the environment
1269	! --------------------------------------
1270	tvK = VirtualTemperature ( tK (k), w (k) )
1271
1272	!~ Buoyancy at this level
1273	! ---------------------
1274	dTvK ( k ) = ptvK - tvK
1275	buoy ( k ) = g * ( dtvK ( k ) / tvK )
1276
1277	END IF
1278
1279	!~ Chirp
1280	! -----
1281	! WRITE ( ,'(I15,6F15.3)' )k,p(k)/100.,ptK,pw1000.,ptvK,tvK,buoy(k)
1282
1283	END DO
1284
1285	!~ Add up the buoyancies, find the LFC
1286	! -----------------------------------
1287	flag = .false.
1288	lfclev = -1
1289	nbuoy = REAL ( 0 )
1290	pbuoy = REAL ( 0 )
1291	DO k = sfc + 1, nz
1292	IF ( tK (k) < 253.15 ) EXIT
1293	CAPE = CAPE + MAX ( buoy (k), 0.0 ) * ( hgt (k) - hgt (k-1) )
1294	CIN = CIN + MIN ( buoy (k), 0.0 ) * ( hgt (k) - hgt (k-1) )
1295
1296	!~ If we've already passed the LFC
1297	! -------------------------------
1298	IF ( flag .and. buoy (k) > REAL (0) ) THEN
1299	pbuoy = pbuoy + buoy (k)
1300	END IF
1301
1302	!~ We are buoyant now - passed the LFC
1303	! -----------------------------------
1304	IF ( .not. flag .and. buoy (k) > REAL (0) .and. p (k) < plcl ) THEN
1305	flag = .true.
1306	pbuoy = pbuoy + buoy (k)
1307	lfclev = k
1308	END IF
1309
1310	!~ If we think we've passed the LFC, but encounter a negative layer
1311	!~ start adding it up.
1312	! ----------------------------------------------------------------
1313	IF ( flag .and. buoy (k) < REAL (0) ) THEN
1314	nbuoy = nbuoy + buoy (k)
1315
1316	!~ If the accumulated negative buoyancy is greater than the
1317	!~ positive buoyancy, then we are capped off. Got to go higher
1318	!~ to find the LFC. Reset positive and negative buoyancy summations
1319	! ----------------------------------------------------------------
1320	IF ( ABS (nbuoy) > pbuoy ) THEN
1321	flag = .false.
1322	nbuoy = REAL ( 0 )
1323	pbuoy = REAL ( 0 )
1324	lfclev = -1
1325	END IF
1326	END IF
1327
1328	END DO
1329
1330	!~ Calculate lifted index by interpolating difference between
1331	!~ parcel and ambient Tv to 500mb.
1332	! ----------------------------------------------------------
1333	DO k = sfc + 1, nz
1334
1335	pm = 50000.
1336	pu = p ( k )
1337	pd = p ( k - 1 )
1338
1339	!~ If we're already above 500mb just set lifted index to 0.
1340	!~ --------------------------------------------------------
1341	IF ( pd .le. pm ) THEN
1342	lidx = zeroRK
1343	EXIT
1344
1345	ELSEIF ( pu .le. pm .and. pd .gt. pm) THEN
1346
1347	!~ Found trapping pressure: up, middle, down.
1348	!~ We are doing first order interpolation.
1349	! ------------------------------------------
1350	lidxu = -dTvK ( k ) * ( pu / 100000. ) ** (Rd/Cp)
1351	lidxd = -dTvK ( k-1 ) * ( pd / 100000. ) ** (Rd/Cp)
1352	lidx = ( lidxu * (pm-pd) + lidxd * (pu-pm) ) / (pu-pd)
1353	EXIT
1354
1355	ENDIF
1356
1357	END DO
1358
1359	!~ Assuming the the LFC is at a pressure level for now
1360	! ---------------------------------------------------
1361	IF ( lfclev > zeroRK ) THEN
1362	PLFC = p ( lfclev )
1363	ZLFC = hgt ( lfclev )
1364	END IF
1365
1366	IF ( PLFC /= PLFC .OR. PLFC < zeroRK ) THEN
1367	PLFC = -oneRK
1368	ZLFC = -oneRK
1369	END IF
1370
1371	IF ( CAPE /= CAPE ) cape = zeroRK
1372
1373	IF ( CIN /= CIN ) cin = zeroRK
1374
1375	!~ Chirp
1376	! -----
1377	! WRITE ( , ) ' CAPE: ', cape, ' CIN: ', cin
1378	! WRITE ( , ) ' LFC: ', ZLFC, ' PLFC: ', PLFC
1379	! WRITE ( , ) ''
1380	! WRITE ( , ) ' Exiting buoyancy.'
1381	! WRITE ( , ) ' ==================================== '
1382	! WRITE ( , ) ''
1383
1384	RETURN
1385
1386	END FUNCTION var_cape_afwa1D
1387
1388	! ---- END modified from module_diag_afwa.F ---- !
1389
1390	SUBROUTINE compute_cape_afwa4D(ta, hur, press, zg, hgt, cape, cin, zlfc, plfc, li, parcelmethod, &
1391	d1, d2, d3, d4)
1392	! Subroutine to use WRF phys/module_diag_afwa.F `buyoancy' subroutine to compute CAPE, CIN, ZLFC, PLFC, LI
1393
1394	IMPLICIT NONE
1395
1396	INTEGER, INTENT(in) :: d1, d2, d3, d4, parcelmethod
1397	REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in) :: ta, hur, press, zg
1398	REAL(r_k), DIMENSION(d1,d2), INTENT(in) :: hgt
1399	REAL(r_k), DIMENSION(d1,d2,d4), INTENT(out) :: cape, cin, zlfc, plfc, li
1400
1401	! Local
1402	INTEGER :: i, j, it
1403	INTEGER :: ofunc
1404
1405	!!!!!!! Variables
1406	! ta: air temperature [K]
1407	! hur: relative humidity [%]
1408	! press: air pressure [Pa]
1409	! zg: geopotential height [gpm]
1410	! hgt: topographical height [m]
1411	! cape: Convective available potential energy [Jkg-1]
1412	! cin: Convective inhibition [Jkg-1]
1413	! zlfc: height at the Level of free convection [m]
1414	! plfc: pressure at the Level of free convection [Pa]
1415	! li: lifted index [1]
1416	! parcelmethod:
1417	! Most Unstable = 1 (default)
1418	! Mean layer = 2
1419	! Surface based = 3
1420
1421	fname = 'compute_cape_afwa4D'
1422
1423	DO i=1, d1
1424	DO j=1, d2
1425	DO it=1, d4
1426	ofunc = var_cape_afwa1D(d3, ta(i,j,:,it), hur(i,j,:,it), press(i,j,:,it), zg(i,j,:,it), &
1427	1, cape(i,j,it), cin(i,j,it), zlfc(i,j,it), plfc(i,j,it), li(i,j,it), parcelmethod)
1428	zlfc(i,j,it) = zlfc(i,j,it)/g - hgt(i,j)
1429	END DO
1430	END DO
1431	END DO
1432
1433	RETURN
1434
1435	END SUBROUTINE compute_cape_afwa4D
1436
1437	END MODULE module_ForDiagnostics

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

Download in other formats:

Original Format