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/r8395_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

source: NEMO/branches/UKMO/r8395_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90 @ 11361

Last change on this file since 11361 was 11361, checked in by jcastill, 5 years ago

First version of the new observation branch - it compiles, but has not been tested

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