New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 9023 for branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

Ignore:
Timestamp:
2017-12-13T18:08:50+01:00 (6 years ago)
Author:
timgraham
Message:

Merged METO_MERCATOR branch and resolved all conflicts in OPA_SRC

Location:
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/OBS
Files:
12 edited
1 copied

Legend:

Unmodified
Added
Removed
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90

    r9019 r9023  
    4646   LOGICAL, PUBLIC :: ln_diaobs   !: Logical switch for the obs operator 
    4747   LOGICAL :: ln_sstnight         !: Logical switch for night mean SST obs 
    48     
     48   LOGICAL :: ln_sla_fp_indegs    !: T=> SLA obs footprint size specified in degrees, F=> in metres 
     49   LOGICAL :: ln_sst_fp_indegs    !: T=> SST obs footprint size specified in degrees, F=> in metres 
     50   LOGICAL :: ln_sss_fp_indegs    !: T=> SSS obs footprint size specified in degrees, F=> in metres 
     51   LOGICAL :: ln_sic_fp_indegs    !: T=> sea-ice obs footprint size specified in degrees, F=> in metres 
     52 
     53   REAL(wp) :: rn_sla_avglamscl !: E/W diameter of SLA observation footprint (metres) 
     54   REAL(wp) :: rn_sla_avgphiscl !: N/S diameter of SLA observation footprint (metres) 
     55   REAL(wp) :: rn_sst_avglamscl !: E/W diameter of SST observation footprint (metres) 
     56   REAL(wp) :: rn_sst_avgphiscl !: N/S diameter of SST observation footprint (metres) 
     57   REAL(wp) :: rn_sss_avglamscl !: E/W diameter of SSS observation footprint (metres) 
     58   REAL(wp) :: rn_sss_avgphiscl !: N/S diameter of SSS observation footprint (metres) 
     59   REAL(wp) :: rn_sic_avglamscl !: E/W diameter of sea-ice observation footprint (metres) 
     60   REAL(wp) :: rn_sic_avgphiscl !: N/S diameter of sea-ice observation footprint (metres) 
     61 
    4962   INTEGER :: nn_1dint       !: Vertical interpolation method 
    50    INTEGER :: nn_2dint       !: Horizontal interpolation method 
     63   INTEGER :: nn_2dint       !: Default horizontal interpolation method 
     64   INTEGER :: nn_2dint_sla   !: SLA horizontal interpolation method  
     65   INTEGER :: nn_2dint_sst   !: SST horizontal interpolation method  
     66   INTEGER :: nn_2dint_sss   !: SSS horizontal interpolation method  
     67   INTEGER :: nn_2dint_sic   !: Seaice horizontal interpolation method  
    5168   INTEGER, DIMENSION(imaxavtypes) ::   nn_profdavtypes   !: Profile data types representing a daily average 
    5269   INTEGER :: nproftypes     !: Number of profile obs types 
    5370   INTEGER :: nsurftypes     !: Number of surface obs types 
    54    INTEGER, DIMENSION(:), ALLOCATABLE ::   nvarsprof, nvarssurf   !: Number of profile & surface variables 
    55    INTEGER, DIMENSION(:), ALLOCATABLE ::   nextrprof, nextrsurf   !: Number of profile & surface extra variables 
    56    INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:) ::   sstbias_type   !: SST bias type     
    57    TYPE(obs_surf), PUBLIC, POINTER, DIMENSION(:) ::   surfdata, surfdataqc   !: Initial surface data before & after quality control 
    58    TYPE(obs_prof), PUBLIC, POINTER, DIMENSION(:) ::   profdata, profdataqc   !: Initial profile data before & after quality control 
     71   INTEGER, DIMENSION(:), ALLOCATABLE :: & 
     72      & nvarsprof, &         !: Number of profile variables 
     73      & nvarssurf            !: Number of surface variables 
     74   INTEGER, DIMENSION(:), ALLOCATABLE :: & 
     75      & nextrprof, &         !: Number of profile extra variables 
     76      & nextrsurf            !: Number of surface extra variables 
     77   INTEGER, DIMENSION(:), ALLOCATABLE :: & 
     78      & n2dintsurf           !: Interpolation option for surface variables 
     79   REAL(wp), DIMENSION(:), ALLOCATABLE :: & 
     80      & zavglamscl, &        !: E/W diameter of averaging footprint for surface variables 
     81      & zavgphiscl           !: N/S diameter of averaging footprint for surface variables 
     82   LOGICAL, DIMENSION(:), ALLOCATABLE :: & 
     83      & lfpindegs, &         !: T=> surface obs footprint size specified in degrees, F=> in metres 
     84      & llnightav            !: Logical for calculating night-time averages 
     85 
     86   TYPE(obs_surf), PUBLIC, POINTER, DIMENSION(:) :: & 
     87      & surfdata, &          !: Initial surface data 
     88      & surfdataqc           !: Surface data after quality control 
     89   TYPE(obs_prof), PUBLIC, POINTER, DIMENSION(:) :: & 
     90      & profdata, &          !: Initial profile data 
     91      & profdataqc           !: Profile data after quality control 
    5992 
    6093   CHARACTER(len=6), PUBLIC, DIMENSION(:), ALLOCATABLE ::   cobstypesprof, cobstypessurf   !: Profile & surface obs types 
     
    95128         & cn_profbfiles, &      ! T/S profile input filenames 
    96129         & cn_sstfbfiles, &      ! Sea surface temperature input filenames 
     130         & cn_sssfbfiles, &      ! Sea surface salinity input filenames 
    97131         & cn_slafbfiles, &      ! Sea level anomaly input filenames 
    98132         & cn_sicfbfiles, &      ! Seaice concentration input filenames 
    99133         & cn_velfbfiles, &      ! Velocity profile input filenames 
    100          & cn_sstbias_files      ! SST bias input filenames 
     134         & cn_sstbiasfiles      ! SST bias input filenames 
    101135      CHARACTER(LEN=128) :: & 
    102136         & cn_altbiasfile        ! Altimeter bias input filename 
     
    109143      LOGICAL :: ln_sla          ! Logical switch for sea level anomalies  
    110144      LOGICAL :: ln_sst          ! Logical switch for sea surface temperature 
     145      LOGICAL :: ln_sss          ! Logical switch for sea surface salinity 
    111146      LOGICAL :: ln_sic          ! Logical switch for sea ice concentration 
    112147      LOGICAL :: ln_vel3d        ! Logical switch for velocity (u,v) obs 
     
    116151      LOGICAL :: ln_ignmis       ! Logical switch for ignoring missing files 
    117152      LOGICAL :: ln_s_at_t       ! Logical switch to compute model S at T obs 
     153      LOGICAL :: ln_bound_reject ! Logical to remove obs near boundaries in LAMs. 
    118154      LOGICAL :: llvar1          ! Logical for profile variable 1 
    119155      LOGICAL :: llvar2          ! Logical for profile variable 1 
    120       LOGICAL :: llnightav       ! Logical for calculating night-time averages 
    121156      LOGICAL, DIMENSION(jpmaxnfiles) :: lmask ! Used for finding number of sstbias files 
    122157 
     
    134169 
    135170      NAMELIST/namobs/ln_diaobs, ln_t3d, ln_s3d, ln_sla,              & 
    136          &            ln_sst, ln_sic, ln_vel3d,                       & 
    137          &            ln_altbias, ln_nea, ln_grid_global,             & 
    138          &            ln_grid_search_lookup,                          & 
    139          &            ln_ignmis, ln_s_at_t, ln_sstnight,              & 
     171         &            ln_sst, ln_sic, ln_sss, ln_vel3d,               & 
     172         &            ln_altbias, ln_sstbias, ln_nea,                 & 
     173         &            ln_grid_global, ln_grid_search_lookup,          & 
     174         &            ln_ignmis, ln_s_at_t, ln_bound_reject,          & 
     175         &            ln_sstnight,                                    & 
     176         &            ln_sla_fp_indegs, ln_sst_fp_indegs,             & 
     177         &            ln_sss_fp_indegs, ln_sic_fp_indegs,             & 
    140178         &            cn_profbfiles, cn_slafbfiles,                   & 
    141179         &            cn_sstfbfiles, cn_sicfbfiles,                   & 
    142          &            cn_velfbfiles, cn_altbiasfile,                  & 
     180         &            cn_velfbfiles, cn_sssfbfiles,                   & 
     181         &            cn_sstbiasfiles, cn_altbiasfile,                & 
    143182         &            cn_gridsearchfile, rn_gridsearchres,            & 
    144          &            rn_dobsini, rn_dobsend, nn_1dint, nn_2dint,     & 
     183         &            rn_dobsini, rn_dobsend,                         & 
     184         &            rn_sla_avglamscl, rn_sla_avgphiscl,             & 
     185         &            rn_sst_avglamscl, rn_sst_avgphiscl,             & 
     186         &            rn_sss_avglamscl, rn_sss_avgphiscl,             & 
     187         &            rn_sic_avglamscl, rn_sic_avgphiscl,             & 
     188         &            nn_1dint, nn_2dint,                             & 
     189         &            nn_2dint_sla, nn_2dint_sst,                     & 
     190         &            nn_2dint_sss, nn_2dint_sic,                     & 
    145191         &            nn_msshc, rn_mdtcorr, rn_mdtcutoff,             & 
    146          &            nn_profdavtypes, ln_sstbias, cn_sstbias_files 
     192         &            nn_profdavtypes 
    147193 
    148194      INTEGER :: jnumsstbias 
     
    157203      ! Read namelist parameters 
    158204      !----------------------------------------------------------------------- 
    159        
    160       !Initalise all values in namelist arrays 
    161       ALLOCATE(sstbias_type(jpmaxnfiles)) 
    162205      ! Some namelist arrays need initialising 
    163206      cn_profbfiles(:) = '' 
     
    166209      cn_sicfbfiles(:) = '' 
    167210      cn_velfbfiles(:) = '' 
    168       cn_sstbias_files(:) = '' 
     211      cn_sssfbfiles(:)    = '' 
     212      cn_sstbiasfiles(:) = '' 
    169213      nn_profdavtypes(:) = -1 
    170214 
     
    187231         RETURN 
    188232      ENDIF 
    189        
    190       !----------------------------------------------------------------------- 
    191       ! Set up list of observation types to be used 
    192       ! and the files associated with each type 
    193       !----------------------------------------------------------------------- 
    194  
    195       nproftypes = COUNT( (/ln_t3d .OR. ln_s3d, ln_vel3d /) ) 
    196       nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic /) ) 
    197  
    198       IF (ln_sstbias) THEN  
    199          lmask(:) = .FALSE.  
    200          WHERE (cn_sstbias_files(:) /= '') lmask(:) = .TRUE.  
    201          jnumsstbias = COUNT(lmask)  
    202          lmask(:) = .FALSE.  
    203       ENDIF       
    204  
    205       IF ( nproftypes == 0 .AND. nsurftypes == 0 ) THEN 
    206          IF(lwp) WRITE(numout,cform_war) 
    207          IF(lwp) WRITE(numout,*) ' ln_diaobs is set to true, but all obs operator logical flags', & 
    208             &                    ' ln_t3d, ln_s3d, ln_sla, ln_sst, ln_sic, ln_vel3d', & 
    209             &                    ' are set to .FALSE. so turning off calls to dia_obs' 
    210          nwarn = nwarn + 1 
    211          ln_diaobs = .FALSE. 
    212          RETURN 
    213       ENDIF 
    214  
    215       IF ( nproftypes > 0 ) THEN 
    216  
    217          ALLOCATE( cobstypesprof(nproftypes) ) 
    218          ALLOCATE( ifilesprof(nproftypes) ) 
    219          ALLOCATE( clproffiles(nproftypes,jpmaxnfiles) ) 
    220  
    221          jtype = 0 
    222          IF (ln_t3d .OR. ln_s3d) THEN 
    223             jtype = jtype + 1 
    224             clproffiles(jtype,:) = cn_profbfiles(:) 
    225             cobstypesprof(jtype) = 'prof  ' 
    226             ifilesprof(jtype) = 0 
    227             DO jfile = 1, jpmaxnfiles 
    228                IF ( trim(clproffiles(jtype,jfile)) /= '' ) & 
    229                   ifilesprof(jtype) = ifilesprof(jtype) + 1 
    230             END DO 
    231          ENDIF 
    232          IF (ln_vel3d) THEN 
    233             jtype = jtype + 1 
    234             clproffiles(jtype,:) = cn_velfbfiles(:) 
    235             cobstypesprof(jtype) = 'vel   ' 
    236             ifilesprof(jtype) = 0 
    237             DO jfile = 1, jpmaxnfiles 
    238                IF ( trim(clproffiles(jtype,jfile)) /= '' ) & 
    239                   ifilesprof(jtype) = ifilesprof(jtype) + 1 
    240             END DO 
    241          ENDIF 
    242  
    243       ENDIF 
    244  
    245       IF ( nsurftypes > 0 ) THEN 
    246  
    247          ALLOCATE( cobstypessurf(nsurftypes) ) 
    248          ALLOCATE( ifilessurf(nsurftypes) ) 
    249          ALLOCATE( clsurffiles(nsurftypes, jpmaxnfiles) ) 
    250  
    251          jtype = 0 
    252          IF (ln_sla) THEN 
    253             jtype = jtype + 1 
    254             clsurffiles(jtype,:) = cn_slafbfiles(:) 
    255             cobstypessurf(jtype) = 'sla   ' 
    256             ifilessurf(jtype) = 0 
    257             DO jfile = 1, jpmaxnfiles 
    258                IF ( trim(clsurffiles(jtype,jfile)) /= '' ) & 
    259                   ifilessurf(jtype) = ifilessurf(jtype) + 1 
    260             END DO 
    261          ENDIF 
    262          IF (ln_sst) THEN 
    263             jtype = jtype + 1 
    264             clsurffiles(jtype,:) = cn_sstfbfiles(:) 
    265             cobstypessurf(jtype) = 'sst   ' 
    266             ifilessurf(jtype) = 0 
    267             DO jfile = 1, jpmaxnfiles 
    268                IF ( trim(clsurffiles(jtype,jfile)) /= '' ) & 
    269                   ifilessurf(jtype) = ifilessurf(jtype) + 1 
    270             END DO 
    271          ENDIF 
    272 #if defined key_lim3 
    273          IF (ln_sic) THEN 
    274             jtype = jtype + 1 
    275             clsurffiles(jtype,:) = cn_sicfbfiles(:) 
    276             cobstypessurf(jtype) = 'sic   ' 
    277             ifilessurf(jtype) = 0 
    278             DO jfile = 1, jpmaxnfiles 
    279                IF ( trim(clsurffiles(jtype,jfile)) /= '' ) & 
    280                   ifilessurf(jtype) = ifilessurf(jtype) + 1 
    281             END DO 
    282          ENDIF 
    283 #endif 
    284  
    285       ENDIF 
    286  
    287       !Write namelist settings to stdout 
     233 
    288234      IF(lwp) THEN 
    289235         WRITE(numout,*) 
     
    297243         WRITE(numout,*) '             Logical switch for Sea Ice observations                  ln_sic = ', ln_sic 
    298244         WRITE(numout,*) '             Logical switch for velocity observations               ln_vel3d = ', ln_vel3d 
    299          WRITE(numout,*) '             Global distribution of observations              ln_grid_global = ',ln_grid_global 
    300          WRITE(numout,*) '             Logical switch for SST bias correction         ln_sstbias = ', ln_sstbias  
    301          WRITE(numout,*) '             Logical switch for obs grid search lookup ln_grid_search_lookup = ',ln_grid_search_lookup 
     245         WRITE(numout,*) '             Logical switch for SSS observations                      ln_sss = ', ln_sss 
     246         WRITE(numout,*) '             Global distribution of observations              ln_grid_global = ', ln_grid_global 
     247         WRITE(numout,*) '             Logical switch for obs grid search lookup ln_grid_search_lookup = ', ln_grid_search_lookup 
    302248         IF (ln_grid_search_lookup) & 
    303249            WRITE(numout,*) '             Grid search lookup file header                cn_gridsearchfile = ', cn_gridsearchfile 
     
    307253         WRITE(numout,*) '             Type of horizontal interpolation method                nn_2dint = ', nn_2dint 
    308254         WRITE(numout,*) '             Rejection of observations near land switch               ln_nea = ', ln_nea 
     255         WRITE(numout,*) '             Rejection of obs near open bdys                 ln_bound_reject = ', ln_bound_reject 
    309256         WRITE(numout,*) '             MSSH correction scheme                                 nn_msshc = ', nn_msshc 
    310257         WRITE(numout,*) '             MDT  correction                                      rn_mdtcorr = ', rn_mdtcorr 
    311258         WRITE(numout,*) '             MDT cutoff for computed correction                 rn_mdtcutoff = ', rn_mdtcutoff 
    312259         WRITE(numout,*) '             Logical switch for alt bias                          ln_altbias = ', ln_altbias 
     260         WRITE(numout,*) '             Logical switch for sst bias                          ln_sstbias = ', ln_sstbias 
    313261         WRITE(numout,*) '             Logical switch for ignoring missing files             ln_ignmis = ', ln_ignmis 
    314262         WRITE(numout,*) '             Daily average types                             nn_profdavtypes = ', nn_profdavtypes 
    315263         WRITE(numout,*) '             Logical switch for night-time SST obs               ln_sstnight = ', ln_sstnight 
    316          WRITE(numout,*) '          Number of profile obs types: ',nproftypes 
    317  
    318          IF ( nproftypes > 0 ) THEN 
    319             DO jtype = 1, nproftypes 
    320                DO jfile = 1, ifilesprof(jtype) 
    321                   WRITE(numout,'(1X,2A)') '             '//cobstypesprof(jtype)//' input observation file names  = ', & 
    322                      TRIM(clproffiles(jtype,jfile)) 
    323                END DO 
    324             END DO 
     264      ENDIF 
     265      !----------------------------------------------------------------------- 
     266      ! Set up list of observation types to be used 
     267      ! and the files associated with each type 
     268      !----------------------------------------------------------------------- 
     269 
     270      nproftypes = COUNT( (/ln_t3d .OR. ln_s3d, ln_vel3d /) ) 
     271      nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic, ln_sss /) ) 
     272 
     273      IF (ln_sstbias) THEN  
     274         lmask(:) = .FALSE.  
     275         WHERE (cn_sstbiasfiles(:) /= '') lmask(:) = .TRUE.  
     276         jnumsstbias = COUNT(lmask)  
     277         lmask(:) = .FALSE.  
     278      ENDIF       
     279 
     280      IF ( nproftypes == 0 .AND. nsurftypes == 0 ) THEN 
     281         IF(lwp) WRITE(numout,cform_war) 
     282         IF(lwp) WRITE(numout,*) ' ln_diaobs is set to true, but all obs operator logical flags', & 
     283            &                    ' ln_t3d, ln_s3d, ln_sla, ln_sst, ln_sic, ln_vel3d', & 
     284            &                    ' are set to .FALSE. so turning off calls to dia_obs' 
     285         nwarn = nwarn + 1 
     286         ln_diaobs = .FALSE. 
     287         RETURN 
     288      ENDIF 
     289 
     290      IF ( nproftypes > 0 ) THEN 
     291 
     292         ALLOCATE( cobstypesprof(nproftypes) ) 
     293         ALLOCATE( ifilesprof(nproftypes) ) 
     294         ALLOCATE( clproffiles(nproftypes,jpmaxnfiles) ) 
     295 
     296         jtype = 0 
     297         IF (ln_t3d .OR. ln_s3d) THEN 
     298            jtype = jtype + 1 
     299            CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'prof  ', & 
     300               &                   cn_profbfiles, ifilesprof, cobstypesprof, clproffiles ) 
    325301         ENDIF 
    326  
    327          WRITE(numout,*)'          Number of surface obs types: ',nsurftypes 
    328          IF ( nsurftypes > 0 ) THEN 
    329             DO jtype = 1, nsurftypes 
    330                DO jfile = 1, ifilessurf(jtype) 
    331                   WRITE(numout,'(1X,2A)') '             '//cobstypessurf(jtype)//' input observation file names  = ', & 
    332                      TRIM(clsurffiles(jtype,jfile)) 
    333                END DO 
    334             END DO 
     302         IF (ln_vel3d) THEN 
     303            jtype = jtype + 1 
     304            CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'vel   ', & 
     305               &                   cn_velfbfiles, ifilesprof, cobstypesprof, clproffiles ) 
    335306         ENDIF 
    336          WRITE(numout,*) '~~~~~~~~~~~~' 
    337  
    338       ENDIF 
     307 
     308      ENDIF 
     309 
     310      IF ( nsurftypes > 0 ) THEN 
     311 
     312         ALLOCATE( cobstypessurf(nsurftypes) ) 
     313         ALLOCATE( ifilessurf(nsurftypes) ) 
     314         ALLOCATE( clsurffiles(nsurftypes, jpmaxnfiles) ) 
     315         ALLOCATE(n2dintsurf(nsurftypes)) 
     316         ALLOCATE(zavglamscl(nsurftypes)) 
     317         ALLOCATE(zavgphiscl(nsurftypes)) 
     318         ALLOCATE(lfpindegs(nsurftypes)) 
     319         ALLOCATE(llnightav(nsurftypes)) 
     320 
     321         jtype = 0 
     322         IF (ln_sla) THEN 
     323            jtype = jtype + 1 
     324            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sla   ', & 
     325               &                   cn_slafbfiles, ifilessurf, cobstypessurf, clsurffiles ) 
     326            CALL obs_setinterpopts( nsurftypes, jtype, 'sla   ',      & 
     327               &                  nn_2dint, nn_2dint_sla,             & 
     328               &                  rn_sla_avglamscl, rn_sla_avgphiscl, & 
     329               &                  ln_sla_fp_indegs, .FALSE.,          & 
     330               &                  n2dintsurf, zavglamscl, zavgphiscl, & 
     331               &                  lfpindegs, llnightav ) 
     332         ENDIF 
     333         IF (ln_sst) THEN 
     334            jtype = jtype + 1 
     335            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sst   ', & 
     336               &                   cn_sstfbfiles, ifilessurf, cobstypessurf, clsurffiles ) 
     337            CALL obs_setinterpopts( nsurftypes, jtype, 'sst   ',      & 
     338               &                  nn_2dint, nn_2dint_sst,             & 
     339               &                  rn_sst_avglamscl, rn_sst_avgphiscl, & 
     340               &                  ln_sst_fp_indegs, ln_sstnight,      & 
     341               &                  n2dintsurf, zavglamscl, zavgphiscl, & 
     342               &                  lfpindegs, llnightav ) 
     343         ENDIF 
     344#if defined key_lim3 || defined key_cice 
     345         IF (ln_sic) THEN 
     346            jtype = jtype + 1 
     347            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sic   ', & 
     348               &                   cn_sicfbfiles, ifilessurf, cobstypessurf, clsurffiles ) 
     349            CALL obs_setinterpopts( nsurftypes, jtype, 'sic   ',      & 
     350               &                  nn_2dint, nn_2dint_sic,             & 
     351               &                  rn_sic_avglamscl, rn_sic_avgphiscl, & 
     352               &                  ln_sic_fp_indegs, .FALSE.,          & 
     353               &                  n2dintsurf, zavglamscl, zavgphiscl, & 
     354               &                  lfpindegs, llnightav ) 
     355         ENDIF 
     356#endif 
     357         IF (ln_sss) THEN 
     358            jtype = jtype + 1 
     359            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sss   ', & 
     360               &                   cn_sssfbfiles, ifilessurf, cobstypessurf, clsurffiles ) 
     361            CALL obs_setinterpopts( nsurftypes, jtype, 'sss   ',      & 
     362               &                  nn_2dint, nn_2dint_sss,             & 
     363               &                  rn_sss_avglamscl, rn_sss_avgphiscl, & 
     364               &                  ln_sss_fp_indegs, .FALSE.,          & 
     365               &                  n2dintsurf, zavglamscl, zavgphiscl, & 
     366               &                  lfpindegs, llnightav ) 
     367         ENDIF 
     368 
     369      ENDIF 
     370 
     371 
    339372 
    340373      !----------------------------------------------------------------------- 
     
    356389      ENDIF 
    357390 
    358       IF ( ( nn_2dint < 0 ) .OR. ( nn_2dint > 4 ) ) THEN 
     391      IF ( ( nn_2dint < 0 ) .OR. ( nn_2dint > 6 ) ) THEN 
    359392         CALL ctl_stop(' Choice of horizontal (2D) interpolation method', & 
    360393            &                    ' is not available') 
     
    421454               &               jpi, jpj, jpk, & 
    422455               &               zmask1, zglam1, zgphi1, zmask2, zglam2, zgphi2,  & 
    423                &               ln_nea, kdailyavtypes = nn_profdavtypes ) 
     456               &               ln_nea, ln_bound_reject, & 
     457               &               kdailyavtypes = nn_profdavtypes ) 
    424458 
    425459         END DO 
     
    440474            nvarssurf(jtype) = 1 
    441475            nextrsurf(jtype) = 0 
    442             llnightav = .FALSE. 
     476            llnightav(jtype) = .FALSE. 
    443477            IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) nextrsurf(jtype) = 2 
    444             IF ( TRIM(cobstypessurf(jtype)) == 'sst' ) llnightav = ln_sstnight 
     478            IF ( TRIM(cobstypessurf(jtype)) == 'sst' ) llnightav(jtype) = ln_sstnight 
    445479 
    446480            !Read in surface obs types 
     
    448482               &               clsurffiles(jtype,1:ifilessurf(jtype)), & 
    449483               &               nvarssurf(jtype), nextrsurf(jtype), nitend-nit000+2, & 
    450                &               rn_dobsini, rn_dobsend, ln_ignmis, .FALSE., llnightav ) 
    451           
    452           
    453             CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea ) 
     484               &               rn_dobsini, rn_dobsend, ln_ignmis, .FALSE., llnightav(jtype) ) 
     485 
     486            CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea, ln_bound_reject ) 
    454487 
    455488            IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN 
    456                CALL obs_rea_mdt( surfdataqc(jtype), nn_2dint ) 
    457                IF ( ln_altbias ) CALL obs_rea_altbias ( surfdataqc(jtype), nn_2dint, cn_altbiasfile ) 
     489               CALL obs_rea_mdt( surfdataqc(jtype), n2dintsurf(jtype) ) 
     490               IF ( ln_altbias ) & 
     491                  & CALL obs_rea_altbias ( surfdataqc(jtype), n2dintsurf(jtype), cn_altbiasfile ) 
    458492            ENDIF 
     493 
    459494            IF ( TRIM(cobstypessurf(jtype)) == 'sst' .AND. ln_sstbias ) THEN 
    460                !Read in bias field and correct SST. 
    461                IF ( jnumsstbias == 0 ) CALL ctl_stop("ln_sstbias set,"// & 
    462                                                      "  but no bias"// & 
    463                                                      " files to read in")    
    464                   CALL obs_app_sstbias( surfdataqc(jtype), nn_2dint, & 
    465                                         jnumsstbias, cn_sstbias_files(1:jnumsstbias) ) 
     495               jnumsstbias = 0 
     496               DO jfile = 1, jpmaxnfiles 
     497                  IF ( TRIM(cn_sstbiasfiles(jfile)) /= '' ) & 
     498                     &  jnumsstbias = jnumsstbias + 1 
     499               END DO 
     500               IF ( jnumsstbias == 0 ) THEN 
     501                  CALL ctl_stop("ln_sstbias set but no bias files to read in")     
     502               ENDIF 
     503 
     504               CALL obs_app_sstbias( surfdataqc(jtype), n2dintsurf(jtype), &  
     505                  &                  jnumsstbias, cn_sstbiasfiles(1:jnumsstbias) )  
     506 
    466507            ENDIF 
    467508         END DO 
     
    512553      USE ice    , ONLY : at_i                ! LIM3 Ice model variables 
    513554#endif 
     555#if defined key_cice 
     556      USE sbc_oce, ONLY : fr_i     ! ice fraction 
     557#endif 
     558 
    514559      IMPLICIT NONE 
    515560 
     
    528573         & zprofmask2              ! Mask associated with zprofvar2 
    529574      REAL(wp), POINTER, DIMENSION(:,:) :: & 
    530          & zsurfvar                ! Model values equivalent to surface ob. 
     575         & zsurfvar, &             ! Model values equivalent to surface ob. 
     576         & zsurfmask               ! Mask associated with surface variable 
    531577      REAL(wp), POINTER, DIMENSION(:,:) :: & 
    532578         & zglam1,    &            ! Model longitudes for prof variable 1 
     
    534580         & zgphi1,    &            ! Model latitudes for prof variable 1 
    535581         & zgphi2                  ! Model latitudes for prof variable 2 
    536 #if ! defined key_lim3 
    537       REAL(wp), POINTER, DIMENSION(:,:) :: at_i 
    538 #endif 
    539       LOGICAL :: llnightav        ! Logical for calculating night-time average 
    540582 
    541583      !Allocate local work arrays 
     
    545587      CALL wrk_alloc( jpi, jpj, jpk, zprofmask2 ) 
    546588      CALL wrk_alloc( jpi, jpj, zsurfvar ) 
     589      CALL wrk_alloc( jpi, jpj, zsurfmask ) 
    547590      CALL wrk_alloc( jpi, jpj, zglam1 ) 
    548591      CALL wrk_alloc( jpi, jpj, zglam2 ) 
    549592      CALL wrk_alloc( jpi, jpj, zgphi1 ) 
    550593      CALL wrk_alloc( jpi, jpj, zgphi2 ) 
    551 #if ! defined key_lim3 
    552       CALL wrk_alloc(jpi,jpj,at_i)  
    553 #endif 
    554594      !----------------------------------------------------------------------- 
    555595 
     
    562602      idaystp = NINT( rday / rdt ) 
    563603 
    564       !----------------------------------------------------------------------- 
    565       ! No LIM => at_i == 0.0_wp 
    566       !----------------------------------------------------------------------- 
    567 #if ! defined key_lim3 
    568       at_i(:,:) = 0.0_wp 
    569 #endif 
    570604      !----------------------------------------------------------------------- 
    571605      ! Call the profile and surface observation operators 
     
    595629               zgphi1(:,:) = gphiu(:,:) 
    596630               zgphi2(:,:) = gphiv(:,:) 
     631            CASE DEFAULT 
     632               CALL ctl_stop( 'Unknown profile observation type '//TRIM(cobstypesprof(jtype))//' in dia_obs' ) 
    597633            END SELECT 
    598634 
    599             IF( ln_zco .OR. ln_zps ) THEN  
    600                CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk,  & 
    601                   &               nit000, idaystp,                         & 
    602                   &               zprofvar1, zprofvar2,                    & 
    603                   &               gdept_1d, zprofmask1, zprofmask2,        & 
    604                   &               zglam1, zglam2, zgphi1, zgphi2,          & 
    605                   &               nn_1dint, nn_2dint,                      & 
    606                   &               kdailyavtypes = nn_profdavtypes ) 
    607             ELSE IF(TRIM(cobstypesprof(jtype)) == 'prof') THEN 
    608                !TG - THIS NEEDS MODIFICATION TO MATCH SIMPLIFICATION 
    609                CALL obs_pro_sco_opt( profdataqc(jtype),                    &  
    610                   &              kstp, jpi, jpj, jpk, nit000, idaystp,   &  
    611                   &              zprofvar1, zprofvar2,                   &  
    612                   &              gdept_n(:,:,:), gdepw_n(:,:,:),           & 
    613                   &              tmask, nn_1dint, nn_2dint,              &  
    614                   &              kdailyavtypes = nn_profdavtypes )  
    615             ELSE 
    616                CALL ctl_stop('DIA_OBS: Generalised vertical interpolation not'// & 
    617                          'yet working for velocity data (turn off velocity observations') 
    618             ENDIF 
     635            CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk,  & 
     636               &               nit000, idaystp,                         & 
     637               &               zprofvar1, zprofvar2,                    & 
     638               &               gdept_n(:,:,:), gdepw_n(:,:,:),            &  
     639               &               zprofmask1, zprofmask2,                  & 
     640               &               zglam1, zglam2, zgphi1, zgphi2,          & 
     641               &               nn_1dint, nn_2dint,                      & 
     642               &               kdailyavtypes = nn_profdavtypes ) 
    619643 
    620644         END DO 
     
    625649 
    626650         DO jtype = 1, nsurftypes 
     651 
     652            !Defaults which might be changed 
     653            zsurfmask(:,:) = tmask(:,:,1) 
    627654 
    628655            SELECT CASE ( TRIM(cobstypessurf(jtype)) ) 
    629656            CASE('sst') 
    630657               zsurfvar(:,:) = tsn(:,:,1,jp_tem) 
    631                llnightav = ln_sstnight 
    632658            CASE('sla') 
    633659               zsurfvar(:,:) = sshn(:,:) 
    634                llnightav = .FALSE. 
    635 #if defined key_lim3 
     660            CASE('sss') 
     661               zsurfvar(:,:) = tsn(:,:,1,jp_sal) 
    636662            CASE('sic') 
    637663               IF ( kstp == 0 ) THEN 
     
    646672                  CYCLE 
    647673               ELSE 
    648                   zsurfvar(:,:) = at_i(:,:) 
     674#if defined key_cice 
     675                  zsurfvar(:,:) = fr_i(:,:) 
     676#elif defined key_lim2 || defined key_lim3 
     677                  zsurfvar(:,:) = 1._wp - frld(:,:) 
     678#else 
     679               CALL ctl_stop( ' Trying to run sea-ice observation operator', & 
     680                  &           ' but no sea-ice model appears to have been defined' ) 
     681#endif 
    649682               ENDIF 
    650683 
    651                llnightav = .FALSE. 
    652 #endif 
    653684            END SELECT 
    654685 
    655686            CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj,       & 
    656                &               nit000, idaystp, zsurfvar, tmask(:,:,1), & 
    657                &               nn_2dint, llnightav ) 
     687               &               nit000, idaystp, zsurfvar, zsurfmask,    & 
     688               &               n2dintsurf(jtype), llnightav(jtype),     & 
     689               &               zavglamscl(jtype), zavgphiscl(jtype),     & 
     690               &               lfpindegs(jtype) ) 
    658691 
    659692         END DO 
     
    666699      CALL wrk_dealloc( jpi, jpj, jpk, zprofmask2 ) 
    667700      CALL wrk_dealloc( jpi, jpj, zsurfvar ) 
     701      CALL wrk_dealloc( jpi, jpj, zsurfmask ) 
    668702      CALL wrk_dealloc( jpi, jpj, zglam1 ) 
    669703      CALL wrk_dealloc( jpi, jpj, zglam2 ) 
    670704      CALL wrk_dealloc( jpi, jpj, zgphi1 ) 
    671705      CALL wrk_dealloc( jpi, jpj, zgphi2 ) 
    672 #if ! defined key_lim3 
    673       CALL wrk_dealloc(jpi,jpj,at_i) 
    674 #endif 
    675706 
    676707   END SUBROUTINE dia_obs 
     
    789820 
    790821      IF ( nsurftypes > 0 ) & 
    791          &   DEALLOCATE( cobstypessurf, surfdata, surfdataqc, nvarssurf, nextrsurf ) 
     822         &   DEALLOCATE( cobstypessurf, surfdata, surfdataqc, nvarssurf, nextrsurf, & 
     823         &               n2dintsurf, zavglamscl, zavgphiscl, lfpindegs, llnightav ) 
    792824 
    793825   END SUBROUTINE dia_obs_dealloc 
     
    938970   END SUBROUTINE fin_date 
    939971    
     972    SUBROUTINE obs_settypefiles( ntypes, jpmaxnfiles, jtype, ctypein, & 
     973       &                         cfilestype, ifiles, cobstypes, cfiles ) 
     974 
     975    INTEGER, INTENT(IN) :: ntypes      ! Total number of obs types 
     976    INTEGER, INTENT(IN) :: jpmaxnfiles ! Maximum number of files allowed for each type 
     977    INTEGER, INTENT(IN) :: jtype       ! Index of the current type of obs 
     978    INTEGER, DIMENSION(ntypes), INTENT(INOUT) :: & 
     979       &                   ifiles      ! Out appended number of files for this type 
     980 
     981    CHARACTER(len=6), INTENT(IN) :: ctypein  
     982    CHARACTER(len=128), DIMENSION(jpmaxnfiles), INTENT(IN) :: & 
     983       &                   cfilestype  ! In list of files for this obs type 
     984    CHARACTER(len=6), DIMENSION(ntypes), INTENT(INOUT) :: & 
     985       &                   cobstypes   ! Out appended list of obs types 
     986    CHARACTER(len=128), DIMENSION(ntypes, jpmaxnfiles), INTENT(INOUT) :: & 
     987       &                   cfiles      ! Out appended list of files for all types 
     988 
     989    !Local variables 
     990    INTEGER :: jfile 
     991 
     992    cfiles(jtype,:) = cfilestype(:) 
     993    cobstypes(jtype) = ctypein 
     994    ifiles(jtype) = 0 
     995    DO jfile = 1, jpmaxnfiles 
     996       IF ( trim(cfiles(jtype,jfile)) /= '' ) & 
     997                 ifiles(jtype) = ifiles(jtype) + 1 
     998    END DO 
     999 
     1000    IF ( ifiles(jtype) == 0 ) THEN 
     1001         CALL ctl_stop( 'Logical for observation type '//TRIM(ctypein)//   & 
     1002            &           ' set to true but no files available to read' ) 
     1003    ENDIF 
     1004 
     1005    IF(lwp) THEN     
     1006       WRITE(numout,*) '             '//cobstypes(jtype)//' input observation file names:' 
     1007       DO jfile = 1, ifiles(jtype) 
     1008          WRITE(numout,*) '                '//TRIM(cfiles(jtype,jfile)) 
     1009       END DO 
     1010    ENDIF 
     1011 
     1012    END SUBROUTINE obs_settypefiles 
     1013 
     1014    SUBROUTINE obs_setinterpopts( ntypes, jtype, ctypein,             & 
     1015               &                  n2dint_default, n2dint_type,        & 
     1016               &                  zavglamscl_type, zavgphiscl_type,   & 
     1017               &                  lfp_indegs_type, lavnight_type,     & 
     1018               &                  n2dint, zavglamscl, zavgphiscl,     & 
     1019               &                  lfpindegs, lavnight ) 
     1020 
     1021    INTEGER, INTENT(IN)  :: ntypes             ! Total number of obs types 
     1022    INTEGER, INTENT(IN)  :: jtype              ! Index of the current type of obs 
     1023    INTEGER, INTENT(IN)  :: n2dint_default     ! Default option for interpolation type 
     1024    INTEGER, INTENT(IN)  :: n2dint_type        ! Option for interpolation type 
     1025    REAL(wp), INTENT(IN) :: & 
     1026       &                    zavglamscl_type, & !E/W diameter of obs footprint for this type 
     1027       &                    zavgphiscl_type    !N/S diameter of obs footprint for this type 
     1028    LOGICAL, INTENT(IN)  :: lfp_indegs_type    !T=> footprint in degrees, F=> in metres 
     1029    LOGICAL, INTENT(IN)  :: lavnight_type      !T=> obs represent night time average 
     1030    CHARACTER(len=6), INTENT(IN) :: ctypein  
     1031 
     1032    INTEGER, DIMENSION(ntypes), INTENT(INOUT) :: & 
     1033       &                    n2dint  
     1034    REAL(wp), DIMENSION(ntypes), INTENT(INOUT) :: & 
     1035       &                    zavglamscl, zavgphiscl 
     1036    LOGICAL, DIMENSION(ntypes), INTENT(INOUT) :: & 
     1037       &                    lfpindegs, lavnight 
     1038 
     1039    lavnight(jtype) = lavnight_type 
     1040 
     1041    IF ( (n2dint_type >= 1) .AND. (n2dint_type <= 6) ) THEN 
     1042       n2dint(jtype) = n2dint_type 
     1043    ELSE 
     1044       n2dint(jtype) = n2dint_default 
     1045    ENDIF 
     1046 
     1047    ! For averaging observation footprints set options for size of footprint  
     1048    IF ( (n2dint(jtype) > 4) .AND. (n2dint(jtype) <= 6) ) THEN 
     1049       IF ( zavglamscl_type > 0._wp ) THEN 
     1050          zavglamscl(jtype) = zavglamscl_type 
     1051       ELSE 
     1052          CALL ctl_stop( 'Incorrect value set for averaging footprint '// & 
     1053                         'scale (zavglamscl) for observation type '//TRIM(ctypein) )       
     1054       ENDIF 
     1055 
     1056       IF ( zavgphiscl_type > 0._wp ) THEN 
     1057          zavgphiscl(jtype) = zavgphiscl_type 
     1058       ELSE 
     1059          CALL ctl_stop( 'Incorrect value set for averaging footprint '// & 
     1060                         'scale (zavgphiscl) for observation type '//TRIM(ctypein) )       
     1061       ENDIF 
     1062 
     1063       lfpindegs(jtype) = lfp_indegs_type  
     1064 
     1065    ENDIF 
     1066 
     1067    ! Write out info  
     1068    IF(lwp) THEN 
     1069       IF ( n2dint(jtype) <= 4 ) THEN 
     1070          WRITE(numout,*) '             '//TRIM(ctypein)// & 
     1071             &            ' model counterparts will be interpolated horizontally' 
     1072       ELSE IF ( n2dint(jtype) <= 6 ) THEN 
     1073          WRITE(numout,*) '             '//TRIM(ctypein)// & 
     1074             &            ' model counterparts will be averaged horizontally' 
     1075          WRITE(numout,*) '             '//'    with E/W scale: ',zavglamscl(jtype) 
     1076          WRITE(numout,*) '             '//'    with N/S scale: ',zavgphiscl(jtype) 
     1077          IF ( lfpindegs(jtype) ) THEN 
     1078              WRITE(numout,*) '             '//'    (in degrees)' 
     1079          ELSE 
     1080              WRITE(numout,*) '             '//'    (in metres)' 
     1081          ENDIF 
     1082       ENDIF 
     1083    ENDIF 
     1084 
     1085    END SUBROUTINE obs_setinterpopts 
     1086 
    9401087END MODULE diaobs 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/OBS/obs_mpp.F90

    r6140 r9023  
    142142      !! 
    143143 
    144       iobsp=kobsp 
     144      iobsp(:)=kobsp(:) 
    145145 
    146146      WHERE( iobsp(:) == -1 ) 
     
    148148      END WHERE 
    149149 
    150       iobsp=-1*iobsp 
     150      iobsp(:)=-1*iobsp(:) 
    151151 
    152152      CALL obs_mpp_max_integer( iobsp, kno ) 
    153153 
    154       kobsp=-1*iobsp 
     154      kobsp(:)=-1*iobsp(:) 
    155155 
    156156      isum=0 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90

    r7646 r9023  
    99   !!   obs_prof_opt :    Compute the model counterpart of profile data 
    1010   !!   obs_surf_opt :    Compute the model counterpart of surface data 
    11    !!   obs_pro_sco_opt: Compute the model counterpart of temperature and  
    12    !!                    salinity observations from profiles in generalised  
    13    !!                    vertical coordinates  
    1411   !!---------------------------------------------------------------------- 
    1512 
     
    2219      & obs_int_h2d, & 
    2320      & obs_int_h2d_init 
     21   USE obs_averg_h2d, ONLY : &    ! Horizontal averaging to the obs footprint 
     22      & obs_avg_h2d, & 
     23      & obs_avg_h2d_init, & 
     24      & obs_max_fpsize 
    2425   USE obs_inter_z1d, ONLY : &    ! Vertical interpolation to the obs pt 
    2526      & obs_int_z1d,    & 
    2627      & obs_int_z1d_spl 
    27    USE obs_const,  ONLY :     & 
    28       & obfillflt      ! Fillvalue    
     28   USE obs_const,  ONLY :    &    ! Obs fill value 
     29      & obfillflt 
    2930   USE dom_oce,       ONLY : & 
    30       & glamt, glamu, glamv, & 
    31       & gphit, gphiu, gphiv, &  
    32       & gdept_n, gdept_0  
    33    USE lib_mpp,       ONLY : & 
     31      & glamt, glamf, & 
     32      & gphit, gphif 
     33   USE lib_mpp,       ONLY : &    ! Warning and stopping routines 
    3434      & ctl_warn, ctl_stop 
     35   USE sbcdcy,        ONLY : &    ! For calculation of where it is night-time 
     36      & sbc_dcy, nday_qsr 
    3537   USE obs_grid,      ONLY : &  
    3638      & obs_level_search      
    37    USE sbcdcy,        ONLY : &    ! For calculation of where it is night-time 
    38       & sbc_dcy, nday_qsr 
    3939 
    4040   IMPLICIT NONE 
     
    4444 
    4545   PUBLIC obs_prof_opt, &  ! Compute the model counterpart of profile obs 
    46       &   obs_pro_sco_opt, &  ! Compute the model counterpart of profile observations  
    4746      &   obs_surf_opt     ! Compute the model counterpart of surface obs 
    4847 
     
    5857CONTAINS 
    5958 
     59 
    6060   SUBROUTINE obs_prof_opt( prodatqc, kt, kpi, kpj, kpk,          & 
    6161      &                     kit000, kdaystp,                      & 
    62       &                     pvar1, pvar2, pgdept, pmask1, pmask2, & 
     62      &                     pvar1, pvar2, pgdept, pgdepw,         & 
     63      &                     pmask1, pmask2,                       &   
    6364      &                     plam1, plam2, pphi1, pphi2,           & 
    6465      &                     k1dint, k2dint, kdailyavtypes ) 
     
    111112      !!      ! 07-03 (K. Mogensen) General handling of profiles 
    112113      !!      ! 15-02 (M. Martin) Combined routine for all profile types 
     114      !!      ! 17-02 (M. Martin) Include generalised vertical coordinate changes 
    113115      !!----------------------------------------------------------------------- 
    114116 
     
    140142         & pphi1,    &               ! Model latitudes for variable 1 
    141143         & pphi2                     ! Model latitudes for variable 2 
    142       REAL(KIND=wp), INTENT(IN), DIMENSION(kpk) :: & 
    143          & pgdept                    ! Model array of depth levels 
     144      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &  
     145         & pgdept, &                 ! Model array of depth T levels  
     146         & pgdepw                    ! Model array of depth W levels  
    144147      INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 
    145148         & kdailyavtypes             ! Types for daily averages 
     
    156159      INTEGER ::   iend 
    157160      INTEGER ::   iobs 
     161      INTEGER ::   iin, ijn, ikn, ik   ! looping indices over interpolation nodes  
     162      INTEGER ::   inum_obs 
    158163      INTEGER, DIMENSION(imaxavtypes) :: & 
    159164         & idailyavtypes 
     
    163168         & igrdj1, & 
    164169         & igrdj2 
     170      INTEGER, ALLOCATABLE, DIMENSION(:) :: iv_indic 
     171 
    165172      REAL(KIND=wp) :: zlam 
    166173      REAL(KIND=wp) :: zphi 
     
    171178         & zobsk,    & 
    172179         & zobs2k 
    173       REAL(KIND=wp), DIMENSION(2,2,kpk) :: & 
     180      REAL(KIND=wp), DIMENSION(2,2,1) :: & 
    174181         & zweig1, & 
    175          & zweig2 
     182         & zweig2, & 
     183         & zweig 
    176184      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: & 
    177185         & zmask1, & 
    178186         & zmask2, & 
    179          & zint1, & 
    180          & zint2, & 
    181          & zinm1, & 
    182          & zinm2 
     187         & zint1,  & 
     188         & zint2,  & 
     189         & zinm1,  & 
     190         & zinm2,  & 
     191         & zgdept, &  
     192         & zgdepw 
    183193      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 
    184194         & zglam1, & 
     
    186196         & zgphi1, & 
    187197         & zgphi2 
     198      REAL(KIND=wp), DIMENSION(1) :: zmsk_1, zmsk_2    
     199      REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner 
     200 
    188201      LOGICAL :: ld_dailyav 
    189202 
     
    266279         & zmask1(2,2,kpk,ipro),  & 
    267280         & zmask2(2,2,kpk,ipro),  & 
    268          & zint1(2,2,kpk,ipro),  & 
    269          & zint2(2,2,kpk,ipro)   & 
     281         & zint1(2,2,kpk,ipro),   & 
     282         & zint2(2,2,kpk,ipro),   & 
     283         & zgdept(2,2,kpk,ipro),  &  
     284         & zgdepw(2,2,kpk,ipro)   &  
    270285         & ) 
    271286 
     
    290305      END DO 
    291306 
     307      ! Initialise depth arrays 
     308      zgdept(:,:,:,:) = 0.0 
     309      zgdepw(:,:,:,:) = 0.0 
     310 
    292311      CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi1, igrdj1, plam1, zglam1 ) 
    293312      CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi1, igrdj1, pphi1, zgphi1 ) 
     
    300319      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, pvar2,   zint2 ) 
    301320 
     321      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pgdept, zgdept )  
     322      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pgdepw, zgdepw )  
     323 
    302324      ! At the end of the day also get interpolated means 
    303325      IF ( ld_dailyav .AND. idayend == 0 ) THEN 
     
    314336 
    315337      ENDIF 
     338 
     339      ! Return if no observations to process  
     340      ! Has to be done after comm commands to ensure processors  
     341      ! stay in sync  
     342      IF ( ipro == 0 ) RETURN  
    316343 
    317344      DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro 
     
    339366         zphi = prodatqc%rphi(jobs) 
    340367 
    341          ! Horizontal weights and vertical mask 
    342  
     368         ! Horizontal weights  
     369         ! Masked values are calculated later.   
    343370         IF ( prodatqc%npvend(jobs,1) > 0 ) THEN 
    344371 
    345             CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi,     & 
     372            CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,     & 
    346373               &                   zglam1(:,:,iobs), zgphi1(:,:,iobs), & 
    347                &                   zmask1(:,:,:,iobs), zweig1, zobsmask1 ) 
     374               &                   zmask1(:,:,1,iobs), zweig1, zmsk_1 ) 
    348375 
    349376         ENDIF 
     
    351378         IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 
    352379 
    353             CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi,     & 
     380            CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,     & 
    354381               &                   zglam2(:,:,iobs), zgphi2(:,:,iobs), & 
    355                &                   zmask2(:,:,:,iobs), zweig2, zobsmask2 ) 
     382               &                   zmask2(:,:,1,iobs), zweig2, zmsk_2 ) 
    356383  
    357384         ENDIF 
     
    365392               IF ( idayend == 0 )  THEN 
    366393                  ! Daily averaged data 
    367                   CALL obs_int_h2d( kpk, kpk,      & 
    368                      &              zweig1, zinm1(:,:,:,iobs), zobsk ) 
    369  
    370                ENDIF 
    371  
    372             ELSE  
    373  
    374                ! Point data 
    375                CALL obs_int_h2d( kpk, kpk,      & 
    376                   &              zweig1, zint1(:,:,:,iobs), zobsk ) 
    377  
    378             ENDIF 
    379  
    380             !------------------------------------------------------------- 
    381             ! Compute vertical second-derivative of the interpolating  
    382             ! polynomial at obs points 
    383             !------------------------------------------------------------- 
    384  
    385             IF ( k1dint == 1 ) THEN 
    386                CALL obs_int_z1d_spl( kpk, zobsk, zobs2k,   & 
    387                   &                  pgdept, zobsmask1 ) 
    388             ENDIF 
    389  
    390             !----------------------------------------------------------------- 
    391             !  Vertical interpolation to the observation point 
    392             !----------------------------------------------------------------- 
    393             ista = prodatqc%npvsta(jobs,1) 
    394             iend = prodatqc%npvend(jobs,1) 
    395             CALL obs_int_z1d( kpk,                & 
    396                & prodatqc%var(1)%mvk(ista:iend),  & 
    397                & k1dint, iend - ista + 1,         & 
    398                & prodatqc%var(1)%vdep(ista:iend), & 
    399                & zobsk, zobs2k,                   & 
    400                & prodatqc%var(1)%vmod(ista:iend), & 
    401                & pgdept, zobsmask1 ) 
    402  
    403          ENDIF 
    404  
     394 
     395                  ! vertically interpolate all 4 corners  
     396                  ista = prodatqc%npvsta(jobs,1)  
     397                  iend = prodatqc%npvend(jobs,1)  
     398                  inum_obs = iend - ista + 1  
     399                  ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs))  
     400 
     401                  DO iin=1,2  
     402                     DO ijn=1,2  
     403 
     404                        IF ( k1dint == 1 ) THEN  
     405                           CALL obs_int_z1d_spl( kpk, &  
     406                              &     zinm1(iin,ijn,:,iobs), &  
     407                              &     zobs2k, zgdept(iin,ijn,:,iobs), &  
     408                              &     zmask1(iin,ijn,:,iobs))  
     409                        ENDIF  
     410        
     411                        CALL obs_level_search(kpk, &  
     412                           &    zgdept(iin,ijn,:,iobs), &  
     413                           &    inum_obs, prodatqc%var(1)%vdep(ista:iend), &  
     414                           &    iv_indic)  
     415 
     416                        CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &  
     417                           &    prodatqc%var(1)%vdep(ista:iend), &  
     418                           &    zinm1(iin,ijn,:,iobs), &  
     419                           &    zobs2k, interp_corner(iin,ijn,:), &  
     420                           &    zgdept(iin,ijn,:,iobs), &  
     421                           &    zmask1(iin,ijn,:,iobs))  
     422        
     423                     ENDDO  
     424                  ENDDO  
     425 
     426               ENDIF !idayend 
     427 
     428            ELSE    
     429 
     430               ! Point data  
     431      
     432               ! vertically interpolate all 4 corners  
     433               ista = prodatqc%npvsta(jobs,1)  
     434               iend = prodatqc%npvend(jobs,1)  
     435               inum_obs = iend - ista + 1  
     436               ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs))  
     437               DO iin=1,2   
     438                  DO ijn=1,2  
     439                     
     440                     IF ( k1dint == 1 ) THEN  
     441                        CALL obs_int_z1d_spl( kpk, &  
     442                           &    zint1(iin,ijn,:,iobs),&  
     443                           &    zobs2k, zgdept(iin,ijn,:,iobs), &  
     444                           &    zmask1(iin,ijn,:,iobs))  
     445   
     446                     ENDIF  
     447        
     448                     CALL obs_level_search(kpk, &  
     449                         &        zgdept(iin,ijn,:,iobs),&  
     450                         &        inum_obs, prodatqc%var(1)%vdep(ista:iend), &  
     451                         &        iv_indic)  
     452 
     453                     CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs,     &  
     454                         &          prodatqc%var(1)%vdep(ista:iend),     &  
     455                         &          zint1(iin,ijn,:,iobs),            &  
     456                         &          zobs2k,interp_corner(iin,ijn,:), &  
     457                         &          zgdept(iin,ijn,:,iobs),         &  
     458                         &          zmask1(iin,ijn,:,iobs) )       
     459          
     460                  ENDDO  
     461               ENDDO  
     462              
     463            ENDIF  
     464 
     465            !-------------------------------------------------------------  
     466            ! Compute the horizontal interpolation for every profile level  
     467            !-------------------------------------------------------------  
     468              
     469            DO ikn=1,inum_obs  
     470               iend=ista+ikn-1 
     471                   
     472               zweig(:,:,1) = 0._wp  
     473    
     474               ! This code forces the horizontal weights to be   
     475               ! zero IF the observation is below the bottom of the   
     476               ! corners of the interpolation nodes, Or if it is in   
     477               ! the mask. This is important for observations near   
     478               ! steep bathymetry  
     479               DO iin=1,2  
     480                  DO ijn=1,2  
     481      
     482                     depth_loop1: DO ik=kpk,2,-1  
     483                        IF(zmask1(iin,ijn,ik-1,iobs ) > 0.9 )THEN    
     484                             
     485                           zweig(iin,ijn,1) = &   
     486                              & zweig1(iin,ijn,1) * &  
     487                              & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) &  
     488                              &  - prodatqc%var(1)%vdep(iend)),0._wp)  
     489                             
     490                           EXIT depth_loop1  
     491 
     492                        ENDIF  
     493 
     494                     ENDDO depth_loop1  
     495      
     496                  ENDDO  
     497               ENDDO  
     498    
     499               CALL obs_int_h2d( 1, 1, zweig, interp_corner(:,:,ikn), &  
     500                  &              prodatqc%var(1)%vmod(iend:iend) )  
     501 
     502                  ! Set QC flag for any observations found below the bottom 
     503                  ! needed as the check here is more strict than that in obs_prep 
     504               IF (sum(zweig) == 0.0_wp) prodatqc%var(1)%nvqc(iend:iend)=4 
     505  
     506            ENDDO  
     507  
     508            DEALLOCATE(interp_corner,iv_indic)  
     509           
     510         ENDIF  
     511 
     512         ! For the second variable 
    405513         IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 
    406514 
     
    410518 
    411519               IF ( idayend == 0 )  THEN 
    412  
    413520                  ! Daily averaged data 
    414                   CALL obs_int_h2d( kpk, kpk,      & 
    415                      &              zweig2, zinm2(:,:,:,iobs), zobsk ) 
    416  
    417                ENDIF 
    418  
    419             ELSE 
    420  
    421                ! Point data 
    422                CALL obs_int_h2d( kpk, kpk,      & 
    423                   &              zweig2, zint2(:,:,:,iobs), zobsk ) 
    424  
    425             ENDIF 
    426  
    427  
    428             !------------------------------------------------------------- 
    429             ! Compute vertical second-derivative of the interpolating  
    430             ! polynomial at obs points 
    431             !------------------------------------------------------------- 
    432  
    433             IF ( k1dint == 1 ) THEN 
    434                CALL obs_int_z1d_spl( kpk, zobsk, zobs2k, & 
    435                   &                  pgdept, zobsmask2 ) 
    436             ENDIF 
    437  
    438             !---------------------------------------------------------------- 
    439             !  Vertical interpolation to the observation point 
    440             !---------------------------------------------------------------- 
    441             ista = prodatqc%npvsta(jobs,2) 
    442             iend = prodatqc%npvend(jobs,2) 
    443             CALL obs_int_z1d( kpk, & 
    444                & prodatqc%var(2)%mvk(ista:iend),& 
    445                & k1dint, iend - ista + 1, & 
    446                & prodatqc%var(2)%vdep(ista:iend),& 
    447                & zobsk, zobs2k, & 
    448                & prodatqc%var(2)%vmod(ista:iend),& 
    449                & pgdept, zobsmask2 ) 
    450  
    451          ENDIF 
    452  
    453       END DO 
     521 
     522                  ! vertically interpolate all 4 corners  
     523                  ista = prodatqc%npvsta(jobs,2)  
     524                  iend = prodatqc%npvend(jobs,2)  
     525                  inum_obs = iend - ista + 1  
     526                  ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs))  
     527 
     528                  DO iin=1,2  
     529                     DO ijn=1,2  
     530 
     531                        IF ( k1dint == 1 ) THEN  
     532                           CALL obs_int_z1d_spl( kpk, &  
     533                              &     zinm2(iin,ijn,:,iobs), &  
     534                              &     zobs2k, zgdept(iin,ijn,:,iobs), &  
     535                              &     zmask2(iin,ijn,:,iobs))  
     536                        ENDIF  
     537        
     538                        CALL obs_level_search(kpk, &  
     539                           &    zgdept(iin,ijn,:,iobs), &  
     540                           &    inum_obs, prodatqc%var(2)%vdep(ista:iend), &  
     541                           &    iv_indic)  
     542 
     543                        CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &  
     544                           &    prodatqc%var(2)%vdep(ista:iend), &  
     545                           &    zinm2(iin,ijn,:,iobs), &  
     546                           &    zobs2k, interp_corner(iin,ijn,:), &  
     547                           &    zgdept(iin,ijn,:,iobs), &  
     548                           &    zmask2(iin,ijn,:,iobs))  
     549        
     550                     ENDDO  
     551                  ENDDO  
     552 
     553               ENDIF !idayend 
     554 
     555            ELSE    
     556 
     557               ! Point data  
     558      
     559               ! vertically interpolate all 4 corners  
     560               ista = prodatqc%npvsta(jobs,2)  
     561               iend = prodatqc%npvend(jobs,2)  
     562               inum_obs = iend - ista + 1  
     563               ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs))  
     564               DO iin=1,2   
     565                  DO ijn=1,2  
     566                     
     567                     IF ( k1dint == 1 ) THEN  
     568                        CALL obs_int_z1d_spl( kpk, &  
     569                           &    zint2(iin,ijn,:,iobs),&  
     570                           &    zobs2k, zgdept(iin,ijn,:,iobs), &  
     571                           &    zmask2(iin,ijn,:,iobs))  
     572   
     573                     ENDIF  
     574        
     575                     CALL obs_level_search(kpk, &  
     576                         &        zgdept(iin,ijn,:,iobs),&  
     577                         &        inum_obs, prodatqc%var(2)%vdep(ista:iend), &  
     578                         &        iv_indic)  
     579 
     580                     CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs,     &  
     581                         &          prodatqc%var(2)%vdep(ista:iend),     &  
     582                         &          zint2(iin,ijn,:,iobs),            &  
     583                         &          zobs2k,interp_corner(iin,ijn,:), &  
     584                         &          zgdept(iin,ijn,:,iobs),         &  
     585                         &          zmask2(iin,ijn,:,iobs) )       
     586          
     587                  ENDDO  
     588               ENDDO  
     589              
     590            ENDIF  
     591 
     592            !-------------------------------------------------------------  
     593            ! Compute the horizontal interpolation for every profile level  
     594            !-------------------------------------------------------------  
     595              
     596            DO ikn=1,inum_obs  
     597               iend=ista+ikn-1 
     598                   
     599               zweig(:,:,1) = 0._wp  
     600    
     601               ! This code forces the horizontal weights to be   
     602               ! zero IF the observation is below the bottom of the   
     603               ! corners of the interpolation nodes, Or if it is in   
     604               ! the mask. This is important for observations near   
     605               ! steep bathymetry  
     606               DO iin=1,2  
     607                  DO ijn=1,2  
     608      
     609                     depth_loop2: DO ik=kpk,2,-1  
     610                        IF(zmask2(iin,ijn,ik-1,iobs ) > 0.9 )THEN    
     611                             
     612                           zweig(iin,ijn,1) = &   
     613                              & zweig2(iin,ijn,1) * &  
     614                              & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) &  
     615                              &  - prodatqc%var(2)%vdep(iend)),0._wp)  
     616                             
     617                           EXIT depth_loop2  
     618 
     619                        ENDIF  
     620 
     621                     ENDDO depth_loop2  
     622      
     623                  ENDDO  
     624               ENDDO  
     625    
     626               CALL obs_int_h2d( 1, 1, zweig, interp_corner(:,:,ikn), &  
     627                  &              prodatqc%var(2)%vmod(iend:iend) )  
     628 
     629                  ! Set QC flag for any observations found below the bottom 
     630                  ! needed as the check here is more strict than that in obs_prep 
     631               IF (sum(zweig) == 0.0_wp) prodatqc%var(2)%nvqc(iend:iend)=4 
     632  
     633            ENDDO  
     634  
     635            DEALLOCATE(interp_corner,iv_indic)  
     636           
     637         ENDIF  
     638 
     639      ENDDO 
    454640 
    455641      ! Deallocate the data for interpolation 
     
    466652         & zmask2, & 
    467653         & zint1,  & 
    468          & zint2   & 
     654         & zint2,  & 
     655         & zgdept, & 
     656         & zgdepw  & 
    469657         & ) 
    470658 
     
    481669   END SUBROUTINE obs_prof_opt 
    482670 
    483    SUBROUTINE obs_pro_sco_opt( prodatqc, kt, kpi, kpj, kpk, kit000, kdaystp, &  
    484       &                    ptn, psn, pgdept, pgdepw, ptmask, k1dint, k2dint, &  
    485       &                    kdailyavtypes )  
    486       !!-----------------------------------------------------------------------  
    487       !!  
    488       !!                     ***  ROUTINE obs_pro_opt  ***  
    489       !!  
    490       !! ** Purpose : Compute the model counterpart of profiles  
    491       !!              data by interpolating from the model grid to the   
    492       !!              observation point. Generalised vertical coordinate version  
    493       !!  
    494       !! ** Method  : Linearly interpolate to each observation point using   
    495       !!              the model values at the corners of the surrounding grid box.  
    496       !!  
    497       !!          First, model values on the model grid are interpolated vertically to the  
    498       !!          Depths of the profile observations.  Two vertical interpolation schemes are  
    499       !!          available:  
    500       !!          - linear       (k1dint = 0)  
    501       !!          - Cubic spline (k1dint = 1)     
    502       !!  
    503       !!  
    504       !!         Secondly the interpolated values are interpolated horizontally to the   
    505       !!         obs (lon, lat) point.  
    506       !!         Several horizontal interpolation schemes are available:  
    507       !!        - distance-weighted (great circle) (k2dint = 0)  
    508       !!        - distance-weighted (small angle)  (k2dint = 1)  
    509       !!        - bilinear (geographical grid)     (k2dint = 2)  
    510       !!        - bilinear (quadrilateral grid)    (k2dint = 3)  
    511       !!        - polynomial (quadrilateral grid)  (k2dint = 4)  
    512       !!  
    513       !!    For the cubic spline the 2nd derivative of the interpolating   
    514       !!    polynomial is computed before entering the vertical interpolation   
    515       !!    routine.  
    516       !!  
    517       !!    For ENACT moored buoy data (e.g., TAO), the model equivalent is  
    518       !!    a daily mean model temperature field. So, we first compute  
    519       !!    the mean, then interpolate only at the end of the day.  
    520       !!  
    521       !!    This is the procedure to be used with generalised vertical model   
    522       !!    coordinates (ie s-coordinates. It is ~4x slower than the equivalent  
    523       !!    horizontal then vertical interpolation algorithm, but can deal with situations  
    524       !!    where the model levels are not flat.  
    525       !!    ONLY PERFORMED if ln_sco=.TRUE.   
    526       !!        
    527       !!    Note: the in situ temperature observations must be converted  
    528       !!    to potential temperature (the model variable) prior to  
    529       !!    assimilation.   
    530       !!??????????????????????????????????????????????????????????????  
    531       !!    INCLUDE POTENTIAL TEMP -> IN SITU TEMP IN OBS OPERATOR???  
    532       !!??????????????????????????????????????????????????????????????  
    533       !!  
    534       !! ** Action  :  
    535       !!  
    536       !! History :  
    537       !!      ! 2014-08 (J. While) Adapted from obs_pro_opt to handel generalised  
    538       !!                           vertical coordinates 
    539       !!-----------------------------------------------------------------------  
    540     
    541       !! * Modules used  
    542       USE obs_profiles_def   ! Definition of storage space for profile obs.  
    543         
    544       IMPLICIT NONE  
    545   
    546       !! * Arguments  
    547       TYPE(obs_prof), INTENT(INOUT) :: prodatqc   ! Subset of profile data not failing screening  
    548       INTEGER, INTENT(IN) :: kt        ! Time step  
    549       INTEGER, INTENT(IN) :: kpi       ! Model grid parameters  
    550       INTEGER, INTENT(IN) :: kpj  
    551       INTEGER, INTENT(IN) :: kpk  
    552       INTEGER, INTENT(IN) :: kit000    ! Number of the first time step   
    553                                        !   (kit000-1 = restart time)  
    554       INTEGER, INTENT(IN) :: k1dint    ! Vertical interpolation type (see header)  
    555       INTEGER, INTENT(IN) :: k2dint    ! Horizontal interpolation type (see header)  
    556       INTEGER, INTENT(IN) :: kdaystp   ! Number of time steps per day                      
    557       REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &  
    558          & ptn,    &    ! Model temperature field  
    559          & psn,    &    ! Model salinity field  
    560          & ptmask       ! Land-sea mask  
    561       REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &  
    562          & pgdept,  &    ! Model array of depth T levels     
    563          & pgdepw       ! Model array of depth W levels  
    564       INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &  
    565          & kdailyavtypes   ! Types for daily averages  
    566        
    567       !! * Local declarations  
    568       INTEGER ::   ji  
    569       INTEGER ::   jj  
    570       INTEGER ::   jk  
    571       INTEGER ::   iico, ijco  
    572       INTEGER ::   jobs  
    573       INTEGER ::   inrc  
    574       INTEGER ::   ipro  
    575       INTEGER ::   idayend  
    576       INTEGER ::   ista  
    577       INTEGER ::   iend  
    578       INTEGER ::   iobs  
    579       INTEGER ::   iin, ijn, ikn, ik   ! looping indices over interpolation nodes  
    580       INTEGER, DIMENSION(imaxavtypes) :: &  
    581          & idailyavtypes  
    582       INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &  
    583          & igrdi, &  
    584          & igrdj  
    585       INTEGER :: &  
    586          & inum_obs 
    587       INTEGER, ALLOCATABLE, DIMENSION(:) :: iv_indic     
    588       REAL(KIND=wp) :: zlam  
    589       REAL(KIND=wp) :: zphi  
    590       REAL(KIND=wp) :: zdaystp  
    591       REAL(KIND=wp), DIMENSION(kpk) :: &  
    592          & zobsmask, &  
    593          & zobsk,    &  
    594          & zobs2k  
    595       REAL(KIND=wp), DIMENSION(2,2,1) :: &  
    596          & zweig, &  
    597          & l_zweig  
    598       REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &  
    599          & zmask, &  
    600          & zintt, &  
    601          & zints, &  
    602          & zinmt, &  
    603          & zgdept,&  
    604          & zgdepw,&  
    605          & zinms  
    606       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &  
    607          & zglam, &  
    608          & zgphi     
    609       REAL(KIND=wp), DIMENSION(1) :: zmsk_1        
    610       REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner        
    611   
    612       !------------------------------------------------------------------------  
    613       ! Local initialization   
    614       !------------------------------------------------------------------------  
    615       ! ... Record and data counters  
    616       inrc = kt - kit000 + 2  
    617       ipro = prodatqc%npstp(inrc)  
    618    
    619       ! Daily average types  
    620       IF ( PRESENT(kdailyavtypes) ) THEN  
    621          idailyavtypes(:) = kdailyavtypes(:)  
    622       ELSE  
    623          idailyavtypes(:) = -1  
    624       ENDIF  
    625   
    626       ! Initialize daily mean for first time-step  
    627       idayend = MOD( kt - kit000 + 1, kdaystp )  
    628   
    629       ! Added kt == 0 test to catch restart case   
    630       IF ( idayend == 1 .OR. kt == 0) THEN  
    631            
    632          IF (lwp) WRITE(numout,*) 'Reset prodatqc%vdmean on time-step: ',kt  
    633          DO jk = 1, jpk  
    634             DO jj = 1, jpj  
    635                DO ji = 1, jpi  
    636                   prodatqc%vdmean(ji,jj,jk,1) = 0.0  
    637                   prodatqc%vdmean(ji,jj,jk,2) = 0.0  
    638                END DO  
    639             END DO  
    640          END DO  
    641         
    642       ENDIF  
    643         
    644       DO jk = 1, jpk  
    645          DO jj = 1, jpj  
    646             DO ji = 1, jpi  
    647                ! Increment the temperature field for computing daily mean  
    648                prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &  
    649                &                        + ptn(ji,jj,jk)  
    650                ! Increment the salinity field for computing daily mean  
    651                prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &  
    652                &                        + psn(ji,jj,jk)  
    653             END DO  
    654          END DO  
    655       END DO  
    656      
    657       ! Compute the daily mean at the end of day  
    658       zdaystp = 1.0 / REAL( kdaystp )  
    659       IF ( idayend == 0 ) THEN  
    660          DO jk = 1, jpk  
    661             DO jj = 1, jpj  
    662                DO ji = 1, jpi  
    663                   prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) &  
    664                   &                        * zdaystp  
    665                   prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) &  
    666                   &                           * zdaystp  
    667                END DO  
    668             END DO  
    669          END DO  
    670       ENDIF  
    671   
    672       ! Get the data for interpolation  
    673       ALLOCATE( &  
    674          & igrdi(2,2,ipro),      &  
    675          & igrdj(2,2,ipro),      &  
    676          & zglam(2,2,ipro),      &  
    677          & zgphi(2,2,ipro),      &  
    678          & zmask(2,2,kpk,ipro),  &  
    679          & zintt(2,2,kpk,ipro),  &  
    680          & zints(2,2,kpk,ipro),  &  
    681          & zgdept(2,2,kpk,ipro), &  
    682          & zgdepw(2,2,kpk,ipro)  &  
    683          & )  
    684   
    685       DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro  
    686          iobs = jobs - prodatqc%nprofup  
    687          igrdi(1,1,iobs) = prodatqc%mi(jobs,1)-1  
    688          igrdj(1,1,iobs) = prodatqc%mj(jobs,1)-1  
    689          igrdi(1,2,iobs) = prodatqc%mi(jobs,1)-1  
    690          igrdj(1,2,iobs) = prodatqc%mj(jobs,1)  
    691          igrdi(2,1,iobs) = prodatqc%mi(jobs,1)  
    692          igrdj(2,1,iobs) = prodatqc%mj(jobs,1)-1  
    693          igrdi(2,2,iobs) = prodatqc%mi(jobs,1)  
    694          igrdj(2,2,iobs) = prodatqc%mj(jobs,1)  
    695       END DO  
    696       
    697       ! Initialise depth arrays 
    698       zgdept = 0.0 
    699       zgdepw = 0.0 
    700   
    701       CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi, igrdj, glamt, zglam )  
    702       CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi, igrdj, gphit, zgphi )  
    703       CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, ptmask,zmask )  
    704       CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, ptn,   zintt )  
    705       CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, psn,   zints )  
    706       CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pgdept(:,:,:), &  
    707         &                     zgdept )  
    708       CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pgdepw(:,:,:), &  
    709         &                     zgdepw )  
    710   
    711       ! At the end of the day also get interpolated means  
    712       IF ( idayend == 0 ) THEN  
    713   
    714          ALLOCATE( &  
    715             & zinmt(2,2,kpk,ipro),  &  
    716             & zinms(2,2,kpk,ipro)   &  
    717             & )  
    718   
    719          CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, &  
    720             &                  prodatqc%vdmean(:,:,:,1), zinmt )  
    721          CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, &  
    722             &                  prodatqc%vdmean(:,:,:,2), zinms )  
    723   
    724       ENDIF  
    725         
    726       ! Return if no observations to process  
    727       ! Has to be done after comm commands to ensure processors  
    728       ! stay in sync  
    729       IF ( ipro == 0 ) RETURN  
    730   
    731       DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro  
    732      
    733          iobs = jobs - prodatqc%nprofup  
    734      
    735          IF ( kt /= prodatqc%mstp(jobs) ) THEN  
    736               
    737             IF(lwp) THEN  
    738                WRITE(numout,*)  
    739                WRITE(numout,*) ' E R R O R : Observation',              &  
    740                   &            ' time step is not consistent with the', &  
    741                   &            ' model time step'  
    742                WRITE(numout,*) ' ========='  
    743                WRITE(numout,*)  
    744                WRITE(numout,*) ' Record  = ', jobs,                    &  
    745                   &            ' kt      = ', kt,                      &  
    746                   &            ' mstp    = ', prodatqc%mstp(jobs), &  
    747                   &            ' ntyp    = ', prodatqc%ntyp(jobs)  
    748             ENDIF  
    749             CALL ctl_stop( 'obs_pro_opt', 'Inconsistent time' )  
    750          ENDIF  
    751            
    752          zlam = prodatqc%rlam(jobs)  
    753          zphi = prodatqc%rphi(jobs)  
    754            
    755          ! Horizontal weights  
    756          ! Only calculated once, for both T and S.  
    757          ! Masked values are calculated later.   
    758   
    759          IF ( ( prodatqc%npvend(jobs,1) > 0 ) .OR. &  
    760             & ( prodatqc%npvend(jobs,2) > 0 ) ) THEN  
    761   
    762             CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,     &  
    763                &                   zglam(:,:,iobs), zgphi(:,:,iobs), &  
    764                &                   zmask(:,:,1,iobs), zweig, zmsk_1 )  
    765   
    766          ENDIF  
    767           
    768          ! IF zmsk_1 = 0; then ob is on land  
    769          IF (zmsk_1(1) < 0.1) THEN  
    770             WRITE(numout,*) 'WARNING (obs_oper) :- profile found within landmask'  
    771     
    772          ELSE   
    773               
    774             ! Temperature  
    775               
    776             IF ( prodatqc%npvend(jobs,1) > 0 ) THEN   
    777      
    778                zobsk(:) = obfillflt  
    779      
    780                IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN  
    781      
    782                   IF ( idayend == 0 )  THEN  
    783                     
    784                      ! Daily averaged moored buoy (MRB) data  
    785                     
    786                      ! vertically interpolate all 4 corners  
    787                      ista = prodatqc%npvsta(jobs,1)  
    788                      iend = prodatqc%npvend(jobs,1)  
    789                      inum_obs = iend - ista + 1  
    790                      ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs))  
    791        
    792                      DO iin=1,2  
    793                         DO ijn=1,2  
    794                                         
    795                                         
    796             
    797                            IF ( k1dint == 1 ) THEN  
    798                               CALL obs_int_z1d_spl( kpk, &  
    799                                  &     zinmt(iin,ijn,:,iobs), &  
    800                                  &     zobs2k, zgdept(iin,ijn,:,iobs), &  
    801                                  &     zmask(iin,ijn,:,iobs))  
    802                            ENDIF  
    803         
    804                            CALL obs_level_search(kpk, &  
    805                               &    zgdept(iin,ijn,:,iobs), &  
    806                               &    inum_obs, prodatqc%var(1)%vdep(ista:iend), &  
    807                               &    iv_indic)  
    808                            CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &  
    809                               &    prodatqc%var(1)%vdep(ista:iend), &  
    810                               &    zinmt(iin,ijn,:,iobs), &  
    811                               &    zobs2k, interp_corner(iin,ijn,:), &  
    812                               &    zgdept(iin,ijn,:,iobs), &  
    813                               &    zmask(iin,ijn,:,iobs))  
    814         
    815                         ENDDO  
    816                      ENDDO  
    817                     
    818                     
    819                   ELSE  
    820                  
    821                      CALL ctl_stop( ' A nonzero' //     &  
    822                         &           ' number of profile T BUOY data should' // &  
    823                         &           ' only occur at the end of a given day' )  
    824      
    825                   ENDIF  
    826           
    827                ELSE   
    828                  
    829                   ! Point data  
    830       
    831                   ! vertically interpolate all 4 corners  
    832                   ista = prodatqc%npvsta(jobs,1)  
    833                   iend = prodatqc%npvend(jobs,1)  
    834                   inum_obs = iend - ista + 1  
    835                   ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs))  
    836                   DO iin=1,2   
    837                      DO ijn=1,2  
    838                                      
    839                                      
    840                         IF ( k1dint == 1 ) THEN  
    841                            CALL obs_int_z1d_spl( kpk, &  
    842                               &    zintt(iin,ijn,:,iobs),&  
    843                               &    zobs2k, zgdept(iin,ijn,:,iobs), &  
    844                               &    zmask(iin,ijn,:,iobs))  
    845    
    846                         ENDIF  
    847         
    848                         CALL obs_level_search(kpk, &  
    849                             &        zgdept(iin,ijn,:,iobs),&  
    850                             &        inum_obs, prodatqc%var(1)%vdep(ista:iend), &  
    851                             &         iv_indic)  
    852                         CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs,     &  
    853                             &          prodatqc%var(1)%vdep(ista:iend),     &  
    854                             &          zintt(iin,ijn,:,iobs),            &  
    855                             &          zobs2k,interp_corner(iin,ijn,:), &  
    856                             &          zgdept(iin,ijn,:,iobs),         &  
    857                             &          zmask(iin,ijn,:,iobs) )       
    858           
    859                      ENDDO  
    860                   ENDDO  
    861               
    862                ENDIF  
    863         
    864                !-------------------------------------------------------------  
    865                ! Compute the horizontal interpolation for every profile level  
    866                !-------------------------------------------------------------  
    867               
    868                DO ikn=1,inum_obs  
    869                   iend=ista+ikn-1  
    870  
    871                   l_zweig(:,:,1) = 0._wp  
    872  
    873                   ! This code forces the horizontal weights to be   
    874                   ! zero IF the observation is below the bottom of the   
    875                   ! corners of the interpolation nodes, Or if it is in   
    876                   ! the mask. This is important for observations are near   
    877                   ! steep bathymetry  
    878                   DO iin=1,2  
    879                      DO ijn=1,2  
    880       
    881                         depth_loop1: DO ik=kpk,2,-1  
    882                            IF(zmask(iin,ijn,ik-1,iobs ) > 0.9 )THEN    
    883                              
    884                               l_zweig(iin,ijn,1) = &   
    885                                  & zweig(iin,ijn,1) * &  
    886                                  & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) &  
    887                                  &  - prodatqc%var(1)%vdep(iend)),0._wp)  
    888                              
    889                               EXIT depth_loop1  
    890                            ENDIF  
    891                         ENDDO depth_loop1  
    892       
    893                      ENDDO  
    894                   ENDDO  
    895     
    896                   CALL obs_int_h2d( 1, 1, l_zweig, interp_corner(:,:,ikn), &  
    897                   &          prodatqc%var(1)%vmod(iend:iend) )  
    898   
    899                ENDDO  
    900   
    901   
    902                DEALLOCATE(interp_corner,iv_indic)  
    903            
    904             ENDIF  
    905         
    906   
    907             ! Salinity   
    908            
    909             IF ( prodatqc%npvend(jobs,2) > 0 ) THEN   
    910      
    911                zobsk(:) = obfillflt  
    912      
    913                IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN  
    914      
    915                   IF ( idayend == 0 )  THEN  
    916                     
    917                      ! Daily averaged moored buoy (MRB) data  
    918                     
    919                      ! vertically interpolate all 4 corners  
    920                      ista = prodatqc%npvsta(iobs,2)  
    921                      iend = prodatqc%npvend(iobs,2)  
    922                      inum_obs = iend - ista + 1  
    923                      ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs))  
    924        
    925                      DO iin=1,2  
    926                         DO ijn=1,2  
    927                                         
    928                                         
    929             
    930                            IF ( k1dint == 1 ) THEN  
    931                               CALL obs_int_z1d_spl( kpk, &  
    932                                  &     zinms(iin,ijn,:,iobs), &  
    933                                  &     zobs2k, zgdept(iin,ijn,:,iobs), &  
    934                                  &     zmask(iin,ijn,:,iobs))  
    935                            ENDIF  
    936         
    937                            CALL obs_level_search(kpk, &  
    938                               &    zgdept(iin,ijn,:,iobs), &  
    939                               &    inum_obs, prodatqc%var(2)%vdep(ista:iend), &  
    940                               &    iv_indic)  
    941                            CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &  
    942                               &    prodatqc%var(2)%vdep(ista:iend), &  
    943                               &    zinms(iin,ijn,:,iobs), &  
    944                               &    zobs2k, interp_corner(iin,ijn,:), &  
    945                               &    zgdept(iin,ijn,:,iobs), &  
    946                               &    zmask(iin,ijn,:,iobs))  
    947         
    948                         ENDDO  
    949                      ENDDO  
    950                     
    951                     
    952                   ELSE  
    953                  
    954                      CALL ctl_stop( ' A nonzero' //     &  
    955                         &           ' number of profile T BUOY data should' // &  
    956                         &           ' only occur at the end of a given day' )  
    957      
    958                   ENDIF  
    959           
    960                ELSE   
    961                  
    962                   ! Point data  
    963       
    964                   ! vertically interpolate all 4 corners  
    965                   ista = prodatqc%npvsta(jobs,2)  
    966                   iend = prodatqc%npvend(jobs,2)  
    967                   inum_obs = iend - ista + 1  
    968                   ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs))  
    969                     
    970                   DO iin=1,2      
    971                      DO ijn=1,2   
    972                                   
    973                                   
    974                         IF ( k1dint == 1 ) THEN  
    975                            CALL obs_int_z1d_spl( kpk, &  
    976                               &    zints(iin,ijn,:,iobs),&  
    977                               &    zobs2k, zgdept(iin,ijn,:,iobs), &  
    978                               &    zmask(iin,ijn,:,iobs))  
    979    
    980                         ENDIF  
    981         
    982                         CALL obs_level_search(kpk, &  
    983                            &        zgdept(iin,ijn,:,iobs),&  
    984                            &        inum_obs, prodatqc%var(2)%vdep(ista:iend), &  
    985                            &         iv_indic)  
    986                         CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs,  &  
    987                            &          prodatqc%var(2)%vdep(ista:iend),     &  
    988                            &          zints(iin,ijn,:,iobs),               &  
    989                            &          zobs2k,interp_corner(iin,ijn,:),     &  
    990                            &          zgdept(iin,ijn,:,iobs),              &  
    991                            &          zmask(iin,ijn,:,iobs) )       
    992           
    993                      ENDDO  
    994                   ENDDO  
    995               
    996                ENDIF  
    997         
    998                !-------------------------------------------------------------  
    999                ! Compute the horizontal interpolation for every profile level  
    1000                !-------------------------------------------------------------  
    1001               
    1002                DO ikn=1,inum_obs  
    1003                   iend=ista+ikn-1  
    1004  
    1005                   l_zweig(:,:,1) = 0._wp 
    1006     
    1007                   ! This code forces the horizontal weights to be   
    1008                   ! zero IF the observation is below the bottom of the   
    1009                   ! corners of the interpolation nodes, Or if it is in   
    1010                   ! the mask. This is important for observations are near   
    1011                   ! steep bathymetry  
    1012                   DO iin=1,2  
    1013                      DO ijn=1,2  
    1014       
    1015                         depth_loop2: DO ik=kpk,2,-1  
    1016                            IF(zmask(iin,ijn,ik-1,iobs ) > 0.9 )THEN    
    1017                              
    1018                               l_zweig(iin,ijn,1) = &   
    1019                                  &  zweig(iin,ijn,1) * &  
    1020                                  &  MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) &  
    1021                                  &  - prodatqc%var(2)%vdep(iend)),0._wp)  
    1022                              
    1023                               EXIT depth_loop2  
    1024                            ENDIF  
    1025                         ENDDO depth_loop2  
    1026       
    1027                      ENDDO  
    1028                   ENDDO  
    1029     
    1030                   CALL obs_int_h2d( 1, 1, l_zweig, interp_corner(:,:,ikn), &  
    1031                   &          prodatqc%var(2)%vmod(iend:iend) )  
    1032   
    1033                ENDDO  
    1034   
    1035   
    1036                DEALLOCATE(interp_corner,iv_indic)  
    1037            
    1038             ENDIF  
    1039            
    1040          ENDIF  
    1041         
    1042       END DO  
    1043       
    1044       ! Deallocate the data for interpolation  
    1045       DEALLOCATE( &  
    1046          & igrdi, &  
    1047          & igrdj, &  
    1048          & zglam, &  
    1049          & zgphi, &  
    1050          & zmask, &  
    1051          & zintt, &  
    1052          & zints, &  
    1053          & zgdept,& 
    1054          & zgdepw & 
    1055          & )  
    1056       ! At the end of the day also get interpolated means  
    1057       IF ( idayend == 0 ) THEN  
    1058          DEALLOCATE( &  
    1059             & zinmt,  &  
    1060             & zinms   &  
    1061             & )  
    1062       ENDIF  
    1063      
    1064       prodatqc%nprofup = prodatqc%nprofup + ipro   
    1065         
    1066    END SUBROUTINE obs_pro_sco_opt  
    1067   
    1068    SUBROUTINE obs_surf_opt( surfdataqc, kt, kpi, kpj,         & 
    1069       &                    kit000, kdaystp, psurf, psurfmask, & 
    1070       &                    k2dint, ldnightav ) 
     671   SUBROUTINE obs_surf_opt( surfdataqc, kt, kpi, kpj,            & 
     672      &                     kit000, kdaystp, psurf, psurfmask,   & 
     673      &                     k2dint, ldnightav, plamscl, pphiscl, & 
     674      &                     lindegrees ) 
    1071675 
    1072676      !!----------------------------------------------------------------------- 
     
    1090694      !!        - polynomial (quadrilateral grid)  (k2dint = 4) 
    1091695      !! 
     696      !!    Two horizontal averaging schemes are also available: 
     697      !!        - weighted radial footprint        (k2dint = 5) 
     698      !!        - weighted rectangular footprint   (k2dint = 6) 
     699      !! 
    1092700      !! 
    1093701      !! ** Action  : 
     
    1096704      !!      ! 07-03 (A. Weaver) 
    1097705      !!      ! 15-02 (M. Martin) Combined routine for surface types 
     706      !!      ! 17-03 (M. Martin) Added horizontal averaging options 
    1098707      !!----------------------------------------------------------------------- 
    1099708 
     
    1117726         & psurfmask                   ! Land-sea mask 
    1118727      LOGICAL, INTENT(IN) :: ldnightav ! Logical for averaging night-time data 
     728      REAL(KIND=wp), INTENT(IN) :: & 
     729         & plamscl, &                  ! Diameter in metres of obs footprint in E/W, N/S directions 
     730         & pphiscl                     ! This is the full width (rather than half-width) 
     731      LOGICAL, INTENT(IN) :: & 
     732         & lindegrees                  ! T=> plamscl and pphiscl are specified in degrees, F=> in metres 
    1119733 
    1120734      !! * Local declarations 
     
    1125739      INTEGER :: isurf 
    1126740      INTEGER :: iobs 
     741      INTEGER :: imaxifp, imaxjfp 
     742      INTEGER :: imodi, imodj 
    1127743      INTEGER :: idayend 
    1128744      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 
    1129          & igrdi, & 
    1130          & igrdj 
     745         & igrdi,   & 
     746         & igrdj,   & 
     747         & igrdip1, & 
     748         & igrdjp1 
    1131749      INTEGER, DIMENSION(:,:), SAVE, ALLOCATABLE :: & 
    1132750         & icount_night,      & 
     
    1136754      REAL(wp), DIMENSION(1) :: zext, zobsmask 
    1137755      REAL(wp) :: zdaystp 
    1138       REAL(wp), DIMENSION(2,2,1) :: & 
    1139          & zweig 
    1140756      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 
     757         & zweig,  & 
    1141758         & zmask,  & 
    1142759         & zsurf,  & 
    1143760         & zsurfm, & 
     761         & zsurftmp, & 
    1144762         & zglam,  & 
    1145          & zgphi 
     763         & zgphi,  & 
     764         & zglamf, & 
     765         & zgphif 
     766 
    1146767      REAL(wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: & 
    1147768         & zintmp,  & 
     
    1155776      inrc = kt - kit000 + 2 
    1156777      isurf = surfdataqc%nsstp(inrc) 
     778 
     779      ! Work out the maximum footprint size for the  
     780      ! interpolation/averaging in model grid-points - has to be even. 
     781 
     782      CALL obs_max_fpsize( k2dint, plamscl, pphiscl, lindegrees, psurfmask, imaxifp, imaxjfp ) 
     783 
    1157784 
    1158785      IF ( ldnightav ) THEN 
     
    1221848 
    1222849      ALLOCATE( & 
    1223          & igrdi(2,2,isurf), & 
    1224          & igrdj(2,2,isurf), & 
    1225          & zglam(2,2,isurf), & 
    1226          & zgphi(2,2,isurf), & 
    1227          & zmask(2,2,isurf), & 
    1228          & zsurf(2,2,isurf)  & 
     850         & zweig(imaxifp,imaxjfp,1),      & 
     851         & igrdi(imaxifp,imaxjfp,isurf), & 
     852         & igrdj(imaxifp,imaxjfp,isurf), & 
     853         & zglam(imaxifp,imaxjfp,isurf), & 
     854         & zgphi(imaxifp,imaxjfp,isurf), & 
     855         & zmask(imaxifp,imaxjfp,isurf), & 
     856         & zsurf(imaxifp,imaxjfp,isurf), & 
     857         & zsurftmp(imaxifp,imaxjfp,isurf),  & 
     858         & zglamf(imaxifp+1,imaxjfp+1,isurf), & 
     859         & zgphif(imaxifp+1,imaxjfp+1,isurf), & 
     860         & igrdip1(imaxifp+1,imaxjfp+1,isurf), & 
     861         & igrdjp1(imaxifp+1,imaxjfp+1,isurf) & 
    1229862         & ) 
    1230863 
    1231864      DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf 
    1232865         iobs = jobs - surfdataqc%nsurfup 
    1233          igrdi(1,1,iobs) = surfdataqc%mi(jobs)-1 
    1234          igrdj(1,1,iobs) = surfdataqc%mj(jobs)-1 
    1235          igrdi(1,2,iobs) = surfdataqc%mi(jobs)-1 
    1236          igrdj(1,2,iobs) = surfdataqc%mj(jobs) 
    1237          igrdi(2,1,iobs) = surfdataqc%mi(jobs) 
    1238          igrdj(2,1,iobs) = surfdataqc%mj(jobs)-1 
    1239          igrdi(2,2,iobs) = surfdataqc%mi(jobs) 
    1240          igrdj(2,2,iobs) = surfdataqc%mj(jobs) 
     866         DO ji = 0, imaxifp 
     867            imodi = surfdataqc%mi(jobs) - int(imaxifp/2) + ji - 1 
     868             
     869            !Deal with wrap around in longitude 
     870            IF ( imodi < 1      ) imodi = imodi + jpiglo 
     871            IF ( imodi > jpiglo ) imodi = imodi - jpiglo 
     872             
     873            DO jj = 0, imaxjfp 
     874               imodj = surfdataqc%mj(jobs) - int(imaxjfp/2) + jj - 1 
     875               !If model values are out of the domain to the north/south then 
     876               !set them to be the edge of the domain 
     877               IF ( imodj < 1      ) imodj = 1 
     878               IF ( imodj > jpjglo ) imodj = jpjglo 
     879 
     880               igrdip1(ji+1,jj+1,iobs) = imodi 
     881               igrdjp1(ji+1,jj+1,iobs) = imodj 
     882                
     883               IF ( ji >= 1 .AND. jj >= 1 ) THEN 
     884                  igrdi(ji,jj,iobs) = imodi 
     885                  igrdj(ji,jj,iobs) = imodj 
     886               ENDIF 
     887                
     888            END DO 
     889         END DO 
    1241890      END DO 
    1242891 
    1243       CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, & 
     892      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 
    1244893         &                  igrdi, igrdj, glamt, zglam ) 
    1245       CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, & 
     894      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 
    1246895         &                  igrdi, igrdj, gphit, zgphi ) 
    1247       CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, & 
     896      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 
    1248897         &                  igrdi, igrdj, psurfmask, zmask ) 
    1249       CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, & 
     898      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 
    1250899         &                  igrdi, igrdj, psurf, zsurf ) 
     900      CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, & 
     901         &                  igrdip1, igrdjp1, glamf, zglamf ) 
     902      CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, & 
     903         &                  igrdip1, igrdjp1, gphif, zgphif ) 
    1251904 
    1252905      ! At the end of the day get interpolated means 
    1253       IF (ldnightav ) THEN 
    1254          IF ( idayend == 0 ) THEN 
    1255  
    1256             ALLOCATE( & 
    1257                & zsurfm(2,2,isurf)  & 
    1258                & ) 
    1259  
    1260             CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, igrdi, igrdj, & 
    1261                &               surfdataqc%vdmean(:,:), zsurfm ) 
    1262  
    1263          ENDIF 
     906      IF ( idayend == 0 .AND. ldnightav ) THEN 
     907 
     908         ALLOCATE( & 
     909            & zsurfm(imaxifp,imaxjfp,isurf)  & 
     910            & ) 
     911 
     912         CALL obs_int_comm_2d( imaxifp,imaxjfp, isurf, kpi, kpj, igrdi, igrdj, & 
     913            &               surfdataqc%vdmean(:,:), zsurfm ) 
     914 
    1264915      ENDIF 
    1265916 
     
    1290941         zphi = surfdataqc%rphi(jobs) 
    1291942 
    1292          ! Get weights to interpolate the model value to the observation point 
    1293          CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         & 
    1294             &                   zglam(:,:,iobs), zgphi(:,:,iobs), & 
    1295             &                   zmask(:,:,iobs), zweig, zobsmask ) 
    1296  
    1297          ! Interpolate the model field to the observation point 
    1298943         IF ( ldnightav .AND. idayend == 0 ) THEN 
    1299944            ! Night-time averaged data 
    1300             CALL obs_int_h2d( 1, 1, zweig, zsurfm(:,:,iobs), zext ) 
     945            zsurftmp(:,:,iobs) = zsurfm(:,:,iobs) 
    1301946         ELSE 
    1302             CALL obs_int_h2d( 1, 1, zweig, zsurf(:,:,iobs),  zext ) 
     947            zsurftmp(:,:,iobs) = zsurf(:,:,iobs) 
     948         ENDIF 
     949 
     950         IF ( k2dint <= 4 ) THEN 
     951 
     952            ! Get weights to interpolate the model value to the observation point 
     953            CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         & 
     954               &                   zglam(:,:,iobs), zgphi(:,:,iobs), & 
     955               &                   zmask(:,:,iobs), zweig, zobsmask ) 
     956 
     957            ! Interpolate the model value to the observation point  
     958            CALL obs_int_h2d( 1, 1, zweig, zsurftmp(:,:,iobs), zext ) 
     959 
     960         ELSE 
     961 
     962            ! Get weights to average the model SLA to the observation footprint 
     963            CALL obs_avg_h2d_init( 1, 1, imaxifp, imaxjfp, k2dint, zlam,  zphi, & 
     964               &                   zglam(:,:,iobs), zgphi(:,:,iobs), & 
     965               &                   zglamf(:,:,iobs), zgphif(:,:,iobs), & 
     966               &                   zmask(:,:,iobs), plamscl, pphiscl, & 
     967               &                   lindegrees, zweig, zobsmask ) 
     968 
     969            ! Average the model SST to the observation footprint 
     970            CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, & 
     971               &              zweig, zsurftmp(:,:,iobs),  zext ) 
     972 
    1303973         ENDIF 
    1304974 
     
    1310980            surfdataqc%rmod(jobs,1) = zext(1) 
    1311981         ENDIF 
     982          
     983         IF ( zext(1) == obfillflt ) THEN 
     984            ! If the observation value is a fill value, set QC flag to bad 
     985            surfdataqc%nqc(jobs) = 4 
     986         ENDIF 
    1312987 
    1313988      END DO 
     
    1315990      ! Deallocate the data for interpolation 
    1316991      DEALLOCATE( & 
     992         & zweig, & 
    1317993         & igrdi, & 
    1318994         & igrdj, & 
     
    1320996         & zgphi, & 
    1321997         & zmask, & 
    1322          & zsurf  & 
     998         & zsurf, & 
     999         & zsurftmp, & 
     1000         & zglamf, & 
     1001         & zgphif, & 
     1002         & igrdip1,& 
     1003         & igrdjp1 & 
    13231004         & ) 
    13241005 
    13251006      ! At the end of the day also deallocate night-time mean array 
    1326       IF ( ldnightav ) THEN 
    1327          IF ( idayend == 0 ) THEN 
    1328             DEALLOCATE( & 
    1329                & zsurfm  & 
    1330                & ) 
    1331          ENDIF 
     1007      IF ( idayend == 0 .AND. ldnightav ) THEN 
     1008         DEALLOCATE( & 
     1009            & zsurfm  & 
     1010            & ) 
    13321011      ENDIF 
    13331012 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90

    r7646 r9023  
    2323   USE obs_oper           ! Observation operators 
    2424   USE lib_mpp, ONLY :   ctl_warn, ctl_stop 
     25   USE bdy_oce, ONLY : &        ! Boundary information 
     26      idx_bdy, nb_bdy, ln_bdy 
    2527 
    2628   IMPLICIT NONE 
     
    4042CONTAINS 
    4143 
    42    SUBROUTINE obs_pre_surf( surfdata, surfdataqc, ld_nea ) 
     44   SUBROUTINE obs_pre_surf( surfdata, surfdataqc, ld_nea, ld_bound_reject, & 
     45                            kqc_cutoff ) 
    4346      !!---------------------------------------------------------------------- 
    4447      !!                    ***  ROUTINE obs_pre_sla  *** 
     
    5760      !!        !  2015-02  (M. Martin) Combined routine for surface types. 
    5861      !!---------------------------------------------------------------------- 
     62      !! * Modules used 
    5963      USE par_oce             ! Ocean parameters 
    6064      USE dom_oce, ONLY       :   glamt, gphit, tmask, nproc   ! Geographical information 
     
    6367      TYPE(obs_surf), INTENT(INOUT) :: surfdataqc   ! Subset of surface data not failing screening 
    6468      LOGICAL, INTENT(IN) :: ld_nea         ! Switch for rejecting observation near land 
    65       ! 
     69      LOGICAL, INTENT(IN) :: ld_bound_reject       ! Switch for rejecting obs near the boundary 
     70      INTEGER, INTENT(IN), OPTIONAL :: kqc_cutoff   ! cut off for QC value 
     71      !! * Local declarations 
     72      INTEGER :: iqc_cutoff = 255   ! cut off for QC value 
    6673      INTEGER :: iyea0        ! Initial date 
    6774      INTEGER :: imon0        !  - (year, month, day, hour, minute) 
     
    7683      INTEGER :: inlasobs     !  - close to land 
    7784      INTEGER :: igrdobs      !  - fail the grid search 
     85      INTEGER :: ibdysobs     !  - close to open boundary 
    7886                              ! Global counters for observations that 
    7987      INTEGER :: iotdobsmpp     !  - outside time domain 
     
    8290      INTEGER :: inlasobsmpp    !  - close to land 
    8391      INTEGER :: igrdobsmpp     !  - fail the grid search 
    84       LOGICAL, DIMENSION(:), ALLOCATABLE ::   llvalid            ! SLA data selection 
     92      INTEGER :: ibdysobsmpp  !  - close to open boundary 
     93      LOGICAL, DIMENSION(:), ALLOCATABLE :: &  
     94         & llvalid            ! SLA data selection 
    8595      INTEGER :: jobs         ! Obs. loop variable 
    8696      INTEGER :: jstp         ! Time loop variable 
     
    107117      ilansobs = 0 
    108118      inlasobs = 0 
     119      ibdysobs = 0  
     120 
     121      ! Set QC cutoff to optional value if provided 
     122      IF ( PRESENT(kqc_cutoff) ) iqc_cutoff=kqc_cutoff 
    109123 
    110124      ! ----------------------------------------------------------------------- 
     
    140154         &                 tmask(:,:,1), surfdata%nqc,  & 
    141155         &                 iosdsobs,     ilansobs,     & 
    142          &                 inlasobs,     ld_nea        ) 
     156         &                 inlasobs,     ld_nea,       & 
     157         &                 ibdysobs,     ld_bound_reject, & 
     158         &                 iqc_cutoff                     ) 
    143159 
    144160      CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp ) 
    145161      CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp ) 
    146162      CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp ) 
     163      CALL obs_mpp_sum_integer( ibdysobs, ibdysobsmpp ) 
    147164 
    148165      ! ----------------------------------------------------------------------- 
     
    155172      ALLOCATE( llvalid(surfdata%nsurf) ) 
    156173       
    157       ! We want all data which has qc flags <= 10 
    158  
    159       llvalid(:)  = ( surfdata%nqc(:)  <= 10 ) 
     174      ! We want all data which has qc flags <= iqc_cutoff 
     175 
     176      llvalid(:)  = ( surfdata%nqc(:)  <= iqc_cutoff ) 
    160177 
    161178      ! The actual copying 
     
    190207               &            inlasobsmpp 
    191208         ENDIF 
     209         WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data near open boundary (removed) = ', & 
     210            &            ibdysobsmpp   
    192211         WRITE(numout,*) ' '//surfdataqc%cvars(1)//' data accepted                             = ', & 
    193212            &            surfdataqc%nsurfmpp 
     
    225244      &                     kpi, kpj, kpk, & 
    226245      &                     zmask1, pglam1, pgphi1, zmask2, pglam2, pgphi2,  & 
    227       &                     ld_nea, kdailyavtypes ) 
     246      &                     ld_nea, ld_bound_reject, kdailyavtypes,  kqc_cutoff ) 
    228247 
    229248!!---------------------------------------------------------------------- 
     
    241260      !! 
    242261      !!---------------------------------------------------------------------- 
    243       USE par_oce           ! Ocean parameters 
    244       USE dom_oce, ONLY :   gdept_1d, nproc   ! Geographical information 
     262      !! * Modules used 
     263      USE par_oce             ! Ocean parameters 
     264      USE dom_oce, ONLY : &   ! Geographical information 
     265         & gdept_1d,             & 
     266         & nproc 
    245267 
    246268      !! * Arguments 
     
    250272      LOGICAL, INTENT(IN) :: ld_var2 
    251273      LOGICAL, INTENT(IN) :: ld_nea               ! Switch for rejecting observation near land 
     274      LOGICAL, INTENT(IN) :: ld_bound_reject      ! Switch for rejecting observations near the boundary 
    252275      INTEGER, INTENT(IN) :: kpi, kpj, kpk        ! Local domain sizes 
    253276      INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 
     
    261284         & pgphi1, & 
    262285         & pgphi2 
     286      INTEGER, INTENT(IN), OPTIONAL :: kqc_cutoff   ! cut off for QC value 
    263287 
    264288      !! * Local declarations 
     289      INTEGER :: iqc_cutoff = 255   ! cut off for QC value 
    265290      INTEGER :: iyea0        ! Initial date 
    266291      INTEGER :: imon0        !  - (year, month, day, hour, minute) 
     
    277302      INTEGER :: inlav1obs    !  - close to land (variable 1) 
    278303      INTEGER :: inlav2obs    !  - close to land (variable 2) 
     304      INTEGER :: ibdyv1obs    !  - boundary (variable 1)  
     305      INTEGER :: ibdyv2obs    !  - boundary (variable 2)       
    279306      INTEGER :: igrdobs      !  - fail the grid search 
    280307      INTEGER :: iuvchku      !  - reject u if v rejected and vice versa 
     
    288315      INTEGER :: inlav1obsmpp !  - close to land (variable 1) 
    289316      INTEGER :: inlav2obsmpp !  - close to land (variable 2) 
     317      INTEGER :: ibdyv1obsmpp !  - boundary (variable 1)  
     318      INTEGER :: ibdyv2obsmpp !  - boundary (variable 2)       
    290319      INTEGER :: igrdobsmpp   !  - fail the grid search 
    291320      INTEGER :: iuvchkumpp   !  - reject var1 if var2 rejected and vice versa 
     
    322351      inlav1obs = 0 
    323352      inlav2obs = 0 
    324       iuvchku  = 0 
    325       iuvchkv = 0 
     353      ibdyv1obs = 0 
     354      ibdyv2obs = 0 
     355      iuvchku   = 0 
     356      iuvchkv   = 0 
     357 
     358 
     359      ! Set QC cutoff to optional value if provided 
     360      IF ( PRESENT(kqc_cutoff) ) iqc_cutoff=kqc_cutoff 
    326361 
    327362      ! ----------------------------------------------------------------------- 
     
    335370            &              profdata%nday,    profdata%nhou, profdata%nmin, & 
    336371            &              profdata%ntyp,    profdata%nqc,  profdata%mstp, & 
    337             &              iotdobs, kdailyavtypes = kdailyavtypes ) 
     372            &              iotdobs, kdailyavtypes = kdailyavtypes,         & 
     373            &              kqc_cutoff = iqc_cutoff ) 
    338374      ELSE 
    339375         CALL obs_coo_tim_prof( icycle, & 
     
    342378            &              profdata%nday,    profdata%nhou, profdata%nmin, & 
    343379            &              profdata%ntyp,    profdata%nqc,  profdata%mstp, & 
    344             &              iotdobs ) 
     380            &              iotdobs,          kqc_cutoff = iqc_cutoff ) 
    345381      ENDIF 
    346382 
     
    359395 
    360396      ! ----------------------------------------------------------------------- 
    361       ! Reject all observations for profiles with nqc > 10 
    362       ! ----------------------------------------------------------------------- 
    363  
    364       CALL obs_pro_rej( profdata ) 
     397      ! Reject all observations for profiles with nqc > iqc_cutoff 
     398      ! ----------------------------------------------------------------------- 
     399 
     400      CALL obs_pro_rej( profdata, kqc_cutoff = iqc_cutoff ) 
    365401 
    366402      ! ----------------------------------------------------------------------- 
     
    381417         &                 gdept_1d,              zmask1,               & 
    382418         &                 profdata%nqc,          profdata%var(1)%nvqc, & 
    383          &                 iosdv1obs,              ilanv1obs,           & 
    384          &                 inlav1obs,              ld_nea                ) 
     419         &                 iosdv1obs,             ilanv1obs,            & 
     420         &                 inlav1obs,             ld_nea,               & 
     421         &                 ibdyv1obs,             ld_bound_reject,      & 
     422         &                 iqc_cutoff       ) 
    385423 
    386424      CALL obs_mpp_sum_integer( iosdv1obs, iosdv1obsmpp ) 
    387425      CALL obs_mpp_sum_integer( ilanv1obs, ilanv1obsmpp ) 
    388426      CALL obs_mpp_sum_integer( inlav1obs, inlav1obsmpp ) 
     427      CALL obs_mpp_sum_integer( ibdyv1obs, ibdyv1obsmpp ) 
    389428 
    390429      ! Variable 2 
     
    400439         &                 gdept_1d,              zmask2,               & 
    401440         &                 profdata%nqc,          profdata%var(2)%nvqc, & 
    402          &                 iosdv2obs,              ilanv2obs,           & 
    403          &                 inlav2obs,              ld_nea                ) 
     441         &                 iosdv2obs,             ilanv2obs,            & 
     442         &                 inlav2obs,             ld_nea,               & 
     443         &                 ibdyv2obs,             ld_bound_reject,      & 
     444         &                 iqc_cutoff       ) 
    404445 
    405446      CALL obs_mpp_sum_integer( iosdv2obs, iosdv2obsmpp ) 
    406447      CALL obs_mpp_sum_integer( ilanv2obs, ilanv2obsmpp ) 
    407448      CALL obs_mpp_sum_integer( inlav2obs, inlav2obsmpp ) 
     449      CALL obs_mpp_sum_integer( ibdyv2obs, ibdyv2obsmpp ) 
    408450 
    409451      ! ----------------------------------------------------------------------- 
     
    412454 
    413455      IF ( TRIM(profdata%cvars(1)) == 'UVEL' ) THEN 
    414          CALL obs_uv_rej( profdata, iuvchku, iuvchkv ) 
     456         CALL obs_uv_rej( profdata, iuvchku, iuvchkv, iqc_cutoff ) 
    415457         CALL obs_mpp_sum_integer( iuvchku, iuvchkumpp ) 
    416458         CALL obs_mpp_sum_integer( iuvchkv, iuvchkvmpp ) 
     
    429471      END DO 
    430472 
    431       ! We want all data which has qc flags = 0 
    432  
    433       llvalid%luse(:) = ( profdata%nqc(:)  <= 10 ) 
     473      ! We want all data which has qc flags <= iqc_cutoff 
     474 
     475      llvalid%luse(:) = ( profdata%nqc(:)  <= iqc_cutoff ) 
    434476      DO jvar = 1,profdata%nvar 
    435          llvvalid(jvar)%luse(:) = ( profdata%var(jvar)%nvqc(:) <= 10 ) 
     477         llvvalid(jvar)%luse(:) = ( profdata%var(jvar)%nvqc(:) <= iqc_cutoff ) 
    436478      END DO 
    437479 
     
    475517               &            iuvchku 
    476518         ENDIF 
     519         WRITE(numout,*) ' Remaining '//prodatqc%cvars(1)//' data near open boundary (removed) = ',& 
     520               &            ibdyv1obsmpp 
    477521         WRITE(numout,*) ' '//prodatqc%cvars(1)//' data accepted                             = ', & 
    478522            &            prodatqc%nvprotmpp(1) 
     
    492536               &            iuvchkv 
    493537         ENDIF 
     538         WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data near open boundary (removed) = ',& 
     539               &            ibdyv2obsmpp 
    494540         WRITE(numout,*) ' '//prodatqc%cvars(2)//' data accepted                             = ', & 
    495541            &            prodatqc%nvprotmpp(2) 
     
    644690            &        .AND. ( kobsmin(jobs) <= kmin0 ) ) ) THEN 
    645691            kobsstp(jobs) = -1 
    646             kobsqc(jobs)  = kobsqc(jobs) + 11 
     692            kobsqc(jobs)  = IBSET(kobsqc(jobs),13) 
    647693            kotdobs       = kotdobs + 1 
    648694            CYCLE 
     
    695741         IF ( ( kobsstp(jobs) < ( nit000 - 1 ) ) & 
    696742            & .OR.( kobsstp(jobs) > nitend ) ) THEN 
    697             kobsqc(jobs) = kobsqc(jobs) + 12 
     743            kobsqc(jobs) = IBSET(kobsqc(jobs),13) 
    698744            kotdobs = kotdobs + 1 
    699745            CYCLE 
     
    739785      &                    kobsno,                                        & 
    740786      &                    kobsyea, kobsmon, kobsday, kobshou, kobsmin,   & 
    741       &                    ktyp,    kobsqc,  kobsstp, kotdobs, kdailyavtypes ) 
     787      &                    ktyp,    kobsqc,  kobsstp, kotdobs, kdailyavtypes, & 
     788      &                    kqc_cutoff ) 
    742789      !!---------------------------------------------------------------------- 
    743790      !!                    ***  ROUTINE obs_coo_tim *** 
     
    783830      INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 
    784831         & kdailyavtypes    ! Types for daily averages 
     832      INTEGER, OPTIONAL, INTENT(IN) :: kqc_cutoff           ! QC cutoff value 
     833 
    785834      !! * Local declarations 
    786835      INTEGER :: jobs 
     836      INTEGER :: iqc_cutoff=255 
    787837 
    788838      !----------------------------------------------------------------------- 
     
    803853         DO jobs = 1, kobsno 
    804854             
    805             IF ( kobsqc(jobs) <= 10 ) THEN 
     855            IF ( kobsqc(jobs) <= iqc_cutoff ) THEN 
    806856                
    807857               IF ( ( kobsstp(jobs) == (nit000 - 1) ).AND.& 
    808858                  & ( ANY (kdailyavtypes(:) == ktyp(jobs)) ) ) THEN 
    809                   kobsqc(jobs) = kobsqc(jobs) + 14 
     859                  kobsqc(jobs) = IBSET(kobsqc(jobs),13) 
    810860                  kotdobs      = kotdobs + 1 
    811861                  CYCLE 
     
    850900      DO jobs = 1, kobsno 
    851901         IF ( ( kobsi(jobs) <= 0 ) .AND. ( kobsj(jobs) <= 0 ) ) THEN 
    852             kobsqc(jobs) = kobsqc(jobs) + 18 
     902            kobsqc(jobs) = IBSET(kobsqc(jobs),12) 
    853903            kgrdobs = kgrdobs + 1 
    854904         ENDIF 
     
    861911      &                       plam,   pphi,    pmask,            & 
    862912      &                       kobsqc, kosdobs, klanobs,          & 
    863       &                       knlaobs,ld_nea                     ) 
     913      &                       knlaobs,ld_nea,                    & 
     914      &                       kbdyobs,ld_bound_reject,           & 
     915      &                       kqc_cutoff                         ) 
    864916      !!---------------------------------------------------------------------- 
    865917      !!                    ***  ROUTINE obs_coo_spc_2d  *** 
     
    894946      INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: & 
    895947         & kobsqc             ! Observation quality control 
    896       INTEGER, INTENT(INOUT) :: kosdobs   ! Observations outside space domain 
    897       INTEGER, INTENT(INOUT) :: klanobs   ! Observations within a model land cell 
    898       INTEGER, INTENT(INOUT) :: knlaobs   ! Observations near land 
    899       LOGICAL, INTENT(IN) :: ld_nea       ! Flag observations near land 
     948      INTEGER, INTENT(INOUT) :: kosdobs          ! Observations outside space domain 
     949      INTEGER, INTENT(INOUT) :: klanobs          ! Observations within a model land cell 
     950      INTEGER, INTENT(INOUT) :: knlaobs          ! Observations near land 
     951      INTEGER, INTENT(INOUT) :: kbdyobs          ! Observations near boundary 
     952      LOGICAL, INTENT(IN)    :: ld_nea           ! Flag observations near land 
     953      LOGICAL, INTENT(IN)    :: ld_bound_reject  ! Flag observations near open boundary  
     954      INTEGER, INTENT(IN)    :: kqc_cutoff       ! Cutoff QC value 
     955 
    900956      !! * Local declarations 
    901957      REAL(KIND=wp), DIMENSION(2,2,kobsno) :: & 
    902958         & zgmsk              ! Grid mask 
     959 
     960      REAL(KIND=wp), DIMENSION(2,2,kobsno) :: & 
     961         & zbmsk              ! Boundary mask 
     962      REAL(KIND=wp), DIMENSION(jpi,jpj) :: zbdymask 
    903963      REAL(KIND=wp), DIMENSION(2,2,kobsno) :: & 
    904964         & zglam, &           ! Model longitude at grid points 
     
    917977         ! For invalid points use 2,2 
    918978 
    919          IF ( kobsqc(jobs) >= 10 ) THEN 
     979         IF ( kobsqc(jobs) >= kqc_cutoff ) THEN 
    920980 
    921981            igrdi(1,1,jobs) = 1 
     
    9421002 
    9431003      END DO 
     1004 
     1005      IF (ln_bdy) THEN 
     1006        ! Create a mask grid points in boundary rim 
     1007        IF (ld_bound_reject) THEN 
     1008           zbdymask(:,:) = 1.0_wp 
     1009           DO ji = 1, nb_bdy 
     1010              DO jj = 1, idx_bdy(ji)%nblen(1) 
     1011                 zbdymask(idx_bdy(ji)%nbi(jj,1),idx_bdy(ji)%nbj(jj,1)) = 0.0_wp 
     1012              ENDDO 
     1013           ENDDO 
     1014 
     1015           CALL obs_int_comm_2d( 2, 2, kobsno, kpi, kpj, igrdi, igrdj, zbdymask, zbmsk ) 
     1016        ENDIF 
     1017      ENDIF 
     1018 
    9441019       
    9451020      CALL obs_int_comm_2d( 2, 2, kobsno, kpi, kpj, igrdi, igrdj, pmask, zgmsk ) 
     
    9501025 
    9511026         ! Skip bad observations 
    952          IF ( kobsqc(jobs) >= 10 ) CYCLE 
     1027         IF ( kobsqc(jobs) >= kqc_cutoff ) CYCLE 
    9531028 
    9541029         ! Flag if the observation falls outside the model spatial domain 
     
    9571032            &  .OR. ( pobsphi(jobs) <  -90. ) & 
    9581033            &  .OR. ( pobsphi(jobs) >   90. ) ) THEN 
    959             kobsqc(jobs) = kobsqc(jobs) + 11 
     1034            kobsqc(jobs) = IBSET(kobsqc(jobs),11) 
    9601035            kosdobs = kosdobs + 1 
    9611036            CYCLE 
     
    9641039         ! Flag if the observation falls with a model land cell 
    9651040         IF ( SUM( zgmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN 
    966             kobsqc(jobs) = kobsqc(jobs)  + 12 
     1041            kobsqc(jobs) = IBSET(kobsqc(jobs),10) 
    9671042            klanobs = klanobs + 1 
    9681043            CYCLE 
     
    9781053               IF ( ( ABS( zgphi(ji,jj,jobs) - pobsphi(jobs) ) < 1.0e-6_wp ) & 
    9791054                  & .AND. & 
    980                   & ( ABS( zglam(ji,jj,jobs) - pobslam(jobs) ) < 1.0e-6_wp ) & 
    981                   & ) THEN 
     1055                  & ( ABS( MOD( zglam(ji,jj,jobs) - pobslam(jobs),360.0) ) & 
     1056                  & < 1.0e-6_wp ) ) THEN 
    9821057                  lgridobs = .TRUE. 
    9831058                  iig = ji 
     
    9921067         IF (lgridobs) THEN 
    9931068            IF (zgmsk(iig,ijg,jobs) == 0.0_wp ) THEN 
    994                kobsqc(jobs) = kobsqc(jobs) + 12 
     1069               kobsqc(jobs) = IBSET(kobsqc(jobs),10) 
    9951070               klanobs = klanobs + 1 
    9961071               CYCLE 
     
    10001075         ! Flag if the observation falls is close to land 
    10011076         IF ( MINVAL( zgmsk(1:2,1:2,jobs) ) == 0.0_wp) THEN 
    1002             IF (ld_nea) kobsqc(jobs) = kobsqc(jobs) + 14 
    10031077            knlaobs = knlaobs + 1 
    1004             CYCLE 
     1078            IF (ld_nea) THEN 
     1079               kobsqc(jobs) = IBSET(kobsqc(jobs),9) 
     1080               CYCLE 
     1081            ENDIF 
     1082         ENDIF 
     1083 
     1084         IF (ln_bdy) THEN 
     1085         ! Flag if the observation falls close to the boundary rim 
     1086           IF (ld_bound_reject) THEN 
     1087              IF ( MINVAL( zbmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN 
     1088                 kobsqc(jobs) = IBSET(kobsqc(jobs),8) 
     1089                 kbdyobs = kbdyobs + 1 
     1090                 CYCLE 
     1091              ENDIF 
     1092              ! for observations on the grid... 
     1093              IF (lgridobs) THEN 
     1094                 IF (zbmsk(iig,ijg,jobs) == 0.0_wp ) THEN 
     1095                    kobsqc(jobs) = IBSET(kobsqc(jobs),8) 
     1096                    kbdyobs = kbdyobs + 1 
     1097                    CYCLE 
     1098                 ENDIF 
     1099              ENDIF 
     1100            ENDIF 
    10051101         ENDIF 
    10061102             
     
    10151111      &                       plam,    pphi,    pdep,    pmask, & 
    10161112      &                       kpobsqc, kobsqc,  kosdobs,        & 
    1017       &                       klanobs, knlaobs, ld_nea          ) 
     1113      &                       klanobs, knlaobs, ld_nea,         & 
     1114      &                       kbdyobs, ld_bound_reject,         & 
     1115      &                       kqc_cutoff                        ) 
    10181116      !!---------------------------------------------------------------------- 
    10191117      !!                    ***  ROUTINE obs_coo_spc_3d  *** 
     
    10771175      INTEGER, INTENT(INOUT) :: klanobs     ! Observations within a model land cell 
    10781176      INTEGER, INTENT(INOUT) :: knlaobs     ! Observations near land 
     1177      INTEGER, INTENT(INOUT) :: kbdyobs     ! Observations near boundary 
    10791178      LOGICAL, INTENT(IN) :: ld_nea         ! Flag observations near land 
     1179      LOGICAL, INTENT(IN) :: ld_bound_reject  ! Flag observations near open boundary 
     1180      INTEGER, INTENT(IN) :: kqc_cutoff     ! Cutoff QC value 
     1181 
    10801182      !! * Local declarations 
    10811183      REAL(KIND=wp), DIMENSION(2,2,kpk,kprofno) :: & 
    10821184         & zgmsk              ! Grid mask 
     1185      REAL(KIND=wp), DIMENSION(2,2,kprofno) :: & 
     1186         & zbmsk              ! Boundary mask 
     1187      REAL(KIND=wp), DIMENSION(jpi,jpj) :: zbdymask 
    10831188      REAL(KIND=wp), DIMENSION(2,2,kpk,kprofno) :: & 
    10841189         & zgdepw 
     
    11001205         ! For invalid points use 2,2 
    11011206 
    1102          IF ( kpobsqc(jobs) >= 10 ) THEN 
     1207         IF ( kpobsqc(jobs) >= kqc_cutoff ) THEN 
    11031208             
    11041209            igrdi(1,1,jobs) = 1 
     
    11251230          
    11261231      END DO 
     1232 
     1233      IF (ln_bdy) THEN 
     1234        ! Create a mask grid points in boundary rim 
     1235        IF (ld_bound_reject) THEN            
     1236           zbdymask(:,:) = 1.0_wp 
     1237           DO ji = 1, nb_bdy 
     1238              DO jj = 1, idx_bdy(ji)%nblen(1) 
     1239                 zbdymask(idx_bdy(ji)%nbi(jj,1),idx_bdy(ji)%nbj(jj,1)) = 0.0_wp 
     1240              ENDDO 
     1241           ENDDO 
     1242        ENDIF 
     1243    
     1244        CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, zbdymask, zbmsk ) 
     1245      ENDIF 
    11271246       
    11281247      CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, pmask, zgmsk ) 
    11291248      CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, plam, zglam ) 
    11301249      CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, pphi, zgphi ) 
    1131       IF ( .NOT.( ln_zps .OR. ln_zco ) ) THEN 
    1132         ! Need to know the bathy depth for each observation for sco 
    1133         CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, gdepw_n(:,:,:), & 
     1250      CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, gdepw_n(:,:,:), & 
    11341251        &                     zgdepw ) 
    1135       ENDIF 
    11361252 
    11371253      DO jobs = 1, kprofno 
    11381254 
    11391255         ! Skip bad profiles 
    1140          IF ( kpobsqc(jobs) >= 10 ) CYCLE 
     1256         IF ( kpobsqc(jobs) >= kqc_cutoff ) CYCLE 
    11411257 
    11421258         ! Check if this observation is on a grid point 
     
    11491265               IF ( ( ABS( zgphi(ji,jj,jobs) - pobsphi(jobs) ) < 1.0e-6_wp ) & 
    11501266                  & .AND. & 
    1151                   & ( ABS( zglam(ji,jj,jobs) - pobslam(jobs) ) < 1.0e-6_wp ) & 
     1267                  & ( ABS( MOD( zglam(ji,jj,jobs) - pobslam(jobs),360.0) ) < 1.0e-6_wp ) & 
    11521268                  & ) THEN 
    11531269                  lgridobs = .TRUE. 
     
    11761292               &  .OR. ( pobsdep(jobsp) < 0.0          )       & 
    11771293               &  .OR. ( pobsdep(jobsp) > gdepw_1d(kpk)) ) THEN 
    1178                kobsqc(jobsp) = kobsqc(jobsp) + 11 
     1294               kobsqc(jobsp) = IBSET(kobsqc(jobsp),11) 
    11791295               kosdobs = kosdobs + 1 
    11801296               CYCLE 
    11811297            ENDIF 
    11821298 
    1183             ! To check if an observations falls within land there are two cases: 
    1184             ! 1: z-coordibnates, where the check uses the mask 
    1185             ! 2: terrain following (eg s-coordinates),  
    1186             !    where we use the depth of the bottom cell to mask observations 
     1299            ! To check if an observations falls within land: 
    11871300              
    1188             IF( ln_zps .OR. ln_zco ) THEN !(CASE 1) 
     1301            ! Flag if the observation is deeper than the bathymetry 
     1302            ! Or if it is within the mask 
     1303            IF ( ALL( zgdepw(1:2,1:2,kpk,jobs) < pobsdep(jobsp) ) & 
     1304               &     .OR. & 
     1305               &  ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) & 
     1306               &  == 0.0_wp) ) THEN 
     1307               kobsqc(jobsp) = IBSET(kobsqc(jobsp),10) 
     1308               klanobs = klanobs + 1 
     1309               CYCLE 
     1310            ENDIF 
    11891311                
    1190                ! Flag if the observation falls with a model land cell 
    1191                IF ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) & 
    1192                   &  == 0.0_wp ) THEN 
    1193                   kobsqc(jobsp) = kobsqc(jobsp) + 12 
    1194                   klanobs = klanobs + 1 
    1195                   CYCLE 
    1196                ENDIF 
    1197               
    1198                ! Flag if the observation is close to land 
    1199               IF ( MINVAL( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) == & 
    1200                   &  0.0_wp) THEN 
    1201                   knlaobs = knlaobs + 1 
    1202                   IF (ld_nea) THEN    
    1203                      kobsqc(jobsp) = kobsqc(jobsp) + 14  
    1204                   ENDIF  
    1205                ENDIF 
    1206               
    1207             ELSE ! Case 2 
    1208   
    1209                ! Flag if the observation is deeper than the bathymetry 
    1210                ! Or if it is within the mask 
    1211                IF ( ALL( zgdepw(1:2,1:2,kpk,jobs) < pobsdep(jobsp) ) & 
    1212                   &     .OR. & 
    1213                   &  ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) & 
    1214                   &  == 0.0_wp) ) THEN 
    1215                   kobsqc(jobsp) = kobsqc(jobsp) + 12 
    1216                   klanobs = klanobs + 1 
    1217                   CYCLE 
    1218                ENDIF 
    1219                 
    1220                ! Flag if the observation is close to land 
    1221                IF ( ll_next_to_land ) THEN 
    1222                   knlaobs = knlaobs + 1 
    1223                   IF (ld_nea) THEN    
    1224                      kobsqc(jobsp) = kobsqc(jobsp) + 14  
    1225                   ENDIF  
    1226                ENDIF 
     1312            ! Flag if the observation is close to land 
     1313            IF ( ll_next_to_land ) THEN 
     1314               knlaobs = knlaobs + 1 
     1315               IF (ld_nea) THEN    
     1316                  kobsqc(jobsp) = IBSET(kobsqc(jobsp),10) 
     1317               ENDIF  
    12271318            ENDIF 
    12281319             
     
    12321323            IF (lgridobs) THEN 
    12331324               IF (zgmsk(iig,ijg,kobsk(jobsp)-1,jobs) == 0.0_wp ) THEN 
    1234                   kobsqc(jobsp) = kobsqc(jobsp) + 12 
     1325                  kobsqc(jobsp) = IBSET(kobsqc(jobs),10) 
    12351326                  klanobs = klanobs + 1 
    12361327                  CYCLE 
     
    12501341            ENDIF 
    12511342             
     1343            IF (ln_bdy) THEN 
     1344               ! Flag if the observation falls close to the boundary rim 
     1345               IF (ld_bound_reject) THEN 
     1346                  IF ( MINVAL( zbmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN 
     1347                     kobsqc(jobsp) = IBSET(kobsqc(jobs),8) 
     1348                     kbdyobs = kbdyobs + 1 
     1349                     CYCLE 
     1350                  ENDIF 
     1351                  ! for observations on the grid... 
     1352                  IF (lgridobs) THEN 
     1353                     IF (zbmsk(iig,ijg,jobs) == 0.0_wp ) THEN 
     1354                        kobsqc(jobsp) = IBSET(kobsqc(jobs),8) 
     1355                        kbdyobs = kbdyobs + 1 
     1356                        CYCLE 
     1357                     ENDIF 
     1358                  ENDIF 
     1359               ENDIF 
     1360            ENDIF 
     1361             
    12521362         END DO 
    12531363      END DO 
     
    12551365   END SUBROUTINE obs_coo_spc_3d 
    12561366 
    1257    SUBROUTINE obs_pro_rej( profdata ) 
     1367   SUBROUTINE obs_pro_rej( profdata, kqc_cutoff ) 
    12581368      !!---------------------------------------------------------------------- 
    12591369      !!                    ***  ROUTINE obs_pro_rej *** 
     
    12731383      !! * Arguments 
    12741384      TYPE(obs_prof), INTENT(INOUT) :: profdata     ! Profile data 
     1385      INTEGER, INTENT(IN) :: kqc_cutoff             ! QC cutoff value 
     1386 
    12751387      !! * Local declarations 
    12761388      INTEGER :: jprof 
     
    12821394      DO jprof = 1, profdata%nprof 
    12831395 
    1284          IF ( profdata%nqc(jprof) > 10 ) THEN 
     1396         IF ( profdata%nqc(jprof) > kqc_cutoff ) THEN 
    12851397             
    12861398            DO jvar = 1, profdata%nvar 
     
    12901402                   
    12911403                  profdata%var(jvar)%nvqc(jobs) = & 
    1292                      & profdata%var(jvar)%nvqc(jobs) + 26 
     1404                     & IBSET(profdata%var(jvar)%nvqc(jobs),14) 
    12931405 
    12941406               END DO 
     
    13021414   END SUBROUTINE obs_pro_rej 
    13031415 
    1304    SUBROUTINE obs_uv_rej( profdata, knumu, knumv ) 
     1416   SUBROUTINE obs_uv_rej( profdata, knumu, knumv, kqc_cutoff ) 
    13051417      !!---------------------------------------------------------------------- 
    13061418      !!                    ***  ROUTINE obs_uv_rej *** 
     
    13221434      INTEGER, INTENT(INOUT) :: knumu             ! Number of u rejected 
    13231435      INTEGER, INTENT(INOUT) :: knumv             ! Number of v rejected 
     1436      INTEGER, INTENT(IN) :: kqc_cutoff           ! QC cutoff value 
     1437 
    13241438      !! * Local declarations 
    13251439      INTEGER :: jprof 
     
    13411455         DO jobs = profdata%npvsta(jprof,1), profdata%npvend(jprof,1) 
    13421456             
    1343             IF ( ( profdata%var(1)%nvqc(jobs) > 10 ) .AND. & 
    1344                & ( profdata%var(2)%nvqc(jobs) <= 10) ) THEN 
    1345                profdata%var(2)%nvqc(jobs) = profdata%var(2)%nvqc(jobs) + 42 
     1457            IF ( ( profdata%var(1)%nvqc(jobs) >  kqc_cutoff ) .AND. & 
     1458               & ( profdata%var(2)%nvqc(jobs) <=  kqc_cutoff) ) THEN 
     1459               profdata%var(2)%nvqc(jobs) = IBSET(profdata%var(1)%nvqc(jobs),15) 
    13461460               knumv = knumv + 1 
    13471461            ENDIF 
    1348             IF ( ( profdata%var(2)%nvqc(jobs) > 10 ) .AND. & 
    1349                & ( profdata%var(1)%nvqc(jobs) <= 10) ) THEN 
    1350                profdata%var(1)%nvqc(jobs) = profdata%var(1)%nvqc(jobs) + 42 
     1462            IF ( ( profdata%var(2)%nvqc(jobs) >  kqc_cutoff ) .AND. & 
     1463               & ( profdata%var(1)%nvqc(jobs) <=  kqc_cutoff) ) THEN 
     1464               profdata%var(1)%nvqc(jobs) = IBSET(profdata%var(1)%nvqc(jobs),15) 
    13511465               knumu = knumu + 1 
    13521466            ENDIF 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/OBS/obs_profiles_def.F90

    r6140 r9023  
    104104      ! Bookkeeping arrays with sizes equal to number of variables 
    105105 
    106       CHARACTER(len=6), POINTER, DIMENSION(:) :: & 
     106      CHARACTER(len=8), POINTER, DIMENSION(:) :: & 
    107107         & cvars          !: Variable names 
    108108 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_prof.F90

    r6140 r9023  
    8787      CHARACTER(LEN=15), PARAMETER :: cpname='obs_rea_prof' 
    8888      CHARACTER(len=8) :: clrefdate 
    89       CHARACTER(len=6), DIMENSION(:), ALLOCATABLE :: clvars 
     89      CHARACTER(len=8), DIMENSION(:), ALLOCATABLE :: clvars 
    9090      INTEGER :: jvar 
    9191      INTEGER :: ji 
     
    307307            inowin = 0 
    308308            DO ji = 1, inpfiles(jj)%nobs 
    309                IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE 
    310                IF ( ( inpfiles(jj)%ivqc(ji,1) > 2 ) .AND. & 
    311                   & ( inpfiles(jj)%ivqc(ji,2) > 2 )) CYCLE 
     309               IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE 
     310               IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. & 
     311                  & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE 
    312312               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. & 
    313313                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN 
     
    325325            inowin = 0 
    326326            DO ji = 1, inpfiles(jj)%nobs 
    327                IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE 
    328                IF ( ( inpfiles(jj)%ivqc(ji,1) > 2 ) .AND. & 
    329                   & ( inpfiles(jj)%ivqc(ji,2) > 2 )) CYCLE 
     327               IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE 
     328               IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. & 
     329                  & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE 
    330330               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. & 
    331331                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN 
     
    351351            inowin = 0 
    352352            DO ji = 1, inpfiles(jj)%nobs 
    353                IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE 
    354                IF ( ( inpfiles(jj)%ivqc(ji,1) > 2 ) .AND. & 
    355                   & ( inpfiles(jj)%ivqc(ji,2) > 2 )) CYCLE 
     353               IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE 
     354               IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. & 
     355                  & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE 
    356356               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. & 
    357357                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN 
     
    373373 
    374374            DO ji = 1, inpfiles(jj)%nobs 
    375                IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE 
    376                IF ( ( inpfiles(jj)%ivqc(ji,1) > 2 ) .AND. & 
    377                   & ( inpfiles(jj)%ivqc(ji,2) > 2 )) CYCLE 
     375               IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE 
     376               IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. & 
     377                  & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE 
    378378               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. & 
    379379                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN 
     
    388388                        IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) & 
    389389                           & CYCLE 
    390                         IF ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. & 
    391                            & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN 
     390                        IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,1),2) .AND. & 
     391                           & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN 
    392392                           ivar1t0 = ivar1t0 + 1 
    393393                        ENDIF 
     
    398398                        IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) & 
    399399                           & CYCLE 
    400                         IF ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. & 
    401                            & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN 
     400                        IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,2),2) .AND. & 
     401                           & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN 
    402402                           ivar2t0 = ivar2t0 + 1 
    403403                        ENDIF 
     
    407407                     IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) & 
    408408                        & CYCLE 
    409                      IF ( ( ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. & 
    410                         &     ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. & 
    411                         &     ldvar1 ) .OR. & 
    412                         & ( ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. & 
    413                         &     ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. & 
     409                     IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,1),2) .AND. & 
     410                        &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. & 
     411                        &    ldvar1 ) .OR. & 
     412                        & ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,2),2) .AND. & 
     413                        &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. & 
    414414                        &     ldvar2 ) ) THEN 
    415415                        ip3dt = ip3dt + 1 
     
    437437      DO jj = 1, inobf 
    438438         DO ji = 1, inpfiles(jj)%nobs 
    439             IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE 
    440             IF ( ( inpfiles(jj)%ivqc(ji,1) > 2 ) .AND. & 
    441                & ( inpfiles(jj)%ivqc(ji,2) > 2 )) CYCLE 
     439            IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE 
     440            IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. & 
     441               & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE 
    442442            IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. & 
    443443               & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN 
     
    452452      DO jj = 1, inobf 
    453453         DO ji = 1, inpfiles(jj)%nobs 
    454             IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE 
    455             IF ( ( inpfiles(jj)%ivqc(ji,1) > 2 ) .AND. & 
    456                & ( inpfiles(jj)%ivqc(ji,2) > 2 )) CYCLE 
     454            IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE 
     455            IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. & 
     456               & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE 
    457457            IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. & 
    458458               & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN 
     
    501501         ji = iprofidx(iindx(jk)) 
    502502 
    503          IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE 
    504          IF ( ( inpfiles(jj)%ivqc(ji,1) > 2 ) .AND. & 
    505             & ( inpfiles(jj)%ivqc(ji,2) > 2 )) CYCLE 
     503            IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE 
     504            IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. & 
     505               & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE 
    506506 
    507507         IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND.  & 
     
    518518            IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE 
    519519 
    520             IF ( ( inpfiles(jj)%ivqc(ji,1) > 2 ) .AND. & 
    521                & ( inpfiles(jj)%ivqc(ji,2) > 2 ) ) CYCLE 
     520            IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE 
     521            IF ( BTEST(inpfiles(jj)%ivqc(ji,1),2) .AND. & 
     522               & BTEST(inpfiles(jj)%ivqc(ji,2),2) ) CYCLE 
    522523 
    523524            loop_prof : DO ij = 1, inpfiles(jj)%nlev 
     
    526527                  & CYCLE 
    527528 
    528                IF ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. & 
    529                   & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN 
     529               IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,1),2) .AND. & 
     530                  & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN 
    530531 
    531532                  llvalprof = .TRUE.  
     
    534535               ENDIF 
    535536 
    536                IF ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. & 
    537                   & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN 
     537               IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,2),2) .AND. & 
     538                  & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN 
    538539 
    539540                  llvalprof = .TRUE.  
     
    615616                  IF (ldsatt) THEN 
    616617 
    617                      IF ( ( ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. & 
    618                         &     ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. & 
    619                         &     ldvar1 ) .OR. & 
    620                         & ( ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. & 
    621                         &     ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. & 
    622                         &     ldvar2 ) ) THEN 
     618                     IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,1),2) .AND. & 
     619                        &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. & 
     620                        &    ldvar1 ) .OR. & 
     621                        & ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,2),2) .AND. & 
     622                        &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. & 
     623                        &   ldvar2 ) ) THEN 
    623624                        ip3dt = ip3dt + 1 
    624625                     ELSE 
     
    628629                  ENDIF 
    629630 
    630                   IF ( ( ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. & 
    631                      &     ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. & 
    632                      &       ldvar1 ) .OR. ldsatt ) THEN 
     631                  IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,1),2) .AND. & 
     632                    &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. & 
     633                    &    ldvar1 ) .OR. ldsatt ) THEN 
    633634 
    634635                     IF (ldsatt) THEN 
     
    661662 
    662663                     ! Profile var1 value 
    663                      IF ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. & 
    664                         & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN 
     664                     IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,1),2) .AND. & 
     665                        & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN 
    665666                        profdata%var(1)%vobs(ivar1t) = & 
    666667                           &                inpfiles(jj)%pob(ij,ji,1) 
     
    692693                  ENDIF 
    693694 
    694                   IF ( ( ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. & 
    695                      &   ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) .AND. & 
     695                  IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,2),2) .AND. & 
     696                     &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2)    .AND. & 
    696697                     &   ldvar2 ) .OR. ldsatt ) THEN 
    697698 
     
    725726 
    726727                     ! Profile var2 value 
    727                      IF ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. & 
    728                         & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN 
     728                  IF (  ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,2),2) ) .AND. & 
     729                    &   ( .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2)    ) ) THEN 
    729730                        profdata%var(2)%vobs(ivar2t) = & 
    730731                           &                inpfiles(jj)%pob(ij,ji,2) 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_surf.F90

    r6140 r9023  
    7777      CHARACTER(LEN=11), PARAMETER :: cpname='obs_rea_surf' 
    7878      CHARACTER(len=8) :: clrefdate 
    79       CHARACTER(len=6), DIMENSION(:), ALLOCATABLE :: clvars 
     79      CHARACTER(len=8), DIMENSION(:), ALLOCATABLE :: clvars 
    8080      INTEGER :: ji 
    8181      INTEGER :: jj 
     
    172172 
    173173            !------------------------------------------------------------------ 
    174             !  Read the profile file into inpfiles 
     174            !  Read the surface file into inpfiles 
    175175            !------------------------------------------------------------------ 
    176176            CALL init_obfbdata( inpfiles(jj) ) 
     
    296296                  ENDIF 
    297297                  llvalprof = .FALSE. 
    298                   IF ( ( inpfiles(jj)%ivlqc(1,ji,1) == 1 ) .OR. & 
    299                      & ( inpfiles(jj)%ivlqc(1,ji,1) == 2 ) ) THEN 
     298                  IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(1,ji,1),2) ) THEN 
    300299                     iobs = iobs + 1 
    301300                  ENDIF 
     
    370369            ! Set observation information 
    371370 
    372             IF ( ( inpfiles(jj)%ivlqc(1,ji,1) == 1 ) .OR. & 
    373                & ( inpfiles(jj)%ivlqc(1,ji,1) == 2 ) ) THEN 
     371            IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(1,ji,1),2) ) THEN 
    374372 
    375373               iobs = iobs + 1 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/OBS/obs_readmdt.F90

    r6140 r9023  
    154154 
    155155! mark any masked data with a QC flag 
    156          IF( zobsmask(1) == 0 )   sladata%nqc(jobs) = 11 
     156         IF( zobsmask(1) == 0 )   sladata%nqc(jobs) = IBSET(sladata%nqc(jobs),15) 
    157157 
    158158         END DO 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/OBS/obs_sstbias.F90

    r6140 r9023  
    11MODULE obs_sstbias 
    22   !!====================================================================== 
    3    !!                       ***  MODULE obs_readsstbias  *** 
    4    !! Observation diagnostics: Read the bias for SLA data 
     3   !!                       ***  MODULE obs_sstbias  *** 
     4   !! Observation diagnostics: Read the bias for SST data 
    55   !!====================================================================== 
    66   !!---------------------------------------------------------------------- 
    7    !!   obs_rea_sstbias : Driver for reading altimeter bias 
     7   !!   obs_app_sstbias : Driver for reading and applying the SST bias 
    88   !!---------------------------------------------------------------------- 
    99   !! * Modules used    
     
    139139               cl_bias_files(jtype) ) 
    140140            ! Get the SST bias data 
    141             CALL iom_get( numsstbias, jpdom_data, 'sst', z_sstbias_2d(:,:), 1 ) 
     141            CALL iom_get( numsstbias, jpdom_data, 'tn', z_sstbias_2d(:,:), 1 ) 
    142142            z_sstbias(:,:,jtype) = z_sstbias_2d(:,:)        
    143143            ! Close the file 
     
    190190               igrdj_tmp(:,:,jt) = igrdj(:,:,jobs) 
    191191               zglam_tmp(:,:,jt) = zglam(:,:,jobs) 
    192                zgphi_tmp(:,:,jt) = zgphi(:,:,jobs) 
    193192               zgphi_tmp(:,:,jt) = zgphi(:,:,jobs) 
    194193               zmask_tmp(:,:,jt) = zmask(:,:,jobs) 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/OBS/obs_surf_def.F90

    r6140 r9023  
    5050      INTEGER :: npj 
    5151      INTEGER :: nsurfup    !: Observation counter used in obs_oper 
     52      INTEGER :: nrec       !: Number of surface observation records in window 
    5253 
    5354      ! Arrays with size equal to the number of surface observations 
     
    5657         & mi,   &        !: i-th grid coord. for interpolating to surface observation 
    5758         & mj,   &        !: j-th grid coord. for interpolating to surface observation 
     59         & mt,   &        !: time record number for gridded data 
    5860         & nsidx,&        !: Surface observation number 
    5961         & nsfil,&        !: Surface observation number in file 
     
    6769         & ntyp           !: Type of surface observation product 
    6870 
    69       CHARACTER(len=6), POINTER, DIMENSION(:) :: & 
     71      CHARACTER(len=8), POINTER, DIMENSION(:) :: & 
    7072         & cvars          !: Variable names 
    7173 
     
    9395         & nsstpmpp       !: Global number of surface observations per time step 
    9496 
     97      ! Arrays with size equal to the number of observation records in the window 
     98      INTEGER, POINTER, DIMENSION(:) :: & 
     99         & mrecstp   ! Time step of the records 
     100 
    95101      ! Arrays used to store source indices when  
    96102      ! compressing obs_surf derived types 
     
    100106      INTEGER, POINTER, DIMENSION(:) :: & 
    101107         & nsind          !: Source indices of surface data in compressed data 
     108 
     109      ! Is this a gridded product? 
     110      
     111      LOGICAL :: lgrid 
    102112 
    103113   END TYPE obs_surf 
     
    160170         & surf%mi(ksurf),      & 
    161171         & surf%mj(ksurf),      & 
     172         & surf%mt(ksurf),      & 
    162173         & surf%nsidx(ksurf),   & 
    163174         & surf%nsfil(ksurf),   & 
     
    176187         & ) 
    177188 
     189      surf%mt(:) = -1 
     190 
    178191 
    179192      ! Allocate arrays of number of surface data size * number of variables 
     
    190203         & ) 
    191204 
     205      surf%rext(:,:) = 0.0_wp  
     206 
    192207      ! Allocate arrays of number of time step size 
    193208 
     
    217232 
    218233      surf%nsurfup     = 0 
     234       
     235      ! Not gridded by default 
     236           
     237      surf%lgrid       = .FALSE. 
    219238               
    220239   END SUBROUTINE obs_surf_alloc 
     
    242261         & surf%mi,      & 
    243262         & surf%mj,      & 
     263         & surf%mt,      & 
    244264         & surf%nsidx,   & 
    245265         & surf%nsfil,   & 
     
    370390            newsurf%mi(insurf)    = surf%mi(ji) 
    371391            newsurf%mj(insurf)    = surf%mj(ji) 
     392            newsurf%mt(insurf)    = surf%mt(ji) 
    372393            newsurf%nsidx(insurf) = surf%nsidx(ji) 
    373394            newsurf%nsfil(insurf) = surf%nsfil(ji) 
     
    414435      newsurf%nstp     = surf%nstp 
    415436      newsurf%cvars(:) = surf%cvars(:) 
     437       
     438      ! Set gridded stuff 
     439       
     440      newsurf%mt(insurf)    = surf%mt(ji) 
    416441  
    417442      ! Deallocate temporary data 
     
    454479         oldsurf%mi(jj)    = surf%mi(ji) 
    455480         oldsurf%mj(jj)    = surf%mj(ji) 
     481         oldsurf%mt(jj)    = surf%mt(ji) 
    456482         oldsurf%nsidx(jj) = surf%nsidx(ji) 
    457483         oldsurf%nsfil(jj) = surf%nsfil(ji) 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/OBS/obs_write.F90

    r6140 r9023  
    88   !!   obs_wri_prof   : Write profile observations in feedback format 
    99   !!   obs_wri_surf   : Write surface observations in feedback format 
    10    !!   obs_wri_stats : Print basic statistics on the data being written out 
     10   !!   obs_wri_stats  : Print basic statistics on the data being written out 
    1111   !!---------------------------------------------------------------------- 
    1212 
     
    8383      TYPE(obfbdata) :: fbdata 
    8484      CHARACTER(LEN=40) :: clfname 
    85       CHARACTER(LEN=6) :: clfiletype 
     85      CHARACTER(LEN=10) :: clfiletype 
    8686      INTEGER :: ilevel 
    8787      INTEGER :: jvar 
     
    196196         fbdata%ivqc(jo,:)    = profdata%ivqc(jo,:) 
    197197         fbdata%ivqcf(:,jo,:) = profdata%ivqcf(:,jo,:) 
    198          IF ( profdata%nqc(jo) > 10 ) THEN 
    199             fbdata%ioqc(jo)    = 4 
     198         IF ( profdata%nqc(jo) > 255 ) THEN 
     199            fbdata%ioqc(jo)    = IBSET(profdata%nqc(jo),2) 
    200200            fbdata%ioqcf(1,jo) = profdata%nqcf(1,jo) 
    201             fbdata%ioqcf(2,jo) = profdata%nqc(jo) - 10 
     201            fbdata%ioqcf(2,jo) = profdata%nqc(jo) 
    202202         ELSE 
    203203            fbdata%ioqc(jo)    = profdata%nqc(jo) 
     
    236236               fbdata%idqc(ik,jo)        = profdata%var(jvar)%idqc(jk) 
    237237               fbdata%idqcf(:,ik,jo)     = profdata%var(jvar)%idqcf(:,jk) 
    238                IF ( profdata%var(jvar)%nvqc(jk) > 10 ) THEN 
    239                   fbdata%ivlqc(ik,jo,jvar) = 4 
     238               IF ( profdata%var(jvar)%nvqc(jk) > 255 ) THEN 
     239                  fbdata%ivlqc(ik,jo,jvar) = IBSET(profdata%var(jvar)%nvqc(jk),2) 
    240240                  fbdata%ivlqcf(1,ik,jo,jvar) = profdata%var(jvar)%nvqcf(1,jk) 
    241                   fbdata%ivlqcf(2,ik,jo,jvar) = profdata%var(jvar)%nvqc(jk) - 10 
     241!$AGRIF_DO_NOT_TREAT 
     242                  fbdata%ivlqcf(2,ik,jo,jvar) = IAND(profdata%var(jvar)%nvqc(jk),b'0000000011111111') 
     243!$AGRIF_END_DO_NOT_TREAT 
    242244               ELSE 
    243245                  fbdata%ivlqc(ik,jo,jvar) = profdata%var(jvar)%nvqc(jk) 
     
    320322      TYPE(obfbdata) :: fbdata 
    321323      CHARACTER(LEN=40) :: clfname         ! netCDF filename 
    322       CHARACTER(LEN=6) :: clfiletype 
     324      CHARACTER(LEN=10) :: clfiletype 
    323325      CHARACTER(LEN=12), PARAMETER :: cpname = 'obs_wri_surf' 
    324326      INTEGER :: jo 
     
    395397         END DO 
    396398 
    397       CASE('ICECON') 
     399      CASE('ICECONC') 
    398400 
    399401         CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, & 
     
    418420         END DO 
    419421 
     422      CASE('SSS') 
     423 
     424         CALL alloc_obfbdata( fbdata, 1, surfdata%nsurf, 1, & 
     425            &                 1 + iadd, iext, .TRUE. ) 
     426 
     427         clfiletype = 'sssfb' 
     428         fbdata%cname(1)      = surfdata%cvars(1) 
     429         fbdata%coblong(1)    = 'Sea surface salinity' 
     430         fbdata%cobunit(1)    = 'psu' 
     431         DO je = 1, iext 
     432            fbdata%cextname(je) = pext%cdname(je) 
     433            fbdata%cextlong(je) = pext%cdlong(je,1) 
     434            fbdata%cextunit(je) = pext%cdunit(je,1) 
     435         END DO 
     436         fbdata%caddlong(1,1) = 'Model interpolated SSS' 
     437         fbdata%caddunit(1,1) = 'psu' 
     438         fbdata%cgrid(1)      = 'T' 
     439         DO ja = 1, iadd 
     440            fbdata%caddname(1+ja) = padd%cdname(ja) 
     441            fbdata%caddlong(1+ja,1) = padd%cdlong(ja,1) 
     442            fbdata%caddunit(1+ja,1) = padd%cdunit(ja,1) 
     443         END DO 
     444 
     445      CASE DEFAULT 
     446 
     447         CALL ctl_stop( 'Unknown observation type '//TRIM(surfdata%cvars(1))//' in obs_wri_surf' ) 
     448 
    420449      END SELECT 
    421450 
     
    439468         fbdata%ivqc(jo,:)    = 0 
    440469         fbdata%ivqcf(:,jo,:) = 0 
    441          IF ( surfdata%nqc(jo) > 10 ) THEN 
     470         IF ( surfdata%nqc(jo) > 255 ) THEN 
    442471            fbdata%ioqc(jo)    = 4 
    443472            fbdata%ioqcf(1,jo) = 0 
    444             fbdata%ioqcf(2,jo) = surfdata%nqc(jo) - 10 
     473!$AGRIF_DO_NOT_TREAT 
     474            fbdata%ioqcf(2,jo) = IAND(surfdata%nqc(jo),b'0000000011111111') 
     475!$AGRIF_END_DO_NOT_TREAT 
    445476         ELSE 
    446477            fbdata%ioqc(jo)    = surfdata%nqc(jo) 
     
    474505         fbdata%idqc(1,jo)     = 0 
    475506         fbdata%idqcf(:,1,jo)  = 0 
    476          IF ( surfdata%nqc(jo) > 10 ) THEN 
     507         IF ( surfdata%nqc(jo) > 255 ) THEN 
    477508            fbdata%ivqc(jo,1)       = 4 
    478509            fbdata%ivlqc(1,jo,1)    = 4 
    479510            fbdata%ivlqcf(1,1,jo,1) = 0 
    480             fbdata%ivlqcf(2,1,jo,1) = surfdata%nqc(jo) - 10 
     511!$AGRIF_DO_NOT_TREAT 
     512            fbdata%ivlqcf(2,1,jo,1) = IAND(surfdata%nqc(jo),b'0000000011111111') 
     513!$AGRIF_END_DO_NOT_TREAT 
    481514         ELSE 
    482515            fbdata%ivqc(jo,1)       = surfdata%nqc(jo) 
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/OBS/obsinter_h2d.h90

    r2474 r9023  
    12401240         & zdum,  & 
    12411241         & zaamax 
    1242         
     1242       
     1243      imax = -1  
    12431244      ! Main computation 
    12441245      pflt = 1.0_wp 
Note: See TracChangeset for help on using the changeset viewer.