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_oper.F90 in trunk/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

source: trunk/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90 @ 2576

Last change on this file since 2576 was 2576, checked in by djlea, 13 years ago

Bug fix to velocity observation operator

  • Property svn:keywords set to Id
File size: 48.3 KB
Line 
1MODULE obs_oper
2   !!======================================================================
3   !!                       ***  MODULE obs_oper  ***
4   !! Observation diagnostics: Observation operators for various observation
5   !!                          types
6   !!======================================================================
7
8   !!----------------------------------------------------------------------
9   !!   obs_pro_opt :    Compute the model counterpart of temperature and
10   !!                    salinity observations from profiles
11   !!   obs_sla_opt :    Compute the model counterpart of sea level anomaly
12   !!                    observations
13   !!   obs_sst_opt :    Compute the model counterpart of sea surface temperature
14   !!                    observations
15   !!   obs_sss_opt :    Compute the model counterpart of sea surface salinity
16   !!                    observations
17   !!   obs_seaice_opt : Compute the model counterpart of sea ice concentration
18   !!                    observations
19   !!
20   !!   obs_vel_opt :    Compute the model counterpart of zonal and meridional
21   !!                    components of velocity from observations.
22   !!----------------------------------------------------------------------
23
24   !! * Modules used   
25   USE par_kind, ONLY : &         ! Precision variables
26      & wp
27   USE in_out_manager             ! I/O manager
28   USE obs_inter_sup              ! Interpolation support
29   USE obs_inter_h2d, ONLY : &    ! Horizontal interpolation to the observation pt
30      & obs_int_h2d, &
31      & obs_int_h2d_init
32   USE obs_inter_z1d, ONLY : &    ! Vertical interpolation to the observation pt
33      & obs_int_z1d,    &
34      & obs_int_z1d_spl
35   USE obs_const,  ONLY :     &
36      & obfillflt      ! Fillvalue   
37   USE dom_oce,       ONLY : &
38      & glamt, glamu, glamv, &
39      & gphit, gphiu, gphiv
40
41   IMPLICIT NONE
42
43   !! * Routine accessibility
44   PRIVATE
45
46   PUBLIC obs_pro_opt, &  ! Compute the model counterpart of profile observations
47      &   obs_sla_opt, &  ! Compute the model counterpart of SLA observations
48      &   obs_sst_opt, &  ! Compute the model counterpart of SST observations
49      &   obs_sss_opt, &  ! Compute the model counterpart of SSS observations
50      &   obs_seaice_opt, &
51      &   obs_vel_opt     ! Compute the model counterpart of velocity profile data
52
53   INTEGER, PARAMETER, PUBLIC :: imaxavtypes = 20 ! Max number of daily avgd obs types
54
55   !!----------------------------------------------------------------------
56   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
57   !! $Id$
58   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
59   !!----------------------------------------------------------------------
60
61CONTAINS
62
63   SUBROUTINE obs_pro_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, &
64      &                    ptn, psn, pgdept, ptmask, k1dint, k2dint, &
65      &                    kdailyavtypes )
66      !!-----------------------------------------------------------------------
67      !!
68      !!                     ***  ROUTINE obs_pro_opt  ***
69      !!
70      !! ** Purpose : Compute the model counterpart of profiles
71      !!              data by interpolating from the model grid to the
72      !!              observation point.
73      !!
74      !! ** Method  : Linearly interpolate to each observation point using
75      !!              the model values at the corners of the surrounding grid box.
76      !!
77      !!    First, a vertical profile of horizontally interpolated model
78      !!    now temperatures is computed at the obs (lon, lat) point.
79      !!    Several horizontal interpolation schemes are available:
80      !!        - distance-weighted (great circle) (k2dint = 0)
81      !!        - distance-weighted (small angle)  (k2dint = 1)
82      !!        - bilinear (geographical grid)     (k2dint = 2)
83      !!        - bilinear (quadrilateral grid)    (k2dint = 3)
84      !!        - polynomial (quadrilateral grid)  (k2dint = 4)
85      !!
86      !!    Next, the vertical temperature profile is interpolated to the
87      !!    data depth points. Two vertical interpolation schemes are
88      !!    available:
89      !!        - linear       (k1dint = 0)
90      !!        - Cubic spline (k1dint = 1)
91      !!
92      !!    For the cubic spline the 2nd derivative of the interpolating
93      !!    polynomial is computed before entering the vertical interpolation
94      !!    routine.
95      !!
96      !!    For ENACT moored buoy data (e.g., TAO), the model equivalent is
97      !!    a daily mean model temperature field. So, we first compute
98      !!    the mean, then interpolate only at the end of the day.
99      !!
100      !!    Note: the in situ temperature observations must be converted
101      !!    to potential temperature (the model variable) prior to
102      !!    assimilation.
103      !!??????????????????????????????????????????????????????????????
104      !!    INCLUDE POTENTIAL TEMP -> IN SITU TEMP IN OBS OPERATOR???
105      !!??????????????????????????????????????????????????????????????
106      !!
107      !! ** Action  :
108      !!
109      !! History :
110      !!      ! 97-11 (A. Weaver, S. Ricci, N. Daget)
111      !!      ! 06-03 (G. Smith) NEMOVAR migration
112      !!      ! 06-10 (A. Weaver) Cleanup
113      !!      ! 07-01 (K. Mogensen) Merge of temperature and salinity
114      !!      ! 07-03 (K. Mogensen) General handling of profiles
115      !!-----------------------------------------------------------------------
116 
117      !! * Modules used
118      USE obs_profiles_def ! Definition of storage space for profile obs.
119
120      IMPLICIT NONE
121
122      !! * Arguments
123      TYPE(obs_prof), INTENT(INOUT) :: prodatqc  ! Subset of profile data not failing screening
124      INTEGER, INTENT(IN) :: kt        ! Time step
125      INTEGER, INTENT(IN) :: kpi       ! Model grid parameters
126      INTEGER, INTENT(IN) :: kpj
127      INTEGER, INTENT(IN) :: kpk
128      INTEGER, INTENT(IN) :: kit000    ! Number of the first time step
129                                       !   (kit000-1 = restart time)
130      INTEGER, INTENT(IN) :: k1dint    ! Vertical interpolation type (see header)
131      INTEGER, INTENT(IN) :: k2dint    ! Horizontal interpolation type (see header)
132      INTEGER, INTENT(IN) :: kdaystp   ! Number of time steps per day                   
133      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &
134         & ptn,    &    ! Model temperature field
135         & psn,    &    ! Model salinity field
136         & ptmask       ! Land-sea mask
137      REAL(KIND=wp), INTENT(IN), DIMENSION(kpk) :: &
138         & pgdept       ! Model array of depth levels
139      INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &
140         & kdailyavtypes! Types for daily averages
141      !! * Local declarations
142      INTEGER ::   ji
143      INTEGER ::   jj
144      INTEGER ::   jk
145      INTEGER ::   jobs
146      INTEGER ::   inrc
147      INTEGER ::   ipro
148      INTEGER ::   idayend
149      INTEGER ::   ista
150      INTEGER ::   iend
151      INTEGER ::   iobs
152      INTEGER, DIMENSION(imaxavtypes) :: &
153         & idailyavtypes
154      REAL(KIND=wp) :: zlam
155      REAL(KIND=wp) :: zphi
156      REAL(KIND=wp) :: zdaystp
157      REAL(KIND=wp), DIMENSION(kpk) :: &
158         & zobsmask, &
159         & zobsk,    &
160         & zobs2k
161      REAL(KIND=wp), DIMENSION(2,2,kpk) :: &
162         & zweig
163      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
164         & zmask, &
165         & zintt, &
166         & zints, &
167         & zinmt, &
168         & zinms
169      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
170         & zglam, &
171         & zgphi
172      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
173         & igrdi, &
174         & igrdj
175
176      !------------------------------------------------------------------------
177      ! Local initialization
178      !------------------------------------------------------------------------
179      ! ... Record and data counters
180      inrc = kt - kit000 + 2
181      ipro = prodatqc%npstp(inrc)
182 
183      ! Daily average types
184      IF ( PRESENT(kdailyavtypes) ) THEN
185         idailyavtypes(:) = kdailyavtypes(:)
186      ELSE
187         idailyavtypes(:) = -1
188      ENDIF
189
190      ! Initialize daily mean for first timestep
191      idayend = MOD( kt - kit000 + 1, kdaystp )
192
193      ! Added kt == 0 test to catch restart case
194      IF ( idayend == 1 .OR. kt == 0) THEN
195         IF (lwp) WRITE(numout,*) 'Reset prodatqc%vdmean on time-step: ',kt
196         DO jk = 1, jpk
197            DO jj = 1, jpj
198               DO ji = 1, jpi
199                  prodatqc%vdmean(ji,jj,jk,1) = 0.0
200                  prodatqc%vdmean(ji,jj,jk,2) = 0.0
201               END DO
202            END DO
203         END DO
204      ENDIF
205
206      DO jk = 1, jpk
207         DO jj = 1, jpj
208            DO ji = 1, jpi
209               ! Increment the temperature field for computing daily mean
210               prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &
211                  &                        + ptn(ji,jj,jk)
212               ! Increment the salinity field for computing daily mean
213               prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &
214                  &                        + psn(ji,jj,jk)
215            END DO
216         END DO
217      END DO
218   
219      ! Compute the daily mean at the end of day
220      zdaystp = 1.0 / REAL( kdaystp )
221      IF ( idayend == 0 ) THEN
222         DO jk = 1, jpk
223            DO jj = 1, jpj
224               DO ji = 1, jpi
225                  prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &
226                     &                        * zdaystp
227                  prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &
228                  &                           * zdaystp
229               END DO
230            END DO
231         END DO
232      ENDIF
233
234      ! Get the data for interpolation
235      ALLOCATE( &
236         & igrdi(2,2,ipro),      &
237         & igrdj(2,2,ipro),      &
238         & zglam(2,2,ipro),      &
239         & zgphi(2,2,ipro),      &
240         & zmask(2,2,kpk,ipro),  &
241         & zintt(2,2,kpk,ipro),  &
242         & zints(2,2,kpk,ipro)   &
243         & )
244
245      DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
246         iobs = jobs - prodatqc%nprofup
247         igrdi(1,1,iobs) = prodatqc%mi(jobs,1)-1
248         igrdj(1,1,iobs) = prodatqc%mj(jobs,1)-1
249         igrdi(1,2,iobs) = prodatqc%mi(jobs,1)-1
250         igrdj(1,2,iobs) = prodatqc%mj(jobs,1)
251         igrdi(2,1,iobs) = prodatqc%mi(jobs,1)
252         igrdj(2,1,iobs) = prodatqc%mj(jobs,1)-1
253         igrdi(2,2,iobs) = prodatqc%mi(jobs,1)
254         igrdj(2,2,iobs) = prodatqc%mj(jobs,1)
255      END DO
256
257      CALL obs_int_comm_2d( 2, 2, ipro, igrdi, igrdj, glamt, zglam )
258      CALL obs_int_comm_2d( 2, 2, ipro, igrdi, igrdj, gphit, zgphi )
259      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, ptmask,zmask )
260      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, ptn,   zintt )
261      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, psn,   zints )
262
263      ! At the end of the day also get interpolated means
264      IF ( idayend == 0 ) THEN
265
266         ALLOCATE( &
267            & zinmt(2,2,kpk,ipro),  &
268            & zinms(2,2,kpk,ipro)   &
269            & )
270
271         CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, &
272            &                  prodatqc%vdmean(:,:,:,1), zinmt )
273         CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, &
274            &                  prodatqc%vdmean(:,:,:,2), zinms )
275
276      ENDIF
277
278      DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
279
280         iobs = jobs - prodatqc%nprofup
281
282         IF ( kt /= prodatqc%mstp(jobs) ) THEN
283           
284            IF(lwp) THEN
285               WRITE(numout,*)
286               WRITE(numout,*) ' E R R O R : Observation',              &
287                  &            ' time step is not consistent with the', &
288                  &            ' model time step'
289               WRITE(numout,*) ' ========='
290               WRITE(numout,*)
291               WRITE(numout,*) ' Record  = ', jobs,                    &
292                  &            ' kt      = ', kt,                      &
293                  &            ' mstp    = ', prodatqc%mstp(jobs), &
294                  &            ' ntyp    = ', prodatqc%ntyp(jobs)
295            ENDIF
296            CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' )
297         ENDIF
298         
299         zlam = prodatqc%rlam(jobs)
300         zphi = prodatqc%rphi(jobs)
301         
302         ! Horizontal weights and vertical mask
303
304         IF ( ( prodatqc%npvend(jobs,1) > 0 ) .OR. &
305            & ( prodatqc%npvend(jobs,2) > 0 ) ) THEN
306
307            CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi,     &
308               &                   zglam(:,:,iobs), zgphi(:,:,iobs), &
309               &                   zmask(:,:,:,iobs), zweig, zobsmask )
310
311         ENDIF
312
313         IF ( prodatqc%npvend(jobs,1) > 0 ) THEN
314
315            zobsk(:) = obfillflt
316
317       IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN
318
319               IF ( idayend == 0 )  THEN
320                 
321                  ! Daily averaged moored buoy (MRB) data
322                 
323                  CALL obs_int_h2d( kpk, kpk,      &
324                     &              zweig, zinmt(:,:,:,iobs), zobsk )
325                 
326                 
327               ELSE
328               
329                  CALL ctl_stop( ' A nonzero' //     &
330                     &           ' number of profile T BUOY data should' // &
331                     &           ' only occur at the end of a given day' )
332
333               ENDIF
334         
335            ELSE 
336               
337               ! Point data
338
339               CALL obs_int_h2d( kpk, kpk,      &
340                  &              zweig, zintt(:,:,:,iobs), zobsk )
341
342            ENDIF
343
344            !-------------------------------------------------------------
345            ! Compute vertical second-derivative of the interpolating
346            ! polynomial at obs points
347            !-------------------------------------------------------------
348           
349            IF ( k1dint == 1 ) THEN
350               CALL obs_int_z1d_spl( kpk, zobsk, zobs2k,   &
351                  &                  pgdept, zobsmask )
352            ENDIF
353           
354            !-----------------------------------------------------------------
355            !  Vertical interpolation to the observation point
356            !-----------------------------------------------------------------
357            ista = prodatqc%npvsta(jobs,1)
358            iend = prodatqc%npvend(jobs,1)
359            CALL obs_int_z1d( kpk,                &
360               & prodatqc%var(1)%mvk(ista:iend),  &
361               & k1dint, iend - ista + 1,         &
362               & prodatqc%var(1)%vdep(ista:iend), &
363               & zobsk, zobs2k,                   &
364               & prodatqc%var(1)%vmod(ista:iend), &
365               & pgdept, zobsmask )
366
367         ENDIF
368
369         IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
370
371            zobsk(:) = obfillflt
372
373            IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN
374
375               IF ( idayend == 0 )  THEN
376
377                  ! Daily averaged moored buoy (MRB) data
378                 
379                  CALL obs_int_h2d( kpk, kpk,      &
380                     &              zweig, zinms(:,:,:,iobs), zobsk )
381                 
382               ELSE
383
384                  CALL ctl_stop( ' A nonzero' //     &
385                     &           ' number of profile S BUOY data should' // &
386                     &           ' only occur at the end of a given day' )
387
388               ENDIF
389
390            ELSE
391               
392               ! Point data
393
394               CALL obs_int_h2d( kpk, kpk,      &
395                  &              zweig, zints(:,:,:,iobs), zobsk )
396
397            ENDIF
398
399
400            !-------------------------------------------------------------
401            ! Compute vertical second-derivative of the interpolating
402            ! polynomial at obs points
403            !-------------------------------------------------------------
404           
405            IF ( k1dint == 1 ) THEN
406               CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, &
407                  &                  pgdept, zobsmask )
408            ENDIF
409           
410            !----------------------------------------------------------------
411            !  Vertical interpolation to the observation point
412            !----------------------------------------------------------------
413            ista = prodatqc%npvsta(jobs,2)
414            iend = prodatqc%npvend(jobs,2)
415            CALL obs_int_z1d( kpk, &
416               & prodatqc%var(2)%mvk(ista:iend),&
417               & k1dint, iend - ista + 1, &
418               & prodatqc%var(2)%vdep(ista:iend),&
419               & zobsk, zobs2k, &
420               & prodatqc%var(2)%vmod(ista:iend),&
421               & pgdept, zobsmask )
422
423         ENDIF
424
425      END DO
426 
427      ! Deallocate the data for interpolation
428      DEALLOCATE( &
429         & igrdi, &
430         & igrdj, &
431         & zglam, &
432         & zgphi, &
433         & zmask, &
434         & zintt, &
435         & zints  &
436         & )
437      ! At the end of the day also get interpolated means
438      IF ( idayend == 0 ) THEN
439         DEALLOCATE( &
440            & zinmt,  &
441            & zinms   &
442            & )
443      ENDIF
444
445      prodatqc%nprofup = prodatqc%nprofup + ipro 
446     
447   END SUBROUTINE obs_pro_opt
448
449   SUBROUTINE obs_sla_opt( sladatqc, kt, kpi, kpj, kit000, &
450      &                    psshn, psshmask, k2dint )
451      !!-----------------------------------------------------------------------
452      !!
453      !!                     ***  ROUTINE obs_sla_opt  ***
454      !!
455      !! ** Purpose : Compute the model counterpart of sea level anomaly
456      !!              data by interpolating from the model grid to the
457      !!              observation point.
458      !!
459      !! ** Method  : Linearly interpolate to each observation point using
460      !!              the model values at the corners of the surrounding grid box.
461      !!
462      !!    The now model SSH is first computed at the obs (lon, lat) point.
463      !!
464      !!    Several horizontal interpolation schemes are available:
465      !!        - distance-weighted (great circle) (k2dint = 0)
466      !!        - distance-weighted (small angle)  (k2dint = 1)
467      !!        - bilinear (geographical grid)     (k2dint = 2)
468      !!        - bilinear (quadrilateral grid)    (k2dint = 3)
469      !!        - polynomial (quadrilateral grid)  (k2dint = 4)
470      !! 
471      !!    The sea level anomaly at the observation points is then computed
472      !!    by removing a mean dynamic topography (defined at the obs. point).
473      !!
474      !! ** Action  :
475      !!
476      !! History :
477      !!      ! 07-03 (A. Weaver)
478      !!-----------------------------------------------------------------------
479 
480      !! * Modules used
481      USE obs_surf_def  ! Definition of storage space for surface observations
482
483      IMPLICIT NONE
484
485      !! * Arguments
486      TYPE(obs_surf), INTENT(INOUT) :: sladatqc     ! Subset of surface data not failing screening
487      INTEGER, INTENT(IN) :: kt      ! Time step
488      INTEGER, INTENT(IN) :: kpi     ! Model grid parameters
489      INTEGER, INTENT(IN) :: kpj
490      INTEGER, INTENT(IN) :: kit000   ! Number of the first time step
491                                      !   (kit000-1 = restart time)
492      INTEGER, INTENT(IN) :: k2dint   ! Horizontal interpolation type (see header)
493      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
494         & psshn,  &    ! Model SSH field
495         & psshmask     ! Land-sea mask
496         
497      !! * Local declarations
498      INTEGER :: ji
499      INTEGER :: jj
500      INTEGER :: jobs
501      INTEGER :: inrc
502      INTEGER :: isla
503      INTEGER :: iobs
504      REAL(KIND=wp) :: zlam
505      REAL(KIND=wp) :: zphi
506      REAL(KIND=wp) :: zext(1), zobsmask(1)
507      REAL(kind=wp), DIMENSION(2,2,1) :: &
508         & zweig
509      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
510         & zmask, &
511         & zsshl, &
512         & zglam, &
513         & zgphi
514      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
515         & igrdi, &
516         & igrdj
517
518      !------------------------------------------------------------------------
519      ! Local initialization
520      !------------------------------------------------------------------------
521      ! ... Record and data counters
522      inrc = kt - kit000 + 2
523      isla = sladatqc%nsstp(inrc)
524
525      ! Get the data for interpolation
526
527      ALLOCATE( &
528         & igrdi(2,2,isla), &
529         & igrdj(2,2,isla), &
530         & zglam(2,2,isla), &
531         & zgphi(2,2,isla), &
532         & zmask(2,2,isla), &
533         & zsshl(2,2,isla)  &
534         & )
535     
536      DO jobs = sladatqc%nsurfup + 1, sladatqc%nsurfup + isla
537         iobs = jobs - sladatqc%nsurfup
538         igrdi(1,1,iobs) = sladatqc%mi(jobs)-1
539         igrdj(1,1,iobs) = sladatqc%mj(jobs)-1
540         igrdi(1,2,iobs) = sladatqc%mi(jobs)-1
541         igrdj(1,2,iobs) = sladatqc%mj(jobs)
542         igrdi(2,1,iobs) = sladatqc%mi(jobs)
543         igrdj(2,1,iobs) = sladatqc%mj(jobs)-1
544         igrdi(2,2,iobs) = sladatqc%mi(jobs)
545         igrdj(2,2,iobs) = sladatqc%mj(jobs)
546      END DO
547
548      CALL obs_int_comm_2d( 2, 2, isla, &
549         &                  igrdi, igrdj, glamt, zglam )
550      CALL obs_int_comm_2d( 2, 2, isla, &
551         &                  igrdi, igrdj, gphit, zgphi )
552      CALL obs_int_comm_2d( 2, 2, isla, &
553         &                  igrdi, igrdj, psshmask, zmask )
554      CALL obs_int_comm_2d( 2, 2, isla, &
555         &                  igrdi, igrdj, psshn, zsshl )
556
557      ! Loop over observations
558
559      DO jobs = sladatqc%nsurfup + 1, sladatqc%nsurfup + isla
560
561         iobs = jobs - sladatqc%nsurfup
562
563         IF ( kt /= sladatqc%mstp(jobs) ) THEN
564           
565            IF(lwp) THEN
566               WRITE(numout,*)
567               WRITE(numout,*) ' E R R O R : Observation',              &
568                  &            ' time step is not consistent with the', &
569                  &            ' model time step'
570               WRITE(numout,*) ' ========='
571               WRITE(numout,*)
572               WRITE(numout,*) ' Record  = ', jobs,                &
573                  &            ' kt      = ', kt,                  &
574                  &            ' mstp    = ', sladatqc%mstp(jobs), &
575                  &            ' ntyp    = ', sladatqc%ntyp(jobs)
576            ENDIF
577            CALL ctl_stop( 'obs_sla_opt', 'Inconsistent time' )
578           
579         ENDIF
580         
581         zlam = sladatqc%rlam(jobs)
582         zphi = sladatqc%rphi(jobs)
583
584         ! Get weights to interpolate the model SSH to the observation point
585         CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         &
586            &                   zglam(:,:,iobs), zgphi(:,:,iobs), &
587            &                   zmask(:,:,iobs), zweig, zobsmask )
588         
589
590         ! Interpolate the model SSH to the observation point
591         CALL obs_int_h2d( 1, 1,      &
592            &              zweig, zsshl(:,:,iobs),  zext )
593         
594         sladatqc%rext(jobs,1) = zext(1)
595         ! ... Remove the MDT at the observation point
596         sladatqc%rmod(jobs,1) = sladatqc%rext(jobs,1) - sladatqc%rext(jobs,2)
597
598      END DO
599
600      ! Deallocate the data for interpolation
601      DEALLOCATE( &
602         & igrdi, &
603         & igrdj, &
604         & zglam, &
605         & zgphi, &
606         & zmask, &
607         & zsshl  &
608         & )
609
610      sladatqc%nsurfup = sladatqc%nsurfup + isla
611
612   END SUBROUTINE obs_sla_opt
613
614   SUBROUTINE obs_sst_opt( sstdatqc, kt, kpi, kpj, kit000, &
615      &                    psstn, psstmask, k2dint )
616
617      !!-----------------------------------------------------------------------
618      !!
619      !!                     ***  ROUTINE obs_sst_opt  ***
620      !!
621      !! ** Purpose : Compute the model counterpart of surface temperature
622      !!              data by interpolating from the model grid to the
623      !!              observation point.
624      !!
625      !! ** Method  : Linearly interpolate to each observation point using
626      !!              the model values at the corners of the surrounding grid box.
627      !!
628      !!    The now model SST is first computed at the obs (lon, lat) point.
629      !!
630      !!    Several horizontal interpolation schemes are available:
631      !!        - distance-weighted (great circle) (k2dint = 0)
632      !!        - distance-weighted (small angle)  (k2dint = 1)
633      !!        - bilinear (geographical grid)     (k2dint = 2)
634      !!        - bilinear (quadrilateral grid)    (k2dint = 3)
635      !!        - polynomial (quadrilateral grid)  (k2dint = 4)
636      !!
637      !!
638      !! ** Action  :
639      !!
640      !! History :
641      !!        !  07-07  (S. Ricci ) : Original
642      !!     
643      !!-----------------------------------------------------------------------
644
645      !! * Modules used
646      USE obs_surf_def  ! Definition of storage space for surface observations
647
648      IMPLICIT NONE
649
650      !! * Arguments
651      TYPE(obs_surf), INTENT(INOUT) :: &
652         & sstdatqc     ! Subset of surface data not failing screening
653      INTEGER, INTENT(IN) :: kt        ! Time step
654      INTEGER, INTENT(IN) :: kpi       ! Model grid parameters
655      INTEGER, INTENT(IN) :: kpj
656      INTEGER, INTENT(IN) :: kit000    ! Number of the first time step
657                                       !   (kit000-1 = restart time)
658      INTEGER, INTENT(IN) :: k2dint    ! Horizontal interpolation type (see header)
659      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
660         & psstn,  &    ! Model SST field
661         & psstmask     ! Land-sea mask
662         
663      !! * Local declarations
664      INTEGER :: ji
665      INTEGER :: jj
666      INTEGER :: jobs
667      INTEGER :: inrc
668      INTEGER :: isst
669      INTEGER :: iobs
670      REAL(KIND=wp) :: zlam
671      REAL(KIND=wp) :: zphi
672      REAL(KIND=wp) :: zext(1), zobsmask(1)
673      REAL(kind=wp), DIMENSION(2,2,1) :: &
674         & zweig
675      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
676         & zmask, &
677         & zsstl, &
678         & zglam, &
679         & zgphi
680      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
681         & igrdi, &
682         & igrdj
683
684      !-----------------------------------------------------------------------
685      ! Local initialization
686      !-----------------------------------------------------------------------
687      ! ... Record and data counters
688      inrc = kt - kit000 + 2
689      isst = sstdatqc%nsstp(inrc)
690
691      ! Get the data for interpolation
692     
693      ALLOCATE( &
694         & igrdi(2,2,isst), &
695         & igrdj(2,2,isst), &
696         & zglam(2,2,isst), &
697         & zgphi(2,2,isst), &
698         & zmask(2,2,isst), &
699         & zsstl(2,2,isst)  &
700         & )
701     
702      DO jobs = sstdatqc%nsurfup + 1, sstdatqc%nsurfup + isst
703         iobs = jobs - sstdatqc%nsurfup
704         igrdi(1,1,iobs) = sstdatqc%mi(jobs)-1
705         igrdj(1,1,iobs) = sstdatqc%mj(jobs)-1
706         igrdi(1,2,iobs) = sstdatqc%mi(jobs)-1
707         igrdj(1,2,iobs) = sstdatqc%mj(jobs)
708         igrdi(2,1,iobs) = sstdatqc%mi(jobs)
709         igrdj(2,1,iobs) = sstdatqc%mj(jobs)-1
710         igrdi(2,2,iobs) = sstdatqc%mi(jobs)
711         igrdj(2,2,iobs) = sstdatqc%mj(jobs)
712      END DO
713     
714      CALL obs_int_comm_2d( 2, 2, isst, &
715         &                  igrdi, igrdj, glamt, zglam )
716      CALL obs_int_comm_2d( 2, 2, isst, &
717         &                  igrdi, igrdj, gphit, zgphi )
718      CALL obs_int_comm_2d( 2, 2, isst, &
719         &                  igrdi, igrdj, psstmask, zmask )
720      CALL obs_int_comm_2d( 2, 2, isst, &
721         &                  igrdi, igrdj, psstn, zsstl )
722     
723      ! Loop over observations
724
725      DO jobs = sstdatqc%nsurfup + 1, sstdatqc%nsurfup + isst
726         
727         iobs = jobs - sstdatqc%nsurfup
728         
729         IF ( kt /= sstdatqc%mstp(jobs) ) THEN
730           
731            IF(lwp) THEN
732               WRITE(numout,*)
733               WRITE(numout,*) ' E R R O R : Observation',              &
734                  &            ' time step is not consistent with the', &
735                  &            ' model time step'
736               WRITE(numout,*) ' ========='
737               WRITE(numout,*)
738               WRITE(numout,*) ' Record  = ', jobs,                &
739                  &            ' kt      = ', kt,                  &
740                  &            ' mstp    = ', sstdatqc%mstp(jobs), &
741                  &            ' ntyp    = ', sstdatqc%ntyp(jobs)
742            ENDIF
743            CALL ctl_stop( 'obs_sst_opt', 'Inconsistent time' )
744           
745         ENDIF
746         
747         zlam = sstdatqc%rlam(jobs)
748         zphi = sstdatqc%rphi(jobs)
749         
750         ! Get weights to interpolate the model SST to the observation point
751         CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         &
752            &                   zglam(:,:,iobs), zgphi(:,:,iobs), &
753            &                   zmask(:,:,iobs), zweig, zobsmask )
754           
755         ! Interpolate the model SST to the observation point
756         CALL obs_int_h2d( 1, 1,      &
757            &              zweig, zsstl(:,:,iobs),  zext )
758         
759         sstdatqc%rmod(jobs,1) = zext(1)
760         
761      END DO
762     
763      ! Deallocate the data for interpolation
764      DEALLOCATE( &
765         & igrdi, &
766         & igrdj, &
767         & zglam, &
768         & zgphi, &
769         & zmask, &
770         & zsstl  &
771         & )
772     
773      sstdatqc%nsurfup = sstdatqc%nsurfup + isst
774
775   END SUBROUTINE obs_sst_opt
776
777   SUBROUTINE obs_sss_opt
778      !!-----------------------------------------------------------------------
779      !!
780      !!                     ***  ROUTINE obs_sss_opt  ***
781      !!
782      !! ** Purpose : Compute the model counterpart of sea surface salinity
783      !!              data by interpolating from the model grid to the
784      !!              observation point.
785      !!
786      !! ** Method  :
787      !!
788      !! ** Action  :
789      !!
790      !! History :
791      !!      ! ??-??
792      !!-----------------------------------------------------------------------
793
794      IMPLICIT NONE
795
796   END SUBROUTINE obs_sss_opt
797
798   SUBROUTINE obs_seaice_opt( seaicedatqc, kt, kpi, kpj, kit000, &
799      &                    pseaicen, pseaicemask, k2dint )
800
801      !!-----------------------------------------------------------------------
802      !!
803      !!                     ***  ROUTINE obs_seaice_opt  ***
804      !!
805      !! ** Purpose : Compute the model counterpart of surface temperature
806      !!              data by interpolating from the model grid to the
807      !!              observation point.
808      !!
809      !! ** Method  : Linearly interpolate to each observation point using
810      !!              the model values at the corners of the surrounding grid box.
811      !!
812      !!    The now model sea ice is first computed at the obs (lon, lat) point.
813      !!
814      !!    Several horizontal interpolation schemes are available:
815      !!        - distance-weighted (great circle) (k2dint = 0)
816      !!        - distance-weighted (small angle)  (k2dint = 1)
817      !!        - bilinear (geographical grid)     (k2dint = 2)
818      !!        - bilinear (quadrilateral grid)    (k2dint = 3)
819      !!        - polynomial (quadrilateral grid)  (k2dint = 4)
820      !!
821      !!
822      !! ** Action  :
823      !!
824      !! History :
825      !!        !  07-07  (S. Ricci ) : Original
826      !!     
827      !!-----------------------------------------------------------------------
828
829      !! * Modules used
830      USE obs_surf_def  ! Definition of storage space for surface observations
831
832      IMPLICIT NONE
833
834      !! * Arguments
835      TYPE(obs_surf), INTENT(INOUT) :: seaicedatqc     ! Subset of surface data not failing screening
836      INTEGER, INTENT(IN) :: kt       ! Time step
837      INTEGER, INTENT(IN) :: kpi      ! Model grid parameters
838      INTEGER, INTENT(IN) :: kpj
839      INTEGER, INTENT(IN) :: kit000   ! Number of the first time step
840                                      !   (kit000-1 = restart time)
841      INTEGER, INTENT(IN) :: k2dint   ! Horizontal interpolation type (see header)
842      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
843         & pseaicen,  &    ! Model sea ice field
844         & pseaicemask     ! Land-sea mask
845         
846      !! * Local declarations
847      INTEGER :: ji
848      INTEGER :: jj
849      INTEGER :: jobs
850      INTEGER :: inrc
851      INTEGER :: iseaice
852      INTEGER :: iobs
853       
854      REAL(KIND=wp) :: zlam
855      REAL(KIND=wp) :: zphi
856      REAL(KIND=wp) :: zext(1), zobsmask(1)
857      REAL(kind=wp), DIMENSION(2,2,1) :: &
858         & zweig
859      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
860         & zmask, &
861         & zseaicel, &
862         & zglam, &
863         & zgphi
864      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
865         & igrdi, &
866         & igrdj
867
868      !------------------------------------------------------------------------
869      ! Local initialization
870      !------------------------------------------------------------------------
871      ! ... Record and data counters
872      inrc = kt - kit000 + 2
873      iseaice = seaicedatqc%nsstp(inrc)
874
875      ! Get the data for interpolation
876     
877      ALLOCATE( &
878         & igrdi(2,2,iseaice), &
879         & igrdj(2,2,iseaice), &
880         & zglam(2,2,iseaice), &
881         & zgphi(2,2,iseaice), &
882         & zmask(2,2,iseaice), &
883         & zseaicel(2,2,iseaice)  &
884         & )
885     
886      DO jobs = seaicedatqc%nsurfup + 1, seaicedatqc%nsurfup + iseaice
887         iobs = jobs - seaicedatqc%nsurfup
888         igrdi(1,1,iobs) = seaicedatqc%mi(jobs)-1
889         igrdj(1,1,iobs) = seaicedatqc%mj(jobs)-1
890         igrdi(1,2,iobs) = seaicedatqc%mi(jobs)-1
891         igrdj(1,2,iobs) = seaicedatqc%mj(jobs)
892         igrdi(2,1,iobs) = seaicedatqc%mi(jobs)
893         igrdj(2,1,iobs) = seaicedatqc%mj(jobs)-1
894         igrdi(2,2,iobs) = seaicedatqc%mi(jobs)
895         igrdj(2,2,iobs) = seaicedatqc%mj(jobs)
896      END DO
897     
898      CALL obs_int_comm_2d( 2, 2, iseaice, &
899         &                  igrdi, igrdj, glamt, zglam )
900      CALL obs_int_comm_2d( 2, 2, iseaice, &
901         &                  igrdi, igrdj, gphit, zgphi )
902      CALL obs_int_comm_2d( 2, 2, iseaice, &
903         &                  igrdi, igrdj, pseaicemask, zmask )
904      CALL obs_int_comm_2d( 2, 2, iseaice, &
905         &                  igrdi, igrdj, pseaicen, zseaicel )
906     
907      DO jobs = seaicedatqc%nsurfup + 1, seaicedatqc%nsurfup + iseaice
908         
909         iobs = jobs - seaicedatqc%nsurfup
910         
911         IF ( kt /= seaicedatqc%mstp(jobs) ) THEN
912           
913            IF(lwp) THEN
914               WRITE(numout,*)
915               WRITE(numout,*) ' E R R O R : Observation',              &
916                  &            ' time step is not consistent with the', &
917                  &            ' model time step'
918               WRITE(numout,*) ' ========='
919               WRITE(numout,*)
920               WRITE(numout,*) ' Record  = ', jobs,                &
921                  &            ' kt      = ', kt,                  &
922                  &            ' mstp    = ', seaicedatqc%mstp(jobs), &
923                  &            ' ntyp    = ', seaicedatqc%ntyp(jobs)
924            ENDIF
925            CALL ctl_stop( 'obs_seaice_opt', 'Inconsistent time' )
926           
927         ENDIF
928         
929         zlam = seaicedatqc%rlam(jobs)
930         zphi = seaicedatqc%rphi(jobs)
931         
932         ! Get weights to interpolate the model sea ice to the observation point
933         CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         &
934            &                   zglam(:,:,iobs), zgphi(:,:,iobs), &
935            &                   zmask(:,:,iobs), zweig, zobsmask )
936         
937         ! ... Interpolate the model sea ice to the observation point
938         CALL obs_int_h2d( 1, 1,      &
939            &              zweig, zseaicel(:,:,iobs),  zext )
940         
941         seaicedatqc%rmod(jobs,1) = zext(1)
942         
943      END DO
944     
945      ! Deallocate the data for interpolation
946      DEALLOCATE( &
947         & igrdi,    &
948         & igrdj,    &
949         & zglam,    &
950         & zgphi,    &
951         & zmask,    &
952         & zseaicel  &
953         & )
954     
955      seaicedatqc%nsurfup = seaicedatqc%nsurfup + iseaice
956
957   END SUBROUTINE obs_seaice_opt
958
959   SUBROUTINE obs_vel_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, &
960      &                    pun, pvn, pgdept, pumask, pvmask, k1dint, k2dint, &
961      &                    ld_dailyav )
962      !!-----------------------------------------------------------------------
963      !!
964      !!                     ***  ROUTINE obs_vel_opt  ***
965      !!
966      !! ** Purpose : Compute the model counterpart of velocity profile
967      !!              data by interpolating from the model grid to the
968      !!              observation point.
969      !!
970      !! ** Method  : Linearly interpolate zonal and meridional components of velocity
971      !!              to each observation point using the model values at the corners of
972      !!              the surrounding grid box. The model velocity components are on a
973      !!              staggered C- grid.
974      !!
975      !!    For velocity data from the TAO array, the model equivalent is
976      !!    a daily mean velocity field. So, we first compute
977      !!    the mean, then interpolate only at the end of the day.
978      !!
979      !! ** Action  :
980      !!
981      !! History :
982      !!    ! 07-03 (K. Mogensen)      : Temperature and Salinity profiles
983      !!    ! 08-10 (Maria Valdivieso) : Velocity component (U,V) profiles
984      !!-----------------------------------------------------------------------
985   
986      !! * Modules used
987      USE obs_profiles_def ! Definition of storage space for profile obs.
988
989      IMPLICIT NONE
990
991      !! * Arguments
992      TYPE(obs_prof), INTENT(INOUT) :: &
993         & prodatqc        ! Subset of profile data not failing screening
994      INTEGER, INTENT(IN) :: kt        ! Time step
995      INTEGER, INTENT(IN) :: kpi       ! Model grid parameters
996      INTEGER, INTENT(IN) :: kpj
997      INTEGER, INTENT(IN) :: kpk 
998      INTEGER, INTENT(IN) :: kit000    ! Number of the first time step
999                                       !   (kit000-1 = restart time)
1000      INTEGER, INTENT(IN) :: k1dint    ! Vertical interpolation type (see header)
1001      INTEGER, INTENT(IN) :: k2dint    ! Horizontal interpolation type (see header)
1002      INTEGER, INTENT(IN) :: kdaystp   ! Number of time steps per day                   
1003      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &
1004         & pun,    &    ! Model zonal component of velocity
1005         & pvn,    &    ! Model meridional component of velocity
1006         & pumask, &    ! Land-sea mask
1007         & pvmask       ! Land-sea mask
1008      REAL(KIND=wp), INTENT(IN), DIMENSION(kpk) :: &
1009         & pgdept       ! Model array of depth levels
1010      LOGICAL, INTENT(IN) :: ld_dailyav
1011         
1012      !! * Local declarations
1013      INTEGER :: ji
1014      INTEGER :: jj
1015      INTEGER :: jk
1016      INTEGER :: jobs
1017      INTEGER :: inrc
1018      INTEGER :: ipro
1019      INTEGER :: idayend
1020      INTEGER :: ista
1021      INTEGER :: iend
1022      INTEGER :: iobs
1023      INTEGER, DIMENSION(imaxavtypes) :: &
1024         & idailyavtypes
1025      REAL(KIND=wp) :: zlam
1026      REAL(KIND=wp) :: zphi
1027      REAL(KIND=wp) :: zdaystp
1028      REAL(KIND=wp), DIMENSION(kpk) :: &
1029         & zobsmasku, &
1030         & zobsmaskv, &
1031         & zobsmask,  &
1032         & zobsk,     &
1033         & zobs2k
1034      REAL(KIND=wp), DIMENSION(2,2,kpk) :: &
1035         & zweigu,zweigv
1036      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
1037         & zumask, zvmask, &
1038         & zintu, &
1039         & zintv, &
1040         & zinmu, &
1041         & zinmv
1042      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
1043         & zglamu, zglamv, &
1044         & zgphiu, zgphiv
1045      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
1046         & igrdiu, &
1047         & igrdju, &
1048         & igrdiv, &
1049         & igrdjv
1050
1051      !------------------------------------------------------------------------
1052      ! Local initialization
1053      !------------------------------------------------------------------------
1054      ! ... Record and data counters
1055      inrc = kt - kit000 + 2
1056      ipro = prodatqc%npstp(inrc)
1057
1058      ! Initialize daily mean for first timestep
1059      idayend = MOD( kt - kit000 + 1, kdaystp )
1060
1061      ! Added kt == 0 test to catch restart case
1062      IF ( idayend == 1 .OR. kt == 0) THEN
1063         IF (lwp) WRITE(numout,*) 'Reset prodatqc%vdmean on time-step: ',kt
1064         prodatqc%vdmean(:,:,:,1) = 0.0
1065         prodatqc%vdmean(:,:,:,2) = 0.0
1066      ENDIF
1067
1068      ! Increment the zonal velocity field for computing daily mean
1069      prodatqc%vdmean(:,:,:,1) = prodatqc%vdmean(:,:,:,1) + pun(:,:,:)
1070      ! Increment the meridional velocity field for computing daily mean
1071      prodatqc%vdmean(:,:,:,2) = prodatqc%vdmean(:,:,:,2) + pvn(:,:,:)
1072   
1073      ! Compute the daily mean at the end of day
1074      zdaystp = 1.0 / REAL( kdaystp )
1075      IF ( idayend == 0 ) THEN
1076         prodatqc%vdmean(:,:,:,1) = prodatqc%vdmean(:,:,:,1) * zdaystp
1077         prodatqc%vdmean(:,:,:,2) = prodatqc%vdmean(:,:,:,2) * zdaystp
1078      ENDIF
1079
1080      ! Get the data for interpolation
1081      ALLOCATE( &
1082         & igrdiu(2,2,ipro),      &
1083         & igrdju(2,2,ipro),      &
1084         & igrdiv(2,2,ipro),      &
1085         & igrdjv(2,2,ipro),      &
1086         & zglamu(2,2,ipro), zglamv(2,2,ipro), &
1087         & zgphiu(2,2,ipro), zgphiv(2,2,ipro), &
1088         & zumask(2,2,kpk,ipro), zvmask(2,2,kpk,ipro), &
1089         & zintu(2,2,kpk,ipro),  &
1090         & zintv(2,2,kpk,ipro)   &
1091         & )
1092
1093      DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
1094         iobs = jobs - prodatqc%nprofup
1095         igrdiu(1,1,iobs) = prodatqc%mi(jobs,1)-1
1096         igrdju(1,1,iobs) = prodatqc%mj(jobs,1)-1
1097         igrdiu(1,2,iobs) = prodatqc%mi(jobs,1)-1
1098         igrdju(1,2,iobs) = prodatqc%mj(jobs,1)
1099         igrdiu(2,1,iobs) = prodatqc%mi(jobs,1)
1100         igrdju(2,1,iobs) = prodatqc%mj(jobs,1)-1
1101         igrdiu(2,2,iobs) = prodatqc%mi(jobs,1)
1102         igrdju(2,2,iobs) = prodatqc%mj(jobs,1)
1103         igrdiv(1,1,iobs) = prodatqc%mi(jobs,2)-1
1104         igrdjv(1,1,iobs) = prodatqc%mj(jobs,2)-1
1105         igrdiv(1,2,iobs) = prodatqc%mi(jobs,2)-1
1106         igrdjv(1,2,iobs) = prodatqc%mj(jobs,2)
1107         igrdiv(2,1,iobs) = prodatqc%mi(jobs,2)
1108         igrdjv(2,1,iobs) = prodatqc%mj(jobs,2)-1
1109         igrdiv(2,2,iobs) = prodatqc%mi(jobs,2)
1110         igrdjv(2,2,iobs) = prodatqc%mj(jobs,2)
1111      END DO
1112
1113      CALL obs_int_comm_2d( 2, 2, ipro, igrdiu, igrdju, glamu, zglamu )
1114      CALL obs_int_comm_2d( 2, 2, ipro, igrdiu, igrdju, gphiu, zgphiu )
1115      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, pumask, zumask )
1116      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, pun, zintu )
1117
1118      CALL obs_int_comm_2d( 2, 2, ipro, igrdiv, igrdjv, glamv, zglamv )
1119      CALL obs_int_comm_2d( 2, 2, ipro, igrdiv, igrdjv, gphiv, zgphiv )
1120      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, pvmask, zvmask )
1121      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, pvn, zintv )
1122
1123      ! At the end of the day also get interpolated means
1124      IF ( idayend == 0 ) THEN
1125
1126         ALLOCATE( &
1127            & zinmu(2,2,kpk,ipro),  &
1128            & zinmv(2,2,kpk,ipro)   &
1129            & )
1130
1131         CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, &
1132            &                  prodatqc%vdmean(:,:,:,1), zinmu )
1133         CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, &
1134            &                  prodatqc%vdmean(:,:,:,2), zinmv )
1135
1136      ENDIF
1137
1138! loop over observations
1139
1140      DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
1141
1142         iobs = jobs - prodatqc%nprofup
1143
1144         IF ( kt /= prodatqc%mstp(jobs) ) THEN
1145           
1146            IF(lwp) THEN
1147               WRITE(numout,*)
1148               WRITE(numout,*) ' E R R O R : Observation',              &
1149                  &            ' time step is not consistent with the', &
1150                  &            ' model time step'
1151               WRITE(numout,*) ' ========='
1152               WRITE(numout,*)
1153               WRITE(numout,*) ' Record  = ', jobs,                    &
1154                  &            ' kt      = ', kt,                      &
1155                  &            ' mstp    = ', prodatqc%mstp(jobs), &
1156                  &            ' ntyp    = ', prodatqc%ntyp(jobs)
1157            ENDIF
1158            CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' )
1159         ENDIF
1160         
1161         zlam = prodatqc%rlam(jobs)
1162         zphi = prodatqc%rphi(jobs)
1163
1164         ! Initialize observation masks
1165
1166         zobsmasku(:) = 0.0
1167         zobsmaskv(:) = 0.0
1168         
1169         ! Horizontal weights and vertical mask
1170
1171         IF  ( prodatqc%npvend(jobs,1) > 0 ) THEN
1172
1173            CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi,     &
1174               &                   zglamu(:,:,iobs), zgphiu(:,:,iobs), &
1175               &                   zumask(:,:,:,iobs), zweigu, zobsmasku )
1176
1177         ENDIF
1178
1179         
1180         IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
1181
1182            CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi,     &
1183               &                   zglamv(:,:,iobs), zgphiv(:,:,iobs), &
1184               &                   zvmask(:,:,:,iobs), zweigv, zobsmasku )
1185
1186         ENDIF
1187
1188         ! Ensure that the vertical mask on u and v are consistent.
1189
1190         zobsmask(:) = MIN( zobsmasku(:), zobsmaskv(:) )
1191
1192         IF ( prodatqc%npvend(jobs,1) > 0 ) THEN
1193
1194            zobsk(:) = obfillflt
1195
1196       IF ( ld_dailyav ) THEN
1197
1198               IF ( idayend == 0 )  THEN
1199                 
1200                  ! Daily averaged data
1201                 
1202                  CALL obs_int_h2d( kpk, kpk,      &
1203                     &              zweigu, zinmu(:,:,:,iobs), zobsk )
1204                 
1205                 
1206               ELSE
1207               
1208                  CALL ctl_stop( ' A nonzero' //     &
1209                     &           ' number of U profile data should' // &
1210                     &           ' only occur at the end of a given day' )
1211
1212               ENDIF
1213         
1214            ELSE 
1215               
1216               ! Point data
1217
1218               CALL obs_int_h2d( kpk, kpk,      &
1219                  &              zweigu, zintu(:,:,:,iobs), zobsk )
1220
1221            ENDIF
1222
1223            !-------------------------------------------------------------
1224            ! Compute vertical second-derivative of the interpolating
1225            ! polynomial at obs points
1226            !-------------------------------------------------------------
1227           
1228            IF ( k1dint == 1 ) THEN
1229               CALL obs_int_z1d_spl( kpk, zobsk, zobs2k,   &
1230                  &                  pgdept, zobsmask )
1231            ENDIF
1232           
1233            !-----------------------------------------------------------------
1234            !  Vertical interpolation to the observation point
1235            !-----------------------------------------------------------------
1236            ista = prodatqc%npvsta(jobs,1)
1237            iend = prodatqc%npvend(jobs,1)
1238            CALL obs_int_z1d( kpk,                &
1239               & prodatqc%var(1)%mvk(ista:iend),  &
1240               & k1dint, iend - ista + 1,         &
1241               & prodatqc%var(1)%vdep(ista:iend), &
1242               & zobsk, zobs2k,                   &
1243               & prodatqc%var(1)%vmod(ista:iend), &
1244               & pgdept, zobsmask )
1245
1246         ENDIF
1247
1248         IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
1249
1250            zobsk(:) = obfillflt
1251
1252            IF ( ld_dailyav ) THEN
1253
1254               IF ( idayend == 0 )  THEN
1255
1256                  ! Daily averaged data
1257                 
1258                  CALL obs_int_h2d( kpk, kpk,      &
1259                     &              zweigv, zinmv(:,:,:,iobs), zobsk )
1260                 
1261               ELSE
1262
1263                  CALL ctl_stop( ' A nonzero' //     &
1264                     &           ' number of V profile data should' // &
1265                     &           ' only occur at the end of a given day' )
1266
1267               ENDIF
1268
1269            ELSE
1270               
1271               ! Point data
1272
1273               CALL obs_int_h2d( kpk, kpk,      &
1274                  &              zweigv, zintv(:,:,:,iobs), zobsk )
1275
1276            ENDIF
1277
1278
1279            !-------------------------------------------------------------
1280            ! Compute vertical second-derivative of the interpolating
1281            ! polynomial at obs points
1282            !-------------------------------------------------------------
1283           
1284            IF ( k1dint == 1 ) THEN
1285               CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, &
1286                  &                  pgdept, zobsmask )
1287            ENDIF
1288           
1289            !----------------------------------------------------------------
1290            !  Vertical interpolation to the observation point
1291            !----------------------------------------------------------------
1292            ista = prodatqc%npvsta(jobs,2)
1293            iend = prodatqc%npvend(jobs,2)
1294            CALL obs_int_z1d( kpk, &
1295               & prodatqc%var(2)%mvk(ista:iend),&
1296               & k1dint, iend - ista + 1, &
1297               & prodatqc%var(2)%vdep(ista:iend),&
1298               & zobsk, zobs2k, &
1299               & prodatqc%var(2)%vmod(ista:iend),&
1300               & pgdept, zobsmask )
1301
1302         ENDIF
1303
1304      END DO
1305 
1306      ! Deallocate the data for interpolation
1307      DEALLOCATE( &
1308         & igrdiu, &
1309         & igrdju, &
1310         & igrdiv, &
1311         & igrdjv, &
1312         & zglamu, zglamv, &
1313         & zgphiu, zgphiv, &
1314         & zumask, zvmask, &
1315         & zintu, &
1316         & zintv  &
1317         & )
1318      ! At the end of the day also get interpolated means
1319      IF ( idayend == 0 ) THEN
1320         DEALLOCATE( &
1321            & zinmu,  &
1322            & zinmv   &
1323            & )
1324      ENDIF
1325
1326      prodatqc%nprofup = prodatqc%nprofup + ipro 
1327     
1328   END SUBROUTINE obs_vel_opt
1329
1330END MODULE obs_oper
1331
Note: See TracBrowser for help on using the repository browser.