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

source: branches/UKMO/dev_r5785_SSS_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90 @ 7837

Last change on this file since 7837 was 7837, checked in by mattmartin, 7 years ago

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