Context Navigation

module_ForDiagnostics.f90 @ 2617

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

Adding:

`tws': Wet Bulb temperature after Stull, 2011

File size: 50.4 KB

Rev	Line
[770]	1	!! Fortran version of different diagnostics
	2	! L. Fita. LMD May 2016
[772]	3	! gfortran module_generic.o module_ForDiagnosticsVars.o -c module_ForDiagnostics.F90
	4	!
	5	! f2py -m module_ForDiagnostics --f90exec=/usr/bin/gfortran-4.7 -c module_generic.F90 module_ForDiagnosticsVars.F90 module_ForDiagnostics.F90
	6
[770]	7	MODULE module_ForDiagnostics
	8
[1608]	9	USE module_definitions
[770]	10	USE module_generic
[2332]	11	USE module_scientific
[772]	12	USE module_ForDiagnosticsVars
[770]	13
[772]	14	CONTAINS
[770]	15
[772]	16	!!!!!!! Calculations
[1769]	17	! compute_cape_afwa4D: Subroutine to use WRF phys/module_diag_afwa.F `buyoancy' subroutine to compute
	18	! CAPE, CIN, ZLFC, PLFC, LI
[2274]	19	! compute_cellbnds: Subroutine to compute cellboundaries using wind-staggered lon, lats as
	20	! intersection of their related parallels and meridians
[2277]	21	! compute_cellbndsreg: Subroutine to compute cellboundaries using lon, lat from a reglar lon/lat
	22	! projection as intersection of their related parallels and meridians
[1758]	23	! compute_cllmh4D3: Computation of low, medium and high cloudiness from a 4D CLDFRA and pressure being
	24	! 3rd dimension the z-dim
	25	! compute_cllmh3D3: Computation of low, medium and high cloudiness from a 3D CLDFRA and pressure being
	26	! 3rd dimension the z-dim
[772]	27	! compute_cllmh: Computation of low, medium and high cloudiness
	28	! compute_clt4D3: Computation of total cloudiness from a 4D CLDFRA being 3rd dimension the z-dim
	29	! compute_clt3D3: Computation of total cloudiness from a 3D CLDFRA being 3rd dimension the z-dim
	30	! compute_clt: Computation of total cloudiness
[1908]	31	! compute_fog_K84: Computation of fog and visibility following Kunkel, (1984)
	32	! compute_fog_RUC: Computation of fog and visibility following RUC method Smirnova, (2000)
[1909]	33	! compute_fog_FRAML50: fog and visibility following Gultepe and Milbrandt, (2010)
[2209]	34	! compute_massvertint1D: Subroutine to vertically integrate a 1D variable in eta vertical coordinates
[1795]	35	! compute_psl_ecmwf: Compute sea level pressure using ECMWF method following Mats Hamrud and Philippe Courtier [Pa]
[2209]	36	! compute_range_faces: Subroutine to compute faces [uphill, valleys, downhill] of a mountain range along a given face
[2387]	37	! compute_tws_RK[1/2/3/4]: Subroutine to compute Wet Bulb temperature of 1/2/3/4D series of values
[1769]	38	! compute_vertint1D: Subroutine to vertically integrate a 1D variable in any vertical coordinates
	39	! compute_zint4D: Subroutine to vertically integrate a 4D variable in any vertical coordinates
[1773]	40	! compute_zmla_generic4D: Subroutine to compute pbl-height following a generic method
[1776]	41	! compute_zwind4D: Subroutine to compute extrapolate the wind at a given height following the 'power law' methodology
[1784]	42	! compute_zwind_log4D: Subroutine to compute extrapolate the wind at a given height following the 'logarithmic law' methodology
[1783]	43	! compute_zwindMCO3D: Subroutine to compute extrapolate the wind at a given height following the 'power law' methodolog
[1773]	44
[772]	45	!!!
	46	! Calculations
	47	!!!
[770]	48
[772]	49	SUBROUTINE compute_cllmh4D2(cldfra4D, pres4D, cllmh4D2, d1, d2, d3, d4)
	50	! Subroutine to compute the low, medium and high cloudiness following 'newmicro.F90' from LMDZ from a 4D CLDFRA and pressure
	51	! where zdim is the 2nd dimension (thus, cldfra4D(d1,d2,d3,d4) --> cllmh(3,d1,d3,d4) 1: low, 2: medium, 3: high
	52	! It should be properly done via an 'INTERFACE', but...
[770]	53
[772]	54	IMPLICIT NONE
	55
[1141]	56	INTEGER, INTENT(in) :: d1, d2, d3, d4
[772]	57	REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in) :: cldfra4D, pres4D
	58	REAL(r_k), DIMENSION(3,d1,d3,d4), INTENT(out) :: cllmh4D2
	59
	60	! Local
	61	INTEGER :: i,j,k, zdim, Ndim
	62
[770]	63	!!!!!!! Variables
[772]	64	! cldfra4D: 4D cloud fraction values [1]
	65	! pres4D: 4D pressure values [Pa]
	66	! Ndim: number of dimensions of the input data
	67	! d[1-4]: dimensions of 'cldfra4D'
	68	! zdim: number of the vertical-dimension within the matrix
	69	! cltlmh4D2: low, medium, high cloudiness for the 4D cldfra and d2 being 'zdim'
[770]	70
[772]	71	fname = 'compute_cllmh4D2'
	72	zdim = 2
	73	Ndim = 4
[770]	74
[772]	75	DO i=1, d1
	76	DO j=1, d3
	77	DO k=1, d4
	78	cllmh4D2(:,i,j,k) = var_cllmh(cldfra4D(i,:,j,k), pres4D(i,:,j,k), d2)
	79	END DO
	80	END DO
	81	END DO
	82
	83	RETURN
[770]	84
[772]	85	END SUBROUTINE compute_cllmh4D2
[770]	86
[772]	87	SUBROUTINE compute_cllmh3D1(cldfra3D, pres3D, cllmh3D1, d1, d2, d3)
	88	! Subroutine to compute the low, medium and high cloudiness following 'newmicro.F90' from LMDZ from a 3D CLDFRA and pressure
	89	! where zdim is the 1st dimension (thus, cldfra3D(d1,d2,d3) --> cllmh(3,d2,d3) 1: low, 2: medium, 3: high
	90	! It should be properly done via an 'INTERFACE', but...
[770]	91
[772]	92	IMPLICIT NONE
[770]	93
[1141]	94	INTEGER, INTENT(in) :: d1, d2, d3
[772]	95	REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in) :: cldfra3D, pres3D
	96	REAL(r_k), DIMENSION(3,d2,d3), INTENT(out) :: cllmh3D1
	97
	98	! Local
	99	INTEGER :: i,j,k, zdim, Ndim
	100
	101	!!!!!!! Variables
	102	! cldfra3D: 3D cloud fraction values [1]
	103	! pres3D: 3D pressure values [Pa]
	104	! Ndim: number of dimensions of the input data
	105	! d[1-3]: dimensions of 'cldfra3D'
	106	! zdim: number of the vertical-dimension within the matrix
	107	! cltlmh3D1: low, medium, high cloudiness for the 3D cldfra and d1 being 'zdim'
	108
	109	fname = 'compute_cllmh3D1'
	110	zdim = 1
	111	Ndim = 3
	112
	113	DO i=1, d1
	114	DO j=1, d2
	115	cllmh3D1(:,i,j) = var_cllmh(cldfra3D(:,i,j), pres3D(:,i,j), d1)
	116	END DO
	117	END DO
	118
	119	RETURN
	120
	121	END SUBROUTINE compute_cllmh3D1
	122
	123	SUBROUTINE compute_cllmh(cldfra1D, cldfra2D, cldfra3D, cldfra4D, pres1D, pres2D, pres3D, pres4D, &
	124	Ndim, zdim, cllmh1D, cllmh2D1, cllmh2D2, cllmh3D1, cllmh3D2, cllmh3D3, cllmh4D1, cllmh4D2, &
	125	cllmh4D3, cllmh4D4, d1, d2, d3, d4)
	126	! Subroutine to compute the low, medium and high cloudiness following 'newmicro.F90' from LMDZ
	127
[770]	128	IMPLICIT NONE
	129
[1141]	130	INTEGER, INTENT(in) :: Ndim, d1, d2, d3, d4, zdim
[772]	131	REAL(r_k), DIMENSION(d1), OPTIONAL, INTENT(in) :: cldfra1D, pres1D
	132	REAL(r_k), DIMENSION(d1,d2), OPTIONAL, INTENT(in) :: cldfra2D, pres2D
	133	REAL(r_k), DIMENSION(d1,d2,d3), OPTIONAL, INTENT(in) :: cldfra3D, pres3D
	134	REAL(r_k), DIMENSION(d1,d2,d3,d4), OPTIONAL, &
	135	INTENT(in) :: cldfra4D, pres4D
	136	REAL(r_k), DIMENSION(3), OPTIONAL, INTENT(out) :: cllmh1D
	137	REAL(r_k), DIMENSION(d1,3), OPTIONAL, INTENT(out) :: cllmh2D1
	138	REAL(r_k), DIMENSION(d2,3), OPTIONAL, INTENT(out) :: cllmh2D2
	139	REAL(r_k), DIMENSION(d2,d3,3), OPTIONAL, INTENT(out) :: cllmh3D1
	140	REAL(r_k), DIMENSION(d1,d3,3), OPTIONAL, INTENT(out) :: cllmh3D2
	141	REAL(r_k), DIMENSION(d1,d2,3), OPTIONAL, INTENT(out) :: cllmh3D3
	142	REAL(r_k), DIMENSION(d2,d3,d4,3), OPTIONAL, &
	143	INTENT(out) :: cllmh4D1
	144	REAL(r_k), DIMENSION(d1,d3,d4,3), OPTIONAL, &
	145	INTENT(out) :: cllmh4D2
	146	REAL(r_k), DIMENSION(d1,d2,d4,3), OPTIONAL, &
	147	INTENT(out) :: cllmh4D3
	148	REAL(r_k), DIMENSION(d1,d2,d3,3), OPTIONAL, &
	149	INTENT(out) :: cllmh4D4
