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

source: branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90 @ 6301

Last change on this file since 6301 was 6301, checked in by djlea, 6 years ago

Update to QC out observations which are caught by the strict bathymetry check in obs_oper. Reverting some changes to obs_prep which are no longer needed.

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