Context Navigation

module_ForInterpolate.F90 @ 1337

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

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

Download in other formats:

Original Format