source: branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90 @ 10276

Last change on this file since 10276 was 10276, checked in by emmafiedler, 3 years ago

Freeboard assimilation updates

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