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

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

First wrking version which produces correct seeming results.
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),      &
416                  &      iv_indic(inum_obs) ) 
417               IF ( prodatqc%lclim ) ALLOCATE( interp_corner_clim(2,2,inum_obs) )                 
418               DO iin=1,2 
419                  DO ijn=1,2 
420                   
421                     IF ( k1dint == 1 ) THEN
422                        CALL obs_int_z1d_spl( kpk, & 
423                           &    zint(iin,ijn,:,iobs),& 
424                           &    zobs2k, zgdept(iin,ijn,:,iobs), & 
425                           &    zmask(iin,ijn,:,iobs)) 
426
427                        IF ( prodatqc%lclim ) THEN
428                           CALL obs_int_z1d_spl( kpk, & 
429                              &    zclim(iin,ijn,:,iobs),& 
430                              &    zclm2k, zgdept(iin,ijn,:,iobs), & 
431                              &    zmask(iin,ijn,:,iobs)) 
432                        ENDIF
433 
434                     ENDIF
435       
436                     CALL obs_level_search(kpk, & 
437                         &        zgdept(iin,ijn,:,iobs),& 
438                         &        inum_obs, prodatqc%var(kvar)%vdep(ista:iend), & 
439                         &        iv_indic) 
440
441                     CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs,     & 
442                         &          prodatqc%var(kvar)%vdep(ista:iend),     & 
443                         &          zint(iin,ijn,:,iobs),            & 
444                         &          zobs2k,interp_corner(iin,ijn,:), & 
445                         &          zgdept(iin,ijn,:,iobs),         & 
446                         &          zmask(iin,ijn,:,iobs) )     
447
448                     IF ( prodatqc%lclim ) THEN
449                        CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs,     & 
450                            &          prodatqc%var(kvar)%vdep(ista:iend),     & 
451                            &          zclim(iin,ijn,:,iobs),            & 
452                            &          zclm2k,interp_corner_clim(iin,ijn,:), & 
453                            &          zgdept(iin,ijn,:,iobs),         & 
454                            &          zmask(iin,ijn,:,iobs) )   
455                     ENDIF   
456         
457                  ENDDO 
458               ENDDO 
459             
460            ENDIF 
461
462            !-------------------------------------------------------------
463            ! Compute the horizontal interpolation for every profile level
464            !-------------------------------------------------------------
465             
466            DO ikn=1,inum_obs 
467               iend=ista+ikn-1
468                 
469               zweig(:,:,1) = 0._wp 
470   
471               ! This code forces the horizontal weights to be 
472               ! zero IF the observation is below the bottom of the 
473               ! corners of the interpolation nodes, Or if it is in 
474               ! the mask. This is important for observations near 
475               ! steep bathymetry
476               DO iin=1,2 
477                  DO ijn=1,2 
478     
479                     depth_loop: DO ik=kpk,2,-1 
480                        IF(zmask(iin,ijn,ik-1,iobs ) > 0.9 )THEN   
481                           
482                           zweig(iin,ijn,1) = & 
483                              & zweig1(iin,ijn,1) * & 
484                              & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) & 
485                              &  - prodatqc%var(kvar)%vdep(iend)),0._wp) 
486                           
487                           EXIT depth_loop
488
489                        ENDIF
490
491                     ENDDO depth_loop
492     
493                  ENDDO 
494               ENDDO 
495   
496               CALL obs_int_h2d( 1, 1, zweig, interp_corner(:,:,ikn), & 
497                  &              prodatqc%var(kvar)%vmod(iend:iend) ) 
498
499               IF ( prodatqc%lclim ) THEN
500                  CALL obs_int_h2d( 1, 1, zweig, interp_corner_clim(:,:,ikn), & 
501                     &              prodatqc%var(kvar)%vclm(iend:iend) )
502               ENDIF
503
504                  ! Set QC flag for any observations found below the bottom
505                  ! needed as the check here is more strict than that in obs_prep
506               IF (sum(zweig) == 0.0_wp) prodatqc%var(kvar)%nvqc(iend:iend)=4
507 
508            ENDDO 
509 
510            DEALLOCATE(interp_corner,iv_indic) 
511            IF ( prodatqc%lclim ) DEALLOCATE( interp_corner_clim )         
512             
513         ENDIF
514
515      ENDDO
516
517      ! Deallocate the data for interpolation
518      DEALLOCATE(  &
519         & igrdi,  &
520         & igrdj,  &
521         & zglam,  &
522         & zgphi,  &
523         & zmask,  &
524         & zint,   &
525         & zgdept, &
526         & zgdepw  &
527         & )
528
529      IF ( prodatqc%lclim ) DEALLOCATE( zclim )
530     
531      ! At the end of the day also get interpolated means
532      IF ( ld_dailyav .AND. idayend == 0 ) THEN
533         DEALLOCATE( zinm )
534      ENDIF
535
536      IF ( kvar == prodatqc%nvar ) THEN
537         prodatqc%nprofup = prodatqc%nprofup + ipro 
538      ENDIF
539
540   END SUBROUTINE obs_prof_opt
541
542   SUBROUTINE obs_surf_opt( surfdataqc, kt, kpi, kpj,            &
543      &                     kit000, kdaystp, psurf, pclim, psurfmask,   &
544      &                     k2dint, ldnightav, plamscl, pphiscl, &
545      &                     lindegrees )
546
547      !!-----------------------------------------------------------------------
548      !!
549      !!                     ***  ROUTINE obs_surf_opt  ***
550      !!
551      !! ** Purpose : Compute the model counterpart of surface
552      !!              data by interpolating from the model grid to the
553      !!              observation point.
554      !!
555      !! ** Method  : Linearly interpolate to each observation point using
556      !!              the model values at the corners of the surrounding grid box.
557      !!
558      !!    The new model value is first computed at the obs (lon, lat) point.
559      !!
560      !!    Several horizontal interpolation schemes are available:
561      !!        - distance-weighted (great circle) (k2dint = 0)
562      !!        - distance-weighted (small angle)  (k2dint = 1)
563      !!        - bilinear (geographical grid)     (k2dint = 2)
564      !!        - bilinear (quadrilateral grid)    (k2dint = 3)
565      !!        - polynomial (quadrilateral grid)  (k2dint = 4)
566      !!
567      !!    Two horizontal averaging schemes are also available:
568      !!        - weighted radial footprint        (k2dint = 5)
569      !!        - weighted rectangular footprint   (k2dint = 6)
570      !!
571      !!
572      !! ** Action  :
573      !!
574      !! History :
575      !!      ! 07-03 (A. Weaver)
576      !!      ! 15-02 (M. Martin) Combined routine for surface types
577      !!      ! 17-03 (M. Martin) Added horizontal averaging options
578      !!-----------------------------------------------------------------------
579
580      !! * Modules used
581      USE obs_surf_def  ! Definition of storage space for surface observations
582
583      IMPLICIT NONE
584
585      !! * Arguments
586      TYPE(obs_surf), INTENT(INOUT) :: &
587         & surfdataqc                  ! Subset of surface data passing QC
588      INTEGER, INTENT(IN) :: kt        ! Time step
589      INTEGER, INTENT(IN) :: kpi       ! Model grid parameters
590      INTEGER, INTENT(IN) :: kpj
591      INTEGER, INTENT(IN) :: kit000    ! Number of the first time step
592                                       !   (kit000-1 = restart time)
593      INTEGER, INTENT(IN) :: kdaystp   ! Number of time steps per day
594      INTEGER, INTENT(IN) :: k2dint    ! Horizontal interpolation type (see header)
595      REAL(wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
596         & psurf,  &                   ! Model surface field
597         & pclim,  &                   ! Climatological surface field         
598         & psurfmask                   ! Land-sea mask
599      LOGICAL, INTENT(IN) :: ldnightav ! Logical for averaging night-time data
600      REAL(KIND=wp), INTENT(IN) :: &
601         & plamscl, &                  ! Diameter in metres of obs footprint in E/W, N/S directions
602         & pphiscl                     ! This is the full width (rather than half-width)
603      LOGICAL, INTENT(IN) :: &
604         & lindegrees                  ! T=> plamscl and pphiscl are specified in degrees, F=> in metres
605
606      !! * Local declarations
607      INTEGER :: ji
608      INTEGER :: jj
609      INTEGER :: jobs
610      INTEGER :: inrc
611      INTEGER :: isurf
612      INTEGER :: iobs
613      INTEGER :: imaxifp, imaxjfp
614      INTEGER :: imodi, imodj
615      INTEGER :: idayend
616      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
617         & igrdi,   &
618         & igrdj,   &
619         & igrdip1, &
620         & igrdjp1
621      INTEGER, DIMENSION(:,:), SAVE, ALLOCATABLE :: &
622         & icount_night,      &
623         & imask_night
624      REAL(wp) :: zlam
625      REAL(wp) :: zphi
626      REAL(wp), DIMENSION(1) :: zext, zobsmask, zclm
627      REAL(wp) :: zdaystp
628      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
629         & zweig,  &
630         & zmask,  &
631         & zsurf,  &
632         & zsurfm, &
633         & zsurftmp, &
634         & zclim,  &
635         & zglam,  &
636         & zgphi,  &
637         & zglamf, &
638         & zgphif
639
640      REAL(wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: &
641         & zintmp,  &
642         & zouttmp, &
643         & zmeanday    ! to compute model sst in region of 24h daylight (pole)
644         
645      !------------------------------------------------------------------------
646      ! Local initialization
647      !------------------------------------------------------------------------
648      ! Record and data counters
649      inrc = kt - kit000 + 2
650      isurf = surfdataqc%nsstp(inrc)
651
652      ! Work out the maximum footprint size for the
653      ! interpolation/averaging in model grid-points - has to be even.
654
655      CALL obs_max_fpsize( k2dint, plamscl, pphiscl, lindegrees, psurfmask, imaxifp, imaxjfp )
656
657
658      IF ( ldnightav ) THEN
659
660      ! Initialize array for night mean
661         IF ( kt == 0 ) THEN
662            ALLOCATE ( icount_night(kpi,kpj) )
663            ALLOCATE ( imask_night(kpi,kpj) )
664            ALLOCATE ( zintmp(kpi,kpj) )
665            ALLOCATE ( zouttmp(kpi,kpj) )
666            ALLOCATE ( zmeanday(kpi,kpj) )
667            nday_qsr = -1   ! initialisation flag for nbc_dcy
668         ENDIF
669
670         ! Night-time means are calculated for night-time values over timesteps:
671         !  [1 <= kt <= kdaystp], [kdaystp+1 <= kt <= 2*kdaystp], .....
672         idayend = MOD( kt - kit000 + 1, kdaystp )
673
674         ! Initialize night-time mean for first timestep of the day
675         IF ( idayend == 1 .OR. kt == 0 ) THEN
676            DO jj = 1, jpj
677               DO ji = 1, jpi
678                  surfdataqc%vdmean(ji,jj) = 0.0
679                  zmeanday(ji,jj) = 0.0
680                  icount_night(ji,jj) = 0
681               END DO
682            END DO
683         ENDIF
684
685         zintmp(:,:) = 0.0
686         zouttmp(:,:) = sbc_dcy( zintmp(:,:), .TRUE. )
687         imask_night(:,:) = INT( zouttmp(:,:) )
688
689         DO jj = 1, jpj
690            DO ji = 1, jpi
691               ! Increment the temperature field for computing night mean and counter
692               surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj)  &
693                      &                    + psurf(ji,jj) * REAL( imask_night(ji,jj) )
694               zmeanday(ji,jj)          = zmeanday(ji,jj) + psurf(ji,jj)
695               icount_night(ji,jj)      = icount_night(ji,jj) + imask_night(ji,jj)
696            END DO
697         END DO
698
699         ! Compute the night-time mean at the end of the day
700         zdaystp = 1.0 / REAL( kdaystp )
701         IF ( idayend == 0 ) THEN
702            IF (lwp) WRITE(numout,*) 'Calculating surfdataqc%vdmean on time-step: ',kt
703            DO jj = 1, jpj
704               DO ji = 1, jpi
705                  ! Test if "no night" point
706                  IF ( icount_night(ji,jj) > 0 ) THEN
707                     surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj) &
708                       &                        / REAL( icount_night(ji,jj) )
709                  ELSE
710                     !At locations where there is no night (e.g. poles),
711                     ! calculate daily mean instead of night-time mean.
712                     surfdataqc%vdmean(ji,jj) = zmeanday(ji,jj) * zdaystp
713                  ENDIF
714               END DO
715            END DO
716         ENDIF
717
718      ENDIF
719
720      ! Get the data for interpolation
721
722      ALLOCATE( &
723         & zweig(imaxifp,imaxjfp,1),      &
724         & igrdi(imaxifp,imaxjfp,isurf), &
725         & igrdj(imaxifp,imaxjfp,isurf), &
726         & zglam(imaxifp,imaxjfp,isurf), &
727         & zgphi(imaxifp,imaxjfp,isurf), &
728         & zmask(imaxifp,imaxjfp,isurf), &
729         & zsurf(imaxifp,imaxjfp,isurf), &
730         & zsurftmp(imaxifp,imaxjfp,isurf),  &
731         & zglamf(imaxifp+1,imaxjfp+1,isurf), &
732         & zgphif(imaxifp+1,imaxjfp+1,isurf), &
733         & igrdip1(imaxifp+1,imaxjfp+1,isurf), &
734         & igrdjp1(imaxifp+1,imaxjfp+1,isurf) &
735         & )
736
737      IF ( surfdataqc%lclim ) ALLOCATE( zclim(imaxifp,imaxjfp,isurf) )
738
739      DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf
740         iobs = jobs - surfdataqc%nsurfup
741         DO ji = 0, imaxifp
742            imodi = surfdataqc%mi(jobs) - int(imaxifp/2) + ji - 1
743           
744            !Deal with wrap around in longitude
745            IF ( imodi < 1      ) imodi = imodi + jpiglo
746            IF ( imodi > jpiglo ) imodi = imodi - jpiglo
747           
748            DO jj = 0, imaxjfp
749               imodj = surfdataqc%mj(jobs) - int(imaxjfp/2) + jj - 1
750               !If model values are out of the domain to the north/south then
751               !set them to be the edge of the domain
752               IF ( imodj < 1      ) imodj = 1
753               IF ( imodj > jpjglo ) imodj = jpjglo
754
755               igrdip1(ji+1,jj+1,iobs) = imodi
756               igrdjp1(ji+1,jj+1,iobs) = imodj
757               
758               IF ( ji >= 1 .AND. jj >= 1 ) THEN
759                  igrdi(ji,jj,iobs) = imodi
760                  igrdj(ji,jj,iobs) = imodj
761               ENDIF
762               
763            END DO
764         END DO
765      END DO
766
767      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
768         &                  igrdi, igrdj, glamt, zglam )
769      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
770         &                  igrdi, igrdj, gphit, zgphi )
771      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
772         &                  igrdi, igrdj, psurfmask, zmask )
773      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
774         &                  igrdi, igrdj, psurf, zsurf )
775      CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, &
776         &                  igrdip1, igrdjp1, glamf, zglamf )
777      CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, &
778         &                  igrdip1, igrdjp1, gphif, zgphif )
779
780      IF ( surfdataqc%lclim ) THEN
781         CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
782            &                  igrdi, igrdj, pclim, zclim )
783      ENDIF
784
785      ! At the end of the day get interpolated means
786      IF ( idayend == 0 .AND. ldnightav ) THEN
787
788         ALLOCATE( &
789            & zsurfm(imaxifp,imaxjfp,isurf)  &
790            & )
791
792         CALL obs_int_comm_2d( imaxifp,imaxjfp, isurf, kpi, kpj, igrdi, igrdj, &
793            &               surfdataqc%vdmean(:,:), zsurfm )
794
795      ENDIF
796
797      ! Loop over observations
798      DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf
799
800         iobs = jobs - surfdataqc%nsurfup
801
802         IF ( kt /= surfdataqc%mstp(jobs) ) THEN
803
804            IF(lwp) THEN
805               WRITE(numout,*)
806               WRITE(numout,*) ' E R R O R : Observation',              &
807                  &            ' time step is not consistent with the', &
808                  &            ' model time step'
809               WRITE(numout,*) ' ========='
810               WRITE(numout,*)
811               WRITE(numout,*) ' Record  = ', jobs,                &
812                  &            ' kt      = ', kt,                  &
813                  &            ' mstp    = ', surfdataqc%mstp(jobs), &
814                  &            ' ntyp    = ', surfdataqc%ntyp(jobs)
815            ENDIF
816            CALL ctl_stop( 'obs_surf_opt', 'Inconsistent time' )
817
818         ENDIF
819
820         zlam = surfdataqc%rlam(jobs)
821         zphi = surfdataqc%rphi(jobs)
822
823         IF ( ldnightav .AND. idayend == 0 ) THEN
824            ! Night-time averaged data
825            zsurftmp(:,:,iobs) = zsurfm(:,:,iobs)
826         ELSE
827            zsurftmp(:,:,iobs) = zsurf(:,:,iobs)
828         ENDIF
829
830         IF ( k2dint <= 4 ) THEN
831
832            ! Get weights to interpolate the model value to the observation point
833            CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         &
834               &                   zglam(:,:,iobs), zgphi(:,:,iobs), &
835               &                   zmask(:,:,iobs), zweig, zobsmask )
836
837            ! Interpolate the model value to the observation point
838            CALL obs_int_h2d( 1, 1, zweig, zsurftmp(:,:,iobs), zext )
839
840            IF ( surfdataqc%lclim ) THEN 
841               CALL obs_int_h2d( 1, 1, zweig, zclim(:,:,iobs), zclm )
842            ENDIF
843
844
845         ELSE
846
847            ! Get weights to average the model SLA to the observation footprint
848            CALL obs_avg_h2d_init( 1, 1, imaxifp, imaxjfp, k2dint, zlam,  zphi, &
849               &                   zglam(:,:,iobs), zgphi(:,:,iobs), &
850               &                   zglamf(:,:,iobs), zgphif(:,:,iobs), &
851               &                   zmask(:,:,iobs), plamscl, pphiscl, &
852               &                   lindegrees, zweig, zobsmask )
853
854            ! Average the model SST to the observation footprint
855            CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, &
856               &              zweig, zsurftmp(:,:,iobs),  zext )
857
858            IF ( surfdataqc%lclim ) THEN 
859               CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, &
860                  &              zweig, zsurftmp(:,:,iobs),  zclm )
861            ENDIF
862
863         ENDIF
864
865         IF ( TRIM(surfdataqc%cvars(1)) == 'SLA' .AND. surfdataqc%nextra == 2 ) THEN
866            ! ... Remove the MDT from the SSH at the observation point to get the SLA
867            surfdataqc%rext(jobs,1) = zext(1)
868            surfdataqc%rmod(jobs,1) = surfdataqc%rext(jobs,1) - surfdataqc%rext(jobs,2)
869         ELSE
870            surfdataqc%rmod(jobs,1) = zext(1)
871         ENDIF
872         
873         IF ( surfdataqc%lclim ) surfdataqc%rclm(jobs,1) = zclm(1)
874         
875         IF ( zext(1) == obfillflt ) THEN
876            ! If the observation value is a fill value, set QC flag to bad
877            surfdataqc%nqc(jobs) = 4
878         ENDIF
879
880      END DO
881
882      ! Deallocate the data for interpolation
883      DEALLOCATE( &
884         & zweig, &
885         & igrdi, &
886         & igrdj, &
887         & zglam, &
888         & zgphi, &
889         & zmask, &
890         & zsurf, &
891         & zsurftmp, &
892         & zglamf, &
893         & zgphif, &
894         & igrdip1,&
895         & igrdjp1 &
896         & )
897
898      IF ( surfdataqc%lclim ) DEALLOCATE( zclim )
899
900      ! At the end of the day also deallocate night-time mean array
901      IF ( idayend == 0 .AND. ldnightav ) THEN
902         DEALLOCATE( &
903            & zsurfm  &
904            & )
905      ENDIF
906
907      surfdataqc%nsurfup = surfdataqc%nsurfup + isurf
908
909   END SUBROUTINE obs_surf_opt
910
911END MODULE obs_oper
Note: See TracBrowser for help on using the repository browser.