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

source: branches/UKMO/dev_r5518_obs_oper_update_addclim/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90 @ 11460

Last change on this file since 11460 was 11460, checked in by mattmartin, 5 years ago

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