Context Navigation

module_ForDiagnosticsVars.f90 @ 1784

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

Adding:

`WRFzwind_log': extrapolate the wind at a given height following the 'logarithmic law' methodology

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

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

Download in other formats:

Original Format