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 branches/UKMO/dev_r5518_obs_oper_update_bgc3d/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

source: branches/UKMO/dev_r5518_obs_oper_update_bgc3d/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90 @ 9216

Last change on this file since 9216 was 9216, checked in by dford, 6 years ago

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