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

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

source: branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90 @ 5682

Last change on this file since 5682 was 5682, checked in by mattmartin, 9 years ago

OBS simplification changes committed to branch after running SETTE tests to make sure we get the same results as the trunk for ORCA2_LIM_OBS.
  • Property svn:keywords set to Id
File size: 27.1 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   !!----------------------------------------------------------------------
12
13   !! * Modules used
14   USE par_kind, ONLY : &         ! Precision variables
15      & wp
16   USE in_out_manager             ! I/O manager
17   USE obs_inter_sup              ! Interpolation support
18   USE obs_inter_h2d, ONLY : &    ! Horizontal interpolation to the obs pt
19      & obs_int_h2d, &
20      & obs_int_h2d_init
21   USE obs_inter_z1d, ONLY : &    ! Vertical interpolation to the obs pt
22      & obs_int_z1d,    &
23      & obs_int_z1d_spl
24   USE obs_const,  ONLY :    &    ! Obs fill value
25      & obfillflt
26   USE dom_oce,       ONLY : &    ! Model lats/lons
27      & glamt, &
28      & gphit
29   USE lib_mpp,       ONLY : &    ! Warning and stopping routines
30      & ctl_warn, ctl_stop
31   USE sbcdcy,        ONLY : &    ! For calculation of where it is night-time
32      & sbc_dcy, nday_qsr
33
34   IMPLICIT NONE
35
36   !! * Routine accessibility
37   PRIVATE
38
39   PUBLIC obs_prof_opt, &  ! Compute the model counterpart of profile obs
40      &   obs_surf_opt     ! Compute the model counterpart of surface obs
41
42   INTEGER, PARAMETER, PUBLIC :: &
43      & imaxavtypes = 20   ! Max number of daily avgd obs types
44
45   !!----------------------------------------------------------------------
46   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
47   !! $Id$
48   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
49   !!----------------------------------------------------------------------
50
51CONTAINS
52
53   SUBROUTINE obs_prof_opt( prodatqc, kt, kpi, kpj, kpk,          &
54      &                     kit000, kdaystp,                      &
55      &                     pvar1, pvar2, pgdept, pmask1, pmask2, &
56      &                     plam1, plam2, pphi1, pphi2,           &
57      &                     k1dint, k2dint, kdailyavtypes )
58
59      !!-----------------------------------------------------------------------
60      !!
61      !!                     ***  ROUTINE obs_pro_opt  ***
62      !!
63      !! ** Purpose : Compute the model counterpart of profiles
64      !!              data by interpolating from the model grid to the
65      !!              observation point.
66      !!
67      !! ** Method  : Linearly interpolate to each observation point using
68      !!              the model values at the corners of the surrounding grid box.
69      !!
70      !!    First, a vertical profile of horizontally interpolated model
71      !!    now values is computed at the obs (lon, lat) point.
72      !!    Several horizontal interpolation schemes are available:
73      !!        - distance-weighted (great circle) (k2dint = 0)
74      !!        - distance-weighted (small angle)  (k2dint = 1)
75      !!        - bilinear (geographical grid)     (k2dint = 2)
76      !!        - bilinear (quadrilateral grid)    (k2dint = 3)
77      !!        - polynomial (quadrilateral grid)  (k2dint = 4)
78      !!
79      !!    Next, the vertical profile is interpolated to the
80      !!    data depth points. Two vertical interpolation schemes are
81      !!    available:
82      !!        - linear       (k1dint = 0)
83      !!        - Cubic spline (k1dint = 1)
84      !!
85      !!    For the cubic spline the 2nd derivative of the interpolating
86      !!    polynomial is computed before entering the vertical interpolation
87      !!    routine.
88      !!
89      !!    If the logical is switched on, the model equivalent is
90      !!    a daily mean model temperature field. So, we first compute
91      !!    the mean, then interpolate only at the end of the day.
92      !!
93      !!    Note: in situ temperature observations must be converted
94      !!    to potential temperature (the model variable) prior to
95      !!    assimilation.
96      !!
97      !! ** Action  :
98      !!
99      !! History :
100      !!      ! 97-11 (A. Weaver, S. Ricci, N. Daget)
101      !!      ! 06-03 (G. Smith) NEMOVAR migration
102      !!      ! 06-10 (A. Weaver) Cleanup
103      !!      ! 07-01 (K. Mogensen) Merge of temperature and salinity
104      !!      ! 07-03 (K. Mogensen) General handling of profiles
105      !!      ! 15-02 (M. Martin) Combined routine for all profile types
106      !!-----------------------------------------------------------------------
107
108      !! * Modules used
109      USE obs_profiles_def ! Definition of storage space for profile obs.
110
111      IMPLICIT NONE
112
113      !! * Arguments
114      TYPE(obs_prof), INTENT(INOUT) :: &
115         & prodatqc                  ! Subset of profile data passing QC
116      INTEGER, INTENT(IN) :: kt      ! Time step
117      INTEGER, INTENT(IN) :: kpi     ! Model grid parameters
118      INTEGER, INTENT(IN) :: kpj
119      INTEGER, INTENT(IN) :: kpk
120      INTEGER, INTENT(IN) :: kit000  ! Number of the first time step
121                                     !   (kit000-1 = restart time)
122      INTEGER, INTENT(IN) :: k1dint  ! Vertical interpolation type (see header)
123      INTEGER, INTENT(IN) :: k2dint  ! Horizontal interpolation type (see header)
124      INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day
125      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &
126         & pvar1,    &               ! Model field 1
127         & pvar2,    &               ! Model field 2
128         & pmask1,   &               ! Land-sea mask 1
129         & pmask2                    ! Land-sea mask 2
130      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
131         & plam1,    &               ! Model longitudes for variable 1
132         & plam2,    &               ! Model longitudes for variable 2
133         & pphi1,    &               ! Model latitudes for variable 1
134         & pphi2                     ! Model latitudes for variable 2
135      REAL(KIND=wp), INTENT(IN), DIMENSION(kpk) :: &
136         & pgdept                    ! Model array of depth levels
137      INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &
138         & kdailyavtypes             ! Types for daily averages
139
140      !! * Local declarations
141      INTEGER ::   ji
142      INTEGER ::   jj
143      INTEGER ::   jk
144      INTEGER ::   jobs
145      INTEGER ::   inrc
146      INTEGER ::   ipro
147      INTEGER ::   idayend
148      INTEGER ::   ista
149      INTEGER ::   iend
150      INTEGER ::   iobs
151      INTEGER, DIMENSION(imaxavtypes) :: &
152         & idailyavtypes
153      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
154         & igrdi1, &
155         & igrdi2, &
156         & igrdj1, &
157         & igrdj2
158      REAL(KIND=wp) :: zlam
159      REAL(KIND=wp) :: zphi
160      REAL(KIND=wp) :: zdaystp
161      REAL(KIND=wp), DIMENSION(kpk) :: &
162         & zobsmask1, &
163         & zobsmask2, &
164         & zobsk,    &
165         & zobs2k
166      REAL(KIND=wp), DIMENSION(2,2,kpk) :: &
167         & zweig1, &
168         & zweig2
169      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
170         & zmask1, &
171         & zmask2, &
172         & zint1, &
173         & zint2, &
174         & zinm1, &
175         & zinm2
176      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
177         & zglam1, &
178         & zglam2, &
179         & zgphi1, &
180         & zgphi2
181      LOGICAL :: ld_dailyav
182
183      !------------------------------------------------------------------------
184      ! Local initialization
185      !------------------------------------------------------------------------
186      ! Record and data counters
187      inrc = kt - kit000 + 2
188      ipro = prodatqc%npstp(inrc)
189
190      ! Daily average types
191      ld_dailyav = .FALSE.
192      IF ( PRESENT(kdailyavtypes) ) THEN
193         idailyavtypes(:) = kdailyavtypes(:)
194         IF ( ANY (idailyavtypes(:) /= -1) ) ld_dailyav = .TRUE.
195      ELSE
196         idailyavtypes(:) = -1
197      ENDIF
198
199      ! Daily means are calculated for values over timesteps:
200      !  [1 <= kt <= kdaystp], [kdaystp+1 <= kt <= 2*kdaystp], ...
201      idayend = MOD( kt - kit000 + 1, kdaystp )
202
203      IF ( ld_dailyav ) THEN
204
205         ! Initialize daily mean for first timestep of the day
206         IF ( idayend == 1 .OR. kt == 0 ) THEN
207            DO jk = 1, jpk
208               DO jj = 1, jpj
209                  DO ji = 1, jpi
210                     prodatqc%vdmean(ji,jj,jk,1) = 0.0
211                     prodatqc%vdmean(ji,jj,jk,2) = 0.0
212                  END DO
213               END DO
214            END DO
215         ENDIF
216
217         DO jk = 1, jpk
218            DO jj = 1, jpj
219               DO ji = 1, jpi
220                  ! Increment field 1 for computing daily mean
221                  prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &
222                     &                        + pvar1(ji,jj,jk)
223                  ! Increment field 2 for computing daily mean
224                  prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &
225                     &                        + pvar2(ji,jj,jk)
226               END DO
227            END DO
228         END DO
229
230         ! Compute the daily mean at the end of day
231         zdaystp = 1.0 / REAL( kdaystp )
232         IF ( idayend == 0 ) THEN
233            IF (lwp) WRITE(numout,*) 'Calculating prodatqc%vdmean on time-step: ',kt
234            CALL FLUSH(numout)
235            DO jk = 1, jpk
236               DO jj = 1, jpj
237                  DO ji = 1, jpi
238                     prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &
239                        &                        * zdaystp
240                     prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &
241                        &                        * zdaystp
242                  END DO
243               END DO
244            END DO
245         ENDIF
246
247      ENDIF
248
249      ! Get the data for interpolation
250      ALLOCATE( &
251         & igrdi1(2,2,ipro),      &
252         & igrdi2(2,2,ipro),      &
253         & igrdj1(2,2,ipro),      &
254         & igrdj2(2,2,ipro),      &
255         & zglam1(2,2,ipro),      &
256         & zglam2(2,2,ipro),      &
257         & zgphi1(2,2,ipro),      &
258         & zgphi2(2,2,ipro),      &
259         & zmask1(2,2,kpk,ipro),  &
260         & zmask2(2,2,kpk,ipro),  &
261         & zint1(2,2,kpk,ipro),  &
262         & zint2(2,2,kpk,ipro)   &
263         & )
264
265      DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
266         iobs = jobs - prodatqc%nprofup
267         igrdi1(1,1,iobs) = prodatqc%mi(jobs,1)-1
268         igrdj1(1,1,iobs) = prodatqc%mj(jobs,1)-1
269         igrdi1(1,2,iobs) = prodatqc%mi(jobs,1)-1
270         igrdj1(1,2,iobs) = prodatqc%mj(jobs,1)
271         igrdi1(2,1,iobs) = prodatqc%mi(jobs,1)
272         igrdj1(2,1,iobs) = prodatqc%mj(jobs,1)-1
273         igrdi1(2,2,iobs) = prodatqc%mi(jobs,1)
274         igrdj1(2,2,iobs) = prodatqc%mj(jobs,1)
275         igrdi2(1,1,iobs) = prodatqc%mi(jobs,2)-1
276         igrdj2(1,1,iobs) = prodatqc%mj(jobs,2)-1
277         igrdi2(1,2,iobs) = prodatqc%mi(jobs,2)-1
278         igrdj2(1,2,iobs) = prodatqc%mj(jobs,2)
279         igrdi2(2,1,iobs) = prodatqc%mi(jobs,2)
280         igrdj2(2,1,iobs) = prodatqc%mj(jobs,2)-1
281         igrdi2(2,2,iobs) = prodatqc%mi(jobs,2)
282         igrdj2(2,2,iobs) = prodatqc%mj(jobs,2)
283      END DO
284
285      CALL obs_int_comm_2d( 2, 2, ipro, igrdi1, igrdj1, plam1, zglam1 )
286      CALL obs_int_comm_2d( 2, 2, ipro, igrdi1, igrdj1, pphi1, zgphi1 )
287      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi1, igrdj1, pmask1, zmask1 )
288      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi1, igrdj1, pvar1,   zint1 )
289     
290      CALL obs_int_comm_2d( 2, 2, ipro, igrdi2, igrdj2, plam2, zglam2 )
291      CALL obs_int_comm_2d( 2, 2, ipro, igrdi2, igrdj2, pphi2, zgphi2 )
292      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi2, igrdj2, pmask2, zmask2 )
293      CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi2, igrdj2, pvar2,   zint2 )
294
295      ! At the end of the day also get interpolated means
296      IF ( ld_dailyav .AND. idayend == 0 ) THEN
297
298         ALLOCATE( &
299            & zinm1(2,2,kpk,ipro),  &
300            & zinm2(2,2,kpk,ipro)   &
301            & )
302
303         CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi1, igrdj1, &
304            &                  prodatqc%vdmean(:,:,:,1), zinm1 )
305         CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi2, igrdj2, &
306            &                  prodatqc%vdmean(:,:,:,2), zinm2 )
307
308      ENDIF
309
310      DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
311
312         iobs = jobs - prodatqc%nprofup
313
314         IF ( kt /= prodatqc%mstp(jobs) ) THEN
315
316            IF(lwp) THEN
317               WRITE(numout,*)
318               WRITE(numout,*) ' E R R O R : Observation',              &
319                  &            ' time step is not consistent with the', &
320                  &            ' model time step'
321               WRITE(numout,*) ' ========='
322               WRITE(numout,*)
323               WRITE(numout,*) ' Record  = ', jobs,                    &
324                  &            ' kt      = ', kt,                      &
325                  &            ' mstp    = ', prodatqc%mstp(jobs), &
326                  &            ' ntyp    = ', prodatqc%ntyp(jobs)
327            ENDIF
328            CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' )
329         ENDIF
330
331         zlam = prodatqc%rlam(jobs)
332         zphi = prodatqc%rphi(jobs)
333
334         ! Horizontal weights and vertical mask
335
336         IF ( prodatqc%npvend(jobs,1) > 0 ) THEN
337
338            CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi,     &
339               &                   zglam1(:,:,iobs), zgphi1(:,:,iobs), &
340               &                   zmask1(:,:,:,iobs), zweig1, zobsmask1 )
341
342         ENDIF
343
344         IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
345
346            CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi,     &
347               &                   zglam2(:,:,iobs), zgphi2(:,:,iobs), &
348               &                   zmask2(:,:,:,iobs), zweig2, zobsmask2 )
349 
350         ENDIF
351
352         IF ( prodatqc%npvend(jobs,1) > 0 ) THEN
353
354            zobsk(:) = obfillflt
355
356            IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN
357
358               IF ( idayend == 0 )  THEN
359                  ! Daily averaged data
360                  CALL obs_int_h2d( kpk, kpk,      &
361                     &              zweig1, zinm1(:,:,:,iobs), zobsk )
362
363               ENDIF
364
365            ELSE 
366
367               ! Point data
368               CALL obs_int_h2d( kpk, kpk,      &
369                  &              zweig1, zint1(:,:,:,iobs), zobsk )
370
371            ENDIF
372
373            !-------------------------------------------------------------
374            ! Compute vertical second-derivative of the interpolating
375            ! polynomial at obs points
376            !-------------------------------------------------------------
377
378            IF ( k1dint == 1 ) THEN
379               CALL obs_int_z1d_spl( kpk, zobsk, zobs2k,   &
380                  &                  pgdept, zobsmask1 )
381            ENDIF
382
383            !-----------------------------------------------------------------
384            !  Vertical interpolation to the observation point
385            !-----------------------------------------------------------------
386            ista = prodatqc%npvsta(jobs,1)
387            iend = prodatqc%npvend(jobs,1)
388            CALL obs_int_z1d( kpk,                &
389               & prodatqc%var(1)%mvk(ista:iend),  &
390               & k1dint, iend - ista + 1,         &
391               & prodatqc%var(1)%vdep(ista:iend), &
392               & zobsk, zobs2k,                   &
393               & prodatqc%var(1)%vmod(ista:iend), &
394               & pgdept, zobsmask1 )
395
396         ENDIF
397
398         IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
399
400            zobsk(:) = obfillflt
401
402            IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN
403
404               IF ( idayend == 0 )  THEN
405
406                  ! Daily averaged data
407                  CALL obs_int_h2d( kpk, kpk,      &
408                     &              zweig2, zinm2(:,:,:,iobs), zobsk )
409
410               ENDIF
411
412            ELSE
413
414               ! Point data
415               CALL obs_int_h2d( kpk, kpk,      &
416                  &              zweig2, zint2(:,:,:,iobs), zobsk )
417
418            ENDIF
419
420
421            !-------------------------------------------------------------
422            ! Compute vertical second-derivative of the interpolating
423            ! polynomial at obs points
424            !-------------------------------------------------------------
425
426            IF ( k1dint == 1 ) THEN
427               CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, &
428                  &                  pgdept, zobsmask2 )
429            ENDIF
430
431            !----------------------------------------------------------------
432            !  Vertical interpolation to the observation point
433            !----------------------------------------------------------------
434            ista = prodatqc%npvsta(jobs,2)
435            iend = prodatqc%npvend(jobs,2)
436            CALL obs_int_z1d( kpk, &
437               & prodatqc%var(2)%mvk(ista:iend),&
438               & k1dint, iend - ista + 1, &
439               & prodatqc%var(2)%vdep(ista:iend),&
440               & zobsk, zobs2k, &
441               & prodatqc%var(2)%vmod(ista:iend),&
442               & pgdept, zobsmask2 )
443
444         ENDIF
445
446      END DO
447
448      ! Deallocate the data for interpolation
449      DEALLOCATE( &
450         & igrdi1, &
451         & igrdi2, &
452         & igrdj1, &
453         & igrdj2, &
454         & zglam1, &
455         & zglam2, &
456         & zgphi1, &
457         & zgphi2, &
458         & zmask1, &
459         & zmask2, &
460         & zint1,  &
461         & zint2   &
462         & )
463
464      ! At the end of the day also get interpolated means
465      IF ( ld_dailyav .AND. idayend == 0 ) THEN
466         DEALLOCATE( &
467            & zinm1,  &
468            & zinm2   &
469            & )
470      ENDIF
471
472      prodatqc%nprofup = prodatqc%nprofup + ipro 
473
474   END SUBROUTINE obs_prof_opt
475
476   SUBROUTINE obs_surf_opt( surfdataqc, kt, kpi, kpj,         &
477      &                    kit000, kdaystp, psurf, psurfmask, &
478      &                    k2dint, ldnightav )
479
480      !!-----------------------------------------------------------------------
481      !!
482      !!                     ***  ROUTINE obs_surf_opt  ***
483      !!
484      !! ** Purpose : Compute the model counterpart of surface
485      !!              data by interpolating from the model grid to the
486      !!              observation point.
487      !!
488      !! ** Method  : Linearly interpolate to each observation point using
489      !!              the model values at the corners of the surrounding grid box.
490      !!
491      !!    The new model value is first computed at the obs (lon, lat) point.
492      !!
493      !!    Several horizontal interpolation schemes are available:
494      !!        - distance-weighted (great circle) (k2dint = 0)
495      !!        - distance-weighted (small angle)  (k2dint = 1)
496      !!        - bilinear (geographical grid)     (k2dint = 2)
497      !!        - bilinear (quadrilateral grid)    (k2dint = 3)
498      !!        - polynomial (quadrilateral grid)  (k2dint = 4)
499      !!
500      !!
501      !! ** Action  :
502      !!
503      !! History :
504      !!      ! 07-03 (A. Weaver)
505      !!      ! 15-02 (M. Martin) Combined routine for surface types
506      !!-----------------------------------------------------------------------
507
508      !! * Modules used
509      USE obs_surf_def  ! Definition of storage space for surface observations
510
511      IMPLICIT NONE
512
513      !! * Arguments
514      TYPE(obs_surf), INTENT(INOUT) :: &
515         & surfdataqc                  ! Subset of surface data passing QC
516      INTEGER, INTENT(IN) :: kt        ! Time step
517      INTEGER, INTENT(IN) :: kpi       ! Model grid parameters
518      INTEGER, INTENT(IN) :: kpj
519      INTEGER, INTENT(IN) :: kit000    ! Number of the first time step
520                                       !   (kit000-1 = restart time)
521      INTEGER, INTENT(IN) :: kdaystp   ! Number of time steps per day
522      INTEGER, INTENT(IN) :: k2dint    ! Horizontal interpolation type (see header)
523      REAL(wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
524         & psurf,  &                   ! Model surface field
525         & psurfmask                   ! Land-sea mask
526      LOGICAL, INTENT(IN) :: ldnightav ! Logical for averaging night-time data
527
528      !! * Local declarations
529      INTEGER :: ji
530      INTEGER :: jj
531      INTEGER :: jobs
532      INTEGER :: inrc
533      INTEGER :: isurf
534      INTEGER :: iobs
535      INTEGER :: idayend
536      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
537         & igrdi, &
538         & igrdj
539      INTEGER, DIMENSION(:,:), SAVE, ALLOCATABLE :: &
540         & icount_night,      &
541         & imask_night
542      REAL(wp) :: zlam
543      REAL(wp) :: zphi
544      REAL(wp), DIMENSION(1) :: zext, zobsmask
545      REAL(wp) :: zdaystp
546      REAL(wp), DIMENSION(2,2,1) :: &
547         & zweig
548      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
549         & zmask,  &
550         & zsurf,  &
551         & zsurfm, &
552         & zglam,  &
553         & zgphi
554      REAL(wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: &
555         & zintmp,  &
556         & zouttmp, &
557         & zmeanday    ! to compute model sst in region of 24h daylight (pole)
558
559      !------------------------------------------------------------------------
560      ! Local initialization
561      !------------------------------------------------------------------------
562      ! Record and data counters
563      inrc = kt - kit000 + 2
564      isurf = surfdataqc%nsstp(inrc)
565
566      IF ( ldnightav ) THEN
567
568      ! Initialize array for night mean
569         IF ( kt == 0 ) THEN
570            ALLOCATE ( icount_night(kpi,kpj) )
571            ALLOCATE ( imask_night(kpi,kpj) )
572            ALLOCATE ( zintmp(kpi,kpj) )
573            ALLOCATE ( zouttmp(kpi,kpj) )
574            ALLOCATE ( zmeanday(kpi,kpj) )
575            nday_qsr = -1   ! initialisation flag for nbc_dcy
576         ENDIF
577
578         ! Night-time means are calculated for night-time values over timesteps:
579         !  [1 <= kt <= kdaystp], [kdaystp+1 <= kt <= 2*kdaystp], .....
580         idayend = MOD( kt - kit000 + 1, kdaystp )
581
582         ! Initialize night-time mean for first timestep of the day
583         IF ( idayend == 1 .OR. kt == 0 ) THEN
584            DO jj = 1, jpj
585               DO ji = 1, jpi
586                  surfdataqc%vdmean(ji,jj) = 0.0
587                  zmeanday(ji,jj) = 0.0
588                  icount_night(ji,jj) = 0
589               END DO
590            END DO
591         ENDIF
592
593         zintmp(:,:) = 0.0
594         zouttmp(:,:) = sbc_dcy( zintmp(:,:), .TRUE. )
595         imask_night(:,:) = INT( zouttmp(:,:) )
596
597         DO jj = 1, jpj
598            DO ji = 1, jpi
599               ! Increment the temperature field for computing night mean and counter
600               surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj)  &
601                      &                    + psurf(ji,jj) * REAL( imask_night(ji,jj) )
602               zmeanday(ji,jj)          = zmeanday(ji,jj) + psurf(ji,jj)
603               icount_night(ji,jj)      = icount_night(ji,jj) + imask_night(ji,jj)
604            END DO
605         END DO
606
607         ! Compute the night-time mean at the end of the day
608         zdaystp = 1.0 / REAL( kdaystp )
609         IF ( idayend == 0 ) THEN
610            IF (lwp) WRITE(numout,*) 'Calculating surfdataqc%vdmean on time-step: ',kt
611            DO jj = 1, jpj
612               DO ji = 1, jpi
613                  ! Test if "no night" point
614                  IF ( icount_night(ji,jj) > 0 ) THEN
615                     surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj) &
616                       &                        / REAL( icount_night(ji,jj) )
617                  ELSE
618                     !At locations where there is no night (e.g. poles),
619                     ! calculate daily mean instead of night-time mean.
620                     surfdataqc%vdmean(ji,jj) = zmeanday(ji,jj) * zdaystp
621                  ENDIF
622               END DO
623            END DO
624         ENDIF
625
626      ENDIF
627
628      ! Get the data for interpolation
629
630      ALLOCATE( &
631         & igrdi(2,2,isurf), &
632         & igrdj(2,2,isurf), &
633         & zglam(2,2,isurf), &
634         & zgphi(2,2,isurf), &
635         & zmask(2,2,isurf), &
636         & zsurf(2,2,isurf)  &
637         & )
638
639      DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf
640         iobs = jobs - surfdataqc%nsurfup
641         igrdi(1,1,iobs) = surfdataqc%mi(jobs)-1
642         igrdj(1,1,iobs) = surfdataqc%mj(jobs)-1
643         igrdi(1,2,iobs) = surfdataqc%mi(jobs)-1
644         igrdj(1,2,iobs) = surfdataqc%mj(jobs)
645         igrdi(2,1,iobs) = surfdataqc%mi(jobs)
646         igrdj(2,1,iobs) = surfdataqc%mj(jobs)-1
647         igrdi(2,2,iobs) = surfdataqc%mi(jobs)
648         igrdj(2,2,iobs) = surfdataqc%mj(jobs)
649      END DO
650
651      CALL obs_int_comm_2d( 2, 2, isurf, &
652         &                  igrdi, igrdj, glamt, zglam )
653      CALL obs_int_comm_2d( 2, 2, isurf, &
654         &                  igrdi, igrdj, gphit, zgphi )
655      CALL obs_int_comm_2d( 2, 2, isurf, &
656         &                  igrdi, igrdj, psurfmask, zmask )
657      CALL obs_int_comm_2d( 2, 2, isurf, &
658         &                  igrdi, igrdj, psurf, zsurf )
659
660      ! At the end of the day get interpolated means
661      IF ( idayend == 0 .AND. ldnightav ) THEN
662
663         ALLOCATE( &
664            & zsurfm(2,2,isurf)  &
665            & )
666
667         CALL obs_int_comm_2d( 2, 2, isurf, igrdi, igrdj, &
668            &               surfdataqc%vdmean(:,:), zsurfm )
669
670      ENDIF
671
672      ! Loop over observations
673      DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf
674
675         iobs = jobs - surfdataqc%nsurfup
676
677         IF ( kt /= surfdataqc%mstp(jobs) ) THEN
678
679            IF(lwp) THEN
680               WRITE(numout,*)
681               WRITE(numout,*) ' E R R O R : Observation',              &
682                  &            ' time step is not consistent with the', &
683                  &            ' model time step'
684               WRITE(numout,*) ' ========='
685               WRITE(numout,*)
686               WRITE(numout,*) ' Record  = ', jobs,                &
687                  &            ' kt      = ', kt,                  &
688                  &            ' mstp    = ', surfdataqc%mstp(jobs), &
689                  &            ' ntyp    = ', surfdataqc%ntyp(jobs)
690            ENDIF
691            CALL ctl_stop( 'obs_surf_opt', 'Inconsistent time' )
692
693         ENDIF
694
695         zlam = surfdataqc%rlam(jobs)
696         zphi = surfdataqc%rphi(jobs)
697
698         ! Get weights to interpolate the model value to the observation point
699         CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         &
700            &                   zglam(:,:,iobs), zgphi(:,:,iobs), &
701            &                   zmask(:,:,iobs), zweig, zobsmask )
702
703         ! Interpolate the model field to the observation point
704         IF ( ldnightav .AND. idayend == 0 ) THEN
705            ! Night-time averaged data
706            CALL obs_int_h2d( 1, 1, zweig, zsurfm(:,:,iobs), zext )
707         ELSE
708            CALL obs_int_h2d( 1, 1, zweig, zsurf(:,:,iobs),  zext )
709         ENDIF
710
711         IF ( TRIM(surfdataqc%cvars(1)) == 'SLA' .AND. surfdataqc%nextra == 2 ) THEN
712            ! ... Remove the MDT from the SSH at the observation point to get the SLA
713            surfdataqc%rext(jobs,1) = zext(1)
714            surfdataqc%rmod(jobs,1) = surfdataqc%rext(jobs,1) - surfdataqc%rext(jobs,2)
715         ELSE
716            surfdataqc%rmod(jobs,1) = zext(1)
717         ENDIF
718
719      END DO
720
721      ! Deallocate the data for interpolation
722      DEALLOCATE( &
723         & igrdi, &
724         & igrdj, &
725         & zglam, &
726         & zgphi, &
727         & zmask, &
728         & zsurf  &
729         & )
730
731      ! At the end of the day also deallocate night-time mean array
732      IF ( idayend == 0 .AND. ldnightav ) THEN
733         DEALLOCATE( &
734            & zsurfm  &
735            & )
736      ENDIF
737
738      surfdataqc%nsurfup = surfdataqc%nsurfup + isurf
739
740   END SUBROUTINE obs_surf_opt
741
742END MODULE obs_oper
Note: See TracBrowser for help on using the repository browser.