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/trunk/src/OCE/OBS – NEMO

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

Last change on this file since 10068 was 10068, checked in by nicolasmartin, 6 years ago

First part of modifications to have a common default header : fix typos and SVN keywords properties

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