Context Navigation

module_ForDiagnosticsVars.f90 @ 2210

Last change on this file since 2210 was 2209, checked in by lfita, 7 years ago
Adding `range_faces'
File size: 57.0 KB

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

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

Download in other formats:

Original Format