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

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

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