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

Last change on this file since 6016 was 6016, checked in by kingr, 5 years ago

#1642 bug fixes for general vertical coord obsoper

File size: 75.5 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               ENDDO 
888 
889 
890               DEALLOCATE(interp_corner,iv_indic) 
891         
892            ENDIF 
893       
894 
895            ! Salinity 
896         
897            IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 
898   
899               zobsk(:) = obfillflt 
900   
901               IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN
902   
903                  IF ( idayend == 0 )  THEN 
904                   
905                     ! Daily averaged moored buoy (MRB) data
906                   
907                     ! vertically interpolate all 4 corners
908                     ista = prodatqc%npvsta(jobs,2) 
909                     iend = prodatqc%npvend(jobs,2) 
910                     inum_obs = iend - ista + 1 
911                     ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs)) 
912     
913                     DO iin=1,2 
914                        DO ijn=1,2 
915                                       
916                                       
917           
918                           IF ( k1dint == 1 ) THEN
919                              CALL obs_int_z1d_spl( kpk, & 
920                                 &     zinms(iin,ijn,:,iobs), & 
921                                 &     zobs2k, zgdept(iin,ijn,:,iobs), & 
922                                 &     zmask(iin,ijn,:,iobs)) 
923                           ENDIF
924       
925                           CALL obs_level_search(kpk, & 
926                              &    zgdept(iin,ijn,:,iobs), & 
927                              &    inum_obs, prodatqc%var(2)%vdep(ista:iend), & 
928                              &    iv_indic) 
929                           CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 
930                              &    prodatqc%var(2)%vdep(ista:iend), & 
931                              &    zinms(iin,ijn,:,iobs), & 
932                              &    zobs2k, interp_corner(iin,ijn,:), & 
933                              &    zgdept(iin,ijn,:,iobs), & 
934                              &    zmask(iin,ijn,:,iobs)) 
935       
936                        ENDDO 
937                     ENDDO 
938                   
939                   
940                  ELSE
941               
942                     CALL ctl_stop( ' A nonzero' //     & 
943                        &           ' number of profile T BUOY data should' // & 
944                        &           ' only occur at the end of a given day' ) 
945   
946                  ENDIF
947         
948               ELSE 
949               
950                  ! Point data
951     
952                  ! vertically interpolate all 4 corners
953                  ista = prodatqc%npvsta(jobs,2) 
954                  iend = prodatqc%npvend(jobs,2) 
955                  inum_obs = iend - ista + 1 
956                  ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs)) 
957                   
958                  DO iin=1,2     
959                     DO ijn=1,2 
960                                 
961                                 
962                        IF ( k1dint == 1 ) THEN
963                           CALL obs_int_z1d_spl( kpk, & 
964                              &    zints(iin,ijn,:,iobs),& 
965                              &    zobs2k, zgdept(iin,ijn,:,iobs), & 
966                              &    zmask(iin,ijn,:,iobs)) 
967 
968                        ENDIF
969       
970                        CALL obs_level_search(kpk, & 
971                           &        zgdept(iin,ijn,:,iobs),& 
972                           &        inum_obs, prodatqc%var(2)%vdep(ista:iend), & 
973                           &         iv_indic) 
974                        CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs,  & 
975                           &          prodatqc%var(2)%vdep(ista:iend),     & 
976                           &          zints(iin,ijn,:,iobs),               & 
977                           &          zobs2k,interp_corner(iin,ijn,:),     & 
978                           &          zgdept(iin,ijn,:,iobs),              & 
979                           &          zmask(iin,ijn,:,iobs) )     
980         
981                     ENDDO 
982                  ENDDO 
983             
984               ENDIF 
985       
986               !-------------------------------------------------------------
987               ! Compute the horizontal interpolation for every profile level
988               !-------------------------------------------------------------
989             
990               DO ikn=1,inum_obs 
991                  iend=ista+ikn-1
992                 
993                  l_zweig(:,:,1) = 0._wp 
994   
995                  ! This code forces the horizontal weights to be 
996                  ! zero IF the observation is below the bottom of the 
997                  ! corners of the interpolation nodes, Or if it is in 
998                  ! the mask. This is important for observations are near 
999                  ! steep bathymetry
1000                  DO iin=1,2 
1001                     DO ijn=1,2 
1002     
1003                        depth_loop2: DO ik=kpk,2,-1 
1004                           IF(zmask(iin,ijn,ik-1,iobs ) > 0.9 )THEN   
1005                           
1006                              l_zweig(iin,ijn,1) = & 
1007                                 &  zweig(iin,ijn,1) * & 
1008                                 &  MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) & 
1009                                 &  - prodatqc%var(2)%vdep(iend)),0._wp) 
1010                           
1011                              EXIT depth_loop2 
1012                           ENDIF
1013                        ENDDO depth_loop2 
1014     
1015                     ENDDO 
1016                  ENDDO 
1017   
1018                  CALL obs_int_h2d( 1, 1, l_zweig, interp_corner(:,:,ikn), & 
1019                  &          prodatqc%var(2)%vmod(iend:iend) ) 
1020 
1021               ENDDO 
1022 
1023 
1024               DEALLOCATE(interp_corner,iv_indic) 
1025         
1026            ENDIF
1027         
1028         ENDIF
1029       
1030      END DO 
1031     
1032      ! Deallocate the data for interpolation
1033      DEALLOCATE( & 
1034         & igrdi, & 
1035         & igrdj, & 
1036         & zglam, & 
1037         & zgphi, & 
1038         & zmask, & 
1039         & zintt, & 
1040         & zints, & 
1041         & zgdept,& 
1042         & zgdepw & 
1043         & ) 
1044      ! At the end of the day also get interpolated means
1045      IF ( idayend == 0 ) THEN
1046         DEALLOCATE( & 
1047            & zinmt,  & 
1048            & zinms   & 
1049            & ) 
1050      ENDIF
1051   
1052      prodatqc%nprofup = prodatqc%nprofup + ipro 
1053       
1054   END SUBROUTINE obs_pro_sco_opt 
1055 
1056   SUBROUTINE obs_sla_opt( sladatqc, kt, kpi, kpj, kit000, &
1057      &                    psshn, psshmask, k2dint )
1058      !!-----------------------------------------------------------------------
1059      !!
1060      !!                     ***  ROUTINE obs_sla_opt  ***
1061      !!
1062      !! ** Purpose : Compute the model counterpart of sea level anomaly
1063      !!              data by interpolating from the model grid to the
1064      !!              observation point.
1065      !!
1066      !! ** Method  : Linearly interpolate to each observation point using
1067      !!              the model values at the corners of the surrounding grid box.
1068      !!
1069      !!    The now model SSH is first computed at the obs (lon, lat) point.
1070      !!
1071      !!    Several horizontal interpolation schemes are available:
1072      !!        - distance-weighted (great circle) (k2dint = 0)
1073      !!        - distance-weighted (small angle)  (k2dint = 1)
1074      !!        - bilinear (geographical grid)     (k2dint = 2)
1075      !!        - bilinear (quadrilateral grid)    (k2dint = 3)
1076      !!        - polynomial (quadrilateral grid)  (k2dint = 4)
1077      !! 
1078      !!    The sea level anomaly at the observation points is then computed
1079      !!    by removing a mean dynamic topography (defined at the obs. point).
1080      !!
1081      !! ** Action  :
1082      !!
1083      !! History :
1084      !!      ! 07-03 (A. Weaver)
1085      !!-----------------------------------------------------------------------
1086 
1087      !! * Modules used
1088      USE obs_surf_def  ! Definition of storage space for surface observations
1089
1090      IMPLICIT NONE
1091
1092      !! * Arguments
1093      TYPE(obs_surf), INTENT(INOUT) :: sladatqc     ! Subset of surface data not failing screening
1094      INTEGER, INTENT(IN) :: kt      ! Time step
1095      INTEGER, INTENT(IN) :: kpi     ! Model grid parameters
1096      INTEGER, INTENT(IN) :: kpj
1097      INTEGER, INTENT(IN) :: kit000   ! Number of the first time step
1098                                      !   (kit000-1 = restart time)
1099      INTEGER, INTENT(IN) :: k2dint   ! Horizontal interpolation type (see header)
1100      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
1101         & psshn,  &    ! Model SSH field
1102         & psshmask     ! Land-sea mask
1103         
1104      !! * Local declarations
1105      INTEGER :: ji
1106      INTEGER :: jj
1107      INTEGER :: jobs
1108      INTEGER :: inrc
1109      INTEGER :: isla
1110      INTEGER :: iobs
1111      REAL(KIND=wp) :: zlam
1112      REAL(KIND=wp) :: zphi
1113      REAL(KIND=wp) :: zext(1), zobsmask(1)
1114      REAL(kind=wp), DIMENSION(2,2,1) :: &
1115         & zweig
1116      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
1117         & zmask, &
1118         & zsshl, &
1119         & zglam, &
1120         & zgphi
1121      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
1122         & igrdi, &
1123         & igrdj
1124
1125      !------------------------------------------------------------------------
1126      ! Local initialization
1127      !------------------------------------------------------------------------
1128      ! ... Record and data counters
1129      inrc = kt - kit000 + 2
1130      isla = sladatqc%nsstp(inrc)
1131
1132      ! Get the data for interpolation
1133
1134      ALLOCATE( &
1135         & igrdi(2,2,isla), &
1136         & igrdj(2,2,isla), &
1137         & zglam(2,2,isla), &
1138         & zgphi(2,2,isla), &
1139         & zmask(2,2,isla), &
1140         & zsshl(2,2,isla)  &
1141         & )
1142     
1143      DO jobs = sladatqc%nsurfup + 1, sladatqc%nsurfup + isla
1144         iobs = jobs - sladatqc%nsurfup
1145         igrdi(1,1,iobs) = sladatqc%mi(jobs)-1
1146         igrdj(1,1,iobs) = sladatqc%mj(jobs)-1
1147         igrdi(1,2,iobs) = sladatqc%mi(jobs)-1
1148         igrdj(1,2,iobs) = sladatqc%mj(jobs)
1149         igrdi(2,1,iobs) = sladatqc%mi(jobs)
1150         igrdj(2,1,iobs) = sladatqc%mj(jobs)-1
1151         igrdi(2,2,iobs) = sladatqc%mi(jobs)
1152         igrdj(2,2,iobs) = sladatqc%mj(jobs)
1153      END DO
1154
1155      CALL obs_int_comm_2d( 2, 2, isla, &
1156         &                  igrdi, igrdj, glamt, zglam )
1157      CALL obs_int_comm_2d( 2, 2, isla, &
1158         &                  igrdi, igrdj, gphit, zgphi )
1159      CALL obs_int_comm_2d( 2, 2, isla, &
1160         &                  igrdi, igrdj, psshmask, zmask )
1161      CALL obs_int_comm_2d( 2, 2, isla, &
1162         &                  igrdi, igrdj, psshn, zsshl )
1163
1164      ! Loop over observations
1165
1166      DO jobs = sladatqc%nsurfup + 1, sladatqc%nsurfup + isla
1167
1168         iobs = jobs - sladatqc%nsurfup
1169
1170         IF ( kt /= sladatqc%mstp(jobs) ) THEN
1171           
1172            IF(lwp) THEN
1173               WRITE(numout,*)
1174               WRITE(numout,*) ' E R R O R : Observation',              &
1175                  &            ' time step is not consistent with the', &
1176                  &            ' model time step'
1177               WRITE(numout,*) ' ========='
1178               WRITE(numout,*)
1179               WRITE(numout,*) ' Record  = ', jobs,                &
1180                  &            ' kt      = ', kt,                  &
1181                  &            ' mstp    = ', sladatqc%mstp(jobs), &
1182                  &            ' ntyp    = ', sladatqc%ntyp(jobs)
1183            ENDIF
1184            CALL ctl_stop( 'obs_sla_opt', 'Inconsistent time' )
1185           
1186         ENDIF
1187         
1188         zlam = sladatqc%rlam(jobs)
1189         zphi = sladatqc%rphi(jobs)
1190
1191         ! Get weights to interpolate the model SSH to the observation point
1192         CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         &
1193            &                   zglam(:,:,iobs), zgphi(:,:,iobs), &
1194            &                   zmask(:,:,iobs), zweig, zobsmask )
1195         
1196
1197         ! Interpolate the model SSH to the observation point
1198         CALL obs_int_h2d( 1, 1,      &
1199            &              zweig, zsshl(:,:,iobs),  zext )
1200         
1201         sladatqc%rext(jobs,1) = zext(1)
1202         ! ... Remove the MDT at the observation point
1203         sladatqc%rmod(jobs,1) = sladatqc%rext(jobs,1) - sladatqc%rext(jobs,2)
1204
1205      END DO
1206
1207      ! Deallocate the data for interpolation
1208      DEALLOCATE( &
1209         & igrdi, &
1210         & igrdj, &
1211         & zglam, &
1212         & zgphi, &
1213         & zmask, &
1214         & zsshl  &
1215         & )
1216
1217      sladatqc%nsurfup = sladatqc%nsurfup + isla
1218
1219   END SUBROUTINE obs_sla_opt
1220
1221   SUBROUTINE obs_sst_opt( sstdatqc, kt, kpi, kpj, kit000, kdaystp, &
1222      &                    psstn, psstmask, k2dint, ld_nightav )
1223      !!-----------------------------------------------------------------------
1224      !!
1225      !!                     ***  ROUTINE obs_sst_opt  ***
1226      !!
1227      !! ** Purpose : Compute the model counterpart of surface temperature
1228      !!              data by interpolating from the model grid to the
1229      !!              observation point.
1230      !!
1231      !! ** Method  : Linearly interpolate to each observation point using
1232      !!              the model values at the corners of the surrounding grid box.
1233      !!
1234      !!    The now model SST is first computed at the obs (lon, lat) point.
1235      !!
1236      !!    Several horizontal interpolation schemes are available:
1237      !!        - distance-weighted (great circle) (k2dint = 0)
1238      !!        - distance-weighted (small angle)  (k2dint = 1)
1239      !!        - bilinear (geographical grid)     (k2dint = 2)
1240      !!        - bilinear (quadrilateral grid)    (k2dint = 3)
1241      !!        - polynomial (quadrilateral grid)  (k2dint = 4)
1242      !!
1243      !!
1244      !! ** Action  :
1245      !!
1246      !! History :
1247      !!        !  07-07  (S. Ricci ) : Original
1248      !!     
1249      !!-----------------------------------------------------------------------
1250
1251      !! * Modules used
1252      USE obs_surf_def  ! Definition of storage space for surface observations
1253      USE sbcdcy
1254
1255      IMPLICIT NONE
1256
1257      !! * Arguments
1258      TYPE(obs_surf), INTENT(INOUT) :: &
1259         & sstdatqc     ! Subset of surface data not failing screening
1260      INTEGER, INTENT(IN) :: kt        ! Time step
1261      INTEGER, INTENT(IN) :: kpi       ! Model grid parameters
1262      INTEGER, INTENT(IN) :: kpj
1263      INTEGER, INTENT(IN) :: kit000    ! Number of the first time step
1264                                       !   (kit000-1 = restart time)
1265      INTEGER, INTENT(IN) :: k2dint    ! Horizontal interpolation type (see header)
1266      INTEGER, INTENT(IN) :: kdaystp   ! Number of time steps per day 
1267      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
1268         & psstn,  &    ! Model SST field
1269         & psstmask     ! Land-sea mask
1270
1271      !! * Local declarations
1272      INTEGER :: ji
1273      INTEGER :: jj
1274      INTEGER :: jobs
1275      INTEGER :: inrc
1276      INTEGER :: isst
1277      INTEGER :: iobs
1278      INTEGER :: idayend
1279      REAL(KIND=wp) :: zlam
1280      REAL(KIND=wp) :: zphi
1281      REAL(KIND=wp) :: zext(1), zobsmask(1)
1282      REAL(KIND=wp) :: zdaystp
1283      INTEGER, DIMENSION(:,:), SAVE, ALLOCATABLE :: &
1284         & icount_sstnight,      &
1285         & imask_night
1286      REAL(kind=wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: &
1287         & zintmp, &
1288         & zouttmp, & 
1289         & zmeanday    ! to compute model sst in region of 24h daylight (pole)
1290      REAL(kind=wp), DIMENSION(2,2,1) :: &
1291         & zweig
1292      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
1293         & zmask, &
1294         & zsstl, &
1295         & zsstm, &
1296         & zglam, &
1297         & zgphi
1298      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
1299         & igrdi, &
1300         & igrdj
1301      LOGICAL, INTENT(IN) :: ld_nightav
1302
1303      !-----------------------------------------------------------------------
1304      ! Local initialization
1305      !-----------------------------------------------------------------------
1306      ! ... Record and data counters
1307      inrc = kt - kit000 + 2
1308      isst = sstdatqc%nsstp(inrc)
1309
1310      IF ( ld_nightav ) THEN
1311
1312      ! Initialize array for night mean
1313
1314      IF ( kt .EQ. 0 ) THEN
1315         ALLOCATE ( icount_sstnight(kpi,kpj) )
1316         ALLOCATE ( imask_night(kpi,kpj) )
1317         ALLOCATE ( zintmp(kpi,kpj) )
1318         ALLOCATE ( zouttmp(kpi,kpj) )
1319         ALLOCATE ( zmeanday(kpi,kpj) )
1320         nday_qsr = -1   ! initialisation flag for nbc_dcy
1321      ENDIF
1322
1323      ! Initialize daily mean for first timestep
1324      idayend = MOD( kt - kit000 + 1, kdaystp )
1325
1326      ! Added kt == 0 test to catch restart case
1327      IF ( idayend == 1 .OR. kt == 0) THEN
1328         IF (lwp) WRITE(numout,*) 'Reset sstdatqc%vdmean on time-step: ',kt
1329         DO jj = 1, jpj
1330            DO ji = 1, jpi
1331               sstdatqc%vdmean(ji,jj) = 0.0
1332               zmeanday(ji,jj) = 0.0
1333               icount_sstnight(ji,jj) = 0
1334            END DO
1335         END DO
1336      ENDIF
1337
1338      zintmp(:,:) = 0.0
1339      zouttmp(:,:) = sbc_dcy( zintmp(:,:), .TRUE. )
1340      imask_night(:,:) = INT( zouttmp(:,:) )
1341
1342      DO jj = 1, jpj
1343         DO ji = 1, jpi
1344            ! Increment the temperature field for computing night mean and counter
1345            sstdatqc%vdmean(ji,jj) = sstdatqc%vdmean(ji,jj)  &
1346                   &                        + psstn(ji,jj)*imask_night(ji,jj)
1347            zmeanday(ji,jj)        = zmeanday(ji,jj) + psstn(ji,jj)
1348            icount_sstnight(ji,jj) = icount_sstnight(ji,jj) + imask_night(ji,jj)
1349         END DO
1350      END DO
1351   
1352      ! Compute the daily mean at the end of day
1353
1354      zdaystp = 1.0 / REAL( kdaystp )
1355
1356      IF ( idayend == 0 ) THEN
1357         DO jj = 1, jpj
1358            DO ji = 1, jpi
1359               ! Test if "no night" point
1360               IF ( icount_sstnight(ji,jj) .NE. 0 ) THEN
1361                  sstdatqc%vdmean(ji,jj) = sstdatqc%vdmean(ji,jj) &
1362                    &                        / icount_sstnight(ji,jj) 
1363               ELSE
1364                  sstdatqc%vdmean(ji,jj) = zmeanday(ji,jj) * zdaystp
1365               ENDIF
1366            END DO
1367         END DO
1368      ENDIF
1369
1370      ENDIF
1371
1372      ! Get the data for interpolation
1373     
1374      ALLOCATE( &
1375         & igrdi(2,2,isst), &
1376         & igrdj(2,2,isst), &
1377         & zglam(2,2,isst), &
1378         & zgphi(2,2,isst), &
1379         & zmask(2,2,isst), &
1380         & zsstl(2,2,isst)  &
1381         & )
1382     
1383      DO jobs = sstdatqc%nsurfup + 1, sstdatqc%nsurfup + isst
1384         iobs = jobs - sstdatqc%nsurfup
1385         igrdi(1,1,iobs) = sstdatqc%mi(jobs)-1
1386         igrdj(1,1,iobs) = sstdatqc%mj(jobs)-1
1387         igrdi(1,2,iobs) = sstdatqc%mi(jobs)-1
1388         igrdj(1,2,iobs) = sstdatqc%mj(jobs)
1389         igrdi(2,1,iobs) = sstdatqc%mi(jobs)
1390         igrdj(2,1,iobs) = sstdatqc%mj(jobs)-1
1391         igrdi(2,2,iobs) = sstdatqc%mi(jobs)
1392         igrdj(2,2,iobs) = sstdatqc%mj(jobs)
1393      END DO
1394     
1395      CALL obs_int_comm_2d( 2, 2, isst, &
1396         &                  igrdi, igrdj, glamt, zglam )
1397      CALL obs_int_comm_2d( 2, 2, isst, &
1398         &                  igrdi, igrdj, gphit, zgphi )
1399      CALL obs_int_comm_2d( 2, 2, isst, &
1400         &                  igrdi, igrdj, psstmask, zmask )
1401      CALL obs_int_comm_2d( 2, 2, isst, &
1402         &                  igrdi, igrdj, psstn, zsstl )
1403
1404      ! At the end of the day get interpolated means
1405      IF ( idayend == 0 .AND. ld_nightav ) THEN
1406
1407         ALLOCATE( &
1408            & zsstm(2,2,isst)  &
1409            & )
1410
1411         CALL obs_int_comm_2d( 2, 2, isst, igrdi, igrdj, &
1412            &               sstdatqc%vdmean(:,:), zsstm )
1413
1414      ENDIF
1415
1416      ! Loop over observations
1417
1418      DO jobs = sstdatqc%nsurfup + 1, sstdatqc%nsurfup + isst
1419         
1420         iobs = jobs - sstdatqc%nsurfup
1421         
1422         IF ( kt /= sstdatqc%mstp(jobs) ) THEN
1423           
1424            IF(lwp) THEN
1425               WRITE(numout,*)
1426               WRITE(numout,*) ' E R R O R : Observation',              &
1427                  &            ' time step is not consistent with the', &
1428                  &            ' model time step'
1429               WRITE(numout,*) ' ========='
1430               WRITE(numout,*)
1431               WRITE(numout,*) ' Record  = ', jobs,                &
1432                  &            ' kt      = ', kt,                  &
1433                  &            ' mstp    = ', sstdatqc%mstp(jobs), &
1434                  &            ' ntyp    = ', sstdatqc%ntyp(jobs)
1435            ENDIF
1436            CALL ctl_stop( 'obs_sst_opt', 'Inconsistent time' )
1437           
1438         ENDIF
1439         
1440         zlam = sstdatqc%rlam(jobs)
1441         zphi = sstdatqc%rphi(jobs)
1442         
1443         ! Get weights to interpolate the model SST to the observation point
1444         CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         &
1445            &                   zglam(:,:,iobs), zgphi(:,:,iobs), &
1446            &                   zmask(:,:,iobs), zweig, zobsmask )
1447           
1448         ! Interpolate the model SST to the observation point
1449
1450         IF ( ld_nightav ) THEN
1451
1452           IF ( idayend == 0 )  THEN
1453               ! Daily averaged/diurnal cycle of SST  data
1454               CALL obs_int_h2d( 1, 1,      & 
1455                     &              zweig, zsstm(:,:,iobs), zext )
1456            ELSE
1457               CALL ctl_stop( ' ld_nightav is set to true: a nonzero' //     &
1458                     &           ' number of night SST data should' // &
1459                     &           ' only occur at the end of a given day' )
1460            ENDIF
1461
1462         ELSE
1463
1464            CALL obs_int_h2d( 1, 1,      &
1465            &              zweig, zsstl(:,:,iobs),  zext )
1466
1467         ENDIF
1468         sstdatqc%rmod(jobs,1) = zext(1)
1469         
1470      END DO
1471     
1472      ! Deallocate the data for interpolation
1473      DEALLOCATE( &
1474         & igrdi, &
1475         & igrdj, &
1476         & zglam, &
1477         & zgphi, &
1478         & zmask, &
1479         & zsstl  &
1480         & )
1481
1482      ! At the end of the day also get interpolated means
1483      IF ( idayend == 0 .AND. ld_nightav ) THEN
1484         DEALLOCATE( &
1485            & zsstm  &
1486            & )
1487      ENDIF
1488     
1489      sstdatqc%nsurfup = sstdatqc%nsurfup + isst
1490
1491   END SUBROUTINE obs_sst_opt
1492
1493   SUBROUTINE obs_sss_opt
1494      !!-----------------------------------------------------------------------
1495      !!
1496      !!                     ***  ROUTINE obs_sss_opt  ***
1497      !!
1498      !! ** Purpose : Compute the model counterpart of sea surface salinity
1499      !!              data by interpolating from the model grid to the
1500      !!              observation point.
1501      !!
1502      !! ** Method  :
1503      !!
1504      !! ** Action  :
1505      !!
1506      !! History :
1507      !!      ! ??-??
1508      !!-----------------------------------------------------------------------
1509
1510      IMPLICIT NONE
1511
1512   END SUBROUTINE obs_sss_opt
1513
1514   SUBROUTINE obs_seaice_opt( seaicedatqc, kt, kpi, kpj, kit000, &
1515      &                    pseaicen, pseaicemask, k2dint )
1516
1517      !!-----------------------------------------------------------------------
1518      !!
1519      !!                     ***  ROUTINE obs_seaice_opt  ***
1520      !!
1521      !! ** Purpose : Compute the model counterpart of surface temperature
1522      !!              data by interpolating from the model grid to the
1523      !!              observation point.
1524      !!
1525      !! ** Method  : Linearly interpolate to each observation point using
1526      !!              the model values at the corners of the surrounding grid box.
1527      !!
1528      !!    The now model sea ice is first computed at the obs (lon, lat) point.
1529      !!
1530      !!    Several horizontal interpolation schemes are available:
1531      !!        - distance-weighted (great circle) (k2dint = 0)
1532      !!        - distance-weighted (small angle)  (k2dint = 1)
1533      !!        - bilinear (geographical grid)     (k2dint = 2)
1534      !!        - bilinear (quadrilateral grid)    (k2dint = 3)
1535      !!        - polynomial (quadrilateral grid)  (k2dint = 4)
1536      !!
1537      !!
1538      !! ** Action  :
1539      !!
1540      !! History :
1541      !!        !  07-07  (S. Ricci ) : Original
1542      !!     
1543      !!-----------------------------------------------------------------------
1544
1545      !! * Modules used
1546      USE obs_surf_def  ! Definition of storage space for surface observations
1547
1548      IMPLICIT NONE
1549
1550      !! * Arguments
1551      TYPE(obs_surf), INTENT(INOUT) :: seaicedatqc     ! Subset of surface data not failing screening
1552      INTEGER, INTENT(IN) :: kt       ! Time step
1553      INTEGER, INTENT(IN) :: kpi      ! Model grid parameters
1554      INTEGER, INTENT(IN) :: kpj
1555      INTEGER, INTENT(IN) :: kit000   ! Number of the first time step
1556                                      !   (kit000-1 = restart time)
1557      INTEGER, INTENT(IN) :: k2dint   ! Horizontal interpolation type (see header)
1558      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
1559         & pseaicen,  &    ! Model sea ice field
1560         & pseaicemask     ! Land-sea mask
1561         
1562      !! * Local declarations
1563      INTEGER :: ji
1564      INTEGER :: jj
1565      INTEGER :: jobs
1566      INTEGER :: inrc
1567      INTEGER :: iseaice
1568      INTEGER :: iobs
1569       
1570      REAL(KIND=wp) :: zlam
1571      REAL(KIND=wp) :: zphi
1572      REAL(KIND=wp) :: zext(1), zobsmask(1)
1573      REAL(kind=wp), DIMENSION(2,2,1) :: &
1574         & zweig
1575      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
1576         & zmask, &
1577         & zseaicel, &
1578         & zglam, &
1579         & zgphi
1580      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
1581         & igrdi, &
1582         & igrdj
1583
1584      !------------------------------------------------------------------------
1585      ! Local initialization
1586      !------------------------------------------------------------------------
1587      ! ... Record and data counters
1588      inrc = kt - kit000 + 2
1589      iseaice = seaicedatqc%nsstp(inrc)
1590
1591      ! Get the data for interpolation
1592     
1593      ALLOCATE( &
1594         & igrdi(2,2,iseaice), &
1595         & igrdj(2,2,iseaice), &
1596         & zglam(2,2,iseaice), &
1597         & zgphi(2,2,iseaice), &
1598         & zmask(2,2,iseaice), &
1599         & zseaicel(2,2,iseaice)  &
1600         & )
1601     
1602      DO jobs = seaicedatqc%nsurfup + 1, seaicedatqc%nsurfup + iseaice
1603         iobs = jobs - seaicedatqc%nsurfup
1604         igrdi(1,1,iobs) = seaicedatqc%mi(jobs)-1
1605         igrdj(1,1,iobs) = seaicedatqc%mj(jobs)-1
1606         igrdi(1,2,iobs) = seaicedatqc%mi(jobs)-1
1607         igrdj(1,2,iobs) = seaicedatqc%mj(jobs)
1608         igrdi(2,1,iobs) = seaicedatqc%mi(jobs)
1609         igrdj(2,1,iobs) = seaicedatqc%mj(jobs)-1
1610         igrdi(2,2,iobs) = seaicedatqc%mi(jobs)
1611         igrdj(2,2,iobs) = seaicedatqc%mj(jobs)
1612      END DO
1613     
1614      CALL obs_int_comm_2d( 2, 2, iseaice, &
1615         &                  igrdi, igrdj, glamt, zglam )
1616      CALL obs_int_comm_2d( 2, 2, iseaice, &
1617         &                  igrdi, igrdj, gphit, zgphi )
1618      CALL obs_int_comm_2d( 2, 2, iseaice, &
1619         &                  igrdi, igrdj, pseaicemask, zmask )
1620      CALL obs_int_comm_2d( 2, 2, iseaice, &
1621         &                  igrdi, igrdj, pseaicen, zseaicel )
1622     
1623      DO jobs = seaicedatqc%nsurfup + 1, seaicedatqc%nsurfup + iseaice
1624         
1625         iobs = jobs - seaicedatqc%nsurfup
1626         
1627         IF ( kt /= seaicedatqc%mstp(jobs) ) THEN
1628           
1629            IF(lwp) THEN
1630               WRITE(numout,*)
1631               WRITE(numout,*) ' E R R O R : Observation',              &
1632                  &            ' time step is not consistent with the', &
1633                  &            ' model time step'
1634               WRITE(numout,*) ' ========='
1635               WRITE(numout,*)
1636               WRITE(numout,*) ' Record  = ', jobs,                &
1637                  &            ' kt      = ', kt,                  &
1638                  &            ' mstp    = ', seaicedatqc%mstp(jobs), &
1639                  &            ' ntyp    = ', seaicedatqc%ntyp(jobs)
1640            ENDIF
1641            CALL ctl_stop( 'obs_seaice_opt', 'Inconsistent time' )
1642           
1643         ENDIF
1644         
1645         zlam = seaicedatqc%rlam(jobs)
1646         zphi = seaicedatqc%rphi(jobs)
1647         
1648         ! Get weights to interpolate the model sea ice to the observation point
1649         CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         &
1650            &                   zglam(:,:,iobs), zgphi(:,:,iobs), &
1651            &                   zmask(:,:,iobs), zweig, zobsmask )
1652         
1653         ! ... Interpolate the model sea ice to the observation point
1654         CALL obs_int_h2d( 1, 1,      &
1655            &              zweig, zseaicel(:,:,iobs),  zext )
1656         
1657         seaicedatqc%rmod(jobs,1) = zext(1)
1658         
1659      END DO
1660     
1661      ! Deallocate the data for interpolation
1662      DEALLOCATE( &
1663         & igrdi,    &
1664         & igrdj,    &
1665         & zglam,    &
1666         & zgphi,    &
1667         & zmask,    &
1668         & zseaicel  &
1669         & )
1670     
1671      seaicedatqc%nsurfup = seaicedatqc%nsurfup + iseaice
1672
1673   END SUBROUTINE obs_seaice_opt
1674
1675   SUBROUTINE obs_vel_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, &
1676      &                    pun, pvn, pgdept, pumask, pvmask, k1dint, k2dint, &
1677      &                    ld_dailyav )
1678      !!-----------------------------------------------------------------------
1679      !!
1680      !!                     ***  ROUTINE obs_vel_opt  ***
1681      !!
1682      !! ** Purpose : Compute the model counterpart of velocity profile
1683      !!              data by interpolating from the model grid to the
1684      !!              observation point.
1685      !!
1686      !! ** Method  : Linearly interpolate zonal and meridional components of velocity
1687      !!              to each observation point using the model values at the corners of
1688      !!              the surrounding grid box. The model velocity components are on a
1689      !!              staggered C- grid.
1690      !!
1691      !!    For velocity data from the TAO array, the model equivalent is
1692      !!    a daily mean velocity field. So, we first compute
1693      !!    the mean, then interpolate only at the end of the day.
1694      !!
1695      !! ** Action  :
1696      !!
1697      !! History :
1698      !!    ! 07-03 (K. Mogensen)      : Temperature and Salinity profiles
1699      !!    ! 08-10 (Maria Valdivieso) : Velocity component (U,V) profiles
1700      !!-----------------------------------------------------------------------
1701   
1702      !! * Modules used
1703      USE obs_profiles_def ! Definition of storage space for profile obs.
1704
1705      IMPLICIT NONE
1706
1707      !! * Arguments
1708      TYPE(obs_prof), INTENT(INOUT) :: &
1709         & prodatqc        ! Subset of profile data not failing screening
1710      INTEGER, INTENT(IN) :: kt        ! Time step
1711      INTEGER, INTENT(IN) :: kpi       ! Model grid parameters
1712      INTEGER, INTENT(IN) :: kpj
1713      INTEGER, INTENT(IN) :: kpk 
1714      INTEGER, INTENT(IN) :: kit000    ! Number of the first time step
1715                                       !   (kit000-1 = restart time)
1716      INTEGER, INTENT(IN) :: k1dint    ! Vertical interpolation type (see header)
1717      INTEGER, INTENT(IN) :: k2dint    ! Horizontal interpolation type (see header)
1718      INTEGER, INTENT(IN) :: kdaystp   ! Number of time steps per day                   
1719      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &
1720         & pun,    &    ! Model zonal component of velocity
1721         & pvn,    &    ! Model meridional component of velocity
1722         & pumask, &    ! Land-sea mask
1723         & pvmask       ! Land-sea mask
1724      REAL(KIND=wp), INTENT(IN), DIMENSION(kpk) :: &
1725         & pgdept       ! Model array of depth levels
1726      LOGICAL, INTENT(IN) :: ld_dailyav
1727         
1728      !! * Local declarations
1729      INTEGER :: ji
1730      INTEGER :: jj
1731      INTEGER :: jk
1732      INTEGER :: jobs
1733      INTEGER :: inrc
1734      INTEGER :: ipro
1735      INTEGER :: idayend
1736      INTEGER :: ista
1737      INTEGER :: iend
1738      INTEGER :: iobs
1739      INTEGER, DIMENSION(imaxavtypes) :: &
1740         & idailyavtypes
1741      REAL(KIND=wp) :: zlam
1742      REAL(KIND=wp) :: zphi
1743      REAL(KIND=wp) :: zdaystp
1744      REAL(KIND=wp), DIMENSION(kpk) :: &
1745         & zobsmasku, &
1746         & zobsmaskv, &
1747         & zobsmask,  &
1748         & zobsk,     &
1749         & zobs2k
1750      REAL(KIND=wp), DIMENSION(2,2,kpk) :: &
1751         & zweigu,zweigv
1752      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
1753         & zumask, zvmask, &
1754         & zintu, &
1755         & zintv, &
1756         & zinmu, &
1757         & zinmv
1758      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
1759         & zglamu, zglamv, &
1760         & zgphiu, zgphiv
1761      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
1762         & igrdiu, &
1763         & igrdju, &
1764         & igrdiv, &
1765         & igrdjv
1766
1767      !------------------------------------------------------------------------
1768      ! Local initialization
1769      !------------------------------------------------------------------------
1770      ! ... Record and data counters
1771      inrc = kt - kit000 + 2
1772      ipro = prodatqc%npstp(inrc)
1773
1774      ! Initialize daily mean for first timestep
1775      idayend = MOD( kt - kit000 + 1, kdaystp )
1776
1777      ! Added kt == 0 test to catch restart case
1778      IF ( idayend == 1 .OR. kt == 0) THEN
1779         IF (lwp) WRITE(numout,*) 'Reset prodatqc%vdmean on time-step: ',kt
1780         prodatqc%vdmean(:,:,:,1) = 0.0
1781         prodatqc%vdmean(:,:,:,2) = 0.0
1782      ENDIF
1783
1784      ! Increment the zonal velocity field for computing daily mean
1785      prodatqc%vdmean(:,:,:,1) = prodatqc%vdmean(:,:,:,1) + pun(:,:,:)
1786      ! Increment the meridional velocity field for computing daily mean
1787      prodatqc%vdmean(:,:,:,2) = prodatqc%vdmean(:,:,:,2) + pvn(:,:,:)
1788   
1789      ! Compute the daily mean at the end of day
1790      zdaystp = 1.0 / REAL( kdaystp )
1791      IF ( idayend == 0 ) THEN
1792         prodatqc%vdmean(:,:,:,1) = prodatqc%vdmean(:,:,:,1) * zdaystp
1793         prodatqc%vdmean(:,:,:,2) = prodatqc%vdmean(:,:,:,2) * zdaystp
1794      ENDIF
1795
1796      ! Get the data for interpolation
1797      ALLOCATE( &
1798         & igrdiu(2,2,ipro),      &
1799         & igrdju(2,2,ipro),      &
1800         & igrdiv(2,2,ipro),      &
1801         & igrdjv(2,2,ipro),      &
1802         & zglamu(2,2,ipro), zglamv(2,2,ipro), &
1803         & zgphiu(2,2,ipro), zgphiv(2,2,ipro), &
1804         & zumask(2,2,kpk,ipro), zvmask(2,2,kpk,ipro), &
1805         & zintu(2,2,kpk,ipro),  &
1806         & zintv(2,2,kpk,ipro)   &
1807         & )
1808
1809      DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
1810         iobs = jobs - prodatqc%nprofup
1811         igrdiu(1,1,iobs) = prodatqc%mi(jobs,1)-1
1812         igrdju(1,1,iobs) = prodatqc%mj(jobs,1)-1
1813         igrdiu(1,2,iobs) = prodatqc%mi(jobs,1)-1
1814         igrdju(1,2,iobs) = prodatqc%mj(jobs,1)
1815         igrdiu(2,1,iobs) = prodatqc%mi(jobs,1)
1816         igrdju(2,1,iobs) = prodatqc%mj(jobs,1)-1
1817         igrdiu(2,2,iobs) = prodatqc%mi(jobs,1)
1818         igrdju(2,2,iobs) = prodatqc%mj(jobs,1)
1819         igrdiv(1,1,iobs) = prodatqc%mi(jobs,2)-1
1820         igrdjv(1,1,iobs) = prodatqc%mj(jobs,2)-1
1821         igrdiv(1,2,iobs) = prodatqc%mi(jobs,2)-1
1822         igrdjv(1,2,iobs) = prodatqc%mj(jobs,2)
1823         igrdiv(2,1,iobs) = prodatqc%mi(jobs,2)
1824         igrdjv(2,1,iobs) = prodatqc%mj(jobs,2)-1
1825         igrdiv(2,2,iobs) = prodatqc%mi(jobs,2)
1826         igrdjv(2,2,iobs) = prodatqc%mj(jobs,2)
1827      END DO
1828
1829      CALL obs_int_comm_2d( 2, 2, ipro, igrdiu, igrdju, glamu, zglamu )
1830      CALL obs_int_comm_2d( 2, 2, ipro, igrdiu, igrdju, gphiu, zgphiu )
1831      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, pumask, zumask )
1832      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, pun, zintu )
1833
1834      CALL obs_int_comm_2d( 2, 2, ipro, igrdiv, igrdjv, glamv, zglamv )
1835      CALL obs_int_comm_2d( 2, 2, ipro, igrdiv, igrdjv, gphiv, zgphiv )
1836      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, pvmask, zvmask )
1837      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, pvn, zintv )
1838
1839      ! At the end of the day also get interpolated means
1840      IF ( idayend == 0 ) THEN
1841
1842         ALLOCATE( &
1843            & zinmu(2,2,kpk,ipro),  &
1844            & zinmv(2,2,kpk,ipro)   &
1845            & )
1846
1847         CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, &
1848            &                  prodatqc%vdmean(:,:,:,1), zinmu )
1849         CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, &
1850            &                  prodatqc%vdmean(:,:,:,2), zinmv )
1851
1852      ENDIF
1853
1854! loop over observations
1855
1856      DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
1857
1858         iobs = jobs - prodatqc%nprofup
1859
1860         IF ( kt /= prodatqc%mstp(jobs) ) THEN
1861           
1862            IF(lwp) THEN
1863               WRITE(numout,*)
1864               WRITE(numout,*) ' E R R O R : Observation',              &
1865                  &            ' time step is not consistent with the', &
1866                  &            ' model time step'
1867               WRITE(numout,*) ' ========='
1868               WRITE(numout,*)
1869               WRITE(numout,*) ' Record  = ', jobs,                    &
1870                  &            ' kt      = ', kt,                      &
1871                  &            ' mstp    = ', prodatqc%mstp(jobs), &
1872                  &            ' ntyp    = ', prodatqc%ntyp(jobs)
1873            ENDIF
1874            CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' )
1875         ENDIF
1876         
1877         zlam = prodatqc%rlam(jobs)
1878         zphi = prodatqc%rphi(jobs)
1879
1880         ! Initialize observation masks
1881
1882         zobsmasku(:) = 0.0
1883         zobsmaskv(:) = 0.0
1884         
1885         ! Horizontal weights and vertical mask
1886
1887         IF  ( prodatqc%npvend(jobs,1) > 0 ) THEN
1888
1889            CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi,     &
1890               &                   zglamu(:,:,iobs), zgphiu(:,:,iobs), &
1891               &                   zumask(:,:,:,iobs), zweigu, zobsmasku )
1892
1893         ENDIF
1894
1895         
1896         IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
1897
1898            CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi,     &
1899               &                   zglamv(:,:,iobs), zgphiv(:,:,iobs), &
1900               &                   zvmask(:,:,:,iobs), zweigv, zobsmaskv )
1901
1902         ENDIF
1903
1904         ! Ensure that the vertical mask on u and v are consistent.
1905
1906         zobsmask(:) = MIN( zobsmasku(:), zobsmaskv(:) )
1907
1908         IF ( prodatqc%npvend(jobs,1) > 0 ) THEN
1909
1910            zobsk(:) = obfillflt
1911
1912       IF ( ld_dailyav ) THEN
1913
1914               IF ( idayend == 0 )  THEN
1915                 
1916                  ! Daily averaged data
1917                 
1918                  CALL obs_int_h2d( kpk, kpk,      &
1919                     &              zweigu, zinmu(:,:,:,iobs), zobsk )
1920                 
1921                 
1922               ELSE
1923               
1924                  CALL ctl_stop( ' A nonzero' //     &
1925                     &           ' number of U profile data should' // &
1926                     &           ' only occur at the end of a given day' )
1927
1928               ENDIF
1929         
1930            ELSE 
1931               
1932               ! Point data
1933
1934               CALL obs_int_h2d( kpk, kpk,      &
1935                  &              zweigu, zintu(:,:,:,iobs), zobsk )
1936
1937            ENDIF
1938
1939            !-------------------------------------------------------------
1940            ! Compute vertical second-derivative of the interpolating
1941            ! polynomial at obs points
1942            !-------------------------------------------------------------
1943           
1944            IF ( k1dint == 1 ) THEN
1945               CALL obs_int_z1d_spl( kpk, zobsk, zobs2k,   &
1946                  &                  pgdept, zobsmask )
1947            ENDIF
1948           
1949            !-----------------------------------------------------------------
1950            !  Vertical interpolation to the observation point
1951            !-----------------------------------------------------------------
1952            ista = prodatqc%npvsta(jobs,1)
1953            iend = prodatqc%npvend(jobs,1)
1954            CALL obs_int_z1d( kpk,                &
1955               & prodatqc%var(1)%mvk(ista:iend),  &
1956               & k1dint, iend - ista + 1,         &
1957               & prodatqc%var(1)%vdep(ista:iend), &
1958               & zobsk, zobs2k,                   &
1959               & prodatqc%var(1)%vmod(ista:iend), &
1960               & pgdept, zobsmask )
1961
1962         ENDIF
1963
1964         IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
1965
1966            zobsk(:) = obfillflt
1967
1968            IF ( ld_dailyav ) THEN
1969
1970               IF ( idayend == 0 )  THEN
1971
1972                  ! Daily averaged data
1973                 
1974                  CALL obs_int_h2d( kpk, kpk,      &
1975                     &              zweigv, zinmv(:,:,:,iobs), zobsk )
1976                 
1977               ELSE
1978
1979                  CALL ctl_stop( ' A nonzero' //     &
1980                     &           ' number of V profile data should' // &
1981                     &           ' only occur at the end of a given day' )
1982
1983               ENDIF
1984
1985            ELSE
1986               
1987               ! Point data
1988
1989               CALL obs_int_h2d( kpk, kpk,      &
1990                  &              zweigv, zintv(:,:,:,iobs), zobsk )
1991
1992            ENDIF
1993
1994
1995            !-------------------------------------------------------------
1996            ! Compute vertical second-derivative of the interpolating
1997            ! polynomial at obs points
1998            !-------------------------------------------------------------
1999           
2000            IF ( k1dint == 1 ) THEN
2001               CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, &
2002                  &                  pgdept, zobsmask )
2003            ENDIF
2004           
2005            !----------------------------------------------------------------
2006            !  Vertical interpolation to the observation point
2007            !----------------------------------------------------------------
2008            ista = prodatqc%npvsta(jobs,2)
2009            iend = prodatqc%npvend(jobs,2)
2010            CALL obs_int_z1d( kpk, &
2011               & prodatqc%var(2)%mvk(ista:iend),&
2012               & k1dint, iend - ista + 1, &
2013               & prodatqc%var(2)%vdep(ista:iend),&
2014               & zobsk, zobs2k, &
2015               & prodatqc%var(2)%vmod(ista:iend),&
2016               & pgdept, zobsmask )
2017
2018         ENDIF
2019
2020      END DO
2021 
2022      ! Deallocate the data for interpolation
2023      DEALLOCATE( &
2024         & igrdiu, &
2025         & igrdju, &
2026         & igrdiv, &
2027         & igrdjv, &
2028         & zglamu, zglamv, &
2029         & zgphiu, zgphiv, &
2030         & zumask, zvmask, &
2031         & zintu, &
2032         & zintv  &
2033         & )
2034      ! At the end of the day also get interpolated means
2035      IF ( idayend == 0 ) THEN
2036         DEALLOCATE( &
2037            & zinmu,  &
2038            & zinmv   &
2039            & )
2040      ENDIF
2041
2042      prodatqc%nprofup = prodatqc%nprofup + ipro 
2043     
2044   END SUBROUTINE obs_vel_opt
2045
2046END MODULE obs_oper
2047
Note: See TracBrowser for help on using the repository browser.