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 @ 10246

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

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