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

Last change on this file since 2287 was 2287, checked in by smasson, 14 years ago

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