source: trunk/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90 @ 2715

Last change on this file since 2715 was 2715, checked in by rblod, 10 years ago

First attempt to put dynamic allocation on the trunk

  • Property svn:keywords set to Id
File size: 48.4 KB
Line 
1MODULE obs_oper
2   !!======================================================================
3   !!                       ***  MODULE obs_oper  ***
4   !! Observation diagnostics: Observation operators for various observation
5   !!                          types
6   !!======================================================================
7
8   !!----------------------------------------------------------------------
9   !!   obs_pro_opt :    Compute the model counterpart of temperature and
10   !!                    salinity observations from profiles
11   !!   obs_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, &
617      &                    psstn, psstmask, k2dint )
618
619      !!-----------------------------------------------------------------------
620      !!
621      !!                     ***  ROUTINE obs_sst_opt  ***
622      !!
623      !! ** Purpose : Compute the model counterpart of surface temperature
624      !!              data by interpolating from the model grid to the
625      !!              observation point.
626      !!
627      !! ** Method  : Linearly interpolate to each observation point using
628      !!              the model values at the corners of the surrounding grid box.
629      !!
630      !!    The now model SST is first computed at the obs (lon, lat) point.
631      !!
632      !!    Several horizontal interpolation schemes are available:
633      !!        - distance-weighted (great circle) (k2dint = 0)
634      !!        - distance-weighted (small angle)  (k2dint = 1)
635      !!        - bilinear (geographical grid)     (k2dint = 2)
636      !!        - bilinear (quadrilateral grid)    (k2dint = 3)
637      !!        - polynomial (quadrilateral grid)  (k2dint = 4)
638      !!
639      !!
640      !! ** Action  :
641      !!
642      !! History :
643      !!        !  07-07  (S. Ricci ) : Original
644      !!     
645      !!-----------------------------------------------------------------------
646
647      !! * Modules used
648      USE obs_surf_def  ! Definition of storage space for surface observations
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      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
662         & psstn,  &    ! Model SST field
663         & psstmask     ! Land-sea mask
664         
665      !! * Local declarations
666      INTEGER :: ji
667      INTEGER :: jj
668      INTEGER :: jobs
669      INTEGER :: inrc
670      INTEGER :: isst
671      INTEGER :: iobs
672      REAL(KIND=wp) :: zlam
673      REAL(KIND=wp) :: zphi
674      REAL(KIND=wp) :: zext(1), zobsmask(1)
675      REAL(kind=wp), DIMENSION(2,2,1) :: &
676         & zweig
677      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
678         & zmask, &
679         & zsstl, &
680         & zglam, &
681         & zgphi
682      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
683         & igrdi, &
684         & igrdj
685
686      !-----------------------------------------------------------------------
687      ! Local initialization
688      !-----------------------------------------------------------------------
689      ! ... Record and data counters
690      inrc = kt - kit000 + 2
691      isst = sstdatqc%nsstp(inrc)
692
693      ! Get the data for interpolation
694     
695      ALLOCATE( &
696         & igrdi(2,2,isst), &
697         & igrdj(2,2,isst), &
698         & zglam(2,2,isst), &
699         & zgphi(2,2,isst), &
700         & zmask(2,2,isst), &
701         & zsstl(2,2,isst)  &
702         & )
703     
704      DO jobs = sstdatqc%nsurfup + 1, sstdatqc%nsurfup + isst
705         iobs = jobs - sstdatqc%nsurfup
706         igrdi(1,1,iobs) = sstdatqc%mi(jobs)-1
707         igrdj(1,1,iobs) = sstdatqc%mj(jobs)-1
708         igrdi(1,2,iobs) = sstdatqc%mi(jobs)-1
709         igrdj(1,2,iobs) = sstdatqc%mj(jobs)
710         igrdi(2,1,iobs) = sstdatqc%mi(jobs)
711         igrdj(2,1,iobs) = sstdatqc%mj(jobs)-1
712         igrdi(2,2,iobs) = sstdatqc%mi(jobs)
713         igrdj(2,2,iobs) = sstdatqc%mj(jobs)
714      END DO
715     
716      CALL obs_int_comm_2d( 2, 2, isst, &
717         &                  igrdi, igrdj, glamt, zglam )
718      CALL obs_int_comm_2d( 2, 2, isst, &
719         &                  igrdi, igrdj, gphit, zgphi )
720      CALL obs_int_comm_2d( 2, 2, isst, &
721         &                  igrdi, igrdj, psstmask, zmask )
722      CALL obs_int_comm_2d( 2, 2, isst, &
723         &                  igrdi, igrdj, psstn, zsstl )
724     
725      ! Loop over observations
726
727      DO jobs = sstdatqc%nsurfup + 1, sstdatqc%nsurfup + isst
728         
729         iobs = jobs - sstdatqc%nsurfup
730         
731         IF ( kt /= sstdatqc%mstp(jobs) ) THEN
732           
733            IF(lwp) THEN
734               WRITE(numout,*)
735               WRITE(numout,*) ' E R R O R : Observation',              &
736                  &            ' time step is not consistent with the', &
737                  &            ' model time step'
738               WRITE(numout,*) ' ========='
739               WRITE(numout,*)
740               WRITE(numout,*) ' Record  = ', jobs,                &
741                  &            ' kt      = ', kt,                  &
742                  &            ' mstp    = ', sstdatqc%mstp(jobs), &
743                  &            ' ntyp    = ', sstdatqc%ntyp(jobs)
744            ENDIF
745            CALL ctl_stop( 'obs_sst_opt', 'Inconsistent time' )
746           
747         ENDIF
748         
749         zlam = sstdatqc%rlam(jobs)
750         zphi = sstdatqc%rphi(jobs)
751         
752         ! Get weights to interpolate the model SST to the observation point
753         CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         &
754            &                   zglam(:,:,iobs), zgphi(:,:,iobs), &
755            &                   zmask(:,:,iobs), zweig, zobsmask )
756           
757         ! Interpolate the model SST to the observation point
758         CALL obs_int_h2d( 1, 1,      &
759            &              zweig, zsstl(:,:,iobs),  zext )
760         
761         sstdatqc%rmod(jobs,1) = zext(1)
762         
763      END DO
764     
765      ! Deallocate the data for interpolation
766      DEALLOCATE( &
767         & igrdi, &
768         & igrdj, &
769         & zglam, &
770         & zgphi, &
771         & zmask, &
772         & zsstl  &
773         & )
774     
775      sstdatqc%nsurfup = sstdatqc%nsurfup + isst
776
777   END SUBROUTINE obs_sst_opt
778
779   SUBROUTINE obs_sss_opt
780      !!-----------------------------------------------------------------------
781      !!
782      !!                     ***  ROUTINE obs_sss_opt  ***
783      !!
784      !! ** Purpose : Compute the model counterpart of sea surface salinity
785      !!              data by interpolating from the model grid to the
786      !!              observation point.
787      !!
788      !! ** Method  :
789      !!
790      !! ** Action  :
791      !!
792      !! History :
793      !!      ! ??-??
794      !!-----------------------------------------------------------------------
795
796      IMPLICIT NONE
797
798   END SUBROUTINE obs_sss_opt
799
800   SUBROUTINE obs_seaice_opt( seaicedatqc, kt, kpi, kpj, kit000, &
801      &                    pseaicen, pseaicemask, k2dint )
802
803      !!-----------------------------------------------------------------------
804      !!
805      !!                     ***  ROUTINE obs_seaice_opt  ***
806      !!
807      !! ** Purpose : Compute the model counterpart of surface temperature
808      !!              data by interpolating from the model grid to the
809      !!              observation point.
810      !!
811      !! ** Method  : Linearly interpolate to each observation point using
812      !!              the model values at the corners of the surrounding grid box.
813      !!
814      !!    The now model sea ice is first computed at the obs (lon, lat) point.
815      !!
816      !!    Several horizontal interpolation schemes are available:
817      !!        - distance-weighted (great circle) (k2dint = 0)
818      !!        - distance-weighted (small angle)  (k2dint = 1)
819      !!        - bilinear (geographical grid)     (k2dint = 2)
820      !!        - bilinear (quadrilateral grid)    (k2dint = 3)
821      !!        - polynomial (quadrilateral grid)  (k2dint = 4)
822      !!
823      !!
824      !! ** Action  :
825      !!
826      !! History :
827      !!        !  07-07  (S. Ricci ) : Original
828      !!     
829      !!-----------------------------------------------------------------------
830
831      !! * Modules used
832      USE obs_surf_def  ! Definition of storage space for surface observations
833
834      IMPLICIT NONE
835
836      !! * Arguments
837      TYPE(obs_surf), INTENT(INOUT) :: seaicedatqc     ! Subset of surface data not failing screening
838      INTEGER, INTENT(IN) :: kt       ! Time step
839      INTEGER, INTENT(IN) :: kpi      ! Model grid parameters
840      INTEGER, INTENT(IN) :: kpj
841      INTEGER, INTENT(IN) :: kit000   ! Number of the first time step
842                                      !   (kit000-1 = restart time)
843      INTEGER, INTENT(IN) :: k2dint   ! Horizontal interpolation type (see header)
844      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
845         & pseaicen,  &    ! Model sea ice field
846         & pseaicemask     ! Land-sea mask
847         
848      !! * Local declarations
849      INTEGER :: ji
850      INTEGER :: jj
851      INTEGER :: jobs
852      INTEGER :: inrc
853      INTEGER :: iseaice
854      INTEGER :: iobs
855       
856      REAL(KIND=wp) :: zlam
857      REAL(KIND=wp) :: zphi
858      REAL(KIND=wp) :: zext(1), zobsmask(1)
859      REAL(kind=wp), DIMENSION(2,2,1) :: &
860         & zweig
861      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
862         & zmask, &
863         & zseaicel, &
864         & zglam, &
865         & zgphi
866      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
867         & igrdi, &
868         & igrdj
869
870      !------------------------------------------------------------------------
871      ! Local initialization
872      !------------------------------------------------------------------------
873      ! ... Record and data counters
874      inrc = kt - kit000 + 2
875      iseaice = seaicedatqc%nsstp(inrc)
876
877      ! Get the data for interpolation
878     
879      ALLOCATE( &
880         & igrdi(2,2,iseaice), &
881         & igrdj(2,2,iseaice), &
882         & zglam(2,2,iseaice), &
883         & zgphi(2,2,iseaice), &
884         & zmask(2,2,iseaice), &
885         & zseaicel(2,2,iseaice)  &
886         & )
887     
888      DO jobs = seaicedatqc%nsurfup + 1, seaicedatqc%nsurfup + iseaice
889         iobs = jobs - seaicedatqc%nsurfup
890         igrdi(1,1,iobs) = seaicedatqc%mi(jobs)-1
891         igrdj(1,1,iobs) = seaicedatqc%mj(jobs)-1
892         igrdi(1,2,iobs) = seaicedatqc%mi(jobs)-1
893         igrdj(1,2,iobs) = seaicedatqc%mj(jobs)
894         igrdi(2,1,iobs) = seaicedatqc%mi(jobs)
895         igrdj(2,1,iobs) = seaicedatqc%mj(jobs)-1
896         igrdi(2,2,iobs) = seaicedatqc%mi(jobs)
897         igrdj(2,2,iobs) = seaicedatqc%mj(jobs)
898      END DO
899     
900      CALL obs_int_comm_2d( 2, 2, iseaice, &
901         &                  igrdi, igrdj, glamt, zglam )
902      CALL obs_int_comm_2d( 2, 2, iseaice, &
903         &                  igrdi, igrdj, gphit, zgphi )
904      CALL obs_int_comm_2d( 2, 2, iseaice, &
905         &                  igrdi, igrdj, pseaicemask, zmask )
906      CALL obs_int_comm_2d( 2, 2, iseaice, &
907         &                  igrdi, igrdj, pseaicen, zseaicel )
908     
909      DO jobs = seaicedatqc%nsurfup + 1, seaicedatqc%nsurfup + iseaice
910         
911         iobs = jobs - seaicedatqc%nsurfup
912         
913         IF ( kt /= seaicedatqc%mstp(jobs) ) THEN
914           
915            IF(lwp) THEN
916               WRITE(numout,*)
917               WRITE(numout,*) ' E R R O R : Observation',              &
918                  &            ' time step is not consistent with the', &
919                  &            ' model time step'
920               WRITE(numout,*) ' ========='
921               WRITE(numout,*)
922               WRITE(numout,*) ' Record  = ', jobs,                &
923                  &            ' kt      = ', kt,                  &
924                  &            ' mstp    = ', seaicedatqc%mstp(jobs), &
925                  &            ' ntyp    = ', seaicedatqc%ntyp(jobs)
926            ENDIF
927            CALL ctl_stop( 'obs_seaice_opt', 'Inconsistent time' )
928           
929         ENDIF
930         
931         zlam = seaicedatqc%rlam(jobs)
932         zphi = seaicedatqc%rphi(jobs)
933         
934         ! Get weights to interpolate the model sea ice to the observation point
935         CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         &
936            &                   zglam(:,:,iobs), zgphi(:,:,iobs), &
937            &                   zmask(:,:,iobs), zweig, zobsmask )
938         
939         ! ... Interpolate the model sea ice to the observation point
940         CALL obs_int_h2d( 1, 1,      &
941            &              zweig, zseaicel(:,:,iobs),  zext )
942         
943         seaicedatqc%rmod(jobs,1) = zext(1)
944         
945      END DO
946     
947      ! Deallocate the data for interpolation
948      DEALLOCATE( &
949         & igrdi,    &
950         & igrdj,    &
951         & zglam,    &
952         & zgphi,    &
953         & zmask,    &
954         & zseaicel  &
955         & )
956     
957      seaicedatqc%nsurfup = seaicedatqc%nsurfup + iseaice
958
959   END SUBROUTINE obs_seaice_opt
960
961   SUBROUTINE obs_vel_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, &
962      &                    pun, pvn, pgdept, pumask, pvmask, k1dint, k2dint, &
963      &                    ld_dailyav )
964      !!-----------------------------------------------------------------------
965      !!
966      !!                     ***  ROUTINE obs_vel_opt  ***
967      !!
968      !! ** Purpose : Compute the model counterpart of velocity profile
969      !!              data by interpolating from the model grid to the
970      !!              observation point.
971      !!
972      !! ** Method  : Linearly interpolate zonal and meridional components of velocity
973      !!              to each observation point using the model values at the corners of
974      !!              the surrounding grid box. The model velocity components are on a
975      !!              staggered C- grid.
976      !!
977      !!    For velocity data from the TAO array, the model equivalent is
978      !!    a daily mean velocity field. So, we first compute
979      !!    the mean, then interpolate only at the end of the day.
980      !!
981      !! ** Action  :
982      !!
983      !! History :
984      !!    ! 07-03 (K. Mogensen)      : Temperature and Salinity profiles
985      !!    ! 08-10 (Maria Valdivieso) : Velocity component (U,V) profiles
986      !!-----------------------------------------------------------------------
987   
988      !! * Modules used
989      USE obs_profiles_def ! Definition of storage space for profile obs.
990
991      IMPLICIT NONE
992
993      !! * Arguments
994      TYPE(obs_prof), INTENT(INOUT) :: &
995         & prodatqc        ! Subset of profile data not failing screening
996      INTEGER, INTENT(IN) :: kt        ! Time step
997      INTEGER, INTENT(IN) :: kpi       ! Model grid parameters
998      INTEGER, INTENT(IN) :: kpj
999      INTEGER, INTENT(IN) :: kpk 
1000      INTEGER, INTENT(IN) :: kit000    ! Number of the first time step
1001                                       !   (kit000-1 = restart time)
1002      INTEGER, INTENT(IN) :: k1dint    ! Vertical interpolation type (see header)
1003      INTEGER, INTENT(IN) :: k2dint    ! Horizontal interpolation type (see header)
1004      INTEGER, INTENT(IN) :: kdaystp   ! Number of time steps per day                   
1005      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &
1006         & pun,    &    ! Model zonal component of velocity
1007         & pvn,    &    ! Model meridional component of velocity
1008         & pumask, &    ! Land-sea mask
1009         & pvmask       ! Land-sea mask
1010      REAL(KIND=wp), INTENT(IN), DIMENSION(kpk) :: &
1011         & pgdept       ! Model array of depth levels
1012      LOGICAL, INTENT(IN) :: ld_dailyav
1013         
1014      !! * Local declarations
1015      INTEGER :: ji
1016      INTEGER :: jj
1017      INTEGER :: jk
1018      INTEGER :: jobs
1019      INTEGER :: inrc
1020      INTEGER :: ipro
1021      INTEGER :: idayend
1022      INTEGER :: ista
1023      INTEGER :: iend
1024      INTEGER :: iobs
1025      INTEGER, DIMENSION(imaxavtypes) :: &
1026         & idailyavtypes
1027      REAL(KIND=wp) :: zlam
1028      REAL(KIND=wp) :: zphi
1029      REAL(KIND=wp) :: zdaystp
1030      REAL(KIND=wp), DIMENSION(kpk) :: &
1031         & zobsmasku, &
1032         & zobsmaskv, &
1033         & zobsmask,  &
1034         & zobsk,     &
1035         & zobs2k
1036      REAL(KIND=wp), DIMENSION(2,2,kpk) :: &
1037         & zweigu,zweigv
1038      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
1039         & zumask, zvmask, &
1040         & zintu, &
1041         & zintv, &
1042         & zinmu, &
1043         & zinmv
1044      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
1045         & zglamu, zglamv, &
1046         & zgphiu, zgphiv
1047      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
1048         & igrdiu, &
1049         & igrdju, &
1050         & igrdiv, &
1051         & igrdjv
1052
1053      !------------------------------------------------------------------------
1054      ! Local initialization
1055      !------------------------------------------------------------------------
1056      ! ... Record and data counters
1057      inrc = kt - kit000 + 2
1058      ipro = prodatqc%npstp(inrc)
1059
1060      ! Initialize daily mean for first timestep
1061      idayend = MOD( kt - kit000 + 1, kdaystp )
1062
1063      ! Added kt == 0 test to catch restart case
1064      IF ( idayend == 1 .OR. kt == 0) THEN
1065         IF (lwp) WRITE(numout,*) 'Reset prodatqc%vdmean on time-step: ',kt
1066         prodatqc%vdmean(:,:,:,1) = 0.0
1067         prodatqc%vdmean(:,:,:,2) = 0.0
1068      ENDIF
1069
1070      ! Increment the zonal velocity field for computing daily mean
1071      prodatqc%vdmean(:,:,:,1) = prodatqc%vdmean(:,:,:,1) + pun(:,:,:)
1072      ! Increment the meridional velocity field for computing daily mean
1073      prodatqc%vdmean(:,:,:,2) = prodatqc%vdmean(:,:,:,2) + pvn(:,:,:)
1074   
1075      ! Compute the daily mean at the end of day
1076      zdaystp = 1.0 / REAL( kdaystp )
1077      IF ( idayend == 0 ) THEN
1078         prodatqc%vdmean(:,:,:,1) = prodatqc%vdmean(:,:,:,1) * zdaystp
1079         prodatqc%vdmean(:,:,:,2) = prodatqc%vdmean(:,:,:,2) * zdaystp
1080      ENDIF
1081
1082      ! Get the data for interpolation
1083      ALLOCATE( &
1084         & igrdiu(2,2,ipro),      &
1085         & igrdju(2,2,ipro),      &
1086         & igrdiv(2,2,ipro),      &
1087         & igrdjv(2,2,ipro),      &
1088         & zglamu(2,2,ipro), zglamv(2,2,ipro), &
1089         & zgphiu(2,2,ipro), zgphiv(2,2,ipro), &
1090         & zumask(2,2,kpk,ipro), zvmask(2,2,kpk,ipro), &
1091         & zintu(2,2,kpk,ipro),  &
1092         & zintv(2,2,kpk,ipro)   &
1093         & )
1094
1095      DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
1096         iobs = jobs - prodatqc%nprofup
1097         igrdiu(1,1,iobs) = prodatqc%mi(jobs,1)-1
1098         igrdju(1,1,iobs) = prodatqc%mj(jobs,1)-1
1099         igrdiu(1,2,iobs) = prodatqc%mi(jobs,1)-1
1100         igrdju(1,2,iobs) = prodatqc%mj(jobs,1)
1101         igrdiu(2,1,iobs) = prodatqc%mi(jobs,1)
1102         igrdju(2,1,iobs) = prodatqc%mj(jobs,1)-1
1103         igrdiu(2,2,iobs) = prodatqc%mi(jobs,1)
1104         igrdju(2,2,iobs) = prodatqc%mj(jobs,1)
1105         igrdiv(1,1,iobs) = prodatqc%mi(jobs,2)-1
1106         igrdjv(1,1,iobs) = prodatqc%mj(jobs,2)-1
1107         igrdiv(1,2,iobs) = prodatqc%mi(jobs,2)-1
1108         igrdjv(1,2,iobs) = prodatqc%mj(jobs,2)
1109         igrdiv(2,1,iobs) = prodatqc%mi(jobs,2)
1110         igrdjv(2,1,iobs) = prodatqc%mj(jobs,2)-1
1111         igrdiv(2,2,iobs) = prodatqc%mi(jobs,2)
1112         igrdjv(2,2,iobs) = prodatqc%mj(jobs,2)
1113      END DO
1114
1115      CALL obs_int_comm_2d( 2, 2, ipro, igrdiu, igrdju, glamu, zglamu )
1116      CALL obs_int_comm_2d( 2, 2, ipro, igrdiu, igrdju, gphiu, zgphiu )
1117      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, pumask, zumask )
1118      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, pun, zintu )
1119
1120      CALL obs_int_comm_2d( 2, 2, ipro, igrdiv, igrdjv, glamv, zglamv )
1121      CALL obs_int_comm_2d( 2, 2, ipro, igrdiv, igrdjv, gphiv, zgphiv )
1122      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, pvmask, zvmask )
1123      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, pvn, zintv )
1124
1125      ! At the end of the day also get interpolated means
1126      IF ( idayend == 0 ) THEN
1127
1128         ALLOCATE( &
1129            & zinmu(2,2,kpk,ipro),  &
1130            & zinmv(2,2,kpk,ipro)   &
1131            & )
1132
1133         CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, &
1134            &                  prodatqc%vdmean(:,:,:,1), zinmu )
1135         CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, &
1136            &                  prodatqc%vdmean(:,:,:,2), zinmv )
1137
1138      ENDIF
1139
1140! loop over observations
1141
1142      DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
1143
1144         iobs = jobs - prodatqc%nprofup
1145
1146         IF ( kt /= prodatqc%mstp(jobs) ) THEN
1147           
1148            IF(lwp) THEN
1149               WRITE(numout,*)
1150               WRITE(numout,*) ' E R R O R : Observation',              &
1151                  &            ' time step is not consistent with the', &
1152                  &            ' model time step'
1153               WRITE(numout,*) ' ========='
1154               WRITE(numout,*)
1155               WRITE(numout,*) ' Record  = ', jobs,                    &
1156                  &            ' kt      = ', kt,                      &
1157                  &            ' mstp    = ', prodatqc%mstp(jobs), &
1158                  &            ' ntyp    = ', prodatqc%ntyp(jobs)
1159            ENDIF
1160            CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' )
1161         ENDIF
1162         
1163         zlam = prodatqc%rlam(jobs)
1164         zphi = prodatqc%rphi(jobs)
1165
1166         ! Initialize observation masks
1167
1168         zobsmasku(:) = 0.0
1169         zobsmaskv(:) = 0.0
1170         
1171         ! Horizontal weights and vertical mask
1172
1173         IF  ( prodatqc%npvend(jobs,1) > 0 ) THEN
1174
1175            CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi,     &
1176               &                   zglamu(:,:,iobs), zgphiu(:,:,iobs), &
1177               &                   zumask(:,:,:,iobs), zweigu, zobsmasku )
1178
1179         ENDIF
1180
1181         
1182         IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
1183
1184            CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi,     &
1185               &                   zglamv(:,:,iobs), zgphiv(:,:,iobs), &
1186               &                   zvmask(:,:,:,iobs), zweigv, zobsmasku )
1187
1188         ENDIF
1189
1190         ! Ensure that the vertical mask on u and v are consistent.
1191
1192         zobsmask(:) = MIN( zobsmasku(:), zobsmaskv(:) )
1193
1194         IF ( prodatqc%npvend(jobs,1) > 0 ) THEN
1195
1196            zobsk(:) = obfillflt
1197
1198       IF ( ld_dailyav ) THEN
1199
1200               IF ( idayend == 0 )  THEN
1201                 
1202                  ! Daily averaged data
1203                 
1204                  CALL obs_int_h2d( kpk, kpk,      &
1205                     &              zweigu, zinmu(:,:,:,iobs), zobsk )
1206                 
1207                 
1208               ELSE
1209               
1210                  CALL ctl_stop( ' A nonzero' //     &
1211                     &           ' number of U profile data should' // &
1212                     &           ' only occur at the end of a given day' )
1213
1214               ENDIF
1215         
1216            ELSE 
1217               
1218               ! Point data
1219
1220               CALL obs_int_h2d( kpk, kpk,      &
1221                  &              zweigu, zintu(:,:,:,iobs), zobsk )
1222
1223            ENDIF
1224
1225            !-------------------------------------------------------------
1226            ! Compute vertical second-derivative of the interpolating
1227            ! polynomial at obs points
1228            !-------------------------------------------------------------
1229           
1230            IF ( k1dint == 1 ) THEN
1231               CALL obs_int_z1d_spl( kpk, zobsk, zobs2k,   &
1232                  &                  pgdept, zobsmask )
1233            ENDIF
1234           
1235            !-----------------------------------------------------------------
1236            !  Vertical interpolation to the observation point
1237            !-----------------------------------------------------------------
1238            ista = prodatqc%npvsta(jobs,1)
1239            iend = prodatqc%npvend(jobs,1)
1240            CALL obs_int_z1d( kpk,                &
1241               & prodatqc%var(1)%mvk(ista:iend),  &
1242               & k1dint, iend - ista + 1,         &
1243               & prodatqc%var(1)%vdep(ista:iend), &
1244               & zobsk, zobs2k,                   &
1245               & prodatqc%var(1)%vmod(ista:iend), &
1246               & pgdept, zobsmask )
1247
1248         ENDIF
1249
1250         IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
1251
1252            zobsk(:) = obfillflt
1253
1254            IF ( ld_dailyav ) THEN
1255
1256               IF ( idayend == 0 )  THEN
1257
1258                  ! Daily averaged data
1259                 
1260                  CALL obs_int_h2d( kpk, kpk,      &
1261                     &              zweigv, zinmv(:,:,:,iobs), zobsk )
1262                 
1263               ELSE
1264
1265                  CALL ctl_stop( ' A nonzero' //     &
1266                     &           ' number of V profile data should' // &
1267                     &           ' only occur at the end of a given day' )
1268
1269               ENDIF
1270
1271            ELSE
1272               
1273               ! Point data
1274
1275               CALL obs_int_h2d( kpk, kpk,      &
1276                  &              zweigv, zintv(:,:,:,iobs), zobsk )
1277
1278            ENDIF
1279
1280
1281            !-------------------------------------------------------------
1282            ! Compute vertical second-derivative of the interpolating
1283            ! polynomial at obs points
1284            !-------------------------------------------------------------
1285           
1286            IF ( k1dint == 1 ) THEN
1287               CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, &
1288                  &                  pgdept, zobsmask )
1289            ENDIF
1290           
1291            !----------------------------------------------------------------
1292            !  Vertical interpolation to the observation point
1293            !----------------------------------------------------------------
1294            ista = prodatqc%npvsta(jobs,2)
1295            iend = prodatqc%npvend(jobs,2)
1296            CALL obs_int_z1d( kpk, &
1297               & prodatqc%var(2)%mvk(ista:iend),&
1298               & k1dint, iend - ista + 1, &
1299               & prodatqc%var(2)%vdep(ista:iend),&
1300               & zobsk, zobs2k, &
1301               & prodatqc%var(2)%vmod(ista:iend),&
1302               & pgdept, zobsmask )
1303
1304         ENDIF
1305
1306      END DO
1307 
1308      ! Deallocate the data for interpolation
1309      DEALLOCATE( &
1310         & igrdiu, &
1311         & igrdju, &
1312         & igrdiv, &
1313         & igrdjv, &
1314         & zglamu, zglamv, &
1315         & zgphiu, zgphiv, &
1316         & zumask, zvmask, &
1317         & zintu, &
1318         & zintv  &
1319         & )
1320      ! At the end of the day also get interpolated means
1321      IF ( idayend == 0 ) THEN
1322         DEALLOCATE( &
1323            & zinmu,  &
1324            & zinmv   &
1325            & )
1326      ENDIF
1327
1328      prodatqc%nprofup = prodatqc%nprofup + ipro 
1329     
1330   END SUBROUTINE obs_vel_opt
1331
1332END MODULE obs_oper
1333
Note: See TracBrowser for help on using the repository browser.