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

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

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90 @ 9023

Last change on this file since 9023 was 9023, checked in by timgraham, 6 years ago

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