Context Navigation

module_ForInterpolate.F90 @ 929

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

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

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

Download in other formats:

Original Format