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_rot_vel.F90 in branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

source: branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_rot_vel.F90 @ 13393

Last change on this file since 13393 was 13393, checked in by mattmartin, 4 years ago

Include surface velocity observation operator. See ticket ​https://code.metoffice.gov.uk/trac/utils/ticket/376.
File size: 15.6 KB
Line 
1MODULE obs_rot_vel
2   !!======================================================================
3   !!                       ***  MODULE obs_rot_vel  ***
4   !! Observation diagnostics: Read the velocity profile observations
5   !!======================================================================
6
7   !!----------------------------------------------------------------------
8   !!   obs_rotvel_pro :  Rotate profile velocity data into N-S,E-W directions
9   !!   obs_rotvel_surf : Rotate surface velocity data into N-S,E-W directions   
10   !!----------------------------------------------------------------------
11   !! * Modules used   
12   USE wrk_nemo                 ! Memory Allocation
13   USE par_kind                 ! Precision variables
14   USE par_oce                  ! Ocean parameters
15   USE in_out_manager           ! I/O manager
16   USE dom_oce                  ! Ocean space and time domain variables
17   USE obs_grid                 ! Grid search
18   USE obs_utils                ! For error handling
19   USE obs_profiles_def         ! Profile definitions
20   USE obs_surf_def             ! Surface definitions   
21   USE obs_inter_h2d            ! Horizontal interpolation
22   USE obs_inter_sup            ! MPP support routines for interpolation
23   USE geo2ocean                ! Rotation of vectors
24   USE obs_fbm                  ! Feedback definitions
25
26   IMPLICIT NONE
27
28   !! * Routine accessibility
29   PRIVATE
30
31   PUBLIC obs_rotvel_pro, &     ! Rotate the profile velocity observations
32      &   obs_rotvel_surf       ! Rotate the surface velocity observations
33     
34   !!----------------------------------------------------------------------
35   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
36   !! $Id$
37   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
38   !!----------------------------------------------------------------------
39
40CONTAINS
41
42   SUBROUTINE obs_rotvel_pro( profdata, k2dint, pu, pv )
43      !!---------------------------------------------------------------------
44      !!
45      !!                   *** ROUTINE obs_rea_pro_dri ***
46      !!
47      !! ** Purpose : Rotate velocity data into N-S,E-W directorions
48      !!
49      !! ** Method  : Interpolation of geo2ocean coefficients on U,V grid
50      !!              to observation point followed by a similar computations
51      !!              as in geo2ocean.
52      !!
53      !! ** Action  : Review if there is a better way to do this.
54      !!
55      !! References :
56      !!
57      !! History : 
58      !!      ! :  2009-02 (K. Mogensen) : New routine
59      !!----------------------------------------------------------------------
60      !! * Modules used
61      !! * Arguments
62      TYPE(obs_prof), INTENT(INOUT) :: profdata    ! Profile data to be read
63      INTEGER, INTENT(IN) :: k2dint     ! Horizontal interpolation methed
64      REAL(wp), DIMENSION(*) :: &
65         & pu, &
66         & pv
67      !! * Local declarations
68      REAL(wp), DIMENSION(2,2,1) :: zweig
69      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
70         & zmasku, &
71         & zmaskv, &
72         & zcoslu, &
73         & zsinlu, &
74         & zcoslv, &
75         & zsinlv, &
76         & zglamu, &
77         & zgphiu, &
78         & zglamv, &
79         & zgphiv
80      REAL(wp), DIMENSION(1) :: &
81         & zsinu, &
82         & zcosu, &
83         & zsinv, &
84         & zcosv
85      REAL(wp) :: zsin
86      REAL(wp) :: zcos
87      REAL(wp), DIMENSION(1) :: zobsmask
88      REAL(wp), POINTER, DIMENSION(:,:) :: zsingu,zcosgu,zsingv,zcosgv
89      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
90         & igrdiu, &
91         & igrdju, &
92         & igrdiv, &
93         & igrdjv
94      INTEGER :: ji
95      INTEGER :: jk
96
97      CALL wrk_alloc(jpi,jpj,zsingu,zcosgu,zsingv,zcosgv) 
98
99      !-----------------------------------------------------------------------
100      ! Allocate data for message parsing and interpolation
101      !-----------------------------------------------------------------------
102
103      ALLOCATE( &
104         & igrdiu(2,2,profdata%nprof), &
105         & igrdju(2,2,profdata%nprof), &
106         & zglamu(2,2,profdata%nprof), &
107         & zgphiu(2,2,profdata%nprof), &
108         & zmasku(2,2,profdata%nprof), &
109         & zcoslu(2,2,profdata%nprof), &
110         & zsinlu(2,2,profdata%nprof), &
111         & igrdiv(2,2,profdata%nprof), &
112         & igrdjv(2,2,profdata%nprof), &
113         & zglamv(2,2,profdata%nprof), &
114         & zgphiv(2,2,profdata%nprof), &
115         & zmaskv(2,2,profdata%nprof), &
116         & zcoslv(2,2,profdata%nprof), &
117         & zsinlv(2,2,profdata%nprof)  &
118         & )
119
120      !-----------------------------------------------------------------------
121      ! Receive the angles on the U and V grids.
122      !-----------------------------------------------------------------------
123
124      CALL obs_rot( zsingu, zcosgu, zsingv, zcosgv )
125
126      DO ji = 1, profdata%nprof
127         igrdiu(1,1,ji) = profdata%mi(ji,1)-1
128         igrdju(1,1,ji) = profdata%mj(ji,1)-1
129         igrdiu(1,2,ji) = profdata%mi(ji,1)-1
130         igrdju(1,2,ji) = profdata%mj(ji,1)
131         igrdiu(2,1,ji) = profdata%mi(ji,1)
132         igrdju(2,1,ji) = profdata%mj(ji,1)-1
133         igrdiu(2,2,ji) = profdata%mi(ji,1)
134         igrdju(2,2,ji) = profdata%mj(ji,1)
135         igrdiv(1,1,ji) = profdata%mi(ji,2)-1
136         igrdjv(1,1,ji) = profdata%mj(ji,2)-1
137         igrdiv(1,2,ji) = profdata%mi(ji,2)-1
138         igrdjv(1,2,ji) = profdata%mj(ji,2)
139         igrdiv(2,1,ji) = profdata%mi(ji,2)
140         igrdjv(2,1,ji) = profdata%mj(ji,2)-1
141         igrdiv(2,2,ji) = profdata%mi(ji,2)
142         igrdjv(2,2,ji) = profdata%mj(ji,2)
143      END DO
144
145      CALL obs_int_comm_2d( 2, 2, profdata%nprof, jpi, jpj, igrdiu, igrdju, &
146         &                  glamu, zglamu )
147      CALL obs_int_comm_2d( 2, 2, profdata%nprof, jpi, jpj, igrdiu, igrdju, &
148         &                  gphiu, zgphiu )
149      CALL obs_int_comm_2d( 2, 2, profdata%nprof, jpi, jpj, igrdiu, igrdju, &
150         &                  umask(:,:,1), zmasku )
151      CALL obs_int_comm_2d( 2, 2, profdata%nprof, jpi, jpj, igrdiu, igrdju, &
152         &                  zsingu, zsinlu )
153      CALL obs_int_comm_2d( 2, 2, profdata%nprof, jpi, jpj, igrdiu, igrdju, &
154         &                  zcosgu, zcoslu )
155      CALL obs_int_comm_2d( 2, 2, profdata%nprof, jpi, jpj, igrdiv, igrdjv, &
156         &                  glamv, zglamv )
157      CALL obs_int_comm_2d( 2, 2, profdata%nprof, jpi, jpj, igrdiv, igrdjv, &
158         &                  gphiv, zgphiv )
159      CALL obs_int_comm_2d( 2, 2, profdata%nprof, jpi, jpj, igrdiv, igrdjv, &
160         &                  vmask(:,:,1), zmaskv )
161      CALL obs_int_comm_2d( 2, 2, profdata%nprof, jpi, jpj, igrdiv, igrdjv, &
162         &                  zsingv, zsinlv )
163      CALL obs_int_comm_2d( 2, 2, profdata%nprof, jpi, jpj, igrdiv, igrdjv, &
164         &                  zcosgv, zcoslv )
165
166      DO ji = 1, profdata%nprof
167           
168         CALL obs_int_h2d_init( 1, 1, k2dint, &
169            &                   profdata%rlam(ji), profdata%rphi(ji), &
170            &                   zglamu(:,:,ji), zgphiu(:,:,ji), &
171            &                   zmasku(:,:,ji), zweig, zobsmask )
172         
173         CALL obs_int_h2d( 1, 1, zweig, zsinlu(:,:,ji),  zsinu )
174
175         CALL obs_int_h2d( 1, 1, zweig, zcoslu(:,:,ji),  zcosu )
176
177         CALL obs_int_h2d_init( 1, 1, k2dint, &
178            &                   profdata%rlam(ji), profdata%rphi(ji), &
179            &                   zglamv(:,:,ji), zgphiv(:,:,ji), &
180            &                   zmaskv(:,:,ji), zweig, zobsmask )
181         
182         CALL obs_int_h2d( 1, 1, zweig, zsinlv(:,:,ji),  zsinv )
183
184         CALL obs_int_h2d( 1, 1, zweig, zcoslv(:,:,ji),  zcosv )
185
186         ! Assume that the angle at observation point is the
187         ! mean of u and v cosines/sines
188
189         zcos = 0.5_wp * ( zcosu(1) + zcosv(1) )
190         zsin = 0.5_wp * ( zsinu(1) + zsinv(1) )
191         
192         IF ( ( profdata%npvsta(ji,1) /= profdata%npvsta(ji,2) ) .OR. &
193            & ( profdata%npvend(ji,1) /= profdata%npvend(ji,2) ) ) THEN
194            CALL fatal_error( 'Different number of U and V observations '// &
195               'in a profile in obs_rotvel', __LINE__ )
196         ENDIF
197
198         DO jk = profdata%npvsta(ji,1), profdata%npvend(ji,1)
199            IF ( ( profdata%var(1)%vmod(jk) /= fbrmdi ) .AND. &
200               & ( profdata%var(2)%vmod(jk) /= fbrmdi ) ) THEN
201               pu(jk) = profdata%var(1)%vmod(jk) * zcos - &
202                  &     profdata%var(2)%vmod(jk) * zsin
203               pv(jk) = profdata%var(2)%vmod(jk) * zcos + &
204                  &     profdata%var(1)%vmod(jk) * zsin
205            ELSE
206               pu(jk) = fbrmdi
207               pv(jk) = fbrmdi
208            ENDIF
209
210         END DO
211
212      END DO
213     
214      DEALLOCATE( &
215         & igrdiu, &
216         & igrdju, &
217         & zglamu, &
218         & zgphiu, &
219         & zmasku, &
220         & zcoslu, &
221         & zsinlu, &
222         & igrdiv, &
223         & igrdjv, &
224         & zglamv, &
225         & zgphiv, &
226         & zmaskv, &
227         & zcoslv, &
228         & zsinlv  &
229         & )
230
231      CALL wrk_dealloc(jpi,jpj,zsingu,zcosgu,zsingv,zcosgv) 
232
233   END SUBROUTINE obs_rotvel_pro
234
235   SUBROUTINE obs_rotvel_surf( surfdata, k2dint, pu, pv )
236      !!---------------------------------------------------------------------
237      !!
238      !!                   *** ROUTINE obs_rotvel_surf ***
239      !!
240      !! ** Purpose : Rotate surface velocity data into N-S,E-W directorions
241      !!
242      !! ** Method  : Interpolation of geo2ocean coefficients on U,V grid
243      !!              to observation point followed by a similar computations
244      !!              as in geo2ocean.
245      !!
246      !! ** Action  : Review if there is a better way to do this.
247      !!
248      !! References :
249      !!
250      !! History : 
251      !!      ! :  2009-02 (K. Mogensen) : New routine
252      !!----------------------------------------------------------------------
253      !! * Modules used
254      !! * Arguments
255      TYPE(obs_surf), INTENT(INOUT) :: surfdata    ! Surface data to be read
256      INTEGER, INTENT(IN) :: k2dint     ! Horizontal interpolation methed
257      REAL(wp), DIMENSION(*) :: &
258         & pu, &
259         & pv
260      !! * Local declarations
261      REAL(wp), DIMENSION(2,2,1) :: zweig
262      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
263         & zmasku, &
264         & zmaskv, &
265         & zcoslu, &
266         & zsinlu, &
267         & zcoslv, &
268         & zsinlv, &
269         & zglamu, &
270         & zgphiu, &
271         & zglamv, &
272         & zgphiv
273      REAL(wp), DIMENSION(1) :: &
274         & zsinu, &
275         & zcosu, &
276         & zsinv, &
277         & zcosv
278      REAL(wp) :: zsin
279      REAL(wp) :: zcos
280      REAL(wp), DIMENSION(1) :: zobsmask
281      REAL(wp), POINTER, DIMENSION(:,:) :: zsingu,zcosgu,zsingv,zcosgv
282      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
283         & igrdiu, &
284         & igrdju, &
285         & igrdiv, &
286         & igrdjv
287      INTEGER :: ji
288      INTEGER :: jk
289
290      CALL wrk_alloc(jpi,jpj,zsingu,zcosgu,zsingv,zcosgv) 
291
292      !-----------------------------------------------------------------------
293      ! Allocate data for message parsing and interpolation
294      !-----------------------------------------------------------------------
295
296      ALLOCATE( &
297         & igrdiu(2,2,surfdata%nsurf), &
298         & igrdju(2,2,surfdata%nsurf), &
299         & zglamu(2,2,surfdata%nsurf), &
300         & zgphiu(2,2,surfdata%nsurf), &
301         & zmasku(2,2,surfdata%nsurf), &
302         & zcoslu(2,2,surfdata%nsurf), &
303         & zsinlu(2,2,surfdata%nsurf), &
304         & igrdiv(2,2,surfdata%nsurf), &
305         & igrdjv(2,2,surfdata%nsurf), &
306         & zglamv(2,2,surfdata%nsurf), &
307         & zgphiv(2,2,surfdata%nsurf), &
308         & zmaskv(2,2,surfdata%nsurf), &
309         & zcoslv(2,2,surfdata%nsurf), &
310         & zsinlv(2,2,surfdata%nsurf)  &
311         & )
312
313      !-----------------------------------------------------------------------
314      ! Receive the angles on the U and V grids.
315      !-----------------------------------------------------------------------
316
317      CALL obs_rot( zsingu, zcosgu, zsingv, zcosgv )
318
319      DO ji = 1, surfdata%nsurf
320         igrdiu(1,1,ji) = surfdata%mi(ji,1)-1
321         igrdju(1,1,ji) = surfdata%mj(ji,1)-1
322         igrdiu(1,2,ji) = surfdata%mi(ji,1)-1
323         igrdju(1,2,ji) = surfdata%mj(ji,1)
324         igrdiu(2,1,ji) = surfdata%mi(ji,1)
325         igrdju(2,1,ji) = surfdata%mj(ji,1)-1
326         igrdiu(2,2,ji) = surfdata%mi(ji,1)
327         igrdju(2,2,ji) = surfdata%mj(ji,1)
328         igrdiv(1,1,ji) = surfdata%mi(ji,2)-1
329         igrdjv(1,1,ji) = surfdata%mj(ji,2)-1
330         igrdiv(1,2,ji) = surfdata%mi(ji,2)-1
331         igrdjv(1,2,ji) = surfdata%mj(ji,2)
332         igrdiv(2,1,ji) = surfdata%mi(ji,2)
333         igrdjv(2,1,ji) = surfdata%mj(ji,2)-1
334         igrdiv(2,2,ji) = surfdata%mi(ji,2)
335         igrdjv(2,2,ji) = surfdata%mj(ji,2)
336      END DO
337
338      CALL obs_int_comm_2d( 2, 2, surfdata%nsurf, jpi, jpj, igrdiu, igrdju, &
339         &                  glamu, zglamu )
340      CALL obs_int_comm_2d( 2, 2, surfdata%nsurf, jpi, jpj, igrdiu, igrdju, &
341         &                  gphiu, zgphiu )
342      CALL obs_int_comm_2d( 2, 2, surfdata%nsurf, jpi, jpj, igrdiu, igrdju, &
343         &                  umask(:,:,1), zmasku )
344      CALL obs_int_comm_2d( 2, 2, surfdata%nsurf, jpi, jpj, igrdiu, igrdju, &
345         &                  zsingu, zsinlu )
346      CALL obs_int_comm_2d( 2, 2, surfdata%nsurf, jpi, jpj, igrdiu, igrdju, &
347         &                  zcosgu, zcoslu )
348      CALL obs_int_comm_2d( 2, 2, surfdata%nsurf, jpi, jpj, igrdiv, igrdjv, &
349         &                  glamv, zglamv )
350      CALL obs_int_comm_2d( 2, 2, surfdata%nsurf, jpi, jpj, igrdiv, igrdjv, &
351         &                  gphiv, zgphiv )
352      CALL obs_int_comm_2d( 2, 2, surfdata%nsurf, jpi, jpj, igrdiv, igrdjv, &
353         &                  vmask(:,:,1), zmaskv )
354      CALL obs_int_comm_2d( 2, 2, surfdata%nsurf, jpi, jpj, igrdiv, igrdjv, &
355         &                  zsingv, zsinlv )
356      CALL obs_int_comm_2d( 2, 2, surfdata%nsurf, jpi, jpj, igrdiv, igrdjv, &
357         &                  zcosgv, zcoslv )
358
359      DO ji = 1, surfdata%nsurf
360           
361         CALL obs_int_h2d_init( 1, 1, k2dint, &
362            &                   surfdata%rlam(ji), surfdata%rphi(ji), &
363            &                   zglamu(:,:,ji), zgphiu(:,:,ji), &
364            &                   zmasku(:,:,ji), zweig, zobsmask )
365         
366         CALL obs_int_h2d( 1, 1, zweig, zsinlu(:,:,ji),  zsinu )
367
368         CALL obs_int_h2d( 1, 1, zweig, zcoslu(:,:,ji),  zcosu )
369
370         CALL obs_int_h2d_init( 1, 1, k2dint, &
371            &                   surfdata%rlam(ji), surfdata%rphi(ji), &
372            &                   zglamv(:,:,ji), zgphiv(:,:,ji), &
373            &                   zmaskv(:,:,ji), zweig, zobsmask )
374         
375         CALL obs_int_h2d( 1, 1, zweig, zsinlv(:,:,ji),  zsinv )
376
377         CALL obs_int_h2d( 1, 1, zweig, zcoslv(:,:,ji),  zcosv )
378
379         ! Assume that the angle at observation point is the
380         ! mean of u and v cosines/sines
381
382         zcos = 0.5_wp * ( zcosu(1) + zcosv(1) )
383         zsin = 0.5_wp * ( zsinu(1) + zsinv(1) )
384
385         IF ( ( surfdata%rmod(ji,1) /= fbrmdi ) .AND. &
386            & ( surfdata%rmod(ji,2) /= fbrmdi ) ) THEN
387            pu(ji) = surfdata%rmod(ji,1) * zcos - &
388               &     surfdata%rmod(ji,2) * zsin
389            pv(ji) = surfdata%rmod(ji,2) * zcos + &
390               &     surfdata%rmod(ji,1) * zsin
391         ELSE
392            pu(ji) = fbrmdi
393            pv(ji) = fbrmdi
394         ENDIF
395
396
397      END DO
398     
399      DEALLOCATE( &
400         & igrdiu, &
401         & igrdju, &
402         & zglamu, &
403         & zgphiu, &
404         & zmasku, &
405         & zcoslu, &
406         & zsinlu, &
407         & igrdiv, &
408         & igrdjv, &
409         & zglamv, &
410         & zgphiv, &
411         & zmaskv, &
412         & zcoslv, &
413         & zsinlv  &
414         & )
415
416      CALL wrk_dealloc(jpi,jpj,zsingu,zcosgu,zsingv,zcosgv) 
417
418   END SUBROUTINE obs_rotvel_surf
419
420END MODULE obs_rot_vel
Note: See TracBrowser for help on using the repository browser.