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

source: branches/2017/dev_r8657_UKMO_OBSoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90 @ 8667

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

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