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

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

source: branches/UKMO/dev_r5518_obs_oper_update_bgc3d/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90 @ 9186

Last change on this file since 9186 was 9186, checked in by dford, 6 years ago

Initial implementation of 3D biogeochemistry observation operator.

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