source: branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90 @ 5838

Last change on this file since 5838 was 5838, checked in by timgraham, 5 years ago

Cleared svn keywords to allow use with fcm make at vn3.6

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