Context Navigation

module_ForInterpolate.F90 @ 1554

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

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

Download in other formats:

Original Format