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

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

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