Context Navigation

source: branches/publications/ORCHIDEE_IPSLCM5A2.1.r5307/ORCHIDEE/src_global/interpweight.f90 @ 7683

Last change on this file since 7683 was 3858, checked in by josefine.ghattas, 8 years ago
Bug corrections for interpolation of lai, the bug was introduced in the trunk rev [3831]. Following corrections correspond to changeset [3853] and [3856] on the branche ORCHIDEE-DYNAMICO, done by L Fita.
File size: 231.5 KB

Line
1	MODULE interpweight
2	!=================================================================================================================================
3	! MODULE : interpweight
4	!
5	! CONTACT : orchidee-help _at_ ipsl.jussieu.fr
6	!
7	! LICENCE : IPSL (2006)
8	! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
9	!
10	!>\BRIEF ! Specific module for the interpolation in ORCHIDEE
11	!!
12	!!\n DESCRIPTION: L. Fita, October 2016
13	!! This module is used in the re-organization done in the section of the code
14	!! responsible of the interpolation of the morphology datat (PFT, albedo, LAI,
15	!! ...) from different files done by L. Fita in the mark of the HEAT project.
16	!! Previous version of the interpolation of these files, proceed by three
17	!! basic steps common for all the files:
18	!! 1.: Modification of values from file
19	!! 2.: Creation of land/sea mask using values from file
20	!! 3.: Interpolate values
21	!! This module servs to generalize all these steps and construct a common
22	!! interface independently of the file. Each step is generalize and the
23	!! different processes almost specific for each file in the previous form
24	!! are introduced as different methodologies. Thus now it provides:
25	!! 1. Modification of initial values: Removing that wrong values from the file
26	!\| `interpweight_modifying_input[1/2/3/4]D': subroutines to modify 1D, 2D, 3D
27	!! or 4D input data from file, using different methods
28	!! 2. Masking input: Compute âon flyâ the values of the land/sea mask taking
29	!! the values in the file
30	!! `interpweight_masking_input[1/2/3/4]D': subroutines to compute the mask
31	!! using data from input file using different methods: nomask, mbelow,
32	!! mabove, msumrange, var
33	!! 3. Area weights: Get coincident areas from input file to the target
34	!! projection using aggregate_p
35	!! No changes on it
36	!! 4. Interpolate: Use obtained areas and interpolate values at each grid
37	!! point at the same format as it will be provided by XIOS
38	!! `interpweight_provide_fractions[1/2/3/4]D': Perform interpolation for
39	!! fraction/category data
40	!! `interpweight_provide_interpolation[2/4]D': Perform interpolation for
41	!! continuous data
42	!! `variableusetypes': Variable to provide the values along the additional
43	!! dimension to perform the interpolation. (e.g.: in case the 13 pfts
44	!! are: 1,3,6,18,23,34,35,39,48,...)
45	!! 5. A new variable is added which explains the âavailabilityâ of data to
46	!! perform the interpolation (0, 1). When it is negative it means that there
47	!! was no data for that grid point
48	!! availability = SUM(area_source)/Area_target_grid_point
49	!! This module will disappear once interpolation will be done throuhguot XIOS
50	!!
51	!! More details at:
52	!! https://forge.ipsl.jussieu.fr/orchidee/wiki/DevelopmentActivities/ORCHIDEE-DYNAMICO
53	!!
54	!! RECENT CHANGE(S): None
55	!!
56	!! REFERENCE(S): None
57	!!
58	!! SVN :
59	!! $HeadURL: $
60	!! $Date: $
61	!! $Revision: $
62	!! \n
63	!_================================================================================================================================
64
65	USE ioipsl_para
66	USE interpol_help
67	USE grid
68
69	IMPLICIT NONE
70
71	PRIVATE
72	PUBLIC interpweight_1D, interpweight_2D, interpweight_3D, interpweight_4D, interpweight_2Dcont, &
73	interpweight_4Dcont, interpweight_get_varNdims_file, interpweight_RangeR, &
74	interpweight_get_var2dims_file, interpweight_get_var3dims_file, interpweight_get_var4dims_file, &
75	interpweight_Index2DLonLat, interpweight_ValVecR
76
77	CONTAINS
78
79	!! ================================================================================================================================
80	!! SUBROUTINE : interpweight_1D
81	!!
82	!>\BRIEF ! Interpolate any input file to the grid of the model for a 1D-variable
83	!!
84	!! DESCRIPTION : (definitions, functional, design, flags):
85	!!
86	!! RECENT CHANGE(S): None
87	!!
88	!! MAIN OUTPUT VARIABLE(S): outvar1D, aoutvar
89	!!
90	!! REFERENCE(S) : None
91	!!
92	!! FLOWCHART : None
93	!! \n
94	!_ ================================================================================================================================
95	SUBROUTINE interpweight_1D(nbpt, Nvariabletypes, variabletypes, lalo, resolution, neighbours, &
96	contfrac, filename, varname, inlonname, inlatname, varmin, varmax, noneg, masktype, &
97	maskvalues, maskvarname, dim1, dim2, initime, typefrac, &
98	maxresollon, maxresollat, outvar1D, aoutvar)
99
100	USE ioipsl
101	USE constantes_var
102	USE mod_orchidee_para
103
104	IMPLICIT NONE
105
106	!! 0. Variables and parameter declaration
107	!
108	!! 0.1 Input variables
109	!
110	INTEGER(i_std), INTENT(in) :: nbpt !! Number of points for which the data
111	!! needs to be interpolated
112	INTEGER(i_std), INTENT(in) :: Nvariabletypes !! Number of types of the variable
113	REAL(r_std), DIMENSION(Nvariabletypes), & !! Vector of values of the types
114	INTENT(in) :: variabletypes !! (-1, not used)
115	REAL(r_std), DIMENSION(nbpt,2), INTENT(in) :: lalo !! Vector of latitude and longitudes
116	!! (beware of the order = 1 : latitude,
117	!! 2 : longitude)
118	REAL(r_std), DIMENSION(nbpt,2), INTENT(in) :: resolution !! The size in km of each grid-box in X and Y
119	!
120	INTEGER(i_std), DIMENSION(nbpt,NbNeighb), INTENT(in) :: neighbours !! Vector of neighbours for each grid point
121	!! 1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW)
122
123	REAL(r_std), DIMENSION(nbpt), INTENT(in) :: contfrac !! Fraction of land in each grid box.
124	CHARACTER(LEN=80), INTENT(in) :: filename !! name of the input map to read
125
126	CHARACTER(LEN=80), INTENT(in) :: varname !! name of the variable to interpolate
127	CHARACTER(LEN=80), INTENT(in) :: inlonname, inlatname !! names of the longitude and latitude
128	REAL(r_std), INTENT(inout) :: varmin, varmax !! min/max values to use for the
129	!! renormalization
130	! Transformations to apply to the original read values from the file
131	LOGICAL, INTENT(in) :: noneg !! no negative values
132	CHARACTER(LEN=50), INTENT(in) :: masktype !! Type of masking
133	!! 'nomask': no-mask is applied
134	!! 'mbelow': take values below maskvalues(1)
135	!! 'mabove': take values above maskvalues(1)
136	!! 'msumrange': take values within 2 ranges;
137	!! maskvalues(2) <= SUM(vals(k)) <= maskvalues(1)
138	!! maskvalues(1) < SUM(vals(k)) <= maskvalues(3)
139	!! (normalized by maskvalues(3))
140	!! 'var': mask values are taken from a
141	!! variable (>0)
142	REAL(r_std), DIMENSION(3), INTENT(in) :: maskvalues !! values to use to mask (according to
143	!! `masktype')
144	CHARACTER(LEN=250), INTENT(in) :: maskvarname !! name of the variable to use to mask
145	INTEGER(i_std), INTENT(in) :: dim1, dim2 !! 4D complementary dimension sizes
146	CHARACTER(LEN=50), INTENT(in) :: typefrac !! Type of fraction retrieval:
147	!! 'default': standard
148
149	INTEGER(i_std), INTENT(in) :: initime !! Initial time to take (for input data
150	!! with time values, -1 not used)
151	REAL(r_std), INTENT(in) :: maxresollon,maxresollat !! lon,lat maximum resolutions (in m)
152	!! (-un, not used)
153	!! 0.2 Modified variables
154	!
155	!! 0.3 output variables
156	!
157	REAL(r_std), DIMENSION(nbpt), INTENT(out) :: outvar1D !! 1D output variable and re-dimensioned
158	REAL(r_std), DIMENSION(nbpt), INTENT(out) :: aoutvar !! availability of input data to
159	!! interpolate output variable
160	!! (on the nbpt space)
161	!
162	!! 0.4 Local variables
163	!
164	INTEGER(i_std) :: iml, jml, lml, tml, fid, ib, ip, jp, it, jj, jv
165	INTEGER :: inNdims !! Number of dimensions of the input data
166	INTEGER :: inNLldims !! Number of dimensions of the lon/lat input data
167	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lat_in, lon_in !! latitude and
168	!! longitude, extract from input map
169	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: sub_area !! the area of the fine grid in the
170	!! model grid ???
171	!! cf src_global/interpol_help.f90,
172	!! line 377, called "areaoverlap"
173	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: sub_index !! the indexes from the grid boxes from the
174	!! data that go into the
175	!! model's boxes
176	!! cf src_global/interpol_help.f90,
177	!! line 300, called "ip"
178
179	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: invar1D !! 1D input values
180	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: invar2D !! 2D input values
181	REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:) :: invar3D !! 3D input values
182	REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:,:) :: invar4D !! 4D input values
183
184	REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:) :: resol_in
185	INTEGER(i_std), ALLOCATABLE, DIMENSION(:) :: mask
186	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: maskvar
187	INTEGER(i_std) :: nix, njx
188	INTEGER(i_std) :: idi, nbvmax !! nbvmax : number of maximum vegetation map
189	!! points in the GCM grid ; idi : its counter
190	CHARACTER(LEN=30) :: callsign !! Allows to specify which variable is beeing
191	!! treated
192	LOGICAL :: ok_interpol !! optionnal return of aggregate_2d
193	INTEGER :: ALLOC_ERR
194	CHARACTER(LEN=50) :: S1
195	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: variableusetypes
196	CHARACTER(LEN=256) :: msg
197	INTEGER, DIMENSION(:), ALLOCATABLE :: indims
198
199	! Debug vars
200	INTEGER :: nb_coord, nb_var, nb_gat
201
202
203	!_ ================================================================================================================================
204	! Allocating input data
205
206	IF (printlev >= 3) WRITE(numout,*)"In interpweight_1D filename: ",TRIM(filename), &
207	" varname:'" // TRIM(varname) // "'"
208
209	IF (is_root_prc) THEN
210	IF (printlev >=5) THEN
211	WRITE(numout,*) "Entering 'interpweight_1D'. Debug mode."
212	WRITE(numout,*)'(/," --> fliodmpf")'
213	CALL fliodmpf (TRIM(filename))
214	WRITE(numout,*)'(/," --> flioopfd")'
215	ENDIF
216	CALL flioopfd (TRIM(filename),fid,nb_dim=nb_coord,nb_var=nb_var,nb_gat=nb_gat)
217	IF (printlev >=5) THEN
218	WRITE (numout,'(" Number of coordinate in the file : ",I2)') nb_coord
219	WRITE (numout,'(" Number of variables in the file : ",I2)') nb_var
220	WRITE (numout,'(" Number of global attributes in the file : ",I2)') nb_gat
221	ENDIF
222	ENDIF
223
224	! Getting shape of input variable from input file
225	inNdims = interpweight_get_varNdims_file(filename, varname)
226	IF (is_root_prc) CALL flininfo(filename, iml, jml, lml, tml, fid)
227
228	InputDims: SELECT CASE (inNdims)
229	CASE (1)
230	CALL bcast(iml)
231
232	ALLOCATE(invar1D(iml), STAT=ALLOC_ERR)
233	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_1D',"Problem in allocation of variable '" // &
234	TRIM(varname) ,'','')
235	IF (is_root_prc) CALL flinget(fid, TRIM(varname), iml, 0, 0, 0, 1, 1, invar1D)
236	IF (printlev >= 5) THEN
237	WRITE(numout,*)' interpweight_1D input data dimensions _______'
238	WRITE(numout,*)' iml:',iml
239	END IF
240
241	CASE DEFAULT
242	WRITE (S1,'(I3)') inNdims
243	CALL ipslerr_p(3,'interpweight_1D',"Input number of dimensions: '" // TRIM(S1) // "' not ready !!", &
244	'','')
245
246	END SELECT InputDims
247	CALL bcast(invar1D)
248
249	IF (printlev >= 5) WRITE(numout,*)' interpweight_1D getting variable ended!'
250
251	! Getting longitudes/latitudes from input file
252	inNLldims = interpweight_get_varNdims_file(filename, inlonname)
253
254	IF (inNLldims == 1) THEN
255	! Allocate memory for latitudes
256	ALLOCATE(lon_in(iml), STAT=ALLOC_ERR)
257	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_1D','Problem in allocation of variable lon_in','','')
258
259	ALLOCATE(lat_in(iml), STAT=ALLOC_ERR)
260	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_1D','Problem in allocation of variable lat_in','','')
261
262	IF (is_root_prc) THEN
263	CALL flinget(fid, TRIM(inlonname), iml, 0, 0, 0, 1, 1, lon_in)
264	CALL flinget(fid, TRIM(inlatname), iml, 0, 0, 0, 1, 1, lat_in)
265	ENDIF
266
267	END IF
268	CALL bcast(lon_in)
269	CALL bcast(lat_in)
270
271	IF (noneg) THEN
272	IF (is_root_prc) THEN
273	SELECT CASE (inNdims)
274	CASE(1)
275	CALL interpweight_modifying_input1D(iml, 0, 0, 0, 1, 1, zero, invar1D)
276	END SELECT
277	END IF
278	END IF
279	IF (printlev >= 5) WRITE(numout,*)' interpweight_1D modifying input ended!'
280
281	IF (is_root_prc) CALL flinclo(fid)
282
283	!
284	! Consider all points a priori
285	!
286	ALLOCATE(mask(iml), STAT=ALLOC_ERR)
287	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_1D','Problem in allocation of variable mask','','')
288	mask = zero
289
290	IF (TRIM(masktype) == 'var') THEN
291	ALLOCATE(maskvar(iml), STAT=ALLOC_ERR)
292	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_1D','Problem in allocation of variable maskvar','','')
293	maskvar = zero
294	! Reads the mask from the file
295	IF (is_root_prc) THEN
296	! ATTTENTION!! One can not use flinget in this way. It always spect to have a REAL variable!
297	CALL flinget(fid, TRIM(maskvarname), iml, 0, 0, 0, 1, 1, maskvar)
298	WHERE (maskvar > zero)
299	mask = un
300	END WHERE
301	END IF
302	CALL bcast(mask)
303	ELSE
304	SELECT CASE (inNdims)
305	CASE(1)
306	CALL interpweight_masking_input1D(iml, 0, 0, 0, 0, 0, masktype, maskvalues, invar1D, &
307	mask)
308	END SELECT
309	END IF
310	IF (printlev >= 5) WRITE(numout,*)' interpweight_1D masking input ended!'
311
312	! Computation of nbvmax
313	! In 'condveg_soilalb, slowproc_update' nbvmax=200
314	!
315	! The number of maximum vegetation map points in the GCM grid is estimated.
316	! Some lmargin is taken.
317	!
318	IF (is_root_prc) THEN
319	IF ( (maxresollon /= -un) .AND. (maxresollat /= -un) ) THEN
320	! Assuming regular lon/lat projection for the input data but it is not caluclated assuming a maximum &
321	! lon,lat resolution given as input paramater
322	nix=INT(MAXVAL(resolution(:,1))/maxresollon)+2
323	njx=INT(MAXVAL(resolution(:,2))/maxresollat)+2
324	ELSE
325	! Assuming regular lon/lat projection for the input data
326	CALL interpweight_calc_resolution_in(lon_in, lat_in, iml, jml, mincos, R_Earth, pi, resol_in)
327	nix=INT(MAXVAL(resolution(:,1))/MAXVAL(resol_in(:,:,1)))+2
328	njx=INT(MAXVAL(resolution(:,2))/MAXVAL(resol_in(:,:,2)))+2
329	END IF
330	nbvmax = nix*njx
331	IF (nbvmax < 200) nbvmax = 200
332	END IF
333	CALL bcast(nbvmax)
334	CALL bcast(nix)
335	CALL bcast(njx)
336
337	!
338	callsign = TRIM(varname) // ' map'
339	!
340	ok_interpol = .FALSE.
341	DO WHILE ( .NOT. ok_interpol )
342	IF (printlev>=3) WRITE(numout,*) "Projection arrays for ",callsign," : "
343	IF (printlev>=3) WRITE(numout,*) "nbvmax = ",nbvmax, nix, njx
344
345
346	ALLOCATE(sub_index(nbpt, nbvmax), STAT=ALLOC_ERR)
347	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_1D','Problem in allocation of variable sub_index','','')
348
349	sub_index(:,:)=zero
350
351	ALLOCATE(sub_area(nbpt, nbvmax), STAT=ALLOC_ERR)
352	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_1D','Problem in allocation of variable sub_area','','')
353
354	sub_area(:,:)=zero
355	IF (printlev>=3) WRITE(numout,*) 'Input data range LON:', MINVAL(lon_in), MAXVAL(lon_in)
356	IF (printlev>=3) WRITE(numout,*) 'Input data range LAT:', MINVAL(lat_in), MAXVAL(lat_in)
357
358	!
359	CALL aggregate_p(nbpt, lalo, neighbours, resolution, contfrac, &
360	& iml, lon_in, lat_in, maxresollon, maxresollat, callsign, &
361	& nbvmax, sub_index, sub_area, ok_interpol)
362
363	!
364	IF ( .NOT. ok_interpol ) THEN
365	DEALLOCATE(sub_area)
366	DEALLOCATE(sub_index)
367	nbvmax = nbvmax * 2
368	ENDIF
369	ENDDO
370	IF (printlev>=4) WRITE (numout,*) ' interpweight_1D: read/allocate OK'
371
372	! Getting variables thresholds
373	ALLOCATE(variableusetypes(Nvariabletypes), STAT=ALLOC_ERR)
374	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_1D','Problem in allocation of variableusetypes','','')
375
376	IF (variabletypes(1) == -un) THEN
377	variableusetypes = (/ (ip*1., ip=1,Nvariabletypes) /)
378	IF (printlev>=4) WRITE(numout,*) 'interpweight_1D: Note: Using default ',&
379	'equivalence between dimension in the file and dimension in the variable, for file ', filename
380	varmin = 1._r_std
381	varmax = Nvariabletypes*1._r_std
382	ELSE
383	variableusetypes = variabletypes
384	END IF
385
386	outvar1D = zero
387
388	! Providing the fractions of each type for each interpolated grid point
389	CALL interpweight_provide_fractions1D(typefrac, nbpt, Nvariabletypes, variableusetypes, &
390	iml, jml, lml, tml, nbvmax, zero, invar1D, sub_area, sub_index, varmax, varmin, &
391	deux, lalo, outvar1D, aoutvar)
392	IF (printlev>=5) WRITE (numout,*)' interpweight_1D end of interp_provide_fractions1D'
393
394	IF (printlev>=3) THEN
395	IF (COUNT(aoutvar == -1) >= AINT(nbpt*0.4)) THEN
396	WRITE(numout,*) ' interpweight_1D total number of points: ', nbpt
397	WRITE(numout,*) ' interpweight_1D interpolated: ', COUNT(aoutvar /= -1), &
398	' non-interpolated: ', COUNT(aoutvar == -1)
399	msg = "Less than 60% of points '" // TRIM(filename) // "' variable '" // TRIM(varname) // &
400	"' have been succesfully interpolated !!"
401	WRITE(numout,*) ' interpweight_1D: ' // TRIM(msg)
402	ELSE
403	msg = TRIM(filename) // "' variable: '" // TRIM(varname) // "' interpolation Done _______"
404	WRITE(numout,*) " interpweight_1D: '" // TRIM(msg)
405	WRITE(numout,*) ' interpweight_1D total number of points: ', nbpt
406	WRITE(numout,*) ' interpweight_1D interpolated: ', COUNT(aoutvar /= -1), &
407	' non-interpolated: ', COUNT(aoutvar == -1)
408	END IF
409	END IF
410
411	! Scaling to 1 the quality variable
412	DO ip=1, nbpt
413	IF (aoutvar(ip) /= -1) aoutvar(ip) = aoutvar(ip)/(resolution(ip,1)*resolution(ip,2))/contfrac(ip)
414	END DO
415
416	IF (ALLOCATED(lon_in)) DEALLOCATE(lon_in)
417	IF (ALLOCATED(lat_in)) DEALLOCATE(lat_in)
418	IF (ALLOCATED(invar1D)) DEALLOCATE(invar1D)
419	DEALLOCATE(mask)
420	IF (ALLOCATED(maskvar)) DEALLOCATE(maskvar)
421	DEALLOCATE(sub_area)
422	DEALLOCATE(sub_index)
423
424	RETURN
425
426	END SUBROUTINE interpweight_1D
427
428	!! ================================================================================================================================
429	!! SUBROUTINE : interpweight_2D
430	!!
431	!>\BRIEF ! Interpolate any input file to the grid of the model for a 2D-variable
432	!!
433	!! DESCRIPTION : (definitions, functional, design, flags):
434	!!
435	!! RECENT CHANGE(S): None
436	!!
437	!! MAIN OUTPUT VARIABLE(S): outvar2D, aoutvar
438	!!
439	!! REFERENCE(S) : None
440	!!
441	!! FLOWCHART : None
442	!! \n
443	!_ ================================================================================================================================
444	SUBROUTINE interpweight_2D(nbpt, Nvariabletypes, variabletypes, lalo, resolution, neighbours, &
445	contfrac, filename, varname, inlonname, inlatname, varmin, varmax, noneg, masktype, &
446	maskvalues, maskvarname, dim1, dim2, initime, typefrac, &
447	maxresollon, maxresollat, outvar2D, aoutvar)
448
449	USE ioipsl
450	USE constantes_var
451	USE mod_orchidee_para
452
453	IMPLICIT NONE
454
455	!! 0. Variables and parameter declaration
456	!
457	!! 0.1 Input variables
458	!
459	INTEGER(i_std), INTENT(in) :: nbpt !! Number of points for which the data
460	!! needs to be interpolated
461	INTEGER(i_std), INTENT(in) :: Nvariabletypes !! Number of types of the variable
462	REAL(r_std), DIMENSION(Nvariabletypes), & !! Vector of values of the types
463	INTENT(in) :: variabletypes !! (-1, not used)
464	REAL(r_std), DIMENSION(nbpt,2), INTENT(in) :: lalo !! Vector of latitude and longitudes
465	!! (beware of the order = 1 : latitude,
466	!! 2 : longitude)
467	REAL(r_std), DIMENSION(nbpt,2), INTENT(in) :: resolution !! The size in km of each grid-box in X and Y
468	!
469	INTEGER(i_std), DIMENSION(nbpt,NbNeighb), INTENT(in) :: neighbours !! Vector of neighbours for each grid point
470	!! 1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW)
471
472	REAL(r_std), DIMENSION(nbpt), INTENT(in) :: contfrac !! Fraction of land in each grid box.
473	CHARACTER(LEN=80), INTENT(in) :: filename !! name of the input map to read
474
475	CHARACTER(LEN=80), INTENT(in) :: varname !! name of the variable to interpolate
476	CHARACTER(LEN=80), INTENT(in) :: inlonname, inlatname !! names of the longitude and latitude
477	REAL(r_std), INTENT(inout) :: varmin, varmax !! min/max values to use for the
478	!! renormalization
479	! Transformations to apply to the original read values from the file
480	LOGICAL, INTENT(in) :: noneg !! no negative values
481	CHARACTER(LEN=50), INTENT(in) :: masktype !! Type of masking
482	!! 'nomask': no-mask is applied
483	!! 'mbelow': take values below maskvalues(1)
484	!! 'mabove': take values above maskvalues(1)
485	!! 'msumrange': take values within 2 ranges;
486	!! maskvalues(2) <= SUM(vals(k)) <= maskvalues(1)
487	!! maskvalues(1) < SUM(vals(k)) <= maskvalues(3)
488	!! (normalized by maskvalues(3))
489	!! 'var': mask values are taken from a
490	!! variable (>0)
491	REAL(r_std), DIMENSION(3), INTENT(in) :: maskvalues !! values to use to mask (according to
492	!! `masktype')
493	CHARACTER(LEN=250), INTENT(in) :: maskvarname !! name of the variable to use to mask
494	CHARACTER(LEN=50), INTENT(in) :: typefrac !! Type of fraction retrieval:
495	!! 'default': standard
496	INTEGER(i_std), INTENT(in) :: dim1, dim2 !! 4D complementary dimension sizes
497	INTEGER(i_std), INTENT(in) :: initime !! Initial time to take (for input data
498	!! with time values, -1 not used)
499	REAL(r_std), INTENT(in) :: maxresollon,maxresollat !! lon,lat maximum resolutions (in m)
500	!! 0.2 Modified variables
501	!
502	! 0.3 OUTPUT
503	!
504	! REAL(r_std), INTENT(out) :: laimap(nbpt,nvm,12) !! lai read variable and re-dimensioned
505	REAL(r_std), DIMENSION(nbpt), INTENT(out) :: outvar2D !! 2D output variable and re-dimensioned
506	REAL(r_std), DIMENSION(nbpt), INTENT(out) :: aoutvar !! availability of input data to
507	!! interpolate output variable
508	!! (on the nbpt space)
509	!
510	! 0.4 LOCAL
511	!
512	INTEGER(i_std) :: iml, jml, lml, tml, fid, ib, ip, jp, it, jj, jv
513	INTEGER :: inNdims !! Number of dimensions of the input data
514	INTEGER :: inNLldims !! Number of dimensions of the lon/lat input data
515	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lat_in, lon_in !! latitude and
516	!! longitude, extract from input map
517	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lat, lon !! en 2D ???
518	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: sub_area !! the area of the fine grid in the
519	!! model grid ???
520	!! cf src_global/interpol_help.f90,
521	!! line 377, called "areaoverlap"
522	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:,:) :: sub_index !! the indexes from the grid boxes from the
523	!! data that go into the
524	!! model's boxes
525	!! cf src_global/interpol_help.f90,
526	!! line 300, called "ip"
527
528	! REAL(r_std), ALLOCATABLE, DIMENSION(:) :: invar1D !! 1D input values
529	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: invar2D !! 2D input values
530	REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:) :: invar3D !! 3D input values
531	REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:,:) :: invar4D !! 4D input values
532
533	REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:) :: resol_in
534	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: mask
535	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: maskvar
536	INTEGER(i_std) :: nix, njx
537	INTEGER(i_std) :: idi, nbvmax !! nbvmax : number of maximum vegetation map
538	!! points in the GCM grid ; idi : its counter
539	CHARACTER(LEN=30) :: callsign !! Allows to specify which variable is beeing
540	!! treated
541	LOGICAL :: ok_interpol !! optionnal return of aggregate_2d
542	!
543	INTEGER :: ALLOC_ERR
544	CHARACTER(LEN=50) :: S1
545	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: variableusetypes
546	CHARACTER(LEN=256) :: msg
547	INTEGER, DIMENSION(:), ALLOCATABLE :: indims
548
549	! Debug vars
550	INTEGER :: nb_coord, nb_var, nb_gat
551
552
553	!_ ================================================================================================================================
554	! Allocating input data
555
556	IF (printlev >= 3) WRITE(numout,*)"In interpweight_2D filename: ",TRIM(filename), &
557	" varname:'" // TRIM(varname) // "'"
558	!
559	IF (is_root_prc) THEN
560	IF (printlev >=5) THEN
561	WRITE(numout,*) "Entering 'interpweight_2D'. Debug mode."
562	WRITE(numout,*)'(/," --> fliodmpf")'
563	CALL fliodmpf (TRIM(filename))
564	WRITE(numout,*)'(/," --> flioopfd")'
565	ENDIF
566	CALL flioopfd (TRIM(filename),fid,nb_dim=nb_coord,nb_var=nb_var,nb_gat=nb_gat)
567	IF (printlev >=5) THEN
568	WRITE (numout,'(" Number of coordinate in the file : ",I2)') nb_coord
569	WRITE (numout,'(" Number of variables in the file : ",I2)') nb_var
570	WRITE (numout,'(" Number of global attributes in the file : ",I2)') nb_gat
571	ENDIF
572	ENDIF
573
574	! Getting shape of input variable from input file
575	inNdims = interpweight_get_varNdims_file(filename, varname)
576	IF (is_root_prc) CALL flininfo(filename, iml, jml, lml, tml, fid)
577
578	! Read variable (varname) from input file
579	InputDims: SELECT CASE (inNdims)
580	CASE (2)
581	CALL bcast(iml)
582	CALL bcast(jml)
583
584	ALLOCATE(invar2D(iml,jml), STAT=ALLOC_ERR)
585	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2D',"Problem in allocation of variable '" // &
586	TRIM(varname) ,'','')
587	IF (is_root_prc) CALL flinget(fid, TRIM(varname), iml, jml, 0, 0, 1, 1, invar2D)
588	IF (printlev >= 5) THEN
589	WRITE(numout,*)' interpweight_2D input data dimensions: '
590	WRITE(numout,*)' iml, jml:',iml, jml
591	END IF
592
593	CASE (3)
594	ALLOCATE(indims(3), STAT=ALLOC_ERR)
595	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2D','Problem in allocation of indims','','')
596	indims = interpweight_get_var3dims_file(filename, varname)
597
598	lml = indims(3)
599	CALL bcast(iml)
600	CALL bcast(jml)
601	CALL bcast(lml)
602
603	ALLOCATE(invar3D(iml,jml, lml), STAT=ALLOC_ERR)
604	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2D',"Problem in allocation of variable '" // &
605	TRIM(varname) ,'','')
606	IF (is_root_prc) CALL flinget(fid, TRIM(varname), iml, jml, lml, 0, 1, 1, invar3D)
607	IF (printlev >= 5) THEN
608	WRITE(numout,*)' interpweight_2D input data dimensions:'
609	WRITE(numout,*)' iml, jml, lml:',iml, jml, lml
610	END IF
611
612	CASE (4)
613	ALLOCATE(indims(4), STAT=ALLOC_ERR)
614	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2D','Problem in allocation of indims','','')
615	indims = interpweight_get_var4dims_file(filename, varname)
616	lml = indims(3)
617	CALL bcast(iml)
618	CALL bcast(jml)
619	CALL bcast(lml)
620	! Taken from `slowproc_update'
621	IF (initime /= -1) THEN
622	inNdims = 3
623	IF (printlev >= 3) WRITE(numout,*) &
624	'interpweight_2D: taking into account the initial reading time:',initime
625	IF (initime <= 0) tml = 0
626
627	CALL bcast(tml)
628	ALLOCATE(invar3D(iml,jml,lml), STAT=ALLOC_ERR)
629	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2D',"Problem in allocation of variable '" // &
630	TRIM(varname) ,'','')
631	IF ( initime > 0 ) THEN
632	IF (is_root_prc) THEN
633	IF (initime <= tml) THEN
634	CALL flinget(fid, TRIM(varname), iml, jml, lml, 1, initime, 1, invar3D)
635	ELSE
636	CALL flinget(fid, TRIM(varname), iml, jml, lml, 1, tml, 1, invar3D)
637	ENDIF
638	ENDIF
639	ELSE
640	IF (is_root_prc) THEN
641	CALL flinget(fid, TRIM(varname), iml, jml, lml, 1, 1, 1, invar3D)
642	ENDIF
643	ENDIF
644	ELSE
645	CALL bcast(tml)
646	ALLOCATE(invar4D(iml,jml,lml,tml), STAT=ALLOC_ERR)
647	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2D',"Problem in allocation of variable '" // &
648	TRIM(varname) ,'','')
649	IF (is_root_prc) CALL flinget(fid, TRIM(varname), iml, jml, lml, tml, 1, tml, invar4D)
650	END IF
651
652	IF (printlev >= 5) THEN
653	WRITE(numout,*)' interpweight_2D input data dimensions'
654	WRITE(numout,*)' iml, jml, lml, tml :',iml, jml, lml, tml
655	END IF
656
657	CASE DEFAULT
658	WRITE (S1,'(I3)') inNdims
659	CALL ipslerr_p(3,'interpweight_2D',"Input number of dimensions: '" // TRIM(S1) // "' not ready !!", &
660	'','')
661
662	END SELECT InputDims
663	CALL bcast(invar2D)
664
665	IF (printlev >= 5) WRITE(numout,*)' interpweight_2D getting variable ended!'
666
667	! Getting longitudes/latitudes from input file
668	inNLldims = interpweight_get_varNdims_file(filename, inlonname)
669
670	IF (inNLldims == 1) THEN
671	! Allocate memory for latitudes
672	ALLOCATE(lon_in(iml), STAT=ALLOC_ERR)
673	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2D','Problem in allocation of variable lon_in','','')
674
675	ALLOCATE(lat_in(jml), STAT=ALLOC_ERR)
676	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2D','Problem in allocation of variable lat_in','','')
677
678	IF (is_root_prc) THEN
679	CALL flinget(fid, TRIM(inlonname), iml, 0, 0, 0, 1, 1, lon_in)
680	CALL flinget(fid, TRIM(inlatname), jml, 0, 0, 0, 1, 1, lat_in)
681	ENDIF
682	CALL bcast(lon_in)
683	CALL bcast(lat_in)
684
685	ALLOCATE(lon(iml,jml), STAT=ALLOC_ERR)
686	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2D','Problem in allocation of variable lon','','')
687	ALLOCATE(lat(iml,jml), STAT=ALLOC_ERR)
688	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2D','Problem in allocation of variable lat','','')
689
690	DO ip=1,iml
691	lat(ip,:) = lat_in(:)
692	ENDDO
693	DO jp=1,jml
694	lon(:,jp) = lon_in(:)
695	ENDDO
696	ELSE
697	! Allocate memory for longitude
698	ALLOCATE(lon(iml,jml), STAT=ALLOC_ERR)
699	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2D','Pb in allocation for lon','','')
700	! Allocate memory for latitudes
701	ALLOCATE(lat(iml,jml), STAT=ALLOC_ERR)
702	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2D','Pb in allocation for lat','','')
703
704	IF (is_root_prc) THEN
705	CALL flinget(fid, TRIM(inlonname), iml, jml, 0, 0, 1, 1, lon)
706	CALL flinget(fid, TRIM(inlatname), iml, jml, 0, 0, 1, 1, lat)
707	ENDIF
708	END IF
709	CALL bcast(lon)
710	CALL bcast(lat)
711
712	IF (is_root_prc) THEN
713	ALLOCATE(resol_in(iml,jml,2), STAT=ALLOC_ERR)
714	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2D','Problem in allocation of variable resol_in','','')
715	END IF
716
717	IF (noneg) THEN
718	IF (is_root_prc) THEN
719	SELECT CASE (inNdims)
720	CASE(2)
721	CALL interpweight_modifying_input2D(iml, jml, 0, 0, 1, 1, zero, invar2D)
722	CASE(3)
723	CALL interpweight_modifying_input3D(iml, jml, lml, 0, 1, 1, zero, invar3D)
724	CASE(4)
725	CALL interpweight_modifying_input4D(iml, jml, lml, tml, 1, tml, zero, invar4D)
726	END SELECT
727	END IF
728	END IF
729	IF (printlev >= 5) WRITE(numout,*)' interpweight_2D modifying input ended!'
730
731	!
732	! Consider all points a priori
733	!
734	ALLOCATE(mask(iml,jml), STAT=ALLOC_ERR)
735	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2D','Problem in allocation of variable mask','','')
736	mask = zero
737
738	IF (TRIM(masktype) == 'var') THEN
739	ALLOCATE(maskvar(iml,jml), STAT=ALLOC_ERR)
740	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2D','Problem in allocation of variable maskvar','','')
741	maskvar(:,:) = 0
742	! Reads the mask from the file
743	IF (is_root_prc) THEN
744	! ATTTENTION!! One can not use flinget in this way. It always spect to have a REAL variable!
745	CALL flinget(fid, TRIM(maskvarname), iml, jml, 0, 0, 1, 1, maskvar)
746	WHERE (maskvar > zero)
747	mask = un
748	END WHERE
749	END IF
750	CALL bcast(mask)
751	ELSE
752	! Create land-sea mask based on values from the variable according to maskvalues and masktype
753	SELECT CASE (inNdims)
754	CASE(2)
755	CALL interpweight_masking_input2D(iml, jml, 0, 0, 0, 0, masktype, maskvalues, invar2D, &
756	mask)
757	CASE(3)
758	CALL interpweight_masking_input3D(iml, jml, 0, 0, lml, 0, masktype, maskvalues, invar3D, &
759	mask)
760	CASE(4)
761	CALL interpweight_masking_input4D(iml, jml, 0, 0, lml, tml, masktype, maskvalues, invar4D, &
762	mask)
763	END SELECT
764	END IF
765	IF (printlev >= 5) WRITE(numout,*)' interpweight_2D masking input ended!'
766
767	IF (is_root_prc) CALL flinclo(fid)
768
769	! Assuming regular lon/lat projection for the input data
770	IF (is_root_prc) THEN
771	CALL interpweight_calc_resolution_in(lon, lat, iml, jml, mincos, R_Earth, pi, resol_in)
772	END IF
773
774	! Computation of nbvmax
775	! In 'condveg_soilalb, slowproc_update' nbvmax=200
776	!
777	! The number of maximum vegetation map points in the GCM grid is estimated.
778	! Some lmargin is taken.
779	!
780	IF (is_root_prc) THEN
781	IF ( (maxresollon /= -un) .AND. (maxresollat /= -un) ) THEN
782	! Assuming regular lon/lat projection for the input data but it is not caluclated assuming a maximum &
783	! lon,lat resolution given as input paramater
784	nix=INT(MAXVAL(resolution(:,1))/maxresollon)+2
785	njx=INT(MAXVAL(resolution(:,2))/maxresollat)+2
786	ELSE
787	! Assuming regular lon/lat projection for the input data
788	CALL interpweight_calc_resolution_in(lon, lat, iml, jml, mincos, R_Earth, pi, resol_in)
789	nix=INT(MAXVAL(resolution(:,1))/MAXVAL(resol_in(:,:,1)))+2
790	njx=INT(MAXVAL(resolution(:,2))/MAXVAL(resol_in(:,:,2)))+2
791	END IF
792	nbvmax = nix*njx
793	IF (nbvmax < 200) nbvmax = 200
794	END IF
795	CALL bcast(nbvmax)
796	CALL bcast(nix)
797	CALL bcast(njx)
798
799	!
800	callsign = TRIM(varname) // ' map'
801	!
802	ok_interpol = .FALSE.
803	DO WHILE ( .NOT. ok_interpol )
804	IF (printlev>=3) WRITE(numout,*) "Projection arrays for ",callsign," : "
805	IF (printlev>=3) WRITE(numout,*) "nbvmax = ",nbvmax, nix, njx
806
807
808	ALLOCATE(sub_index(nbpt, nbvmax, 2), STAT=ALLOC_ERR)
809	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2D','Problem in allocation of variable sub_index','','')
810
811	sub_index(:,:,:)=zero
812
813	ALLOCATE(sub_area(nbpt, nbvmax), STAT=ALLOC_ERR)
814	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2D','Problem in allocation of variable sub_area','','')
815
816	sub_area(:,:)=zero
817	IF (printlev>=3) WRITE(numout,*) 'Input data range LON:', MINVAL(lon), MAXVAL(lon)
818	IF (printlev>=3) WRITE(numout,*) 'Input data range LAT:', MINVAL(lat), MAXVAL(lat)
819
820	!
821	CALL aggregate_p(nbpt, lalo, neighbours, resolution, contfrac, &
822	& iml, jml, lon, lat, mask, callsign, &
823	& nbvmax, sub_index, sub_area, ok_interpol)
824
825	!
826	IF ( .NOT. ok_interpol ) THEN
827	DEALLOCATE(sub_area)
828	DEALLOCATE(sub_index)
829	nbvmax = nbvmax * 2
830	ENDIF
831	ENDDO
832	IF (printlev>=4) WRITE (numout,*) ' interpweight_2D: read/allocate OK'
833
834	! Getting variables thresholds
835	ALLOCATE(variableusetypes(Nvariabletypes), STAT=ALLOC_ERR)
836	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2D','Problem in allocation of variableusetypes','','')
837
838
839	IF (variabletypes(1) == -un) THEN
840	variableusetypes = (/ (ip*1., ip=1,Nvariabletypes) /)
841	IF (printlev>=4) THEN
842	WRITE(numout,*)' interpweight_2D: Note: Using default equivalence between dimension ',&
843	'in the file and dimension in the variable, for file ', filename
844	END IF
845	varmin = 1._r_std
846	varmax = Nvariabletypes*1._r_std
847	ELSE
848	variableusetypes = variabletypes
849	END IF
850
851	outvar2D = zero
852
853	! Providing the fractions of each type for each interpolated grid point
854	CALL interpweight_provide_fractions2D(typefrac, nbpt, Nvariabletypes, variableusetypes, &
855	iml, jml, lml, tml, nbvmax, zero, invar2D, sub_area, sub_index, varmax, varmin, &
856	deux, lalo, outvar2D, aoutvar)
857	IF (printlev>=5) WRITE (numout,*)' interpweight_2D end of provide_fractions2D'
858
859	IF (printlev>=3) THEN
860	IF (COUNT(aoutvar == -1) >= AINT(nbpt*0.4)) THEN
861	WRITE(numout,*) ' interpweight_2D total number of points: ', nbpt
862	WRITE(numout,*) ' interpweight_2D interpolated: ', COUNT(aoutvar /= -1), &
863	' non-interpolated: ', COUNT(aoutvar == -1)
864	msg = "Less than 60% of points '" // TRIM(filename) // "' variable '" // TRIM(varname) // &
865	"' have been succesfully interpolated !!"
866	WRITE(numout,*) ' interpweight_2D: ' // TRIM(msg)
867	ELSE
868	msg = TRIM(filename) // "' variable: '" // TRIM(varname) // "' interpolation Done _______"
869	WRITE(numout,*) "interpweight_2D: '" // TRIM(msg)
870	WRITE(numout,*) ' Total number of points: ', nbpt
871	WRITE(numout,*) ' Interpolated: ', COUNT(aoutvar /= -1), &
872	' Non-interpolated: ', COUNT(aoutvar == -1)
873	END IF
874	END IF
875
876	! Scaling to 1 the quality variable
877	DO ip=1, nbpt
878	IF (aoutvar(ip) /= -1) aoutvar(ip) = aoutvar(ip)/(resolution(ip,1)*resolution(ip,2))/contfrac(ip)
879	END DO
880
881	IF (ALLOCATED(lon_in)) DEALLOCATE(lon_in)
882	IF (ALLOCATED(lat_in)) DEALLOCATE(lat_in)
883	DEALLOCATE(lon)
884	DEALLOCATE(lat)
885	IF (ALLOCATED(invar2D)) DEALLOCATE(invar2D)
886	IF (ALLOCATED(invar3D)) DEALLOCATE(invar3D)
887	IF (ALLOCATED(invar4D)) DEALLOCATE(invar4D)
888	DEALLOCATE(mask)
889	IF (ALLOCATED(maskvar)) DEALLOCATE(maskvar)
890	DEALLOCATE(sub_area)
891	DEALLOCATE(sub_index)
892	IF (ALLOCATED(resol_in)) DEALLOCATE(resol_in)
893
894	RETURN
895
896	END SUBROUTINE interpweight_2D
897
898	!! ================================================================================================================================
899	!! SUBROUTINE : interpweight_3D
900	!!
901	!>\BRIEF ! Interpolate any input file to the grid of the model for a 3D-variable
902	!!
903	!! DESCRIPTION : (definitions, functional, design, flags):
904	!!
905	!! RECENT CHANGE(S): None
906	!!
907	!! MAIN OUTPUT VARIABLE(S): outvar3D, aoutvar
908	!!
909	!! REFERENCE(S) : None
910	!!
911	!! FLOWCHART : None
912	!! \n
913	!_ ================================================================================================================================
914	SUBROUTINE interpweight_3D(nbpt, Nvariabletypes, variabletypes, lalo, resolution, neighbours, &
915	contfrac, filename, varname, inlonname, inlatname, varmin, varmax, noneg, masktype, &
916	maskvalues, maskvarname, dim1, dim2, initime, typefrac, &
917	maxresollon, maxresollat, outvar3D, aoutvar)
918
919	USE ioipsl
920	USE constantes_var
921	USE mod_orchidee_para
922
923	IMPLICIT NONE
924
925	!! 0. Variables and parameter declaration
926	!
927	! 0.1 INPUT
928	!
929	INTEGER(i_std), INTENT(in) :: nbpt !! Number of points for which the data
930	!! needs to be interpolated
931	INTEGER(i_std), INTENT(in) :: Nvariabletypes !! Number of types of the variable
932	REAL(r_std), DIMENSION(Nvariabletypes), & !! Vector of values of the types
933	INTENT(in) :: variabletypes !! (-1, not used)
934	REAL(r_std), DIMENSION(nbpt,2), INTENT(in) :: lalo !! Vector of latitude and longitudes
935	!! (beware of the order = 1 : latitude,
936	!! 2 : longitude)
937	REAL(r_std), DIMENSION(nbpt,2), INTENT(in) :: resolution !! The size in km of each grid-box in X and Y
938	INTEGER(i_std), DIMENSION(nbpt,NbNeighb), INTENT(in) :: neighbours !! Vector of neighbours for each grid point
939	!! 1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW)
940
941	REAL(r_std), DIMENSION(nbpt), INTENT(in) :: contfrac !! Fraction of land in each grid box.
942	CHARACTER(LEN=80), INTENT(in) :: filename !! name of the input map to read
943
944	CHARACTER(LEN=80), INTENT(in) :: varname !! name of the variable to interpolate
945	CHARACTER(LEN=80), INTENT(in) :: inlonname, inlatname !! names of the longitude and latitude
946	!! in the input file
947	! Transformations to apply to the original read values from the file
948	LOGICAL, INTENT(in) :: noneg !! no negative values
949	CHARACTER(LEN=50), INTENT(in) :: masktype !! Type of masking
950	!! 'nomask': no-mask is applied
951	!! 'mbelow': take values below maskvalues(1)
952	!! 'mabove': take values above maskvalues(1)
953	!! 'msumrange': take values within 2 ranges;
954	!! maskvalues(2) <= SUM(vals(k)) <= maskvalues(1)
955	!! maskvalues(1) < SUM(vals(k)) <= maskvalues(3)
956	!! (normalized by maskvalues(3))
957	!! 'var': mask values are taken from a
958	!! variable (>0)
959	REAL(r_std), DIMENSION(3), INTENT(in) :: maskvalues !! values to use to mask (according to
960	!! `masktype')
961	CHARACTER(LEN=250), INTENT(in) :: maskvarname !! name of the variable to use to mask
962	INTEGER, INTENT(in) :: dim1, dim2 !! 3/4D dimensions of the output variable
963	REAL(r_std), DIMENSION(dim1), INTENT(inout) :: varmin, varmax !! min/max values to use for the renormalization
964	CHARACTER(LEN=50), INTENT(in) :: typefrac !! Type of fraction retrieval:
965	!! 'default': standard
966	INTEGER(i_std), INTENT(in) :: initime !! Initial time to take (for input data
967	!! with time values, -1 not used)
968	REAL(r_std), INTENT(in) :: maxresollon,maxresollat !! lon,lat maximum resolutions (in m)
969	!! (-un, not used)
970	!! 0.2 Modified variables
971	!
972	! 0.3 OUTPUT
973	!
974	REAL(r_std), DIMENSION(nbpt,dim1), INTENT(out) :: outvar3D !! 3D output variable and re-dimensioned
975	! aovar: availability of input data to interpolate output variable (on the nbpt space)
976	REAL(r_std), DIMENSION(nbpt), INTENT(out) :: aoutvar !! availability of input data to
977	!! interpolate output variable
978	!! (on the nbpt space)
979	!
980	! 0.4 Local variables
981	!
982	INTEGER(i_std) :: iml, jml, lml, tml, fid, ib, ip, jp, it, jj, jv
983	INTEGER :: inNdims !! Number of dimensions of the input data
984	INTEGER :: inNLldims !! Number of dimensions of the lon/lat input data
985	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lat_in, lon_in !! latitude and
986	!! longitude, extract from input map
987	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lat, lon !! en 2D ???
988	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: sub_area !! the area of the fine grid in the
989	!! model grid ???
990	!! cf src_global/interpol_help.f90,
991	!! line 377, called "areaoverlap"
992	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:,:) :: sub_index !! the indexes from the grid boxes from the
993	!! data that go into the
994	!! model's boxes
995	!! cf src_global/interpol_help.f90,
996	!! line 300, called "ip"
997
998	! REAL(r_std), ALLOCATABLE, DIMENSION(:) :: invar1D !! 1D input values
999	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: invar2D !! 2D input values
1000	REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:) :: invar3D !! 3D input values
1001	REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:,:) :: invar4D !! 4D input values
1002
1003	REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:) :: resol_in
1004	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: mask
1005	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: maskvar
1006	INTEGER(i_std) :: nix, njx
1007	INTEGER(i_std) :: idi, nbvmax !! nbvmax : number of maximum vegetation map
1008	!! points in the GCM grid ; idi : its counter
1009	CHARACTER(LEN=30) :: callsign !! Allows to specify which variable is beeing
1010	!! treated
1011	LOGICAL :: ok_interpol !! optionnal return of aggregate_2d
1012	INTEGER :: ALLOC_ERR
1013
1014	CHARACTER(LEN=50) :: S1
1015	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: variableusetypes
1016	CHARACTER(LEN=256) :: msg
1017
1018	! Debug vars
1019	INTEGER :: nb_coord, nb_var, nb_gat
1020	CHARACTER(LEN=5) :: iS
1021	CHARACTER(LEN=50) :: suminS
1022
1023	INTEGER, ALLOCATABLE, DIMENSION(:) :: indims
1024
1025
1026	!_ ================================================================================================================================
1027	IF (printlev >= 3) WRITE(numout, *)"In interpweight_3D filename: ",TRIM(filename), &
1028	" varname:'" // TRIM(varname) // "'"
1029	! Allocating input data
1030	IF (is_root_prc) THEN
1031	IF (printlev >=5) THEN
1032	WRITE(numout,*) "Entering 'interpweight_3D'. Debug mode."
1033	WRITE(numout,*)'(/," --> fliodmpf")'
1034	CALL fliodmpf (TRIM(filename))
1035	WRITE(numout,*)'(/," --> flioopfd")'
1036	ENDIF
1037	CALL flioopfd (TRIM(filename),fid,nb_dim=nb_coord,nb_var=nb_var,nb_gat=nb_gat)
1038	IF (printlev >=5) THEN
1039	WRITE (numout,'(" Number of coordinate in the file : ",I2)') nb_coord
1040	WRITE (numout,'(" Number of variables in the file : ",I2)') nb_var
1041	WRITE (numout,'(" Number of global attributes in the file : ",I2)') nb_gat
1042	ENDIF
1043	ENDIF
1044
1045	! Getting shape of input variable from input file
1046	inNdims = interpweight_get_varNdims_file(filename, varname)
1047	IF (is_root_prc) CALL flininfo(filename, iml, jml, lml, tml, fid)
1048
1049	InputDims: SELECT CASE (inNdims)
1050	CASE (2)
1051	CALL bcast(iml)
1052	CALL bcast(jml)
1053
1054	ALLOCATE(invar2D(iml,jml), STAT=ALLOC_ERR)
1055	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_3D',"Problem in allocation of variable '" // &
1056	TRIM(varname) ,'','')
1057	IF (is_root_prc) CALL flinget(fid, TRIM(varname), iml, jml, 0, 0, 1, 1, invar2D)
1058	IF (printlev >= 5) THEN
1059	WRITE(numout,*)' interpweight_3D input data dimensions'
1060	WRITE(numout,*)' iml, jml:',iml, jml
1061	END IF
1062
1063	CASE (3)
1064	ALLOCATE(indims(3), STAT=ALLOC_ERR)
1065	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_3D','Problem in allocation of indims','','')
1066	indims = interpweight_get_var3dims_file(filename, varname)
1067
1068	lml = indims(3)
1069	CALL bcast(iml)
1070	CALL bcast(jml)
1071	CALL bcast(lml)
1072
1073	ALLOCATE(invar3D(iml,jml, lml), STAT=ALLOC_ERR)
1074	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_3D',"Problem in allocation of variable '" // &
1075	TRIM(varname) ,'','')
1076	IF (is_root_prc) CALL flinget(fid, TRIM(varname), iml, jml, lml, 0, 1, 1, invar3D)
1077	IF (printlev >= 5) THEN
1078	WRITE(numout,*)' interpweight_3D input data dimensions'
1079	WRITE(numout,*)' iml, jml, lml :',iml, jml, lml
1080	END IF
1081
1082	CASE (4)
1083	ALLOCATE(indims(4), STAT=ALLOC_ERR)
1084	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_3D','Problem in allocation of indims','','')
1085
1086	indims = interpweight_get_var4dims_file(filename, varname)
1087	lml = indims(3)
1088	CALL bcast(iml)
1089	CALL bcast(jml)
1090	CALL bcast(lml)
1091	! Taken from `slowproc_update'
1092	IF (initime /= -1) THEN
1093	inNdims = 3
1094	IF (printlev >= 3) WRITE (numout,*) &
1095	'interpweight_3D: taking into account the initial reading time:',initime
1096	IF (initime <= 0) tml = 0
1097
1098	CALL bcast(tml)
1099	ALLOCATE(invar3D(iml,jml,lml), STAT=ALLOC_ERR)
1100	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_3D','Problem in allocation of variable invar3D','','')
1101	IF ( initime > 0 ) THEN
1102	IF (is_root_prc) THEN
1103	IF (initime <= tml) THEN
1104	CALL flinget(fid, TRIM(varname), iml, jml, lml, 1, initime, 1, invar3D)
1105	ELSE
1106	CALL flinget(fid, TRIM(varname), iml, jml, lml, 1, tml, 1, invar3D)
1107	ENDIF
1108	ENDIF
1109	ELSE
1110	IF (is_root_prc) THEN
1111	CALL flinget(fid, TRIM(varname), iml, jml, lml, 1, 1, 1, invar3D)
1112	ENDIF
1113	ENDIF
1114	ELSE
1115	CALL bcast(tml)
1116	ALLOCATE(invar4D(iml,jml,lml,tml), STAT=ALLOC_ERR)
1117	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_3D','Problem in allocation of variable invar4D','','')
1118	IF (is_root_prc) CALL flinget(fid, TRIM(varname), iml, jml, lml, tml, 1, tml, invar4D)
1119	END IF
1120
1121	IF (printlev >= 5) THEN
1122	WRITE(numout,*)' interpweight_3D input data dimensions'
1123	WRITE(numout,*)' iml, jml, lml, tml:',iml, jml, lml, tml
1124	END IF
1125
1126	CASE DEFAULT
1127	WRITE (S1,'(I3)') inNdims
1128	CALL ipslerr_p(3,'interpweight_3D',"Input number of dimensions: '" // TRIM(S1) // "' not ready !!", &
1129	'','')
1130
1131	END SELECT InputDims
1132	CALL bcast(invar3D)
1133
1134	IF (printlev >= 5) WRITE(numout,*)' interpweight_3D getting variable ended!'
1135	! Test lecture
1136	DO ip=1,lml
1137	IF ( SUM(invar3D(:,:,ip)) < 0.00001 ) THEN
1138	WRITE(iS, '(i5)')ip
1139	WRITE(suminS, '(f20.10)')SUM(invar3d(:,:,ip))
1140	msg = 'Wrong lecture of data for type ' // TRIM(iS) // ' all values are zero !!'
1141	CALL ipslerr_p(3,'interpweight_3D',TRIM(msg),'','')
1142	END IF
1143	END DO
1144
1145	! Getting longitudes/latitudes from input file
1146	inNLldims = interpweight_get_varNdims_file(filename, inlonname)
1147	IF (inNLldims == 1) THEN
1148	! Allocate memory for latitudes
1149	ALLOCATE(lon_in(iml), STAT=ALLOC_ERR)
1150	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_3D','Problem in allocation of variable lon_in','','')
1151
1152	ALLOCATE(lat_in(jml), STAT=ALLOC_ERR)
1153	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_3D','Problem in allocation of variable lat_in','','')
1154
1155	IF (is_root_prc) THEN
1156	CALL flininspect(fid)
1157	CALL flinget(fid, TRIM(inlonname), iml, 0, 0, 0, 1, 1, lon_in)
1158	CALL flinget(fid, TRIM(inlatname), jml, 0, 0, 0, 1, 1, lat_in)
1159	ENDIF
1160	CALL bcast(lon_in)
1161	CALL bcast(lat_in)
1162
1163	ALLOCATE(lon(iml,jml), STAT=ALLOC_ERR)
1164	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_3D','Problem in allocation of variable lon','','')
1165	ALLOCATE(lat(iml,jml), STAT=ALLOC_ERR)
1166	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_3D','Problem in allocation of variable lat','','')
1167
1168	DO ip=1,iml
1169	lat(ip,:) = lat_in(:)
1170	ENDDO
1171	DO jp=1,jml
1172	lon(:,jp) = lon_in(:)
1173	ENDDO
1174	ELSE
1175	! Allocate memory for longitude
1176	ALLOCATE(lon(iml,jml), STAT=ALLOC_ERR)
1177	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_3D','Pb in allocation for lon','','')
1178	! Allocate memory for latitudes
1179	ALLOCATE(lat(iml,jml), STAT=ALLOC_ERR)
1180	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_3D','Pb in allocation for lat','','')
1181
1182	IF (is_root_prc) THEN
1183	CALL flinget(fid, TRIM(inlonname), iml, jml, 0, 0, 1, 1, lon)
1184	CALL flinget(fid, TRIM(inlatname), iml, jml, 0, 0, 1, 1, lat)
1185	ENDIF
1186	END IF
1187	CALL bcast(lon)
1188	CALL bcast(lat)
1189
1190	ALLOCATE(resol_in(iml,jml,2), STAT=ALLOC_ERR)
1191	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_3D','Problem in allocation of variable resol_in','','')
1192
1193	IF (noneg) THEN
1194	IF (is_root_prc) THEN
1195	SELECT CASE (inNdims)
1196	CASE(2)
1197	CALL interpweight_modifying_input2D(iml, jml, 0, 0, 1, 1, zero, invar2D)
1198	CASE(3)
1199	CALL interpweight_modifying_input3D(iml, jml, lml, 0, 1, 1, zero, invar3D)
1200	CASE(4)
1201	CALL interpweight_modifying_input4D(iml, jml, lml, tml, 1, tml, zero, invar4D)
1202	END SELECT
1203	END IF
1204	END IF
1205	IF (printlev >= 5) WRITE(numout,*)' interpweight_3D modifying input ended!'
1206
1207	ALLOCATE(mask(iml,jml), STAT=ALLOC_ERR)
1208	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_3D','Problem in allocation of variable mask','','')
1209	mask(:,:) = 0
1210
1211	IF (TRIM(masktype) == 'var') THEN
1212	ALLOCATE(maskvar(iml,jml), STAT=ALLOC_ERR)
1213	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_3D','Problem in allocation of variable maskvar','','')
1214	maskvar(:,:) = 0.
1215	! Reads the mask from the file
1216	IF (is_root_prc) THEN
1217	! ATTTENTION!! One can not use flinget in this way. It always spect to have a REAL variable!
1218	CALL flinget(fid, TRIM(maskvarname), iml, jml, 0, 0, 1, 1, maskvar)
1219	WHERE (maskvar > zero)
1220	mask = un
1221	END WHERE
1222	END IF
1223	CALL bcast(mask)
1224	ELSE
1225	SELECT CASE (inNdims)
1226	CASE(2)
1227	CALL interpweight_masking_input2D(iml, jml, 0, 0, 0, 0, masktype, maskvalues, invar2D, &
1228	mask)
1229	CASE(3)
1230	CALL interpweight_masking_input3D(iml, jml, 0, 0, lml, 0, masktype, maskvalues, invar3D, &
1231	mask)
1232	CASE(4)
1233	CALL interpweight_masking_input4D(iml, jml, 0, 0, lml, tml, masktype, maskvalues, invar4D, &
1234	mask)
1235	END SELECT
1236	END IF
1237	IF (printlev >= 5) WRITE(numout,*)' interpweight_3D masking input ended!'
1238
1239	IF (is_root_prc) CALL flinclo(fid)
1240
1241	! Computation of nbvmax
1242	! In 'condveg_soilalb, slowproc_update' nbvmax=200
1243	!
1244	! The number of maximum vegetation map points in the GCM grid is estimated.
1245	! Some lmargin is taken.
1246	!
1247	IF (is_root_prc) THEN
1248	IF ( (maxresollon /= -un) .AND. (maxresollat /= -un) ) THEN
1249	! Assuming regular lon/lat projection for the input data but it is not caluclated assuming a maximum &
1250	! lon,lat resolution given as input paramater
1251	nix=INT(MAXVAL(resolution(:,1))/maxresollon)+2
1252	njx=INT(MAXVAL(resolution(:,2))/maxresollat)+2
1253	ELSE
1254	! Assuming regular lon/lat projection for the input data
1255	CALL interpweight_calc_resolution_in(lon, lat, iml, jml, mincos, R_Earth, pi, resol_in)
1256	nix=INT(MAXVAL(resolution(:,1))/MAXVAL(resol_in(:,:,1)))+2
1257	njx=INT(MAXVAL(resolution(:,2))/MAXVAL(resol_in(:,:,2)))+2
1258	END IF
1259	nbvmax = nix*njx
1260	IF (nbvmax < 200) nbvmax = 200
1261	END IF
1262	CALL bcast(nbvmax)
1263	CALL bcast(nix)
1264	CALL bcast(njx)
1265
1266	!
1267	callsign = TRIM(varname) // ' map'
1268	!
1269	ok_interpol = .FALSE.
1270	DO WHILE ( .NOT. ok_interpol )
1271	IF (printlev>=3) WRITE(numout,*) "Projection arrays for ",callsign," : "
1272	IF (printlev>=3) WRITE(numout,*) "nbvmax = ",nbvmax, nix, njx
1273
1274	ALLOCATE(sub_index(nbpt, nbvmax, 2), STAT=ALLOC_ERR)
1275	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_3D','Problem in allocation of variable sub_index','','')
1276
1277	sub_index(:,:,:)=zero
1278
1279	ALLOCATE(sub_area(nbpt, nbvmax), STAT=ALLOC_ERR)
1280	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_3D','Problem in allocation of variable sub_area','','')
1281
1282	sub_area(:,:)=zero
1283
1284	IF (printlev >= 3) THEN
1285	WRITE(numout,*) 'interpweight_3D:'
1286	WRITE(numout,*) ' Carteveg range LON:', MINVAL(lon), MAXVAL(lon)
1287	WRITE(numout,*) ' Carteveg range LAT:', MINVAL(lat), MAXVAL(lat)
1288	END IF
1289
1290	CALL aggregate_p(nbpt, lalo, neighbours, resolution, contfrac, &
1291	& iml, jml, lon, lat, mask, callsign, &
1292	& nbvmax, sub_index, sub_area, ok_interpol)
1293
1294	IF ( .NOT. ok_interpol ) THEN
1295	DEALLOCATE(sub_area)
1296	DEALLOCATE(sub_index)
1297	nbvmax = nbvmax * 2
1298	ENDIF
1299	ENDDO
1300
1301	! Getting variables thresholds
1302	ALLOCATE(variableusetypes(Nvariabletypes),STAT=ALLOC_ERR)
1303	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_3D','Problem in allocation of variable variableusetypes','','')
1304
1305	IF (variabletypes(1) == -un) THEN
1306	variableusetypes = (/ (ip*1., ip=1,Nvariabletypes) /)
1307	IF (printlev>=4) WRITE(numout,*)'interpweight_3D: Note: Using default equivalence between ',&
1308	'dimension in the file and dimension in the variable, for file ', filename
1309	varmin = 1._r_std
1310	varmax = Nvariabletypes*1._r_std
1311	ELSE
1312	variableusetypes = variabletypes
1313	END IF
1314
1315	outvar3D = zero
1316
1317	! Providing the fractions of each type for each interpolated grid point
1318	CALL interpweight_provide_fractions3D(typefrac, nbpt, Nvariabletypes, variableusetypes, &
1319	iml, jml, lml, tml, nbvmax, zero, invar3D, sub_area, sub_index, varmax, varmin, &
1320	deux, lalo, outvar3D, aoutvar)
1321	IF (printlev >= 5) WRITE(numout,*)'interpweight_3D end of interpweight_provide_fractions3D'
1322
1323	IF (printlev>=2) THEN
1324	IF (COUNT(aoutvar == -1) >= AINT(nbpt*0.4)) THEN
1325	WRITE(numout,*) 'interpweight_3D total number of points: ', nbpt
1326	WRITE(numout,*) ' Interpolated: ', COUNT(aoutvar /= -1), &
1327	' Non-interpolated: ', COUNT(aoutvar == -1)
1328	msg = "Less than 60% of points '" // TRIM(filename) // "' variable '" // TRIM(varname) // &
1329	"' have been succesfully interpolated on the current process !!"
1330	WRITE(numout,*)' ' // TRIM(msg)
1331	ELSE
1332	msg = TRIM(filename) // "' variable: '" // TRIM(varname) // "' interpolation Done"
1333	WRITE(numout,*) "interpweight_3D: '" // TRIM(msg)
1334	WRITE(numout,*) ' Total number of points: ', nbpt
1335	WRITE(numout,*) ' Interpolated: ', COUNT(aoutvar /= -1), &
1336	' Non-interpolated: ', COUNT(aoutvar == -1)
1337	END IF
1338	END IF
1339
1340	! Scaling to 1 the quality variable
1341	DO ip=1, nbpt
1342	IF (aoutvar(ip) /= -1) aoutvar(ip) = aoutvar(ip)/(resolution(ip,1)*resolution(ip,2))/contfrac(ip)
1343	END DO
1344
1345	IF (ALLOCATED(lat_in)) DEALLOCATE(lat_in)
1346	IF (ALLOCATED(lon_in)) DEALLOCATE(lon_in)
1347	DEALLOCATE(lon)
1348	DEALLOCATE(lat)
1349	IF (ALLOCATED(invar2D)) DEALLOCATE(invar2D)
1350	IF (ALLOCATED(invar3D)) DEALLOCATE(invar3D)
1351	IF (ALLOCATED(invar4D)) DEALLOCATE(invar4D)
1352	DEALLOCATE(mask)
1353	IF (ALLOCATED(maskvar)) DEALLOCATE(maskvar)
1354	DEALLOCATE(sub_area)
1355	DEALLOCATE(sub_index)
1356	IF (ALLOCATED(resol_in)) DEALLOCATE(resol_in)
1357
1358	RETURN
1359
1360	END SUBROUTINE interpweight_3D
1361
1362	!! ================================================================================================================================
1363	!! SUBROUTINE : interpweight_4D
1364	!!
1365	!>\BRIEF ! Interpolate any input file to the grid of the model for a 4D-variable
1366	!!
1367	!! DESCRIPTION : (definitions, functional, design, flags):
1368	!!
1369	!! RECENT CHANGE(S): None
1370	!!
1371	!! MAIN OUTPUT VARIABLE(S): outvar4D, aoutvar
1372	!!
1373	!! REFERENCE(S) : None
1374	!!
1375	!! FLOWCHART : None
1376	!! \n
1377	!_ ================================================================================================================================
1378	SUBROUTINE interpweight_4D(nbpt, Nvariabletypes, variabletypes, lalo, resolution, neighbours, &
1379	contfrac, filename, varname, inlonname, inlatname, varmin, varmax, noneg, masktype, &
1380	maskvalues, maskvarname, dim1, dim2, initime, typefrac, &
1381	maxresollon, maxresollat, outvar4D, aoutvar)
1382
1383	USE ioipsl
1384	USE constantes_var
1385	USE mod_orchidee_para
1386
1387	IMPLICIT NONE
1388
1389	!! 0. Variables and parameter declaration
1390	!
1391	!! 0.1 Input variables
1392	!
1393	INTEGER(i_std), INTENT(in) :: nbpt !! Number of points for which the data
1394	!! needs to be interpolated
1395	INTEGER(i_std), INTENT(in) :: Nvariabletypes !! Number of types of the variable
1396	REAL(r_std), DIMENSION(Nvariabletypes), & !! Vector of values of the types
1397	INTENT(in) :: variabletypes !! (-1, not used)
1398	REAL(r_std), DIMENSION(nbpt,2), INTENT(in) :: lalo !! Vector of latitude and longitudes
1399	!! (beware of the order = 1 : latitude,
1400	!! 2 : longitude)
1401	REAL(r_std), DIMENSION(nbpt,2), INTENT(in) :: resolution !! The size in km of each grid-box in X and Y
1402	INTEGER(i_std), DIMENSION(nbpt,NbNeighb), INTENT(in) :: neighbours !! Vector of neighbours for each grid point
1403	!! 1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW)
1404
1405	REAL(r_std), DIMENSION(nbpt), INTENT(in) :: contfrac !! Fraction of land in each grid box.
1406	CHARACTER(LEN=80), INTENT(in) :: filename !! name of the input map to read
1407
1408	CHARACTER(LEN=80), INTENT(in) :: varname !! name of the variable to interpolate
1409	CHARACTER(LEN=80), INTENT(in) :: inlonname, inlatname !! names of the longitude and latitude
1410	! Transformations to apply to the original read values from the file
1411	LOGICAL, INTENT(in) :: noneg !! no negative values
1412	CHARACTER(LEN=50), INTENT(in) :: masktype !! Type of masking
1413	!! 'nomask': no-mask is applied
1414	!! 'mbelow': take values below maskvalues(1)
1415	!! 'mabove': take values above maskvalues(1)
1416	!! 'msumrange': take values within 2 ranges;
1417	!! maskvalues(2) <= SUM(vals(k)) <= maskvalues(1)
1418	!! maskvalues(1) < SUM(vals(k)) <= maskvalues(3)
1419	!! (normalized by maskedvalues(3))
1420	!! 'var': mask values are taken from a
1421	!! variable (>0)
1422	REAL(r_std), DIMENSION(3), INTENT(in) :: maskvalues !! values to use to mask (according to
1423	!! `masktype')
1424	CHARACTER(LEN=250), INTENT(in) :: maskvarname !! name of the variable to use to mask
1425	INTEGER, INTENT(in) :: dim1, dim2 !! 3/4D dimensions of the output variable
1426	!! dim1 = 1 fpr outvar2D
1427	!! in the input file
1428	REAL(r_std), DIMENSION(dim1), INTENT(inout) :: varmin, varmax !! min/max values to use for the renormalization
1429	CHARACTER(LEN=50), INTENT(in) :: typefrac !! Type of fractions retrieval:
1430	!! 'default': standard
1431
1432	INTEGER(i_std), INTENT(in) :: initime !! Initial time to take (for input data
1433	!! with time values, -1 not used)
1434	REAL(r_std), INTENT(in) :: maxresollon,maxresollat !! lon,lat maximum resolutions (in m)
1435	!! (-un, not used)
1436	!! 0.2 Modified variables
1437	!
1438	!! 0.3 Output variables
1439	!
1440	! REAL(r_std), INTENT(out) :: laimap(nbpt,nvm,12) !! lai read variable and re-dimensioned
1441	REAL(r_std), DIMENSION(nbpt,dim1,dim2), INTENT(out) :: outvar4D !! 4D output variable and re-dimensioned
1442	REAL(r_std), DIMENSION(nbpt), INTENT(out) :: aoutvar !! availability of input data to
1443	!! interpolate output variable
1444	!! (on the nbpt space)
1445	!
1446	!! 0.4 Local variables
1447	!
1448	INTEGER(i_std) :: iml, jml, lml, tml, fid, ib, ip, jp, it, jj, jv
1449	INTEGER :: inNdims !! Number of dimensions of the input data
1450	INTEGER :: inNLldims !! Number of dimensions of the lon/lat input data
1451	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lat_in, lon_in !! latitude and
1452	!! longitude, extract from input map
1453	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lat, lon !! en 2D ???
1454	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: sub_area !! the area of the fine grid in the
1455	!! model grid ???
1456	!! cf src_global/interpol_help.f90,
1457	!! line 377, called "areaoverlap"
1458	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:,:) :: sub_index !! the indexes from the grid boxes from the
1459	!! data that go into the
1460	!! model's boxes
1461	!! cf src_global/interpol_help.f90,
1462	!! line 300, called "ip"
1463
1464	INTEGER :: ALLOC_ERR
1465	! REAL(r_std), ALLOCATABLE, DIMENSION(:) :: invar1D !! 1D input values
1466	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: invar2D !! 2D input values
1467	REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:) :: invar3D !! 3D input values
1468	REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:,:) :: invar4D !! 4D input values
1469
1470	REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:) :: resol_in
1471	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: mask
1472	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: maskvar
1473	INTEGER(i_std) :: nix, njx
1474	INTEGER(i_std) :: idi, nbvmax !! nbvmax : number of maximum vegetation map
1475	!! points in the GCM grid ; idi : its counter
1476	CHARACTER(LEN=30) :: callsign !! Allows to specify which variable is beeing
1477	!! treated
1478	LOGICAL :: ok_interpol !! optionnal return of aggregate_2d
1479	CHARACTER(LEN=50) :: S1
1480	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: variableusetypes
1481	CHARACTER(LEN=256) :: msg
1482
1483	! Debug vars
1484	INTEGER :: nb_coord, nb_var, nb_gat
1485
1486	INTEGER, ALLOCATABLE, DIMENSION(:) :: indims
1487
1488
1489	!_ ================================================================================================================================
1490	IF (printlev >= 3) WRITE(numout, *)"In interpweight_4D filename: ",TRIM(filename), &
1491	" varname:'" // TRIM(varname) // "'"
1492
1493	! Allocating input data
1494	IF (is_root_prc) THEN
1495	IF (printlev >=5) THEN
1496	WRITE(numout,*) "Entering 'interpweight_4D'. Debug mode."
1497	WRITE(numout,*)'(/," --> fliodmpf")'
1498	CALL fliodmpf (TRIM(filename))
1499	WRITE(numout,*)'(/," --> flioopfd")'
1500	ENDIF
1501	CALL flioopfd (TRIM(filename),fid,nb_dim=nb_coord,nb_var=nb_var,nb_gat=nb_gat)
1502	IF (printlev >=5) THEN
1503	WRITE (numout,'(" Number of coordinate in the file : ",I2)') nb_coord
1504	WRITE (numout,'(" Number of variables in the file : ",I2)') nb_var
1505	WRITE (numout,'(" Number of global attributes in the file : ",I2)') nb_gat
1506	ENDIF
1507	ENDIF
1508
1509	! Getting shape of input variable from input file
1510	inNdims = interpweight_get_varNdims_file(filename, varname)
1511	IF (is_root_prc) CALL flininfo(filename, iml, jml, lml, tml, fid)
1512
1513	InputDims: SELECT CASE (inNdims)
1514	CASE (2)
1515	CALL bcast(iml)
1516	CALL bcast(jml)
1517
1518	ALLOCATE(invar2D(iml,jml), STAT=ALLOC_ERR)
1519	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4D',"Problem in allocation of variable '" // &
1520	TRIM(varname) ,'','')
1521	IF (is_root_prc) CALL flinget(fid, TRIM(varname), iml, jml, 0, 0, 1, 1, invar2D)
1522	IF (printlev >= 5) THEN
1523	WRITE(numout,*)' interpweight_4D input data dimensions _______'
1524	WRITE(numout,*)' iml, jml:',iml, jml
1525	END IF
1526
1527	CASE (3)
1528	ALLOCATE(indims(3), STAT=ALLOC_ERR)
1529	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4D','Problem in allocation of indims','','')
1530	indims = interpweight_get_var3dims_file(filename, varname)
1531	lml = indims(3)
1532	CALL bcast(iml)
1533	CALL bcast(jml)
1534	CALL bcast(lml)
1535
1536	ALLOCATE(invar3D(iml,jml, lml), STAT=ALLOC_ERR)
1537	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4D',"Problem in allocation of variable '" // &
1538	TRIM(varname) ,'','')
1539	IF (is_root_prc) CALL flinget(fid, TRIM(varname), iml, jml, lml, 0, 1, 1, invar3D)
1540	IF (printlev >= 5) THEN
1541	WRITE(numout,*)' interpweight_4D input data dimensions:'
1542	WRITE(numout,*)' iml, jml, lml:',iml, jml, lml
1543	END IF
1544
1545	CASE (4)
1546	ALLOCATE(indims(4), STAT=ALLOC_ERR)
1547	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4D','Problem in allocation of indims','','')
1548	indims = interpweight_get_var4dims_file(filename, varname)
1549	lml = indims(3)
1550	CALL bcast(iml)
1551	CALL bcast(jml)
1552	CALL bcast(lml)
1553	! Taken from `slowproc_update'
1554	IF (initime /= -1) THEN
1555	inNdims = 3
1556	IF (printlev >= 3) WRITE(numout,*) &
1557	'interpweight_4D: taking into account the initial reading time:',initime
1558	IF (initime <= 0) tml = 0
1559
1560	CALL bcast(tml)
1561	ALLOCATE(invar3D(iml,jml,lml), STAT=ALLOC_ERR)
1562	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4D',"Problem in allocation of variable '" // &
1563	TRIM(varname) ,'','')
1564	IF ( initime > 0 ) THEN
1565	IF (is_root_prc) THEN
1566	IF (initime <= tml) THEN
1567	CALL flinget(fid, TRIM(varname), iml, jml, lml, 1, initime, 1, invar3D)
1568	ELSE
1569	CALL flinget(fid, TRIM(varname), iml, jml, lml, 1, tml, 1, invar3D)
1570	ENDIF
1571	ENDIF
1572	ELSE
1573	IF (is_root_prc) THEN
1574	CALL flinget(fid, TRIM(varname), iml, jml, lml, 1, 1, 1, invar3D)
1575	ENDIF
1576	ENDIF
1577	ELSE
1578	CALL bcast(tml)
1579	ALLOCATE(invar4D(iml,jml,lml,tml), STAT=ALLOC_ERR)
1580	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4D',"Problem in allocation of variable '" // &
1581	TRIM(varname) ,'','')
1582	IF (is_root_prc) CALL flinget(fid, TRIM(varname), iml, jml, lml, tml, 1, tml, invar4D)
1583	END IF
1584
1585	IF (printlev >= 5) THEN
1586	WRITE(numout,*)' interpweight_4D input data dimensions'
1587	WRITE(numout,*)' iml, jml, lml, tml:',iml, jml, lml, tml
1588	END IF
1589
1590	CASE DEFAULT
1591	WRITE (S1,'(I3)') inNdims
1592	CALL ipslerr_p(3,'interpweight_4D',"Input number of dimensions: '" // TRIM(S1) // "' not ready !!", &
1593	'','')
1594
1595	END SELECT InputDims
1596	CALL bcast(invar4D)
1597
1598	! Getting longitudes/latitudes from input file
1599	inNLldims = interpweight_get_varNdims_file(filename, inlonname)
1600	IF (inNLldims == 1) THEN
1601	! Allocate memory for latitudes
1602	ALLOCATE(lon_in(iml), STAT=ALLOC_ERR)
1603	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4D','Problem in allocation of variable lon_in','','')
1604
1605	ALLOCATE(lat_in(jml), STAT=ALLOC_ERR)
1606	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4D','Problem in allocation of variable lat_in','','')
1607
1608	IF (is_root_prc) THEN
1609	CALL flinget(fid, TRIM(inlonname), iml, 0, 0, 0, 1, 1, lon_in)
1610	CALL flinget(fid, TRIM(inlatname), jml, 0, 0, 0, 1, 1, lat_in)
1611	ENDIF
1612	CALL bcast(lon_in)
1613	CALL bcast(lat_in)
1614
1615	ALLOCATE(lon(iml,jml), STAT=ALLOC_ERR)
1616	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4D','Problem in allocation of variable lon','','')
1617	ALLOCATE(lat(iml,jml), STAT=ALLOC_ERR)
1618	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4D','Problem in allocation of variable lat','','')
1619
1620	DO ip=1,iml
1621	lat(ip,:) = lat_in(:)
1622	ENDDO
1623	DO jp=1,jml
1624	lon(:,jp) = lon_in(:)
1625	ENDDO
1626	ELSE
1627	! Allocate memory for longitude
1628	ALLOCATE(lon(iml,jml), STAT=ALLOC_ERR)
1629	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4D','Pb in allocation for lon','','')
1630	! Allocate memory for latitudes
1631	ALLOCATE(lat(iml,jml), STAT=ALLOC_ERR)
1632	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4D','Pb in allocation for lat','','')
1633
1634	IF (is_root_prc) THEN
1635	CALL flinget(fid, TRIM(inlonname), iml, jml, 0, 0, 1, 1, lon)
1636	CALL flinget(fid, TRIM(inlatname), iml, jml, 0, 0, 1, 1, lat)
1637	ENDIF
1638	END IF
1639	CALL bcast(lon)
1640	CALL bcast(lat)
1641
1642	IF (is_root_prc) THEN
1643	ALLOCATE(resol_in(iml,jml,2), STAT=ALLOC_ERR)
1644	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4D','Problem in allocation of variable resol_in','','')
1645	END IF
1646
1647	IF (noneg) THEN
1648	IF (is_root_prc) THEN
1649	SELECT CASE (inNdims)
1650	CASE(2)
1651	CALL interpweight_modifying_input2D(iml, jml, 0, 0, 1, 1, zero, invar2D)
1652	CASE(3)
1653	CALL interpweight_modifying_input3D(iml, jml, lml, 0, 1, 1, zero, invar3D)
1654	CASE(4)
1655	CALL interpweight_modifying_input4D(iml, jml, lml, tml, 1, tml, zero, invar4D)
1656	END SELECT
1657	END IF
1658	END IF
1659
1660	IF (printlev >= 5) WRITE(numout,*)' interpweight_4D modifying input ended!'
1661
1662	ALLOCATE(mask(iml,jml), STAT=ALLOC_ERR)
1663	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4D','Problem in allocation of variable mask','','')
1664
1665	IF (TRIM(masktype) == 'var') THEN
1666	ALLOCATE(maskvar(iml,jml), STAT=ALLOC_ERR)
1667	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4D','Problem in allocation of variable maskvar','','')
1668	maskvar(:,:) = 0
1669	! Reads the mask from the file
1670	IF (is_root_prc) THEN
1671	! ATTTENTION!! One can not use flinget in this way. It always spect to have a REAL variable!
1672	CALL flinget(fid, TRIM(maskvarname), iml, jml, 0, 0, 1, 1, maskvar)
1673	WHERE (maskvar > zero)
1674	mask = un
1675	END WHERE
1676	END IF
1677	CALL bcast(mask)
1678	ELSE
1679	SELECT CASE (inNdims)
1680	CASE(2)
1681	CALL interpweight_masking_input2D(iml, jml, 0, 0, 0, 0, masktype, maskvalues, invar2D, &
1682	mask)
1683	CASE(3)
1684	CALL interpweight_masking_input3D(iml, jml, 0, 0, lml, 0, masktype, maskvalues, invar3D, &
1685	mask)
1686	CASE(4)
1687	CALL interpweight_masking_input4D(iml, jml, 0, 0, lml, tml, masktype, maskvalues, invar4D, &
1688	mask)
1689	END SELECT
1690	END IF
1691	IF (printlev >= 5) WRITE(numout,*)' interpweight_4D masking input ended!'
1692
1693	IF (is_root_prc) CALL flinclo(fid)
1694
1695	! Computation of nbvmax
1696	! In 'condveg_soilalb, slowproc_update' nbvmax=200
1697	!
1698	! The number of maximum vegetation map points in the GCM grid is estimated.
1699	! Some lmargin is taken.
1700	!
1701	IF (is_root_prc) THEN
1702	IF ( (maxresollon /= -un) .AND. (maxresollat /= -un) ) THEN
1703	! Assuming regular lon/lat projection for the input data but it is not caluclated assuming a maximum &
1704	! lon,lat resolution given as input paramater
1705	nix=INT(MAXVAL(resolution(:,1))/maxresollon)+2
1706	njx=INT(MAXVAL(resolution(:,2))/maxresollat)+2
1707	ELSE
1708	! Assuming regular lon/lat projection for the input data
1709	CALL interpweight_calc_resolution_in(lon, lat, iml, jml, mincos, R_Earth, pi, resol_in)
1710	nix=INT(MAXVAL(resolution(:,1))/MAXVAL(resol_in(:,:,1)))+2
1711	njx=INT(MAXVAL(resolution(:,2))/MAXVAL(resol_in(:,:,2)))+2
1712	END IF
1713	nbvmax = nix*njx
1714	IF (nbvmax < 200) nbvmax = 200
1715	END IF
1716	CALL bcast(nbvmax)
1717	CALL bcast(nix)
1718	CALL bcast(njx)
1719
1720	callsign = TRIM(varname) // ' map'
1721
1722	ok_interpol = .FALSE.
1723	DO WHILE ( .NOT. ok_interpol )
1724	IF (printlev>=3) WRITE(numout,*) "Projection arrays for ",callsign," : "
1725	IF (printlev>=3) WRITE(numout,*) "nbvmax = ",nbvmax, nix, njx
1726
1727	ALLOCATE(sub_index(nbpt, nbvmax, 2), STAT=ALLOC_ERR)
1728	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4D','Problem in allocation of variable sub_index','','')
1729
1730	sub_index(:,:,:)=zero
1731
1732	ALLOCATE(sub_area(nbpt, nbvmax), STAT=ALLOC_ERR)
1733	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4D','Problem in allocation of variable sub_area','','')
1734
1735	sub_area(:,:)=zero
1736	CALL aggregate_p(nbpt, lalo, neighbours, resolution, contfrac, &
1737	& iml, jml, lon, lat, mask, callsign, &
1738	& nbvmax, sub_index, sub_area, ok_interpol)
1739
1740	IF ( .NOT. ok_interpol ) THEN
1741	DEALLOCATE(sub_area)
1742	DEALLOCATE(sub_index)
1743	nbvmax = nbvmax * 2
1744	ENDIF
1745	ENDDO
1746
1747	! Getting variables thresholds
1748	ALLOCATE(variableusetypes(Nvariabletypes), STAT=ALLOC_ERR)
1749	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4D','Problem in allocation of variable variableusetypesr','','')
1750
1751	IF (variabletypes(1) == -un) THEN
1752	variableusetypes = (/ (ip*1., ip=1,Nvariabletypes) /)
1753	IF (printlev>=4) WRITE(numout,*)' interpweight_4D: Note: Using default ',&
1754	'equivalence between dimension in the file and dimension in the variable, for file ', filename
1755	varmin = 1._r_std
1756	varmax = Nvariabletypes*1._r_std
1757	ELSE
1758	variableusetypes = variabletypes
1759	END IF
1760
1761	outvar4D = zero
1762
1763	! Doing the interpolation as function of the type fo the values
1764	CALL interpweight_provide_fractions4D(typefrac, nbpt, Nvariabletypes, variableusetypes, &
1765	iml, jml, lml, tml, nbvmax, zero, invar4D, sub_area, sub_index, varmax, varmin, &
1766	deux, lalo, outvar4D, aoutvar)
1767	IF (printlev >= 5) WRITE(numout,*)' interpweight_4D end of interpweight_provide_fractions4D'
1768
1769	IF (printlev>=2) THEN
1770	IF (COUNT(aoutvar == -1) >= AINT(nbpt*0.4)) THEN
1771	WRITE(numout,*) 'interpweight_4D total number of points: ', nbpt
1772	WRITE(numout,*) ' Interpolated: ', COUNT(aoutvar /= -1), &
1773	' Non-interpolated: ', COUNT(aoutvar == -1)
1774	msg = "Less than 60% of points '" // TRIM(filename) // "' variable '" // TRIM(varname) // &
1775	"' have been succesfully interpolated !!"
1776	WRITE(numout,*)' ' // TRIM(msg)
1777	ELSE
1778	msg = TRIM(filename) // "' variable: '" // TRIM(varname) // "' interpolation Done _______"
1779	WRITE(numout,*) "interpweight_4D: '" // TRIM(msg)
1780	WRITE(numout,*) ' Total number of points: ', nbpt
1781	WRITE(numout,*) ' Interpweight_4D interpolated: ', COUNT(aoutvar /= -1), &
1782	' Non-interpolated: ', COUNT(aoutvar == -1)
1783	END IF
1784	END IF
1785
1786	! Scaling to 1 the quality variable
1787	DO ip=1, nbpt
1788	IF (aoutvar(ip) /= -1) aoutvar(ip) = aoutvar(ip)/(resolution(ip,1)*resolution(ip,2))/contfrac(ip)
1789	END DO
1790
1791	IF (ALLOCATED(lat_in)) DEALLOCATE(lat_in)
1792	IF (ALLOCATED(lon_in)) DEALLOCATE(lon_in)
1793	DEALLOCATE(lon)
1794	DEALLOCATE(lat)
1795	IF (ALLOCATED(invar2D)) DEALLOCATE(invar2D)
1796	IF (ALLOCATED(invar3D)) DEALLOCATE(invar3D)
1797	IF (ALLOCATED(invar4D)) DEALLOCATE(invar4D)
1798	DEALLOCATE(mask)
1799	IF (ALLOCATED(maskvar)) DEALLOCATE(maskvar)
1800	DEALLOCATE(sub_area)
1801	DEALLOCATE(sub_index)
1802	IF (ALLOCATED(resol_in)) DEALLOCATE(resol_in)
1803
1804	RETURN
1805
1806	END SUBROUTINE interpweight_4D
1807
1808	!! ================================================================================================================================
1809	!! SUBROUTINE : interpweight_2Dcont
1810	!!
1811	!>\BRIEF ! Interpolate any input file to the grid of the model for a continuos 2D-variable
1812	!!
1813	!! DESCRIPTION : (definitions, functional, design, flags):
1814	!!
1815	!! RECENT CHANGE(S): None
1816	!!
1817	!! MAIN OUTPUT VARIABLE(S): outvar2D, aoutvar
1818	!!
1819	!! REFERENCE(S) : None
1820	!!
1821	!! FLOWCHART : None
1822	!! \n
1823	!_ ================================================================================================================================
1824	SUBROUTINE interpweight_2Dcont(nbpt, dim1, dim2, lalo, resolution, neighbours, &
1825	contfrac, filename, varname, inlonname, inlatname, varmin, varmax, noneg, masktype, &
1826	maskvalues, maskvarname, initime, typefrac, defaultvalue, defaultNOvalue, &
1827	outvar2D, aoutvar)
1828
1829	USE ioipsl
1830	USE constantes_var
1831	USE mod_orchidee_para
1832
1833	IMPLICIT NONE
1834
1835	!! 0. Variables and parameter declaration
1836	!
1837	!! 0.1 Input variables
1838	!
1839	INTEGER(i_std), INTENT(in) :: nbpt !! Number of points for which the data
1840	!! needs to be interpolated
1841	INTEGER(i_std), INTENT(in) :: dim1, dim2 !! 4D complementary dimension sizes
1842	REAL(r_std), DIMENSION(nbpt,2), INTENT(in) :: lalo !! Vector of latitude and longitudes
1843	!! (beware of the order = 1 : latitude,
1844	!! 2 : longitude)
1845	REAL(r_std), DIMENSION(nbpt,2), INTENT(in) :: resolution !! The size in km of each grid-box in X and Y
1846	INTEGER(i_std), DIMENSION(nbpt,NbNeighb), INTENT(in) :: neighbours !! Vector of neighbours for each grid point
1847	!! 1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW)
1848
1849	REAL(r_std), DIMENSION(nbpt), INTENT(in) :: contfrac !! Fraction of land in each grid box.
1850	CHARACTER(LEN=80), INTENT(in) :: filename !! name of the input map to read
1851
1852	CHARACTER(LEN=80), INTENT(in) :: varname !! name of the variable to interpolate
1853	CHARACTER(LEN=80), INTENT(in) :: inlonname, inlatname !! names of the longitude and latitude
1854	REAL(r_std), INTENT(in) :: varmin, varmax !! min/max values to use for the
1855	!! renormalization
1856	REAL(r_std), INTENT(in) :: defaultvalue !! default interpolated value
1857	REAL(r_std), INTENT(in) :: defaultNOvalue !! default interpolated NO-value
1858
1859	! Transformations to apply to the original read values from the file
1860	LOGICAL, INTENT(in) :: noneg !! no negative values
1861	CHARACTER(LEN=50), INTENT(in) :: masktype !! Type of masking
1862	!! 'nomask': no-mask is applied
1863	!! 'mbelow': take values below maskvalues(1)
1864	!! 'mabove': take values above maskvalues(1)
1865	!! 'msumrange': take values within 2 ranges;
1866	!! maskvalues(2) <= SUM(vals(k)) <= maskvalues(1)
1867	!! maskvalues(1) < SUM(vals(k)) <= maskvalues(3)
1868	!! (normalized by maskedvalues(3))
1869	!! 'var': mask values are taken from a
1870	!! variable (1/0)
1871	REAL(r_std), DIMENSION(3), INTENT(in) :: maskvalues !! values to use to mask (according to
1872	!! `masktype')
1873	CHARACTER(LEN=250), INTENT(in) :: maskvarname !! name of the variable to use to mask
1874	!! (according to `masktype')
1875	CHARACTER(LEN=50), INTENT(in) :: typefrac !! Type of fraction retrieval:
1876	!! 'default': standard
1877
1878	INTEGER(i_std), INTENT(in) :: initime !! Initial time to take (for input data
1879	!! with time values, -1: not used)
1880	!! 0.2 Modified variables
1881	!
1882	!! 0.3 Output variables
1883	!
1884	REAL(r_std), DIMENSION(nbpt), INTENT(out) :: outvar2D !! 2D output variable and re-dimensioned
1885	REAL(r_std), DIMENSION(nbpt), INTENT(out) :: aoutvar !! availability of input data to
1886	!! interpolate output variable
1887	!! (on the nbpt space)
1888	!
1889	!! 0.4 Local variables
1890	!
1891	INTEGER(i_std) :: iml, jml, lml, tml, fid, ib, ip, jp, it, jj, jv
1892	INTEGER :: inNdims !! Number of dimensions of the input data
1893	INTEGER :: inNLldims !! Number of dimensions of the lon/lat input data
1894	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lat_in, lon_in !! latitude and
1895	!! longitude, extract from input map
1896	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lat, lon !! en 2D ???
1897	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: sub_area !! the area of the fine grid in the
1898	!! model grid ???
1899	!! cf src_global/interpol_help.f90,
1900	!! line 377, called "areaoverlap"
1901	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:,:) :: sub_index !! the indexes from the grid boxes from the
1902	!! data that go into the
1903	!! model's boxes
1904	!! cf src_global/interpol_help.f90,
1905	!! line 300, called "ip"
1906
1907	! REAL(r_std), ALLOCATABLE, DIMENSION(:) :: invar1D !! 1D input values
1908	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: invar2D !! 2D input values
1909	REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:) :: invar3D !! 3D input values
1910	REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:,:) :: invar4D !! 4D input values
1911
1912	REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:) :: resol_in
1913	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: mask
1914	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: maskvar
1915	INTEGER(i_std) :: nix, njx
1916	INTEGER(i_std) :: idi, nbvmax !! nbvmax : number of maximum vegetation map
1917	!! points in the GCM grid ; idi : its counter
1918	CHARACTER(LEN=30) :: callsign !! Allows to specify which variable is beeing
1919	!! treated
1920	LOGICAL :: ok_interpol !! optionnal return of aggregate_2d
1921	INTEGER :: ALLOC_ERR
1922
1923	CHARACTER(LEN=50) :: S1
1924	CHARACTER(LEN=250) :: msg
1925
1926
1927	! Debug vars
1928	INTEGER :: nb_coord, nb_var, nb_gat
1929
1930	!_ ================================================================================================================================
1931	IF (printlev >= 3) WRITE(numout,*)"In interpweight_2Dcont filename: ",TRIM(filename), &
1932	" varname:'" // TRIM(varname) // "'"
1933	! Allocating input data
1934	IF (is_root_prc) THEN
1935	IF (printlev >=5) THEN
1936	WRITE(numout,*) "Entering 'interpweight_2Dcont'. Debug mode."
1937	WRITE(numout,*)'(/," --> fliodmpf")'
1938	CALL fliodmpf (TRIM(filename))
1939	WRITE(numout,*)'(/," --> flioopfd")'
1940	ENDIF
1941	CALL flioopfd (TRIM(filename),fid,nb_dim=nb_coord,nb_var=nb_var,nb_gat=nb_gat)
1942	IF (printlev >=5) THEN
1943	WRITE (numout,'(" Number of coordinate in the file : ",I2)') nb_coord
1944	WRITE (numout,'(" Number of variables in the file : ",I2)') nb_var
1945	WRITE (numout,'(" Number of global attributes in the file : ",I2)') nb_gat
1946	ENDIF
1947	ENDIF
1948
1949	! Getting shape of input variable from input file
1950	inNdims = interpweight_get_varNdims_file(filename, varname)
1951	IF (is_root_prc) CALL flininfo(filename, iml, jml, lml, tml, fid)
1952
1953	InputDims: SELECT CASE (inNdims)
1954	CASE (2)
1955	CALL bcast(iml)
1956	CALL bcast(jml)
1957
1958	ALLOCATE(invar2D(iml,jml), STAT=ALLOC_ERR)
1959	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2Dcont',"Problem in allocation of variable '" // &
1960	TRIM(varname) ,'','')
1961	IF (is_root_prc) CALL flinget(fid, TRIM(varname), iml, jml, 0, 0, 1, 1, invar2D)
1962	IF (printlev >= 5) THEN
1963	WRITE(numout,*)' interpweight_2Dcont input data dimensions _______'
1964	WRITE(numout,*)' iml, jml:',iml, jml
1965	END IF
1966
1967	CASE (3)
1968	CALL bcast(iml)
1969	CALL bcast(jml)
1970	CALL bcast(lml)
1971
1972	ALLOCATE(invar3D(iml,jml, lml), STAT=ALLOC_ERR)
1973	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2Dcont',"Problem in allocation of variable '" // &
1974	TRIM(varname) ,'','')
1975	IF (is_root_prc) CALL flinget(fid, TRIM(varname), iml, jml, lml, 0, 1, 1, invar3D)
1976	IF (printlev >= 5) THEN
1977	WRITE(numout,*)' interpweight_2Dcont input data dimensions:'
1978	WRITE(numout,*)' iml, jml, lml:',iml, jml, lml
1979	END IF
1980
1981	CASE (4)
1982	CALL bcast(iml)
1983	CALL bcast(jml)
1984	CALL bcast(lml)
1985	! Taken from `slowproc_update'
1986	IF (initime /= -1) THEN
1987	inNdims = 3
1988	IF (printlev >= 3) WRITE(numout,*) &
1989	'interpweight_2Dcont: taking into account the initial reading time:',initime
1990	IF (initime <= 0) tml = 0
1991
1992	CALL bcast(tml)
1993	ALLOCATE(invar3D(iml,jml,lml), STAT=ALLOC_ERR)
1994	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2Dcont',"Problem in allocation of variable '" // &
1995	TRIM(varname) ,'','')
1996	IF ( initime > 0 ) THEN
1997	IF (is_root_prc) THEN
1998	IF (initime <= tml) THEN
1999	CALL flinget(fid, TRIM(varname), iml, jml, lml, 1, initime, 1, invar3D)
2000	ELSE
2001	CALL flinget(fid, TRIM(varname), iml, jml, lml, 1, tml, 1, invar3D)
2002	ENDIF
2003	ENDIF
2004	ELSE
2005	IF (is_root_prc) THEN
2006	CALL flinget(fid, TRIM(varname), iml, jml, lml, 1, 1, 1, invar3D)
2007	ENDIF
2008	ENDIF
2009	ELSE
2010	CALL bcast(tml)
2011	ALLOCATE(invar4D(iml,jml,lml,tml), STAT=ALLOC_ERR)
2012	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2Dcont',"Problem in allocation of variable '" // &
2013	TRIM(varname) ,'','')
2014	IF (is_root_prc) CALL flinget(fid, TRIM(varname), iml, jml, lml, tml, 1, tml, invar4D)
2015	END IF
2016
2017	IF (printlev >= 5) THEN
2018	WRITE(numout,*)' interpweight_2Dcont input data dimensions:'
2019	WRITE(numout,*)' iml, jml, lml, tml :',iml, jml, lml, tml
2020	END IF
2021
2022	CASE DEFAULT
2023	WRITE (S1,'(I3)') inNdims
2024	CALL ipslerr_p(3,'interpweight_2Dcont',"Input number of dimensions: '" // TRIM(S1) // "' not ready !!", &
2025	'','')
2026
2027	END SELECT InputDims
2028	CALL bcast(invar2D)
2029
2030	IF (printlev >= 5) WRITE(numout,*)' interpweight_2Dcont getting variable ended!'
2031
2032	! Getting longitudes/latitudes from input file
2033	inNLldims = interpweight_get_varNdims_file(filename, inlonname)
2034
2035	IF (inNLldims == 1) THEN
2036	! Allocate memory for latitudes
2037	ALLOCATE(lon_in(iml), STAT=ALLOC_ERR)
2038	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2Dcont','Problem in allocation of variable lon_in','','')
2039
2040	ALLOCATE(lat_in(jml), STAT=ALLOC_ERR)
2041	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2Dcont','Problem in allocation of variable lat_in','','')
2042
2043	IF (is_root_prc) THEN
2044	CALL flinget(fid, TRIM(inlonname), iml, 0, 0, 0, 1, 1, lon_in)
2045	CALL flinget(fid, TRIM(inlatname), jml, 0, 0, 0, 1, 1, lat_in)
2046	ENDIF
2047	CALL bcast(lon_in)
2048	CALL bcast(lat_in)
2049
2050	ALLOCATE(lon(iml,jml), STAT=ALLOC_ERR)
2051	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2Dcont','Problem in allocation of variable lon','','')
2052	ALLOCATE(lat(iml,jml), STAT=ALLOC_ERR)
2053	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2Dcont','Problem in allocation of variable lat','','')
2054
2055	DO ip=1,iml
2056	lat(ip,:) = lat_in(:)
2057	ENDDO
2058	DO jp=1,jml
2059	lon(:,jp) = lon_in(:)
2060	ENDDO
2061	ELSE
2062	! Allocate memory for longitude
2063	ALLOCATE(lon(iml,jml), STAT=ALLOC_ERR)
2064	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2Dcont','Pb in allocation for lon','','')
2065	! Allocate memory for latitudes
2066	ALLOCATE(lat(iml,jml), STAT=ALLOC_ERR)
2067	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2Dcont','Pb in allocation for lat','','')
2068
2069	IF (is_root_prc) THEN
2070	CALL flinget(fid, TRIM(inlonname), iml, jml, 0, 0, 1, 1, lon)
2071	CALL flinget(fid, TRIM(inlatname), iml, jml, 0, 0, 1, 1, lat)
2072	ENDIF
2073	END IF
2074	CALL bcast(lon)
2075	CALL bcast(lat)
2076
2077	IF (is_root_prc) THEN
2078	ALLOCATE(resol_in(iml,jml,2), STAT=ALLOC_ERR)
2079	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2Dcont','Problem in allocation of variable resol_in','','')
2080	END IF
2081
2082	IF (noneg) THEN
2083	IF (is_root_prc) THEN
2084	SELECT CASE (inNdims)
2085	CASE(2)
2086	CALL interpweight_modifying_input2D(iml, jml, 0, 0, 1, 1, zero, invar2D)
2087	CASE(3)
2088	CALL interpweight_modifying_input3D(iml, jml, lml, 0, 1, 1, zero, invar3D)
2089	CASE(4)
2090	CALL interpweight_modifying_input4D(iml, jml, lml, tml, 1, tml, zero, invar4D)
2091	END SELECT
2092	END IF
2093	END IF
2094	IF (printlev >= 5) WRITE(numout,*)' interpweight_2Dcont modifying input ended!'
2095
2096	!
2097	! Consider all points a priori
2098	!
2099	ALLOCATE(mask(iml,jml), STAT=ALLOC_ERR)
2100	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2Dcont','Problem in allocation of variable mask','','')
2101	mask(:,:) = 0
2102
2103	IF (TRIM(masktype) == 'var') THEN
2104	ALLOCATE(maskvar(iml,jml), STAT=ALLOC_ERR)
2105	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2Dcont','Problem in allocation of variable maskvar','','')
2106	maskvar(:,:) = 0
2107	! Reads the mask from the file
2108	IF (is_root_prc) THEN
2109	! ATTTENTION!! One can not use flinget in this way. It always spect to have a REAL variable!
2110	CALL flinget(fid, TRIM(maskvarname), iml, jml, 0, 0, 1, 1, maskvar)
2111	WHERE (maskvar > zero)
2112	mask = un
2113	END WHERE
2114	END IF
2115	CALL bcast(mask)
2116	ELSE
2117	SELECT CASE (inNdims)
2118	CASE(2)
2119	CALL interpweight_masking_input2D(iml, jml, 0, 0, 0, 0, masktype, maskvalues, invar2D, &
2120	mask)
2121	CASE(3)
2122	CALL interpweight_masking_input3D(iml, jml, 0, 0, lml, 0, masktype, maskvalues, invar3D, &
2123	mask)
2124	CASE(4)
2125	CALL interpweight_masking_input4D(iml, jml, 0, 0, lml, tml, masktype, maskvalues, invar4D, &
2126	mask)
2127	END SELECT
2128	END IF
2129	IF (printlev >= 5) WRITE(numout,*)' interpweight_2Dcont masking input ended!'
2130	IF (is_root_prc) CALL flinclo(fid)
2131
2132	IF (printlev>=5) WRITE(numout,*)" interpweight_2Dcont '" // TRIM(typefrac) // "' interpolation"
2133
2134	! Assuming regular lon/lat projection for the input data
2135	IF (is_root_prc) THEN
2136	CALL interpweight_calc_resolution_in(lon, lat, iml, jml, mincos, R_Earth, pi, resol_in)
2137	END IF
2138
2139	! Computation of nbvmax
2140	! In 'condveg_soilalb, slowproc_update' nbvmax=200
2141
2142	!
2143	! The number of maximum vegetation map points in the GCM grid is estimated.
2144	! Some lmargin is taken.
2145	!
2146	IF (is_root_prc) THEN
2147	nix=INT(MAXVAL(resolution(:,1))/MAXVAL(resol_in(:,:,1)))+2
2148	njx=INT(MAXVAL(resolution(:,2))/MAXVAL(resol_in(:,:,2)))+2
2149	nbvmax = nix*njx
2150	ENDIF
2151	CALL bcast(nbvmax)
2152	CALL bcast(nix)
2153	CALL bcast(njx)
2154
2155	callsign = TRIM(varname) // ' map'
2156
2157	ok_interpol = .FALSE.
2158	DO WHILE ( .NOT. ok_interpol )
2159	IF (printlev>=3) WRITE(numout,*) "Projection arrays for ",callsign," : "
2160	IF (printlev>=3) WRITE(numout,*) "nbvmax = ",nbvmax, nix, njx
2161
2162	ALLOCATE(sub_index(nbpt, nbvmax, 2), STAT=ALLOC_ERR)
2163	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2Dcont','Problem in allocation of variable sub_index','','')
2164
2165	sub_index(:,:,:)=0
2166
2167	ALLOCATE(sub_area(nbpt, nbvmax), STAT=ALLOC_ERR)
2168	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_2Dcont','Problem in allocation of variable sub_area','','')
2169
2170	sub_area(:,:)=zero
2171	IF (printlev>=3) WRITE(numout,*) 'Input data range LON:', MINVAL(lon), MAXVAL(lon)
2172	IF (printlev>=3) WRITE(numout,*) 'Input data range LAT:', MINVAL(lat), MAXVAL(lat)
2173
2174	CALL aggregate_p(nbpt, lalo, neighbours, resolution, contfrac, &
2175	& iml, jml, lon, lat, mask, callsign, &
2176	& nbvmax, sub_index, sub_area, ok_interpol)
2177
2178	IF ( .NOT. ok_interpol ) THEN
2179	DEALLOCATE(sub_area)
2180	DEALLOCATE(sub_index)
2181	nbvmax = nbvmax * 2
2182	ENDIF
2183	ENDDO
2184	IF (printlev>=4) WRITE (numout,*) ' interpweight_2Dcont: read/allocate OK'
2185
2186	outvar2D = zero
2187	! Providing the interpolated grid point
2188
2189	CALL interpweight_provide_interpolation2D(typefrac, nbpt, 0, 0, iml, jml, lml, tml, &
2190	nbvmax, zero, invar2D, sub_area, sub_index, varmax, varmin, un, deux, lalo, &
2191	defaultvalue, defaultNOvalue, outvar2D, aoutvar)
2192	IF (printlev >= 5) WRITE(numout,*)' interpweight_2Dcont end of interpweight_provide_interpolation2D'
2193
2194	IF (printlev>=2) THEN
2195	IF (COUNT(aoutvar == -1) >= AINT(nbpt*0.4)) THEN
2196	WRITE(numout,*) 'interpweight_2Dcont total number of points: ', nbpt
2197	WRITE(numout,*) ' Interpolated: ', COUNT(aoutvar /= -1), &
2198	' non-interpolated: ', COUNT(aoutvar == -1)
2199	msg = "Less than 60% of points '" // TRIM(filename) // "' variable '" // TRIM(varname) // &
2200	"' have been succesfully interpolated !!"
2201	WRITE(numout,*) ' ' // TRIM(msg)
2202	ELSE
2203	msg = TRIM(filename) // "' variable: '" // TRIM(varname) // "' interpolation Done _______"
2204	WRITE(numout,*) "interpweight_2Dcont: '" // TRIM(msg)
2205	WRITE(numout,*) ' Total number of points: ', nbpt
2206	WRITE(numout,*) ' Interpolated: ', COUNT(aoutvar /= -1), &
2207	' Non-interpolated: ', COUNT(aoutvar == -1)
2208	END IF
2209	END IF
2210
2211	! Scaling to 1 the quality variable
2212	DO ip=1, nbpt
2213	IF (aoutvar(ip) /= -1) aoutvar(ip) = aoutvar(ip)/(resolution(ip,1)*resolution(ip,2))/contfrac(ip)
2214	END DO
2215
2216	IF (printlev>=5) WRITE(numout,*) ' interpweight_2Dcont: Interpolation Done'
2217
2218	IF (ALLOCATED(lat_in)) DEALLOCATE(lat_in)
2219	IF (ALLOCATED(lon_in)) DEALLOCATE(lon_in)
2220	DEALLOCATE(lon)
2221	DEALLOCATE(lat)
2222	IF (ALLOCATED(invar2D)) DEALLOCATE(invar2D)
2223	IF (ALLOCATED(invar3D)) DEALLOCATE(invar3D)
2224	IF (ALLOCATED(invar4D)) DEALLOCATE(invar4D)
2225	DEALLOCATE(mask)
2226	IF (ALLOCATED(maskvar)) DEALLOCATE(maskvar)
2227	IF (ALLOCATED(invar4D)) DEALLOCATE(sub_area)
2228	IF (ALLOCATED(invar4D)) DEALLOCATE(sub_index)
2229	IF (ALLOCATED(resol_in)) DEALLOCATE(resol_in)
2230
2231	RETURN
2232
2233	END SUBROUTINE interpweight_2Dcont
2234
2235	!! ================================================================================================================================
2236	!! SUBROUTINE : interpweight_4Dcont
2237	!!
2238	!>\BRIEF ! Interpolate any input file to the grid of the model for a continuos 4D-variable
2239	!!
2240	!! DESCRIPTION : (definitions, functional, design, flags):
2241	!!
2242	!! RECENT CHANGE(S): None
2243	!!
2244	!! MAIN OUTPUT VARIABLE(S): outvar4D, aoutvar
2245	!!
2246	!! REFERENCE(S) : None
2247	!!
2248	!! FLOWCHART : None
2249	!! \n
2250	!_ ================================================================================================================================
2251	SUBROUTINE interpweight_4Dcont(nbpt, dim1, dim2, lalo, resolution, neighbours, &
2252	contfrac, filename, varname, inlonname, inlatname, varmin, varmax, noneg, masktype, &
2253	maskvalues, maskvarname, initime, typefrac, defaultvalue, defaultNOvalue, &
2254	outvar4D, aoutvar)
2255
2256	USE ioipsl
2257	USE constantes_var
2258	USE mod_orchidee_para
2259
2260	IMPLICIT NONE
2261
2262	!! 0. Variables and parameter declaration
2263	!
2264	! 0.1 INPUT
2265	!
2266	INTEGER(i_std), INTENT(in) :: nbpt !! Number of points for which the data
2267	!! needs to be interpolated
2268	INTEGER(i_std), INTENT(in) :: dim1, dim2 !! 4D complementary dimension sizes
2269	REAL(r_std), DIMENSION(nbpt,2), INTENT(in) :: lalo !! Vector of latitude and longitudes
2270	!! (beware of the order = 1 : latitude,
2271	!! 2 : longitude)
2272	REAL(r_std), DIMENSION(nbpt,2), INTENT(in) :: resolution !! The size in km of each grid-box in X and Y
2273	INTEGER(i_std), DIMENSION(nbpt,NbNeighb), INTENT(in) :: neighbours !! Vector of neighbours for each grid point
2274	!! 1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW)
2275
2276	REAL(r_std), DIMENSION(nbpt), INTENT(in) :: contfrac !! Fraction of land in each grid box.
2277	CHARACTER(LEN=80), INTENT(in) :: filename !! name of the input map to read
2278
2279	CHARACTER(LEN=80), INTENT(in) :: varname !! name of the variable to interpolate
2280	CHARACTER(LEN=80), INTENT(in) :: inlonname, inlatname !! names of the longitude and latitude
2281	REAL(r_std), INTENT(in) :: varmin, varmax !! min/max values to use for the
2282	!! renormalization
2283	REAL(r_std), DIMENSION(dim1,dim2), INTENT(in) :: defaultvalue !! default interpolated value
2284	REAL(r_std), INTENT(in) :: defaultNOvalue !! default interpolated NO-value
2285
2286	! Transformations to apply to the original read values from the file
2287	LOGICAL, INTENT(in) :: noneg !! no negative values
2288	CHARACTER(LEN=50), INTENT(in) :: masktype !! Type of masking
2289	!! 'nomask': no-mask is applied
2290	!! 'mbelow': take values below maskvalues(1)
2291	!! 'mabove': take values above maskvalues(1)
2292	!! 'msumrange': take values within 2 ranges;
2293	!! maskvalues(2) <= SUM(vals(k)) <= maskvalues(1)
2294	!! maskvalues(1) < SUM(vals(k)) <= maskvalues(3)
2295	!! (normalized by maskedvalues(3))
2296	!! 'var': mask values are taken from a
2297	!! variable (1/0)
2298	REAL(r_std), DIMENSION(3), INTENT(in) :: maskvalues !! values to use to mask (according to
2299	!! `masktype')
2300	CHARACTER(LEN=250), INTENT(in) :: maskvarname !! name of the variable to use to mask
2301	!! (according to `masktype')
2302	CHARACTER(LEN=50), INTENT(in) :: typefrac !! Type of fraction retrieval:
2303	!! 'default': standard
2304
2305	INTEGER(i_std), INTENT(in) :: initime !! Initial time to take (for input data
2306	!! with time values, -1: not used)
2307	!! 0.2 Modified variables
2308	!
2309	!! 0.3 Output variables
2310	!
2311	REAL(r_std), DIMENSION(nbpt,dim1,dim2), INTENT(out) :: outvar4D !! 2D output variable and re-dimensioned
2312	REAL(r_std), DIMENSION(nbpt), INTENT(out) :: aoutvar !! availability of input data to
2313	!! interpolate output variable
2314	!! (on the nbpt space)
2315	!
2316	!! 0.4 Local variables
2317	!
2318	INTEGER(i_std) :: iml, jml, lml, tml, fid, ib, ip, jp, it, jj, jv
2319	INTEGER :: inNdims !! Number of dimensions of the input data
2320	INTEGER :: inNLldims !! Number of dimensions of the lon/lat input data
2321	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lat_in, lon_in !! latitude and
2322	!! longitude, extract from input map
2323	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lat, lon !! en 2D ???
2324	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: sub_area !! the area of the fine grid in the
2325	!! model grid ???
2326	!! cf src_global/interpol_help.f90,
2327	!! line 377, called "areaoverlap"
2328	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:,:) :: sub_index !! the indexes from the grid boxes from the
2329	!! data that go into the
2330	!! model's boxes
2331	!! cf src_global/interpol_help.f90,
2332	!! line 300, called "ip"
2333
2334	! REAL(r_std), ALLOCATABLE, DIMENSION(:) :: invar1D !! 1D input values
2335	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: invar2D !! 2D input values
2336	REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:) :: invar3D !! 3D input values
2337	REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:,:) :: invar4D !! 4D input values
2338
2339	REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:) :: resol_in
2340	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: maskvar
2341	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: mask
2342	INTEGER(i_std) :: nix, njx
2343	INTEGER(i_std) :: idi, nbvmax !! nbvmax : number of maximum vegetation map
2344	!! points in the GCM grid ; idi : its counter
2345	CHARACTER(LEN=30) :: callsign !! Allows to specify which variable is beeing
2346	!! treated
2347	LOGICAL :: ok_interpol !! optionnal return of aggregate_2d
2348	INTEGER :: ALLOC_ERR
2349
2350	CHARACTER(LEN=50) :: S1
2351	INTEGER, ALLOCATABLE, DIMENSION(:) :: indims
2352	CHARACTER(LEN=250) :: msg
2353
2354	! Debug vars
2355	INTEGER :: nb_coord, nb_var, nb_gat
2356	INTEGER :: nblL
2357	INTEGER, DIMENSION(2) :: lLpt
2358
2359	!_ ================================================================================================================================
2360	IF (printlev >= 3) WRITE(numout, *)"In interpweight_4Dcont filename: '",TRIM(filename), &
2361	" varname:'" // TRIM(varname) // "'"
2362	! Allocating input data
2363	IF (is_root_prc) THEN
2364	IF (printlev >=5) THEN
2365	WRITE(numout,*) "Entering 'interpweight_4Dcont'. Debug mode."
2366	WRITE(numout,*)'(/," --> fliodmpf")'
2367	CALL fliodmpf (TRIM(filename))
2368	WRITE(numout,*)'(/," --> flioopfd")'
2369	ENDIF
2370	CALL flioopfd (TRIM(filename),fid,nb_dim=nb_coord,nb_var=nb_var,nb_gat=nb_gat)
2371	IF (printlev >=5) THEN
2372	WRITE (numout,'(" Number of coordinate in the file : ",I2)') nb_coord
2373	WRITE (numout,'(" Number of variables in the file : ",I2)') nb_var
2374	WRITE (numout,'(" Number of global attributes in the file : ",I2)') nb_gat
2375	ENDIF
2376	ENDIF
2377
2378	! Getting shape of input variable from input file
2379	inNdims = interpweight_get_varNdims_file(filename, varname)
2380	IF (is_root_prc) CALL flininfo(filename, iml, jml, lml, tml, fid)
2381
2382	InputDims: SELECT CASE (inNdims)
2383	CASE (2)
2384	CALL bcast(iml)
2385	CALL bcast(jml)
2386
2387	ALLOCATE(invar2D(iml,jml), STAT=ALLOC_ERR)
2388	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4Dcont',"Problem in allocation of variable '" // &
2389	TRIM(varname) ,'','')
2390	IF (is_root_prc) CALL flinget(fid, TRIM(varname), iml, jml, 0, 0, 1, 1, invar2D)
2391	IF (printlev >= 5) THEN
2392	WRITE(numout,*)' interpweight_4Dcont input data dimensions:'
2393	WRITE(numout,*)' iml, jml:',iml, jml
2394	END IF
2395
2396	CASE (3)
2397	ALLOCATE(indims(3),STAT=ALLOC_ERR)
2398	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4Dcont','Problem in allocation of variable indims','','')
2399	indims = interpweight_get_var3dims_file(filename, varname)
2400	lml = indims(3)
2401	CALL bcast(iml)
2402	CALL bcast(jml)
2403	CALL bcast(lml)
2404
2405	ALLOCATE(invar3D(iml,jml, lml), STAT=ALLOC_ERR)
2406	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4Dcont',"Problem in allocation of variable '" // &
2407	TRIM(varname) ,'','')
2408	IF (is_root_prc) CALL flinget(fid, TRIM(varname), iml, jml, lml, 0, 1, 1, invar3D)
2409	IF (printlev >= 5) THEN
2410	WRITE(numout,*)' interpweight_4Dcont input data dimensions:'
2411	WRITE(numout,*)' iml, jml, lml:',iml, jml, lml
2412	END IF
2413
2414	CASE (4)
2415	ALLOCATE(indims(4), STAT=ALLOC_ERR)
2416	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4Dcont','Problem in allocation of indims','','')
2417	indims = interpweight_get_var4dims_file(filename, varname)
2418	lml = indims(3)
2419	CALL bcast(iml)
2420	CALL bcast(jml)
2421	CALL bcast(lml)
2422	! Taken from `slowproc_update'
2423	IF (initime /= -1) THEN
2424	inNdims = 3
2425	IF (printlev >= 3) WRITE(numout,*) &
2426	'interpweight_4Dcont: taking into account the initial reading time:',initime
2427	IF (initime <= 0) tml = 0
2428
2429	CALL bcast(tml)
2430	ALLOCATE(invar3D(iml,jml,lml), STAT=ALLOC_ERR)
2431	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4Dcont',"Problem in allocation of variable '" // &
2432	TRIM(varname) ,'','')
2433	IF ( initime > 0 ) THEN
2434	IF (is_root_prc) THEN
2435	IF (initime <= tml) THEN
2436	CALL flinget(fid, TRIM(varname), iml, jml, lml, 1, initime, 1, invar3D)
2437	ELSE
2438	CALL flinget(fid, TRIM(varname), iml, jml, lml, 1, tml, 1, invar3D)
2439	ENDIF
2440	ENDIF
2441	ELSE
2442	IF (is_root_prc) THEN
2443	CALL flinget(fid, TRIM(varname), iml, jml, lml, 1, 1, 1, invar3D)
2444	ENDIF
2445	ENDIF
2446	ELSE
2447	CALL bcast(tml)
2448	ALLOCATE(invar4D(iml,jml,lml,tml), STAT=ALLOC_ERR)
2449	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4Dcont',"Problem in allocation of variable '" // &
2450	TRIM(varname) ,'','')
2451	IF (is_root_prc) CALL flinget(fid, TRIM(varname), iml, jml, lml, tml, 1, tml, invar4D)
2452	END IF
2453
2454	IF (printlev >= 5) THEN
2455	WRITE(numout,*)'interpweight_4Dcont input data dimensions:'
2456	WRITE(numout,*)' iml, jml, lml, tml:',iml, jml, lml, tml
2457	END IF
2458
2459	CASE DEFAULT
2460	WRITE (S1,'(I3)') inNdims
2461	CALL ipslerr_p(3,'interpweight_4Dcont',"Input number of dimensions: '" // S1 // "' not ready !!",'','')
2462
2463	END SELECT InputDims
2464	CALL bcast(invar4D)
2465
2466	IF (printlev >= 5) WRITE(numout,*)' interpweight_4Dcont getting variable ended!'
2467
2468	! Getting longitudes/latitudes from input file
2469	inNLldims = interpweight_get_varNdims_file(filename, inlonname)
2470
2471	IF (inNLldims == 1) THEN
2472	! Allocate memory for latitudes
2473	ALLOCATE(lon_in(iml), STAT=ALLOC_ERR)
2474	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4Dcont','Problem in allocation of variable lon_in','','')
2475
2476	ALLOCATE(lat_in(jml), STAT=ALLOC_ERR)
2477	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4Dcont','Problem in allocation of variable lat_in','','')
2478
2479	IF (is_root_prc) THEN
2480	CALL flinget(fid, TRIM(inlonname), iml, 0, 0, 0, 1, 1, lon_in)
2481	CALL flinget(fid, TRIM(inlatname), jml, 0, 0, 0, 1, 1, lat_in)
2482	ENDIF
2483	CALL bcast(lon_in)
2484	CALL bcast(lat_in)
2485
2486	ALLOCATE(lon(iml,jml), STAT=ALLOC_ERR)
2487	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4Dcont','Problem in allocation of variable lon','','')
2488	ALLOCATE(lat(iml,jml), STAT=ALLOC_ERR)
2489	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4Dcont','Problem in allocation of variable lat','','')
2490
2491	DO ip=1,iml
2492	lat(ip,:) = lat_in(:)
2493	ENDDO
2494	DO jp=1,jml
2495	lon(:,jp) = lon_in(:)
2496	ENDDO
2497	ELSE
2498	! Allocate memory for longitude
2499	ALLOCATE(lon(iml,jml), STAT=ALLOC_ERR)
2500	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4Dcont','Pb in allocation for lon','','')
2501	! Allocate memory for latitudes
2502	ALLOCATE(lat(iml,jml), STAT=ALLOC_ERR)
2503	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4Dcont','Pb in allocation for lat','','')
2504
2505	IF (is_root_prc) THEN
2506	CALL flinget(fid, TRIM(inlonname), iml, jml, 0, 0, 1, 1, lon)
2507	CALL flinget(fid, TRIM(inlatname), iml, jml, 0, 0, 1, 1, lat)
2508	ENDIF
2509	END IF
2510	CALL bcast(lon)
2511	CALL bcast(lat)
2512
2513	IF (is_root_prc) THEN
2514	ALLOCATE(resol_in(iml,jml,2), STAT=ALLOC_ERR)
2515	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4Dcont','Problem in allocation of variable resol_in','','')
2516	END IF
2517
2518	IF (noneg) THEN
2519	IF (is_root_prc) THEN
2520	SELECT CASE (inNdims)
2521	CASE(2)
2522	CALL interpweight_modifying_input2D(iml, jml, 0, 0, 0, 0, zero, invar2D)
2523	CASE(3)
2524	CALL interpweight_modifying_input3D(iml, jml, lml, 0, 0, 0, zero, invar3D)
2525	CASE(4)
2526	CALL interpweight_modifying_input4D(iml, jml, lml, tml, 1, tml, zero, invar4D)
2527	END SELECT
2528	END IF
2529	END IF
2530	CALL bcast(invar4D)
2531	IF (printlev >= 5) WRITE(numout,*)' interpweight_4Dcont modifying input ended!'
2532
2533	!
2534	! Consider all points a priori
2535	!
2536	ALLOCATE(mask(iml,jml), STAT=ALLOC_ERR)
2537	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4Dcont','Problem in allocation of variable mask','','')
2538	mask(:,:) = zero
2539
2540	IF (TRIM(masktype) == 'var') THEN
2541	ALLOCATE(maskvar(iml,jml), STAT=ALLOC_ERR)
2542	maskvar(:,:) = zero
2543	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4Dcont','Problem in allocation of variable maskvar','','')
2544	! Reads the mask from the file
2545	IF (is_root_prc) THEN
2546	! ATTTENTION!! One can not use flinget in this way. It always spect to have a REAL variable!
2547	CALL flinget(fid, TRIM(maskvarname), iml, jml, 0, 0, 1, 1, maskvar)
2548	WHERE (maskvar > zero)
2549	mask = un
2550	ENDWHERE
2551	END IF
2552	CALL bcast(mask)
2553	ELSE
2554	SELECT CASE (inNdims)
2555	CASE(2)
2556	CALL interpweight_masking_input2D(iml, jml, 0, 0, 0, 0, masktype, maskvalues, invar2D, &
2557	mask)
2558	CASE(3)
2559	CALL interpweight_masking_input3D(iml, jml, 0, 0, lml, 0, masktype, maskvalues, invar3D, &
2560	mask)
2561	CASE(4)
2562	CALL interpweight_masking_input4D(iml, jml, 0, 0, lml, tml, masktype, maskvalues, invar4D, &
2563	mask)
2564	END SELECT
2565	END IF
2566
2567	IF (is_root_prc) CALL flinclo(fid)
2568
2569	! Assuming regular lon/lat projection for the input data
2570	IF (is_root_prc) THEN
2571	CALL interpweight_calc_resolution_in(lon, lat, iml, jml, mincos, R_Earth, pi, resol_in)
2572	END IF
2573
2574	! Computation of nbvmax
2575	! In 'condveg_soilalb, slowproc_update' nbvmax=200
2576
2577	!
2578	! The number of maximum vegetation map points in the GCM grid is estimated.
2579	! Some lmargin is taken.
2580	!
2581	IF (is_root_prc) THEN
2582	nix=INT(MAXVAL(resolution(:,1))/MAXVAL(resol_in(:,:,1)))+2
2583	njx=INT(MAXVAL(resolution(:,2))/MAXVAL(resol_in(:,:,2)))+2
2584	nbvmax = nix*njx
2585	ENDIF
2586	CALL bcast(nbvmax)
2587	CALL bcast(nix)
2588	CALL bcast(njx)
2589
2590	!
2591	callsign = TRIM(varname) // ' map'
2592	!
2593	ok_interpol = .FALSE.
2594	DO WHILE ( .NOT. ok_interpol )
2595	IF (printlev>=3) WRITE(numout,*) "Projection arrays for ",callsign," : "
2596	IF (printlev>=3) WRITE(numout,*) "nbvmax = ",nbvmax, nix, njx
2597
2598
2599	ALLOCATE(sub_index(nbpt, nbvmax, 2), STAT=ALLOC_ERR)
2600	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4Dcont','Problem in allocation of variable sub_index','','')
2601
2602	sub_index(:,:,:)=0
2603
2604	ALLOCATE(sub_area(nbpt, nbvmax), STAT=ALLOC_ERR)
2605	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_4Dcont','Problem in allocation of variable sub_area','','')
2606
2607	sub_area(:,:)=zero
2608	IF (printlev>=3) WRITE(numout,*) 'Input data range LON:', MINVAL(lon), MAXVAL(lon)
2609	IF (printlev>=3) WRITE(numout,*) 'Input data range LAT:', MINVAL(lat), MAXVAL(lat)
2610
2611	CALL aggregate_p(nbpt, lalo, neighbours, resolution, contfrac, &
2612	& iml, jml, lon, lat, mask, callsign, &
2613	& nbvmax, sub_index, sub_area, ok_interpol)
2614
2615	IF ( .NOT. ok_interpol ) THEN
2616	DEALLOCATE(sub_area)
2617	DEALLOCATE(sub_index)
2618	nbvmax = nbvmax * 2
2619	ENDIF
2620	ENDDO
2621	IF (printlev>=4) WRITE (numout,*) ' interpweight_4Dcont: read/allocate OK'
2622
2623	outvar4D = zero
2624
2625	! Providing the interpolated grid point
2626	CALL interpweight_provide_interpolation4D(typefrac, nbpt, dim1, dim2, iml, jml, lml, tml, &
2627	nbvmax, zero, invar4D, sub_area, sub_index, varmax, varmin, un, deux, lalo, &
2628	defaultvalue, defaultNOvalue, outvar4D, aoutvar)
2629	IF (printlev >= 5) WRITE(numout,*)' interpweight_4Dcont end of interpweight_provide_interpolation4D'
2630
2631	IF (printlev>=2) THEN
2632	IF (COUNT(aoutvar == -1) >= AINT(nbpt*0.4)) THEN
2633	WRITE(numout,*) 'interpweight_4Dcont total number of points: ', nbpt
2634	WRITE(numout,*) ' Interpolated: ', COUNT(aoutvar /= -1), &
2635	' non-interpolated: ', COUNT(aoutvar == -1)
2636	msg = "Less than 60% of points '" // TRIM(filename) // "' variable '" // TRIM(varname) // &
2637	"' have been succesfully interpolated !!"
2638	WRITE(numout,*) ' ' // TRIM(msg)
2639	ELSE
2640	msg = TRIM(filename) // "' variable: '" // TRIM(varname) // "' interpolation Done _______"
2641	WRITE(numout,*) "interpweight_4Dcont: '" // TRIM(msg)
2642	WRITE(numout,*) ' Total number of points: ', nbpt
2643	WRITE(numout,*) ' Interpolated: ', COUNT(aoutvar /= -1), &
2644	' Non-interpolated: ', COUNT(aoutvar == -1)
2645	END IF
2646	END IF
2647
2648	! Scaling to 1 the quality variable
2649	DO ip=1, nbpt
2650	IF (aoutvar(ip) /= -1) aoutvar(ip) = aoutvar(ip)/(resolution(ip,1)*resolution(ip,2))/contfrac(ip)
2651	END DO
2652
2653	IF (printlev>=5) WRITE(numout,*) ' interpweight_4Dcont: Interpolation Done'
2654
2655	IF (ALLOCATED(lat_in)) DEALLOCATE(lat_in)
2656	IF (ALLOCATED(lon_in)) DEALLOCATE(lon_in)
2657	DEALLOCATE(lon)
2658	DEALLOCATE(lat)
2659	IF (ALLOCATED(invar2D)) DEALLOCATE(invar2D)
2660	IF (ALLOCATED(invar3D)) DEALLOCATE(invar3D)
2661	IF (ALLOCATED(invar4D)) DEALLOCATE(invar4D)
2662	DEALLOCATE(mask)
2663	IF (ALLOCATED(maskvar)) DEALLOCATE(maskvar)
2664	DEALLOCATE(sub_area)
2665	DEALLOCATE(sub_index)
2666	IF (ALLOCATED(resol_in)) DEALLOCATE(resol_in)
2667
2668	RETURN
2669
2670	END SUBROUTINE interpweight_4Dcont
2671
2672	!! ================================================================================================================================
2673	!! SUBROUTINE : interpweight_calc_resolution_in
2674	!!
2675	!>\BRIEF ! Subroutine to compute the resolution of the input data
2676	!!
2677	!! DESCRIPTION : (definitions, functional, design, flags):
2678	!!
2679	!! RECENT CHANGE(S): None
2680	!!
2681	!! MAIN OUTPUT VARIABLE(S): resolin
2682	!!
2683	!! REFERENCE(S) : None
2684	!!
2685	!! FLOWCHART : None
2686	!! \n
2687	!_ ================================================================================================================================
2688	SUBROUTINE interpweight_calc_resolution_in(lons, lats, dx, dy, mcos, REarth, piv, resolin)
2689
2690	USE constantes_var
2691
2692	IMPLICIT NONE
2693
2694	!! 0. Variables and parameter declaration
2695	!
2696	!! 0.1 Input variables
2697	!
2698	INTEGER(i_std), INTENT(in) :: dx, dy !! Resolution input data
2699	REAL(r_std), INTENT(in) :: mcos !! cosinus
2700	REAL(r_std), INTENT(in) :: REarth !! Earth Radius
2701	REAL(r_std), INTENT(in) :: piv !! pi
2702	REAL(r_std), DIMENSION(dx,dy), INTENT(in) :: lons, lats !! longitude and latitudes (degrees)
2703	!! 0.2 Modified variables
2704	!
2705	!! 0.3 Output variables
2706	!
2707	REAL(r_std), DIMENSION(dx,dy,2), INTENT(out) :: resolin !! distance between grid points (km)
2708	!
2709	!! 0.4 Local variables
2710	!
2711	INTEGER :: ip,jp
2712	REAL(r_std) :: coslat
2713
2714	DO ip=1,dx
2715	DO jp=1,dy
2716	!
2717	! Resolution in longitude
2718	!
2719	coslat = MAX( COS( lats(ip,jp) * piv/180. ), mcos )
2720	IF ( ip .EQ. 1 ) THEN
2721	resolin(ip,jp,1) = ABS( lons(ip+1,jp) - lons(ip,jp) ) * piv/180. * REarth * coslat
2722	ELSEIF ( ip .EQ. dx ) THEN
2723	resolin(ip,jp,1) = ABS( lons(ip,jp) - lons(ip-1,jp) ) * piv/180. * REarth * coslat
2724	ELSE
2725	resolin(ip,jp,1) = ABS( lons(ip+1,jp) - lons(ip-1,jp) )/2. * piv/180. * REarth * coslat
2726	ENDIF
2727	!
2728	! Resolution in latitude
2729	!
2730	IF ( jp .EQ. 1 ) THEN
2731	resolin(ip,jp,2) = ABS( lats(ip,jp) - lats(ip,jp+1) ) * piv/180. * REarth
2732	ELSEIF ( jp .EQ. dy ) THEN
2733	resolin(ip,jp,2) = ABS( lats(ip,jp-1) - lats(ip,jp) ) * piv/180. * REarth
2734	ELSE
2735	resolin(ip,jp,2) = ABS( lats(ip,jp-1) - lats(ip,jp+1) )/2. * piv/180. * REarth
2736	ENDIF
2737	!
2738	ENDDO
2739	ENDDO
2740
2741	END SUBROUTINE interpweight_calc_resolution_in
2742
2743	!! ================================================================================================================================
2744	!! SUBROUTINE : interpweight_modifying_input1D
2745	!!
2746	!>\BRIEF ! modify the initial 1D values from the given file
2747	!!
2748	!! DESCRIPTION : (definitions, functional, design, flags):
2749	!!
2750	!! RECENT CHANGE(S): None
2751	!!
2752	!! MAIN OUTPUT VARIABLE(S): ivar
2753	!!
2754	!! REFERENCE(S) : None
2755	!!
2756	!! FLOWCHART : None
2757	!! \n
2758	!_ ================================================================================================================================
2759	SUBROUTINE interpweight_modifying_input1D(dim1, dim2, dim3, dim4, stime, etime, zeroval, ivar)
2760
2761	USE constantes_var
2762	USE ioipsl_para
2763
2764	IMPLICIT NONE
2765
2766	!! 0. Variables and parameter declaration
2767	!! 0.1 Input variables
2768	INTEGER, INTENT(in) :: dim1, dim2, dim3, dim4 !! size of the dimensions to get
2769	INTEGER, INTENT(in) :: stime, etime !! starting/ending time step to get
2770	REAL(r_std), INTENT(in) :: zeroval !! zero value
2771	!! 0.2 Modified variables
2772	REAL(r_std), DIMENSION(dim1), INTENT(inout) :: ivar !! modified input data
2773
2774
2775	WHERE (ivar(:) < zeroval )
2776	ivar(:) = zeroval
2777	END WHERE
2778
2779	END SUBROUTINE interpweight_modifying_input1D
2780
2781	!! ================================================================================================================================
2782	!! SUBROUTINE : interpweight_modifying_input2D
2783	!!
2784	!>\BRIEF ! modify the initial 2D values from the given file
2785	!!
2786	!! DESCRIPTION : (definitions, functional, design, flags):
2787	!!
2788	!! RECENT CHANGE(S): None
2789	!!
2790	!! MAIN OUTPUT VARIABLE(S): ivar
2791	!!
2792	!! REFERENCE(S) : None
2793	!!
2794	!! FLOWCHART : None
2795	!! \n
2796	!_ ================================================================================================================================
2797	SUBROUTINE interpweight_modifying_input2D(dim1, dim2, dim3, dim4, stime, etime, zeroval, ivar)
2798
2799	USE constantes_var
2800	USE ioipsl_para
2801
2802	IMPLICIT NONE
2803
2804	!! 0. Variables and parameter declaration
2805	!! 0.1 Input variables
2806	INTEGER, INTENT(in) :: dim1, dim2, & !! size of the dimensions to get
2807	dim3, dim4
2808	INTEGER, INTENT(in) :: stime, etime !! starting/ending time step to get
2809	REAL(r_std), INTENT(in) :: zeroval !! zero value
2810	!! 0.2 Modified variables
2811	REAL(r_std), DIMENSION(dim1,dim2), INTENT(inout) :: ivar !! modified input data
2812
2813	WHERE (ivar(:,:) < zeroval )
2814	ivar(:,:) = zeroval
2815	END WHERE
2816
2817	END SUBROUTINE interpweight_modifying_input2D
2818
2819	!! ================================================================================================================================
2820	!! SUBROUTINE : interpweight_modifying_input3D
2821	!!
2822	!>\BRIEF ! modify the initial 3D values from the given file
2823	!!
2824	!! DESCRIPTION : (definitions, functional, design, flags):
2825	!!
2826	!! RECENT CHANGE(S): None
2827	!!
2828	!! MAIN OUTPUT VARIABLE(S): ivar
2829	!!
2830	!! REFERENCE(S) : None
2831	!!
2832	!! FLOWCHART : None
2833	!! \n
2834	!_ ================================================================================================================================
2835	SUBROUTINE interpweight_modifying_input3D(dim1, dim2, dim3, dim4, stime, etime, zeroval, ivar)
2836
2837	USE constantes_var
2838	USE ioipsl_para
2839
2840	IMPLICIT NONE
2841
2842	!! 0. Variables and parameter declaration
2843	!! 0.1 Input variables
2844	INTEGER, INTENT(in) :: dim1, dim2, & !! size of the dimensions to get
2845	dim3, dim4
2846	INTEGER, INTENT(in) :: stime, etime !! starting/ending time step to get
2847	REAL(r_std), INTENT(in) :: zeroval !! zero value
2848	!! 0.2 Modified variables
2849	REAL(r_std), DIMENSION(dim1,dim2,dim3), INTENT(inout):: ivar !! modified input data
2850
2851	WHERE (ivar(:,:,:) < zeroval )
2852	ivar(:,:,:) = zeroval
2853	END WHERE
2854
2855	END SUBROUTINE interpweight_modifying_input3D
2856
2857	!! ================================================================================================================================
2858	!! SUBROUTINE : interpweight_modifying_input4D
2859	!!
2860	!>\BRIEF ! modify the initial 4D values from the given file
2861	!!
2862	!! DESCRIPTION : (definitions, functional, design, flags):
2863	!!
2864	!! RECENT CHANGE(S): None
2865	!!
2866	!! MAIN OUTPUT VARIABLE(S): ivar
2867	!!
2868	!! REFERENCE(S) : None
2869	!!
2870	!! FLOWCHART : None
2871	!! \n
2872	!_ ================================================================================================================================
2873	SUBROUTINE interpweight_modifying_input4D(dim1, dim2, dim3, dim4, stime, etime, zeroval, ivar)
2874
2875	USE constantes_var
2876	USE ioipsl_para
2877
2878	IMPLICIT NONE
2879
2880	!! 0. Variables and parameter declaration
2881	!! 0.1 Input variables
2882	INTEGER, INTENT(in) :: dim1, dim2, & !! size of the dimensions to get
2883	dim3, dim4
2884	INTEGER, INTENT(in) :: stime, etime !! starting/ending time step to get
2885	REAL(r_std), INTENT(in) :: zeroval !! zero value
2886	!! 0.2 Modified variables
2887	REAL(r_std), DIMENSION(dim1,dim2,dim3,dim4), INTENT(inout) :: ivar !! modified input data
2888
2889	WHERE (ivar(:,:,:,:) < zeroval )
2890	ivar(:,:,:,:) = zeroval
2891	END WHERE
2892
2893	END SUBROUTINE interpweight_modifying_input4D
2894
2895	!! ================================================================================================================================
2896	!! SUBROUTINE : interpweight_interpweight_1D
2897	!!
2898	!>\BRIEF ! mask 1D input values
2899	!!
2900	!! DESCRIPTION : (definitions, functional, design, flags):
2901	!!
2902	!! RECENT CHANGE(S): None
2903	!!
2904	!! MAIN OUTPUT VARIABLE(S): msk
2905	!!
2906	!! REFERENCE(S) : None
2907	!!
2908	!! FLOWCHART : None
2909	!! \n
2910	!_ ================================================================================================================================
2911	SUBROUTINE interpweight_masking_input1D(dx, dy, d1, d2, d3, d4, mtype, mvalues, var, &
2912	msk)
2913
2914	USE constantes_var
2915
2916	IMPLICIT NONE
2917
2918	!! 0. Variables and parameter declaration
2919	!! 0.1 Input variables
2920	INTEGER, INTENT(in) :: dx, dy, d1, & !! shapes of the input variables
2921	d2, d3, d4
2922	CHARACTER(LEN=50), INTENT(in) :: mtype !! Type of masking
2923	!! 'nomask': no-mask is applied
2924	!! 'mbelow': take values below values(1)
2925	!! 'mabove': take values above values(1)
2926	!! 'msumrange': take values within 2 ranges;
2927	!! values(2) <= SUM(vals(k)) <= values(1)
2928	!! values(1) < SUM(vals(k)) <= values(3):
2929	!! renormalize
2930	!! 'var': mask values are taken from a
2931	!! variable (>0)
2932	REAL(r_std), DIMENSION(3), INTENT(in) :: mvalues !! values to use to mask (according to `mtype')
2933	!! 0.2 Modified variables
2934	REAL(r_std), DIMENSION(dx), INTENT(inout) :: var !! variable from which provide the mask
2935	INTEGER(i_std), DIMENSION(dx), INTENT(inout) :: msk !! mask to retrieve
2936	!! 0.3 Output variables
2937	!! 0.4 Local variables
2938	INTEGER :: i,j,k,l
2939
2940	SELECT CASE (mtype)
2941	CASE('nomask')
2942	msk=un
2943	CASE('mbelow')
2944	! e.g.:
2945	! Exclude the points where there is never a LAI value. It is probably
2946	! an ocean point.
2947	!
2948	DO i=1,dx
2949	IF ( var(i) < mvalues(1) ) msk(i) = un
2950	END DO
2951	CASE('mabove')
2952	DO i=1,dx
2953	IF ( var(i) > mvalues(1) ) msk(i) = un
2954	ENDDO
2955	CASE('msumrange')
2956	CALL ipslerr_p(3, 'interpweight_masking_input1D',"'msumrange' no sens on a 1D variable",'','')
2957	CASE DEFAULT
2958	CALL ipslerr_p(3, 'interpweight_masking_input1D',"mask with '" // TRIM(mtype) // "' not ready !!",'','')
2959	END SELECT
2960
2961	END SUBROUTINE interpweight_masking_input1D
2962
2963	!! ================================================================================================================================
2964	!! SUBROUTINE : interpweight_masking_input2D
2965	!!
2966	!>\BRIEF ! mask 2D input values
2967	!!
2968	!! DESCRIPTION : (definitions, functional, design, flags):
2969	!!
2970	!! RECENT CHANGE(S): None
2971	!!
2972	!! MAIN OUTPUT VARIABLE(S): msk
2973	!!
2974	!! REFERENCE(S) : None
2975	!!
2976	!! FLOWCHART : None
2977	!! \n
2978	!_ ================================================================================================================================
2979	SUBROUTINE interpweight_masking_input2D(dx, dy, d1, d2, d3, d4, mtype, mvalues, var, &
2980	msk)
2981
2982	USE constantes_var
2983
2984	IMPLICIT NONE
2985
2986	!! 0. Variables and parameter declaration
2987	!! 0.1 Input variables
2988	INTEGER, INTENT(in) :: dx, dy, d1, & !! shapes of the input variables
2989	d2, d3, d4
2990	CHARACTER(LEN=50), INTENT(in) :: mtype !! Type of masking
2991	!! 'nomask': no-mask is applied
2992	!! 'mbelow': take values below values(1)
2993	!! 'mabove': take values above values(1)
2994	!! 'msumrange': take values within 2 ranges;
2995	!! values(2) <= SUM(vals(k)) <= values(1)
2996	!! values(1) < SUM(vals(k)) <= values(3):
2997	!! renormalize
2998	!! 'var': mask values are taken from a
2999	!! variable (>0)
3000	REAL(r_std), DIMENSION(3), INTENT(in) :: mvalues !! values to use to mask (according to `mtype')
3001	!! 0.2 Modified variables
3002	REAL(r_std), DIMENSION(dx,dy), INTENT(inout) :: var !! variable from which provide the mask
3003	INTEGER(i_std), DIMENSION(dx,dy), INTENT(inout) :: msk !! mask to retrieve
3004	!! 0.3 Output variables
3005	!! 0.4 Local variables
3006	INTEGER :: i,j,k,l
3007	REAL(r_std) :: sumval
3008
3009
3010	SELECT CASE (mtype)
3011	CASE('nomask')
3012	msk=un
3013	CASE('mbelow')
3014	! e.g.:
3015	! Exclude the points where there is never a LAI value. It is probably
3016	! an ocean point.
3017	!
3018	DO i=1,dx
3019	DO j=1,dy
3020	IF ( var(i,j) < mvalues(1) ) msk(i,j) = un
3021	END DO
3022	END DO
3023	CASE('mabove')
3024	DO i=1,dx
3025	DO j=1,dy
3026	IF ( var(i,j) > mvalues(1) ) msk(i,j) = un
3027	ENDDO
3028	ENDDO
3029	CASE('msumrange')
3030	IF (printlev>=3) WRITE(numout,*) &
3031	'interpweight_masking_input2D: masking with mask sum range. Interval:', mvalues
3032	DO i=1,dx
3033	DO j=1,dy
3034	sumval=SUM(var(i,:))
3035	IF ( sumval .GE. mvalues(2) .AND. sumval .LE. mvalues(1) ) THEN
3036	msk(i,j) = un
3037	ELSEIF ( sumval .GT. mvalues(1) .AND. sumval .LE. mvalues(3) ) THEN
3038	! normalization
3039	var(i,j) = var(i,j) / sumval
3040	msk(i,j) = un
3041	ENDIF
3042	ENDDO
3043	ENDDO
3044	CASE DEFAULT
3045	CALL ipslerr_p(3, 'interpweight_masking_input2D',"mask with '" // TRIM(mtype) // "' not ready !!",'','')
3046	END SELECT
3047
3048	END SUBROUTINE interpweight_masking_input2D
3049
3050	!! ================================================================================================================================
3051	!! SUBROUTINE : interpweight_masking_input3D
3052	!!
3053	!>\BRIEF ! mask 3D input values
3054	!!
3055	!! DESCRIPTION : (definitions, functional, design, flags):
3056	!!
3057	!! RECENT CHANGE(S): None
3058	!!
3059	!! MAIN OUTPUT VARIABLE(S): msk
3060	!!
3061	!! REFERENCE(S) : None
3062	!!
3063	!! FLOWCHART : None
3064	!! \n
3065	!_ ================================================================================================================================
3066	SUBROUTINE interpweight_masking_input3D(dx, dy, d1, d2, d3, d4, mtype, mvalues, var, &
3067	msk)
3068
3069	USE constantes_var
3070
3071	IMPLICIT NONE
3072
3073	!! 0. Variables and parameter declaration
3074	!! 0.1 Input variables
3075	INTEGER, INTENT(in) :: dx, dy, d1, & !! shapes of the input variables
3076	d2, d3, d4
3077	CHARACTER(LEN=50), INTENT(in) :: mtype !! Type of masking
3078	!! 'nomask': no-mask is applied
3079	!! 'mbelow': take values below values(1)
3080	!! 'mabove': take values above values(1)
3081	!! 'msumrange': take values within 2 ranges;
3082	!! values(2) <= SUM(vals(k)) <= values(1)
3083	!! values(1) < SUM(vals(k)) <= values(3):
3084	!! renormalize
3085	!! 'var': mask values are taken from a
3086	!! variable (>0)
3087	REAL(r_std), DIMENSION(3), INTENT(in) :: mvalues !! values to use to mask (according to `mtype')
3088	!! 0.2 Modified variables
3089	REAL(r_std), DIMENSION(dx,dy,d3), INTENT(inout) :: var !! variable from which provide the mask
3090	INTEGER(i_std), DIMENSION(dx,dy), INTENT(inout) :: msk !! mask to retrieve
3091	!! 0.3 Output variables
3092	!! 0.4 Local variables
3093	INTEGER :: i,j,k,l
3094	REAL(r_std) :: sumval
3095
3096
3097	SELECT CASE (mtype)
3098	CASE('nomask')
3099	msk=un
3100	CASE('mbelow')
3101	! e.g.:
3102	! Exclude the points where there is never a LAI value. It is probably
3103	! an ocean point.
3104	!
3105	DO i=1,dx
3106	DO j=1,dy
3107	IF ( ANY(var(i,j,:) < mvalues(1)) ) msk(i,j) = un
3108	END DO
3109	END DO
3110	CASE('mabove')
3111	DO i=1,dx
3112	DO j=1,dy
3113	IF ( ANY(var(i,j,:) > mvalues(1)) ) msk(i,j) = un
3114	ENDDO
3115	ENDDO
3116	CASE('msumrange')
3117	DO i=1,dx
3118	DO j=1,dy
3119	sumval=SUM(var(i,j,:))
3120	IF ( sumval .GE. mvalues(2) .AND. sumval .LE. mvalues(1)) THEN
3121	msk(i,j) = un
3122	ELSEIF ( sumval .GT. mvalues(1) .AND. sumval .LE. mvalues(3)) THEN
3123	! normalization
3124	var(i,j,:) = var(i,j,:) / sumval
3125	msk(i,j) = un
3126	ENDIF
3127	ENDDO
3128	ENDDO
3129	CASE DEFAULT
3130	CALL ipslerr_p(3, 'interpweight_masking_input3D',"mask with '" // TRIM(mtype) // "' not ready !!",'','')
3131	END SELECT
3132	END SUBROUTINE interpweight_masking_input3D
3133
3134	!! ================================================================================================================================
3135	!! SUBROUTINE : interpweight_masking_input4D
3136	!!
3137	!>\BRIEF ! mask 4D input values
3138	!!
3139	!! DESCRIPTION : (definitions, functional, design, flags):
3140	!!
3141	!! RECENT CHANGE(S): None
3142	!!
3143	!! MAIN OUTPUT VARIABLE(S): msk
3144	!!
3145	!! REFERENCE(S) : None
3146	!!
3147	!! FLOWCHART : None
3148	!! \n
3149	!_ ================================================================================================================================
3150	SUBROUTINE interpweight_masking_input4D(dx, dy, d1, d2, d3, d4, mtype, mvalues, var, &
3151	msk)
3152
3153	USE constantes_var
3154
3155	IMPLICIT NONE
3156
3157	!! 0. Variables and parameter declaration
3158	!! 0.1 Input variables
3159	INTEGER, INTENT(in) :: dx, dy, d1, & !! shapes of the input variables
3160	d2, d3, d4
3161	CHARACTER(LEN=50), INTENT(in) :: mtype !! Type of masking
3162	!! 'nomask': no-mask is applied
3163	!! 'mbelow': take values below values(1)
3164	!! 'mabove': take values above values(1)
3165	!! 'msumrange': take values within 2 ranges;
3166	!! values(2) <= SUM(vals(k)) <= values(1)
3167	!! values(1) < SUM(vals(k)) <= values(3):
3168	!! renormalize
3169	!! 'var': mask values are taken from a
3170	!! variable (>0)
3171	REAL(r_std), DIMENSION(3), INTENT(in) :: mvalues !! values to use to mask (according to `mtype')
3172	!! 0.2 Modified variables
3173	REAL(r_std), DIMENSION(dx,dy,d3,d4), INTENT(inout) :: var !! variable from which provide the mask
3174	INTEGER(i_std), DIMENSION(dx,dy), INTENT(inout) :: msk !! mask to retrieve
3175	!! 0.3 Output variables
3176	!! 0.4 Local variables
3177	INTEGER :: i,j,k,l
3178	REAL(r_std) :: sumval
3179
3180
3181	SELECT CASE (mtype)
3182	CASE('nomask')
3183	msk=un
3184	CASE('mbelow')
3185	! e.g.:
3186	! Exclude the points where there is never a LAI value. It is probably
3187	! an ocean point.
3188	!
3189	DO i=1,dx
3190	DO j=1,dy
3191	IF ( ANY(var(i,j,:,:) < mvalues(1)) ) msk(i,j) = un
3192	END DO
3193	END DO
3194	CASE('mabove')
3195	DO i=1,dx
3196	DO j=1,dy
3197	IF ( ANY(var(i,j,:,:) > mvalues(1)) ) msk(i,j) = un
3198	ENDDO
3199	ENDDO
3200	CASE('msumrange')
3201	IF (printlev>=3) WRITE(numout,*) &
3202	'interpweight_masking_input4D: masking with mask sum range. Interval:', mvalues
3203	DO i=1,dx
3204	DO j=1,dy
3205	DO l=1,d4
3206	sumval=SUM(var(i,j,:,l))
3207	IF ( sumval .GE. mvalues(2) .AND. sumval .LE. mvalues(1) .AND. msk(i,j) /= un) THEN
3208	msk(i,j) = un
3209	ELSEIF ( sumval .GT. mvalues(1) .AND. sumval .LE. mvalues(3)) THEN
3210	! normalization
3211	var(i,j,:,l) = var(i,j,:,l) / sumval
3212	msk(i,j) = un
3213	ENDIF
3214	ENDDO
3215	ENDDO
3216	ENDDO
3217	CASE DEFAULT
3218	CALL ipslerr_p(3, 'interpweight_masking_input4D',"mask with '" // TRIM(mtype) // "' not ready !!",'','')
3219	END SELECT
3220
3221	END SUBROUTINE interpweight_masking_input4D
3222
3223	!! ================================================================================================================================
3224	!! SUBROUTINE : interpweight_provide_fractions1D
3225	!!
3226	!>\BRIEF ! provide fractions from a 1D incoming variable
3227	!!
3228	!! DESCRIPTION : (definitions, functional, design, flags):
3229	!!
3230	!! RECENT CHANGE(S): None
3231	!!
3232	!! MAIN OUTPUT VARIABLE(S): ovar, aovar
3233	!!
3234	!! REFERENCE(S) : None
3235	!!
3236	!! FLOWCHART : None
3237	!! \n
3238	!_ ================================================================================================================================
3239	SUBROUTINE interpweight_provide_fractions1D(tint, npt, Ntypes, valstypes, &
3240	dx, dy, d3, d4, nbmax, zeroval, ivar1D, sarea, sindex, vmax, vmin, &
3241	two, latlon, ovar, aovar)
3242
3243	USE constantes_var
3244
3245	IMPLICIT NONE
3246
3247	!! 0. Variables and parameter declaration
3248	!! 0.1 Input variables
3249	CHARACTER(LEN=50), INTENT(in) :: tint !! tint: type of interpolation
3250	!! 'XYKindTime': Input values are kinds
3251	!! of something with a temporal
3252	!! evolution on the dx*dy matrix
3253	INTEGER, INTENT(in) :: npt !! total number of grid points
3254	INTEGER, INTENT(in) :: dx,dy,d3,d4 !! shape of the input values
3255	INTEGER, INTENT(in) :: nbmax !! maximum input grid cells to use
3256	!! for in 1 grid point
3257	REAL(r_std), INTENT(in) :: zeroval, two !! zero and 2. real values
3258	INTEGER, INTENT(in) :: Ntypes !! number of types
3259	REAL(r_std), DIMENSION(Ntypes), INTENT(in) :: valstypes !! value for each type
3260	REAL(r_std), DIMENSION(dx), INTENT(in) :: ivar1D !! values of the input data
3261	REAL(r_std), DIMENSION(npt,nbmax), INTENT(in) :: sarea !! overlaping area of a given input
3262	!! grid point
3263	INTEGER(i_std), DIMENSION(npt,nbmax), INTENT(in) :: sindex !! index of the overlaping input
3264	!! grid point
3265	REAL(r_std), DIMENSION(npt,2), INTENT(in) :: latlon !! longitude and latitude of the destiny
3266	!! grid points
3267	!! 0.2 Modified variables
3268	REAL(r_std), INTENT(inout) :: vmax, vmin !! thresholds of the values
3269	!! 0.3 Output variables
3270	REAL(r_std), DIMENSION(npt,Ntypes), INTENT(out) :: ovar !! output variable (on the nbpt space)
3271	REAL(r_std), DIMENSION(npt), INTENT(out) :: aovar !! availability of input data to
3272	!! interpolate output variable (on
3273	!! the nbpt space)
3274	!! >0.: As the total sum of the
3275	!! space fraction used from the
3276	!! source
3277	!! -1: Did not find any required values
3278	!! from input data
3279	!! 0.4 Local variables
3280	INTEGER :: i,j,k,ij,ib,idi,jj,ip,jp,jv,it
3281	CHARACTER(LEN=250) :: msg
3282
3283
3284	! Initialization
3285	!!
3286	ovar = zeroval
3287
3288	SELECT CASE (tint)
3289	CASE ('default')
3290	IF (printlev>=4) WRITE(numout,*)' interpweight_provide_fractions1D Fractions following default1D method'
3291	IF (ALL(ivar1D == zeroval)) THEN
3292	msg = "' wrong 1D input variable"
3293	CALL ipslerr_p(3,'interpweight_provide_fractions1D',"'" // TRIM(tint) // TRIM(msg),'','')
3294	END IF
3295	DO ib=1,npt
3296	aovar(ib) = zeroval
3297	idi = COUNT(sarea(ib,:) > zeroval)
3298	IF ( idi > 0 ) THEN
3299	DO jj=1,idi
3300	ip = sindex(ib,jj)
3301	DO it = 1, Ntypes
3302	IF (ivar1d(ip) == valstypes(it)*un) THEN
3303	ovar(ib,it) = ovar(ib,it) + sarea(ib,jj)
3304	END IF
3305	END DO
3306	aovar(ib) = aovar(ib) + sarea(ib,jj)
3307	ENDDO
3308	!
3309	! Normalize
3310	!
3311	ovar(ib,:) = ovar(ib,:)/aovar(ib)
3312	!
3313	! Quality variable
3314	ELSE
3315	! No points found.
3316	aovar(ib) = -un
3317	ovar(ib,INT((vmax+vmin)/two)) = un
3318	ENDIF
3319	IF (SUM(ovar(ib,:)) > 1.0000001) THEN
3320	IF (printlev>=3) WRITE(numout,*) ' interpweight_provide_fractions1D On point ', ib, &
3321	' location: ', latlon(ib,2), latlon(ib,1),' total fractions=', SUM(ovar(ib,:))
3322	IF (printlev>=3) THEN
3323	WRITE(numout,*) ' type fraction: '
3324	DO i=1, Ntypes
3325	WRITE(numout,*) i, ovar(ib,i)
3326	END DO
3327	END IF
3328	msg='total of fractions above 1.!'
3329	CALL ipslerr_p(3,'interpweight_provide_fractions1D',TRIM(msg),'','')
3330	END IF
3331	END DO
3332
3333	IF (printlev >= 3 .AND. ANY(aovar == -un)) THEN
3334	WRITE(numout,*) 'interpweight_provide_fractions1D: '
3335	WRITE(numout,*) 'Some points on the model grid did not have any points on the source grid for interpolation.'
3336	WRITE(numout,*) 'These points are initialized by putting the average VAR ', INT((vmax+vmin)/two), ' to 1.'
3337	WRITE(numout,*) 'The points are (ib, lon, lat) :'
3338	DO ib=1, npt
3339	IF (aovar(ib) == -un) WRITE(numout,*) ib, latlon(ib,2), latlon(ib,1)
3340	END DO
3341	END IF
3342	CASE DEFAULT
3343	CALL ipslerr_p(3,'interpweight_provide_fractions1D',"Inerpolation type '" // TRIM(tint) // "' not ready!",'','')
3344	END SELECT
3345
3346	END SUBROUTINE interpweight_provide_fractions1D
3347
3348	!! ================================================================================================================================
3349	!! SUBROUTINE : interpweight_provide_fractions2D
3350	!!
3351	!>\BRIEF ! provide fractions from a 2D incoming variable
3352	!!
3353	!! DESCRIPTION : (definitions, functional, design, flags):
3354	!!
3355	!! RECENT CHANGE(S): None
3356	!!
3357	!! MAIN OUTPUT VARIABLE(S): ovar, aovar
3358	!!
3359	!! REFERENCE(S) : None
3360	!!
3361	!! FLOWCHART : None
3362	!! \n
3363	!_ ================================================================================================================================
3364	SUBROUTINE interpweight_provide_fractions2D(tint, npt, Ntypes, valstypes, &
3365	dx, dy, d3, d4, nbmax, zeroval, ivar2D, sarea, sindex, vmax, vmin, &
3366	two, latlon, ovar, aovar)
3367
3368	USE constantes_var
3369
3370	IMPLICIT NONE
3371
3372	!! 0. Variables and parameter declaration
3373	!! 0.1 Input variables
3374	CHARACTER(LEN=50), INTENT(in) :: tint !! tint: type of interpolation
3375	!! 'XYKindTime': Input values are kinds
3376	!! of something with a temporal
3377	!! evolution on the dx*dy matrix
3378	INTEGER, INTENT(in) :: npt !! total number of grid points
3379	INTEGER, INTENT(in) :: dx,dy,d3,d4 !! shape of the input values
3380	INTEGER, INTENT(in) :: nbmax !! maximum input grid cells to use
3381	!! for in 1 grid point
3382	REAL(r_std), INTENT(in) :: zeroval, two !! zero and 2. real values
3383	INTEGER, INTENT(in) :: Ntypes !! number of types
3384	REAL(r_std), DIMENSION(Ntypes), INTENT(in) :: valstypes !! value for each type
3385	REAL(r_std), DIMENSION(dx,dy), INTENT(in) :: ivar2D !! values of the input data
3386	REAL(r_std), DIMENSION(npt,nbmax), INTENT(in) :: sarea !! overlaping area of a given input
3387	!! grid point
3388	INTEGER(i_std), DIMENSION(npt,nbmax,2), INTENT(in) :: sindex !! index of the overlaping input
3389	!! grid point
3390	REAL(r_std), DIMENSION(npt,2), INTENT(in) :: latlon !! longitude and latitude of the destiny
3391	!! grid points
3392	!! 0.2 Modified variables
3393	REAL(r_std), INTENT(inout) :: vmax, vmin !! thresholds of the values
3394	!! 0.3 Output variables
3395	REAL(r_std), DIMENSION(npt,Ntypes), INTENT(out) :: ovar !! output variable (on the nbpt space)
3396	REAL(r_std), DIMENSION(npt), INTENT(out) :: aovar !! availability of input data to
3397	!! interpolate output variable (on
3398	!! the nbpt space)
3399	!! >0.: As the total sum of the
3400	!! space fraction used from the
3401	!! source
3402	!! -1: Found any required values
3403	!! from input data
3404	!! 0.4 Local variables
3405	INTEGER :: i,j,k,ij,ib,idi,jj,ip,jp,jv,it
3406	CHARACTER(LEN=250) :: msg
3407
3408
3409	! Initialization
3410	!!
3411	ovar = zeroval
3412
3413	SELECT CASE (tint)
3414	CASE ('default')
3415	IF (printlev>=3) WRITE(numout,*)' interpweight_provide_fractions2D Fractions following default2D method'
3416	IF (ALL(ivar2D == zeroval)) THEN
3417	msg = "' wrong 2D input variable"
3418	CALL ipslerr_p(3,'interpweight_provide_fractions2D',"'" // TRIM(tint) // TRIM(msg),'','')
3419	END IF
3420	DO ib=1,npt
3421	aovar(ib) = zeroval
3422	idi = COUNT(sarea(ib,:) > zeroval)
3423	IF ( idi > 0 ) THEN
3424	DO jj=1,idi
3425	ip = sindex(ib,jj,1)
3426	jp = sindex(ib,jj,2)
3427	DO it = 1, Ntypes
3428	IF (ivar2d(ip,jp) == valstypes(it)*un) THEN
3429	ovar(ib,it) = ovar(ib,it) + sarea(ib,jj)
3430	END IF
3431	END DO
3432	aovar(ib) = aovar(ib) + sarea(ib,jj)
3433	ENDDO
3434	!
3435	! Normalize
3436	!
3437	ovar(ib,:) = ovar(ib,:)/aovar(ib)
3438	!
3439	! Quality variable
3440	ELSE
3441	! No points on the source grid were found for the current point.
3442	! We provide here a default value.
3443	aovar(ib) = -un
3444	ovar(ib,INT((vmax+vmin)/two)) = un
3445	ENDIF
3446
3447	IF (SUM(ovar(ib,:)) > 1.0000001) THEN
3448	WRITE(numout,*) ' interpweight_provide_fractions2D On point ', ib, ' location: ', latlon(ib,2), &
3449	latlon(ib,1),' total fractions=', SUM(ovar(ib,:))
3450	WRITE(numout,*) ' type fraction _______'
3451	DO i=1, Ntypes
3452	WRITE(numout,*) i, ovar(ib,i)
3453	END DO
3454	msg='total of fractions above 1.!'
3455	CALL ipslerr_p(3,'interpweight_provide_fractions2D',TRIM(msg),'','')
3456	END IF
3457	END DO
3458
3459	IF (printlev >= 3 .AND. ANY(aovar == -un)) THEN
3460	WRITE(numout,*) 'interpweight_provide_fractions2D: '
3461	WRITE(numout,*) 'Some points on the model grid did not have any points on the source grid for interpolation.'
3462	WRITE(numout,*) 'These points are initialized by putting the average VAR ', INT((vmax+vmin)/two), ' to 1.'
3463	WRITE(numout,*) 'The points are (ib, lon, lat) :'
3464	DO ib=1, npt
3465	IF (aovar(ib) == -un) WRITE(numout,*) ib, latlon(ib,2), latlon(ib,1)
3466	END DO
3467	END IF
3468
3469	CASE DEFAULT
3470	CALL ipslerr_p(3,'interpweight_provide_fractions2D',"Inerpolation type '" // TRIM(tint) // "' not ready!",'','')
3471	END SELECT
3472
3473	END SUBROUTINE interpweight_provide_fractions2D
3474
3475	!! ================================================================================================================================
3476	!! SUBROUTINE : interpweight_provide_fractions3D
3477	!!
3478	!>\BRIEF ! provide fractions from a 3D incoming variable
3479	!!
3480	!! DESCRIPTION : (definitions, functional, design, flags):
3481	!!
3482	!! RECENT CHANGE(S): None
3483	!!
3484	!! MAIN OUTPUT VARIABLE(S): ovar, aovar
3485	!!
3486	!! REFERENCE(S) : None
3487	!!
3488	!! FLOWCHART : None
3489	!! \n
3490	!_ ================================================================================================================================
3491	SUBROUTINE interpweight_provide_fractions3D(tint, npt, Ntypes, valstypes, &
3492	dx, dy, d3, d4, nbmax, zeroval, ivar3D, sarea, sindex, vmax, vmin, &
3493	two, latlon, ovar, aovar)
3494
3495	USE constantes_var
3496
3497	IMPLICIT NONE
3498
3499	!! 0. Variables and parameter declaration
3500	!! 0.1 Input variables
3501	CHARACTER(LEN=50), INTENT(in) :: tint !! tint: type of interpolation
3502	!! 'XYKindTime': Input values are kinds
3503	!! of something with a temporal
3504	!! evolution on the dx*dy matrix
3505	INTEGER, INTENT(in) :: npt !! total number of grid points
3506	INTEGER, INTENT(in) :: dx,dy,d3,d4 !! shape of the input values
3507	INTEGER, INTENT(in) :: nbmax !! maximum input grid cells to use
3508	!! for in 1 grid point
3509	REAL(r_std), INTENT(in) :: zeroval, two !! zero and 2. real values
3510	INTEGER, INTENT(in) :: Ntypes !! number of types
3511	REAL(r_std), DIMENSION(Ntypes), INTENT(in) :: valstypes !! value for each type
3512	REAL(r_std), DIMENSION(dx,dy,d3), INTENT(in) :: ivar3D !! values of the input data
3513	REAL(r_std), DIMENSION(npt,nbmax), INTENT(in) :: sarea !! overlaping area of a given input
3514	!! grid point
3515	INTEGER(i_std), DIMENSION(npt,nbmax,2), INTENT(in) :: sindex !! index of the overlaping input
3516	!! grid point
3517	REAL(r_std), DIMENSION(npt,2), INTENT(in) :: latlon !! longitude and latitude of the destiny
3518	!! grid points
3519	!! 0.2 Modified variables
3520	REAL(r_std), DIMENSION(Ntypes), INTENT(inout) :: vmax, vmin !! thresholds of the values
3521	!! 0.3 Output variables
3522	REAL(r_std), DIMENSION(npt,Ntypes), INTENT(out) :: ovar !! output variable (on the nbpt space)
3523	REAL(r_std), DIMENSION(npt), INTENT(out) :: aovar !! availability of input data to
3524	!! interpolate output variable (on
3525	!! the nbpt space)
3526	!! >0.: As the total sum of the
3527	!! space fraction used from the
3528	!! source
3529	!! -1: Found any required values
3530	!! from input data
3531	!! 0.4 Local variables
3532	INTEGER :: i,j,k,ij,ib,idi,jj,ip,jp,jv,it
3533	CHARACTER(LEN=250) :: msg
3534	INTEGER :: Nzeros
3535	INTEGER :: ALLOC_ERR
3536	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: subfrac
3537	REAL(r_std) :: sumpt
3538
3539
3540	! Initialization
3541	!!
3542	ovar = zeroval
3543
3544	IF (Ntypes /= d3) THEN
3545	msg = 'Different number of types than third dimension in the input data!!'
3546	CALL ipslerr_p(3,'interpweight_provide_fractions3D',TRIM(msg),'','')
3547	END IF
3548
3549	Nzeros = 0
3550	SELECT CASE (tint)
3551	CASE ('default')
3552	IF (printlev>=3) WRITE(numout,*)'Fractions following default3D method'
3553	IF (ALL(ivar3D == zeroval)) THEN
3554	msg = "' wrong 3D input variable"
3555	CALL ipslerr_p(3,'interpweight_provide_fractions3D',"'" // TRIM(tint) // TRIM(msg),'','')
3556	END IF
3557	DO ib=1,npt
3558	aovar(ib) = zeroval
3559	idi = COUNT(sarea(ib,:) > zeroval)
3560	IF ( idi > 0 ) THEN
3561	sumpt = zeroval
3562	DO jj=1,idi
3563	ip = sindex(ib,jj,1)
3564	jp = sindex(ib,jj,2)
3565	DO it = 1, Ntypes
3566	! Do not get that areas with zero input value and with missing value
3567	IF (ivar3D(ip,jp,it) > zeroval .AND. ivar3D(ip,jp,it) < two) THEN
3568	ovar(ib,it) = ovar(ib,it) + ivar3D(ip,jp,it)*sarea(ib,jj)
3569	sumpt = sumpt + ivar3D(ip,jp,it)*sarea(ib,jj)
3570	END IF
3571	END DO
3572	aovar(ib) = aovar(ib) + sarea(ib,jj)
3573	ENDDO
3574	!
3575	! Normalize
3576	!
3577	ovar(ib,:) = ovar(ib,:)/aovar(ib)
3578	!
3579	! Quality variable
3580	ELSE
3581	Nzeros = Nzeros + 1
3582	aovar(ib) = -un
3583
3584	DO it = 1, d3
3585	ovar(ib,INT((vmax(it)+vmin(it))/two)) = un
3586	END DO
3587	ENDIF
3588	! Tests
3589	!
3590	IF (SUM(ovar(ib,:)) > 1.0000001) THEN
3591	idi = COUNT(sarea(ib,:) > zeroval)
3592	ALLOCATE(subfrac(idi), STAT=ALLOC_ERR)
3593	IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'interpweight_provide_fractions3D','Problem in allocation of subfrac','','')
3594	IF (printlev>=3) THEN
3595	WRITE(numout,*) 'interpweight_provide_fractions3D On point ', ib, &
3596	' location: ', latlon(ib,2), latlon(ib,1),' total fractions=', SUM(ovar(ib,:)), &
3597	' total points from invar: ', idi, ' total invar area: ', aovar(ib)
3598	WRITE(numout,*) ' type ivar3D jj sarea fraction: '
3599
3600	DO it = 1, Ntypes
3601	DO jj=1,idi
3602	ip = sindex(ib,jj,1)
3603	jp = sindex(ib,jj,2)
3604	IF (ivar3D(ip,jp,it) > zeroval) THEN
3605	WRITE(numout,*) it, ivar3D(ip,jp,it), jj, sarea(ib,jj), ovar(ib,it)
3606	END IF
3607	END DO
3608	END DO
3609	END IF
3610	msg='total of fractions above 1.!'
3611	CALL ipslerr_p(3,'interpweight_provide_fractions3D',TRIM(msg),'','')
3612	END IF
3613
3614	END DO
3615
3616	IF (printlev >= 3 .AND. ANY(aovar == -un)) THEN
3617	WRITE(numout,*) 'interpweight_provide_fractions3D: '
3618	WRITE(numout,*) 'Some points on the model grid did not have any points on the source grid for interpolation.'
3619	WRITE(numout,*) 'These points are initialized by putting the average VAR ', INT((vmax+vmin)/two), ' to 1.'
3620	WRITE(numout,*) 'The points are (ib, lon, lat) :'
3621	DO ib=1, npt
3622	IF (aovar(ib) == -un) WRITE(numout,*) ib, latlon(ib,2), latlon(ib,1)
3623	END DO
3624	END IF
3625
3626	IF (COUNT(aovar == -1) == npt) THEN
3627	msg='No grid points could be interpolated from the source grid for the current process'
3628	CALL ipslerr_p(2,'interpweight_provide_fractions3D',TRIM(msg),'','')
3629	END IF
3630
3631	IF (printlev >=3) WRITE(numout,*) 'interpweight_provide_fractions3D: Nzeros=', &
3632	Nzeros,' aovar=', COUNT(aovar == -1), ' nbpt=', npt
3633
3634	IF (COUNT(aovar == -1) /= Nzeros) THEN
3635	msg='Something went wrong COUNT(aovar == -1) /= Nzeros'
3636	CALL ipslerr_p(3,'interpweight_provide_fractions3D',TRIM(msg),'','')
3637	END IF
3638
3639	CASE DEFAULT
3640	CALL ipslerr_p(3,'interpweight_provide_fractions3D',"Inerpolation type '" // TRIM(tint) // "' not ready!",'','')
3641	END SELECT
3642
3643	END SUBROUTINE interpweight_provide_fractions3D
3644
3645	!! ================================================================================================================================
3646	!! SUBROUTINE : interpweight_provide_fractions4D
3647	!!
3648	!>\BRIEF ! provide fractions from a 4D incoming variable
3649	!!
3650	!! DESCRIPTION : (definitions, functional, design, flags):
3651	!!
3652	!! RECENT CHANGE(S): None
3653	!!
3654	!! MAIN OUTPUT VARIABLE(S): ovar, aovar
3655	!!
3656	!! REFERENCE(S) : None
3657	!!
3658	!! FLOWCHART : None
3659	!! \n
3660	!_ ================================================================================================================================
3661	SUBROUTINE interpweight_provide_fractions4D(tint, npt, Ntypes, valstypes, &
3662	dx, dy, d3, d4, nbmax, zeroval, ivar4D, sarea, sindex, vmax, vmin, &
3663	two, latlon, ovar, aovar)
3664
3665	USE constantes_var
3666
3667	IMPLICIT NONE
3668
3669	!! 0. Variables and parameter declaration
3670	!! 0.1 Input variables
3671	CHARACTER(LEN=50), INTENT(in) :: tint !! tint: type of interpolation
3672	!! 'XYKindTime': Input values are kinds
3673	!! of something with a temporal
3674	!! evolution on the dx*dy matrix
3675	INTEGER, INTENT(in) :: npt !! total number of grid points
3676	INTEGER, INTENT(in) :: dx,dy,d3,d4 !! shape of the input values
3677	INTEGER, INTENT(in) :: nbmax !! maximum input grid cells to use
3678	!! for in 1 grid point
3679	REAL(r_std), INTENT(in) :: zeroval, two !! zero and 2. real values
3680	INTEGER, INTENT(in) :: Ntypes !! number of types
3681	REAL(r_std), DIMENSION(Ntypes), INTENT(in) :: valstypes !! value for each type
3682	REAL(r_std), DIMENSION(dx,dy,d3,d4), INTENT(in) :: ivar4D !! values of the input data
3683	REAL(r_std), DIMENSION(npt,nbmax), INTENT(in) :: sarea !! overlaping area of a given input
3684	!! grid point
3685	INTEGER(i_std), DIMENSION(npt,nbmax,2), INTENT(in) :: sindex !! index of the overlaping input
3686	!! grid point
3687	REAL(r_std), DIMENSION(npt,2), INTENT(in) :: latlon !! longitude and latitude of the destiny
3688	!! grid points
3689	!! 0.2 Modified variables
3690	REAL(r_std), DIMENSION(Ntypes), INTENT(inout) :: vmax, vmin !! thresholds of the values
3691	!! 0.3 Output variables
3692	REAL(r_std), DIMENSION(npt,Ntypes,d4), INTENT(out) :: ovar !! output variable (on the nbpt space)
3693	REAL(r_std), DIMENSION(npt), INTENT(out) :: aovar !! availability of input data to
3694	!! interpolate output variable (on
3695	!! the nbpt space)
3696	!! >0.: As the total sum of the
3697	!! space fraction used from the
3698	!! source
3699	!! -1: Found any required values
3700	!! from input data
3701	!! 0.4 Local variables
3702	INTEGER :: i,j,k,ij,ib,idi,jj,ip,jp,jv,it,itt
3703	CHARACTER(LEN=250) :: msg
3704
3705
3706	! Initialization
3707	!!
3708	ovar = zeroval
3709
3710	IF (Ntypes /= d3) THEN
3711	msg = 'Different number of types than third dimension in the input data!!'
3712	CALL ipslerr_p(3,'interpweight_provide_fractions4D',TRIM(msg),'','')
3713	END IF
3714
3715	SELECT CASE (tint)
3716	CASE ('default')
3717	IF (printlev>=3) WRITE(numout,*)'Fractions following default4D method'
3718	IF (ALL(ivar4D == zeroval)) THEN
3719	msg = "' wrong 4D input variable"
3720	CALL ipslerr_p(3,'interpweight_provide_fractions4D',"'" // TRIM(tint) // TRIM(msg),'','')
3721	END IF
3722	DO ib=1,npt
3723	aovar(ib) = zeroval
3724	idi = COUNT(sarea(ib,:) > zeroval)
3725	IF ( idi > 0 ) THEN
3726	DO jj=1,idi
3727	ip = sindex(ib,jj,1)
3728	jp = sindex(ib,jj,2)
3729	DO it = 1, Ntypes
3730	DO itt = 1, d4
3731	! Do not get that areas with zero input value and with missing value
3732	IF (ivar4D(ip,jp,it,itt) > zeroval .AND. ivar4D(ip,jp,it,itt) < 20.) THEN
3733	ovar(ib,it,itt) = ovar(ib,it,itt) + ivar4D(ip,jp,it,itt)*sarea(ib,jj)
3734	END IF
3735	END DO
3736	END DO
3737	aovar(ib) = aovar(ib) + sarea(ib,jj)
3738	ENDDO
3739	!
3740	! Normalize
3741	!
3742	ovar(ib,:,:) = ovar(ib,:,:)/aovar(ib)
3743	!
3744	! Quality variable
3745	ELSE
3746	! No points found. Set the average of input threashold (vmax, vmin) as value.
3747	aovar(ib) = -un
3748	IF (printlev >= 3) WRITE(numout,*) 'interpweight_provide_fractions4D: ',&
3749	'No points found for interpolating of point ib=',ib,&
3750	' with location lon, lat = ',latlon(ib,2), latlon(ib,1)
3751	DO it=1, Ntypes
3752	DO itt=1, d4
3753	ovar(ib,INT((vmax(it)+vmin(it))/two),itt) = un
3754	END DO
3755	END DO
3756	ENDIF
3757	! DO itt = 1, d4
3758	! IF (SUM(ovar(ib,:,itt)) > 1.0000001) THEN
3759	! IF (printlev>=3) THEN
3760	! WRITE(numout,*) ' interpweight_provide_fractions4D On point ', ib, &
3761	! ' location: ', latlon(ib,2), latlon(ib,1),' itt=', itt, &
3762	! ' total fractions=', SUM(ovar(ib,:,itt))
3763	! WRITE(numout,*) ' type fraction:'
3764	! DO i=1, Ntypes
3765	! WRITE(numout,*) i, ovar(ib,:,itt)
3766	! END DO
3767	! END IF
3768	! msg='total of fractions above 1.!'
3769	! CALL ipslerr_p(3,'interpweight_provide_fractions4D',TRIM(msg),'','')
3770	! END IF
3771	! END DO
3772	END DO
3773
3774	CASE DEFAULT
3775	CALL ipslerr_p(3,'interpweight_provide_fractions4D',"Inerpolation type '" // TRIM(tint) // "' not ready!",'','')
3776	END SELECT
3777
3778	END SUBROUTINE interpweight_provide_fractions4D
3779
3780	!! ================================================================================================================================
3781	!! SUBROUTINE : interpweight_provide_interpolation2D
3782	!!
3783	!>\BRIEF ! perform the interpolation for a continuos 2D field
3784	!!
3785	!! DESCRIPTION : (definitions, functional, design, flags):
3786	!!
3787	!! RECENT CHANGE(S): None
3788	!!
3789	!! MAIN OUTPUT VARIABLE(S): ovar, aovar
3790	!!
3791	!! REFERENCE(S) : None
3792	!!
3793	!! FLOWCHART : None
3794	!! \n
3795	!_ ================================================================================================================================
3796	SUBROUTINE interpweight_provide_interpolation2D(tint, npt, d1, d2, dx, dy, d3, d4, &
3797	nbmax, zeroval, ivar2D, sarea, sindex, vmax, vmin, one, two, latlon, &
3798	defaultval, defaultNOval, ovar, aovar)
3799
3800	USE constantes_var
3801
3802	IMPLICIT NONE
3803
3804	!! 0. Variables and parameter declaration
3805	!! 0.1 Input variables
3806	CHARACTER(LEN=50), INTENT(in) :: tint !! tint: type of interpolation
3807	!! 'XYKindTime': Input values are kinds
3808	!! of something with a temporal
3809	!! evolution on the dx*dy matrix
3810	INTEGER, INTENT(in) :: npt !! total number of grid points
3811	INTEGER, INTENT(in) :: d1, d2, d3, & !! shape of the input values
3812	d4, dx, dy
3813	INTEGER, INTENT(in) :: nbmax !! maximum input grid cells to use
3814	!! for in 1 grid point
3815	REAL(r_std), INTENT(in) :: zeroval,one,two !! zero 1. and 2. real values
3816	REAL(r_std), DIMENSION(dx,dy), INTENT(in) :: ivar2D !! values of the input data
3817	REAL(r_std), DIMENSION(npt,nbmax), INTENT(in) :: sarea !! overlaping area of a given input
3818	!! grid point
3819	INTEGER(i_std), DIMENSION(npt,nbmax,2), INTENT(in) :: sindex !! index of the overlaping input
3820	!! grid point
3821	REAL(r_std), DIMENSION(npt,2), INTENT(in) :: latlon !! longitude and latitude of the destiny
3822	!! grid points
3823	REAL(r_std), INTENT(in) :: defaultval !! default interpolated value, only usef if tint='default'
3824	REAL(r_std), INTENT(in) :: defaultNOval !! default interpolated NO value, only used if tint='slopecalc'
3825	!! 0.2 Modified variables
3826	REAL(r_std), INTENT(in) :: vmax, vmin !! thresholds of the values
3827	!! 0.3 Output variables
3828	REAL(r_std), DIMENSION(npt), INTENT(out) :: ovar !! output variable (on the nbpt space)
3829	REAL(r_std), DIMENSION(npt), INTENT(out) :: aovar !! availability of input data to
3830	!! interpolate output variable (on
3831	!! the nbpt space)
3832	!! >0.: As the total sum of the
3833	!! space fraction used from the
3834	!! source
3835	!! -1: Found any required values
3836	!! from input data
3837	!! 0.4 Local variables
3838	INTEGER :: i,j,k,ij,ib,idi,jj,ip,jp,jv,it
3839	INTEGER :: idi_last
3840	REAL(r_std) :: int1D
3841	CHARACTER(LEN=250) :: msg
3842
3843
3844	! Initialization
3845	!!
3846	ovar = zeroval
3847
3848	SELECT CASE (tint)
3849	CASE ('default')
3850	IF (printlev>=3) WRITE(numout,*)'Fractions following default method'
3851	IF (ALL(ivar2D == zeroval)) THEN
3852	msg = "' wrong 2D input variable"
3853	CALL ipslerr_p(3,'interpweight_provide_interpolation2D',"'" // TRIM(tint) // TRIM(msg),'','')
3854	END IF
3855
3856	DO ib = 1, npt
3857	!-
3858	!- Reinfiltration coefficient due to the slope: Calculation with parameteres maxlope_ro
3859	!-
3860	int1D = zeroval
3861
3862	! Initialize last index to the highest possible
3863	idi_last=nbmax
3864	aovar(ib) = zero
3865	DO idi=1, nbmax
3866	! Leave the do loop if all sub areas are treated, sub_area <= 0
3867	IF ( sarea(ib,idi) <= zeroval ) THEN
3868	! Set last index to the last one used
3869	idi_last=idi-1
3870	! Exit do loop
3871	EXIT
3872	END IF
3873
3874	ip = sindex(ib,idi,1)
3875	jp = sindex(ib,idi,2)
3876
3877	int1D = int1D + ivar2D(ip,jp) * sarea(ib,idi)
3878	aovar(ib) = aovar(ib) + sarea(ib,idi)
3879	ENDDO
3880
3881	IF ( idi_last >= 1 ) THEN
3882	ovar(ib) = int1D / aovar(ib)
3883	ELSE
3884	ovar(ib) = defaultval
3885	!
3886	! Quality variable
3887	aovar(ib) = -1.
3888	ENDIF
3889	ENDDO
3890	CASE ('slopecalc')
3891	IF (printlev>=3) WRITE(numout,*)'Fractions following slopecalc method'
3892	IF (ALL(ivar2D == zeroval)) THEN
3893	msg = "' requires a 2D input variable"
3894	CALL ipslerr_p(3,'interpweight_provide_interpolation2D',"'" // TRIM(tint) // TRIM(msg),'','')
3895	END IF
3896
3897	DO ib = 1, npt
3898	!-
3899	!- Reinfiltration coefficient due to the slope: Calculation with parameteres maxlope_ro
3900	!-
3901	int1D = zeroval
3902
3903	! Initialize last index to the highest possible
3904	idi_last=nbmax
3905	aovar(ib) = zero
3906	DO idi=1, nbmax
3907	! Leave the do loop if all sub areas are treated, sub_area <= 0
3908	IF ( sarea(ib,idi) <= zeroval ) THEN
3909	! Set last index to the last one used
3910	idi_last=idi-1
3911	! Exit do loop
3912	EXIT
3913	END IF
3914
3915	ip = sindex(ib,idi,1)
3916	jp = sindex(ib,idi,2)
3917
3918	int1D = int1D + MIN(ivar2D(ip,jp)/defaultNOval,un) * sarea(ib,idi)
3919	aovar(ib) = aovar(ib) + sarea(ib,idi)
3920	ENDDO
3921
3922	IF ( idi_last >= 1 ) THEN
3923	ovar(ib) = un - int1D / aovar(ib)
3924	ELSE
3925	ovar(ib) = defaultval
3926	!
3927	! Quality variable
3928	aovar(ib) = -1.
3929	ENDIF
3930	ENDDO
3931
3932	CASE DEFAULT
3933	CALL ipslerr_p(3,'interpweight_provide_interpolation2D',"Inerpolation type '" // TRIM(tint) // "' not ready!",'','')
3934	END SELECT
3935
3936	END SUBROUTINE interpweight_provide_interpolation2D
3937
3938	!! ================================================================================================================================
3939	!! SUBROUTINE : interpweight_provide_interpolation4D
3940	!!
3941	!>\BRIEF ! perform the interpolation for a continuos 4D field
3942	!!
3943	!! DESCRIPTION : (definitions, functional, design, flags):
3944	!!
3945	!! RECENT CHANGE(S): None
3946	!!
3947	!! MAIN OUTPUT VARIABLE(S): ovar, aovar
3948	!!
3949	!! REFERENCE(S) : None
3950	!!
3951	!! FLOWCHART : None
3952	!! \n
3953	!_ ================================================================================================================================
3954	SUBROUTINE interpweight_provide_interpolation4D(tint, npt, d1, d2, dx, dy, d3, d4, &
3955	nbmax, zeroval, ivar4D, sarea, sindex, vmax, vmin, one, two, latlon, &
3956	defaultval, defaultNOval, ovar, aovar)
3957
3958	USE constantes_var
3959
3960	IMPLICIT NONE
3961
3962	!! 0. Variables and parameter declaration
3963	!! 0.1 Input variables
3964	CHARACTER(LEN=50), INTENT(in) :: tint !! tint: type of interpolation
3965	!! 'XYKindTime': Input values are kinds
3966	!! of something with a temporal
3967	!! evolution on the dx*dy matrix
3968	INTEGER, INTENT(in) :: npt !! total number of grid points
3969	INTEGER, INTENT(in) :: d1, d2, d3, & !! shape of the input values
3970	d4, dx, dy
3971	INTEGER, INTENT(in) :: nbmax !! maximum input grid cells to use
3972	!! for in 1 grid point
3973	REAL(r_std), INTENT(in) :: zeroval,two,one !! zero 1. and 2. real values
3974	REAL(r_std), DIMENSION(dx,dy,d3,d4), INTENT(in) :: ivar4D !! values of the input data
3975	REAL(r_std), DIMENSION(npt,nbmax), INTENT(in) :: sarea !! overlaping area of a given input
3976	!! grid point
3977	INTEGER(i_std), DIMENSION(npt,nbmax,2), INTENT(in) :: sindex !! index of the overlaping input
3978	!! grid point
3979	REAL(r_std), DIMENSION(npt,2), INTENT(in) :: latlon !! longitude and latitude of the destiny
3980	!! grid points
3981	REAL(r_std), DIMENSION(d1,d2), INTENT(in) :: defaultval !! default interpolated value
3982	REAL(r_std), INTENT(in) :: defaultNOval !! default interpolated NO value
3983	!! 0.2 Modified variables
3984	REAL(r_std), INTENT(in) :: vmax, vmin !! thresholds of the values
3985	!! 0.3 Output variables
3986	REAL(r_std), DIMENSION(npt,d1,d2), INTENT(out) :: ovar !! output variable (on the nbpt space)
3987	REAL(r_std), DIMENSION(npt), INTENT(out) :: aovar !! availability of input data to
3988	!! interpolate output variable (on
3989	!! the nbpt space)
3990	!! >0.: As the total sum of the
3991	!! space fraction used from the
3992	!! source
3993	!! -1: Found any required values
3994	!! from input data
3995	!! 0.4 Local variables
3996	INTEGER :: i,j,k,ij,ib,idi,jj,ip,jp,jv,it
3997	INTEGER :: ilf, fopt
3998	REAL(r_std) :: totarea
3999	CHARACTER(LEN=250) :: msg
4000	! Debug
4001	INTEGER :: nbptlL
4002	REAL(r_std), DIMENSION(npt,d1,d2) :: ovartst
4003
4004
4005	! Initialization
4006	!!
4007	ovar = zero
4008	ovartst = zero
4009
4010	SELECT CASE (tint)
4011	CASE ('default')
4012	IF (printlev>=3) WRITE(numout,*)'Fractions following default method'
4013	IF ( ALL(ivar4D == zeroval) ) THEN
4014	msg = "' wrong 4D input variable"
4015	CALL ipslerr_p(3,'provide_interpolation4D',"'" // TRIM(tint) // TRIM(msg),'','')
4016	END IF
4017
4018	DO ib = 1, npt
4019
4020	fopt = COUNT(sarea(ib,:) > zero)
4021	aovar(ib) = zero
4022
4023	IF ( fopt > 0 ) THEN
4024	!-
4025	!- Reinfiltration coefficient
4026	!-
4027	totarea = zeroval
4028
4029	DO ilf=1, fopt
4030	! Leave the do loop if all sub areas are treated, sub_area <= 0
4031	ip = sindex(ib,ilf,1)
4032	jp = sindex(ib,ilf,2)
4033
4034	DO jv=1,d1
4035	DO it=1,d2
4036	ovar(ib,jv,it) = ovar(ib,jv,it) + ivar4D(ip,jp,jv,it)*sarea(ib,ilf)
4037	END DO
4038	END DO
4039	totarea = totarea + sarea(ib,ilf)
4040	ENDDO
4041	! Normalize
4042	ovar(ib,:,:) = ovar(ib,:,:)/totarea
4043	aovar(ib) = totarea
4044	ELSE
4045	! Set defalut value for points where the interpolation fail
4046	IF (printlev>=3) THEN
4047	WRITE(numout,*) 'provide_interpolation4D: no point found !!'
4048	WRITE(numout,*) 'On point ', ib, ' no points were found for interpolation data. ' // &
4049	'Default values are used.'
4050	WRITE(numout,*) 'Location : ', latlon(ib,2), latlon(ib,1)
4051	END IF
4052	ovar(ib,ivis,:) = defaultval(ivis,:)
4053	ovar(ib,inir,:) = defaultval(inir,:)
4054	aovar(ib) = -un
4055	END IF
4056	ENDDO
4057
4058	CASE DEFAULT
4059	CALL ipslerr_p(3,'provide_interpolation4D',"Inerpolation type '" // TRIM(tint) // "' not ready!",'','')
4060	END SELECT
4061
4062	END SUBROUTINE interpweight_provide_interpolation4D
4063
4064
4065	!!================================================================================================================================
4066	!! FUNCTION : interpweight_get_var2dims_file
4067	!!
4068	!>\BRIEF ! Function to get the dimensions of a given 2D variable inside a file
4069	!!
4070	!! DESCRIPTION : None
4071	!!
4072	!! RECENT CHANGE(S): None
4073	!!
4074	!! RETURN VALUE : interpweight_get_var2dims_file
4075	!!
4076	!! REFERENCE(S) : See above, module description.
4077	!!
4078	!! FLOWCHART : None
4079	!! \n
4080	!_================================================================================================================================
4081	FUNCTION interpweight_get_var2dims_file(filename, varname)
4082
4083	USE netcdf
4084
4085	IMPLICIT NONE
4086
4087	!! 0. Variables and parameter declaration
4088	!! 0.1 Input variables
4089	CHARACTER(LEN=*), INTENT(in) :: filename !! filename: name of the file to open
4090	CHARACTER(LEN=*), INTENT(in) :: varname !! varname: name of the variable
4091	!! 0.2 Modified variables
4092	!! 0.3 Output variables
4093	! Following: http://stackoverflow.com/questions/3828094/function-returning-an-array-in-fortran
4094	INTEGER, DIMENSION(2) :: interpweight_get_var2dims_file
4095	!! 0.4 Local variables
4096	INTEGER :: nid, vid, Ndims
4097	INTEGER :: rcode
4098	INTEGER, DIMENSION(2) :: dimsid
4099	CHARACTER(LEN=250) :: msg
4100
4101
4102	IF (LEN_TRIM(varname) < 1) THEN
4103	msg = " interpweight_get_var2dims_file: any variable name '" // TRIM(varname) // "' was provided !!"
4104	CALL ipslerr_p(3,'interpweight_get_var2dims_file',TRIM(msg), '', '')
4105	END IF
4106
4107	IF (is_root_prc) THEN
4108	! Open file for read only mode
4109	rcode = nf90_open(TRIM(filename), NF90_NOWRITE, nid)
4110	IF (rcode /= NF90_NOERR) CALL ipslerr_p(3, 'interpweight_get_var2dims_file', 'Error in ng90_open', &
4111	TRIM(nf90_strerror(rcode)),'')
4112
4113	rcode = nf90_inq_varid(nid, varname, vid)
4114	IF (rcode /= NF90_NOERR) CALL ipslerr_p(3, 'interpweight_get_var2dims_file', 'Error in ng90_inq_varid', &
4115	TRIM(nf90_strerror(rcode)),'')
4116
4117	rcode = nf90_inquire_variable(nid, vid, NDIMS = Ndims)
4118	IF (rcode /= NF90_NOERR) CALL ipslerr_p(3, 'interpweight_get_var2dims_file', 'Error in ng90_inquire_variable 1', &
4119	TRIM(nf90_strerror(rcode)),'')
4120
4121	IF (Ndims /= 2) THEN
4122	msg = "variable '" // TRIM(varname) // "' has not 2 dimensions!!"
4123	CALL ipslerr_p(3,'interpweight_get_var2dims_file',TRIM(msg), '', '')
4124	END IF
4125
4126	rcode = nf90_inquire_variable(nid, vid, DIMIDS = dimsid)
4127	IF (rcode /= NF90_NOERR) CALL ipslerr_p(3, 'interpweight_get_var2dims_file', 'Error in ng90_inquire_variable 2', &
4128	TRIM(nf90_strerror(rcode)),'')
4129
4130	rcode = nf90_inquire_dimension(nid, dimsid(1), LEN = interpweight_get_var2dims_file(1))
4131	IF (rcode /= NF90_NOERR) CALL ipslerr_p(3, 'interpweight_get_var2dims_file', 'Error in ng90_inquire_dimension 1', &
4132	TRIM(nf90_strerror(rcode)),'')
4133
4134	rcode = nf90_inquire_dimension(nid, dimsid(2), LEN = interpweight_get_var2dims_file(2))
4135	IF (rcode /= NF90_NOERR) CALL ipslerr_p(3, 'interpweight_get_var2dims_file', 'Error in ng90_inquire_dimension 2', &
4136	TRIM(nf90_strerror(rcode)),'')
4137
4138	rcode = NF90_CLOSE(nid)
4139	IF (rcode /= NF90_NOERR) CALL ipslerr_p(3, 'interpweight_get_var2dims_file', 'Error in ng90_close', &
4140	TRIM(nf90_strerror(rcode)),'')
4141	END IF
4142
4143	CALL bcast(interpweight_get_var2dims_file)
4144
4145	END FUNCTION interpweight_get_var2dims_file
4146
4147	!!================================================================================================================================
4148	!! FUNCTION : interpweight_get_var3dims_file
4149	!!
4150	!>\BRIEF ! Function to get the dimensions of a given 3D variable inside a file
4151	!!
4152	!! DESCRIPTION : None
4153	!!
4154	!! RECENT CHANGE(S): None
4155	!!
4156	!! RETURN VALUE : interpweight_get_var3dims_file
4157	!!
4158	!! REFERENCE(S) : See above, module description.
4159	!!
4160	!! FLOWCHART : None
4161	!! \n
4162	!_================================================================================================================================
4163	FUNCTION interpweight_get_var3dims_file(filename, varname)
4164
4165	USE netcdf
4166
4167	IMPLICIT NONE
4168
4169	!! 0. Variables and parameter declaration
4170	!! 0.1 Input variables
4171	CHARACTER(LEN=*), INTENT(in) :: filename !! filename: name of the file to open
4172	CHARACTER(LEN=*), INTENT(in) :: varname !! varname: name of the variable
4173	!! 0.2 Modified variables
4174	!! 0.3 Output variables
4175	INTEGER, DIMENSION(3) :: interpweight_get_var3dims_file
4176	!! 0.4 Local variables
4177	INTEGER :: nid, vid, Ndims
4178	INTEGER :: rcode
4179	INTEGER, DIMENSION(3) :: dimsid
4180	CHARACTER(LEN=250) :: msg
4181
4182
4183	IF (LEN_TRIM(varname) < 1) THEN
4184	msg = " interpweight_get_var3dims_file: any variable name '" // TRIM(varname) // "' was provided !!"
4185	CALL ipslerr_p(3,'interpweight_get_var3dims_file',TRIM(msg), '', '')
4186	END IF
4187
4188	IF (is_root_prc) THEN
4189	! Open file for read only mode
4190	rcode = nf90_open(TRIM(filename), NF90_NOWRITE, nid)
4191	IF (rcode /= NF90_NOERR) CALL ipslerr_p(3,'interpweight_get_var3dims_file', 'Error in ng90_open', &
4192	TRIM(nf90_strerror(rcode)),'')
4193
4194	rcode = nf90_inq_varid(nid, varname, vid)
4195	IF (rcode /= NF90_NOERR) CALL ipslerr_p(3,'interpweight_get_var3dims_file', 'Error in ng90_inq_varid', &
4196	TRIM(nf90_strerror(rcode)),'')
4197
4198	rcode = nf90_inquire_variable(nid, vid, NDIMS = Ndims)
4199	IF (rcode /= NF90_NOERR) CALL ipslerr_p(3,'interpweight_get_var3dims_file', 'Error in ng90_inquire_variable 1', &
4200	TRIM(nf90_strerror(rcode)),'')
4201
4202	IF (Ndims /= 3) THEN
4203	msg = "variable '" // TRIM(varname) // "' has not 3 dimensions!!"
4204	CALL ipslerr_p(3,'interpweight_get_var3dims_file',TRIM(msg), '', '')
4205	END IF
4206
4207	rcode = nf90_inquire_variable(nid, vid, DIMIDS = dimsid)
4208	IF (rcode /= NF90_NOERR) CALL ipslerr_p(3,'interpweight_get_var3dims_file', 'Error in ng90_inquire_variable 2', &
4209	TRIM(nf90_strerror(rcode)),'')
4210
4211	rcode = nf90_inquire_dimension(nid, dimsid(1), LEN = interpweight_get_var3dims_file(1))
4212	IF (rcode /= NF90_NOERR) CALL ipslerr_p(3,'interpweight_get_var3dims_file', 'Error in ng90_dimesion 1', &
4213	TRIM(nf90_strerror(rcode)),'')
4214
4215	rcode = nf90_inquire_dimension(nid, dimsid(2), LEN = interpweight_get_var3dims_file(2))
4216	IF (rcode /= NF90_NOERR) CALL ipslerr_p(3,'interpweight_get_var3dims_file', 'Error in ng90_dimension 2', &
4217	TRIM(nf90_strerror(rcode)),'')
4218
4219	rcode = nf90_inquire_dimension(nid, dimsid(3), LEN = interpweight_get_var3dims_file(3))
4220	IF (rcode /= NF90_NOERR) CALL ipslerr_p(3,'interpweight_get_var3dims_file', 'Error in ng90_dimension 3', &
4221	TRIM(nf90_strerror(rcode)),'')
4222
4223	rcode = NF90_CLOSE(nid)
4224	IF (rcode /= NF90_NOERR) CALL ipslerr_p(3,'interpweight_get_var3dims_file', 'Error in ng90_close', &
4225	TRIM(nf90_strerror(rcode)),'')
4226	END IF
4227	CALL bcast(interpweight_get_var3dims_file)
4228
4229	END FUNCTION interpweight_get_var3dims_file
4230
4231	!!================================================================================================================================
4232	!! FUNCTION : interpweight_get_var4dims_file
4233	!!
4234	!>\BRIEF ! Function to get the dimensions of a given 4D variable inside a file
4235	!!
4236	!! DESCRIPTION : None
4237	!!
4238	!! RECENT CHANGE(S): None
4239	!!
4240	!! RETURN VALUE : interpweight_get_var4dims_file
4241	!!
4242	!! REFERENCE(S) : See above, module description.
4243	!!
4244	!! FLOWCHART : None
4245	!! \n
4246	!_================================================================================================================================
4247	FUNCTION interpweight_get_var4dims_file(filename, varname)
4248
4249	USE netcdf
4250
4251	IMPLICIT NONE
4252
4253	!! 0. Variables and parameter declaration
4254	!! 0.1 Input variables
4255	CHARACTER(LEN=*), INTENT(in) :: filename !! filename: name of the file to open
4256	CHARACTER(LEN=*), INTENT(in) :: varname !! varname: name of the variable
4257	!! 0.2 Modified variables
4258	!! 0.3 Output variables
4259	INTEGER, DIMENSION(4) :: interpweight_get_var4dims_file
4260	!! 0.4 Local variables
4261	INTEGER :: nid, vid, Ndims
4262	INTEGER :: rcode
4263	INTEGER, DIMENSION(4) :: dimsid
4264	CHARACTER(LEN=250) :: msg
4265
4266	IF (LEN_TRIM(varname) < 1) THEN
4267	msg = " interpweight_get_var4dims_file: any variable name '" // TRIM(varname) // "' was provided !!"
4268	CALL ipslerr_p(3,'interpweight_get_var4dims_file',TRIM(msg), '', '')
4269	END IF
4270
4271	IF (is_root_prc) THEN
4272	! Open file for read only mode
4273	rcode = nf90_open(TRIM(filename), NF90_NOWRITE, nid)
4274	IF (rcode /= NF90_NOERR) CALL ipslerr_p(3,'interpweight_get_var4dims_file','Problem with nf90_open', &
4275	TRIM(nf90_strerror(rcode)),'')
4276
4277	rcode = nf90_inq_varid(nid, varname, vid)
4278	IF (rcode /= NF90_NOERR) CALL ipslerr_p(3,'interpweight_get_var4dims_file','Problem with nf90_inq_varid', &
4279	TRIM(nf90_strerror(rcode)),'')
4280
4281	rcode = nf90_inquire_variable(nid, vid, NDIMS = Ndims)
4282	IF (rcode /= NF90_NOERR) CALL ipslerr_p(3,'interpweight_get_var4dims_file','Problem with nf90_inquire_variable 1', &
4283	TRIM(nf90_strerror(rcode)),'')
4284
4285	IF (Ndims /= 4) THEN
4286	msg = "variable '" // TRIM(varname) // "' has not 4 dimensions!!"
4287	CALL ipslerr_p(3,'interpweight_get_var4dims_file',TRIM(msg), '', '')
4288	END IF
4289
4290	rcode = nf90_inquire_variable(nid, vid, DIMIDS = dimsid)
4291	IF (rcode /= NF90_NOERR) CALL ipslerr_p(3,'interpweight_get_var4dims_file','Problem with nf90_inquire_variable 2', &
4292	TRIM(nf90_strerror(rcode)),'')
4293
4294	rcode = nf90_inquire_dimension(nid, dimsid(1), LEN = interpweight_get_var4dims_file(1))
4295	IF (rcode /= NF90_NOERR) CALL ipslerr_p(3,'interpweight_get_var4dims_file','Problem with nf90_inquire_dimension 1', &
4296	TRIM(nf90_strerror(rcode)),'')
4297
4298	rcode = nf90_inquire_dimension(nid, dimsid(2), LEN = interpweight_get_var4dims_file(2))
4299	IF (rcode /= NF90_NOERR) CALL ipslerr_p(3,'interpweight_get_var4dims_file','Problem with nf90_inquire_dimension 2', &
4300	TRIM(nf90_strerror(rcode)),'')
4301
4302	rcode = nf90_inquire_dimension(nid, dimsid(3), LEN = interpweight_get_var4dims_file(3))
4303	IF (rcode /= NF90_NOERR) CALL ipslerr_p(3,'interpweight_get_var4dims_file','Problem with nf90_inquire_dimension 3', &
4304	TRIM(nf90_strerror(rcode)),'')
4305
4306	rcode = nf90_inquire_dimension(nid, dimsid(4), LEN = interpweight_get_var4dims_file(4))
4307	IF (rcode /= NF90_NOERR) CALL ipslerr_p(3,'interpweight_get_var4dims_file','Problem with nf90_inquire_dimension 4', &
4308	TRIM(nf90_strerror(rcode)),'')
4309
4310	rcode = NF90_CLOSE(nid)
4311	IF (rcode /= NF90_NOERR) CALL ipslerr_p(3,'interpweight_get_var4dims_file','Problem with nf90_close', &
4312	TRIM(nf90_strerror(rcode)),'')
4313	END iF
4314	CALL bcast(interpweight_get_var4dims_file)
4315
4316	END FUNCTION interpweight_get_var4dims_file
4317
4318	!!================================================================================================================================
4319	!! FUNCTION : interpweight_get_varNdims_file
4320	!!
4321	!>\BRIEF ! Function to get the rank(number of dimensions) of a given variable inside a file
4322	!!
4323	!! DESCRIPTION : None
4324	!!
4325	!! RECENT CHANGE(S): None
4326	!!
4327	!! RETURN VALUE : interpweight_get_varNdims_file
4328	!!
4329	!! REFERENCE(S) : See above, module description.
4330	!!
4331	!! FLOWCHART : None
4332	!! \n
4333	!_================================================================================================================================
4334	INTEGER FUNCTION interpweight_get_varNdims_file(filename, varname)
4335
4336	USE netcdf
4337
4338	IMPLICIT NONE
4339
4340	!! 0. Variables and parameter declaration
4341	!! 0.1 Input variables
4342	CHARACTER(LEN=*), INTENT(in) :: filename !! filename: name of the file to open
4343	CHARACTER(LEN=*), INTENT(in) :: varname !! varname: name of the variable
4344	!! 0.2 Modified variables
4345	!! 0.3 Output variables
4346	!! 0.4 Local variables
4347	INTEGER :: nid, vid
4348	INTEGER :: rcode
4349	CHARACTER(LEN=256) :: msg
4350
4351
4352	IF (LEN_TRIM(varname) < 1) THEN
4353	msg = " interpweight_get_varNdims_file: any variable name '" // TRIM(varname) // "' was provided !!"
4354	CALL ipslerr_p(3,'interpweight_get_varNdims_file',TRIM(msg), '', '')
4355	END IF
4356
4357	IF (is_root_prc) THEN
4358
4359	! Open file for read only mode
4360	rcode = nf90_open(TRIM(filename), NF90_NOWRITE, nid)
4361	IF (rcode /= NF90_NOERR) CALL ipslerr_p(3,'interpweight_get_varNdims_file','Problem with nf90_open', &
4362	TRIM(nf90_strerror(rcode)), '')
4363
4364	rcode = nf90_inq_varid(nid, varname, vid)
4365	IF (rcode /= NF90_NOERR) CALL ipslerr_p(3,'interpweight_get_varNdims_file','Problem with nf90_inq_varid', &
4366	TRIM(nf90_strerror(rcode)), '')
4367
4368	rcode = nf90_inquire_variable(nid, vid, NDIMS = interpweight_get_varNdims_file)
4369	IF (rcode /= NF90_NOERR) CALL ipslerr_p(3,'interpweight_get_varNdims_file','Problem with nf90_inquire_variable', &
4370	TRIM(nf90_strerror(rcode)), '')
4371
4372	rcode = NF90_CLOSE(nid)
4373	IF (rcode /= NF90_NOERR) CALL ipslerr_p(3,'interpweight_get_varNdims_file','Problem with nf90_close', &
4374	TRIM(nf90_strerror(rcode)), '')
4375	END IF
4376	CALL bcast(interpweight_get_varNdims_file)
4377
4378	END FUNCTION interpweight_get_varNdims_file
4379
4380	!!================================================================================================================================
4381	!! FUNCTION : isin_file
4382	!!
4383	!>\BRIEF ! Function to tell if a given variable is inside a file
4384	!!
4385	!! DESCRIPTION : None
4386	!!
4387	!! RECENT CHANGE(S): None
4388	!!
4389	!! RETURN VALUE : isin_file
4390	!!
4391	!! REFERENCE(S) : None
4392	!!
4393	!! FLOWCHART : None
4394	!! \n
4395	!_================================================================================================================================
4396	LOGICAL FUNCTION isin_file(filename, varname)
4397
4398	USE netcdf
4399
4400	IMPLICIT NONE
4401
4402	!! 0. Variables and parameter declaration
4403	!! 0.1 Input variables
4404	CHARACTER(LEN=*), INTENT(in) :: filename !! filename: name of the file to open
4405	CHARACTER(LEN=*), INTENT(in) :: varname !! varname: name of the variable
4406	!! 0.2 Modified variables
4407	!! 0.3 Output variables
4408	!! 0.4 Local variables
4409	INTEGER :: nid, vid, Ndims, Nvars
4410	INTEGER :: iv, rcode
4411	CHARACTER(LEN=1000) :: varinfile
4412	CHARACTER(LEN=250) :: msg
4413
4414
4415	IF (is_root_prc) THEN
4416	rcode = nf90_open(TRIM(filename), NF90_NOWRITE, nid)
4417	IF (rcode /= NF90_NOERR) CALL ipslerr_p(3,'isin_file','Problem with nf90_open',TRIM(nf90_strerror(rcode)),'')
4418
4419	rcode = nf90_inquire(nid, Ndims, Nvars)
4420	IF (rcode /= NF90_NOERR) CALL ipslerr_p(3,'isin_file','Problem with nf90_inquire',TRIM(nf90_strerror(rcode)),'')
4421
4422	DO iv=1, Nvars
4423	rcode = nf90_inquire_variable(nid, iv, name=varinfile)
4424	IF (rcode /= NF90_NOERR) CALL ipslerr_p(3,'isin_file','Problem with nf90_inquire_variable',TRIM(nf90_strerror(rcode)),'')
4425
4426	IF (TRIM(varinfile) == TRIM(varname)) THEN
4427	isin_file = .TRUE.
4428	EXIT
4429	ELSE
4430	isin_file = .FALSE.
4431	END IF
4432	END DO
4433
4434	rcode = NF90_CLOSE(nid)
4435	IF (rcode /= NF90_NOERR) CALL ipslerr_p(3,'isin_file','Problem with nf90_close',TRIM(nf90_strerror(rcode)),'')
4436	END IF
4437
4438	CALL bcast(isin_file)
4439
4440	END FUNCTION isin_file
4441
4442	!!!!!!! !!!!!! !!!!! !!!! !!! !! !
4443	! Generic functions no netCDF
4444
4445	!!================================================================================================================================
4446	!! FUNCTION : interpweight_Index1DArrayI
4447	!!
4448	!>\BRIEF ! Function to provide the first index of a given value inside a 1D intger array
4449	!!
4450	!! DESCRIPTION : None
4451	!!
4452	!! RECENT CHANGE(S): None
4453	!!
4454	!! RETURN VALUE : interpweight_Index1DArrayI
4455	!!
4456	!! REFERENCE(S) : See above, module description.
4457	!!
4458	!! FLOWCHART : None
4459	!! \n
4460	!_================================================================================================================================
4461	INTEGER FUNCTION interpweight_Index1DArrayI(array1D, d1, val)
4462
4463	IMPLICIT NONE
4464
4465	!! 0. Variables and parameter declaration
4466	!! 0.1 Input variables
4467	INTEGER, INTENT(in) :: d1 !! size of the array
4468	INTEGER, INTENT(in) :: val !! value to search
4469	INTEGER, DIMENSION(d1), INTENT(in) :: array1D !! values of the array
4470	!! 0.2 Modified variables
4471	!! 0.3 Output variables
4472	!! 0.4 Local variables
4473	INTEGER :: i
4474
4475
4476	interpweight_Index1DArrayI = -1
4477
4478	DO i=1,d1
4479	IF (array1d(i) == val) THEN
4480	interpweight_Index1DArrayI = i
4481	EXIT
4482	END IF
4483	END DO
4484
4485	END FUNCTION interpweight_Index1DArrayI
4486
4487	!!================================================================================================================================
4488	!! FUNCTION : interpweight_Index1DArrayR
4489	!!
4490	!>\BRIEF ! Function to provide the first index of a given value inside a 1D real array
4491	!!
4492	!! DESCRIPTION : None
4493	!!
4494	!! RECENT CHANGE(S): None
4495	!!
4496	!! RETURN VALUE : interpweight_Index1DArrayR
4497	!!
4498	!! REFERENCE(S) : See above, module description.
4499	!!
4500	!! FLOWCHART : None
4501	!! \n
4502	!_================================================================================================================================
4503	INTEGER FUNCTION interpweight_Index1DArrayR(array1D, d1, val)
4504
4505	IMPLICIT NONE
4506
4507	!! 0. Variables and parameter declaration
4508	!! 0.1 Input variables
4509	INTEGER, INTENT(in) :: d1 !! size of the array
4510	REAL, INTENT(in) :: val !! value to search
4511	REAL, DIMENSION(d1), INTENT(in) :: array1D !! values of the array
4512	!! 0.2 Modified variables
4513	!! 0.3 Output variables
4514	!! 0.4 Local variables
4515	INTEGER :: i
4516
4517	interpweight_Index1DArrayR = -1
4518
4519	DO i=1,d1
4520	IF (array1d(i) == val) THEN
4521	interpweight_Index1DArrayR = i
4522	EXIT
4523	END IF
4524	END DO
4525
4526	END FUNCTION interpweight_Index1DArrayR
4527
4528	!!================================================================================================================================
4529	!! FUNCTION : interpweight_Index2DArrayI
4530	!!
4531	!>\BRIEF ! Function to provide the first index of a given value inside a 2D integer array
4532	!!
4533	!! DESCRIPTION : None
4534	!!
4535	!! RECENT CHANGE(S): None
4536	!!
4537	!! RETURN VALUE : interpweight_Index2DArrayI
4538	!!
4539	!! REFERENCE(S) : See above, module description.
4540	!!
4541	!! FLOWCHART : None
4542	!! \n
4543	!_================================================================================================================================
4544	FUNCTION interpweight_Index2DArrayI(array2D, d1, d2, val)
4545
4546	IMPLICIT NONE
4547
4548	!! 0. Variables and parameter declaration
4549	!! 0.1 Input variables
4550	INTEGER, INTENT(in) :: d1, d2 !! size of the array
4551	INTEGER, INTENT(in) :: val !! value to search
4552	INTEGER, DIMENSION(d1,d2), INTENT(in) :: array2D !! values of the array
4553	!! 0.2 Modified variables
4554	!! 0.3 Output variables
4555	INTEGER, DIMENSION(2) :: interpweight_Index2DArrayI
4556	!! 0.4 Local variables
4557	INTEGER :: i, j
4558
4559
4560	interpweight_Index2DArrayI = -1
4561
4562	DO i=1,d1
4563	DO j=1,d2
4564	IF (array2d(i,j) == val) THEN
4565	interpweight_Index2DArrayI(1) = i
4566	interpweight_Index2DArrayI(2) = j
4567	EXIT
4568	END IF
4569	END DO
4570	END DO
4571
4572	END FUNCTION interpweight_Index2DArrayI
4573
4574	!!================================================================================================================================
4575	!! FUNCTION : interpweight_Index2DArrayR
4576	!!
4577	!>\BRIEF ! Function to provide the first index of a given value inside a 2D real array
4578	!!
4579	!! DESCRIPTION : None
4580	!!
4581	!! RECENT CHANGE(S): None
4582	!!
4583	!! RETURN VALUE : bm_vol
4584	!!
4585	!! REFERENCE(S) : See above, module description.
4586	!!
4587	!! FLOWCHART : None
4588	!! \n
4589	!_================================================================================================================================
4590	FUNCTION interpweight_Index2DArrayR(array2D, d1, d2, val)
4591
4592	IMPLICIT NONE
4593
4594	!! 0. Variables and parameter declaration
4595	!! 0.1 Input variables
4596	INTEGER, INTENT(in) :: d1, d2 !! size of the array
4597	REAL, INTENT(in) :: val !! value to search
4598	REAL, DIMENSION(d1,d2), INTENT(in) :: array2D !! values of the array
4599	!! 0.2 Modified variables
4600	!! 0.3 Output variables
4601	INTEGER, DIMENSION(2) :: interpweight_Index2DArrayR
4602	!! 0.4 Local variables
4603	INTEGER :: i, j
4604
4605	interpweight_Index2DArrayR = -1
4606
4607	DO i=1,d1
4608	DO j=1,d2
4609	IF (array2d(i,j) == val) THEN
4610	interpweight_Index2DArrayR(1) = i
4611	interpweight_Index2DArrayR(2) = j
4612	EXIT
4613	END IF
4614	END DO
4615	END DO
4616
4617	END FUNCTION interpweight_Index2DArrayR
4618
4619	!!================================================================================================================================
4620	!! FUNCTION : interpweight_Index1DLonLat
4621	!!
4622	!>\BRIEF ! Function to provide the first index of a given pair of
4623	!! longitude and latitude from a set of lon, lat 1D arrays
4624	!!
4625	!! DESCRIPTION : None
4626	!!
4627	!! RECENT CHANGE(S): None
4628	!!
4629	!! RETURN VALUE : interpweight_Index1DLonLat
4630	!!
4631	!! REFERENCE(S) : See above, module description.
4632	!!
4633	!! FLOWCHART : None
4634	!! \n
4635	!_================================================================================================================================
4636	INTEGER FUNCTION interpweight_Index1DLonLat(lon, lat, d1, lonval, latval)
4637
4638	USE constantes_var
4639
4640	IMPLICIT NONE
4641
4642	!! 0. Variables and parameter declaration
4643	!! 0.1 Input variables
4644	INTEGER, INTENT(in) :: d1 !! size of the array
4645	REAL(r_std), INTENT(in) :: lonval,latval !! longitude and latitude value
4646	REAL(r_std), DIMENSION(d1), INTENT(in) :: lon, lat !! longitude and latitudes to search in
4647	!! 0.2 Modified variables
4648	!! 0.3 Output variables
4649	!! 0.4 Local variables
4650	INTEGER :: i
4651	REAL(r_std) :: mindist
4652	REAL(r_std), DIMENSION(d1) :: dist
4653
4654
4655	interpweight_Index1DLonLat = -1
4656
4657	dist = SQRT((lon - lonval)2. + (lat - latval)2.)
4658	! mindist = MINVAL(ABS(dist))
4659	mindist = 0.
4660
4661	DO i=1,d1
4662	IF (dist(i) == mindist) THEN
4663	interpweight_Index1DLonLat = i
4664	EXIT
4665	END IF
4666	END DO
4667
4668	END FUNCTION interpweight_Index1DLonLat
4669
4670	!!================================================================================================================================
4671	!! FUNCTION : interpweight_Index2DLonLat
4672	!!
4673	!>\BRIEF Function to provide the first index of a given pair of longitude and latitude
4674	!! from a set of lon, lat 2D arrays
4675	!!
4676	!!
4677	!! DESCRIPTION : None
4678	!!
4679	!! RECENT CHANGE(S): None
4680	!!
4681	!! RETURN VALUE : interpweight_Index2DLonLat
4682	!!
4683	!! REFERENCE(S) : See above, module description.
4684	!!
4685	!! FLOWCHART : None
4686	!! \n
4687	!_================================================================================================================================
4688	FUNCTION interpweight_Index2DLonLat(lon, lat, d1, d2, lonval, latval)
4689
4690	USE constantes_var
4691
4692	IMPLICIT NONE
4693
4694	!! 0. Variables and parameter declaration
4695	!! 0.1 Input variables
4696	INTEGER, INTENT(in) :: d1, d2 !! size of the array
4697	REAL(r_std), INTENT(in) :: lonval,latval !! longitude and latitude value
4698	REAL(r_std), DIMENSION(d1,d2), INTENT(in) :: lon, lat !! longitude and latitudes to search in
4699	!! 0.2 Modified variables
4700	!! 0.3 Output variables
4701	INTEGER, DIMENSION(2) :: interpweight_Index2DLonLat
4702	!! 0.4 Local variables
4703	INTEGER :: i, j
4704	REAL(r_std) :: mindist
4705	REAL(r_std), DIMENSION(d1,d2) :: dist
4706
4707
4708	interpweight_Index2DLonLat = -1
4709
4710	dist = SQRT((lon - lonval)2. + (lat - latval)2.)
4711	mindist = 0.
4712
4713	DO i=1,d1
4714	DO j=1,d2
4715	IF (dist(i,j) == mindist) THEN
4716	interpweight_Index2DLonLat(1) = i
4717	interpweight_Index2DLonLat(2) = j
4718	EXIT
4719	END IF
4720	END DO
4721	END DO
4722
4723	END FUNCTION interpweight_Index2DLonLat
4724
4725	!!================================================================================================================================
4726	!! FUNCTION : interpweight_RangeI
4727	!!
4728	!>\BRIEF ! Function to provide a range of d1 values from 'iniv' to
4729	!! 'endv', of integer values in a vector
4730	!!
4731	!! DESCRIPTION : None
4732	!!
4733	!! RECENT CHANGE(S): None
4734	!!
4735	!! RETURN VALUE : interpweight_RangeI
4736	!!
4737	!! REFERENCE(S) : See above, module description.
4738	!!
4739	!! FLOWCHART : None
4740	!! \n
4741	!_================================================================================================================================
4742	FUNCTION interpweight_RangeI(d1, iniv, endv)
4743
4744	IMPLICIT NONE
4745
4746	!! 0. Variables and parameter declaration
4747	!! 0.1 Input variables
4748	INTEGER, INTENT(in) :: d1 !! number of values
4749	INTEGER, INTENT(in) :: iniv !! initial value
4750	INTEGER, INTENT(in) :: endv !! end value
4751	!! 0.2 Modified variables
4752	!! 0.3 Output variables
4753	INTEGER, DIMENSION(d1) :: interpweight_RangeI
4754	!! 0.4 Local variables
4755	INTEGER :: i, intv
4756
4757	intv = (endv - iniv) / (d1*1 - 1)
4758
4759	interpweight_RangeI(1) = iniv
4760	DO i=2,d1
4761	interpweight_RangeI(i) = interpweight_RangeI(i-1) + intv
4762	END DO
4763
4764	END FUNCTION interpweight_RangeI
4765
4766	!!================================================================================================================================
4767	!! FUNCTION : interpweight_RangeR
4768	!!
4769	!>\BRIEF ! Function to provide a range of d1 values from 'iniv' to
4770	!! 'endv', of real values in a vector
4771	!!
4772	!! DESCRIPTION : None
4773	!!
4774	!! RECENT CHANGE(S): None
4775	!!
4776	!! RETURN VALUE : interpweight_RangeR
4777	!!
4778	!! REFERENCE(S) : See above, module description.
4779	!!
4780	!! FLOWCHART : None
4781	!! \n
4782	!_================================================================================================================================
4783	FUNCTION interpweight_RangeR(d1, iniv, endv)
4784
4785	IMPLICIT NONE
4786
4787	!! 0. Variables and parameter declaration
4788	!! 0.1 Input variables
4789	INTEGER, INTENT(in) :: d1 !! number of values
4790	REAL, INTENT(in) :: iniv !! initial value
4791	REAL, INTENT(in) :: endv !! end value
4792	!! 0.2 Modified variables
4793	!! 0.3 Output variables
4794	REAL, DIMENSION(d1) :: interpweight_RangeR
4795	!! 0.4 Local variables
4796	INTEGER :: i
4797	REAL :: intv
4798
4799	intv = (endv - iniv) / (d1*1. - 1.)
4800
4801	interpweight_RangeR(1) = iniv
4802	DO i=2,d1
4803	interpweight_RangeR(i) = interpweight_RangeR(i-1) + intv
4804	END DO
4805
4806	END FUNCTION interpweight_RangeR
4807
4808	!!================================================================================================================================
4809	!! FUNCTION : interpweight_ValVecI
4810	!!
4811	!>\BRIEF ! Function to provide the number of times and where that a
4812	!! given value 'oper' on a vector of integers
4813	!!
4814	!! DESCRIPTION : None
4815	!!
4816	!! RECENT CHANGE(S): None
4817	!!
4818	!! RETURN VALUE : interpweight_ValVecI
4819	!!
4820	!! REFERENCE(S) : See above, module description.
4821	!!
4822	!! FLOWCHART : None
4823	!! \n
4824	!_================================================================================================================================
4825	FUNCTION interpweight_ValVecI(vec, d1, val, oper)
4826
4827	IMPLICIT NONE
4828
4829	!! 0. Variables and parameter declaration
4830	!! 0.1 Input variables
4831	INTEGER, INTENT(in) :: d1 !! Number of values
4832	INTEGER, DIMENSION(d1), INTENT(in) :: vec !! vector of integer values
4833	INTEGER, INTENT(in) :: val !! value to search
4834	CHARACTER(LEN=*), INTENT(in) :: oper !! operation:
4835	!! 'eq': equal
4836	!! 'ge': greater or equal
4837	!! 'le': less or equal
4838	!! 'ne': not equal
4839	!! 0.2 Modified variables
4840	!! 0.3 Output variables
4841	INTEGER, DIMENSION(d1) :: interpweight_ValVecI !! Vector with positions
4842	!! (NtimesVal, pos1, pos2, ..., posNtimesVal, ...)
4843	!! NtimesVal=number of correspondencies with 'oper'
4844	!! 0.4 Local variables
4845	INTEGER :: i, NtimesVal
4846	CHARACTER(LEN=50) :: errormsg
4847
4848	errormsg = 'ERROR -- error -- ERROR -- error'
4849
4850	NtimesVal = 0
4851	DO i=1,d1
4852	SELECT CASE(oper)
4853	CASE ('eq')
4854	IF (vec(i) == val) THEN
4855	NtimesVal = NtimesVal + 1
4856	IF (NtimesVal < d1) interpweight_ValVecI(NtimesVal + 1) = i
4857	END IF
4858	CASE ('ge')
4859	IF (vec(i) >= val) THEN
4860	NtimesVal = NtimesVal + 1
4861	IF (NtimesVal < d1) interpweight_ValVecI(NtimesVal + 1) = i
4862	END IF
4863	CASE ('le')
4864	IF (vec(i) <= val) THEN
4865	NtimesVal = NtimesVal + 1
4866	IF (NtimesVal < d1) interpweight_ValVecI(NtimesVal + 1) = i
4867	END IF
4868	CASE ('neq')
4869	IF (vec(i) /= val) THEN
4870	NtimesVal = NtimesVal + 1
4871	IF (NtimesVal < d1) interpweight_ValVecI(NtimesVal + 1) = i
4872	END IF
4873	CASE DEFAULT
4874	WRITE(numout,*)TRIM(errormsg)
4875	WRITE(numout,*)" interpweight_ValVecI: operation '" // TRIM(oper) // "' not ready!!"
4876	WRITE(numout,*)" only available: 'eq', 'ge', 'le', 'neq'"
4877	CALL ipslerr_p (3, 'interpweight_ValVecI', 'wrong operation', 'provide another one', '')
4878	END SELECT
4879	END DO
4880	interpweight_ValVecI(1) = NtimesVal
4881
4882
4883	END FUNCTION interpweight_ValVecI
4884
4885	!!================================================================================================================================
4886	!! FUNCTION : interpweight_ValVecR
4887	!!
4888	!>\BRIEF Function to provide the number of times and where that a
4889	!! given value 'oper' on a vector of reals
4890	!!
4891	!! DESCRIPTION : None
4892	!!
4893	!! RECENT CHANGE(S): None
4894	!!
4895	!! RETURN VALUE : interpweight_ValVecR
4896	!!
4897	!! REFERENCE(S) : See above, module description.
4898	!!
4899	!! FLOWCHART : None
4900	!! \n
4901	!_================================================================================================================================
4902	FUNCTION interpweight_ValVecR(vec, d1, val, oper)
4903
4904	IMPLICIT NONE
4905
4906	!! 0. Variables and parameter declaration
4907	!! 0.1 Input variables
4908	INTEGER, INTENT(in) :: d1 !! Number of values
4909	REAL, DIMENSION(d1), INTENT(in) :: vec !! vector of integer values
4910	REAL, INTENT(in) :: val !! value to search
4911	CHARACTER(LEN=*), INTENT(in) :: oper !! operation:
4912	!! 'eq': equal
4913	!! 'ge': greater or equal
4914	!! 'le': less or equal
4915	!! 'neq': not equal
4916	!! 0.2 Modified variables
4917	!! 0.3 Output variables
4918	INTEGER, DIMENSION(d1) :: interpweight_ValVecR !! Vector with positions
4919	!! (NtimesVal, pos1, pos2, ..., posNtimesVal, ...)
4920	!! NtimesVal=number of correspondencies with 'oper'
4921	!! 0.4 Local variables
4922	INTEGER :: i, NtimesVal
4923	CHARACTER(LEN=50) :: errormsg
4924
4925	errormsg = 'ERROR -- error -- ERROR -- error'
4926
4927	NtimesVal = 0
4928	DO i=1,d1
4929	SELECT CASE(oper)
4930	CASE ('eq')
4931	IF (vec(i) == val) THEN
4932	NtimesVal = NtimesVal + 1
4933	IF (NtimesVal < d1) interpweight_ValVecR(NtimesVal + 1) = i
4934	END IF
4935	CASE ('ge')
4936	IF (vec(i) >= val) THEN
4937	NtimesVal = NtimesVal + 1
4938	IF (NtimesVal < d1) interpweight_ValVecR(NtimesVal + 1) = i
4939	END IF
4940	CASE ('le')
4941	IF (vec(i) <= val) THEN
4942	NtimesVal = NtimesVal + 1
4943	IF (NtimesVal < d1) interpweight_ValVecR(NtimesVal + 1) = i
4944	END IF
4945	CASE ('neq')
4946	IF (vec(i) /= val) THEN
4947	NtimesVal = NtimesVal + 1
4948	IF (NtimesVal < d1) interpweight_ValVecR(NtimesVal + 1) = i
4949	END IF
4950	CASE DEFAULT
4951	WRITE(numout,*)TRIM(errormsg)
4952	WRITE(numout,*)" interpweight_ValVecR: operation '" // TRIM(oper) // "' not ready!!"
4953	WRITE(numout,*)" only available: 'eq', 'ge', 'le', 'neq'"
4954	CALL ipslerr_p (3, 'interpweight_ValVecR', 'wrong operation', 'provide another one', '')
4955	END SELECT
4956	END DO
4957	interpweight_ValVecR(1) = NtimesVal
4958
4959
4960	END FUNCTION interpweight_ValVecR
4961
4962	!!================================================================================================================================
4963	!! FUNCTION : interpweight_From2DmatTo1Dvec
4964	!!
4965	!>\BRIEF ! Suroutine to provide value on a transformation from a 2D
4966	!! matrix to a 1D vector
4967	!! e.g.: from T2(i,j) to T2M(ijklon)
4968	!!
4969	!! DESCRIPTION : None
4970	!!
4971	!! RECENT CHANGE(S): None
4972	!!
4973	!! RETURN VALUE : interpweight_From2DmatTo1Dvec
4974	!!
4975	!! REFERENCE(S) : See above, module description.
4976	!!
4977	!! FLOWCHART : None
4978	!! \n
4979	!_================================================================================================================================
4980	INTEGER FUNCTION interpweight_From2DmatTo1Dvec(ix,iy,dx,dy)
4981
4982	IMPLICIT NONE
4983
4984	!! 0. Variables and parameter declaration
4985	!! 0.1 Input variables
4986	INTEGER, INTENT(in) :: ix,iy !! location within the 2D matrix
4987	INTEGER, INTENT(in) :: dx,dy !! dimension of the 2D matrix
4988
4989	IF (ix > dx) THEN
4990	WRITE(numout,*) "Error: "
4991	WRITE(numout,*)" interpweight_From2DmatTo1Dvec: 'ix' too large ix > dx ", ix, '>', dx, '!!'
4992	CALL ipslerr_p (3, 'interpweight_From2DmatTo1Dvec', 'ix value too big', 'out of bounds', '')
4993	END IF
4994
4995	IF (iy > dy) THEN
4996	WRITE(numout,*) "Error: "
4997	WRITE(numout,*)" interpweight_From2DmatTo1Dvec: 'iy' too large iy > dy ", iy, '>', dy, '!!'
4998	CALL ipslerr_p (3, 'interpweight_From2DmatTo1Dvec', 'iy value too big', 'out of bounds', '')
4999	END IF
5000
5001	interpweight_From2DmatTo1Dvec = (ix-1)*dy + iy
5002
5003	END FUNCTION interpweight_From2DmatTo1Dvec
5004
5005	!!================================================================================================================================
5006	!! FUNCTION : interpweight_From1DvecTo2Dmat
5007	!!
5008	!>\BRIEF ! Suroutine to provide value on a transformation from a 1D
5009	!! vector to a 2D matrix
5010	!! e.g.: from T2M(ijklon) to T2(i,j)
5011	!!
5012	!! DESCRIPTION : None
5013	!!
5014	!! RECENT CHANGE(S): None
5015	!!
5016	!! RETURN VALUE : interpweight_From1DvecTo2Dmat
5017	!!
5018	!! REFERENCE(S) : See above, module description.
5019	!!
5020	!! FLOWCHART : None
5021	!! \n
5022	!_================================================================================================================================
5023	FUNCTION interpweight_From1DvecTo2Dmat(ivec,dx,dy)
5024
5025	IMPLICIT NONE
5026
5027	!! 0. Variables and parameter declaration
5028	!! 0.1 Input variables
5029	INTEGER, INTENT(in) :: ivec !! position within the 1D vec
5030	INTEGER, INTENT(in) :: dx,dy !! dimension of the 2D matrix
5031	!! 0.2 Modified variables
5032	!! 0.3 Output variables
5033	INTEGER, DIMENSION(2) :: interpweight_From1DvecTo2Dmat
5034
5035
5036	IF (ivec > dx*dy) THEN
5037	WRITE(numout,*) "Error: "
5038	WRITE(numout,)" interpweight_From1DvecTo2Dmat: 'ivec' too large ivec > dxdy ", ivec, '>', dx*dy, '!!'
5039	CALL ipslerr_p (3, 'interpweight_From1DvecTo2Dmat', "'ivec' too large", 'out of bounds!!', '')
5040	END IF
5041
5042	interpweight_From1DvecTo2Dmat(1) = INT((ivec-1)/dy) + 1
5043	interpweight_From1DvecTo2Dmat(2) = ivec - (interpweight_From1DvecTo2Dmat(1)-1)*dy
5044
5045	END FUNCTION interpweight_From1DvecTo2Dmat
5046
5047	END MODULE interpweight

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

Download in other formats: