Context Navigation

module_ForInterpolate.F90 @ 1476

Last change on this file since 1476 was 1355, checked in by lfita, 9 years ago
Fixing some issues of the `reprojection' modules Generalizing the use of `module_generic' Definition of `errormsg' and` warnmsg' in 'module_generic.F90'
File size: 64.6 KB

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

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

Download in other formats:

Original Format