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

Last change on this file since 11932 was 11932, checked in by dcarneir, 4 years ago

Changes in OBS and SBC routines for sea ice thickness data assimilation
File size: 36.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#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 )
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      REAL(wp), INTENT(IN), DIMENSION(kpi,kpj) :: &
599         & psurf,  &                   ! Model surface field
600         & pclim,  &                   ! Climatological surface field         
601         & psurfmask                   ! Land-sea mask
602      LOGICAL, INTENT(IN) :: ldnightav ! Logical for averaging night-time data
603      REAL(KIND=wp), INTENT(IN) :: &
604         & plamscl, &                  ! Diameter in metres of obs footprint in E/W, N/S directions
605         & pphiscl                     ! This is the full width (rather than half-width)
606      LOGICAL, INTENT(IN) :: &
607         & lindegrees                  ! T=> plamscl and pphiscl are specified in degrees, F=> in metres
608
609      !! * Local declarations
610      INTEGER :: ji
611      INTEGER :: jj
612      INTEGER :: jobs
613      INTEGER :: inrc
614      INTEGER :: isurf
615      INTEGER :: iobs
616      INTEGER :: imaxifp, imaxjfp
617      INTEGER :: imodi, imodj
618      INTEGER :: idayend
619      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
620         & igrdi,   &
621         & igrdj,   &
622         & igrdip1, &
623         & igrdjp1
624      INTEGER, DIMENSION(:,:), SAVE, ALLOCATABLE :: &
625         & icount_night,      &
626         & imask_night
627      REAL(wp) :: zlam
628      REAL(wp) :: zphi
629      REAL(wp), DIMENSION(1) :: zext, zobsmask, zclm
630      REAL(wp) :: zdaystp
631      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
632         & zweig,  &
633         & zmask,  &
634         & zsurf,  &
635         & zsurfm, &
636         & zsurftmp, &
637         & zclim,  &
638         & zglam,  &
639         & zgphi,  &
640         & zglamf, &
641         & zgphif
642
643      REAL(wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: &
644         & zintmp,  &
645         & zouttmp, &
646         & zmeanday    ! to compute model sst in region of 24h daylight (pole)
647         
648      !------------------------------------------------------------------------
649      ! Local initialization
650      !------------------------------------------------------------------------
651      ! Record and data counters
652      inrc = kt - kit000 + 2
653      isurf = surfdataqc%nsstp(inrc)
654
655      ! Work out the maximum footprint size for the
656      ! interpolation/averaging in model grid-points - has to be even.
657
658      CALL obs_max_fpsize( k2dint, plamscl, pphiscl, lindegrees, psurfmask, imaxifp, imaxjfp )
659
660
661      IF ( ldnightav ) THEN
662
663      ! Initialize array for night mean
664         IF ( kt == 0 ) THEN
665            ALLOCATE ( icount_night(kpi,kpj) )
666            ALLOCATE ( imask_night(kpi,kpj) )
667            ALLOCATE ( zintmp(kpi,kpj) )
668            ALLOCATE ( zouttmp(kpi,kpj) )
669            ALLOCATE ( zmeanday(kpi,kpj) )
670            nday_qsr = -1   ! initialisation flag for nbc_dcy
671         ENDIF
672
673         ! Night-time means are calculated for night-time values over timesteps:
674         !  [1 <= kt <= kdaystp], [kdaystp+1 <= kt <= 2*kdaystp], .....
675         idayend = MOD( kt - kit000 + 1, kdaystp )
676
677         ! Initialize night-time mean for first timestep of the day
678         IF ( idayend == 1 .OR. kt == 0 ) THEN
679            DO jj = 1, jpj
680               DO ji = 1, jpi
681                  surfdataqc%vdmean(ji,jj) = 0.0
682                  zmeanday(ji,jj) = 0.0
683                  icount_night(ji,jj) = 0
684               END DO
685            END DO
686         ENDIF
687
688         zintmp(:,:) = 0.0
689         zouttmp(:,:) = sbc_dcy( zintmp(:,:), .TRUE. )
690         imask_night(:,:) = INT( zouttmp(:,:) )
691
692         DO jj = 1, jpj
693            DO ji = 1, jpi
694               ! Increment the temperature field for computing night mean and counter
695               surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj)  &
696                      &                    + psurf(ji,jj) * REAL( imask_night(ji,jj) )
697               zmeanday(ji,jj)          = zmeanday(ji,jj) + psurf(ji,jj)
698               icount_night(ji,jj)      = icount_night(ji,jj) + imask_night(ji,jj)
699            END DO
700         END DO
701
702         ! Compute the night-time mean at the end of the day
703         zdaystp = 1.0 / REAL( kdaystp )
704         IF ( idayend == 0 ) THEN
705            IF (lwp) WRITE(numout,*) 'Calculating surfdataqc%vdmean on time-step: ',kt
706            DO jj = 1, jpj
707               DO ji = 1, jpi
708                  ! Test if "no night" point
709                  IF ( icount_night(ji,jj) > 0 ) THEN
710                     surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj) &
711                       &                        / REAL( icount_night(ji,jj) )
712                  ELSE
713                     !At locations where there is no night (e.g. poles),
714                     ! calculate daily mean instead of night-time mean.
715                     surfdataqc%vdmean(ji,jj) = zmeanday(ji,jj) * zdaystp
716                  ENDIF
717               END DO
718            END DO
719         ENDIF
720
721      ENDIF
722
723      ! Get the data for interpolation
724
725      ALLOCATE( &
726         & zweig(imaxifp,imaxjfp,1),      &
727         & igrdi(imaxifp,imaxjfp,isurf), &
728         & igrdj(imaxifp,imaxjfp,isurf), &
729         & zglam(imaxifp,imaxjfp,isurf), &
730         & zgphi(imaxifp,imaxjfp,isurf), &
731         & zmask(imaxifp,imaxjfp,isurf), &
732         & zsurf(imaxifp,imaxjfp,isurf), &
733         & zsurftmp(imaxifp,imaxjfp,isurf),  &
734         & zglamf(imaxifp+1,imaxjfp+1,isurf), &
735         & zgphif(imaxifp+1,imaxjfp+1,isurf), &
736         & igrdip1(imaxifp+1,imaxjfp+1,isurf), &
737         & igrdjp1(imaxifp+1,imaxjfp+1,isurf) &
738         & )
739
740      IF ( surfdataqc%lclim ) ALLOCATE( zclim(imaxifp,imaxjfp,isurf) )
741
742      DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf
743         iobs = jobs - surfdataqc%nsurfup
744         DO ji = 0, imaxifp
745            imodi = surfdataqc%mi(jobs) - int(imaxifp/2) + ji - 1
746           
747            !Deal with wrap around in longitude
748            IF ( imodi < 1      ) imodi = imodi + jpiglo
749            IF ( imodi > jpiglo ) imodi = imodi - jpiglo
750           
751            DO jj = 0, imaxjfp
752               imodj = surfdataqc%mj(jobs) - int(imaxjfp/2) + jj - 1
753               !If model values are out of the domain to the north/south then
754               !set them to be the edge of the domain
755               IF ( imodj < 1      ) imodj = 1
756               IF ( imodj > jpjglo ) imodj = jpjglo
757
758               igrdip1(ji+1,jj+1,iobs) = imodi
759               igrdjp1(ji+1,jj+1,iobs) = imodj
760               
761               IF ( ji >= 1 .AND. jj >= 1 ) THEN
762                  igrdi(ji,jj,iobs) = imodi
763                  igrdj(ji,jj,iobs) = imodj
764               ENDIF
765               
766            END DO
767         END DO
768      END DO
769
770      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
771         &                  igrdi, igrdj, glamt, zglam )
772      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
773         &                  igrdi, igrdj, gphit, zgphi )
774      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
775         &                  igrdi, igrdj, psurfmask, zmask )
776      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
777         &                  igrdi, igrdj, psurf, zsurf )
778
779      IF ( k2dint > 4 ) THEN         
780         CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, &
781            &                  igrdip1, igrdjp1, glamf, zglamf )
782         CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, &
783            &                  igrdip1, igrdjp1, gphif, zgphif )
784      ENDIF
785     
786      IF ( surfdataqc%lclim ) THEN
787         CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, &
788            &                  igrdi, igrdj, pclim, zclim )
789      ENDIF
790
791      ! At the end of the day get interpolated means
792      IF ( idayend == 0 .AND. ldnightav ) THEN
793
794         ALLOCATE( &
795            & zsurfm(imaxifp,imaxjfp,isurf)  &
796            & )
797
798         CALL obs_int_comm_2d( imaxifp,imaxjfp, isurf, kpi, kpj, igrdi, igrdj, &
799            &               surfdataqc%vdmean(:,:), zsurfm )
800
801      ENDIF
802
803      ! Loop over observations
804      DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf
805
806         iobs = jobs - surfdataqc%nsurfup
807
808         IF ( kt /= surfdataqc%mstp(jobs) ) THEN
809
810            IF(lwp) THEN
811               WRITE(numout,*)
812               WRITE(numout,*) ' E R R O R : Observation',              &
813                  &            ' time step is not consistent with the', &
814                  &            ' model time step'
815               WRITE(numout,*) ' ========='
816               WRITE(numout,*)
817               WRITE(numout,*) ' Record  = ', jobs,                &
818                  &            ' kt      = ', kt,                  &
819                  &            ' mstp    = ', surfdataqc%mstp(jobs), &
820                  &            ' ntyp    = ', surfdataqc%ntyp(jobs)
821            ENDIF
822            CALL ctl_stop( 'obs_surf_opt', 'Inconsistent time' )
823
824         ENDIF
825
826         zlam = surfdataqc%rlam(jobs)
827         zphi = surfdataqc%rphi(jobs)
828
829         IF ( ldnightav .AND. idayend == 0 ) THEN
830            ! Night-time averaged data
831            zsurftmp(:,:,iobs) = zsurfm(:,:,iobs)
832         ELSE
833            zsurftmp(:,:,iobs) = zsurf(:,:,iobs)
834         ENDIF
835
836         IF ( k2dint <= 4 ) THEN
837
838            ! Get weights to interpolate the model value to the observation point
839            CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         &
840               &                   zglam(:,:,iobs), zgphi(:,:,iobs), &
841               &                   zmask(:,:,iobs), zweig, zobsmask )
842
843            ! Interpolate the model value to the observation point
844            CALL obs_int_h2d( 1, 1, zweig, zsurftmp(:,:,iobs), zext )
845
846            IF ( surfdataqc%lclim ) THEN 
847               CALL obs_int_h2d( 1, 1, zweig, zclim(:,:,iobs), zclm )
848            ENDIF
849
850
851         ELSE
852
853            ! Get weights to average the model SLA to the observation footprint
854            CALL obs_avg_h2d_init( 1, 1, imaxifp, imaxjfp, k2dint, zlam,  zphi, &
855               &                   zglam(:,:,iobs), zgphi(:,:,iobs), &
856               &                   zglamf(:,:,iobs), zgphif(:,:,iobs), &
857               &                   zmask(:,:,iobs), plamscl, pphiscl, &
858               &                   lindegrees, zweig, zobsmask )
859
860            ! Average the model SST to the observation footprint
861            CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, &
862               &              zweig, zsurftmp(:,:,iobs),  zext )
863
864            IF ( surfdataqc%lclim ) THEN 
865               CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, &
866                  &              zweig, zclim(:,:,iobs),  zclm )
867            ENDIF
868
869         ENDIF
870
871         IF ( TRIM(surfdataqc%cvars(1)) == 'SLA' .AND. surfdataqc%nextra == 2 ) THEN
872            ! ... Remove the MDT from the SSH at the observation point to get the SLA
873            surfdataqc%rext(jobs,1) = zext(1)
874            surfdataqc%rmod(jobs,1) = surfdataqc%rext(jobs,1) - surfdataqc%rext(jobs,2)
875         ELSE
876            surfdataqc%rmod(jobs,1) = zext(1)
877         ENDIF
878         
879         IF ( surfdataqc%lclim ) surfdataqc%rclm(jobs,1) = zclm(1)
880
881         IF ( TRIM(surfdataqc%cvars(1)) == 'FBD' ) THEN
882            ! 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))
883            surfdataqc%rext(jobs,1) = surfdataqc%robs(jobs,1) 
884            surfdataqc%robs(jobs,1) = surfdataqc%rext(jobs,1) + 0.25*surfdataqc%rext(jobs,2)
885            ! If the corrected freeboard observation is outside -0.3 to 3.0 m (CPOM) then set the QC flag to bad
886            IF ((surfdataqc%robs(jobs,1) < -0.3) .OR. (surfdataqc%robs(jobs,1) > 3.0)) THEN
887               surfdataqc%nqc(jobs) = 4
888            ENDIF           
889            ! Convert corrected freeboard to ice thickness following Tilling et al. (2016)
890            surfdataqc%robs(jobs,1) = (surfdataqc%robs(jobs,1)*rhow + surfdataqc%rext(jobs,2)*rhos)/ &
891                                      (rhow - rhoi)
892            ! Flag any negative ice thickness values as bad
893            IF (surfdataqc%robs(jobs,1) < 0.0) THEN
894               surfdataqc%nqc(jobs) = 4
895            ENDIF                                     
896         ENDIF
897         
898         IF ( zext(1) == obfillflt ) THEN
899            ! If the observation value is a fill value, set QC flag to bad
900            surfdataqc%nqc(jobs) = 4
901         ENDIF
902
903      END DO
904
905      ! Deallocate the data for interpolation
906      DEALLOCATE( &
907         & zweig, &
908         & igrdi, &
909         & igrdj, &
910         & zglam, &
911         & zgphi, &
912         & zmask, &
913         & zsurf, &
914         & zsurftmp, &
915         & zglamf, &
916         & zgphif, &
917         & igrdip1,&
918         & igrdjp1 &
919         & )
920
921      IF ( surfdataqc%lclim ) DEALLOCATE( zclim )
922
923      ! At the end of the day also deallocate night-time mean array
924      IF ( idayend == 0 .AND. ldnightav ) THEN
925         DEALLOCATE( &
926            & zsurfm  &
927            & )
928      ENDIF
929
930      surfdataqc%nsurfup = surfdataqc%nsurfup + isurf
931
932   END SUBROUTINE obs_surf_opt
933
934END MODULE obs_oper
Note: See TracBrowser for help on using the repository browser.