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 @ 15487

Last change on this file since 15487 was 15487, checked in by dford, 8 months ago

Remove some comments and add some spaces.

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: ', &
681               &                     kt, ' with weight: ', zmeanstp
682            DO jj = 1, jpj
683               DO ji = 1, jpi
684                  surfdataqc%vdmean(ji,jj,kvar) = surfdataqc%vdmean(ji,jj,kvar) &
685                     &                            * zmeanstp
686               END DO
687            END DO
688         ENDIF
689      ENDIF
690
691      IF ( ldnightav ) THEN
692
693         ! Initialize array for night mean
694         IF ( kt == 0 ) THEN
695            ALLOCATE ( icount_night(kpi,kpj) )
696            ALLOCATE ( imask_night(kpi,kpj) )
697            ALLOCATE ( zintmp(kpi,kpj) )
698            ALLOCATE ( zouttmp(kpi,kpj) )
699            ALLOCATE ( zmeanday(kpi,kpj) )
700            nday_qsr = -1   ! initialisation flag for nbc_dcy
701         ENDIF
702
703         ! Night-time means are calculated for night-time values over timesteps:
704         !  [1 <= kt <= kdaystp], [kdaystp+1 <= kt <= 2*kdaystp], .....
705         idayend = MOD( kt - kit000 + 1, kdaystp )
706
707         ! Initialize night-time mean for first timestep of the day
708         IF ( idayend == 1 .OR. kt == 0 ) THEN
709            DO jj = 1, jpj
710               DO ji = 1, jpi
711                  surfdataqc%vdmean(ji,jj,:) = 0.0
712                  zmeanday(ji,jj) = 0.0
713                  icount_night(ji,jj) = 0
714               END DO
715            END DO
716         ENDIF
717
718         zintmp(:,:) = 0.0
719         zouttmp(:,:) = sbc_dcy( zintmp(:,:), .TRUE. )
720         imask_night(:,:) = INT( zouttmp(:,:) )
721
722         DO jj = 1, jpj
723            DO ji = 1, jpi
724               ! Increment the model field for computing night mean and counter
725               surfdataqc%vdmean(ji,jj,kvar) = surfdataqc%vdmean(ji,jj,kvar)  &
726                      &                        + psurf(ji,jj) * REAL( imask_night(ji,jj) )
727               zmeanday(ji,jj)          = zmeanday(ji,jj) + psurf(ji,jj)
728               icount_night(ji,jj)      = icount_night(ji,jj) + imask_night(ji,jj)
729            END DO
730         END DO
731
732         ! Compute the night-time mean at the end of the day
733         zdaystp = 1.0 / REAL( kdaystp )
734         IF ( idayend == 0 ) THEN
735            IF (lwp) WRITE(numout,*) 'Calculating surfdataqc%vdmean on time-step: ', kt
736            DO jj = 1, jpj
737               DO ji = 1, jpi
738                  ! Test if "no night" point
739                  IF ( icount_night(ji,jj) > 0 ) THEN
740                     surfdataqc%vdmean(ji,jj,kvar) = surfdataqc%vdmean(ji,jj,kvar) &
741                       &                             / REAL( icount_night(ji,jj) )
742                  ELSE
743                     !At locations where there is no night (e.g. poles),
744                     ! calculate daily mean instead of night-time mean.
745                     surfdataqc%vdmean(ji,jj,kvar) = zmeanday(ji,jj) * zdaystp
746                  ENDIF
747               END DO
748            END DO
749         ENDIF
750
751      ENDIF
752
753      ! Get the data for interpolation
754
755      ALLOCATE( &
756         & zweig(imaxifp,imaxjfp,1),      &
757         & igrdi(imaxifp,imaxjfp,isurf), &
758         & igrdj(imaxifp,imaxjfp,isurf), &
759         & zglam(imaxifp,imaxjfp,isurf), &
760         & zgphi(imaxifp,imaxjfp,isurf), &
761         & zmask(imaxifp,imaxjfp,isurf), &
762         & zsurf(imaxifp,imaxjfp,isurf), &
763         & zsurftmp(imaxifp,imaxjfp,isurf) &
764         & )
765
766      IF ( k2dint > 4 ) THEN
767         ALLOCATE( &
768            & zglamf(imaxifp+1,imaxjfp+1,isurf),  &
769            & zgphif(imaxifp+1,imaxjfp+1,isurf),  &
770            & igrdip1(imaxifp+1,imaxjfp+1,isurf), &
771            & igrdjp1(imaxifp+1,imaxjfp+1,isurf)  &
772            & )
773      ENDIF
774
775      IF ( ldclim ) THEN
776         ALLOCATE( zclim(imaxifp,imaxjfp,isurf) )
777      ENDIF
778
779      DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf
780         iobs = jobs - surfdataqc%nsurfup
781         DO ji = 0, imaxifp
782            imodi = surfdataqc%mi(jobs,kvar) - int(imaxifp/2) + ji - 1
783            !
784            !Deal with wrap around in longitude
785            IF ( imodi < 1      ) imodi = imodi + jpiglo
786            IF ( imodi > jpiglo ) imodi = imodi - jpiglo
787            !
788            DO jj = 0, imaxjfp
789               imodj = surfdataqc%mj(jobs,kvar) - int(imaxjfp/2) + jj - 1
790               !If model values are out of the domain to the north/south then
791               !set them to be the edge of the domain
792               IF ( imodj < 1      ) imodj = 1
793               IF ( imodj > jpjglo ) imodj = jpjglo
794               !
795               IF ( k2dint > 4 ) THEN
796                  igrdip1(ji+1,jj+1,iobs) = imodi
797                  igrdjp1(ji+1,jj+1,iobs) = imodj
798               ENDIF
799               !
800               IF ( ji >= 1 .AND. jj >= 1 ) THEN
801                  igrdi(ji,jj,iobs) = imodi
802                  igrdj(ji,jj,iobs) = imodj
803               ENDIF
804               !
805            END DO
806         END DO
807      END DO
808
809      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
810         &                  igrdi, igrdj, glamt, zglam )
811      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
812         &                  igrdi, igrdj, gphit, zgphi )
813      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
814         &                  igrdi, igrdj, psurfmask, zmask )
815
816      ! At the end of the averaging period get interpolated means
817      IF ( ldtime_mean ) THEN
818         IF ( imeanend == 0 ) THEN
819            ALLOCATE( zsurfm(imaxifp,imaxjfp,isurf) )
820            IF (lwp) WRITE(numout,*)' Interpolating the time mean values on time step: ', kt
821            CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
822               &                  igrdi, igrdj, surfdataqc%vdmean(:,:,kvar), zsurfm )
823         ENDIF
824      ELSE
825         CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
826            &                  igrdi, igrdj, psurf, zsurf )
827      ENDIF
828
829      IF ( k2dint > 4 ) THEN
830         CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, &
831            &                  igrdip1, igrdjp1, glamf, zglamf )
832         CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, &
833            &                  igrdip1, igrdjp1, gphif, zgphif )
834      ENDIF
835     
836      IF ( ldclim ) THEN
837         CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
838            &                  igrdi, igrdj, pclim, zclim )
839      ENDIF
840
841      ! At the end of the day get interpolated means
842      IF ( idayend == 0 .AND. ldnightav ) THEN
843
844         ALLOCATE( &
845            & zsurfm(imaxifp,imaxjfp,isurf)  &
846            & )
847
848         CALL obs_int_comm_2d( imaxifp,imaxjfp, isurf, kpi, kpj, igrdi, igrdj, &
849            &               surfdataqc%vdmean(:,:,kvar), zsurfm )
850
851      ENDIF
852
853      ! Loop over observations
854      DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf
855
856         iobs = jobs - surfdataqc%nsurfup
857
858         IF ( kt /= surfdataqc%mstp(jobs) ) THEN
859
860            IF(lwp) THEN
861               WRITE(numout,*)
862               WRITE(numout,*) ' E R R O R : Observation',              &
863                  &            ' time step is not consistent with the', &
864                  &            ' model time step'
865               WRITE(numout,*) ' ========='
866               WRITE(numout,*)
867               WRITE(numout,*) ' Record  = ', jobs,                &
868                  &            ' kt      = ', kt,                  &
869                  &            ' mstp    = ', surfdataqc%mstp(jobs), &
870                  &            ' ntyp    = ', surfdataqc%ntyp(jobs)
871            ENDIF
872            CALL ctl_stop( 'obs_surf_opt', 'Inconsistent time' )
873
874         ENDIF
875
876         zlam = surfdataqc%rlam(jobs)
877         zphi = surfdataqc%rphi(jobs)
878
879         IF ( ( ldnightav .AND. idayend == 0 ) .OR. (ldtime_mean .AND. imeanend == 0) ) THEN
880            ! Night-time or N=kmeanstp timestep averaged data
881            zsurftmp(:,:,iobs) = zsurfm(:,:,iobs)
882         ELSE
883            zsurftmp(:,:,iobs) = zsurf(:,:,iobs)
884         ENDIF
885
886         IF ( ( .NOT. ldtime_mean ) .OR. ( ldtime_mean .AND. imeanend == 0) ) THEN
887
888            IF ( k2dint <= 4 ) THEN
889
890               ! Get weights to interpolate the model value to the observation point
891               CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         &
892                  &                   zglam(:,:,iobs), zgphi(:,:,iobs), &
893                  &                   zmask(:,:,iobs), zweig, zobsmask )
894
895               ! Interpolate the model value to the observation point
896               CALL obs_int_h2d( 1, 1, zweig, zsurftmp(:,:,iobs), zext )
897
898               IF ( ldclim ) THEN
899                  CALL obs_int_h2d( 1, 1, zweig, zclim(:,:,iobs), zclm )
900               ENDIF
901
902            ELSE
903
904               ! Get weights to average the model field to the observation footprint
905               CALL obs_avg_h2d_init( 1, 1, imaxifp, imaxjfp, k2dint, zlam,  zphi, &
906                  &                   zglam(:,:,iobs), zgphi(:,:,iobs), &
907                  &                   zglamf(:,:,iobs), zgphif(:,:,iobs), &
908                  &                   zmask(:,:,iobs), plamscl, pphiscl, &
909                  &                   lindegrees, zweig, zobsmask )
910
911               ! Average the model field to the observation footprint
912               CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, &
913                  &              zweig, zsurftmp(:,:,iobs),  zext )
914
915               IF ( ldclim ) THEN
916                  CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, &
917                     &              zweig, zclim(:,:,iobs),  zclm )
918               ENDIF
919
920            ENDIF
921
922            IF ( TRIM(surfdataqc%cvars(kvar)) == cobsname_sla .AND. PRESENT(kssh) .AND. PRESENT(kmdt) ) THEN
923               ! ... Remove the MDT from the SSH at the observation point to get the SLA
924               surfdataqc%radd(jobs,kssh,kvar) = zext(1)
925               surfdataqc%rmod(jobs,kvar) = surfdataqc%radd(jobs,kssh,kvar) - surfdataqc%rext(jobs,kmdt)
926#if defined key_si3 || defined key_cice
927            ELSE IF ( TRIM(surfdataqc%cvars(kvar)) == cobsname_fbd .AND. PRESENT(kfbd) .AND. PRESENT(ksnow) ) THEN
928               surfdataqc%rmod(jobs,kvar) = zext(1)
929               ! Convert radar freeboard to true freeboard
930               ! (add 1/4 snow depth; 1/4 based on ratio of speed of light in vacuum
931               !  compared to snow (3.0e8 vs 2.4e8 m/s))
932               surfdataqc%radd(jobs,kfbd,kvar) = surfdataqc%robs(jobs,kvar)
933               surfdataqc%robs(jobs,kvar) = surfdataqc%radd(jobs,kfbd,kvar) + 0.25*surfdataqc%rext(jobs,ksnow)
934               ! If the corrected freeboard observation is outside -0.3 to 3.0 m (CPOM) then set the QC flag to bad
935               IF ((surfdataqc%robs(jobs,kvar) < -0.3) .OR. (surfdataqc%robs(jobs,kvar) > 3.0)) THEN
936                  surfdataqc%nqc(jobs) = 4
937               ENDIF           
938               ! Convert corrected freeboard to ice thickness following Tilling et al. (2016)
939               surfdataqc%robs(jobs,kvar) = (surfdataqc%robs(jobs,kvar)*rhow + surfdataqc%rext(jobs,ksnow)*rhos)/ &
940                  &                         (rhow - rhoi)
941#endif
942            ELSE
943               surfdataqc%rmod(jobs,kvar) = zext(1)
944            ENDIF
945
946            IF ( ldclim ) THEN
947               surfdataqc%radd(jobs,kclim,kvar) = zclm(1)
948            ENDIF
949
950            IF ( zext(1) == obfillflt ) THEN
951               ! If the observation value is a fill value, set QC flag to bad
952               surfdataqc%nqc(jobs) = 4
953            ENDIF
954
955         ENDIF
956
957      END DO
958
959      ! Deallocate the data for interpolation
960      DEALLOCATE( &
961         & zweig, &
962         & igrdi, &
963         & igrdj, &
964         & zglam, &
965         & zgphi, &
966         & zmask, &
967         & zsurf, &
968         & zsurftmp &
969         & )
970
971      IF ( k2dint > 4 ) THEN
972         DEALLOCATE( &     
973            & zglamf, &
974            & zgphif, &
975            & igrdip1,&
976            & igrdjp1 &
977            & )
978      ENDIF
979           
980      IF ( ldclim ) THEN
981         DEALLOCATE( zclim )
982      ENDIF
983
984      ! At the end of the day also deallocate time mean array
985      IF ( ( idayend == 0 .AND. ldnightav ) .OR. ( imeanend == 0 .AND. ldtime_mean ) ) THEN
986         DEALLOCATE( &
987            & zsurfm  &
988            & )
989      ENDIF
990      !
991      IF ( kvar == surfdataqc%nvar ) THEN
992         surfdataqc%nsurfup = surfdataqc%nsurfup + isurf
993      ENDIF
994      !
995   END SUBROUTINE obs_surf_opt
996
997   !!======================================================================
998END MODULE obs_oper
Note: See TracBrowser for help on using the repository browser.