Context Navigation

interpweight.f90 @ 8005

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

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

Context Navigation

source: branches/publications/ORCHIDEE_gmd-2018-57/src_global/interpweight.f90 @ 8005

Download in other formats: