Context Navigation

module_ForDiagnosticsVars.f90 @ 1776

Last change on this file since 1776 was 1776, checked in by lfita, 7 years ago
Adding `zwind': wind at a given height following power-law methodology
File size: 39.0 KB

Rev	Line
[772]	1	!! Fortran version of different diagnostics
	2	! L. Fita. LMD May 2016
	3	! gfortran module_generic.o -c module_ForDiagnostics.F90
	4	!
	5	! f2py -m module_ForDiagnostics --f90exec=/usr/bin/gfortran-4.7 -c module_generic.F90 module_ForDiagnostics.F90
	6
	7	MODULE module_ForDiagnosticsVars
	8
[1608]	9	USE module_definitions
[772]	10	USE module_generic
	11
	12	IMPLICIT NONE
	13
	14	CONTAINS
	15
	16	!!!!!!! Variables
[1769]	17	! compute_psl_ptarget4d2: Compute sea level pressure using a target pressure. Similar to the Benjamin
	18	! and Miller (1990). Method found in p_interp.F90
	19	! compute_tv4d: 4D calculation of virtual temperaure
	20	! SaturationMixingRatio: WRF's AFWA method to compute the saturation mixing ratio
	21	! The2T: WRF's AFWA method to compute the temperature at any pressure level along a saturation adiabat
	22	! by iteratively solving for it from the parcel thetae.
	23	! Theta: WRF's AFWA method to compute potential temperature
	24	! Thetae: WRF's AFWA method to compute equivalent potential temperature
	25	! TLCL: WRF's AFWA method to compute the temperature of a parcel of air would have if lifed dry
	26	! adiabatically to it's lifting condensation level (lcl)
	27	! var_cape_afwa1D: WRF's AFWA method to compute cape, cin, fclp, fclz and li
[772]	28	! var_cllmh: low, medium, high-cloud [0,1]
	29	! var_clt: total cloudiness [0,1]
[1773]	30	! var_zmla_generic: Subroutine to compute pbl-height following a generic method
[1776]	31	! var_zwind: Subroutine to extrapolate the wind at a given height following the 'power law' methodology
[1769]	32	! VirtualTemp1D: Function for 1D calculation of virtual temperaure
	33	! VirtualTemperature: WRF's AFWA method to compute virtual temperature
[772]	34
	35	!!!!!!! Calculations
	36	! compute_clt: Computation of total cloudiness
	37
	38	!!!
	39	! Variables
	40	!!!
	41
	42	FUNCTION var_cllmh(clfra, p, dz)
	43	! Function to compute cllmh on a 1D column 1: low-cloud; 2: medium-cloud; 3: high-cloud [1]
	44
	45	IMPLICIT NONE
	46
	47	INTEGER, INTENT(in) :: dz
	48	REAL(r_k), DIMENSION(dz), INTENT(in) :: clfra, p
	49	REAL(r_k), DIMENSION(3) :: var_cllmh
	50
	51	! Local
	52	INTEGER :: iz
	53	REAL(r_k) :: zclearl, zcloudl, zclearm, zcloudm, &
	54	zclearh, zcloudh
	55
	56	!!!!!!! Variables
	57	! clfra: cloudfraction as 1D verical-column [1]
	58	! p: pressure values of the column
	59	fname = 'var_cllmh'
	60
[1608]	61	zclearl = oneRK
	62	zcloudl = zeroRK
	63	zclearm = oneRK
	64	zcloudm = zeroRK
	65	zclearh = oneRK
	66	zcloudh = zeroRK
[772]	67
[1608]	68	var_cllmh = oneRK
[772]	69
	70	DO iz=1, dz
	71	IF (p(iz) < prmhc) THEN
[1608]	72	var_cllmh(3) = var_cllmh(3)*(oneRK-MAX(clfra(iz),zcloudh))/(oneRK-MIN(zcloudh,oneRK-ZEPSEC))
[772]	73	zcloudh = clfra(iz)
	74	ELSE IF ( (p(iz) >= prmhc) .AND. (p(iz) < prmlc) ) THEN
[1608]	75	var_cllmh(2) = var_cllmh(2)*(oneRK-MAX(clfra(iz),zcloudm))/(oneRK-MIN(zcloudm,oneRK-ZEPSEC))
[772]	76	zcloudm = clfra(iz)
	77	ELSE IF (p(iz) >= prmlc) THEN
[1608]	78	var_cllmh(1) = var_cllmh(1)*(oneRK-MAX(clfra(iz),zcloudl))/(oneRK-MIN(zcloudl,oneRK-ZEPSEC))
[772]	79	zcloudl = clfra(iz)
	80	ELSE
	81	PRINT *,' ' // TRIM(fname) // ': This is weird, pressure:', p(iz), ' Pa fails out!!'
	82	PRINT *,' from high, low cloud pressures:', prmhc, ' ,', prmlc,' Pa at z-level:', iz
	83	PRINT *,' p_high > p:', prmhc,'> ',p(iz),' Pa'
	84	PRINT *,' p_low > p >= p_high:', prmlc,'> ',p(iz),' >=', prmhc,' Pa'
	85	PRINT *,' p_low >= p:', prmlc,'>= ',p(iz),' Pa'
	86	STOP
	87	END IF
	88	END DO
	89
[1608]	90	var_cllmh = oneRK - var_cllmh
[772]	91
	92	RETURN
	93
	94	END FUNCTION var_cllmh
	95
	96	REAL(r_k) FUNCTION var_clt(clfra, dz)
	97	! Function to compute the total cloud following 'newmicro.F90' from LMDZ using 1D vertical
	98	! column values
	99
	100	IMPLICIT NONE
	101
[1141]	102	INTEGER, INTENT(in) :: dz
[772]	103	REAL(r_k), DIMENSION(dz), INTENT(in) :: clfra
	104	! Local
	105	INTEGER :: iz
	106	REAL(r_k) :: zclear, zcloud
[1608]	107
[772]	108	!!!!!!! Variables
	109	! cfra: 1-column cloud fraction values
	110
	111	fname = 'var_clt'
	112
[1608]	113	zclear = oneRK
	114	zcloud = zeroRK
[772]	115
	116	DO iz=1,dz
[1608]	117	zclear = zclear*(oneRK-MAX(clfra(iz),zcloud))/(oneRK-MIN(zcloud,1.-ZEPSEC))
	118	var_clt = oneRK - zclear
[772]	119	zcloud = clfra(iz)
	120	END DO
	121
	122	RETURN
	123
	124	END FUNCTION var_clt
	125
[1769]	126
	127	SUBROUTINE compute_psl_ptarget4d2(press, ps, hgt, ta, qv, ptarget, psl, d1, d2, d3, d4)
	128	! Subroutine to compute sea level pressure using a target pressure. Similar to the Benjamin
	129	! and Miller (1990). Method found in p_interp.F90
	130
	131	IMPLICIT NONE
	132
	133	INTEGER, INTENT(in) :: d1, d2, d3, d4
	134	REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in) :: press, ta, qv
	135	REAL(r_k), DIMENSION(d1,d2,d4), INTENT(in) :: ps
	136	REAL(r_k), DIMENSION(d1,d2), INTENT(in) :: hgt
	137	REAL(r_k), INTENT(in) :: ptarget
	138	REAL(r_k), DIMENSION(d1,d2,d4), INTENT(out) :: psl
	139
	140	! Local
	141	INTEGER :: i, j, it
	142	INTEGER :: kin
	143	INTEGER :: kupper
	144	REAL(r_k) :: dpmin, dp, tbotextrap, &
	145	tvbotextrap, virtual
	146	! Exponential related to standard atmosphere lapse rate r_d*gammav/g
	147	REAL(r_k), PARAMETER :: expon=r_d*gammav/grav
	148
	149	!!!!!!! Variables
	150	! press: Atmospheric pressure [Pa]
	151	! ps: surface pressure [Pa]
	152	! hgt: surface height
	153	! ta: temperature [K]
	154	! qv: water vapor mixing ratio
	155	! dz: number of vertical levels
	156	! psl: sea-level pressure
	157
	158	fname = 'compute_psl_ptarget4d2'
	159
	160	! Minimal distance between pressures [Pa]
	161	dpmin=1.e4
	162	psl=0.
	163
	164	DO i=1,d1
	165	DO j=1,d2
	166	IF (hgt(i,j) /= 0.) THEN
	167	DO it=1,d4
	168
	169	! target pressure to be used for the extrapolation [Pa] (defined in namelist.input)
	170	! ptarget = 70000. default value
	171
	172	! We are below both the ground and the lowest data level.
	173
	174	! First, find the model level that is closest to a "target" pressure
	175	! level, where the "target" pressure is delta-p less that the local
	176	! value of a horizontally smoothed surface pressure field. We use
	177	! delta-p = 150 hPa here. A standard lapse rate temperature profile
	178	! passing through the temperature at this model level will be used
	179	! to define the temperature profile below ground. This is similar
	180	! to the Benjamin and Miller (1990) method, using
	181	! 700 hPa everywhere for the "target" pressure.
	182
	183	kupper = 0
	184	loop_kIN: DO kin=d3,1,-1
	185	kupper = kin
	186	dp=abs( press(i,j,kin,it) - ptarget )
	187	IF (dp .GT. dpmin) EXIT loop_kIN
	188	dpmin=min(dpmin,dp)
	189	ENDDO loop_kIN
	190
	191	tbotextrap=ta(i,j,kupper,it)(ps(i,j,it)/ptarget)*expon
	192	tvbotextrap=virtualTemp1D(tbotextrap,qv(i,j,kupper,it))
	193
	194	psl(i,j,it) = ps(i,j,it)((tvbotextrap+gammavhgt(i,j))/tvbotextrap)**(1/expon)
	195	END DO
	196	ELSE
	197	psl(i,j,:) = ps(i,j,:)
	198	END IF
	199	END DO
	200	END DO
	201
	202	RETURN
	203
	204	END SUBROUTINE compute_psl_ptarget4d2
	205
	206	SUBROUTINE compute_tv4d(ta,qv,tv,d1,d2,d3,d4)
	207	! 4D calculation of virtual temperaure
	208
	209	IMPLICIT NONE
	210
	211	INTEGER, INTENT(in) :: d1, d2, d3, d4
	212	REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in) :: ta, qv
	213	REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(out) :: tv
	214
	215	! Variables
	216	! ta: temperature [K]
	217	! qv: mixing ratio [kgkg-1]
	218	! tv: virtual temperature
	219
	220	tv = ta*(oneRK+(qv/epsilonv))/(oneRK+qv)
	221
	222	END SUBROUTINE compute_tv4d
	223
	224	FUNCTION VirtualTemp1D (ta,qv) result (tv)
	225	! 1D calculation of virtual temperaure
	226
	227	IMPLICIT NONE
	228
	229	REAL(r_k), INTENT(in) :: ta, qv
	230	REAL(r_k) :: tv
	231
	232	! Variables
	233	! ta: temperature [K]
	234	! qv: mixing ratio [kgkg-1]
	235
	236	tv = ta*(oneRK+(qv/epsilonv))/(oneRK+qv)
	237
	238	END FUNCTION VirtualTemp1D
	239
	240	! ---- BEGIN modified from module_diag_afwa.F ---- !
	241
	242	!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
	243	!~
	244	!~ Name:
	245	!~ Theta
	246	!~
	247	!~ Description:
	248	!~ This function calculates potential temperature as defined by
	249	!~ Poisson's equation, given temperature and pressure ( hPa ).
	250	!~
	251	!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
	252	FUNCTION Theta ( t, p )
	253
	254	IMPLICIT NONE
	255
	256	!~ Variable declaration
	257	! --------------------
	258	REAL(r_k), INTENT ( IN ) :: t
	259	REAL(r_k), INTENT ( IN ) :: p
	260	REAL(r_k) :: theta
	261
	262	! Using WRF values
	263	!REAL :: Rd ! Dry gas constant
	264	!REAL :: Cp ! Specific heat of dry air at constant pressure
	265	!REAL :: p00 ! Standard pressure ( 1000 hPa )
	266	REAL(r_k) :: Rd, p00
	267
	268	!Rd = 287.04
	269	!Cp = 1004.67
	270	!p00 = 1000.00
	271
	272	Rd = r_d
	273	p00 = p1000mb/100.
	274
	275	!~ Poisson's equation
	276	! ------------------
	277	theta = t * ( (p00/p)**(Rd/Cp) )
	278
	279	END FUNCTION Theta
	280
	281	!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
	282	!~
	283	!~ Name:
	284	!~ Thetae
	285	!~
	286	!~ Description:
	287	!~ This function returns equivalent potential temperature using the
	288	!~ method described in Bolton 1980, Monthly Weather Review, equation 43.
	289	!~
	290	!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
	291	FUNCTION Thetae ( tK, p, rh, mixr )
	292
	293	IMPLICIT NONE
	294
	295	!~ Variable Declarations
	296	! ---------------------
	297	REAL(r_k) :: tK ! Temperature ( K )
	298	REAL(r_k) :: p ! Pressure ( hPa )
	299	REAL(r_k) :: rh ! Relative humidity
	300	REAL(r_k) :: mixr ! Mixing Ratio ( kg kg^-1)
	301	REAL(r_k) :: te ! Equivalent temperature ( K )
	302	REAL(r_k) :: thetae ! Equivalent potential temperature
	303
	304	! Using WRF values
	305	!REAL, PARAMETER :: R = 287.04 ! Universal gas constant (J/deg kg)
	306	!REAL, PARAMETER :: P0 = 1000.0 ! Standard pressure at surface (hPa)
	307	REAL(r_k) :: R, p00, Lv
	308	!REAL, PARAMETER :: lv = 2.54(10*6) ! Latent heat of vaporization
	309	! (J kg^-1)
	310	!REAL, PARAMETER :: cp = 1004.67 ! Specific heat of dry air constant
	311	! at pressure (J/deg kg)
	312	REAL(r_k) :: tlc ! LCL temperature
	313
	314	R = r_d
	315	p00 = p1000mb/100.
	316	lv = XLV
	317
	318	!~ Calculate the temperature of the LCL
	319	! ------------------------------------
	320	tlc = TLCL ( tK, rh )
	321
	322	!~ Calculate theta-e
	323	! -----------------
	324	thetae = (tK * (p00/p)*( (R/Cp)(1.- ( (.28E-3)mixr1000.) ) ) )* &
	325	exp( (((3.376/tlc)-.00254))*&
	326	(mixr1000.(1.+(.81E-3)mixr1000.)) )
	327
	328	END FUNCTION Thetae
	329
	330	!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
	331	!~
	332	!~ Name:
	333	!~ The2T.f90
	334	!~
	335	!~ Description:
	336	!~ This function returns the temperature at any pressure level along a
	337	!~ saturation adiabat by iteratively solving for it from the parcel
	338	!~ thetae.
	339	!~
	340	!~ Dependencies:
	341	!~ function thetae.f90
	342	!~
	343	!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
	344	FUNCTION The2T ( thetaeK, pres, flag ) result ( tparcel )
	345
	346	IMPLICIT NONE
	347
	348	!~ Variable Declaration
	349	! --------------------
	350	REAL(r_k), INTENT ( IN ) :: thetaeK
	351	REAL(r_k), INTENT ( IN ) :: pres
	352	LOGICAL, INTENT ( INOUT ) :: flag
	353	REAL(r_k) :: tparcel
	354
	355	REAL(r_k) :: thetaK
	356	REAL(r_k) :: tovtheta
	357	REAL(r_k) :: tcheck
	358	REAL(r_k) :: svpr, svpr2
	359	REAL(r_k) :: smixr, smixr2
	360	REAL(r_k) :: thetae_check, thetae_check2
	361	REAL(r_k) :: tguess_2, correction
	362
	363	LOGICAL :: found
	364	INTEGER :: iter
	365
	366	! Using WRF values
	367	!REAL :: R ! Dry gas constant
	368	!REAL :: Cp ! Specific heat for dry air
	369	!REAL :: kappa ! Rd / Cp
	370	!REAL :: Lv ! Latent heat of vaporization at 0 deg. C
	371	REAL(r_k) :: R, kappa, Lv
	372
	373	R = r_d
	374	Lv = XLV
	375	!R = 287.04
	376	!Cp = 1004.67
	377	Kappa = R/Cp
	378	!Lv = 2.500E+6
	379
	380	!~ Make initial guess for temperature of the parcel
	381	! ------------------------------------------------
	382	tovtheta = (pres/100000.0)**(r/cp)
	383	tparcel = thetaeK/exp(lv.012/(cp295.))*tovtheta
	384
	385	iter = 1
	386	found = .false.
	387	flag = .false.
	388
	389	DO
	390	IF ( iter > 105 ) EXIT
	391
	392	tguess_2 = tparcel + REAL ( 1 )
	393
	394	svpr = 6.122 * exp ( (17.67*(tparcel-273.15)) / (tparcel-29.66) )
	395	smixr = ( 0.622*svpr ) / ( (pres/100.0)-svpr )
	396	svpr2 = 6.122 * exp ( (17.67*(tguess_2-273.15)) / (tguess_2-29.66) )
	397	smixr2 = ( 0.622*svpr2 ) / ( (pres/100.0)-svpr2 )
	398
	399	! ------------------------------------------------------------------ ~!
	400	!~ When this function was orinially written, the final parcel ~!
	401	!~ temperature check was based off of the parcel temperature and ~!
	402	!~ not the theta-e it produced. As there are multiple temperature- ~!
	403	!~ mixing ratio combinations that can produce a single theta-e value, ~!
	404	!~ we change the check to be based off of the resultant theta-e ~!
	405	!~ value. This seems to be the most accurate way of backing out ~!
	406	!~ temperature from theta-e. ~!
	407	!~ ~!
	408	!~ Rentschler, April 2010 ~!
	409	! ------------------------------------------------------------------ !
	410
	411	!~ Old way...
	412	!thetaK = thetaeK / EXP (lv * smixr /(cp*tparcel) )
	413	!tcheck = thetaK * tovtheta
	414
	415	!~ New way
	416	thetae_check = Thetae ( tparcel, pres/100., 100., smixr )
	417	thetae_check2 = Thetae ( tguess_2, pres/100., 100., smixr2 )
	418
	419	!~ Whew doggies - that there is some accuracy...
	420	!IF ( ABS (tparcel-tcheck) < .05) THEN
	421	IF ( ABS (thetaeK-thetae_check) < .001) THEN
	422	found = .true.
	423	flag = .true.
	424	EXIT
	425	END IF
	426
	427	!~ Old
	428	!tparcel = tparcel + (tcheck - tparcel)*.3
	429
	430	!~ New
	431	correction = ( thetaeK-thetae_check ) / ( thetae_check2-thetae_check )
	432	tparcel = tparcel + correction
	433
	434	iter = iter + 1
	435	END DO
	436
	437	!IF ( .not. found ) THEN
	438	! print*, "Warning! Thetae to temperature calculation did not converge!"
	439	! print*, "Thetae ", thetaeK, "Pressure ", pres
	440	!END IF
	441
	442	END FUNCTION The2T
	443
	444	!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
	445	!~
	446	!~ Name:
	447	!~ VirtualTemperature
	448	!~
	449	!~ Description:
	450	!~ This function returns virtual temperature given temperature ( K )
	451	!~ and mixing ratio.
	452	!~
	453	!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
	454	FUNCTION VirtualTemperature ( tK, w ) result ( Tv )
	455
	456	IMPLICIT NONE
	457
	458	!~ Variable declaration
	459	real(r_k), intent ( in ) :: tK !~ Temperature
	460	real(r_k), intent ( in ) :: w !~ Mixing ratio ( kg kg^-1 )
	461	real(r_k) :: Tv !~ Virtual temperature
	462
	463	Tv = tK * ( 1.0 + (w/0.622) ) / ( 1.0 + w )
	464
	465	END FUNCTION VirtualTemperature
	466
	467	!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
	468	!~
	469	!~ Name:
	470	!~ SaturationMixingRatio
	471	!~
	472	!~ Description:
	473	!~ This function calculates saturation mixing ratio given the
	474	!~ temperature ( K ) and the ambient pressure ( Pa ). Uses
	475	!~ approximation of saturation vapor pressure.
	476	!~
	477	!~ References:
	478	!~ Bolton (1980), Monthly Weather Review, pg. 1047, Eq. 10
	479	!~
	480	!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
	481	FUNCTION SaturationMixingRatio ( tK, p ) result ( ws )
	482
	483	IMPLICIT NONE
	484
	485	REAL(r_k), INTENT ( IN ) :: tK
	486	REAL(r_k), INTENT ( IN ) :: p
	487	REAL(r_k) :: ws
	488
	489	REAL(r_k) :: es
	490
	491	es = 6.122 * exp ( (17.67*(tK-273.15))/ (tK-29.66) )
	492	ws = ( 0.622*es ) / ( (p/100.0)-es )
	493
	494	END FUNCTION SaturationMixingRatio
	495
	496	!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
	497	!~
	498	!~ Name:
	499	!~ tlcl
	500	!~
	501	!~ Description:
	502	!~ This function calculates the temperature of a parcel of air would have
	503	!~ if lifed dry adiabatically to it's lifting condensation level (lcl).
	504	!~
	505	!~ References:
	506	!~ Bolton (1980), Monthly Weather Review, pg. 1048, Eq. 22
	507	!~
	508	!!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!!!
	509	FUNCTION TLCL ( tk, rh )
	510
	511	IMPLICIT NONE
	512
	513	REAL(r_k), INTENT ( IN ) :: tK !~ Temperature ( K )
	514	REAL(r_k), INTENT ( IN ) :: rh !~ Relative Humidity ( % )
	515	REAL(r_k) :: tlcl
	516
	517	REAL(r_k) :: denom, term1, term2
	518
	519	term1 = 1.0 / ( tK - 55.0 )
	520	!! Lluis
	521	! IF ( rh > REAL (0) ) THEN
	522	IF ( rh > zeroRK ) THEN
	523	term2 = ( LOG (rh/100.0) / 2840.0 )
	524	ELSE
	525	term2 = ( LOG (0.001/oneRK) / 2840.0 )
	526	END IF
	527	denom = term1 - term2
	528	!! Lluis
	529	! tlcl = ( 1.0 / denom ) + REAL ( 55 )
	530	tlcl = ( oneRK / denom ) + 55*oneRK
	531
	532	END FUNCTION TLCL
	533
	534	FUNCTION var_cape_afwa1D(nz, tk, rhv, p, hgt, sfc, cape, cin, zlfc, plfc, lidx, parcel) RESULT (ostat)
	535	! Function to compute cape on a 1D column following implementation in phys/module_diag_afwa.F
	536
	537	IMPLICIT NONE
	538
	539	INTEGER, INTENT(in) :: nz, sfc
	540	REAL(r_k), DIMENSION(nz), INTENT(in) :: tk, rhv, p, hgt
	541	REAL(r_k), INTENT(out) :: cape, cin, zlfc, plfc, lidx
	542	INTEGER :: ostat
	543	INTEGER, INTENT(in) :: parcel
	544
	545	! Local
	546	!~ Derived profile variables
	547	! -------------------------
	548	REAL(r_k), DIMENSION(nz) :: rh, ws, w, dTvK, buoy
	549	REAL(r_k) :: tlclK, plcl, nbuoy, pbuoy
	550
	551	!~ Source parcel information
	552	! -------------------------
	553	REAL(r_k) :: srctK, srcrh, srcws, srcw, srcp, &
	554	srctheta, srcthetaeK
	555	INTEGER :: srclev
	556	REAL(r_k) :: spdiff
	557
	558	!~ Parcel variables
	559	! ----------------
	560	REAL(r_k) :: ptK, ptvK, tvK, pw
	561
	562	!~ Other utility variables
	563	! -----------------------
	564	INTEGER :: i, j, k
	565	INTEGER :: lfclev
	566	INTEGER :: prcl
	567	INTEGER :: mlev
	568	INTEGER :: lyrcnt
	569	LOGICAL :: flag
	570	LOGICAL :: wflag
	571	REAL(r_k) :: freeze
	572	REAL(r_k) :: pdiff
	573	REAL(r_k) :: pm, pu, pd
	574	REAL(r_k) :: lidxu
	575	REAL(r_k) :: lidxd
	576
	577	REAL(r_k), PARAMETER :: Rd = r_d
	578	REAL(r_k), PARAMETER :: RUNDEF = -9.999E30
	579
	580	!!!!!!! Variables
	581	! nz: Number of vertical levels
	582	! sfc: Surface level in the profile
	583	! tk: Temperature profile [K]
	584	! rhv: Relative Humidity profile [1]
	585	! rh: Relative Humidity profile [%]
	586	! p: Pressure profile [Pa]
	587	! hgt: Geopotential height profile [gpm]
	588	! cape: CAPE [Jkg-1]
	589	! cin: CIN [Jkg-1]
	590	! zlfc: LFC Height [gpm]
	591	! plfc: LFC Pressure [Pa]
	592	! lidx: Lifted index
	593	! FROM: https://en.wikipedia.org/wiki/Lifted_index
	594	! lidx >= 6: Very Stable Conditions
	595	! 6 > lidx > 1: Stable Conditions, Thunderstorms Not Likely
	596	! 0 > lidx > -2: Slightly Unstable, Thunderstorms Possible, With Lifting Mechanism (i.e., cold front, daytime heating, ...)
	597	! -2 > lidx > -6: Unstable, Thunderstorms Likely, Some Severe With Lifting Mechanism
	598	! -6 > lidx: Very Unstable, Severe Thunderstorms Likely With Lifting Mechanism
	599	! ostat: Function return status (Nonzero is bad)
	600	! parcel:
	601	! Most Unstable = 1 (default)
	602	! Mean layer = 2
	603	! Surface based = 3
	604	!~ Derived profile variables
	605	! -------------------------
	606	! ws: Saturation mixing ratio
	607	! w: Mixing ratio
	608	! dTvK: Parcel / ambient Tv difference
	609	! buoy: Buoyancy
	610	! tlclK: LCL temperature [K]
	611	! plcl: LCL pressure [Pa]
	612	! nbuoy: Negative buoyancy
	613	! pbuoy: Positive buoyancy
	614
	615	!~ Source parcel information
	616	! -------------------------
	617	! srctK: Source parcel temperature [K]
	618	! srcrh: Source parcel rh [%]
	619	! srcws: Source parcel sat. mixing ratio
	620	! srcw: Source parcel mixing ratio
	621	! srcp: Source parcel pressure [Pa]
	622	! srctheta: Source parcel theta [K]
	623	! srcthetaeK: Source parcel theta-e [K]
	624	! srclev: Level of the source parcel
	625	! spdiff: Pressure difference
	626
	627	!~ Parcel variables
	628	! ----------------
	629	! ptK: Parcel temperature [K]
	630	! ptvK: Parcel virtual temperature [K]
	631	! tvK: Ambient virtual temperature [K]
	632	! pw: Parcel mixing ratio
	633
	634	!~ Other utility variables
	635	! -----------------------
	636	! lfclev: Level of LFC
	637	! prcl: Internal parcel type indicator
	638	! mlev: Level for ML calculation
	639	! lyrcnt: Number of layers in mean layer
	640	! flag: Dummy flag
	641	! wflag: Saturation flag
	642	! freeze: Water loading multiplier
	643	! pdiff: Pressure difference between levs
	644	! pm, pu, pd: Middle, upper, lower pressures
	645	! lidxu: Lifted index at upper level
	646	! lidxd: Lifted index at lower level
	647
	648	fname = 'var_cape_afwa'
	649
	650	!~ Initialize variables
	651	! --------------------
	652	rh = rhv*100.
	653	ostat = 0
	654	CAPE = zeroRK
	655	CIN = zeroRK
	656	ZLFC = RUNDEF
	657	PLFC = RUNDEF
	658
	659	!~ Look for submitted parcel definition
	660	!~ 1 = Most unstable
	661	!~ 2 = Mean layer
	662	!~ 3 = Surface based
	663	! -------------------------------------
	664	IF ( parcel > 3 .or. parcel < 1 ) THEN
	665	prcl = 1
	666	ELSE
	667	prcl = parcel
	668	END IF
	669
	670	!~ Initalize our parcel to be (sort of) surface based. Because of
	671	!~ issues we've been observing in the WRF model, specifically with
	672	!~ excessive surface moisture values at the surface, using a true
	673	!~ surface based parcel is resulting a more unstable environment
	674	!~ than is actually occuring. To address this, our surface parcel
	675	!~ is now going to be defined as the parcel between 25-50 hPa
	676	!~ above the surface. UPDATE - now that this routine is in WRF,
	677	!~ going to trust surface info. GAC 20140415
	678	! ----------------------------------------------------------------
	679
	680	!~ Compute mixing ratio values for the layer
	681	! -----------------------------------------
	682	DO k = sfc, nz
	683	ws ( k ) = SaturationMixingRatio ( tK(k), p(k) )
	684	w ( k ) = ( rh(k)/100.0 ) * ws ( k )
	685	END DO
	686
	687	srclev = sfc
	688	srctK = tK ( sfc )
	689	srcrh = rh ( sfc )
	690	srcp = p ( sfc )
	691	srcws = ws ( sfc )
	692	srcw = w ( sfc )
	693	srctheta = Theta ( tK(sfc), p(sfc)/100.0 )
	694
	695	!~ Compute the profile mixing ratio. If the parcel is the MU parcel,
	696	!~ define our parcel to be the most unstable parcel within the lowest
	697	!~ 180 mb.
	698	! -------------------------------------------------------------------
	699	mlev = sfc + 1
	700	DO k = sfc + 1, nz
	701
	702	!~ Identify the last layer within 100 hPa of the surface
	703	! -----------------------------------------------------
	704	pdiff = ( p (sfc) - p (k) ) / REAL ( 100 )
	705	IF ( pdiff <= REAL (100) ) mlev = k
	706
	707	!~ If we've made it past the lowest 180 hPa, exit the loop
	708	! -------------------------------------------------------
	709	IF ( pdiff >= REAL (180) ) EXIT
	710
	711	IF ( prcl == 1 ) THEN
	712	!IF ( (p(k) > 70000.0) .and. (w(k) > srcw) ) THEN
	713	IF ( (w(k) > srcw) ) THEN
	714	srctheta = Theta ( tK(k), p(k)/100.0 )
	715	srcw = w ( k )
	716	srclev = k
	717	srctK = tK ( k )
	718	srcrh = rh ( k )
	719	srcp = p ( k )
	720	END IF
	721	END IF
	722
	723	END DO
	724
	725	!~ If we want the mean layer parcel, compute the mean values in the
	726	!~ lowest 100 hPa.
	727	! ----------------------------------------------------------------
	728	lyrcnt = mlev - sfc + 1
	729	IF ( prcl == 2 ) THEN
	730
	731	srclev = sfc
	732	srctK = SUM ( tK (sfc:mlev) ) / REAL ( lyrcnt )
	733	srcw = SUM ( w (sfc:mlev) ) / REAL ( lyrcnt )
	734	srcrh = SUM ( rh (sfc:mlev) ) / REAL ( lyrcnt )
	735	srcp = SUM ( p (sfc:mlev) ) / REAL ( lyrcnt )
	736	srctheta = Theta ( srctK, srcp/100. )
	737
	738	END IF
	739
	740	srcthetaeK = Thetae ( srctK, srcp/100.0, srcrh, srcw )
	741
	742	!~ Calculate temperature and pressure of the LCL
	743	! ---------------------------------------------
	744	tlclK = TLCL ( tK(srclev), rh(srclev) )
	745	plcl = p(srclev) * ( (tlclK/tK(srclev))**(Cp/Rd) )
	746
	747	!~ Now lift the parcel
	748	! -------------------
	749
	750	buoy = REAL ( 0 )
	751	pw = srcw
	752	wflag = .false.
	753	DO k = srclev, nz
	754	IF ( p (k) <= plcl ) THEN
	755
	756	!~ The first level after we pass the LCL, we're still going to
	757	!~ lift the parcel dry adiabatically, as we haven't added the
	758	!~ the required code to switch between the dry adiabatic and moist
	759	!~ adiabatic cooling. Since the dry version results in a greater
	760	!~ temperature loss, doing that for the first step so we don't over
	761	!~ guesstimate the instability.
	762	! ----------------------------------------------------------------
	763
	764	IF ( wflag ) THEN
	765	flag = .false.
	766
	767	!~ Above the LCL, our parcel is now undergoing moist adiabatic
	768	!~ cooling. Because of the latent heating being undergone as
	769	!~ the parcel rises above the LFC, must iterative solve for the
	770	!~ parcel temperature using equivalant potential temperature,
	771	!~ which is conserved during both dry adiabatic and
	772	!~ pseudoadiabatic displacements.
	773	! --------------------------------------------------------------
	774	ptK = The2T ( srcthetaeK, p(k), flag )
	775
	776	!~ Calculate the parcel mixing ratio, which is now changing
	777	!~ as we condense moisture out of the parcel, and is equivalent
	778	!~ to the saturation mixing ratio, since we are, in theory, at
	779	!~ saturation.
	780	! ------------------------------------------------------------
	781	pw = SaturationMixingRatio ( ptK, p(k) )
	782
	783	!~ Now we can calculate the virtual temperature of the parcel
	784	!~ and the surrounding environment to assess the buoyancy.
	785	! ----------------------------------------------------------
	786	ptvK = VirtualTemperature ( ptK, pw )
	787	tvK = VirtualTemperature ( tK (k), w (k) )
	788
	789	!~ Modification to account for water loading
	790	! -----------------------------------------
	791	freeze = 0.033 * ( 263.15 - pTvK )
	792	IF ( freeze > 1.0 ) freeze = 1.0
	793	IF ( freeze < 0.0 ) freeze = 0.0
	794
	795	!~ Approximate how much of the water vapor has condensed out
	796	!~ of the parcel at this level
	797	! ---------------------------------------------------------
	798	freeze = freeze * 333700.0 * ( srcw - pw ) / 1005.7
	799
	800	pTvK = pTvK - pTvK * ( srcw - pw ) + freeze
	801	dTvK ( k ) = ptvK - tvK
	802	buoy ( k ) = g * ( dTvK ( k ) / tvK )
	803
	804	ELSE
	805
	806	!~ Since the theta remains constant whilst undergoing dry
	807	!~ adiabatic processes, can back out the parcel temperature
	808	!~ from potential temperature below the LCL
	809	! --------------------------------------------------------
	810	ptK = srctheta / ( 100000.0/p(k) )**(Rd/Cp)
	811
	812	!~ Grab the parcel virtual temperture, can use the source
	813	!~ mixing ratio since we are undergoing dry adiabatic cooling
	814	! ----------------------------------------------------------
	815	ptvK = VirtualTemperature ( ptK, srcw )
	816
	817	!~ Virtual temperature of the environment
	818	! --------------------------------------
	819	tvK = VirtualTemperature ( tK (k), w (k) )
	820
	821	!~ Buoyancy at this level
	822	! ----------------------
	823	dTvK ( k ) = ptvK - tvK
	824	buoy ( k ) = g * ( dtvK ( k ) / tvK )
	825
	826	wflag = .true.
	827
	828	END IF
	829
	830	ELSE
	831
	832	!~ Since the theta remains constant whilst undergoing dry
	833	!~ adiabatic processes, can back out the parcel temperature
	834	!~ from potential temperature below the LCL
	835	! --------------------------------------------------------
	836	ptK = srctheta / ( 100000.0/p(k) )**(Rd/Cp)
	837
	838	!~ Grab the parcel virtual temperture, can use the source
	839	!~ mixing ratio since we are undergoing dry adiabatic cooling
	840	! ----------------------------------------------------------
	841	ptvK = VirtualTemperature ( ptK, srcw )
	842
	843	!~ Virtual temperature of the environment
	844	! --------------------------------------
	845	tvK = VirtualTemperature ( tK (k), w (k) )
	846
	847	!~ Buoyancy at this level
	848	! ---------------------
	849	dTvK ( k ) = ptvK - tvK
	850	buoy ( k ) = g * ( dtvK ( k ) / tvK )
	851
	852	END IF
	853
	854	!~ Chirp
	855	! -----
	856	! WRITE ( ,'(I15,6F15.3)' )k,p(k)/100.,ptK,pw1000.,ptvK,tvK,buoy(k)
	857
	858	END DO
	859
	860	!~ Add up the buoyancies, find the LFC
	861	! -----------------------------------
	862	flag = .false.
	863	lfclev = -1
	864	nbuoy = REAL ( 0 )
	865	pbuoy = REAL ( 0 )
	866	DO k = sfc + 1, nz
	867	IF ( tK (k) < 253.15 ) EXIT
	868	CAPE = CAPE + MAX ( buoy (k), 0.0 ) * ( hgt (k) - hgt (k-1) )
	869	CIN = CIN + MIN ( buoy (k), 0.0 ) * ( hgt (k) - hgt (k-1) )
	870
	871	!~ If we've already passed the LFC
	872	! -------------------------------
	873	IF ( flag .and. buoy (k) > REAL (0) ) THEN
	874	pbuoy = pbuoy + buoy (k)
	875	END IF
	876
	877	!~ We are buoyant now - passed the LFC
	878	! -----------------------------------
	879	IF ( .not. flag .and. buoy (k) > REAL (0) .and. p (k) < plcl ) THEN
	880	flag = .true.
	881	pbuoy = pbuoy + buoy (k)
	882	lfclev = k
	883	END IF
	884
	885	!~ If we think we've passed the LFC, but encounter a negative layer
	886	!~ start adding it up.
	887	! ----------------------------------------------------------------
	888	IF ( flag .and. buoy (k) < REAL (0) ) THEN
	889	nbuoy = nbuoy + buoy (k)
	890
	891	!~ If the accumulated negative buoyancy is greater than the
	892	!~ positive buoyancy, then we are capped off. Got to go higher
	893	!~ to find the LFC. Reset positive and negative buoyancy summations
	894	! ----------------------------------------------------------------
	895	IF ( ABS (nbuoy) > pbuoy ) THEN
	896	flag = .false.
	897	nbuoy = REAL ( 0 )
	898	pbuoy = REAL ( 0 )
	899	lfclev = -1
	900	END IF
	901	END IF
	902
	903	END DO
	904
	905	!~ Calculate lifted index by interpolating difference between
	906	!~ parcel and ambient Tv to 500mb.
	907	! ----------------------------------------------------------
	908	DO k = sfc + 1, nz
	909
	910	pm = 50000.
	911	pu = p ( k )
	912	pd = p ( k - 1 )
	913
	914	!~ If we're already above 500mb just set lifted index to 0.
	915	!~ --------------------------------------------------------
	916	IF ( pd .le. pm ) THEN
	917	lidx = zeroRK
	918	EXIT
	919
	920	ELSEIF ( pu .le. pm .and. pd .gt. pm) THEN
	921
	922	!~ Found trapping pressure: up, middle, down.
	923	!~ We are doing first order interpolation.
	924	! ------------------------------------------
	925	lidxu = -dTvK ( k ) * ( pu / 100000. ) ** (Rd/Cp)
	926	lidxd = -dTvK ( k-1 ) * ( pd / 100000. ) ** (Rd/Cp)
	927	lidx = ( lidxu * (pm-pd) + lidxd * (pu-pm) ) / (pu-pd)
	928	EXIT
	929
	930	ENDIF
	931
	932	END DO
	933
	934	!~ Assuming the the LFC is at a pressure level for now
	935	! ---------------------------------------------------
	936	IF ( lfclev > zeroRK ) THEN
	937	PLFC = p ( lfclev )
	938	ZLFC = hgt ( lfclev )
	939	END IF
	940
	941	IF ( PLFC /= PLFC .OR. PLFC < zeroRK ) THEN
	942	PLFC = -oneRK
	943	ZLFC = -oneRK
	944	END IF
	945
	946	IF ( CAPE /= CAPE ) cape = zeroRK
	947
	948	IF ( CIN /= CIN ) cin = zeroRK
	949
	950	!~ Chirp
	951	! -----
	952	! WRITE ( , ) ' CAPE: ', cape, ' CIN: ', cin
	953	! WRITE ( , ) ' LFC: ', ZLFC, ' PLFC: ', PLFC
	954	! WRITE ( , ) ''
	955	! WRITE ( , ) ' Exiting buoyancy.'
	956	! WRITE ( , ) ' ==================================== '
	957	! WRITE ( , ) ''
	958
	959	RETURN
	960
	961	END FUNCTION var_cape_afwa1D
	962
	963	! ---- END modified from module_diag_afwa.F ---- !
	964
[1773]	965	SUBROUTINE var_zmla_generic(dz, qv, tpot, z, topo, zmla)
	966	! Subroutine to compute pbl-height following a generic method
	967	! from Nielsen-Gammon et al., 2008 J. Appl. Meteor. Clim.
	968	! applied also in Garcia-Diez et al., 2013, QJRMS
	969	! where
	970	! "The technique identifies the ML height as a threshold increase of potential temperature from
	971	! its minimum value within the boundary layer."
	972	! here applied similarly to Garcia-Diez et al. where
	973	! zmla = "...first level where potential temperature exceeds the minimum potential temperature
	974	! reached in the mixed layer by more than 1.5 K"
[1769]	975
[1773]	976	IMPLICIT NONE
	977
	978	INTEGER, INTENT(in) :: dz
	979	REAL(r_k), DIMENSION(dz), INTENT(in) :: qv, tpot, z
	980	REAL(r_k), INTENT(in) :: topo
	981	REAL(r_k), INTENT(out) :: zmla
	982
	983	! Local
	984	INTEGER :: i
	985	INTEGER :: mldlev, bllev
	986	REAL(r_k) :: dqvar, tpotmin, refdt
	987
	988	!!!!!!! Variables
	989	! qv: water vapour mixing ratio
	990	! tpot: potential temperature [K]
	991	! z: height above sea level [m]
	992	! topo: topographic height [m]
	993	! zmla: boundary layer height [m]
	994
	995	fname = 'var_zmla_generic'
	996
	997	! Pecentage of difference of mixing ratio used to determine Mixed layer depth
	998	dqvar = 0.1
	999
	1000	! MLD = Mixed layer with no substantial variation of mixing ratio /\qv < 10% ?
	1001	!PRINT *,' Mixed layer mixing ratios qv[1] lev qv[lev] dqvar% _______'
	1002	DO mldlev = 2, dz
	1003	IF (ABS(qv(mldlev)-qv(1))/qv(1) > dqvar ) EXIT
	1004	! PRINT *,qv(1), mldlev, qv(mldlev), ABS(qv(mldlev)-qv(1))/qv(1)
	1005	END DO
	1006
	1007	! Looking for the minimum potential temperature within the MLD [tpotmin = min(tpot)_0^MLD]
	1008	tpotmin = MINVAL(tpot(1:mldlev))
	1009
	1010	! Change in temperature to determine boundary layer height
	1011	refdt = 1.5
	1012
	1013	! Determine the first level where tpot > tpotmin + 1.5 K
	1014	!PRINT *,' Mixed layer tpotmin lev tpotmin[lev] dtpot _______'
	1015	DO bllev = 1, dz
	1016	IF (ABS(tpot(bllev)-tpotmin) > refdt ) EXIT
	1017	! PRINT *,tpotmin, bllev, tpot(bllev), ABS(tpot(bllev)-tpotmin)
	1018	END DO
	1019
	1020	!PRINT *,' height end MLD:', z(mldlev)
	1021	!PRINT *,' pbl height:', z(bllev)
	1022
	1023	zmla = z(bllev) - topo
	1024
	1025	RETURN
	1026
	1027	END SUBROUTINE var_zmla_generic
	1028
[1776]	1029	SUBROUTINE var_zwind(d1, u, v, z, u10, v10, sa, ca, topo, newz, unewz, vnewz)
	1030	! Subroutine to extrapolate the wind at a given height following the 'power law' methodology
	1031	! wss[newz] = wss[z1](newz/z1)*alpha
	1032	! alpha = (ln(wss[z2])-ln(wss[z1]))/(ln(z2)-ln(z1))
	1033	! AFTER: Phd Thesis:
	1034	! Benedicte Jourdier. Ressource eolienne en France metropolitaine : methodes dâevaluation du
	1035	! potentiel, variabilite et tendances. Climatologie. Ecole Doctorale Polytechnique, 2015. French
	1036	!
	1037	IMPLICIT NONE
	1038
	1039	INTEGER, INTENT(in) :: d1
	1040	REAL(r_k), DIMENSION(d1), INTENT(in) :: u,v,z
	1041	REAL(r_k), INTENT(in) :: u10, v10, topo, sa, ca, newz
	1042	REAL(r_k), INTENT(out) :: unewz, vnewz
	1043
	1044	! Local
	1045	INTEGER :: inear
	1046	REAL(r_k) :: zaground
	1047	REAL(r_k), DIMENSION(2) :: v1, v2, zz, alpha, uvnewz
	1048
	1049	!!!!!!! Variables
	1050	! u,v: vertical wind components [ms-1]
	1051	! z: height above surface [m]
	1052	! u10,v10: 10-m wind components [ms-1]
	1053	! topo: topographical height [m]
	1054	! sa, ca: local sine and cosine of map rotation [1.]
	1055	! newz: desired height above grpund of extrapolation
	1056	! unewz,vnewz: Wind compoonents at the given height [ms-1]
	1057
	1058	fname = 'var_zwind'
	1059
	1060	PRINT *,' ilev zaground newz z[ilev+1] z[ilev+2] _______'
	1061	IF (z(1) < newz ) THEN
	1062	DO inear = 1,d1-2
	1063	zaground = z(inear+2)
	1064	PRINT *, inear, z(inear), newz, z(inear+1), z(inear+2)
	1065	IF ( zaground >= newz) EXIT
	1066	END DO
	1067	ELSE
	1068	PRINT *, 1, z(1), newz, z(2), z(3), ' z(1) > newz'
	1069	inear = d1 - 2
	1070	END IF
	1071
	1072	IF (inear == d1-2) THEN
	1073	! No vertical pair of levels is below newz, using 10m wind as first value and the first level as the second
	1074	v1(1) = u10
	1075	v1(2) = v10
	1076	v2(1) = u(1)
	1077	v2(2) = v(1)
	1078	zz(1) = 10.
	1079	zz(2) = z(1)
	1080	ELSE
	1081	v1(1) = u(inear)
	1082	v1(2) = v(inear)
	1083	v2(1) = u(inear+1)
	1084	v2(2) = v(inear+1)
	1085	zz(1) = z(inear) - topo
	1086	zz(2) = z(inear+1) - topo
	1087	END IF
	1088
	1089	! Computing for each component
	1090	alpha = (LOG(ABS(v2))-LOG(ABS(v1)))/(LOG(zz(2))-LOG(zz(1)))
	1091	PRINT *,' Computing with v1:', v1, ' ms-1 v2:', v2, ' ms-1'
	1092	PRINT *,' z1:', zz(1), 'm z2:', zz(2), ' m'
	1093	PRINT *,' alhpa u:', alpha(1), ' alpha 2:', alpha(2)
	1094
	1095	uvnewz = v1(newz/zz(1))*alpha
	1096	! Earth-rotation
	1097	unewz = uvnewz(1)ca - uvnewz(2)sa
	1098	vnewz = uvnewz(1)sa + uvnewz(2)ca
	1099
	1100	PRINT *,' result vz:', uvnewz
	1101
	1102	!STOP
	1103
	1104	RETURN
	1105
	1106	END SUBROUTINE var_zwind
	1107
[772]	1108	END MODULE module_ForDiagnosticsVars

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

Download in other formats:

Original Format