Context Navigation

module_ForInterpolate.F90 @ 1030

Last change on this file since 1030 was 900, checked in by lfita, 8 years ago
Adding `pinterp': Function to vertically interpolate using subroutines from the p_interp.F90 NCAR program
File size: 44.9 KB

Line
1	! Module to interpolate values from a giving projection and pressure interpolation
2	! To be included in a python
3	! f2py -m module_ForInterpolate --f90exec=/usr/bin/gfortran-4.7 -c module_generic.F90 module_ForInterpolate.F90 >& run_f2py.log
4	MODULE module_ForInterpolate
5
6	CONTAINS
7
8	SUBROUTINE CoarselonlatFind(dx, dy, ilon, ilat, nxlon, nxlat, fraclon, fraclat, lonv, latv, per, &
9	Nperx, Npery, ilonlat, mindiffLl)
10	! Function to search a given value from a coarser version of the data
11
12	USE module_generic
13
14	IMPLICIT NONE
15
16	INTEGER, PARAMETER :: r_k = KIND(1.d0)
17	INTEGER, INTENT(in) :: dx, dy
18	REAL(r_k), DIMENSION(dx,dy), INTENT(in) :: ilon, ilat
19	REAL(r_k), DIMENSION(Nperx,Npery), INTENT(in) :: fraclon, fraclat
20	REAL(r_k), INTENT(in) :: lonv, latv, per
21	REAL(r_k), DIMENSION(2), INTENT(in) :: nxlon, nxlat
22	INTEGER, INTENT(in) :: Nperx, Npery
23	INTEGER, DIMENSION(2), INTENT(out) :: ilonlat
24	REAL(r_k), INTENT(out) :: mindiffLl
25	! Local
26	REAL(r_k), DIMENSION(Nperx,Npery) :: difffraclonlat
27	REAL(r_k) :: mindifffracLl
28	INTEGER, DIMENSION(2) :: ilonlatfrac
29	INTEGER :: ixbeg, ixend, iybeg, iyend
30	INTEGER :: fracx, fracy
31	REAL(r_k) :: fraclonv, fraclatv
32	REAL(r_k), ALLOCATABLE, DIMENSION(:,:) :: difflonlat, lon, lat
33	CHARACTER(LEN=50) :: fname
34
35	! Variables
36	! ilon, ilat: original 2D matrices with the longitudes and the latitudes
37	! lonv, latv: longitude and latitude to find
38	! nxlon, nxlat: minimum and maximum longitude and latitude of the target lon,lat
39	! per: fraction of the whole domain (as percentage)
40	! Nper[x/y]: period (as fraction over 1) of the fractions of the original grid to use to explore
41	! fraclon, fraclat: longitude and latitude fractional matricies to perform the first guess
42
43	fname = 'CoarselonlatFind'
44
45	IF (lonv < nxlon(1) .OR. lonv > nxlon(2)) THEN
46	PRINT *, TRIM(ErrWarnMsg('err'))
47	PRINT *,' ' // TRIM(fname) // ': longitude outside data range!!'
48	PRINT *,' given value:', lonv,' outside (',nxlon(1),' ,',nxlon(2),' )'
49	STOP
50	END IF
51	IF (latv < nxlat(1) .OR. latv > nxlat(2)) THEN
52	PRINT *, TRIM(ErrWarnMsg('err'))
53	PRINT *,' ' // TRIM(fname) // ': latitude outside data range!!'
54	PRINT *,' given value:', latv,' outside (',nxlat(1),' ,',nxlat(2),' )'
55	STOP
56	END IF
57
58	! Initializing variables
59	ixbeg = 0
60	ixend = 0
61	iybeg = 0
62	iyend = 0
63
64	fracx = int(dx*per)
65	fracy = int(dy*per)
66
67	! PRINT *,'fraclon _______'
68	! PRINT *,fraclon
69
70	! PRINT *,'fraclat _______'
71	! PRINT *,fraclat
72
73	! Fraction point
74	difffraclonlat = SQRT((fraclon-lonv)2. + (fraclat-latv)2.)
75	mindifffracLl = MINVAL(difffraclonlat)
76	ilonlatfrac = index2DArrayR(difffraclonlat, Nperx, Npery, mindifffracLl)
77
78	! PRINT *, 'mindifffracLl:', mindifffracLl, ' ilonlatfrac:', ilonlatfrac
79	! PRINT *, 'frac lon, lat:', fraclon(ilonlatfrac(1),ilonlatfrac(2)), ' ,', &
80	! fraclat(ilonlatfrac(1),ilonlatfrac(2))
81	! PRINT *, 'values lon, lat:', lonv, latv
82
83	! Providing fraction range
84	fraclonv = fraclon(ilonlatfrac(1),ilonlatfrac(2))
85	fraclatv = fraclat(ilonlatfrac(1),ilonlatfrac(2))
86
87	IF (fraclonv >= lonv .AND. fraclatv >= latv) THEN
88	IF (ilonlatfrac(1) > 0) THEN
89	ixbeg = (ilonlatfrac(1)-1)*fracx
90	ixend = ilonlatfrac(1)*fracx+1
91	ELSE
92	ixbeg = 0
93	ixend = fracx+1
94	END IF
95	IF (ilonlatfrac(2) > 0) THEN
96	iybeg = (ilonlatfrac(2)-1)*fracy
97	iyend = ilonlatfrac(2)*fracy+1
98	ELSE
99	iybeg = 0
100	iyend = fracy+1
101	END IF
102	ELSE IF (fraclonv < lonv .AND. fraclatv >= latv) THEN
103	IF (ilonlatfrac(1) < Nperx) THEN
104	IF (ilonlatfrac(1) /= 0) THEN
105	ixbeg = (ilonlatfrac(1)-1)*fracx
106	ixend = ilonlatfrac(1)*fracx+1
107	ELSE
108	ixbeg = 0
109	ixend = fracx+1
110	END IF
111	ELSE
112	ixbeg = Nperx*fracx
113	ixend = dx+1
114	END IF
115	IF (ilonlatfrac(2) > 0) THEN
116	iybeg = (ilonlatfrac(2)-1)*fracy
117	iyend = ilonlatfrac(2)*fracy+1
118	ELSE
119	iybeg = 0
120	iyend = fracy+1
121	END IF
122	ELSE IF (fraclonv < lonv .AND. fraclatv < latv) THEN
123	IF (ilonlatfrac(1) < Nperx) THEN
124	IF (ilonlatfrac(1) /= 0) THEN
125	ixbeg = (ilonlatfrac(1)-1)*fracx
126	ixend = ilonlatfrac(1)*fracx+1
127	ELSE
128	ixbeg = 0
129	ixend = fracx+1
130	END IF
131	ELSE
132	ixbeg = Nperx*fracx
133	ixend = dx+1
134	ENDIF
135	IF (ilonlatfrac(2) < Npery) THEN
136	IF (ilonlatfrac(2) /= 0) THEN
137	iybeg = (ilonlatfrac(2)-1)*fracy
138	iyend = ilonlatfrac(2)*fracy+1
139	ELSE
140	iybeg = 0
141	iyend = fracy+1
142	END IF
143	ELSE
144	iybeg = Npery*fracy
145	iyend = dy+1
146	END IF
147	ELSE IF (fraclonv >= lonv .AND. fraclatv < latv) THEN
148	IF (ilonlatfrac(1) > 0) THEN
149	ixbeg = (ilonlatfrac(1)-1)*fracx
150	ixend = ilonlatfrac(1)*fracx+1
151	ELSE
152	ixbeg = 0
153	ixend = fracx+1
154	END IF
155	IF (ilonlatfrac(2) < Npery) THEN
156	IF (ilonlatfrac(2) /= 0) THEN
157	iybeg = (ilonlatfrac(2)-1)*fracy
158	iyend = ilonlatfrac(2)*fracy+1
159	ELSE
160	iybeg = 0
161	iyend = fracy+1
162	END IF
163	ELSE
164	iybeg = Npery*fracy
165	iyend = dy+1
166	END IF
167	END IF
168
169	IF (ALLOCATED(lon)) DEALLOCATE(lon)
170	ALLOCATE(lon(ixend-ixbeg+1, iyend-iybeg+1))
171	IF (ALLOCATED(lat)) DEALLOCATE(lat)
172	ALLOCATE(lat(ixend-ixbeg+1, iyend-iybeg+1))
173	IF (ALLOCATED(difflonlat)) DEALLOCATE(difflonlat)
174	ALLOCATE(difflonlat(ixend-ixbeg+1, iyend-iybeg+1))
175
176	lon = ilon(ixbeg:ixend,iybeg:iyend)
177	lat = ilat(ixbeg:ixend,iybeg:iyend)
178
179	! print *,'lon _______'
180	! print *,lon
181	! print *,'lat _______'
182	! print *,lat
183
184	! Find point
185	difflonlat = SQRT((lon-lonv)2. + (lat-latv)2.)
186	mindiffLl = MINVAL(difflonlat)
187	ilonlat = index2DArrayR(difflonlat, ixend-ixbeg+1, iyend-iybeg+1, mindiffLl)
188
189	ilonlat(1) = ilonlat(1) + ixbeg
190	ilonlat(2) = ilonlat(2) + iybeg
191
192	! PRINT *,'mindiffLl:', mindiffLl, ' ilatlon:', ilatlon
193	! PRINT *,'lon, lat:', lon(ilonlat(1),ilonlat(2)), ' ,', lat(ilonlat(1),ilonlat(2))
194
195	RETURN
196
197	END SUBROUTINE CoarselonlatFind
198
199	SUBROUTINE CoarselonlatFindExact(dx, dy, ilon, ilat, nxlon, nxlat, fracx, fracy, fraclon, fraclat, &
200	iv, lonv, latv, per, Nperx, Npery, mindiff, ilonlat, mindiffLl)
201	! Function to search a given value from a coarser version of the data
202
203	USE module_generic
204
205	IMPLICIT NONE
206
207	INTEGER, PARAMETER :: r_k = KIND(1.d0)
208	INTEGER, INTENT(in) :: dx, dy, iv
209	REAL(r_k), DIMENSION(dx,dy), INTENT(in) :: ilon, ilat
210	INTEGER, INTENT(in) :: fracx, fracy
211	REAL(r_k), DIMENSION(Nperx,Npery), INTENT(in) :: fraclon, fraclat
212	REAL(r_k), INTENT(in) :: lonv, latv, per, mindiff
213	REAL(r_k), DIMENSION(2), INTENT(in) :: nxlon, nxlat
214	INTEGER, INTENT(in) :: Nperx, Npery
215	INTEGER, DIMENSION(2), INTENT(out) :: ilonlat
216	REAL(r_k), INTENT(out) :: mindiffLl
217	! Local
218	INTEGER :: i
219	REAL(r_k), DIMENSION(Nperx,Npery) :: difffraclonlat
220	REAL(r_k) :: mindifffracLl
221	INTEGER, DIMENSION(2) :: ilonlatfrac
222	INTEGER :: ixbeg, ixend, iybeg, iyend
223	REAL(r_k) :: fraclonv, fraclatv
224	REAL(r_k), ALLOCATABLE, DIMENSION(:,:) :: difflonlat, lon, lat
225	CHARACTER(LEN=50) :: fname
226
227	! Variables
228	! ilon, ilat: original 2D matrices with the longitudes and the latitudes
229	! lonv, latv: longitude and latitude to find
230	! iv: point in the input data
231	! nxlon, nxlat: minimum and maximum longitude and latitude of the target lon,lat
232	! per: fraction of the whole domain (as percentage)
233	! Nper[x/y]: period (as fraction over 1) of the fractions of the original grid to use to explore
234	! frac[x/y]: Number of grid points for each fraction
235	! fraclon, fraclat: longitude and latitude fractional matricies to perform the first guess
236	! mindiff: authorized minimal distance between input and interpolated point
237	! ilonlat: grid point on the total lon,lat matrix
238	! mindiffLl: distance between input and interpolated point
239
240	fname = 'CoarselonlatFindExact'
241
242	IF (lonv < nxlon(1) .OR. lonv > nxlon(2)) THEN
243	PRINT *, TRIM(ErrWarnMsg('err'))
244	PRINT *,' ' // TRIM(fname) // ': longitude outside data range!!'
245	PRINT *,' given value:', lonv,' outside (',nxlon(1),' ,',nxlon(2),' )'
246	STOP
247	END IF
248	IF (latv < nxlat(1) .OR. latv > nxlat(2)) THEN
249	PRINT *, TRIM(ErrWarnMsg('err'))
250	PRINT *,' ' // TRIM(fname) // ': latitude outside data range!!'
251	PRINT *,' given value:', latv,' outside (',nxlat(1),' ,',nxlat(2),' )'
252	STOP
253	END IF
254
255	! Initializing variables
256	ixbeg = 0
257	ixend = 0
258	iybeg = 0
259	iyend = 0
260
261	! Fraction point
262	difffraclonlat = SQRT((fraclon-lonv)2. + (fraclat-latv)2.)
263	mindifffracLl = MINVAL(difffraclonlat)
264	ilonlatfrac = index2DArrayR(difffraclonlat, Nperx, Npery, mindifffracLl)
265
266	! PRINT *, 'mindifffracLl:', mindifffracLl, ' ilonlatfrac:', ilonlatfrac
267	! PRINT *, 'frac lon, lat:', fraclon(ilonlatfrac(1),ilonlatfrac(2)), ' ,', &
268	! fraclat(ilonlatfrac(1),ilonlatfrac(2))
269	! PRINT *, 'values lon, lat:', lonv, latv
270
271	! Providing fraction range
272	fraclonv = fraclon(ilonlatfrac(1),ilonlatfrac(2))
273	fraclatv = fraclat(ilonlatfrac(1),ilonlatfrac(2))
274
275	IF (fraclonv >= lonv .AND. fraclatv >= latv) THEN
276	PRINT *,'Lluis!',fraclonv, '>=', lonv,'&', fraclatv, '>=', latv
277	IF (ilonlatfrac(1) > 1) THEN
278	ixbeg = (ilonlatfrac(1)-1)*fracx
279	ixend = ilonlatfrac(1)*fracx+1
280	ELSE
281	PRINT *,'Lluis 2'
282	ixbeg = 1
283	ixend = fracx+1
284	END IF
285	IF (ilonlatfrac(2) > 1) THEN
286	iybeg = (ilonlatfrac(2)-1)*fracy
287	iyend = ilonlatfrac(2)*fracy+1
288	ELSE
289	iybeg = 1
290	iyend = fracy+1
291	END IF
292	ELSE IF (fraclonv < lonv .AND. fraclatv >= latv) THEN
293	PRINT *,'Lluis!',fraclonv, '<', lonv,'&', fraclatv, '>=', latv
294	IF (ilonlatfrac(1) < Nperx) THEN
295	PRINT *,'Lluis 2'
296	IF (ilonlatfrac(1) /= 1) THEN
297	ixbeg = (ilonlatfrac(1)-1)*fracx
298	ixend = ilonlatfrac(1)*fracx+1
299	ELSE
300	ixbeg = 1
301	ixend = fracx+1
302	END IF
303	ELSE
304	ixbeg = Nperx*fracx
305	ixend = dx+1
306	END IF
307	IF (ilonlatfrac(2) > 1) THEN
308	iybeg = (ilonlatfrac(2)-1)*fracy
309	iyend = ilonlatfrac(2)*fracy+1
310	ELSE
311	iybeg = 1
312	iyend = fracy+1
313	END IF
314	ELSE IF (fraclonv < lonv .AND. fraclatv < latv) THEN
315	PRINT *,'Lluis!',fraclonv, '<', lonv,'&', fraclatv, '<', latv
316	IF (ilonlatfrac(1) < Nperx) THEN
317	IF (ilonlatfrac(1) /= 1) THEN
318	ixbeg = (ilonlatfrac(1)-1)*fracx
319	ixend = ilonlatfrac(1)*fracx+1
320	ELSE
321	ixbeg = 1
322	ixend = fracx+1
323	END IF
324	ELSE
325	ixbeg = Nperx*fracx
326	ixend = dx+1
327	ENDIF
328	IF (ilonlatfrac(2) < Npery) THEN
329	IF (ilonlatfrac(2) /= 1) THEN
330	iybeg = (ilonlatfrac(2)-1)*fracy
331	iyend = ilonlatfrac(2)*fracy+1
332	ELSE
333	iybeg = 1
334	iyend = fracy+1
335	END IF
336	ELSE
337	iybeg = Npery*fracy
338	iyend = dy+1
339	END IF
340	ELSE IF (fraclonv >= lonv .AND. fraclatv < latv) THEN
341	PRINT *,'Llui!',fraclonv, '>=', lonv,'&', fraclatv, '<', latv
342	IF (ilonlatfrac(1) > 1) THEN
343	ixbeg = (ilonlatfrac(1)-1)*fracx
344	ixend = ilonlatfrac(1)*fracx+1
345	ELSE
346	ixbeg = 1
347	ixend = fracx+1
348	END IF
349	IF (ilonlatfrac(2) < Npery) THEN
350	IF (ilonlatfrac(2) /= 1) THEN
351	iybeg = (ilonlatfrac(2)-1)*fracy
352	iyend = ilonlatfrac(2)*fracy+1
353	ELSE
354	iybeg = 1
355	iyend = fracy+1
356	END IF
357	ELSE
358	iybeg = Npery*fracy
359	iyend = dy+1
360	END IF
361	END IF
362
363	IF (ALLOCATED(lon)) DEALLOCATE(lon)
364	ALLOCATE(lon(ixend-ixbeg+1, iyend-iybeg+1))
365	IF (ALLOCATED(lat)) DEALLOCATE(lat)
366	ALLOCATE(lat(ixend-ixbeg+1, iyend-iybeg+1))
367	IF (ALLOCATED(difflonlat)) DEALLOCATE(difflonlat)
368	ALLOCATE(difflonlat(ixend-ixbeg+1, iyend-iybeg+1))
369
370	lon = ilon(ixbeg:ixend,iybeg:iyend)
371	lat = ilat(ixbeg:ixend,iybeg:iyend)
372
373	! print *,'lon _______'
374	! print *,lon
375	! print *,'lat _______'
376	! print *,lat
377
378	! Find point
379	difflonlat = SQRT((lon-lonv)2. + (lat-latv)2.)
380	mindiffLl = MINVAL(difflonlat)
381
382	IF (mindiffLl > mindiff) THEN
383	difflonlat = SQRT((lon-lonv)2. + (lat-latv)2.)
384	mindiffLl = MINVAL(difflonlat)
385	END IF
386
387	IF (mindiffLl > mindiff) THEN
388	PRINT *,TRIM(ErrWarnMsg('err'))
389	PRINT *,' ' // TRIM(fname) // ': not equivalent point closer than:',mindiff,' found!!'
390	PRINT *,' at input point iv:', iv,' lon/lat:', lonv,', ',latv,' distance:',mindiffLl
391	PRINT *,' Fraction values _______ (',Nperx,', ',Npery ,')'
392	PRINT *,' fraclon'
393	DO i=1, Nperx
394	PRINT *,' ',fraclon(i,:)
395	END DO
396	PRINT *,' fraclat'
397	DO i=1, Nperx
398	PRINT *,' ',fraclat(i,:)
399	END DO
400	PRINT *,' frac lon, lat:', fraclon(ilonlatfrac(1),ilonlatfrac(2)), ' ,', &
401	fraclat(ilonlatfrac(1),ilonlatfrac(2))
402	PRINT *,' mindifffracLl:', mindifffracLl, ' ilonlatfrac:', ilonlatfrac
403	PRINT *,' Coarse values _______'
404	PRINT *,' indices. x:', ixbeg, ', ', ixend, ' y:', iybeg, ', ', iyend
405	PRINT *,' lon range:', '(',ilon(ixbeg,iybeg),', ',ilon(ixend,iyend),')'
406	PRINT *,' lat range:', '(',ilat(ixbeg,iybeg),', ',ilat(ixend,iyend),')'
407	PRINT *,' lon', UBOUND(lon)
408	DO i=1, ixend-ixbeg+1
409	PRINT *,' ',lon(i,:)
410	END DO
411	PRINT *,' lat', UBOUND(lat)
412	DO i=1, ixend-ixbeg+1
413	PRINT *,' ',lat(i,:)
414	END DO
415	STOP
416	END IF
417
418	ilonlat = index2DArrayR(difflonlat, ixend-ixbeg+1, iyend-iybeg+1, mindiffLl)
419
420	ilonlat(1) = ilonlat(1) + ixbeg
421	ilonlat(2) = ilonlat(2) + iybeg
422
423	! PRINT *,'mindiffLl:', mindiffLl, ' ilatlon:', ilatlon
424	! PRINT *,'lon, lat:', lon(ilonlat(1),ilonlat(2)), ' ,', lat(ilonlat(1),ilonlat(2))
425
426	RETURN
427
428	END SUBROUTINE CoarselonlatFindExact
429
430	SUBROUTINE lonlatFind(dx, dy, ilon, ilat, nxlon, nxlat, lonv, latv, ilonlat, mindiffLl)
431	! Function to search a given value from a coarser version of the data
432
433	USE module_generic
434
435	IMPLICIT NONE
436
437	INTEGER, PARAMETER :: r_k = KIND(1.d0)
438	INTEGER, INTENT(in) :: dx, dy
439	REAL(r_k), DIMENSION(dx,dy), INTENT(in) :: ilon, ilat
440	REAL(r_k), INTENT(in) :: lonv, latv
441	REAL(r_k), DIMENSION(2), INTENT(in) :: nxlon, nxlat
442	INTEGER, DIMENSION(2), INTENT(out) :: ilonlat
443	REAL(r_k), INTENT(out) :: mindiffLl
444	! Local
445	REAL(r_k), DIMENSION(dx,dy) :: difflonlat
446	CHARACTER(LEN=50) :: fname
447
448	! Variables
449	! ilon, ilat: original 2D matrices with the longitudes and the latitudes
450	! lonv, latv: longitude and latitude to find
451	! nxlon, nxlat: minimum and maximum longitude and latitude of the target lon,lat
452
453	fname = 'lonlatFind'
454
455	IF (lonv < nxlon(1) .OR. lonv > nxlon(2)) THEN
456	PRINT *, TRIM(ErrWarnMsg('err'))
457	PRINT *,' ' // TRIM(fname) // ': longitude outside data range!!'
458	PRINT *,' given value:', lonv,' outside (',nxlon(1),' ,',nxlon(2),' )'
459	STOP
460	END IF
461	IF (latv < nxlat(1) .OR. latv > nxlat(2)) THEN
462	PRINT *, TRIM(ErrWarnMsg('err'))
463	PRINT *,' ' // TRIM(fname) // ': latitude outside data range!!'
464	PRINT *,' given value:', latv,' outside (',nxlat(1),' ,',nxlat(2),' )'
465	STOP
466	END IF
467
468	! Find point
469	difflonlat = SQRT((ilon-lonv)2. + (ilat-latv)2.)
470	mindiffLl = MINVAL(difflonlat)
471	ilonlat = index2DArrayR(difflonlat, dx, dy, mindiffLl)
472
473	! PRINT *,'mindiffLl:', mindiffLl, ' ilatlon:', ilatlon
474	! PRINT *,'lon, lat:', lon(ilonlat(1),ilonlat(2)), ' ,', lat(ilonlat(1),ilonlat(2))
475
476	RETURN
477
478	END SUBROUTINE lonlatFind
479
480	SUBROUTINE CoarseInterpolate(projlon, projlat, lonvs, latvs, percen, mindiff, inpt, ilonlat, &
481	mindiffLl, dimx, dimy, Ninpts)
482	! Subroutine which finds the closest grid point within a projection throughout a first guest
483	! approche from percentages of the whole domain
484
485	USE module_generic
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	USE module_generic
590
591	IMPLICIT NONE
592
593	INTEGER, PARAMETER :: r_k = KIND(1.d0)
594	INTEGER, INTENT(in) :: dimx, dimy
595	REAL(r_k), DIMENSION(dimx,dimy), INTENT(in) :: projlon, projlat
596	INTEGER, INTENT(in) :: Ninpts
597	REAL(r_k), DIMENSION(Ninpts), INTENT(in) :: ivar, lonvs, latvs
598	REAL(r_k), INTENT(in) :: mindiff, percen
599	REAL(r_k), DIMENSION(dimx,dimy), INTENT(out) :: newvar
600	INTEGER, DIMENSION(dimx,dimy), INTENT(out) :: newvarin
601	INTEGER, DIMENSION(Ninpts), INTENT(out) :: newvarinpt
602	REAL(r_k), DIMENSION(Ninpts), INTENT(out) :: newvarindiff
603
604	! Local
605	INTEGER :: iv,i,j
606	INTEGER :: ierr
607	INTEGER, DIMENSION(2) :: ilonlat
608	REAL(r_k) :: mindiffLl
609	INTEGER :: Ninpts1
610	REAL(r_k), DIMENSION(2) :: extremelon, extremelat
611	REAL(r_k), ALLOCATABLE, DIMENSION(:,:) :: fractionlon, fractionlat
612	INTEGER :: dfracdx, dfracdy, fracdx, fracdy
613	CHARACTER(LEN=50) :: fname
614
615	!!!!!!! Variables
616	! dimx, dimy: dimension length of the target interpolation
617	! proj[lon/lat]: longitudes and latitudes of the target interpolation
618	! Ninpts: number of points to interpolate
619	! [lon/lat]vs: longitudes and latitudes of the points to interpolate
620	! mindiff: minimal accepted distance to the target point
621	! percen: size (as percentage of the total domain) of the first guess portions to provide the first gues
622	! ivar: values to localize in the target projection
623	! newvar: localisation of the [lon/lat]vs point in the target projection
624	! newvarin: number of point from the input data
625	! newvarinpt: integer value indicating if the value has been already located (0: no, 1: yes)
626	! newvarindiff: distance of point from the input data to the closest target point
627	! ncid: netCDF output file id
628
629	fname = 'CoarseInterpolateExact'
630	Ninpts1 = Ninpts/100
631
632	extremelon = (/ MINVAL(projlon), MAXVAL(projlon) /)
633	extremelat = (/ MINVAL(projlat), MAXVAL(projlat) /)
634
635	PRINT *,' ' // TRIM(fname) //' total space:', dimx, ', ', dimy, ' %', percen
636
637	dfracdx = INT(1./percen+1)
638	dfracdy = INT(1./percen+1)
639	fracdx = INT(dimx*percen)
640	fracdy = INT(dimy*percen)
641	PRINT *,' ' // TRIM(fname) //' fraction:', dfracdx, ', ', dfracdy, ' freq:', fracdx,', ',fracdy
642
643	IF (ALLOCATED(fractionlon)) DEALLOCATE(fractionlon)
644	ALLOCATE(fractionlon(dfracdx, dfracdy), STAT=ierr)
645	IF (ierr /= 0) THEN
646	PRINT *,TRIM(ErrWarnMsg('err'))
647	PRINT *,' ' // TRIM(fname) //": problem allocating 'fractionlon' !!"
648	STOP
649	END IF
650	IF (ALLOCATED(fractionlat)) DEALLOCATE(fractionlat)
651	ALLOCATE(fractionlat(dfracdx, dfracdy), STAT=ierr)
652	IF (ierr /= 0) THEN
653	PRINT *,TRIM(ErrWarnMsg('err'))
654	PRINT *,' ' // TRIM(fname) //": problem allocating 'fractionlat' !!"
655	STOP
656	END IF
657
658	DO i=1,dfracdx
659	DO j=1,dfracdy
660	fractionlon(i,j) = projlon(fracdx(i-1)+1,fracdy(j-1)+1)
661	fractionlat(i,j) = projlat(fracdx(i-1)+1,fracdy(j-1)+1)
662	! PRINT ,'i,j:',i,', ',j,' frac ij:',fracdx(i-1),', ',fracdy*(j-1),' lonlat:', fractionlon(i,j),&
663	! ', ',fractionlat(i,j)
664	END DO
665	END DO
666
667	! PRINT *,' ' // TRIM(fname) // ' fractions of:'
668	! PRINT *,' lon _______ (',dfracdx,', ',dfracdy,')'
669	! DO i=1,dfracdx
670	! PRINT *,fractionlon(i,:)
671	! END DO
672	! PRINT *,' lat_______'
673	! DO i=1,dfracdx
674	! PRINT *,fractionlat(i,:)
675	! END DO
676
677	DO iv=1,Ninpts
678	IF (newvarinpt(iv) == 0) THEN
679	CALL CoarselonlatFindExact(dimx, dimy, projlon, projlat, extremelon, extremelat, fracdx, fracdy,&
680	fractionlon, fractionlat, iv, lonvs(iv), latvs(iv), percen, dfracdx, dfracdy, mindiff, &
681	ilonlat, mindiffLl)
682
683	IF (mindiffLl >= mindiff) THEN
684	! percendone(iv,Ninpts,0.5,'done:')
685
686	IF (ilonlat(1) >= 0 .AND. ilonlat(1) >= 0) THEN
687	newvar(ilonlat(1),ilonlat(2)) = ivar(iv)
688	newvarin(ilonlat(1),ilonlat(2)) = iv
689	newvarinpt(iv) = 1
690	newvarindiff(iv) = mindiffLl
691	! PRINT *,'Lluis iv:', newvarin(ilonlat(1),ilonlat(2)), ' localized:', newvarinpt(iv), &
692	! ' values:', newvar(ilonlat(1),ilonlat(2)), ' invalues:', ivar(iv), ' mindist:', &
693	! newvarindiff(iv), ' point:',ilonlat
694	ELSE
695	PRINT *,TRIM(ErrWarnMsg('err'))
696	PRINT *,' ' // TRIM(fname) // ': point iv:', iv, ' at', lonvs(iv), ' ,', latvs(iv), &
697	' not relocated !!'
698	PRINT *,' mindiffl:', mindiffLl, ' ilon:', ilonlat(1), ' ilat:', ilonlat(2)
699	STOP
700	END IF
701
702	! IF (MOD(iv,Ninpts1) == 0) newnc.sync()
703	ELSE
704	PRINT *,TRIM(ErrWarnMsg('err'))
705	PRINT *,' ' // TRIM(fname) // ': for point #', iv,' lon,lat in incomplet map:', lonvs(iv), &
706	' ,', latvs(iv), ' there is not a set of lon,lat in the completed map closer than: ', &
707	mindiff, ' !!'
708	PRINT *,' found minimum difference:', mindiffLl
709	STOP
710	END IF
711	END IF
712	END DO
713
714	END SUBROUTINE CoarseInterpolateExact
715
716	SUBROUTINE Interpolate(projlon, projlat, lonvs, latvs, mindiff, inpt, diffs, ilonlat, dimx, dimy, &
717	Ninpts)
718	! Subroutine which finds the closest grid point within a projection
719
720	USE module_generic
721
722	IMPLICIT NONE
723
724	INTEGER, PARAMETER :: r_k = KIND(1.d0)
725	INTEGER, INTENT(in) :: dimx, dimy
726	REAL(r_k), DIMENSION(dimx,dimy), INTENT(in) :: projlon, projlat
727	INTEGER, INTENT(in) :: Ninpts
728	REAL(r_k), DIMENSION(Ninpts), INTENT(in) :: lonvs, latvs
729	REAL(r_k), INTENT(in) :: mindiff
730	INTEGER, DIMENSION(Ninpts), INTENT(inout) :: inpt
731	REAL(r_k), DIMENSION(Ninpts), INTENT(out) :: diffs
732	INTEGER, DIMENSION(Ninpts,2), INTENT(out) :: ilonlat
733
734	! Local
735	INTEGER :: iv
736	REAL(r_k) :: mindiffLl
737	INTEGER :: Ninpts1
738	REAL(r_k), DIMENSION(dimx,dimy) :: difflonlat
739	REAL(r_k), DIMENSION(2) :: extremelon, extremelat
740	CHARACTER(LEN=50) :: fname
741
742	!!!!!!! Variables
743	! dimx, dimy: dimension length of the target interpolation
744	! proj[lon/lat]: longitudes and latitudes of the target interpolation
745	! Ninpts: number of points to interpolate
746	! [lon/lat]vs: longitudes and latitudes of the points to interpolate
747	! mindiff: minimal accepted distance to the target point
748	! inpt: whether the point has already been localized
749	! diffs: distance of point from the input data to the closest target point
750	! ilonlat: longitude and latitude of the point
751	! ncid: netCDF output file id
752
753	fname = 'Interpolate'
754	Ninpts1 = Ninpts/100
755
756	extremelon = (/ MINVAL(projlon), MAXVAL(projlon) /)
757	extremelat = (/ MINVAL(projlat), MAXVAL(projlat) /)
758
759	DO iv=1,Ninpts
760	IF (inpt(iv) <= 0) THEN
761	! Not using the subroutine, not efficient!
762	! CALL lonlatFind(dimx, dimy, projlon, projlat, extremelon, extremelat, lonvs(iv), latvs(iv), &
763	! ilonlat, mindiffLl)
764
765	IF (lonvs(iv) < extremelon(1) .OR. lonvs(iv) > extremelon(2)) THEN
766	PRINT *, TRIM(ErrWarnMsg('err'))
767	PRINT *,' ' // TRIM(fname) // ': longitude outside data range!!'
768	PRINT *,' given value:', lonvs(iv),' outside (',extremelon(1),' ,',extremelon(2),' )'
769	STOP
770	END IF
771	IF (latvs(iv) < extremelat(1) .OR. latvs(iv) > extremelat(2)) THEN
772	PRINT *, TRIM(ErrWarnMsg('err'))
773	PRINT *,' ' // TRIM(fname) // ': latitude outside data range!!'
774	PRINT *,' given value:', latvs(iv),' outside (',extremelat(1),' ,',extremelat(2),' )'
775	STOP
776	END IF
777
778	! Find point
779	difflonlat = SQRT((projlon-lonvs(iv))2. + (projlat-latvs(iv))2.)
780	mindiffLl = MINVAL(difflonlat)
781	ilonlat(iv,:) = index2DArrayR(difflonlat, dimx, dimy, mindiffLl)
782
783	IF (mindiffLl <= mindiff) THEN
784	! percendone(iv,Ninpts,0.5,'done:')
785
786	IF (ilonlat(iv,1) >= 0 .AND. ilonlat(iv,2) >= 0) THEN
787	diffs(iv) = mindiffLl
788	inpt(iv) = 1
789	! PRINT *,'Lluis iv:', newvarin(ilonlat(1),ilonlat(2)), ' localized:', newvarinpt(iv), &
790	! ' values:', newvar(ilonlat(1),ilonlat(2)), ' invalues:', ivar(iv), ' mindist:', &
791	! newvarindiff(iv), ' point:',ilonlat
792	ELSE
793	PRINT *,TRIM(ErrWarnMsg('err'))
794	PRINT *,' ' // TRIM(fname) // ': point iv:', iv, ' at', lonvs(iv), ' ,', latvs(iv), &
795	' not relocated !!'
796	PRINT *,' mindiffl:', mindiffLl, ' ilon:', ilonlat(iv,1), ' ilat:', ilonlat(iv,2)
797	STOP
798	END IF
799
800	! IF (MOD(iv,Ninpts1) == 0) newnc.sync()
801	ELSE
802	! Because doing boxes and Goode is not conitnuos, we should jump this error message
803	PRINT *,TRIM(ErrWarnMsg('err'))
804	PRINT *,' ' // TRIM(fname) // ': for point #', iv,' lon,lat in incomplet map:', lonvs(iv), &
805	' ,', latvs(iv), ' there is not a set of lon,lat in the completed map closer than: ', &
806	mindiff, ' !!'
807	PRINT *,' found minimum difference:', mindiffLl
808	STOP
809	END IF
810	END IF
811	END DO
812
813	END SUBROUTINE Interpolate
814
815	SUBROUTINE Interpolate1DLl(projlon, projlat, lonvs, latvs, mindiff, inpt, diffs, ilonlat, dimx, dimy, &
816	Ninpts)
817	! Subroutine which finds the closest grid point within a projection with 1D longitudes and latitudes
818
819	USE module_generic
820
821	IMPLICIT NONE
822
823	INTEGER, PARAMETER :: r_k = KIND(1.d0)
824	INTEGER, INTENT(in) :: dimx, dimy
825	REAL(r_k), DIMENSION(dimx), INTENT(in) :: projlon
826	REAL(r_k), DIMENSION(dimy), INTENT(in) :: projlat
827	INTEGER, INTENT(in) :: Ninpts
828	REAL(r_k), DIMENSION(Ninpts), INTENT(in) :: lonvs, latvs
829	REAL(r_k), INTENT(in) :: mindiff
830	INTEGER, DIMENSION(Ninpts), INTENT(inout) :: inpt
831	REAL(r_k), DIMENSION(Ninpts), INTENT(out) :: diffs
832	INTEGER, DIMENSION(Ninpts,2), INTENT(out) :: ilonlat
833
834	! Local
835	INTEGER :: iv
836	REAL(r_k) :: mindifflo, mindiffLa, mindiffLl
837	INTEGER :: Ninpts1
838	REAL(r_k), DIMENSION(dimx) :: difflon
839	REAL(r_k), DIMENSION(dimy) :: difflat
840	REAL(r_k), DIMENSION(2) :: extremelon, extremelat
841	CHARACTER(LEN=50) :: fname
842
843	!!!!!!! Variables
844	! dimx, dimy: dimension length of the target interpolation
845	! proj[lon/lat]: longitudes and latitudes of the target interpolation
846	! Ninpts: number of points to interpolate
847	! [lon/lat]vs: longitudes and latitudes of the points to interpolate
848	! mindiff: minimal accepted distance to the target point
849	! inpt: whether the point has already been localized
850	! diffs: distance of point from the input data to the closest target point
851	! ilonlat: longitude and latitude of the point
852	! ncid: netCDF output file id
853
854	fname = 'Interpolate1DLl'
855	Ninpts1 = Ninpts/100
856
857	extremelon = (/ MINVAL(projlon), MAXVAL(projlon) /)
858	extremelat = (/ MINVAL(projlat), MAXVAL(projlat) /)
859
860	DO iv=1,Ninpts
861	IF (inpt(iv) <= 0) THEN
862	! Not using the subroutine, not efficient!
863	! CALL lonlatFind(dimx, dimy, projlon, projlat, extremelon, extremelat, lonvs(iv), latvs(iv), &
864	! ilonlat, mindiffLl)
865
866	IF (lonvs(iv) < extremelon(1) .OR. lonvs(iv) > extremelon(2)) THEN
867	PRINT *, TRIM(ErrWarnMsg('err'))
868	PRINT *,' ' // TRIM(fname) // ': longitude outside data range!!'
869	PRINT *,' given value:', lonvs(iv),' outside (',extremelon(1),' ,',extremelon(2),' )'
870	STOP
871	END IF
872	IF (latvs(iv) < extremelat(1) .OR. latvs(iv) > extremelat(2)) THEN
873	PRINT *, TRIM(ErrWarnMsg('err'))
874	PRINT *,' ' // TRIM(fname) // ': latitude outside data range!!'
875	PRINT *,' given value:', latvs(iv),' outside (',extremelat(1),' ,',extremelat(2),' )'
876	STOP
877	END IF
878
879	! Find point
880	difflon = SQRT((projlon-lonvs(iv))**2.)
881	difflat = SQRT((projlat-latvs(iv))**2.)
882	mindifflo = MINVAL(difflon)
883	mindiffLa = MINVAL(difflat)
884	mindifflL = SQRT(mindifflomindifflo + mindiffLamindiffLa)
885	ilonlat(iv,1) = index1DArrayR(difflon, dimx, mindifflo)
886	ilonlat(iv,2) = index1DArrayR(difflat, dimy, mindiffLa)
887	! PRINT *,' Lluis: iv',iv,' lonvs:', lonvs(iv),' latvs:',latvs(iv)
888	! PRINT *,' Lluis: mindifflo:', mindifflo,' ilonlat(1):',ilonlat(iv,1)
889	! PRINT *,' Lluis: mindiffLa:', mindiffLa,' ilonlat(2):',ilonlat(iv,2)
890
891
892	IF (mindiffLl <= mindiff) THEN
893	! percendone(iv,Ninpts,0.5,'done:')
894
895	IF (ilonlat(iv,1) >= 1 .AND. ilonlat(iv,2) >= 1) THEN
896	diffs(iv) = mindiffLl
897	inpt(iv) = 1
898	! PRINT *,'Lluis iv:', newvarin(ilonlat(1),ilonlat(2)), ' localized:', newvarinpt(iv), &
899	! ' values:', newvar(ilonlat(1),ilonlat(2)), ' invalues:', ivar(iv), ' mindist:', &
900	! newvarindiff(iv), ' point:',ilonlat
901	ELSE
902	PRINT *,TRIM(ErrWarnMsg('err'))
903	PRINT *,' ' // TRIM(fname) // ': point iv:', iv, ' at', lonvs(iv), ' ,', latvs(iv), &
904	' not relocated !!'
905	PRINT *,' mindiffl:', mindiffLl, ' ilon:', ilonlat(iv,1), ' ilat:', ilonlat(iv,2)
906	STOP
907	END IF
908
909	! IF (MOD(iv,Ninpts1) == 0) newnc.sync()
910	ELSE
911	! Because doing boxes and Goode is not conitnuos, we should jump this error message
912	PRINT *,TRIM(ErrWarnMsg('err'))
913	PRINT *,' ' // TRIM(fname) // ': for point #', iv,' lon,lat in incomplet map:', lonvs(iv), &
914	' ,', latvs(iv), ' there is not a set of lon,lat in the completed map closer than: ', &
915	mindiff, ' !!'
916	PRINT *,' found minimum difference:', mindiffLl
917	STOP
918	END IF
919	END IF
920	END DO
921
922	END SUBROUTINE Interpolate1DLl
923
924
925	SUBROUTINE interp (data_in, pres_field, interp_levels, psfc, ter, tk, qv, LINLOG, extrapolate, &
926	GEOPT, MISSING, data_out, ix, iy, iz, it, num_metgrid_levels)
927	! Interpolation subroutine from the p_interp.F90 NCAR program
928	! Program to read wrfout data and interpolate to pressure levels
929	! The program reads namelist.pinterp
930	! November 2007 - Cindy Bruyere
931	!
932	INTEGER, INTENT(IN) :: ix, iy, iz, it
933	INTEGER, INTENT(IN) :: num_metgrid_levels, LINLOG
934	REAL, DIMENSION(ix,iy,iz,it), INTENT(IN) :: data_in, pres_field, tk, qv
935	REAL, DIMENSION(ix,iy,it), INTENT(IN) :: psfc
936	REAL, DIMENSION(ix,iy), INTENT(IN) :: ter
937	REAL, DIMENSION(num_metgrid_levels), INTENT(IN) :: interp_levels
938	INTEGER, INTENT(IN) :: extrapolate
939	REAL, INTENT(IN) :: MISSING
940	LOGICAL, INTENT(IN) :: GEOPT
941	REAL, DIMENSION(ix,iy,num_metgrid_levels,it), &
942	INTENT(OUT) :: data_out
943
944	! Local
945	INTEGER :: i, j, itt, k, kk, kin
946	REAL, DIMENSION(num_metgrid_levels) :: data_out1D
947	REAL, DIMENSION(iz) :: data_in1D, pres_field1D
948	REAL, DIMENSION(ix, iy, num_metgrid_levels, it) :: N
949	REAL :: sumA, sumN, AVE_geopt
950
951	!!!!!!! Variables
952	! data_out: interpolated field
953	! data_in: field to interpolate
954	! pres_field: pressure field [Pa]
955	! interp_levels: pressure levels to interpolate [hPa]
956	! psfc: surface pressure [Pa]
957	! ter: terrein height [m]
958	! tk: temperature [K]
959	! qv: mositure mizing ratio [kg/kg]
960	! i[x/y/z/t]: size of the matrices
961	! num_metgrid_levels: number of pressure values to interpolate
962	! LINLOG: if abs(linlog)=1 use linear interp in pressure
963	! if abs(linlog)=2 linear interp in ln(pressure)
964	! extrapolate: whether to set to missing value below/above model ground and top (0), or extrapolate (1)
965	! GEOPT: Wether the file is the geopotential file or not
966	! MISSING: Missing value
967
968	N = 1.0
969
970	expon=287.04*.0065/9.81
971
972	do itt = 1, it
973	do j = 1, iy
974	do i = 1, ix
975	data_in1D(:) = data_in(i,j,:,itt)
976	pres_field1D(:) = pres_field(i,j,:,itt)
977	CALL int1D (data_out1D, data_in1D, pres_field1D, interp_levels, iz, num_metgrid_levels, LINLOG, MISSING)
978	data_out(i,j,:,itt) = data_out1D(:)
979	end do
980	end do
981	end do
982
983
984	! Fill in missing values
985	IF ( extrapolate == 0 ) RETURN !! no extrapolation - we are out of here
986
987	! First find where about 400 hPa is located
988	kk = 0
989	find_kk : do k = 1, num_metgrid_levels
990	kk = k
991	if ( interp_levels(k) <= 40000. ) exit find_kk
992	end do find_kk
993
994
995	IF ( GEOPT ) THEN !! geopt is treated different below ground
996
997	do itt = 1, it
998	do k = 1, kk
999	do j = 1, iy
1000	do i = 1, ix
1001	IF ( data_out(i,j,k,itt) == MISSING .AND. interp_levels(k) < psfc(i,j,itt) ) THEN
1002
1003	! We are below the first model level, but above the ground
1004
1005	data_out(i,j,k,itt) = ((interp_levels(k) - pres_field(i,j,1,itt))ter(i,j)9.81 + &
1006	(psfc(i,j,itt) - interp_levels(k))*data_in(i,j,1,itt) ) / &
1007	(psfc(i,j,itt) - pres_field(i,j,1,itt))
1008
1009	ELSEIF ( data_out(i,j,k,itt) == MISSING ) THEN
1010
1011	! We are below both the ground and the lowest data level.
1012
1013	! First, find the model level that is closest to a "target" pressure
1014	! level, where the "target" pressure is delta-p less that the local
1015	! value of a horizontally smoothed surface pressure field. We use
1016	! delta-p = 150 hPa here. A standard lapse rate temperature profile
1017	! passing through the temperature at this model level will be used
1018	! to define the temperature profile below ground. This is similar
1019	! to the Benjamin and Miller (1990) method, except that for
1020	! simplicity, they used 700 hPa everywhere for the "target" pressure.
1021	! Code similar to what is implemented in RIP4
1022
1023	ptarget = (psfc(i,j,itt)*.01) - 150.
1024	dpmin=1.e4
1025	kupper = 0
1026	loop_kIN : do kin=iz,1,-1
1027	kupper = kin
1028	dp=abs( (pres_field(i,j,kin,itt)*.01) - ptarget )
1029	if (dp.gt.dpmin) exit loop_kIN
1030	dpmin=min(dpmin,dp)
1031	enddo loop_kIN
1032
1033	pbot=max(pres_field(i,j,1,itt),psfc(i,j,itt))
1034	zbot=min(data_in(i,j,1,itt)/9.81,ter(i,j))
1035
1036	tbotextrap=tk(i,j,kupper,itt)(pbot/pres_field(i,j,kupper,itt))*expon
1037	tvbotextrap=virtual(tbotextrap,qv(i,j,1,itt))
1038
1039	data_out(i,j,k,itt) = (zbot+tvbotextrap/.0065(1.-(interp_levels(k)/pbot)expon))9.81
1040
1041	ENDIF
1042	enddo
1043	enddo
1044	enddo
1045	enddo
1046
1047
1048	!!! Code for filling missing data with an average - we don't want to do this
1049	!!do itt = 1, it
1050	!!loop_levels : do k = 1, num_metgrid_levels
1051	!!sumA = SUM(data_out(:,:,k,itt), MASK = data_out(:,:,k,itt) /= MISSING)
1052	!!sumN = SUM(N(:,:,k,itt), MASK = data_out(:,:,k,itt) /= MISSING)
1053	!!IF ( sumN == 0. ) CYCLE loop_levels
1054	!!AVE_geopt = sumA/sumN
1055	!!WHERE ( data_out(:,:,k,itt) == MISSING )
1056	!!data_out(:,:,k,itt) = AVE_geopt
1057	!!END WHERE
1058	!!end do loop_levels
1059	!!end do
1060
1061	END IF
1062
1063	!!! All other fields and geopt at higher levels come here
1064	do itt = 1, it
1065	do j = 1, iy
1066	do i = 1, ix
1067	do k = 1, kk
1068	if ( data_out(i,j,k,itt) == MISSING ) data_out(i,j,k,itt) = data_in(i,j,1,itt)
1069	end do
1070	do k = kk+1, num_metgrid_levels
1071	if ( data_out(i,j,k,itt) == MISSING ) data_out(i,j,k,itt) = data_in(i,j,iz,itt)
1072	end do
1073	end do
1074	end do
1075	end do
1076
1077	END SUBROUTINE interp
1078
1079	SUBROUTINE int1D(xxout, xxin, ppin, ppout, npin, npout, LINLOG, MISSING)
1080
1081	! Modified from int2p - NCL code
1082	! routine to interpolate from one set of pressure levels
1083	! . to another set using linear or ln(p) interpolation
1084	!
1085	! NCL: xout = int2p (pin,xin,pout,linlog)
1086	! This code was originally written for a specific purpose.
1087	! . Several features were added for incorporation into NCL's
1088	! . function suite including linear extrapolation.
1089	!
1090	! nomenclature:
1091	!
1092	! . ppin - input pressure levels. The pin can be
1093	! . be in ascending or descending order
1094	! . xxin - data at corresponding input pressure levels
1095	! . npin - number of input pressure levels >= 2
1096	! . ppout - output pressure levels (input by user)
1097	! . same (ascending or descending) order as pin
1098	! . xxout - data at corresponding output pressure levels
1099	! . npout - number of output pressure levels
1100	! . linlog - if abs(linlog)=1 use linear interp in pressure
1101	! . if abs(linlog)=2 linear interp in ln(pressure)
1102	! . missing- missing data code.
1103
1104	! ! input types
1105	INTEGER :: npin,npout,linlog,ier
1106	real :: ppin(npin),xxin(npin),ppout(npout)
1107	real :: MISSING
1108	logical :: AVERAGE
1109	! ! output
1110	real :: xxout(npout)
1111	INTEGER :: j1,np,nl,nin,nlmax,nplvl
1112	INTEGER :: nlsave,np1,no1,n1,n2,nlstrt
1113	real :: slope,pa,pb,pc
1114
1115	! automatic arrays
1116	real :: pin(npin),xin(npin),p(npin),x(npin)
1117	real :: pout(npout),xout(npout)
1118
1119
1120	xxout = MISSING
1121	pout = ppout
1122	p = ppin
1123	x = xxin
1124	nlmax = npin
1125
1126	! exact p-level matches
1127	nlstrt = 1
1128	nlsave = 1
1129	do np = 1,npout
1130	xout(np) = MISSING
1131	do nl = nlstrt,nlmax
1132	if (pout(np).eq.p(nl)) then
1133	xout(np) = x(nl)
1134	nlsave = nl + 1
1135	go to 10
1136	end if
1137	end do
1138	10 nlstrt = nlsave
1139	end do
1140
1141	if (LINLOG.eq.1) then
1142	do np = 1,npout
1143	do nl = 1,nlmax - 1
1144	if (pout(np).lt.p(nl) .and. pout(np).gt.p(nl+1)) then
1145	slope = (x(nl)-x(nl+1))/ (p(nl)-p(nl+1))
1146	xout(np) = x(nl+1) + slope* (pout(np)-p(nl+1))
1147	end if
1148	end do
1149	end do
1150	elseif (LINLOG.eq.2) then
1151	do np = 1,npout
1152	do nl = 1,nlmax - 1
1153	if (pout(np).lt.p(nl) .and. pout(np).gt.p(nl+1)) then
1154	pa = log(p(nl))
1155	pb = log(pout(np))
1156	! special case: in case someone inadvertently enter p=0.
1157	if (p(nl+1).gt.0.d0) then
1158	pc = log(p(nl+1))
1159	else
1160	pc = log(1.d-4)
1161	end if
1162
1163	slope = (x(nl)-x(nl+1))/ (pa-pc)
1164	xout(np) = x(nl+1) + slope* (pb-pc)
1165	end if
1166	end do
1167	end do
1168	end if
1169
1170
1171	! place results in the return array;
1172	xxout = xout
1173
1174	END SUBROUTINE int1D
1175
1176	FUNCTION virtual (tmp,rmix)
1177	! This function returns virtual temperature in K, given temperature
1178	! in K and mixing ratio in kg/kg.
1179
1180	real :: tmp, rmix, virtual
1181
1182	virtual=tmp(0.622+rmix)/(0.622(1.+rmix))
1183
1184	END FUNCTION virtual
1185
1186	END MODULE module_ForInterpolate

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

Download in other formats:

Original Format