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 NEMO/branches/UKMO/NEMO_4.0.4_generic_obs/src/OCE/OBS – NEMO

source: NEMO/branches/UKMO/NEMO_4.0.4_generic_obs/src/OCE/OBS/obs_oper.F90 @ 15144

Last change on this file since 15144 was 15144, checked in by dford, 3 years ago

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