Context Navigation

module_ForDiagnostics.f90 @ 2299

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

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

Download in other formats:

Original Format