Newer
Older
1
2
3
4
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
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
MODULE obs_inter_sup
!!=====================================================================
!! *** MODULE obs_inter_sup ***
!! Observation diagnostics: Support for interpolation
!!=====================================================================
!!----------------------------------------------------------------------
!! obs_int_comm_3d : Get 3D interpolation stencil
!! obs_int_comm_2d : Get 2D interpolation stencil
!!---------------------------------------------------------------------
!! * Modules used
USE par_kind ! Precision variables
USE dom_oce ! Domain variables
USE mpp_map ! Map of processor points
USE lib_mpp ! MPP stuff
USE obs_mpp ! MPP stuff for observations
USE obs_grid ! Grid tools
USE in_out_manager ! I/O stuff
IMPLICIT NONE
!! * Routine accessibility
PRIVATE
PUBLIC obs_int_comm_3d, & ! Get 3D interpolation stencil
& obs_int_comm_2d ! Get 2D interpolation stencil
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: obs_inter_sup.F90 10068 2018-08-28 14:09:04Z nicolasmartin $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE obs_int_comm_3d( kptsi, kptsj, kobs, kpi, kpj, kpk, kgrdi, kgrdj, &
& pval, pgval, kproc )
!!----------------------------------------------------------------------
!! *** ROUTINE obs_int_comm_3d ***
!!
!! ** Purpose : Get 3D interpolation stencil
!!
!! ** Method : Either on-demand communication with
!! obs_int_comm_3d_global
!! or local memory with
!! obs_int_comm_3D_local
!! depending on ln_global_grid
!!
!! ** Action :
!!
!! History :
!! ! 08-02 (K. Mogensen) Original code
!!----------------------------------------------------------------------
!! * Arguments
INTEGER, INTENT(IN) :: kptsi ! Number of i horizontal points per stencil
INTEGER, INTENT(IN) :: kptsj ! Number of j horizontal points per stencil
INTEGER, INTENT(IN) :: kobs ! Local number of observations
INTEGER, INTENT(IN) :: kpi ! Number of points in i direction
INTEGER, INTENT(IN) :: kpj ! Number of points in j direction
INTEGER, INTENT(IN) :: kpk ! Number of levels
INTEGER, DIMENSION(kptsi,kptsj,kobs), INTENT(IN) :: &
& kgrdi, & ! i,j indicies for each stencil
& kgrdj
INTEGER, OPTIONAL, DIMENSION(kptsi,kptsj,kobs), INTENT(IN) :: &
& kproc ! Precomputed processor for each i,j,iobs points
REAL(KIND=wp), DIMENSION(kpi,kpj,kpk), INTENT(IN) ::&
& pval ! Local 3D array to extract data from
REAL(KIND=wp), DIMENSION(kptsi,kptsj,kpk,kobs), INTENT(OUT) ::&
& pgval ! Stencil at each point
!! * Local declarations
IF (ln_grid_global) THEN
IF (PRESENT(kproc)) THEN
CALL obs_int_comm_3d_global( kptsi, kptsj, kobs, kpi, kpj, kpk, kgrdi, &
& kgrdj, pval, pgval, kproc=kproc )
ELSE
CALL obs_int_comm_3d_global( kptsi, kptsj, kobs, kpi, kpj, kpk, kgrdi, &
& kgrdj, pval, pgval )
ENDIF
ELSE
CALL obs_int_comm_3d_local( kptsi, kptsj, kobs, kpi, kpj, kpk, kgrdi, kgrdj, &
& pval, pgval )
ENDIF
END SUBROUTINE obs_int_comm_3d
SUBROUTINE obs_int_comm_2d( kptsi, kptsj, kobs, kpi, kpj, kgrdi, kgrdj, pval, pgval, &
& kproc )
!!----------------------------------------------------------------------
!! *** ROUTINE obs_int_comm_2d ***
!!
!! ** Purpose : Get 2D interpolation stencil
!!
!! ** Method : Call to obs_int_comm_3d
!!
!! ** Action :
!!
!! History :
!! ! 08-02 (K. Mogensen) Original code
!!----------------------------------------------------------------------
!!
!! * Arguments
INTEGER, INTENT(IN) :: kptsi ! Number of i horizontal points per stencil
INTEGER, INTENT(IN) :: kptsj ! Number of j horizontal points per stencil
INTEGER, INTENT(IN) :: kobs ! Local number of observations
INTEGER, INTENT(IN) :: kpi ! Number of model grid points in i direction
INTEGER, INTENT(IN) :: kpj ! Number of model grid points in j direction
INTEGER, DIMENSION(kptsi,kptsj,kobs), INTENT(IN) :: &
& kgrdi, & ! i,j indicies for each stencil
& kgrdj
INTEGER, OPTIONAL, DIMENSION(kptsi,kptsj,kobs), INTENT(IN) :: &
& kproc ! Precomputed processor for each i,j,iobs points
REAL(KIND=wp), DIMENSION(kpi,kpj), INTENT(IN) ::&
& pval ! Local 3D array to extra data from
REAL(KIND=wp), DIMENSION(kptsi,kptsj,kobs), INTENT(OUT) ::&
& pgval ! Stencil at each point
!! * Local declarations
REAL(KIND=wp), DIMENSION(jpi,jpj,1) :: zval
REAL(KIND=wp), DIMENSION(kptsi,kptsj,1,kobs) ::&
& zgval
! Set up local "3D" buffer
zval(:,:,1) = pval(:,:)
! Call the 3D version
IF (PRESENT(kproc)) THEN
CALL obs_int_comm_3d( kptsi, kptsj, kobs, kpi, kpj, 1, kgrdi, kgrdj, zval, &
& zgval, kproc=kproc )
ELSE
CALL obs_int_comm_3d( kptsi, kptsj, kobs, kpi, kpj, 1, kgrdi, kgrdj, zval, &
& zgval )
ENDIF
! Copy "3D" data back to 2D
pgval(:,:,:) = zgval(:,:,1,:)
END SUBROUTINE obs_int_comm_2d
SUBROUTINE obs_int_comm_3d_global( kptsi, kptsj, kobs, kpi, kpj, kpk, kgrdi, kgrdj, &
& pval, pgval, kproc )
!!----------------------------------------------------------------------
!! *** ROUTINE obs_int_comm_3d_global ***
!!
!! ** Purpose : Get 3D interpolation stencil (global version)
!!
!! ** Method : On-demand communication where each processor send its
!! list of (i,j) of points to all processors and receive
!! the corresponding values
!!
!! ** Action :
!!
!! History :
!! ! 08-02 (K. Mogensen) Original code
!!----------------------------------------------------------------------
!! * Arguments
INTEGER, INTENT(IN) :: kptsi ! Number of i horizontal points per stencil
INTEGER, INTENT(IN) :: kptsj ! Number of j horizontal points per stencil
INTEGER, INTENT(IN) :: kobs ! Local number of observations
INTEGER, INTENT(IN) :: kpi ! Number of model points in i direction
INTEGER, INTENT(IN) :: kpj ! Number of model points in j direction
INTEGER, INTENT(IN) :: kpk ! Number of levels
INTEGER, DIMENSION(kptsi,kptsj,kobs), INTENT(IN) :: &
& kgrdi, & ! i,j indicies for each stencil
& kgrdj
INTEGER, OPTIONAL, DIMENSION(kptsi,kptsj,kobs), INTENT(IN) :: &
& kproc ! Precomputed processor for each i,j,iobs points
REAL(KIND=wp), DIMENSION(kpi,kpj,kpk), INTENT(IN) ::&
& pval ! Local 3D array to extract data from
REAL(KIND=wp), DIMENSION(kptsi,kptsj,kpk,kobs), INTENT(OUT) ::&
& pgval ! Stencil at each point
!! * Local declarations
REAL(KIND=wp), DIMENSION(:,:), ALLOCATABLE :: &
& zsend, &
& zrecv
INTEGER, DIMENSION(:), ALLOCATABLE :: &
& igrdij_send, &
& igrdij_recv
INTEGER, DIMENSION(kptsi,kptsj,kobs) :: &
& iorder, &
& iproc
INTEGER :: nplocal(jpnij)
INTEGER :: npglobal(jpnij)
INTEGER :: ji
INTEGER :: jj
INTEGER :: jk
INTEGER :: jp
INTEGER :: jobs
INTEGER :: it
INTEGER :: itot
INTEGER :: ii
INTEGER :: ij
! Check valid points
IF ( ( MAXVAL(kgrdi) > jpiglo ) .OR. ( MINVAL(kgrdi) < 1 ) .OR. &
& ( MAXVAL(kgrdj) > jpjglo ) .OR. ( MINVAL(kgrdj) < 1 ) ) THEN
CALL ctl_stop( 'Error in obs_int_comm_3d_global', &
& 'Point outside global domain' )
ENDIF
! Count number of points on each processors
nplocal(:) = 0
IF (PRESENT(kproc)) THEN
iproc(:,:,:) = kproc(:,:,:)
DO jobs = 1, kobs
DO jj = 1, kptsj
DO ji = 1, kptsi
nplocal(iproc(ji,jj,jobs)) = nplocal(iproc(ji,jj,jobs)) + 1
END DO
END DO
END DO
ELSE
DO jobs = 1, kobs
DO jj = 1, kptsj
DO ji = 1, kptsi
iproc(ji,jj,jobs) = mppmap(kgrdi(ji,jj,jobs),&
& kgrdj(ji,jj,jobs))
nplocal(iproc(ji,jj,jobs)) = nplocal(iproc(ji,jj,jobs)) + 1
END DO
END DO
END DO
ENDIF
! Send local number of points and receive points on current domain
CALL mpp_alltoall_int( 1, nplocal, npglobal )
! Allocate message parsing workspace
itot = SUM(npglobal)
ALLOCATE( &
& igrdij_send(kptsi*kptsj*kobs*2), &
& igrdij_recv(itot*2), &
& zsend(kpk,itot), &
& zrecv(kpk,kptsi*kptsj*kobs) &
& )
! Pack buffers for list of points
it = 0
DO jp = 1, jpnij
DO jobs = 1, kobs
DO jj = 1, kptsj
DO ji = 1, kptsi
IF ( iproc(ji,jj,jobs) == jp ) THEN
it = it + 1
iorder(ji,jj,jobs) = it
igrdij_send(2*it-1) = kgrdi(ji,jj,jobs)
igrdij_send(2*it ) = kgrdj(ji,jj,jobs)
ENDIF
END DO
END DO
END DO
END DO
! Send and recieve buffers for list of points
CALL mpp_alltoallv_int( igrdij_send, kptsi*kptsj*kobs*2, nplocal(:)*2, &
& igrdij_recv, itot*2, npglobal(:)*2 )
! Pack interpolation data to be sent
DO ji = 1, itot
ii = mi1(igrdij_recv(2*ji-1),nn_hls)
ij = mj1(igrdij_recv(2*ji ),nn_hls)
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
DO jk = 1, kpk
zsend(jk,ji) = pval(ii,ij,jk)
END DO
END DO
! Re-adjust sizes
nplocal(:) = kpk*nplocal(:)
npglobal(:) = kpk*npglobal(:)
! Send and receive data for interpolation stencil
CALL mpp_alltoallv_real( zsend, kpk*itot, npglobal, &
& zrecv, kpk*kptsi*kptsj*kobs, nplocal )
! Copy the received data into output data structure
DO jobs = 1, kobs
DO jj = 1, kptsj
DO ji = 1, kptsi
it = iorder(ji,jj,jobs)
DO jk = 1, kpk
pgval(ji,jj,jk,jobs) = zrecv(jk,it)
END DO
END DO
END DO
END DO
! Deallocate message parsing workspace
DEALLOCATE( &
& igrdij_send, &
& igrdij_recv, &
& zsend, &
& zrecv &
& )
END SUBROUTINE obs_int_comm_3d_global
SUBROUTINE obs_int_comm_3d_local( kptsi, kptsj, kobs, kpi, kpj, kpk, kgrdi, kgrdj, &
& pval, pgval )
!!----------------------------------------------------------------------
!! *** ROUTINE obs_int_comm_3d_global ***
!!
!! ** Purpose : Get 3D interpolation stencil (global version)
!!
!! ** Method : On-demand communication where each processor send its
!! list of (i,j) of points to all processors and receive
!! the corresponding values
!!
!! ** Action :
!!
!! History :
!! ! 08-02 (K. Mogensen) Original code
!!----------------------------------------------------------------------
!! * Arguments
INTEGER, INTENT(IN) :: kptsi ! Number of i horizontal points per stencil
INTEGER, INTENT(IN) :: kptsj ! Number of j horizontal points per stencil
INTEGER, INTENT(IN) :: kobs ! Local number of observations
INTEGER, INTENT(IN) :: kpi ! Number of model points in i direction
INTEGER, INTENT(IN) :: kpj ! Number of model points in j direction
INTEGER, INTENT(IN) :: kpk ! Number of levels
INTEGER, DIMENSION(kptsi,kptsj,kobs), INTENT(IN) :: &
& kgrdi, & ! i,j indicies for each stencil
& kgrdj
REAL(KIND=wp), DIMENSION(kpi,kpj,kpk), INTENT(IN) ::&
& pval ! Local 3D array to extract data from
REAL(KIND=wp), DIMENSION(kptsi,kptsj,kpk,kobs), INTENT(OUT) ::&
& pgval ! Stencil at each point
!! * Local declarations
INTEGER :: ji
INTEGER :: jj
INTEGER :: jk
INTEGER :: jobs
! Check valid points
IF ( ( MAXVAL(kgrdi) > jpi ) .OR. ( MINVAL(kgrdi) < 1 ) .OR. &
& ( MAXVAL(kgrdj) > jpj ) .OR. ( MINVAL(kgrdj) < 1 ) ) THEN
CALL ctl_stop( 'Error in obs_int_comm_3d_local', &
& 'Point outside local domain' )
ENDIF
! Copy local data
DO jobs = 1, kobs
DO jj = 1, kptsj
DO ji = 1, kptsi
DO jk = 1, kpk
pgval(ji,jj,jk,jobs) = &
& pval(kgrdi(ji,jj,jobs),kgrdj(ji,jj,jobs),jk)
END DO
END DO
END DO
END DO
END SUBROUTINE obs_int_comm_3d_local
END MODULE obs_inter_sup