Context Navigation

module_ForDiagnosticsVars.f90 @ 1879

Last change on this file since 1879 was 1804, checked in by lfita, 8 years ago

Adding:

`evspsblpot': potential evapotranspiration following Penman-Monteith formulation implemented in ORCHIDEE

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

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

Download in other formats:

Original Format