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

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

source: NEMO/branches/UKMO/NEMO_4.0.4_generic_obs/src/OCE/OBS/obs_oper.F90 @ 15285

Last change on this file since 15285 was 15285, checked in by dford, 9 months ago

Some bug fixes.

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