[770]	150
	151	! Local
[772]	152	INTEGER :: i,j,k
[770]	153
	154	!!!!!!! Variables
[772]	155	! cldfra[1-4]D: cloud fraction values [1]
	156	! pres[1-4]D: pressure values [Pa]
	157	! Ndim: number of dimensions of the input data
	158	! d[1-4]: dimensions of 'cldfra'
	159	! zdim: number of the vertical-dimension within the matrix
	160	! cllmh1D: low, medium and high cloudiness for the 1D cldfra
	161	! cllmh2D1: low, medium and high cloudiness for the 2D cldfra and d1 being 'zdim'
	162	! cllmh2D2: low, medium and high cloudiness for the 2D cldfra and d2 being 'zdim'
	163	! cllmh3D1: low, medium and high cloudiness for the 3D cldfra and d1 being 'zdim'
	164	! cllmh3D2: low, medium and high cloudiness for the 3D cldfra and d2 being 'zdim'
	165	! cllmh3D3: low, medium and high cloudiness for the 3D cldfra and d3 being 'zdim'
	166	! cllmh4D1: low, medium and high cloudiness for the 4D cldfra and d1 being 'zdim'
	167	! cllmh4D2: low, medium and high cloudiness for the 4D cldfra and d2 being 'zdim'
	168	! cllmh4D3: low, medium and high cloudiness for the 4D cldfra and d3 being 'zdim'
	169	! cllmh4D4: low, medium and high cloudiness for the 4D cldfra and d4 being 'zdim'
[770]	170
[772]	171	fname = 'compute_cllmh'
[770]	172
[772]	173	SELECT CASE (Ndim)
	174	CASE (1)
	175	cllmh1D = var_cllmh(cldfra1D, pres1D, d1)
	176	CASE (2)
	177	IF (zdim == 1) THEN
	178	DO i=1, d2
	179	cllmh2D1(i,:) = var_cllmh(cldfra2D(:,i), pres2D(:,i), d1)
	180	END DO
	181	ELSE IF (zdim == 2) THEN
	182	DO i=1, d1
	183	cllmh2D2(i,:) = var_cllmh(cldfra2D(:,i), pres2D(i,:), d2)
	184	END DO
	185	ELSE
	186	PRINT *,TRIM(ErrWarnMsg('err'))
	187	PRINT *,' ' // TRIM(fname) // ': wrong zdim:', zdim,' for Ndim=', Ndim, ' !!'
	188	PRINT *,' accepted values: 1,2'
	189	STOP
	190	END IF
	191	CASE (3)
	192	IF (zdim == 1) THEN
	193	DO i=1, d2
	194	DO j=1, d3
	195	cllmh3D1(i,j,:) = var_cllmh(cldfra3D(:,i,j), pres3D(:,i,j), d1)
	196	END DO
	197	END DO
	198	ELSE IF (zdim == 2) THEN
	199	DO i=1, d1
	200	DO j=1, d3
	201	cllmh3D2(i,j,:) = var_cllmh(cldfra3D(i,:,j), pres3D(i,:,j), d2)
	202	END DO
	203	END DO
	204	ELSE IF (zdim == 3) THEN
	205	DO i=1, d1
	206	DO j=1, d2
	207	cllmh3D3(i,j,:) = var_cllmh(cldfra3D(i,j,:), pres3D(i,j,:), d3)
	208	END DO
	209	END DO
	210	ELSE
	211	PRINT *,TRIM(ErrWarnMsg('err'))
	212	PRINT *,' ' // TRIM(fname) // ': wrong zdim:', zdim,' for Ndim=', Ndim, ' !!'
	213	PRINT *,' accepted values: 1,2,3'
	214	STOP
	215	END IF
	216	CASE (4)
	217	IF (zdim == 1) THEN
	218	DO i=1, d2
	219	DO j=1, d3
	220	DO k=1, d4
	221	cllmh4D1(i,j,k,:) = var_cllmh(cldfra4D(:,i,j,k), pres4D(:,i,j,k), d1)
	222	END DO
	223	END DO
	224	END DO
	225	ELSE IF (zdim == 2) THEN
	226	DO i=1, d1
	227	DO j=1, d3
	228	DO k=1, d4
	229	cllmh4D2(i,j,k,:) = var_cllmh(cldfra4D(i,:,j,k), pres4D(i,:,j,k), d2)
	230	END DO
	231	END DO
	232	END DO
	233	ELSE IF (zdim == 3) THEN
	234	DO i=1, d2
	235	DO j=1, d3
	236	DO k=1, d4
	237	cllmh4D3(i,j,k,:) = var_cllmh(cldfra4D(i,j,:,k), pres4D(i,j,:,k), d3)
	238	END DO
	239	END DO
	240	END DO
	241	ELSE IF (zdim == 4) THEN
	242	DO i=1, d1
	243	DO j=1, d2
	244	DO k=1, d3
	245	cllmh4D4(i,j,k,:) = var_cllmh(cldfra4D(i,j,k,:), pres4D(i,j,k,:), d4)
	246	END DO
	247	END DO
	248	END DO
	249	ELSE
	250	PRINT *,TRIM(ErrWarnMsg('err'))
	251	PRINT *,' ' // TRIM(fname) // ': wrong zdim:', zdim,' for Ndim=', Ndim, ' !!'
	252	PRINT *,' accepted values: 1,2,3,4'
	253	STOP
	254	END IF
	255	CASE DEFAULT
	256	PRINT *,TRIM(ErrWarnMsg('err'))
	257	PRINT *,' ' // TRIM(fname) // ': Ndim:', Ndim,' not ready !!'
	258	STOP
	259	END SELECT
[770]	260
	261	RETURN
	262
[772]	263	END SUBROUTINE compute_cllmh
[770]	264
[772]	265	SUBROUTINE compute_clt4D2(cldfra4D, clt4D2, d1, d2, d3, d4)
	266	! Subroutine to compute the total cloudiness following 'newmicro.F90' from LMDZ from a 4D CLDFRA
	267	! where zdim is the 2nd dimension (thus, cldfra4D(d1,d2,d3,d4) --> clt(d1,d3,d4)
	268	! It should be properly done via an 'INTERFACE', but...
[770]	269
	270	IMPLICIT NONE
	271
[1141]	272	INTEGER, INTENT(in) :: d1, d2, d3, d4
[772]	273	REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in) :: cldfra4D
	274	REAL(r_k), DIMENSION(d1,d3,d4), INTENT(out) :: clt4D2
	275
[770]	276	! Local
[772]	277	INTEGER :: i,j,k, zdim, Ndim
	278
[770]	279	!!!!!!! Variables
[772]	280	! cldfra4D: 4D cloud fraction values [1]
	281	! Ndim: number of dimensions of the input data
	282	! d[1-4]: dimensions of 'cldfra4D'
	283	! zdim: number of the vertical-dimension within the matrix
	284	! clt4D2: total cloudiness for the 4D cldfra and d2 being 'zdim'
[770]	285
[772]	286	fname = 'compute_clt4D2'
	287	zdim = 2
	288	Ndim = 4
[770]	289
[772]	290	DO i=1, d1
	291	DO j=1, d3
	292	DO k=1, d4
	293	clt4D2(i,j,k) = var_clt(cldfra4D(i,:,j,k), d2)
	294	END DO
	295	END DO
[770]	296	END DO
[772]	297
	298	RETURN
[770]	299
[772]	300	END SUBROUTINE compute_clt4D2
[770]	301
[772]	302	SUBROUTINE compute_clt3D1(cldfra3D, clt3D1, d1, d2, d3)
	303	! Subroutine to compute the total cloudiness following 'newmicro.F90' from LMDZ from a 3D CLDFRA
	304	! where zdim is the 1st dimension (thus, cldfra4D(d1,d2,d3) --> clt(d2,d3)
	305	! It should be properly done via an 'INTERFACE', but...
[770]	306
[772]	307	IMPLICIT NONE
[770]	308
[1141]	309	INTEGER, INTENT(in) :: d1, d2, d3
[772]	310	REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in) :: cldfra3D
	311	REAL(r_k), DIMENSION(d2,d3), INTENT(out) :: clt3D1
