Context Navigation

module_ForDiagnosticsVars.f90 @ 2370

Last change on this file since 2370 was 2223, checked in by lfita, 7 years ago

Adding to 'range_faces':

-2:|downhill;-1:|downhill-valley;|0:|peak|;1:|uphill-valley;2:|uphill
filtering outside any range

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

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

Download in other formats:

Original Format