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/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/OPA_SRC/OBS/obs_inter_sup.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 7 years ago

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