Context Navigation

module_ForDiagnosticsVars.f90 @ 1908

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

Adding:

`fog_K84': Computation of fog and visibility following Kunkel, (1984)
`fog_RUC': Computation of fog and visibility following RUC method Smirnova, (2000)

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

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

Download in other formats:

Original Format