Ignore:
Timestamp:
2019-07-01T12:44:06+02:00 (15 months ago)
Author:
jcastill
Message:

Copy of branch branches/UKMO/dev_r5518_obs_oper_update@11130 without namelist_ref changes to allow merging with coupling and biogeochemistry branches

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90

    r4245 r11202  
    77 
    88   !!---------------------------------------------------------------------- 
    9    !!   obs_pro_opt :    Compute the model counterpart of temperature and 
    10    !!                    salinity observations from profiles 
    11    !!   obs_sla_opt :    Compute the model counterpart of sea level anomaly 
    12    !!                    observations 
    13    !!   obs_sst_opt :    Compute the model counterpart of sea surface temperature 
    14    !!                    observations 
    15    !!   obs_sss_opt :    Compute the model counterpart of sea surface salinity 
    16    !!                    observations 
    17    !!   obs_seaice_opt : Compute the model counterpart of sea ice concentration 
    18    !!                    observations 
    19    !! 
    20    !!   obs_vel_opt :    Compute the model counterpart of zonal and meridional 
    21    !!                    components of velocity from observations. 
     9   !!   obs_prof_opt :    Compute the model counterpart of profile data 
     10   !!   obs_surf_opt :    Compute the model counterpart of surface data 
    2211   !!---------------------------------------------------------------------- 
    2312 
    24    !! * Modules used    
     13   !! * Modules used 
    2514   USE par_kind, ONLY : &         ! Precision variables 
    2615      & wp 
    2716   USE in_out_manager             ! I/O manager 
    2817   USE obs_inter_sup              ! Interpolation support 
    29    USE obs_inter_h2d, ONLY : &    ! Horizontal interpolation to the observation pt 
     18   USE obs_inter_h2d, ONLY : &    ! Horizontal interpolation to the obs pt 
    3019      & obs_int_h2d, & 
    3120      & obs_int_h2d_init 
    32    USE obs_inter_z1d, ONLY : &    ! Vertical interpolation to the observation pt 
     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 
    3326      & obs_int_z1d,    & 
    3427      & obs_int_z1d_spl 
    35    USE obs_const,  ONLY :     & 
    36       & obfillflt      ! Fillvalue    
     28   USE obs_const,  ONLY :    &    ! Obs fill value 
     29      & obfillflt 
    3730   USE dom_oce,       ONLY : & 
    38       & glamt, glamu, glamv, & 
    39       & gphit, gphiu, gphiv 
    40    USE lib_mpp,       ONLY : & 
     31      & glamt, glamf, & 
     32      & gphit, gphif 
     33   USE lib_mpp,       ONLY : &    ! Warning and stopping routines 
    4134      & 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      
    4239 
    4340   IMPLICIT NONE 
     
    4643   PRIVATE 
    4744 
    48    PUBLIC obs_pro_opt, &  ! Compute the model counterpart of profile observations 
    49       &   obs_sla_opt, &  ! Compute the model counterpart of SLA observations 
    50       &   obs_sst_opt, &  ! Compute the model counterpart of SST observations 
    51       &   obs_sss_opt, &  ! Compute the model counterpart of SSS observations 
    52       &   obs_seaice_opt, & 
    53       &   obs_vel_opt     ! Compute the model counterpart of velocity profile data 
    54  
    55    INTEGER, PARAMETER, PUBLIC :: imaxavtypes = 20 ! Max number of daily avgd obs types 
     45   PUBLIC obs_prof_opt, &  ! Compute the model counterpart of profile obs 
     46      &   obs_surf_opt     ! Compute the model counterpart of surface obs 
     47 
     48   INTEGER, PARAMETER, PUBLIC :: & 
     49      & imaxavtypes = 20   ! Max number of daily avgd obs types 
    5650 
    5751   !!---------------------------------------------------------------------- 
     
    6155   !!---------------------------------------------------------------------- 
    6256 
     57   !! * Substitutions  
     58#  include "domzgr_substitute.h90"  
    6359CONTAINS 
    6460 
    65    SUBROUTINE obs_pro_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, & 
    66       &                    ptn, psn, pgdept, ptmask, k1dint, k2dint, & 
    67       &                    kdailyavtypes ) 
     61 
     62   SUBROUTINE obs_prof_opt( prodatqc, kt, kpi, kpj, kpk, & 
     63      &                     kit000, kdaystp, kvar,       & 
     64      &                     pvar, pgdept, pgdepw,        & 
     65      &                     pmask,                       &   
     66      &                     plam, pphi,                  & 
     67      &                     k1dint, k2dint, kdailyavtypes ) 
     68 
    6869      !!----------------------------------------------------------------------- 
    6970      !! 
     
    7879      !! 
    7980      !!    First, a vertical profile of horizontally interpolated model 
    80       !!    now temperatures is computed at the obs (lon, lat) point. 
     81      !!    now values is computed at the obs (lon, lat) point. 
    8182      !!    Several horizontal interpolation schemes are available: 
    8283      !!        - distance-weighted (great circle) (k2dint = 0) 
     
    8687      !!        - polynomial (quadrilateral grid)  (k2dint = 4) 
    8788      !! 
    88       !!    Next, the vertical temperature profile is interpolated to the 
     89      !!    Next, the vertical profile is interpolated to the 
    8990      !!    data depth points. Two vertical interpolation schemes are 
    9091      !!    available: 
     
    9697      !!    routine. 
    9798      !! 
    98       !!    For ENACT moored buoy data (e.g., TAO), the model equivalent is 
     99      !!    If the logical is switched on, the model equivalent is 
    99100      !!    a daily mean model temperature field. So, we first compute 
    100101      !!    the mean, then interpolate only at the end of the day. 
    101102      !! 
    102       !!    Note: the in situ temperature observations must be converted 
     103      !!    Note: in situ temperature observations must be converted 
    103104      !!    to potential temperature (the model variable) prior to 
    104105      !!    assimilation.  
    105       !!?????????????????????????????????????????????????????????????? 
    106       !!    INCLUDE POTENTIAL TEMP -> IN SITU TEMP IN OBS OPERATOR??? 
    107       !!?????????????????????????????????????????????????????????????? 
    108106      !! 
    109107      !! ** Action  : 
     
    115113      !!      ! 07-01 (K. Mogensen) Merge of temperature and salinity 
    116114      !!      ! 07-03 (K. Mogensen) General handling of profiles 
     115      !!      ! 15-02 (M. Martin) Combined routine for all profile types 
     116      !!      ! 17-02 (M. Martin) Include generalised vertical coordinate changes 
    117117      !!----------------------------------------------------------------------- 
    118    
     118 
    119119      !! * Modules used 
    120120      USE obs_profiles_def ! Definition of storage space for profile obs. 
     
    123123 
    124124      !! * Arguments 
    125       TYPE(obs_prof), INTENT(INOUT) :: prodatqc  ! Subset of profile data not failing screening 
    126       INTEGER, INTENT(IN) :: kt        ! Time step 
    127       INTEGER, INTENT(IN) :: kpi       ! Model grid parameters 
     125      TYPE(obs_prof), INTENT(INOUT) :: & 
     126         & prodatqc                  ! Subset of profile data passing QC 
     127      INTEGER, INTENT(IN) :: kt      ! Time step 
     128      INTEGER, INTENT(IN) :: kpi     ! Model grid parameters 
    128129      INTEGER, INTENT(IN) :: kpj 
    129130      INTEGER, INTENT(IN) :: kpk 
    130       INTEGER, INTENT(IN) :: kit000    ! Number of the first time step  
    131                                        !   (kit000-1 = restart time) 
    132       INTEGER, INTENT(IN) :: k1dint    ! Vertical interpolation type (see header) 
    133       INTEGER, INTENT(IN) :: k2dint    ! Horizontal interpolation type (see header) 
    134       INTEGER, INTENT(IN) :: kdaystp   ! Number of time steps per day                     
     131      INTEGER, INTENT(IN) :: kit000  ! Number of the first time step 
     132                                     !   (kit000-1 = restart time) 
     133      INTEGER, INTENT(IN) :: k1dint  ! Vertical interpolation type (see header) 
     134      INTEGER, INTENT(IN) :: k2dint  ! Horizontal interpolation type (see header) 
     135      INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day 
     136      INTEGER, INTENT(IN) :: kvar    ! Number of variable in prodatqc 
    135137      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 
    136          & ptn,    &    ! Model temperature field 
    137          & psn,    &    ! Model salinity field 
    138          & ptmask       ! Land-sea mask 
    139       REAL(KIND=wp), INTENT(IN), DIMENSION(kpk) :: & 
    140          & pgdept       ! Model array of depth levels 
     138         & pvar,   &                 ! Model field for variable 
     139         & pmask                     ! Land-sea mask for variable 
     140      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 
     141         & plam,   &                 ! Model longitudes for variable 
     142         & pphi                      ! Model latitudes for variable 
     143      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &  
     144         & pgdept, &                 ! Model array of depth T levels  
     145         & pgdepw                    ! Model array of depth W levels  
    141146      INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 
    142          & kdailyavtypes! Types for daily averages 
     147         & kdailyavtypes             ! Types for daily averages 
     148 
    143149      !! * Local declarations 
    144150      INTEGER ::   ji 
     
    152158      INTEGER ::   iend 
    153159      INTEGER ::   iobs 
     160      INTEGER ::   iin, ijn, ikn, ik   ! looping indices over interpolation nodes  
     161      INTEGER ::   inum_obs 
    154162      INTEGER, DIMENSION(imaxavtypes) :: & 
    155163         & idailyavtypes 
     164      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 
     165         & igrdi, & 
     166         & igrdj 
     167      INTEGER, ALLOCATABLE, DIMENSION(:) :: iv_indic 
     168 
    156169      REAL(KIND=wp) :: zlam 
    157170      REAL(KIND=wp) :: zphi 
    158171      REAL(KIND=wp) :: zdaystp 
    159172      REAL(KIND=wp), DIMENSION(kpk) :: & 
    160          & zobsmask, & 
    161173         & zobsk,    & 
    162174         & zobs2k 
    163       REAL(KIND=wp), DIMENSION(2,2,kpk) :: & 
     175      REAL(KIND=wp), DIMENSION(2,2,1) :: & 
     176         & zweig1, & 
    164177         & zweig 
    165178      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: & 
    166          & zmask, & 
    167          & zintt, & 
    168          & zints, & 
    169          & zinmt, & 
    170          & zinms 
     179         & zmask,  & 
     180         & zint & 
     181         & zinm,  & 
     182         & zgdept, &  
     183         & zgdepw 
    171184      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 
    172          & zglam, & 
     185         & zglam,  & 
    173186         & zgphi 
    174       INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 
    175          & igrdi, & 
    176          & igrdj 
     187      REAL(KIND=wp), DIMENSION(1) :: zmsk 
     188      REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner 
     189 
     190      LOGICAL :: ld_dailyav 
    177191 
    178192      !------------------------------------------------------------------------ 
    179193      ! Local initialization  
    180194      !------------------------------------------------------------------------ 
    181       ! ... Record and data counters 
     195      ! Record and data counters 
    182196      inrc = kt - kit000 + 2 
    183197      ipro = prodatqc%npstp(inrc) 
    184   
     198 
    185199      ! Daily average types 
     200      ld_dailyav = .FALSE. 
    186201      IF ( PRESENT(kdailyavtypes) ) THEN 
    187202         idailyavtypes(:) = kdailyavtypes(:) 
     203         IF ( ANY (idailyavtypes(:) /= -1) ) ld_dailyav = .TRUE. 
    188204      ELSE 
    189205         idailyavtypes(:) = -1 
    190206      ENDIF 
    191207 
    192       ! Initialize daily mean for first timestep 
     208      ! Daily means are calculated for values over timesteps: 
     209      !  [1 <= kt <= kdaystp], [kdaystp+1 <= kt <= 2*kdaystp], ... 
    193210      idayend = MOD( kt - kit000 + 1, kdaystp ) 
    194211 
    195       ! Added kt == 0 test to catch restart case  
    196       IF ( idayend == 1 .OR. kt == 0) THEN 
    197          IF (lwp) WRITE(numout,*) 'Reset prodatqc%vdmean on time-step: ',kt 
     212      IF ( ld_dailyav ) THEN 
     213 
     214         ! Initialize daily mean for first timestep of the day 
     215         IF ( idayend == 1 .OR. kt == 0 ) THEN 
     216            DO jk = 1, jpk 
     217               DO jj = 1, jpj 
     218                  DO ji = 1, jpi 
     219                     prodatqc%vdmean(ji,jj,jk,kvar) = 0.0 
     220                  END DO 
     221               END DO 
     222            END DO 
     223         ENDIF 
     224 
    198225         DO jk = 1, jpk 
    199226            DO jj = 1, jpj 
    200227               DO ji = 1, jpi 
    201                   prodatqc%vdmean(ji,jj,jk,1) = 0.0 
    202                   prodatqc%vdmean(ji,jj,jk,2) = 0.0 
     228                  ! Increment field for computing daily mean 
     229                  prodatqc%vdmean(ji,jj,jk,kvar) = prodatqc%vdmean(ji,jj,jk,kvar) & 
     230                     &                           + pvar(ji,jj,jk) 
    203231               END DO 
    204232            END DO 
    205233         END DO 
    206       ENDIF 
    207  
    208       DO jk = 1, jpk 
    209          DO jj = 1, jpj 
    210             DO ji = 1, jpi 
    211                ! Increment the temperature field for computing daily mean 
    212                prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) & 
    213                   &                        + ptn(ji,jj,jk) 
    214                ! Increment the salinity field for computing daily mean 
    215                prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 
    216                   &                        + psn(ji,jj,jk) 
    217             END DO 
    218          END DO 
    219       END DO 
    220     
    221       ! Compute the daily mean at the end of day 
    222       zdaystp = 1.0 / REAL( kdaystp ) 
    223       IF ( idayend == 0 ) THEN 
    224          DO jk = 1, jpk 
    225             DO jj = 1, jpj 
    226                DO ji = 1, jpi 
    227                   prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) & 
    228                      &                        * zdaystp 
    229                   prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 
    230                   &                           * zdaystp 
     234 
     235         ! Compute the daily mean at the end of day 
     236         zdaystp = 1.0 / REAL( kdaystp ) 
     237         IF ( idayend == 0 ) THEN 
     238            IF (lwp) WRITE(numout,*) 'Calculating prodatqc%vdmean on time-step: ',kt 
     239            CALL FLUSH(numout) 
     240            DO jk = 1, jpk 
     241               DO jj = 1, jpj 
     242                  DO ji = 1, jpi 
     243                     prodatqc%vdmean(ji,jj,jk,kvar) = prodatqc%vdmean(ji,jj,jk,kvar) & 
     244                        &                           * zdaystp 
     245                  END DO 
    231246               END DO 
    232247            END DO 
    233          END DO 
     248         ENDIF 
     249 
    234250      ENDIF 
    235251 
    236252      ! Get the data for interpolation 
    237253      ALLOCATE( & 
    238          & igrdi(2,2,ipro),      & 
    239          & igrdj(2,2,ipro),      & 
    240          & zglam(2,2,ipro),      & 
    241          & zgphi(2,2,ipro),      & 
    242          & zmask(2,2,kpk,ipro),  & 
    243          & zintt(2,2,kpk,ipro),  & 
    244          & zints(2,2,kpk,ipro)   & 
     254         & igrdi(2,2,ipro),       & 
     255         & igrdj(2,2,ipro),       & 
     256         & zglam(2,2,ipro),       & 
     257         & zgphi(2,2,ipro),       & 
     258         & zmask(2,2,kpk,ipro),   & 
     259         & zint(2,2,kpk,ipro),    & 
     260         & zgdept(2,2,kpk,ipro),  &  
     261         & zgdepw(2,2,kpk,ipro)   &  
    245262         & ) 
    246263 
    247264      DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro 
    248265         iobs = jobs - prodatqc%nprofup 
    249          igrdi(1,1,iobs) = prodatqc%mi(jobs,1)-1 
    250          igrdj(1,1,iobs) = prodatqc%mj(jobs,1)-1 
    251          igrdi(1,2,iobs) = prodatqc%mi(jobs,1)-1 
    252          igrdj(1,2,iobs) = prodatqc%mj(jobs,1) 
    253          igrdi(2,1,iobs) = prodatqc%mi(jobs,1) 
    254          igrdj(2,1,iobs) = prodatqc%mj(jobs,1)-1 
    255          igrdi(2,2,iobs) = prodatqc%mi(jobs,1) 
    256          igrdj(2,2,iobs) = prodatqc%mj(jobs,1) 
     266         igrdi(1,1,iobs) = prodatqc%mi(jobs,kvar)-1 
     267         igrdj(1,1,iobs) = prodatqc%mj(jobs,kvar)-1 
     268         igrdi(1,2,iobs) = prodatqc%mi(jobs,kvar)-1 
     269         igrdj(1,2,iobs) = prodatqc%mj(jobs,kvar) 
     270         igrdi(2,1,iobs) = prodatqc%mi(jobs,kvar) 
     271         igrdj(2,1,iobs) = prodatqc%mj(jobs,kvar)-1 
     272         igrdi(2,2,iobs) = prodatqc%mi(jobs,kvar) 
     273         igrdj(2,2,iobs) = prodatqc%mj(jobs,kvar) 
    257274      END DO 
    258275 
    259       CALL obs_int_comm_2d( 2, 2, ipro, igrdi, igrdj, glamt, zglam ) 
    260       CALL obs_int_comm_2d( 2, 2, ipro, igrdi, igrdj, gphit, zgphi ) 
    261       CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, ptmask,zmask ) 
    262       CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, ptn,   zintt ) 
    263       CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, psn,   zints ) 
     276      ! Initialise depth arrays 
     277      zgdept(:,:,:,:) = 0.0 
     278      zgdepw(:,:,:,:) = 0.0 
     279 
     280      CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi, igrdj, plam, zglam ) 
     281      CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi, igrdj, pphi, zgphi ) 
     282      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pmask, zmask ) 
     283      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pvar,   zint ) 
     284 
     285      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pgdept, zgdept )  
     286      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pgdepw, zgdepw )  
    264287 
    265288      ! At the end of the day also get interpolated means 
    266       IF ( idayend == 0 ) THEN 
    267  
    268          ALLOCATE( & 
    269             & zinmt(2,2,kpk,ipro),  & 
    270             & zinms(2,2,kpk,ipro)   & 
    271             & ) 
    272  
    273          CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, & 
    274             &                  prodatqc%vdmean(:,:,:,1), zinmt ) 
    275          CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdi, igrdj, & 
    276             &                  prodatqc%vdmean(:,:,:,2), zinms ) 
     289      IF ( ld_dailyav .AND. idayend == 0 ) THEN 
     290 
     291         ALLOCATE( zinm(2,2,kpk,ipro) ) 
     292 
     293         CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, & 
     294            &                  prodatqc%vdmean(:,:,:,kvar), zinm ) 
    277295 
    278296      ENDIF 
    279297 
     298      ! Return if no observations to process  
     299      ! Has to be done after comm commands to ensure processors  
     300      ! stay in sync  
     301      IF ( ipro == 0 ) RETURN  
     302 
    280303      DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro 
    281304 
     
    283306 
    284307         IF ( kt /= prodatqc%mstp(jobs) ) THEN 
    285              
     308 
    286309            IF(lwp) THEN 
    287310               WRITE(numout,*) 
     
    298321            CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' ) 
    299322         ENDIF 
    300           
     323 
    301324         zlam = prodatqc%rlam(jobs) 
    302325         zphi = prodatqc%rphi(jobs) 
     326 
     327         ! Horizontal weights  
     328         ! Masked values are calculated later.   
     329         IF ( prodatqc%npvend(jobs,kvar) > 0 ) THEN 
     330 
     331            CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,     & 
     332               &                   zglam(:,:,iobs), zgphi(:,:,iobs), & 
     333               &                   zmask(:,:,1,iobs), zweig1, zmsk ) 
     334 
     335         ENDIF 
     336 
     337         IF ( prodatqc%npvend(jobs,kvar) > 0 ) THEN 
     338 
     339            zobsk(:) = obfillflt 
     340 
     341            IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN 
     342 
     343               IF ( idayend == 0 )  THEN 
     344                  ! Daily averaged data 
     345 
     346                  ! vertically interpolate all 4 corners  
     347                  ista = prodatqc%npvsta(jobs,kvar)  
     348                  iend = prodatqc%npvend(jobs,kvar)  
     349                  inum_obs = iend - ista + 1  
     350                  ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs))  
     351 
     352                  DO iin=1,2  
     353                     DO ijn=1,2  
     354 
     355                        IF ( k1dint == 1 ) THEN  
     356                           CALL obs_int_z1d_spl( kpk, &  
     357                              &     zinm(iin,ijn,:,iobs), &  
     358                              &     zobs2k, zgdept(iin,ijn,:,iobs), &  
     359                              &     zmask(iin,ijn,:,iobs))  
     360                        ENDIF  
     361        
     362                        CALL obs_level_search(kpk, &  
     363                           &    zgdept(iin,ijn,:,iobs), &  
     364                           &    inum_obs, prodatqc%var(kvar)%vdep(ista:iend), &  
     365                           &    iv_indic)  
     366 
     367                        CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &  
     368                           &    prodatqc%var(kvar)%vdep(ista:iend), &  
     369                           &    zinm(iin,ijn,:,iobs), &  
     370                           &    zobs2k, interp_corner(iin,ijn,:), &  
     371                           &    zgdept(iin,ijn,:,iobs), &  
     372                           &    zmask(iin,ijn,:,iobs))  
     373        
     374                     ENDDO  
     375                  ENDDO  
     376 
     377               ENDIF !idayend 
     378 
     379            ELSE    
     380 
     381               ! Point data  
     382      
     383               ! vertically interpolate all 4 corners  
     384               ista = prodatqc%npvsta(jobs,kvar)  
     385               iend = prodatqc%npvend(jobs,kvar)  
     386               inum_obs = iend - ista + 1  
     387               ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs))  
     388               DO iin=1,2   
     389                  DO ijn=1,2  
     390                     
     391                     IF ( k1dint == 1 ) THEN  
     392                        CALL obs_int_z1d_spl( kpk, &  
     393                           &    zint(iin,ijn,:,iobs),&  
     394                           &    zobs2k, zgdept(iin,ijn,:,iobs), &  
     395                           &    zmask(iin,ijn,:,iobs))  
     396   
     397                     ENDIF  
     398        
     399                     CALL obs_level_search(kpk, &  
     400                         &        zgdept(iin,ijn,:,iobs),&  
     401                         &        inum_obs, prodatqc%var(kvar)%vdep(ista:iend), &  
     402                         &        iv_indic)  
     403 
     404                     CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs,     &  
     405                         &          prodatqc%var(kvar)%vdep(ista:iend),     &  
     406                         &          zint(iin,ijn,:,iobs),            &  
     407                         &          zobs2k,interp_corner(iin,ijn,:), &  
     408                         &          zgdept(iin,ijn,:,iobs),         &  
     409                         &          zmask(iin,ijn,:,iobs) )       
    303410          
    304          ! Horizontal weights and vertical mask 
    305  
    306          IF ( ( prodatqc%npvend(jobs,1) > 0 ) .OR. & 
    307             & ( prodatqc%npvend(jobs,2) > 0 ) ) THEN 
    308  
    309             CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi,     & 
    310                &                   zglam(:,:,iobs), zgphi(:,:,iobs), & 
    311                &                   zmask(:,:,:,iobs), zweig, zobsmask ) 
    312  
    313          ENDIF 
    314  
    315          IF ( prodatqc%npvend(jobs,1) > 0 ) THEN 
    316  
    317             zobsk(:) = obfillflt 
    318  
    319        IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN 
    320  
    321                IF ( idayend == 0 )  THEN 
     411                  ENDDO  
     412               ENDDO  
     413              
     414            ENDIF  
     415 
     416            !-------------------------------------------------------------  
     417            ! Compute the horizontal interpolation for every profile level  
     418            !-------------------------------------------------------------  
     419              
     420            DO ikn=1,inum_obs  
     421               iend=ista+ikn-1 
    322422                   
    323                   ! Daily averaged moored buoy (MRB) data 
    324                    
    325                   CALL obs_int_h2d( kpk, kpk,      & 
    326                      &              zweig, zinmt(:,:,:,iobs), zobsk ) 
    327                    
    328                    
    329                ELSE 
    330                 
    331                   CALL ctl_stop( ' A nonzero' //     & 
    332                      &           ' number of profile T BUOY data should' // & 
    333                      &           ' only occur at the end of a given day' ) 
    334  
    335                ENDIF 
    336            
    337             ELSE  
    338                 
    339                ! Point data 
    340  
    341                CALL obs_int_h2d( kpk, kpk,      & 
    342                   &              zweig, zintt(:,:,:,iobs), zobsk ) 
    343  
    344             ENDIF 
    345  
    346             !------------------------------------------------------------- 
    347             ! Compute vertical second-derivative of the interpolating  
    348             ! polynomial at obs points 
    349             !------------------------------------------------------------- 
    350              
    351             IF ( k1dint == 1 ) THEN 
    352                CALL obs_int_z1d_spl( kpk, zobsk, zobs2k,   & 
    353                   &                  pgdept, zobsmask ) 
    354             ENDIF 
    355              
    356             !----------------------------------------------------------------- 
    357             !  Vertical interpolation to the observation point 
    358             !----------------------------------------------------------------- 
    359             ista = prodatqc%npvsta(jobs,1) 
    360             iend = prodatqc%npvend(jobs,1) 
    361             CALL obs_int_z1d( kpk,                & 
    362                & prodatqc%var(1)%mvk(ista:iend),  & 
    363                & k1dint, iend - ista + 1,         & 
    364                & prodatqc%var(1)%vdep(ista:iend), & 
    365                & zobsk, zobs2k,                   & 
    366                & prodatqc%var(1)%vmod(ista:iend), & 
    367                & pgdept, zobsmask ) 
    368  
    369          ENDIF 
    370  
    371          IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 
    372  
    373             zobsk(:) = obfillflt 
    374  
    375             IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN 
    376  
    377                IF ( idayend == 0 )  THEN 
    378  
    379                   ! Daily averaged moored buoy (MRB) data 
    380                    
    381                   CALL obs_int_h2d( kpk, kpk,      & 
    382                      &              zweig, zinms(:,:,:,iobs), zobsk ) 
    383                    
    384                ELSE 
    385  
    386                   CALL ctl_stop( ' A nonzero' //     & 
    387                      &           ' number of profile S BUOY data should' // & 
    388                      &           ' only occur at the end of a given day' ) 
    389  
    390                ENDIF 
    391  
    392             ELSE 
    393                 
    394                ! Point data 
    395  
    396                CALL obs_int_h2d( kpk, kpk,      & 
    397                   &              zweig, zints(:,:,:,iobs), zobsk ) 
    398  
    399             ENDIF 
    400  
    401  
    402             !------------------------------------------------------------- 
    403             ! Compute vertical second-derivative of the interpolating  
    404             ! polynomial at obs points 
    405             !------------------------------------------------------------- 
    406              
    407             IF ( k1dint == 1 ) THEN 
    408                CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, & 
    409                   &                  pgdept, zobsmask ) 
    410             ENDIF 
    411              
    412             !---------------------------------------------------------------- 
    413             !  Vertical interpolation to the observation point 
    414             !---------------------------------------------------------------- 
    415             ista = prodatqc%npvsta(jobs,2) 
    416             iend = prodatqc%npvend(jobs,2) 
    417             CALL obs_int_z1d( kpk, & 
    418                & prodatqc%var(2)%mvk(ista:iend),& 
    419                & k1dint, iend - ista + 1, & 
    420                & prodatqc%var(2)%vdep(ista:iend),& 
    421                & zobsk, zobs2k, & 
    422                & prodatqc%var(2)%vmod(ista:iend),& 
    423                & pgdept, zobsmask ) 
    424  
    425          ENDIF 
    426  
    427       END DO 
     423               zweig(:,:,1) = 0._wp  
     424    
     425               ! This code forces the horizontal weights to be   
     426               ! zero IF the observation is below the bottom of the   
     427               ! corners of the interpolation nodes, Or if it is in   
     428               ! the mask. This is important for observations near   
     429               ! steep bathymetry  
     430               DO iin=1,2  
     431                  DO ijn=1,2  
     432      
     433                     depth_loop: DO ik=kpk,2,-1  
     434                        IF(zmask(iin,ijn,ik-1,iobs ) > 0.9 )THEN    
     435                             
     436                           zweig(iin,ijn,1) = &   
     437                              & zweig1(iin,ijn,1) * &  
     438                              & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) &  
     439                              &  - prodatqc%var(kvar)%vdep(iend)),0._wp)  
     440                             
     441                           EXIT depth_loop 
     442 
     443                        ENDIF  
     444 
     445                     ENDDO depth_loop 
     446      
     447                  ENDDO  
     448               ENDDO  
     449    
     450               CALL obs_int_h2d( 1, 1, zweig, interp_corner(:,:,ikn), &  
     451                  &              prodatqc%var(kvar)%vmod(iend:iend) )  
     452 
     453                  ! Set QC flag for any observations found below the bottom 
     454                  ! needed as the check here is more strict than that in obs_prep 
     455               IF (sum(zweig) == 0.0_wp) prodatqc%var(kvar)%nvqc(iend:iend)=4 
    428456  
     457            ENDDO  
     458  
     459            DEALLOCATE(interp_corner,iv_indic)  
     460           
     461         ENDIF 
     462 
     463      ENDDO 
     464 
    429465      ! Deallocate the data for interpolation 
    430       DEALLOCATE( & 
    431          & igrdi, & 
    432          & igrdj, & 
    433          & zglam, & 
    434          & zgphi, & 
    435          & zmask, & 
    436          & zintt, & 
    437          & zints  & 
     466      DEALLOCATE(  & 
     467         & igrdi,  & 
     468         & igrdj,  & 
     469         & zglam,  & 
     470         & zgphi,  & 
     471         & zmask,  & 
     472         & zint,   & 
     473         & zgdept, & 
     474         & zgdepw  & 
    438475         & ) 
     476 
    439477      ! At the end of the day also get interpolated means 
    440       IF ( idayend == 0 ) THEN 
    441          DEALLOCATE( & 
    442             & zinmt,  & 
    443             & zinms   & 
    444             & ) 
     478      IF ( ld_dailyav .AND. idayend == 0 ) THEN 
     479         DEALLOCATE( zinm ) 
    445480      ENDIF 
    446481 
    447       prodatqc%nprofup = prodatqc%nprofup + ipro  
    448        
    449    END SUBROUTINE obs_pro_opt 
    450  
    451    SUBROUTINE obs_sla_opt( sladatqc, kt, kpi, kpj, kit000, & 
    452       &                    psshn, psshmask, k2dint ) 
     482      IF ( kvar == prodatqc%nvar ) THEN 
     483         prodatqc%nprofup = prodatqc%nprofup + ipro  
     484      ENDIF 
     485 
     486   END SUBROUTINE obs_prof_opt 
     487 
     488   SUBROUTINE obs_surf_opt( surfdataqc, kt, kpi, kpj,            & 
     489      &                     kit000, kdaystp, psurf, psurfmask,   & 
     490      &                     k2dint, ldnightav, plamscl, pphiscl, & 
     491      &                     lindegrees ) 
     492 
    453493      !!----------------------------------------------------------------------- 
    454494      !! 
    455       !!                     ***  ROUTINE obs_sla_opt  *** 
    456       !! 
    457       !! ** Purpose : Compute the model counterpart of sea level anomaly 
     495      !!                     ***  ROUTINE obs_surf_opt  *** 
     496      !! 
     497      !! ** Purpose : Compute the model counterpart of surface 
    458498      !!              data by interpolating from the model grid to the  
    459499      !!              observation point. 
     
    462502      !!              the model values at the corners of the surrounding grid box. 
    463503      !! 
    464       !!    The now model SSH is first computed at the obs (lon, lat) point. 
     504      !!    The new model value is first computed at the obs (lon, lat) point. 
    465505      !! 
    466506      !!    Several horizontal interpolation schemes are available: 
     
    470510      !!        - bilinear (quadrilateral grid)    (k2dint = 3) 
    471511      !!        - polynomial (quadrilateral grid)  (k2dint = 4) 
    472       !!   
    473       !!    The sea level anomaly at the observation points is then computed  
    474       !!    by removing a mean dynamic topography (defined at the obs. point). 
     512      !! 
     513      !!    Two horizontal averaging schemes are also available: 
     514      !!        - weighted radial footprint        (k2dint = 5) 
     515      !!        - weighted rectangular footprint   (k2dint = 6) 
     516      !! 
    475517      !! 
    476518      !! ** Action  : 
     
    478520      !! History : 
    479521      !!      ! 07-03 (A. Weaver) 
     522      !!      ! 15-02 (M. Martin) Combined routine for surface types 
     523      !!      ! 17-03 (M. Martin) Added horizontal averaging options 
    480524      !!----------------------------------------------------------------------- 
    481    
     525 
    482526      !! * Modules used 
    483527      USE obs_surf_def  ! Definition of storage space for surface observations 
     
    486530 
    487531      !! * Arguments 
    488       TYPE(obs_surf), INTENT(INOUT) :: sladatqc     ! Subset of surface data not failing screening 
    489       INTEGER, INTENT(IN) :: kt      ! Time step 
    490       INTEGER, INTENT(IN) :: kpi     ! Model grid parameters 
     532      TYPE(obs_surf), INTENT(INOUT) :: & 
     533         & surfdataqc                  ! Subset of surface data passing QC 
     534      INTEGER, INTENT(IN) :: kt        ! Time step 
     535      INTEGER, INTENT(IN) :: kpi       ! Model grid parameters 
    491536      INTEGER, INTENT(IN) :: kpj 
    492       INTEGER, INTENT(IN) :: kit000   ! Number of the first time step  
    493                                       !   (kit000-1 = restart time) 
    494       INTEGER, INTENT(IN) :: k2dint   ! Horizontal interpolation type (see header) 
    495       REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 
    496          & psshn,  &    ! Model SSH field 
    497          & psshmask     ! Land-sea mask 
    498           
     537      INTEGER, INTENT(IN) :: kit000    ! Number of the first time step  
     538                                       !   (kit000-1 = restart time) 
     539      INTEGER, INTENT(IN) :: kdaystp   ! Number of time steps per day 
     540      INTEGER, INTENT(IN) :: k2dint    ! Horizontal interpolation type (see header) 
     541      REAL(wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 
     542         & psurf,  &                   ! Model surface field 
     543         & psurfmask                   ! Land-sea mask 
     544      LOGICAL, INTENT(IN) :: ldnightav ! Logical for averaging night-time data 
     545      REAL(KIND=wp), INTENT(IN) :: & 
     546         & plamscl, &                  ! Diameter in metres of obs footprint in E/W, N/S directions 
     547         & pphiscl                     ! This is the full width (rather than half-width) 
     548      LOGICAL, INTENT(IN) :: & 
     549         & lindegrees                  ! T=> plamscl and pphiscl are specified in degrees, F=> in metres 
     550 
    499551      !! * Local declarations 
    500552      INTEGER :: ji 
     
    502554      INTEGER :: jobs 
    503555      INTEGER :: inrc 
    504       INTEGER :: isla 
     556      INTEGER :: isurf 
    505557      INTEGER :: iobs 
    506       REAL(KIND=wp) :: zlam 
    507       REAL(KIND=wp) :: zphi 
    508       REAL(KIND=wp) :: zext(1), zobsmask(1) 
    509       REAL(kind=wp), DIMENSION(2,2,1) :: & 
    510          & zweig 
     558      INTEGER :: imaxifp, imaxjfp 
     559      INTEGER :: imodi, imodj 
     560      INTEGER :: idayend 
     561      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 
     562         & igrdi,   & 
     563         & igrdj,   & 
     564         & igrdip1, & 
     565         & igrdjp1 
     566      INTEGER, DIMENSION(:,:), SAVE, ALLOCATABLE :: & 
     567         & icount_night,      & 
     568         & imask_night 
     569      REAL(wp) :: zlam 
     570      REAL(wp) :: zphi 
     571      REAL(wp), DIMENSION(1) :: zext, zobsmask 
     572      REAL(wp) :: zdaystp 
    511573      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 
    512          & zmask, & 
    513          & zsshl, & 
    514          & zglam, & 
    515          & zgphi 
    516       INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 
    517          & igrdi, & 
    518          & igrdj 
     574         & zweig,  & 
     575         & zmask,  & 
     576         & zsurf,  & 
     577         & zsurfm, & 
     578         & zsurftmp, & 
     579         & zglam,  & 
     580         & zgphi,  & 
     581         & zglamf, & 
     582         & zgphif 
     583 
     584      REAL(wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: & 
     585         & zintmp,  & 
     586         & zouttmp, & 
     587         & zmeanday    ! to compute model sst in region of 24h daylight (pole) 
    519588 
    520589      !------------------------------------------------------------------------ 
    521590      ! Local initialization  
    522591      !------------------------------------------------------------------------ 
    523       ! ... Record and data counters 
     592      ! Record and data counters 
    524593      inrc = kt - kit000 + 2 
    525       isla = sladatqc%nsstp(inrc) 
     594      isurf = surfdataqc%nsstp(inrc) 
     595 
     596      ! Work out the maximum footprint size for the  
     597      ! interpolation/averaging in model grid-points - has to be even. 
     598 
     599      CALL obs_max_fpsize( k2dint, plamscl, pphiscl, lindegrees, psurfmask, imaxifp, imaxjfp ) 
     600 
     601 
     602      IF ( ldnightav ) THEN 
     603 
     604      ! Initialize array for night mean 
     605         IF ( kt == 0 ) THEN 
     606            ALLOCATE ( icount_night(kpi,kpj) ) 
     607            ALLOCATE ( imask_night(kpi,kpj) ) 
     608            ALLOCATE ( zintmp(kpi,kpj) ) 
     609            ALLOCATE ( zouttmp(kpi,kpj) ) 
     610            ALLOCATE ( zmeanday(kpi,kpj) ) 
     611            nday_qsr = -1   ! initialisation flag for nbc_dcy 
     612         ENDIF 
     613 
     614         ! Night-time means are calculated for night-time values over timesteps: 
     615         !  [1 <= kt <= kdaystp], [kdaystp+1 <= kt <= 2*kdaystp], ..... 
     616         idayend = MOD( kt - kit000 + 1, kdaystp ) 
     617 
     618         ! Initialize night-time mean for first timestep of the day 
     619         IF ( idayend == 1 .OR. kt == 0 ) THEN 
     620            DO jj = 1, jpj 
     621               DO ji = 1, jpi 
     622                  surfdataqc%vdmean(ji,jj) = 0.0 
     623                  zmeanday(ji,jj) = 0.0 
     624                  icount_night(ji,jj) = 0 
     625               END DO 
     626            END DO 
     627         ENDIF 
     628 
     629         zintmp(:,:) = 0.0 
     630         zouttmp(:,:) = sbc_dcy( zintmp(:,:), .TRUE. ) 
     631         imask_night(:,:) = INT( zouttmp(:,:) ) 
     632 
     633         DO jj = 1, jpj 
     634            DO ji = 1, jpi 
     635               ! Increment the temperature field for computing night mean and counter 
     636               surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj)  & 
     637                      &                    + psurf(ji,jj) * REAL( imask_night(ji,jj) ) 
     638               zmeanday(ji,jj)          = zmeanday(ji,jj) + psurf(ji,jj) 
     639               icount_night(ji,jj)      = icount_night(ji,jj) + imask_night(ji,jj) 
     640            END DO 
     641         END DO 
     642 
     643         ! Compute the night-time mean at the end of the day 
     644         zdaystp = 1.0 / REAL( kdaystp ) 
     645         IF ( idayend == 0 ) THEN 
     646            IF (lwp) WRITE(numout,*) 'Calculating surfdataqc%vdmean on time-step: ',kt 
     647            DO jj = 1, jpj 
     648               DO ji = 1, jpi 
     649                  ! Test if "no night" point 
     650                  IF ( icount_night(ji,jj) > 0 ) THEN 
     651                     surfdataqc%vdmean(ji,jj) = surfdataqc%vdmean(ji,jj) & 
     652                       &                        / REAL( icount_night(ji,jj) ) 
     653                  ELSE 
     654                     !At locations where there is no night (e.g. poles), 
     655                     ! calculate daily mean instead of night-time mean. 
     656                     surfdataqc%vdmean(ji,jj) = zmeanday(ji,jj) * zdaystp 
     657                  ENDIF 
     658               END DO 
     659            END DO 
     660         ENDIF 
     661 
     662      ENDIF 
    526663 
    527664      ! Get the data for interpolation 
    528665 
    529666      ALLOCATE( & 
    530          & igrdi(2,2,isla), & 
    531          & igrdj(2,2,isla), & 
    532          & zglam(2,2,isla), & 
    533          & zgphi(2,2,isla), & 
    534          & zmask(2,2,isla), & 
    535          & zsshl(2,2,isla)  & 
     667         & zweig(imaxifp,imaxjfp,1),      & 
     668         & igrdi(imaxifp,imaxjfp,isurf), & 
     669         & igrdj(imaxifp,imaxjfp,isurf), & 
     670         & zglam(imaxifp,imaxjfp,isurf), & 
     671         & zgphi(imaxifp,imaxjfp,isurf), & 
     672         & zmask(imaxifp,imaxjfp,isurf), & 
     673         & zsurf(imaxifp,imaxjfp,isurf), & 
     674         & zsurftmp(imaxifp,imaxjfp,isurf),  & 
     675         & zglamf(imaxifp+1,imaxjfp+1,isurf), & 
     676         & zgphif(imaxifp+1,imaxjfp+1,isurf), & 
     677         & igrdip1(imaxifp+1,imaxjfp+1,isurf), & 
     678         & igrdjp1(imaxifp+1,imaxjfp+1,isurf) & 
    536679         & ) 
    537        
    538       DO jobs = sladatqc%nsurfup + 1, sladatqc%nsurfup + isla 
    539          iobs = jobs - sladatqc%nsurfup 
    540          igrdi(1,1,iobs) = sladatqc%mi(jobs)-1 
    541          igrdj(1,1,iobs) = sladatqc%mj(jobs)-1 
    542          igrdi(1,2,iobs) = sladatqc%mi(jobs)-1 
    543          igrdj(1,2,iobs) = sladatqc%mj(jobs) 
    544          igrdi(2,1,iobs) = sladatqc%mi(jobs) 
    545          igrdj(2,1,iobs) = sladatqc%mj(jobs)-1 
    546          igrdi(2,2,iobs) = sladatqc%mi(jobs) 
    547          igrdj(2,2,iobs) = sladatqc%mj(jobs) 
     680 
     681      DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf 
     682         iobs = jobs - surfdataqc%nsurfup 
     683         DO ji = 0, imaxifp 
     684            imodi = surfdataqc%mi(jobs) - int(imaxifp/2) + ji - 1 
     685             
     686            !Deal with wrap around in longitude 
     687            IF ( imodi < 1      ) imodi = imodi + jpiglo 
     688            IF ( imodi > jpiglo ) imodi = imodi - jpiglo 
     689             
     690            DO jj = 0, imaxjfp 
     691               imodj = surfdataqc%mj(jobs) - int(imaxjfp/2) + jj - 1 
     692               !If model values are out of the domain to the north/south then 
     693               !set them to be the edge of the domain 
     694               IF ( imodj < 1      ) imodj = 1 
     695               IF ( imodj > jpjglo ) imodj = jpjglo 
     696 
     697               igrdip1(ji+1,jj+1,iobs) = imodi 
     698               igrdjp1(ji+1,jj+1,iobs) = imodj 
     699                
     700               IF ( ji >= 1 .AND. jj >= 1 ) THEN 
     701                  igrdi(ji,jj,iobs) = imodi 
     702                  igrdj(ji,jj,iobs) = imodj 
     703               ENDIF 
     704                
     705            END DO 
     706         END DO 
    548707      END DO 
    549708 
    550       CALL obs_int_comm_2d( 2, 2, isla, & 
     709      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 
    551710         &                  igrdi, igrdj, glamt, zglam ) 
    552       CALL obs_int_comm_2d( 2, 2, isla, & 
     711      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 
    553712         &                  igrdi, igrdj, gphit, zgphi ) 
    554       CALL obs_int_comm_2d( 2, 2, isla, & 
    555          &                  igrdi, igrdj, psshmask, zmask ) 
    556       CALL obs_int_comm_2d( 2, 2, isla, & 
    557          &                  igrdi, igrdj, psshn, zsshl ) 
     713      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 
     714         &                  igrdi, igrdj, psurfmask, zmask ) 
     715      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 
     716         &                  igrdi, igrdj, psurf, zsurf ) 
     717      CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, & 
     718         &                  igrdip1, igrdjp1, glamf, zglamf ) 
     719      CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, & 
     720         &                  igrdip1, igrdjp1, gphif, zgphif ) 
     721 
     722      ! At the end of the day get interpolated means 
     723      IF ( idayend == 0 .AND. ldnightav ) THEN 
     724 
     725         ALLOCATE( & 
     726            & zsurfm(imaxifp,imaxjfp,isurf)  & 
     727            & ) 
     728 
     729         CALL obs_int_comm_2d( imaxifp,imaxjfp, isurf, kpi, kpj, igrdi, igrdj, & 
     730            &               surfdataqc%vdmean(:,:), zsurfm ) 
     731 
     732      ENDIF 
    558733 
    559734      ! Loop over observations 
    560  
    561       DO jobs = sladatqc%nsurfup + 1, sladatqc%nsurfup + isla 
    562  
    563          iobs = jobs - sladatqc%nsurfup 
    564  
    565          IF ( kt /= sladatqc%mstp(jobs) ) THEN 
    566              
     735      DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf 
     736 
     737         iobs = jobs - surfdataqc%nsurfup 
     738 
     739         IF ( kt /= surfdataqc%mstp(jobs) ) THEN 
     740 
    567741            IF(lwp) THEN 
    568742               WRITE(numout,*) 
     
    574748               WRITE(numout,*) ' Record  = ', jobs,                & 
    575749                  &            ' kt      = ', kt,                  & 
    576                   &            ' mstp    = ', sladatqc%mstp(jobs), & 
    577                   &            ' ntyp    = ', sladatqc%ntyp(jobs) 
     750                  &            ' mstp    = ', surfdataqc%mstp(jobs), & 
     751                  &            ' ntyp    = ', surfdataqc%ntyp(jobs) 
    578752            ENDIF 
    579             CALL ctl_stop( 'obs_sla_opt', 'Inconsistent time' ) 
    580              
     753            CALL ctl_stop( 'obs_surf_opt', 'Inconsistent time' ) 
     754 
     755         ENDIF 
     756 
     757         zlam = surfdataqc%rlam(jobs) 
     758         zphi = surfdataqc%rphi(jobs) 
     759 
     760         IF ( ldnightav .AND. idayend == 0 ) THEN 
     761            ! Night-time averaged data 
     762            zsurftmp(:,:,iobs) = zsurfm(:,:,iobs) 
     763         ELSE 
     764            zsurftmp(:,:,iobs) = zsurf(:,:,iobs) 
     765         ENDIF 
     766 
     767         IF ( k2dint <= 4 ) THEN 
     768 
     769            ! Get weights to interpolate the model value to the observation point 
     770            CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         & 
     771               &                   zglam(:,:,iobs), zgphi(:,:,iobs), & 
     772               &                   zmask(:,:,iobs), zweig, zobsmask ) 
     773 
     774            ! Interpolate the model value to the observation point  
     775            CALL obs_int_h2d( 1, 1, zweig, zsurftmp(:,:,iobs), zext ) 
     776 
     777         ELSE 
     778 
     779            ! Get weights to average the model SLA to the observation footprint 
     780            CALL obs_avg_h2d_init( 1, 1, imaxifp, imaxjfp, k2dint, zlam,  zphi, & 
     781               &                   zglam(:,:,iobs), zgphi(:,:,iobs), & 
     782               &                   zglamf(:,:,iobs), zgphif(:,:,iobs), & 
     783               &                   zmask(:,:,iobs), plamscl, pphiscl, & 
     784               &                   lindegrees, zweig, zobsmask ) 
     785 
     786            ! Average the model SST to the observation footprint 
     787            CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, & 
     788               &              zweig, zsurftmp(:,:,iobs),  zext ) 
     789 
     790         ENDIF 
     791 
     792         IF ( TRIM(surfdataqc%cvars(1)) == 'SLA' .AND. surfdataqc%nextra == 2 ) THEN 
     793            ! ... Remove the MDT from the SSH at the observation point to get the SLA 
     794            surfdataqc%rext(jobs,1) = zext(1) 
     795            surfdataqc%rmod(jobs,1) = surfdataqc%rext(jobs,1) - surfdataqc%rext(jobs,2) 
     796         ELSE 
     797            surfdataqc%rmod(jobs,1) = zext(1) 
    581798         ENDIF 
    582799          
    583          zlam = sladatqc%rlam(jobs) 
    584          zphi = sladatqc%rphi(jobs) 
    585  
    586          ! Get weights to interpolate the model SSH to the observation point 
    587          CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         & 
    588             &                   zglam(:,:,iobs), zgphi(:,:,iobs), & 
    589             &                   zmask(:,:,iobs), zweig, zobsmask ) 
    590           
    591  
    592          ! Interpolate the model SSH to the observation point 
    593          CALL obs_int_h2d( 1, 1,      & 
    594             &              zweig, zsshl(:,:,iobs),  zext ) 
    595           
    596          sladatqc%rext(jobs,1) = zext(1) 
    597          ! ... Remove the MDT at the observation point 
    598          sladatqc%rmod(jobs,1) = sladatqc%rext(jobs,1) - sladatqc%rext(jobs,2) 
     800         IF ( zext(1) == obfillflt ) THEN 
     801            ! If the observation value is a fill value, set QC flag to bad 
     802            surfdataqc%nqc(jobs) = 4 
     803         ENDIF 
    599804 
    600805      END DO 
     
    602807      ! Deallocate the data for interpolation 
    603808      DEALLOCATE( & 
     809         & zweig, & 
    604810         & igrdi, & 
    605811         & igrdj, & 
     
    607813         & zgphi, & 
    608814         & zmask, & 
    609          & zsshl  & 
     815         & zsurf, & 
     816         & zsurftmp, & 
     817         & zglamf, & 
     818         & zgphif, & 
     819         & igrdip1,& 
     820         & igrdjp1 & 
    610821         & ) 
    611822 
    612       sladatqc%nsurfup = sladatqc%nsurfup + isla 
    613  
    614    END SUBROUTINE obs_sla_opt 
    615  
    616    SUBROUTINE obs_sst_opt( sstdatqc, kt, kpi, kpj, kit000, kdaystp, & 
    617       &                    psstn, psstmask, k2dint, ld_nightav ) 
    618       !!----------------------------------------------------------------------- 
    619       !! 
    620       !!                     ***  ROUTINE obs_sst_opt  *** 
    621       !! 
    622       !! ** Purpose : Compute the model counterpart of surface temperature 
    623       !!              data by interpolating from the model grid to the  
    624       !!              observation point. 
    625       !! 
    626       !! ** Method  : Linearly interpolate to each observation point using  
    627       !!              the model values at the corners of the surrounding grid box. 
    628       !! 
    629       !!    The now model SST is first computed at the obs (lon, lat) point. 
    630       !! 
    631       !!    Several horizontal interpolation schemes are available: 
    632       !!        - distance-weighted (great circle) (k2dint = 0) 
    633       !!        - distance-weighted (small angle)  (k2dint = 1) 
    634       !!        - bilinear (geographical grid)     (k2dint = 2) 
    635       !!        - bilinear (quadrilateral grid)    (k2dint = 3) 
    636       !!        - polynomial (quadrilateral grid)  (k2dint = 4) 
    637       !! 
    638       !! 
    639       !! ** Action  : 
    640       !! 
    641       !! History : 
    642       !!        !  07-07  (S. Ricci ) : Original 
    643       !!       
    644       !!----------------------------------------------------------------------- 
    645  
    646       !! * Modules used 
    647       USE obs_surf_def  ! Definition of storage space for surface observations 
    648       USE sbcdcy 
    649  
    650       IMPLICIT NONE 
    651  
    652       !! * Arguments 
    653       TYPE(obs_surf), INTENT(INOUT) :: & 
    654          & sstdatqc     ! Subset of surface data not failing screening 
    655       INTEGER, INTENT(IN) :: kt        ! Time step 
    656       INTEGER, INTENT(IN) :: kpi       ! Model grid parameters 
    657       INTEGER, INTENT(IN) :: kpj 
    658       INTEGER, INTENT(IN) :: kit000    ! Number of the first time step  
    659                                        !   (kit000-1 = restart time) 
    660       INTEGER, INTENT(IN) :: k2dint    ! Horizontal interpolation type (see header) 
    661       INTEGER, INTENT(IN) :: kdaystp   ! Number of time steps per day   
    662       REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 
    663          & psstn,  &    ! Model SST field 
    664          & psstmask     ! Land-sea mask 
    665  
    666       !! * Local declarations 
    667       INTEGER :: ji 
    668       INTEGER :: jj 
    669       INTEGER :: jobs 
    670       INTEGER :: inrc 
    671       INTEGER :: isst 
    672       INTEGER :: iobs 
    673       INTEGER :: idayend 
    674       REAL(KIND=wp) :: zlam 
    675       REAL(KIND=wp) :: zphi 
    676       REAL(KIND=wp) :: zext(1), zobsmask(1) 
    677       REAL(KIND=wp) :: zdaystp 
    678       INTEGER, DIMENSION(:,:), SAVE, ALLOCATABLE :: & 
    679          & icount_sstnight,      & 
    680          & imask_night 
    681       REAL(kind=wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: & 
    682          & zintmp, & 
    683          & zouttmp, &  
    684          & zmeanday    ! to compute model sst in region of 24h daylight (pole) 
    685       REAL(kind=wp), DIMENSION(2,2,1) :: & 
    686          & zweig 
    687       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 
    688          & zmask, & 
    689          & zsstl, & 
    690          & zsstm, & 
    691          & zglam, & 
    692          & zgphi 
    693       INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 
    694          & igrdi, & 
    695          & igrdj 
    696       LOGICAL, INTENT(IN) :: ld_nightav 
    697  
    698       !----------------------------------------------------------------------- 
    699       ! Local initialization  
    700       !----------------------------------------------------------------------- 
    701       ! ... Record and data counters 
    702       inrc = kt - kit000 + 2 
    703       isst = sstdatqc%nsstp(inrc) 
    704  
    705       IF ( ld_nightav ) THEN 
    706  
    707       ! Initialize array for night mean 
    708  
    709       IF ( kt .EQ. 0 ) THEN 
    710          ALLOCATE ( icount_sstnight(kpi,kpj) ) 
    711          ALLOCATE ( imask_night(kpi,kpj) ) 
    712          ALLOCATE ( zintmp(kpi,kpj) ) 
    713          ALLOCATE ( zouttmp(kpi,kpj) ) 
    714          ALLOCATE ( zmeanday(kpi,kpj) ) 
    715          nday_qsr = -1   ! initialisation flag for nbc_dcy 
    716       ENDIF 
    717  
    718       ! Initialize daily mean for first timestep 
    719       idayend = MOD( kt - kit000 + 1, kdaystp ) 
    720  
    721       ! Added kt == 0 test to catch restart case  
    722       IF ( idayend == 1 .OR. kt == 0) THEN 
    723          IF (lwp) WRITE(numout,*) 'Reset sstdatqc%vdmean on time-step: ',kt 
    724          DO jj = 1, jpj 
    725             DO ji = 1, jpi 
    726                sstdatqc%vdmean(ji,jj) = 0.0 
    727                zmeanday(ji,jj) = 0.0 
    728                icount_sstnight(ji,jj) = 0 
    729             END DO 
    730          END DO 
    731       ENDIF 
    732  
    733       zintmp(:,:) = 0.0 
    734       zouttmp(:,:) = sbc_dcy( zintmp(:,:), .TRUE. ) 
    735       imask_night(:,:) = INT( zouttmp(:,:) ) 
    736  
    737       DO jj = 1, jpj 
    738          DO ji = 1, jpi 
    739             ! Increment the temperature field for computing night mean and counter 
    740             sstdatqc%vdmean(ji,jj) = sstdatqc%vdmean(ji,jj)  & 
    741                    &                        + psstn(ji,jj)*imask_night(ji,jj) 
    742             zmeanday(ji,jj)        = zmeanday(ji,jj) + psstn(ji,jj) 
    743             icount_sstnight(ji,jj) = icount_sstnight(ji,jj) + imask_night(ji,jj) 
    744          END DO 
    745       END DO 
    746     
    747       ! Compute the daily mean at the end of day 
    748  
    749       zdaystp = 1.0 / REAL( kdaystp ) 
    750  
    751       IF ( idayend == 0 ) THEN  
    752          DO jj = 1, jpj 
    753             DO ji = 1, jpi 
    754                ! Test if "no night" point 
    755                IF ( icount_sstnight(ji,jj) .NE. 0 ) THEN 
    756                   sstdatqc%vdmean(ji,jj) = sstdatqc%vdmean(ji,jj) & 
    757                     &                        / icount_sstnight(ji,jj)  
    758                ELSE 
    759                   sstdatqc%vdmean(ji,jj) = zmeanday(ji,jj) * zdaystp 
    760                ENDIF 
    761             END DO 
    762          END DO 
    763       ENDIF 
    764  
    765       ENDIF 
    766  
    767       ! Get the data for interpolation 
    768        
    769       ALLOCATE( & 
    770          & igrdi(2,2,isst), & 
    771          & igrdj(2,2,isst), & 
    772          & zglam(2,2,isst), & 
    773          & zgphi(2,2,isst), & 
    774          & zmask(2,2,isst), & 
    775          & zsstl(2,2,isst)  & 
    776          & ) 
    777        
    778       DO jobs = sstdatqc%nsurfup + 1, sstdatqc%nsurfup + isst 
    779          iobs = jobs - sstdatqc%nsurfup 
    780          igrdi(1,1,iobs) = sstdatqc%mi(jobs)-1 
    781          igrdj(1,1,iobs) = sstdatqc%mj(jobs)-1 
    782          igrdi(1,2,iobs) = sstdatqc%mi(jobs)-1 
    783          igrdj(1,2,iobs) = sstdatqc%mj(jobs) 
    784          igrdi(2,1,iobs) = sstdatqc%mi(jobs) 
    785          igrdj(2,1,iobs) = sstdatqc%mj(jobs)-1 
    786          igrdi(2,2,iobs) = sstdatqc%mi(jobs) 
    787          igrdj(2,2,iobs) = sstdatqc%mj(jobs) 
    788       END DO 
    789        
    790       CALL obs_int_comm_2d( 2, 2, isst, & 
    791          &                  igrdi, igrdj, glamt, zglam ) 
    792       CALL obs_int_comm_2d( 2, 2, isst, & 
    793          &                  igrdi, igrdj, gphit, zgphi ) 
    794       CALL obs_int_comm_2d( 2, 2, isst, & 
    795          &                  igrdi, igrdj, psstmask, zmask ) 
    796       CALL obs_int_comm_2d( 2, 2, isst, & 
    797          &                  igrdi, igrdj, psstn, zsstl ) 
    798  
    799       ! At the end of the day get interpolated means 
    800       IF ( idayend == 0 .AND. ld_nightav ) THEN 
    801  
    802          ALLOCATE( & 
    803             & zsstm(2,2,isst)  & 
    804             & ) 
    805  
    806          CALL obs_int_comm_2d( 2, 2, isst, igrdi, igrdj, & 
    807             &               sstdatqc%vdmean(:,:), zsstm ) 
    808  
    809       ENDIF 
    810  
    811       ! Loop over observations 
    812  
    813       DO jobs = sstdatqc%nsurfup + 1, sstdatqc%nsurfup + isst 
    814           
    815          iobs = jobs - sstdatqc%nsurfup 
    816           
    817          IF ( kt /= sstdatqc%mstp(jobs) ) THEN 
    818              
    819             IF(lwp) THEN 
    820                WRITE(numout,*) 
    821                WRITE(numout,*) ' E R R O R : Observation',              & 
    822                   &            ' time step is not consistent with the', & 
    823                   &            ' model time step' 
    824                WRITE(numout,*) ' =========' 
    825                WRITE(numout,*) 
    826                WRITE(numout,*) ' Record  = ', jobs,                & 
    827                   &            ' kt      = ', kt,                  & 
    828                   &            ' mstp    = ', sstdatqc%mstp(jobs), & 
    829                   &            ' ntyp    = ', sstdatqc%ntyp(jobs) 
    830             ENDIF 
    831             CALL ctl_stop( 'obs_sst_opt', 'Inconsistent time' ) 
    832              
    833          ENDIF 
    834           
    835          zlam = sstdatqc%rlam(jobs) 
    836          zphi = sstdatqc%rphi(jobs) 
    837           
    838          ! Get weights to interpolate the model SST 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 SST to the observation point  
    844  
    845          IF ( ld_nightav ) THEN 
    846  
    847            IF ( idayend == 0 )  THEN 
    848                ! Daily averaged/diurnal cycle of SST  data 
    849                CALL obs_int_h2d( 1, 1,      &  
    850                      &              zweig, zsstm(:,:,iobs), zext ) 
    851             ELSE  
    852                CALL ctl_stop( ' ld_nightav is set to true: a nonzero' //     & 
    853                      &           ' number of night SST data should' // & 
    854                      &           ' only occur at the end of a given day' ) 
    855             ENDIF 
    856  
    857          ELSE 
    858  
    859             CALL obs_int_h2d( 1, 1,      & 
    860             &              zweig, zsstl(:,:,iobs),  zext ) 
    861  
    862          ENDIF 
    863          sstdatqc%rmod(jobs,1) = zext(1) 
    864           
    865       END DO 
    866        
    867       ! Deallocate the data for interpolation 
    868       DEALLOCATE( & 
    869          & igrdi, & 
    870          & igrdj, & 
    871          & zglam, & 
    872          & zgphi, & 
    873          & zmask, & 
    874          & zsstl  & 
    875          & ) 
    876  
    877       ! At the end of the day also get interpolated means 
    878       IF ( idayend == 0 .AND. ld_nightav ) THEN 
     823      ! At the end of the day also deallocate night-time mean array 
     824      IF ( idayend == 0 .AND. ldnightav ) THEN 
    879825         DEALLOCATE( & 
    880             & zsstm  & 
     826            & zsurfm  & 
    881827            & ) 
    882828      ENDIF 
    883        
    884       sstdatqc%nsurfup = sstdatqc%nsurfup + isst 
    885  
    886    END SUBROUTINE obs_sst_opt 
    887  
    888    SUBROUTINE obs_sss_opt 
    889       !!----------------------------------------------------------------------- 
    890       !! 
    891       !!                     ***  ROUTINE obs_sss_opt  *** 
    892       !! 
    893       !! ** Purpose : Compute the model counterpart of sea surface salinity 
    894       !!              data by interpolating from the model grid to the  
    895       !!              observation point. 
    896       !! 
    897       !! ** Method  :  
    898       !! 
    899       !! ** Action  : 
    900       !! 
    901       !! History : 
    902       !!      ! ??-??  
    903       !!----------------------------------------------------------------------- 
    904  
    905       IMPLICIT NONE 
    906  
    907    END SUBROUTINE obs_sss_opt 
    908  
    909    SUBROUTINE obs_seaice_opt( seaicedatqc, kt, kpi, kpj, kit000, & 
    910       &                    pseaicen, pseaicemask, k2dint ) 
    911  
    912       !!----------------------------------------------------------------------- 
    913       !! 
    914       !!                     ***  ROUTINE obs_seaice_opt  *** 
    915       !! 
    916       !! ** Purpose : Compute the model counterpart of surface temperature 
    917       !!              data by interpolating from the model grid to the  
    918       !!              observation point. 
    919       !! 
    920       !! ** Method  : Linearly interpolate to each observation point using  
    921       !!              the model values at the corners of the surrounding grid box. 
    922       !! 
    923       !!    The now model sea ice is first computed at the obs (lon, lat) point. 
    924       !! 
    925       !!    Several horizontal interpolation schemes are available: 
    926       !!        - distance-weighted (great circle) (k2dint = 0) 
    927       !!        - distance-weighted (small angle)  (k2dint = 1) 
    928       !!        - bilinear (geographical grid)     (k2dint = 2) 
    929       !!        - bilinear (quadrilateral grid)    (k2dint = 3) 
    930       !!        - polynomial (quadrilateral grid)  (k2dint = 4) 
    931       !! 
    932       !! 
    933       !! ** Action  : 
    934       !! 
    935       !! History : 
    936       !!        !  07-07  (S. Ricci ) : Original 
    937       !!       
    938       !!----------------------------------------------------------------------- 
    939  
    940       !! * Modules used 
    941       USE obs_surf_def  ! Definition of storage space for surface observations 
    942  
    943       IMPLICIT NONE 
    944  
    945       !! * Arguments 
    946       TYPE(obs_surf), INTENT(INOUT) :: seaicedatqc     ! Subset of surface data not failing screening 
    947       INTEGER, INTENT(IN) :: kt       ! Time step 
    948       INTEGER, INTENT(IN) :: kpi      ! Model grid parameters 
    949       INTEGER, INTENT(IN) :: kpj 
    950       INTEGER, INTENT(IN) :: kit000   ! Number of the first time step  
    951                                       !   (kit000-1 = restart time) 
    952       INTEGER, INTENT(IN) :: k2dint   ! Horizontal interpolation type (see header) 
    953       REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 
    954          & pseaicen,  &    ! Model sea ice field 
    955          & pseaicemask     ! Land-sea mask 
    956           
    957       !! * Local declarations 
    958       INTEGER :: ji 
    959       INTEGER :: jj 
    960       INTEGER :: jobs 
    961       INTEGER :: inrc 
    962       INTEGER :: iseaice 
    963       INTEGER :: iobs 
    964         
    965       REAL(KIND=wp) :: zlam 
    966       REAL(KIND=wp) :: zphi 
    967       REAL(KIND=wp) :: zext(1), zobsmask(1) 
    968       REAL(kind=wp), DIMENSION(2,2,1) :: & 
    969          & zweig 
    970       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 
    971          & zmask, & 
    972          & zseaicel, & 
    973          & zglam, & 
    974          & zgphi 
    975       INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 
    976          & igrdi, & 
    977          & igrdj 
    978  
    979       !------------------------------------------------------------------------ 
    980       ! Local initialization  
    981       !------------------------------------------------------------------------ 
    982       ! ... Record and data counters 
    983       inrc = kt - kit000 + 2 
    984       iseaice = seaicedatqc%nsstp(inrc) 
    985  
    986       ! Get the data for interpolation 
    987        
    988       ALLOCATE( & 
    989          & igrdi(2,2,iseaice), & 
    990          & igrdj(2,2,iseaice), & 
    991          & zglam(2,2,iseaice), & 
    992          & zgphi(2,2,iseaice), & 
    993          & zmask(2,2,iseaice), & 
    994          & zseaicel(2,2,iseaice)  & 
    995          & ) 
    996        
    997       DO jobs = seaicedatqc%nsurfup + 1, seaicedatqc%nsurfup + iseaice 
    998          iobs = jobs - seaicedatqc%nsurfup 
    999          igrdi(1,1,iobs) = seaicedatqc%mi(jobs)-1 
    1000          igrdj(1,1,iobs) = seaicedatqc%mj(jobs)-1 
    1001          igrdi(1,2,iobs) = seaicedatqc%mi(jobs)-1 
    1002          igrdj(1,2,iobs) = seaicedatqc%mj(jobs) 
    1003          igrdi(2,1,iobs) = seaicedatqc%mi(jobs) 
    1004          igrdj(2,1,iobs) = seaicedatqc%mj(jobs)-1 
    1005          igrdi(2,2,iobs) = seaicedatqc%mi(jobs) 
    1006          igrdj(2,2,iobs) = seaicedatqc%mj(jobs) 
    1007       END DO 
    1008        
    1009       CALL obs_int_comm_2d( 2, 2, iseaice, & 
    1010          &                  igrdi, igrdj, glamt, zglam ) 
    1011       CALL obs_int_comm_2d( 2, 2, iseaice, & 
    1012          &                  igrdi, igrdj, gphit, zgphi ) 
    1013       CALL obs_int_comm_2d( 2, 2, iseaice, & 
    1014          &                  igrdi, igrdj, pseaicemask, zmask ) 
    1015       CALL obs_int_comm_2d( 2, 2, iseaice, & 
    1016          &                  igrdi, igrdj, pseaicen, zseaicel ) 
    1017        
    1018       DO jobs = seaicedatqc%nsurfup + 1, seaicedatqc%nsurfup + iseaice 
    1019           
    1020          iobs = jobs - seaicedatqc%nsurfup 
    1021           
    1022          IF ( kt /= seaicedatqc%mstp(jobs) ) THEN 
    1023              
    1024             IF(lwp) THEN 
    1025                WRITE(numout,*) 
    1026                WRITE(numout,*) ' E R R O R : Observation',              & 
    1027                   &            ' time step is not consistent with the', & 
    1028                   &            ' model time step' 
    1029                WRITE(numout,*) ' =========' 
    1030                WRITE(numout,*) 
    1031                WRITE(numout,*) ' Record  = ', jobs,                & 
    1032                   &            ' kt      = ', kt,                  & 
    1033                   &            ' mstp    = ', seaicedatqc%mstp(jobs), & 
    1034                   &            ' ntyp    = ', seaicedatqc%ntyp(jobs) 
    1035             ENDIF 
    1036             CALL ctl_stop( 'obs_seaice_opt', 'Inconsistent time' ) 
    1037              
    1038          ENDIF 
    1039           
    1040          zlam = seaicedatqc%rlam(jobs) 
    1041          zphi = seaicedatqc%rphi(jobs) 
    1042           
    1043          ! Get weights to interpolate the model sea ice to the observation point 
    1044          CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         & 
    1045             &                   zglam(:,:,iobs), zgphi(:,:,iobs), & 
    1046             &                   zmask(:,:,iobs), zweig, zobsmask ) 
    1047           
    1048          ! ... Interpolate the model sea ice to the observation point 
    1049          CALL obs_int_h2d( 1, 1,      & 
    1050             &              zweig, zseaicel(:,:,iobs),  zext ) 
    1051           
    1052          seaicedatqc%rmod(jobs,1) = zext(1) 
    1053           
    1054       END DO 
    1055        
    1056       ! Deallocate the data for interpolation 
    1057       DEALLOCATE( & 
    1058          & igrdi,    & 
    1059          & igrdj,    & 
    1060          & zglam,    & 
    1061          & zgphi,    & 
    1062          & zmask,    & 
    1063          & zseaicel  & 
    1064          & ) 
    1065        
    1066       seaicedatqc%nsurfup = seaicedatqc%nsurfup + iseaice 
    1067  
    1068    END SUBROUTINE obs_seaice_opt 
    1069  
    1070    SUBROUTINE obs_vel_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, & 
    1071       &                    pun, pvn, pgdept, pumask, pvmask, k1dint, k2dint, & 
    1072       &                    ld_dailyav ) 
    1073       !!----------------------------------------------------------------------- 
    1074       !! 
    1075       !!                     ***  ROUTINE obs_vel_opt  *** 
    1076       !! 
    1077       !! ** Purpose : Compute the model counterpart of velocity profile 
    1078       !!              data by interpolating from the model grid to the  
    1079       !!              observation point. 
    1080       !! 
    1081       !! ** Method  : Linearly interpolate zonal and meridional components of velocity  
    1082       !!              to each observation point using the model values at the corners of  
    1083       !!              the surrounding grid box. The model velocity components are on a  
    1084       !!              staggered C- grid. 
    1085       !! 
    1086       !!    For velocity data from the TAO array, the model equivalent is 
    1087       !!    a daily mean velocity field. So, we first compute 
    1088       !!    the mean, then interpolate only at the end of the day. 
    1089       !! 
    1090       !! ** Action  : 
    1091       !! 
    1092       !! History : 
    1093       !!    ! 07-03 (K. Mogensen)      : Temperature and Salinity profiles 
    1094       !!    ! 08-10 (Maria Valdivieso) : Velocity component (U,V) profiles 
    1095       !!----------------------------------------------------------------------- 
    1096      
    1097       !! * Modules used 
    1098       USE obs_profiles_def ! Definition of storage space for profile obs. 
    1099  
    1100       IMPLICIT NONE 
    1101  
    1102       !! * Arguments 
    1103       TYPE(obs_prof), INTENT(INOUT) :: & 
    1104          & prodatqc        ! Subset of profile data not failing screening 
    1105       INTEGER, INTENT(IN) :: kt        ! Time step 
    1106       INTEGER, INTENT(IN) :: kpi       ! Model grid parameters 
    1107       INTEGER, INTENT(IN) :: kpj 
    1108       INTEGER, INTENT(IN) :: kpk  
    1109       INTEGER, INTENT(IN) :: kit000    ! Number of the first time step  
    1110                                        !   (kit000-1 = restart time) 
    1111       INTEGER, INTENT(IN) :: k1dint    ! Vertical interpolation type (see header) 
    1112       INTEGER, INTENT(IN) :: k2dint    ! Horizontal interpolation type (see header) 
    1113       INTEGER, INTENT(IN) :: kdaystp   ! Number of time steps per day                     
    1114       REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 
    1115          & pun,    &    ! Model zonal component of velocity 
    1116          & pvn,    &    ! Model meridional component of velocity 
    1117          & pumask, &    ! Land-sea mask 
    1118          & pvmask       ! Land-sea mask 
    1119       REAL(KIND=wp), INTENT(IN), DIMENSION(kpk) :: & 
    1120          & pgdept       ! Model array of depth levels 
    1121       LOGICAL, INTENT(IN) :: ld_dailyav 
    1122           
    1123       !! * Local declarations 
    1124       INTEGER :: ji 
    1125       INTEGER :: jj 
    1126       INTEGER :: jk 
    1127       INTEGER :: jobs 
    1128       INTEGER :: inrc 
    1129       INTEGER :: ipro 
    1130       INTEGER :: idayend 
    1131       INTEGER :: ista 
    1132       INTEGER :: iend 
    1133       INTEGER :: iobs 
    1134       INTEGER, DIMENSION(imaxavtypes) :: & 
    1135          & idailyavtypes 
    1136       REAL(KIND=wp) :: zlam 
    1137       REAL(KIND=wp) :: zphi 
    1138       REAL(KIND=wp) :: zdaystp 
    1139       REAL(KIND=wp), DIMENSION(kpk) :: & 
    1140          & zobsmasku, & 
    1141          & zobsmaskv, & 
    1142          & zobsmask,  & 
    1143          & zobsk,     & 
    1144          & zobs2k 
    1145       REAL(KIND=wp), DIMENSION(2,2,kpk) :: & 
    1146          & zweigu,zweigv 
    1147       REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: & 
    1148          & zumask, zvmask, & 
    1149          & zintu, & 
    1150          & zintv, & 
    1151          & zinmu, & 
    1152          & zinmv 
    1153       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 
    1154          & zglamu, zglamv, & 
    1155          & zgphiu, zgphiv 
    1156       INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 
    1157          & igrdiu, & 
    1158          & igrdju, & 
    1159          & igrdiv, & 
    1160          & igrdjv 
    1161  
    1162       !------------------------------------------------------------------------ 
    1163       ! Local initialization  
    1164       !------------------------------------------------------------------------ 
    1165       ! ... Record and data counters 
    1166       inrc = kt - kit000 + 2 
    1167       ipro = prodatqc%npstp(inrc) 
    1168  
    1169       ! Initialize daily mean for first timestep 
    1170       idayend = MOD( kt - kit000 + 1, kdaystp ) 
    1171  
    1172       ! Added kt == 0 test to catch restart case  
    1173       IF ( idayend == 1 .OR. kt == 0) THEN 
    1174          IF (lwp) WRITE(numout,*) 'Reset prodatqc%vdmean on time-step: ',kt 
    1175          prodatqc%vdmean(:,:,:,1) = 0.0 
    1176          prodatqc%vdmean(:,:,:,2) = 0.0 
    1177       ENDIF 
    1178  
    1179       ! Increment the zonal velocity field for computing daily mean 
    1180       prodatqc%vdmean(:,:,:,1) = prodatqc%vdmean(:,:,:,1) + pun(:,:,:) 
    1181       ! Increment the meridional velocity field for computing daily mean 
    1182       prodatqc%vdmean(:,:,:,2) = prodatqc%vdmean(:,:,:,2) + pvn(:,:,:) 
    1183     
    1184       ! Compute the daily mean at the end of day 
    1185       zdaystp = 1.0 / REAL( kdaystp ) 
    1186       IF ( idayend == 0 ) THEN 
    1187          prodatqc%vdmean(:,:,:,1) = prodatqc%vdmean(:,:,:,1) * zdaystp 
    1188          prodatqc%vdmean(:,:,:,2) = prodatqc%vdmean(:,:,:,2) * zdaystp 
    1189       ENDIF 
    1190  
    1191       ! Get the data for interpolation 
    1192       ALLOCATE( & 
    1193          & igrdiu(2,2,ipro),      & 
    1194          & igrdju(2,2,ipro),      & 
    1195          & igrdiv(2,2,ipro),      & 
    1196          & igrdjv(2,2,ipro),      & 
    1197          & zglamu(2,2,ipro), zglamv(2,2,ipro), & 
    1198          & zgphiu(2,2,ipro), zgphiv(2,2,ipro), & 
    1199          & zumask(2,2,kpk,ipro), zvmask(2,2,kpk,ipro), & 
    1200          & zintu(2,2,kpk,ipro),  & 
    1201          & zintv(2,2,kpk,ipro)   & 
    1202          & ) 
    1203  
    1204       DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro 
    1205          iobs = jobs - prodatqc%nprofup 
    1206          igrdiu(1,1,iobs) = prodatqc%mi(jobs,1)-1 
    1207          igrdju(1,1,iobs) = prodatqc%mj(jobs,1)-1 
    1208          igrdiu(1,2,iobs) = prodatqc%mi(jobs,1)-1 
    1209          igrdju(1,2,iobs) = prodatqc%mj(jobs,1) 
    1210          igrdiu(2,1,iobs) = prodatqc%mi(jobs,1) 
    1211          igrdju(2,1,iobs) = prodatqc%mj(jobs,1)-1 
    1212          igrdiu(2,2,iobs) = prodatqc%mi(jobs,1) 
    1213          igrdju(2,2,iobs) = prodatqc%mj(jobs,1) 
    1214          igrdiv(1,1,iobs) = prodatqc%mi(jobs,2)-1 
    1215          igrdjv(1,1,iobs) = prodatqc%mj(jobs,2)-1 
    1216          igrdiv(1,2,iobs) = prodatqc%mi(jobs,2)-1 
    1217          igrdjv(1,2,iobs) = prodatqc%mj(jobs,2) 
    1218          igrdiv(2,1,iobs) = prodatqc%mi(jobs,2) 
    1219          igrdjv(2,1,iobs) = prodatqc%mj(jobs,2)-1 
    1220          igrdiv(2,2,iobs) = prodatqc%mi(jobs,2) 
    1221          igrdjv(2,2,iobs) = prodatqc%mj(jobs,2) 
    1222       END DO 
    1223  
    1224       CALL obs_int_comm_2d( 2, 2, ipro, igrdiu, igrdju, glamu, zglamu ) 
    1225       CALL obs_int_comm_2d( 2, 2, ipro, igrdiu, igrdju, gphiu, zgphiu ) 
    1226       CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, pumask, zumask ) 
    1227       CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, pun, zintu ) 
    1228  
    1229       CALL obs_int_comm_2d( 2, 2, ipro, igrdiv, igrdjv, glamv, zglamv ) 
    1230       CALL obs_int_comm_2d( 2, 2, ipro, igrdiv, igrdjv, gphiv, zgphiv ) 
    1231       CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, pvmask, zvmask ) 
    1232       CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, pvn, zintv ) 
    1233  
    1234       ! At the end of the day also get interpolated means 
    1235       IF ( idayend == 0 ) THEN 
    1236  
    1237          ALLOCATE( & 
    1238             & zinmu(2,2,kpk,ipro),  & 
    1239             & zinmv(2,2,kpk,ipro)   & 
    1240             & ) 
    1241  
    1242          CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiu, igrdju, & 
    1243             &                  prodatqc%vdmean(:,:,:,1), zinmu ) 
    1244          CALL obs_int_comm_3d( 2, 2, ipro, kpk, igrdiv, igrdjv, & 
    1245             &                  prodatqc%vdmean(:,:,:,2), zinmv ) 
    1246  
    1247       ENDIF 
    1248  
    1249 ! loop over observations 
    1250  
    1251       DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro 
    1252  
    1253          iobs = jobs - prodatqc%nprofup 
    1254  
    1255          IF ( kt /= prodatqc%mstp(jobs) ) THEN 
    1256              
    1257             IF(lwp) THEN 
    1258                WRITE(numout,*) 
    1259                WRITE(numout,*) ' E R R O R : Observation',              & 
    1260                   &            ' time step is not consistent with the', & 
    1261                   &            ' model time step' 
    1262                WRITE(numout,*) ' =========' 
    1263                WRITE(numout,*) 
    1264                WRITE(numout,*) ' Record  = ', jobs,                    & 
    1265                   &            ' kt      = ', kt,                      & 
    1266                   &            ' mstp    = ', prodatqc%mstp(jobs), & 
    1267                   &            ' ntyp    = ', prodatqc%ntyp(jobs) 
    1268             ENDIF 
    1269             CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' ) 
    1270          ENDIF 
    1271           
    1272          zlam = prodatqc%rlam(jobs) 
    1273          zphi = prodatqc%rphi(jobs) 
    1274  
    1275          ! Initialize observation masks 
    1276  
    1277          zobsmasku(:) = 0.0 
    1278          zobsmaskv(:) = 0.0 
    1279           
    1280          ! Horizontal weights and vertical mask 
    1281  
    1282          IF  ( prodatqc%npvend(jobs,1) > 0 ) THEN 
    1283  
    1284             CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi,     & 
    1285                &                   zglamu(:,:,iobs), zgphiu(:,:,iobs), & 
    1286                &                   zumask(:,:,:,iobs), zweigu, zobsmasku ) 
    1287  
    1288          ENDIF 
    1289  
    1290           
    1291          IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 
    1292  
    1293             CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi,     & 
    1294                &                   zglamv(:,:,iobs), zgphiv(:,:,iobs), & 
    1295                &                   zvmask(:,:,:,iobs), zweigv, zobsmasku ) 
    1296  
    1297          ENDIF 
    1298  
    1299          ! Ensure that the vertical mask on u and v are consistent. 
    1300  
    1301          zobsmask(:) = MIN( zobsmasku(:), zobsmaskv(:) ) 
    1302  
    1303          IF ( prodatqc%npvend(jobs,1) > 0 ) THEN 
    1304  
    1305             zobsk(:) = obfillflt 
    1306  
    1307        IF ( ld_dailyav ) THEN 
    1308  
    1309                IF ( idayend == 0 )  THEN 
    1310                    
    1311                   ! Daily averaged data 
    1312                    
    1313                   CALL obs_int_h2d( kpk, kpk,      & 
    1314                      &              zweigu, zinmu(:,:,:,iobs), zobsk ) 
    1315                    
    1316                    
    1317                ELSE 
    1318                 
    1319                   CALL ctl_stop( ' A nonzero' //     & 
    1320                      &           ' number of U profile data should' // & 
    1321                      &           ' only occur at the end of a given day' ) 
    1322  
    1323                ENDIF 
    1324            
    1325             ELSE  
    1326                 
    1327                ! Point data 
    1328  
    1329                CALL obs_int_h2d( kpk, kpk,      & 
    1330                   &              zweigu, zintu(:,:,:,iobs), zobsk ) 
    1331  
    1332             ENDIF 
    1333  
    1334             !------------------------------------------------------------- 
    1335             ! Compute vertical second-derivative of the interpolating  
    1336             ! polynomial at obs points 
    1337             !------------------------------------------------------------- 
    1338              
    1339             IF ( k1dint == 1 ) THEN 
    1340                CALL obs_int_z1d_spl( kpk, zobsk, zobs2k,   & 
    1341                   &                  pgdept, zobsmask ) 
    1342             ENDIF 
    1343              
    1344             !----------------------------------------------------------------- 
    1345             !  Vertical interpolation to the observation point 
    1346             !----------------------------------------------------------------- 
    1347             ista = prodatqc%npvsta(jobs,1) 
    1348             iend = prodatqc%npvend(jobs,1) 
    1349             CALL obs_int_z1d( kpk,                & 
    1350                & prodatqc%var(1)%mvk(ista:iend),  & 
    1351                & k1dint, iend - ista + 1,         & 
    1352                & prodatqc%var(1)%vdep(ista:iend), & 
    1353                & zobsk, zobs2k,                   & 
    1354                & prodatqc%var(1)%vmod(ista:iend), & 
    1355                & pgdept, zobsmask ) 
    1356  
    1357          ENDIF 
    1358  
    1359          IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 
    1360  
    1361             zobsk(:) = obfillflt 
    1362  
    1363             IF ( ld_dailyav ) THEN 
    1364  
    1365                IF ( idayend == 0 )  THEN 
    1366  
    1367                   ! Daily averaged data 
    1368                    
    1369                   CALL obs_int_h2d( kpk, kpk,      & 
    1370                      &              zweigv, zinmv(:,:,:,iobs), zobsk ) 
    1371                    
    1372                ELSE 
    1373  
    1374                   CALL ctl_stop( ' A nonzero' //     & 
    1375                      &           ' number of V profile data should' // & 
    1376                      &           ' only occur at the end of a given day' ) 
    1377  
    1378                ENDIF 
    1379  
    1380             ELSE 
    1381                 
    1382                ! Point data 
    1383  
    1384                CALL obs_int_h2d( kpk, kpk,      & 
    1385                   &              zweigv, zintv(:,:,:,iobs), zobsk ) 
    1386  
    1387             ENDIF 
    1388  
    1389  
    1390             !------------------------------------------------------------- 
    1391             ! Compute vertical second-derivative of the interpolating  
    1392             ! polynomial at obs points 
    1393             !------------------------------------------------------------- 
    1394              
    1395             IF ( k1dint == 1 ) THEN 
    1396                CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, & 
    1397                   &                  pgdept, zobsmask ) 
    1398             ENDIF 
    1399              
    1400             !---------------------------------------------------------------- 
    1401             !  Vertical interpolation to the observation point 
    1402             !---------------------------------------------------------------- 
    1403             ista = prodatqc%npvsta(jobs,2) 
    1404             iend = prodatqc%npvend(jobs,2) 
    1405             CALL obs_int_z1d( kpk, & 
    1406                & prodatqc%var(2)%mvk(ista:iend),& 
    1407                & k1dint, iend - ista + 1, & 
    1408                & prodatqc%var(2)%vdep(ista:iend),& 
    1409                & zobsk, zobs2k, & 
    1410                & prodatqc%var(2)%vmod(ista:iend),& 
    1411                & pgdept, zobsmask ) 
    1412  
    1413          ENDIF 
    1414  
    1415       END DO 
    1416   
    1417       ! Deallocate the data for interpolation 
    1418       DEALLOCATE( & 
    1419          & igrdiu, & 
    1420          & igrdju, & 
    1421          & igrdiv, & 
    1422          & igrdjv, & 
    1423          & zglamu, zglamv, & 
    1424          & zgphiu, zgphiv, & 
    1425          & zumask, zvmask, & 
    1426          & zintu, & 
    1427          & zintv  & 
    1428          & ) 
    1429       ! At the end of the day also get interpolated means 
    1430       IF ( idayend == 0 ) THEN 
    1431          DEALLOCATE( & 
    1432             & zinmu,  & 
    1433             & zinmv   & 
    1434             & ) 
    1435       ENDIF 
    1436  
    1437       prodatqc%nprofup = prodatqc%nprofup + ipro  
    1438        
    1439    END SUBROUTINE obs_vel_opt 
     829 
     830      surfdataqc%nsurfup = surfdataqc%nsurfup + isurf 
     831 
     832   END SUBROUTINE obs_surf_opt 
    1440833 
    1441834END MODULE obs_oper 
    1442  
Note: See TracChangeset for help on using the changeset viewer.