Context Navigation

module_ForDiagnosticsVars.f90 @ 2644

Last change on this file since 2644 was 2642, checked in by lfita, 6 years ago

Adding:

`compute_front_R04d3': Subroutine to compute presence of a front following Rodrigues et al.(2004)
`var_front_R04': Subroutine to compute presence of a front following Rodrigues et al.(2004)

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

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

Download in other formats:

Original Format