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 branches/UKMO/dev_r4650_general_vert_coord_obsoper_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

source: branches/UKMO/dev_r4650_general_vert_coord_obsoper_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90 @ 6855

Last change on this file since 6855 was 6855, checked in by dford, 8 years ago

Initial implementation of observation operator for SPM.

File size: 87.9 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_pro_sco_opt: Compute the model counterpart of temperature and
12   !!                    salinity observations from profiles in generalised
13   !!                    vertical coordinates
14   !!   obs_sla_opt :    Compute the model counterpart of sea level anomaly
15   !!                    observations
16   !!   obs_sst_opt :    Compute the model counterpart of sea surface temperature
17   !!                    observations
18   !!   obs_sss_opt :    Compute the model counterpart of sea surface salinity
19   !!                    observations
20   !!   obs_seaice_opt : Compute the model counterpart of sea ice concentration
21   !!                    observations
22   !!
23   !!   obs_vel_opt :    Compute the model counterpart of zonal and meridional
24   !!                    components of velocity from observations.
25   !!   obs_logchl_opt : Compute the model counterpart of log10(chlorophyll)
26   !!                    observations
27   !!   obs_spm_opt :    Compute the model counterpart of spm
28   !!                    observations
29   !!----------------------------------------------------------------------
30
31   !! * Modules used   
32   USE par_kind, ONLY : &         ! Precision variables
33      & wp
34   USE in_out_manager             ! I/O manager
35   USE obs_inter_sup              ! Interpolation support
36   USE obs_inter_h2d, ONLY : &    ! Horizontal interpolation to the observation pt
37      & obs_int_h2d, &
38      & obs_int_h2d_init
39   USE obs_inter_z1d, ONLY : &    ! Vertical interpolation to the observation pt
40      & obs_int_z1d,    &
41      & obs_int_z1d_spl
42   USE obs_const,  ONLY :     &
43      & obfillflt      ! Fillvalue   
44   USE dom_oce,       ONLY : &
45      & glamt, glamu, glamv, &
46      & gphit, gphiu, gphiv, & 
47#if defined key_vvl 
48      & gdept_n 
49#else 
50      & gdept_0 
51#endif 
52   USE lib_mpp,       ONLY : &
53      & ctl_warn, ctl_stop
54   USE obs_grid,      ONLY : & 
55      & obs_level_search     
56     
57   IMPLICIT NONE
58
59   !! * Routine accessibility
60   PRIVATE
61
62   PUBLIC obs_pro_opt, &  ! Compute the model counterpart of profile observations
63      &   obs_pro_sco_opt, &  ! Compute the model counterpart of profile observations
64                              ! in generalised vertical coordinates
65      &   obs_sla_opt, &  ! Compute the model counterpart of SLA observations
66      &   obs_sst_opt, &  ! Compute the model counterpart of SST observations
67      &   obs_sss_opt, &  ! Compute the model counterpart of SSS observations
68      &   obs_seaice_opt, &
69      &   obs_vel_opt, &  ! Compute the model counterpart of velocity profile data
70      &   obs_logchl_opt, & ! Compute the model counterpart of logchl data
71      &   obs_spm_opt     ! Compute the model counterpart of spm data
72
73   INTEGER, PARAMETER, PUBLIC :: imaxavtypes = 20 ! Max number of daily avgd obs types
74
75   !!----------------------------------------------------------------------
76   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
77   !! $Id$
78   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
79   !!----------------------------------------------------------------------
80
81   !! * Substitutions
82#  include "domzgr_substitute.h90" 
83CONTAINS
84
85   SUBROUTINE obs_pro_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, &
86      &                    ptn, psn, pgdept, ptmask, k1dint, k2dint, &
87      &                    kdailyavtypes )
88      !!-----------------------------------------------------------------------
89      !!
90      !!                     ***  ROUTINE obs_pro_opt  ***
91      !!
92      !! ** Purpose : Compute the model counterpart of profiles
93      !!              data by interpolating from the model grid to the
94      !!              observation point.
95      !!
96      !! ** Method  : Linearly interpolate to each observation point using
97      !!              the model values at the corners of the surrounding grid box.
98      !!
99      !!    First, a vertical profile of horizontally interpolated model
100      !!    now temperatures is computed at the obs (lon, lat) point.
101      !!    Several horizontal interpolation schemes are available:
102      !!        - distance-weighted (great circle) (k2dint = 0)
103      !!        - distance-weighted (small angle)  (k2dint = 1)
104      !!        - bilinear (geographical grid)     (k2dint = 2)
105      !!        - bilinear (quadrilateral grid)    (k2dint = 3)
106      !!        - polynomial (quadrilateral grid)  (k2dint = 4)
107      !!
108      !!    Next, the vertical temperature profile is interpolated to the
109      !!    data depth points. Two vertical interpolation schemes are
110      !!    available:
111      !!        - linear       (k1dint = 0)
112      !!        - Cubic spline (k1dint = 1)
113      !!
114      !!    For the cubic spline the 2nd derivative of the interpolating
115      !!    polynomial is computed before entering the vertical interpolation
116      !!    routine.
117      !!
118      !!    For ENACT moored buoy data (e.g., TAO), the model equivalent is
119      !!    a daily mean model temperature field. So, we first compute
120      !!    the mean, then interpolate only at the end of the day.
121      !!
122      !!    Note: the in situ temperature observations must be converted
123      !!    to potential temperature (the model variable) prior to
124      !!    assimilation.
125      !!??????????????????????????????????????????????????????????????
126      !!    INCLUDE POTENTIAL TEMP -> IN SITU TEMP IN OBS OPERATOR???
127      !!??????????????????????????????????????????????????????????????
128      !!
129      !! ** Action  :
130      !!
131      !! History :
132      !!      ! 97-11 (A. Weaver, S. Ricci, N. Daget)
133      !!      ! 06-03 (G. Smith) NEMOVAR migration
134      !!      ! 06-10 (A. Weaver) Cleanup
135      !!      ! 07-01 (K. Mogensen) Merge of temperature and salinity
136      !!      ! 07-03 (K. Mogensen) General handling of profiles
137      !!-----------------------------------------------------------------------
138 
139      !! * Modules used
140      USE obs_profiles_def ! Definition of storage space for profile obs.
141
142      IMPLICIT NONE
143
144      !! * Arguments
145      TYPE(obs_prof), INTENT(INOUT) :: prodatqc  ! Subset of profile data not failing screening
146      INTEGER, INTENT(IN) :: kt        ! Time step
147      INTEGER, INTENT(IN) :: kpi       ! Model grid parameters
148      INTEGER, INTENT(IN) :: kpj
149      INTEGER, INTENT(IN) :: kpk
150      INTEGER, INTENT(IN) :: kit000    ! Number of the first time step
151                                       !   (kit000-1 = restart time)
152      INTEGER, INTENT(IN) :: k1dint    ! Vertical interpolation type (see header)
153      INTEGER, INTENT(IN) :: k2dint    ! Horizontal interpolation type (see header)
154      INTEGER, INTENT(IN) :: kdaystp   ! Number of time steps per day                   
155      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &
156         & ptn,    &    ! Model temperature field
157         & psn,    &    ! Model salinity field
158         & ptmask       ! Land-sea mask
159      REAL(KIND=wp), INTENT(IN), DIMENSION(kpk) :: &
160         & pgdept       ! Model array of depth levels
161      INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &
162         & kdailyavtypes! Types for daily averages
163      !! * Local declarations
164      INTEGER ::   ji
165      INTEGER ::   jj
166      INTEGER ::   jk
167      INTEGER ::   jobs
168      INTEGER ::   inrc
169      INTEGER ::   ipro
170      INTEGER ::   idayend
171      INTEGER ::   ista
172      INTEGER ::   iend
173      INTEGER ::   iobs
174      INTEGER, DIMENSION(imaxavtypes) :: &
175         & idailyavtypes
176      REAL(KIND=wp) :: zlam
177      REAL(KIND=wp) :: zphi
178      REAL(KIND=wp) :: zdaystp
179      REAL(KIND=wp), DIMENSION(kpk) :: &
180         & zobsmask, &
181         & zobsk,    &
182         & zobs2k
183      REAL(KIND=wp), DIMENSION(2,2,kpk) :: &
184         & zweig
185      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
186         & zmask, &
187         & zintt, &
188         & zints, &
189         & zinmt, &
190         & zinms
191      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
192         & zglam, &
193         & zgphi
194      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
195         & igrdi, &
196         & igrdj
197
198      !------------------------------------------------------------------------
199      ! Local initialization
200      !------------------------------------------------------------------------
201      ! ... Record and data counters
202      inrc = kt - kit000 + 2
203      ipro = prodatqc%npstp(inrc)
204 
205      ! Daily average types
206      IF ( PRESENT(kdailyavtypes) ) THEN
207         idailyavtypes(:) = kdailyavtypes(:)
208      ELSE
209         idailyavtypes(:) = -1
210      ENDIF
211
212      ! Initialize daily mean for first timestep
213      idayend = MOD( kt - kit000 + 1, kdaystp )
214
215      ! Added kt == 0 test to catch restart case
216      IF ( idayend == 1 .OR. kt == 0) THEN
217         IF (lwp) WRITE(numout,*) 'Reset prodatqc%vdmean on time-step: ',kt
218         DO jk = 1, jpk
219            DO jj = 1, jpj
220               DO ji = 1, jpi
221                  prodatqc%vdmean(ji,jj,jk,1) = 0.0
222                  prodatqc%vdmean(ji,jj,jk,2) = 0.0
223               END DO
224            END DO
225         END DO
226      ENDIF
227
228      DO jk = 1, jpk
229         DO jj = 1, jpj
230            DO ji = 1, jpi
231               ! Increment the temperature field for computing daily mean
232               prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &
233                  &                        + ptn(ji,jj,jk)
234               ! Increment the salinity field for computing daily mean
235               prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &
236                  &                        + psn(ji,jj,jk)
237            END DO
238         END DO
239      END DO
240   
241      ! Compute the daily mean at the end of day
242      zdaystp = 1.0 / REAL( kdaystp )
243      IF ( idayend == 0 ) THEN
244         DO jk = 1, jpk
245            DO jj = 1, jpj
246               DO ji = 1, jpi
247                  prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &
248                     &                        * zdaystp
249                  prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &
250                  &                           * zdaystp
251               END DO
252            END DO
253         END DO
254      ENDIF
255
256      ! Get the data for interpolation
257      ALLOCATE( &
258         & igrdi(2,2,ipro),      &
259         & igrdj(2,2,ipro),      &
260         & zglam(2,2,ipro),      &
261         & zgphi(2,2,ipro),      &
262         & zmask(2,2,kpk,ipro),  &
263         & zintt(2,2,kpk,ipro),  &
264         & zints(2,2,kpk,ipro)   &
265         & )
266
267      DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
268         iobs = jobs - prodatqc%nprofup
269         igrdi(1,1,iobs) = prodatqc%mi(jobs,1)-1
270         igrdj(1,1,iobs) = prodatqc%mj(jobs,1)-1
271         igrdi(1,2,iobs) = prodatqc%mi(jobs,1)-1
272         igrdj(1,2,iobs) = prodatqc%mj(jobs,1)
273         igrdi(2,1,iobs) = prodatqc%mi(jobs,1)
274         igrdj(2,1,iobs) = prodatqc%mj(jobs,1)-1
275         igrdi(2,2,iobs) = prodatqc%mi(jobs,1)
276         igrdj(2,2,iobs) = prodatqc%mj(jobs,1)
277      END DO
278
279      CALL obs_int_comm_2d( 2, 2, ipro, igrdi, igrdj, glamt, zglam )
280      CALL obs_int_comm_2d( 2, 2, ipro, igrdi, igrdj, gphit, zgphi )
281      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, ptmask,zmask )
282      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, ptn,   zintt )
283      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, psn,   zints )
284
285      ! At the end of the day also get interpolated means
286      IF ( idayend == 0 ) THEN
287
288         ALLOCATE( &
289            & zinmt(2,2,kpk,ipro),  &
290            & zinms(2,2,kpk,ipro)   &
291            & )
292
293         CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, &
294            &                  prodatqc%vdmean(:,:,:,1), zinmt )
295         CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, &
296            &                  prodatqc%vdmean(:,:,:,2), zinms )
297
298      ENDIF
299
300      DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
301
302         iobs = jobs - prodatqc%nprofup
303
304         IF ( kt /= prodatqc%mstp(jobs) ) THEN
305           
306            IF(lwp) THEN
307               WRITE(numout,*)
308               WRITE(numout,*) ' E R R O R : Observation',              &
309                  &            ' time step is not consistent with the', &
310                  &            ' model time step'
311               WRITE(numout,*) ' ========='
312               WRITE(numout,*)
313               WRITE(numout,*) ' Record  = ', jobs,                    &
314                  &            ' kt      = ', kt,                      &
315                  &            ' mstp    = ', prodatqc%mstp(jobs), &
316                  &            ' ntyp    = ', prodatqc%ntyp(jobs)
317            ENDIF
318            CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' )
319         ENDIF
320         
321         zlam = prodatqc%rlam(jobs)
322         zphi = prodatqc%rphi(jobs)
323         
324         ! Horizontal weights and vertical mask
325
326         IF ( ( prodatqc%npvend(jobs,1) > 0 ) .OR. &
327            & ( prodatqc%npvend(jobs,2) > 0 ) ) THEN
328
329            CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi,     &
330               &                   zglam(:,:,iobs), zgphi(:,:,iobs), &
331               &                   zmask(:,:,:,iobs), zweig, zobsmask )
332
333         ENDIF
334
335         IF ( prodatqc%npvend(jobs,1) > 0 ) THEN
336
337            zobsk(:) = obfillflt
338
339       IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN
340
341               IF ( idayend == 0 )  THEN
342                 
343                  ! Daily averaged moored buoy (MRB) data
344                 
345                  CALL obs_int_h2d( kpk, kpk,      &
346                     &              zweig, zinmt(:,:,:,iobs), zobsk )
347                 
348                 
349               ELSE
350               
351                  CALL ctl_stop( ' A nonzero' //     &
352                     &           ' number of profile T BUOY data should' // &
353                     &           ' only occur at the end of a given day' )
354
355               ENDIF
356         
357            ELSE 
358               
359               ! Point data
360
361               CALL obs_int_h2d( kpk, kpk,      &
362                  &              zweig, zintt(:,:,:,iobs), zobsk )
363
364            ENDIF
365
366            !-------------------------------------------------------------
367            ! Compute vertical second-derivative of the interpolating
368            ! polynomial at obs points
369            !-------------------------------------------------------------
370           
371            IF ( k1dint == 1 ) THEN
372               CALL obs_int_z1d_spl( kpk, zobsk, zobs2k,   &
373                  &                  pgdept, zobsmask )
374            ENDIF
375           
376            !-----------------------------------------------------------------
377            !  Vertical interpolation to the observation point
378            !-----------------------------------------------------------------
379            ista = prodatqc%npvsta(jobs,1)
380            iend = prodatqc%npvend(jobs,1)
381            CALL obs_int_z1d( kpk,                &
382               & prodatqc%var(1)%mvk(ista:iend),  &
383               & k1dint, iend - ista + 1,         &
384               & prodatqc%var(1)%vdep(ista:iend), &
385               & zobsk, zobs2k,                   &
386               & prodatqc%var(1)%vmod(ista:iend), &
387               & pgdept, zobsmask )
388
389         ENDIF
390
391         IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
392
393            zobsk(:) = obfillflt
394
395            IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN
396
397               IF ( idayend == 0 )  THEN
398
399                  ! Daily averaged moored buoy (MRB) data
400                 
401                  CALL obs_int_h2d( kpk, kpk,      &
402                     &              zweig, zinms(:,:,:,iobs), zobsk )
403                 
404               ELSE
405
406                  CALL ctl_stop( ' A nonzero' //     &
407                     &           ' number of profile S BUOY data should' // &
408                     &           ' only occur at the end of a given day' )
409
410               ENDIF
411
412            ELSE
413               
414               ! Point data
415
416               CALL obs_int_h2d( kpk, kpk,      &
417                  &              zweig, zints(:,:,:,iobs), zobsk )
418
419            ENDIF
420
421
422            !-------------------------------------------------------------
423            ! Compute vertical second-derivative of the interpolating
424            ! polynomial at obs points
425            !-------------------------------------------------------------
426           
427            IF ( k1dint == 1 ) THEN
428               CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, &
429                  &                  pgdept, zobsmask )
430            ENDIF
431           
432            !----------------------------------------------------------------
433            !  Vertical interpolation to the observation point
434            !----------------------------------------------------------------
435            ista = prodatqc%npvsta(jobs,2)
436            iend = prodatqc%npvend(jobs,2)
437            CALL obs_int_z1d( kpk, &
438               & prodatqc%var(2)%mvk(ista:iend),&
439               & k1dint, iend - ista + 1, &
440               & prodatqc%var(2)%vdep(ista:iend),&
441               & zobsk, zobs2k, &
442               & prodatqc%var(2)%vmod(ista:iend),&
443               & pgdept, zobsmask )
444
445         ENDIF
446
447      END DO
448 
449      ! Deallocate the data for interpolation
450      DEALLOCATE( &
451         & igrdi, &
452         & igrdj, &
453         & zglam, &
454         & zgphi, &
455         & zmask, &
456         & zintt, &
457         & zints  &
458         & )
459      ! At the end of the day also get interpolated means
460      IF ( idayend == 0 ) THEN
461         DEALLOCATE( &
462            & zinmt,  &
463            & zinms   &
464            & )
465      ENDIF
466
467      prodatqc%nprofup = prodatqc%nprofup + ipro 
468     
469   END SUBROUTINE obs_pro_opt
470
471   SUBROUTINE obs_pro_sco_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, & 
472      &                    ptn, psn, pgdept, pgdepw, ptmask, k1dint, k2dint, & 
473      &                    kdailyavtypes ) 
474      !!-----------------------------------------------------------------------
475      !!
476      !!                     ***  ROUTINE obs_pro_opt  ***
477      !!
478      !! ** Purpose : Compute the model counterpart of profiles
479      !!              data by interpolating from the model grid to the 
480      !!              observation point. Generalised vertical coordinate version
481      !!
482      !! ** Method  : Linearly interpolate to each observation point using 
483      !!              the model values at the corners of the surrounding grid box.
484      !!
485      !!          First, model values on the model grid are interpolated vertically to the
486      !!          Depths of the profile observations.  Two vertical interpolation schemes are
487      !!          available:
488      !!          - linear       (k1dint = 0)
489      !!          - Cubic spline (k1dint = 1)   
490      !!
491      !!
492      !!         Secondly the interpolated values are interpolated horizontally to the 
493      !!         obs (lon, lat) point.
494      !!         Several horizontal interpolation schemes are available:
495      !!        - distance-weighted (great circle) (k2dint = 0)
496      !!        - distance-weighted (small angle)  (k2dint = 1)
497      !!        - bilinear (geographical grid)     (k2dint = 2)
498      !!        - bilinear (quadrilateral grid)    (k2dint = 3)
499      !!        - polynomial (quadrilateral grid)  (k2dint = 4)
500      !!
501      !!    For the cubic spline the 2nd derivative of the interpolating 
502      !!    polynomial is computed before entering the vertical interpolation 
503      !!    routine.
504      !!
505      !!    For ENACT moored buoy data (e.g., TAO), the model equivalent is
506      !!    a daily mean model temperature field. So, we first compute
507      !!    the mean, then interpolate only at the end of the day.
508      !!
509      !!    This is the procedure to be used with generalised vertical model 
510      !!    coordinates (ie s-coordinates. It is ~4x slower than the equivalent
511      !!    horizontal then vertical interpolation algorithm, but can deal with situations
512      !!    where the model levels are not flat.
513      !!    ONLY PERFORMED if ln_sco=.TRUE. 
514      !!       
515      !!    Note: the in situ temperature observations must be converted
516      !!    to potential temperature (the model variable) prior to
517      !!    assimilation. 
518      !!??????????????????????????????????????????????????????????????
519      !!    INCLUDE POTENTIAL TEMP -> IN SITU TEMP IN OBS OPERATOR???
520      !!??????????????????????????????????????????????????????????????
521      !!
522      !! ** Action  :
523      !!
524      !! History :
525      !!      ! 2014-08 (J. While) Adapted from obs_pro_opt to handel generalised
526      !!                           vertical coordinates
527      !!-----------------------------------------------------------------------
528   
529      !! * Modules used
530      USE obs_profiles_def   ! Definition of storage space for profile obs.
531      USE dom_oce,  ONLY : & 
532#if defined key_vvl 
533      &   gdepw_n 
534#else 
535      &   gdepw_0 
536#endif
537       
538      IMPLICIT NONE 
539 
540      !! * Arguments
541      TYPE(obs_prof), INTENT(INOUT) :: prodatqc   ! Subset of profile data not failing screening
542      INTEGER, INTENT(IN) :: kt        ! Time step
543      INTEGER, INTENT(IN) :: kpi       ! Model grid parameters
544      INTEGER, INTENT(IN) :: kpj 
545      INTEGER, INTENT(IN) :: kpk 
546      INTEGER, INTENT(IN) :: kit000    ! Number of the first time step 
547                                       !   (kit000-1 = restart time)
548      INTEGER, INTENT(IN) :: k1dint    ! Vertical interpolation type (see header)
549      INTEGER, INTENT(IN) :: k2dint    ! Horizontal interpolation type (see header)
550      INTEGER, INTENT(IN) :: kdaystp   ! Number of time steps per day                     
551      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 
552         & ptn,    &    ! Model temperature field
553         & psn,    &    ! Model salinity field
554         & ptmask       ! Land-sea mask
555      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 
556         & pgdept, &       ! Model array of depth T levels
557         & pgdepw       ! Model array of depth W levels
558      INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 
559         & kdailyavtypes   ! Types for daily averages
560     
561      !! * Local declarations
562      INTEGER ::   ji 
563      INTEGER ::   jj 
564      INTEGER ::   jk 
565      INTEGER ::   iico, ijco 
566      INTEGER ::   jobs 
567      INTEGER ::   inrc 
568      INTEGER ::   ipro 
569      INTEGER ::   idayend 
570      INTEGER ::   ista 
571      INTEGER ::   iend 
572      INTEGER ::   iobs 
573      INTEGER ::   iin, ijn, ikn, ik   ! looping indices over interpolation nodes
574      INTEGER, DIMENSION(imaxavtypes) :: & 
575         & idailyavtypes 
576      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 
577         & igrdi, & 
578         & igrdj 
579      INTEGER :: & 
580         & inum_obs
581      INTEGER, ALLOCATABLE, DIMENSION(:) :: iv_indic   
582      REAL(KIND=wp) :: zlam 
583      REAL(KIND=wp) :: zphi 
584      REAL(KIND=wp) :: zdaystp 
585      REAL(KIND=wp), DIMENSION(kpk) :: & 
586         & zobsmask, & 
587         & zobsk,    & 
588         & zobs2k 
589      REAL(KIND=wp), DIMENSION(2,2,1) :: & 
590         & zweig, & 
591         & l_zweig 
592      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: & 
593         & zmask, & 
594         & zintt, & 
595         & zints, & 
596         & zinmt, & 
597         & zgdept,& 
598         & zgdepw,& 
599         & zinms 
600      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 
601         & zglam, & 
602         & zgphi   
603      REAL(KIND=wp), DIMENSION(1) :: zmsk_1       
604      REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner       
605 
606      !------------------------------------------------------------------------
607      ! Local initialization 
608      !------------------------------------------------------------------------
609      ! ... Record and data counters
610      inrc = kt - kit000 + 2 
611      ipro = prodatqc%npstp(inrc) 
612 
613      ! Daily average types
614      IF ( PRESENT(kdailyavtypes) ) THEN
615         idailyavtypes(:) = kdailyavtypes(:) 
616      ELSE
617         idailyavtypes(:) = -1 
618      ENDIF 
619 
620      ! Initialize daily mean for first time-step
621      idayend = MOD( kt - kit000 + 1, kdaystp ) 
622 
623      ! Added kt == 0 test to catch restart case 
624      IF ( idayend == 1 .OR. kt == 0) THEN
625         
626         IF (lwp) WRITE(numout,*) 'Reset prodatqc%vdmean on time-step: ',kt 
627         DO jk = 1, jpk 
628            DO jj = 1, jpj 
629               DO ji = 1, jpi 
630                  prodatqc%vdmean(ji,jj,jk,1) = 0.0 
631                  prodatqc%vdmean(ji,jj,jk,2) = 0.0 
632               END DO
633            END DO
634         END DO
635       
636      ENDIF
637       
638      DO jk = 1, jpk 
639         DO jj = 1, jpj 
640            DO ji = 1, jpi 
641               ! Increment the temperature field for computing daily mean
642               prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) & 
643               &                        + ptn(ji,jj,jk) 
644               ! Increment the salinity field for computing daily mean
645               prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 
646               &                        + psn(ji,jj,jk) 
647            END DO
648         END DO
649      END DO 
650   
651      ! Compute the daily mean at the end of day
652      zdaystp = 1.0 / REAL( kdaystp ) 
653      IF ( idayend == 0 ) THEN
654         DO jk = 1, jpk 
655            DO jj = 1, jpj 
656               DO ji = 1, jpi 
657                  prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) & 
658                  &                        * zdaystp 
659                  prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 
660                  &                           * zdaystp 
661               END DO
662            END DO
663         END DO
664      ENDIF 
665 
666      ! Get the data for interpolation
667      ALLOCATE( & 
668         & igrdi(2,2,ipro),      & 
669         & igrdj(2,2,ipro),      & 
670         & zglam(2,2,ipro),      & 
671         & zgphi(2,2,ipro),      & 
672         & zmask(2,2,kpk,ipro),  & 
673         & zintt(2,2,kpk,ipro),  & 
674         & zints(2,2,kpk,ipro),  & 
675         & zgdept(2,2,kpk,ipro), & 
676         & zgdepw(2,2,kpk,ipro)  & 
677         & ) 
678 
679      DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro 
680         iobs = jobs - prodatqc%nprofup 
681         igrdi(1,1,iobs) = prodatqc%mi(jobs,1)-1 
682         igrdj(1,1,iobs) = prodatqc%mj(jobs,1)-1 
683         igrdi(1,2,iobs) = prodatqc%mi(jobs,1)-1 
684         igrdj(1,2,iobs) = prodatqc%mj(jobs,1) 
685         igrdi(2,1,iobs) = prodatqc%mi(jobs,1) 
686         igrdj(2,1,iobs) = prodatqc%mj(jobs,1)-1 
687         igrdi(2,2,iobs) = prodatqc%mi(jobs,1) 
688         igrdj(2,2,iobs) = prodatqc%mj(jobs,1) 
689      END DO 
690 
691      ! Initiialise depth arrays
692      zgdept = 0.0
693      zgdepw = 0.0
694     
695      CALL obs_int_comm_2d( 2, 2, ipro, igrdi, igrdj, glamt, zglam ) 
696      CALL obs_int_comm_2d( 2, 2, ipro, igrdi, igrdj, gphit, zgphi ) 
697      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, ptmask,zmask ) 
698      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, ptn,   zintt ) 
699      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, psn,   zints ) 
700      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, pgdept(:,:,:), & 
701        &                     zgdept ) 
702      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, pgdepw(:,:,:), & 
703        &                     zgdepw ) 
704 
705      ! At the end of the day also get interpolated means
706      IF ( idayend == 0 ) THEN
707 
708         ALLOCATE( & 
709            & zinmt(2,2,kpk,ipro),  & 
710            & zinms(2,2,kpk,ipro)   & 
711            & ) 
712 
713         CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, & 
714            &                  prodatqc%vdmean(:,:,:,1), zinmt ) 
715         CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, & 
716            &                  prodatqc%vdmean(:,:,:,2), zinms ) 
717 
718      ENDIF 
719       
720      ! Return if no observations to process
721      ! Has to be done after comm commands to ensure processors
722      ! stay in sync
723      IF ( ipro == 0 ) RETURN
724 
725      DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro 
726   
727         iobs = jobs - prodatqc%nprofup 
728   
729         IF ( kt /= prodatqc%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    = ', prodatqc%mstp(jobs), & 
741                  &            ' ntyp    = ', prodatqc%ntyp(jobs) 
742            ENDIF
743            CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' ) 
744         ENDIF
745         
746         zlam = prodatqc%rlam(jobs) 
747         zphi = prodatqc%rphi(jobs) 
748         
749         ! Horizontal weights
750         ! Only calculated once, for both T and S.
751         ! Masked values are calculated later. 
752 
753         IF ( ( prodatqc%npvend(jobs,1) > 0 ) .OR. & 
754            & ( prodatqc%npvend(jobs,2) > 0 ) ) THEN
755 
756            CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,     & 
757               &                   zglam(:,:,iobs), zgphi(:,:,iobs), & 
758               &                   zmask(:,:,1,iobs), zweig, zmsk_1 ) 
759 
760         ENDIF 
761         
762         ! IF zmsk_1 = 0; then ob is on land
763         IF (zmsk_1(1) < 0.1) THEN
764            WRITE(numout,*) 'WARNING (obs_oper) :- profile found within landmask' 
765   
766         ELSE 
767             
768            ! Temperature
769             
770            IF ( prodatqc%npvend(jobs,1) > 0 ) THEN 
771   
772               zobsk(:) = obfillflt 
773   
774               IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN
775   
776                  IF ( idayend == 0 )  THEN 
777                   
778                     ! Daily averaged moored buoy (MRB) data
779                   
780                     ! vertically interpolate all 4 corners
781                     ista = prodatqc%npvsta(jobs,1) 
782                     iend = prodatqc%npvend(jobs,1) 
783                     inum_obs = iend - ista + 1 
784                     ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs)) 
785     
786                     DO iin=1,2 
787                        DO ijn=1,2 
788                                       
789                                       
790           
791                           IF ( k1dint == 1 ) THEN
792                              CALL obs_int_z1d_spl( kpk, & 
793                                 &     zinmt(iin,ijn,:,iobs), & 
794                                 &     zobs2k, zgdept(iin,ijn,:,iobs), & 
795                                 &     zmask(iin,ijn,:,iobs)) 
796                           ENDIF
797       
798                           CALL obs_level_search(kpk, & 
799                              &    zgdept(iin,ijn,:,iobs), & 
800                              &    inum_obs, prodatqc%var(1)%vdep(ista:iend), & 
801                              &    iv_indic) 
802                           CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 
803                              &    prodatqc%var(1)%vdep(ista:iend), & 
804                              &    zinmt(iin,ijn,:,iobs), & 
805                              &    zobs2k, interp_corner(iin,ijn,:), & 
806                              &    zgdept(iin,ijn,:,iobs), & 
807                              &    zmask(iin,ijn,:,iobs)) 
808       
809                        ENDDO 
810                     ENDDO 
811                   
812                   
813                  ELSE
814               
815                     CALL ctl_stop( ' A nonzero' //     & 
816                        &           ' number of profile T BUOY data should' // & 
817                        &           ' only occur at the end of a given day' ) 
818   
819                  ENDIF
820         
821               ELSE 
822               
823                  ! Point data
824     
825                  ! vertically interpolate all 4 corners
826                  ista = prodatqc%npvsta(jobs,1) 
827                  iend = prodatqc%npvend(jobs,1) 
828                  inum_obs = iend - ista + 1 
829                  ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs)) 
830                  DO iin=1,2 
831                     DO ijn=1,2 
832                                   
833                                   
834                        IF ( k1dint == 1 ) THEN
835                           CALL obs_int_z1d_spl( kpk, & 
836                              &    zintt(iin,ijn,:,iobs),& 
837                              &    zobs2k, zgdept(iin,ijn,:,iobs), & 
838                              &    zmask(iin,ijn,:,iobs)) 
839 
840                        ENDIF
841       
842                        CALL obs_level_search(kpk, & 
843                            &        zgdept(iin,ijn,:,iobs),& 
844                            &        inum_obs, prodatqc%var(1)%vdep(ista:iend), & 
845                            &         iv_indic) 
846                        CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs,     & 
847                            &          prodatqc%var(1)%vdep(ista:iend),     & 
848                            &          zintt(iin,ijn,:,iobs),            & 
849                            &          zobs2k,interp_corner(iin,ijn,:), & 
850                            &          zgdept(iin,ijn,:,iobs),         & 
851                            &          zmask(iin,ijn,:,iobs) )     
852         
853                     ENDDO 
854                  ENDDO 
855             
856               ENDIF 
857       
858               !-------------------------------------------------------------
859               ! Compute the horizontal interpolation for every profile level
860               !-------------------------------------------------------------
861             
862               DO ikn=1,inum_obs 
863                  iend=ista+ikn-1
864                 
865                  l_zweig(:,:,1) = 0._wp 
866   
867                  ! This code forces the horizontal weights to be 
868                  ! zero IF the observation is below the bottom of the 
869                  ! corners of the interpolation nodes, Or if it is in 
870                  ! the mask. This is important for observations are near 
871                  ! steep bathymetry
872                  DO iin=1,2 
873                     DO ijn=1,2 
874     
875                        depth_loop1: DO ik=kpk,2,-1 
876                           IF(zmask(iin,ijn,ik-1,iobs ) > 0.9 )THEN   
877                           
878                              l_zweig(iin,ijn,1) = & 
879                                 & zweig(iin,ijn,1) * & 
880                                 & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) & 
881                                 &  - prodatqc%var(1)%vdep(iend)),0._wp) 
882                           
883                              EXIT depth_loop1 
884                           ENDIF
885                        ENDDO depth_loop1 
886     
887                     ENDDO 
888                  ENDDO 
889   
890                  CALL obs_int_h2d( 1, 1, l_zweig, interp_corner(:,:,ikn), & 
891                  &          prodatqc%var(1)%vmod(iend:iend) ) 
892
893                  ! Set QC flag for any observations found below the bottom
894                  ! needed as the check here is more strict than that in obs_prep
895                  IF (sum(l_zweig) == 0.0_wp) prodatqc%var(1)%nvqc(iend:iend)=4
896 
897               ENDDO 
898 
899 
900               DEALLOCATE(interp_corner,iv_indic) 
901         
902            ENDIF 
903       
904 
905            ! Salinity 
906         
907            IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 
908   
909               zobsk(:) = obfillflt 
910   
911               IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN
912   
913                  IF ( idayend == 0 )  THEN 
914                   
915                     ! Daily averaged moored buoy (MRB) data
916                   
917                     ! vertically interpolate all 4 corners
918                     ista = prodatqc%npvsta(jobs,2) 
919                     iend = prodatqc%npvend(jobs,2) 
920                     inum_obs = iend - ista + 1 
921                     ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs)) 
922     
923                     DO iin=1,2 
924                        DO ijn=1,2 
925                                       
926                                       
927           
928                           IF ( k1dint == 1 ) THEN
929                              CALL obs_int_z1d_spl( kpk, & 
930                                 &     zinms(iin,ijn,:,iobs), & 
931                                 &     zobs2k, zgdept(iin,ijn,:,iobs), & 
932                                 &     zmask(iin,ijn,:,iobs)) 
933                           ENDIF
934       
935                           CALL obs_level_search(kpk, & 
936                              &    zgdept(iin,ijn,:,iobs), & 
937                              &    inum_obs, prodatqc%var(2)%vdep(ista:iend), & 
938                              &    iv_indic) 
939                           CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 
940                              &    prodatqc%var(2)%vdep(ista:iend), & 
941                              &    zinms(iin,ijn,:,iobs), & 
942                              &    zobs2k, interp_corner(iin,ijn,:), & 
943                              &    zgdept(iin,ijn,:,iobs), & 
944                              &    zmask(iin,ijn,:,iobs)) 
945       
946                        ENDDO 
947                     ENDDO 
948                   
949                   
950                  ELSE
951               
952                     CALL ctl_stop( ' A nonzero' //     & 
953                        &           ' number of profile T BUOY data should' // & 
954                        &           ' only occur at the end of a given day' ) 
955   
956                  ENDIF
957         
958               ELSE 
959               
960                  ! Point data
961     
962                  ! vertically interpolate all 4 corners
963                  ista = prodatqc%npvsta(jobs,2) 
964                  iend = prodatqc%npvend(jobs,2) 
965                  inum_obs = iend - ista + 1 
966                  ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs)) 
967                   
968                  DO iin=1,2     
969                     DO ijn=1,2 
970                                 
971                                 
972                        IF ( k1dint == 1 ) THEN
973                           CALL obs_int_z1d_spl( kpk, & 
974                              &    zints(iin,ijn,:,iobs),& 
975                              &    zobs2k, zgdept(iin,ijn,:,iobs), & 
976                              &    zmask(iin,ijn,:,iobs)) 
977 
978                        ENDIF
979       
980                        CALL obs_level_search(kpk, & 
981                           &        zgdept(iin,ijn,:,iobs),& 
982                           &        inum_obs, prodatqc%var(2)%vdep(ista:iend), & 
983                           &         iv_indic) 
984                        CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs,  & 
985                           &          prodatqc%var(2)%vdep(ista:iend),     & 
986                           &          zints(iin,ijn,:,iobs),               & 
987                           &          zobs2k,interp_corner(iin,ijn,:),     & 
988                           &          zgdept(iin,ijn,:,iobs),              & 
989                           &          zmask(iin,ijn,:,iobs) )     
990         
991                     ENDDO 
992                  ENDDO 
993             
994               ENDIF 
995       
996               !-------------------------------------------------------------
997               ! Compute the horizontal interpolation for every profile level
998               !-------------------------------------------------------------
999             
1000               DO ikn=1,inum_obs 
1001                  iend=ista+ikn-1
1002                 
1003                  l_zweig(:,:,1) = 0._wp 
1004   
1005                  ! This code forces the horizontal weights to be 
1006                  ! zero IF the observation is below the bottom of the 
1007                  ! corners of the interpolation nodes, Or if it is in 
1008                  ! the mask. This is important for observations are near 
1009                  ! steep bathymetry
1010                  DO iin=1,2 
1011                     DO ijn=1,2 
1012     
1013                        depth_loop2: DO ik=kpk,2,-1 
1014                           IF(zmask(iin,ijn,ik-1,iobs ) > 0.9 )THEN   
1015                           
1016                              l_zweig(iin,ijn,1) = & 
1017                                 &  zweig(iin,ijn,1) * & 
1018                                 &  MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) & 
1019                                 &  - prodatqc%var(2)%vdep(iend)),0._wp) 
1020                           
1021                              EXIT depth_loop2 
1022                           ENDIF
1023                        ENDDO depth_loop2 
1024     
1025                     ENDDO 
1026                  ENDDO 
1027   
1028                  CALL obs_int_h2d( 1, 1, l_zweig, interp_corner(:,:,ikn), & 
1029                  &          prodatqc%var(2)%vmod(iend:iend) ) 
1030
1031                  ! Set QC flag for any observations found below the bottom
1032                  ! needed as the check here is more strict than that in obs_prep
1033                  IF (sum(l_zweig) == 0.0_wp) prodatqc%var(2)%nvqc(iend:iend)=4
1034 
1035               ENDDO 
1036 
1037 
1038               DEALLOCATE(interp_corner,iv_indic) 
1039         
1040            ENDIF
1041         
1042         ENDIF
1043       
1044      END DO 
1045     
1046      ! Deallocate the data for interpolation
1047      DEALLOCATE( & 
1048         & igrdi, & 
1049         & igrdj, & 
1050         & zglam, & 
1051         & zgphi, & 
1052         & zmask, & 
1053         & zintt, & 
1054         & zints, & 
1055         & zgdept,& 
1056         & zgdepw & 
1057         & ) 
1058      ! At the end of the day also get interpolated means
1059      IF ( idayend == 0 ) THEN
1060         DEALLOCATE( & 
1061            & zinmt,  & 
1062            & zinms   & 
1063            & ) 
1064      ENDIF
1065   
1066      prodatqc%nprofup = prodatqc%nprofup + ipro 
1067       
1068   END SUBROUTINE obs_pro_sco_opt 
1069 
1070   SUBROUTINE obs_sla_opt( sladatqc, kt, kpi, kpj, kit000, &
1071      &                    psshn, psshmask, k2dint )
1072      !!-----------------------------------------------------------------------
1073      !!
1074      !!                     ***  ROUTINE obs_sla_opt  ***
1075      !!
1076      !! ** Purpose : Compute the model counterpart of sea level anomaly
1077      !!              data by interpolating from the model grid to the
1078      !!              observation point.
1079      !!
1080      !! ** Method  : Linearly interpolate to each observation point using
1081      !!              the model values at the corners of the surrounding grid box.
1082      !!
1083      !!    The now model SSH is first computed at the obs (lon, lat) point.
1084      !!
1085      !!    Several horizontal interpolation schemes are available:
1086      !!        - distance-weighted (great circle) (k2dint = 0)
1087      !!        - distance-weighted (small angle)  (k2dint = 1)
1088      !!        - bilinear (geographical grid)     (k2dint = 2)
1089      !!        - bilinear (quadrilateral grid)    (k2dint = 3)
1090      !!        - polynomial (quadrilateral grid)  (k2dint = 4)
1091      !! 
1092      !!    The sea level anomaly at the observation points is then computed
1093      !!    by removing a mean dynamic topography (defined at the obs. point).
1094      !!
1095      !! ** Action  :
1096      !!
1097      !! History :
1098      !!      ! 07-03 (A. Weaver)
1099      !!-----------------------------------------------------------------------
1100 
1101      !! * Modules used
1102      USE obs_surf_def  ! Definition of storage space for surface observations
1103
1104      IMPLICIT NONE
1105
1106      !! * Arguments
1107      TYPE(obs_surf), INTENT(INOUT) :: sladatqc     ! Subset of surface data not failing screening
1108      INTEGER, INTENT(IN) :: kt      ! Time step
1109      INTEGER, INTENT(IN) :: kpi     ! Model grid parameters
1110      INTEGER, INTENT(IN) :: kpj
1111      INTEGER, INTENT(IN) :: kit000   ! Number of the first time step
1112                                      !   (kit000-1 = restart time)
1113      INTEGER, INTENT(IN) :: k2dint   ! Horizontal interpolation type (see header)
1114      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
1115         & psshn,  &    ! Model SSH field
1116         & psshmask     ! Land-sea mask
1117         
1118      !! * Local declarations
1119      INTEGER :: ji
1120      INTEGER :: jj
1121      INTEGER :: jobs
1122      INTEGER :: inrc
1123      INTEGER :: isla
1124      INTEGER :: iobs
1125      REAL(KIND=wp) :: zlam
1126      REAL(KIND=wp) :: zphi
1127      REAL(KIND=wp) :: zext(1), zobsmask(1)
1128      REAL(kind=wp), DIMENSION(2,2,1) :: &
1129         & zweig
1130      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
1131         & zmask, &
1132         & zsshl, &
1133         & zglam, &
1134         & zgphi
1135      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
1136         & igrdi, &
1137         & igrdj
1138
1139      !------------------------------------------------------------------------
1140      ! Local initialization
1141      !------------------------------------------------------------------------
1142      ! ... Record and data counters
1143      inrc = kt - kit000 + 2
1144      isla = sladatqc%nsstp(inrc)
1145
1146      ! Get the data for interpolation
1147
1148      ALLOCATE( &
1149         & igrdi(2,2,isla), &
1150         & igrdj(2,2,isla), &
1151         & zglam(2,2,isla), &
1152         & zgphi(2,2,isla), &
1153         & zmask(2,2,isla), &
1154         & zsshl(2,2,isla)  &
1155         & )
1156     
1157      DO jobs = sladatqc%nsurfup + 1, sladatqc%nsurfup + isla
1158         iobs = jobs - sladatqc%nsurfup
1159         igrdi(1,1,iobs) = sladatqc%mi(jobs)-1
1160         igrdj(1,1,iobs) = sladatqc%mj(jobs)-1
1161         igrdi(1,2,iobs) = sladatqc%mi(jobs)-1
1162         igrdj(1,2,iobs) = sladatqc%mj(jobs)
1163         igrdi(2,1,iobs) = sladatqc%mi(jobs)
1164         igrdj(2,1,iobs) = sladatqc%mj(jobs)-1
1165         igrdi(2,2,iobs) = sladatqc%mi(jobs)
1166         igrdj(2,2,iobs) = sladatqc%mj(jobs)
1167      END DO
1168
1169      CALL obs_int_comm_2d( 2, 2, isla, &
1170         &                  igrdi, igrdj, glamt, zglam )
1171      CALL obs_int_comm_2d( 2, 2, isla, &
1172         &                  igrdi, igrdj, gphit, zgphi )
1173      CALL obs_int_comm_2d( 2, 2, isla, &
1174         &                  igrdi, igrdj, psshmask, zmask )
1175      CALL obs_int_comm_2d( 2, 2, isla, &
1176         &                  igrdi, igrdj, psshn, zsshl )
1177
1178      ! Loop over observations
1179
1180      DO jobs = sladatqc%nsurfup + 1, sladatqc%nsurfup + isla
1181
1182         iobs = jobs - sladatqc%nsurfup
1183
1184         IF ( kt /= sladatqc%mstp(jobs) ) THEN
1185           
1186            IF(lwp) THEN
1187               WRITE(numout,*)
1188               WRITE(numout,*) ' E R R O R : Observation',              &
1189                  &            ' time step is not consistent with the', &
1190                  &            ' model time step'
1191               WRITE(numout,*) ' ========='
1192               WRITE(numout,*)
1193               WRITE(numout,*) ' Record  = ', jobs,                &
1194                  &            ' kt      = ', kt,                  &
1195                  &            ' mstp    = ', sladatqc%mstp(jobs), &
1196                  &            ' ntyp    = ', sladatqc%ntyp(jobs)
1197            ENDIF
1198            CALL ctl_stop( 'obs_sla_opt', 'Inconsistent time' )
1199           
1200         ENDIF
1201         
1202         zlam = sladatqc%rlam(jobs)
1203         zphi = sladatqc%rphi(jobs)
1204
1205         ! Get weights to interpolate the model SSH to the observation point
1206         CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         &
1207            &                   zglam(:,:,iobs), zgphi(:,:,iobs), &
1208            &                   zmask(:,:,iobs), zweig, zobsmask )
1209         
1210
1211         ! Interpolate the model SSH to the observation point
1212         CALL obs_int_h2d( 1, 1,      &
1213            &              zweig, zsshl(:,:,iobs),  zext )
1214         
1215         sladatqc%rext(jobs,1) = zext(1)
1216         ! ... Remove the MDT at the observation point
1217         sladatqc%rmod(jobs,1) = sladatqc%rext(jobs,1) - sladatqc%rext(jobs,2)
1218
1219      END DO
1220
1221      ! Deallocate the data for interpolation
1222      DEALLOCATE( &
1223         & igrdi, &
1224         & igrdj, &
1225         & zglam, &
1226         & zgphi, &
1227         & zmask, &
1228         & zsshl  &
1229         & )
1230
1231      sladatqc%nsurfup = sladatqc%nsurfup + isla
1232
1233   END SUBROUTINE obs_sla_opt
1234
1235   SUBROUTINE obs_sst_opt( sstdatqc, kt, kpi, kpj, kit000, kdaystp, &
1236      &                    psstn, psstmask, k2dint, ld_nightav )
1237      !!-----------------------------------------------------------------------
1238      !!
1239      !!                     ***  ROUTINE obs_sst_opt  ***
1240      !!
1241      !! ** Purpose : Compute the model counterpart of surface temperature
1242      !!              data by interpolating from the model grid to the
1243      !!              observation point.
1244      !!
1245      !! ** Method  : Linearly interpolate to each observation point using
1246      !!              the model values at the corners of the surrounding grid box.
1247      !!
1248      !!    The now model SST is first computed at the obs (lon, lat) point.
1249      !!
1250      !!    Several horizontal interpolation schemes are available:
1251      !!        - distance-weighted (great circle) (k2dint = 0)
1252      !!        - distance-weighted (small angle)  (k2dint = 1)
1253      !!        - bilinear (geographical grid)     (k2dint = 2)
1254      !!        - bilinear (quadrilateral grid)    (k2dint = 3)
1255      !!        - polynomial (quadrilateral grid)  (k2dint = 4)
1256      !!
1257      !!
1258      !! ** Action  :
1259      !!
1260      !! History :
1261      !!        !  07-07  (S. Ricci ) : Original
1262      !!     
1263      !!-----------------------------------------------------------------------
1264
1265      !! * Modules used
1266      USE obs_surf_def  ! Definition of storage space for surface observations
1267      USE sbcdcy
1268
1269      IMPLICIT NONE
1270
1271      !! * Arguments
1272      TYPE(obs_surf), INTENT(INOUT) :: &
1273         & sstdatqc     ! Subset of surface data not failing screening
1274      INTEGER, INTENT(IN) :: kt        ! Time step
1275      INTEGER, INTENT(IN) :: kpi       ! Model grid parameters
1276      INTEGER, INTENT(IN) :: kpj
1277      INTEGER, INTENT(IN) :: kit000    ! Number of the first time step
1278                                       !   (kit000-1 = restart time)
1279      INTEGER, INTENT(IN) :: k2dint    ! Horizontal interpolation type (see header)
1280      INTEGER, INTENT(IN) :: kdaystp   ! Number of time steps per day 
1281      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
1282         & psstn,  &    ! Model SST field
1283         & psstmask     ! Land-sea mask
1284
1285      !! * Local declarations
1286      INTEGER :: ji
1287      INTEGER :: jj
1288      INTEGER :: jobs
1289      INTEGER :: inrc
1290      INTEGER :: isst
1291      INTEGER :: iobs
1292      INTEGER :: idayend
1293      REAL(KIND=wp) :: zlam
1294      REAL(KIND=wp) :: zphi
1295      REAL(KIND=wp) :: zext(1), zobsmask(1)
1296      REAL(KIND=wp) :: zdaystp
1297      INTEGER, DIMENSION(:,:), SAVE, ALLOCATABLE :: &
1298         & icount_sstnight,      &
1299         & imask_night
1300      REAL(kind=wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: &
1301         & zintmp, &
1302         & zouttmp, & 
1303         & zmeanday    ! to compute model sst in region of 24h daylight (pole)
1304      REAL(kind=wp), DIMENSION(2,2,1) :: &
1305         & zweig
1306      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
1307         & zmask, &
1308         & zsstl, &
1309         & zsstm, &
1310         & zglam, &
1311         & zgphi
1312      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
1313         & igrdi, &
1314         & igrdj
1315      LOGICAL, INTENT(IN) :: ld_nightav
1316
1317      !-----------------------------------------------------------------------
1318      ! Local initialization
1319      !-----------------------------------------------------------------------
1320      ! ... Record and data counters
1321      inrc = kt - kit000 + 2
1322      isst = sstdatqc%nsstp(inrc)
1323
1324      IF ( ld_nightav ) THEN
1325
1326      ! Initialize array for night mean
1327
1328      IF ( kt .EQ. 0 ) THEN
1329         ALLOCATE ( icount_sstnight(kpi,kpj) )
1330         ALLOCATE ( imask_night(kpi,kpj) )
1331         ALLOCATE ( zintmp(kpi,kpj) )
1332         ALLOCATE ( zouttmp(kpi,kpj) )
1333         ALLOCATE ( zmeanday(kpi,kpj) )
1334         nday_qsr = -1   ! initialisation flag for nbc_dcy
1335      ENDIF
1336
1337      ! Initialize daily mean for first timestep
1338      idayend = MOD( kt - kit000 + 1, kdaystp )
1339
1340      ! Added kt == 0 test to catch restart case
1341      IF ( idayend == 1 .OR. kt == 0) THEN
1342         IF (lwp) WRITE(numout,*) 'Reset sstdatqc%vdmean on time-step: ',kt
1343         DO jj = 1, jpj
1344            DO ji = 1, jpi
1345               sstdatqc%vdmean(ji,jj) = 0.0
1346               zmeanday(ji,jj) = 0.0
1347               icount_sstnight(ji,jj) = 0
1348            END DO
1349         END DO
1350      ENDIF
1351
1352      zintmp(:,:) = 0.0
1353      zouttmp(:,:) = sbc_dcy( zintmp(:,:), .TRUE. )
1354      imask_night(:,:) = INT( zouttmp(:,:) )
1355
1356      DO jj = 1, jpj
1357         DO ji = 1, jpi
1358            ! Increment the temperature field for computing night mean and counter
1359            sstdatqc%vdmean(ji,jj) = sstdatqc%vdmean(ji,jj)  &
1360                   &                        + psstn(ji,jj)*imask_night(ji,jj)
1361            zmeanday(ji,jj)        = zmeanday(ji,jj) + psstn(ji,jj)
1362            icount_sstnight(ji,jj) = icount_sstnight(ji,jj) + imask_night(ji,jj)
1363         END DO
1364      END DO
1365   
1366      ! Compute the daily mean at the end of day
1367
1368      zdaystp = 1.0 / REAL( kdaystp )
1369
1370      IF ( idayend == 0 ) THEN
1371         DO jj = 1, jpj
1372            DO ji = 1, jpi
1373               ! Test if "no night" point
1374               IF ( icount_sstnight(ji,jj) .NE. 0 ) THEN
1375                  sstdatqc%vdmean(ji,jj) = sstdatqc%vdmean(ji,jj) &
1376                    &                        / icount_sstnight(ji,jj) 
1377               ELSE
1378                  sstdatqc%vdmean(ji,jj) = zmeanday(ji,jj) * zdaystp
1379               ENDIF
1380            END DO
1381         END DO
1382      ENDIF
1383
1384      ENDIF
1385
1386      ! Get the data for interpolation
1387     
1388      ALLOCATE( &
1389         & igrdi(2,2,isst), &
1390         & igrdj(2,2,isst), &
1391         & zglam(2,2,isst), &
1392         & zgphi(2,2,isst), &
1393         & zmask(2,2,isst), &
1394         & zsstl(2,2,isst)  &
1395         & )
1396     
1397      DO jobs = sstdatqc%nsurfup + 1, sstdatqc%nsurfup + isst
1398         iobs = jobs - sstdatqc%nsurfup
1399         igrdi(1,1,iobs) = sstdatqc%mi(jobs)-1
1400         igrdj(1,1,iobs) = sstdatqc%mj(jobs)-1
1401         igrdi(1,2,iobs) = sstdatqc%mi(jobs)-1
1402         igrdj(1,2,iobs) = sstdatqc%mj(jobs)
1403         igrdi(2,1,iobs) = sstdatqc%mi(jobs)
1404         igrdj(2,1,iobs) = sstdatqc%mj(jobs)-1
1405         igrdi(2,2,iobs) = sstdatqc%mi(jobs)
1406         igrdj(2,2,iobs) = sstdatqc%mj(jobs)
1407      END DO
1408     
1409      CALL obs_int_comm_2d( 2, 2, isst, &
1410         &                  igrdi, igrdj, glamt, zglam )
1411      CALL obs_int_comm_2d( 2, 2, isst, &
1412         &                  igrdi, igrdj, gphit, zgphi )
1413      CALL obs_int_comm_2d( 2, 2, isst, &
1414         &                  igrdi, igrdj, psstmask, zmask )
1415      CALL obs_int_comm_2d( 2, 2, isst, &
1416         &                  igrdi, igrdj, psstn, zsstl )
1417
1418      ! At the end of the day get interpolated means
1419      IF ( idayend == 0 .AND. ld_nightav ) THEN
1420
1421         ALLOCATE( &
1422            & zsstm(2,2,isst)  &
1423            & )
1424
1425         CALL obs_int_comm_2d( 2, 2, isst, igrdi, igrdj, &
1426            &               sstdatqc%vdmean(:,:), zsstm )
1427
1428      ENDIF
1429
1430      ! Loop over observations
1431
1432      DO jobs = sstdatqc%nsurfup + 1, sstdatqc%nsurfup + isst
1433         
1434         iobs = jobs - sstdatqc%nsurfup
1435         
1436         IF ( kt /= sstdatqc%mstp(jobs) ) THEN
1437           
1438            IF(lwp) THEN
1439               WRITE(numout,*)
1440               WRITE(numout,*) ' E R R O R : Observation',              &
1441                  &            ' time step is not consistent with the', &
1442                  &            ' model time step'
1443               WRITE(numout,*) ' ========='
1444               WRITE(numout,*)
1445               WRITE(numout,*) ' Record  = ', jobs,                &
1446                  &            ' kt      = ', kt,                  &
1447                  &            ' mstp    = ', sstdatqc%mstp(jobs), &
1448                  &            ' ntyp    = ', sstdatqc%ntyp(jobs)
1449            ENDIF
1450            CALL ctl_stop( 'obs_sst_opt', 'Inconsistent time' )
1451           
1452         ENDIF
1453         
1454         zlam = sstdatqc%rlam(jobs)
1455         zphi = sstdatqc%rphi(jobs)
1456         
1457         ! Get weights to interpolate the model SST to the observation point
1458         CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         &
1459            &                   zglam(:,:,iobs), zgphi(:,:,iobs), &
1460            &                   zmask(:,:,iobs), zweig, zobsmask )
1461           
1462         ! Interpolate the model SST to the observation point
1463
1464         IF ( ld_nightav ) THEN
1465
1466           IF ( idayend == 0 )  THEN
1467               ! Daily averaged/diurnal cycle of SST  data
1468               CALL obs_int_h2d( 1, 1,      & 
1469                     &              zweig, zsstm(:,:,iobs), zext )
1470            ELSE
1471               CALL ctl_stop( ' ld_nightav is set to true: a nonzero' //     &
1472                     &           ' number of night SST data should' // &
1473                     &           ' only occur at the end of a given day' )
1474            ENDIF
1475
1476         ELSE
1477
1478            CALL obs_int_h2d( 1, 1,      &
1479            &              zweig, zsstl(:,:,iobs),  zext )
1480
1481         ENDIF
1482         sstdatqc%rmod(jobs,1) = zext(1)
1483         
1484      END DO
1485     
1486      ! Deallocate the data for interpolation
1487      DEALLOCATE( &
1488         & igrdi, &
1489         & igrdj, &
1490         & zglam, &
1491         & zgphi, &
1492         & zmask, &
1493         & zsstl  &
1494         & )
1495
1496      ! At the end of the day also get interpolated means
1497      IF ( idayend == 0 .AND. ld_nightav ) THEN
1498         DEALLOCATE( &
1499            & zsstm  &
1500            & )
1501      ENDIF
1502     
1503      sstdatqc%nsurfup = sstdatqc%nsurfup + isst
1504
1505   END SUBROUTINE obs_sst_opt
1506
1507   SUBROUTINE obs_sss_opt
1508      !!-----------------------------------------------------------------------
1509      !!
1510      !!                     ***  ROUTINE obs_sss_opt  ***
1511      !!
1512      !! ** Purpose : Compute the model counterpart of sea surface salinity
1513      !!              data by interpolating from the model grid to the
1514      !!              observation point.
1515      !!
1516      !! ** Method  :
1517      !!
1518      !! ** Action  :
1519      !!
1520      !! History :
1521      !!      ! ??-??
1522      !!-----------------------------------------------------------------------
1523
1524      IMPLICIT NONE
1525
1526   END SUBROUTINE obs_sss_opt
1527
1528   SUBROUTINE obs_seaice_opt( seaicedatqc, kt, kpi, kpj, kit000, &
1529      &                    pseaicen, pseaicemask, k2dint )
1530
1531      !!-----------------------------------------------------------------------
1532      !!
1533      !!                     ***  ROUTINE obs_seaice_opt  ***
1534      !!
1535      !! ** Purpose : Compute the model counterpart of surface temperature
1536      !!              data by interpolating from the model grid to the
1537      !!              observation point.
1538      !!
1539      !! ** Method  : Linearly interpolate to each observation point using
1540      !!              the model values at the corners of the surrounding grid box.
1541      !!
1542      !!    The now model sea ice is first computed at the obs (lon, lat) point.
1543      !!
1544      !!    Several horizontal interpolation schemes are available:
1545      !!        - distance-weighted (great circle) (k2dint = 0)
1546      !!        - distance-weighted (small angle)  (k2dint = 1)
1547      !!        - bilinear (geographical grid)     (k2dint = 2)
1548      !!        - bilinear (quadrilateral grid)    (k2dint = 3)
1549      !!        - polynomial (quadrilateral grid)  (k2dint = 4)
1550      !!
1551      !!
1552      !! ** Action  :
1553      !!
1554      !! History :
1555      !!        !  07-07  (S. Ricci ) : Original
1556      !!     
1557      !!-----------------------------------------------------------------------
1558
1559      !! * Modules used
1560      USE obs_surf_def  ! Definition of storage space for surface observations
1561
1562      IMPLICIT NONE
1563
1564      !! * Arguments
1565      TYPE(obs_surf), INTENT(INOUT) :: seaicedatqc     ! Subset of surface data not failing screening
1566      INTEGER, INTENT(IN) :: kt       ! Time step
1567      INTEGER, INTENT(IN) :: kpi      ! Model grid parameters
1568      INTEGER, INTENT(IN) :: kpj
1569      INTEGER, INTENT(IN) :: kit000   ! Number of the first time step
1570                                      !   (kit000-1 = restart time)
1571      INTEGER, INTENT(IN) :: k2dint   ! Horizontal interpolation type (see header)
1572      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
1573         & pseaicen,  &    ! Model sea ice field
1574         & pseaicemask     ! Land-sea mask
1575         
1576      !! * Local declarations
1577      INTEGER :: ji
1578      INTEGER :: jj
1579      INTEGER :: jobs
1580      INTEGER :: inrc
1581      INTEGER :: iseaice
1582      INTEGER :: iobs
1583       
1584      REAL(KIND=wp) :: zlam
1585      REAL(KIND=wp) :: zphi
1586      REAL(KIND=wp) :: zext(1), zobsmask(1)
1587      REAL(kind=wp), DIMENSION(2,2,1) :: &
1588         & zweig
1589      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
1590         & zmask, &
1591         & zseaicel, &
1592         & zglam, &
1593         & zgphi
1594      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
1595         & igrdi, &
1596         & igrdj
1597
1598      !------------------------------------------------------------------------
1599      ! Local initialization
1600      !------------------------------------------------------------------------
1601      ! ... Record and data counters
1602      inrc = kt - kit000 + 2
1603      iseaice = seaicedatqc%nsstp(inrc)
1604
1605      ! Get the data for interpolation
1606     
1607      ALLOCATE( &
1608         & igrdi(2,2,iseaice), &
1609         & igrdj(2,2,iseaice), &
1610         & zglam(2,2,iseaice), &
1611         & zgphi(2,2,iseaice), &
1612         & zmask(2,2,iseaice), &
1613         & zseaicel(2,2,iseaice)  &
1614         & )
1615     
1616      DO jobs = seaicedatqc%nsurfup + 1, seaicedatqc%nsurfup + iseaice
1617         iobs = jobs - seaicedatqc%nsurfup
1618         igrdi(1,1,iobs) = seaicedatqc%mi(jobs)-1
1619         igrdj(1,1,iobs) = seaicedatqc%mj(jobs)-1
1620         igrdi(1,2,iobs) = seaicedatqc%mi(jobs)-1
1621         igrdj(1,2,iobs) = seaicedatqc%mj(jobs)
1622         igrdi(2,1,iobs) = seaicedatqc%mi(jobs)
1623         igrdj(2,1,iobs) = seaicedatqc%mj(jobs)-1
1624         igrdi(2,2,iobs) = seaicedatqc%mi(jobs)
1625         igrdj(2,2,iobs) = seaicedatqc%mj(jobs)
1626      END DO
1627     
1628      CALL obs_int_comm_2d( 2, 2, iseaice, &
1629         &                  igrdi, igrdj, glamt, zglam )
1630      CALL obs_int_comm_2d( 2, 2, iseaice, &
1631         &                  igrdi, igrdj, gphit, zgphi )
1632      CALL obs_int_comm_2d( 2, 2, iseaice, &
1633         &                  igrdi, igrdj, pseaicemask, zmask )
1634      CALL obs_int_comm_2d( 2, 2, iseaice, &
1635         &                  igrdi, igrdj, pseaicen, zseaicel )
1636     
1637      DO jobs = seaicedatqc%nsurfup + 1, seaicedatqc%nsurfup + iseaice
1638         
1639         iobs = jobs - seaicedatqc%nsurfup
1640         
1641         IF ( kt /= seaicedatqc%mstp(jobs) ) THEN
1642           
1643            IF(lwp) THEN
1644               WRITE(numout,*)
1645               WRITE(numout,*) ' E R R O R : Observation',              &
1646                  &            ' time step is not consistent with the', &
1647                  &            ' model time step'
1648               WRITE(numout,*) ' ========='
1649               WRITE(numout,*)
1650               WRITE(numout,*) ' Record  = ', jobs,                &
1651                  &            ' kt      = ', kt,                  &
1652                  &            ' mstp    = ', seaicedatqc%mstp(jobs), &
1653                  &            ' ntyp    = ', seaicedatqc%ntyp(jobs)
1654            ENDIF
1655            CALL ctl_stop( 'obs_seaice_opt', 'Inconsistent time' )
1656           
1657         ENDIF
1658         
1659         zlam = seaicedatqc%rlam(jobs)
1660         zphi = seaicedatqc%rphi(jobs)
1661         
1662         ! Get weights to interpolate the model sea ice to the observation point
1663         CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         &
1664            &                   zglam(:,:,iobs), zgphi(:,:,iobs), &
1665            &                   zmask(:,:,iobs), zweig, zobsmask )
1666         
1667         ! ... Interpolate the model sea ice to the observation point
1668         CALL obs_int_h2d( 1, 1,      &
1669            &              zweig, zseaicel(:,:,iobs),  zext )
1670         
1671         seaicedatqc%rmod(jobs,1) = zext(1)
1672         
1673      END DO
1674     
1675      ! Deallocate the data for interpolation
1676      DEALLOCATE( &
1677         & igrdi,    &
1678         & igrdj,    &
1679         & zglam,    &
1680         & zgphi,    &
1681         & zmask,    &
1682         & zseaicel  &
1683         & )
1684     
1685      seaicedatqc%nsurfup = seaicedatqc%nsurfup + iseaice
1686
1687   END SUBROUTINE obs_seaice_opt
1688
1689   SUBROUTINE obs_vel_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, &
1690      &                    pun, pvn, pgdept, pumask, pvmask, k1dint, k2dint, &
1691      &                    ld_dailyav )
1692      !!-----------------------------------------------------------------------
1693      !!
1694      !!                     ***  ROUTINE obs_vel_opt  ***
1695      !!
1696      !! ** Purpose : Compute the model counterpart of velocity profile
1697      !!              data by interpolating from the model grid to the
1698      !!              observation point.
1699      !!
1700      !! ** Method  : Linearly interpolate zonal and meridional components of velocity
1701      !!              to each observation point using the model values at the corners of
1702      !!              the surrounding grid box. The model velocity components are on a
1703      !!              staggered C- grid.
1704      !!
1705      !!    For velocity data from the TAO array, the model equivalent is
1706      !!    a daily mean velocity field. So, we first compute
1707      !!    the mean, then interpolate only at the end of the day.
1708      !!
1709      !! ** Action  :
1710      !!
1711      !! History :
1712      !!    ! 07-03 (K. Mogensen)      : Temperature and Salinity profiles
1713      !!    ! 08-10 (Maria Valdivieso) : Velocity component (U,V) profiles
1714      !!-----------------------------------------------------------------------
1715   
1716      !! * Modules used
1717      USE obs_profiles_def ! Definition of storage space for profile obs.
1718
1719      IMPLICIT NONE
1720
1721      !! * Arguments
1722      TYPE(obs_prof), INTENT(INOUT) :: &
1723         & prodatqc        ! Subset of profile data not failing screening
1724      INTEGER, INTENT(IN) :: kt        ! Time step
1725      INTEGER, INTENT(IN) :: kpi       ! Model grid parameters
1726      INTEGER, INTENT(IN) :: kpj
1727      INTEGER, INTENT(IN) :: kpk 
1728      INTEGER, INTENT(IN) :: kit000    ! Number of the first time step
1729                                       !   (kit000-1 = restart time)
1730      INTEGER, INTENT(IN) :: k1dint    ! Vertical interpolation type (see header)
1731      INTEGER, INTENT(IN) :: k2dint    ! Horizontal interpolation type (see header)
1732      INTEGER, INTENT(IN) :: kdaystp   ! Number of time steps per day                   
1733      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &
1734         & pun,    &    ! Model zonal component of velocity
1735         & pvn,    &    ! Model meridional component of velocity
1736         & pumask, &    ! Land-sea mask
1737         & pvmask       ! Land-sea mask
1738      REAL(KIND=wp), INTENT(IN), DIMENSION(kpk) :: &
1739         & pgdept       ! Model array of depth levels
1740      LOGICAL, INTENT(IN) :: ld_dailyav
1741         
1742      !! * Local declarations
1743      INTEGER :: ji
1744      INTEGER :: jj
1745      INTEGER :: jk
1746      INTEGER :: jobs
1747      INTEGER :: inrc
1748      INTEGER :: ipro
1749      INTEGER :: idayend
1750      INTEGER :: ista
1751      INTEGER :: iend
1752      INTEGER :: iobs
1753      INTEGER, DIMENSION(imaxavtypes) :: &
1754         & idailyavtypes
1755      REAL(KIND=wp) :: zlam
1756      REAL(KIND=wp) :: zphi
1757      REAL(KIND=wp) :: zdaystp
1758      REAL(KIND=wp), DIMENSION(kpk) :: &
1759         & zobsmasku, &
1760         & zobsmaskv, &
1761         & zobsmask,  &
1762         & zobsk,     &
1763         & zobs2k
1764      REAL(KIND=wp), DIMENSION(2,2,kpk) :: &
1765         & zweigu,zweigv
1766      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
1767         & zumask, zvmask, &
1768         & zintu, &
1769         & zintv, &
1770         & zinmu, &
1771         & zinmv
1772      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
1773         & zglamu, zglamv, &
1774         & zgphiu, zgphiv
1775      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
1776         & igrdiu, &
1777         & igrdju, &
1778         & igrdiv, &
1779         & igrdjv
1780
1781      !------------------------------------------------------------------------
1782      ! Local initialization
1783      !------------------------------------------------------------------------
1784      ! ... Record and data counters
1785      inrc = kt - kit000 + 2
1786      ipro = prodatqc%npstp(inrc)
1787
1788      ! Initialize daily mean for first timestep
1789      idayend = MOD( kt - kit000 + 1, kdaystp )
1790
1791      ! Added kt == 0 test to catch restart case
1792      IF ( idayend == 1 .OR. kt == 0) THEN
1793         IF (lwp) WRITE(numout,*) 'Reset prodatqc%vdmean on time-step: ',kt
1794         prodatqc%vdmean(:,:,:,1) = 0.0
1795         prodatqc%vdmean(:,:,:,2) = 0.0
1796      ENDIF
1797
1798      ! Increment the zonal velocity field for computing daily mean
1799      prodatqc%vdmean(:,:,:,1) = prodatqc%vdmean(:,:,:,1) + pun(:,:,:)
1800      ! Increment the meridional velocity field for computing daily mean
1801      prodatqc%vdmean(:,:,:,2) = prodatqc%vdmean(:,:,:,2) + pvn(:,:,:)
1802   
1803      ! Compute the daily mean at the end of day
1804      zdaystp = 1.0 / REAL( kdaystp )
1805      IF ( idayend == 0 ) THEN
1806         prodatqc%vdmean(:,:,:,1) = prodatqc%vdmean(:,:,:,1) * zdaystp
1807         prodatqc%vdmean(:,:,:,2) = prodatqc%vdmean(:,:,:,2) * zdaystp
1808      ENDIF
1809
1810      ! Get the data for interpolation
1811      ALLOCATE( &
1812         & igrdiu(2,2,ipro),      &
1813         & igrdju(2,2,ipro),      &
1814         & igrdiv(2,2,ipro),      &
1815         & igrdjv(2,2,ipro),      &
1816         & zglamu(2,2,ipro), zglamv(2,2,ipro), &
1817         & zgphiu(2,2,ipro), zgphiv(2,2,ipro), &
1818         & zumask(2,2,kpk,ipro), zvmask(2,2,kpk,ipro), &
1819         & zintu(2,2,kpk,ipro),  &
1820         & zintv(2,2,kpk,ipro)   &
1821         & )
1822
1823      DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
1824         iobs = jobs - prodatqc%nprofup
1825         igrdiu(1,1,iobs) = prodatqc%mi(jobs,1)-1
1826         igrdju(1,1,iobs) = prodatqc%mj(jobs,1)-1
1827         igrdiu(1,2,iobs) = prodatqc%mi(jobs,1)-1
1828         igrdju(1,2,iobs) = prodatqc%mj(jobs,1)
1829         igrdiu(2,1,iobs) = prodatqc%mi(jobs,1)
1830         igrdju(2,1,iobs) = prodatqc%mj(jobs,1)-1
1831         igrdiu(2,2,iobs) = prodatqc%mi(jobs,1)
1832         igrdju(2,2,iobs) = prodatqc%mj(jobs,1)
1833         igrdiv(1,1,iobs) = prodatqc%mi(jobs,2)-1
1834         igrdjv(1,1,iobs) = prodatqc%mj(jobs,2)-1
1835         igrdiv(1,2,iobs) = prodatqc%mi(jobs,2)-1
1836         igrdjv(1,2,iobs) = prodatqc%mj(jobs,2)
1837         igrdiv(2,1,iobs) = prodatqc%mi(jobs,2)
1838         igrdjv(2,1,iobs) = prodatqc%mj(jobs,2)-1
1839         igrdiv(2,2,iobs) = prodatqc%mi(jobs,2)
1840         igrdjv(2,2,iobs) = prodatqc%mj(jobs,2)
1841      END DO
1842
1843      CALL obs_int_comm_2d( 2, 2, ipro, igrdiu, igrdju, glamu, zglamu )
1844      CALL obs_int_comm_2d( 2, 2, ipro, igrdiu, igrdju, gphiu, zgphiu )
1845      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, pumask, zumask )
1846      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, pun, zintu )
1847
1848      CALL obs_int_comm_2d( 2, 2, ipro, igrdiv, igrdjv, glamv, zglamv )
1849      CALL obs_int_comm_2d( 2, 2, ipro, igrdiv, igrdjv, gphiv, zgphiv )
1850      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, pvmask, zvmask )
1851      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, pvn, zintv )
1852
1853      ! At the end of the day also get interpolated means
1854      IF ( idayend == 0 ) THEN
1855
1856         ALLOCATE( &
1857            & zinmu(2,2,kpk,ipro),  &
1858            & zinmv(2,2,kpk,ipro)   &
1859            & )
1860
1861         CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, &
1862            &                  prodatqc%vdmean(:,:,:,1), zinmu )
1863         CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, &
1864            &                  prodatqc%vdmean(:,:,:,2), zinmv )
1865
1866      ENDIF
1867
1868! loop over observations
1869
1870      DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
1871
1872         iobs = jobs - prodatqc%nprofup
1873
1874         IF ( kt /= prodatqc%mstp(jobs) ) THEN
1875           
1876            IF(lwp) THEN
1877               WRITE(numout,*)
1878               WRITE(numout,*) ' E R R O R : Observation',              &
1879                  &            ' time step is not consistent with the', &
1880                  &            ' model time step'
1881               WRITE(numout,*) ' ========='
1882               WRITE(numout,*)
1883               WRITE(numout,*) ' Record  = ', jobs,                    &
1884                  &            ' kt      = ', kt,                      &
1885                  &            ' mstp    = ', prodatqc%mstp(jobs), &
1886                  &            ' ntyp    = ', prodatqc%ntyp(jobs)
1887            ENDIF
1888            CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' )
1889         ENDIF
1890         
1891         zlam = prodatqc%rlam(jobs)
1892         zphi = prodatqc%rphi(jobs)
1893
1894         ! Initialize observation masks
1895
1896         zobsmasku(:) = 0.0
1897         zobsmaskv(:) = 0.0
1898         
1899         ! Horizontal weights and vertical mask
1900
1901         IF  ( prodatqc%npvend(jobs,1) > 0 ) THEN
1902
1903            CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi,     &
1904               &                   zglamu(:,:,iobs), zgphiu(:,:,iobs), &
1905               &                   zumask(:,:,:,iobs), zweigu, zobsmasku )
1906
1907         ENDIF
1908
1909         
1910         IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
1911
1912            CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi,     &
1913               &                   zglamv(:,:,iobs), zgphiv(:,:,iobs), &
1914               &                   zvmask(:,:,:,iobs), zweigv, zobsmaskv )
1915
1916         ENDIF
1917
1918         ! Ensure that the vertical mask on u and v are consistent.
1919
1920         zobsmask(:) = MIN( zobsmasku(:), zobsmaskv(:) )
1921
1922         IF ( prodatqc%npvend(jobs,1) > 0 ) THEN
1923
1924            zobsk(:) = obfillflt
1925
1926       IF ( ld_dailyav ) THEN
1927
1928               IF ( idayend == 0 )  THEN
1929                 
1930                  ! Daily averaged data
1931                 
1932                  CALL obs_int_h2d( kpk, kpk,      &
1933                     &              zweigu, zinmu(:,:,:,iobs), zobsk )
1934                 
1935                 
1936               ELSE
1937               
1938                  CALL ctl_stop( ' A nonzero' //     &
1939                     &           ' number of U profile data should' // &
1940                     &           ' only occur at the end of a given day' )
1941
1942               ENDIF
1943         
1944            ELSE 
1945               
1946               ! Point data
1947
1948               CALL obs_int_h2d( kpk, kpk,      &
1949                  &              zweigu, zintu(:,:,:,iobs), zobsk )
1950
1951            ENDIF
1952
1953            !-------------------------------------------------------------
1954            ! Compute vertical second-derivative of the interpolating
1955            ! polynomial at obs points
1956            !-------------------------------------------------------------
1957           
1958            IF ( k1dint == 1 ) THEN
1959               CALL obs_int_z1d_spl( kpk, zobsk, zobs2k,   &
1960                  &                  pgdept, zobsmask )
1961            ENDIF
1962           
1963            !-----------------------------------------------------------------
1964            !  Vertical interpolation to the observation point
1965            !-----------------------------------------------------------------
1966            ista = prodatqc%npvsta(jobs,1)
1967            iend = prodatqc%npvend(jobs,1)
1968            CALL obs_int_z1d( kpk,                &
1969               & prodatqc%var(1)%mvk(ista:iend),  &
1970               & k1dint, iend - ista + 1,         &
1971               & prodatqc%var(1)%vdep(ista:iend), &
1972               & zobsk, zobs2k,                   &
1973               & prodatqc%var(1)%vmod(ista:iend), &
1974               & pgdept, zobsmask )
1975
1976         ENDIF
1977
1978         IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
1979
1980            zobsk(:) = obfillflt
1981
1982            IF ( ld_dailyav ) THEN
1983
1984               IF ( idayend == 0 )  THEN
1985
1986                  ! Daily averaged data
1987                 
1988                  CALL obs_int_h2d( kpk, kpk,      &
1989                     &              zweigv, zinmv(:,:,:,iobs), zobsk )
1990                 
1991               ELSE
1992
1993                  CALL ctl_stop( ' A nonzero' //     &
1994                     &           ' number of V profile data should' // &
1995                     &           ' only occur at the end of a given day' )
1996
1997               ENDIF
1998
1999            ELSE
2000               
2001               ! Point data
2002
2003               CALL obs_int_h2d( kpk, kpk,      &
2004                  &              zweigv, zintv(:,:,:,iobs), zobsk )
2005
2006            ENDIF
2007
2008
2009            !-------------------------------------------------------------
2010            ! Compute vertical second-derivative of the interpolating
2011            ! polynomial at obs points
2012            !-------------------------------------------------------------
2013           
2014            IF ( k1dint == 1 ) THEN
2015               CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, &
2016                  &                  pgdept, zobsmask )
2017            ENDIF
2018           
2019            !----------------------------------------------------------------
2020            !  Vertical interpolation to the observation point
2021            !----------------------------------------------------------------
2022            ista = prodatqc%npvsta(jobs,2)
2023            iend = prodatqc%npvend(jobs,2)
2024            CALL obs_int_z1d( kpk, &
2025               & prodatqc%var(2)%mvk(ista:iend),&
2026               & k1dint, iend - ista + 1, &
2027               & prodatqc%var(2)%vdep(ista:iend),&
2028               & zobsk, zobs2k, &
2029               & prodatqc%var(2)%vmod(ista:iend),&
2030               & pgdept, zobsmask )
2031
2032         ENDIF
2033
2034      END DO
2035 
2036      ! Deallocate the data for interpolation
2037      DEALLOCATE( &
2038         & igrdiu, &
2039         & igrdju, &
2040         & igrdiv, &
2041         & igrdjv, &
2042         & zglamu, zglamv, &
2043         & zgphiu, zgphiv, &
2044         & zumask, zvmask, &
2045         & zintu, &
2046         & zintv  &
2047         & )
2048      ! At the end of the day also get interpolated means
2049      IF ( idayend == 0 ) THEN
2050         DEALLOCATE( &
2051            & zinmu,  &
2052            & zinmv   &
2053            & )
2054      ENDIF
2055
2056      prodatqc%nprofup = prodatqc%nprofup + ipro 
2057     
2058   END SUBROUTINE obs_vel_opt
2059
2060   SUBROUTINE obs_logchl_opt( logchldatqc, kt, kpi, kpj, kit000, &
2061      &                    plogchln, plogchlmask, k2dint )
2062
2063      !!-----------------------------------------------------------------------
2064      !!
2065      !!                     ***  ROUTINE obs_logchl_opt  ***
2066      !!
2067      !! ** Purpose : Compute the model counterpart of logchl
2068      !!              data by interpolating from the model grid to the
2069      !!              observation point.
2070      !!
2071      !! ** Method  : Linearly interpolate to each observation point using
2072      !!              the model values at the corners of the surrounding grid box.
2073      !!
2074      !!    The now model logchl is first computed at the obs (lon, lat) point.
2075      !!
2076      !!    Several horizontal interpolation schemes are available:
2077      !!        - distance-weighted (great circle) (k2dint = 0)
2078      !!        - distance-weighted (small angle)  (k2dint = 1)
2079      !!        - bilinear (geographical grid)     (k2dint = 2)
2080      !!        - bilinear (quadrilateral grid)    (k2dint = 3)
2081      !!        - polynomial (quadrilateral grid)  (k2dint = 4)
2082      !!
2083      !!
2084      !! ** Action  :
2085      !!
2086      !! History :
2087      !!     
2088      !!-----------------------------------------------------------------------
2089
2090      !! * Modules used
2091      USE obs_surf_def  ! Definition of storage space for surface observations
2092
2093      IMPLICIT NONE
2094
2095      !! * Arguments
2096      TYPE(obs_surf), INTENT(INOUT) :: logchldatqc     ! Subset of surface data not failing screening
2097      INTEGER, INTENT(IN) :: kt       ! Time step
2098      INTEGER, INTENT(IN) :: kpi      ! Model grid parameters
2099      INTEGER, INTENT(IN) :: kpj
2100      INTEGER, INTENT(IN) :: kit000   ! Number of the first time step
2101                                      !   (kit000-1 = restart time)
2102      INTEGER, INTENT(IN) :: k2dint   ! Horizontal interpolation type (see header)
2103      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
2104         & plogchln,  &    ! Model logchl field
2105         & plogchlmask     ! Land-sea mask
2106         
2107      !! * Local declarations
2108      INTEGER :: ji
2109      INTEGER :: jj
2110      INTEGER :: jobs
2111      INTEGER :: inrc
2112      INTEGER :: ilogchl
2113      INTEGER :: iobs
2114       
2115      REAL(KIND=wp) :: zlam
2116      REAL(KIND=wp) :: zphi
2117      REAL(KIND=wp) :: zext(1), zobsmask(1)
2118      REAL(kind=wp), DIMENSION(2,2,1) :: &
2119         & zweig
2120      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
2121         & zmask, &
2122         & zlogchll, &
2123         & zglam, &
2124         & zgphi
2125      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
2126         & igrdi, &
2127         & igrdj
2128
2129      !------------------------------------------------------------------------
2130      ! Local initialization
2131      !------------------------------------------------------------------------
2132      ! ... Record and data counters
2133      inrc = kt - kit000 + 2
2134      ilogchl = logchldatqc%nsstp(inrc)
2135
2136      ! Get the data for interpolation
2137     
2138      ALLOCATE( &
2139         & igrdi(2,2,ilogchl), &
2140         & igrdj(2,2,ilogchl), &
2141         & zglam(2,2,ilogchl), &
2142         & zgphi(2,2,ilogchl), &
2143         & zmask(2,2,ilogchl), &
2144         & zlogchll(2,2,ilogchl)  &
2145         & )
2146     
2147      DO jobs = logchldatqc%nsurfup + 1, logchldatqc%nsurfup + ilogchl
2148         iobs = jobs - logchldatqc%nsurfup
2149         igrdi(1,1,iobs) = logchldatqc%mi(jobs)-1
2150         igrdj(1,1,iobs) = logchldatqc%mj(jobs)-1
2151         igrdi(1,2,iobs) = logchldatqc%mi(jobs)-1
2152         igrdj(1,2,iobs) = logchldatqc%mj(jobs)
2153         igrdi(2,1,iobs) = logchldatqc%mi(jobs)
2154         igrdj(2,1,iobs) = logchldatqc%mj(jobs)-1
2155         igrdi(2,2,iobs) = logchldatqc%mi(jobs)
2156         igrdj(2,2,iobs) = logchldatqc%mj(jobs)
2157      END DO
2158     
2159      CALL obs_int_comm_2d( 2, 2, ilogchl, &
2160         &                  igrdi, igrdj, glamt, zglam )
2161      CALL obs_int_comm_2d( 2, 2, ilogchl, &
2162         &                  igrdi, igrdj, gphit, zgphi )
2163      CALL obs_int_comm_2d( 2, 2, ilogchl, &
2164         &                  igrdi, igrdj, plogchlmask, zmask )
2165      CALL obs_int_comm_2d( 2, 2, ilogchl, &
2166         &                  igrdi, igrdj, plogchln, zlogchll )
2167     
2168      DO jobs = logchldatqc%nsurfup + 1, logchldatqc%nsurfup + ilogchl
2169         
2170         iobs = jobs - logchldatqc%nsurfup
2171         
2172         IF ( kt /= logchldatqc%mstp(jobs) ) THEN
2173           
2174            IF(lwp) THEN
2175               WRITE(numout,*)
2176               WRITE(numout,*) ' E R R O R : Observation',              &
2177                  &            ' time step is not consistent with the', &
2178                  &            ' model time step'
2179               WRITE(numout,*) ' ========='
2180               WRITE(numout,*)
2181               WRITE(numout,*) ' Record  = ', jobs,                &
2182                  &            ' kt      = ', kt,                  &
2183                  &            ' mstp    = ', logchldatqc%mstp(jobs), &
2184                  &            ' ntyp    = ', logchldatqc%ntyp(jobs)
2185            ENDIF
2186            CALL ctl_stop( 'obs_logchl_opt', 'Inconsistent time' )
2187           
2188         ENDIF
2189         
2190         zlam = logchldatqc%rlam(jobs)
2191         zphi = logchldatqc%rphi(jobs)
2192         
2193         ! Get weights to interpolate the model logchl to the observation point
2194         CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         &
2195            &                   zglam(:,:,iobs), zgphi(:,:,iobs), &
2196            &                   zmask(:,:,iobs), zweig, zobsmask )
2197         
2198         ! ... Interpolate the model logchl to the observation point
2199         CALL obs_int_h2d( 1, 1,      &
2200            &              zweig, zlogchll(:,:,iobs),  zext )
2201         
2202         logchldatqc%rmod(jobs,1) = zext(1)
2203         
2204      END DO
2205     
2206      ! Deallocate the data for interpolation
2207      DEALLOCATE( &
2208         & igrdi,    &
2209         & igrdj,    &
2210         & zglam,    &
2211         & zgphi,    &
2212         & zmask,    &
2213         & zlogchll  &
2214         & )
2215     
2216      logchldatqc%nsurfup = logchldatqc%nsurfup + ilogchl
2217
2218   END SUBROUTINE obs_logchl_opt
2219
2220   SUBROUTINE obs_spm_opt( spmdatqc, kt, kpi, kpj, kit000, &
2221      &                    pspmn, pspmmask, k2dint )
2222
2223      !!-----------------------------------------------------------------------
2224      !!
2225      !!                     ***  ROUTINE obs_spm_opt  ***
2226      !!
2227      !! ** Purpose : Compute the model counterpart of spm
2228      !!              data by interpolating from the model grid to the
2229      !!              observation point.
2230      !!
2231      !! ** Method  : Linearly interpolate to each observation point using
2232      !!              the model values at the corners of the surrounding grid box.
2233      !!
2234      !!    The now model spm is first computed at the obs (lon, lat) point.
2235      !!
2236      !!    Several horizontal interpolation schemes are available:
2237      !!        - distance-weighted (great circle) (k2dint = 0)
2238      !!        - distance-weighted (small angle)  (k2dint = 1)
2239      !!        - bilinear (geographical grid)     (k2dint = 2)
2240      !!        - bilinear (quadrilateral grid)    (k2dint = 3)
2241      !!        - polynomial (quadrilateral grid)  (k2dint = 4)
2242      !!
2243      !!
2244      !! ** Action  :
2245      !!
2246      !! History :
2247      !!     
2248      !!-----------------------------------------------------------------------
2249
2250      !! * Modules used
2251      USE obs_surf_def  ! Definition of storage space for surface observations
2252
2253      IMPLICIT NONE
2254
2255      !! * Arguments
2256      TYPE(obs_surf), INTENT(INOUT) :: spmdatqc     ! Subset of surface data not failing screening
2257      INTEGER, INTENT(IN) :: kt       ! Time step
2258      INTEGER, INTENT(IN) :: kpi      ! Model grid parameters
2259      INTEGER, INTENT(IN) :: kpj
2260      INTEGER, INTENT(IN) :: kit000   ! Number of the first time step
2261                                      !   (kit000-1 = restart time)
2262      INTEGER, INTENT(IN) :: k2dint   ! Horizontal interpolation type (see header)
2263      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
2264         & pspmn,  &    ! Model spm field
2265         & pspmmask     ! Land-sea mask
2266         
2267      !! * Local declarations
2268      INTEGER :: ji
2269      INTEGER :: jj
2270      INTEGER :: jobs
2271      INTEGER :: inrc
2272      INTEGER :: ispm
2273      INTEGER :: iobs
2274       
2275      REAL(KIND=wp) :: zlam
2276      REAL(KIND=wp) :: zphi
2277      REAL(KIND=wp) :: zext(1), zobsmask(1)
2278      REAL(kind=wp), DIMENSION(2,2,1) :: &
2279         & zweig
2280      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
2281         & zmask, &
2282         & zspml, &
2283         & zglam, &
2284         & zgphi
2285      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
2286         & igrdi, &
2287         & igrdj
2288
2289      !------------------------------------------------------------------------
2290      ! Local initialization
2291      !------------------------------------------------------------------------
2292      ! ... Record and data counters
2293      inrc = kt - kit000 + 2
2294      ispm = spmdatqc%nsstp(inrc)
2295
2296      ! Get the data for interpolation
2297     
2298      ALLOCATE( &
2299         & igrdi(2,2,ispm), &
2300         & igrdj(2,2,ispm), &
2301         & zglam(2,2,ispm), &
2302         & zgphi(2,2,ispm), &
2303         & zmask(2,2,ispm), &
2304         & zspml(2,2,ispm)  &
2305         & )
2306     
2307      DO jobs = spmdatqc%nsurfup + 1, spmdatqc%nsurfup + ispm
2308         iobs = jobs - spmdatqc%nsurfup
2309         igrdi(1,1,iobs) = spmdatqc%mi(jobs)-1
2310         igrdj(1,1,iobs) = spmdatqc%mj(jobs)-1
2311         igrdi(1,2,iobs) = spmdatqc%mi(jobs)-1
2312         igrdj(1,2,iobs) = spmdatqc%mj(jobs)
2313         igrdi(2,1,iobs) = spmdatqc%mi(jobs)
2314         igrdj(2,1,iobs) = spmdatqc%mj(jobs)-1
2315         igrdi(2,2,iobs) = spmdatqc%mi(jobs)
2316         igrdj(2,2,iobs) = spmdatqc%mj(jobs)
2317      END DO
2318     
2319      CALL obs_int_comm_2d( 2, 2, ispm, &
2320         &                  igrdi, igrdj, glamt, zglam )
2321      CALL obs_int_comm_2d( 2, 2, ispm, &
2322         &                  igrdi, igrdj, gphit, zgphi )
2323      CALL obs_int_comm_2d( 2, 2, ispm, &
2324         &                  igrdi, igrdj, pspmmask, zmask )
2325      CALL obs_int_comm_2d( 2, 2, ispm, &
2326         &                  igrdi, igrdj, pspmn, zspml )
2327     
2328      DO jobs = spmdatqc%nsurfup + 1, spmdatqc%nsurfup + ispm
2329         
2330         iobs = jobs - spmdatqc%nsurfup
2331         
2332         IF ( kt /= spmdatqc%mstp(jobs) ) THEN
2333           
2334            IF(lwp) THEN
2335               WRITE(numout,*)
2336               WRITE(numout,*) ' E R R O R : Observation',              &
2337                  &            ' time step is not consistent with the', &
2338                  &            ' model time step'
2339               WRITE(numout,*) ' ========='
2340               WRITE(numout,*)
2341               WRITE(numout,*) ' Record  = ', jobs,                &
2342                  &            ' kt      = ', kt,                  &
2343                  &            ' mstp    = ', spmdatqc%mstp(jobs), &
2344                  &            ' ntyp    = ', spmdatqc%ntyp(jobs)
2345            ENDIF
2346            CALL ctl_stop( 'obs_spm_opt', 'Inconsistent time' )
2347           
2348         ENDIF
2349         
2350         zlam = spmdatqc%rlam(jobs)
2351         zphi = spmdatqc%rphi(jobs)
2352         
2353         ! Get weights to interpolate the model spm to the observation point
2354         CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         &
2355            &                   zglam(:,:,iobs), zgphi(:,:,iobs), &
2356            &                   zmask(:,:,iobs), zweig, zobsmask )
2357         
2358         ! ... Interpolate the model spm to the observation point
2359         CALL obs_int_h2d( 1, 1,      &
2360            &              zweig, zspml(:,:,iobs),  zext )
2361         
2362         spmdatqc%rmod(jobs,1) = zext(1)
2363         
2364      END DO
2365     
2366      ! Deallocate the data for interpolation
2367      DEALLOCATE( &
2368         & igrdi,    &
2369         & igrdj,    &
2370         & zglam,    &
2371         & zgphi,    &
2372         & zmask,    &
2373         & zspml  &
2374         & )
2375     
2376      spmdatqc%nsurfup = spmdatqc%nsurfup + ispm
2377
2378   END SUBROUTINE obs_spm_opt
2379
2380END MODULE obs_oper
2381
Note: See TracBrowser for help on using the repository browser.