Context Navigation

module_ForDiagnosticsVars.f90 @ 2470

Last change on this file since 2470 was 2387, checked in by lfita, 6 years ago

Adding:

`tws': Wet Bulb temperature after Stull, 2011

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

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

Download in other formats:

Original Format