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

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

Change to loop over obs_prof_opt one variable at a time.
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   USE dom_oce                  ! Ocean space and time domain variables 
40
41   IMPLICIT NONE
42
43   !! * Routine accessibility
44   PRIVATE
45
46   PUBLIC obs_prof_opt, &  ! Compute the model counterpart of profile obs
47      &   obs_surf_opt     ! Compute the model counterpart of surface obs
48
49   INTEGER, PARAMETER, PUBLIC :: &
50      & imaxavtypes = 20   ! Max number of daily avgd obs types
51
52   !!----------------------------------------------------------------------
53   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
54   !! $Id$
55   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
56   !!----------------------------------------------------------------------
57
58   !! * Substitutions
59#  include "domzgr_substitute.h90" 
60CONTAINS
61
62
63   SUBROUTINE obs_prof_opt( prodatqc, kt, kpi, kpj, kpk, &
64      &                     kit000, kdaystp, kvar,       &
65      &                     pvar, pgdept, pgdepw,        &
66      &                     pmask,                       & 
67      &                     plam, pphi,                  &
68      &                     k1dint, k2dint, kdailyavtypes )
69
70      !!-----------------------------------------------------------------------
71      !!
72      !!                     ***  ROUTINE obs_pro_opt  ***
73      !!
74      !! ** Purpose : Compute the model counterpart of profiles
75      !!              data by interpolating from the model grid to the
76      !!              observation point.
77      !!
78      !! ** Method  : Linearly interpolate to each observation point using
79      !!              the model values at the corners of the surrounding grid box.
80      !!
81      !!    First, a vertical profile of horizontally interpolated model
82      !!    now values is computed at the obs (lon, lat) point.
83      !!    Several horizontal interpolation schemes are available:
84      !!        - distance-weighted (great circle) (k2dint = 0)
85      !!        - distance-weighted (small angle)  (k2dint = 1)
86      !!        - bilinear (geographical grid)     (k2dint = 2)
87      !!        - bilinear (quadrilateral grid)    (k2dint = 3)
88      !!        - polynomial (quadrilateral grid)  (k2dint = 4)
89      !!
90      !!    Next, the vertical profile is interpolated to the
91      !!    data depth points. Two vertical interpolation schemes are
92      !!    available:
93      !!        - linear       (k1dint = 0)
94      !!        - Cubic spline (k1dint = 1)
95      !!
96      !!    For the cubic spline the 2nd derivative of the interpolating
97      !!    polynomial is computed before entering the vertical interpolation
98      !!    routine.
99      !!
100      !!    If the logical is switched on, the model equivalent is
101      !!    a daily mean model temperature field. So, we first compute
102      !!    the mean, then interpolate only at the end of the day.
103      !!
104      !!    Note: in situ temperature observations must be converted
105      !!    to potential temperature (the model variable) prior to
106      !!    assimilation.
107      !!
108      !! ** Action  :
109      !!
110      !! History :
111      !!      ! 97-11 (A. Weaver, S. Ricci, N. Daget)
112      !!      ! 06-03 (G. Smith) NEMOVAR migration
113      !!      ! 06-10 (A. Weaver) Cleanup
114      !!      ! 07-01 (K. Mogensen) Merge of temperature and salinity
115      !!      ! 07-03 (K. Mogensen) General handling of profiles
116      !!      ! 15-02 (M. Martin) Combined routine for all profile types
117      !!      ! 17-02 (M. Martin) Include generalised vertical coordinate changes
118      !!-----------------------------------------------------------------------
119
120      !! * Modules used
121      USE obs_profiles_def ! Definition of storage space for profile obs.
122
123      IMPLICIT NONE
124
125      !! * Arguments
126      TYPE(obs_prof), INTENT(INOUT) :: &
127         & prodatqc                  ! Subset of profile data passing QC
128      INTEGER, INTENT(IN) :: kt      ! Time step
129      INTEGER, INTENT(IN) :: kpi     ! Model grid parameters
130      INTEGER, INTENT(IN) :: kpj
131      INTEGER, INTENT(IN) :: kpk
132      INTEGER, INTENT(IN) :: kit000  ! Number of the first time step
133                                     !   (kit000-1 = restart time)
134      INTEGER, INTENT(IN) :: k1dint  ! Vertical interpolation type (see header)
135      INTEGER, INTENT(IN) :: k2dint  ! Horizontal interpolation type (see header)
136      INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day
137      INTEGER, INTENT(IN) :: kvar    ! Number of variable in prodatqc
138      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &
139         & pvar,   &                 ! Model field for variable
140         & pmask                     ! Land-sea mask for variable
141      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
142         & plam,   &                 ! Model longitudes for variable
143         & pphi                      ! Model latitudes for variable
144      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 
145         & pgdept, &                 ! Model array of depth T levels
146         & pgdepw                    ! Model array of depth W levels
147      INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &
148         & kdailyavtypes             ! Types for daily averages
149
150      !! * Local declarations
151      INTEGER ::   ji
152      INTEGER ::   jj
153      INTEGER ::   jk
154      INTEGER ::   jobs
155      INTEGER ::   inrc
156      INTEGER ::   ipro
157      INTEGER ::   idayend
158      INTEGER ::   ista
159      INTEGER ::   iend
160      INTEGER ::   iobs
161      INTEGER ::   iin, ijn, ikn, ik   ! looping indices over interpolation nodes
162      INTEGER ::   inum_obs
163      INTEGER, DIMENSION(imaxavtypes) :: &
164         & idailyavtypes
165      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
166         & igrdi, &
167         & igrdj
168      INTEGER, ALLOCATABLE, DIMENSION(:) :: iv_indic
169
170      REAL(KIND=wp) :: zlam
171      REAL(KIND=wp) :: zphi
172      REAL(KIND=wp) :: zdaystp
173      REAL(KIND=wp), DIMENSION(kpk) :: &
174         & zobsk,    &
175         & zobs2k
176      REAL(KIND=wp), DIMENSION(2,2,1) :: &
177         & zweig1, &
178         & zweig
179      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
180         & zmask,  &
181         & zint,   &
182         & zinm,   &
183         & zgdept, & 
184         & zgdepw
185      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
186         & zglam,  &
187         & zgphi
188      REAL(KIND=wp), DIMENSION(1) :: zmsk
189      REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner
190
191      LOGICAL :: ld_dailyav
192
193      !------------------------------------------------------------------------
194      ! Local initialization
195      !------------------------------------------------------------------------
196      ! Record and data counters
197      inrc = kt - kit000 + 2
198      ipro = prodatqc%npstp(inrc)
199
200      ! Daily average types
201      ld_dailyav = .FALSE.
202      IF ( PRESENT(kdailyavtypes) ) THEN
203         idailyavtypes(:) = kdailyavtypes(:)
204         IF ( ANY (idailyavtypes(:) /= -1) ) ld_dailyav = .TRUE.
205      ELSE
206         idailyavtypes(:) = -1
207      ENDIF
208
209      ! Daily means are calculated for values over timesteps:
210      !  [1 <= kt <= kdaystp], [kdaystp+1 <= kt <= 2*kdaystp], ...
211      idayend = MOD( kt - kit000 + 1, kdaystp )
212
213      IF ( ld_dailyav ) THEN
214
215         ! Initialize daily mean for first timestep of the day
216         IF ( idayend == 1 .OR. kt == 0 ) THEN
217            DO jk = 1, jpk
218               DO jj = 1, jpj
219                  DO ji = 1, jpi
220                     prodatqc%vdmean(ji,jj,jk,kvar) = 0.0
221                  END DO
222               END DO
223            END DO
224         ENDIF
225
226         DO jk = 1, jpk
227            DO jj = 1, jpj
228               DO ji = 1, jpi
229                  ! Increment field for computing daily mean
230                  prodatqc%vdmean(ji,jj,jk,kvar) = prodatqc%vdmean(ji,jj,jk,kvar) &
231                     &                           + pvar(ji,jj,jk)
232               END DO
233            END DO
234         END DO
235
236         ! Compute the daily mean at the end of day
237         zdaystp = 1.0 / REAL( kdaystp )
238         IF ( idayend == 0 ) THEN
239            IF (lwp) WRITE(numout,*) 'Calculating prodatqc%vdmean on time-step: ',kt
240            CALL FLUSH(numout)
241            DO jk = 1, jpk
242               DO jj = 1, jpj
243                  DO ji = 1, jpi
244                     prodatqc%vdmean(ji,jj,jk,kvar) = prodatqc%vdmean(ji,jj,jk,kvar) &
245                        &                           * zdaystp
246                  END DO
247               END DO
248            END DO
249         ENDIF
250
251      ENDIF
252
253      ! Get the data for interpolation
254      ALLOCATE( &
255         & igrdi(2,2,ipro),       &
256         & igrdj(2,2,ipro),       &
257         & zglam(2,2,ipro),       &
258         & zgphi(2,2,ipro),       &
259         & zmask(2,2,kpk,ipro),   &
260         & zint(2,2,kpk,ipro),    &
261         & zgdept(2,2,kpk,ipro),  & 
262         & zgdepw(2,2,kpk,ipro)   & 
263         & )
264
265      DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
266         iobs = jobs - prodatqc%nprofup
267         igrdi(1,1,iobs) = prodatqc%mi(jobs,kvar)-1
268         igrdj(1,1,iobs) = prodatqc%mj(jobs,kvar)-1
269         igrdi(1,2,iobs) = prodatqc%mi(jobs,kvar)-1
270         igrdj(1,2,iobs) = prodatqc%mj(jobs,kvar)
271         igrdi(2,1,iobs) = prodatqc%mi(jobs,kvar)
272         igrdj(2,1,iobs) = prodatqc%mj(jobs,kvar)-1
273         igrdi(2,2,iobs) = prodatqc%mi(jobs,kvar)
274         igrdj(2,2,iobs) = prodatqc%mj(jobs,kvar)
275      END DO
276
277      ! Initialise depth arrays
278      zgdept(:,:,:,:) = 0.0
279      zgdepw(:,:,:,:) = 0.0
280
281      CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi, igrdj, plam, zglam )
282      CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi, igrdj, pphi, zgphi )
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,1) 
348                  iend = prodatqc%npvend(jobs,1) 
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.