Newer
Older
!SURFEX_LIC Copyright 1994-2014 Meteo-France
!SURFEX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
!SURFEX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
!SURFEX_LIC for details. version 1.
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
! #########
SUBROUTINE INTERPOL_NPTS(HPROGRAM,KLUOUT,KNPTS,KCODE,PX,PY,PFIELD)
! #########################################################
!
!!**** *INTERPOL_NPTS* interpolates with ###ine f77 programs a 2D field
!! from all grid points valid values
!!
!! PURPOSE
!! -------
!!
!! The points are all on only one grid (defined with the coordinates
!! of all the points). The code to apply for each point is:
!!
!! KCODE>0 : data point (with field valid for interpolation)
!! KCODE=-1: point to ignore
!! KCODE=0 : point to interpolate
!!
!!
!!
!! METHOD
!! ------
!!
!! EXTERNAL
!! --------
!!
!! IMPLICIT ARGUMENTS
!! ------------------
!!
!!
!!
!! REFERENCE
!! ---------
!!
!! AUTHOR
!! ------
!!
!! V. Masson Meteo-France
!!
!! MODIFICATION
!! ------------
!!
!! Original 03/2004
!! Modification
!----------------------------------------------------------------------------
!
!* 0. DECLARATION
! -----------
!
USE MODD_SURF_PAR, ONLY : XUNDEF
USE MODD_SURF_ATM_GRID_n, ONLY : CGRID, XGRID_PAR, NGRID_PAR, NNEAR
USE MODD_SURF_ATM_n, ONLY : NSIZE_FULL, NDIM_FULL
!
USE MODI_GET_INTERP_HALO
USE MODI_GET_NEAR_MESHES
!
USE YOMHOOK ,ONLY : LHOOK, DR_HOOK
USE PARKIND1 ,ONLY : JPRB
!
IMPLICIT NONE
!
!* 0.1 Declaration of arguments
! ------------------------
!
CHARACTER(LEN=6), INTENT(IN) :: HPROGRAM ! host program
INTEGER, INTENT(IN) :: KLUOUT ! output listing
INTEGER, INTENT(IN) :: KNPTS ! number of points to interpolate with
INTEGER,DIMENSION(:), INTENT(INOUT) :: KCODE ! code for each point
! >0 point used for interpolation
! 0 point to interpolate
! -1 point not used
! -2 point not used
! ! -3 if spline is no computed
! ! for this point
REAL, DIMENSION(:), INTENT(IN) :: PX ! x of each grid mesh.
REAL, DIMENSION(:), INTENT(IN) :: PY ! y of each grid mesh.
REAL, DIMENSION(:,:),INTENT(INOUT) :: PFIELD ! pgd field on grid mesh.
!
!* 0.2 Declaration of local variables
! ------------------------------
!
INTEGER :: IL ! number of points
INTEGER :: JD ! data point index
INTEGER :: JS ! loop counter on data points
INTEGER :: JL ! loop counter on points to initialize
INTEGER :: JP, JPP ! loops counter on KNPTS
REAL :: ZDIST ! square distance between two interpolating and interpolated points
REAL, DIMENSION(0:KNPTS) :: ZNDIST ! 3 nearest square distances
REAL, DIMENSION(0:KNPTS,SIZE(PFIELD,2)) :: ZNVAL ! 3 corresponding field values
REAL, DIMENSION(SIZE(PFIELD,2)) :: ZSUM
!
INTEGER :: INEAR_NBR ! number of points to scan
INTEGER :: JLIST ! loop counter on points to interpolate
INTEGER :: ICOUNT ! counter
INTEGER :: INPTS
INTEGER :: ISCAN ! number of points to scan
INTEGER, DIMENSION(:), ALLOCATABLE :: IINDEX ! list of index to scan
INTEGER :: IHALO ! halo available
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!-------------------------------------------------------------------------------
!
IF (LHOOK) CALL DR_HOOK('INTERPOL_NPTS',0,ZHOOK_HANDLE)
IL = SIZE(PFIELD,1)
!
CALL GET_INTERP_HALO(HPROGRAM,CGRID,IHALO)
!
INEAR_NBR = (2*IHALO+1)**2
!
!
ALLOCATE(IINDEX(IL))
IINDEX(:) = 0
!
!
IF (.NOT.ASSOCIATED(NNEAR)) THEN
ALLOCATE(NNEAR(IL,INEAR_NBR))
NNEAR(:,:) = 0
CALL GET_NEAR_MESHES(CGRID,NGRID_PAR,NSIZE_FULL,XGRID_PAR,INEAR_NBR,NNEAR)
ENDIF
!
DO JL=1,IL
IF (KCODE(JL)/=0) CYCLE
ZNDIST (1:KNPTS) = 1.E20
ZNDIST (0) = 0.
ZNVAL(0:KNPTS,:) = 0.
!
ICOUNT = 0
DO JD=1,INEAR_NBR
IF (NNEAR(JL,JD)>0) THEN
IF (KCODE(NNEAR(JL,JD))>0) THEN
ICOUNT = ICOUNT+1
IINDEX(ICOUNT) = NNEAR(JL,JD)
END IF
END IF
END DO
!
IF (ICOUNT>=1) THEN
ISCAN = ICOUNT
INPTS = MIN(ICOUNT,KNPTS)
ELSE
KCODE(JL) = -4
CYCLE
END IF
!
!
DO JS=1,ISCAN
!
JD = IINDEX(JS)
!
ZDIST= ( ( PX(JD)-PX(JL) ) ** 2 ) + ( ( PY(JD)-PY(JL) ) ** 2 )
!
IF ( ZDIST>ZNDIST(INPTS) ) CYCLE
!
DO JP = INPTS,1,-1
!
IF ( ZDIST>ZNDIST(JP-1) ) THEN
!
IF ( JP<INPTS ) THEN
DO JPP = INPTS,JP+1,-1
ZNDIST(JPP) = ZNDIST(JPP-1)
ZNVAL(JPP,:) = ZNVAL(JPP-1,:)
ENDDO
ENDIF
!
ZNDIST(JP) = ZDIST
ZNVAL(JP,:) = PFIELD(JD,:)
!
EXIT
!
ENDIF
!
ENDDO
!
ENDDO
!
ZNDIST(:) = SQRT(ZNDIST(:))
!
PFIELD(JL,:) = 0.
ZSUM(:) = 0.
DO JP = 1, INPTS
PFIELD(JL,:) = PFIELD(JL,:) + ZNVAL(JP,:)/ZNDIST(JP)
ZSUM(:) = ZSUM(:) + 1./ZNDIST(JP)
ENDDO
PFIELD(JL,:) = PFIELD(JL,:) / ZSUM(:)
!
END DO
!
DEALLOCATE(IINDEX )
IF (LHOOK) CALL DR_HOOK('INTERPOL_NPTS',1,ZHOOK_HANDLE)
!-------------------------------------------------------------------------------
!
END SUBROUTINE INTERPOL_NPTS