[770]	312
[772]	313	! Local
	314	INTEGER :: i,j,k, zdim, Ndim
	315
	316	!!!!!!! Variables
	317	! cldfra3D: 3D cloud fraction values [1]
	318	! Ndim: number of dimensions of the input data
	319	! d[1-3]: dimensions of 'cldfra3D'
	320	! zdim: number of the vertical-dimension within the matrix
	321	! clt3D1: total cloudiness for the 3D cldfra and d1 being 'zdim'
	322
	323	fname = 'compute_clt3D1'
	324	zdim = 1
	325	Ndim = 3
	326
	327	DO i=1, d2
	328	DO j=1, d3
	329	clt3D1(i,j) = var_clt(cldfra3D(:,i,j), d1)
	330	END DO
	331	END DO
	332
	333	RETURN
	334
	335	END SUBROUTINE compute_clt3D1
	336
	337	SUBROUTINE compute_clt(cldfra1D, cldfra2D, cldfra3D, cldfra4D, Ndim, zdim, clt1D, clt2D1, clt2D2, &
	338	clt3D1, clt3D2, clt3D3, clt4D1, clt4D2, clt4D3, clt4D4, d1, d2, d3, d4)
	339	! Subroutine to compute the total cloudiness following 'newmicro.F90' from LMDZ
	340
[770]	341	IMPLICIT NONE
	342
[1141]	343	INTEGER, INTENT(in) :: Ndim, d1, d2, d3, d4, zdim
[770]	344	REAL(r_k), DIMENSION(d1), OPTIONAL, INTENT(in) :: cldfra1D
	345	REAL(r_k), DIMENSION(d1,d2), OPTIONAL, INTENT(in) :: cldfra2D
	346	REAL(r_k), DIMENSION(d1,d2,d3), OPTIONAL, INTENT(in) :: cldfra3D
	347	REAL(r_k), DIMENSION(d1,d2,d3,d4), OPTIONAL, &
	348	INTENT(in) :: cldfra4D
	349	REAL(r_k), OPTIONAL, INTENT(out) :: clt1D
	350	REAL(r_k), DIMENSION(d1), OPTIONAL, INTENT(out) :: clt2D1
	351	REAL(r_k), DIMENSION(d2), OPTIONAL, INTENT(out) :: clt2D2
	352	REAL(r_k), DIMENSION(d2,d3), OPTIONAL, INTENT(out) :: clt3D1
	353	REAL(r_k), DIMENSION(d1,d3), OPTIONAL, INTENT(out) :: clt3D2
	354	REAL(r_k), DIMENSION(d1,d2), OPTIONAL, INTENT(out) :: clt3D3
	355	REAL(r_k), DIMENSION(d2,d3,d4), OPTIONAL,INTENT(out) :: clt4D1
	356	REAL(r_k), DIMENSION(d1,d3,d4), OPTIONAL,INTENT(out) :: clt4D2
	357	REAL(r_k), DIMENSION(d1,d2,d4), OPTIONAL,INTENT(out) :: clt4D3
	358	REAL(r_k), DIMENSION(d1,d2,d3), OPTIONAL,INTENT(out) :: clt4D4
	359
	360	! Local
	361	INTEGER :: i,j,k
	362
	363	!!!!!!! Variables
	364	! cldfra[1-4]D: cloud fraction values [1]
	365	! Ndim: number of dimensions of the input data
	366	! d[1-4]: dimensions of 'cldfra'
	367	! zdim: number of the vertical-dimension within the matrix
	368	! clt1D: total cloudiness for the 1D cldfra
	369	! clt2D1: total cloudiness for the 2D cldfra and d1 being 'zdim'
	370	! clt2D2: total cloudiness for the 2D cldfra and d2 being 'zdim'
	371	! clt3D1: total cloudiness for the 3D cldfra and d1 being 'zdim'
	372	! clt3D2: total cloudiness for the 3D cldfra and d2 being 'zdim'
	373	! clt3D3: total cloudiness for the 3D cldfra and d3 being 'zdim'
	374	! clt4D1: total cloudiness for the 4D cldfra and d1 being 'zdim'
	375	! clt4D2: total cloudiness for the 4D cldfra and d2 being 'zdim'
	376	! clt4D3: total cloudiness for the 4D cldfra and d3 being 'zdim'
	377	! clt4D4: total cloudiness for the 4D cldfra and d4 being 'zdim'
	378
	379	fname = 'compute_clt'
	380
	381	SELECT CASE (Ndim)
	382	CASE (1)
	383	clt1D = var_clt(cldfra1D, d1)
	384	CASE (2)
	385	IF (zdim == 1) THEN
	386	DO i=1, d2
	387	clt2D1(i) = var_clt(cldfra2D(:,i), d1)
	388	END DO
	389	ELSE IF (zdim == 2) THEN
	390	DO i=1, d1
	391	clt2D2(i) = var_clt(cldfra2D(:,i), d2)
	392	END DO
	393	ELSE
	394	PRINT *,TRIM(ErrWarnMsg('err'))
	395	PRINT *,' ' // TRIM(fname) // ': wrong zdim:', zdim,' for Ndim=', Ndim, ' !!'
	396	PRINT *,' accepted values: 1,2'
	397	STOP
	398	END IF
	399	CASE (3)
	400	IF (zdim == 1) THEN
	401	DO i=1, d2
	402	DO j=1, d3
	403	clt3D1(i,j) = var_clt(cldfra3D(:,i,j), d1)
	404	END DO
	405	END DO
	406	ELSE IF (zdim == 2) THEN
	407	DO i=1, d1
	408	DO j=1, d3
	409	clt3D2(i,j) = var_clt(cldfra3D(i,:,j), d2)
	410	END DO
	411	END DO
	412	ELSE IF (zdim == 3) THEN
	413	DO i=1, d1
	414	DO j=1, d2
	415	clt3D3(i,j) = var_clt(cldfra3D(i,j,:), d3)
	416	END DO
	417	END DO
	418	ELSE
	419	PRINT *,TRIM(ErrWarnMsg('err'))
	420	PRINT *,' ' // TRIM(fname) // ': wrong zdim:', zdim,' for Ndim=', Ndim, ' !!'
	421	PRINT *,' accepted values: 1,2,3'
	422	STOP
	423	END IF
	424	CASE (4)
	425	IF (zdim == 1) THEN
	426	DO i=1, d2
	427	DO j=1, d3
	428	DO k=1, d4
	429	clt4D1(i,j,k) = var_clt(cldfra4D(:,i,j,k), d1)
	430	END DO
	431	END DO
	432	END DO
	433	ELSE IF (zdim == 2) THEN
	434	DO i=1, d1
	435	DO j=1, d3
	436	DO k=1, d4
	437	clt4D2(i,j,k) = var_clt(cldfra4D(i,:,j,k), d2)
	438	END DO
	439	END DO
	440	END DO
	441	ELSE IF (zdim == 3) THEN
	442	DO i=1, d2
	443	DO j=1, d3
	444	DO k=1, d4
	445	clt4D3(i,j,k) = var_clt(cldfra4D(i,j,:,k), d3)
	446	END DO
	447	END DO
	448	END DO
	449	ELSE IF (zdim == 4) THEN
	450	DO i=1, d1
	451	DO j=1, d2
	452	DO k=1, d3
	453	clt4D4(i,j,k) = var_clt(cldfra4D(i,j,k,:), d4)
	454	END DO
	455	END DO
	456	END DO
	457	ELSE
	458	PRINT *,TRIM(ErrWarnMsg('err'))
	459	PRINT *,' ' // TRIM(fname) // ': wrong zdim:', zdim,' for Ndim=', Ndim, ' !!'
	460	PRINT *,' accepted values: 1,2,3,4'
	461	STOP
	462	END IF
	463	CASE DEFAULT
	464	PRINT *,TRIM(ErrWarnMsg('err'))
	465	PRINT *,' ' // TRIM(fname) // ': Ndim:', Ndim,' not ready !!'
	466	STOP
	467	END SELECT
	468
	469	RETURN
	470
	471	END SUBROUTINE compute_clt
	472
