Context Navigation

module_ForDiagnostics.f90 @ 2734

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

Adding:

'frontogenesis': Calculation of the frontogenesis diagnostics

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

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

Download in other formats:

Original Format