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

source: branches/dev_1784_OBS/NEMO/OPA_SRC/OBS/obs_inter_sup.F90 @ 2001

Last change on this file since 2001 was 2001, checked in by djlea, 14 years ago

Adding observation operator code

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