[1762]	473	SUBROUTINE compute_massvertint1D(var, mutot, dz, deta, integral)
	474	! Subroutine to vertically integrate a 1D variable in eta vertical coordinates
	475
	476	IMPLICIT NONE
	477
	478	INTEGER, INTENT(in) :: dz
	479	REAL(r_k), INTENT(in) :: mutot
	480	REAL(r_k), DIMENSION(dz), INTENT(in) :: var, deta
	481	REAL(r_k), INTENT(out) :: integral
	482
	483	! Local
	484	INTEGER :: k
	485
	486	!!!!!!! Variables
	487	! var: vertical variable to integrate (assuming kgkg-1)
	488	! mutot: total dry-air mass in column
	489	! dz: vertical dimension
	490	! deta: eta-levels difference between full eta-layers
	491
	492	fname = 'compute_massvertint1D'
	493
	494	! integral=0.
	495	! DO k=1,dz
	496	! integral = integral + var(k)*deta(k)
	497	! END DO
	498	integral = SUM(var*deta)
	499
	500	integral=integral*mutot/g
	501
	502	RETURN
	503
	504	END SUBROUTINE compute_massvertint1D
	505
	506	SUBROUTINE compute_zint4D(var4D, dlev, zweight, d1, d2, d3, d4, int3D)
	507	! Subroutine to vertically integrate a 4D variable in any vertical coordinates
	508
	509	IMPLICIT NONE
	510
	511	INTEGER, INTENT(in) :: d1,d2,d3,d4
	512	REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in) :: var4D, dlev, zweight
	513	REAL(r_k), DIMENSION(d1,d2,d4), INTENT(out) :: int3D
	514
	515	! Local
	516	INTEGER :: i,j,l
	517
	518	!!!!!!! Variables
	519	! var4D: vertical variable to integrate
	520	! dlev: height of layers
	521	! zweight: weight for each level to be applied (=1. for no effect)
	522
	523	fname = 'compute_zint4D'
	524
	525	DO i=1,d1
	526	DO j=1,d2
	527	DO l=1,d4
	528	CALL compute_vertint1D(var4D(i,j,:,l),d3, dlev(i,j,:,l), zweight(i,j,:,l), &
	529	int3D(i,j,l))
	530	END DO
	531	END DO
	532	END DO
	533
	534	RETURN
	535
	536	END SUBROUTINE compute_zint4D
	537
	538	SUBROUTINE compute_vertint1D(var, dz, deta, zweight, integral)
	539	! Subroutine to vertically integrate a 1D variable in any vertical coordinates
	540
	541	IMPLICIT NONE
	542
	543	INTEGER, INTENT(in) :: dz
	544	REAL(r_k), DIMENSION(dz), INTENT(in) :: var, deta, zweight
	545	REAL(r_k), INTENT(out) :: integral
	546
	547	! Local
	548	INTEGER :: k
	549
	550	!!!!!!! Variables
	551	! var: vertical variable to integrate
	552	! dz: vertical dimension
	553	! deta: eta-levels difference between layers
	554	! zweight: weight for each level to be applied (=1. for no effect)
	555
	556	fname = 'compute_vertint1D'
	557
	558	! integral=0.
	559	! DO k=1,dz
	560	! integral = integral + var(k)*deta(k)
	561	! END DO
	562	integral = SUM(vardetazweight)
	563
	564	RETURN
	565
	566	END SUBROUTINE compute_vertint1D
	567
