Context Navigation

interpweight.f90 @ 6940

Last change on this file since 6940 was 6940, checked in by jinfeng.chang, 4 years ago
add missing files for ORCHIDEE-GMv3.2
File size: 231.8 KB

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

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

Context Navigation

source: branches/publications/ORCHIDEE-GMv3.2/ORCHIDEE/src_global/interpweight.f90 @ 6940

Download in other formats: