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

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

Commit changes which are in one of Jennie's orca12 branches.
File size: 34.8 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
775      IF ( k2dint > 4 ) THEN         
776         CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, &
777            &                  igrdip1, igrdjp1, glamf, zglamf )
778         CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, &
779            &                  igrdip1, igrdjp1, gphif, zgphif )
780      ENDIF
781     
782      IF ( surfdataqc%lclim ) THEN
783         CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
784            &                  igrdi, igrdj, pclim, zclim )
785      ENDIF
786
787      ! At the end of the day get interpolated means
788      IF ( idayend == 0 .AND. ldnightav ) THEN
789
790         ALLOCATE( &
791            & zsurfm(imaxifp,imaxjfp,isurf)  &
792            & )
793
794         CALL obs_int_comm_2d( imaxifp,imaxjfp, isurf, kpi, kpj, igrdi, igrdj, &
795            &               surfdataqc%vdmean(:,:), zsurfm )
796
797      ENDIF
798
799      ! Loop over observations
800      DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf
801
802         iobs = jobs - surfdataqc%nsurfup
803
804         IF ( kt /= surfdataqc%mstp(jobs) ) THEN
805
806            IF(lwp) THEN
807               WRITE(numout,*)
808               WRITE(numout,*) ' E R R O R : Observation',              &
809                  &            ' time step is not consistent with the', &
810                  &            ' model time step'
811               WRITE(numout,*) ' ========='
812               WRITE(numout,*)
813               WRITE(numout,*) ' Record  = ', jobs,                &
814                  &            ' kt      = ', kt,                  &
815                  &            ' mstp    = ', surfdataqc%mstp(jobs), &
816                  &            ' ntyp    = ', surfdataqc%ntyp(jobs)
817            ENDIF
818            CALL ctl_stop( 'obs_surf_opt', 'Inconsistent time' )
819
820         ENDIF
821
822         zlam = surfdataqc%rlam(jobs)
823         zphi = surfdataqc%rphi(jobs)
824
825         IF ( ldnightav .AND. idayend == 0 ) THEN
826            ! Night-time averaged data
827            zsurftmp(:,:,iobs) = zsurfm(:,:,iobs)
828         ELSE
829            zsurftmp(:,:,iobs) = zsurf(:,:,iobs)
830         ENDIF
831
832         IF ( k2dint <= 4 ) THEN
833
834            ! Get weights to interpolate the model value to the observation point
835            CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         &
836               &                   zglam(:,:,iobs), zgphi(:,:,iobs), &
837               &                   zmask(:,:,iobs), zweig, zobsmask )
838
839            ! Interpolate the model value to the observation point
840            CALL obs_int_h2d( 1, 1, zweig, zsurftmp(:,:,iobs), zext )
841
842            IF ( surfdataqc%lclim ) THEN 
843               CALL obs_int_h2d( 1, 1, zweig, zclim(:,:,iobs), zclm )
844            ENDIF
845
846
847         ELSE
848
849            ! Get weights to average the model SLA to the observation footprint
850            CALL obs_avg_h2d_init( 1, 1, imaxifp, imaxjfp, k2dint, zlam,  zphi, &
851               &                   zglam(:,:,iobs), zgphi(:,:,iobs), &
852               &                   zglamf(:,:,iobs), zgphif(:,:,iobs), &
853               &                   zmask(:,:,iobs), plamscl, pphiscl, &
854               &                   lindegrees, zweig, zobsmask )
855
856            ! Average the model SST to the observation footprint
857            CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, &
858               &              zweig, zsurftmp(:,:,iobs),  zext )
859
860            IF ( surfdataqc%lclim ) THEN 
861               CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, &
862                  &              zweig, zclim(:,:,iobs),  zclm )
863            ENDIF
864
865         ENDIF
866
867         IF ( TRIM(surfdataqc%cvars(1)) == 'SLA' .AND. surfdataqc%nextra == 2 ) THEN
868            ! ... Remove the MDT from the SSH at the observation point to get the SLA
869            surfdataqc%rext(jobs,1) = zext(1)
870            surfdataqc%rmod(jobs,1) = surfdataqc%rext(jobs,1) - surfdataqc%rext(jobs,2)
871         ELSE
872            surfdataqc%rmod(jobs,1) = zext(1)
873         ENDIF
874         
875         IF ( surfdataqc%lclim ) surfdataqc%rclm(jobs,1) = zclm(1)
876         
877         IF ( zext(1) == obfillflt ) THEN
878            ! If the observation value is a fill value, set QC flag to bad
879            surfdataqc%nqc(jobs) = 4
880         ENDIF
881
882      END DO
883
884      ! Deallocate the data for interpolation
885      DEALLOCATE( &
886         & zweig, &
887         & igrdi, &
888         & igrdj, &
889         & zglam, &
890         & zgphi, &
891         & zmask, &
892         & zsurf, &
893         & zsurftmp, &
894         & zglamf, &
895         & zgphif, &
896         & igrdip1,&
897         & igrdjp1 &
898         & )
899
900      IF ( surfdataqc%lclim ) DEALLOCATE( zclim )
901
902      ! At the end of the day also deallocate night-time mean array
903      IF ( idayend == 0 .AND. ldnightav ) THEN
904         DEALLOCATE( &
905            & zsurfm  &
906            & )
907      ENDIF
908
909      surfdataqc%nsurfup = surfdataqc%nsurfup + isurf
910
911   END SUBROUTINE obs_surf_opt
912
913END MODULE obs_oper
Note: See TracBrowser for help on using the repository browser.