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/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

source: branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBS/obs_inter_sup.F90 @ 2797

Last change on this file since 2797 was 2715, checked in by rblod, 13 years ago

First attempt to put dynamic allocation on the trunk
  • Property svn:keywords set to Id
File size: 13.6 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, 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) :: kpk       ! Number of levels
59      INTEGER, DIMENSION(kptsi,kptsj,kobs), INTENT(IN) :: &
60         & kgrdi, &         ! i,j indicies for each stencil
61         & kgrdj
62      INTEGER, OPTIONAL, DIMENSION(kptsi,kptsj,kobs), INTENT(IN) :: &
63         & kproc            ! Precomputed processor for each i,j,iobs points
64      REAL(KIND=wp), DIMENSION(jpi,jpj,kpk), INTENT(IN) ::&
65         & pval             ! Local 3D array to extract data from
66      REAL(KIND=wp), DIMENSION(kptsi,kptsj,kpk,kobs), INTENT(OUT) ::&
67         & pgval            ! Stencil at each point
68      !! * Local declarations
69     
70      IF (ln_grid_global) THEN
71         
72         IF (PRESENT(kproc)) THEN
73
74            CALL obs_int_comm_3d_global( kptsi, kptsj, kobs, kpk, kgrdi, &
75               &                         kgrdj, pval, pgval, kproc=kproc )
76
77         ELSE
78
79            CALL obs_int_comm_3d_global( kptsi, kptsj, kobs, kpk, kgrdi, &
80               &                         kgrdj, pval, pgval )
81
82         ENDIF
83
84      ELSE
85
86         CALL obs_int_comm_3d_local( kptsi, kptsj, kobs, kpk, kgrdi, kgrdj, &
87            &                        pval, pgval )
88
89      ENDIF
90
91   END SUBROUTINE obs_int_comm_3d
92
93   SUBROUTINE obs_int_comm_2d( kptsi, kptsj, kobs, kgrdi, kgrdj, pval, pgval, &
94      &                        kproc )
95      !!----------------------------------------------------------------------
96      !!                    ***  ROUTINE obs_int_comm_2d  ***
97      !!         
98      !! ** Purpose : Get 2D interpolation stencil
99      !!
100      !! ** Method  : Call to obs_int_comm_3d
101      !!
102      !! ** Action  :
103      !!
104      !! History :
105      !!        !  08-02  (K. Mogensen)  Original code
106      !!----------------------------------------------------------------------
107      USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
108      USE wrk_nemo, ONLY: wrk_3d_1
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, DIMENSION(kptsi,kptsj,kobs), INTENT(IN) :: &
115         & kgrdi, &         ! i,j indicies for each stencil
116         & kgrdj
117      INTEGER, OPTIONAL, DIMENSION(kptsi,kptsj,kobs), INTENT(IN) :: &
118         & kproc            ! Precomputed processor for each i,j,iobs points
119      REAL(KIND=wp), DIMENSION(jpi,jpj), INTENT(IN) ::&
120         & pval             ! Local 3D array to extra data from
121      REAL(KIND=wp), DIMENSION(kptsi,kptsj,kobs), INTENT(OUT) ::&
122         & pgval            ! Stencil at each point
123      !! * Local declarations
124      REAL(KIND=wp), POINTER, DIMENSION(:,:,:) :: &
125         & zval
126      REAL(KIND=wp), DIMENSION(kptsi,kptsj,1,kobs) ::&
127         & zgval 
128
129      ! Check workspace array and set-up pointer
130      IF(wrk_in_use(3, 1))THEN
131         CALL ctl_stop('obs_int_comm_2d : requested workspace array unavailable.')
132         RETURN
133      END IF
134      zval => wrk_3d_1(:,:,1:1)
135
136      ! Set up local "3D" buffer
137
138      zval(:,:,1) = pval(:,:)
139
140      ! Call the 3D version
141
142      IF (PRESENT(kproc)) THEN
143
144         CALL obs_int_comm_3d( kptsi, kptsj, kobs, 1, kgrdi, kgrdj, zval, &
145            &                  zgval, kproc=kproc )
146      ELSE
147
148         CALL obs_int_comm_3d( kptsi, kptsj, kobs, 1, kgrdi, kgrdj, zval, &
149            &                  zgval )
150
151      ENDIF
152
153      ! Copy "3D" data back to 2D
154
155      pgval(:,:,:) = zgval(:,:,1,:)
156
157      ! 'Release' workspace array back to pool
158      IF(wrk_not_released(3, 1))THEN
159         CALL ctl_stop('obs_int_comm_2d : failed to release workspace array.')
160      END IF
161
162   END SUBROUTINE obs_int_comm_2d
163
164   SUBROUTINE obs_int_comm_3d_global( kptsi, kptsj, kobs, kpk, kgrdi, kgrdj, &
165      &                               pval, pgval, kproc )
166      !!----------------------------------------------------------------------
167      !!                    ***  ROUTINE obs_int_comm_3d_global  ***
168      !!         
169      !! ** Purpose : Get 3D interpolation stencil (global version)
170      !!
171      !! ** Method  : On-demand communication where each processor send its
172      !!              list of (i,j) of points to all processors and receive
173      !!              the corresponding values
174      !!
175      !! ** Action  :
176      !!
177      !! History :
178      !!        !  08-02  (K. Mogensen)  Original code
179      !!----------------------------------------------------------------------
180      !! * Arguments
181      INTEGER, INTENT(IN) :: kptsi     ! Number of i horizontal points per stencil
182      INTEGER, INTENT(IN) :: kptsj     ! Number of j horizontal points per stencil
183      INTEGER, INTENT(IN) :: kobs      ! Local number of observations
184      INTEGER, INTENT(IN) :: kpk       ! Number of levels
185      INTEGER, DIMENSION(kptsi,kptsj,kobs), INTENT(IN) :: &
186         & kgrdi, &         ! i,j indicies for each stencil
187         & kgrdj
188      INTEGER, OPTIONAL, DIMENSION(kptsi,kptsj,kobs), INTENT(IN) :: &
189         & kproc            ! Precomputed processor for each i,j,iobs points
190      REAL(KIND=wp), DIMENSION(jpi,jpj,kpk), INTENT(IN) ::&
191         & pval             ! Local 3D array to extract data from
192      REAL(KIND=wp), DIMENSION(kptsi,kptsj,kpk,kobs), INTENT(OUT) ::&
193         & pgval            ! Stencil at each point
194      !! * Local declarations
195      REAL(KIND=wp), DIMENSION(:,:), ALLOCATABLE :: &
196         & zsend, &
197         & zrecv
198      INTEGER, DIMENSION(:), ALLOCATABLE :: &
199         & igrdij_send, &
200         & igrdij_recv
201      INTEGER, DIMENSION(kptsi,kptsj,kobs) :: &
202         & iorder, &
203         & iproc
204      INTEGER :: nplocal(jpnij)
205      INTEGER :: npglobal(jpnij)
206      INTEGER :: ji
207      INTEGER :: jj
208      INTEGER :: jk
209      INTEGER :: jp
210      INTEGER :: jobs
211      INTEGER :: it
212      INTEGER :: itot
213      INTEGER :: ii
214      INTEGER :: ij
215
216      ! Check valid points
217     
218      IF ( ( MAXVAL(kgrdi) > jpiglo ) .OR. ( MINVAL(kgrdi) < 1 ) .OR. &
219         & ( MAXVAL(kgrdj) > jpjglo ) .OR. ( MINVAL(kgrdj) < 1 ) ) THEN
220         
221         CALL ctl_stop( 'Error in obs_int_comm_3d_global', &
222            &           'Point outside global domain' )
223         
224      ENDIF
225
226      ! Count number of points on each processors
227
228      nplocal(:) = 0
229      IF (PRESENT(kproc)) THEN
230         iproc(:,:,:) = kproc(:,:,:)
231         DO jobs = 1, kobs
232            DO jj = 1, kptsj
233               DO ji = 1, kptsi
234                  nplocal(iproc(ji,jj,jobs)) = nplocal(iproc(ji,jj,jobs)) + 1
235               END DO
236            END DO
237         END DO
238      ELSE
239         DO jobs = 1, kobs
240            DO jj = 1, kptsj
241               DO ji = 1, kptsi
242                  iproc(ji,jj,jobs) = mppmap(kgrdi(ji,jj,jobs),&
243                     &                       kgrdj(ji,jj,jobs))
244                  nplocal(iproc(ji,jj,jobs)) = nplocal(iproc(ji,jj,jobs)) + 1
245               END DO
246            END DO
247         END DO
248      ENDIF
249
250      ! Send local number of points and receive points on current domain
251
252      CALL mpp_alltoall_int( 1, nplocal, npglobal )
253
254      ! Allocate message parsing workspace
255
256      itot = SUM(npglobal)
257
258      ALLOCATE( &
259         & igrdij_send(kptsi*kptsj*kobs*2), &
260         & igrdij_recv(itot*2),             &
261         & zsend(kpk,itot),                 &
262         & zrecv(kpk,kptsi*kptsj*kobs)      &
263         & )
264
265      ! Pack buffers for list of points
266
267      it = 0
268      DO jp = 1, jpnij 
269         DO jobs = 1, kobs
270            DO jj = 1, kptsj
271               DO ji = 1, kptsi
272                  IF ( iproc(ji,jj,jobs) == jp ) THEN
273                     it = it + 1
274                     iorder(ji,jj,jobs) = it
275                     igrdij_send(2*it-1) = kgrdi(ji,jj,jobs)
276                     igrdij_send(2*it  ) = kgrdj(ji,jj,jobs)
277                  ENDIF
278               END DO
279            END DO
280         END DO
281      END DO
282
283      ! Send and recieve buffers for list of points
284
285      CALL mpp_alltoallv_int( igrdij_send, kptsi*kptsj*kobs*2, nplocal(:)*2, &
286         &                    igrdij_recv, itot*2, npglobal(:)*2 )
287
288      ! Pack interpolation data to be sent
289
290      DO ji = 1, itot
291         ii = mi1(igrdij_recv(2*ji-1))
292         ij = mj1(igrdij_recv(2*ji))
293         DO jk = 1, kpk
294            zsend(jk,ji) = pval(ii,ij,jk)
295         END DO
296      END DO
297
298      ! Re-adjust sizes
299
300      nplocal(:)  = kpk*nplocal(:)
301      npglobal(:) = kpk*npglobal(:)
302
303
304      ! Send and receive data for interpolation stencil
305
306      CALL mpp_alltoallv_real( zsend, kpk*itot,             npglobal, &
307         &                     zrecv, kpk*kptsi*kptsj*kobs, nplocal )
308
309      ! Copy the received data into output data structure
310     
311      DO jobs = 1, kobs
312         DO jj = 1, kptsj
313            DO ji = 1, kptsi
314               it = iorder(ji,jj,jobs) 
315               DO jk = 1, kpk
316                  pgval(ji,jj,jk,jobs) = zrecv(jk,it)
317               END DO
318            END DO
319         END DO
320      END DO
321
322      ! Deallocate message parsing workspace
323
324      DEALLOCATE( &
325         & igrdij_send, &
326         & igrdij_recv, &
327         & zsend,       &
328         & zrecv        &
329         & )
330
331   END SUBROUTINE obs_int_comm_3d_global
332   
333   SUBROUTINE obs_int_comm_3d_local( kptsi, kptsj, kobs, kpk, kgrdi, kgrdj, &
334      &                              pval, pgval )
335      !!----------------------------------------------------------------------
336      !!                    ***  ROUTINE obs_int_comm_3d_global  ***
337      !!         
338      !! ** Purpose : Get 3D interpolation stencil (global version)
339      !!
340      !! ** Method  : On-demand communication where each processor send its
341      !!              list of (i,j) of points to all processors and receive
342      !!              the corresponding values
343      !!
344      !! ** Action  :
345      !!
346      !! History :
347      !!        !  08-02  (K. Mogensen)  Original code
348      !!----------------------------------------------------------------------
349      !! * Arguments
350      INTEGER, INTENT(IN) :: kptsi        ! Number of i horizontal points per stencil
351      INTEGER, INTENT(IN) :: kptsj        ! Number of j horizontal points per stencil
352      INTEGER, INTENT(IN) :: kobs         ! Local number of observations
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(jpi,jpj,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.