Context Navigation

module_ForDiagnosticsVars.f90 @ 2196

Last change on this file since 2196 was 1909, checked in by lfita, 7 years ago

Adding:

`fog_FRAML50': fog and visibility following Gultepe and Milbrandt, (2010)
`var_hus': relative humidity using August-Roche-Magnus approximation [1]

Fixing: fog_K84' and fog_RUC'

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

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

Download in other formats:

Original Format