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

source: branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/OBS/obs_inter_sup.F90 @ 3183

Last change on this file since 3183 was 3183, checked in by davestorkey, 13 years ago

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