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/NEMOGCM/NEMO/OCE_SRC/OBS – NEMO

source: NEMO/trunk/NEMOGCM/NEMO/OCE_SRC/OBS/obs_oper.F90 @ 9594

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

Global renaming for core routines (./NEMO)

  • Folders
    • LIM_SRC_3 -> ICE_SRC
    • OPA_SRC -> OCE_SRC
  • CPP key: key_lim3 -> key_si3
  • Modules, (sub)routines and variables names
    • MPI: mpi_comm_opa -> mpi_comm_oce, MPI_COMM_OPA -> MPI_COMM_OCE, mpi_init_opa -> mpi_init_oce
    • AGRIF: agrif_opa_* -> agrif_oce_*, agrif_lim3_* -> agrif_si3_* and few more
    • TOP-PISCES: p.zlim -> p.zice, namp.zlim -> namp.zice
  • Comments
    • NEMO/OPA -> NEMO/OCE
    • ESIM|LIM3 -> SI3
  • Property svn:keywords set to Id
File size: 38.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 licence (NEMOGCM/NEMO_CeCILL.txt)
37   !!----------------------------------------------------------------------
38CONTAINS
39
40   SUBROUTINE obs_prof_opt( prodatqc, kt, kpi, kpj, kpk,          &
41      &                     kit000, kdaystp,                      &
42      &                     pvar1, pvar2, pgdept, pgdepw,         &
43      &                     pmask1, pmask2,                       & 
44      &                     plam1, plam2, pphi1, pphi2,           &
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      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
111
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
123      INTEGER ::   iin, ijn, ikn, ik   ! looping indices over interpolation nodes
124      INTEGER ::   inum_obs
125      INTEGER, DIMENSION(imaxavtypes) :: &
126         & idailyavtypes
127      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
128         & igrdi1, &
129         & igrdi2, &
130         & igrdj1, &
131         & igrdj2
132      INTEGER, ALLOCATABLE, DIMENSION(:) :: iv_indic
133
134      REAL(KIND=wp) :: zlam
135      REAL(KIND=wp) :: zphi
136      REAL(KIND=wp) :: zdaystp
137      REAL(KIND=wp), DIMENSION(kpk) :: &
138         & zobsmask1, &
139         & zobsmask2, &
140         & zobsk,    &
141         & zobs2k
142      REAL(KIND=wp), DIMENSION(2,2,1) :: &
143         & zweig1, &
144         & zweig2, &
145         & zweig
146      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
147         & zmask1, &
148         & zmask2, &
149         & zint1,  &
150         & zint2,  &
151         & zinm1,  &
152         & zinm2,  &
153         & zgdept, & 
154         & zgdepw
155      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
156         & zglam1, &
157         & zglam2, &
158         & zgphi1, &
159         & zgphi2
160      REAL(KIND=wp), DIMENSION(1) :: zmsk_1, zmsk_2   
161      REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner
162
163      LOGICAL :: ld_dailyav
164
165      !------------------------------------------------------------------------
166      ! Local initialization
167      !------------------------------------------------------------------------
168      ! Record and data counters
169      inrc = kt - kit000 + 2
170      ipro = prodatqc%npstp(inrc)
171
172      ! Daily average types
173      ld_dailyav = .FALSE.
174      IF ( PRESENT(kdailyavtypes) ) THEN
175         idailyavtypes(:) = kdailyavtypes(:)
176         IF ( ANY (idailyavtypes(:) /= -1) ) ld_dailyav = .TRUE.
177      ELSE
178         idailyavtypes(:) = -1
179      ENDIF
180
181      ! Daily means are calculated for values over timesteps:
182      !  [1 <= kt <= kdaystp], [kdaystp+1 <= kt <= 2*kdaystp], ...
183      idayend = MOD( kt - kit000 + 1, kdaystp )
184
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
195               END DO
196            END DO
197         ENDIF
198
199         DO jk = 1, jpk
200            DO jj = 1, jpj
201               DO ji = 1, jpi
202                  ! Increment field 1 for computing daily mean
203                  prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &
204                     &                        + pvar1(ji,jj,jk)
205                  ! Increment field 2 for computing daily mean
206                  prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &
207                     &                        + pvar2(ji,jj,jk)
208               END DO
209            END DO
210         END DO
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
229      ENDIF
230
231      ! Get the data for interpolation
232      ALLOCATE( &
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),  &
243         & zint1(2,2,kpk,ipro),   &
244         & zint2(2,2,kpk,ipro),   &
245         & zgdept(2,2,kpk,ipro),  & 
246         & zgdepw(2,2,kpk,ipro)   & 
247         & )
248
249      DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
250         iobs = jobs - prodatqc%nprofup
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)
267      END DO
268
269      ! Initialise depth arrays
270      zgdept(:,:,:,:) = 0.0
271      zgdepw(:,:,:,:) = 0.0
272
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 )
282
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
286      ! At the end of the day also get interpolated means
287      IF ( ld_dailyav .AND. idayend == 0 ) THEN
288
289         ALLOCATE( &
290            & zinm1(2,2,kpk,ipro),  &
291            & zinm2(2,2,kpk,ipro)   &
292            & )
293
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 )
298
299      ENDIF
300
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
306      DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
307
308         iobs = jobs - prodatqc%nprofup
309
310         IF ( kt /= prodatqc%mstp(jobs) ) THEN
311
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
326
327         zlam = prodatqc%rlam(jobs)
328         zphi = prodatqc%rphi(jobs)
329
330         ! Horizontal weights
331         ! Masked values are calculated later. 
332         IF ( prodatqc%npvend(jobs,1) > 0 ) THEN
333
334            CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,     &
335               &                   zglam1(:,:,iobs), zgphi1(:,:,iobs), &
336               &                   zmask1(:,:,1,iobs), zweig1, zmsk_1 )
337
338         ENDIF
339
340         IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
341
342            CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,     &
343               &                   zglam2(:,:,iobs), zgphi2(:,:,iobs), &
344               &                   zmask2(:,:,1,iobs), zweig2, zmsk_2 )
345 
346         ENDIF
347
348         IF ( prodatqc%npvend(jobs,1) > 0 ) THEN
349
350            zobsk(:) = obfillflt
351
352            IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN
353
354               IF ( idayend == 0 )  THEN
355                  ! Daily averaged data
356
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)) 
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                              &     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) 
377
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 
387
388               ENDIF !idayend
389
390            ELSE   
391
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) 
414
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 
426
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
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
482                  ! Daily averaged data
483
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)) 
489
490                  DO iin=1,2 
491                     DO ijn=1,2 
492
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) 
504
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 
514
515               ENDIF !idayend
516
517            ELSE   
518
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) 
541
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 
553
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 
580
581                        ENDIF
582
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
603      ! Deallocate the data for interpolation
604      DEALLOCATE( &
605         & igrdi1, &
606         & igrdi2, &
607         & igrdj1, &
608         & igrdj2, &
609         & zglam1, &
610         & zglam2, &
611         & zgphi1, &
612         & zgphi2, &
613         & zmask1, &
614         & zmask2, &
615         & zint1,  &
616         & zint2,  &
617         & zgdept, &
618         & zgdepw  &
619         & )
620
621      ! At the end of the day also get interpolated means
622      IF ( ld_dailyav .AND. idayend == 0 ) THEN
623         DEALLOCATE( &
624            & zinm1,  &
625            & zinm2   &
626            & )
627      ENDIF
628
629      prodatqc%nprofup = prodatqc%nprofup + ipro 
630
631   END SUBROUTINE obs_prof_opt
632
633   SUBROUTINE obs_surf_opt( surfdataqc, kt, kpi, kpj,            &
634      &                     kit000, kdaystp, psurf, psurfmask,   &
635      &                     k2dint, ldnightav, plamscl, pphiscl, &
636      &                     lindegrees )
637
638      !!-----------------------------------------------------------------------
639      !!
640      !!                     ***  ROUTINE obs_surf_opt  ***
641      !!
642      !! ** Purpose : Compute the model counterpart of surface
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      !!
649      !!    The new model value is first computed at the obs (lon, lat) point.
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      !!
658      !!    Two horizontal averaging schemes are also available:
659      !!        - weighted radial footprint        (k2dint = 5)
660      !!        - weighted rectangular footprint   (k2dint = 6)
661      !!
662      !!
663      !! ** Action  :
664      !!
665      !! History :
666      !!      ! 07-03 (A. Weaver)
667      !!      ! 15-02 (M. Martin) Combined routine for surface types
668      !!      ! 17-03 (M. Martin) Added horizontal averaging options
669      !!-----------------------------------------------------------------------
670      USE obs_surf_def  ! Definition of storage space for surface observations
671
672      IMPLICIT NONE
673
674      TYPE(obs_surf), INTENT(INOUT) :: &
675         & surfdataqc                  ! Subset of surface data passing QC
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)
681      INTEGER, INTENT(IN) :: kdaystp   ! Number of time steps per day
682      INTEGER, INTENT(IN) :: k2dint    ! Horizontal interpolation type (see header)
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
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
692
693      !! * Local declarations
694      INTEGER :: ji
695      INTEGER :: jj
696      INTEGER :: jobs
697      INTEGER :: inrc
698      INTEGER :: isurf
699      INTEGER :: iobs
700      INTEGER :: imaxifp, imaxjfp
701      INTEGER :: imodi, imodj
702      INTEGER :: idayend
703      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
704         & igrdi,   &
705         & igrdj,   &
706         & igrdip1, &
707         & igrdjp1
708      INTEGER, DIMENSION(:,:), SAVE, ALLOCATABLE :: &
709         & icount_night,      &
710         & imask_night
711      REAL(wp) :: zlam
712      REAL(wp) :: zphi
713      REAL(wp), DIMENSION(1) :: zext, zobsmask
714      REAL(wp) :: zdaystp
715      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
716         & zweig,  &
717         & zmask,  &
718         & zsurf,  &
719         & zsurfm, &
720         & zsurftmp, &
721         & zglam,  &
722         & zgphi,  &
723         & zglamf, &
724         & zgphif
725
726      REAL(wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: &
727         & zintmp,  &
728         & zouttmp, &
729         & zmeanday    ! to compute model sst in region of 24h daylight (pole)
730
731      !------------------------------------------------------------------------
732      ! Local initialization
733      !------------------------------------------------------------------------
734      ! Record and data counters
735      inrc = kt - kit000 + 2
736      isurf = surfdataqc%nsstp(inrc)
737
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
744      IF ( ldnightav ) THEN
745
746      ! Initialize array for night mean
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
755
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 )
759
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
768            END DO
769         ENDIF
770
771         zintmp(:,:) = 0.0
772         zouttmp(:,:) = sbc_dcy( zintmp(:,:), .TRUE. )
773         imask_night(:,:) = INT( zouttmp(:,:) )
774
775         DO jj = 1, jpj
776            DO ji = 1, jpi
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)
782            END DO
783         END DO
784
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
804      ENDIF
805
806      ! Get the data for interpolation
807
808      ALLOCATE( &
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) &
821         & )
822
823      DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf
824         iobs = jobs - surfdataqc%nsurfup
825         DO ji = 0, imaxifp
826            imodi = surfdataqc%mi(jobs) - int(imaxifp/2) + ji - 1
827            !
828            !Deal with wrap around in longitude
829            IF ( imodi < 1      ) imodi = imodi + jpiglo
830            IF ( imodi > jpiglo ) imodi = imodi - jpiglo
831            !
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
838               !
839               igrdip1(ji+1,jj+1,iobs) = imodi
840               igrdjp1(ji+1,jj+1,iobs) = imodj
841               !
842               IF ( ji >= 1 .AND. jj >= 1 ) THEN
843                  igrdi(ji,jj,iobs) = imodi
844                  igrdj(ji,jj,iobs) = imodj
845               ENDIF
846               !
847            END DO
848         END DO
849      END DO
850
851      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
852         &                  igrdi, igrdj, glamt, zglam )
853      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
854         &                  igrdi, igrdj, gphit, zgphi )
855      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
856         &                  igrdi, igrdj, psurfmask, zmask )
857      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
858         &                  igrdi, igrdj, psurf, zsurf )
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 )
863
864      ! At the end of the day get interpolated means
865      IF ( idayend == 0 .AND. ldnightav ) THEN
866
867         ALLOCATE( &
868            & zsurfm(imaxifp,imaxjfp,isurf)  &
869            & )
870
871         CALL obs_int_comm_2d( imaxifp,imaxjfp, isurf, kpi, kpj, igrdi, igrdj, &
872            &               surfdataqc%vdmean(:,:), zsurfm )
873
874      ENDIF
875
876      ! Loop over observations
877      DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf
878
879         iobs = jobs - surfdataqc%nsurfup
880
881         IF ( kt /= surfdataqc%mstp(jobs) ) THEN
882
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,                  &
892                  &            ' mstp    = ', surfdataqc%mstp(jobs), &
893                  &            ' ntyp    = ', surfdataqc%ntyp(jobs)
894            ENDIF
895            CALL ctl_stop( 'obs_surf_opt', 'Inconsistent time' )
896
897         ENDIF
898
899         zlam = surfdataqc%rlam(jobs)
900         zphi = surfdataqc%rphi(jobs)
901
902         IF ( ldnightav .AND. idayend == 0 ) THEN
903            ! Night-time averaged data
904            zsurftmp(:,:,iobs) = zsurfm(:,:,iobs)
905         ELSE
906            zsurftmp(:,:,iobs) = zsurf(:,:,iobs)
907         ENDIF
908
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
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)
938         ELSE
939            surfdataqc%rmod(jobs,1) = zext(1)
940         ENDIF
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
946
947      END DO
948
949      ! Deallocate the data for interpolation
950      DEALLOCATE( &
951         & zweig, &
952         & igrdi, &
953         & igrdj, &
954         & zglam, &
955         & zgphi, &
956         & zmask, &
957         & zsurf, &
958         & zsurftmp, &
959         & zglamf, &
960         & zgphif, &
961         & igrdip1,&
962         & igrdjp1 &
963         & )
964
965      ! At the end of the day also deallocate night-time mean array
966      IF ( idayend == 0 .AND. ldnightav ) THEN
967         DEALLOCATE( &
968            & zsurfm  &
969            & )
970      ENDIF
971      !
972      surfdataqc%nsurfup = surfdataqc%nsurfup + isurf
973      !
974   END SUBROUTINE obs_surf_opt
975
976   !!======================================================================
977END MODULE obs_oper
Note: See TracBrowser for help on using the repository browser.