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

source: branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90 @ 10246

Last change on this file since 10246 was 10246, checked in by kingr, 6 years ago

Cleared SVN keywords.

File size: 51.8 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_sla_opt :    Compute the model counterpart of sea level anomaly
12   !!                    observations
13   !!   obs_sst_opt :    Compute the model counterpart of sea surface temperature
14   !!                    observations
15   !!   obs_sss_opt :    Compute the model counterpart of sea surface salinity
16   !!                    observations
17   !!   obs_seaice_opt : Compute the model counterpart of sea ice concentration
18   !!                    observations
19   !!
20   !!   obs_vel_opt :    Compute the model counterpart of zonal and meridional
21   !!                    components of velocity from observations.
22   !!----------------------------------------------------------------------
23
24   !! * Modules used   
25   USE par_kind, ONLY : &         ! Precision variables
26      & wp
27   USE in_out_manager             ! I/O manager
28   USE obs_inter_sup              ! Interpolation support
29   USE obs_inter_h2d, ONLY : &    ! Horizontal interpolation to the observation pt
30      & obs_int_h2d, &
31      & obs_int_h2d_init
32   USE obs_inter_z1d, ONLY : &    ! Vertical interpolation to the observation pt
33      & obs_int_z1d,    &
34      & obs_int_z1d_spl
35   USE obs_const,  ONLY :     &
36      & obfillflt      ! Fillvalue   
37   USE dom_oce,       ONLY : &
38      & glamt, glamu, glamv, &
39      & gphit, gphiu, gphiv
40   USE lib_mpp,       ONLY : &
41      & ctl_warn, ctl_stop
42
43   IMPLICIT NONE
44
45   !! * Routine accessibility
46   PRIVATE
47
48   PUBLIC obs_pro_opt, &  ! Compute the model counterpart of profile observations
49      &   obs_sla_opt, &  ! Compute the model counterpart of SLA observations
50      &   obs_sst_opt, &  ! Compute the model counterpart of SST observations
51      &   obs_sss_opt, &  ! Compute the model counterpart of SSS observations
52      &   obs_seaice_opt, &
53      &   obs_vel_opt     ! Compute the model counterpart of velocity profile data
54
55   INTEGER, PARAMETER, PUBLIC :: imaxavtypes = 20 ! Max number of daily avgd obs types
56
57   !!----------------------------------------------------------------------
58   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
59   !! $Id$
60   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
61   !!----------------------------------------------------------------------
62
63CONTAINS
64
65   SUBROUTINE obs_pro_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, &
66      &                    ptn, psn, pgdept, ptmask, k1dint, k2dint, &
67      &                    kdailyavtypes )
68      !!-----------------------------------------------------------------------
69      !!
70      !!                     ***  ROUTINE obs_pro_opt  ***
71      !!
72      !! ** Purpose : Compute the model counterpart of profiles
73      !!              data by interpolating from the model grid to the
74      !!              observation point.
75      !!
76      !! ** Method  : Linearly interpolate to each observation point using
77      !!              the model values at the corners of the surrounding grid box.
78      !!
79      !!    First, a vertical profile of horizontally interpolated model
80      !!    now temperatures is computed at the obs (lon, lat) point.
81      !!    Several horizontal interpolation schemes are available:
82      !!        - distance-weighted (great circle) (k2dint = 0)
83      !!        - distance-weighted (small angle)  (k2dint = 1)
84      !!        - bilinear (geographical grid)     (k2dint = 2)
85      !!        - bilinear (quadrilateral grid)    (k2dint = 3)
86      !!        - polynomial (quadrilateral grid)  (k2dint = 4)
87      !!
88      !!    Next, the vertical temperature profile is interpolated to the
89      !!    data depth points. Two vertical interpolation schemes are
90      !!    available:
91      !!        - linear       (k1dint = 0)
92      !!        - Cubic spline (k1dint = 1)
93      !!
94      !!    For the cubic spline the 2nd derivative of the interpolating
95      !!    polynomial is computed before entering the vertical interpolation
96      !!    routine.
97      !!
98      !!    For ENACT moored buoy data (e.g., TAO), the model equivalent is
99      !!    a daily mean model temperature field. So, we first compute
100      !!    the mean, then interpolate only at the end of the day.
101      !!
102      !!    Note: the in situ temperature observations must be converted
103      !!    to potential temperature (the model variable) prior to
104      !!    assimilation.
105      !!??????????????????????????????????????????????????????????????
106      !!    INCLUDE POTENTIAL TEMP -> IN SITU TEMP IN OBS OPERATOR???
107      !!??????????????????????????????????????????????????????????????
108      !!
109      !! ** Action  :
110      !!
111      !! History :
112      !!      ! 97-11 (A. Weaver, S. Ricci, N. Daget)
113      !!      ! 06-03 (G. Smith) NEMOVAR migration
114      !!      ! 06-10 (A. Weaver) Cleanup
115      !!      ! 07-01 (K. Mogensen) Merge of temperature and salinity
116      !!      ! 07-03 (K. Mogensen) General handling of profiles
117      !!-----------------------------------------------------------------------
118 
119      !! * Modules used
120      USE obs_profiles_def ! Definition of storage space for profile obs.
121
122      IMPLICIT NONE
123
124      !! * Arguments
125      TYPE(obs_prof), INTENT(INOUT) :: prodatqc  ! Subset of profile data not failing screening
126      INTEGER, INTENT(IN) :: kt        ! Time step
127      INTEGER, INTENT(IN) :: kpi       ! Model grid parameters
128      INTEGER, INTENT(IN) :: kpj
129      INTEGER, INTENT(IN) :: kpk
130      INTEGER, INTENT(IN) :: kit000    ! Number of the first time step
131                                       !   (kit000-1 = restart time)
132      INTEGER, INTENT(IN) :: k1dint    ! Vertical interpolation type (see header)
133      INTEGER, INTENT(IN) :: k2dint    ! Horizontal interpolation type (see header)
134      INTEGER, INTENT(IN) :: kdaystp   ! Number of time steps per day                   
135      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &
136         & ptn,    &    ! Model temperature field
137         & psn,    &    ! Model salinity field
138         & ptmask       ! Land-sea mask
139      REAL(KIND=wp), INTENT(IN), DIMENSION(kpk) :: &
140         & pgdept       ! Model array of depth levels
141      INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &
142         & kdailyavtypes! Types for daily averages
143      !! * Local declarations
144      INTEGER ::   ji
145      INTEGER ::   jj
146      INTEGER ::   jk
147      INTEGER ::   jobs
148      INTEGER ::   inrc
149      INTEGER ::   ipro
150      INTEGER ::   idayend
151      INTEGER ::   ista
152      INTEGER ::   iend
153      INTEGER ::   iobs
154      INTEGER, DIMENSION(imaxavtypes) :: &
155         & idailyavtypes
156      REAL(KIND=wp) :: zlam
157      REAL(KIND=wp) :: zphi
158      REAL(KIND=wp) :: zdaystp
159      REAL(KIND=wp), DIMENSION(kpk) :: &
160         & zobsmask, &
161         & zobsk,    &
162         & zobs2k
163      REAL(KIND=wp), DIMENSION(2,2,kpk) :: &
164         & zweig
165      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
166         & zmask, &
167         & zintt, &
168         & zints, &
169         & zinmt, &
170         & zinms
171      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
172         & zglam, &
173         & zgphi
174      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
175         & igrdi, &
176         & igrdj
177
178      !------------------------------------------------------------------------
179      ! Local initialization
180      !------------------------------------------------------------------------
181      ! ... Record and data counters
182      inrc = kt - kit000 + 2
183      ipro = prodatqc%npstp(inrc)
184 
185      ! Daily average types
186      IF ( PRESENT(kdailyavtypes) ) THEN
187         idailyavtypes(:) = kdailyavtypes(:)
188      ELSE
189         idailyavtypes(:) = -1
190      ENDIF
191
192      ! Initialize daily mean for first timestep
193      idayend = MOD( kt - kit000 + 1, kdaystp )
194
195      ! Added kt == 0 test to catch restart case
196      IF ( idayend == 1 .OR. kt == 0) THEN
197         IF (lwp) WRITE(numout,*) 'Reset prodatqc%vdmean on time-step: ',kt
198         DO jk = 1, jpk
199            DO jj = 1, jpj
200               DO ji = 1, jpi
201                  prodatqc%vdmean(ji,jj,jk,1) = 0.0
202                  prodatqc%vdmean(ji,jj,jk,2) = 0.0
203               END DO
204            END DO
205         END DO
206      ENDIF
207
208      DO jk = 1, jpk
209         DO jj = 1, jpj
210            DO ji = 1, jpi
211               ! Increment the temperature field for computing daily mean
212               prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &
213                  &                        + ptn(ji,jj,jk)
214               ! Increment the salinity field for computing daily mean
215               prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &
216                  &                        + psn(ji,jj,jk)
217            END DO
218         END DO
219      END DO
220   
221      ! Compute the daily mean at the end of day
222      zdaystp = 1.0 / REAL( kdaystp )
223      IF ( idayend == 0 ) THEN
224         DO jk = 1, jpk
225            DO jj = 1, jpj
226               DO ji = 1, jpi
227                  prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &
228                     &                        * zdaystp
229                  prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &
230                  &                           * zdaystp
231               END DO
232            END DO
233         END DO
234      ENDIF
235
236      ! Get the data for interpolation
237      ALLOCATE( &
238         & igrdi(2,2,ipro),      &
239         & igrdj(2,2,ipro),      &
240         & zglam(2,2,ipro),      &
241         & zgphi(2,2,ipro),      &
242         & zmask(2,2,kpk,ipro),  &
243         & zintt(2,2,kpk,ipro),  &
244         & zints(2,2,kpk,ipro)   &
245         & )
246
247      DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
248         iobs = jobs - prodatqc%nprofup
249         igrdi(1,1,iobs) = prodatqc%mi(jobs,1)-1
250         igrdj(1,1,iobs) = prodatqc%mj(jobs,1)-1
251         igrdi(1,2,iobs) = prodatqc%mi(jobs,1)-1
252         igrdj(1,2,iobs) = prodatqc%mj(jobs,1)
253         igrdi(2,1,iobs) = prodatqc%mi(jobs,1)
254         igrdj(2,1,iobs) = prodatqc%mj(jobs,1)-1
255         igrdi(2,2,iobs) = prodatqc%mi(jobs,1)
256         igrdj(2,2,iobs) = prodatqc%mj(jobs,1)
257      END DO
258
259      CALL obs_int_comm_2d( 2, 2, ipro, igrdi, igrdj, glamt, zglam )
260      CALL obs_int_comm_2d( 2, 2, ipro, igrdi, igrdj, gphit, zgphi )
261      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, ptmask,zmask )
262      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, ptn,   zintt )
263      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, psn,   zints )
264
265      ! At the end of the day also get interpolated means
266      IF ( idayend == 0 ) THEN
267
268         ALLOCATE( &
269            & zinmt(2,2,kpk,ipro),  &
270            & zinms(2,2,kpk,ipro)   &
271            & )
272
273         CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, &
274            &                  prodatqc%vdmean(:,:,:,1), zinmt )
275         CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, &
276            &                  prodatqc%vdmean(:,:,:,2), zinms )
277
278      ENDIF
279
280      DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
281
282         iobs = jobs - prodatqc%nprofup
283
284         IF ( kt /= prodatqc%mstp(jobs) ) THEN
285           
286            IF(lwp) THEN
287               WRITE(numout,*)
288               WRITE(numout,*) ' E R R O R : Observation',              &
289                  &            ' time step is not consistent with the', &
290                  &            ' model time step'
291               WRITE(numout,*) ' ========='
292               WRITE(numout,*)
293               WRITE(numout,*) ' Record  = ', jobs,                    &
294                  &            ' kt      = ', kt,                      &
295                  &            ' mstp    = ', prodatqc%mstp(jobs), &
296                  &            ' ntyp    = ', prodatqc%ntyp(jobs)
297            ENDIF
298            CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' )
299         ENDIF
300         
301         zlam = prodatqc%rlam(jobs)
302         zphi = prodatqc%rphi(jobs)
303         
304         ! Horizontal weights and vertical mask
305
306         IF ( ( prodatqc%npvend(jobs,1) > 0 ) .OR. &
307            & ( prodatqc%npvend(jobs,2) > 0 ) ) THEN
308
309            CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi,     &
310               &                   zglam(:,:,iobs), zgphi(:,:,iobs), &
311               &                   zmask(:,:,:,iobs), zweig, zobsmask )
312
313         ENDIF
314
315         IF ( prodatqc%npvend(jobs,1) > 0 ) THEN
316
317            zobsk(:) = obfillflt
318
319       IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN
320
321               IF ( idayend == 0 )  THEN
322                 
323                  ! Daily averaged moored buoy (MRB) data
324                 
325                  CALL obs_int_h2d( kpk, kpk,      &
326                     &              zweig, zinmt(:,:,:,iobs), zobsk )
327                 
328                 
329               ELSE
330               
331                  CALL ctl_stop( ' A nonzero' //     &
332                     &           ' number of profile T BUOY data should' // &
333                     &           ' only occur at the end of a given day' )
334
335               ENDIF
336         
337            ELSE 
338               
339               ! Point data
340
341               CALL obs_int_h2d( kpk, kpk,      &
342                  &              zweig, zintt(:,:,:,iobs), zobsk )
343
344            ENDIF
345
346            !-------------------------------------------------------------
347            ! Compute vertical second-derivative of the interpolating
348            ! polynomial at obs points
349            !-------------------------------------------------------------
350           
351            IF ( k1dint == 1 ) THEN
352               CALL obs_int_z1d_spl( kpk, zobsk, zobs2k,   &
353                  &                  pgdept, zobsmask )
354            ENDIF
355           
356            !-----------------------------------------------------------------
357            !  Vertical interpolation to the observation point
358            !-----------------------------------------------------------------
359            ista = prodatqc%npvsta(jobs,1)
360            iend = prodatqc%npvend(jobs,1)
361            CALL obs_int_z1d( kpk,                &
362               & prodatqc%var(1)%mvk(ista:iend),  &
363               & k1dint, iend - ista + 1,         &
364               & prodatqc%var(1)%vdep(ista:iend), &
365               & zobsk, zobs2k,                   &
366               & prodatqc%var(1)%vmod(ista:iend), &
367               & pgdept, zobsmask )
368
369         ENDIF
370
371         IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
372
373            zobsk(:) = obfillflt
374
375            IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN
376
377               IF ( idayend == 0 )  THEN
378
379                  ! Daily averaged moored buoy (MRB) data
380                 
381                  CALL obs_int_h2d( kpk, kpk,      &
382                     &              zweig, zinms(:,:,:,iobs), zobsk )
383                 
384               ELSE
385
386                  CALL ctl_stop( ' A nonzero' //     &
387                     &           ' number of profile S BUOY data should' // &
388                     &           ' only occur at the end of a given day' )
389
390               ENDIF
391
392            ELSE
393               
394               ! Point data
395
396               CALL obs_int_h2d( kpk, kpk,      &
397                  &              zweig, zints(:,:,:,iobs), zobsk )
398
399            ENDIF
400
401
402            !-------------------------------------------------------------
403            ! Compute vertical second-derivative of the interpolating
404            ! polynomial at obs points
405            !-------------------------------------------------------------
406           
407            IF ( k1dint == 1 ) THEN
408               CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, &
409                  &                  pgdept, zobsmask )
410            ENDIF
411           
412            !----------------------------------------------------------------
413            !  Vertical interpolation to the observation point
414            !----------------------------------------------------------------
415            ista = prodatqc%npvsta(jobs,2)
416            iend = prodatqc%npvend(jobs,2)
417            CALL obs_int_z1d( kpk, &
418               & prodatqc%var(2)%mvk(ista:iend),&
419               & k1dint, iend - ista + 1, &
420               & prodatqc%var(2)%vdep(ista:iend),&
421               & zobsk, zobs2k, &
422               & prodatqc%var(2)%vmod(ista:iend),&
423               & pgdept, zobsmask )
424
425         ENDIF
426
427      END DO
428 
429      ! Deallocate the data for interpolation
430      DEALLOCATE( &
431         & igrdi, &
432         & igrdj, &
433         & zglam, &
434         & zgphi, &
435         & zmask, &
436         & zintt, &
437         & zints  &
438         & )
439      ! At the end of the day also get interpolated means
440      IF ( idayend == 0 ) THEN
441         DEALLOCATE( &
442            & zinmt,  &
443            & zinms   &
444            & )
445      ENDIF
446
447      prodatqc%nprofup = prodatqc%nprofup + ipro 
448     
449   END SUBROUTINE obs_pro_opt
450
451   SUBROUTINE obs_sla_opt( sladatqc, kt, kpi, kpj, kit000, &
452      &                    psshn, psshmask, k2dint )
453      !!-----------------------------------------------------------------------
454      !!
455      !!                     ***  ROUTINE obs_sla_opt  ***
456      !!
457      !! ** Purpose : Compute the model counterpart of sea level anomaly
458      !!              data by interpolating from the model grid to the
459      !!              observation point.
460      !!
461      !! ** Method  : Linearly interpolate to each observation point using
462      !!              the model values at the corners of the surrounding grid box.
463      !!
464      !!    The now model SSH is first computed at the obs (lon, lat) point.
465      !!
466      !!    Several horizontal interpolation schemes are available:
467      !!        - distance-weighted (great circle) (k2dint = 0)
468      !!        - distance-weighted (small angle)  (k2dint = 1)
469      !!        - bilinear (geographical grid)     (k2dint = 2)
470      !!        - bilinear (quadrilateral grid)    (k2dint = 3)
471      !!        - polynomial (quadrilateral grid)  (k2dint = 4)
472      !! 
473      !!    The sea level anomaly at the observation points is then computed
474      !!    by removing a mean dynamic topography (defined at the obs. point).
475      !!
476      !! ** Action  :
477      !!
478      !! History :
479      !!      ! 07-03 (A. Weaver)
480      !!-----------------------------------------------------------------------
481 
482      !! * Modules used
483      USE obs_surf_def  ! Definition of storage space for surface observations
484
485      IMPLICIT NONE
486
487      !! * Arguments
488      TYPE(obs_surf), INTENT(INOUT) :: sladatqc     ! Subset of surface data not failing screening
489      INTEGER, INTENT(IN) :: kt      ! Time step
490      INTEGER, INTENT(IN) :: kpi     ! Model grid parameters
491      INTEGER, INTENT(IN) :: kpj
492      INTEGER, INTENT(IN) :: kit000   ! Number of the first time step
493                                      !   (kit000-1 = restart time)
494      INTEGER, INTENT(IN) :: k2dint   ! Horizontal interpolation type (see header)
495      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
496         & psshn,  &    ! Model SSH field
497         & psshmask     ! Land-sea mask
498         
499      !! * Local declarations
500      INTEGER :: ji
501      INTEGER :: jj
502      INTEGER :: jobs
503      INTEGER :: inrc
504      INTEGER :: isla
505      INTEGER :: iobs
506      REAL(KIND=wp) :: zlam
507      REAL(KIND=wp) :: zphi
508      REAL(KIND=wp) :: zext(1), zobsmask(1)
509      REAL(kind=wp), DIMENSION(2,2,1) :: &
510         & zweig
511      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
512         & zmask, &
513         & zsshl, &
514         & zglam, &
515         & zgphi
516      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
517         & igrdi, &
518         & igrdj
519
520      !------------------------------------------------------------------------
521      ! Local initialization
522      !------------------------------------------------------------------------
523      ! ... Record and data counters
524      inrc = kt - kit000 + 2
525      isla = sladatqc%nsstp(inrc)
526
527      ! Get the data for interpolation
528
529      ALLOCATE( &
530         & igrdi(2,2,isla), &
531         & igrdj(2,2,isla), &
532         & zglam(2,2,isla), &
533         & zgphi(2,2,isla), &
534         & zmask(2,2,isla), &
535         & zsshl(2,2,isla)  &
536         & )
537     
538      DO jobs = sladatqc%nsurfup + 1, sladatqc%nsurfup + isla
539         iobs = jobs - sladatqc%nsurfup
540         igrdi(1,1,iobs) = sladatqc%mi(jobs)-1
541         igrdj(1,1,iobs) = sladatqc%mj(jobs)-1
542         igrdi(1,2,iobs) = sladatqc%mi(jobs)-1
543         igrdj(1,2,iobs) = sladatqc%mj(jobs)
544         igrdi(2,1,iobs) = sladatqc%mi(jobs)
545         igrdj(2,1,iobs) = sladatqc%mj(jobs)-1
546         igrdi(2,2,iobs) = sladatqc%mi(jobs)
547         igrdj(2,2,iobs) = sladatqc%mj(jobs)
548      END DO
549
550      CALL obs_int_comm_2d( 2, 2, isla, &
551         &                  igrdi, igrdj, glamt, zglam )
552      CALL obs_int_comm_2d( 2, 2, isla, &
553         &                  igrdi, igrdj, gphit, zgphi )
554      CALL obs_int_comm_2d( 2, 2, isla, &
555         &                  igrdi, igrdj, psshmask, zmask )
556      CALL obs_int_comm_2d( 2, 2, isla, &
557         &                  igrdi, igrdj, psshn, zsshl )
558
559      ! Loop over observations
560
561      DO jobs = sladatqc%nsurfup + 1, sladatqc%nsurfup + isla
562
563         iobs = jobs - sladatqc%nsurfup
564
565         IF ( kt /= sladatqc%mstp(jobs) ) THEN
566           
567            IF(lwp) THEN
568               WRITE(numout,*)
569               WRITE(numout,*) ' E R R O R : Observation',              &
570                  &            ' time step is not consistent with the', &
571                  &            ' model time step'
572               WRITE(numout,*) ' ========='
573               WRITE(numout,*)
574               WRITE(numout,*) ' Record  = ', jobs,                &
575                  &            ' kt      = ', kt,                  &
576                  &            ' mstp    = ', sladatqc%mstp(jobs), &
577                  &            ' ntyp    = ', sladatqc%ntyp(jobs)
578            ENDIF
579            CALL ctl_stop( 'obs_sla_opt', 'Inconsistent time' )
580           
581         ENDIF
582         
583         zlam = sladatqc%rlam(jobs)
584         zphi = sladatqc%rphi(jobs)
585
586         ! Get weights to interpolate the model SSH to the observation point
587         CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         &
588            &                   zglam(:,:,iobs), zgphi(:,:,iobs), &
589            &                   zmask(:,:,iobs), zweig, zobsmask )
590         
591
592         ! Interpolate the model SSH to the observation point
593         CALL obs_int_h2d( 1, 1,      &
594            &              zweig, zsshl(:,:,iobs),  zext )
595         
596         sladatqc%rext(jobs,1) = zext(1)
597         ! ... Remove the MDT at the observation point
598         sladatqc%rmod(jobs,1) = sladatqc%rext(jobs,1) - sladatqc%rext(jobs,2)
599
600      END DO
601
602      ! Deallocate the data for interpolation
603      DEALLOCATE( &
604         & igrdi, &
605         & igrdj, &
606         & zglam, &
607         & zgphi, &
608         & zmask, &
609         & zsshl  &
610         & )
611
612      sladatqc%nsurfup = sladatqc%nsurfup + isla
613
614   END SUBROUTINE obs_sla_opt
615
616   SUBROUTINE obs_sst_opt( sstdatqc, kt, kpi, kpj, kit000, kdaystp, &
617      &                    psstn, psstmask, k2dint, ld_nightav )
618      !!-----------------------------------------------------------------------
619      !!
620      !!                     ***  ROUTINE obs_sst_opt  ***
621      !!
622      !! ** Purpose : Compute the model counterpart of surface temperature
623      !!              data by interpolating from the model grid to the
624      !!              observation point.
625      !!
626      !! ** Method  : Linearly interpolate to each observation point using
627      !!              the model values at the corners of the surrounding grid box.
628      !!
629      !!    The now model SST is first computed at the obs (lon, lat) point.
630      !!
631      !!    Several horizontal interpolation schemes are available:
632      !!        - distance-weighted (great circle) (k2dint = 0)
633      !!        - distance-weighted (small angle)  (k2dint = 1)
634      !!        - bilinear (geographical grid)     (k2dint = 2)
635      !!        - bilinear (quadrilateral grid)    (k2dint = 3)
636      !!        - polynomial (quadrilateral grid)  (k2dint = 4)
637      !!
638      !!
639      !! ** Action  :
640      !!
641      !! History :
642      !!        !  07-07  (S. Ricci ) : Original
643      !!     
644      !!-----------------------------------------------------------------------
645
646      !! * Modules used
647      USE obs_surf_def  ! Definition of storage space for surface observations
648      USE sbcdcy
649
650      IMPLICIT NONE
651
652      !! * Arguments
653      TYPE(obs_surf), INTENT(INOUT) :: &
654         & sstdatqc     ! Subset of surface data not failing screening
655      INTEGER, INTENT(IN) :: kt        ! Time step
656      INTEGER, INTENT(IN) :: kpi       ! Model grid parameters
657      INTEGER, INTENT(IN) :: kpj
658      INTEGER, INTENT(IN) :: kit000    ! Number of the first time step
659                                       !   (kit000-1 = restart time)
660      INTEGER, INTENT(IN) :: k2dint    ! Horizontal interpolation type (see header)
661      INTEGER, INTENT(IN) :: kdaystp   ! Number of time steps per day 
662      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
663         & psstn,  &    ! Model SST field
664         & psstmask     ! Land-sea mask
665
666      !! * Local declarations
667      INTEGER :: ji
668      INTEGER :: jj
669      INTEGER :: jobs
670      INTEGER :: inrc
671      INTEGER :: isst
672      INTEGER :: iobs
673      INTEGER :: idayend
674      REAL(KIND=wp) :: zlam
675      REAL(KIND=wp) :: zphi
676      REAL(KIND=wp) :: zext(1), zobsmask(1)
677      REAL(KIND=wp) :: zdaystp
678      INTEGER, DIMENSION(:,:), SAVE, ALLOCATABLE :: &
679         & icount_sstnight,      &
680         & imask_night
681      REAL(kind=wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: &
682         & zintmp, &
683         & zouttmp, & 
684         & zmeanday    ! to compute model sst in region of 24h daylight (pole)
685      REAL(kind=wp), DIMENSION(2,2,1) :: &
686         & zweig
687      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
688         & zmask, &
689         & zsstl, &
690         & zsstm, &
691         & zglam, &
692         & zgphi
693      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
694         & igrdi, &
695         & igrdj
696      LOGICAL, INTENT(IN) :: ld_nightav
697
698      !-----------------------------------------------------------------------
699      ! Local initialization
700      !-----------------------------------------------------------------------
701      ! ... Record and data counters
702      inrc = kt - kit000 + 2
703      isst = sstdatqc%nsstp(inrc)
704
705      IF ( ld_nightav ) THEN
706
707      ! Initialize array for night mean
708
709      IF ( kt .EQ. 0 ) THEN
710         ALLOCATE ( icount_sstnight(kpi,kpj) )
711         ALLOCATE ( imask_night(kpi,kpj) )
712         ALLOCATE ( zintmp(kpi,kpj) )
713         ALLOCATE ( zouttmp(kpi,kpj) )
714         ALLOCATE ( zmeanday(kpi,kpj) )
715         nday_qsr = -1   ! initialisation flag for nbc_dcy
716      ENDIF
717
718      ! Initialize daily mean for first timestep
719      idayend = MOD( kt - kit000 + 1, kdaystp )
720
721      ! Added kt == 0 test to catch restart case
722      IF ( idayend == 1 .OR. kt == 0) THEN
723         IF (lwp) WRITE(numout,*) 'Reset sstdatqc%vdmean on time-step: ',kt
724         DO jj = 1, jpj
725            DO ji = 1, jpi
726               sstdatqc%vdmean(ji,jj) = 0.0
727               zmeanday(ji,jj) = 0.0
728               icount_sstnight(ji,jj) = 0
729            END DO
730         END DO
731      ENDIF
732
733      zintmp(:,:) = 0.0
734      zouttmp(:,:) = sbc_dcy( zintmp(:,:), .TRUE. )
735      imask_night(:,:) = INT( zouttmp(:,:) )
736
737      DO jj = 1, jpj
738         DO ji = 1, jpi
739            ! Increment the temperature field for computing night mean and counter
740            sstdatqc%vdmean(ji,jj) = sstdatqc%vdmean(ji,jj)  &
741                   &                        + psstn(ji,jj)*imask_night(ji,jj)
742            zmeanday(ji,jj)        = zmeanday(ji,jj) + psstn(ji,jj)
743            icount_sstnight(ji,jj) = icount_sstnight(ji,jj) + imask_night(ji,jj)
744         END DO
745      END DO
746   
747      ! Compute the daily mean at the end of day
748
749      zdaystp = 1.0 / REAL( kdaystp )
750
751      IF ( idayend == 0 ) THEN
752         DO jj = 1, jpj
753            DO ji = 1, jpi
754               ! Test if "no night" point
755               IF ( icount_sstnight(ji,jj) .NE. 0 ) THEN
756                  sstdatqc%vdmean(ji,jj) = sstdatqc%vdmean(ji,jj) &
757                    &                        / icount_sstnight(ji,jj) 
758               ELSE
759                  sstdatqc%vdmean(ji,jj) = zmeanday(ji,jj) * zdaystp
760               ENDIF
761            END DO
762         END DO
763      ENDIF
764
765      ENDIF
766
767      ! Get the data for interpolation
768     
769      ALLOCATE( &
770         & igrdi(2,2,isst), &
771         & igrdj(2,2,isst), &
772         & zglam(2,2,isst), &
773         & zgphi(2,2,isst), &
774         & zmask(2,2,isst), &
775         & zsstl(2,2,isst)  &
776         & )
777     
778      DO jobs = sstdatqc%nsurfup + 1, sstdatqc%nsurfup + isst
779         iobs = jobs - sstdatqc%nsurfup
780         igrdi(1,1,iobs) = sstdatqc%mi(jobs)-1
781         igrdj(1,1,iobs) = sstdatqc%mj(jobs)-1
782         igrdi(1,2,iobs) = sstdatqc%mi(jobs)-1
783         igrdj(1,2,iobs) = sstdatqc%mj(jobs)
784         igrdi(2,1,iobs) = sstdatqc%mi(jobs)
785         igrdj(2,1,iobs) = sstdatqc%mj(jobs)-1
786         igrdi(2,2,iobs) = sstdatqc%mi(jobs)
787         igrdj(2,2,iobs) = sstdatqc%mj(jobs)
788      END DO
789     
790      CALL obs_int_comm_2d( 2, 2, isst, &
791         &                  igrdi, igrdj, glamt, zglam )
792      CALL obs_int_comm_2d( 2, 2, isst, &
793         &                  igrdi, igrdj, gphit, zgphi )
794      CALL obs_int_comm_2d( 2, 2, isst, &
795         &                  igrdi, igrdj, psstmask, zmask )
796      CALL obs_int_comm_2d( 2, 2, isst, &
797         &                  igrdi, igrdj, psstn, zsstl )
798
799      ! At the end of the day get interpolated means
800      IF ( idayend == 0 .AND. ld_nightav ) THEN
801
802         ALLOCATE( &
803            & zsstm(2,2,isst)  &
804            & )
805
806         CALL obs_int_comm_2d( 2, 2, isst, igrdi, igrdj, &
807            &               sstdatqc%vdmean(:,:), zsstm )
808
809      ENDIF
810
811      ! Loop over observations
812
813      DO jobs = sstdatqc%nsurfup + 1, sstdatqc%nsurfup + isst
814         
815         iobs = jobs - sstdatqc%nsurfup
816         
817         IF ( kt /= sstdatqc%mstp(jobs) ) THEN
818           
819            IF(lwp) THEN
820               WRITE(numout,*)
821               WRITE(numout,*) ' E R R O R : Observation',              &
822                  &            ' time step is not consistent with the', &
823                  &            ' model time step'
824               WRITE(numout,*) ' ========='
825               WRITE(numout,*)
826               WRITE(numout,*) ' Record  = ', jobs,                &
827                  &            ' kt      = ', kt,                  &
828                  &            ' mstp    = ', sstdatqc%mstp(jobs), &
829                  &            ' ntyp    = ', sstdatqc%ntyp(jobs)
830            ENDIF
831            CALL ctl_stop( 'obs_sst_opt', 'Inconsistent time' )
832           
833         ENDIF
834         
835         zlam = sstdatqc%rlam(jobs)
836         zphi = sstdatqc%rphi(jobs)
837         
838         ! Get weights to interpolate the model SST to the observation point
839         CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         &
840            &                   zglam(:,:,iobs), zgphi(:,:,iobs), &
841            &                   zmask(:,:,iobs), zweig, zobsmask )
842           
843         ! Interpolate the model SST to the observation point
844
845         IF ( ld_nightav ) THEN
846
847           IF ( idayend == 0 )  THEN
848               ! Daily averaged/diurnal cycle of SST  data
849               CALL obs_int_h2d( 1, 1,      & 
850                     &              zweig, zsstm(:,:,iobs), zext )
851            ELSE
852               CALL ctl_stop( ' ld_nightav is set to true: a nonzero' //     &
853                     &           ' number of night SST data should' // &
854                     &           ' only occur at the end of a given day' )
855            ENDIF
856
857         ELSE
858
859            CALL obs_int_h2d( 1, 1,      &
860            &              zweig, zsstl(:,:,iobs),  zext )
861
862         ENDIF
863         sstdatqc%rmod(jobs,1) = zext(1)
864         
865      END DO
866     
867      ! Deallocate the data for interpolation
868      DEALLOCATE( &
869         & igrdi, &
870         & igrdj, &
871         & zglam, &
872         & zgphi, &
873         & zmask, &
874         & zsstl  &
875         & )
876
877      ! At the end of the day also get interpolated means
878      IF ( idayend == 0 .AND. ld_nightav ) THEN
879         DEALLOCATE( &
880            & zsstm  &
881            & )
882      ENDIF
883     
884      sstdatqc%nsurfup = sstdatqc%nsurfup + isst
885
886   END SUBROUTINE obs_sst_opt
887
888   SUBROUTINE obs_sss_opt
889      !!-----------------------------------------------------------------------
890      !!
891      !!                     ***  ROUTINE obs_sss_opt  ***
892      !!
893      !! ** Purpose : Compute the model counterpart of sea surface salinity
894      !!              data by interpolating from the model grid to the
895      !!              observation point.
896      !!
897      !! ** Method  :
898      !!
899      !! ** Action  :
900      !!
901      !! History :
902      !!      ! ??-??
903      !!-----------------------------------------------------------------------
904
905      IMPLICIT NONE
906
907   END SUBROUTINE obs_sss_opt
908
909   SUBROUTINE obs_seaice_opt( seaicedatqc, kt, kpi, kpj, kit000, &
910      &                    pseaicen, pseaicemask, k2dint )
911
912      !!-----------------------------------------------------------------------
913      !!
914      !!                     ***  ROUTINE obs_seaice_opt  ***
915      !!
916      !! ** Purpose : Compute the model counterpart of surface temperature
917      !!              data by interpolating from the model grid to the
918      !!              observation point.
919      !!
920      !! ** Method  : Linearly interpolate to each observation point using
921      !!              the model values at the corners of the surrounding grid box.
922      !!
923      !!    The now model sea ice is first computed at the obs (lon, lat) point.
924      !!
925      !!    Several horizontal interpolation schemes are available:
926      !!        - distance-weighted (great circle) (k2dint = 0)
927      !!        - distance-weighted (small angle)  (k2dint = 1)
928      !!        - bilinear (geographical grid)     (k2dint = 2)
929      !!        - bilinear (quadrilateral grid)    (k2dint = 3)
930      !!        - polynomial (quadrilateral grid)  (k2dint = 4)
931      !!
932      !!
933      !! ** Action  :
934      !!
935      !! History :
936      !!        !  07-07  (S. Ricci ) : Original
937      !!     
938      !!-----------------------------------------------------------------------
939
940      !! * Modules used
941      USE obs_surf_def  ! Definition of storage space for surface observations
942
943      IMPLICIT NONE
944
945      !! * Arguments
946      TYPE(obs_surf), INTENT(INOUT) :: seaicedatqc     ! Subset of surface data not failing screening
947      INTEGER, INTENT(IN) :: kt       ! Time step
948      INTEGER, INTENT(IN) :: kpi      ! Model grid parameters
949      INTEGER, INTENT(IN) :: kpj
950      INTEGER, INTENT(IN) :: kit000   ! Number of the first time step
951                                      !   (kit000-1 = restart time)
952      INTEGER, INTENT(IN) :: k2dint   ! Horizontal interpolation type (see header)
953      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
954         & pseaicen,  &    ! Model sea ice field
955         & pseaicemask     ! Land-sea mask
956         
957      !! * Local declarations
958      INTEGER :: ji
959      INTEGER :: jj
960      INTEGER :: jobs
961      INTEGER :: inrc
962      INTEGER :: iseaice
963      INTEGER :: iobs
964       
965      REAL(KIND=wp) :: zlam
966      REAL(KIND=wp) :: zphi
967      REAL(KIND=wp) :: zext(1), zobsmask(1)
968      REAL(kind=wp), DIMENSION(2,2,1) :: &
969         & zweig
970      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
971         & zmask, &
972         & zseaicel, &
973         & zglam, &
974         & zgphi
975      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
976         & igrdi, &
977         & igrdj
978
979      !------------------------------------------------------------------------
980      ! Local initialization
981      !------------------------------------------------------------------------
982      ! ... Record and data counters
983      inrc = kt - kit000 + 2
984      iseaice = seaicedatqc%nsstp(inrc)
985
986      ! Get the data for interpolation
987     
988      ALLOCATE( &
989         & igrdi(2,2,iseaice), &
990         & igrdj(2,2,iseaice), &
991         & zglam(2,2,iseaice), &
992         & zgphi(2,2,iseaice), &
993         & zmask(2,2,iseaice), &
994         & zseaicel(2,2,iseaice)  &
995         & )
996     
997      DO jobs = seaicedatqc%nsurfup + 1, seaicedatqc%nsurfup + iseaice
998         iobs = jobs - seaicedatqc%nsurfup
999         igrdi(1,1,iobs) = seaicedatqc%mi(jobs)-1
1000         igrdj(1,1,iobs) = seaicedatqc%mj(jobs)-1
1001         igrdi(1,2,iobs) = seaicedatqc%mi(jobs)-1
1002         igrdj(1,2,iobs) = seaicedatqc%mj(jobs)
1003         igrdi(2,1,iobs) = seaicedatqc%mi(jobs)
1004         igrdj(2,1,iobs) = seaicedatqc%mj(jobs)-1
1005         igrdi(2,2,iobs) = seaicedatqc%mi(jobs)
1006         igrdj(2,2,iobs) = seaicedatqc%mj(jobs)
1007      END DO
1008     
1009      CALL obs_int_comm_2d( 2, 2, iseaice, &
1010         &                  igrdi, igrdj, glamt, zglam )
1011      CALL obs_int_comm_2d( 2, 2, iseaice, &
1012         &                  igrdi, igrdj, gphit, zgphi )
1013      CALL obs_int_comm_2d( 2, 2, iseaice, &
1014         &                  igrdi, igrdj, pseaicemask, zmask )
1015      CALL obs_int_comm_2d( 2, 2, iseaice, &
1016         &                  igrdi, igrdj, pseaicen, zseaicel )
1017     
1018      DO jobs = seaicedatqc%nsurfup + 1, seaicedatqc%nsurfup + iseaice
1019         
1020         iobs = jobs - seaicedatqc%nsurfup
1021         
1022         IF ( kt /= seaicedatqc%mstp(jobs) ) THEN
1023           
1024            IF(lwp) THEN
1025               WRITE(numout,*)
1026               WRITE(numout,*) ' E R R O R : Observation',              &
1027                  &            ' time step is not consistent with the', &
1028                  &            ' model time step'
1029               WRITE(numout,*) ' ========='
1030               WRITE(numout,*)
1031               WRITE(numout,*) ' Record  = ', jobs,                &
1032                  &            ' kt      = ', kt,                  &
1033                  &            ' mstp    = ', seaicedatqc%mstp(jobs), &
1034                  &            ' ntyp    = ', seaicedatqc%ntyp(jobs)
1035            ENDIF
1036            CALL ctl_stop( 'obs_seaice_opt', 'Inconsistent time' )
1037           
1038         ENDIF
1039         
1040         zlam = seaicedatqc%rlam(jobs)
1041         zphi = seaicedatqc%rphi(jobs)
1042         
1043         ! Get weights to interpolate the model sea ice to the observation point
1044         CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         &
1045            &                   zglam(:,:,iobs), zgphi(:,:,iobs), &
1046            &                   zmask(:,:,iobs), zweig, zobsmask )
1047         
1048         ! ... Interpolate the model sea ice to the observation point
1049         CALL obs_int_h2d( 1, 1,      &
1050            &              zweig, zseaicel(:,:,iobs),  zext )
1051         
1052         seaicedatqc%rmod(jobs,1) = zext(1)
1053         
1054      END DO
1055     
1056      ! Deallocate the data for interpolation
1057      DEALLOCATE( &
1058         & igrdi,    &
1059         & igrdj,    &
1060         & zglam,    &
1061         & zgphi,    &
1062         & zmask,    &
1063         & zseaicel  &
1064         & )
1065     
1066      seaicedatqc%nsurfup = seaicedatqc%nsurfup + iseaice
1067
1068   END SUBROUTINE obs_seaice_opt
1069
1070   SUBROUTINE obs_vel_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, &
1071      &                    pun, pvn, pgdept, pumask, pvmask, k1dint, k2dint, &
1072      &                    ld_dailyav )
1073      !!-----------------------------------------------------------------------
1074      !!
1075      !!                     ***  ROUTINE obs_vel_opt  ***
1076      !!
1077      !! ** Purpose : Compute the model counterpart of velocity profile
1078      !!              data by interpolating from the model grid to the
1079      !!              observation point.
1080      !!
1081      !! ** Method  : Linearly interpolate zonal and meridional components of velocity
1082      !!              to each observation point using the model values at the corners of
1083      !!              the surrounding grid box. The model velocity components are on a
1084      !!              staggered C- grid.
1085      !!
1086      !!    For velocity data from the TAO array, the model equivalent is
1087      !!    a daily mean velocity field. So, we first compute
1088      !!    the mean, then interpolate only at the end of the day.
1089      !!
1090      !! ** Action  :
1091      !!
1092      !! History :
1093      !!    ! 07-03 (K. Mogensen)      : Temperature and Salinity profiles
1094      !!    ! 08-10 (Maria Valdivieso) : Velocity component (U,V) profiles
1095      !!-----------------------------------------------------------------------
1096   
1097      !! * Modules used
1098      USE obs_profiles_def ! Definition of storage space for profile obs.
1099
1100      IMPLICIT NONE
1101
1102      !! * Arguments
1103      TYPE(obs_prof), INTENT(INOUT) :: &
1104         & prodatqc        ! Subset of profile data not failing screening
1105      INTEGER, INTENT(IN) :: kt        ! Time step
1106      INTEGER, INTENT(IN) :: kpi       ! Model grid parameters
1107      INTEGER, INTENT(IN) :: kpj
1108      INTEGER, INTENT(IN) :: kpk 
1109      INTEGER, INTENT(IN) :: kit000    ! Number of the first time step
1110                                       !   (kit000-1 = restart time)
1111      INTEGER, INTENT(IN) :: k1dint    ! Vertical interpolation type (see header)
1112      INTEGER, INTENT(IN) :: k2dint    ! Horizontal interpolation type (see header)
1113      INTEGER, INTENT(IN) :: kdaystp   ! Number of time steps per day                   
1114      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &
1115         & pun,    &    ! Model zonal component of velocity
1116         & pvn,    &    ! Model meridional component of velocity
1117         & pumask, &    ! Land-sea mask
1118         & pvmask       ! Land-sea mask
1119      REAL(KIND=wp), INTENT(IN), DIMENSION(kpk) :: &
1120         & pgdept       ! Model array of depth levels
1121      LOGICAL, INTENT(IN) :: ld_dailyav
1122         
1123      !! * Local declarations
1124      INTEGER :: ji
1125      INTEGER :: jj
1126      INTEGER :: jk
1127      INTEGER :: jobs
1128      INTEGER :: inrc
1129      INTEGER :: ipro
1130      INTEGER :: idayend
1131      INTEGER :: ista
1132      INTEGER :: iend
1133      INTEGER :: iobs
1134      INTEGER, DIMENSION(imaxavtypes) :: &
1135         & idailyavtypes
1136      REAL(KIND=wp) :: zlam
1137      REAL(KIND=wp) :: zphi
1138      REAL(KIND=wp) :: zdaystp
1139      REAL(KIND=wp), DIMENSION(kpk) :: &
1140         & zobsmasku, &
1141         & zobsmaskv, &
1142         & zobsmask,  &
1143         & zobsk,     &
1144         & zobs2k
1145      REAL(KIND=wp), DIMENSION(2,2,kpk) :: &
1146         & zweigu,zweigv
1147      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
1148         & zumask, zvmask, &
1149         & zintu, &
1150         & zintv, &
1151         & zinmu, &
1152         & zinmv
1153      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
1154         & zglamu, zglamv, &
1155         & zgphiu, zgphiv
1156      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
1157         & igrdiu, &
1158         & igrdju, &
1159         & igrdiv, &
1160         & igrdjv
1161
1162      !------------------------------------------------------------------------
1163      ! Local initialization
1164      !------------------------------------------------------------------------
1165      ! ... Record and data counters
1166      inrc = kt - kit000 + 2
1167      ipro = prodatqc%npstp(inrc)
1168
1169      ! Initialize daily mean for first timestep
1170      idayend = MOD( kt - kit000 + 1, kdaystp )
1171
1172      ! Added kt == 0 test to catch restart case
1173      IF ( idayend == 1 .OR. kt == 0) THEN
1174         IF (lwp) WRITE(numout,*) 'Reset prodatqc%vdmean on time-step: ',kt
1175         prodatqc%vdmean(:,:,:,1) = 0.0
1176         prodatqc%vdmean(:,:,:,2) = 0.0
1177      ENDIF
1178
1179      ! Increment the zonal velocity field for computing daily mean
1180      prodatqc%vdmean(:,:,:,1) = prodatqc%vdmean(:,:,:,1) + pun(:,:,:)
1181      ! Increment the meridional velocity field for computing daily mean
1182      prodatqc%vdmean(:,:,:,2) = prodatqc%vdmean(:,:,:,2) + pvn(:,:,:)
1183   
1184      ! Compute the daily mean at the end of day
1185      zdaystp = 1.0 / REAL( kdaystp )
1186      IF ( idayend == 0 ) THEN
1187         prodatqc%vdmean(:,:,:,1) = prodatqc%vdmean(:,:,:,1) * zdaystp
1188         prodatqc%vdmean(:,:,:,2) = prodatqc%vdmean(:,:,:,2) * zdaystp
1189      ENDIF
1190
1191      ! Get the data for interpolation
1192      ALLOCATE( &
1193         & igrdiu(2,2,ipro),      &
1194         & igrdju(2,2,ipro),      &
1195         & igrdiv(2,2,ipro),      &
1196         & igrdjv(2,2,ipro),      &
1197         & zglamu(2,2,ipro), zglamv(2,2,ipro), &
1198         & zgphiu(2,2,ipro), zgphiv(2,2,ipro), &
1199         & zumask(2,2,kpk,ipro), zvmask(2,2,kpk,ipro), &
1200         & zintu(2,2,kpk,ipro),  &
1201         & zintv(2,2,kpk,ipro)   &
1202         & )
1203
1204      DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
1205         iobs = jobs - prodatqc%nprofup
1206         igrdiu(1,1,iobs) = prodatqc%mi(jobs,1)-1
1207         igrdju(1,1,iobs) = prodatqc%mj(jobs,1)-1
1208         igrdiu(1,2,iobs) = prodatqc%mi(jobs,1)-1
1209         igrdju(1,2,iobs) = prodatqc%mj(jobs,1)
1210         igrdiu(2,1,iobs) = prodatqc%mi(jobs,1)
1211         igrdju(2,1,iobs) = prodatqc%mj(jobs,1)-1
1212         igrdiu(2,2,iobs) = prodatqc%mi(jobs,1)
1213         igrdju(2,2,iobs) = prodatqc%mj(jobs,1)
1214         igrdiv(1,1,iobs) = prodatqc%mi(jobs,2)-1
1215         igrdjv(1,1,iobs) = prodatqc%mj(jobs,2)-1
1216         igrdiv(1,2,iobs) = prodatqc%mi(jobs,2)-1
1217         igrdjv(1,2,iobs) = prodatqc%mj(jobs,2)
1218         igrdiv(2,1,iobs) = prodatqc%mi(jobs,2)
1219         igrdjv(2,1,iobs) = prodatqc%mj(jobs,2)-1
1220         igrdiv(2,2,iobs) = prodatqc%mi(jobs,2)
1221         igrdjv(2,2,iobs) = prodatqc%mj(jobs,2)
1222      END DO
1223
1224      CALL obs_int_comm_2d( 2, 2, ipro, igrdiu, igrdju, glamu, zglamu )
1225      CALL obs_int_comm_2d( 2, 2, ipro, igrdiu, igrdju, gphiu, zgphiu )
1226      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, pumask, zumask )
1227      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, pun, zintu )
1228
1229      CALL obs_int_comm_2d( 2, 2, ipro, igrdiv, igrdjv, glamv, zglamv )
1230      CALL obs_int_comm_2d( 2, 2, ipro, igrdiv, igrdjv, gphiv, zgphiv )
1231      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, pvmask, zvmask )
1232      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, pvn, zintv )
1233
1234      ! At the end of the day also get interpolated means
1235      IF ( idayend == 0 ) THEN
1236
1237         ALLOCATE( &
1238            & zinmu(2,2,kpk,ipro),  &
1239            & zinmv(2,2,kpk,ipro)   &
1240            & )
1241
1242         CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, &
1243            &                  prodatqc%vdmean(:,:,:,1), zinmu )
1244         CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, &
1245            &                  prodatqc%vdmean(:,:,:,2), zinmv )
1246
1247      ENDIF
1248
1249! loop over observations
1250
1251      DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
1252
1253         iobs = jobs - prodatqc%nprofup
1254
1255         IF ( kt /= prodatqc%mstp(jobs) ) THEN
1256           
1257            IF(lwp) THEN
1258               WRITE(numout,*)
1259               WRITE(numout,*) ' E R R O R : Observation',              &
1260                  &            ' time step is not consistent with the', &
1261                  &            ' model time step'
1262               WRITE(numout,*) ' ========='
1263               WRITE(numout,*)
1264               WRITE(numout,*) ' Record  = ', jobs,                    &
1265                  &            ' kt      = ', kt,                      &
1266                  &            ' mstp    = ', prodatqc%mstp(jobs), &
1267                  &            ' ntyp    = ', prodatqc%ntyp(jobs)
1268            ENDIF
1269            CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' )
1270         ENDIF
1271         
1272         zlam = prodatqc%rlam(jobs)
1273         zphi = prodatqc%rphi(jobs)
1274
1275         ! Initialize observation masks
1276
1277         zobsmasku(:) = 0.0
1278         zobsmaskv(:) = 0.0
1279         
1280         ! Horizontal weights and vertical mask
1281
1282         IF  ( prodatqc%npvend(jobs,1) > 0 ) THEN
1283
1284            CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi,     &
1285               &                   zglamu(:,:,iobs), zgphiu(:,:,iobs), &
1286               &                   zumask(:,:,:,iobs), zweigu, zobsmasku )
1287
1288         ENDIF
1289
1290         
1291         IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
1292
1293            CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi,     &
1294               &                   zglamv(:,:,iobs), zgphiv(:,:,iobs), &
1295               &                   zvmask(:,:,:,iobs), zweigv, zobsmasku )
1296
1297         ENDIF
1298
1299         ! Ensure that the vertical mask on u and v are consistent.
1300
1301         zobsmask(:) = MIN( zobsmasku(:), zobsmaskv(:) )
1302
1303         IF ( prodatqc%npvend(jobs,1) > 0 ) THEN
1304
1305            zobsk(:) = obfillflt
1306
1307       IF ( ld_dailyav ) THEN
1308
1309               IF ( idayend == 0 )  THEN
1310                 
1311                  ! Daily averaged data
1312                 
1313                  CALL obs_int_h2d( kpk, kpk,      &
1314                     &              zweigu, zinmu(:,:,:,iobs), zobsk )
1315                 
1316                 
1317               ELSE
1318               
1319                  CALL ctl_stop( ' A nonzero' //     &
1320                     &           ' number of U profile data should' // &
1321                     &           ' only occur at the end of a given day' )
1322
1323               ENDIF
1324         
1325            ELSE 
1326               
1327               ! Point data
1328
1329               CALL obs_int_h2d( kpk, kpk,      &
1330                  &              zweigu, zintu(:,:,:,iobs), zobsk )
1331
1332            ENDIF
1333
1334            !-------------------------------------------------------------
1335            ! Compute vertical second-derivative of the interpolating
1336            ! polynomial at obs points
1337            !-------------------------------------------------------------
1338           
1339            IF ( k1dint == 1 ) THEN
1340               CALL obs_int_z1d_spl( kpk, zobsk, zobs2k,   &
1341                  &                  pgdept, zobsmask )
1342            ENDIF
1343           
1344            !-----------------------------------------------------------------
1345            !  Vertical interpolation to the observation point
1346            !-----------------------------------------------------------------
1347            ista = prodatqc%npvsta(jobs,1)
1348            iend = prodatqc%npvend(jobs,1)
1349            CALL obs_int_z1d( kpk,                &
1350               & prodatqc%var(1)%mvk(ista:iend),  &
1351               & k1dint, iend - ista + 1,         &
1352               & prodatqc%var(1)%vdep(ista:iend), &
1353               & zobsk, zobs2k,                   &
1354               & prodatqc%var(1)%vmod(ista:iend), &
1355               & pgdept, zobsmask )
1356
1357         ENDIF
1358
1359         IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
1360
1361            zobsk(:) = obfillflt
1362
1363            IF ( ld_dailyav ) THEN
1364
1365               IF ( idayend == 0 )  THEN
1366
1367                  ! Daily averaged data
1368                 
1369                  CALL obs_int_h2d( kpk, kpk,      &
1370                     &              zweigv, zinmv(:,:,:,iobs), zobsk )
1371                 
1372               ELSE
1373
1374                  CALL ctl_stop( ' A nonzero' //     &
1375                     &           ' number of V profile data should' // &
1376                     &           ' only occur at the end of a given day' )
1377
1378               ENDIF
1379
1380            ELSE
1381               
1382               ! Point data
1383
1384               CALL obs_int_h2d( kpk, kpk,      &
1385                  &              zweigv, zintv(:,:,:,iobs), zobsk )
1386
1387            ENDIF
1388
1389
1390            !-------------------------------------------------------------
1391            ! Compute vertical second-derivative of the interpolating
1392            ! polynomial at obs points
1393            !-------------------------------------------------------------
1394           
1395            IF ( k1dint == 1 ) THEN
1396               CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, &
1397                  &                  pgdept, zobsmask )
1398            ENDIF
1399           
1400            !----------------------------------------------------------------
1401            !  Vertical interpolation to the observation point
1402            !----------------------------------------------------------------
1403            ista = prodatqc%npvsta(jobs,2)
1404            iend = prodatqc%npvend(jobs,2)
1405            CALL obs_int_z1d( kpk, &
1406               & prodatqc%var(2)%mvk(ista:iend),&
1407               & k1dint, iend - ista + 1, &
1408               & prodatqc%var(2)%vdep(ista:iend),&
1409               & zobsk, zobs2k, &
1410               & prodatqc%var(2)%vmod(ista:iend),&
1411               & pgdept, zobsmask )
1412
1413         ENDIF
1414
1415      END DO
1416 
1417      ! Deallocate the data for interpolation
1418      DEALLOCATE( &
1419         & igrdiu, &
1420         & igrdju, &
1421         & igrdiv, &
1422         & igrdjv, &
1423         & zglamu, zglamv, &
1424         & zgphiu, zgphiv, &
1425         & zumask, zvmask, &
1426         & zintu, &
1427         & zintv  &
1428         & )
1429      ! At the end of the day also get interpolated means
1430      IF ( idayend == 0 ) THEN
1431         DEALLOCATE( &
1432            & zinmu,  &
1433            & zinmv   &
1434            & )
1435      ENDIF
1436
1437      prodatqc%nprofup = prodatqc%nprofup + ipro 
1438     
1439   END SUBROUTINE obs_vel_opt
1440
1441END MODULE obs_oper
1442
Note: See TracBrowser for help on using the repository browser.