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

source: branches/UKMO/dev_r5518_obs_oper_update_sit/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90 @ 12380

Last change on this file since 12380 was 12380, checked in by dcarneir, 2 years ago

Updating obs_oper with obs_oper version in the apps trunk

File size: 38.9 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#if defined key_cice
40   USE ice_constants, ONLY : &    ! For conversion from sea ice freeboard to thickness
41      & rhos, rhoi, rhow
42#endif
43
44   IMPLICIT NONE
45
46   !! * Routine accessibility
47   PRIVATE
48
49   PUBLIC obs_prof_opt, &  ! Compute the model counterpart of profile obs
50      &   obs_surf_opt     ! Compute the model counterpart of surface obs
51
52   INTEGER, PARAMETER, PUBLIC :: &
53      & imaxavtypes = 20   ! Max number of daily avgd obs types
54
55   !!----------------------------------------------------------------------
56   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
57   !! $Id$
58   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
59   !!----------------------------------------------------------------------
60
61   !! * Substitutions
62#  include "domzgr_substitute.h90" 
63CONTAINS
64
65
66   SUBROUTINE obs_prof_opt( prodatqc, kt, kpi, kpj, kpk, &
67      &                     kit000, kdaystp, kvar,       &
68      &                     pvar, pclim,                 &
69      &                     pgdept, pgdepw, pmask,       & 
70      &                     plam, pphi,                  &
71      &                     k1dint, k2dint, kdailyavtypes )
72
73      !!-----------------------------------------------------------------------
74      !!
75      !!                     ***  ROUTINE obs_pro_opt  ***
76      !!
77      !! ** Purpose : Compute the model counterpart of profiles
78      !!              data by interpolating from the model grid to the
79      !!              observation point.
80      !!
81      !! ** Method  : Linearly interpolate to each observation point using
82      !!              the model values at the corners of the surrounding grid box.
83      !!
84      !!    First, a vertical profile of horizontally interpolated model
85      !!    now values is computed at the obs (lon, lat) point.
86      !!    Several horizontal interpolation schemes are available:
87      !!        - distance-weighted (great circle) (k2dint = 0)
88      !!        - distance-weighted (small angle)  (k2dint = 1)
89      !!        - bilinear (geographical grid)     (k2dint = 2)
90      !!        - bilinear (quadrilateral grid)    (k2dint = 3)
91      !!        - polynomial (quadrilateral grid)  (k2dint = 4)
92      !!
93      !!    Next, the vertical profile is interpolated to the
94      !!    data depth points. Two vertical interpolation schemes are
95      !!    available:
96      !!        - linear       (k1dint = 0)
97      !!        - Cubic spline (k1dint = 1)
98      !!
99      !!    For the cubic spline the 2nd derivative of the interpolating
100      !!    polynomial is computed before entering the vertical interpolation
101      !!    routine.
102      !!
103      !!    If the logical is switched on, the model equivalent is
104      !!    a daily mean model temperature field. So, we first compute
105      !!    the mean, then interpolate only at the end of the day.
106      !!
107      !!    Note: in situ temperature observations must be converted
108      !!    to potential temperature (the model variable) prior to
109      !!    assimilation.
110      !!
111      !! ** Action  :
112      !!
113      !! History :
114      !!      ! 97-11 (A. Weaver, S. Ricci, N. Daget)
115      !!      ! 06-03 (G. Smith) NEMOVAR migration
116      !!      ! 06-10 (A. Weaver) Cleanup
117      !!      ! 07-01 (K. Mogensen) Merge of temperature and salinity
118      !!      ! 07-03 (K. Mogensen) General handling of profiles
119      !!      ! 15-02 (M. Martin) Combined routine for all profile types
120      !!      ! 17-02 (M. Martin) Include generalised vertical coordinate changes
121      !!-----------------------------------------------------------------------
122
123      !! * Modules used
124      USE obs_profiles_def ! Definition of storage space for profile obs.
125
126      IMPLICIT NONE
127
128      !! * Arguments
129      TYPE(obs_prof), INTENT(INOUT) :: &
130         & prodatqc                  ! Subset of profile data passing QC
131      INTEGER, INTENT(IN) :: kt      ! Time step
132      INTEGER, INTENT(IN) :: kpi     ! Model grid parameters
133      INTEGER, INTENT(IN) :: kpj
134      INTEGER, INTENT(IN) :: kpk
135      INTEGER, INTENT(IN) :: kit000  ! Number of the first time step
136                                     !   (kit000-1 = restart time)
137      INTEGER, INTENT(IN) :: k1dint  ! Vertical interpolation type (see header)
138      INTEGER, INTENT(IN) :: k2dint  ! Horizontal interpolation type (see header)
139      INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day
140      INTEGER, INTENT(IN) :: kvar    ! Number of variable in prodatqc
141      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &
142         & pvar,   &                 ! Model field for variable
143         & pclim,  &                 ! Climatology field for variable         
144         & pmask                     ! Land-sea mask for variable
145      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
146         & plam,   &                 ! Model longitudes for variable
147         & pphi                      ! Model latitudes for variable
148      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 
149         & pgdept, &                 ! Model array of depth T levels
150         & pgdepw                    ! Model array of depth W levels
151      INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &
152         & kdailyavtypes             ! Types for daily averages
153
154      !! * Local declarations
155      INTEGER ::   ji
156      INTEGER ::   jj
157      INTEGER ::   jk
158      INTEGER ::   jobs
159      INTEGER ::   inrc
160      INTEGER ::   ipro
161      INTEGER ::   idayend
162      INTEGER ::   ista
163      INTEGER ::   iend
164      INTEGER ::   iobs
165      INTEGER ::   iin, ijn, ikn, ik   ! looping indices over interpolation nodes
166      INTEGER ::   inum_obs
167      INTEGER, DIMENSION(imaxavtypes) :: &
168         & idailyavtypes
169      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
170         & igrdi, &
171         & igrdj
172      INTEGER, ALLOCATABLE, DIMENSION(:) :: iv_indic
173
174      REAL(KIND=wp) :: zlam
175      REAL(KIND=wp) :: zphi
176      REAL(KIND=wp) :: zdaystp
177      REAL(KIND=wp), DIMENSION(kpk) :: &
178         & zobsk,    &
179         & zobs2k,   &
180         & zclm2k         
181      REAL(KIND=wp), DIMENSION(2,2,1) :: &
182         & zweig1, &
183         & zweig
184      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
185         & zmask,  &
186         & zclim,  &         
187         & zint,   &
188         & zinm,   &
189         & zgdept, & 
190         & zgdepw
191      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
192         & zglam,  &
193         & zgphi
194      REAL(KIND=wp), DIMENSION(1) :: zmsk
195      REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner
196      REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner_clim
197     
198      LOGICAL :: ld_dailyav
199
200      !------------------------------------------------------------------------
201      ! Local initialization
202      !------------------------------------------------------------------------
203      ! Record and data counters
204      inrc = kt - kit000 + 2
205      ipro = prodatqc%npstp(inrc)
206
207      ! Daily average types
208      ld_dailyav = .FALSE.
209      IF ( PRESENT(kdailyavtypes) ) THEN
210         idailyavtypes(:) = kdailyavtypes(:)
211         IF ( ANY (idailyavtypes(:) /= -1) ) ld_dailyav = .TRUE.
212      ELSE
213         idailyavtypes(:) = -1
214      ENDIF
215
216      ! Daily means are calculated for values over timesteps:
217      !  [1 <= kt <= kdaystp], [kdaystp+1 <= kt <= 2*kdaystp], ...
218      idayend = MOD( kt - kit000 + 1, kdaystp )
219
220      IF ( ld_dailyav ) THEN
221
222         ! Initialize daily mean for first timestep of the day
223         IF ( idayend == 1 .OR. kt == 0 ) THEN
224            DO jk = 1, jpk
225               DO jj = 1, jpj
226                  DO ji = 1, jpi
227                     prodatqc%vdmean(ji,jj,jk,kvar) = 0.0
228                  END DO
229               END DO
230            END DO
231         ENDIF
232
233         DO jk = 1, jpk
234            DO jj = 1, jpj
235               DO ji = 1, jpi
236                  ! Increment field for computing daily mean
237                  prodatqc%vdmean(ji,jj,jk,kvar) = prodatqc%vdmean(ji,jj,jk,kvar) &
238                     &                           + pvar(ji,jj,jk)
239               END DO
240            END DO
241         END DO
242
243         ! Compute the daily mean at the end of day
244         zdaystp = 1.0 / REAL( kdaystp )
245         IF ( idayend == 0 ) THEN
246            IF (lwp) WRITE(numout,*) 'Calculating prodatqc%vdmean on time-step: ',kt
247            CALL FLUSH(numout)
248            DO jk = 1, jpk
249               DO jj = 1, jpj
250                  DO ji = 1, jpi
251                     prodatqc%vdmean(ji,jj,jk,kvar) = prodatqc%vdmean(ji,jj,jk,kvar) &
252                        &                           * zdaystp
253                  END DO
254               END DO
255            END DO
256         ENDIF
257
258      ENDIF
259
260      ! Get the data for interpolation
261      ALLOCATE( &
262         & igrdi(2,2,ipro),       &
263         & igrdj(2,2,ipro),       &
264         & zglam(2,2,ipro),       &
265         & zgphi(2,2,ipro),       &
266         & zmask(2,2,kpk,ipro),   &
267         & zint(2,2,kpk,ipro),    &
268         & zgdept(2,2,kpk,ipro),  & 
269         & zgdepw(2,2,kpk,ipro)   & 
270         & )
271
272      IF ( prodatqc%lclim ) ALLOCATE( zclim(2,2,kpk,ipro) )
273
274      DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
275         iobs = jobs - prodatqc%nprofup
276         igrdi(1,1,iobs) = prodatqc%mi(jobs,kvar)-1
277         igrdj(1,1,iobs) = prodatqc%mj(jobs,kvar)-1
278         igrdi(1,2,iobs) = prodatqc%mi(jobs,kvar)-1
279         igrdj(1,2,iobs) = prodatqc%mj(jobs,kvar)
280         igrdi(2,1,iobs) = prodatqc%mi(jobs,kvar)
281         igrdj(2,1,iobs) = prodatqc%mj(jobs,kvar)-1
282         igrdi(2,2,iobs) = prodatqc%mi(jobs,kvar)
283         igrdj(2,2,iobs) = prodatqc%mj(jobs,kvar)
284      END DO
285
286      ! Initialise depth arrays
287      zgdept(:,:,:,:) = 0.0
288      zgdepw(:,:,:,:) = 0.0
289
290      CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi, igrdj, plam, zglam )
291      CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi, igrdj, pphi, zgphi )
292      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pmask, zmask )
293      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pvar,   zint )
294
295      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pgdept, zgdept ) 
296      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pgdepw, zgdepw ) 
297
298      IF ( prodatqc%lclim ) THEN
299         CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pclim, zclim )           
300      ENDIF 
301     
302      ! At the end of the day also get interpolated means
303      IF ( ld_dailyav .AND. idayend == 0 ) THEN
304
305         ALLOCATE( zinm(2,2,kpk,ipro) )
306
307         CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, &
308            &                  prodatqc%vdmean(:,:,:,kvar), zinm )
309
310      ENDIF
311
312      ! Return if no observations to process
313      ! Has to be done after comm commands to ensure processors
314      ! stay in sync
315      IF ( ipro == 0 ) RETURN
316
317      DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro
318
319         iobs = jobs - prodatqc%nprofup
320
321         IF ( kt /= prodatqc%mstp(jobs) ) THEN
322
323            IF(lwp) THEN
324               WRITE(numout,*)
325               WRITE(numout,*) ' E R R O R : Observation',              &
326                  &            ' time step is not consistent with the', &
327                  &            ' model time step'
328               WRITE(numout,*) ' ========='
329               WRITE(numout,*)
330               WRITE(numout,*) ' Record  = ', jobs,                    &
331                  &            ' kt      = ', kt,                      &
332                  &            ' mstp    = ', prodatqc%mstp(jobs), &
333                  &            ' ntyp    = ', prodatqc%ntyp(jobs)
334            ENDIF
335            CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' )
336         ENDIF
337
338         zlam = prodatqc%rlam(jobs)
339         zphi = prodatqc%rphi(jobs)
340
341         ! Horizontal weights
342         ! Masked values are calculated later. 
343         IF ( prodatqc%npvend(jobs,kvar) > 0 ) THEN
344
345            CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,     &
346               &                   zglam(:,:,iobs), zgphi(:,:,iobs), &
347               &                   zmask(:,:,1,iobs), zweig1, zmsk )
348
349         ENDIF
350
351         IF ( prodatqc%npvend(jobs,kvar) > 0 ) THEN
352
353            zobsk(:) = obfillflt
354
355            IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN
356
357               IF ( idayend == 0 )  THEN
358                  ! Daily averaged data
359
360                  ! vertically interpolate all 4 corners
361                  ista = prodatqc%npvsta(jobs,kvar) 
362                  iend = prodatqc%npvend(jobs,kvar) 
363                  inum_obs = iend - ista + 1 
364                  ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs)) 
365                  IF ( prodatqc%lclim ) ALLOCATE( interp_corner_clim(2,2,inum_obs) )
366                 
367                  DO iin=1,2 
368                     DO ijn=1,2 
369
370                        IF ( k1dint == 1 ) THEN
371                           CALL obs_int_z1d_spl( kpk, & 
372                              &     zinm(iin,ijn,:,iobs), & 
373                              &     zobs2k, zgdept(iin,ijn,:,iobs), & 
374                              &     zmask(iin,ijn,:,iobs)) 
375
376                           IF ( prodatqc%lclim ) THEN
377                              CALL obs_int_z1d_spl( kpk, & 
378                                 &     zclim(iin,ijn,:,iobs), & 
379                                 &     zclm2k, zgdept(iin,ijn,:,iobs), & 
380                                 &     zmask(iin,ijn,:,iobs)) 
381                           ENDIF
382
383                        ENDIF
384       
385                        CALL obs_level_search(kpk, & 
386                           &    zgdept(iin,ijn,:,iobs), & 
387                           &    inum_obs, prodatqc%var(kvar)%vdep(ista:iend), & 
388                           &    iv_indic) 
389
390                        CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 
391                           &    prodatqc%var(kvar)%vdep(ista:iend), & 
392                           &    zinm(iin,ijn,:,iobs), & 
393                           &    zobs2k, interp_corner(iin,ijn,:), & 
394                           &    zgdept(iin,ijn,:,iobs), & 
395                           &    zmask(iin,ijn,:,iobs)) 
396
397                        IF ( prodatqc%lclim ) THEN
398                           CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, & 
399                              &    prodatqc%var(kvar)%vdep(ista:iend), & 
400                              &    zclim(iin,ijn,:,iobs), & 
401                              &    zclm2k, interp_corner_clim(iin,ijn,:), & 
402                              &    zgdept(iin,ijn,:,iobs), & 
403                              &    zmask(iin,ijn,:,iobs)) 
404                        ENDIF
405                       
406                     ENDDO 
407                  ENDDO 
408
409               ENDIF !idayend
410
411            ELSE   
412
413               ! Point data
414     
415               ! vertically interpolate all 4 corners
416               ista = prodatqc%npvsta(jobs,kvar) 
417               iend = prodatqc%npvend(jobs,kvar) 
418               inum_obs = iend - ista + 1 
419               ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs)) 
420               IF ( prodatqc%lclim ) ALLOCATE( interp_corner_clim(2,2,inum_obs) )                 
421               DO iin=1,2 
422                  DO ijn=1,2 
423                   
424                     IF ( k1dint == 1 ) THEN
425                        CALL obs_int_z1d_spl( kpk, & 
426                           &    zint(iin,ijn,:,iobs),& 
427                           &    zobs2k, zgdept(iin,ijn,:,iobs), & 
428                           &    zmask(iin,ijn,:,iobs)) 
429
430                        IF ( prodatqc%lclim ) THEN
431                           CALL obs_int_z1d_spl( kpk, & 
432                              &    zclim(iin,ijn,:,iobs),& 
433                              &    zclm2k, zgdept(iin,ijn,:,iobs), & 
434                              &    zmask(iin,ijn,:,iobs)) 
435                        ENDIF
436 
437                     ENDIF
438       
439                     CALL obs_level_search(kpk, & 
440                         &        zgdept(iin,ijn,:,iobs),& 
441                         &        inum_obs, prodatqc%var(kvar)%vdep(ista:iend), & 
442                         &        iv_indic) 
443
444                     CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs,     & 
445                         &          prodatqc%var(kvar)%vdep(ista:iend),     & 
446                         &          zint(iin,ijn,:,iobs),            & 
447                         &          zobs2k,interp_corner(iin,ijn,:), & 
448                         &          zgdept(iin,ijn,:,iobs),         & 
449                         &          zmask(iin,ijn,:,iobs) )     
450
451                     IF ( prodatqc%lclim ) THEN
452                        CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs,     & 
453                            &          prodatqc%var(kvar)%vdep(ista:iend),     & 
454                            &          zclim(iin,ijn,:,iobs),            & 
455                            &          zclm2k,interp_corner_clim(iin,ijn,:), & 
456                            &          zgdept(iin,ijn,:,iobs),         & 
457                            &          zmask(iin,ijn,:,iobs) )   
458                     ENDIF   
459         
460                  ENDDO 
461               ENDDO 
462             
463            ENDIF 
464
465            !-------------------------------------------------------------
466            ! Compute the horizontal interpolation for every profile level
467            !-------------------------------------------------------------
468             
469            DO ikn=1,inum_obs 
470               iend=ista+ikn-1
471                 
472               zweig(:,:,1) = 0._wp 
473   
474               ! This code forces the horizontal weights to be 
475               ! zero IF the observation is below the bottom of the 
476               ! corners of the interpolation nodes, Or if it is in 
477               ! the mask. This is important for observations near 
478               ! steep bathymetry
479               DO iin=1,2 
480                  DO ijn=1,2 
481     
482                     depth_loop: DO ik=kpk,2,-1 
483                        IF(zmask(iin,ijn,ik-1,iobs ) > 0.9 )THEN   
484                           
485                           zweig(iin,ijn,1) = & 
486                              & zweig1(iin,ijn,1) * & 
487                              & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) & 
488                              &  - prodatqc%var(kvar)%vdep(iend)),0._wp) 
489                           
490                           EXIT depth_loop
491
492                        ENDIF
493
494                     ENDDO depth_loop
495     
496                  ENDDO 
497               ENDDO 
498   
499               CALL obs_int_h2d( 1, 1, zweig, interp_corner(:,:,ikn), & 
500                  &              prodatqc%var(kvar)%vmod(iend:iend) ) 
501
502               IF ( prodatqc%lclim ) THEN
503                  CALL obs_int_h2d( 1, 1, zweig, interp_corner_clim(:,:,ikn), & 
504                     &              prodatqc%var(kvar)%vclm(iend:iend) )
505               ENDIF
506
507                  ! Set QC flag for any observations found below the bottom
508                  ! needed as the check here is more strict than that in obs_prep
509               IF (sum(zweig) == 0.0_wp) prodatqc%var(kvar)%nvqc(iend:iend)=4
510 
511            ENDDO 
512 
513            DEALLOCATE(interp_corner,iv_indic) 
514            IF ( prodatqc%lclim ) DEALLOCATE( interp_corner_clim )         
515             
516         ENDIF
517
518      ENDDO
519
520      ! Deallocate the data for interpolation
521      DEALLOCATE(  &
522         & igrdi,  &
523         & igrdj,  &
524         & zglam,  &
525         & zgphi,  &
526         & zmask,  &
527         & zint,   &
528         & zgdept, &
529         & zgdepw  &
530         & )
531
532      IF ( prodatqc%lclim ) DEALLOCATE( zclim )
533     
534      ! At the end of the day also get interpolated means
535      IF ( ld_dailyav .AND. idayend == 0 ) THEN
536         DEALLOCATE( zinm )
537      ENDIF
538
539      IF ( kvar == prodatqc%nvar ) THEN
540         prodatqc%nprofup = prodatqc%nprofup + ipro 
541      ENDIF
542
543   END SUBROUTINE obs_prof_opt
544
545   SUBROUTINE obs_surf_opt( surfdataqc, kt, kpi, kpj,            &
546      &                     kit000, kdaystp, psurf, pclim, psurfmask,   &
547      &                     k2dint, ldnightav, plamscl, pphiscl, &
548      &                     lindegrees, kmeanstp )
549
550      !!-----------------------------------------------------------------------
551      !!
552      !!                     ***  ROUTINE obs_surf_opt  ***
553      !!
554      !! ** Purpose : Compute the model counterpart of surface
555      !!              data by interpolating from the model grid to the
556      !!              observation point.
557      !!
558      !! ** Method  : Linearly interpolate to each observation point using
559      !!              the model values at the corners of the surrounding grid box.
560      !!
561      !!    The new model value is first computed at the obs (lon, lat) point.
562      !!
563      !!    Several horizontal interpolation schemes are available:
564      !!        - distance-weighted (great circle) (k2dint = 0)
565      !!        - distance-weighted (small angle)  (k2dint = 1)
566      !!        - bilinear (geographical grid)     (k2dint = 2)
567      !!        - bilinear (quadrilateral grid)    (k2dint = 3)
568      !!        - polynomial (quadrilateral grid)  (k2dint = 4)
569      !!
570      !!    Two horizontal averaging schemes are also available:
571      !!        - weighted radial footprint        (k2dint = 5)
572      !!        - weighted rectangular footprint   (k2dint = 6)
573      !!
574      !!
575      !! ** Action  :
576      !!
577      !! History :
578      !!      ! 07-03 (A. Weaver)
579      !!      ! 15-02 (M. Martin) Combined routine for surface types
580      !!      ! 17-03 (M. Martin) Added horizontal averaging options
581      !!-----------------------------------------------------------------------
582
583      !! * Modules used
584      USE obs_surf_def  ! Definition of storage space for surface observations
585
586      IMPLICIT NONE
587
588      !! * Arguments
589      TYPE(obs_surf), INTENT(INOUT) :: &
590         & surfdataqc                  ! Subset of surface data passing QC
591      INTEGER, INTENT(IN) :: kt        ! Time step
592      INTEGER, INTENT(IN) :: kpi       ! Model grid parameters
593      INTEGER, INTENT(IN) :: kpj
594      INTEGER, INTENT(IN) :: kit000    ! Number of the first time step
595                                       !   (kit000-1 = restart time)
596      INTEGER, INTENT(IN) :: kdaystp   ! Number of time steps per day
597      INTEGER, INTENT(IN) :: k2dint    ! Horizontal interpolation type (see header)
598      INTEGER, INTENT(IN), OPTIONAL :: &
599                             kmeanstp  ! Number of time steps for the time meaning
600                                       ! Averaging is triggered if present and greater than one                   
601      REAL(wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
602         & psurf,  &                   ! Model surface field
603         & pclim,  &                   ! Climatological surface field         
604         & psurfmask                   ! Land-sea mask
605      LOGICAL, INTENT(IN) :: ldnightav ! Logical for averaging night-time data
606      REAL(KIND=wp), INTENT(IN) :: &
607         & plamscl, &                  ! Diameter in metres of obs footprint in E/W, N/S directions
608         & pphiscl                     ! This is the full width (rather than half-width)
609      LOGICAL, INTENT(IN) :: &
610         & lindegrees                  ! T=> plamscl and pphiscl are specified in degrees, F=> in metres
611
612      !! * Local declarations
613      INTEGER :: ji
614      INTEGER :: jj
615      INTEGER :: jobs
616      INTEGER :: inrc
617      INTEGER :: isurf
618      INTEGER :: iobs
619      INTEGER :: imaxifp, imaxjfp
620      INTEGER :: imodi, imodj
621      INTEGER :: idayend
622      INTEGER :: imeanend
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) :: zmeanstp
636      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
637         & zweig,  &
638         & zmask,  &
639         & zsurf,  &
640         & zsurfm, &
641         & zsurftmp, &
642         & zclim,  &
643         & zglam,  &
644         & zgphi,  &
645         & zglamf, &
646         & zgphif
647
648      REAL(wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: &
649         & zintmp,  &
650         & zouttmp, &
651         & zmeanday    ! to compute model sst in region of 24h daylight (pole)
652
653      LOGICAL :: l_timemean
654         
655      !------------------------------------------------------------------------
656      ! Local initialization
657      !------------------------------------------------------------------------
658      ! Record and data counters
659      inrc = kt - kit000 + 2
660      isurf = surfdataqc%nsstp(inrc)
661
662      l_timemean = .FALSE.
663      IF ( PRESENT( kmeanstp ) ) THEN
664         IF ( kmeanstp > 1 ) l_timemean = .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 ( l_timemean ) THEN
674         ! Initialize time mean for first timestep
675         imeanend = MOD( kt - kit000 + 1, kmeanstp )
676         IF (lwp) WRITE(numout,*) 'Obs time mean ', kt, kit000, kmeanstp, imeanend
677
678         ! Added kt == 0 test to catch restart case
679         IF ( ( imeanend == 1 ) .OR. ( kt == 0 ) ) THEN
680            IF (lwp) WRITE(numout,*) 'Reset surfdataqc%vdmean on time-step: ',kt
681            DO jj = 1, jpj
682               DO ji = 1, jpi
683                  surfdataqc%vdmean(ji,jj) = 0.0
684               END DO
685            END DO
686         ENDIF
687
688         ! On each time-step, increment the field for computing time mean
689         IF (lwp) WRITE(numout,*)'Accumulating surfdataqc%vdmean on time-step: ',kt
690         DO jj = 1, jpj
691            DO ji = 1, jpi
692               surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj) &
693                  &                        + psurf(ji,jj)
694            END DO
695         END DO
696
697         ! Compute the time mean at the end of time period
698         IF ( imeanend == 0 ) THEN
699            zmeanstp = 1.0 / REAL( kmeanstp )
700            IF (lwp) WRITE(numout,*)'Calculating surfdataqc%vdmean time mean on time-step: ',kt,' with weight: ',zmeanstp
701            DO jj = 1, jpj
702               DO ji = 1, jpi
703                  surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj) &
704                     &                       * zmeanstp
705               END DO
706            END DO
707         ENDIF
708      ENDIF !l_timemean
709
710
711      IF ( ldnightav ) THEN
712
713      ! Initialize array for night mean
714         IF ( kt == 0 ) THEN
715            ALLOCATE ( icount_night(kpi,kpj) )
716            ALLOCATE ( imask_night(kpi,kpj) )
717            ALLOCATE ( zintmp(kpi,kpj) )
718            ALLOCATE ( zouttmp(kpi,kpj) )
719            ALLOCATE ( zmeanday(kpi,kpj) )
720            nday_qsr = -1   ! initialisation flag for nbc_dcy
721         ENDIF
722
723         ! Night-time means are calculated for night-time values over timesteps:
724         !  [1 <= kt <= kdaystp], [kdaystp+1 <= kt <= 2*kdaystp], .....
725         idayend = MOD( kt - kit000 + 1, kdaystp )
726
727         ! Initialize night-time mean for first timestep of the day
728         IF ( idayend == 1 .OR. kt == 0 ) THEN
729            DO jj = 1, jpj
730               DO ji = 1, jpi
731                  surfdataqc%vdmean(ji,jj) = 0.0
732                  zmeanday(ji,jj) = 0.0
733                  icount_night(ji,jj) = 0
734               END DO
735            END DO
736         ENDIF
737
738         zintmp(:,:) = 0.0
739         zouttmp(:,:) = sbc_dcy( zintmp(:,:), .TRUE. )
740         imask_night(:,:) = INT( zouttmp(:,:) )
741
742         DO jj = 1, jpj
743            DO ji = 1, jpi
744               ! Increment the temperature field for computing night mean and counter
745               surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj)  &
746                      &                    + psurf(ji,jj) * REAL( imask_night(ji,jj) )
747               zmeanday(ji,jj)          = zmeanday(ji,jj) + psurf(ji,jj)
748               icount_night(ji,jj)      = icount_night(ji,jj) + imask_night(ji,jj)
749            END DO
750         END DO
751
752         ! Compute the night-time mean at the end of the day
753         zdaystp = 1.0 / REAL( kdaystp )
754         IF ( idayend == 0 ) THEN
755            IF (lwp) WRITE(numout,*) 'Calculating surfdataqc%vdmean on time-step: ',kt
756            DO jj = 1, jpj
757               DO ji = 1, jpi
758                  ! Test if "no night" point
759                  IF ( icount_night(ji,jj) > 0 ) THEN
760                     surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj) &
761                       &                        / REAL( icount_night(ji,jj) )
762                  ELSE
763                     !At locations where there is no night (e.g. poles),
764                     ! calculate daily mean instead of night-time mean.
765                     surfdataqc%vdmean(ji,jj) = zmeanday(ji,jj) * zdaystp
766                  ENDIF
767               END DO
768            END DO
769         ENDIF
770
771      ENDIF
772
773      ! Get the data for interpolation
774
775      ALLOCATE( &
776         & zweig(imaxifp,imaxjfp,1),      &
777         & igrdi(imaxifp,imaxjfp,isurf), &
778         & igrdj(imaxifp,imaxjfp,isurf), &
779         & zglam(imaxifp,imaxjfp,isurf), &
780         & zgphi(imaxifp,imaxjfp,isurf), &
781         & zmask(imaxifp,imaxjfp,isurf), &
782         & zsurf(imaxifp,imaxjfp,isurf), &
783         & zsurftmp(imaxifp,imaxjfp,isurf),  &
784         & zglamf(imaxifp+1,imaxjfp+1,isurf), &
785         & zgphif(imaxifp+1,imaxjfp+1,isurf), &
786         & igrdip1(imaxifp+1,imaxjfp+1,isurf), &
787         & igrdjp1(imaxifp+1,imaxjfp+1,isurf) &
788         & )
789
790      IF ( surfdataqc%lclim ) ALLOCATE( zclim(imaxifp,imaxjfp,isurf) )
791
792      DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf
793         iobs = jobs - surfdataqc%nsurfup
794         DO ji = 0, imaxifp
795            imodi = surfdataqc%mi(jobs) - int(imaxifp/2) + ji - 1
796           
797            !Deal with wrap around in longitude
798            IF ( imodi < 1      ) imodi = imodi + jpiglo
799            IF ( imodi > jpiglo ) imodi = imodi - jpiglo
800           
801            DO jj = 0, imaxjfp
802               imodj = surfdataqc%mj(jobs) - int(imaxjfp/2) + jj - 1
803               !If model values are out of the domain to the north/south then
804               !set them to be the edge of the domain
805               IF ( imodj < 1      ) imodj = 1
806               IF ( imodj > jpjglo ) imodj = jpjglo
807
808               igrdip1(ji+1,jj+1,iobs) = imodi
809               igrdjp1(ji+1,jj+1,iobs) = imodj
810               
811               IF ( ji >= 1 .AND. jj >= 1 ) THEN
812                  igrdi(ji,jj,iobs) = imodi
813                  igrdj(ji,jj,iobs) = imodj
814               ENDIF
815               
816            END DO
817         END DO
818      END DO
819
820      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
821         &                  igrdi, igrdj, glamt, zglam )
822      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
823         &                  igrdi, igrdj, gphit, zgphi )
824      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
825         &                  igrdi, igrdj, psurfmask, zmask )
826
827      ! At the end of the averaging period get interpolated means
828      IF ( l_timemean ) THEN
829         IF ( imeanend == 0 ) THEN
830            ALLOCATE( zsurfm(imaxifp,imaxjfp,isurf) )
831            IF (lwp) WRITE(numout,*)' Interpolating the time mean values on time step: ',kt
832            CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
833               &                  igrdi, igrdj, surfdataqc%vdmean(:,:), zsurfm )
834         ENDIF
835      ELSE
836         CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
837            &                  igrdi, igrdj, psurf, zsurf )
838      ENDIF
839
840      IF ( k2dint > 4 ) THEN         
841         CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, &
842            &                  igrdip1, igrdjp1, glamf, zglamf )
843         CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, &
844            &                  igrdip1, igrdjp1, gphif, zgphif )
845      ENDIF
846     
847      IF ( surfdataqc%lclim ) THEN
848         CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
849            &                  igrdi, igrdj, pclim, zclim )
850      ENDIF
851
852      ! At the end of the day get interpolated means
853      IF ( idayend == 0 .AND. ldnightav ) THEN
854
855         ALLOCATE( &
856            & zsurfm(imaxifp,imaxjfp,isurf)  &
857            & )
858
859         CALL obs_int_comm_2d( imaxifp,imaxjfp, isurf, kpi, kpj, igrdi, igrdj, &
860            &               surfdataqc%vdmean(:,:), zsurfm )
861
862      ENDIF
863
864      ! Loop over observations
865      DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf
866
867         iobs = jobs - surfdataqc%nsurfup
868
869         IF ( kt /= surfdataqc%mstp(jobs) ) THEN
870
871            IF(lwp) THEN
872               WRITE(numout,*)
873               WRITE(numout,*) ' E R R O R : Observation',              &
874                  &            ' time step is not consistent with the', &
875                  &            ' model time step'
876               WRITE(numout,*) ' ========='
877               WRITE(numout,*)
878               WRITE(numout,*) ' Record  = ', jobs,                &
879                  &            ' kt      = ', kt,                  &
880                  &            ' mstp    = ', surfdataqc%mstp(jobs), &
881                  &            ' ntyp    = ', surfdataqc%ntyp(jobs)
882            ENDIF
883            CALL ctl_stop( 'obs_surf_opt', 'Inconsistent time' )
884
885         ENDIF
886
887         zlam = surfdataqc%rlam(jobs)
888         zphi = surfdataqc%rphi(jobs)
889
890         IF (( ldnightav .AND. idayend == 0 ) .OR. (l_timemean .AND. imeanend == 0)) THEN
891            ! Night-time or N=kmeanstp timestep averaged data
892            zsurftmp(:,:,iobs) = zsurfm(:,:,iobs)
893         ELSE
894            zsurftmp(:,:,iobs) = zsurf(:,:,iobs)
895         ENDIF
896
897         IF (   ( .NOT. l_timemean ) .OR. & 
898             &  ( l_timemean .AND. imeanend == 0) ) THEN
899            IF ( k2dint <= 4 ) THEN
900
901               ! Get weights to interpolate the model value to the observation point
902               CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         &
903                  &                   zglam(:,:,iobs), zgphi(:,:,iobs), &
904                  &                   zmask(:,:,iobs), zweig, zobsmask )
905
906               ! Interpolate the model value to the observation point
907               CALL obs_int_h2d( 1, 1, zweig, zsurftmp(:,:,iobs), zext )
908
909               IF ( surfdataqc%lclim ) THEN 
910                  CALL obs_int_h2d( 1, 1, zweig, zclim(:,:,iobs), zclm )
911               ENDIF
912
913
914            ELSE
915
916               ! Get weights to average the model SLA to the observation footprint
917               CALL obs_avg_h2d_init( 1, 1, imaxifp, imaxjfp, k2dint, zlam,  zphi, &
918                  &                   zglam(:,:,iobs), zgphi(:,:,iobs), &
919                  &                   zglamf(:,:,iobs), zgphif(:,:,iobs), &
920                  &                   zmask(:,:,iobs), plamscl, pphiscl, &
921                  &                   lindegrees, zweig, zobsmask )
922
923               ! Average the model SST to the observation footprint
924               CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, &
925                  &              zweig, zsurftmp(:,:,iobs),  zext )
926
927               IF ( surfdataqc%lclim ) THEN 
928                  CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, &
929                     &              zweig, zclim(:,:,iobs),  zclm )
930               ENDIF
931
932            ENDIF
933
934            IF ( TRIM(surfdataqc%cvars(1)) == 'SLA' .AND. surfdataqc%nextra == 2 ) THEN
935               ! ... Remove the MDT from the SSH at the observation point to get the SLA
936               surfdataqc%rext(jobs,1) = zext(1)
937               surfdataqc%rmod(jobs,1) = surfdataqc%rext(jobs,1) - surfdataqc%rext(jobs,2)
938            ELSE
939               surfdataqc%rmod(jobs,1) = zext(1)
940            ENDIF
941
942            IF ( surfdataqc%lclim ) surfdataqc%rclm(jobs,1) = zclm(1)
943
944            IF ( zext(1) == obfillflt ) THEN
945               ! If the observation value is a fill value, set QC flag to bad
946               surfdataqc%nqc(jobs) = 4
947            ENDIF         
948         ENDIF
949         
950         IF ( surfdataqc%lclim ) surfdataqc%rclm(jobs,1) = zclm(1)
951
952         IF ( TRIM(surfdataqc%cvars(1)) == 'FBD' ) THEN
953            ! Convert radar freeboard to true freeboard (add 1/4 snow depth; 1/4 based on ratio of speed of light in vacuum compared to snow (3.0e8 vs 2.4e8 m/s))
954            surfdataqc%rext(jobs,1) = surfdataqc%robs(jobs,1) 
955            surfdataqc%robs(jobs,1) = surfdataqc%rext(jobs,1) + 0.25*surfdataqc%rext(jobs,2)
956            ! If the corrected freeboard observation is outside -0.3 to 3.0 m (CPOM) then set the QC flag to bad
957            IF ((surfdataqc%robs(jobs,1) < -0.3) .OR. (surfdataqc%robs(jobs,1) > 3.0)) THEN
958               surfdataqc%nqc(jobs) = 4
959            ENDIF           
960            ! Convert corrected freeboard to ice thickness following Tilling et al. (2016)
961            surfdataqc%robs(jobs,1) = (surfdataqc%robs(jobs,1)*rhow + surfdataqc%rext(jobs,2)*rhos)/ &
962                                      (rhow - rhoi)
963            ! Flag any negative ice thickness values as bad
964            IF (surfdataqc%robs(jobs,1) < 0.0) THEN
965               surfdataqc%nqc(jobs) = 4
966            ENDIF                                     
967         ENDIF
968         
969         IF ( zext(1) == obfillflt ) THEN
970            ! If the observation value is a fill value, set QC flag to bad
971            surfdataqc%nqc(jobs) = 4
972         ENDIF
973         
974      END DO
975
976      ! Deallocate the data for interpolation
977      DEALLOCATE( &
978         & zweig, &
979         & igrdi, &
980         & igrdj, &
981         & zglam, &
982         & zgphi, &
983         & zmask, &
984         & zsurf, &
985         & zsurftmp, &
986         & zglamf, &
987         & zgphif, &
988         & igrdip1,&
989         & igrdjp1 &
990         & )
991
992      IF ( surfdataqc%lclim ) DEALLOCATE( zclim )
993
994      ! At the end of the day also deallocate night-time mean array
995      IF (( idayend == 0 .AND. ldnightav ) .OR. ( imeanend == 0 .AND. l_timemean )) THEN
996         DEALLOCATE( &
997            & zsurfm  &
998            & )
999      ENDIF
1000
1001      surfdataqc%nsurfup = surfdataqc%nsurfup + isurf
1002
1003   END SUBROUTINE obs_surf_opt
1004
1005END MODULE obs_oper
Note: See TracBrowser for help on using the repository browser.