Context Navigation

module_ForDiagnosticsVars.f90 @ 2601

Last change on this file since 2601 was 2387, checked in by lfita, 6 years ago

Adding:

`tws': Wet Bulb temperature after Stull, 2011

File size: 66.7 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	USE module_scientific
12
13	IMPLICIT NONE
14
15	CONTAINS
16
17	!!!!!!! Variables
18	! Cdrag_0: Fuction to compute a first order generic approximation of the drag coefficient
19	! compute_psl_ptarget4d2: Compute sea level pressure using a target pressure. Similar to the Benjamin
20	! and Miller (1990). Method found in p_interp.F90
21	! compute_tv4d: 4D calculation of virtual temperaure
22	! SaturationMixingRatio: WRF's AFWA method to compute the saturation mixing ratio
23	! The2T: WRF's AFWA method to compute the temperature at any pressure level along a saturation adiabat
24	! by iteratively solving for it from the parcel thetae.
25	! Theta: WRF's AFWA method to compute potential temperature
26	! Thetae: WRF's AFWA method to compute equivalent potential temperature
27	! TLCL: WRF's AFWA method to compute the temperature of a parcel of air would have if lifed dry
28	! adiabatically to it's lifting condensation level (lcl)
29	! var_cape_afwa1D: WRF's AFWA method to compute cape, cin, fclp, fclz and li
30	! var_cllmh: low, medium, high-cloud [0,1]
31	! var_clt: total cloudiness [0,1]
32	! var_hur: relative humidity using August-Roche-Magnus approximation [1]
33	! var_fog_K84: fog and visibility following Kunkel, (1984)
34	! var_fog_RUC: fog and visibility following RUC method Smirnova, (2000)
35	! var_fog_FRAML50: fog and visibility following Gultepe and Milbrandt, (2010)
36	! var_potevap_orPM: potential evapotranspiration following Penman-Monteith formulation implemented in ORCHIDEE
37	! var_psl_ecmwf: sea level pressure using ECMWF method following Mats Hamrud and Philippe Courtier [Pa]
38	! var_range_faces: compute faces [uphill, valleys, downhill] of a monuntain range along a face
39	! var_rh: Subroutine to compute relative humidity following 'Tetens' equation (T,P) ...'
40	! var_tws_S11: Function to compute wet bulb temperature after Stull, 2011
41	! var_zmla_generic: Subroutine to compute pbl-height following a generic method
42	! var_zwind: Subroutine to extrapolate the wind at a given height following the 'power law' methodology
43	! var_zwind_log: Subroutine to extrapolate the wind at a given height following the 'logarithmic law' methodology
44	! var_zwind_MOtheor: Subroutine of wind extrapolation following Moin-Obukhov theory
45	! VirtualTemp1D: Function for 1D calculation of virtual temperaure
46	! VirtualTemperature: WRF's AFWA method to compute virtual temperature
47	! stabfunc_businger: Fucntion of the stability function after Businger et al. (1971)
48
49	!!!!!!! Calculations
50	! compute_clt: Computation of total cloudiness
51
52	!!!
53	! Variables
54	!!!
55
56	FUNCTION var_cllmh(clfra, p, dz)
57	! Function to compute cllmh on a 1D column 1: low-cloud; 2: medium-cloud; 3: high-cloud [1]
58
59	IMPLICIT NONE
60
61	INTEGER, INTENT(in) :: dz
62	REAL(r_k), DIMENSION(dz), INTENT(in) :: clfra, p
63	REAL(r_k), DIMENSION(3) :: var_cllmh
64
65	! Local
66	INTEGER :: iz
67	REAL(r_k) :: zclearl, zcloudl, zclearm, zcloudm, &
68	zclearh, zcloudh
69
70	!!!!!!! Variables
71	! clfra: cloudfraction as 1D verical-column [1]
72	! p: pressure values of the column
73	fname = 'var_cllmh'
74
75	zclearl = oneRK
76	zcloudl = zeroRK
77	zclearm = oneRK
78	zcloudm = zeroRK
79	zclearh = oneRK
80	zcloudh = zeroRK
81
82	var_cllmh = oneRK
83
84	DO iz=1, dz
85	IF (p(iz) < prmhc) THEN
86	var_cllmh(3) = var_cllmh(3)*(oneRK-MAX(clfra(iz),zcloudh))/(oneRK-MIN(zcloudh,oneRK-ZEPSEC))
87	zcloudh = clfra(iz)
88	ELSE IF ( (p(iz) >= prmhc) .AND. (p(iz) < prmlc) ) THEN
89	var_cllmh(2) = var_cllmh(2)*(oneRK-MAX(clfra(iz),zcloudm))/(oneRK-MIN(zcloudm,oneRK-ZEPSEC))
90	zcloudm = clfra(iz)
91	ELSE IF (p(iz) >= prmlc) THEN
92	var_cllmh(1) = var_cllmh(1)*(oneRK-MAX(clfra(iz),zcloudl))/(oneRK-MIN(zcloudl,oneRK-ZEPSEC))
93	zcloudl = clfra(iz)
94	ELSE
95	PRINT *,' ' // TRIM(fname) // ': This is weird, pressure:', p(iz), ' Pa fails out!!'
96	PRINT *,' from high, low cloud pressures:', prmhc, ' ,', prmlc,' Pa at z-level:', iz
97	PRINT *,' p_high > p:', prmhc,'> ',p(iz),' Pa'
98	PRINT *,' p_low > p >= p_high:', prmlc,'> ',p(iz),' >=', prmhc,' Pa'
99	PRINT *,' p_low >= p:', prmlc,'>= ',p(iz),' Pa'
100	STOP
101	END IF
102	END DO
103
104	var_cllmh = oneRK - var_cllmh
105
106	RETURN
107
108	END FUNCTION var_cllmh
109
110	REAL(r_k) FUNCTION var_clt(clfra, dz)
111	! Function to compute the total cloud following 'newmicro.F90' from LMDZ using 1D vertical
112	! column values
113
114	IMPLICIT NONE
115
116	INTEGER, INTENT(in) :: dz
117	REAL(r_k), DIMENSION(dz), INTENT(in) :: clfra
118	! Local
119	INTEGER :: iz
120	REAL(r_k) :: zclear, zcloud
121
122	!!!!!!! Variables
123	! cfra: 1-column cloud fraction values
124
125	fname = 'var_clt'
126
127	zclear = oneRK
128	zcloud = zeroRK
129
130	DO iz=1,dz
131	zclear = zclear*(oneRK-MAX(clfra(iz),zcloud))/(oneRK-MIN(zcloud,1.-ZEPSEC))
132	var_clt = oneRK - zclear
133	zcloud = clfra(iz)
134	END DO
135
136	RETURN
137
138	END FUNCTION var_clt
139
140	SUBROUTINE var_psl_ecmwf(PRPRESS, hgt, PTB, PRESBH, PRESBF, psl)
141	! Subroutine to compute sea level pressure using ECMWF method following Mats Hamrud and Philippe Courtier
142	! method found in LMDZ in phylmd/pppmer.F90 in combination with phylmd/ctsar.F90
143
144	! IMPLICIT ARGUMENTS : CONSTANTS FROM YOMCST,YOMGEM,YOMSTA.
145	! --------------------
146
147	IMPLICIT NONE
148
149	REAL, INTENT(in) :: PRPRESS, hgt, PTB, PRESBH, PRESBF
150	REAL, INTENT(out) :: psl
151
152	! Local
153	REAL :: ghgt, PTSTAR, PT0, ZTSTAR
154	REAL :: ZALPHA, POROG
155	REAL :: ZDTDZSG, ZOROG, ZT0, ZTX, ZTY, ZX, ZY, ZY2
156	REAL, PARAMETER :: RDTDZ1 = -gammav
157
158	!!!!!!! Variables
159	! PRPRESS: Surface pressure [Pa]
160	! hgt: Terrain height [m]
161	! PTB: Temperature first half-level [K]
162	! PRESBH: Pressure first half-level [Pa]
163	! PRESBF: Pressure second full-level [Pa]
164	! psl: sea-level pressure
165
166	fname = 'var_psl_ecmwf'
167
168	! Height by gravity
169	POROG = hgt*g
170
171	!* 1. COMPUTES SURFACE TEMPERATURE
172	!* THEN STANDARD SURFACE TEMPERATURE.
173
174	ZDTDZSG=-RDTDZ1/g
175	ZALPHA=ZDTDZSG*r_d
176
177	PTSTAR=PTB(1.0+ZALPHA(PRESBH/PRESBF-1.0))
178	PT0=PTSTAR+ZDTDZSG*POROG
179
180	!* 2. POST-PROCESS MSL PRESSURE.
181	! --------------------------
182
183	!* 2.1 COMPUTATION OF MODIFIED ALPHA AND TSTAR.
184
185	ZTX=290.5
186	ZTY=255.0
187
188	IF (PTSTAR < ZTY) THEN
189	ZTSTAR=0.5*(ZTY+PTSTAR)
190	ELSEIF (PTSTAR < ZTX) THEN
191	ZTSTAR=PTSTAR
192	ELSE
193	ZTSTAR=0.5*(ZTX+PTSTAR)
194	ENDIF
195
196	ZT0=ZTSTAR+ZDTDZSG*POROG
197	IF (ZTX > ZTSTAR .AND. ZT0 > ZTX) THEN
198	ZT0=ZTX
199	ELSEIF (ZTX <= ZTSTAR .AND. ZT0 > ZTSTAR) THEN
200	ZT0=ZTSTAR
201	ELSE
202	ZT0=PT0
203	ENDIF
204
205	ZOROG=SIGN(MAX(1.0,ABS(POROG)),POROG)
206	ZALPHA=r_d*(ZT0-ZTSTAR)/ZOROG
207
208	!* 2.2 COMPUTATION OF MSL PRESSURE.
209
210	IF (ABS(POROG) >= 0.001) THEN
211	ZX=POROG/(r_d*ZTSTAR)
212	ZY=ZALPHA*ZX
213	ZY2=ZY*ZY
214
215	psl=PRPRESSEXP(ZX(1.0-0.5ZY+1.0/3.ZY2))
216	ELSE
217	psl=PRPRESS
218	ENDIF
219
220	RETURN
221
222	END SUBROUTINE var_psl_ecmwf
223
224	SUBROUTINE compute_psl_ptarget4d2(press, ps, hgt, ta, qv, ptarget, psl, d1, d2, d3, d4)
225	! Subroutine to compute sea level pressure using a target pressure. Similar to the Benjamin
226	! and Miller (1990). Method found in p_interp.F90
227
228	IMPLICIT NONE
229
230	INTEGER, INTENT(in) :: d1, d2, d3, d4
231	REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in) :: press, ta, qv
232	REAL(r_k), DIMENSION(d1,d2,d4), INTENT(in) :: ps
233	REAL(r_k), DIMENSION(d1,d2), INTENT(in) :: hgt
234	REAL(r_k), INTENT(in) :: ptarget
235	REAL(r_k), DIMENSION(d1,d2,d4), INTENT(out) :: psl
236
237	! Local
238	INTEGER :: i, j, it
239	INTEGER :: kin
240	INTEGER :: kupper
241	REAL(r_k) :: dpmin, dp, tbotextrap, &
242	tvbotextrap, virtual
243	! Exponential related to standard atmosphere lapse rate r_d*gammav/g
244	REAL(r_k), PARAMETER :: expon=r_d*gammav/grav
245
246	!!!!!!! Variables
247	! press: Atmospheric pressure [Pa]
248	! ps: surface pressure [Pa]
249	! hgt: surface height
250	! ta: temperature [K]
251	! qv: water vapor mixing ratio
252	! dz: number of vertical levels
253	! psl: sea-level pressure
254
255	fname = 'compute_psl_ptarget4d2'
256
257	! Minimal distance between pressures [Pa]
258	dpmin=1.e4
259	psl=0.
260
261	DO i=1,d1
262	DO j=1,d2
263	IF (hgt(i,j) /= 0.) THEN
264	DO it=1,d4
265
266	! target pressure to be used for the extrapolation [Pa] (defined in namelist.input)
267	! ptarget = 70000. default value
268
269	! We are below both the ground and the lowest data level.
270
271	! First, find the model level that is closest to a "target" pressure
272	! level, where the "target" pressure is delta-p less that the local
273	! value of a horizontally smoothed surface pressure field. We use
274	! delta-p = 150 hPa here. A standard lapse rate temperature profile
275	! passing through the temperature at this model level will be used
276	! to define the temperature profile below ground. This is similar
277	! to the Benjamin and Miller (1990) method, using
278	! 700 hPa everywhere for the "target" pressure.
279
280	kupper = 0
281	loop_kIN: DO kin=d3,1,-1
282	kupper = kin
283	dp=abs( press(i,j,kin,it) - ptarget )
284	IF (dp .GT. dpmin) EXIT loop_kIN
285	dpmin=min(dpmin,dp)
286	ENDDO loop_kIN
287
288	tbotextrap=ta(i,j,kupper,it)(ps(i,j,it)/ptarget)*expon
289	tvbotextrap=virtualTemp1D(tbotextrap,qv(i,j,kupper,it))
290
291	psl(i,j,it) = ps(i,j,it)((tvbotextrap+gammavhgt(i,j))/tvbotextrap)**(1/expon)
292	END DO
293	ELSE
294	psl(i,j,:) = ps(i,j,:)
295	END IF
296	END DO
297	END DO
298
299	RETURN
300
301	END SUBROUTINE compute_psl_ptarget4d2
302
303	SUBROUTINE compute_tv4d(ta,qv,tv,d1,d2,d3,d4)
304	! 4D calculation of virtual temperaure
305
306	IMPLICIT NONE
307
308	INTEGER, INTENT(in) :: d1, d2, d3, d4
309	REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in) :: ta, qv
310	REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(out) :: tv
311
312	! Variables
313	! ta: temperature [K]
314	! qv: mixing ratio [kgkg-1]
315	! tv: virtual temperature
316
317	tv = ta*(oneRK+(qv/epsilonv))/(oneRK+qv)
318
319	END SUBROUTINE compute_tv4d
320
321	FUNCTION VirtualTemp1D (ta,qv) result (tv)
322	! 1D calculation of virtual temperaure
323
324	IMPLICIT NONE
325
326	REAL(r_k), INTENT(in) :: ta, qv
327	REAL(r_k) :: tv
328
329	! Variables
330	! ta: temperature [K]
331	! qv: mixing ratio [kgkg-1]
332
333	tv = ta*(oneRK+(qv/epsilonv))/(oneRK+qv)
334
335	END FUNCTION VirtualTemp1D
336
337	! ---- BEGIN modified from module_diag_afwa.F ---- !
338
339	!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
340	!~
341	!~ Name:
342	!~ Theta
343	!~
344	!~ Description:
345	!~ This function calculates potential temperature as defined by
346	!~ Poisson's equation, given temperature and pressure ( hPa ).
347	!~
348	!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
349	FUNCTION Theta ( t, p )
350
351	IMPLICIT NONE
352
353	!~ Variable declaration
354	! --------------------
355	REAL(r_k), INTENT ( IN ) :: t
356	REAL(r_k), INTENT ( IN ) :: p
357	REAL(r_k) :: theta
358
359	! Using WRF values
360	!REAL :: Rd ! Dry gas constant
361	!REAL :: Cp ! Specific heat of dry air at constant pressure
362	!REAL :: p00 ! Standard pressure ( 1000 hPa )
363	REAL(r_k) :: Rd, p00
364
365	!Rd = 287.04
366	!Cp = 1004.67
367	!p00 = 1000.00
368
369	Rd = r_d
370	p00 = p1000mb/100.
371
372	!~ Poisson's equation
373	! ------------------
374	theta = t * ( (p00/p)**(Rd/Cp) )
375
376	END FUNCTION Theta
377
378	!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
379	!~
380	!~ Name:
381	!~ Thetae
382	!~
383	!~ Description:
384	!~ This function returns equivalent potential temperature using the
385	!~ method described in Bolton 1980, Monthly Weather Review, equation 43.
386	!~
387	!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
388	FUNCTION Thetae ( tK, p, rh, mixr )
389
390	IMPLICIT NONE
391
392	!~ Variable Declarations
393	! ---------------------
394	REAL(r_k) :: tK ! Temperature ( K )
395	REAL(r_k) :: p ! Pressure ( hPa )
396	REAL(r_k) :: rh ! Relative humidity
397	REAL(r_k) :: mixr ! Mixing Ratio ( kg kg^-1)
398	REAL(r_k) :: te ! Equivalent temperature ( K )
399	REAL(r_k) :: thetae ! Equivalent potential temperature
400
401	! Using WRF values
402	!REAL, PARAMETER :: R = 287.04 ! Universal gas constant (J/deg kg)
403	!REAL, PARAMETER :: P0 = 1000.0 ! Standard pressure at surface (hPa)
404	REAL(r_k) :: R, p00, Lv
405	!REAL, PARAMETER :: lv = 2.54(10*6) ! Latent heat of vaporization
406	! (J kg^-1)
407	!REAL, PARAMETER :: cp = 1004.67 ! Specific heat of dry air constant
408	! at pressure (J/deg kg)
409	REAL(r_k) :: tlc ! LCL temperature
410
411	R = r_d
412	p00 = p1000mb/100.
413	lv = XLV
414
415	!~ Calculate the temperature of the LCL
416	! ------------------------------------
417	tlc = TLCL ( tK, rh )
418
419	!~ Calculate theta-e
420	! -----------------
421	thetae = (tK * (p00/p)*( (R/Cp)(1.- ( (.28E-3)mixr1000.) ) ) )* &
422	exp( (((3.376/tlc)-.00254))*&
423	(mixr1000.(1.+(.81E-3)mixr1000.)) )
424
425	END FUNCTION Thetae
426
427	!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
428	!~
429	!~ Name:
430	!~ The2T.f90
431	!~
432	!~ Description:
433	!~ This function returns the temperature at any pressure level along a
434	!~ saturation adiabat by iteratively solving for it from the parcel
435	!~ thetae.
436	!~
437	!~ Dependencies:
438	!~ function thetae.f90
439	!~
440	!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
441	FUNCTION The2T ( thetaeK, pres, flag ) result ( tparcel )
442
443	IMPLICIT NONE
444
445	!~ Variable Declaration
446	! --------------------
447	REAL(r_k), INTENT ( IN ) :: thetaeK
448	REAL(r_k), INTENT ( IN ) :: pres
449	LOGICAL, INTENT ( INOUT ) :: flag
450	REAL(r_k) :: tparcel
451
452	REAL(r_k) :: thetaK
453	REAL(r_k) :: tovtheta
454	REAL(r_k) :: tcheck
455	REAL(r_k) :: svpr, svpr2
456	REAL(r_k) :: smixr, smixr2
457	REAL(r_k) :: thetae_check, thetae_check2
458	REAL(r_k) :: tguess_2, correction
459
460	LOGICAL :: found
461	INTEGER :: iter
462
463	! Using WRF values
464	!REAL :: R ! Dry gas constant
465	!REAL :: Cp ! Specific heat for dry air
466	!REAL :: kappa ! Rd / Cp
467	!REAL :: Lv ! Latent heat of vaporization at 0 deg. C
468	REAL(r_k) :: R, kappa, Lv
469
470	R = r_d
471	Lv = XLV
472	!R = 287.04
473	!Cp = 1004.67
474	Kappa = R/Cp
475	!Lv = 2.500E+6
476
477	!~ Make initial guess for temperature of the parcel
478	! ------------------------------------------------
479	tovtheta = (pres/100000.0)**(r/cp)
480	tparcel = thetaeK/exp(lv.012/(cp295.))*tovtheta
481
482	iter = 1
483	found = .false.
484	flag = .false.
485
486	DO
487	IF ( iter > 105 ) EXIT
488
489	tguess_2 = tparcel + REAL ( 1 )
490
491	svpr = 6.122 * exp ( (17.67*(tparcel-273.15)) / (tparcel-29.66) )
492	smixr = ( 0.622*svpr ) / ( (pres/100.0)-svpr )
493	svpr2 = 6.122 * exp ( (17.67*(tguess_2-273.15)) / (tguess_2-29.66) )
494	smixr2 = ( 0.622*svpr2 ) / ( (pres/100.0)-svpr2 )
495
496	! ------------------------------------------------------------------ ~!
497	!~ When this function was orinially written, the final parcel ~!
498	!~ temperature check was based off of the parcel temperature and ~!
499	!~ not the theta-e it produced. As there are multiple temperature- ~!
500	!~ mixing ratio combinations that can produce a single theta-e value, ~!
501	!~ we change the check to be based off of the resultant theta-e ~!
502	!~ value. This seems to be the most accurate way of backing out ~!
503	!~ temperature from theta-e. ~!
504	!~ ~!
505	!~ Rentschler, April 2010 ~!
506	! ------------------------------------------------------------------ !
507
508	!~ Old way...
509	!thetaK = thetaeK / EXP (lv * smixr /(cp*tparcel) )
510	!tcheck = thetaK * tovtheta
511
512	!~ New way
513	thetae_check = Thetae ( tparcel, pres/100., 100., smixr )
514	thetae_check2 = Thetae ( tguess_2, pres/100., 100., smixr2 )
515
516	!~ Whew doggies - that there is some accuracy...
517	!IF ( ABS (tparcel-tcheck) < .05) THEN
518	IF ( ABS (thetaeK-thetae_check) < .001) THEN
519	found = .true.
520	flag = .true.
521	EXIT
522	END IF
523
524	!~ Old
525	!tparcel = tparcel + (tcheck - tparcel)*.3
526
527	!~ New
528	correction = ( thetaeK-thetae_check ) / ( thetae_check2-thetae_check )
529	tparcel = tparcel + correction
530
531	iter = iter + 1
532	END DO
533
534	!IF ( .not. found ) THEN
535	! print*, "Warning! Thetae to temperature calculation did not converge!"
536	! print*, "Thetae ", thetaeK, "Pressure ", pres
537	!END IF
538
539	END FUNCTION The2T
540
541	!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
542	!~
543	!~ Name:
544	!~ VirtualTemperature
545	!~
546	!~ Description:
547	!~ This function returns virtual temperature given temperature ( K )
548	!~ and mixing ratio.
549	!~
550	!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
551	FUNCTION VirtualTemperature ( tK, w ) result ( Tv )
552
553	IMPLICIT NONE
554
555	!~ Variable declaration
556	real(r_k), intent ( in ) :: tK !~ Temperature
557	real(r_k), intent ( in ) :: w !~ Mixing ratio ( kg kg^-1 )
558	real(r_k) :: Tv !~ Virtual temperature
559
560	Tv = tK * ( 1.0 + (w/0.622) ) / ( 1.0 + w )
561
562	END FUNCTION VirtualTemperature
563
564	!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
565	!~
566	!~ Name:
567	!~ SaturationMixingRatio
568	!~
569	!~ Description:
570	!~ This function calculates saturation mixing ratio given the
571	!~ temperature ( K ) and the ambient pressure ( Pa ). Uses
572	!~ approximation of saturation vapor pressure.
573	!~
574	!~ References:
575	!~ Bolton (1980), Monthly Weather Review, pg. 1047, Eq. 10
576	!~
577	!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
578	FUNCTION SaturationMixingRatio ( tK, p ) result ( ws )
579
580	IMPLICIT NONE
581
582	REAL(r_k), INTENT ( IN ) :: tK
583	REAL(r_k), INTENT ( IN ) :: p
584	REAL(r_k) :: ws
585
586	REAL(r_k) :: es
587
588	es = 6.122 * exp ( (17.67*(tK-273.15))/ (tK-29.66) )
589	ws = ( 0.622*es ) / ( (p/100.0)-es )
590
591	END FUNCTION SaturationMixingRatio
592
593	!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
594	!~
595	!~ Name:
596	!~ tlcl
597	!~
598	!~ Description:
599	!~ This function calculates the temperature of a parcel of air would have
600	!~ if lifed dry adiabatically to it's lifting condensation level (lcl).
601	!~
602	!~ References:
603	!~ Bolton (1980), Monthly Weather Review, pg. 1048, Eq. 22
604	!~
605	!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
606	FUNCTION TLCL ( tk, rh )
607
608	IMPLICIT NONE
609
610	REAL(r_k), INTENT ( IN ) :: tK !~ Temperature ( K )
611	REAL(r_k), INTENT ( IN ) :: rh !~ Relative Humidity ( % )
612	REAL(r_k) :: tlcl
613
614	REAL(r_k) :: denom, term1, term2
615
616	term1 = 1.0 / ( tK - 55.0 )
617	!! Lluis
618	! IF ( rh > REAL (0) ) THEN
619	IF ( rh > zeroRK ) THEN
620	term2 = ( LOG (rh/100.0) / 2840.0 )
621	ELSE
622	term2 = ( LOG (0.001/oneRK) / 2840.0 )
623	END IF
624	denom = term1 - term2
625	!! Lluis
626	! tlcl = ( 1.0 / denom ) + REAL ( 55 )
627	tlcl = ( oneRK / denom ) + 55*oneRK
628
629	END FUNCTION TLCL
630
631	FUNCTION var_cape_afwa1D(nz, tk, rhv, p, hgt, sfc, cape, cin, zlfc, plfc, lidx, parcel) RESULT (ostat)
632	! Function to compute cape on a 1D column following implementation in phys/module_diag_afwa.F
633
634	IMPLICIT NONE
635
636	INTEGER, INTENT(in) :: nz, sfc
637	REAL(r_k), DIMENSION(nz), INTENT(in) :: tk, rhv, p, hgt
638	REAL(r_k), INTENT(out) :: cape, cin, zlfc, plfc, lidx
639	INTEGER :: ostat
640	INTEGER, INTENT(in) :: parcel
641
642	! Local
643	!~ Derived profile variables
644	! -------------------------
645	REAL(r_k), DIMENSION(nz) :: rh, ws, w, dTvK, buoy
646	REAL(r_k) :: tlclK, plcl, nbuoy, pbuoy
647
648	!~ Source parcel information
649	! -------------------------
650	REAL(r_k) :: srctK, srcrh, srcws, srcw, srcp, &
651	srctheta, srcthetaeK
652	INTEGER :: srclev
653	REAL(r_k) :: spdiff
654
655	!~ Parcel variables
656	! ----------------
657	REAL(r_k) :: ptK, ptvK, tvK, pw
658
659	!~ Other utility variables
660	! -----------------------
661	INTEGER :: i, j, k
662	INTEGER :: lfclev
663	INTEGER :: prcl
664	INTEGER :: mlev
665	INTEGER :: lyrcnt
666	LOGICAL :: flag
667	LOGICAL :: wflag
668	REAL(r_k) :: freeze
669	REAL(r_k) :: pdiff
670	REAL(r_k) :: pm, pu, pd
671	REAL(r_k) :: lidxu
672	REAL(r_k) :: lidxd
673
674	REAL(r_k), PARAMETER :: Rd = r_d
675	REAL(r_k), PARAMETER :: RUNDEF = -9.999E30
676
677	!!!!!!! Variables
678	! nz: Number of vertical levels
679	! sfc: Surface level in the profile
680	! tk: Temperature profile [K]
681	! rhv: Relative Humidity profile [1]
682	! rh: Relative Humidity profile [%]
683	! p: Pressure profile [Pa]
684	! hgt: Geopotential height profile [gpm]
685	! cape: CAPE [Jkg-1]
686	! cin: CIN [Jkg-1]
687	! zlfc: LFC Height [gpm]
688	! plfc: LFC Pressure [Pa]
689	! lidx: Lifted index
690	! FROM: https://en.wikipedia.org/wiki/Lifted_index
691	! lidx >= 6: Very Stable Conditions
692	! 6 > lidx > 1: Stable Conditions, Thunderstorms Not Likely
693	! 0 > lidx > -2: Slightly Unstable, Thunderstorms Possible, With Lifting Mechanism (i.e., cold front, daytime heating, ...)
694	! -2 > lidx > -6: Unstable, Thunderstorms Likely, Some Severe With Lifting Mechanism
695	! -6 > lidx: Very Unstable, Severe Thunderstorms Likely With Lifting Mechanism
696	! ostat: Function return status (Nonzero is bad)
697	! parcel:
698	! Most Unstable = 1 (default)
699	! Mean layer = 2
700	! Surface based = 3
701	!~ Derived profile variables
702	! -------------------------
703	! ws: Saturation mixing ratio
704	! w: Mixing ratio
705	! dTvK: Parcel / ambient Tv difference
706	! buoy: Buoyancy
707	! tlclK: LCL temperature [K]
708	! plcl: LCL pressure [Pa]
709	! nbuoy: Negative buoyancy
710	! pbuoy: Positive buoyancy
711
712	!~ Source parcel information
713	! -------------------------
714	! srctK: Source parcel temperature [K]
715	! srcrh: Source parcel rh [%]
716	! srcws: Source parcel sat. mixing ratio
717	! srcw: Source parcel mixing ratio
718	! srcp: Source parcel pressure [Pa]
719	! srctheta: Source parcel theta [K]
720	! srcthetaeK: Source parcel theta-e [K]
721	! srclev: Level of the source parcel
722	! spdiff: Pressure difference
723
724	!~ Parcel variables
725	! ----------------
726	! ptK: Parcel temperature [K]
727	! ptvK: Parcel virtual temperature [K]
728	! tvK: Ambient virtual temperature [K]
729	! pw: Parcel mixing ratio
730
731	!~ Other utility variables
732	! -----------------------
733	! lfclev: Level of LFC
734	! prcl: Internal parcel type indicator
735	! mlev: Level for ML calculation
736	! lyrcnt: Number of layers in mean layer
737	! flag: Dummy flag
738	! wflag: Saturation flag
739	! freeze: Water loading multiplier
740	! pdiff: Pressure difference between levs
741	! pm, pu, pd: Middle, upper, lower pressures
742	! lidxu: Lifted index at upper level
743	! lidxd: Lifted index at lower level
744
745	fname = 'var_cape_afwa'
746
747	!~ Initialize variables
748	! --------------------
749	rh = rhv*100.
750	ostat = 0
751	CAPE = zeroRK
752	CIN = zeroRK
753	ZLFC = RUNDEF
754	PLFC = RUNDEF
755
756	!~ Look for submitted parcel definition
757	!~ 1 = Most unstable
758	!~ 2 = Mean layer
759	!~ 3 = Surface based
760	! -------------------------------------
761	IF ( parcel > 3 .or. parcel < 1 ) THEN
762	prcl = 1
763	ELSE
764	prcl = parcel
765	END IF
766
767	!~ Initalize our parcel to be (sort of) surface based. Because of
768	!~ issues we've been observing in the WRF model, specifically with
769	!~ excessive surface moisture values at the surface, using a true
770	!~ surface based parcel is resulting a more unstable environment
771	!~ than is actually occuring. To address this, our surface parcel
772	!~ is now going to be defined as the parcel between 25-50 hPa
773	!~ above the surface. UPDATE - now that this routine is in WRF,
774	!~ going to trust surface info. GAC 20140415
775	! ----------------------------------------------------------------
776
777	!~ Compute mixing ratio values for the layer
778	! -----------------------------------------
779	DO k = sfc, nz
780	ws ( k ) = SaturationMixingRatio ( tK(k), p(k) )
781	w ( k ) = ( rh(k)/100.0 ) * ws ( k )
782	END DO
783
784	srclev = sfc
785	srctK = tK ( sfc )
786	srcrh = rh ( sfc )
787	srcp = p ( sfc )
788	srcws = ws ( sfc )
789	srcw = w ( sfc )
790	srctheta = Theta ( tK(sfc), p(sfc)/100.0 )
791
792	!~ Compute the profile mixing ratio. If the parcel is the MU parcel,
793	!~ define our parcel to be the most unstable parcel within the lowest
794	!~ 180 mb.
795	! -------------------------------------------------------------------
796	mlev = sfc + 1
797	DO k = sfc + 1, nz
798
799	!~ Identify the last layer within 100 hPa of the surface
800	! -----------------------------------------------------
801	pdiff = ( p (sfc) - p (k) ) / REAL ( 100 )
802	IF ( pdiff <= REAL (100) ) mlev = k
803
804	!~ If we've made it past the lowest 180 hPa, exit the loop
805	! -------------------------------------------------------
806	IF ( pdiff >= REAL (180) ) EXIT
807
808	IF ( prcl == 1 ) THEN
809	!IF ( (p(k) > 70000.0) .and. (w(k) > srcw) ) THEN
810	IF ( (w(k) > srcw) ) THEN
811	srctheta = Theta ( tK(k), p(k)/100.0 )
812	srcw = w ( k )
813	srclev = k
814	srctK = tK ( k )
815	srcrh = rh ( k )
816	srcp = p ( k )
817	END IF
818	END IF
819
820	END DO
821
822	!~ If we want the mean layer parcel, compute the mean values in the
823	!~ lowest 100 hPa.
824	! ----------------------------------------------------------------
825	lyrcnt = mlev - sfc + 1
826	IF ( prcl == 2 ) THEN
827
828	srclev = sfc
829	srctK = SUM ( tK (sfc:mlev) ) / REAL ( lyrcnt )
830	srcw = SUM ( w (sfc:mlev) ) / REAL ( lyrcnt )
831	srcrh = SUM ( rh (sfc:mlev) ) / REAL ( lyrcnt )
832	srcp = SUM ( p (sfc:mlev) ) / REAL ( lyrcnt )
833	srctheta = Theta ( srctK, srcp/100. )
834
835	END IF
836
837	srcthetaeK = Thetae ( srctK, srcp/100.0, srcrh, srcw )
838
839	!~ Calculate temperature and pressure of the LCL
840	! ---------------------------------------------
841	tlclK = TLCL ( tK(srclev), rh(srclev) )
842	plcl = p(srclev) * ( (tlclK/tK(srclev))**(Cp/Rd) )
843
844	!~ Now lift the parcel
845	! -------------------
846
847	buoy = REAL ( 0 )
848	pw = srcw
849	wflag = .false.
850	DO k = srclev, nz
851	IF ( p (k) <= plcl ) THEN
852
853	!~ The first level after we pass the LCL, we're still going to
854	!~ lift the parcel dry adiabatically, as we haven't added the
855	!~ the required code to switch between the dry adiabatic and moist
856	!~ adiabatic cooling. Since the dry version results in a greater
857	!~ temperature loss, doing that for the first step so we don't over
858	!~ guesstimate the instability.
859	! ----------------------------------------------------------------
860
861	IF ( wflag ) THEN
862	flag = .false.
863
864	!~ Above the LCL, our parcel is now undergoing moist adiabatic
865	!~ cooling. Because of the latent heating being undergone as
866	!~ the parcel rises above the LFC, must iterative solve for the
867	!~ parcel temperature using equivalant potential temperature,
868	!~ which is conserved during both dry adiabatic and
869	!~ pseudoadiabatic displacements.
870	! --------------------------------------------------------------
871	ptK = The2T ( srcthetaeK, p(k), flag )
872
873	!~ Calculate the parcel mixing ratio, which is now changing
874	!~ as we condense moisture out of the parcel, and is equivalent
875	!~ to the saturation mixing ratio, since we are, in theory, at
876	!~ saturation.
877	! ------------------------------------------------------------
878	pw = SaturationMixingRatio ( ptK, p(k) )
879
880	!~ Now we can calculate the virtual temperature of the parcel
881	!~ and the surrounding environment to assess the buoyancy.
882	! ----------------------------------------------------------
883	ptvK = VirtualTemperature ( ptK, pw )
884	tvK = VirtualTemperature ( tK (k), w (k) )
885
886	!~ Modification to account for water loading
887	! -----------------------------------------
888	freeze = 0.033 * ( 263.15 - pTvK )
889	IF ( freeze > 1.0 ) freeze = 1.0
890	IF ( freeze < 0.0 ) freeze = 0.0
891
892	!~ Approximate how much of the water vapor has condensed out
893	!~ of the parcel at this level
894	! ---------------------------------------------------------
895	freeze = freeze * 333700.0 * ( srcw - pw ) / 1005.7
896
897	pTvK = pTvK - pTvK * ( srcw - pw ) + freeze
898	dTvK ( k ) = ptvK - tvK
899	buoy ( k ) = g * ( dTvK ( k ) / tvK )
900
901	ELSE
902
903	!~ Since the theta remains constant whilst undergoing dry
904	!~ adiabatic processes, can back out the parcel temperature
905	!~ from potential temperature below the LCL
906	! --------------------------------------------------------
907	ptK = srctheta / ( 100000.0/p(k) )**(Rd/Cp)
908
909	!~ Grab the parcel virtual temperture, can use the source
910	!~ mixing ratio since we are undergoing dry adiabatic cooling
911	! ----------------------------------------------------------
912	ptvK = VirtualTemperature ( ptK, srcw )
913
914	!~ Virtual temperature of the environment
915	! --------------------------------------
916	tvK = VirtualTemperature ( tK (k), w (k) )
917
918	!~ Buoyancy at this level
919	! ----------------------
920	dTvK ( k ) = ptvK - tvK
921	buoy ( k ) = g * ( dtvK ( k ) / tvK )
922
923	wflag = .true.
924
925	END IF
926
927	ELSE
928
929	!~ Since the theta remains constant whilst undergoing dry
930	!~ adiabatic processes, can back out the parcel temperature
931	!~ from potential temperature below the LCL
932	! --------------------------------------------------------
933	ptK = srctheta / ( 100000.0/p(k) )**(Rd/Cp)
934
935	!~ Grab the parcel virtual temperture, can use the source
936	!~ mixing ratio since we are undergoing dry adiabatic cooling
937	! ----------------------------------------------------------
938	ptvK = VirtualTemperature ( ptK, srcw )
939
940	!~ Virtual temperature of the environment
941	! --------------------------------------
942	tvK = VirtualTemperature ( tK (k), w (k) )
943
944	!~ Buoyancy at this level
945	! ---------------------
946	dTvK ( k ) = ptvK - tvK
947	buoy ( k ) = g * ( dtvK ( k ) / tvK )
948
949	END IF
950
951	!~ Chirp
952	! -----
953	! WRITE ( ,'(I15,6F15.3)' )k,p(k)/100.,ptK,pw1000.,ptvK,tvK,buoy(k)
954
955	END DO
956
957	!~ Add up the buoyancies, find the LFC
958	! -----------------------------------
959	flag = .false.
960	lfclev = -1
961	nbuoy = REAL ( 0 )
962	pbuoy = REAL ( 0 )
963	DO k = sfc + 1, nz
964	IF ( tK (k) < 253.15 ) EXIT
965	CAPE = CAPE + MAX ( buoy (k), 0.0 ) * ( hgt (k) - hgt (k-1) )
966	CIN = CIN + MIN ( buoy (k), 0.0 ) * ( hgt (k) - hgt (k-1) )
967
968	!~ If we've already passed the LFC
969	! -------------------------------
970	IF ( flag .and. buoy (k) > REAL (0) ) THEN
971	pbuoy = pbuoy + buoy (k)
972	END IF
973
974	!~ We are buoyant now - passed the LFC
975	! -----------------------------------
976	IF ( .not. flag .and. buoy (k) > REAL (0) .and. p (k) < plcl ) THEN
977	flag = .true.
978	pbuoy = pbuoy + buoy (k)
979	lfclev = k
980	END IF
981
982	!~ If we think we've passed the LFC, but encounter a negative layer
983	!~ start adding it up.
984	! ----------------------------------------------------------------
985	IF ( flag .and. buoy (k) < REAL (0) ) THEN
986	nbuoy = nbuoy + buoy (k)
987
988	!~ If the accumulated negative buoyancy is greater than the
989	!~ positive buoyancy, then we are capped off. Got to go higher
990	!~ to find the LFC. Reset positive and negative buoyancy summations
991	! ----------------------------------------------------------------
992	IF ( ABS (nbuoy) > pbuoy ) THEN
993	flag = .false.
994	nbuoy = REAL ( 0 )
995	pbuoy = REAL ( 0 )
996	lfclev = -1
997	END IF
998	END IF
999
1000	END DO
1001
1002	!~ Calculate lifted index by interpolating difference between
1003	!~ parcel and ambient Tv to 500mb.
1004	! ----------------------------------------------------------
1005	DO k = sfc + 1, nz
1006
1007	pm = 50000.
1008	pu = p ( k )
1009	pd = p ( k - 1 )
1010
1011	!~ If we're already above 500mb just set lifted index to 0.
1012	!~ --------------------------------------------------------
1013	IF ( pd .le. pm ) THEN
1014	lidx = zeroRK
1015	EXIT
1016
1017	ELSEIF ( pu .le. pm .and. pd .gt. pm) THEN
1018
1019	!~ Found trapping pressure: up, middle, down.
1020	!~ We are doing first order interpolation.
1021	! ------------------------------------------
1022	lidxu = -dTvK ( k ) * ( pu / 100000. ) ** (Rd/Cp)
1023	lidxd = -dTvK ( k-1 ) * ( pd / 100000. ) ** (Rd/Cp)
1024	lidx = ( lidxu * (pm-pd) + lidxd * (pu-pm) ) / (pu-pd)
1025	EXIT
1026
1027	ENDIF
1028
1029	END DO
1030
1031	!~ Assuming the the LFC is at a pressure level for now
1032	! ---------------------------------------------------
1033	IF ( lfclev > zeroRK ) THEN
1034	PLFC = p ( lfclev )
1035	ZLFC = hgt ( lfclev )
1036	END IF
1037
1038	IF ( PLFC /= PLFC .OR. PLFC < zeroRK ) THEN
1039	PLFC = -oneRK
1040	ZLFC = -oneRK
1041	END IF
1042
1043	IF ( CAPE /= CAPE ) cape = zeroRK
1044
1045	IF ( CIN /= CIN ) cin = zeroRK
1046
1047	!~ Chirp
1048	! -----
1049	! WRITE ( , ) ' CAPE: ', cape, ' CIN: ', cin
1050	! WRITE ( , ) ' LFC: ', ZLFC, ' PLFC: ', PLFC
1051	! WRITE ( , ) ''
1052	! WRITE ( , ) ' Exiting buoyancy.'
1053	! WRITE ( , ) ' ==================================== '
1054	! WRITE ( , ) ''
1055
1056	RETURN
1057
1058	END FUNCTION var_cape_afwa1D
1059
1060	! ---- END modified from module_diag_afwa.F ---- !
1061
1062	SUBROUTINE var_zmla_generic(dz, qv, tpot, z, topo, zmla)
1063	! Subroutine to compute pbl-height following a generic method
1064	! from Nielsen-Gammon et al., 2008 J. Appl. Meteor. Clim.
1065	! applied also in Garcia-Diez et al., 2013, QJRMS
1066	! where
1067	! "The technique identifies the ML height as a threshold increase of potential temperature from
1068	! its minimum value within the boundary layer."
1069	! here applied similarly to Garcia-Diez et al. where
1070	! zmla = "...first level where potential temperature exceeds the minimum potential temperature
1071	! reached in the mixed layer by more than 1.5 K"
1072
1073	IMPLICIT NONE
1074
1075	INTEGER, INTENT(in) :: dz
1076	REAL(r_k), DIMENSION(dz), INTENT(in) :: qv, tpot, z
1077	REAL(r_k), INTENT(in) :: topo
1078	REAL(r_k), INTENT(out) :: zmla
1079
1080	! Local
1081	INTEGER :: i
1082	INTEGER :: mldlev, bllev
1083	REAL(r_k) :: dqvar, tpotmin, refdt
1084
1085	!!!!!!! Variables
1086	! qv: water vapour mixing ratio
1087	! tpot: potential temperature [K]
1088	! z: height above sea level [m]
1089	! topo: topographic height [m]
1090	! zmla: boundary layer height [m]
1091
1092	fname = 'var_zmla_generic'
1093
1094	! Pecentage of difference of mixing ratio used to determine Mixed layer depth
1095	dqvar = 0.1
1096
1097	! MLD = Mixed layer with no substantial variation of mixing ratio /\qv < 10% ?
1098	!PRINT *,' Mixed layer mixing ratios qv[1] lev qv[lev] dqvar% _______'
1099	DO mldlev = 2, dz
1100	IF (ABS(qv(mldlev)-qv(1))/qv(1) > dqvar ) EXIT
1101	! PRINT *,qv(1), mldlev, qv(mldlev), ABS(qv(mldlev)-qv(1))/qv(1)
1102	END DO
1103
1104	! Looking for the minimum potential temperature within the MLD [tpotmin = min(tpot)_0^MLD]
1105	tpotmin = MINVAL(tpot(1:mldlev))
1106
1107	! Change in temperature to determine boundary layer height
1108	refdt = 1.5
1109
1110	! Determine the first level where tpot > tpotmin + 1.5 K
1111	!PRINT *,' Mixed layer tpotmin lev tpotmin[lev] dtpot _______'
1112	DO bllev = 1, dz
1113	IF (ABS(tpot(bllev)-tpotmin) > refdt ) EXIT
1114	! PRINT *,tpotmin, bllev, tpot(bllev), ABS(tpot(bllev)-tpotmin)
1115	END DO
1116
1117	!PRINT *,' height end MLD:', z(mldlev)
1118	!PRINT *,' pbl height:', z(bllev)
1119
1120	zmla = z(bllev) - topo
1121
1122	RETURN
1123
1124	END SUBROUTINE var_zmla_generic
1125
1126	SUBROUTINE var_zwind(d1, u, v, z, u10, v10, sa, ca, newz, unewz, vnewz)
1127	! Subroutine to extrapolate the wind at a given height following the 'power law' methodology
1128	! wss[newz] = wss[z1](newz/z1)*alpha
1129	! alpha = (ln(wss[z2])-ln(wss[z1]))/(ln(z2)-ln(z1))
1130	! AFTER: Phd Thesis:
1131	! Benedicte Jourdier. Ressource eolienne en France metropolitaine : methodes dâevaluation du
1132	! potentiel, variabilite et tendances. Climatologie. Ecole Doctorale Polytechnique, 2015. French
1133	!
1134	IMPLICIT NONE
1135
1136	INTEGER, INTENT(in) :: d1
1137	REAL(r_k), DIMENSION(d1), INTENT(in) :: u,v,z
1138	REAL(r_k), INTENT(in) :: u10, v10, sa, ca, newz
1139	REAL(r_k), INTENT(out) :: unewz, vnewz
1140
1141	! Local
1142	INTEGER :: inear
1143	REAL(r_k) :: zaground
1144	REAL(r_k), DIMENSION(2) :: v1, v2, zz, alpha, uvnewz
1145
1146	!!!!!!! Variables
1147	! u,v: vertical wind components [ms-1]
1148	! z: height above surface on half-mass levels [m]
1149	! u10,v10: 10-m wind components [ms-1]
1150	! sa, ca: local sine and cosine of map rotation [1.]
1151	! newz: desired height above grpund of extrapolation
1152	! unewz,vnewz: Wind compoonents at the given height [ms-1]
1153
1154	fname = 'var_zwind'
1155
1156	!PRINT *,' ilev zaground newz z[ilev+1] z[ilev+2] _______'
1157	IF (z(1) < newz ) THEN
1158	DO inear = 1,d1-2
1159	! L. Fita, CIMA. Feb. 2018
1160	!! Choose between extra/inter-polate. Maybe better interpolate?
1161	! Here we extrapolate from two closest lower levels
1162	!zaground = z(inear+2)
1163	zaground = z(inear+1)
1164	!PRINT *, inear, z(inear), newz, z(inear+1), z(inear+2)
1165	IF ( zaground >= newz) EXIT
1166	END DO
1167	ELSE
1168	!PRINT *, 1, z(1), newz, z(2), z(3), ' z(1) > newz'
1169	inear = d1 - 2
1170	END IF
1171
1172	IF (inear == d1-2) THEN
1173	! No vertical pair of levels is below newz, using 10m wind as first value and the first level as the second
1174	v1(1) = u10
1175	v1(2) = v10
1176	v2(1) = u(1)
1177	v2(2) = v(1)
1178	zz(1) = 10.
1179	zz(2) = z(1)
1180	ELSE
1181	v1(1) = u(inear)
1182	v1(2) = v(inear)
1183	v2(1) = u(inear+1)
1184	v2(2) = v(inear+1)
1185	zz(1) = z(inear)
1186	zz(2) = z(inear+1)
1187	END IF
1188
1189	! Computing for each component
1190	alpha = (LOG(ABS(v2))-LOG(ABS(v1)))/(LOG(zz(2))-LOG(zz(1)))
1191	!PRINT *,' Computing with v1:', v1, ' ms-1 v2:', v2, ' ms-1'
1192	!PRINT *,' z1:', zz(1), 'm z2:', zz(2), ' m'
1193	!PRINT *,' alhpa u:', alpha(1), ' alpha 2:', alpha(2)
1194
1195	uvnewz = v1(newz/zz(1))*alpha
1196	! Earth-rotation
1197	unewz = uvnewz(1)ca - uvnewz(2)sa
1198	vnewz = uvnewz(1)sa + uvnewz(2)ca
1199
1200	!PRINT *,' result vz:', uvnewz
1201
1202	!STOP
1203
1204	RETURN
1205
1206	END SUBROUTINE var_zwind
1207
1208	SUBROUTINE var_zwind_log(d1, u, v, z, u10, v10, sa, ca, newz, unewz, vnewz)
1209	! Subroutine to extrapolate the wind at a given height following the 'logarithmic law' methodology
1210	! wsz = wss[z2]*(ln(newz)-ln(z0))(ln(z2)-ln(z0))
1211	! ln(z0) = (ws(z2)ln(z1)-ws(z1)ln(z2))/(ws(z2)-ws(z1))
1212	! AFTER: Phd Thesis:
1213	! Benedicte Jourdier. Ressource eolienne en France metropolitaine : methodes dâevaluation du
1214	! potentiel, variabilite et tendances. Climatologie. Ecole Doctorale Polytechnique, 2015. French
1215	!
1216	IMPLICIT NONE
1217
1218	INTEGER, INTENT(in) :: d1
1219	REAL(r_k), DIMENSION(d1), INTENT(in) :: u,v,z
1220	REAL(r_k), INTENT(in) :: u10, v10, sa, ca, newz
1221	REAL(r_k), INTENT(out) :: unewz, vnewz
1222
1223	! Local
1224	INTEGER :: inear
1225	REAL(r_k) :: zaground
1226	REAL(r_k), DIMENSION(2) :: v1, v2, zz, logz0, uvnewz
1227
1228	!!!!!!! Variables
1229	! u,v: vertical wind components [ms-1]
1230	! z: height above surface on half-mass levels [m]
1231	! u10,v10: 10-m wind components [ms-1]
1232	! sa, ca: local sine and cosine of map rotation [1.]
1233	! newz: desired height above grpund of extrapolation
1234	! unewz,vnewz: Wind compoonents at the given height [ms-1]
1235
1236	fname = 'var_zwind_log'
1237
1238	!PRINT *,' ilev zaground newz z[ilev+1] z[ilev+2] _______'
1239	IF (z(1) < newz ) THEN
1240	DO inear = 1,d1-2
1241	! L. Fita, CIMA. Feb. 2018
1242	!! Choose between extra/inter-polate. Maybe better interpolate?
1243	! Here we extrapolate from two closest lower levels
1244	!zaground = z(inear+2)
1245	zaground = z(inear+1)
1246	!PRINT *, inear, z(inear), newz, z(inear+1), z(inear+2)
1247	IF ( zaground >= newz) EXIT
1248	END DO
1249	ELSE
1250	!PRINT *, 1, z(1), newz, z(2), z(3), ' z(1) > newz'
1251	inear = d1 - 2
1252	END IF
1253
1254	IF (inear == d1-2) THEN
1255	! No vertical pair of levels is below newz, using 10m wind as first value and the first level as the second
1256	v1(1) = u10
1257	v1(2) = v10
1258	v2(1) = u(1)
1259	v2(2) = v(1)
1260	zz(1) = 10.
1261	zz(2) = z(1)
1262	ELSE
1263	v1(1) = u(inear)
1264	v1(2) = v(inear)
1265	v2(1) = u(inear+1)
1266	v2(2) = v(inear+1)
1267	zz(1) = z(inear)
1268	zz(2) = z(inear+1)
1269	END IF
1270
1271	! Computing for each component
1272	logz0 = (v2LOG(zz(1))-v1LOG(zz(2)))/(v2-v1)
1273
1274	uvnewz = v2*(LOG(newz)-logz0)/(LOG(zz(2))-logz0)
1275	! Earth-rotation
1276	unewz = uvnewz(1)ca - uvnewz(2)sa
1277	vnewz = uvnewz(1)sa + uvnewz(2)ca
1278
1279	!PRINT *,' result vz:', uvnewz
1280
1281	!STOP
1282
1283	RETURN
1284
1285	END SUBROUTINE var_zwind_log
1286
1287	SUBROUTINE var_zwind_MOtheor(ust, znt, rmol, u10, v10, sa, ca, newz, uznew, vznew)
1288	! Subroutine of wind extrapolation following Moin-Obukhov theory R. B. Stull, 1988,
1289	! Springer (p376-383)
1290	! Only usefull for newz < 80. m
1291	! Ackonwledgement: M. A. Jimenez, UIB
1292
1293	IMPLICIT NONE
1294
1295	REAL, INTENT(in) :: ust, znt, rmol, u10, v10, sa, ca
1296	REAL, INTENT(in) :: newz
1297	REAL, INTENT(out) :: uznew, vznew
1298
1299	! Local
1300	REAL :: OL
1301	REAL :: stability
1302	REAL :: wsz, alpha
1303	REAL, DIMENSION(2) :: uvnewz
1304
1305	!!!!!!! Variables
1306	! ust: u* in similarity theory [ms-1]
1307	! z0: roughness length [m]
1308	!!! L. Fita, CIMA. Feb. 2018
1309	!! NOT SURE if it should be z0 instead?
1310	! znt: thermal time-varying roughness length [m]
1311	! rmol: inverse of Obukhov length [m-1]
1312	! u10: x-component 10-m wind speed [ms-1]
1313	! v10: y-component 10-m wind speed [ms-1]
1314	! sa, ca: local sine and cosine of map rotation [1.]
1315	!
1316	fname = 'var_zwind_MOtheor'
1317
1318	! Obukhov Length (using the Boussinesq approximation giving Tv from t2)
1319	OL = 1/rmol
1320
1321	! Wind speed at desired height
1322	PRINT *,'ust:', ust, 'znt:', znt, 'OL:', OL
1323
1324	CALL stabfunc_businger(newz,OL,stability)
1325	PRINT *,' z/L:', newz/OL,' stabfunc:', stability, 'log:', LOG(newz/znt), 'log+stability:', LOG(newz/znt) + stability
1326	PRINT *,' ust/karman:', ust/karman
1327
1328	wsz = ust/karman*( LOG(newz/znt) + stability)
1329	PRINT *,' wsz:', wsz
1330
1331	! Without taking into account Ekcman pumping, etc... redistributed by components unsing 10-m wind
1332	! as reference...
1333	alpha = ATAN2(v10,u10)
1334	uvnewz(1) = wsz*COS(alpha)
1335	uvnewz(2) = wsz*SIN(alpha)
1336	PRINT *,' uvnewz:', uvnewz
1337
1338	! Earth-rotation
1339	uznew = uvnewz(1)ca - uvnewz(2)sa
1340	vznew = uvnewz(1)sa + uvnewz(2)ca
1341	PRINT *,' uznew:', uznew, ' vznew:', vznew
1342
1343	RETURN
1344
1345	END SUBROUTINE var_zwind_MOtheor
1346
1347	! L. Fita, CIMA. Feb. 2018
1348	! WRF seems to have problems with my functions, let'suse subroutine instead
1349	!REAL FUNCTION stabfunc_businger(z,L)
1350	SUBROUTINE stabfunc_businger(z,L,stabfunc_busingerv)
1351	! Fucntion of the stability function after Businger et al. (1971), JAS, 28(2), 181â189
1352
1353	IMPLICIT NONE
1354
1355	REAL, INTENT(in) :: z,L
1356	REAL, INTENT(out) :: stabfunc_busingerv
1357
1358	! Local
1359	REAL :: zL, X
1360
1361	!!!!!!! Variables
1362	! z: height [m]
1363	! L: Obukhov length [-]
1364
1365	fname = 'stabfunc_businger'
1366
1367	IF (L /= 0.) THEN
1368	zL = z/L
1369	ELSE
1370	! Neutral
1371	zL = 0.
1372	END IF
1373
1374	IF (zL > 0.) THEN
1375	! Stable case
1376	stabfunc_busingerv = 4.7*z/L
1377	ELSE IF (zL < 0.) THEN
1378	! unstable
1379	X = (1. - 15.z/L)*(0.25)
1380	!stabfunc_busingerv = -2.LOG((1.+X)/2.)-LOG((1.+X2)/2.)+2.ATAN(X)-piconst/2.
1381	stabfunc_busingerv = LOG( ((1.+X*2)/2.)((1.+X)/2.)*2)-2.ATAN(X)+piconst/2.
1382	ELSE
1383	stabfunc_busingerv = 0.
1384	END IF
1385
1386	RETURN
1387
1388	! END FUNCTION stabfunc_businger
1389	END SUBROUTINE stabfunc_businger
1390
1391	REAL(r_k) FUNCTION Cdrag_0(ust,uas,vas)
1392	! Fuction to compute a first order generic approximation of the drag coefficient as
1393	! CD = (ust/wss)**2
1394	! after, Garratt, J.R., 1992.: The Atmospheric Boundary Layer. Cambridge Univ. Press,
1395	! Cambridge, U.K., 316 pp
1396	! Ackonwledgement: M. A. Jimenez, UIB
1397	!
1398	IMPLICIT NONE
1399
1400	REAL(r_k), INTENT(in) :: ust, uas, vas
1401
1402	!!!!!!! Variables
1403	! ust: u* in similarity theory [ms-1]
1404	! uas, vas: x/y-components of wind at 10 m
1405
1406	fname = 'Cdrag_0'
1407
1408	Cdrag_0 = ust2/(uas2+vas**2)
1409
1410	END FUNCTION Cdrag_0
1411
1412	SUBROUTINE var_potevap_orPM(rho1, ust, uas, vas, tas, ps, qv1, potevap)
1413	! Subroutine to compute the potential evapotranspiration following Penman-Monteith formulation
1414	! implemented in ORCHIDEE
1415	! potevap = dtrho1qc*(q2sat-qv1)
1416
1417	IMPLICIT NONE
1418
1419	REAL(r_k), INTENT(in) :: rho1, ust, uas, vas, tas, ps, qv1
1420	REAL(r_k), INTENT(out) :: potevap
1421
1422	! Local
1423	REAL(r_k) :: q2sat, Cd, qc
1424
1425	!!!!!!! Variables
1426	! rho1: atsmophere density at the first layer [kgm-3]
1427	! ust: u* in similarity theory [ms-1]
1428	! uas, vas: x/y-components of 10-m wind [ms-1]
1429	! tas: 2-m atmosphere temperature [K]
1430	! ps: surface pressure [Pa]
1431	! qv1: 1st layer atmospheric mixing ratio [kgkg-1]
1432	! potevap: potential evapo transpiration [kgm-2s-1]
1433	fname = 'var_potevap_orPM'
1434
1435	! q2sat: Saturated air at 2m (can be assumed to be q2 == qsfc?)
1436	q2sat = SaturationMixingRatio(tas, ps)
1437
1438	! Cd: drag coeffiecient
1439	Cd = Cdrag_0(ust, uas, vas)
1440
1441	! qc: surface drag coefficient
1442	qc = SQRT(uas2 + vas2)*Cd
1443
1444	potevap = MAX(zeroRK, rho1qc(q2sat - qv1))
1445
1446	END SUBROUTINE var_potevap_orPM
1447
1448	SUBROUTINE var_fog_K84(qc, qi, fog, vis)
1449	! Computation of fog (vis < 1km) only computed where qcloud, qice /= 0.
1450	! And visibility following Kunkel, B. A., (1984): Parameterization of droplet terminal velocity and
1451	! extinction coefficient in fog models. J. Climate Appl. Meteor., 23, 34â41.
1452
1453	IMPLICIT NONE
1454
1455	REAL(r_k), INTENT(in) :: qc, qi
1456	INTEGER, INTENT(out) :: fog
1457	REAL(r_k), INTENT(out) :: vis
1458
1459	! Local
1460	REAL(r_k) :: visc, visi
1461
1462	!!!!!!! Variables
1463	! qc: cloud mixing ratio [kgkg-1]
1464	! qi, ice mixing ratio [kgkg-1]
1465	! fog: presence of fog (1: yes, 0: no)
1466	! vis: visibility within fog [km]
1467
1468	fname = 'var_fog_K84'
1469
1470	IF (qi > nullv .OR. qc > nullv) THEN
1471	visc = 100000.*oneRK
1472	visi = 100000.*oneRK
1473	! From: Gultepe, 2006, JAM, 45, 1469-1480
1474	IF (qc > nullv) visc = 0.027(qc1000.)**(-0.88)
1475	IF (qi > nullv) visi = 0.024(qi1000.)**(-1.0)
1476	vis = MINVAL((/visc, visi/))
1477	IF (vis <= oneRK) THEN
1478	fog = 1
1479	ELSE
1480	fog = 0
1481	vis = -oneRK
1482	END IF
1483	ELSE
1484	fog = 0
1485	vis = -oneRK
1486	END IF
1487
1488	END SUBROUTINE var_fog_K84
1489
1490	SUBROUTINE var_fog_RUC(qv, ta, pres, fog, vis)
1491	! Computation of fog (vis < 1km) only computed where qcloud, qice /= 0.
1492	! And visibility following RUC method Smirnova, T. G., S. G. Benjamin, and J. M. Brown, 2000: Case
1493	! study verification of RUC/MAPS fog and visibility forecasts. Preprints, 9 th Conference on
1494	! Aviation, Range, and Aerospace Meteorlogy, AMS, Orlando, FL, Sep. 2000. Paper#2.3, 6 pp.
1495
1496	IMPLICIT NONE
1497
1498	REAL(r_k), INTENT(in) :: qv, ta, pres
1499	INTEGER, INTENT(out) :: fog
1500	REAL(r_k), INTENT(out) :: vis
1501
1502	! Local
1503	REAL(r_k) :: rh
1504
1505	!!!!!!! Variables
1506	! qc: cloud mixing ratio [kgkg-1]
1507	! qi, ice mixing ratio [kgkg-1]
1508	! fog: presence of fog (1: yes, 0: no)
1509	! vis: visibility within fog [km]
1510
1511	fname = 'var_fog_RUC'
1512
1513	CALL var_hur(ta, pres, qv, rh)
1514	! Avoiding supersaturation
1515	rh = MINVAL((/1.,rh/))
1516
1517	IF (rh > 0.3) THEN
1518	! From: Gultepe, I., and G. Isaac, 2006: Visbility versus precipitation rate and relative
1519	! humidity. Preprints, 12th Cloud Physics Conf, Madison, WI, Amer. Meteor. Soc., P2.55.
1520	! [Available online at http://ams.confex.com/ams/Madison2006/techprogram/paper_l13177.htm]
1521	vis = 60.EXP(-2.5(rh*100.-15.)/80.)
1522	IF (vis <= oneRK) THEN
1523	fog = 1
1524	ELSE
1525	fog = 0
1526	vis = -oneRK
1527	END IF
1528	ELSE
1529	fog = 0
1530	vis = -oneRK
1531	END IF
1532
1533	END SUBROUTINE var_fog_RUC
1534
1535	SUBROUTINE var_fog_FRAML50(qv, ta, pres, fog, vis)
1536	! Computation of fog (vis < 1km)
1537	! And visibility following Gultepe, I. and J.A. Milbrandt, 2010: Probabilistic Parameterizations
1538	! of Visibility Using Observations of Rain Precipitation Rate, Relative Humidity, and Visibility.
1539	! J. Appl. Meteor. Climatol., 49, 36-46, https://doi.org/10.1175/2009JAMC1927.1
1540	! Interest is focused on a 'general' fog/visibilty approach, thus the fit at 50 % of probability
1541	! is chosen
1542	! Effects from precipitation are not considered
1543
1544	IMPLICIT NONE
1545
1546	REAL(r_k), INTENT(in) :: qv, ta, pres
1547	INTEGER, INTENT(out) :: fog
1548	REAL(r_k), INTENT(out) :: vis
1549
1550	! Local
1551	REAL(r_k) :: rh
1552
1553	!!!!!!! Variables
1554	! qv: mixing ratio in [kgkg-1]
1555	! ta: temperature [K]
1556	! pres: pressure field [Pa]
1557	! rh: relative humidity [1]
1558	! fog: presence of fog (1: yes, 0: no)
1559	! vis: visibility within fog [km]
1560
1561	fname = 'var_fog_FRAML50'
1562
1563	CALL var_hur(ta, pres, qv, rh)
1564	! Avoiding supersaturation
1565	rh = MINVAL((/1.,rh/))
1566
1567	IF (rh > 0.3) THEN
1568	vis = -5.1910.(-10)(rh100.)*5.44+40.10
1569	! Fog definition (vis <= 1. km)
1570	IF (vis <= oneRK) THEN
1571	fog = 1
1572	ELSE
1573	vis = -oneRK
1574	fog = 0
1575	END IF
1576	ELSE
1577	vis = -oneRK
1578	fog = 0
1579	END IF
1580
1581	END SUBROUTINE var_fog_FRAML50
1582
1583	SUBROUTINE var_range_faces(d, lon, lat, hgt, dist, filt, newrange, hvalleyrange, hgtmax, ihgtmax, &
1584	dhgt, Npeaks, ipeaks, Nvalleys, ivalleys, faces0, faces, Nranges, ranges, rangeshgtmax, &
1585	irangeshgtmax)
1586	! Subroutine to compute faces [uphill, valleys, downhill] of a monuntain range along a face
1587
1588	IMPLICIT NONE
1589
1590	INTEGER, INTENT(in) :: d
1591	REAL(r_k), INTENT(in) :: filt, newrange, hvalleyrange
1592	REAL(r_k), DIMENSION(d), INTENT(in) :: lon, lat, hgt, dist
1593	INTEGER, INTENT(out) :: ihgtmax, Npeaks, Nvalleys, Nranges
1594	INTEGER, DIMENSION(d), INTENT(out) :: ipeaks, ivalleys, faces0, faces, &
1595	irangeshgtmax
1596	INTEGER, DIMENSION(2,d), INTENT(out) :: ranges
1597	REAL(r_k), INTENT(out) :: hgtmax
1598	REAL(r_k), DIMENSION(d), INTENT(out) :: dhgt, rangeshgtmax
1599
1600	! Local
1601	INTEGER :: i, j, j1, k, l, m
1602	INTEGER :: iface
1603	INTEGER :: Nfaces, Nfaces1, Npeaks1, Nvalleys1
1604	INTEGER :: fbeg, fend
1605	REAL(r_k) :: dd
1606	INTEGER, DIMENSION(1) :: ihmax
1607	INTEGER, DIMENSION(d) :: ddhgt, Ndhgt
1608	INTEGER, DIMENSION(d) :: faces1, Ndhgt1, ipeaks1, ivalleys1
1609	REAL(r_k), DIMENSION(d) :: dLl, peaks, valleys, peaks1, valleys1
1610	REAL(r_k), DIMENSION(d) :: sortedpeaks
1611	INTEGER, DIMENSION(d) :: isortedpeaks
1612	LOGICAL :: firstvalley, rangewithin, peakwithin, &
1613	valleywithin
1614
1615	!!!!!!! Variables
1616	! lon: longitude [degrees east]
1617	! lat: latitude [degrees north]
1618	! hgt: topograpical height [m]
1619	! dist: distance between grid points [m]
1620	! filt: distance to use to filter the topographic values [m]. used to define:
1621	! the minimum length of a given section
1622	! newrange: distance to start a new mountain range [m]
1623	! hvalleyrange: maximum height of a valley to mark change of range [m]
1624
1625	fname = 'var_range_faces'
1626
1627	!PRINT *, 'heights:', hgt
1628
1629	! Looking for the maximum height
1630	hgtmax = MAXVAL(hgt)
1631	ihmax = MAXLOC(hgt)
1632	ihgtmax = ihmax(1)
1633
1634	!PRINT *, 'height max:', hgtmax, 'location:', ihmax
1635
1636	! Filtering height [NOT needed]
1637	! CALL runmean_F1D(d, hgt, filt, 'progressfill', hgtrunmean)
1638
1639	! range slope
1640	dhgt(1:d) = hgt(2:d) - hgt(1:d-1)
1641	dhgt = dhgt/dist(1:d)
1642
1643	!PRINT *, 'slope:', dhgt
1644
1645	! Initial classification
1646	Npeaks = 0
1647	Nvalleys = 0
1648	DO i=1, d-1
1649	IF (dhgt(i) > 0.) THEN
1650	IF (hgt(i) < hvalleyrange) THEN
1651	faces0(i) = fillValI
1652	ELSE
1653	faces0(i) = 2
1654	END IF
1655	ELSE IF (dhgt(i) < 0.) THEN
1656	IF (hgt(i) < hvalleyrange) THEN
1657	faces0(i) = fillValI
1658	ELSE
1659	faces0(i) = -2
1660	END IF
1661	ELSE
1662	IF (hgt(i) == zeroRK .AND. hgt(i+1) == zeroRK) THEN
1663	faces0(i) = fillValI
1664	ELSE
1665	faces0(i) = 0
1666	END IF
1667	END IF
1668	! peaks
1669	IF (dhgt(i) > 0. .AND. dhgt(i+1) < 0.) THEN
1670	Npeaks = Npeaks + 1
1671	peaks(Npeaks) = hgt(i+1)
1672	ipeaks(Npeaks) = i+1
1673	END IF
1674	! valleys
1675	IF (dhgt(i) < 0. .AND. dhgt(i+1) > 0.) THEN
1676	Nvalleys = Nvalleys + 1
1677	valleys(Nvalleys) = hgt(i+1)
1678	ivalleys(Nvalleys) = i+1
1679	END IF
1680	END DO
1681
1682	!PRINT *, 'faces:', faces0
1683	!PRINT *, Npeaks, ' peaks:', peaks(1:Npeaks), ' ipeak:', ipeaks(1:Npeaks)
1684	!PRINT *, Nvalleys, ' valleys:', valleys(1:Nvalleys), ' ivalley:', ivalleys(1:Nvalleys)
1685
1686	! tendency of the slope
1687	ddhgt(1) = 1
1688	Nfaces = 1
1689	Ndhgt(Nfaces) = 1
1690	DO i=2, d-1
1691	IF (faces0(i) /= faces0(i-1)) THEN
1692	Nfaces = Nfaces + 1
1693	Ndhgt(Nfaces) = 1
1694	ddhgt(i) = ddhgt(i-1) + 1
1695	ELSE
1696	Ndhgt(Nfaces) = Ndhgt(Nfaces) + 1
1697	ddhgt(i) = ddhgt(i-1)
1698	END IF
1699	END DO
1700
1701	!PRINT *, Nfaces, ' length faces:', Ndhgt(1:Nfaces)
1702
1703	! Defining quantitiy of ranges within the face
1704	! sorting peaks within the face and defining ranges as maximum peaks distanced > newrage
1705
1706	!sortedpeaks = peaks
1707	!CALL SortR_K(sortedpeaks, Npeaks)
1708	!isortedpeaks = 0
1709	!DO i=1, Npeaks
1710	! DO j=1, Npeaks
1711	! IF (peaks(j) == sortedpeaks(i)) isortedpeaks(i) = ipeaks(j)
1712	! END DO
1713	!END DO
1714
1715	ranges = -1
1716	Nranges = 0
1717	l = 0
1718	DO i=1,d
1719	IF (hgt(i) >= hvalleyrange) THEN
1720	IF (Nranges == 0) THEN
1721	Nranges = 1
1722	ranges(1,Nranges) = i
1723	ELSE
1724	IF (ranges(2,Nranges) /= -1) THEN
1725	Nranges = Nranges + 1
1726	ranges(1,Nranges) = i
1727	END IF
1728	END IF
1729	ELSE
1730	IF ((ranges(1,Nranges) /= -1) .AND. (ranges(2,Nranges) == -1)) ranges(2,Nranges) = i-1
1731	END IF
1732	END DO
1733	IF (ranges(2,Nranges) == -1) ranges(2,Nranges) = d
1734
1735	! Getting characteristics of ranges
1736	irangeshgtmax = -1
1737	rangeshgtmax = -1000.
1738	k = 1
1739	DO i=1, Nranges
1740	DO j=ranges(1,i), ranges(2,i)
1741	IF (hgt(j) > rangeshgtmax(i)) THEN
1742	irangeshgtmax(i) = j
1743	rangeshgtmax(i) = hgt(j)
1744	END IF
1745	END DO
1746	END DO
1747	!PRINT *, Nranges, ' _______'
1748	!DO i=1, Nranges
1749	! PRINT *,i,':', ranges(:,i),' max:', rangeshgtmax(i), ' ,', irangeshgtmax(i)
1750	!END DO
1751
1752	! Defining valleys as that consecutive grid points below surrounding peaks from the same range and
1753	! below range max
1754	k = 1
1755	IF (Npeaks > 1) THEN
1756	j = 1
1757	j1 = 2
1758	DO i=2, d
1759	IF (i >= ipeaks(j) .AND. i < ipeaks(j1)) THEN
1760	IF (hgt(i) < peaks(j) .AND. hgt(i) < peaks(j1)) faces0(i) = 0
1761	IF (i == ipeaks(j) .AND. hgt(i) < rangeshgtmax(k)) faces0(i) = 0
1762	ELSE IF (i == ipeaks(j1)) THEN
1763	j = j1
1764	j1 = j1 + 1
1765	IF (j1 > Npeaks) EXIT
1766	END IF
1767	END DO
1768	IF (ipeaks(j1) > ranges(2,k)) k = k+1
1769	END IF
1770
1771	! Using defined ranges for the last refinement.
1772	! Uphills after range maximum peak as valley: if dhgt[i] > 0. and i > rangemax --> face[i] = 0
1773	! Downhills before range maximum peak as valley: if dhgt[i] < 0. and i < rangemax --> face[i] = 0
1774	! face values before first range, removed: if dhgt[i] != 0. and i < range[1,1] --> face[i] = fillvalueI
1775	! face values after last range, removed: if dhgt[i] != 0. and i < range[2,Nranges] --> face[i] = fillvalueI
1776	! Uphill valleys as 0.5: if face[i] == 0. and i < rangemax --> face[i] = 0.5
1777	! downhill valleys as 0.5: if face[i] == 0. and i > rangemax --> face[i] = -0.5
1778	! range peaks as 0. face[rangemax] = 0.
1779	DO i=1,Nranges
1780	DO j=ranges(1,i), ranges(2,i)
1781	IF (dhgt(j) > 0. .AND. j > irangeshgtmax(i)) faces0(j) = 0
1782	IF (dhgt(j) < 0. .AND. j < irangeshgtmax(i)) faces0(j) = 0
1783	IF (faces0(j) == 0. .AND. j < irangeshgtmax(i)) faces0(j) = 1
1784	IF (faces0(j) == 0. .AND. j > irangeshgtmax(i)) faces0(j) = -1
1785	END DO
1786	faces0(irangeshgtmax(i)) = 0
1787	END DO
1788	IF (dhgt(j) /= 0. .AND. j < ranges(1,1)) faces0(j) = fillValI
1789	IF (dhgt(j) /= 0. .AND. j > ranges(2,Nranges)) faces0(j) = fillValI
1790
1791	! Removing all that grid points outside any range
1792	DO j=1, d
1793	IF (j < ranges(1,1)) THEN
1794	faces0(j) = fillValI
1795	ELSE IF (j > ranges(2,Nranges)) THEN
1796	faces0(j) = fillValI
1797	ELSE
1798	rangewithin = .FALSE.
1799	DO i=1,Nranges
1800	IF (j >= ranges(1,i) .AND. j <= ranges(2,i)) THEN
1801	rangewithin = .TRUE.
1802	EXIT
1803	END IF
1804	END DO
1805	IF (.NOT.rangewithin) faces0(j) = fillValI
1806	END IF
1807	END DO
1808
1809	! classification using filt.
1810	! On that section of the hill smaller than the selected filter. Replace values by the ones from
1811	! the previous section (except if it is already valley!)
1812	! e.g.: filt= 3
1813	! faces0 = 0, 0, 1, 1, 1, 1, -1, -1, 1, 1, 1, 1, 1, 1, -1, -1, -1, -1 -1, -1, 0, 0
1814	! faces = 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, -1, -1, -1, -1 -1, -1, 0, 0
1815	!!faces1 = faces0
1816	!!Nfaces1 = 1
1817	!!Npeaks1 = 0
1818	!!Nvalleys1 = 0
1819	!!Ndhgt1(1) = Ndhgt(1)
1820	!!DO i=2, Nfaces
1821	!! fbeg = SUM(Ndhgt(1:i-1))
1822	!! fend = fbeg + Ndhgt(i)
1823	!! peakwithin = .FALSE.
1824	!! valleywithin = .FALSE.
1825	!! DO j=1, Npeaks
1826	!! IF (ipeaks(j) == fbeg) THEN
1827	!! k = j
1828	!! peakwithin = .TRUE.
1829	!! EXIT
1830	!! END IF
1831	!! END DO
1832	!! DO j=1, Nvalleys
1833	!! IF (ivalleys(j) == fbeg) THEN
1834	!! m = j
1835	!! valleywithin = .TRUE.
1836	!! EXIT
1837	!! END IF
1838	!! END DO
1839
1840	!! dd = SUM(dist(fbeg:fend))
1841	!! IF (Ndhgt(i) < filt) THEN
1842	!! PRINT *, 'Lluis replacing !!!', i, ':', fbeg, fend, '<>', faces0(fbeg-1)
1843	!! IF (faces0(fbeg) /= 0) faces1(fbeg:fend) = faces0(fbeg-1)
1844	!! Ndhgt1(Nfaces1) = Ndhgt1(Nfaces1) + fend - fbeg + 1
1845	!! peakwithin = .FALSE.
1846	!! valleywithin = .FALSE.
1847	!! ELSE
1848	!! Nfaces1 = Nfaces1 + 1
1849	!! Ndhgt1(Nfaces1) = Ndhgt(i)
1850	!! END IF
1851
1852	!! IF (peakwithin) THEN
1853	!! Npeaks1 = Npeaks1 + 1
1854	!! peaks1(Npeaks1) = peaks(k)
1855	!! END IF
1856	!! IF (valleywithin) THEN
1857	!! Nvalleys1 = Nvalleys1 + 1
1858	!! valleys1(Nvalleys1) = valleys(m)
1859	!! END IF
1860	!!END DO
1861
1862	! Introducing valleys between peaks lower than the highest peak in the section
1863	faces = faces1
1864	IF (Npeaks1 > 1) THEN
1865	DO i=1,Npeaks1,2
1866	fbeg = ipeaks1(i)
1867	fend = ipeaks1(i+1)
1868	faces(fbeg:fend) = 0
1869	END DO
1870	END IF
1871
1872
1873	! Removing above sea points
1874	DO i=1, d-1
1875	IF (hgt(i) == zeroRK .AND. hgt(i+1) == zeroRK) dhgt(i) = fillVal64
1876	END DO
1877
1878	!PRINT *, '# lon[0] lat[1] heights[2] slope[3] faces0[4] faces1[5] faces1[6] ddhgt[7]'
1879	!DO i=1, d
1880	! PRINT *, lon(i), lat(i), hgt(i), dhgt(i), faces0(i), faces1(i), faces(i), ddhgt(i)
1881	!END DO
1882
1883	RETURN
1884
1885	END SUBROUTINE var_range_faces
1886
1887	SUBROUTINE var_hur(t, press, qv, hur)
1888	! Subroutine to compute relative humidity using August-Roche-Magnus approximation [1]
1889
1890	IMPLICIT NONE
1891
1892	REAL, INTENT(in) :: t, press, qv
1893	REAL, INTENT(out) :: hur
1894
1895	! Local
1896	REAL :: tC, es, ws
1897
1898	!!!!!!! Variables
1899	! t: temperature [K]
1900	! press: pressure [Pa]
1901	! q: mixing ratio [kgkg-1]
1902	! dz: vertical extension
1903	! hur: relative humidity [1]
1904
1905	fname = 'var_hur'
1906
1907	! August - Roche - Magnus formula (Avoiding overflow on last level)
1908	tC = t - SVPT0
1909
1910	es = ARM1 * exp(ARM2*tC/(tC+ARM3))
1911	! Saturated mixing ratio
1912	ws = mol_watdryes/(0.01press-es)
1913
1914	! Relative humidity
1915	hur = qv / ws
1916
1917	RETURN
1918
1919	END SUBROUTINE var_hur
1920
1921	REAL(r_k) FUNCTION var_tws_S11(ta0, hur0)
1922	! Function to compute wet bulb temperature using equation after:
1923	! Stull, R. (2011), J. Appl. Meteor. Climatol. 50(11):2267-2269. doi: 10.1175/JAMC-D-11-0143.1
1924
1925	IMPLICIT NONE
1926
1927	REAL(r_k), INTENT(in) :: ta0, hur0
1928
1929	! Local
1930	REAL :: ta, hur
1931
1932	!!!!!!! Variables
1933	! ta0: temperature [K] (does it only make sense if it is at 2 m?)
1934	! hur0: relative humidity [1] (does it only make sense if it is at 2 m?)
1935	! tws: wet bulb temperature [C]
1936
1937	fname = 'var_tws_S11'
1938
1939	ta = ta0 - SVPT0
1940	hur = hur0100.oneRK
1941
1942	var_tws_S11 = ta * ATAN(0.151977*SQRT(hur+8.313659)) + ATAN(ta+hur) - ATAN(hur-1.676331) + &
1943	0.00391838(hur)(1.5)ATAN(0.023101*hur) - 4.686035
1944
1945	RETURN
1946
1947	END FUNCTION var_tws_S11
1948
1949	SUBROUTINE Svar_tws_S11(ta0, hur0, tws)
1950	! Subroutine to compute wet bulb temperature using equation after:
1951	! Stull, R. (2011), J. Appl. Meteor. Climatol. 50(11):2267-2269. doi: 10.1175/JAMC-D-11-0143.1
1952
1953	IMPLICIT NONE
1954
1955	REAL(r_k), INTENT(in) :: ta0, hur0
1956	REAL(r_k), INTENT(out) :: tws
1957
1958	! Local
1959	REAL :: ta, hur
1960
1961	!!!!!!! Variables
1962	! ta0: temperature [K] (does it only make sense if it is at 2 m?)
1963	! hur0: relative humidity [1] (does it only make sense if it is at 2 m?)
1964	! tws: wet bulb temperature [C]
1965
1966	fname = 'var_tws_S11'
1967
1968	ta = ta0 - SVPT0
1969	hur = hur0100.oneRK
1970
1971	tws = ta * ATAN(0.151977*SQRT(hur+8.313659)) + ATAN(ta+hur) - ATAN(hur-1.676331) + &
1972	0.00391838(hur)(1.5)ATAN(0.023101*hur) - 4.686035
1973
1974	RETURN
1975
1976	END SUBROUTINE Svar_tws_S11
1977
1978	END MODULE module_ForDiagnosticsVars

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

Download in other formats:

Original Format