Context Navigation

module_ForInterpolate.f90 @ 2261

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

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

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

Download in other formats:

Original Format