source: branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90 @ 3586

Last change on this file since 3586 was 3586, checked in by cbricaud, 8 years ago

add modification from dev_r3342_MERCATOR7_SST in dev_MERCATOR_2012_rev3555

  • Property svn:keywords set to Id
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         
864         sstdatqc%rmod(jobs,1) = zext(1)
865         
866      END DO
867     
868      ! Deallocate the data for interpolation
869      DEALLOCATE( &
870         & igrdi, &
871         & igrdj, &
872         & zglam, &
873         & zgphi, &
874         & zmask, &
875         & zsstl  &
876         & )
877
878      ! At the end of the day also get interpolated means
879      IF ( idayend == 0 .AND. ld_nightav ) THEN
880         DEALLOCATE( &
881            & zsstm  &
882            & )
883      ENDIF
884     
885      sstdatqc%nsurfup = sstdatqc%nsurfup + isst
886
887   END SUBROUTINE obs_sst_opt
888
889   SUBROUTINE obs_sss_opt
890      !!-----------------------------------------------------------------------
891      !!
892      !!                     ***  ROUTINE obs_sss_opt  ***
893      !!
894      !! ** Purpose : Compute the model counterpart of sea surface salinity
895      !!              data by interpolating from the model grid to the
896      !!              observation point.
897      !!
898      !! ** Method  :
899      !!
900      !! ** Action  :
901      !!
902      !! History :
903      !!      ! ??-??
904      !!-----------------------------------------------------------------------
905
906      IMPLICIT NONE
907
908   END SUBROUTINE obs_sss_opt
909
910   SUBROUTINE obs_seaice_opt( seaicedatqc, kt, kpi, kpj, kit000, &
911      &                    pseaicen, pseaicemask, k2dint )
912
913      !!-----------------------------------------------------------------------
914      !!
915      !!                     ***  ROUTINE obs_seaice_opt  ***
916      !!
917      !! ** Purpose : Compute the model counterpart of surface temperature
918      !!              data by interpolating from the model grid to the
919      !!              observation point.
920      !!
921      !! ** Method  : Linearly interpolate to each observation point using
922      !!              the model values at the corners of the surrounding grid box.
923      !!
924      !!    The now model sea ice is first computed at the obs (lon, lat) point.
925      !!
926      !!    Several horizontal interpolation schemes are available:
927      !!        - distance-weighted (great circle) (k2dint = 0)
928      !!        - distance-weighted (small angle)  (k2dint = 1)
929      !!        - bilinear (geographical grid)     (k2dint = 2)
930      !!        - bilinear (quadrilateral grid)    (k2dint = 3)
931      !!        - polynomial (quadrilateral grid)  (k2dint = 4)
932      !!
933      !!
934      !! ** Action  :
935      !!
936      !! History :
937      !!        !  07-07  (S. Ricci ) : Original
938      !!     
939      !!-----------------------------------------------------------------------
940
941      !! * Modules used
942      USE obs_surf_def  ! Definition of storage space for surface observations
943
944      IMPLICIT NONE
945
946      !! * Arguments
947      TYPE(obs_surf), INTENT(INOUT) :: seaicedatqc     ! Subset of surface data not failing screening
948      INTEGER, INTENT(IN) :: kt       ! Time step
949      INTEGER, INTENT(IN) :: kpi      ! Model grid parameters
950      INTEGER, INTENT(IN) :: kpj
951      INTEGER, INTENT(IN) :: kit000   ! Number of the first time step
952                                      !   (kit000-1 = restart time)
953      INTEGER, INTENT(IN) :: k2dint   ! Horizontal interpolation type (see header)
954      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
955         & pseaicen,  &    ! Model sea ice field
956         & pseaicemask     ! Land-sea mask
957         
958      !! * Local declarations
959      INTEGER :: ji
960      INTEGER :: jj
961      INTEGER :: jobs
962      INTEGER :: inrc
963      INTEGER :: iseaice
964      INTEGER :: iobs
965       
966      REAL(KIND=wp) :: zlam
967      REAL(KIND=wp) :: zphi
968      REAL(KIND=wp) :: zext(1), zobsmask(1)
969      REAL(kind=wp), DIMENSION(2,2,1) :: &
970         & zweig
971      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
972         & zmask, &
973         & zseaicel, &
974         & zglam, &
975         & zgphi
976      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
977         & igrdi, &
978         & igrdj
979
980      !------------------------------------------------------------------------
981      ! Local initialization
982      !------------------------------------------------------------------------
983      ! ... Record and data counters
984      inrc = kt - kit000 + 2
985      iseaice = seaicedatqc%nsstp(inrc)
986
987      ! Get the data for interpolation
988     
989      ALLOCATE( &
990         & igrdi(2,2,iseaice), &
991         & igrdj(2,2,iseaice), &
992         & zglam(2,2,iseaice), &
993         & zgphi(2,2,iseaice), &
994         & zmask(2,2,iseaice), &
995         & zseaicel(2,2,iseaice)  &
996         & )
997     
998      DO jobs = seaicedatqc%nsurfup + 1, seaicedatqc%nsurfup + iseaice
999         iobs = jobs - seaicedatqc%nsurfup
1000         igrdi(1,1,iobs) = seaicedatqc%mi(jobs)-1
1001         igrdj(1,1,iobs) = seaicedatqc%mj(jobs)-1
1002         igrdi(1,2,iobs) = seaicedatqc%mi(jobs)-1
1003         igrdj(1,2,iobs) = seaicedatqc%mj(jobs)
1004         igrdi(2,1,iobs) = seaicedatqc%mi(jobs)
1005         igrdj(2,1,iobs) = seaicedatqc%mj(jobs)-1
1006         igrdi(2,2,iobs) = seaicedatqc%mi(jobs)
1007         igrdj(2,2,iobs) = seaicedatqc%mj(jobs)
1008      END DO
1009     
1010      CALL obs_int_comm_2d( 2, 2, iseaice, &
1011         &                  igrdi, igrdj, glamt, zglam )
1012      CALL obs_int_comm_2d( 2, 2, iseaice, &
1013         &                  igrdi, igrdj, gphit, zgphi )
1014      CALL obs_int_comm_2d( 2, 2, iseaice, &
1015         &                  igrdi, igrdj, pseaicemask, zmask )
1016      CALL obs_int_comm_2d( 2, 2, iseaice, &
1017         &                  igrdi, igrdj, pseaicen, zseaicel )
1018     
1019      DO jobs = seaicedatqc%nsurfup + 1, seaicedatqc%nsurfup + iseaice
1020         
1021         iobs = jobs - seaicedatqc%nsurfup
1022         
1023         IF ( kt /= seaicedatqc%mstp(jobs) ) THEN
1024           
1025            IF(lwp) THEN
1026               WRITE(numout,*)
1027               WRITE(numout,*) ' E R R O R : Observation',              &
1028                  &            ' time step is not consistent with the', &
1029                  &            ' model time step'
1030               WRITE(numout,*) ' ========='
1031               WRITE(numout,*)
1032               WRITE(numout,*) ' Record  = ', jobs,                &
1033                  &            ' kt      = ', kt,                  &
1034                  &            ' mstp    = ', seaicedatqc%mstp(jobs), &
1035                  &            ' ntyp    = ', seaicedatqc%ntyp(jobs)
1036            ENDIF
1037            CALL ctl_stop( 'obs_seaice_opt', 'Inconsistent time' )
1038           
1039         ENDIF
1040         
1041         zlam = seaicedatqc%rlam(jobs)
1042         zphi = seaicedatqc%rphi(jobs)
1043         
1044         ! Get weights to interpolate the model sea ice to the observation point
1045         CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         &
1046            &                   zglam(:,:,iobs), zgphi(:,:,iobs), &
1047            &                   zmask(:,:,iobs), zweig, zobsmask )
1048         
1049         ! ... Interpolate the model sea ice to the observation point
1050         CALL obs_int_h2d( 1, 1,      &
1051            &              zweig, zseaicel(:,:,iobs),  zext )
1052         
1053         seaicedatqc%rmod(jobs,1) = zext(1)
1054         
1055      END DO
1056     
1057      ! Deallocate the data for interpolation
1058      DEALLOCATE( &
1059         & igrdi,    &
1060         & igrdj,    &
1061         & zglam,    &
1062         & zgphi,    &
1063         & zmask,    &
1064         & zseaicel  &
1065         & )
1066     
1067      seaicedatqc%nsurfup = seaicedatqc%nsurfup + iseaice
1068
1069   END SUBROUTINE obs_seaice_opt
1070
1071   SUBROUTINE obs_vel_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, &
1072      &                    pun, pvn, pgdept, pumask, pvmask, k1dint, k2dint, &
1073      &                    ld_dailyav )
1074      !!-----------------------------------------------------------------------
1075      !!
1076      !!                     ***  ROUTINE obs_vel_opt  ***
1077      !!
1078      !! ** Purpose : Compute the model counterpart of velocity profile
1079      !!              data by interpolating from the model grid to the
1080      !!              observation point.
1081      !!
1082      !! ** Method  : Linearly interpolate zonal and meridional components of velocity
1083      !!              to each observation point using the model values at the corners of
1084      !!              the surrounding grid box. The model velocity components are on a
1085      !!              staggered C- grid.
1086      !!
1087      !!    For velocity data from the TAO array, the model equivalent is
1088      !!    a daily mean velocity field. So, we first compute
1089      !!    the mean, then interpolate only at the end of the day.
1090      !!
1091      !! ** Action  :
1092      !!
1093      !! History :
1094      !!    ! 07-03 (K. Mogensen)      : Temperature and Salinity profiles
1095      !!    ! 08-10 (Maria Valdivieso) : Velocity component (U,V) profiles
1096      !!-----------------------------------------------------------------------
1097   
1098      !! * Modules used
1099      USE obs_profiles_def ! Definition of storage space for profile obs.
1100
1101      IMPLICIT NONE
1102
1103      !! * Arguments
1104      TYPE(obs_prof), INTENT(INOUT) :: &
1105         & prodatqc        ! Subset of profile data not failing screening
1106      INTEGER, INTENT(IN) :: kt        ! Time step
1107      INTEGER, INTENT(IN) :: kpi       ! Model grid parameters
1108      INTEGER, INTENT(IN) :: kpj
1109      INTEGER, INTENT(IN) :: kpk 
1110      INTEGER, INTENT(IN) :: kit000    ! Number of the first time step
1111                                       !   (kit000-1 = restart time)
1112      INTEGER, INTENT(IN) :: k1dint    ! Vertical interpolation type (see header)
1113      INTEGER, INTENT(IN) :: k2dint    ! Horizontal interpolation type (see header)
1114      INTEGER, INTENT(IN) :: kdaystp   ! Number of time steps per day                   
1115      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &
1116         & pun,    &    ! Model zonal component of velocity
1117         & pvn,    &    ! Model meridional component of velocity
1118         & pumask, &    ! Land-sea mask
1119         & pvmask       ! Land-sea mask
1120      REAL(KIND=wp), INTENT(IN), DIMENSION(kpk) :: &
1121         & pgdept       ! Model array of depth levels
1122      LOGICAL, INTENT(IN) :: ld_dailyav
1123         
1124      !! * Local declarations
1125      INTEGER :: ji
1126      INTEGER :: jj
1127      INTEGER :: jk
1128      INTEGER :: jobs
1129      INTEGER :: inrc
1130      INTEGER :: ipro
1131      INTEGER :: idayend
1132      INTEGER :: ista
1133      INTEGER :: iend
1134      INTEGER :: iobs
1135      INTEGER, DIMENSION(imaxavtypes) :: &
1136         & idailyavtypes
1137      REAL(KIND=wp) :: zlam
1138      REAL(KIND=wp) :: zphi
1139      REAL(KIND=wp) :: zdaystp
1140      REAL(KIND=wp), DIMENSION(kpk) :: &
1141         & zobsmasku, &
1142         & zobsmaskv, &
1143         & zobsmask,  &
1144         & zobsk,     &
1145         & zobs2k
1146      REAL(KIND=wp), DIMENSION(2,2,kpk) :: &
1147         & zweigu,zweigv
1148      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
1149         & zumask, zvmask, &
1150         & zintu, &
1151         & zintv, &
1152         & zinmu, &
1153         & zinmv
1154      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
1155         & zglamu, zglamv, &
1156         & zgphiu, zgphiv
1157      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
1158         & igrdiu, &
1159         & igrdju, &
1160         & igrdiv, &
1161         & igrdjv
1162
1163      !------------------------------------------------------------------------
1164      ! Local initialization
1165      !------------------------------------------------------------------------
1166      ! ... Record and data counters
1167      inrc = kt - kit000 + 2
1168      ipro = prodatqc%npstp(inrc)
1169
1170      ! Initialize daily mean for first timestep
1171      idayend = MOD( kt - kit000 + 1, kdaystp )
1172
1173      ! Added kt == 0 test to catch restart case
1174      IF ( idayend == 1 .OR. kt == 0) THEN
1175         IF (lwp) WRITE(numout,*) 'Reset prodatqc%vdmean on time-step: ',kt
1176         prodatqc%vdmean(:,:,:,1) = 0.0
1177         prodatqc%vdmean(:,:,:,2) = 0.0
1178      ENDIF
1179
1180      ! Increment the zonal velocity field for computing daily mean
1181      prodatqc%vdmean(:,:,:,1) = prodatqc%vdmean(:,:,:,1) + pun(:,:,:)
1182      ! Increment the meridional velocity field for computing daily mean
1183      prodatqc%vdmean(:,:,:,2) = prodatqc%vdmean(:,:,:,2) + pvn(:,:,:)
1184   
1185      ! Compute the daily mean at the end of day
1186      zdaystp = 1.0 / REAL( kdaystp )
1187      IF ( idayend == 0 ) THEN
1188         prodatqc%vdmean(:,:,:,1) = prodatqc%vdmean(:,:,:,1) * zdaystp
1189         prodatqc%vdmean(:,:,:,2) = prodatqc%vdmean(:,:,:,2) * zdaystp
1190      ENDIF
1191
1192      ! Get the data for interpolation
1193      ALLOCATE( &
1194         & igrdiu(2,2,ipro),      &
1195         & igrdju(2,2,ipro),      &
1196         & igrdiv(2,2,ipro),      &
1197         & igrdjv(2,2,ipro),      &
1198         & zglamu(2,2,ipro), zglamv(2,2,ipro), &
1199         & zgphiu(2,2,ipro), zgphiv(2,2,ipro), &
1200         & zumask(2,2,kpk,ipro), zvmask(2,2,kpk,ipro), &
1201         & zintu(2,2,kpk,ipro),  &
1202         & zintv(2,2,kpk,ipro)   &
1203         & )
1204
1205      DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
1206         iobs = jobs - prodatqc%nprofup
1207         igrdiu(1,1,iobs) = prodatqc%mi(jobs,1)-1
1208         igrdju(1,1,iobs) = prodatqc%mj(jobs,1)-1
1209         igrdiu(1,2,iobs) = prodatqc%mi(jobs,1)-1
1210         igrdju(1,2,iobs) = prodatqc%mj(jobs,1)
1211         igrdiu(2,1,iobs) = prodatqc%mi(jobs,1)
1212         igrdju(2,1,iobs) = prodatqc%mj(jobs,1)-1
1213         igrdiu(2,2,iobs) = prodatqc%mi(jobs,1)
1214         igrdju(2,2,iobs) = prodatqc%mj(jobs,1)
1215         igrdiv(1,1,iobs) = prodatqc%mi(jobs,2)-1
1216         igrdjv(1,1,iobs) = prodatqc%mj(jobs,2)-1
1217         igrdiv(1,2,iobs) = prodatqc%mi(jobs,2)-1
1218         igrdjv(1,2,iobs) = prodatqc%mj(jobs,2)
1219         igrdiv(2,1,iobs) = prodatqc%mi(jobs,2)
1220         igrdjv(2,1,iobs) = prodatqc%mj(jobs,2)-1
1221         igrdiv(2,2,iobs) = prodatqc%mi(jobs,2)
1222         igrdjv(2,2,iobs) = prodatqc%mj(jobs,2)
1223      END DO
1224
1225      CALL obs_int_comm_2d( 2, 2, ipro, igrdiu, igrdju, glamu, zglamu )
1226      CALL obs_int_comm_2d( 2, 2, ipro, igrdiu, igrdju, gphiu, zgphiu )
1227      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, pumask, zumask )
1228      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, pun, zintu )
1229
1230      CALL obs_int_comm_2d( 2, 2, ipro, igrdiv, igrdjv, glamv, zglamv )
1231      CALL obs_int_comm_2d( 2, 2, ipro, igrdiv, igrdjv, gphiv, zgphiv )
1232      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, pvmask, zvmask )
1233      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, pvn, zintv )
1234
1235      ! At the end of the day also get interpolated means
1236      IF ( idayend == 0 ) THEN
1237
1238         ALLOCATE( &
1239            & zinmu(2,2,kpk,ipro),  &
1240            & zinmv(2,2,kpk,ipro)   &
1241            & )
1242
1243         CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, &
1244            &                  prodatqc%vdmean(:,:,:,1), zinmu )
1245         CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, &
1246            &                  prodatqc%vdmean(:,:,:,2), zinmv )
1247
1248      ENDIF
1249
1250! loop over observations
1251
1252      DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
1253
1254         iobs = jobs - prodatqc%nprofup
1255
1256         IF ( kt /= prodatqc%mstp(jobs) ) THEN
1257           
1258            IF(lwp) THEN
1259               WRITE(numout,*)
1260               WRITE(numout,*) ' E R R O R : Observation',              &
1261                  &            ' time step is not consistent with the', &
1262                  &            ' model time step'
1263               WRITE(numout,*) ' ========='
1264               WRITE(numout,*)
1265               WRITE(numout,*) ' Record  = ', jobs,                    &
1266                  &            ' kt      = ', kt,                      &
1267                  &            ' mstp    = ', prodatqc%mstp(jobs), &
1268                  &            ' ntyp    = ', prodatqc%ntyp(jobs)
1269            ENDIF
1270            CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' )
1271         ENDIF
1272         
1273         zlam = prodatqc%rlam(jobs)
1274         zphi = prodatqc%rphi(jobs)
1275
1276         ! Initialize observation masks
1277
1278         zobsmasku(:) = 0.0
1279         zobsmaskv(:) = 0.0
1280         
1281         ! Horizontal weights and vertical mask
1282
1283         IF  ( prodatqc%npvend(jobs,1) > 0 ) THEN
1284
1285            CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi,     &
1286               &                   zglamu(:,:,iobs), zgphiu(:,:,iobs), &
1287               &                   zumask(:,:,:,iobs), zweigu, zobsmasku )
1288
1289         ENDIF
1290
1291         
1292         IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
1293
1294            CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi,     &
1295               &                   zglamv(:,:,iobs), zgphiv(:,:,iobs), &
1296               &                   zvmask(:,:,:,iobs), zweigv, zobsmasku )
1297
1298         ENDIF
1299
1300         ! Ensure that the vertical mask on u and v are consistent.
1301
1302         zobsmask(:) = MIN( zobsmasku(:), zobsmaskv(:) )
1303
1304         IF ( prodatqc%npvend(jobs,1) > 0 ) THEN
1305
1306            zobsk(:) = obfillflt
1307
1308       IF ( ld_dailyav ) THEN
1309
1310               IF ( idayend == 0 )  THEN
1311                 
1312                  ! Daily averaged data
1313                 
1314                  CALL obs_int_h2d( kpk, kpk,      &
1315                     &              zweigu, zinmu(:,:,:,iobs), zobsk )
1316                 
1317                 
1318               ELSE
1319               
1320                  CALL ctl_stop( ' A nonzero' //     &
1321                     &           ' number of U profile data should' // &
1322                     &           ' only occur at the end of a given day' )
1323
1324               ENDIF
1325         
1326            ELSE 
1327               
1328               ! Point data
1329
1330               CALL obs_int_h2d( kpk, kpk,      &
1331                  &              zweigu, zintu(:,:,:,iobs), zobsk )
1332
1333            ENDIF
1334
1335            !-------------------------------------------------------------
1336            ! Compute vertical second-derivative of the interpolating
1337            ! polynomial at obs points
1338            !-------------------------------------------------------------
1339           
1340            IF ( k1dint == 1 ) THEN
1341               CALL obs_int_z1d_spl( kpk, zobsk, zobs2k,   &
1342                  &                  pgdept, zobsmask )
1343            ENDIF
1344           
1345            !-----------------------------------------------------------------
1346            !  Vertical interpolation to the observation point
1347            !-----------------------------------------------------------------
1348            ista = prodatqc%npvsta(jobs,1)
1349            iend = prodatqc%npvend(jobs,1)
1350            CALL obs_int_z1d( kpk,                &
1351               & prodatqc%var(1)%mvk(ista:iend),  &
1352               & k1dint, iend - ista + 1,         &
1353               & prodatqc%var(1)%vdep(ista:iend), &
1354               & zobsk, zobs2k,                   &
1355               & prodatqc%var(1)%vmod(ista:iend), &
1356               & pgdept, zobsmask )
1357
1358         ENDIF
1359
1360         IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
1361
1362            zobsk(:) = obfillflt
1363
1364            IF ( ld_dailyav ) THEN
1365
1366               IF ( idayend == 0 )  THEN
1367
1368                  ! Daily averaged data
1369                 
1370                  CALL obs_int_h2d( kpk, kpk,      &
1371                     &              zweigv, zinmv(:,:,:,iobs), zobsk )
1372                 
1373               ELSE
1374
1375                  CALL ctl_stop( ' A nonzero' //     &
1376                     &           ' number of V profile data should' // &
1377                     &           ' only occur at the end of a given day' )
1378
1379               ENDIF
1380
1381            ELSE
1382               
1383               ! Point data
1384
1385               CALL obs_int_h2d( kpk, kpk,      &
1386                  &              zweigv, zintv(:,:,:,iobs), zobsk )
1387
1388            ENDIF
1389
1390
1391            !-------------------------------------------------------------
1392            ! Compute vertical second-derivative of the interpolating
1393            ! polynomial at obs points
1394            !-------------------------------------------------------------
1395           
1396            IF ( k1dint == 1 ) THEN
1397               CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, &
1398                  &                  pgdept, zobsmask )
1399            ENDIF
1400           
1401            !----------------------------------------------------------------
1402            !  Vertical interpolation to the observation point
1403            !----------------------------------------------------------------
1404            ista = prodatqc%npvsta(jobs,2)
1405            iend = prodatqc%npvend(jobs,2)
1406            CALL obs_int_z1d( kpk, &
1407               & prodatqc%var(2)%mvk(ista:iend),&
1408               & k1dint, iend - ista + 1, &
1409               & prodatqc%var(2)%vdep(ista:iend),&
1410               & zobsk, zobs2k, &
1411               & prodatqc%var(2)%vmod(ista:iend),&
1412               & pgdept, zobsmask )
1413
1414         ENDIF
1415
1416      END DO
1417 
1418      ! Deallocate the data for interpolation
1419      DEALLOCATE( &
1420         & igrdiu, &
1421         & igrdju, &
1422         & igrdiv, &
1423         & igrdjv, &
1424         & zglamu, zglamv, &
1425         & zgphiu, zgphiv, &
1426         & zumask, zvmask, &
1427         & zintu, &
1428         & zintv  &
1429         & )
1430      ! At the end of the day also get interpolated means
1431      IF ( idayend == 0 ) THEN
1432         DEALLOCATE( &
1433            & zinmu,  &
1434            & zinmv   &
1435            & )
1436      ENDIF
1437
1438      prodatqc%nprofup = prodatqc%nprofup + ipro 
1439     
1440   END SUBROUTINE obs_vel_opt
1441
1442END MODULE obs_oper
1443
Note: See TracBrowser for help on using the repository browser.