source: NEMO/trunk/src/OCE/OBS/obs_oper.F90

Last change on this file was 14056, checked in by ayoung, 5 months ago

Adding branch for ticket #2567 to trunk.

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