Context Navigation

module_ForDiagnosticsVars.f90 @ 2196

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

Adding:

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

Fixing: fog_K84' and fog_RUC'

File size: 54.7 KB

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

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

Download in other formats:

Original Format