source: branches/UKMO/dev_r5518_obs_oper_update_medusa/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90 @ 8646

Last change on this file since 8646 was 8646, checked in by dford, 3 years ago

Add comments in response to review.

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
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                     prodatqc%vdmean(ji,jj,jk,2) = 0.0
234                  END DO
235               END DO
236            END DO
237         ENDIF
238
239         DO jk = 1, jpk
240            DO jj = 1, jpj
241               DO ji = 1, jpi
242                  ! Increment field 1 for computing daily mean
243                  prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &
244                     &                        + pvar1(ji,jj,jk)
245                  ! Increment field 2 for computing daily mean
246                  prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &
247                     &                        + pvar2(ji,jj,jk)
248               END DO
249            END DO
250         END DO
251
252         ! Compute the daily mean at the end of day
253         zdaystp = 1.0 / REAL( kdaystp )
254         IF ( idayend == 0 ) THEN
255            IF (lwp) WRITE(numout,*) 'Calculating prodatqc%vdmean on time-step: ',kt
256            CALL FLUSH(numout)
257            DO jk = 1, jpk
258               DO jj = 1, jpj
259                  DO ji = 1, jpi
260                     prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &
261                        &                        * zdaystp
262                     prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &
263                        &                        * zdaystp
264                  END DO
265               END DO
266            END DO
267         ENDIF
268
269      ENDIF
270
271      ! Get the data for interpolation
272      ALLOCATE( &
273         & igrdi1(2,2,ipro),      &
274         & igrdi2(2,2,ipro),      &
275         & igrdj1(2,2,ipro),      &
276         & igrdj2(2,2,ipro),      &
277         & zglam1(2,2,ipro),      &
278         & zglam2(2,2,ipro),      &
279         & zgphi1(2,2,ipro),      &
280         & zgphi2(2,2,ipro),      &
281         & zmask1(2,2,kpk,ipro),  &
282         & zmask2(2,2,kpk,ipro),  &
283         & zint1(2,2,kpk,ipro),   &
284         & zint2(2,2,kpk,ipro),   &
285         & zgdept(2,2,kpk,ipro),  & 
286         & zgdepw(2,2,kpk,ipro)   & 
287         & )
288
289      DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
290         iobs = jobs - prodatqc%nprofup
291         igrdi1(1,1,iobs) = prodatqc%mi(jobs,1)-1
292         igrdj1(1,1,iobs) = prodatqc%mj(jobs,1)-1
293         igrdi1(1,2,iobs) = prodatqc%mi(jobs,1)-1
294         igrdj1(1,2,iobs) = prodatqc%mj(jobs,1)
295         igrdi1(2,1,iobs) = prodatqc%mi(jobs,1)
296         igrdj1(2,1,iobs) = prodatqc%mj(jobs,1)-1
297         igrdi1(2,2,iobs) = prodatqc%mi(jobs,1)
298         igrdj1(2,2,iobs) = prodatqc%mj(jobs,1)
299         igrdi2(1,1,iobs) = prodatqc%mi(jobs,2)-1
300         igrdj2(1,1,iobs) = prodatqc%mj(jobs,2)-1
301         igrdi2(1,2,iobs) = prodatqc%mi(jobs,2)-1
302         igrdj2(1,2,iobs) = prodatqc%mj(jobs,2)
303         igrdi2(2,1,iobs) = prodatqc%mi(jobs,2)
304         igrdj2(2,1,iobs) = prodatqc%mj(jobs,2)-1
305         igrdi2(2,2,iobs) = prodatqc%mi(jobs,2)
306         igrdj2(2,2,iobs) = prodatqc%mj(jobs,2)
307      END DO
308
309      ! Initialise depth arrays
310      zgdept(:,:,:,:) = 0.0
311      zgdepw(:,:,:,:) = 0.0
312
313      CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi1, igrdj1, plam1, zglam1 )
314      CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi1, igrdj1, pphi1, zgphi1 )
315      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pmask1, zmask1 )
316      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pvar1,   zint1 )
317     
318      CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi2, igrdj2, plam2, zglam2 )
319      CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi2, igrdj2, pphi2, zgphi2 )
320      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, pmask2, zmask2 )
321      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, pvar2,   zint2 )
322
323      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pgdept, zgdept ) 
324      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pgdepw, zgdepw ) 
325
326      ! At the end of the day also get interpolated means
327      IF ( ld_dailyav .AND. idayend == 0 ) THEN
328
329         ALLOCATE( &
330            & zinm1(2,2,kpk,ipro),  &
331            & zinm2(2,2,kpk,ipro)   &
332            & )
333
334         CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, &
335            &                  prodatqc%vdmean(:,:,:,1), zinm1 )
336         CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, &
337            &                  prodatqc%vdmean(:,:,:,2), zinm2 )
338
339      ENDIF
340
341      ! Return if no observations to process
342      ! Has to be done after comm commands to ensure processors
343      ! stay in sync
344      IF ( ipro == 0 ) RETURN
345
346      DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
347
348         iobs = jobs - prodatqc%nprofup
349
350         IF ( kt /= prodatqc%mstp(jobs) ) THEN
351
352            IF(lwp) THEN
353               WRITE(numout,*)
354               WRITE(numout,*) ' E R R O R : Observation',              &
355                  &            ' time step is not consistent with the', &
356                  &            ' model time step'
357               WRITE(numout,*) ' ========='
358               WRITE(numout,*)
359               WRITE(numout,*) ' Record  = ', jobs,                    &
360                  &            ' kt      = ', kt,                      &
361                  &            ' mstp    = ', prodatqc%mstp(jobs), &
362                  &            ' ntyp    = ', prodatqc%ntyp(jobs)
363            ENDIF
364            CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' )
365         ENDIF
366
367         zlam = prodatqc%rlam(jobs)
368         zphi = prodatqc%rphi(jobs)
369
370         ! Horizontal weights
371         ! Masked values are calculated later. 
372         IF ( prodatqc%npvend(jobs,1) > 0 ) THEN
373
374            CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,     &
375               &                   zglam1(:,:,iobs), zgphi1(:,:,iobs), &
376               &                   zmask1(:,:,1,iobs), zweig1, zmsk_1 )
377
378         ENDIF
379
380         IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
381
382            CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,     &
383               &                   zglam2(:,:,iobs), zgphi2(:,:,iobs), &
384               &                   zmask2(:,:,1,iobs), zweig2, zmsk_2 )
385 
386         ENDIF
387
388         IF ( prodatqc%npvend(jobs,1) > 0 ) THEN
389
390            zobsk(:) = obfillflt
391
392            IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN
393
394               IF ( idayend == 0 )  THEN
395                  ! Daily averaged data
396
397                  ! vertically interpolate all 4 corners
398                  ista = prodatqc%npvsta(jobs,1) 
399                  iend = prodatqc%npvend(jobs,1) 
400                  inum_obs = iend - ista + 1 
401                  ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs)) 
402
403                  DO iin=1,2 
404                     DO ijn=1,2 
405
406                        IF ( k1dint == 1 ) THEN
407                           CALL obs_int_z1d_spl( kpk, & 
408                              &     zinm1(iin,ijn,:,iobs), & 
409                              &     zobs2k, zgdept(iin,ijn,:,iobs), & 
410                              &     zmask1(iin,ijn,:,iobs)) 
411                        ENDIF
412       
413                        CALL obs_level_search(kpk, & 
414                           &    zgdept(iin,ijn,:,iobs), & 
415                           &    inum_obs, prodatqc%var(1)%vdep(ista:iend), & 
416                           &    iv_indic) 
417
418                        CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 
419                           &    prodatqc%var(1)%vdep(ista:iend), & 
420                           &    zinm1(iin,ijn,:,iobs), & 
421                           &    zobs2k, interp_corner(iin,ijn,:), & 
422                           &    zgdept(iin,ijn,:,iobs), & 
423                           &    zmask1(iin,ijn,:,iobs)) 
424       
425                     ENDDO 
426                  ENDDO 
427
428               ENDIF !idayend
429
430            ELSE   
431
432               ! Point data
433     
434               ! vertically interpolate all 4 corners
435               ista = prodatqc%npvsta(jobs,1) 
436               iend = prodatqc%npvend(jobs,1) 
437               inum_obs = iend - ista + 1 
438               ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs)) 
439               DO iin=1,2 
440                  DO ijn=1,2 
441                   
442                     IF ( k1dint == 1 ) THEN
443                        CALL obs_int_z1d_spl( kpk, & 
444                           &    zint1(iin,ijn,:,iobs),& 
445                           &    zobs2k, zgdept(iin,ijn,:,iobs), & 
446                           &    zmask1(iin,ijn,:,iobs)) 
447 
448                     ENDIF
449       
450                     CALL obs_level_search(kpk, & 
451                         &        zgdept(iin,ijn,:,iobs),& 
452                         &        inum_obs, prodatqc%var(1)%vdep(ista:iend), & 
453                         &        iv_indic) 
454
455                     CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs,     & 
456                         &          prodatqc%var(1)%vdep(ista:iend),     & 
457                         &          zint1(iin,ijn,:,iobs),            & 
458                         &          zobs2k,interp_corner(iin,ijn,:), & 
459                         &          zgdept(iin,ijn,:,iobs),         & 
460                         &          zmask1(iin,ijn,:,iobs) )     
461         
462                  ENDDO 
463               ENDDO 
464             
465            ENDIF 
466
467            !-------------------------------------------------------------
468            ! Compute the horizontal interpolation for every profile level
469            !-------------------------------------------------------------
470             
471            DO ikn=1,inum_obs 
472               iend=ista+ikn-1
473                 
474               zweig(:,:,1) = 0._wp 
475   
476               ! This code forces the horizontal weights to be 
477               ! zero IF the observation is below the bottom of the 
478               ! corners of the interpolation nodes, Or if it is in 
479               ! the mask. This is important for observations near 
480               ! steep bathymetry
481               DO iin=1,2 
482                  DO ijn=1,2 
483     
484                     depth_loop1: DO ik=kpk,2,-1 
485                        IF(zmask1(iin,ijn,ik-1,iobs ) > 0.9 )THEN   
486                           
487                           zweig(iin,ijn,1) = & 
488                              & zweig1(iin,ijn,1) * & 
489                              & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) & 
490                              &  - prodatqc%var(1)%vdep(iend)),0._wp) 
491                           
492                           EXIT depth_loop1 
493
494                        ENDIF
495
496                     ENDDO depth_loop1 
497     
498                  ENDDO 
499               ENDDO 
500   
501               CALL obs_int_h2d( 1, 1, zweig, interp_corner(:,:,ikn), & 
502                  &              prodatqc%var(1)%vmod(iend:iend) ) 
503
504                  ! Set QC flag for any observations found below the bottom
505                  ! needed as the check here is more strict than that in obs_prep
506               IF (sum(zweig) == 0.0_wp) prodatqc%var(1)%nvqc(iend:iend)=4
507 
508            ENDDO 
509 
510            DEALLOCATE(interp_corner,iv_indic) 
511         
512         ENDIF 
513
514         ! For the second variable
515         IF ( prodatqc%npvend(jobs,2) > 0 ) THEN
516
517            zobsk(:) = obfillflt
518
519            IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN
520
521               IF ( idayend == 0 )  THEN
522                  ! Daily averaged data
523
524                  ! vertically interpolate all 4 corners
525                  ista = prodatqc%npvsta(jobs,2) 
526                  iend = prodatqc%npvend(jobs,2) 
527                  inum_obs = iend - ista + 1 
528                  ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs)) 
529
530                  DO iin=1,2 
531                     DO ijn=1,2 
532
533                        IF ( k1dint == 1 ) THEN
534                           CALL obs_int_z1d_spl( kpk, & 
535                              &     zinm2(iin,ijn,:,iobs), & 
536                              &     zobs2k, zgdept(iin,ijn,:,iobs), & 
537                              &     zmask2(iin,ijn,:,iobs)) 
538                        ENDIF
539       
540                        CALL obs_level_search(kpk, & 
541                           &    zgdept(iin,ijn,:,iobs), & 
542                           &    inum_obs, prodatqc%var(2)%vdep(ista:iend), & 
543                           &    iv_indic) 
544
545                        CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 
546                           &    prodatqc%var(2)%vdep(ista:iend), & 
547                           &    zinm2(iin,ijn,:,iobs), & 
548                           &    zobs2k, interp_corner(iin,ijn,:), & 
549                           &    zgdept(iin,ijn,:,iobs), & 
550                           &    zmask2(iin,ijn,:,iobs)) 
551       
552                     ENDDO 
553                  ENDDO 
554
555               ENDIF !idayend
556
557            ELSE   
558
559               ! Point data
560     
561               ! vertically interpolate all 4 corners
562               ista = prodatqc%npvsta(jobs,2) 
563               iend = prodatqc%npvend(jobs,2) 
564               inum_obs = iend - ista + 1 
565               ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs)) 
566               DO iin=1,2 
567                  DO ijn=1,2 
568                   
569                     IF ( k1dint == 1 ) THEN
570                        CALL obs_int_z1d_spl( kpk, & 
571                           &    zint2(iin,ijn,:,iobs),& 
572                           &    zobs2k, zgdept(iin,ijn,:,iobs), & 
573                           &    zmask2(iin,ijn,:,iobs)) 
574 
575                     ENDIF
576       
577                     CALL obs_level_search(kpk, & 
578                         &        zgdept(iin,ijn,:,iobs),& 
579                         &        inum_obs, prodatqc%var(2)%vdep(ista:iend), & 
580                         &        iv_indic) 
581
582                     CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs,     & 
583                         &          prodatqc%var(2)%vdep(ista:iend),     & 
584                         &          zint2(iin,ijn,:,iobs),            & 
585                         &          zobs2k,interp_corner(iin,ijn,:), & 
586                         &          zgdept(iin,ijn,:,iobs),         & 
587                         &          zmask2(iin,ijn,:,iobs) )     
588         
589                  ENDDO 
590               ENDDO 
591             
592            ENDIF 
593
594            !-------------------------------------------------------------
595            ! Compute the horizontal interpolation for every profile level
596            !-------------------------------------------------------------
597             
598            DO ikn=1,inum_obs 
599               iend=ista+ikn-1
600                 
601               zweig(:,:,1) = 0._wp 
602   
603               ! This code forces the horizontal weights to be 
604               ! zero IF the observation is below the bottom of the 
605               ! corners of the interpolation nodes, Or if it is in 
606               ! the mask. This is important for observations near 
607               ! steep bathymetry
608               DO iin=1,2 
609                  DO ijn=1,2 
610     
611                     depth_loop2: DO ik=kpk,2,-1 
612                        IF(zmask2(iin,ijn,ik-1,iobs ) > 0.9 )THEN   
613                           
614                           zweig(iin,ijn,1) = & 
615                              & zweig2(iin,ijn,1) * & 
616                              & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) & 
617                              &  - prodatqc%var(2)%vdep(iend)),0._wp) 
618                           
619                           EXIT depth_loop2 
620
621                        ENDIF
622
623                     ENDDO depth_loop2 
624     
625                  ENDDO 
626               ENDDO 
627   
628               CALL obs_int_h2d( 1, 1, zweig, interp_corner(:,:,ikn), & 
629                  &              prodatqc%var(2)%vmod(iend:iend) ) 
630
631                  ! Set QC flag for any observations found below the bottom
632                  ! needed as the check here is more strict than that in obs_prep
633               IF (sum(zweig) == 0.0_wp) prodatqc%var(2)%nvqc(iend:iend)=4
634 
635            ENDDO 
636 
637            DEALLOCATE(interp_corner,iv_indic) 
638         
639         ENDIF
640
641      ENDDO
642
643      ! Deallocate the data for interpolation
644      DEALLOCATE( &
645         & igrdi1, &
646         & igrdi2, &
647         & igrdj1, &
648         & igrdj2, &
649         & zglam1, &
650         & zglam2, &
651         & zgphi1, &
652         & zgphi2, &
653         & zmask1, &
654         & zmask2, &
655         & zint1,  &
656         & zint2,  &
657         & zgdept, &
658         & zgdepw  &
659         & )
660
661      ! At the end of the day also get interpolated means
662      IF ( ld_dailyav .AND. idayend == 0 ) THEN
663         DEALLOCATE( &
664            & zinm1,  &
665            & zinm2   &
666            & )
667      ENDIF
668
669      prodatqc%nprofup = prodatqc%nprofup + ipro 
670
671   END SUBROUTINE obs_prof_opt
672
673   SUBROUTINE obs_surf_opt( surfdataqc, kt, kpi, kpj,            &
674      &                     kit000, kdaystp, psurf, psurfmask,   &
675      &                     k2dint, ldnightav, plamscl, pphiscl, &
676      &                     lindegrees )
677
678      !!-----------------------------------------------------------------------
679      !!
680      !!                     ***  ROUTINE obs_surf_opt  ***
681      !!
682      !! ** Purpose : Compute the model counterpart of surface
683      !!              data by interpolating from the model grid to the
684      !!              observation point.
685      !!
686      !! ** Method  : Linearly interpolate to each observation point using
687      !!              the model values at the corners of the surrounding grid box.
688      !!
689      !!    The new model value is first computed at the obs (lon, lat) point.
690      !!
691      !!    Several horizontal interpolation schemes are available:
692      !!        - distance-weighted (great circle) (k2dint = 0)
693      !!        - distance-weighted (small angle)  (k2dint = 1)
694      !!        - bilinear (geographical grid)     (k2dint = 2)
695      !!        - bilinear (quadrilateral grid)    (k2dint = 3)
696      !!        - polynomial (quadrilateral grid)  (k2dint = 4)
697      !!
698      !!    Two horizontal averaging schemes are also available:
699      !!        - weighted radial footprint        (k2dint = 5)
700      !!        - weighted rectangular footprint   (k2dint = 6)
701      !!
702      !!
703      !! ** Action  :
704      !!
705      !! History :
706      !!      ! 07-03 (A. Weaver)
707      !!      ! 15-02 (M. Martin) Combined routine for surface types
708      !!      ! 17-03 (M. Martin) Added horizontal averaging options
709      !!-----------------------------------------------------------------------
710
711      !! * Modules used
712      USE obs_surf_def  ! Definition of storage space for surface observations
713
714      IMPLICIT NONE
715
716      !! * Arguments
717      TYPE(obs_surf), INTENT(INOUT) :: &
718         & surfdataqc                  ! Subset of surface data passing QC
719      INTEGER, INTENT(IN) :: kt        ! Time step
720      INTEGER, INTENT(IN) :: kpi       ! Model grid parameters
721      INTEGER, INTENT(IN) :: kpj
722      INTEGER, INTENT(IN) :: kit000    ! Number of the first time step
723                                       !   (kit000-1 = restart time)
724      INTEGER, INTENT(IN) :: kdaystp   ! Number of time steps per day
725      INTEGER, INTENT(IN) :: k2dint    ! Horizontal interpolation type (see header)
726      REAL(wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
727         & psurf,  &                   ! Model surface field
728         & psurfmask                   ! Land-sea mask
729      LOGICAL, INTENT(IN) :: ldnightav ! Logical for averaging night-time data
730      REAL(KIND=wp), INTENT(IN) :: &
731         & plamscl, &                  ! Diameter in metres of obs footprint in E/W, N/S directions
732         & pphiscl                     ! This is the full width (rather than half-width)
733      LOGICAL, INTENT(IN) :: &
734         & lindegrees                  ! T=> plamscl and pphiscl are specified in degrees, F=> in metres
735
736      !! * Local declarations
737      INTEGER :: ji
738      INTEGER :: jj
739      INTEGER :: jobs
740      INTEGER :: inrc
741      INTEGER :: isurf
742      INTEGER :: iobs
743      INTEGER :: imaxifp, imaxjfp
744      INTEGER :: imodi, imodj
745      INTEGER :: idayend
746      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
747         & igrdi,   &
748         & igrdj,   &
749         & igrdip1, &
750         & igrdjp1
751      INTEGER, DIMENSION(:,:), SAVE, ALLOCATABLE :: &
752         & icount_night,      &
753         & imask_night
754      REAL(wp) :: zlam
755      REAL(wp) :: zphi
756      REAL(wp), DIMENSION(1) :: zext, zobsmask
757      REAL(wp) :: zdaystp
758      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
759         & zweig,  &
760         & zmask,  &
761         & zsurf,  &
762         & zsurfm, &
763         & zsurftmp, &
764         & zglam,  &
765         & zgphi,  &
766         & zglamf, &
767         & zgphif
768
769      REAL(wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: &
770         & zintmp,  &
771         & zouttmp, &
772         & zmeanday    ! to compute model sst in region of 24h daylight (pole)
773
774      !------------------------------------------------------------------------
775      ! Local initialization
776      !------------------------------------------------------------------------
777      ! Record and data counters
778      inrc = kt - kit000 + 2
779      isurf = surfdataqc%nsstp(inrc)
780
781      ! Work out the maximum footprint size for the
782      ! interpolation/averaging in model grid-points - has to be even.
783
784      CALL obs_max_fpsize( k2dint, plamscl, pphiscl, lindegrees, psurfmask, imaxifp, imaxjfp )
785
786
787      IF ( ldnightav ) THEN
788
789      ! Initialize array for night mean
790         IF ( kt == 0 ) THEN
791            ALLOCATE ( icount_night(kpi,kpj) )
792            ALLOCATE ( imask_night(kpi,kpj) )
793            ALLOCATE ( zintmp(kpi,kpj) )
794            ALLOCATE ( zouttmp(kpi,kpj) )
795            ALLOCATE ( zmeanday(kpi,kpj) )
796            nday_qsr = -1   ! initialisation flag for nbc_dcy
797         ENDIF
798
799         ! Night-time means are calculated for night-time values over timesteps:
800         !  [1 <= kt <= kdaystp], [kdaystp+1 <= kt <= 2*kdaystp], .....
801         idayend = MOD( kt - kit000 + 1, kdaystp )
802
803         ! Initialize night-time mean for first timestep of the day
804         IF ( idayend == 1 .OR. kt == 0 ) THEN
805            DO jj = 1, jpj
806               DO ji = 1, jpi
807                  surfdataqc%vdmean(ji,jj) = 0.0
808                  zmeanday(ji,jj) = 0.0
809                  icount_night(ji,jj) = 0
810               END DO
811            END DO
812         ENDIF
813
814         zintmp(:,:) = 0.0
815         zouttmp(:,:) = sbc_dcy( zintmp(:,:), .TRUE. )
816         imask_night(:,:) = INT( zouttmp(:,:) )
817
818         DO jj = 1, jpj
819            DO ji = 1, jpi
820               ! Increment the temperature field for computing night mean and counter
821               surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj)  &
822                      &                    + psurf(ji,jj) * REAL( imask_night(ji,jj) )
823               zmeanday(ji,jj)          = zmeanday(ji,jj) + psurf(ji,jj)
824               icount_night(ji,jj)      = icount_night(ji,jj) + imask_night(ji,jj)
825            END DO
826         END DO
827
828         ! Compute the night-time mean at the end of the day
829         zdaystp = 1.0 / REAL( kdaystp )
830         IF ( idayend == 0 ) THEN
831            IF (lwp) WRITE(numout,*) 'Calculating surfdataqc%vdmean on time-step: ',kt
832            DO jj = 1, jpj
833               DO ji = 1, jpi
834                  ! Test if "no night" point
835                  IF ( icount_night(ji,jj) > 0 ) THEN
836                     surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj) &
837                       &                        / REAL( icount_night(ji,jj) )
838                  ELSE
839                     !At locations where there is no night (e.g. poles),
840                     ! calculate daily mean instead of night-time mean.
841                     surfdataqc%vdmean(ji,jj) = zmeanday(ji,jj) * zdaystp
842                  ENDIF
843               END DO
844            END DO
845         ENDIF
846
847      ENDIF
848
849      ! Get the data for interpolation
850
851      ALLOCATE( &
852         & zweig(imaxifp,imaxjfp,1),      &
853         & igrdi(imaxifp,imaxjfp,isurf), &
854         & igrdj(imaxifp,imaxjfp,isurf), &
855         & zglam(imaxifp,imaxjfp,isurf), &
856         & zgphi(imaxifp,imaxjfp,isurf), &
857         & zmask(imaxifp,imaxjfp,isurf), &
858         & zsurf(imaxifp,imaxjfp,isurf), &
859         & zsurftmp(imaxifp,imaxjfp,isurf),  &
860         & zglamf(imaxifp+1,imaxjfp+1,isurf), &
861         & zgphif(imaxifp+1,imaxjfp+1,isurf), &
862         & igrdip1(imaxifp+1,imaxjfp+1,isurf), &
863         & igrdjp1(imaxifp+1,imaxjfp+1,isurf) &
864         & )
865
866      DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf
867         iobs = jobs - surfdataqc%nsurfup
868         DO ji = 0, imaxifp
869            imodi = surfdataqc%mi(jobs) - int(imaxifp/2) + ji - 1
870           
871            !Deal with wrap around in longitude
872            IF ( imodi < 1      ) imodi = imodi + jpiglo
873            IF ( imodi > jpiglo ) imodi = imodi - jpiglo
874           
875            DO jj = 0, imaxjfp
876               imodj = surfdataqc%mj(jobs) - int(imaxjfp/2) + jj - 1
877               !If model values are out of the domain to the north/south then
878               !set them to be the edge of the domain
879               IF ( imodj < 1      ) imodj = 1
880               IF ( imodj > jpjglo ) imodj = jpjglo
881
882               igrdip1(ji+1,jj+1,iobs) = imodi
883               igrdjp1(ji+1,jj+1,iobs) = imodj
884               
885               IF ( ji >= 1 .AND. jj >= 1 ) THEN
886                  igrdi(ji,jj,iobs) = imodi
887                  igrdj(ji,jj,iobs) = imodj
888               ENDIF
889               
890            END DO
891         END DO
892      END DO
893
894      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
895         &                  igrdi, igrdj, glamt, zglam )
896      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
897         &                  igrdi, igrdj, gphit, zgphi )
898      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
899         &                  igrdi, igrdj, psurfmask, zmask )
900      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
901         &                  igrdi, igrdj, psurf, zsurf )
902      CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, &
903         &                  igrdip1, igrdjp1, glamf, zglamf )
904      CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, &
905         &                  igrdip1, igrdjp1, gphif, zgphif )
906
907      ! At the end of the day get interpolated means
908      IF ( idayend == 0 .AND. ldnightav ) THEN
909
910         ALLOCATE( &
911            & zsurfm(imaxifp,imaxjfp,isurf)  &
912            & )
913
914         CALL obs_int_comm_2d( imaxifp,imaxjfp, isurf, kpi, kpj, igrdi, igrdj, &
915            &               surfdataqc%vdmean(:,:), zsurfm )
916
917      ENDIF
918
919      ! Loop over observations
920      DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf
921
922         iobs = jobs - surfdataqc%nsurfup
923
924         IF ( kt /= surfdataqc%mstp(jobs) ) THEN
925
926            IF(lwp) THEN
927               WRITE(numout,*)
928               WRITE(numout,*) ' E R R O R : Observation',              &
929                  &            ' time step is not consistent with the', &
930                  &            ' model time step'
931               WRITE(numout,*) ' ========='
932               WRITE(numout,*)
933               WRITE(numout,*) ' Record  = ', jobs,                &
934                  &            ' kt      = ', kt,                  &
935                  &            ' mstp    = ', surfdataqc%mstp(jobs), &
936                  &            ' ntyp    = ', surfdataqc%ntyp(jobs)
937            ENDIF
938            CALL ctl_stop( 'obs_surf_opt', 'Inconsistent time' )
939
940         ENDIF
941
942         zlam = surfdataqc%rlam(jobs)
943         zphi = surfdataqc%rphi(jobs)
944
945         IF ( ldnightav .AND. idayend == 0 ) THEN
946            ! Night-time averaged data
947            zsurftmp(:,:,iobs) = zsurfm(:,:,iobs)
948         ELSE
949            zsurftmp(:,:,iobs) = zsurf(:,:,iobs)
950         ENDIF
951
952         IF ( k2dint <= 4 ) THEN
953
954            ! Get weights to interpolate the model value to the observation point
955            CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         &
956               &                   zglam(:,:,iobs), zgphi(:,:,iobs), &
957               &                   zmask(:,:,iobs), zweig, zobsmask )
958
959            ! Interpolate the model value to the observation point
960            CALL obs_int_h2d( 1, 1, zweig, zsurftmp(:,:,iobs), zext )
961
962         ELSE
963
964            ! Get weights to average the model SLA to the observation footprint
965            CALL obs_avg_h2d_init( 1, 1, imaxifp, imaxjfp, k2dint, zlam,  zphi, &
966               &                   zglam(:,:,iobs), zgphi(:,:,iobs), &
967               &                   zglamf(:,:,iobs), zgphif(:,:,iobs), &
968               &                   zmask(:,:,iobs), plamscl, pphiscl, &
969               &                   lindegrees, zweig, zobsmask )
970
971            ! Average the model SST to the observation footprint
972            CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, &
973               &              zweig, zsurftmp(:,:,iobs),  zext )
974
975         ENDIF
976
977         IF ( TRIM(surfdataqc%cvars(1)) == 'SLA' .AND. surfdataqc%nextra == 2 ) THEN
978            ! ... Remove the MDT from the SSH at the observation point to get the SLA
979            surfdataqc%rext(jobs,1) = zext(1)
980            surfdataqc%rmod(jobs,1) = surfdataqc%rext(jobs,1) - surfdataqc%rext(jobs,2)
981         ELSE
982            surfdataqc%rmod(jobs,1) = zext(1)
983         ENDIF
984         
985         IF ( zext(1) == obfillflt ) THEN
986            ! If the observation value is a fill value, set QC flag to bad
987            surfdataqc%nqc(jobs) = 4
988         ENDIF
989
990      END DO
991
992      ! Deallocate the data for interpolation
993      DEALLOCATE( &
994         & zweig, &
995         & igrdi, &
996         & igrdj, &
997         & zglam, &
998         & zgphi, &
999         & zmask, &
1000         & zsurf, &
1001         & zsurftmp, &
1002         & zglamf, &
1003         & zgphif, &
1004         & igrdip1,&
1005         & igrdjp1 &
1006         & )
1007
1008      ! At the end of the day also deallocate night-time mean array
1009      IF ( idayend == 0 .AND. ldnightav ) THEN
1010         DEALLOCATE( &
1011            & zsurfm  &
1012            & )
1013      ENDIF
1014
1015      surfdataqc%nsurfup = surfdataqc%nsurfup + isurf
1016
1017   END SUBROUTINE obs_surf_opt
1018
1019END MODULE obs_oper
Note: See TracBrowser for help on using the repository browser.