Context Navigation

module_ForDiagnostics.f90 @ 2783

Last change on this file since 2783 was 2765, checked in by lfita, 6 years ago

Adding:

`1st order 2D horizontal gradient'

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

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

Download in other formats:

Original Format