Context Navigation

module_ForDiagnosticsVars.f90 @ 1802

Last change on this file since 1802 was 1795, checked in by lfita, 7 years ago

Adding:

`psl_ecmwf': sea-level pressure computation following ECMWF

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

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

Download in other formats:

Original Format