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/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

source: branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90 @ 12140

Last change on this file since 12140 was 12140, checked in by kingr, 4 years ago

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