Context Navigation

module_ForInterpolate.f90 @ 2831

Last change on this file since 2831 was 1812, checked in by lfita, 7 years ago
Fixing right interpolation!
File size: 64.4 KB

Line
1	! Module to interpolate values from a giving projection and pressure interpolation
2	! To be included in a python
3	! Follow compilation instructions from Makefile
4	! Content
5	! LlInterpolateProjection: Subroutine which provides the indices for a given interpolation of a projection
6	! var2D_IntProj: Subroutine to interpolate a 2D variable
7	! var3D_IntProj: Subroutine to interpolate a 3D variable
8	! var4D_IntProj: Subroutine to interpolate a 4D variable
9	! var5D_IntProj: Subroutine to interpolate a 5D variable
10	MODULE module_ForInterpolate
11
12	USE module_definitions
13	USE module_generic
14
15	CONTAINS
16
17	SUBROUTINE CoarselonlatFind(dx, dy, ilon, ilat, nxlon, nxlat, fraclon, fraclat, lonv, latv, per, &
18	Nperx, Npery, ilonlat, mindiffLl)
19	! Function to search a given value from a coarser version of the data
20
21	IMPLICIT NONE
22
23	! INTEGER, PARAMETER :: r_k = KIND(1.d0)
24	INTEGER, INTENT(in) :: dx, dy
25	REAL(r_k), DIMENSION(dx,dy), INTENT(in) :: ilon, ilat
26	REAL(r_k), DIMENSION(Nperx,Npery), INTENT(in) :: fraclon, fraclat
27	REAL(r_k), INTENT(in) :: lonv, latv, per
28	REAL(r_k), DIMENSION(2), INTENT(in) :: nxlon, nxlat
29	INTEGER, INTENT(in) :: Nperx, Npery
30	INTEGER, DIMENSION(2), INTENT(out) :: ilonlat
31	REAL(r_k), INTENT(out) :: mindiffLl
32	! Local
33	REAL(r_k), DIMENSION(Nperx,Npery) :: difffraclonlat
34	REAL(r_k) :: mindifffracLl
35	INTEGER, DIMENSION(2) :: ilonlatfrac
36	INTEGER :: ixbeg, ixend, iybeg, iyend
37	INTEGER :: fracx, fracy
38	REAL(r_k) :: fraclonv, fraclatv
39	REAL(r_k), ALLOCATABLE, DIMENSION(:,:) :: difflonlat, lon, lat
40
41	! Variables
42	! ilon, ilat: original 2D matrices with the longitudes and the latitudes
43	! lonv, latv: longitude and latitude to find
44	! nxlon, nxlat: minimum and maximum longitude and latitude of the target lon,lat
45	! per: fraction of the whole domain (as percentage)
46	! Nper[x/y]: period (as fraction over 1) of the fractions of the original grid to use to explore
47	! fraclon, fraclat: longitude and latitude fractional matricies to perform the first guess
48
49	fname = 'CoarselonlatFind'
50
51	IF (lonv < nxlon(1) .OR. lonv > nxlon(2)) THEN
52	PRINT *, TRIM(ErrWarnMsg('err'))
53	PRINT *,' ' // TRIM(fname) // ': longitude outside data range!!'
54	PRINT *,' given value:', lonv,' outside (',nxlon(1),' ,',nxlon(2),' )'
55	STOP
56	END IF
57	IF (latv < nxlat(1) .OR. latv > nxlat(2)) THEN
58	PRINT *, TRIM(ErrWarnMsg('err'))
59	PRINT *,' ' // TRIM(fname) // ': latitude outside data range!!'
60	PRINT *,' given value:', latv,' outside (',nxlat(1),' ,',nxlat(2),' )'
61	STOP
62	END IF
63
64	! Initializing variables
65	ixbeg = 0
66	ixend = 0
67	iybeg = 0
68	iyend = 0
69
70	fracx = int(dx*per)
71	fracy = int(dy*per)
72
73	! PRINT *,'fraclon _______'
74	! PRINT *,fraclon
75
76	! PRINT *,'fraclat _______'
77	! PRINT *,fraclat
78
79	! Fraction point
80	difffraclonlat = SQRT((fraclon-lonv)2. + (fraclat-latv)2.)
81	mindifffracLl = MINVAL(difffraclonlat)
82	ilonlatfrac = index2DArrayR(difffraclonlat, Nperx, Npery, mindifffracLl)
83
84	! PRINT *, 'mindifffracLl:', mindifffracLl, ' ilonlatfrac:', ilonlatfrac
85	! PRINT *, 'frac lon, lat:', fraclon(ilonlatfrac(1),ilonlatfrac(2)), ' ,', &
86	! fraclat(ilonlatfrac(1),ilonlatfrac(2))
87	! PRINT *, 'values lon, lat:', lonv, latv
88
89	! Providing fraction range
90	fraclonv = fraclon(ilonlatfrac(1),ilonlatfrac(2))
91	fraclatv = fraclat(ilonlatfrac(1),ilonlatfrac(2))
92
93	IF (fraclonv >= lonv .AND. fraclatv >= latv) THEN
94	IF (ilonlatfrac(1) > 0) THEN
95	ixbeg = (ilonlatfrac(1)-1)*fracx
96	ixend = ilonlatfrac(1)*fracx+1
97	ELSE
98	ixbeg = 0
99	ixend = fracx+1
100	END IF
101	IF (ilonlatfrac(2) > 0) THEN
102	iybeg = (ilonlatfrac(2)-1)*fracy
103	iyend = ilonlatfrac(2)*fracy+1
104	ELSE
105	iybeg = 0
106	iyend = fracy+1
107	END IF
108	ELSE IF (fraclonv < lonv .AND. fraclatv >= latv) THEN
109	IF (ilonlatfrac(1) < Nperx) THEN
110	IF (ilonlatfrac(1) /= 0) THEN
111	ixbeg = (ilonlatfrac(1)-1)*fracx
112	ixend = ilonlatfrac(1)*fracx+1
113	ELSE
114	ixbeg = 0
115	ixend = fracx+1
116	END IF
117	ELSE
118	ixbeg = Nperx*fracx
119	ixend = dx+1
120	END IF
121	IF (ilonlatfrac(2) > 0) THEN
122	iybeg = (ilonlatfrac(2)-1)*fracy
123	iyend = ilonlatfrac(2)*fracy+1
124	ELSE
125	iybeg = 0
126	iyend = fracy+1
127	END IF
128	ELSE IF (fraclonv < lonv .AND. fraclatv < latv) THEN
129	IF (ilonlatfrac(1) < Nperx) THEN
130	IF (ilonlatfrac(1) /= 0) THEN
131	ixbeg = (ilonlatfrac(1)-1)*fracx
132	ixend = ilonlatfrac(1)*fracx+1
133	ELSE
134	ixbeg = 0
135	ixend = fracx+1
136	END IF
137	ELSE
138	ixbeg = Nperx*fracx
139	ixend = dx+1
140	ENDIF
141	IF (ilonlatfrac(2) < Npery) THEN
142	IF (ilonlatfrac(2) /= 0) THEN
143	iybeg = (ilonlatfrac(2)-1)*fracy
144	iyend = ilonlatfrac(2)*fracy+1
145	ELSE
146	iybeg = 0
147	iyend = fracy+1
148	END IF
149	ELSE
150	iybeg = Npery*fracy
151	iyend = dy+1
152	END IF
153	ELSE IF (fraclonv >= lonv .AND. fraclatv < latv) THEN
154	IF (ilonlatfrac(1) > 0) THEN
155	ixbeg = (ilonlatfrac(1)-1)*fracx
156	ixend = ilonlatfrac(1)*fracx+1
157	ELSE
158	ixbeg = 0
159	ixend = fracx+1
160	END IF
161	IF (ilonlatfrac(2) < Npery) THEN
162	IF (ilonlatfrac(2) /= 0) THEN
163	iybeg = (ilonlatfrac(2)-1)*fracy
164	iyend = ilonlatfrac(2)*fracy+1
165	ELSE
166	iybeg = 0
167	iyend = fracy+1
168	END IF
169	ELSE
170	iybeg = Npery*fracy
171	iyend = dy+1
172	END IF
173	END IF
174
175	IF (ALLOCATED(lon)) DEALLOCATE(lon)
176	ALLOCATE(lon(ixend-ixbeg+1, iyend-iybeg+1))
177	IF (ALLOCATED(lat)) DEALLOCATE(lat)
178	ALLOCATE(lat(ixend-ixbeg+1, iyend-iybeg+1))
179	IF (ALLOCATED(difflonlat)) DEALLOCATE(difflonlat)
180	ALLOCATE(difflonlat(ixend-ixbeg+1, iyend-iybeg+1))
181
182	lon = ilon(ixbeg:ixend,iybeg:iyend)
183	lat = ilat(ixbeg:ixend,iybeg:iyend)
184
185	! print *,'lon _______'
186	! print *,lon
187	! print *,'lat _______'
188	! print *,lat
189
190	! Find point
191	difflonlat = SQRT((lon-lonv)2. + (lat-latv)2.)
192	mindiffLl = MINVAL(difflonlat)
193	ilonlat = index2DArrayR(difflonlat, ixend-ixbeg+1, iyend-iybeg+1, mindiffLl)
194
195	ilonlat(1) = ilonlat(1) + ixbeg
196	ilonlat(2) = ilonlat(2) + iybeg
197
198	! PRINT *,'mindiffLl:', mindiffLl, ' ilatlon:', ilatlon
199	! PRINT *,'lon, lat:', lon(ilonlat(1),ilonlat(2)), ' ,', lat(ilonlat(1),ilonlat(2))
200
201	RETURN
202
203	END SUBROUTINE CoarselonlatFind
204
205	SUBROUTINE CoarselonlatFindExact(dx, dy, ilon, ilat, nxlon, nxlat, fracx, fracy, fraclon, fraclat, &
206	iv, lonv, latv, per, Nperx, Npery, mindiff, ilonlat, mindiffLl)
207	! Function to search a given value from a coarser version of the data
208
209	IMPLICIT NONE
210
211	! INTEGER, PARAMETER :: r_k = KIND(1.d0)
212	INTEGER, INTENT(in) :: dx, dy, iv
213	REAL(r_k), DIMENSION(dx,dy), INTENT(in) :: ilon, ilat
214	INTEGER, INTENT(in) :: fracx, fracy
215	REAL(r_k), DIMENSION(Nperx,Npery), INTENT(in) :: fraclon, fraclat
216	REAL(r_k), INTENT(in) :: lonv, latv, per, mindiff
217	REAL(r_k), DIMENSION(2), INTENT(in) :: nxlon, nxlat
218	INTEGER, INTENT(in) :: Nperx, Npery
219	INTEGER, DIMENSION(2), INTENT(out) :: ilonlat
220	REAL(r_k), INTENT(out) :: mindiffLl
221	! Local
222	INTEGER :: i
223	REAL(r_k), DIMENSION(Nperx,Npery) :: difffraclonlat
224	REAL(r_k) :: mindifffracLl
225	INTEGER, DIMENSION(2) :: ilonlatfrac
226	INTEGER :: ixbeg, ixend, iybeg, iyend
227	REAL(r_k) :: fraclonv, fraclatv
228	REAL(r_k), ALLOCATABLE, DIMENSION(:,:) :: difflonlat, lon, lat
229
230	! Variables
231	! ilon, ilat: original 2D matrices with the longitudes and the latitudes
232	! lonv, latv: longitude and latitude to find
233	! iv: point in the input data
234	! nxlon, nxlat: minimum and maximum longitude and latitude of the target lon,lat
235	! per: fraction of the whole domain (as percentage)
236	! Nper[x/y]: period (as fraction over 1) of the fractions of the original grid to use to explore
237	! frac[x/y]: Number of grid points for each fraction
238	! fraclon, fraclat: longitude and latitude fractional matricies to perform the first guess
239	! mindiff: authorized minimal distance between input and interpolated point
240	! ilonlat: grid point on the total lon,lat matrix
241	! mindiffLl: distance between input and interpolated point
242
243	fname = 'CoarselonlatFindExact'
244
245	IF (lonv < nxlon(1) .OR. lonv > nxlon(2)) THEN
246	PRINT *, TRIM(ErrWarnMsg('err'))
247	PRINT *,' ' // TRIM(fname) // ': longitude outside data range!!'
248	PRINT *,' given value:', lonv,' outside (',nxlon(1),' ,',nxlon(2),' )'
249	STOP
250	END IF
251	IF (latv < nxlat(1) .OR. latv > nxlat(2)) THEN
252	PRINT *, TRIM(ErrWarnMsg('err'))
253	PRINT *,' ' // TRIM(fname) // ': latitude outside data range!!'
254	PRINT *,' given value:', latv,' outside (',nxlat(1),' ,',nxlat(2),' )'
255	STOP
256	END IF
257
258	! Initializing variables
259	ixbeg = 0
260	ixend = 0
261	iybeg = 0
262	iyend = 0
263
264	! Fraction point
265	difffraclonlat = SQRT((fraclon-lonv)2. + (fraclat-latv)2.)
266	mindifffracLl = MINVAL(difffraclonlat)
267	ilonlatfrac = index2DArrayR(difffraclonlat, Nperx, Npery, mindifffracLl)
268
269	! PRINT *, 'mindifffracLl:', mindifffracLl, ' ilonlatfrac:', ilonlatfrac
270	! PRINT *, 'frac lon, lat:', fraclon(ilonlatfrac(1),ilonlatfrac(2)), ' ,', &
271	! fraclat(ilonlatfrac(1),ilonlatfrac(2))
272	! PRINT *, 'values lon, lat:', lonv, latv
273
274	! Providing fraction range
275	fraclonv = fraclon(ilonlatfrac(1),ilonlatfrac(2))
276	fraclatv = fraclat(ilonlatfrac(1),ilonlatfrac(2))
277
278	IF (fraclonv >= lonv .AND. fraclatv >= latv) THEN
279	PRINT *,'Lluis!',fraclonv, '>=', lonv,'&', fraclatv, '>=', latv
280	IF (ilonlatfrac(1) > 1) THEN
281	ixbeg = (ilonlatfrac(1)-1)*fracx
282	ixend = ilonlatfrac(1)*fracx+1
283	ELSE
284	PRINT *,'Lluis 2'
285	ixbeg = 1
286	ixend = fracx+1
287	END IF
288	IF (ilonlatfrac(2) > 1) THEN
289	iybeg = (ilonlatfrac(2)-1)*fracy
290	iyend = ilonlatfrac(2)*fracy+1
291	ELSE
292	iybeg = 1
293	iyend = fracy+1
294	END IF
295	ELSE IF (fraclonv < lonv .AND. fraclatv >= latv) THEN
296	PRINT *,'Lluis!',fraclonv, '<', lonv,'&', fraclatv, '>=', latv
297	IF (ilonlatfrac(1) < Nperx) THEN
298	PRINT *,'Lluis 2'
299	IF (ilonlatfrac(1) /= 1) THEN
300	ixbeg = (ilonlatfrac(1)-1)*fracx
301	ixend = ilonlatfrac(1)*fracx+1
302	ELSE
303	ixbeg = 1
304	ixend = fracx+1
305	END IF
306	ELSE
307	ixbeg = Nperx*fracx
308	ixend = dx+1
309	END IF
310	IF (ilonlatfrac(2) > 1) THEN
311	iybeg = (ilonlatfrac(2)-1)*fracy
312	iyend = ilonlatfrac(2)*fracy+1
313	ELSE
314	iybeg = 1
315	iyend = fracy+1
316	END IF
317	ELSE IF (fraclonv < lonv .AND. fraclatv < latv) THEN
318	PRINT *,'Lluis!',fraclonv, '<', lonv,'&', fraclatv, '<', latv
319	IF (ilonlatfrac(1) < Nperx) THEN
320	IF (ilonlatfrac(1) /= 1) THEN
321	ixbeg = (ilonlatfrac(1)-1)*fracx
322	ixend = ilonlatfrac(1)*fracx+1
323	ELSE
324	ixbeg = 1
325	ixend = fracx+1
326	END IF
327	ELSE
328	ixbeg = Nperx*fracx
329	ixend = dx+1
330	ENDIF
331	IF (ilonlatfrac(2) < Npery) THEN
332	IF (ilonlatfrac(2) /= 1) THEN
333	iybeg = (ilonlatfrac(2)-1)*fracy
334	iyend = ilonlatfrac(2)*fracy+1
335	ELSE
336	iybeg = 1
337	iyend = fracy+1
338	END IF
339	ELSE
340	iybeg = Npery*fracy
341	iyend = dy+1
342	END IF
343	ELSE IF (fraclonv >= lonv .AND. fraclatv < latv) THEN
344	PRINT *,'Llui!',fraclonv, '>=', lonv,'&', fraclatv, '<', latv
345	IF (ilonlatfrac(1) > 1) THEN
346	ixbeg = (ilonlatfrac(1)-1)*fracx
347	ixend = ilonlatfrac(1)*fracx+1
348	ELSE
349	ixbeg = 1
350	ixend = fracx+1
351	END IF
352	IF (ilonlatfrac(2) < Npery) THEN
353	IF (ilonlatfrac(2) /= 1) THEN
354	iybeg = (ilonlatfrac(2)-1)*fracy
355	iyend = ilonlatfrac(2)*fracy+1
356	ELSE
357	iybeg = 1
358	iyend = fracy+1
359	END IF
360	ELSE
361	iybeg = Npery*fracy
362	iyend = dy+1
363	END IF
364	END IF
365
366	IF (ALLOCATED(lon)) DEALLOCATE(lon)
367	ALLOCATE(lon(ixend-ixbeg+1, iyend-iybeg+1))
368	IF (ALLOCATED(lat)) DEALLOCATE(lat)
369	ALLOCATE(lat(ixend-ixbeg+1, iyend-iybeg+1))
370	IF (ALLOCATED(difflonlat)) DEALLOCATE(difflonlat)
371	ALLOCATE(difflonlat(ixend-ixbeg+1, iyend-iybeg+1))
372
373	lon = ilon(ixbeg:ixend,iybeg:iyend)
374	lat = ilat(ixbeg:ixend,iybeg:iyend)
375
376	! print *,'lon _______'
377	! print *,lon
378	! print *,'lat _______'
379	! print *,lat
380
381	! Find point
382	difflonlat = SQRT((lon-lonv)2. + (lat-latv)2.)
383	mindiffLl = MINVAL(difflonlat)
384
385	IF (mindiffLl > mindiff) THEN
386	difflonlat = SQRT((lon-lonv)2. + (lat-latv)2.)
387	mindiffLl = MINVAL(difflonlat)
388	END IF
389
390	IF (mindiffLl > mindiff) THEN
391	PRINT *,TRIM(ErrWarnMsg('err'))
392	PRINT *,' ' // TRIM(fname) // ': not equivalent point closer than:',mindiff,' found!!'
393	PRINT *,' at input point iv:', iv,' lon/lat:', lonv,', ',latv,' distance:',mindiffLl
394	PRINT *,' Fraction values _______ (',Nperx,', ',Npery ,')'
395	PRINT *,' fraclon'
396	DO i=1, Nperx
397	PRINT *,' ',fraclon(i,:)
398	END DO
399	PRINT *,' fraclat'
400	DO i=1, Nperx
401	PRINT *,' ',fraclat(i,:)
402	END DO
403	PRINT *,' frac lon, lat:', fraclon(ilonlatfrac(1),ilonlatfrac(2)), ' ,', &
404	fraclat(ilonlatfrac(1),ilonlatfrac(2))
405	PRINT *,' mindifffracLl:', mindifffracLl, ' ilonlatfrac:', ilonlatfrac
406	PRINT *,' Coarse values _______'
407	PRINT *,' indices. x:', ixbeg, ', ', ixend, ' y:', iybeg, ', ', iyend
408	PRINT *,' lon range:', '(',ilon(ixbeg,iybeg),', ',ilon(ixend,iyend),')'
409	PRINT *,' lat range:', '(',ilat(ixbeg,iybeg),', ',ilat(ixend,iyend),')'
410	PRINT *,' lon', UBOUND(lon)
411	DO i=1, ixend-ixbeg+1
412	PRINT *,' ',lon(i,:)
413	END DO
414	PRINT *,' lat', UBOUND(lat)
415	DO i=1, ixend-ixbeg+1
416	PRINT *,' ',lat(i,:)
417	END DO
418	STOP
419	END IF
420
421	ilonlat = index2DArrayR(difflonlat, ixend-ixbeg+1, iyend-iybeg+1, mindiffLl)
422
423	ilonlat(1) = ilonlat(1) + ixbeg
424	ilonlat(2) = ilonlat(2) + iybeg
425
426	! PRINT *,'mindiffLl:', mindiffLl, ' ilatlon:', ilatlon
427	! PRINT *,'lon, lat:', lon(ilonlat(1),ilonlat(2)), ' ,', lat(ilonlat(1),ilonlat(2))
428
429	RETURN
430
431	END SUBROUTINE CoarselonlatFindExact
432
433	SUBROUTINE lonlatFind(dx, dy, ilon, ilat, nxlon, nxlat, lonv, latv, ilonlat, mindiffLl)
434	! Function to search a given value from a coarser version of the data
435
436	IMPLICIT NONE
437
438	! INTEGER, PARAMETER :: r_k = KIND(1.d0)
439	INTEGER, INTENT(in) :: dx, dy
440	REAL(r_k), DIMENSION(dx,dy), INTENT(in) :: ilon, ilat
441	REAL(r_k), INTENT(in) :: lonv, latv
442	REAL(r_k), DIMENSION(2), INTENT(in) :: nxlon, nxlat
443	INTEGER, DIMENSION(2), INTENT(out) :: ilonlat
444	REAL(r_k), INTENT(out) :: mindiffLl
445	! Local
446	REAL(r_k), DIMENSION(dx,dy) :: difflonlat
447
448	! Variables
449	! ilon, ilat: original 2D matrices with the longitudes and the latitudes
450	! lonv, latv: longitude and latitude to find
451	! nxlon, nxlat: minimum and maximum longitude and latitude of the target lon,lat
452
453	fname = 'lonlatFind'
454
455	IF (lonv < nxlon(1) .OR. lonv > nxlon(2)) THEN
456	PRINT *, TRIM(ErrWarnMsg('err'))
457	PRINT *,' ' // TRIM(fname) // ': longitude outside data range!!'
458	PRINT *,' given value:', lonv,' outside (',nxlon(1),' ,',nxlon(2),' )'
459	STOP
460	END IF
461	IF (latv < nxlat(1) .OR. latv > nxlat(2)) THEN
462	PRINT *, TRIM(ErrWarnMsg('err'))
463	PRINT *,' ' // TRIM(fname) // ': latitude outside data range!!'
464	PRINT *,' given value:', latv,' outside (',nxlat(1),' ,',nxlat(2),' )'
465	STOP
466	END IF
467
468	! Find point
469	difflonlat = SQRT((ilon-lonv)2. + (ilat-latv)2.)
470	mindiffLl = MINVAL(difflonlat)
471	ilonlat = index2DArrayR(difflonlat, dx, dy, mindiffLl)
472
473	! PRINT *,'mindiffLl:', mindiffLl, ' ilatlon:', ilatlon
474	! PRINT *,'lon, lat:', lon(ilonlat(1),ilonlat(2)), ' ,', lat(ilonlat(1),ilonlat(2))
475
476	RETURN
477
478	END SUBROUTINE lonlatFind
479
480	SUBROUTINE CoarseInterpolate(projlon, projlat, lonvs, latvs, percen, mindiff, inpt, ilonlat, &
481	mindiffLl, dimx, dimy, Ninpts)
482	! Subroutine which finds the closest grid point within a projection throughout a first guest
483	! approche from percentages of the whole domain
484
485	IMPLICIT NONE
486
487	! INTEGER, PARAMETER :: r_k = KIND(1.d0)
488	INTEGER, INTENT(in) :: dimx, dimy
489	REAL(r_k), DIMENSION(dimx,dimy), INTENT(in) :: projlon, projlat
490	INTEGER, INTENT(in) :: Ninpts
491	REAL(r_k), DIMENSION(Ninpts), INTENT(in) :: inpt, lonvs, latvs
492	REAL(r_k), INTENT(in) :: mindiff, percen
493	INTEGER, DIMENSION(Ninpts,2), INTENT(out) :: ilonlat
494	REAL(r_k), DIMENSION(Ninpts), INTENT(out) :: mindiffLl
495
496	! Local
497	INTEGER :: iv,i,j
498	INTEGER :: ierr
499	INTEGER :: Ninpts1
500	REAL(r_k), DIMENSION(2) :: extremelon, extremelat
501	REAL(r_k), ALLOCATABLE, DIMENSION(:,:) :: fractionlon, fractionlat
502	INTEGER :: dfracdx, dfracdy, fracdx, fracdy
503
504	!!!!!!! Variables
505	! dimx, dimy: dimension length of the target interpolation
506	! proj[lon/lat]: longitudes and latitudes of the target interpolation
507	! Ninpts: number of points to interpolate
508	! [lon/lat]vs: longitudes and latitudes of the points to interpolate
509	! mindiff: minimal accepted distance to the target point
510	! percen: size (as percentage of the total domain) of the first guess portions to provide the first gues
511	! inpt: whether the point has already been localized (1) or not (0)
512	! ilonlat: Longitude and Latitude of the input points
513	! mindiffLl: minimum difference between target and source longitude/latitude (in degrees)
514
515	fname = 'CoarseInterpolate'
516	Ninpts1 = Ninpts/100
517
518	extremelon = (/ MINVAL(projlon), MAXVAL(projlon) /)
519	extremelat = (/ MINVAL(projlat), MAXVAL(projlat) /)
520
521	PRINT *,' ' // TRIM(fname) //' total space:', dimx, ', ', dimy, ' %', percen
522
523	dfracdx = INT(1./percen+1)
524	dfracdy = INT(1./percen+1)
525	fracdx = INT(dimx*percen)
526	fracdy = INT(dimy*percen)
527	PRINT *,' ' // TRIM(fname) //' fraction:', dfracdx, ', ', dfracdy, ' freq:', fracdx,', ',fracdy
528
529	IF (ALLOCATED(fractionlon)) DEALLOCATE(fractionlon)
530	ALLOCATE(fractionlon(dfracdx, dfracdy), STAT=ierr)
531	IF (ierr /= 0) THEN
532	PRINT *,TRIM(ErrWarnMsg('err'))
533	PRINT *,' ' // TRIM(fname) //": problem allocating 'fractionlon' !!"
534	STOP
535	END IF
536	IF (ALLOCATED(fractionlat)) DEALLOCATE(fractionlat)
537	ALLOCATE(fractionlat(dfracdx, dfracdy), STAT=ierr)
538	IF (ierr /= 0) THEN
539	PRINT *,TRIM(ErrWarnMsg('err'))
540	PRINT *,' ' // TRIM(fname) //": problem allocating 'fractionlat' !!"
541	STOP
542	END IF
543
544	DO i=1,dfracdx
545	DO j=1,dfracdy
546	fractionlon(i,j) = projlon(fracdx(i-1)+1,fracdy(j-1)+1)
547	fractionlat(i,j) = projlat(fracdx(i-1)+1,fracdy(j-1)+1)
548	! PRINT ,'i,j:',i,', ',j,' frac ij:',fracdx(i-1),', ',fracdy*(j-1),' lonlat:', fractionlon(i,j),&
549	! ', ',fractionlat(i,j)
550	END DO
551	END DO
552
553	! PRINT *,' ' // TRIM(fname) // ' fractions of:'
554	! PRINT *,' lon _______ (',dfracdx,', ',dfracdy,')'
555	! DO i=1,dfracdx
556	! PRINT *,fractionlon(i,:)
557	! END DO
558	! PRINT *,' lat_______'
559	! DO i=1,dfracdx
560	! PRINT *,fractionlat(i,:)
561	! END DO
562
563	DO iv=1,Ninpts
564	IF (inpt(iv) == 0) THEN
565	CALL CoarselonlatFind(dimx, dimy, projlon, projlat, extremelon, extremelat, fractionlon, &
566	fractionlat, lonvs(iv), latvs(iv), percen, dfracdx, dfracdy, ilonlat(iv,:), mindiffLl(iv))
567
568	IF ((mindiffLl(iv) <= mindiff) .AND. .NOT.(ilonlat(iv,1) >= 0 .AND. ilonlat(iv,1) >= 0)) THEN
569	PRINT *,TRIM(ErrWarnMsg('err'))
570	PRINT *,' ' // TRIM(fname) // ': point iv:', iv, ' at', lonvs(iv), ' ,', latvs(iv), &
571	' not relocated !!'
572	PRINT *,' mindiffl:', mindiffLl(iv), ' ilon:', ilonlat(iv,1), ' ilat:', ilonlat(iv,2)
573	STOP
574	END IF
575
576	END IF
577	END DO
578
579	END SUBROUTINE CoarseInterpolate
580
581	SUBROUTINE CoarseInterpolateExact(projlon, projlat, lonvs, latvs, percen, mindiff, ivar, newvar, &
582	newvarin, newvarinpt, newvarindiff, dimx, dimy, Ninpts)
583	! Subroutine which finds the closest grid point within a projection throughout a first guest
584	! and then whole domain approche from percentages of the whole domain
585
586	IMPLICIT NONE
587
588	! INTEGER, PARAMETER :: r_k = KIND(1.d0)
589	INTEGER, INTENT(in) :: dimx, dimy
590	REAL(r_k), DIMENSION(dimx,dimy), INTENT(in) :: projlon, projlat
591	INTEGER, INTENT(in) :: Ninpts
592	REAL(r_k), DIMENSION(Ninpts), INTENT(in) :: ivar, lonvs, latvs
593	REAL(r_k), INTENT(in) :: mindiff, percen
594	REAL(r_k), DIMENSION(dimx,dimy), INTENT(out) :: newvar
595	INTEGER, DIMENSION(dimx,dimy), INTENT(out) :: newvarin
596	INTEGER, DIMENSION(Ninpts), INTENT(out) :: newvarinpt
597	REAL(r_k), DIMENSION(Ninpts), INTENT(out) :: newvarindiff
598
599	! Local
600	INTEGER :: iv,i,j
601	INTEGER :: ierr
602	INTEGER, DIMENSION(2) :: ilonlat
603	REAL(r_k) :: mindiffLl
604	INTEGER :: Ninpts1
605	REAL(r_k), DIMENSION(2) :: extremelon, extremelat
606	REAL(r_k), ALLOCATABLE, DIMENSION(:,:) :: fractionlon, fractionlat
607	INTEGER :: dfracdx, dfracdy, fracdx, fracdy
608
609	!!!!!!! Variables
610	! dimx, dimy: dimension length of the target interpolation
611	! proj[lon/lat]: longitudes and latitudes of the target interpolation
612	! Ninpts: number of points to interpolate
613	! [lon/lat]vs: longitudes and latitudes of the points to interpolate
614	! mindiff: minimal accepted distance to the target point
615	! percen: size (as percentage of the total domain) of the first guess portions to provide the first gues
616	! ivar: values to localize in the target projection
617	! newvar: localisation of the [lon/lat]vs point in the target projection
618	! newvarin: number of point from the input data
619	! newvarinpt: integer value indicating if the value has been already located (0: no, 1: yes)
620	! newvarindiff: distance of point from the input data to the closest target point
621	! ncid: netCDF output file id
622
623	fname = 'CoarseInterpolateExact'
624	Ninpts1 = Ninpts/100
625
626	extremelon = (/ MINVAL(projlon), MAXVAL(projlon) /)
627	extremelat = (/ MINVAL(projlat), MAXVAL(projlat) /)
628
629	PRINT *,' ' // TRIM(fname) //' total space:', dimx, ', ', dimy, ' %', percen
630
631	dfracdx = INT(1./percen+1)
632	dfracdy = INT(1./percen+1)
633	fracdx = INT(dimx*percen)
634	fracdy = INT(dimy*percen)
635	PRINT *,' ' // TRIM(fname) //' fraction:', dfracdx, ', ', dfracdy, ' freq:', fracdx,', ',fracdy
636
637	IF (ALLOCATED(fractionlon)) DEALLOCATE(fractionlon)
638	ALLOCATE(fractionlon(dfracdx, dfracdy), STAT=ierr)
639	IF (ierr /= 0) THEN
640	PRINT *,TRIM(ErrWarnMsg('err'))
641	PRINT *,' ' // TRIM(fname) //": problem allocating 'fractionlon' !!"
642	STOP
643	END IF
644	IF (ALLOCATED(fractionlat)) DEALLOCATE(fractionlat)
645	ALLOCATE(fractionlat(dfracdx, dfracdy), STAT=ierr)
646	IF (ierr /= 0) THEN
647	PRINT *,TRIM(ErrWarnMsg('err'))
648	PRINT *,' ' // TRIM(fname) //": problem allocating 'fractionlat' !!"
649	STOP
650	END IF
651
652	DO i=1,dfracdx
653	DO j=1,dfracdy
654	fractionlon(i,j) = projlon(fracdx(i-1)+1,fracdy(j-1)+1)
655	fractionlat(i,j) = projlat(fracdx(i-1)+1,fracdy(j-1)+1)
656	! PRINT ,'i,j:',i,', ',j,' frac ij:',fracdx(i-1),', ',fracdy*(j-1),' lonlat:', fractionlon(i,j),&
657	! ', ',fractionlat(i,j)
658	END DO
659	END DO
660
661	! PRINT *,' ' // TRIM(fname) // ' fractions of:'
662	! PRINT *,' lon _______ (',dfracdx,', ',dfracdy,')'
663	! DO i=1,dfracdx
664	! PRINT *,fractionlon(i,:)
665	! END DO
666	! PRINT *,' lat_______'
667	! DO i=1,dfracdx
668	! PRINT *,fractionlat(i,:)
669	! END DO
670
671	DO iv=1,Ninpts
672	IF (newvarinpt(iv) == 0) THEN
673	CALL CoarselonlatFindExact(dimx, dimy, projlon, projlat, extremelon, extremelat, fracdx, fracdy,&
674	fractionlon, fractionlat, iv, lonvs(iv), latvs(iv), percen, dfracdx, dfracdy, mindiff, &
675	ilonlat, mindiffLl)
676
677	IF (mindiffLl >= mindiff) THEN
678	! percendone(iv,Ninpts,0.5,'done:')
679
680	IF (ilonlat(1) >= 0 .AND. ilonlat(1) >= 0) THEN
681	newvar(ilonlat(1),ilonlat(2)) = ivar(iv)
682	newvarin(ilonlat(1),ilonlat(2)) = iv
683	newvarinpt(iv) = 1
684	newvarindiff(iv) = mindiffLl
685	! PRINT *,'Lluis iv:', newvarin(ilonlat(1),ilonlat(2)), ' localized:', newvarinpt(iv), &
686	! ' values:', newvar(ilonlat(1),ilonlat(2)), ' invalues:', ivar(iv), ' mindist:', &
687	! newvarindiff(iv), ' point:',ilonlat
688	ELSE
689	PRINT *,TRIM(ErrWarnMsg('err'))
690	PRINT *,' ' // TRIM(fname) // ': point iv:', iv, ' at', lonvs(iv), ' ,', latvs(iv), &
691	' not relocated !!'
692	PRINT *,' mindiffl:', mindiffLl, ' ilon:', ilonlat(1), ' ilat:', ilonlat(2)
693	STOP
694	END IF
695
696	! IF (MOD(iv,Ninpts1) == 0) newnc.sync()
697	ELSE
698	PRINT *,TRIM(ErrWarnMsg('err'))
699	PRINT *,' ' // TRIM(fname) // ': for point #', iv,' lon,lat in incomplet map:', lonvs(iv), &
700	' ,', latvs(iv), ' there is not a set of lon,lat in the completed map closer than: ', &
701	mindiff, ' !!'
702	PRINT *,' found minimum difference:', mindiffLl
703	STOP
704	END IF
705	END IF
706	END DO
707
708	END SUBROUTINE CoarseInterpolateExact
709
710	SUBROUTINE LlInterpolateProjection(inlonv, inlatv, projlon, projlat, intkind, outLlw, idimx, idimy, &
711	pdimx, pdimy)
712	! Subroutine which provides the indices for a given interpolation of a projection
713
714	IMPLICIT NONE
715
716	! INTEGER, PARAMETER :: r_k = KIND(1.d0)
717	INTEGER, INTENT(in) :: idimx, idimy, pdimx, pdimy
718	REAL(r_k), DIMENSION(pdimx,pdimy), INTENT(in) :: projlon, projlat
719	REAL(r_k), DIMENSION(idimx,idimy), INTENT(in) :: inlonv, inlatv
720	CHARACTER(LEN=50), INTENT(in) :: intkind
721	REAL(r_k), DIMENSION(3,16,pdimx,pdimy), INTENT(out) :: outLlw
722
723	! Local
724	INTEGER :: i,j,iv, ix, iy
725	REAL(r_k) :: mindiffLl, dist
726	REAL(r_k) :: inclon, inclat, maxdiffprojlonlat, &
727	maxdiffinlonlat
728	REAL(r_k), DIMENSION(idimx,idimy) :: difflonlat
729	REAL(r_k), DIMENSION(idimx,idimy) :: idifflon, idifflat
730	REAL(r_k), DIMENSION(pdimx,pdimy) :: difflon, difflat
731	REAL(r_k), DIMENSION(2) :: extremelon, extremelat, ipos
732	INTEGER, DIMENSION(2) :: iLl
733
734	!!!!!!! Variables
735	! idimx, idimy: dimension length of the input projection
736	! pdimx, pdimy: dimension length of the target projection
737	! in[lon/lat]: longitudes and latitudes of the target projection
738	! proj[lon/lat]: longitudes and latitudes of the target projection
739	! intkind: kind of interpolation
740	! 'npp': nearest neighbourgh
741	! 'dis': weighted distance, grid-output for SW, NW, NE, SE
742	! outLlw: output interpolation result
743	! for point pi,pj; up to 16 different values of
744	! 1st: i-index in input projection
745	! 2nd: j-index in input projection
746	! 3rd: weight for i-index, j-index to use for ponderation (0<1.)
747	fname = 'LlInterpolateProjection'
748
749	extremelon = (/ MINVAL(projlon), MAXVAL(projlon) /)
750	extremelat = (/ MINVAL(projlat), MAXVAL(projlat) /)
751
752	! Maximum distance between grid points in input projection
753	idifflon = 0.
754	idifflat = 0.
755	idifflon(1:idimx-1,:) = inlonv(2:idimx,:)-inlonv(1:idimx-1,:)
756	idifflat(:,1:idimy-1) = inlatv(:,2:idimy)-inlatv(:,1:idimy-1)
757	maxdiffinlonlat = MAXVAL(SQRT(idifflon2. + idifflat2.))
758	! Maximum distance between grid points in target projection
759	difflon = 0.
760	difflat = 0.
761	difflon(1:pdimx-1,:) = projlon(2:pdimx,:)-projlon(1:pdimx-1,:)
762	difflat(:,1:pdimy-1) = projlat(:,2:pdimy)-projlat(:,1:pdimy-1)
763	maxdiffprojlonlat = MAXVAL(SQRT(difflon2. + difflat2.))
764
765	IF (maxdiffinlonlat > maxdiffprojlonlat) THEN
766	PRINT *,TRIM(warnmsg)
767	PRINT *,' ' //TRIM(fname)// '; input resolution: ', maxdiffinlonlat, ' is coarser than target:', &
768	maxdiffprojlonlat, ' !!'
769	END IF
770
771	! Using case outside loop to be more efficient
772	SELECT CASE(TRIM(intkind))
773	CASE ('dis')
774	inclon = inlonv(2,1) - inlonv(1,1)
775	inclat = inlatv(1,2) - inlatv(1,1)
776
777	DO i=1, pdimx
778	DO j=1, pdimy
779	! Find point
780	difflonlat = SQRT((projlon(i,j)-inlonv)2. + (projlat(i,j)-inlatv)2.)
781	mindiffLl = MINVAL(difflonlat)
782	IF ( (mindiffLl > maxdiffprojlonlat) .AND. (mindiffLl > maxdiffinlonlat)) THEN
783	outLlw(3,:,i,j) = 0.
784	outLlw(3,:,i,j) = -1.
785	ELSE
786	! Getting the four surrounding values
787	iLl = index2DArrayR(difflonlat, idimx, idimy, mindiffLl)
788	IF ( (projlon(i,j) < inlonv(iLl(1),iLl(2)) .AND. inclon > 0.) .OR. &
789	(projlon(i,j) > inlonv(iLl(1),iLl(2)) .AND. inclon < 0.) ) THEN
790	outLlw(1,1,i,j) = MAX(iLl(1)-1,1)
791	outLlw(1,2,i,j) = MAX(iLl(1)-1,1)
792	outLlw(1,3,i,j) = MIN(iLl(1),idimx)
793	outLlw(1,4,i,j) = MIN(iLl(1),idimx)
794	ELSE
795	outLlw(1,1,i,j) = MAX(iLl(1),1)
796	outLlw(1,2,i,j) = MAX(iLl(1),1)
797	outLlw(1,3,i,j) = MIN(iLl(1)+1,idimx)
798	outLlw(1,4,i,j) = MIN(iLl(1)+1,idimx)
799	END IF
800	IF ( (projlat(i,j) < inlatv(iLl(2),iLl(2)) .AND. inclat > 0.) .OR. &
801	(projlat(i,j) > inlatv(iLl(2),iLl(2)) .AND. inclat < 0.) ) THEN
802	outLlw(2,1,i,j) = MAX(iLl(2)-1,1)
803	outLlw(2,2,i,j) = MIN(iLl(2),idimy)
804	outLlw(2,3,i,j) = MIN(iLl(2),idimy)
805	outLlw(2,4,i,j) = MAX(iLl(2)-1,1)
806	ELSE
807	outLlw(2,1,i,j) = MAX(iLl(2),1)
808	outLlw(2,2,i,j) = MIN(iLl(2)+1,idimy)
809	outLlw(2,3,i,j) = MIN(iLl(2)+1,idimy)
810	outLlw(2,4,i,j) = MAX(iLl(2),1)
811	END IF
812	! Computing distances
813	!Keep the print for future checks?
814	! IF (MOD(i+j,200) == 0) THEN
815	! PRINT *,'center point:',i,j,'=', iLl,':',projlon(i,j),projlat(i,j),' incs',inclon,' ,',inclat
816	! PRINT *,'SW:', outLlw(1,1,i,j), outLlw(2,1,i,j),':',inlonv(outLlw(1,1,i,j), outLlw(2,1,i,j)),&
817	! inlatv(outLlw(1,1,i,j), outLlw(2,1,i,j))
818	! PRINT *,'NW:', outLlw(1,2,i,j), outLlw(2,2,i,j),':',inlonv(outLlw(1,2,i,j), outLlw(2,2,i,j)),&
819	! inlatv(outLlw(1,2,i,j), outLlw(2,2,i,j))
820	! PRINT *,'NE:', outLlw(1,3,i,j), outLlw(2,3,i,j),':',inlonv(outLlw(1,3,i,j), outLlw(2,3,i,j)),&
821	! inlatv(outLlw(1,3,i,j), outLlw(2,3,i,j))
822	! PRINT *,'SE:', outLlw(1,4,i,j), outLlw(2,4,i,j),':',inlonv(outLlw(1,4,i,j), outLlw(2,4,i,j)),&
823	! inlatv(outLlw(1,4,i,j), outLlw(2,4,i,j))
824	! END IF
825	DO iv=1,4
826	ix = INT(outLlw(1,iv,i,j))
827	iy = INT(outLlw(2,iv,i,j))
828	dist = SQRT( (projlon(i,j)-inlonv(ix,iy))2. + (projlat(i,j)-inlatv(ix,iy))2. )
829	IF ( dist /= 0.) THEN
830	outLlw(3,iv,i,j) = 1./dist
831	ELSE
832	outLlw(3,iv,i,j) = 1.
833	END IF
834	! IF (i+j == 2) PRINT *,'iv:',iv,'dist:',dist,'weight:',outLlw(3,iv,i,j)
835	END DO
836	! Normalizing
837	outLlw(3,:,i,j) = outLlw(3,:,i,j)/(SUM(outLlw(3,:,i,j)))
838	! IF (i+j == 2) PRINT *,'Normalized weights:',outLlw(3,:,i,j),':',SUM(outLlw(3,:,i,j))
839	END IF
840	END DO
841	END DO
842	CASE ('npp')
843	DO i=1, pdimx
844	DO j=1, pdimy
845	! Find point
846	difflonlat = SQRT((projlon(i,j)-inlonv)2. + (projlat(i,j)-inlatv)2.)
847	mindiffLl = MINVAL(difflonlat)
848	ipos = index2DArrayR(difflonlat, idimx, idimy, mindiffLl)
849	outLlw(1,1,i,j) = REAL(ipos(1))
850	outLlw(2,1,i,j) = REAL(ipos(2))
851	! We do not want that values larger that the maximum distance between target grid points
852	! PRINT *,i,j,':',mindiffLl,'maxdiffLl:',maxdiffprojlonlat
853	IF ((mindiffLl .gt. maxdiffprojlonlat) .AND. (mindiffLl > maxdiffinlonlat)) THEN
854	! PRINT *,' ' // TRIM(fname) // ': reprojected minimum distance to nearest grid point:', &
855	! mindiffLl, ' larger than the maximum distance between grid points in target projection!!'
856	outLlw(3,1,i,j) = -1.
857	ELSE
858	outLlw(3,1,i,j) = 1.
859	END IF
860	ix = INT(outLlw(1,1,i,j))
861	iy = INT(outLlw(2,1,i,j))
862	END DO
863	END DO
864	END SELECT
865
866	END SUBROUTINE LlInterpolateProjection
867
868	SUBROUTINE var2D_IntProj(var2Din, inlonv, inlatv, projlon, projlat, intkind, mask, varout, idimx, &
869	idimy, pdimx, pdimy)
870	! Subroutine to interpolate a 2D variable
871
872	IMPLICIT NONE
873
874	INTEGER, INTENT(in) :: idimx, idimy, pdimx, pdimy
875	REAL(r_k), DIMENSION(pdimx,pdimy), INTENT(in) :: projlon, projlat
876	REAL(r_k), DIMENSION(idimx,idimy), INTENT(in) :: inlonv, inlatv
877	CHARACTER(LEN=50), INTENT(in) :: intkind
878	REAL(r_k), DIMENSION(idimx,idimy), INTENT(in) :: var2Din
879	INTEGER, DIMENSION(idimx,idimy), INTENT(in) :: mask
880	REAL(r_k), DIMENSION(pdimx,pdimy), INTENT(out) :: varout
881
882	! Local
883	INTEGER :: i, j, k, iv, ix, iy
884	REAL(r_k) :: w
885	REAL(r_k), DIMENSION(3,16,pdimx,pdimy) :: outLlw
886
887	!!!!!!! Variables
888	! idimx, idimy: dimension length of the input projection
889	! pdimx, pdimy: dimension length of the target projection
890	! in[lon/lat]: longitudes and latitudes of the target projection
891	! proj[lon/lat]: longitudes and latitudes of the target projection
892	! intkind: kind of interpolation
893	! 'npp': nearest neighbourgh
894	! 'dis': weighted distance, grid-output for SW, NW, NE, SE
895	! outLlw: output interpolation result
896	! for point pi,pj; up to 16 different values of
897	! 1st: i-index in input projection
898	! 2nd: j-index in input projection
899	! 3rd: weight for i-index, j-index to use for ponderation (0<1.)
900	! var2Din: 2D variable to interpolate
901	! mask: mask of the intpu values (1: good, 0: none)
902	! varout: variable interpolated on the target projection
903	fname = 'var2D_IntProj'
904
905	CALL LlInterpolateProjection(inlonv, inlatv, projlon, projlat, intkind, outLlw, idimx, idimy, pdimx,&
906	pdimy)
907
908	SELECT CASE (intkind)
909	CASE('dis')
910	DO i=1, pdimx
911	DO j=1, pdimy
912	IF (outLlw(3,1,i,j) == -1.) THEN
913	varout(i,j) = fillVal64
914	ELSE
915	varout(i,j) = 0.
916	DO iv=1, 4
917	ix = INT(outLlw(1,iv,i,j))
918	iy = INT(outLlw(2,iv,i,j))
919	IF (mask(ix,iy) == 1) THEN
920	w = outLlw(3,iv,i,j)
921	varout(i,j) = varout(i,j) + w*var2Din(ix,iy)
922	END IF
923	END DO
924	END IF
925	END DO
926	END DO
927	CASE('npp')
928	DO i=1, pdimx
929	DO j=1, pdimy
930	ix = INT(outLlw(1,1,i,j))
931	iy = INT(outLlw(2,1,i,j))
932	IF ( (outLlw(3,1,i,j) == -1.) .OR. (mask(ix,iy) == 0) ) THEN
933	varout(i,j) = fillVal64
934	ELSE
935	varout(i,j) = var2Din(ix,iy)*outLlw(3,1,i,j)
936	END IF
937	END DO
938	END DO
939	END SELECT
940
941	END SUBROUTINE var2D_IntProj
942
943
944	SUBROUTINE var3D_IntProj(var3Din, inlonv, inlatv, projlon, projlat, intkind, mask, varout, idimx, &
945	idimy, pdimx, pdimy, d3)
946	! Subroutine to interpolate a 3D variable
947
948	IMPLICIT NONE
949
950	! INTEGER, PARAMETER :: r_k = KIND(1.d0)
951	INTEGER, INTENT(in) :: idimx, idimy, pdimx, pdimy, d3
952	REAL(r_k), DIMENSION(pdimx,pdimy), INTENT(in) :: projlon, projlat
953	REAL(r_k), DIMENSION(idimx,idimy), INTENT(in) :: inlonv, inlatv
954	CHARACTER(LEN=50), INTENT(in) :: intkind
955	REAL(r_k), DIMENSION(idimx,idimy,d3), INTENT(in) :: var3Din
956	INTEGER, DIMENSION(idimx,idimy,d3), INTENT(in) :: mask
957	REAL(r_k), DIMENSION(pdimx,pdimy,d3), INTENT(out) :: varout
958
959	! Local
960	INTEGER :: i, j, k, iv, ix, iy
961	REAL(r_k) :: w
962	REAL(r_k), DIMENSION(3,16,pdimx,pdimy) :: outLlw
963
964	!!!!!!! Variables
965	! idimx, idimy: dimension length of the input projection
966	! pdimx, pdimy: dimension length of the target projection
967	! in[lon/lat]: longitudes and latitudes of the target projection
968	! proj[lon/lat]: longitudes and latitudes of the target projection
969	! intkind: kind of interpolation
970	! 'npp': nearest neighbourgh
971	! 'dis': weighted distance, grid-output for SW, NW, NE, SE
972	! outLlw: output interpolation result
973	! for point pi,pj; up to 16 different values of
974	! 1st: i-index in input projection
975	! 2nd: j-index in input projection
976	! 3rd: weight for i-index, j-index to use for ponderation (0<1.)
977	! var3Din: 3D variable to interpolate
978	! mask: mask of the intpu values (1: good, 0: none)
979	! varout: variable interpolated on the target projection
980	fname = 'var3D_IntProj'
981
982	CALL LlInterpolateProjection(inlonv, inlatv, projlon, projlat, intkind, outLlw, idimx, idimy, pdimx,&
983	pdimy)
984
985	SELECT CASE (intkind)
986	CASE('dis')
987	DO i=1, pdimx
988	DO j=1, pdimy
989	IF (ALL(outLlw(3,:,i,j) == -1.)) THEN
990	varout(i,j,:) = fillVal64
991	ELSE
992	DO k=1, d3
993	varout(i,j,k) = 0.
994	DO iv=1, 4
995	ix = INT(outLlw(1,iv,i,j))
996	iy = INT(outLlw(2,iv,i,j))
997	IF (mask(ix,iy,k) == 1) THEN
998	w = outLlw(3,iv,i,j)
999	varout(i,j,k) = varout(i,j,k) + w*var3Din(ix,iy,k)
1000	END IF
1001	END DO
1002	END DO
1003	END IF
1004	END DO
1005	END DO
1006	CASE('npp')
1007	DO i=1, pdimx
1008	DO j=1, pdimy
1009	ix = INT(outLlw(1,1,i,j))
1010	iy = INT(outLlw(2,1,i,j))
1011	IF ( (outLlw(3,1,i,j) == -1.) .OR. (mask(ix,iy,1) == 0) ) THEN
1012	varout(i,j,:) = fillVal64
1013	ELSE
1014	DO k=1, d3
1015	varout(i,j,k) = var3Din(ix,iy,k)*outLlw(3,1,i,j)
1016	END DO
1017	END IF
1018	END DO
1019	END DO
1020	END SELECT
1021
1022	END SUBROUTINE var3D_IntProj
1023
1024	SUBROUTINE var4D_IntProj(var4Din, inlonv, inlatv, projlon, projlat, intkind, mask, varout, idimx, &
1025	idimy, pdimx, pdimy, d3, d4)
1026	! Subroutine to interpolate a 4D variable
1027
1028	IMPLICIT NONE
1029
1030	! INTEGER, PARAMETER :: r_k = KIND(1.d0)
1031	INTEGER, INTENT(in) :: idimx, idimy, pdimx, pdimy, d3, d4
1032	REAL(r_k), DIMENSION(pdimx,pdimy), INTENT(in) :: projlon, projlat
1033	REAL(r_k), DIMENSION(idimx,idimy), INTENT(in) :: inlonv, inlatv
1034	CHARACTER(LEN=50), INTENT(in) :: intkind
1035	REAL(r_k), DIMENSION(idimx,idimy,d3,d4), INTENT(in) :: var4Din
1036	INTEGER, DIMENSION(idimx,idimy,d3,d4), INTENT(in) :: mask
1037	REAL(r_k), DIMENSION(pdimx,pdimy,d3,d4), INTENT(out) :: varout
1038
1039	! Local
1040	INTEGER :: i, j, k, l, iv, ix, iy
1041	REAL(r_k) :: w
1042	REAL(r_k), DIMENSION(3,16,pdimx,pdimy) :: outLlw
1043
1044	!!!!!!! Variables
1045	! idimx, idimy: dimension length of the input projection
1046	! pdimx, pdimy: dimension length of the target projection
1047	! in[lon/lat]: longitudes and latitudes of the target projection
1048	! proj[lon/lat]: longitudes and latitudes of the target projection
1049	! intkind: kind of interpolation
1050	! 'npp': nearest neighbourgh
1051	! 'dis': weighted distance, grid-output for SW, NW, NE, SE
1052	! outLlw: output interpolation result
1053	! for point pi,pj; up to 16 different values of
1054	! 1st: i-index in input projection
1055	! 2nd: j-index in input projection
1056	! 3rd: weight for i-index, j-index to use for ponderation (0<1.)
1057	! var4Din: 4D variable to interpolate
1058	! mask: mask of the intpu values (1: good, 0: none)
1059	! varout: variable interpolated on the target projection
1060	fname = 'var4D_IntProj'
1061
1062	CALL LlInterpolateProjection(inlonv, inlatv, projlon, projlat, intkind, outLlw, idimx, idimy, pdimx,&
1063	pdimy)
1064
1065	SELECT CASE (intkind)
1066	CASE('dis')
1067	DO i=1, pdimx
1068	DO j=1, pdimy
1069	IF (ALL(outLlw(3,:,i,j) == -1.)) THEN
1070	varout(i,j,:,:) = fillVal64
1071	ELSE
1072	DO k=1, d3
1073	DO l=1, d4
1074	varout(i,j,k,l) = 0.
1075	DO iv=1, 4
1076	ix = INT(outLlw(1,iv,i,j))
1077	iy = INT(outLlw(2,iv,i,j))
1078	IF (mask(ix,iy,k,l) == 1) THEN
1079	w = outLlw(3,iv,i,j)
1080	varout(i,j,k,l) = varout(i,j,k,l) + w*var4Din(ix,iy,k,l)
1081	END IF
1082	END DO
1083	END DO
1084	END DO
1085	END IF
1086	END DO
1087	END DO
1088	CASE('npp')
1089	DO i=1, pdimx
1090	DO j=1, pdimy
1091	ix = INT(outLlw(1,1,i,j))
1092	iy = INT(outLlw(2,1,i,j))
1093	IF ( (outLlw(3,1,i,j) == -1.) .OR. (mask(ix,iy,1,1) == 0) ) THEN
1094	varout(i,j,:,:) = fillVal64
1095	ELSE
1096	DO k=1, d3
1097	DO l=1, d4
1098	varout(i,j,k,l) = var4Din(ix,iy,k,l)*outLlw(3,1,i,j)
1099	END DO
1100	END DO
1101	END IF
1102	END DO
1103	END DO
1104	END SELECT
1105
1106	END SUBROUTINE var4D_IntProj
1107
1108	SUBROUTINE var5D_IntProj(var5Din, inlonv, inlatv, projlon, projlat, intkind, mask, varout, idimx, &
1109	idimy, pdimx, pdimy, d3, d4, d5)
1110	! Subroutine to interpolate a 5D variable
1111
1112	IMPLICIT NONE
1113
1114	! INTEGER, PARAMETER :: r_k = KIND(1.d0)
1115	INTEGER, INTENT(in) :: idimx, idimy, pdimx, pdimy, d3, d4, d5
1116	REAL(r_k), DIMENSION(pdimx,pdimy), INTENT(in) :: projlon, projlat
1117	REAL(r_k), DIMENSION(idimx,idimy), INTENT(in) :: inlonv, inlatv
1118	CHARACTER(LEN=50), INTENT(in) :: intkind
1119	REAL(r_k), DIMENSION(idimx,idimy,d3,d4,d5), INTENT(in) :: var5Din
1120	INTEGER, DIMENSION(idimx,idimy,d3,d4,d5), INTENT(in) :: mask
1121	REAL(r_k), DIMENSION(pdimx,pdimy,d3,d4,d5), INTENT(out) :: varout
1122
1123	! Local
1124	INTEGER :: i, j, k, l, m, iv, ix, iy
1125	REAL(r_k) :: w
1126	REAL(r_k), DIMENSION(3,16,pdimx,pdimy) :: outLlw
1127
1128	!!!!!!! Variables
1129	! idimx, idimy: dimension length of the input projection
1130	! pdimx, pdimy: dimension length of the target projection
1131	! in[lon/lat]: longitudes and latitudes of the target projection
1132	! proj[lon/lat]: longitudes and latitudes of the target projection
1133	! intkind: kind of interpolation
1134	! 'npp': nearest neighbourgh
1135	! 'dis': weighted distance, grid-output for SW, NW, NE, SE
1136	! outLlw: output interpolation result
1137	! for point pi,pj; up to 16 different values of
1138	! 1st: i-index in input projection
1139	! 2nd: j-index in input projection
1140	! 3rd: weight for i-index, j-index to use for ponderation (0<1.)
1141	! var5Din: 5D variable to interpolate
1142	! mask: mask of the intpu values (1: good, 0: none)
1143	! varout: variable interpolated on the target projection
1144	fname = 'var5D_IntProj'
1145
1146	CALL LlInterpolateProjection(inlonv, inlatv, projlon, projlat, intkind, outLlw, idimx, idimy, pdimx,&
1147	pdimy)
1148
1149	SELECT CASE (intkind)
1150	CASE('dis')
1151	DO i=1, pdimx
1152	DO j=1, pdimy
1153	IF (ALL(outLlw(3,:,i,j) == -1.)) THEN
1154	varout(i,j,:,:,:) = fillVal64
1155	ELSE
1156	DO k=1, d3
1157	DO l=1, d4
1158	DO m=1, d5
1159	varout(i,j,k,l,m) = 0.
1160	DO iv=1, 4
1161	ix = INT(outLlw(1,iv,i,j))
1162	iy = INT(outLlw(2,iv,i,j))
1163	IF (mask(ix,iy,k,l,m) == 1) THEN
1164	w = outLlw(3,iv,i,j)
1165	varout(i,j,k,l,m) = varout(i,j,k,l,m) + w*var5Din(ix,iy,k,l,m)
1166	END IF
1167	END DO
1168	END DO
1169	END DO
1170	END DO
1171	END IF
1172	END DO
1173	END DO
1174	CASE('npp')
1175	DO i=1, pdimx
1176	DO j=1, pdimy
1177	ix = INT(outLlw(1,1,i,j))
1178	iy = INT(outLlw(2,1,i,j))
1179	IF ( (outLlw(3,1,i,j) == -1.) .OR. (mask(ix,iy,1,1,1) == 0) ) THEN
1180	varout(i,j,:,:,:) = fillVal64
1181	ELSE
1182	DO k=1, d3
1183	DO l=1, d4
1184	DO m=1, d5
1185	varout(i,j,k,l,m) = var5Din(ix,iy,k,l,m)*outLlw(3,1,i,j)
1186	END DO
1187	END DO
1188	END DO
1189	END IF
1190	END DO
1191	END DO
1192	END SELECT
1193
1194	END SUBROUTINE var5D_IntProj
1195
1196	SUBROUTINE Interpolate(projlon, projlat, lonvs, latvs, mindiff, inpt, diffs, ilonlat, dimx, dimy, &
1197	Ninpts)
1198	! Subroutine which finds the closest grid point within a projection
1199
1200	IMPLICIT NONE
1201
1202	! INTEGER, PARAMETER :: r_k = KIND(1.d0)
1203	INTEGER, INTENT(in) :: dimx, dimy
1204	REAL(r_k), DIMENSION(dimx,dimy), INTENT(in) :: projlon, projlat
1205	INTEGER, INTENT(in) :: Ninpts
1206	REAL(r_k), DIMENSION(Ninpts), INTENT(in) :: lonvs, latvs
1207	REAL(r_k), INTENT(in) :: mindiff
1208	INTEGER, DIMENSION(Ninpts), INTENT(inout) :: inpt
1209	REAL(r_k), DIMENSION(Ninpts), INTENT(out) :: diffs
1210	INTEGER, DIMENSION(Ninpts,2), INTENT(out) :: ilonlat
1211
1212	! Local
1213	INTEGER :: iv
1214	REAL(r_k) :: mindiffLl
1215	INTEGER :: Ninpts1
1216	REAL(r_k), DIMENSION(dimx,dimy) :: difflonlat
1217	REAL(r_k), DIMENSION(2) :: extremelon, extremelat
1218
1219	!!!!!!! Variables
1220	! dimx, dimy: dimension length of the target interpolation
1221	! proj[lon/lat]: longitudes and latitudes of the target interpolation
1222	! Ninpts: number of points to interpolate
1223	! [lon/lat]vs: longitudes and latitudes of the points to interpolate
1224	! mindiff: minimal accepted distance to the target point
1225	! inpt: whether the point has already been localized
1226	! diffs: distance of point from the input data to the closest target point
1227	! ilonlat: longitude and latitude of the point
1228	! ncid: netCDF output file id
1229
1230	fname = 'Interpolate'
1231	Ninpts1 = Ninpts/100
1232
1233	extremelon = (/ MINVAL(projlon), MAXVAL(projlon) /)
1234	extremelat = (/ MINVAL(projlat), MAXVAL(projlat) /)
1235
1236	DO iv=1,Ninpts
1237	IF (inpt(iv) <= 0) THEN
1238	! Not using the subroutine, not efficient!
1239	! CALL lonlatFind(dimx, dimy, projlon, projlat, extremelon, extremelat, lonvs(iv), latvs(iv), &
1240	! ilonlat, mindiffLl)
1241
1242	IF (lonvs(iv) < extremelon(1) .OR. lonvs(iv) > extremelon(2)) THEN
1243	PRINT *, TRIM(ErrWarnMsg('err'))
1244	PRINT *,' ' // TRIM(fname) // ': longitude outside data range!!'
1245	PRINT *,' given value:', lonvs(iv),' outside (',extremelon(1),' ,',extremelon(2),' )'
1246	STOP
1247	END IF
1248	IF (latvs(iv) < extremelat(1) .OR. latvs(iv) > extremelat(2)) THEN
1249	PRINT *, TRIM(ErrWarnMsg('err'))
1250	PRINT *,' ' // TRIM(fname) // ': latitude outside data range!!'
1251	PRINT *,' given value:', latvs(iv),' outside (',extremelat(1),' ,',extremelat(2),' )'
1252	STOP
1253	END IF
1254
1255	! Find point
1256	difflonlat = SQRT((projlon-lonvs(iv))2. + (projlat-latvs(iv))2.)
1257	mindiffLl = MINVAL(difflonlat)
1258	ilonlat(iv,:) = index2DArrayR(difflonlat, dimx, dimy, mindiffLl)
1259
1260	IF (mindiffLl <= mindiff) THEN
1261	! percendone(iv,Ninpts,0.5,'done:')
1262
1263	IF (ilonlat(iv,1) >= 0 .AND. ilonlat(iv,2) >= 0) THEN
1264	diffs(iv) = mindiffLl
1265	inpt(iv) = 1
1266	! PRINT *,'Lluis iv:', newvarin(ilonlat(1),ilonlat(2)), ' localized:', newvarinpt(iv), &
1267	! ' values:', newvar(ilonlat(1),ilonlat(2)), ' invalues:', ivar(iv), ' mindist:', &
1268	! newvarindiff(iv), ' point:',ilonlat
1269	ELSE
1270	PRINT *,TRIM(ErrWarnMsg('err'))
1271	PRINT *,' ' // TRIM(fname) // ': point iv:', iv, ' at', lonvs(iv), ' ,', latvs(iv), &
1272	' not relocated !!'
1273	PRINT *,' mindiffl:', mindiffLl, ' ilon:', ilonlat(iv,1), ' ilat:', ilonlat(iv,2)
1274	STOP
1275	END IF
1276
1277	! IF (MOD(iv,Ninpts1) == 0) newnc.sync()
1278	ELSE
1279	! Because doing boxes and Goode is not conitnuos, we should jump this error message
1280	PRINT *,TRIM(ErrWarnMsg('err'))
1281	PRINT *,' ' // TRIM(fname) // ': for point #', iv,' lon,lat in incomplet map:', lonvs(iv), &
1282	' ,', latvs(iv), ' there is not a set of lon,lat in the completed map closer than: ', &
1283	mindiff, ' !!'
1284	PRINT *,' found minimum difference:', mindiffLl
1285	STOP
1286	END IF
1287	END IF
1288	END DO
1289
1290	END SUBROUTINE Interpolate
1291
1292	SUBROUTINE Interpolate1DLl(projlon, projlat, lonvs, latvs, mindiff, inpt, diffs, ilonlat, dimx, dimy, &
1293	Ninpts)
1294	! Subroutine which finds the closest grid point within a projection with 1D longitudes and latitudes
1295
1296	IMPLICIT NONE
1297
1298	! INTEGER, PARAMETER :: r_k = KIND(1.d0)
1299	INTEGER, INTENT(in) :: dimx, dimy
1300	REAL(r_k), DIMENSION(dimx), INTENT(in) :: projlon
1301	REAL(r_k), DIMENSION(dimy), INTENT(in) :: projlat
1302	INTEGER, INTENT(in) :: Ninpts
1303	REAL(r_k), DIMENSION(Ninpts), INTENT(in) :: lonvs, latvs
1304	REAL(r_k), INTENT(in) :: mindiff
1305	INTEGER, DIMENSION(Ninpts), INTENT(inout) :: inpt
1306	REAL(r_k), DIMENSION(Ninpts), INTENT(out) :: diffs
1307	INTEGER, DIMENSION(Ninpts,2), INTENT(out) :: ilonlat
1308
1309	! Local
1310	INTEGER :: iv
1311	REAL(r_k) :: mindifflo, mindiffLa, mindiffLl
1312	INTEGER :: Ninpts1
1313	REAL(r_k), DIMENSION(dimx) :: difflon
1314	REAL(r_k), DIMENSION(dimy) :: difflat
1315	REAL(r_k), DIMENSION(2) :: extremelon, extremelat
1316
1317	!!!!!!! Variables
1318	! dimx, dimy: dimension length of the target interpolation
1319	! proj[lon/lat]: longitudes and latitudes of the target interpolation
1320	! Ninpts: number of points to interpolate
1321	! [lon/lat]vs: longitudes and latitudes of the points to interpolate
1322	! mindiff: minimal accepted distance to the target point
1323	! inpt: whether the point has already been localized
1324	! diffs: distance of point from the input data to the closest target point
1325	! ilonlat: longitude and latitude of the point
1326	! ncid: netCDF output file id
1327
1328	fname = 'Interpolate1DLl'
1329	Ninpts1 = Ninpts/100
1330
1331	extremelon = (/ MINVAL(projlon), MAXVAL(projlon) /)
1332	extremelat = (/ MINVAL(projlat), MAXVAL(projlat) /)
1333
1334	DO iv=1,Ninpts
1335	IF (inpt(iv) <= 0) THEN
1336	! Not using the subroutine, not efficient!
1337	! CALL lonlatFind(dimx, dimy, projlon, projlat, extremelon, extremelat, lonvs(iv), latvs(iv), &
1338	! ilonlat, mindiffLl)
1339
1340	IF (lonvs(iv) < extremelon(1) .OR. lonvs(iv) > extremelon(2)) THEN
1341	PRINT *, TRIM(ErrWarnMsg('err'))
1342	PRINT *,' ' // TRIM(fname) // ': longitude outside data range!!'
1343	PRINT *,' given value:', lonvs(iv),' outside (',extremelon(1),' ,',extremelon(2),' )'
1344	STOP
1345	END IF
1346	IF (latvs(iv) < extremelat(1) .OR. latvs(iv) > extremelat(2)) THEN
1347	PRINT *, TRIM(ErrWarnMsg('err'))
1348	PRINT *,' ' // TRIM(fname) // ': latitude outside data range!!'
1349	PRINT *,' given value:', latvs(iv),' outside (',extremelat(1),' ,',extremelat(2),' )'
1350	STOP
1351	END IF
1352
1353	! Find point
1354	difflon = SQRT((projlon-lonvs(iv))**2.)
1355	difflat = SQRT((projlat-latvs(iv))**2.)
1356	mindifflo = MINVAL(difflon)
1357	mindiffLa = MINVAL(difflat)
1358	mindifflL = SQRT(mindifflomindifflo + mindiffLamindiffLa)
1359	ilonlat(iv,1) = index1DArrayR(difflon, dimx, mindifflo)
1360	ilonlat(iv,2) = index1DArrayR(difflat, dimy, mindiffLa)
1361	! PRINT *,' Lluis: iv',iv,' lonvs:', lonvs(iv),' latvs:',latvs(iv)
1362	! PRINT *,' Lluis: mindifflo:', mindifflo,' ilonlat(1):',ilonlat(iv,1)
1363	! PRINT *,' Lluis: mindiffLa:', mindiffLa,' ilonlat(2):',ilonlat(iv,2)
1364
1365
1366	IF (mindiffLl <= mindiff) THEN
1367	! percendone(iv,Ninpts,0.5,'done:')
1368
1369	IF (ilonlat(iv,1) >= 1 .AND. ilonlat(iv,2) >= 1) THEN
1370	diffs(iv) = mindiffLl
1371	inpt(iv) = 1
1372	! PRINT *,'Lluis iv:', newvarin(ilonlat(1),ilonlat(2)), ' localized:', newvarinpt(iv), &
1373	! ' values:', newvar(ilonlat(1),ilonlat(2)), ' invalues:', ivar(iv), ' mindist:', &
1374	! newvarindiff(iv), ' point:',ilonlat
1375	ELSE
1376	PRINT *,TRIM(ErrWarnMsg('err'))
1377	PRINT *,' ' // TRIM(fname) // ': point iv:', iv, ' at', lonvs(iv), ' ,', latvs(iv), &
1378	' not relocated !!'
1379	PRINT *,' mindiffl:', mindiffLl, ' ilon:', ilonlat(iv,1), ' ilat:', ilonlat(iv,2)
1380	STOP
1381	END IF
1382
1383	! IF (MOD(iv,Ninpts1) == 0) newnc.sync()
1384	ELSE
1385	! Because doing boxes and Goode is not conitnuos, we should jump this error message
1386	PRINT *,TRIM(ErrWarnMsg('err'))
1387	PRINT *,' ' // TRIM(fname) // ': for point #', iv,' lon,lat in incomplet map:', lonvs(iv), &
1388	' ,', latvs(iv), ' there is not a set of lon,lat in the completed map closer than: ', &
1389	mindiff, ' !!'
1390	PRINT *,' found minimum difference:', mindiffLl
1391	STOP
1392	END IF
1393	END IF
1394	END DO
1395
1396	END SUBROUTINE Interpolate1DLl
1397
1398
1399	SUBROUTINE interp (data_in, pres_field, interp_levels, psfc, ter, tk, qv, LINLOG, extrapolate, &
1400	GEOPT, MISSING, data_out, ix, iy, iz, it, num_metgrid_levels)
1401	! Interpolation subroutine from the p_interp.F90 NCAR program
1402	! Program to read wrfout data and interpolate to pressure levels
1403	! The program reads namelist.pinterp
1404	! November 2007 - Cindy Bruyere
1405	!
1406	INTEGER, INTENT(IN) :: ix, iy, iz, it
1407	INTEGER, INTENT(IN) :: num_metgrid_levels, LINLOG
1408	REAL(r_k),DIMENSION(ix,iy,iz,it), INTENT(IN) :: data_in, pres_field, tk, qv
1409	REAL(r_k),DIMENSION(ix,iy,it), INTENT(IN) :: psfc
1410	REAL(r_k),DIMENSION(ix,iy), INTENT(IN) :: ter
1411	REAL(r_k),DIMENSION(num_metgrid_levels), INTENT(IN) :: interp_levels
1412	INTEGER, INTENT(IN) :: extrapolate
1413	REAL(r_k),INTENT(IN) :: MISSING
1414	LOGICAL, INTENT(IN) :: GEOPT
1415	REAL(r_k),DIMENSION(ix,iy,num_metgrid_levels,it), &
1416	INTENT(OUT) :: data_out
1417
1418	! Local
1419	INTEGER :: i, j, itt, k, kk, kin
1420	INTEGER :: kupper
1421	REAL(r_k),DIMENSION(num_metgrid_levels) :: data_out1D
1422	REAL(r_k),DIMENSION(iz) :: data_in1D, pres_field1D
1423	REAL(r_k),DIMENSION(ix, iy, num_metgrid_levels, it) :: N
1424	REAL(r_k) :: sumA, sumN, AVE_geopt
1425	REAL(r_k) :: dp, dpmin, expon
1426	REAL(r_k) :: pbot, ptarget, tbotextrap, tvbotextrap, &
1427	zbot
1428
1429	!!!!!!! Variables
1430	! data_out: interpolated field
1431	! data_in: field to interpolate
1432	! pres_field: pressure field [Pa]
1433	! interp_levels: pressure levels to interpolate [hPa]
1434	! psfc: surface pressure [Pa]
1435	! ter: terrein height [m]
1436	! tk: temperature [K]
1437	! qv: mositure mizing ratio [kg/kg]
1438	! i[x/y/z/t]: size of the matrices
1439	! num_metgrid_levels: number of pressure values to interpolate
1440	! LINLOG: if abs(linlog)=1 use linear interp in pressure
1441	! if abs(linlog)=2 linear interp in ln(pressure)
1442	! extrapolate: whether to set to missing value below/above model ground and top (0), or extrapolate (1)
1443	! GEOPT: Wether the file is the geopotential file or not
1444	! MISSING: Missing value
1445
1446	N = 1.0
1447
1448	expon=287.04*.0065/9.81
1449
1450	do itt = 1, it
1451	do j = 1, iy
1452	do i = 1, ix
1453	data_in1D(:) = data_in(i,j,:,itt)
1454	pres_field1D(:) = pres_field(i,j,:,itt)
1455	CALL int1D (data_out1D, data_in1D, pres_field1D, interp_levels, iz, num_metgrid_levels, LINLOG, MISSING)
1456	data_out(i,j,:,itt) = data_out1D(:)
1457	end do
1458	end do
1459	end do
1460
1461
1462	! Fill in missing values
1463	IF ( extrapolate == 0 ) RETURN !! no extrapolation - we are out of here
1464
1465	! First find where about 400 hPa is located
1466	kk = 0
1467	find_kk : do k = 1, num_metgrid_levels
1468	kk = k
1469	if ( interp_levels(k) <= 40000. ) exit find_kk
1470	end do find_kk
1471
1472
1473	IF ( GEOPT ) THEN !! geopt is treated different below ground
1474
1475	do itt = 1, it
1476	do k = 1, kk
1477	do j = 1, iy
1478	do i = 1, ix
1479	IF ( data_out(i,j,k,itt) == MISSING .AND. interp_levels(k) < psfc(i,j,itt) ) THEN
1480
1481	! We are below the first model level, but above the ground
1482
1483	data_out(i,j,k,itt) = ((interp_levels(k) - pres_field(i,j,1,itt))ter(i,j)9.81 + &
1484	(psfc(i,j,itt) - interp_levels(k))*data_in(i,j,1,itt) ) / &
1485	(psfc(i,j,itt) - pres_field(i,j,1,itt))
1486
1487	ELSEIF ( data_out(i,j,k,itt) == MISSING ) THEN
1488
1489	! We are below both the ground and the lowest data level.
1490
1491	! First, find the model level that is closest to a "target" pressure
1492	! level, where the "target" pressure is delta-p less that the local
1493	! value of a horizontally smoothed surface pressure field. We use
1494	! delta-p = 150 hPa here. A standard lapse rate temperature profile
1495	! passing through the temperature at this model level will be used
1496	! to define the temperature profile below ground. This is similar
1497	! to the Benjamin and Miller (1990) method, except that for
1498	! simplicity, they used 700 hPa everywhere for the "target" pressure.
1499	! Code similar to what is implemented in RIP4
1500
1501	ptarget = (psfc(i,j,itt)*.01) - 150.
1502	dpmin=1.e4
1503	kupper = 0
1504	loop_kIN : do kin=iz,1,-1
1505	kupper = kin
1506	dp=abs( (pres_field(i,j,kin,itt)*.01) - ptarget )
1507	if (dp.gt.dpmin) exit loop_kIN
1508	dpmin=min(dpmin,dp)
1509	enddo loop_kIN
1510
1511	pbot=max(pres_field(i,j,1,itt),psfc(i,j,itt))
1512	zbot=min(data_in(i,j,1,itt)/9.81,ter(i,j))
1513
1514	tbotextrap=tk(i,j,kupper,itt)(pbot/pres_field(i,j,kupper,itt))*expon
1515	tvbotextrap=virtual(tbotextrap,qv(i,j,1,itt))
1516
1517	data_out(i,j,k,itt) = (zbot+tvbotextrap/.0065(1.-(interp_levels(k)/pbot)expon))9.81
1518
1519	ENDIF
1520	enddo
1521	enddo
1522	enddo
1523	enddo
1524
1525
1526	!!! Code for filling missing data with an average - we don't want to do this
1527	!!do itt = 1, it
1528	!!loop_levels : do k = 1, num_metgrid_levels
1529	!!sumA = SUM(data_out(:,:,k,itt), MASK = data_out(:,:,k,itt) /= MISSING)
1530	!!sumN = SUM(N(:,:,k,itt), MASK = data_out(:,:,k,itt) /= MISSING)
1531	!!IF ( sumN == 0. ) CYCLE loop_levels
1532	!!AVE_geopt = sumA/sumN
1533	!!WHERE ( data_out(:,:,k,itt) == MISSING )
1534	!!data_out(:,:,k,itt) = AVE_geopt
1535	!!END WHERE
1536	!!end do loop_levels
1537	!!end do
1538
1539	END IF
1540
1541	!!! All other fields and geopt at higher levels come here
1542	do itt = 1, it
1543	do j = 1, iy
1544	do i = 1, ix
1545	do k = 1, kk
1546	if ( data_out(i,j,k,itt) == MISSING ) data_out(i,j,k,itt) = data_in(i,j,1,itt)
1547	end do
1548	do k = kk+1, num_metgrid_levels
1549	if ( data_out(i,j,k,itt) == MISSING ) data_out(i,j,k,itt) = data_in(i,j,iz,itt)
1550	end do
1551	end do
1552	end do
1553	end do
1554
1555	END SUBROUTINE interp
1556
1557	SUBROUTINE int1D(xxout, xxin, ppin, ppout, npin, npout, LINLOG, MISSING)
1558
1559	! Modified from int2p - NCL code
1560	! routine to interpolate from one set of pressure levels
1561	! . to another set using linear or ln(p) interpolation
1562	!
1563	! NCL: xout = int2p (pin,xin,pout,linlog)
1564	! This code was originally written for a specific purpose.
1565	! . Several features were added for incorporation into NCL's
1566	! . function suite including linear extrapolation.
1567	!
1568	! nomenclature:
1569	!
1570	! . ppin - input pressure levels. The pin can be
1571	! . be in ascending or descending order
1572	! . xxin - data at corresponding input pressure levels
1573	! . npin - number of input pressure levels >= 2
1574	! . ppout - output pressure levels (input by user)
1575	! . same (ascending or descending) order as pin
1576	! . xxout - data at corresponding output pressure levels
1577	! . npout - number of output pressure levels
1578	! . linlog - if abs(linlog)=1 use linear interp in pressure
1579	! . if abs(linlog)=2 linear interp in ln(pressure)
1580	! . missing- missing data code.
1581
1582	! ! input types
1583	INTEGER :: npin,npout,linlog,ier
1584	REAL(r_k) :: ppin(npin),xxin(npin),ppout(npout)
1585	REAL(r_k) :: MISSING
1586	logical :: AVERAGE
1587	! ! output
1588	REAL(r_k) :: xxout(npout)
1589	INTEGER :: j1,np,nl,nin,nlmax,nplvl
1590	INTEGER :: nlsave,np1,no1,n1,n2,nlstrt
1591	REAL(r_k) :: slope,pa,pb,pc
1592
1593	! automatic arrays
1594	REAL(r_k) :: pin(npin),xin(npin),p(npin),x(npin)
1595	REAL(r_k) :: pout(npout),xout(npout)
1596
1597
1598	xxout = MISSING
1599	pout = ppout
1600	p = ppin
1601	x = xxin
1602	nlmax = npin
1603
1604	! exact p-level matches
1605	nlstrt = 1
1606	nlsave = 1
1607	do np = 1,npout
1608	xout(np) = MISSING
1609	do nl = nlstrt,nlmax
1610	if (pout(np).eq.p(nl)) then
1611	xout(np) = x(nl)
1612	nlsave = nl + 1
1613	go to 10
1614	end if
1615	end do
1616	10 nlstrt = nlsave
1617	end do
1618
1619	if (LINLOG.eq.1) then
1620	do np = 1,npout
1621	do nl = 1,nlmax - 1
1622	if (pout(np).lt.p(nl) .and. pout(np).gt.p(nl+1)) then
1623	slope = (x(nl)-x(nl+1))/ (p(nl)-p(nl+1))
1624	xout(np) = x(nl+1) + slope* (pout(np)-p(nl+1))
1625	end if
1626	end do
1627	end do
1628	elseif (LINLOG.eq.2) then
1629	do np = 1,npout
1630	do nl = 1,nlmax - 1
1631	if (pout(np).lt.p(nl) .and. pout(np).gt.p(nl+1)) then
1632	pa = log(p(nl))
1633	pb = log(pout(np))
1634	! special case: in case someone inadvertently enter p=0.
1635	if (p(nl+1).gt.0.d0) then
1636	pc = log(p(nl+1))
1637	else
1638	pc = log(1.d-4)
1639	end if
1640
1641	slope = (x(nl)-x(nl+1))/ (pa-pc)
1642	xout(np) = x(nl+1) + slope* (pb-pc)
1643	end if
1644	end do
1645	end do
1646	end if
1647
1648
1649	! place results in the return array;
1650	xxout = xout
1651
1652	END SUBROUTINE int1D
1653
1654	FUNCTION virtual (tmp,rmix)
1655	! This function returns virtual temperature in K, given temperature
1656	! in K and mixing ratio in kg/kg.
1657
1658	REAL(r_k) :: tmp, rmix, virtual
1659
1660	virtual=tmp(0.622+rmix)/(0.622(1.+rmix))
1661
1662	END FUNCTION virtual
1663
1664	END MODULE module_ForInterpolate

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

Download in other formats:

Original Format