Context Navigation

module_ForDiagnosticsVars.f90 @ 2828

Last change on this file since 2828 was 2675, checked in by lfita, 5 years ago

Adding:

'frontogenesis': Calculation of the frontogenesis diagnostics

File size: 72.2 KB

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

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

Download in other formats:

Original Format