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

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

Merge of branches/2015/dev_merge_2015 back into trunk. Merge excludes NEMOGCM/TOOLS/OBSTOOLS/ for now due to issues with the change of file type. Will sort these manually with further commits.

Branch merged as follows:
In the working copy of branch ran:
svn merge svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk@HEAD
Small conflicts due to bug fixes applied to trunk since the dev_merge_2015 was copied. Bug fixes were applied to the branch as well so these were easy to resolve.
Branch committed at this stage

In working copy run:
svn switch svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk
to switch working copy

Run:
svn merge —reintegrate svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/2015/dev_merge_2015
to merge the branch into the trunk and then commit - no conflicts at this stage.

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