New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
obs_inter_sup.F90 in branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

source: branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/OBS/obs_inter_sup.F90 @ 10247

Last change on this file since 10247 was 10247, checked in by kingr, 6 years ago

Merged dev_r5518_obs_oper_update@10168

File size: 14.0 KB
Line 
1MODULE obs_inter_sup
2   !!=====================================================================
3   !!                       ***  MODULE obs_inter_sup  ***
4   !! Observation diagnostics: Support for interpolation
5   !!=====================================================================
6
7   !!----------------------------------------------------------------------
8   !!   obs_int_comm_3d : Get 3D interpolation stencil
9   !!   obs_int_comm_2d : Get 2D interpolation stencil
10   !!---------------------------------------------------------------------
11   !! * Modules used
12   USE wrk_nemo        ! Memory Allocation
13   USE par_kind        ! Precision variables
14   USE dom_oce         ! Domain variables
15   USE mpp_map         ! Map of processor points
16   USE lib_mpp         ! MPP stuff
17   USE obs_mpp         ! MPP stuff for observations
18   USE obs_grid        ! Grid tools
19   USE in_out_manager  ! I/O stuff
20
21   IMPLICIT NONE
22
23   !! * Routine accessibility
24   PRIVATE
25   
26   PUBLIC obs_int_comm_3d, & ! Get 3D interpolation stencil
27      &   obs_int_comm_2d    ! Get 2D interpolation stencil
28   
29   !!----------------------------------------------------------------------
30   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
31   !! $Id$
32   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
33   !!----------------------------------------------------------------------
34
35CONTAINS
36
37   SUBROUTINE obs_int_comm_3d( kptsi, kptsj, kobs, kpi, kpj, kpk, kgrdi, kgrdj, &
38      &                        pval, pgval, kproc )
39      !!----------------------------------------------------------------------
40      !!                    ***  ROUTINE obs_int_comm_3d  ***
41      !!         
42      !! ** Purpose : Get 3D interpolation stencil
43      !!
44      !! ** Method  : Either on-demand communication with
45      !!              obs_int_comm_3d_global
46      !!              or local memory with
47      !!              obs_int_comm_3D_local
48      !!              depending on ln_global_grid
49      !!
50      !! ** Action  :
51      !!
52      !! History :
53      !!        !  08-02  (K. Mogensen)  Original code
54      !!----------------------------------------------------------------------
55      !! * Arguments
56      INTEGER, INTENT(IN) :: kptsi     ! Number of i horizontal points per stencil
57      INTEGER, INTENT(IN) :: kptsj     ! Number of j horizontal points per stencil
58      INTEGER, INTENT(IN) :: kobs      ! Local number of observations
59      INTEGER, INTENT(IN) :: kpi       ! Number of points in i direction
60      INTEGER, INTENT(IN) :: kpj       ! Number of points in j direction
61      INTEGER, INTENT(IN) :: kpk       ! Number of levels
62      INTEGER, DIMENSION(kptsi,kptsj,kobs), INTENT(IN) :: &
63         & kgrdi, &         ! i,j indicies for each stencil
64         & kgrdj
65      INTEGER, OPTIONAL, DIMENSION(kptsi,kptsj,kobs), INTENT(IN) :: &
66         & kproc            ! Precomputed processor for each i,j,iobs points
67      REAL(KIND=wp), DIMENSION(kpi,kpj,kpk), INTENT(IN) ::&
68         & pval             ! Local 3D array to extract data from
69      REAL(KIND=wp), DIMENSION(kptsi,kptsj,kpk,kobs), INTENT(OUT) ::&
70         & pgval            ! Stencil at each point
71      !! * Local declarations
72     
73      IF (ln_grid_global) THEN
74         
75         IF (PRESENT(kproc)) THEN
76
77            CALL obs_int_comm_3d_global( kptsi, kptsj, kobs, kpi, kpj, kpk, kgrdi, &
78               &                         kgrdj, pval, pgval, kproc=kproc )
79
80         ELSE
81
82            CALL obs_int_comm_3d_global( kptsi, kptsj, kobs, kpi, kpj, kpk, kgrdi, &
83               &                         kgrdj, pval, pgval )
84
85         ENDIF
86
87      ELSE
88
89         CALL obs_int_comm_3d_local( kptsi, kptsj, kobs, kpi, kpj, kpk, kgrdi, kgrdj, &
90            &                        pval, pgval )
91
92      ENDIF
93
94   END SUBROUTINE obs_int_comm_3d
95
96   SUBROUTINE obs_int_comm_2d( kptsi, kptsj, kobs, kpi, kpj, kgrdi, kgrdj, pval, pgval, &
97      &                        kproc )
98      !!----------------------------------------------------------------------
99      !!                    ***  ROUTINE obs_int_comm_2d  ***
100      !!         
101      !! ** Purpose : Get 2D interpolation stencil
102      !!
103      !! ** Method  : Call to obs_int_comm_3d
104      !!
105      !! ** Action  :
106      !!
107      !! History :
108      !!        !  08-02  (K. Mogensen)  Original code
109      !!----------------------------------------------------------------------
110      !!
111      !! * Arguments
112      INTEGER, INTENT(IN) :: kptsi        ! Number of i horizontal points per stencil
113      INTEGER, INTENT(IN) :: kptsj        ! Number of j horizontal points per stencil
114      INTEGER, INTENT(IN) :: kobs          ! Local number of observations
115      INTEGER, INTENT(IN) :: kpi          ! Number of model grid points in i direction
116      INTEGER, INTENT(IN) :: kpj          ! Number of model grid points in j direction
117      INTEGER, DIMENSION(kptsi,kptsj,kobs), INTENT(IN) :: &
118         & kgrdi, &         ! i,j indicies for each stencil
119         & kgrdj
120      INTEGER, OPTIONAL, DIMENSION(kptsi,kptsj,kobs), INTENT(IN) :: &
121         & kproc            ! Precomputed processor for each i,j,iobs points
122      REAL(KIND=wp), DIMENSION(kpi,kpj), INTENT(IN) ::&
123         & pval             ! Local 3D array to extra data from
124      REAL(KIND=wp), DIMENSION(kptsi,kptsj,kobs), INTENT(OUT) ::&
125         & pgval            ! Stencil at each point
126      !! * Local declarations
127      REAL(KIND=wp), POINTER, DIMENSION(:,:,:) ::   zval
128      REAL(KIND=wp), DIMENSION(kptsi,kptsj,1,kobs) ::&
129         & zgval 
130
131      ! Check workspace array and set-up pointer
132      CALL wrk_alloc(jpi,jpj,1,zval) 
133
134      ! Set up local "3D" buffer
135
136      zval(:,:,1) = pval(:,:)
137
138      ! Call the 3D version
139
140      IF (PRESENT(kproc)) THEN
141
142         CALL obs_int_comm_3d( kptsi, kptsj, kobs, kpi, kpj, 1, kgrdi, kgrdj, zval, &
143            &                  zgval, kproc=kproc )
144      ELSE
145
146         CALL obs_int_comm_3d( kptsi, kptsj, kobs, kpi, kpj, 1, kgrdi, kgrdj, zval, &
147            &                  zgval )
148
149      ENDIF
150
151      ! Copy "3D" data back to 2D
152
153      pgval(:,:,:) = zgval(:,:,1,:)
154
155      ! 'Release' workspace array back to pool
156      CALL wrk_dealloc(jpi,jpj,1,zval) 
157
158   END SUBROUTINE obs_int_comm_2d
159
160   SUBROUTINE obs_int_comm_3d_global( kptsi, kptsj, kobs, kpi, kpj, kpk, kgrdi, kgrdj, &
161      &                               pval, pgval, kproc )
162      !!----------------------------------------------------------------------
163      !!                    ***  ROUTINE obs_int_comm_3d_global  ***
164      !!         
165      !! ** Purpose : Get 3D interpolation stencil (global version)
166      !!
167      !! ** Method  : On-demand communication where each processor send its
168      !!              list of (i,j) of points to all processors and receive
169      !!              the corresponding values
170      !!
171      !! ** Action  :
172      !!
173      !! History :
174      !!        !  08-02  (K. Mogensen)  Original code
175      !!----------------------------------------------------------------------
176      !! * Arguments
177      INTEGER, INTENT(IN) :: kptsi     ! Number of i horizontal points per stencil
178      INTEGER, INTENT(IN) :: kptsj     ! Number of j horizontal points per stencil
179      INTEGER, INTENT(IN) :: kobs      ! Local number of observations
180      INTEGER, INTENT(IN) :: kpi       ! Number of model points in i direction
181      INTEGER, INTENT(IN) :: kpj       ! Number of model points in j direction
182      INTEGER, INTENT(IN) :: kpk       ! Number of levels
183      INTEGER, DIMENSION(kptsi,kptsj,kobs), INTENT(IN) :: &
184         & kgrdi, &         ! i,j indicies for each stencil
185         & kgrdj
186      INTEGER, OPTIONAL, DIMENSION(kptsi,kptsj,kobs), INTENT(IN) :: &
187         & kproc            ! Precomputed processor for each i,j,iobs points
188      REAL(KIND=wp), DIMENSION(kpi,kpj,kpk), INTENT(IN) ::&
189         & pval             ! Local 3D array to extract data from
190      REAL(KIND=wp), DIMENSION(kptsi,kptsj,kpk,kobs), INTENT(OUT) ::&
191         & pgval            ! Stencil at each point
192      !! * Local declarations
193      REAL(KIND=wp), DIMENSION(:,:), ALLOCATABLE :: &
194         & zsend, &
195         & zrecv
196      INTEGER, DIMENSION(:), ALLOCATABLE :: &
197         & igrdij_send, &
198         & igrdij_recv
199      INTEGER, DIMENSION(kptsi,kptsj,kobs) :: &
200         & iorder, &
201         & iproc
202      INTEGER :: nplocal(jpnij)
203      INTEGER :: npglobal(jpnij)
204      INTEGER :: ji
205      INTEGER :: jj
206      INTEGER :: jk
207      INTEGER :: jp
208      INTEGER :: jobs
209      INTEGER :: it
210      INTEGER :: itot
211      INTEGER :: ii
212      INTEGER :: ij
213
214      ! Check valid points
215
216      IF ( ( MAXVAL(kgrdi) > jpiglo ) .OR. ( MINVAL(kgrdi) < 1 ) .OR. &
217         & ( MAXVAL(kgrdj) > jpjglo ) .OR. ( MINVAL(kgrdj) < 1 ) ) THEN
218
219         CALL ctl_stop( 'Error in obs_int_comm_3d_global', &
220            &           'Point outside global domain' )
221
222      ENDIF
223
224      ! Count number of points on each processors
225
226      nplocal(:) = 0
227      IF (PRESENT(kproc)) THEN
228         iproc(:,:,:) = kproc(:,:,:)
229         DO jobs = 1, kobs
230            DO jj = 1, kptsj
231               DO ji = 1, kptsi
232                  nplocal(iproc(ji,jj,jobs)) = nplocal(iproc(ji,jj,jobs)) + 1
233               END DO
234            END DO
235         END DO
236      ELSE
237         DO jobs = 1, kobs
238            DO jj = 1, kptsj
239               DO ji = 1, kptsi
240                  iproc(ji,jj,jobs) = mppmap(kgrdi(ji,jj,jobs),&
241                     &                       kgrdj(ji,jj,jobs))
242                  nplocal(iproc(ji,jj,jobs)) = nplocal(iproc(ji,jj,jobs)) + 1
243               END DO
244            END DO
245         END DO
246      ENDIF
247
248      ! Send local number of points and receive points on current domain
249
250      CALL mpp_alltoall_int( 1, nplocal, npglobal )
251
252      ! Allocate message parsing workspace
253
254      itot = SUM(npglobal)
255
256      ALLOCATE( &
257         & igrdij_send(kptsi*kptsj*kobs*2), &
258         & igrdij_recv(itot*2),             &
259         & zsend(kpk,itot),                 &
260         & zrecv(kpk,kptsi*kptsj*kobs)      &
261         & )
262
263      ! Pack buffers for list of points
264
265      it = 0
266      DO jp = 1, jpnij 
267         DO jobs = 1, kobs
268            DO jj = 1, kptsj
269               DO ji = 1, kptsi
270                  IF ( iproc(ji,jj,jobs) == jp ) THEN
271                     it = it + 1
272                     iorder(ji,jj,jobs) = it
273                     igrdij_send(2*it-1) = kgrdi(ji,jj,jobs)
274                     igrdij_send(2*it  ) = kgrdj(ji,jj,jobs)
275                  ENDIF
276               END DO
277            END DO
278         END DO
279      END DO
280
281      ! Send and recieve buffers for list of points
282
283      CALL mpp_alltoallv_int( igrdij_send, kptsi*kptsj*kobs*2, nplocal(:)*2, &
284         &                    igrdij_recv, itot*2, npglobal(:)*2 )
285
286      ! Pack interpolation data to be sent
287
288      DO ji = 1, itot
289         ii = mi1(igrdij_recv(2*ji-1))
290         ij = mj1(igrdij_recv(2*ji))
291         DO jk = 1, kpk
292            zsend(jk,ji) = pval(ii,ij,jk)
293         END DO
294      END DO
295
296      ! Re-adjust sizes
297
298      nplocal(:)  = kpk*nplocal(:)
299      npglobal(:) = kpk*npglobal(:)
300
301
302      ! Send and receive data for interpolation stencil
303
304      CALL mpp_alltoallv_real( zsend, kpk*itot,             npglobal, &
305         &                     zrecv, kpk*kptsi*kptsj*kobs, nplocal )
306
307      ! Copy the received data into output data structure
308     
309      DO jobs = 1, kobs
310         DO jj = 1, kptsj
311            DO ji = 1, kptsi
312               it = iorder(ji,jj,jobs) 
313               DO jk = 1, kpk
314                  pgval(ji,jj,jk,jobs) = zrecv(jk,it)
315               END DO
316            END DO
317         END DO
318      END DO
319
320      ! Deallocate message parsing workspace
321
322      DEALLOCATE( &
323         & igrdij_send, &
324         & igrdij_recv, &
325         & zsend,       &
326         & zrecv        &
327         & )
328
329   END SUBROUTINE obs_int_comm_3d_global
330   
331   SUBROUTINE obs_int_comm_3d_local( kptsi, kptsj, kobs, kpi, kpj, kpk, kgrdi, kgrdj, &
332      &                              pval, pgval )
333      !!----------------------------------------------------------------------
334      !!                    ***  ROUTINE obs_int_comm_3d_global  ***
335      !!         
336      !! ** Purpose : Get 3D interpolation stencil (global version)
337      !!
338      !! ** Method  : On-demand communication where each processor send its
339      !!              list of (i,j) of points to all processors and receive
340      !!              the corresponding values
341      !!
342      !! ** Action  :
343      !!
344      !! History :
345      !!        !  08-02  (K. Mogensen)  Original code
346      !!----------------------------------------------------------------------
347      !! * Arguments
348      INTEGER, INTENT(IN) :: kptsi        ! Number of i horizontal points per stencil
349      INTEGER, INTENT(IN) :: kptsj        ! Number of j horizontal points per stencil
350      INTEGER, INTENT(IN) :: kobs         ! Local number of observations
351      INTEGER, INTENT(IN) :: kpi          ! Number of model points in i direction
352      INTEGER, INTENT(IN) :: kpj          ! Number of model points in j direction
353      INTEGER, INTENT(IN) :: kpk          ! Number of levels
354      INTEGER, DIMENSION(kptsi,kptsj,kobs), INTENT(IN) :: &
355         & kgrdi, &         ! i,j indicies for each stencil
356         & kgrdj
357      REAL(KIND=wp), DIMENSION(kpi,kpj,kpk), INTENT(IN) ::&
358         & pval             ! Local 3D array to extract data from
359      REAL(KIND=wp), DIMENSION(kptsi,kptsj,kpk,kobs), INTENT(OUT) ::&
360         & pgval            ! Stencil at each point
361      !! * Local declarations
362      INTEGER ::  ji
363      INTEGER ::  jj
364      INTEGER ::  jk
365      INTEGER ::  jobs
366
367      ! Check valid points
368
369      IF ( ( MAXVAL(kgrdi) > jpi ) .OR. ( MINVAL(kgrdi) < 1 ) .OR. &
370         & ( MAXVAL(kgrdj) > jpj ) .OR. ( MINVAL(kgrdj) < 1 ) ) THEN
371         
372         CALL ctl_stop( 'Error in obs_int_comm_3d_local', &
373            &           'Point outside local domain' )
374         
375      ENDIF
376
377      ! Copy local data
378
379      DO jobs = 1, kobs
380         DO jj = 1, kptsj
381            DO ji = 1, kptsi
382               DO jk = 1, kpk
383                  pgval(ji,jj,jk,jobs) = &
384                     &            pval(kgrdi(ji,jj,jobs),kgrdj(ji,jj,jobs),jk)
385               END DO
386            END DO
387         END DO
388      END DO
389
390   END SUBROUTINE obs_int_comm_3d_local
391
392END MODULE obs_inter_sup
393 
Note: See TracBrowser for help on using the repository browser.