Context Navigation

module_ForDiagnosticsVars.f90 @ 1777

Last change on this file since 1777 was 1777, checked in by lfita, 7 years ago

Adding:

`WRFz': as height above surface unstaggering WRFgeop
Fixing zwind computation

File size: 38.9 KB

Line
1	!! Fortran version of different diagnostics
2	! L. Fita. LMD May 2016
3	! gfortran module_generic.o -c module_ForDiagnostics.F90
4	!
5	! f2py -m module_ForDiagnostics --f90exec=/usr/bin/gfortran-4.7 -c module_generic.F90 module_ForDiagnostics.F90
6
7	MODULE module_ForDiagnosticsVars
8
9	USE module_definitions
10	USE module_generic
11
12	IMPLICIT NONE
13
14	CONTAINS
15
16	!!!!!!! Variables
17	! compute_psl_ptarget4d2: Compute sea level pressure using a target pressure. Similar to the Benjamin
18	! and Miller (1990). Method found in p_interp.F90
19	! compute_tv4d: 4D calculation of virtual temperaure
20	! SaturationMixingRatio: WRF's AFWA method to compute the saturation mixing ratio
21	! The2T: WRF's AFWA method to compute the temperature at any pressure level along a saturation adiabat
22	! by iteratively solving for it from the parcel thetae.
23	! Theta: WRF's AFWA method to compute potential temperature
24	! Thetae: WRF's AFWA method to compute equivalent potential temperature
25	! TLCL: WRF's AFWA method to compute the temperature of a parcel of air would have if lifed dry
26	! adiabatically to it's lifting condensation level (lcl)
27	! var_cape_afwa1D: WRF's AFWA method to compute cape, cin, fclp, fclz and li
28	! var_cllmh: low, medium, high-cloud [0,1]
29	! var_clt: total cloudiness [0,1]
30	! var_zmla_generic: Subroutine to compute pbl-height following a generic method
31	! var_zwind: Subroutine to extrapolate the wind at a given height following the 'power law' methodology
32	! VirtualTemp1D: Function for 1D calculation of virtual temperaure
33	! VirtualTemperature: WRF's AFWA method to compute virtual temperature
34
35	!!!!!!! Calculations
36	! compute_clt: Computation of total cloudiness
37
38	!!!
39	! Variables
40	!!!
41
42	FUNCTION var_cllmh(clfra, p, dz)
43	! Function to compute cllmh on a 1D column 1: low-cloud; 2: medium-cloud; 3: high-cloud [1]
44
45	IMPLICIT NONE
46
47	INTEGER, INTENT(in) :: dz
48	REAL(r_k), DIMENSION(dz), INTENT(in) :: clfra, p
49	REAL(r_k), DIMENSION(3) :: var_cllmh
50
51	! Local
52	INTEGER :: iz
53	REAL(r_k) :: zclearl, zcloudl, zclearm, zcloudm, &
54	zclearh, zcloudh
55
56	!!!!!!! Variables
57	! clfra: cloudfraction as 1D verical-column [1]
58	! p: pressure values of the column
59	fname = 'var_cllmh'
60
61	zclearl = oneRK
62	zcloudl = zeroRK
63	zclearm = oneRK
64	zcloudm = zeroRK
65	zclearh = oneRK
66	zcloudh = zeroRK
67
68	var_cllmh = oneRK
69
70	DO iz=1, dz
71	IF (p(iz) < prmhc) THEN
72	var_cllmh(3) = var_cllmh(3)*(oneRK-MAX(clfra(iz),zcloudh))/(oneRK-MIN(zcloudh,oneRK-ZEPSEC))
73	zcloudh = clfra(iz)
74	ELSE IF ( (p(iz) >= prmhc) .AND. (p(iz) < prmlc) ) THEN
75	var_cllmh(2) = var_cllmh(2)*(oneRK-MAX(clfra(iz),zcloudm))/(oneRK-MIN(zcloudm,oneRK-ZEPSEC))
76	zcloudm = clfra(iz)
77	ELSE IF (p(iz) >= prmlc) THEN
78	var_cllmh(1) = var_cllmh(1)*(oneRK-MAX(clfra(iz),zcloudl))/(oneRK-MIN(zcloudl,oneRK-ZEPSEC))
79	zcloudl = clfra(iz)
80	ELSE
81	PRINT *,' ' // TRIM(fname) // ': This is weird, pressure:', p(iz), ' Pa fails out!!'
82	PRINT *,' from high, low cloud pressures:', prmhc, ' ,', prmlc,' Pa at z-level:', iz
83	PRINT *,' p_high > p:', prmhc,'> ',p(iz),' Pa'
84	PRINT *,' p_low > p >= p_high:', prmlc,'> ',p(iz),' >=', prmhc,' Pa'
85	PRINT *,' p_low >= p:', prmlc,'>= ',p(iz),' Pa'
86	STOP
87	END IF
88	END DO
89
90	var_cllmh = oneRK - var_cllmh
91
92	RETURN
93
94	END FUNCTION var_cllmh
95
96	REAL(r_k) FUNCTION var_clt(clfra, dz)
97	! Function to compute the total cloud following 'newmicro.F90' from LMDZ using 1D vertical
98	! column values
99
100	IMPLICIT NONE
101
102	INTEGER, INTENT(in) :: dz
103	REAL(r_k), DIMENSION(dz), INTENT(in) :: clfra
104	! Local
105	INTEGER :: iz
106	REAL(r_k) :: zclear, zcloud
107
108	!!!!!!! Variables
109	! cfra: 1-column cloud fraction values
110
111	fname = 'var_clt'
112
113	zclear = oneRK
114	zcloud = zeroRK
115
116	DO iz=1,dz
117	zclear = zclear*(oneRK-MAX(clfra(iz),zcloud))/(oneRK-MIN(zcloud,1.-ZEPSEC))
118	var_clt = oneRK - zclear
119	zcloud = clfra(iz)
120	END DO
121
122	RETURN
123
124	END FUNCTION var_clt
125
126
127	SUBROUTINE compute_psl_ptarget4d2(press, ps, hgt, ta, qv, ptarget, psl, d1, d2, d3, d4)
128	! Subroutine to compute sea level pressure using a target pressure. Similar to the Benjamin
129	! and Miller (1990). Method found in p_interp.F90
130
131	IMPLICIT NONE
132
133	INTEGER, INTENT(in) :: d1, d2, d3, d4
134	REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in) :: press, ta, qv
135	REAL(r_k), DIMENSION(d1,d2,d4), INTENT(in) :: ps
136	REAL(r_k), DIMENSION(d1,d2), INTENT(in) :: hgt
137	REAL(r_k), INTENT(in) :: ptarget
138	REAL(r_k), DIMENSION(d1,d2,d4), INTENT(out) :: psl
139
140	! Local
141	INTEGER :: i, j, it
142	INTEGER :: kin
143	INTEGER :: kupper
144	REAL(r_k) :: dpmin, dp, tbotextrap, &
145	tvbotextrap, virtual
146	! Exponential related to standard atmosphere lapse rate r_d*gammav/g
147	REAL(r_k), PARAMETER :: expon=r_d*gammav/grav
148
149	!!!!!!! Variables
150	! press: Atmospheric pressure [Pa]
151	! ps: surface pressure [Pa]
152	! hgt: surface height
153	! ta: temperature [K]
154	! qv: water vapor mixing ratio
155	! dz: number of vertical levels
156	! psl: sea-level pressure
157
158	fname = 'compute_psl_ptarget4d2'
159
160	! Minimal distance between pressures [Pa]
161	dpmin=1.e4
162	psl=0.
163
164	DO i=1,d1
165	DO j=1,d2
166	IF (hgt(i,j) /= 0.) THEN
167	DO it=1,d4
168
169	! target pressure to be used for the extrapolation [Pa] (defined in namelist.input)
170	! ptarget = 70000. default value
171
172	! We are below both the ground and the lowest data level.
173
174	! First, find the model level that is closest to a "target" pressure
175	! level, where the "target" pressure is delta-p less that the local
176	! value of a horizontally smoothed surface pressure field. We use
177	! delta-p = 150 hPa here. A standard lapse rate temperature profile
178	! passing through the temperature at this model level will be used
179	! to define the temperature profile below ground. This is similar
180	! to the Benjamin and Miller (1990) method, using
181	! 700 hPa everywhere for the "target" pressure.
182
183	kupper = 0
184	loop_kIN: DO kin=d3,1,-1
185	kupper = kin
186	dp=abs( press(i,j,kin,it) - ptarget )
187	IF (dp .GT. dpmin) EXIT loop_kIN
188	dpmin=min(dpmin,dp)
189	ENDDO loop_kIN
190
191	tbotextrap=ta(i,j,kupper,it)(ps(i,j,it)/ptarget)*expon
192	tvbotextrap=virtualTemp1D(tbotextrap,qv(i,j,kupper,it))
193
194	psl(i,j,it) = ps(i,j,it)((tvbotextrap+gammavhgt(i,j))/tvbotextrap)**(1/expon)
195	END DO
196	ELSE
197	psl(i,j,:) = ps(i,j,:)
198	END IF
199	END DO
200	END DO
201
202	RETURN
203
204	END SUBROUTINE compute_psl_ptarget4d2
205
206	SUBROUTINE compute_tv4d(ta,qv,tv,d1,d2,d3,d4)
207	! 4D calculation of virtual temperaure
208
209	IMPLICIT NONE
210
211	INTEGER, INTENT(in) :: d1, d2, d3, d4
212	REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in) :: ta, qv
213	REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(out) :: tv
214
215	! Variables
216	! ta: temperature [K]
217	! qv: mixing ratio [kgkg-1]
218	! tv: virtual temperature
219
220	tv = ta*(oneRK+(qv/epsilonv))/(oneRK+qv)
221
222	END SUBROUTINE compute_tv4d
223
224	FUNCTION VirtualTemp1D (ta,qv) result (tv)
225	! 1D calculation of virtual temperaure
226
227	IMPLICIT NONE
228
229	REAL(r_k), INTENT(in) :: ta, qv
230	REAL(r_k) :: tv
231
232	! Variables
233	! ta: temperature [K]
234	! qv: mixing ratio [kgkg-1]
235
236	tv = ta*(oneRK+(qv/epsilonv))/(oneRK+qv)
237
238	END FUNCTION VirtualTemp1D
239
240	! ---- BEGIN modified from module_diag_afwa.F ---- !
241
242	!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
243	!~
244	!~ Name:
245	!~ Theta
246	!~
247	!~ Description:
248	!~ This function calculates potential temperature as defined by
249	!~ Poisson's equation, given temperature and pressure ( hPa ).
250	!~
251	!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
252	FUNCTION Theta ( t, p )
253
254	IMPLICIT NONE
255
256	!~ Variable declaration
257	! --------------------
258	REAL(r_k), INTENT ( IN ) :: t
259	REAL(r_k), INTENT ( IN ) :: p
260	REAL(r_k) :: theta
261
262	! Using WRF values
263	!REAL :: Rd ! Dry gas constant
264	!REAL :: Cp ! Specific heat of dry air at constant pressure
265	!REAL :: p00 ! Standard pressure ( 1000 hPa )
266	REAL(r_k) :: Rd, p00
267
268	!Rd = 287.04
269	!Cp = 1004.67
270	!p00 = 1000.00
271
272	Rd = r_d
273	p00 = p1000mb/100.
274
275	!~ Poisson's equation
276	! ------------------
277	theta = t * ( (p00/p)**(Rd/Cp) )
278
279	END FUNCTION Theta
280
281	!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
282	!~
283	!~ Name:
284	!~ Thetae
285	!~
286	!~ Description:
287	!~ This function returns equivalent potential temperature using the
288	!~ method described in Bolton 1980, Monthly Weather Review, equation 43.
289	!~
290	!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
291	FUNCTION Thetae ( tK, p, rh, mixr )
292
293	IMPLICIT NONE
294
295	!~ Variable Declarations
296	! ---------------------
297	REAL(r_k) :: tK ! Temperature ( K )
298	REAL(r_k) :: p ! Pressure ( hPa )
299	REAL(r_k) :: rh ! Relative humidity
300	REAL(r_k) :: mixr ! Mixing Ratio ( kg kg^-1)
301	REAL(r_k) :: te ! Equivalent temperature ( K )
302	REAL(r_k) :: thetae ! Equivalent potential temperature
303
304	! Using WRF values
305	!REAL, PARAMETER :: R = 287.04 ! Universal gas constant (J/deg kg)
306	!REAL, PARAMETER :: P0 = 1000.0 ! Standard pressure at surface (hPa)
307	REAL(r_k) :: R, p00, Lv
308	!REAL, PARAMETER :: lv = 2.54(10*6) ! Latent heat of vaporization
309	! (J kg^-1)
310	!REAL, PARAMETER :: cp = 1004.67 ! Specific heat of dry air constant
311	! at pressure (J/deg kg)
312	REAL(r_k) :: tlc ! LCL temperature
313
314	R = r_d
315	p00 = p1000mb/100.
316	lv = XLV
317
318	!~ Calculate the temperature of the LCL
319	! ------------------------------------
320	tlc = TLCL ( tK, rh )
321
322	!~ Calculate theta-e
323	! -----------------
324	thetae = (tK * (p00/p)*( (R/Cp)(1.- ( (.28E-3)mixr1000.) ) ) )* &
325	exp( (((3.376/tlc)-.00254))*&
326	(mixr1000.(1.+(.81E-3)mixr1000.)) )
327
328	END FUNCTION Thetae
329
330	!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
331	!~
332	!~ Name:
333	!~ The2T.f90
334	!~
335	!~ Description:
336	!~ This function returns the temperature at any pressure level along a
337	!~ saturation adiabat by iteratively solving for it from the parcel
338	!~ thetae.
339	!~
340	!~ Dependencies:
341	!~ function thetae.f90
342	!~
343	!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
344	FUNCTION The2T ( thetaeK, pres, flag ) result ( tparcel )
345
346	IMPLICIT NONE
347
348	!~ Variable Declaration
349	! --------------------
350	REAL(r_k), INTENT ( IN ) :: thetaeK
351	REAL(r_k), INTENT ( IN ) :: pres
352	LOGICAL, INTENT ( INOUT ) :: flag
353	REAL(r_k) :: tparcel
354
355	REAL(r_k) :: thetaK
356	REAL(r_k) :: tovtheta
357	REAL(r_k) :: tcheck
358	REAL(r_k) :: svpr, svpr2
359	REAL(r_k) :: smixr, smixr2
360	REAL(r_k) :: thetae_check, thetae_check2
361	REAL(r_k) :: tguess_2, correction
362
363	LOGICAL :: found
364	INTEGER :: iter
365
366	! Using WRF values
367	!REAL :: R ! Dry gas constant
368	!REAL :: Cp ! Specific heat for dry air
369	!REAL :: kappa ! Rd / Cp
370	!REAL :: Lv ! Latent heat of vaporization at 0 deg. C
371	REAL(r_k) :: R, kappa, Lv
372
373	R = r_d
374	Lv = XLV
375	!R = 287.04
376	!Cp = 1004.67
377	Kappa = R/Cp
378	!Lv = 2.500E+6
379
380	!~ Make initial guess for temperature of the parcel
381	! ------------------------------------------------
382	tovtheta = (pres/100000.0)**(r/cp)
383	tparcel = thetaeK/exp(lv.012/(cp295.))*tovtheta
384
385	iter = 1
386	found = .false.
387	flag = .false.
388
389	DO
390	IF ( iter > 105 ) EXIT
391
392	tguess_2 = tparcel + REAL ( 1 )
393
394	svpr = 6.122 * exp ( (17.67*(tparcel-273.15)) / (tparcel-29.66) )
395	smixr = ( 0.622*svpr ) / ( (pres/100.0)-svpr )
396	svpr2 = 6.122 * exp ( (17.67*(tguess_2-273.15)) / (tguess_2-29.66) )
397	smixr2 = ( 0.622*svpr2 ) / ( (pres/100.0)-svpr2 )
398
399	! ------------------------------------------------------------------ ~!
400	!~ When this function was orinially written, the final parcel ~!
401	!~ temperature check was based off of the parcel temperature and ~!
402	!~ not the theta-e it produced. As there are multiple temperature- ~!
403	!~ mixing ratio combinations that can produce a single theta-e value, ~!
404	!~ we change the check to be based off of the resultant theta-e ~!
405	!~ value. This seems to be the most accurate way of backing out ~!
406	!~ temperature from theta-e. ~!
407	!~ ~!
408	!~ Rentschler, April 2010 ~!
409	! ------------------------------------------------------------------ !
410
411	!~ Old way...
412	!thetaK = thetaeK / EXP (lv * smixr /(cp*tparcel) )
413	!tcheck = thetaK * tovtheta
414
415	!~ New way
416	thetae_check = Thetae ( tparcel, pres/100., 100., smixr )
417	thetae_check2 = Thetae ( tguess_2, pres/100., 100., smixr2 )
418
419	!~ Whew doggies - that there is some accuracy...
420	!IF ( ABS (tparcel-tcheck) < .05) THEN
421	IF ( ABS (thetaeK-thetae_check) < .001) THEN
422	found = .true.
423	flag = .true.
424	EXIT
425	END IF
426
427	!~ Old
428	!tparcel = tparcel + (tcheck - tparcel)*.3
429
430	!~ New
431	correction = ( thetaeK-thetae_check ) / ( thetae_check2-thetae_check )
432	tparcel = tparcel + correction
433
434	iter = iter + 1
435	END DO
436
437	!IF ( .not. found ) THEN
438	! print*, "Warning! Thetae to temperature calculation did not converge!"
439	! print*, "Thetae ", thetaeK, "Pressure ", pres
440	!END IF
441
442	END FUNCTION The2T
443
444	!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
445	!~
446	!~ Name:
447	!~ VirtualTemperature
448	!~
449	!~ Description:
450	!~ This function returns virtual temperature given temperature ( K )
451	!~ and mixing ratio.
452	!~
453	!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
454	FUNCTION VirtualTemperature ( tK, w ) result ( Tv )
455
456	IMPLICIT NONE
457
458	!~ Variable declaration
459	real(r_k), intent ( in ) :: tK !~ Temperature
460	real(r_k), intent ( in ) :: w !~ Mixing ratio ( kg kg^-1 )
461	real(r_k) :: Tv !~ Virtual temperature
462
463	Tv = tK * ( 1.0 + (w/0.622) ) / ( 1.0 + w )
464
465	END FUNCTION VirtualTemperature
466
467	!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
468	!~
469	!~ Name:
470	!~ SaturationMixingRatio
471	!~
472	!~ Description:
473	!~ This function calculates saturation mixing ratio given the
474	!~ temperature ( K ) and the ambient pressure ( Pa ). Uses
475	!~ approximation of saturation vapor pressure.
476	!~
477	!~ References:
478	!~ Bolton (1980), Monthly Weather Review, pg. 1047, Eq. 10
479	!~
480	!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
481	FUNCTION SaturationMixingRatio ( tK, p ) result ( ws )
482
483	IMPLICIT NONE
484
485	REAL(r_k), INTENT ( IN ) :: tK
486	REAL(r_k), INTENT ( IN ) :: p
487	REAL(r_k) :: ws
488
489	REAL(r_k) :: es
490
491	es = 6.122 * exp ( (17.67*(tK-273.15))/ (tK-29.66) )
492	ws = ( 0.622*es ) / ( (p/100.0)-es )
493
494	END FUNCTION SaturationMixingRatio
495
496	!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
497	!~
498	!~ Name:
499	!~ tlcl
500	!~
501	!~ Description:
502	!~ This function calculates the temperature of a parcel of air would have
503	!~ if lifed dry adiabatically to it's lifting condensation level (lcl).
504	!~
505	!~ References:
506	!~ Bolton (1980), Monthly Weather Review, pg. 1048, Eq. 22
507	!~
508	!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
509	FUNCTION TLCL ( tk, rh )
510
511	IMPLICIT NONE
512
513	REAL(r_k), INTENT ( IN ) :: tK !~ Temperature ( K )
514	REAL(r_k), INTENT ( IN ) :: rh !~ Relative Humidity ( % )
515	REAL(r_k) :: tlcl
516
517	REAL(r_k) :: denom, term1, term2
518
519	term1 = 1.0 / ( tK - 55.0 )
520	!! Lluis
521	! IF ( rh > REAL (0) ) THEN
522	IF ( rh > zeroRK ) THEN
523	term2 = ( LOG (rh/100.0) / 2840.0 )
524	ELSE
525	term2 = ( LOG (0.001/oneRK) / 2840.0 )
526	END IF
527	denom = term1 - term2
528	!! Lluis
529	! tlcl = ( 1.0 / denom ) + REAL ( 55 )
530	tlcl = ( oneRK / denom ) + 55*oneRK
531
532	END FUNCTION TLCL
533
534	FUNCTION var_cape_afwa1D(nz, tk, rhv, p, hgt, sfc, cape, cin, zlfc, plfc, lidx, parcel) RESULT (ostat)
535	! Function to compute cape on a 1D column following implementation in phys/module_diag_afwa.F
536
537	IMPLICIT NONE
538
539	INTEGER, INTENT(in) :: nz, sfc
540	REAL(r_k), DIMENSION(nz), INTENT(in) :: tk, rhv, p, hgt
541	REAL(r_k), INTENT(out) :: cape, cin, zlfc, plfc, lidx
542	INTEGER :: ostat
543	INTEGER, INTENT(in) :: parcel
544
545	! Local
546	!~ Derived profile variables
547	! -------------------------
548	REAL(r_k), DIMENSION(nz) :: rh, ws, w, dTvK, buoy
549	REAL(r_k) :: tlclK, plcl, nbuoy, pbuoy
550
551	!~ Source parcel information
552	! -------------------------
553	REAL(r_k) :: srctK, srcrh, srcws, srcw, srcp, &
554	srctheta, srcthetaeK
555	INTEGER :: srclev
556	REAL(r_k) :: spdiff
557
558	!~ Parcel variables
559	! ----------------
560	REAL(r_k) :: ptK, ptvK, tvK, pw
561
562	!~ Other utility variables
563	! -----------------------
564	INTEGER :: i, j, k
565	INTEGER :: lfclev
566	INTEGER :: prcl
567	INTEGER :: mlev
568	INTEGER :: lyrcnt
569	LOGICAL :: flag
570	LOGICAL :: wflag
571	REAL(r_k) :: freeze
572	REAL(r_k) :: pdiff
573	REAL(r_k) :: pm, pu, pd
574	REAL(r_k) :: lidxu
575	REAL(r_k) :: lidxd
576
577	REAL(r_k), PARAMETER :: Rd = r_d
578	REAL(r_k), PARAMETER :: RUNDEF = -9.999E30
579
580	!!!!!!! Variables
581	! nz: Number of vertical levels
582	! sfc: Surface level in the profile
583	! tk: Temperature profile [K]
584	! rhv: Relative Humidity profile [1]
585	! rh: Relative Humidity profile [%]
586	! p: Pressure profile [Pa]
587	! hgt: Geopotential height profile [gpm]
588	! cape: CAPE [Jkg-1]
589	! cin: CIN [Jkg-1]
590	! zlfc: LFC Height [gpm]
591	! plfc: LFC Pressure [Pa]
592	! lidx: Lifted index
593	! FROM: https://en.wikipedia.org/wiki/Lifted_index
594	! lidx >= 6: Very Stable Conditions
595	! 6 > lidx > 1: Stable Conditions, Thunderstorms Not Likely
596	! 0 > lidx > -2: Slightly Unstable, Thunderstorms Possible, With Lifting Mechanism (i.e., cold front, daytime heating, ...)
597	! -2 > lidx > -6: Unstable, Thunderstorms Likely, Some Severe With Lifting Mechanism
598	! -6 > lidx: Very Unstable, Severe Thunderstorms Likely With Lifting Mechanism
599	! ostat: Function return status (Nonzero is bad)
600	! parcel:
601	! Most Unstable = 1 (default)
602	! Mean layer = 2
603	! Surface based = 3
604	!~ Derived profile variables
605	! -------------------------
606	! ws: Saturation mixing ratio
607	! w: Mixing ratio
608	! dTvK: Parcel / ambient Tv difference
609	! buoy: Buoyancy
610	! tlclK: LCL temperature [K]
611	! plcl: LCL pressure [Pa]
612	! nbuoy: Negative buoyancy
613	! pbuoy: Positive buoyancy
614
615	!~ Source parcel information
616	! -------------------------
617	! srctK: Source parcel temperature [K]
618	! srcrh: Source parcel rh [%]
619	! srcws: Source parcel sat. mixing ratio
620	! srcw: Source parcel mixing ratio
621	! srcp: Source parcel pressure [Pa]
622	! srctheta: Source parcel theta [K]
623	! srcthetaeK: Source parcel theta-e [K]
624	! srclev: Level of the source parcel
625	! spdiff: Pressure difference
626
627	!~ Parcel variables
628	! ----------------
629	! ptK: Parcel temperature [K]
630	! ptvK: Parcel virtual temperature [K]
631	! tvK: Ambient virtual temperature [K]
632	! pw: Parcel mixing ratio
633
634	!~ Other utility variables
635	! -----------------------
636	! lfclev: Level of LFC
637	! prcl: Internal parcel type indicator
638	! mlev: Level for ML calculation
639	! lyrcnt: Number of layers in mean layer
640	! flag: Dummy flag
641	! wflag: Saturation flag
642	! freeze: Water loading multiplier
643	! pdiff: Pressure difference between levs
644	! pm, pu, pd: Middle, upper, lower pressures
645	! lidxu: Lifted index at upper level
646	! lidxd: Lifted index at lower level
647
648	fname = 'var_cape_afwa'
649
650	!~ Initialize variables
651	! --------------------
652	rh = rhv*100.
653	ostat = 0
654	CAPE = zeroRK
655	CIN = zeroRK
656	ZLFC = RUNDEF
657	PLFC = RUNDEF
658
659	!~ Look for submitted parcel definition
660	!~ 1 = Most unstable
661	!~ 2 = Mean layer
662	!~ 3 = Surface based
663	! -------------------------------------
664	IF ( parcel > 3 .or. parcel < 1 ) THEN
665	prcl = 1
666	ELSE
667	prcl = parcel
668	END IF
669
670	!~ Initalize our parcel to be (sort of) surface based. Because of
671	!~ issues we've been observing in the WRF model, specifically with
672	!~ excessive surface moisture values at the surface, using a true
673	!~ surface based parcel is resulting a more unstable environment
674	!~ than is actually occuring. To address this, our surface parcel
675	!~ is now going to be defined as the parcel between 25-50 hPa
676	!~ above the surface. UPDATE - now that this routine is in WRF,
677	!~ going to trust surface info. GAC 20140415
678	! ----------------------------------------------------------------
679
680	!~ Compute mixing ratio values for the layer
681	! -----------------------------------------
682	DO k = sfc, nz
683	ws ( k ) = SaturationMixingRatio ( tK(k), p(k) )
684	w ( k ) = ( rh(k)/100.0 ) * ws ( k )
685	END DO
686
687	srclev = sfc
688	srctK = tK ( sfc )
689	srcrh = rh ( sfc )
690	srcp = p ( sfc )
691	srcws = ws ( sfc )
692	srcw = w ( sfc )
693	srctheta = Theta ( tK(sfc), p(sfc)/100.0 )
694
695	!~ Compute the profile mixing ratio. If the parcel is the MU parcel,
696	!~ define our parcel to be the most unstable parcel within the lowest
697	!~ 180 mb.
698	! -------------------------------------------------------------------
699	mlev = sfc + 1
700	DO k = sfc + 1, nz
701
702	!~ Identify the last layer within 100 hPa of the surface
703	! -----------------------------------------------------
704	pdiff = ( p (sfc) - p (k) ) / REAL ( 100 )
705	IF ( pdiff <= REAL (100) ) mlev = k
706
707	!~ If we've made it past the lowest 180 hPa, exit the loop
708	! -------------------------------------------------------
709	IF ( pdiff >= REAL (180) ) EXIT
710
711	IF ( prcl == 1 ) THEN
712	!IF ( (p(k) > 70000.0) .and. (w(k) > srcw) ) THEN
713	IF ( (w(k) > srcw) ) THEN
714	srctheta = Theta ( tK(k), p(k)/100.0 )
715	srcw = w ( k )
716	srclev = k
717	srctK = tK ( k )
718	srcrh = rh ( k )
719	srcp = p ( k )
720	END IF
721	END IF
722
723	END DO
724
725	!~ If we want the mean layer parcel, compute the mean values in the
726	!~ lowest 100 hPa.
727	! ----------------------------------------------------------------
728	lyrcnt = mlev - sfc + 1
729	IF ( prcl == 2 ) THEN
730
731	srclev = sfc
732	srctK = SUM ( tK (sfc:mlev) ) / REAL ( lyrcnt )
733	srcw = SUM ( w (sfc:mlev) ) / REAL ( lyrcnt )
734	srcrh = SUM ( rh (sfc:mlev) ) / REAL ( lyrcnt )
735	srcp = SUM ( p (sfc:mlev) ) / REAL ( lyrcnt )
736	srctheta = Theta ( srctK, srcp/100. )
737
738	END IF
739
740	srcthetaeK = Thetae ( srctK, srcp/100.0, srcrh, srcw )
741
742	!~ Calculate temperature and pressure of the LCL
743	! ---------------------------------------------
744	tlclK = TLCL ( tK(srclev), rh(srclev) )
745	plcl = p(srclev) * ( (tlclK/tK(srclev))**(Cp/Rd) )
746
747	!~ Now lift the parcel
748	! -------------------
749
750	buoy = REAL ( 0 )
751	pw = srcw
752	wflag = .false.
753	DO k = srclev, nz
754	IF ( p (k) <= plcl ) THEN
755
756	!~ The first level after we pass the LCL, we're still going to
757	!~ lift the parcel dry adiabatically, as we haven't added the
758	!~ the required code to switch between the dry adiabatic and moist
759	!~ adiabatic cooling. Since the dry version results in a greater
760	!~ temperature loss, doing that for the first step so we don't over
761	!~ guesstimate the instability.
762	! ----------------------------------------------------------------
763
764	IF ( wflag ) THEN
765	flag = .false.
766
767	!~ Above the LCL, our parcel is now undergoing moist adiabatic
768	!~ cooling. Because of the latent heating being undergone as
769	!~ the parcel rises above the LFC, must iterative solve for the
770	!~ parcel temperature using equivalant potential temperature,
771	!~ which is conserved during both dry adiabatic and
772	!~ pseudoadiabatic displacements.
773	! --------------------------------------------------------------
774	ptK = The2T ( srcthetaeK, p(k), flag )
775
776	!~ Calculate the parcel mixing ratio, which is now changing
777	!~ as we condense moisture out of the parcel, and is equivalent
778	!~ to the saturation mixing ratio, since we are, in theory, at
779	!~ saturation.
780	! ------------------------------------------------------------
781	pw = SaturationMixingRatio ( ptK, p(k) )
782
783	!~ Now we can calculate the virtual temperature of the parcel
784	!~ and the surrounding environment to assess the buoyancy.
785	! ----------------------------------------------------------
786	ptvK = VirtualTemperature ( ptK, pw )
787	tvK = VirtualTemperature ( tK (k), w (k) )
788
789	!~ Modification to account for water loading
790	! -----------------------------------------
791	freeze = 0.033 * ( 263.15 - pTvK )
792	IF ( freeze > 1.0 ) freeze = 1.0
793	IF ( freeze < 0.0 ) freeze = 0.0
794
795	!~ Approximate how much of the water vapor has condensed out
796	!~ of the parcel at this level
797	! ---------------------------------------------------------
798	freeze = freeze * 333700.0 * ( srcw - pw ) / 1005.7
799
800	pTvK = pTvK - pTvK * ( srcw - pw ) + freeze
801	dTvK ( k ) = ptvK - tvK
802	buoy ( k ) = g * ( dTvK ( k ) / tvK )
803
804	ELSE
805
806	!~ Since the theta remains constant whilst undergoing dry
807	!~ adiabatic processes, can back out the parcel temperature
808	!~ from potential temperature below the LCL
809	! --------------------------------------------------------
810	ptK = srctheta / ( 100000.0/p(k) )**(Rd/Cp)
811
812	!~ Grab the parcel virtual temperture, can use the source
813	!~ mixing ratio since we are undergoing dry adiabatic cooling
814	! ----------------------------------------------------------
815	ptvK = VirtualTemperature ( ptK, srcw )
816
817	!~ Virtual temperature of the environment
818	! --------------------------------------
819	tvK = VirtualTemperature ( tK (k), w (k) )
820
821	!~ Buoyancy at this level
822	! ----------------------
823	dTvK ( k ) = ptvK - tvK
824	buoy ( k ) = g * ( dtvK ( k ) / tvK )
825
826	wflag = .true.
827
828	END IF
829
830	ELSE
831
832	!~ Since the theta remains constant whilst undergoing dry
833	!~ adiabatic processes, can back out the parcel temperature
834	!~ from potential temperature below the LCL
835	! --------------------------------------------------------
836	ptK = srctheta / ( 100000.0/p(k) )**(Rd/Cp)
837
838	!~ Grab the parcel virtual temperture, can use the source
839	!~ mixing ratio since we are undergoing dry adiabatic cooling
840	! ----------------------------------------------------------
841	ptvK = VirtualTemperature ( ptK, srcw )
842
843	!~ Virtual temperature of the environment
844	! --------------------------------------
845	tvK = VirtualTemperature ( tK (k), w (k) )
846
847	!~ Buoyancy at this level
848	! ---------------------
849	dTvK ( k ) = ptvK - tvK
850	buoy ( k ) = g * ( dtvK ( k ) / tvK )
851
852	END IF
853
854	!~ Chirp
855	! -----
856	! WRITE ( ,'(I15,6F15.3)' )k,p(k)/100.,ptK,pw1000.,ptvK,tvK,buoy(k)
857
858	END DO
859
860	!~ Add up the buoyancies, find the LFC
861	! -----------------------------------
862	flag = .false.
863	lfclev = -1
864	nbuoy = REAL ( 0 )
865	pbuoy = REAL ( 0 )
866	DO k = sfc + 1, nz
867	IF ( tK (k) < 253.15 ) EXIT
868	CAPE = CAPE + MAX ( buoy (k), 0.0 ) * ( hgt (k) - hgt (k-1) )
869	CIN = CIN + MIN ( buoy (k), 0.0 ) * ( hgt (k) - hgt (k-1) )
870
871	!~ If we've already passed the LFC
872	! -------------------------------
873	IF ( flag .and. buoy (k) > REAL (0) ) THEN
874	pbuoy = pbuoy + buoy (k)
875	END IF
876
877	!~ We are buoyant now - passed the LFC
878	! -----------------------------------
879	IF ( .not. flag .and. buoy (k) > REAL (0) .and. p (k) < plcl ) THEN
880	flag = .true.
881	pbuoy = pbuoy + buoy (k)
882	lfclev = k
883	END IF
884
885	!~ If we think we've passed the LFC, but encounter a negative layer
886	!~ start adding it up.
887	! ----------------------------------------------------------------
888	IF ( flag .and. buoy (k) < REAL (0) ) THEN
889	nbuoy = nbuoy + buoy (k)
890
891	!~ If the accumulated negative buoyancy is greater than the
892	!~ positive buoyancy, then we are capped off. Got to go higher
893	!~ to find the LFC. Reset positive and negative buoyancy summations
894	! ----------------------------------------------------------------
895	IF ( ABS (nbuoy) > pbuoy ) THEN
896	flag = .false.
897	nbuoy = REAL ( 0 )
898	pbuoy = REAL ( 0 )
899	lfclev = -1
900	END IF
901	END IF
902
903	END DO
904
905	!~ Calculate lifted index by interpolating difference between
906	!~ parcel and ambient Tv to 500mb.
907	! ----------------------------------------------------------
908	DO k = sfc + 1, nz
909
910	pm = 50000.
911	pu = p ( k )
912	pd = p ( k - 1 )
913
914	!~ If we're already above 500mb just set lifted index to 0.
915	!~ --------------------------------------------------------
916	IF ( pd .le. pm ) THEN
917	lidx = zeroRK
918	EXIT
919
920	ELSEIF ( pu .le. pm .and. pd .gt. pm) THEN
921
922	!~ Found trapping pressure: up, middle, down.
923	!~ We are doing first order interpolation.
924	! ------------------------------------------
925	lidxu = -dTvK ( k ) * ( pu / 100000. ) ** (Rd/Cp)
926	lidxd = -dTvK ( k-1 ) * ( pd / 100000. ) ** (Rd/Cp)
927	lidx = ( lidxu * (pm-pd) + lidxd * (pu-pm) ) / (pu-pd)
928	EXIT
929
930	ENDIF
931
932	END DO
933
934	!~ Assuming the the LFC is at a pressure level for now
935	! ---------------------------------------------------
936	IF ( lfclev > zeroRK ) THEN
937	PLFC = p ( lfclev )
938	ZLFC = hgt ( lfclev )
939	END IF
940
941	IF ( PLFC /= PLFC .OR. PLFC < zeroRK ) THEN
942	PLFC = -oneRK
943	ZLFC = -oneRK
944	END IF
945
946	IF ( CAPE /= CAPE ) cape = zeroRK
947
948	IF ( CIN /= CIN ) cin = zeroRK
949
950	!~ Chirp
951	! -----
952	! WRITE ( , ) ' CAPE: ', cape, ' CIN: ', cin
953	! WRITE ( , ) ' LFC: ', ZLFC, ' PLFC: ', PLFC
954	! WRITE ( , ) ''
955	! WRITE ( , ) ' Exiting buoyancy.'
956	! WRITE ( , ) ' ==================================== '
957	! WRITE ( , ) ''
958
959	RETURN
960
961	END FUNCTION var_cape_afwa1D
962
963	! ---- END modified from module_diag_afwa.F ---- !
964
965	SUBROUTINE var_zmla_generic(dz, qv, tpot, z, topo, zmla)
966	! Subroutine to compute pbl-height following a generic method
967	! from Nielsen-Gammon et al., 2008 J. Appl. Meteor. Clim.
968	! applied also in Garcia-Diez et al., 2013, QJRMS
969	! where
970	! "The technique identifies the ML height as a threshold increase of potential temperature from
971	! its minimum value within the boundary layer."
972	! here applied similarly to Garcia-Diez et al. where
973	! zmla = "...first level where potential temperature exceeds the minimum potential temperature
974	! reached in the mixed layer by more than 1.5 K"
975
976	IMPLICIT NONE
977
978	INTEGER, INTENT(in) :: dz
979	REAL(r_k), DIMENSION(dz), INTENT(in) :: qv, tpot, z
980	REAL(r_k), INTENT(in) :: topo
981	REAL(r_k), INTENT(out) :: zmla
982
983	! Local
984	INTEGER :: i
985	INTEGER :: mldlev, bllev
986	REAL(r_k) :: dqvar, tpotmin, refdt
987
988	!!!!!!! Variables
989	! qv: water vapour mixing ratio
990	! tpot: potential temperature [K]
991	! z: height above sea level [m]
992	! topo: topographic height [m]
993	! zmla: boundary layer height [m]
994
995	fname = 'var_zmla_generic'
996
997	! Pecentage of difference of mixing ratio used to determine Mixed layer depth
998	dqvar = 0.1
999
1000	! MLD = Mixed layer with no substantial variation of mixing ratio /\qv < 10% ?
1001	!PRINT *,' Mixed layer mixing ratios qv[1] lev qv[lev] dqvar% _______'
1002	DO mldlev = 2, dz
1003	IF (ABS(qv(mldlev)-qv(1))/qv(1) > dqvar ) EXIT
1004	! PRINT *,qv(1), mldlev, qv(mldlev), ABS(qv(mldlev)-qv(1))/qv(1)
1005	END DO
1006
1007	! Looking for the minimum potential temperature within the MLD [tpotmin = min(tpot)_0^MLD]
1008	tpotmin = MINVAL(tpot(1:mldlev))
1009
1010	! Change in temperature to determine boundary layer height
1011	refdt = 1.5
1012
1013	! Determine the first level where tpot > tpotmin + 1.5 K
1014	!PRINT *,' Mixed layer tpotmin lev tpotmin[lev] dtpot _______'
1015	DO bllev = 1, dz
1016	IF (ABS(tpot(bllev)-tpotmin) > refdt ) EXIT
1017	! PRINT *,tpotmin, bllev, tpot(bllev), ABS(tpot(bllev)-tpotmin)
1018	END DO
1019
1020	!PRINT *,' height end MLD:', z(mldlev)
1021	!PRINT *,' pbl height:', z(bllev)
1022
1023	zmla = z(bllev) - topo
1024
1025	RETURN
1026
1027	END SUBROUTINE var_zmla_generic
1028
1029	SUBROUTINE var_zwind(d1, u, v, z, u10, v10, sa, ca, newz, unewz, vnewz)
1030	! Subroutine to extrapolate the wind at a given height following the 'power law' methodology
1031	! wss[newz] = wss[z1](newz/z1)*alpha
1032	! alpha = (ln(wss[z2])-ln(wss[z1]))/(ln(z2)-ln(z1))
1033	! AFTER: Phd Thesis:
1034	! Benedicte Jourdier. Ressource eolienne en France metropolitaine : methodes dâevaluation du
1035	! potentiel, variabilite et tendances. Climatologie. Ecole Doctorale Polytechnique, 2015. French
1036	!
1037	IMPLICIT NONE
1038
1039	INTEGER, INTENT(in) :: d1
1040	REAL(r_k), DIMENSION(d1), INTENT(in) :: u,v,z
1041	REAL(r_k), INTENT(in) :: u10, v10, sa, ca, newz
1042	REAL(r_k), INTENT(out) :: unewz, vnewz
1043
1044	! Local
1045	INTEGER :: inear
1046	REAL(r_k) :: zaground
1047	REAL(r_k), DIMENSION(2) :: v1, v2, zz, alpha, uvnewz
1048
1049	!!!!!!! Variables
1050	! u,v: vertical wind components [ms-1]
1051	! z: height above surface on half-mass levels [m]
1052	! u10,v10: 10-m wind components [ms-1]
1053	! sa, ca: local sine and cosine of map rotation [1.]
1054	! newz: desired height above grpund of extrapolation
1055	! unewz,vnewz: Wind compoonents at the given height [ms-1]
1056
1057	fname = 'var_zwind'
1058
1059	!PRINT *,' ilev zaground newz z[ilev+1] z[ilev+2] _______'
1060	IF (z(1) < newz ) THEN
1061	DO inear = 1,d1-2
1062	zaground = z(inear+2)
1063	!PRINT *, inear, z(inear), newz, z(inear+1), z(inear+2)
1064	IF ( zaground >= newz) EXIT
1065	END DO
1066	ELSE
1067	!PRINT *, 1, z(1), newz, z(2), z(3), ' z(1) > newz'
1068	inear = d1 - 2
1069	END IF
1070
1071	IF (inear == d1-2) THEN
1072	! No vertical pair of levels is below newz, using 10m wind as first value and the first level as the second
1073	v1(1) = u10
1074	v1(2) = v10
1075	v2(1) = u(1)
1076	v2(2) = v(1)
1077	zz(1) = 10.
1078	zz(2) = z(1)
1079	ELSE
1080	v1(1) = u(inear)
1081	v1(2) = v(inear)
1082	v2(1) = u(inear+1)
1083	v2(2) = v(inear+1)
1084	zz(1) = z(inear)
1085	zz(2) = z(inear+1)
1086	END IF
1087
1088	! Computing for each component
1089	alpha = (LOG(ABS(v2))-LOG(ABS(v1)))/(LOG(zz(2))-LOG(zz(1)))
1090	!PRINT *,' Computing with v1:', v1, ' ms-1 v2:', v2, ' ms-1'
1091	!PRINT *,' z1:', zz(1), 'm z2:', zz(2), ' m'
1092	!PRINT *,' alhpa u:', alpha(1), ' alpha 2:', alpha(2)
1093
1094	uvnewz = v1(newz/zz(1))*alpha
1095	! Earth-rotation
1096	unewz = uvnewz(1)ca - uvnewz(2)sa
1097	vnewz = uvnewz(1)sa + uvnewz(2)ca
1098
1099	!PRINT *,' result vz:', uvnewz
1100
1101	!STOP
1102
1103	RETURN
1104
1105	END SUBROUTINE var_zwind
1106
1107	END MODULE module_ForDiagnosticsVars

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

Download in other formats:

Original Format