[1759]	568	SUBROUTINE compute_cape_afwa4D(ta, hur, press, zg, hgt, cape, cin, zlfc, plfc, li, parcelmethod, &
	569	d1, d2, d3, d4)
	570	! Subroutine to use WRF phys/module_diag_afwa.F `buyoancy' subroutine to compute CAPE, CIN, ZLFC, PLFC, LI
	571
	572	IMPLICIT NONE
	573
	574	INTEGER, INTENT(in) :: d1, d2, d3, d4, parcelmethod
	575	REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in) :: ta, hur, press, zg
	576	REAL(r_k), DIMENSION(d1,d2), INTENT(in) :: hgt
	577	REAL(r_k), DIMENSION(d1,d2,d4), INTENT(out) :: cape, cin, zlfc, plfc, li
	578
	579	! Local
	580	INTEGER :: i, j, it
	581	INTEGER :: ofunc
	582
	583	!!!!!!! Variables
	584	! ta: air temperature [K]
	585	! hur: relative humidity [%]
	586	! press: air pressure [Pa]
	587	! zg: geopotential height [gpm]
	588	! hgt: topographical height [m]
	589	! cape: Convective available potential energy [Jkg-1]
	590	! cin: Convective inhibition [Jkg-1]
	591	! zlfc: height at the Level of free convection [m]
	592	! plfc: pressure at the Level of free convection [Pa]
	593	! li: lifted index [1]
	594	! parcelmethod:
	595	! Most Unstable = 1 (default)
	596	! Mean layer = 2
	597	! Surface based = 3
	598
	599	fname = 'compute_cape_afwa4D'
	600
	601	DO i=1, d1
	602	DO j=1, d2
	603	DO it=1, d4
	604	ofunc = var_cape_afwa1D(d3, ta(i,j,:,it), hur(i,j,:,it), press(i,j,:,it), zg(i,j,:,it), &
	605	1, cape(i,j,it), cin(i,j,it), zlfc(i,j,it), plfc(i,j,it), li(i,j,it), parcelmethod)
[1833]	606	IF (zlfc(i,j,it) /= -1.) zlfc(i,j,it) = zlfc(i,j,it) - hgt(i,j)
[1759]	607	END DO
	608	END DO
	609	END DO
	610
	611	RETURN
	612
	613	END SUBROUTINE compute_cape_afwa4D
	614
[1795]	615	SUBROUTINE compute_psl_ecmwf(ps, hgt, T, press, unpress, psl, d1, d2, d4)
	616	! Subroutine to compute sea level pressure using ECMWF method following Mats Hamrud and Philippe Courtier [Pa]
	617
	618	IMPLICIT NONE
	619
	620	INTEGER, INTENT(in) :: d1, d2, d4
	621	REAL(r_k), DIMENSION(d1,d2,d4), INTENT(in) :: ps, T, press, unpress
	622	REAL(r_k), DIMENSION(d1,d2), INTENT(in) :: hgt
	623	REAL(r_k), DIMENSION(d1,d2,d4), INTENT(out) :: psl
	624
	625	! Local
	626	INTEGER :: i, j, it
	627
	628	!!!!!!! Variables
	629	! ps: surface pressure [Pa]
	630	! hgt: terrain height [m]
	631	! T: temperature at first half-mass level [K]
	632	! press: pressure at first full levels [Pa]
	633	! unpress: pressure at first mass (half) levels [Pa]
	634	! psl: sea-level pressure [Pa]
	635
	636	fname = 'compute_psl_ecmwf'
	637
	638	DO i=1, d1
	639	DO j=1, d2
	640	DO it=1, d4
	641	CALL var_psl_ecmwf(ps(i,j,it), hgt(i,j), T(i,j,it), unpress(i,j,it), press(i,j,it), &
	642	psl(i,j,it))
	643	END DO
	644	END DO
	645	END DO
	646
	647	RETURN
	648
	649	END SUBROUTINE compute_psl_ecmwf
	650
[1773]	651	SUBROUTINE compute_zmla_generic4D(tpot, qratio, z, hgt, zmla3D, d1, d2, d3, d4)
	652	! Subroutine to compute pbl-height following a generic method
	653	! from Nielsen-Gammon et al., 2008 J. Appl. Meteor. Clim.
	654	! applied also in Garcia-Diez et al., 2013, QJRMS
	655	! where
	656	! "The technique identifies the ML height as a threshold increase of potential temperature from
	657	! its minimum value within the boundary layer."
	658	! here applied similarly to Garcia-Diez et al. where
	659	! zmla = "...first level where potential temperature exceeds the minimum potential temperature
	660	! reached in the mixed layer by more than 1.5 K"
	661
	662	IMPLICIT NONE
	663
	664	INTEGER, INTENT(in) :: d1, d2, d3, d4
	665	REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in) :: tpot, qratio, z
	666	REAL(r_k), DIMENSION(d1,d2), INTENT(in) :: hgt
	667	REAL(r_k), DIMENSION(d1,d2,d4), INTENT(out) :: zmla3D
	668
	669	! Local
	670	INTEGER :: i, j, it
	671
	672	!!!!!!! Variables
	673	! tpot: potential air temperature [K]
	674	! qratio: water vapour mixing ratio [kgkg-1]
	675	! z: height above sea level [m]
	676	! hgt: terrain height [m]
	677	! zmla3D: boundary layer height from surface [m]
	678
	679	fname = 'compute_zmla_generic4D'
	680
	681	DO i=1, d1
	682	DO j=1, d2
	683	DO it=1, d4
	684	CALL var_zmla_generic(d3, qratio(i,j,:,it), tpot(i,j,:,it), z(i,j,:,it), hgt(i,j), &
	685	zmla3D(i,j,it))
	686	END DO
	687	END DO
	688	END DO
	689
	690	RETURN
	691
	692	END SUBROUTINE compute_zmla_generic4D
	693
[1777]	694	SUBROUTINE compute_zwind4D(ua, va, z, uas, vas, sina, cosa, zextrap, uaz, vaz, d1, d2, d3, d4)
[1776]	695	! Subroutine to compute extrapolate the wind at a given height following the 'power law' methodology
[1773]	696
[1776]	697	IMPLICIT NONE
	698
	699	INTEGER, INTENT(in) :: d1, d2, d3, d4
[1777]	700	REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in) :: ua, va, z
[1776]	701	REAL(r_k), DIMENSION(d1,d2,d4), INTENT(in) :: uas, vas
[1777]	702	REAL(r_k), DIMENSION(d1,d2), INTENT(in) :: sina, cosa
[1776]	703	REAL(r_k), INTENT(in) :: zextrap
	704	REAL(r_k), DIMENSION(d1,d2,d4), INTENT(out) :: uaz, vaz
	705
	706	! Local
	707	INTEGER :: i, j, it
	708
	709	!!!!!!! Variables
	710	! tpot: potential air temperature [K]
	711	! qratio: water vapour mixing ratio [kgkg-1]
[1777]	712	! z: height above surface [m]
[1776]	713	! sina, cosa: local sine and cosine of map rotation [1.]
	714	! zmla3D: boundary layer height from surface [m]
	715
	716	fname = 'compute_zwind4D'
	717
	718	DO i=1, d1
	719	DO j=1, d2
	720	DO it=1, d4
[1777]	721	CALL var_zwind(d3, ua(i,j,:,it), va(i,j,:,it), z(i,j,:,it), uas(i,j,it), vas(i,j,it), &
	722	sina(i,j), cosa(i,j), zextrap, uaz(i,j,it), vaz(i,j,it))
[1776]	723	END DO
	724	END DO
	725	END DO
	726
	727	RETURN
	728
	729	END SUBROUTINE compute_zwind4D
	730
[1784]	731	SUBROUTINE compute_zwind_log4D(ua, va, z, uas, vas, sina, cosa, zextrap, uaz, vaz, d1, d2, d3, d4)
	732	! Subroutine to compute extrapolate the wind at a given height following the 'logarithmic law' methodology
	733
	734	IMPLICIT NONE
	735
	736	INTEGER, INTENT(in) :: d1, d2, d3, d4
	737	REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in) :: ua, va, z
	738	REAL(r_k), DIMENSION(d1,d2,d4), INTENT(in) :: uas, vas
	739	REAL(r_k), DIMENSION(d1,d2), INTENT(in) :: sina, cosa
	740	REAL(r_k), INTENT(in) :: zextrap
	741	REAL(r_k), DIMENSION(d1,d2,d4), INTENT(out) :: uaz, vaz
	742
	743	! Local
	744	INTEGER :: i, j, it
	745
	746	!!!!!!! Variables
	747	! tpot: potential air temperature [K]
	748	! qratio: water vapour mixing ratio [kgkg-1]
	749	! z: height above surface [m]
	750	! sina, cosa: local sine and cosine of map rotation [1.]
	751	! zmla3D: boundary layer height from surface [m]
	752
	753	fname = 'compute_zwind_log4D'
	754
	755	DO i=1, d1
	756	DO j=1, d2
	757	DO it=1, d4
	758	CALL var_zwind_log(d3, ua(i,j,:,it), va(i,j,:,it), z(i,j,:,it), uas(i,j,it), vas(i,j,it), &
	759	sina(i,j), cosa(i,j), zextrap, uaz(i,j,it), vaz(i,j,it))
	760	END DO
	761	END DO
	762	END DO
	763
	764	RETURN
	765
	766	END SUBROUTINE compute_zwind_log4D
	767
[1783]	768	SUBROUTINE compute_zwindMO3D(d1, d2, d3, ust, znt, rmol, uas, vas, sina, cosa, newz, uznew, vznew)
	769	! Subroutine to compute extrapolate the wind at a given height following the 'power law' methodology
[1784]	770	! NOTE: only usefull for newz < 80. m
[1783]	771
	772	IMPLICIT NONE
	773
	774	INTEGER, INTENT(in) :: d1, d2, d3
	775	REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in) :: ust, znt, rmol
	776	REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in) :: uas, vas
	777	REAL(r_k), DIMENSION(d1,d2), INTENT(in) :: sina, cosa
	778	REAL(r_k), INTENT(in) :: newz
	779	REAL(r_k), DIMENSION(d1,d2,d3), INTENT(out) :: uznew, vznew
	780
	781	! Local
	782	INTEGER :: i, j, it
	783
	784	!!!!!!! Variables
	785	! ust: u* in similarity theory [ms-1]
	786	! znt: thermal time-varying roughness length [m]
	787	! rmol: Inverse of the Obukhov length [m-1]
	788	! uas: x-component 10-m wind speed [ms-1]
	789	! vas: y-component 10-m wind speed [ms-1]
	790	! sina, cosa: local sine and cosine of map rotation [1.]
	791
	792	fname = 'compute_zwindMO3D'
	793
	794	DO i=1, d1
	795	DO j=1, d2
	796	DO it=1, d3
	797	CALL var_zwind_MOtheor(ust(i,j,it), znt(i,j,it), rmol(i,j,it), uas(i,j,it), vas(i,j,it), &
	798	sina(i,j), cosa(i,j), newz, uznew(i,j,it), vznew(i,j,it))
	799	END DO
	800	END DO
	801	END DO
	802
	803	RETURN
	804
	805	END SUBROUTINE compute_zwindMO3D
	806
[1804]	807	SUBROUTINE compute_potevap_orPM3D(d1, d2, d3, rho1, ust, uas, vas, tas, ps, qv1, potevap)
	808	! Subroutine to compute potential evapotranspiration Penman-Monteith formulation implemented in
[1833]	809	! ORCHIDEE in src_sechiba/enerbil.f90
[1804]	810
	811	IMPLICIT NONE
	812
	813	INTEGER, INTENT(in) :: d1, d2, d3
	814	REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in) :: rho1, ust, uas, vas, tas, ps, qv1
	815	REAL(r_k), DIMENSION(d1,d2,d3), INTENT(out) :: potevap
	816
	817	! Local
	818	INTEGER :: i, j, it
	819
	820	!!!!!!! Variables
	821	! rho1: atsmophere density at the first layer [kgm-3]
	822	! ust: u* in similarity theory [ms-1]
	823	! uas: x-component 10-m wind speed [ms-1]
	824	! vas: y-component 10-m wind speed [ms-1]
	825	! tas: 2-m atmosphere temperature [K]
	826	! ps: surface pressure [Pa]
	827	! qv1: 1st layer atmospheric mixing ratio [kgkg-1]
	828	! potevap: potential evapo transpiration [kgm-2s-1]
	829
	830	fname = 'compute_potevap_orPM3D'
	831
	832	DO i=1, d1
	833	DO j=1, d2
	834	DO it=1, d3
	835	CALL var_potevap_orPM(rho1(i,j,it), ust(i,j,it), uas(i,j,it), vas(i,j,it), tas(i,j,it), &
	836	ps(i,j,it), qv1(i,j,it), potevap(i,j,it))
	837	END DO
	838	END DO
	839	END DO
	840
	841	RETURN
	842
	843	END SUBROUTINE compute_potevap_orPM3D
	844
[1908]	845	SUBROUTINE compute_fog_K84(d1, d2, d3, qc, qi, fog, vis)
	846	! Subroutine to compute fog: qcloud + qice /= 0.
	847	! And visibility following Kunkel, B. A., (1984): Parameterization of droplet terminal velocity and
[2387]	848	! extinction coefficient in fog models. J. Climate Appl. Meteor., 23, 34-41.
[1908]	849
	850	IMPLICIT NONE
	851
	852	INTEGER, INTENT(in) :: d1, d2, d3
	853	REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in) :: qc, qi
	854	INTEGER, DIMENSION(d1,d2,d3), INTENT(out) :: fog
	855	REAL(r_k), DIMENSION(d1,d2,d3), INTENT(out) :: vis
	856
	857	! Local
	858	INTEGER :: i, j, it
	859
	860	!!!!!!! Variables
	861	! qc: cloud mixing ratio [kgkg-1]
	862	! qi, ice mixing ratio [kgkg-1]
	863	! fog: presence of fog (1: yes, 0: no)
	864	! vis: visibility within fog [km]
	865
	866	fname = 'compute_fog_K84'
	867
	868	DO i=1, d1
	869	DO j=1, d2
	870	DO it=1, d3
	871	CALL var_fog_K84(qc(i,j,it), qi(i,j,it), fog(i,j,it), vis(i,j,it))
	872	END DO
	873	END DO
	874	END DO
	875
	876	RETURN
	877
	878	END SUBROUTINE compute_fog_K84
	879
[1909]	880	SUBROUTINE compute_fog_RUC(d1, d2, d3, qv, ta, pres, fog, vis)
[1908]	881	! Subroutine to compute fog: qcloud + qice /= 0.
	882	! And visibility following RUC method Smirnova, T. G., S. G. Benjamin, and J. M. Brown, 2000: Case
	883	! study verification of RUC/MAPS fog and visibility forecasts. Preprints, 9 th Conference on
	884	! Aviation, Range, and Aerospace Meteorlogy, AMS, Orlando, FL, Sep. 2000. Paper#2.3, 6 pp.
	885
	886	IMPLICIT NONE
	887
	888	INTEGER, INTENT(in) :: d1, d2, d3
[1909]	889	REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in) :: qv, ta, pres
[1908]	890	INTEGER, DIMENSION(d1,d2,d3), INTENT(out) :: fog
	891	REAL(r_k), DIMENSION(d1,d2,d3), INTENT(out) :: vis
	892
	893	! Local
	894	INTEGER :: i, j, it
	895
	896	!!!!!!! Variables
[1909]	897	! qv: water vapor mixing ratio [kgkg-1]
	898	! ta: temperature [K]
	899	! pres: pressure [Pa]
[1908]	900	! fog: presence of fog (1: yes, 0: no)
	901	! vis: visibility within fog [km]
	902
[1909]	903	fname = 'compute_fog_RUC'
[1908]	904
	905	DO i=1, d1
	906	DO j=1, d2
	907	DO it=1, d3
[1909]	908	CALL var_fog_RUC(qv(i,j,it), ta(i,j,it), pres(i,j,it), fog(i,j,it), vis(i,j,it))
[1908]	909	END DO
	910	END DO
	911	END DO
	912
	913	RETURN
	914
	915	END SUBROUTINE compute_fog_RUC
	916
[1909]	917	SUBROUTINE compute_fog_FRAML50(d1, d2, d3, qv, ta, pres, fog, vis)
	918	! Subroutine to compute fog (vis < 1 km) and visibility following
	919	! Gultepe, I. and J.A. Milbrandt, 2010: Probabilistic Parameterizations of Visibility Using
	920	! Observations of Rain Precipitation Rate, Relative Humidity, and Visibility. J. Appl. Meteor.
	921	! Climatol., 49, 36-46, https://doi.org/10.1175/2009JAMC1927.1
	922	! Interest is focused on a 'general' fog/visibilty approach, thus the fit at 50 % of probability is
	923	! chosen
	924	! Effects from precipitation are not considered
	925
	926	IMPLICIT NONE
	927
	928	INTEGER, INTENT(in) :: d1, d2, d3
	929	REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in) :: qv, ta, pres
	930	INTEGER, DIMENSION(d1,d2,d3), INTENT(out) :: fog
	931	REAL(r_k), DIMENSION(d1,d2,d3), INTENT(out) :: vis
	932
	933	! Local
	934	INTEGER :: i, j, it
	935
	936	!!!!!!! Variables
	937	! qv: mixing ratio in [kgkg-1]
	938	! ta: temperature [K]
	939	! pres: pressure field [Pa]
	940	! fog: presence of fog (1: yes, 0: no)
	941	! vis: visibility within fog [km]
	942
	943	fname = 'compute_fog_FRAML50'
	944
	945	DO i=1, d1
	946	DO j=1, d2
	947	DO it=1, d3
	948	CALL var_fog_FRAML50(qv(i,j,it), ta(i,j,it), pres(i,j,it), fog(i,j,it), vis(i,j,it))
	949	END DO
	950	END DO
	951	END DO
	952
	953	RETURN
	954
	955	END SUBROUTINE compute_fog_FRAML50
	956
[2260]	957	SUBROUTINE compute_range_faces(d1, d2, lon, lat, hgt, xdist, ydist, dist, face, dsfilt, dsnewrange, &
	958	hvalrng, hgtmax, pthgtmax, derivhgt, peaks, valleys, origfaces, filtfaces, ranges, rangeshgtmax, &
[2223]	959	ptrangeshgtmax)
[2208]	960	! Subroutine to compute faces [uphill, valleys, downhill] of a mountain range along a given face
	961
	962	IMPLICIT NONE
	963
[2215]	964	INTEGER, INTENT(in) :: d1, d2
	965	REAL(r_k), INTENT(in) :: dsfilt, dsnewrange, hvalrng
[2260]	966	REAL(r_k), DIMENSION(d1,d2), INTENT(in) :: lon, lat, hgt, xdist, ydist, dist
[2208]	967	CHARACTER(len=*) :: face
[2214]	968	REAL(r_k), DIMENSION(d1,d2), INTENT(out) :: derivhgt, hgtmax, rangeshgtmax
[2213]	969	INTEGER, DIMENSION(d1,d2), INTENT(out) :: pthgtmax, origfaces, filtfaces, peaks, &
[2223]	970	valleys, ranges, ptrangeshgtmax
[2208]	971	! Local
	972	INTEGER :: i, j
[2214]	973	INTEGER :: pthgtmax1, Npeaks, Nvalleys, Nranges
[2213]	974	REAL(r_k) :: hgtmax1
[2214]	975	INTEGER, DIMENSION(d1) :: ipeaks1, ivalleys1, irangeshgtmax1
	976	INTEGER, DIMENSION(d2) :: ipeaks2, ivalleys2, irangeshgtmax2
	977	REAL(r_k), DIMENSION(d1) :: rangeshgtmax1
	978	REAL(r_k), DIMENSION(d2) :: rangeshgtmax2
	979	INTEGER, DIMENSION(2,d1) :: ranges1
	980	INTEGER, DIMENSION(2,d2) :: ranges2
[2330]	981	INTEGER, DIMENSION(d1,d2) :: iranges
[2332]	982	LOGICAL, DIMENSION(d1,d2) :: Lranges
[2208]	983
	984	!!!!!!! Variables
	985	! lon: longitude [degrees east]
	986	! lat: latitude [degrees north]
	987	! hgt: topograpical height [m]
[2212]	988	! face: which face (axis along which produce slices) to use to compute the faces: WE, SN
[2215]	989	! dsfilt: distance to filter orography smaller scale of it [m]
	990	! dsnewrange: distance to start a new mountain range [m]
	991	! hvalrng: maximum height of a valley to mark change of range [m]
[2213]	992	! hgtmax: maximum height of the face [m]
	993	! pthgtmax: grid point of the maximum height [1]
[2212]	994	! derivhgt: topograpic derivative along axis [m deg-1]
	995	! peaks: peak point
	996	! valleys: valley point
	997	! origfaces: original faces (-1, downhill; 0: valley; 1: uphill)
	998	! filtfaces: filtered faces (-1, downhill; 0: valley; 1: uphill)
[2223]	999	! ranges: number of range
[2214]	1000	! rangeshgtmax: maximum height for each individual range [m]
	1001	! ptrangeshgtmax: grid point maximum height for each individual range [1]
[2208]	1002
	1003	fname = 'compute_range_faces'
	1004
[2212]	1005	peaks = 0
	1006	valleys = 0
[2213]	1007	pthgtmax = 0
[2215]	1008	rangeshgtmax = fillVal64
[2208]	1009	IF (TRIM(face) == 'WE') THEN
	1010	DO j=1, d2
[2217]	1011	!PRINT ,'Lluis:', j-1, '**'
[2260]	1012	CALL var_range_faces(d1, lon(:,j), lat(:,j), hgt(:,j), xdist(:,j), dsfilt, &
	1013	dsnewrange, hvalrng, hgtmax1, pthgtmax1, derivhgt(:,j), Npeaks, ipeaks1, &
	1014	Nvalleys, ivalleys1, origfaces(:,j), filtfaces(:,j), Nranges, ranges1, &
	1015	rangeshgtmax1, irangeshgtmax1)
[2213]	1016	hgtmax(:,j) = hgtmax1
	1017	pthgtmax(pthgtmax1,j) = 1
[2212]	1018	DO i=1, Npeaks
	1019	peaks(ipeaks1(i),j) = 1
	1020	END DO
	1021	DO i=1, Nvalleys
	1022	valleys(ivalleys1(i),j) = 1
	1023	END DO
[2214]	1024	DO i=1, Nranges
[2330]	1025	iranges(ranges1(1,i):ranges1(2,i),j) = i
[2214]	1026	rangeshgtmax(ranges1(1,i):ranges1(2,i),j) = rangeshgtmax1(i)
	1027	ptrangeshgtmax(irangeshgtmax1(i),j) = 1
	1028	END DO
[2208]	1029	END DO
	1030	ELSE IF (TRIM(face) == 'SN') THEN
	1031	DO i=1, d1
[2260]	1032	CALL var_range_faces(d2, lon(i,:), lat(i,:), hgt(i,:), ydist(i,:), dsfilt, &
	1033	dsnewrange, hvalrng, hgtmax1, pthgtmax1, derivhgt(i,:), Npeaks, ipeaks2, &
	1034	Nvalleys, ivalleys2, origfaces(i,:), filtfaces(i,:), Nranges, ranges2, &
	1035	rangeshgtmax2, irangeshgtmax2)
[2213]	1036	hgtmax(i,:) = hgtmax1
	1037	pthgtmax(i,pthgtmax1) = 1
[2212]	1038	DO j=1, Npeaks
	1039	peaks(i,ipeaks2(j)) = 1
	1040	END DO
	1041	DO j=1, Nvalleys
	1042	valleys(i,ivalleys2(j)) = 1
	1043	END DO
[2214]	1044	DO j=1, Nranges
[2330]	1045	iranges(i,ranges2(1,j):ranges2(2,j)) = j
[2214]	1046	rangeshgtmax(i,ranges2(1,j):ranges2(2,j)) = rangeshgtmax2(j)
	1047	ptrangeshgtmax(i,irangeshgtmax2(j)) = 1
	1048	END DO
[2208]	1049	END DO
	1050	ELSE
	1051	PRINT *,TRIM(ErrWarnMsg('err'))
	1052	PRINT *,' ' // TRIM(fname) // ": wrong face: '" // TRIM(face) // "' !!"
	1053	PRINT *,' accepted ones: WE, SN'
	1054	STOP
	1055	END IF
	1056
[2330]	1057	! Homogenizing indices of the ranges
[2341]	1058	CALL continguos_homogene_zones(d1, d2, iranges, Nranges, ranges)
[2345]	1059	WHERE (ranges == -1)
	1060	ranges = fillValI
	1061	END WHERE
[2330]	1062
[2208]	1063	RETURN
	1064
	1065	END SUBROUTINE compute_range_faces
	1066
[2274]	1067	SUBROUTINE compute_cellbnds(dx, dy, sdx, sdy, ulon, ulat, vlon, vlat, xbnds, ybnds)
	1068	! Subroutine to compute cellboundaries using wind-staggered lon, lats as intersection of their related
	1069	! parallels and meridians
	1070
	1071	IMPLICIT NONE
	1072
	1073	INTEGER, INTENT(in) :: dx, dy, sdx, sdy
	1074	REAL(r_k), DIMENSION(sdx, dy), INTENT(in) :: ulon, ulat
	1075	REAL(r_k), DIMENSION(dx, sdy), INTENT(in) :: vlon, vlat
	1076	REAL(r_k), DIMENSION(dx, dy, 4), INTENT(out) :: xbnds, ybnds
	1077
	1078	! Local
	1079	INTEGER :: i,j,iv
	1080	INTEGER :: ix,ex,iy,ey
[2285]	1081	REAL(r_k) :: tmpval1, tmpval2
[2274]	1082	CHARACTER(len=2), DIMENSION(4) :: Svertex
	1083	INTEGER, DIMENSION(4,2,2,2) :: indices
	1084	REAL(r_k), DIMENSION(2) :: ptintsct
	1085	REAL(r_k), DIMENSION(2,2) :: merid, paral
	1086	LOGICAL :: intsct
	1087
	1088	!!!!!!! Variables
	1089	! dx, dy: un-staggered dimensions
	1090	! sdx, sdy: staggered dimensions
	1091	! ulon, ulat: x-wind staggered longitudes and latitudes
	1092	! vlon, vlat: y-wind staggered longitudes and latitudes
	1093	! xbnds, ybnds: x and y cell boundaries
	1094
	1095	fname = 'compute_cellbnds'
	1096
	1097	! Indices to use indices[SW/NW/NE/SE, m/p, x/y, i/e]
	1098	Svertex = (/ 'SW', 'NW', 'NE', 'SE' /)
	1099
	1100	! SW
	1101	indices(1,1,1,1) = 0
	1102	indices(1,1,1,2) = 0
	1103	indices(1,1,2,1) = -1
	1104	indices(1,1,2,2) = 0
	1105	indices(1,2,1,1) = -1
	1106	indices(1,2,1,2) = 0
	1107	indices(1,2,2,1) = -1
	1108	indices(1,2,2,2) = -1
	1109	! NW
	1110	indices(2,1,1,1) = 0
	1111	indices(2,1,1,2) = 0
	1112	indices(2,1,2,1) = 0
	1113	indices(2,1,2,2) = 1
	1114	indices(2,2,1,1) = -1
	1115	indices(2,2,1,2) = 0
	1116	indices(2,2,2,1) = 1
	1117	indices(2,2,2,2) = 1
	1118	! NE
	1119	indices(3,1,1,1) = 1
	1120	indices(3,1,1,2) = 1
	1121	indices(3,1,2,1) = 0
	1122	indices(3,1,2,2) = 1
	1123	indices(3,2,1,1) = 0
	1124	indices(3,2,1,2) = 1
	1125	indices(3,2,2,1) = 1
	1126	indices(3,2,2,2) = 1
	1127	! SE
[2285]	1128	indices(4,1,1,1) = 1
[2274]	1129	indices(4,1,1,2) = 1
	1130	indices(4,1,2,1) = -1
	1131	indices(4,1,2,2) = 0
	1132	indices(4,2,1,1) = 0
	1133	indices(4,2,1,2) = 1
	1134	indices(4,2,2,1) = -1
	1135	indices(4,2,2,2) = -1
	1136
[2285]	1137	DO i=2,dx-1
[2274]	1138	DO j=1,dy
	1139	DO iv=1,4
	1140
	1141	ix = MAX(i+indices(iv,1,1,1),1)
	1142	!ex = MIN(i+indices(iv,1,1,2),dx)
[2285]	1143	ex = MAX(i+indices(iv,1,1,2),1)
[2274]	1144	iy = MAX(j+indices(iv,1,2,1),1)
	1145	ey = MIN(j+indices(iv,1,2,2),dy)
	1146
	1147	merid(1,1) = ulon(ix,iy)
	1148	merid(1,2) = ulat(ix,iy)
	1149	merid(2,1) = ulon(ex,ey)
	1150	merid(2,2) = ulat(ex,ey)
	1151
	1152	ix = MAX(i+indices(iv,2,1,1),1)
	1153	ex = MIN(i+indices(iv,2,1,2),dx)
	1154	iy = MAX(j+indices(iv,2,2,1),1)
	1155	!ey = MIN(i+indices(iv,2,2,2),dy)
[2285]	1156	ey = MAX(j+indices(iv,2,2,2),1)
[2274]	1157	paral(1,1) = vlon(ix,iy)
	1158	paral(1,2) = vlat(ix,iy)
	1159	paral(2,1) = vlon(ex,ey)
	1160	paral(2,2) = vlat(ex,ey)
	1161
	1162	CALL intersection_2Dlines(merid, paral, intsct, ptintsct)
	1163	IF (.NOT.intsct) THEN
[2285]	1164	msg = 'not intersection found for ' // Svertex(iv) // ' vertex'
[2274]	1165	CALL ErrMsg(msg, fname, -1)
	1166	END IF
	1167	xbnds(i,j,iv) = ptintsct(1)
	1168	ybnds(i,j,iv) = ptintsct(2)
[2285]	1169
[2274]	1170	END DO
	1171	END DO
	1172	END DO
	1173
[2285]	1174	! Dealing with the boundary values
	1175	i = 1
	1176	DO j=1,dy
	1177	DO iv=1,4
	1178
	1179	ix = MAX(i+indices(iv,1,1,1),1)
	1180	!ex = MIN(i+indices(iv,1,1,2),dx)
	1181	ex = MAX(i+indices(iv,1,1,2),1)
	1182	iy = MAX(j+indices(iv,1,2,1),1)
	1183	ey = MIN(j+indices(iv,1,2,2),dy)
	1184	merid(1,1) = ulon(ix,iy)
	1185	merid(1,2) = ulat(ix,iy)
	1186	merid(2,1) = ulon(ex,ey)
	1187	merid(2,2) = ulat(ex,ey)
	1188
	1189	ix = MAX(i+indices(iv,2,1,1),1)
	1190	ex = MIN(i+indices(iv,2,1,2),dx)
	1191	iy = MAX(j+indices(iv,2,2,1),1)
	1192	!ey = MIN(i+indices(iv,2,2,2),dy)
	1193	ey = MAX(j+indices(iv,2,2,2),1)
	1194	IF (iv == 1 .OR. iv == 2) THEN
	1195	! Projecting values using dx from next grid point
	1196	tmpval1 = vlon(2,iy)
	1197	paral(2,1) = vlon(ex,ey)
	1198	tmpval2 = tmpval1 - paral(2,1)
	1199	paral(1,1) = paral(2,1) - tmpval2
	1200	tmpval1 = vlat(2,iy)
	1201	paral(2,2) = vlat(ex,ey)
	1202	tmpval2 = tmpval1 - paral(2,2)
	1203	paral(1,2) = paral(2,2) - tmpval2
	1204	ELSE
	1205	paral(1,1) = vlon(ix,iy)
	1206	paral(1,2) = vlat(ix,iy)
	1207	paral(2,1) = vlon(ex,ey)
	1208	paral(2,2) = vlat(ex,ey)
	1209	END IF
	1210
	1211	CALL intersection_2Dlines(merid, paral, intsct, ptintsct)
	1212	IF (.NOT.intsct) THEN
	1213	msg = 'not intersection found for ' // Svertex(iv) // ' vertex'
	1214	CALL ErrMsg(msg, fname, -1)
	1215	END IF
	1216	xbnds(i,j,iv) = ptintsct(1)
	1217	ybnds(i,j,iv) = ptintsct(2)
	1218
	1219	END DO
	1220	END DO
	1221
	1222	i = dx
	1223	DO j=1,dy
	1224	DO iv=1,4
	1225
	1226	ix = MAX(i+indices(iv,1,1,1),1)
	1227	!ex = MIN(i+indices(iv,1,1,2),dx)
	1228	ex = MAX(i+indices(iv,1,1,2),1)
	1229	iy = MAX(j+indices(iv,1,2,1),1)
	1230	ey = MIN(j+indices(iv,1,2,2),dy)
	1231	merid(1,1) = ulon(ix,iy)
	1232	merid(1,2) = ulat(ix,iy)
	1233	merid(2,1) = ulon(ex,ey)
	1234	merid(2,2) = ulat(ex,ey)
	1235
	1236	ix = MAX(i+indices(iv,2,1,1),1)
	1237	ex = MIN(i+indices(iv,2,1,2),dx)
	1238	iy = MAX(j+indices(iv,2,2,1),1)
	1239	!ey = MIN(i+indices(iv,2,2,2),dy)
	1240	ey = MAX(j+indices(iv,2,2,2),1)
	1241	IF (iv == 3 .OR. iv == 4) THEN
	1242	! Projecting values using dx from previous grid point
	1243	tmpval1 = vlon(dx-1,iy)
	1244	paral(2,1) = vlon(ex,ey)
	1245	tmpval2 = tmpval1 - paral(2,1)
	1246	paral(1,1) = paral(2,1) - tmpval2
	1247	tmpval1 = vlat(dx-1,iy)
	1248	paral(2,2) = vlat(ex,ey)
	1249	tmpval2 = tmpval1 - paral(2,2)
	1250	paral(1,2) = paral(2,2) - tmpval2
	1251	ELSE
	1252	paral(1,1) = vlon(ix,iy)
	1253	paral(1,2) = vlat(ix,iy)
	1254	paral(2,1) = vlon(ex,ey)
	1255	paral(2,2) = vlat(ex,ey)
	1256	END IF
	1257
	1258	CALL intersection_2Dlines(merid, paral, intsct, ptintsct)
	1259	IF (.NOT.intsct) THEN
	1260	msg = 'not intersection found for ' // Svertex(iv) // ' vertex'
	1261	CALL ErrMsg(msg, fname, -1)
	1262	END IF
	1263	xbnds(i,j,iv) = ptintsct(1)
	1264	ybnds(i,j,iv) = ptintsct(2)
	1265
	1266	END DO
	1267	END DO
	1268
[2277]	1269	END SUBROUTINE compute_cellbnds
	1270
	1271	SUBROUTINE compute_cellbndsreg(dx, dy, lon, lat, xbnds, ybnds)
	1272	! Subroutine to compute cellboundaries using lon, lat from a reglar lon/lat projection as intersection
	1273	! of their related parallels and meridians
	1274
	1275	IMPLICIT NONE
	1276
	1277	INTEGER, INTENT(in) :: dx, dy
	1278	REAL(r_k), DIMENSION(dx, dy), INTENT(in) :: lon, lat
	1279	REAL(r_k), DIMENSION(dx, dy, 4), INTENT(out) :: xbnds, ybnds
	1280
	1281	! Local
	1282	INTEGER :: i,j,iv
	1283	INTEGER :: ix,ex,iy,ey
	1284	CHARACTER(len=2), DIMENSION(4) :: Svertex
[2278]	1285	INTEGER, DIMENSION(4,2,2) :: indices
[2277]	1286
	1287	!!!!!!! Variables
	1288	! dx, dy: un-staggered dimensions
	1289	! lon, lat: longitudes and latitudes
	1290	! xbnds, ybnds: x and y cell boundaries
	1291
	1292	fname = 'compute_cellbndsreg'
	1293
	1294	! Indices to use indices[SW/NW/NE/SE, m/p, x/y, i/e]
	1295	Svertex = (/ 'SW', 'NW', 'NE', 'SE' /)
	1296
	1297	! SW
[2278]	1298	indices(1,1,1) = -1
	1299	indices(1,1,2) = 0
	1300	indices(1,2,1) = -1
	1301	indices(1,2,2) = 0
[2277]	1302	! NW
[2278]	1303	indices(2,1,1) = -1
	1304	indices(2,1,2) = 0
	1305	indices(2,2,1) = 0
	1306	indices(2,2,2) = 1
[2277]	1307	! NE
[2278]	1308	indices(3,1,1) = 0
	1309	indices(3,1,2) = 1
	1310	indices(3,2,1) = 0
	1311	indices(3,2,2) = 1
[2277]	1312	! SE
[2278]	1313	indices(4,1,1) = 0
	1314	indices(4,1,2) = 1
	1315	indices(4,2,1) = -1
	1316	indices(4,2,2) = 0
[2277]	1317
	1318	DO i=1,dx
	1319	DO j=1,dy
	1320	DO iv=1,4
	1321
[2278]	1322	ix = MAX(i+indices(iv,1,1),1)
	1323	ix = MIN(ix,dx)
	1324	ex = MAX(i+indices(iv,1,2),1)
	1325	ex = MIN(ex,dx)
	1326	iy = MAX(j+indices(iv,2,1),1)
	1327	iy = MIN(iy,dy)
	1328	ey = MAX(j+indices(iv,2,2),1)
	1329	ey = MIN(ey,dy)
[2277]	1330
[2278]	1331	xbnds(i,j,iv) = 0.5*(lon(ix,iy) + lon(ex,ey))
	1332	ybnds(i,j,iv) = 0.5*(lat(ix,iy) + lat(ex,ey))
[2277]	1333
	1334	END DO
	1335	END DO
	1336	END DO
	1337
[2274]	1338	END SUBROUTINE
	1339
[2260]	1340	SUBROUTINE compute_Koeppen_Geiger_climates(dx, dy, dt, pr, tas, climatesI, climatesS, climlegend)
	1341	! Subroutine to compute the Koeppen-Geiger climates after:
	1342	! Kottek et al., Meteorologische Zeitschrift, Vol. 15, No. 3, 259-263
	1343
	1344	IMPLICIT NONE
	1345
	1346	INTEGER, INTENT(in) :: dx, dy, dt
	1347	REAL(r_k), DIMENSION(dx, dy, dt), INTENT(in) :: pr, tas
	1348	INTEGER, DIMENSION(dy, dt), INTENT(out) :: climatesI
	1349	CHARACTER(LEN=3), DIMENSION(dy, dt), INTENT(out) :: climatesS
	1350	CHARACTER(LEN=1000), INTENT(out) :: climlegend
	1351
	1352	! Local
	1353	INTEGER :: i,j
	1354
	1355
	1356
	1357
	1358	END SUBROUTINE compute_Koeppen_Geiger_climates
	1359
[2387]	1360	SUBROUTINE compute_tws_RK1(d1, tas, hurs, tws)
	1361	! Subroutine to compute Wet Bulb temperature of 1D series of values using equation after:
	1362	! Stull, R. (2011), J. Appl. Meteor. Climatol. 50(11):2267-2269. doi: 10.1175/JAMC-D-11-0143.1
	1363
	1364	IMPLICIT NONE
	1365
	1366	INTEGER, INTENT(in) :: d1
	1367	REAL(r_k), DIMENSION(d1), INTENT(in) :: tas, hurs
	1368	REAL(r_k), DIMENSION(d1), INTENT(out) :: tws
	1369
	1370	! Local
	1371	INTEGER :: it
	1372
	1373	!!!!!!! Variables
	1374	! tas: 2-m air temperature [K]
	1375	! hurs: 2-m relative humidity [1]
	1376
	1377	fname = 'compute_tws_RK1'
	1378
	1379	DO it=1, d1
	1380	tws(it) = var_tws_S11(tas(it), hurs(it))
	1381	END DO
	1382
	1383	RETURN
	1384
	1385	END SUBROUTINE compute_tws_RK1
	1386
	1387	SUBROUTINE compute_tws_RK2(d1, d2, tas, hurs, tws)
	1388	! Subroutine to compute Wet Bulb temperature of 2D series of values using equation after:
	1389	! Stull, R. (2011), J. Appl. Meteor. Climatol. 50(11):2267-2269. doi: 10.1175/JAMC-D-11-0143.1
	1390
	1391	IMPLICIT NONE
	1392
	1393	INTEGER, INTENT(in) :: d1, d2
	1394	REAL(r_k), DIMENSION(d1,d2), INTENT(in) :: tas, hurs
	1395	REAL(r_k), DIMENSION(d1,d2), INTENT(out) :: tws
	1396
	1397	! Local
	1398	INTEGER :: i, j
	1399
	1400	!!!!!!! Variables
	1401	! tas: 2-m air temperature [K]
	1402	! hurs: 2-m relative humidity [1]
	1403
	1404	fname = 'compute_tws_RK2'
	1405
	1406	DO i=1, d1
	1407	DO j=1, d2
	1408	tws(i,j) = var_tws_S11(tas(i,j), hurs(i,j))
	1409	END DO
	1410	END DO
	1411
	1412	RETURN
	1413
	1414	END SUBROUTINE compute_tws_RK2
	1415
	1416	SUBROUTINE compute_tws_RK3(d1, d2, d3, tas, hurs, tws)
	1417	! Subroutine to compute Wet Bulb temperature of 3D series of values using equation after:
	1418	! Stull, R. (2011), J. Appl. Meteor. Climatol. 50(11):2267-2269. doi: 10.1175/JAMC-D-11-0143.1
	1419
	1420	IMPLICIT NONE
	1421
	1422	INTEGER, INTENT(in) :: d1, d2, d3
	1423	REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in) :: tas, hurs
	1424	REAL(r_k), DIMENSION(d1,d2,d3), INTENT(out) :: tws
	1425
	1426	! Local
	1427	INTEGER :: i, j, k
	1428
	1429	!!!!!!! Variables
	1430	! tas: 2-m air temperature [K]
	1431	! hurs: 2-m relative humidity [1]
	1432
	1433	fname = 'compute_tws_RK3'
	1434
	1435	DO i=1, d1
	1436	DO j=1, d2
	1437	DO k=1, d3
	1438	tws(i,j,k) = var_tws_S11(tas(i,j,k), hurs(i,j,k))
	1439	END DO
	1440	END DO
	1441	END DO
	1442
	1443	RETURN
	1444
	1445	END SUBROUTINE compute_tws_RK3
	1446
	1447	SUBROUTINE compute_tws_RK4(d1, d2, d3, d4, tas, hurs, tws)
	1448	! Subroutine to compute Wet Bulb temperature of 4D series of values using equation after:
	1449	! Stull, R. (2011), J. Appl. Meteor. Climatol. 50(11):2267-2269. doi: 10.1175/JAMC-D-11-0143.1
	1450
	1451	IMPLICIT NONE
	1452
	1453	INTEGER, INTENT(in) :: d1, d2, d3, d4
	1454	REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in) :: tas, hurs
	1455	REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(out) :: tws
	1456
	1457	! Local
	1458	INTEGER :: i,j,k,l
	1459
	1460	!!!!!!! Variables
	1461	! tas: 2-m air temperature [K]
	1462	! hurs: 2-m relative humidity [1]
	1463
	1464	fname = 'compute_tws_RK4'
	1465
	1466	DO i=1, d1
	1467	DO j=1, d2
	1468	DO k=1, d3
	1469	DO l=1, d4
	1470	tws(i,j,k,l) = var_tws_S11(tas(i,j,k,l), hurs(i,j,k,l))
	1471	END DO
	1472	END DO
	1473	END DO
	1474	END DO
	1475
	1476	RETURN
	1477
	1478	END SUBROUTINE compute_tws_RK4
	1479
[770]	1480	END MODULE module_ForDiagnostics

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

Download in other formats:

Original Format