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

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/OBS/obs_inter_sup.F90 @ 2245

Last change on this file since 2245 was 2128, checked in by rfurner, 14 years ago

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