Context Navigation

module_ForDiagnosticsVars.f90 @ 2341

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

Adding to 'range_faces':

-2:|downhill;-1:|downhill-valley;|0:|peak|;1:|uphill-valley;2:|uphill
filtering outside any range

File size: 64.9 KB

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

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

Download in other formats:

Original Format