Context Navigation

module_ForDiagnosticsVars.f90 @ 2790

Last change on this file since 2790 was 2675, checked in by lfita, 6 years ago

Adding:

'frontogenesis': Calculation of the frontogenesis diagnostics

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

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

Download in other formats:

Original Format