New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
obs_oper.F90 in NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/OBS – NEMO

source: NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/OBS/obs_oper.F90 @ 12928

Last change on this file since 12928 was 12928, checked in by smueller, 4 years ago

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