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

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/OBS/obs_inter_sup.F90 @ 4416

Last change on this file since 4416 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

Commit to bring the rest of the DCSE NEMO development branch
in line with the latest development version. This includes
array index re-ordering of all OPA_SRC/.
  • Property svn:keywords set to Id
File size: 15.1 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   !! * Control permutation of array indices
29#  include "dom_oce_ftrans.h90"
30
31   !!----------------------------------------------------------------------
32   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
33   !! $Id$
34   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
35   !!----------------------------------------------------------------------
36
37CONTAINS
38
39   SUBROUTINE obs_int_comm_3d( kptsi, kptsj, kobs, kpk, kgrdi, kgrdj, &
40      &                        pval, pgval, kproc )
41      !!----------------------------------------------------------------------
42      !!                    ***  ROUTINE obs_int_comm_3d  ***
43      !!         
44      !! ** Purpose : Get 3D interpolation stencil
45      !!
46      !! ** Method  : Either on-demand communication with
47      !!              obs_int_comm_3d_global
48      !!              or local memory with
49      !!              obs_int_comm_3D_local
50      !!              depending on ln_global_grid
51      !!
52      !! ** Action  :
53      !!
54      !! History :
55      !!        !  08-02  (K. Mogensen)  Original code
56      !!----------------------------------------------------------------------
57      !! * Arguments
58      INTEGER, INTENT(IN) :: kptsi     ! Number of i horizontal points per stencil
59      INTEGER, INTENT(IN) :: kptsj     ! Number of j horizontal points per stencil
60      INTEGER, INTENT(IN) :: kobs      ! Local number of observations
61      INTEGER, INTENT(IN) :: kpk       ! Number of levels
62      INTEGER, DIMENSION(kptsi,kptsj,kobs), INTENT(IN) :: &
63         & kgrdi, &         ! i,j indices 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
68!! DCSE_NEMO: This style defeats ftrans
69!      REAL(KIND=wp), DIMENSION(jpi,jpj,kpk), INTENT(IN) ::&
70!        & pval             ! Local 3D array to extract data from
71!FTRANS pval :I :I :z
72      REAL(KIND=wp), INTENT(IN) ::&
73         & pval(jpi,jpj,kpk) ! Local 3D array to extract data from
74      REAL(KIND=wp), DIMENSION(kptsi,kptsj,kpk,kobs), INTENT(OUT) ::&
75         & pgval            ! Stencil at each point
76      !! * Local declarations
77     
78#if defined key_z_first
79      IF ( kpk /= jpk ) THEN
80         CALL ctl_stop( 'Error in obs_int_comm_3d', &
81            &           'index reordering requires that jpk==kpk' )
82      ENDIF
83#endif
84
85      IF (ln_grid_global) THEN
86         
87         IF (PRESENT(kproc)) THEN
88
89            CALL obs_int_comm_3d_global( kptsi, kptsj, kobs, kpk, kgrdi, &
90               &                         kgrdj, pval, pgval, kproc=kproc )
91
92         ELSE
93
94            CALL obs_int_comm_3d_global( kptsi, kptsj, kobs, kpk, kgrdi, &
95               &                         kgrdj, pval, pgval )
96
97         ENDIF
98
99      ELSE
100
101         CALL obs_int_comm_3d_local( kptsi, kptsj, kobs, kpk, kgrdi, kgrdj, &
102            &                        pval, pgval )
103
104      ENDIF
105
106   END SUBROUTINE obs_int_comm_3d
107
108   SUBROUTINE obs_int_comm_2d( kptsi, kptsj, kobs, kgrdi, kgrdj, pval, pgval, &
109      &                        kproc )
110      !!----------------------------------------------------------------------
111      !!                    ***  ROUTINE obs_int_comm_2d  ***
112      !!         
113      !! ** Purpose : Get 2D interpolation stencil
114      !!
115      !! ** Method  : Call to obs_int_comm_3d
116      !!
117      !! ** Action  :
118      !!
119      !! History :
120      !!        !  08-02  (K. Mogensen)  Original code
121      !!----------------------------------------------------------------------
122      USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
123      USE wrk_nemo, ONLY: wrk_3d_1
124
125   !! * Control permutation of array indices
126!FTRANS CLEAR
127#  include "dom_oce_ftrans.h90"
128!FTRANS wrk_3d_1 :I :I :z
129!FTRANS zval :I :I :z
130
131      !!
132      !! * Arguments
133      INTEGER, INTENT(IN) :: kptsi        ! Number of i horizontal points per stencil
134      INTEGER, INTENT(IN) :: kptsj        ! Number of j horizontal points per stencil
135      INTEGER, INTENT(IN) :: kobs          ! Local number of observations
136      INTEGER, DIMENSION(kptsi,kptsj,kobs), INTENT(IN) :: &
137         & kgrdi, &         ! i,j indicies for each stencil
138         & kgrdj
139      INTEGER, OPTIONAL, DIMENSION(kptsi,kptsj,kobs), INTENT(IN) :: &
140         & kproc            ! Precomputed processor for each i,j,iobs points
141      REAL(KIND=wp), DIMENSION(jpi,jpj), INTENT(IN) ::&
142         & pval             ! Local 3D array to extra data from
143      REAL(KIND=wp), DIMENSION(kptsi,kptsj,kobs), INTENT(OUT) ::&
144         & pgval            ! Stencil at each point
145      !! * Local declarations
146      REAL(KIND=wp), POINTER, DIMENSION(:,:,:) :: &
147         & zval
148      REAL(KIND=wp), DIMENSION(kptsi,kptsj,1,kobs) ::&
149         & zgval 
150
151      ! Check workspace array and set-up pointer
152      IF(wrk_in_use(3, 1))THEN
153         CALL ctl_stop('obs_int_comm_2d : requested workspace array unavailable.')
154         RETURN
155      END IF
156
157      zval => wrk_3d_1
158
159      ! Set up local "3D" buffer
160
161      zval(:,:,1) = pval(:,:)
162
163      ! Call the 3D version
164
165!! DCSE_NEMO: this is not going to work with index re-ordering
166!! Really want obs_int_comm_2d to do its own stuff, instead of calling
167!! obs_int_comm_3d !!
168
169      IF (PRESENT(kproc)) THEN
170
171         CALL obs_int_comm_3d( kptsi, kptsj, kobs, 1, kgrdi, kgrdj, zval, &
172            &                  zgval, kproc=kproc )
173      ELSE
174
175         CALL obs_int_comm_3d( kptsi, kptsj, kobs, 1, kgrdi, kgrdj, zval, &
176            &                  zgval )
177
178      ENDIF
179
180      ! Copy "3D" data back to 2D
181
182      pgval(:,:,:) = zgval(:,:,1,:)
183
184      ! 'Release' workspace array back to pool
185      IF(wrk_not_released(3, 1))THEN
186         CALL ctl_stop('obs_int_comm_2d : failed to release workspace array.')
187      END IF
188
189   END SUBROUTINE obs_int_comm_2d
190
191   SUBROUTINE obs_int_comm_3d_global( kptsi, kptsj, kobs, kpk, kgrdi, kgrdj, &
192      &                               pval, pgval, kproc )
193      !!----------------------------------------------------------------------
194      !!                    ***  ROUTINE obs_int_comm_3d_global  ***
195      !!         
196      !! ** Purpose : Get 3D interpolation stencil (global version)
197      !!
198      !! ** Method  : On-demand communication where each processor send its
199      !!              list of (i,j) of points to all processors and receive
200      !!              the corresponding values
201      !!
202      !! ** Action  :
203      !!
204      !! History :
205      !!        !  08-02  (K. Mogensen)  Original code
206      !!----------------------------------------------------------------------
207      !! * Arguments
208      INTEGER, INTENT(IN) :: kptsi     ! Number of i horizontal points per stencil
209      INTEGER, INTENT(IN) :: kptsj     ! Number of j horizontal points per stencil
210      INTEGER, INTENT(IN) :: kobs      ! Local number of observations
211      INTEGER, INTENT(IN) :: kpk       ! Number of levels
212      INTEGER, DIMENSION(kptsi,kptsj,kobs), INTENT(IN) :: &
213         & kgrdi, &         ! i,j indices for each stencil
214         & kgrdj
215      INTEGER, OPTIONAL, DIMENSION(kptsi,kptsj,kobs), INTENT(IN) :: &
216         & kproc            ! Precomputed processor for each i,j,iobs points
217
218   !! * Control permutation of array indices
219!FTRANS CLEAR
220#  include "dom_oce_ftrans.h90"
221
222!! DCSE_NEMO: this style defeats ftrans
223!     REAL(KIND=wp), DIMENSION(jpi,jpj,kpk), INTENT(IN) ::&
224!        & pval             ! Local 3D array to extract data from
225!     REAL(KIND=wp), DIMENSION(kptsi,kptsj,kpk,kobs), INTENT(OUT) ::&
226!        & pgval            ! Stencil at each point
227
228!FTRANS pval :I :I :z
229      REAL(KIND=wp), INTENT(IN) ::&
230         & pval(jpi,jpj,kpk)             ! Local 3D array to extract data from
231
232!FTRANS pgval :I :I :z :
233      REAL(KIND=wp), INTENT(OUT) ::&
234         & pgval(kptsi,kptsj,kpk,kobs)   ! Stencil at each point
235
236      !! * Local declarations
237      REAL(KIND=wp), DIMENSION(:,:), ALLOCATABLE :: &
238         & zsend, &
239         & zrecv
240      INTEGER, DIMENSION(:), ALLOCATABLE :: &
241         & igrdij_send, &
242         & igrdij_recv
243      INTEGER, DIMENSION(kptsi,kptsj,kobs) :: &
244         & iorder, &
245         & iproc
246      INTEGER :: nplocal(jpnij)
247      INTEGER :: npglobal(jpnij)
248      INTEGER :: ji
249      INTEGER :: jj
250      INTEGER :: jk
251      INTEGER :: jp
252      INTEGER :: jobs
253      INTEGER :: it
254      INTEGER :: itot
255      INTEGER :: ii
256      INTEGER :: ij
257
258      ! Check valid points
259     
260      IF ( ( MAXVAL(kgrdi) > jpiglo ) .OR. ( MINVAL(kgrdi) < 1 ) .OR. &
261         & ( MAXVAL(kgrdj) > jpjglo ) .OR. ( MINVAL(kgrdj) < 1 ) ) THEN
262         
263         CALL ctl_stop( 'Error in obs_int_comm_3d_global', &
264            &           'Point outside global domain' )
265         
266      ENDIF
267
268      ! Count number of points on each processors
269
270      nplocal(:) = 0
271      IF (PRESENT(kproc)) THEN
272         iproc(:,:,:) = kproc(:,:,:)
273         DO jobs = 1, kobs
274            DO jj = 1, kptsj
275               DO ji = 1, kptsi
276                  nplocal(iproc(ji,jj,jobs)) = nplocal(iproc(ji,jj,jobs)) + 1
277               END DO
278            END DO
279         END DO
280      ELSE
281         DO jobs = 1, kobs
282            DO jj = 1, kptsj
283               DO ji = 1, kptsi
284                  iproc(ji,jj,jobs) = mppmap(kgrdi(ji,jj,jobs),&
285                     &                       kgrdj(ji,jj,jobs))
286                  nplocal(iproc(ji,jj,jobs)) = nplocal(iproc(ji,jj,jobs)) + 1
287               END DO
288            END DO
289         END DO
290      ENDIF
291
292      ! Send local number of points and receive points on current domain
293
294      CALL mpp_alltoall_int( 1, nplocal, npglobal )
295
296      ! Allocate message parsing workspace
297
298      itot = SUM(npglobal)
299
300      ALLOCATE( &
301         & igrdij_send(kptsi*kptsj*kobs*2), &
302         & igrdij_recv(itot*2),             &
303         & zsend(kpk,itot),                 &
304         & zrecv(kpk,kptsi*kptsj*kobs)      &
305         & )
306
307      ! Pack buffers for list of points
308
309      it = 0
310      DO jp = 1, jpnij 
311         DO jobs = 1, kobs
312            DO jj = 1, kptsj
313               DO ji = 1, kptsi
314                  IF ( iproc(ji,jj,jobs) == jp ) THEN
315                     it = it + 1
316                     iorder(ji,jj,jobs) = it
317                     igrdij_send(2*it-1) = kgrdi(ji,jj,jobs)
318                     igrdij_send(2*it  ) = kgrdj(ji,jj,jobs)
319                  ENDIF
320               END DO
321            END DO
322         END DO
323      END DO
324
325      ! Send and receive buffers for list of points
326
327      CALL mpp_alltoallv_int( igrdij_send, kptsi*kptsj*kobs*2, nplocal(:)*2, &
328         &                    igrdij_recv, itot*2, npglobal(:)*2 )
329
330      ! Pack interpolation data to be sent
331
332      DO ji = 1, itot
333         ii = mi1(igrdij_recv(2*ji-1))
334         ij = mj1(igrdij_recv(2*ji))
335         DO jk = 1, kpk
336            zsend(jk,ji) = pval(ii,ij,jk)
337         END DO
338      END DO
339
340      ! Re-adjust sizes
341
342      nplocal(:)  = kpk*nplocal(:)
343      npglobal(:) = kpk*npglobal(:)
344
345
346      ! Send and receive data for interpolation stencil
347
348      CALL mpp_alltoallv_real( zsend, kpk*itot,             npglobal, &
349         &                     zrecv, kpk*kptsi*kptsj*kobs, nplocal )
350
351      ! Copy the received data into output data structure
352     
353      DO jobs = 1, kobs
354         DO jj = 1, kptsj
355            DO ji = 1, kptsi
356               it = iorder(ji,jj,jobs) 
357               DO jk = 1, kpk
358                  pgval(ji,jj,jk,jobs) = zrecv(jk,it)
359               END DO
360            END DO
361         END DO
362      END DO
363
364      ! Deallocate message passing workspace
365
366      DEALLOCATE( &
367         & igrdij_send, &
368         & igrdij_recv, &
369         & zsend,       &
370         & zrecv        &
371         & )
372
373   END SUBROUTINE obs_int_comm_3d_global
374   
375   SUBROUTINE obs_int_comm_3d_local( kptsi, kptsj, kobs, kpk, kgrdi, kgrdj, &
376      &                              pval, pgval )
377      !!----------------------------------------------------------------------
378      !!                    ***  ROUTINE obs_int_comm_3d_global  ***
379      !!         
380      !! ** Purpose : Get 3D interpolation stencil (global version)
381      !!
382      !! ** Method  : On-demand communication where each processor send its
383      !!              list of (i,j) of points to all processors and receive
384      !!              the corresponding values
385      !!
386      !! ** Action  :
387      !!
388      !! History :
389      !!        !  08-02  (K. Mogensen)  Original code
390      !!----------------------------------------------------------------------
391      !! * Arguments
392      INTEGER, INTENT(IN) :: kptsi        ! Number of i horizontal points per stencil
393      INTEGER, INTENT(IN) :: kptsj        ! Number of j horizontal points per stencil
394      INTEGER, INTENT(IN) :: kobs         ! Local number of observations
395      INTEGER, INTENT(IN) :: kpk          ! Number of levels
396      INTEGER, DIMENSION(kptsi,kptsj,kobs), INTENT(IN) :: &
397         & kgrdi, &         ! i,j indices for each stencil
398         & kgrdj
399
400   !! * Control permutation of array indices
401!FTRANS CLEAR
402#  include "dom_oce_ftrans.h90"
403
404!! DCSE_NEMO: this style defeats ftrans
405!     REAL(KIND=wp), DIMENSION(jpi,jpj,kpk), INTENT(IN) ::&
406!        & pval             ! Local 3D array to extract data from
407!     REAL(KIND=wp), DIMENSION(kptsi,kptsj,kpk,kobs), INTENT(OUT) ::&
408!        & pgval            ! Stencil at each point
409
410!FTRANS pval :I :I :z
411      REAL(KIND=wp), INTENT(IN)  ::&
412         & pval(jpi,jpj,kpk)             ! Local 3D array to extract data from
413!FTRANS pgval :I :I :z :
414      REAL(KIND=wp), INTENT(OUT) ::&
415         & pgval(kptsi,kptsj,kpk,kobs)   ! Stencil at each point
416
417      !! * Local declarations
418      INTEGER ::  ji
419      INTEGER ::  jj
420      INTEGER ::  jk
421      INTEGER ::  jobs
422
423      ! Check valid points
424
425      IF ( ( MAXVAL(kgrdi) > jpi ) .OR. ( MINVAL(kgrdi) < 1 ) .OR. &
426         & ( MAXVAL(kgrdj) > jpj ) .OR. ( MINVAL(kgrdj) < 1 ) ) THEN
427         
428         CALL ctl_stop( 'Error in obs_int_comm_3d_local', &
429            &           'Point outside local domain' )
430         
431      ENDIF
432
433      ! Copy local data
434
435      DO jobs = 1, kobs
436         DO jj = 1, kptsj
437            DO ji = 1, kptsi
438               DO jk = 1, kpk
439                  pgval(ji,jj,jk,jobs) = &
440                     &            pval(kgrdi(ji,jj,jobs),kgrdj(ji,jj,jobs),jk)
441               END DO
442            END DO
443         END DO
444      END DO
445
446   END SUBROUTINE obs_int_comm_3d_local
447
448END MODULE obs_inter_sup
449 
Note: See TracBrowser for help on using the repository browser.