New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
obs_oper.F90 in branches/UKMO/dev_r4650_general_vert_coord_obsoper_removetides/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

source: branches/UKMO/dev_r4650_general_vert_coord_obsoper_removetides/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90 @ 7375

Last change on this file since 7375 was 7375, checked in by mattmartin, 5 years ago

Added code to branch to allow sla model counterpart to be a time-average.

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