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

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

Add option to use time mean background.

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