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 11361 – NEMO

Changeset 11361


Ignore:
Timestamp:
2019-07-29T11:26:23+02:00 (5 years ago)
Author:
jcastill
Message:

First version of the new observation branch - it compiles, but has not been tested

Location:
NEMO/branches/UKMO/r8395_obs_oper_update/NEMOGCM/NEMO/OPA_SRC
Files:
2 added
13 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/UKMO/r8395_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90

    r11350 r11361  
    7171   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   u_bkginc, v_bkginc   !: Increment to the u- & v-components  
    7272   REAL(wp), PUBLIC, DIMENSION(:)    , ALLOCATABLE ::   wgtiau               !: IAU weights for each time step 
    73 #if defined key_asminc 
    74    REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   ssh_iau           !: IAU-weighted sea surface height increment 
    75 #endif 
     73   REAL(wp), PUBLIC, DIMENSION(:,:)  , ALLOCATABLE ::   ssh_iau              !: IAU-weighted sea surface height increment 
    7674   !                                !!! time steps relative to the cycle interval [0,nitend-nit000-1] 
    7775   INTEGER , PUBLIC ::   nitbkg      !: Time step of the background state used in the Jb term 

  • NEMO/branches/UKMO/r8395_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90

    r11350 r11361  
    4646 
    4747   !! * Module variables 
    48    LOGICAL, PUBLIC :: ln_diaobs   !: Logical switch for the obs operator 
    49    LOGICAL :: ln_sstnight         !: Logical switch for night mean SST obs 
    50     
    51    INTEGER :: nn_1dint       !: Vertical interpolation method 
    52    INTEGER :: nn_2dint       !: Horizontal interpolation method 
     48   LOGICAL, PUBLIC :: & 
     49      &       ln_diaobs            !: Logical switch for the obs operator 
     50   LOGICAL :: ln_sstnight          !: Logical switch for night mean SST obs 
     51   LOGICAL :: ln_default_fp_indegs !: T=> Default obs footprint size specified in degrees, F=> in metres 
     52   LOGICAL :: ln_sla_fp_indegs     !: T=>     SLA obs footprint size specified in degrees, F=> in metres 
     53   LOGICAL :: ln_sst_fp_indegs     !: T=>     SST obs footprint size specified in degrees, F=> in metres 
     54   LOGICAL :: ln_sss_fp_indegs     !: T=>     SSS obs footprint size specified in degrees, F=> in metres 
     55   LOGICAL :: ln_sic_fp_indegs     !: T=> sea-ice obs footprint size specified in degrees, F=> in metres 
     56 
     57   REAL(wp) :: rn_default_avglamscl !: Default E/W diameter of observation footprint 
     58   REAL(wp) :: rn_default_avgphiscl !: Default N/S diameter of observation footprint 
     59   REAL(wp) :: rn_sla_avglamscl     !: E/W diameter of SLA observation footprint 
     60   REAL(wp) :: rn_sla_avgphiscl     !: N/S diameter of SLA observation footprint 
     61   REAL(wp) :: rn_sst_avglamscl     !: E/W diameter of SST observation footprint 
     62   REAL(wp) :: rn_sst_avgphiscl     !: N/S diameter of SST observation footprint 
     63   REAL(wp) :: rn_sss_avglamscl     !: E/W diameter of SSS observation footprint 
     64   REAL(wp) :: rn_sss_avgphiscl     !: N/S diameter of SSS observation footprint 
     65   REAL(wp) :: rn_sic_avglamscl     !: E/W diameter of sea-ice observation footprint 
     66   REAL(wp) :: rn_sic_avgphiscl     !: N/S diameter of sea-ice observation footprint 
     67 
     68   INTEGER :: nn_1dint         !: Vertical interpolation method 
     69   INTEGER :: nn_2dint_default !: Default horizontal interpolation method 
     70   INTEGER :: nn_2dint_sla     !: SLA horizontal interpolation method (-1 = default) 
     71   INTEGER :: nn_2dint_sst     !: SST horizontal interpolation method (-1 = default) 
     72   INTEGER :: nn_2dint_sss     !: SSS horizontal interpolation method (-1 = default) 
     73   INTEGER :: nn_2dint_sic     !: Seaice horizontal interpolation method (-1 = default) 
     74  
    5375   INTEGER, DIMENSION(imaxavtypes) :: & 
    5476      & nn_profdavtypes      !: Profile data types representing a daily average 
    … …  
    6284      & nextrsurf            !: Number of surface extra variables 
    6385   INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:) :: sstbias_type !SST bias type     
     86   INTEGER, DIMENSION(:), ALLOCATABLE :: & 
     87      & n2dintsurf           !: Interpolation option for surface variables 
     88   REAL(wp), DIMENSION(:), ALLOCATABLE :: & 
     89      & ravglamscl, &        !: E/W diameter of averaging footprint for surface variables 
     90      & ravgphiscl           !: N/S diameter of averaging footprint for surface variables 
     91   LOGICAL, DIMENSION(:), ALLOCATABLE :: & 
     92      & lfpindegs, &         !: T=> surface obs footprint size specified in degrees, F=> in metres 
     93      & llnightav            !: Logical for calculating night-time averages 
    6494   TYPE(obs_surf), PUBLIC, POINTER, DIMENSION(:) :: & 
    6595      & surfdata, &          !: Initial surface data 
    … …  
    6999      & profdataqc           !: Profile data after quality control 
    70100 
    71    CHARACTER(len=6), PUBLIC, DIMENSION(:), ALLOCATABLE :: & 
     101   CHARACTER(len=8), PUBLIC, DIMENSION(:), ALLOCATABLE :: & 
    72102      & cobstypesprof, &     !: Profile obs types 
    73103      & cobstypessurf        !: Surface obs types 
    … …  
    96126      !!        !  06-10  (A. Weaver) Cleaning and add controls 
    97127      !!        !  07-03  (K. Mogensen) General handling of profiles 
    98       !!        !  14-08  (J.While) Incorporated SST bias correction   
     128      !!        !  14-08  (J.While) Incorporated SST bias correction 
    99129      !!        !  15-02  (M. Martin) Simplification of namelist and code 
    100130      !!---------------------------------------------------------------------- 
    … …  
    112142      INTEGER :: jvar            ! Counter for variables 
    113143      INTEGER :: jfile           ! Counter for files 
     144      INTEGER :: jnumsstbias     ! Number of SST bias files to read and apply 
     145      INTEGER :: n2dint_type     ! Local version of nn_2dint* 
    114146 
    115147      CHARACTER(len=128), DIMENSION(jpmaxnfiles) :: & 
    116          & cn_profbfiles, &      ! T/S profile input filenames 
    117          & cn_sstfbfiles, &      ! Sea surface temperature input filenames 
    118          & cn_slafbfiles, &      ! Sea level anomaly input filenames 
    119          & cn_sicfbfiles, &      ! Seaice concentration input filenames 
    120          & cn_velfbfiles, &      ! Velocity profile input filenames 
    121          & cn_sstbias_files      ! SST bias input filenames 
     148         & cn_profbfiles,      & ! T/S profile input filenames 
     149         & cn_sstfbfiles,      & ! Sea surface temperature input filenames 
     150         & cn_slafbfiles,      & ! Sea level anomaly input filenames 
     151         & cn_sicfbfiles,      & ! Seaice concentration input filenames 
     152         & cn_velfbfiles,      & ! Velocity profile input filenames 
     153         & cn_sssfbfiles,      & ! Sea surface salinity input filenames 
     154         & cn_slchltotfbfiles, & ! Surface total              log10(chlorophyll) input filenames 
     155         & cn_slchldiafbfiles, & ! Surface diatom             log10(chlorophyll) input filenames 
     156         & cn_slchlnonfbfiles, & ! Surface non-diatom         log10(chlorophyll) input filenames 
     157         & cn_slchldinfbfiles, & ! Surface dinoflagellate     log10(chlorophyll) input filenames 
     158         & cn_slchlmicfbfiles, & ! Surface microphytoplankton log10(chlorophyll) input filenames 
     159         & cn_slchlnanfbfiles, & ! Surface nanophytoplankton  log10(chlorophyll) input filenames 
     160         & cn_slchlpicfbfiles, & ! Surface picophytoplankton  log10(chlorophyll) input filenames 
     161         & cn_schltotfbfiles,  & ! Surface total              chlorophyll        input filenames 
     162         & cn_slphytotfbfiles, & ! Surface total      log10(phytoplankton carbon) input filenames 
     163         & cn_slphydiafbfiles, & ! Surface diatom     log10(phytoplankton carbon) input filenames 
     164         & cn_slphynonfbfiles, & ! Surface non-diatom log10(phytoplankton carbon) input filenames 
     165         & cn_sspmfbfiles,     & ! Surface suspended particulate matter input filenames 
     166         & cn_sfco2fbfiles,    & ! Surface fugacity         of carbon dioxide input filenames 
     167         & cn_spco2fbfiles,    & ! Surface partial pressure of carbon dioxide input filenames 
     168         & cn_plchltotfbfiles, & ! Profile total log10(chlorophyll) input filenames 
     169         & cn_pchltotfbfiles,  & ! Profile total chlorophyll input filenames 
     170         & cn_pno3fbfiles,     & ! Profile nitrate input filenames 
     171         & cn_psi4fbfiles,     & ! Profile silicate input filenames 
     172         & cn_ppo4fbfiles,     & ! Profile phosphate input filenames 
     173         & cn_pdicfbfiles,     & ! Profile dissolved inorganic carbon input filenames 
     174         & cn_palkfbfiles,     & ! Profile alkalinity input filenames 
     175         & cn_pphfbfiles,      & ! Profile pH input filenames 
     176         & cn_po2fbfiles,      & ! Profile dissolved oxygen input filenames 
     177         & cn_sstbiasfiles       ! SST bias input filenames 
     178 
    122179      CHARACTER(LEN=128) :: & 
    123180         & cn_altbiasfile        ! Altimeter bias input filename 
    124       CHARACTER(len=128), DIMENSION(:,:), ALLOCATABLE :: & 
    125          & clproffiles, &        ! Profile filenames 
    126          & clsurffiles           ! Surface filenames 
    127181 
    128182      LOGICAL :: ln_t3d          ! Logical switch for temperature profiles 
    … …  
    131185      LOGICAL :: ln_sst          ! Logical switch for sea surface temperature 
    132186      LOGICAL :: ln_sic          ! Logical switch for sea ice concentration 
     187      LOGICAL :: ln_sss          ! Logical switch for sea surface salinity obs 
    133188      LOGICAL :: ln_vel3d        ! Logical switch for velocity (u,v) obs 
     189      LOGICAL :: ln_slchltot     ! Logical switch for surface total              log10(chlorophyll) obs 
     190      LOGICAL :: ln_slchldia     ! Logical switch for surface diatom             log10(chlorophyll) obs 
     191      LOGICAL :: ln_slchlnon     ! Logical switch for surface non-diatom         log10(chlorophyll) obs 
     192      LOGICAL :: ln_slchldin     ! Logical switch for surface dinoflagellate     log10(chlorophyll) obs 
     193      LOGICAL :: ln_slchlmic     ! Logical switch for surface microphytoplankton log10(chlorophyll) obs 
     194      LOGICAL :: ln_slchlnan     ! Logical switch for surface nanophytoplankton  log10(chlorophyll) obs 
     195      LOGICAL :: ln_slchlpic     ! Logical switch for surface picophytoplankton  log10(chlorophyll) obs 
     196      LOGICAL :: ln_schltot      ! Logical switch for surface total              chlorophyll        obs 
     197      LOGICAL :: ln_slphytot     ! Logical switch for surface total      log10(phytoplankton carbon) obs 
     198      LOGICAL :: ln_slphydia     ! Logical switch for surface diatom     log10(phytoplankton carbon) obs 
     199      LOGICAL :: ln_slphynon     ! Logical switch for surface non-diatom log10(phytoplankton carbon) obs 
     200      LOGICAL :: ln_sspm         ! Logical switch for surface suspended particulate matter obs 
     201      LOGICAL :: ln_sfco2        ! Logical switch for surface fugacity         of carbon dioxide obs 
     202      LOGICAL :: ln_spco2        ! Logical switch for surface partial pressure of carbon dioxide obs 
     203      LOGICAL :: ln_plchltot     ! Logical switch for profile total log10(chlorophyll) obs 
     204      LOGICAL :: ln_pchltot      ! Logical switch for profile total chlorophyll obs 
     205      LOGICAL :: ln_pno3         ! Logical switch for profile nitrate obs 
     206      LOGICAL :: ln_psi4         ! Logical switch for profile silicate obs 
     207      LOGICAL :: ln_ppo4         ! Logical switch for profile phosphate obs 
     208      LOGICAL :: ln_pdic         ! Logical switch for profile dissolved inorganic carbon obs 
     209      LOGICAL :: ln_palk         ! Logical switch for profile alkalinity obs 
     210      LOGICAL :: ln_pph          ! Logical switch for profile pH obs 
     211      LOGICAL :: ln_po2          ! Logical switch for profile dissolved oxygen obs 
    134212      LOGICAL :: ln_nea          ! Logical switch to remove obs near land 
    135213      LOGICAL :: ln_altbias      ! Logical switch for altimeter bias 
    136       LOGICAL :: ln_sstbias     !: Logical switch for bias corection of SST  
     214      LOGICAL :: ln_sstbias      ! Logical switch for bias correction of SST 
    137215      LOGICAL :: ln_ignmis       ! Logical switch for ignoring missing files 
    138216      LOGICAL :: ln_s_at_t       ! Logical switch to compute model S at T obs 
    139       LOGICAL :: llvar1          ! Logical for profile variable 1 
    140       LOGICAL :: llvar2          ! Logical for profile variable 1 
    141       LOGICAL :: llnightav       ! Logical for calculating night-time averages 
     217      LOGICAL :: ln_bound_reject ! Logical switch for rejecting obs near the boundary 
    142218      LOGICAL, DIMENSION(jpmaxnfiles) :: lmask ! Used for finding number of sstbias files 
    143219 
    144220      REAL(dp) :: rn_dobsini     ! Obs window start date YYYYMMDD.HHMMSS 
    145221      REAL(dp) :: rn_dobsend     ! Obs window end date   YYYYMMDD.HHMMSS 
    146       REAL(wp), POINTER, DIMENSION(:,:) :: & 
    147          & zglam1, &             ! Model longitudes for profile variable 1 
    148          & zglam2                ! Model longitudes for profile variable 2 
    149       REAL(wp), POINTER, DIMENSION(:,:) :: & 
    150          & zgphi1, &             ! Model latitudes for profile variable 1 
    151          & zgphi2                ! Model latitudes for profile variable 2 
     222 
     223      REAL(wp) :: ztype_avglamscl ! Local version of rn_*_avglamscl 
     224      REAL(wp) :: ztype_avgphiscl ! Local version of rn_*_avgphiscl 
     225 
     226      CHARACTER(len=128), DIMENSION(:,:), ALLOCATABLE :: & 
     227         & clproffiles, &        ! Profile filenames 
     228         & clsurffiles           ! Surface filenames 
     229 
     230      LOGICAL, DIMENSION(:), ALLOCATABLE :: llvar   ! Logical for profile variable read 
     231      LOGICAL :: ltype_fp_indegs ! Local version of ln_*_fp_indegs 
     232      LOGICAL :: ltype_night     ! Local version of ln_sstnight (false for other variables) 
     233 
    152234      REAL(wp), POINTER, DIMENSION(:,:,:) :: & 
    153          & zmask1, &             ! Model land/sea mask associated with variable 1 
    154          & zmask2                ! Model land/sea mask associated with variable 2 
     235         & zglam                 ! Model longitudes for profile variables 
     236      REAL(wp), POINTER, DIMENSION(:,:,:) :: & 
     237         & zgphi                 ! Model latitudes for profile variables 
     238      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: & 
     239         & zmask                 ! Model land/sea mask associated with variables 
     240 
    155241 
    156242      NAMELIST/namobs/ln_diaobs, ln_t3d, ln_s3d, ln_sla,              & 
    157          &            ln_sst, ln_sic, ln_vel3d,                       & 
    158          &            ln_altbias, ln_nea, ln_grid_global,             & 
    159          &            ln_grid_search_lookup,                          & 
    160          &            ln_ignmis, ln_s_at_t, ln_sstnight,              & 
     243         &            ln_sst, ln_sic, ln_sss, ln_vel3d,               & 
     244         &            ln_slchltot, ln_slchldia, ln_slchlnon,          & 
     245         &            ln_slchldin, ln_slchlmic, ln_slchlnan,          & 
     246         &            ln_slchlpic, ln_schltot,                        & 
     247         &            ln_slphytot, ln_slphydia, ln_slphynon,          & 
     248         &            ln_sspm,     ln_sfco2,    ln_spco2,             & 
     249         &            ln_plchltot, ln_pchltot,  ln_pno3,              & 
     250         &            ln_psi4,     ln_ppo4,     ln_pdic,              & 
     251         &            ln_palk,     ln_pph,      ln_po2,               & 
     252         &            ln_altbias, ln_sstbias, ln_nea,                 & 
     253         &            ln_grid_global, ln_grid_search_lookup,          & 
     254         &            ln_ignmis, ln_s_at_t, ln_bound_reject,          & 
     255         &            ln_sstnight, ln_default_fp_indegs,              & 
     256         &            ln_sla_fp_indegs, ln_sst_fp_indegs,             & 
     257         &            ln_sss_fp_indegs, ln_sic_fp_indegs,             & 
    161258         &            cn_profbfiles, cn_slafbfiles,                   & 
    162259         &            cn_sstfbfiles, cn_sicfbfiles,                   & 
    163          &            cn_velfbfiles, cn_altbiasfile,                  & 
     260         &            cn_velfbfiles, cn_sssfbfiles,                   & 
     261         &            cn_slchltotfbfiles, cn_slchldiafbfiles,         & 
     262         &            cn_slchlnonfbfiles, cn_slchldinfbfiles,         & 
     263         &            cn_slchlmicfbfiles, cn_slchlnanfbfiles,         & 
     264         &            cn_slchlpicfbfiles, cn_schltotfbfiles,          & 
     265         &            cn_slphytotfbfiles, cn_slphydiafbfiles,         & 
     266         &            cn_slphynonfbfiles, cn_sspmfbfiles,             & 
     267         &            cn_sfco2fbfiles, cn_spco2fbfiles,               & 
     268         &            cn_plchltotfbfiles, cn_pchltotfbfiles,          & 
     269         &            cn_pno3fbfiles, cn_psi4fbfiles, cn_ppo4fbfiles, & 
     270         &            cn_pdicfbfiles, cn_palkfbfiles, cn_pphfbfiles,  & 
     271         &            cn_po2fbfiles,                                  & 
     272         &            cn_sstbiasfiles, cn_altbiasfile,                & 
    164273         &            cn_gridsearchfile, rn_gridsearchres,            & 
    165          &            rn_dobsini, rn_dobsend, nn_1dint, nn_2dint,     & 
     274         &            rn_dobsini, rn_dobsend,                         & 
     275         &            rn_default_avglamscl, rn_default_avgphiscl,     & 
     276         &            rn_sla_avglamscl, rn_sla_avgphiscl,             & 
     277         &            rn_sst_avglamscl, rn_sst_avgphiscl,             & 
     278         &            rn_sss_avglamscl, rn_sss_avgphiscl,             & 
     279         &            rn_sic_avglamscl, rn_sic_avgphiscl,             & 
     280         &            nn_1dint, nn_2dint_default,                     & 
     281         &            nn_2dint_sla, nn_2dint_sst,                     & 
     282         &            nn_2dint_sss, nn_2dint_sic,                     & 
    166283         &            nn_msshc, rn_mdtcorr, rn_mdtcutoff,             & 
    167          &            nn_profdavtypes, ln_sstbias, cn_sstbias_files 
     284         &            nn_profdavtypes 
    168285 
    169286      INTEGER :: jnumsstbias 
    170       CALL wrk_alloc( jpi, jpj, zglam1 ) 
    171       CALL wrk_alloc( jpi, jpj, zglam2 ) 
    172       CALL wrk_alloc( jpi, jpj, zgphi1 ) 
    173       CALL wrk_alloc( jpi, jpj, zgphi2 ) 
    174       CALL wrk_alloc( jpi, jpj, jpk, zmask1 ) 
    175       CALL wrk_alloc( jpi, jpj, jpk, zmask2 ) 
     287      ALLOCATE(sstbias_type(jpmaxnfiles)) 
    176288 
    177289      !----------------------------------------------------------------------- 
    178290      ! Read namelist parameters 
    179291      !----------------------------------------------------------------------- 
    180        
    181       !Initalise all values in namelist arrays 
    182       ALLOCATE(sstbias_type(jpmaxnfiles)) 
     292 
    183293      ! Some namelist arrays need initialising 
    184       cn_profbfiles(:) = '' 
    185       cn_slafbfiles(:) = '' 
    186       cn_sstfbfiles(:) = '' 
    187       cn_sicfbfiles(:) = '' 
    188       cn_velfbfiles(:) = '' 
    189       cn_sstbias_files(:) = '' 
    190       nn_profdavtypes(:) = -1 
     294      cn_profbfiles(:)      = '' 
     295      cn_slafbfiles(:)      = '' 
     296      cn_sstfbfiles(:)      = '' 
     297      cn_sicfbfiles(:)      = '' 
     298      cn_velfbfiles(:)      = '' 
     299      cn_sssfbfiles(:)      = '' 
     300      cn_slchltotfbfiles(:) = '' 
     301      cn_slchldiafbfiles(:) = '' 
     302      cn_slchlnonfbfiles(:) = '' 
     303      cn_slchldinfbfiles(:) = '' 
     304      cn_slchlmicfbfiles(:) = '' 
     305      cn_slchlnanfbfiles(:) = '' 
     306      cn_slchlpicfbfiles(:) = '' 
     307      cn_schltotfbfiles(:)  = '' 
     308      cn_slphytotfbfiles(:) = '' 
     309      cn_slphydiafbfiles(:) = '' 
     310      cn_slphynonfbfiles(:) = '' 
     311      cn_sspmfbfiles(:)     = '' 
     312      cn_sfco2fbfiles(:)    = '' 
     313      cn_spco2fbfiles(:)    = '' 
     314      cn_plchltotfbfiles(:) = '' 
     315      cn_pchltotfbfiles(:)  = '' 
     316      cn_pno3fbfiles(:)     = '' 
     317      cn_psi4fbfiles(:)     = '' 
     318      cn_ppo4fbfiles(:)     = '' 
     319      cn_pdicfbfiles(:)     = '' 
     320      cn_palkfbfiles(:)     = '' 
     321      cn_pphfbfiles(:)      = '' 
     322      cn_po2fbfiles(:)      = '' 
     323      cn_sstbiasfiles(:)    = '' 
     324      nn_profdavtypes(:)    = -1 
    191325 
    192326      CALL ini_date( rn_dobsini ) 
    … …  
    205339      IF ( .NOT. ln_diaobs ) THEN 
    206340         IF(lwp) WRITE(numout,cform_war) 
    207          IF(lwp) WRITE(numout,*)' ln_diaobs is set to false so not calling dia_obs' 
     341         IF(lwp) WRITE(numout,*)' ln_diaobs is set to false, so not calling dia_obs' 
    208342         RETURN 
    209       ENDIF 
    210        
    211       !----------------------------------------------------------------------- 
    212       ! Set up list of observation types to be used 
    213       ! and the files associated with each type 
    214       !----------------------------------------------------------------------- 
    215  
    216       nproftypes = COUNT( (/ln_t3d .OR. ln_s3d, ln_vel3d /) ) 
    217       nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic /) ) 
    218  
    219       IF (ln_sstbias) THEN  
    220          lmask(:) = .FALSE.  
    221          WHERE (cn_sstbias_files(:) /= '') lmask(:) = .TRUE.  
    222          jnumsstbias = COUNT(lmask)  
    223          lmask(:) = .FALSE.  
    224       ENDIF       
    225  
    226       IF ( nproftypes == 0 .AND. nsurftypes == 0 ) THEN 
    227          IF(lwp) WRITE(numout,cform_war) 
    228          IF(lwp) WRITE(numout,*) ' ln_diaobs is set to true, but all obs operator logical flags', & 
    229             &                    ' ln_t3d, ln_s3d, ln_sla, ln_sst, ln_sic, ln_vel3d', & 
    230             &                    ' are set to .FALSE. so turning off calls to dia_obs' 
    231          nwarn = nwarn + 1 
    232          ln_diaobs = .FALSE. 
    233          RETURN 
    234       ENDIF 
    235  
    236       IF ( nproftypes > 0 ) THEN 
    237  
    238          ALLOCATE( cobstypesprof(nproftypes) ) 
    239          ALLOCATE( ifilesprof(nproftypes) ) 
    240          ALLOCATE( clproffiles(nproftypes,jpmaxnfiles) ) 
    241  
    242          jtype = 0 
    243          IF (ln_t3d .OR. ln_s3d) THEN 
    244             jtype = jtype + 1 
    245             clproffiles(jtype,:) = cn_profbfiles(:) 
    246             cobstypesprof(jtype) = 'prof  ' 
    247             ifilesprof(jtype) = 0 
    248             DO jfile = 1, jpmaxnfiles 
    249                IF ( trim(clproffiles(jtype,jfile)) /= '' ) & 
    250                   ifilesprof(jtype) = ifilesprof(jtype) + 1 
    251             END DO 
    252          ENDIF 
    253          IF (ln_vel3d) THEN 
    254             jtype = jtype + 1 
    255             clproffiles(jtype,:) = cn_velfbfiles(:) 
    256             cobstypesprof(jtype) = 'vel   ' 
    257             ifilesprof(jtype) = 0 
    258             DO jfile = 1, jpmaxnfiles 
    259                IF ( trim(clproffiles(jtype,jfile)) /= '' ) & 
    260                   ifilesprof(jtype) = ifilesprof(jtype) + 1 
    261             END DO 
    262          ENDIF 
    263  
    264       ENDIF 
    265  
    266       IF ( nsurftypes > 0 ) THEN 
    267  
    268          ALLOCATE( cobstypessurf(nsurftypes) ) 
    269          ALLOCATE( ifilessurf(nsurftypes) ) 
    270          ALLOCATE( clsurffiles(nsurftypes, jpmaxnfiles) ) 
    271  
    272          jtype = 0 
    273          IF (ln_sla) THEN 
    274             jtype = jtype + 1 
    275             clsurffiles(jtype,:) = cn_slafbfiles(:) 
    276             cobstypessurf(jtype) = 'sla   ' 
    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          IF (ln_sst) THEN 
    284             jtype = jtype + 1 
    285             clsurffiles(jtype,:) = cn_sstfbfiles(:) 
    286             cobstypessurf(jtype) = 'sst   ' 
    287             ifilessurf(jtype) = 0 
    288             DO jfile = 1, jpmaxnfiles 
    289                IF ( trim(clsurffiles(jtype,jfile)) /= '' ) & 
    290                   ifilessurf(jtype) = ifilessurf(jtype) + 1 
    291             END DO 
    292          ENDIF 
    293 #if defined key_lim2 || defined key_lim3 
    294          IF (ln_sic) THEN 
    295             jtype = jtype + 1 
    296             clsurffiles(jtype,:) = cn_sicfbfiles(:) 
    297             cobstypessurf(jtype) = 'sic   ' 
    298             ifilessurf(jtype) = 0 
    299             DO jfile = 1, jpmaxnfiles 
    300                IF ( trim(clsurffiles(jtype,jfile)) /= '' ) & 
    301                   ifilessurf(jtype) = ifilessurf(jtype) + 1 
    302             END DO 
    303          ENDIF 
    304 #endif 
    305  
    306343      ENDIF 
    307344 
    … …  
    318355         WRITE(numout,*) '             Logical switch for Sea Ice observations                  ln_sic = ', ln_sic 
    319356         WRITE(numout,*) '             Logical switch for velocity observations               ln_vel3d = ', ln_vel3d 
    320          WRITE(numout,*) '             Global distribution of observations              ln_grid_global = ',ln_grid_global 
    321          WRITE(numout,*) '             Logical switch for SST bias correction         ln_sstbias = ', ln_sstbias  
    322          WRITE(numout,*) '             Logical switch for obs grid search lookup ln_grid_search_lookup = ',ln_grid_search_lookup 
     357         WRITE(numout,*) '             Logical switch for SSS observations                      ln_sss = ', ln_sss 
     358         WRITE(numout,*) '             Logical switch for surface total logchl obs         ln_slchltot = ', ln_slchltot 
     359         WRITE(numout,*) '             Logical switch for surface diatom logchl obs        ln_slchldia = ', ln_slchldia 
     360         WRITE(numout,*) '             Logical switch for surface non-diatom logchl obs    ln_slchlnon = ', ln_slchlnon 
     361         WRITE(numout,*) '             Logical switch for surface dino logchl obs          ln_slchldin = ', ln_slchldin 
     362         WRITE(numout,*) '             Logical switch for surface micro logchl obs         ln_slchlmic = ', ln_slchlmic 
     363         WRITE(numout,*) '             Logical switch for surface nano logchl obs          ln_slchlnan = ', ln_slchlnan 
     364         WRITE(numout,*) '             Logical switch for surface pico logchl obs          ln_slchlpic = ', ln_slchlpic 
     365         WRITE(numout,*) '             Logical switch for surface total chl obs             ln_schltot = ', ln_schltot 
     366         WRITE(numout,*) '             Logical switch for surface total log(phyC) obs      ln_slphytot = ', ln_slphytot 
     367         WRITE(numout,*) '             Logical switch for surface diatom log(phyC) obs     ln_slphydia = ', ln_slphydia 
     368         WRITE(numout,*) '             Logical switch for surface non-diatom log(phyC) obs ln_slphynon = ', ln_slphynon 
     369         WRITE(numout,*) '             Logical switch for surface SPM observations             ln_sspm = ', ln_sspm 
     370         WRITE(numout,*) '             Logical switch for surface fCO2 observations           ln_sfco2 = ', ln_sfco2 
     371         WRITE(numout,*) '             Logical switch for surface pCO2 observations           ln_spco2 = ', ln_spco2 
     372         WRITE(numout,*) '             Logical switch for profile total logchl obs         ln_plchltot = ', ln_plchltot 
     373         WRITE(numout,*) '             Logical switch for profile total chl obs             ln_pchltot = ', ln_pchltot 
     374         WRITE(numout,*) '             Logical switch for profile nitrate obs                  ln_pno3 = ', ln_pno3 
     375         WRITE(numout,*) '             Logical switch for profile silicate obs                 ln_psi4 = ', ln_psi4 
     376         WRITE(numout,*) '             Logical switch for profile phosphate obs                ln_ppo4 = ', ln_ppo4 
     377         WRITE(numout,*) '             Logical switch for profile DIC obs                      ln_pdic = ', ln_pdic 
     378         WRITE(numout,*) '             Logical switch for profile alkalinity obs               ln_palk = ', ln_palk 
     379         WRITE(numout,*) '             Logical switch for profile pH obs                        ln_pph = ', ln_pph 
     380         WRITE(numout,*) '             Logical switch for profile oxygen obs                    ln_po2 = ', ln_po2 
     381         WRITE(numout,*) '             Global distribution of observations              ln_grid_global = ', ln_grid_global 
     382         WRITE(numout,*) '             Logical switch for obs grid search lookup ln_grid_search_lookup = ', ln_grid_search_lookup 
    323383         IF (ln_grid_search_lookup) & 
    324384            WRITE(numout,*) '             Grid search lookup file header                cn_gridsearchfile = ', cn_gridsearchfile 
    … …  
    326386         WRITE(numout,*) '             Final date in window YYYYMMDD.HHMMSS                 rn_dobsend = ', rn_dobsend 
    327387         WRITE(numout,*) '             Type of vertical interpolation method                  nn_1dint = ', nn_1dint 
    328          WRITE(numout,*) '             Type of horizontal interpolation method                nn_2dint = ', nn_2dint 
     388         WRITE(numout,*) '             Default horizontal interpolation method        nn_2dint_default = ', nn_2dint_default 
     389         WRITE(numout,*) '             Type of horizontal interpolation method for SLA    nn_2dint_sla = ', nn_2dint_sla 
     390         WRITE(numout,*) '             Type of horizontal interpolation method for SST    nn_2dint_sst = ', nn_2dint_sst 
     391         WRITE(numout,*) '             Type of horizontal interpolation method for SSS    nn_2dint_sss = ', nn_2dint_sss 
     392         WRITE(numout,*) '             Type of horizontal interpolation method for SIC    nn_2dint_sic = ', nn_2dint_sic 
     393         WRITE(numout,*) '             Default E/W diameter of obs footprint      rn_default_avglamscl = ', rn_default_avglamscl 
     394         WRITE(numout,*) '             Default N/S diameter of obs footprint      rn_default_avgphiscl = ', rn_default_avgphiscl 
     395         WRITE(numout,*) '             Default obs footprint in deg [T] or m [F]  ln_default_fp_indegs = ', ln_default_fp_indegs 
     396         WRITE(numout,*) '             SLA E/W diameter of obs footprint              rn_sla_avglamscl = ', rn_sla_avglamscl 
     397         WRITE(numout,*) '             SLA N/S diameter of obs footprint              rn_sla_avgphiscl = ', rn_sla_avgphiscl 
     398         WRITE(numout,*) '             SLA obs footprint in deg [T] or m [F]          ln_sla_fp_indegs = ', ln_sla_fp_indegs 
     399         WRITE(numout,*) '             SST E/W diameter of obs footprint              rn_sst_avglamscl = ', rn_sst_avglamscl 
     400         WRITE(numout,*) '             SST N/S diameter of obs footprint              rn_sst_avgphiscl = ', rn_sst_avgphiscl 
     401         WRITE(numout,*) '             SST obs footprint in deg [T] or m [F]          ln_sst_fp_indegs = ', ln_sst_fp_indegs 
     402         WRITE(numout,*) '             SIC E/W diameter of obs footprint              rn_sic_avglamscl = ', rn_sic_avglamscl 
     403         WRITE(numout,*) '             SIC N/S diameter of obs footprint              rn_sic_avgphiscl = ', rn_sic_avgphiscl 
     404         WRITE(numout,*) '             SIC obs footprint in deg [T] or m [F]          ln_sic_fp_indegs = ', ln_sic_fp_indegs 
    329405         WRITE(numout,*) '             Rejection of observations near land switch               ln_nea = ', ln_nea 
     406         WRITE(numout,*) '             Rejection of obs near open bdys                 ln_bound_reject = ', ln_bound_reject 
    330407         WRITE(numout,*) '             MSSH correction scheme                                 nn_msshc = ', nn_msshc 
    331408         WRITE(numout,*) '             MDT  correction                                      rn_mdtcorr = ', rn_mdtcorr 
    332409         WRITE(numout,*) '             MDT cutoff for computed correction                 rn_mdtcutoff = ', rn_mdtcutoff 
    333410         WRITE(numout,*) '             Logical switch for alt bias                          ln_altbias = ', ln_altbias 
     411         WRITE(numout,*) '             Logical switch for sst bias                          ln_sstbias = ', ln_sstbias 
    334412         WRITE(numout,*) '             Logical switch for ignoring missing files             ln_ignmis = ', ln_ignmis 
    335413         WRITE(numout,*) '             Daily average types                             nn_profdavtypes = ', nn_profdavtypes 
    336414         WRITE(numout,*) '             Logical switch for night-time SST obs               ln_sstnight = ', ln_sstnight 
    337          WRITE(numout,*) '          Number of profile obs types: ',nproftypes 
    338  
    339          IF ( nproftypes > 0 ) THEN 
    340             DO jtype = 1, nproftypes 
    341                DO jfile = 1, ifilesprof(jtype) 
    342                   WRITE(numout,'(1X,2A)') '             '//cobstypesprof(jtype)//' input observation file names  = ', & 
    343                      TRIM(clproffiles(jtype,jfile)) 
    344                END DO 
    345             END DO 
    346          ENDIF 
    347  
    348          WRITE(numout,*)'          Number of surface obs types: ',nsurftypes 
    349          IF ( nsurftypes > 0 ) THEN 
    350             DO jtype = 1, nsurftypes 
    351                DO jfile = 1, ifilessurf(jtype) 
    352                   WRITE(numout,'(1X,2A)') '             '//cobstypessurf(jtype)//' input observation file names  = ', & 
    353                      TRIM(clsurffiles(jtype,jfile)) 
    354                END DO 
    355             END DO 
    356          ENDIF 
    357          WRITE(numout,*) '~~~~~~~~~~~~' 
    358  
    359       ENDIF 
     415      ENDIF 
     416      !----------------------------------------------------------------------- 
     417      ! Set up list of observation types to be used 
     418      ! and the files associated with each type 
     419      !----------------------------------------------------------------------- 
     420 
     421      nproftypes = COUNT( (/ln_t3d .OR. ln_s3d, ln_vel3d, ln_plchltot,          & 
     422         &                  ln_pchltot,  ln_pno3,     ln_psi4,     ln_ppo4,     & 
     423         &                  ln_pdic,     ln_palk,     ln_pph,      ln_po2 /) ) 
     424      nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic, ln_sss,                     & 
     425         &                  ln_slchltot, ln_slchldia, ln_slchlnon, ln_slchldin, & 
     426         &                  ln_slchlmic, ln_slchlnan, ln_slchlpic, ln_schltot,  & 
     427         &                  ln_slphytot, ln_slphydia, ln_slphynon, ln_sspm,     & 
     428         &                  ln_sfco2,    ln_spco2 /) ) 
     429 
     430      IF (ln_sstbias) THEN  
     431         lmask(:) = .FALSE.  
     432         WHERE (cn_sstbiasfiles(:) /= '') lmask(:) = .TRUE.  
     433         jnumsstbias = COUNT(lmask)  
     434         lmask(:) = .FALSE.  
     435      ENDIF       
     436 
     437      IF ( nproftypes == 0 .AND. nsurftypes == 0 ) THEN 
     438         IF(lwp) WRITE(numout,cform_war) 
     439         IF(lwp) WRITE(numout,*) ' ln_diaobs is set to true, but all obs operator logical flags', & 
     440            &                    ' ln_t3d, ln_s3d, ln_sla, ln_sst, ln_sic, ln_vel3d', & 
     441            &                    ' are set to .FALSE. so turning off calls to dia_obs' 
     442         nwarn = nwarn + 1 
     443         ln_diaobs = .FALSE. 
     444         RETURN 
     445      ENDIF 
     446 
     447      IF(lwp) WRITE(numout,*) '          Number of profile obs types: ',nproftypes 
     448      IF ( nproftypes > 0 ) THEN 
     449 
     450         ALLOCATE( cobstypesprof(nproftypes) ) 
     451         ALLOCATE( ifilesprof(nproftypes) ) 
     452         ALLOCATE( clproffiles(nproftypes,jpmaxnfiles) ) 
     453 
     454         jtype = 0 
     455         IF (ln_t3d .OR. ln_s3d) THEN 
     456            jtype = jtype + 1 
     457            cobstypesprof(jtype) = 'prof' 
     458            clproffiles(jtype,:) = cn_profbfiles 
     459         ENDIF 
     460         IF (ln_vel3d) THEN 
     461            jtype = jtype + 1 
     462            cobstypesprof(jtype) =  'vel' 
     463            clproffiles(jtype,:) = cn_velfbfiles 
     464         ENDIF 
     465         IF (ln_plchltot) THEN 
     466            jtype = jtype + 1 
     467            cobstypesprof(jtype) = 'plchltot' 
     468            clproffiles(jtype,:) = cn_plchltotfbfiles 
     469         ENDIF 
     470         IF (ln_pchltot) THEN 
     471            jtype = jtype + 1 
     472            cobstypesprof(jtype) = 'pchltot' 
     473            clproffiles(jtype,:) = cn_pchltotfbfiles 
     474         ENDIF 
     475         IF (ln_pno3) THEN 
     476            jtype = jtype + 1 
     477            cobstypesprof(jtype) = 'pno3' 
     478            clproffiles(jtype,:) = cn_pno3fbfiles 
     479         ENDIF 
     480         IF (ln_psi4) THEN 
     481            jtype = jtype + 1 
     482            cobstypesprof(jtype) = 'psi4' 
     483            clproffiles(jtype,:) = cn_psi4fbfiles 
     484         ENDIF 
     485         IF (ln_ppo4) THEN 
     486            jtype = jtype + 1 
     487            cobstypesprof(jtype) = 'ppo4' 
     488            clproffiles(jtype,:) = cn_ppo4fbfiles 
     489         ENDIF 
     490         IF (ln_pdic) THEN 
     491            jtype = jtype + 1 
     492            cobstypesprof(jtype) = 'pdic' 
     493            clproffiles(jtype,:) = cn_pdicfbfiles 
     494         ENDIF 
     495         IF (ln_palk) THEN 
     496            jtype = jtype + 1 
     497            cobstypesprof(jtype) = 'palk' 
     498            clproffiles(jtype,:) = cn_palkfbfiles 
     499         ENDIF 
     500         IF (ln_pph) THEN 
     501            jtype = jtype + 1 
     502            cobstypesprof(jtype) = 'pph' 
     503            clproffiles(jtype,:) = cn_pphfbfiles 
     504         ENDIF 
     505         IF (ln_po2) THEN 
     506            jtype = jtype + 1 
     507            cobstypesprof(jtype) = 'po2' 
     508            clproffiles(jtype,:) = cn_po2fbfiles 
     509         ENDIF 
     510 
     511         CALL obs_settypefiles( nproftypes, jpmaxnfiles, ifilesprof, cobstypesprof, clproffiles ) 
     512 
     513      ENDIF 
     514 
     515      IF(lwp) WRITE(numout,*)'          Number of surface obs types: ',nsurftypes 
     516      IF ( nsurftypes > 0 ) THEN 
     517 
     518         ALLOCATE( cobstypessurf(nsurftypes) ) 
     519         ALLOCATE( ifilessurf(nsurftypes) ) 
     520         ALLOCATE( clsurffiles(nsurftypes, jpmaxnfiles) ) 
     521         ALLOCATE(n2dintsurf(nsurftypes)) 
     522         ALLOCATE(ravglamscl(nsurftypes)) 
     523         ALLOCATE(ravgphiscl(nsurftypes)) 
     524         ALLOCATE(lfpindegs(nsurftypes)) 
     525         ALLOCATE(llnightav(nsurftypes)) 
     526 
     527         jtype = 0 
     528         IF (ln_sla) THEN 
     529            jtype = jtype + 1 
     530            cobstypessurf(jtype) = 'sla' 
     531            clsurffiles(jtype,:) = cn_slafbfiles 
     532         ENDIF 
     533         IF (ln_sst) THEN 
     534            jtype = jtype + 1 
     535            cobstypessurf(jtype) = 'sst' 
     536            clsurffiles(jtype,:) = cn_sstfbfiles 
     537         ENDIF 
     538#if defined key_lim2 || defined key_lim3 
     539         IF (ln_sic) THEN 
     540            jtype = jtype + 1 
     541            cobstypessurf(jtype) = 'sic' 
     542            clsurffiles(jtype,:) = cn_sicfbfiles 
     543         ENDIF 
     544#endif 
     545         IF (ln_sss) THEN 
     546            jtype = jtype + 1 
     547            cobstypessurf(jtype) = 'sss' 
     548            clsurffiles(jtype,:) = cn_sssfbfiles 
     549         ENDIF 
     550         IF (ln_slchltot) THEN 
     551            jtype = jtype + 1 
     552            cobstypessurf(jtype) = 'slchltot' 
     553            clsurffiles(jtype,:) = cn_slchltotfbfiles 
     554         ENDIF 
     555         IF (ln_slchldia) THEN 
     556            jtype = jtype + 1 
     557            cobstypessurf(jtype) = 'slchldia' 
     558            clsurffiles(jtype,:) = cn_slchldiafbfiles 
     559         ENDIF 
     560         IF (ln_slchlnon) THEN 
     561            jtype = jtype + 1 
     562            cobstypessurf(jtype) = 'slchlnon' 
     563            clsurffiles(jtype,:) = cn_slchlnonfbfiles 
     564         ENDIF 
     565         IF (ln_slchldin) THEN 
     566            jtype = jtype + 1 
     567            cobstypessurf(jtype) = 'slchldin' 
     568            clsurffiles(jtype,:) = cn_slchldinfbfiles 
     569         ENDIF 
     570         IF (ln_slchlmic) THEN 
     571            jtype = jtype + 1 
     572            cobstypessurf(jtype) = 'slchlmic' 
     573            clsurffiles(jtype,:) = cn_slchlmicfbfiles 
     574         ENDIF 
     575         IF (ln_slchlnan) THEN 
     576            jtype = jtype + 1 
     577            cobstypessurf(jtype) = 'slchlnan' 
     578            clsurffiles(jtype,:) = cn_slchlnanfbfiles 
     579         ENDIF 
     580         IF (ln_slchlpic) THEN 
     581            jtype = jtype + 1 
     582            cobstypessurf(jtype) = 'slchlpic' 
     583            clsurffiles(jtype,:) = cn_slchlpicfbfiles 
     584         ENDIF 
     585         IF (ln_schltot) THEN 
     586            jtype = jtype + 1 
     587            cobstypessurf(jtype) = 'schltot' 
     588            clsurffiles(jtype,:) = cn_schltotfbfiles 
     589         ENDIF 
     590         IF (ln_slphytot) THEN 
     591            jtype = jtype + 1 
     592            cobstypessurf(jtype) = 'slphytot' 
     593            clsurffiles(jtype,:) = cn_slphytotfbfiles 
     594         ENDIF 
     595         IF (ln_slphydia) THEN 
     596            jtype = jtype + 1 
     597            cobstypessurf(jtype) = 'slphydia' 
     598            clsurffiles(jtype,:) = cn_slphydiafbfiles 
     599         ENDIF 
     600         IF (ln_slphynon) THEN 
     601            jtype = jtype + 1 
     602            cobstypessurf(jtype) = 'slphynon' 
     603            clsurffiles(jtype,:) = cn_slphynonfbfiles 
     604         ENDIF 
     605         IF (ln_sspm) THEN 
     606            jtype = jtype + 1 
     607            cobstypessurf(jtype) = 'sspm' 
     608            clsurffiles(jtype,:) = cn_sspmfbfiles 
     609         ENDIF 
     610         IF (ln_sfco2) THEN 
     611            jtype = jtype + 1 
     612            cobstypessurf(jtype) = 'sfco2' 
     613            clsurffiles(jtype,:) = cn_sfco2fbfiles 
     614         ENDIF 
     615         IF (ln_spco2) THEN 
     616            jtype = jtype + 1 
     617            cobstypessurf(jtype) = 'spco2' 
     618            clsurffiles(jtype,:) = cn_spco2fbfiles 
     619         ENDIF 
     620 
     621         CALL obs_settypefiles( nsurftypes, jpmaxnfiles, ifilessurf, cobstypessurf, clsurffiles ) 
     622 
     623         DO jtype = 1, nsurftypes 
     624 
     625            IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN 
     626               IF ( nn_2dint_sla == -1 ) THEN 
     627                  n2dint_type  = nn_2dint_default 
     628               ELSE 
     629                  n2dint_type  = nn_2dint_sla 
     630               ENDIF 
     631               ztype_avglamscl = rn_sla_avglamscl 
     632               ztype_avgphiscl = rn_sla_avgphiscl 
     633               ltype_fp_indegs = ln_sla_fp_indegs 
     634               ltype_night     = .FALSE. 
     635            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sst' ) THEN 
     636               IF ( nn_2dint_sst == -1 ) THEN 
     637                  n2dint_type  = nn_2dint_default 
     638               ELSE 
     639                  n2dint_type  = nn_2dint_sst 
     640               ENDIF 
     641               ztype_avglamscl = rn_sst_avglamscl 
     642               ztype_avgphiscl = rn_sst_avgphiscl 
     643               ltype_fp_indegs = ln_sst_fp_indegs 
     644               ltype_night     = ln_sstnight 
     645            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sic' ) THEN 
     646               IF ( nn_2dint_sic == -1 ) THEN 
     647                  n2dint_type  = nn_2dint_default 
     648               ELSE 
     649                  n2dint_type  = nn_2dint_sic 
     650               ENDIF 
     651               ztype_avglamscl = rn_sic_avglamscl 
     652               ztype_avgphiscl = rn_sic_avgphiscl 
     653               ltype_fp_indegs = ln_sic_fp_indegs 
     654               ltype_night     = .FALSE. 
     655            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sss' ) THEN 
     656               IF ( nn_2dint_sss == -1 ) THEN 
     657                  n2dint_type  = nn_2dint_default 
     658               ELSE 
     659                  n2dint_type  = nn_2dint_sss 
     660               ENDIF 
     661               ztype_avglamscl = rn_sss_avglamscl 
     662               ztype_avgphiscl = rn_sss_avgphiscl 
     663               ltype_fp_indegs = ln_sss_fp_indegs 
     664               ltype_night     = .FALSE. 
     665            ELSE 
     666               n2dint_type     = nn_2dint_default 
     667               ztype_avglamscl = rn_default_avglamscl 
     668               ztype_avgphiscl = rn_default_avgphiscl 
     669               ltype_fp_indegs = ln_default_fp_indegs 
     670               ltype_night     = .FALSE. 
     671            ENDIF 
     672             
     673            CALL obs_setinterpopts( nsurftypes, jtype, TRIM(cobstypessurf(jtype)), & 
     674               &                    nn_2dint_default, n2dint_type,                 & 
     675               &                    ztype_avglamscl, ztype_avgphiscl,              & 
     676               &                    ltype_fp_indegs, ltype_night,                  & 
     677               &                    n2dintsurf, ravglamscl, ravgphiscl,            & 
     678               &                    lfpindegs, llnightav ) 
     679 
     680         END DO 
     681 
     682      ENDIF 
     683 
     684      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 
     685 
    360686 
    361687      !----------------------------------------------------------------------- 
    … …  
    377703      ENDIF 
    378704 
    379       IF ( ( nn_2dint < 0 ) .OR. ( nn_2dint > 4 ) ) THEN 
    380          CALL ctl_stop(' Choice of horizontal (2D) interpolation method', & 
     705      IF ( ( nn_2dint_default < 0 ) .OR. ( nn_2dint_default > 6 ) ) THEN 
     706         CALL ctl_stop(' Choice of default horizontal (2D) interpolation method', & 
    381707            &                    ' is not available') 
    382708      ENDIF 
    … …  
    402728         DO jtype = 1, nproftypes 
    403729 
    404             nvarsprof(jtype) = 2 
    405730            IF ( TRIM(cobstypesprof(jtype)) == 'prof' ) THEN 
     731               nvarsprof(jtype) = 2 
    406732               nextrprof(jtype) = 1 
    407                llvar1 = ln_t3d 
    408                llvar2 = ln_s3d 
    409                zglam1 = glamt 
    410                zgphi1 = gphit 
    411                zmask1 = tmask 
    412                zglam2 = glamt 
    413                zgphi2 = gphit 
    414                zmask2 = tmask 
    415             ENDIF 
    416             IF ( TRIM(cobstypesprof(jtype)) == 'vel' )  THEN 
     733               ALLOCATE(llvar(nvarsprof(jtype))) 
     734               CALL wrk_alloc( jpi, jpj,      nvarsprof(jtype), zglam ) 
     735               CALL wrk_alloc( jpi, jpj,      nvarsprof(jtype), zgphi ) 
     736               CALL wrk_alloc( jpi, jpj, jpk, nvarsprof(jtype), zmask ) 
     737               llvar(1)       = ln_t3d 
     738               llvar(2)       = ln_s3d 
     739               zglam(:,:,1)   = glamt(:,:) 
     740               zglam(:,:,2)   = glamt(:,:) 
     741               zgphi(:,:,1)   = gphit(:,:) 
     742               zgphi(:,:,2)   = gphit(:,:) 
     743               zmask(:,:,:,1) = tmask(:,:,:) 
     744               zmask(:,:,:,2) = tmask(:,:,:) 
     745            ELSE IF ( TRIM(cobstypesprof(jtype)) == 'vel' )  THEN 
     746               nvarsprof(jtype) = 2 
    417747               nextrprof(jtype) = 2 
    418                llvar1 = ln_vel3d 
    419                llvar2 = ln_vel3d 
    420                zglam1 = glamu 
    421                zgphi1 = gphiu 
    422                zmask1 = umask 
    423                zglam2 = glamv 
    424                zgphi2 = gphiv 
    425                zmask2 = vmask 
     748               ALLOCATE(llvar(nvarsprof(jtype))) 
     749               CALL wrk_alloc( jpi, jpj,      nvarsprof(jtype), zglam ) 
     750               CALL wrk_alloc( jpi, jpj,      nvarsprof(jtype), zgphi ) 
     751               CALL wrk_alloc( jpi, jpj, jpk, nvarsprof(jtype), zmask ) 
     752               llvar(1)       = ln_vel3d 
     753               llvar(2)       = ln_vel3d 
     754               zglam(:,:,1)   = glamu(:,:) 
     755               zglam(:,:,2)   = glamv(:,:) 
     756               zgphi(:,:,1)   = gphiu(:,:) 
     757               zgphi(:,:,2)   = gphiv(:,:) 
     758               zmask(:,:,:,1) = umask(:,:,:) 
     759               zmask(:,:,:,2) = vmask(:,:,:) 
     760            ELSE 
     761               nvarsprof(jtype) = 1 
     762               nextrprof(jtype) = 0 
     763               ALLOCATE(llvar(nvarsprof(jtype))) 
     764               CALL wrk_alloc( jpi, jpj,      nvarsprof(jtype), zglam ) 
     765               CALL wrk_alloc( jpi, jpj,      nvarsprof(jtype), zgphi ) 
     766               CALL wrk_alloc( jpi, jpj, jpk, nvarsprof(jtype), zmask ) 
     767               llvar(1)       = .TRUE. 
     768               zglam(:,:,1)   = glamt(:,:) 
     769               zgphi(:,:,1)   = gphit(:,:) 
     770               zmask(:,:,:,1) = tmask(:,:,:) 
    426771            ENDIF 
    427772 
    … …  
    430775               &               clproffiles(jtype,1:ifilesprof(jtype)), & 
    431776               &               nvarsprof(jtype), nextrprof(jtype), nitend-nit000+2, & 
    432                &               rn_dobsini, rn_dobsend, llvar1, llvar2, & 
     777               &               rn_dobsini, rn_dobsend, llvar, & 
    433778               &               ln_ignmis, ln_s_at_t, .FALSE., & 
    434779               &               kdailyavtypes = nn_profdavtypes ) 
    … …  
    439784 
    440785            CALL obs_pre_prof( profdata(jtype), profdataqc(jtype), & 
    441                &               llvar1, llvar2, & 
     786               &               llvar, & 
    442787               &               jpi, jpj, jpk, & 
    443                &               zmask1, zglam1, zgphi1, zmask2, zglam2, zgphi2,  & 
    444                &               ln_nea, kdailyavtypes = nn_profdavtypes ) 
     788               &               zmask, zglam, zgphi,  & 
     789               &               ln_nea, ln_bound_reject, & 
     790               &               kdailyavtypes = nn_profdavtypes ) 
     791             
     792            DEALLOCATE( llvar ) 
     793            CALL wrk_dealloc( jpi, jpj,      nvarsprof(jtype), zglam ) 
     794            CALL wrk_dealloc( jpi, jpj,      nvarsprof(jtype), zgphi ) 
     795            CALL wrk_dealloc( jpi, jpj, jpk, nvarsprof(jtype), zmask ) 
    445796 
    446797         END DO 
    … …  
    461812            nvarssurf(jtype) = 1 
    462813            nextrsurf(jtype) = 0 
    463             llnightav = .FALSE. 
    464814            IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) nextrsurf(jtype) = 2 
    465             IF ( TRIM(cobstypessurf(jtype)) == 'sst' ) llnightav = ln_sstnight 
     815            IF ( TRIM(cobstypessurf(jtype)) == 'sst' ) llnightav(jtype) = ln_sstnight 
    466816 
    467817            !Read in surface obs types 
    … …  
    469819               &               clsurffiles(jtype,1:ifilessurf(jtype)), & 
    470820               &               nvarssurf(jtype), nextrsurf(jtype), nitend-nit000+2, & 
    471                &               rn_dobsini, rn_dobsend, ln_ignmis, .FALSE., llnightav ) 
    472           
    473           
    474             CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea ) 
     821               &               rn_dobsini, rn_dobsend, ln_ignmis, .FALSE., llnightav(jtype) ) 
     822 
     823            CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea, ln_bound_reject ) 
    475824 
    476825            IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN 
    477                CALL obs_rea_mdt( surfdataqc(jtype), nn_2dint ) 
    478                IF ( ln_altbias ) CALL obs_rea_altbias ( surfdataqc(jtype), nn_2dint, cn_altbiasfile ) 
     826               CALL obs_rea_mdt( surfdataqc(jtype), n2dintsurf(jtype) ) 
     827               IF ( ln_altbias ) & 
     828                  & CALL obs_rea_altbias ( surfdataqc(jtype), n2dintsurf(jtype), cn_altbiasfile ) 
    479829            ENDIF 
    480830            IF ( TRIM(cobstypessurf(jtype)) == 'sst' .AND. ln_sstbias ) THEN 
    481                !Read in bias field and correct SST. 
    482                IF ( jnumsstbias == 0 ) CALL ctl_stop("ln_sstbias set,"// & 
    483                                                      "  but no bias"// & 
    484                                                      " files to read in")    
    485                   CALL obs_app_sstbias( surfdataqc(jtype), nn_2dint, & 
    486                                         jnumsstbias, cn_sstbias_files(1:jnumsstbias) ) 
     831               jnumsstbias = 0 
     832               DO jfile = 1, jpmaxnfiles 
     833                  IF ( TRIM(cn_sstbiasfiles(jfile)) /= '' ) & 
     834                     &  jnumsstbias = jnumsstbias + 1 
     835               END DO 
     836               IF ( jnumsstbias == 0 ) THEN 
     837                  CALL ctl_stop("ln_sstbias set but no bias files to read in")     
     838               ENDIF 
     839 
     840               CALL obs_app_sstbias( surfdataqc(jtype), n2dintsurf(jtype), &  
     841                  &                  jnumsstbias, cn_sstbiasfiles(1:jnumsstbias) )  
     842 
    487843            ENDIF 
    488844         END DO 
    … …  
    491847 
    492848      ENDIF 
    493  
    494       CALL wrk_dealloc( jpi, jpj, zglam1 ) 
    495       CALL wrk_dealloc( jpi, jpj, zglam2 ) 
    496       CALL wrk_dealloc( jpi, jpj, zgphi1 ) 
    497       CALL wrk_dealloc( jpi, jpj, zgphi2 ) 
    498       CALL wrk_dealloc( jpi, jpj, jpk, zmask1 ) 
    499       CALL wrk_dealloc( jpi, jpj, jpk, zmask2 ) 
    500849 
    501850   END SUBROUTINE dia_obs_init 
    … …  
    526875      !!---------------------------------------------------------------------- 
    527876      !! * Modules used 
    528       USE dom_oce, ONLY : &             ! Ocean space and time domain variables 
    529          & gdept_n,       &       
    530          & gdept_1d       
    531       USE phycst, ONLY : &              ! Physical constants 
    532          & rday                          
    533       USE oce, ONLY : &                 ! Ocean dynamics and tracers variables 
    534          & tsn,  &              
    535          & un, vn, & 
    536          & sshn   
    537877      USE phycst, ONLY : &         ! Physical constants 
     878#if defined key_fabm 
     879         & rt0,          & 
     880#endif 
    538881         & rday 
     882      USE oce, ONLY : &            ! Ocean dynamics and tracers variables 
     883         & tsn,       & 
     884         & un,        & 
     885         & vn,        & 
     886         & sshn 
    539887#if defined  key_lim3 
    540888      USE ice, ONLY : &            ! LIM3 Ice model variables 
    … …  
    545893         & frld 
    546894#endif 
     895#if defined key_cice 
     896      USE sbc_oce, ONLY : fr_i     ! ice fraction 
     897#endif 
     898#if defined key_top 
     899      USE trc, ONLY :  &           ! Biogeochemical state variables 
     900         & trn 
     901#endif 
     902#if defined key_hadocc 
     903      USE par_hadocc               ! HadOCC parameters 
     904      USE trc, ONLY :  & 
     905         & HADOCC_CHL, & 
     906         & HADOCC_FCO2, & 
     907         & HADOCC_PCO2, & 
     908         & HADOCC_FILL_FLT 
     909      USE had_bgc_const, ONLY: c2n_p 
     910#elif defined key_medusa 
     911      USE par_medusa               ! MEDUSA parameters 
     912      USE sms_medusa, ONLY: & 
     913         & xthetapn, & 
     914         & xthetapd 
     915#if defined key_roam 
     916      USE sms_medusa, ONLY: & 
     917         & f2_pco2w, & 
     918         & f2_fco2w, & 
     919         & f3_pH 
     920#endif 
     921#elif defined key_fabm 
     922      USE par_fabm                 ! FABM parameters 
     923      USE fabm, ONLY: & 
     924         & fabm_get_interior_diagnostic_data 
     925#endif 
     926#if defined key_spm 
     927      USE par_spm, ONLY: &         ! Sediment parameters 
     928         & jp_spm 
     929#endif 
     930 
    547931      IMPLICIT NONE 
    548932 
    … …  
    553937      INTEGER :: jtype             ! Data loop variable 
    554938      INTEGER :: jvar              ! Variable number 
    555       INTEGER :: ji, jj            ! Loop counters 
     939      INTEGER :: ji, jj, jk        ! Loop counters 
     940      REAL(wp) :: tiny             ! small number 
     941      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: & 
     942         & zprofvar                ! Model values for variables in a prof ob 
     943      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: & 
     944         & zprofmask               ! Mask associated with zprofvar 
     945      REAL(wp), POINTER, DIMENSION(:,:) :: & 
     946         & zsurfvar, &             ! Model values equivalent to surface ob. 
     947         & zsurfmask               ! Mask associated with surface variable 
    556948      REAL(wp), POINTER, DIMENSION(:,:,:) :: & 
    557          & zprofvar1, &            ! Model values for 1st variable in a prof ob 
    558          & zprofvar2               ! Model values for 2nd variable in a prof ob 
     949         & zglam,    &             ! Model longitudes for prof variables 
     950         & zgphi                   ! Model latitudes for prof variables 
     951      LOGICAL :: llog10            ! Perform log10 transform of variable 
     952#if defined key_fabm 
    559953      REAL(wp), POINTER, DIMENSION(:,:,:) :: & 
    560          & zprofmask1, &           ! Mask associated with zprofvar1 
    561          & zprofmask2              ! Mask associated with zprofvar2 
    562       REAL(wp), POINTER, DIMENSION(:,:) :: & 
    563          & zsurfvar                ! Model values equivalent to surface ob. 
    564       REAL(wp), POINTER, DIMENSION(:,:) :: & 
    565          & zglam1,    &            ! Model longitudes for prof variable 1 
    566          & zglam2,    &            ! Model longitudes for prof variable 2 
    567          & zgphi1,    &            ! Model latitudes for prof variable 1 
    568          & zgphi2                  ! Model latitudes for prof variable 2 
    569 #if ! defined key_lim2 && ! defined key_lim3 
    570       REAL(wp), POINTER, DIMENSION(:,:) :: frld 
    571 #endif 
    572       LOGICAL :: llnightav        ! Logical for calculating night-time average 
    573  
    574       !Allocate local work arrays 
    575       CALL wrk_alloc( jpi, jpj, jpk, zprofvar1 ) 
    576       CALL wrk_alloc( jpi, jpj, jpk, zprofvar2 ) 
    577       CALL wrk_alloc( jpi, jpj, jpk, zprofmask1 ) 
    578       CALL wrk_alloc( jpi, jpj, jpk, zprofmask2 ) 
    579       CALL wrk_alloc( jpi, jpj, zsurfvar ) 
    580       CALL wrk_alloc( jpi, jpj, zglam1 ) 
    581       CALL wrk_alloc( jpi, jpj, zglam2 ) 
    582       CALL wrk_alloc( jpi, jpj, zgphi1 ) 
    583       CALL wrk_alloc( jpi, jpj, zgphi2 ) 
    584 #if ! defined key_lim2 && ! defined key_lim3 
    585       CALL wrk_alloc(jpi,jpj,frld)  
     954         & pco2_3d                 ! 3D pCO2 from FABM 
    586955#endif 
    587956 
    … …  
    595964 
    596965      !----------------------------------------------------------------------- 
    597       ! No LIM => frld == 0.0_wp 
    598       !----------------------------------------------------------------------- 
    599 #if ! defined key_lim2 && ! defined key_lim3 
    600       frld(:,:) = 0.0_wp 
    601 #endif 
    602       !----------------------------------------------------------------------- 
    603966      ! Call the profile and surface observation operators 
    604967      !----------------------------------------------------------------------- 
    … …  
    608971         DO jtype = 1, nproftypes 
    609972 
     973            ! Allocate local work arrays 
     974            CALL wrk_alloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofvar  ) 
     975            CALL wrk_alloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofmask ) 
     976            CALL wrk_alloc( jpi, jpj,      profdataqc(jtype)%nvar, zglam     ) 
     977            CALL wrk_alloc( jpi, jpj,      profdataqc(jtype)%nvar, zgphi     ) 
     978             
     979            ! Defaults which might change 
     980            DO jvar = 1, profdataqc(jtype)%nvar 
     981               zprofmask(:,:,:,jvar) = tmask(:,:,:) 
     982               zglam(:,:,jvar)       = glamt(:,:) 
     983               zgphi(:,:,jvar)       = gphit(:,:) 
     984            END DO 
     985 
    610986            SELECT CASE ( TRIM(cobstypesprof(jtype)) ) 
    611987            CASE('prof') 
    612                zprofvar1(:,:,:) = tsn(:,:,:,jp_tem) 
    613                zprofvar2(:,:,:) = tsn(:,:,:,jp_sal) 
    614                zprofmask1(:,:,:) = tmask(:,:,:) 
    615                zprofmask2(:,:,:) = tmask(:,:,:) 
    616                zglam1(:,:) = glamt(:,:) 
    617                zglam2(:,:) = glamt(:,:) 
    618                zgphi1(:,:) = gphit(:,:) 
    619                zgphi2(:,:) = gphit(:,:) 
     988               zprofvar(:,:,:,1) = tsn(:,:,:,jp_tem) 
     989               zprofvar(:,:,:,2) = tsn(:,:,:,jp_sal) 
     990 
    620991            CASE('vel') 
    621                zprofvar1(:,:,:) = un(:,:,:) 
    622                zprofvar2(:,:,:) = vn(:,:,:) 
    623                zprofmask1(:,:,:) = umask(:,:,:) 
    624                zprofmask2(:,:,:) = vmask(:,:,:) 
    625                zglam1(:,:) = glamu(:,:) 
    626                zglam2(:,:) = glamv(:,:) 
    627                zgphi1(:,:) = gphiu(:,:) 
    628                zgphi2(:,:) = gphiv(:,:) 
     992               zprofvar(:,:,:,1) = un(:,:,:) 
     993               zprofvar(:,:,:,2) = vn(:,:,:) 
     994               zprofmask(:,:,:,1) = umask(:,:,:) 
     995               zprofmask(:,:,:,2) = vmask(:,:,:) 
     996               zglam(:,:,1) = glamu(:,:) 
     997               zglam(:,:,2) = glamv(:,:) 
     998               zgphi(:,:,1) = gphiu(:,:) 
     999               zgphi(:,:,2) = gphiv(:,:) 
     1000 
     1001            CASE('plchltot') 
     1002#if defined key_hadocc 
     1003               ! Chlorophyll from HadOCC 
     1004               zprofvar(:,:,:,1) = HADOCC_CHL(:,:,:) 
     1005#elif defined key_medusa 
     1006               ! Add non-diatom and diatom chlorophyll from MEDUSA 
     1007               zprofvar(:,:,:,1) = trn(:,:,:,jpchn) + trn(:,:,:,jpchd) 
     1008#elif defined key_fabm 
     1009               ! Add all chlorophyll groups from ERSEM 
     1010               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl2) + & 
     1011                  &                trn(:,:,:,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl4) 
     1012#else 
     1013               CALL ctl_stop( ' Trying to run plchltot observation operator', & 
     1014                  &           ' but no biogeochemical model appears to have been defined' ) 
     1015#endif 
     1016               ! Take the log10 where we can, otherwise exclude 
     1017               tiny = 1.0e-20 
     1018               WHERE(zprofvar(:,:,:,:) > tiny .AND. zprofvar(:,:,:,:) /= obfillflt ) 
     1019                  zprofvar(:,:,:,:)  = LOG10(zprofvar(:,:,:,:)) 
     1020               ELSEWHERE 
     1021                  zprofvar(:,:,:,:)  = obfillflt 
     1022                  zprofmask(:,:,:,:) = 0 
     1023               END WHERE 
     1024               ! Mask out model below any excluded values, 
     1025               ! to avoid interpolation issues 
     1026               DO jvar = 1, profdataqc(jtype)%nvar 
     1027                 DO jj = 1, jpj 
     1028                    DO ji = 1, jpi 
     1029                       depth_loop: DO jk = 1, jpk 
     1030                          IF ( zprofmask(ji,jj,jk,jvar) == 0 ) THEN 
     1031                             zprofmask(ji,jj,jk:jpk,jvar) = 0 
     1032                             EXIT depth_loop 
     1033                          ENDIF 
     1034                       END DO depth_loop 
     1035                    END DO 
     1036                 END DO 
     1037              END DO 
     1038 
     1039            CASE('pchltot') 
     1040#if defined key_hadocc 
     1041               ! Chlorophyll from HadOCC 
     1042               zprofvar(:,:,:,1) = HADOCC_CHL(:,:,:) 
     1043#elif defined key_medusa 
     1044               ! Add non-diatom and diatom chlorophyll from MEDUSA 
     1045               zprofvar(:,:,:,1) = trn(:,:,:,jpchn) + trn(:,:,:,jpchd) 
     1046#elif defined key_fabm 
     1047               ! Add all chlorophyll groups from ERSEM 
     1048               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl2) + & 
     1049                  &                trn(:,:,:,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl4) 
     1050#else 
     1051               CALL ctl_stop( ' Trying to run pchltot observation operator', & 
     1052                  &           ' but no biogeochemical model appears to have been defined' ) 
     1053#endif 
     1054 
     1055            CASE('pno3') 
     1056#if defined key_hadocc 
     1057               ! Dissolved inorganic nitrogen from HadOCC 
     1058               zprofvar(:,:,:,1) = trn(:,:,:,jp_had_nut) 
     1059#elif defined key_medusa 
     1060               ! Dissolved inorganic nitrogen from MEDUSA 
     1061               zprofvar(:,:,:,1) = trn(:,:,:,jpdin) 
     1062#elif defined key_fabm 
     1063               ! Nitrate from ERSEM 
     1064               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_n3n) 
     1065#else 
     1066               CALL ctl_stop( ' Trying to run pno3 observation operator', & 
     1067                  &           ' but no biogeochemical model appears to have been defined' ) 
     1068#endif 
     1069 
     1070            CASE('psi4') 
     1071#if defined key_hadocc 
     1072               CALL ctl_stop( ' Trying to run psi4 observation operator', & 
     1073                  &           ' but HadOCC does not simulate silicate' ) 
     1074#elif defined key_medusa 
     1075               ! Silicate from MEDUSA 
     1076               zprofvar(:,:,:,1) = trn(:,:,:,jpsil) 
     1077#elif defined key_fabm 
     1078               ! Silicate from ERSEM 
     1079               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_n5s) 
     1080#else 
     1081               CALL ctl_stop( ' Trying to run psi4 observation operator', & 
     1082                  &           ' but no biogeochemical model appears to have been defined' ) 
     1083#endif 
     1084 
     1085            CASE('ppo4') 
     1086#if defined key_hadocc 
     1087               CALL ctl_stop( ' Trying to run ppo4 observation operator', & 
     1088                  &           ' but HadOCC does not simulate phosphate' ) 
     1089#elif defined key_medusa 
     1090               CALL ctl_stop( ' Trying to run ppo4 observation operator', & 
     1091                  &           ' but MEDUSA does not simulate phosphate' ) 
     1092#elif defined key_fabm 
     1093               ! Phosphate from ERSEM 
     1094               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_n1p) 
     1095#else 
     1096               CALL ctl_stop( ' Trying to run ppo4 observation operator', & 
     1097                  &           ' but no biogeochemical model appears to have been defined' ) 
     1098#endif 
     1099 
     1100            CASE('pdic') 
     1101#if defined key_hadocc 
     1102               ! Dissolved inorganic carbon from HadOCC 
     1103               zprofvar(:,:,:,1) = trn(:,:,:,jp_had_dic) 
     1104#elif defined key_medusa 
     1105               ! Dissolved inorganic carbon from MEDUSA 
     1106               zprofvar(:,:,:,1) = trn(:,:,:,jpdic) 
     1107#elif defined key_fabm 
     1108               ! Dissolved inorganic carbon from ERSEM 
     1109               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_o3c) 
     1110#else 
     1111               CALL ctl_stop( ' Trying to run pdic observation operator', & 
     1112                  &           ' but no biogeochemical model appears to have been defined' ) 
     1113#endif 
     1114 
     1115            CASE('palk') 
     1116#if defined key_hadocc 
     1117               ! Alkalinity from HadOCC 
     1118               zprofvar(:,:,:,1) = trn(:,:,:,jp_had_alk) 
     1119#elif defined key_medusa 
     1120               ! Alkalinity from MEDUSA 
     1121               zprofvar(:,:,:,1) = trn(:,:,:,jpalk) 
     1122#elif defined key_fabm 
     1123               ! Alkalinity from ERSEM 
     1124               zprofvar(:,:,:,1) = fabm_get_interior_diagnostic_data(model, jp_fabm_o3ta) 
     1125#else 
     1126               CALL ctl_stop( ' Trying to run palk observation operator', & 
     1127                  &           ' but no biogeochemical model appears to have been defined' ) 
     1128#endif 
     1129 
     1130            CASE('pph') 
     1131#if defined key_hadocc 
     1132               CALL ctl_stop( ' Trying to run pph observation operator', & 
     1133                  &           ' but HadOCC has no pH diagnostic defined' ) 
     1134#elif defined key_medusa && defined key_roam 
     1135               ! pH from MEDUSA 
     1136               zprofvar(:,:,:,1) = f3_pH(:,:,:) 
     1137#elif defined key_fabm 
     1138               ! pH from ERSEM 
     1139               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_o3ph) 
     1140#else 
     1141               CALL ctl_stop( ' Trying to run pph observation operator', & 
     1142                  &           ' but no biogeochemical model appears to have been defined' ) 
     1143#endif 
     1144 
     1145            CASE('po2') 
     1146#if defined key_hadocc 
     1147               CALL ctl_stop( ' Trying to run po2 observation operator', & 
     1148                  &           ' but HadOCC does not simulate oxygen' ) 
     1149#elif defined key_medusa 
     1150               ! Oxygen from MEDUSA 
     1151               zprofvar(:,:,:,1) = trn(:,:,:,jpoxy) 
     1152#elif defined key_fabm 
     1153               ! Oxygen from ERSEM 
     1154               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_o2o) 
     1155#else 
     1156               CALL ctl_stop( ' Trying to run po2 observation operator', & 
     1157                  &           ' but no biogeochemical model appears to have been defined' ) 
     1158#endif 
     1159 
     1160            CASE DEFAULT 
     1161               CALL ctl_stop( 'Unknown profile observation type '//TRIM(cobstypesprof(jtype))//' in dia_obs' ) 
     1162 
    6291163            END SELECT 
    6301164 
    631             IF( ln_zco .OR. ln_zps ) THEN  
     1165            DO jvar = 1, profdataqc(jtype)%nvar 
    6321166               CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk,  & 
    633                   &               nit000, idaystp,                         & 
    634                   &               zprofvar1, zprofvar2,                    & 
    635                   &               gdept_1d, zprofmask1, zprofmask2,        & 
    636                   &               zglam1, zglam2, zgphi1, zgphi2,          & 
    637                   &               nn_1dint, nn_2dint,                      & 
     1167                  &               nit000, idaystp, jvar,                   & 
     1168                  &               zprofvar(:,:,:,jvar),                    & 
     1169                  &               gdept_n(:,:,:), gdepw_n(:,:,:),          &  
     1170                  &               zprofmask(:,:,:,jvar),                   & 
     1171                  &               zglam(:,:,jvar), zgphi(:,:,jvar),        & 
     1172                  &               nn_1dint, nn_2dint_default,              & 
    6381173                  &               kdailyavtypes = nn_profdavtypes ) 
    639             ELSE IF(TRIM(cobstypesprof(jtype)) == 'prof') THEN 
    640                !TG - THIS NEEDS MODIFICATION TO MATCH SIMPLIFICATION 
    641                CALL obs_pro_sco_opt( profdataqc(jtype),                    &  
    642                   &              kstp, jpi, jpj, jpk, nit000, idaystp,   &  
    643                   &              zprofvar1, zprofvar2,                   &  
    644                   &              gdept_n(:,:,:), gdepw_n(:,:,:),           & 
    645                   &              tmask, nn_1dint, nn_2dint,              &  
    646                   &              kdailyavtypes = nn_profdavtypes )  
    647             ELSE 
    648                CALL ctl_stop('DIA_OBS: Generalised vertical interpolation not'// & 
    649                          'yet working for velocity data (turn off velocity observations') 
    650             ENDIF 
     1174            END DO 
     1175 
     1176            CALL wrk_dealloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofvar  ) 
     1177            CALL wrk_dealloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofmask ) 
     1178            CALL wrk_dealloc( jpi, jpj,      profdataqc(jtype)%nvar, zglam     ) 
     1179            CALL wrk_dealloc( jpi, jpj,      profdataqc(jtype)%nvar, zgphi     ) 
    6511180 
    6521181         END DO 
    … …  
    6561185      IF ( nsurftypes > 0 ) THEN 
    6571186 
     1187         !Allocate local work arrays 
     1188         CALL wrk_alloc( jpi, jpj, zsurfvar ) 
     1189         CALL wrk_alloc( jpi, jpj, zsurfmask ) 
     1190#if defined key_fabm 
     1191         CALL wrk_alloc( jpi, jpj, jpk, pco2_3d ) 
     1192#endif 
     1193 
    6581194         DO jtype = 1, nsurftypes 
     1195 
     1196            !Defaults which might be changed 
     1197            zsurfmask(:,:) = tmask(:,:,1) 
     1198            llog10 = .FALSE. 
    6591199 
    6601200            SELECT CASE ( TRIM(cobstypessurf(jtype)) ) 
    6611201            CASE('sst') 
    6621202               zsurfvar(:,:) = tsn(:,:,1,jp_tem) 
    663                llnightav = ln_sstnight 
     1203               llnightav(jtype) = ln_sstnight 
    6641204            CASE('sla') 
    6651205               zsurfvar(:,:) = sshn(:,:) 
    666                llnightav = .FALSE. 
    667 #if defined key_lim2 || defined key_lim3 
     1206            CASE('sss') 
     1207               zsurfvar(:,:) = tsn(:,:,1,jp_sal) 
    6681208            CASE('sic') 
    6691209               IF ( kstp == 0 ) THEN 
    … …  
    6781218                  CYCLE 
    6791219               ELSE 
     1220#if defined key_cice 
     1221                  zsurfvar(:,:) = fr_i(:,:) 
     1222#elif defined key_lim2 || defined key_lim3 
    6801223                  zsurfvar(:,:) = 1._wp - frld(:,:) 
     1224#else 
     1225               CALL ctl_stop( ' Trying to run sea-ice observation operator', & 
     1226                  &           ' but no sea-ice model appears to have been defined' ) 
     1227#endif 
    6811228               ENDIF 
    682  
    683                llnightav = .FALSE. 
    684 #endif 
     1229            CASE('slchltot') 
     1230#if defined key_hadocc 
     1231               ! Surface chlorophyll from HadOCC 
     1232               zsurfvar(:,:) = HADOCC_CHL(:,:,1) 
     1233#elif defined key_medusa 
     1234               ! Add non-diatom and diatom surface chlorophyll from MEDUSA 
     1235               zsurfvar(:,:) = trn(:,:,1,jpchn) + trn(:,:,1,jpchd) 
     1236#elif defined key_fabm 
     1237               ! Add all surface chlorophyll groups from ERSEM 
     1238               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) + & 
     1239                  &            trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4) 
     1240#else 
     1241               CALL ctl_stop( ' Trying to run slchltot observation operator', & 
     1242                  &           ' but no biogeochemical model appears to have been defined' ) 
     1243#endif 
     1244               llog10 = .TRUE. 
     1245 
     1246            CASE('slchldia') 
     1247#if defined key_hadocc 
     1248               CALL ctl_stop( ' Trying to run slchldia observation operator', & 
     1249                  &           ' but HadOCC does not explicitly simulate diatoms' ) 
     1250#elif defined key_medusa 
     1251               ! Diatom surface chlorophyll from MEDUSA 
     1252               zsurfvar(:,:) = trn(:,:,1,jpchd) 
     1253#elif defined key_fabm 
     1254               ! Diatom surface chlorophyll from ERSEM 
     1255               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) 
     1256#else 
     1257               CALL ctl_stop( ' Trying to run slchldia observation operator', & 
     1258                  &           ' but no biogeochemical model appears to have been defined' ) 
     1259#endif 
     1260               llog10 = .TRUE. 
     1261 
     1262            CASE('slchlnon') 
     1263#if defined key_hadocc 
     1264               CALL ctl_stop( ' Trying to run slchlnon observation operator', & 
     1265                  &           ' but HadOCC does not explicitly simulate non-diatoms' ) 
     1266#elif defined key_medusa 
     1267               ! Non-diatom surface chlorophyll from MEDUSA 
     1268               zsurfvar(:,:) = trn(:,:,1,jpchn) 
     1269#elif defined key_fabm 
     1270               ! Add all non-diatom surface chlorophyll groups from ERSEM 
     1271               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) + & 
     1272                  &            trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4) 
     1273#else 
     1274               CALL ctl_stop( ' Trying to run slchlnon observation operator', & 
     1275                  &           ' but no biogeochemical model appears to have been defined' ) 
     1276#endif 
     1277               llog10 = .TRUE. 
     1278 
     1279            CASE('slchldin') 
     1280#if defined key_hadocc 
     1281               CALL ctl_stop( ' Trying to run slchldin observation operator', & 
     1282                  &           ' but HadOCC does not explicitly simulate dinoflagellates' ) 
     1283#elif defined key_medusa 
     1284               CALL ctl_stop( ' Trying to run slchldin observation operator', & 
     1285                  &           ' but MEDUSA does not explicitly simulate dinoflagellates' ) 
     1286#elif defined key_fabm 
     1287               ! Dinoflagellate surface chlorophyll from ERSEM 
     1288               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl4) 
     1289#else 
     1290               CALL ctl_stop( ' Trying to run slchldin observation operator', & 
     1291                  &           ' but no biogeochemical model appears to have been defined' ) 
     1292#endif 
     1293               llog10 = .TRUE. 
     1294 
     1295            CASE('slchlmic') 
     1296#if defined key_hadocc 
     1297               CALL ctl_stop( ' Trying to run slchlmic observation operator', & 
     1298                  &           ' but HadOCC does not explicitly simulate microphytoplankton' ) 
     1299#elif defined key_medusa 
     1300               CALL ctl_stop( ' Trying to run slchlmic observation operator', & 
     1301                  &           ' but MEDUSA does not explicitly simulate microphytoplankton' ) 
     1302#elif defined key_fabm 
     1303               ! Add diatom and dinoflagellate surface chlorophyll from ERSEM 
     1304               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4) 
     1305#else 
     1306               CALL ctl_stop( ' Trying to run slchlmic observation operator', & 
     1307                  &           ' but no biogeochemical model appears to have been defined' ) 
     1308#endif 
     1309               llog10 = .TRUE. 
     1310 
     1311            CASE('slchlnan') 
     1312#if defined key_hadocc 
     1313               CALL ctl_stop( ' Trying to run slchlnan observation operator', & 
     1314                  &           ' but HadOCC does not explicitly simulate nanophytoplankton' ) 
     1315#elif defined key_medusa 
     1316               CALL ctl_stop( ' Trying to run slchlnan observation operator', & 
     1317                  &           ' but MEDUSA does not explicitly simulate nanophytoplankton' ) 
     1318#elif defined key_fabm 
     1319               ! Nanophytoplankton surface chlorophyll from ERSEM 
     1320               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) 
     1321#else 
     1322               CALL ctl_stop( ' Trying to run slchlnan observation operator', & 
     1323                  &           ' but no biogeochemical model appears to have been defined' ) 
     1324#endif 
     1325               llog10 = .TRUE. 
     1326 
     1327            CASE('slchlpic') 
     1328#if defined key_hadocc 
     1329               CALL ctl_stop( ' Trying to run slchlpic observation operator', & 
     1330                  &           ' but HadOCC does not explicitly simulate picophytoplankton' ) 
     1331#elif defined key_medusa 
     1332               CALL ctl_stop( ' Trying to run slchlpic observation operator', & 
     1333                  &           ' but MEDUSA does not explicitly simulate picophytoplankton' ) 
     1334#elif defined key_fabm 
     1335               ! Picophytoplankton surface chlorophyll from ERSEM 
     1336               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) 
     1337#else 
     1338               CALL ctl_stop( ' Trying to run slchlpic observation operator', & 
     1339                  &           ' but no biogeochemical model appears to have been defined' ) 
     1340#endif 
     1341               llog10 = .TRUE. 
     1342 
     1343            CASE('schltot') 
     1344#if defined key_hadocc 
     1345               ! Surface chlorophyll from HadOCC 
     1346               zsurfvar(:,:) = HADOCC_CHL(:,:,1) 
     1347#elif defined key_medusa 
     1348               ! Add non-diatom and diatom surface chlorophyll from MEDUSA 
     1349               zsurfvar(:,:) = trn(:,:,1,jpchn) + trn(:,:,1,jpchd) 
     1350#elif defined key_fabm 
     1351               ! Add all surface chlorophyll groups from ERSEM 
     1352               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) + & 
     1353                  &            trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4) 
     1354#else 
     1355               CALL ctl_stop( ' Trying to run schltot observation operator', & 
     1356                  &           ' but no biogeochemical model appears to have been defined' ) 
     1357#endif 
     1358 
     1359            CASE('slphytot') 
     1360#if defined key_hadocc 
     1361               ! Surface phytoplankton nitrogen from HadOCC multiplied by C:N ratio 
     1362               zsurfvar(:,:) = trn(:,:,1,jp_had_phy) * c2n_p 
     1363#elif defined key_medusa 
     1364               ! Add non-diatom and diatom surface phytoplankton nitrogen from MEDUSA 
     1365               ! multiplied by C:N ratio for each 
     1366               zsurfvar(:,:) = (trn(:,:,1,jpphn) * xthetapn) + (trn(:,:,1,jpphd) * xthetapd) 
     1367#elif defined key_fabm 
     1368               ! Add all surface phytoplankton carbon groups from ERSEM 
     1369               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_p1c) + trn(:,:,1,jp_fabm_m1+jp_fabm_p2c) + & 
     1370                  &            trn(:,:,1,jp_fabm_m1+jp_fabm_p3c) + trn(:,:,1,jp_fabm_m1+jp_fabm_p4c) 
     1371#else 
     1372               CALL ctl_stop( ' Trying to run slphytot observation operator', & 
     1373                  &           ' but no biogeochemical model appears to have been defined' ) 
     1374#endif 
     1375               llog10 = .TRUE. 
     1376 
     1377            CASE('slphydia') 
     1378#if defined key_hadocc 
     1379               CALL ctl_stop( ' Trying to run slphydia observation operator', & 
     1380                  &           ' but HadOCC does not explicitly simulate diatoms' ) 
     1381#elif defined key_medusa 
     1382               ! Diatom surface phytoplankton nitrogen from MEDUSA multiplied by C:N ratio 
     1383               zsurfvar(:,:) = trn(:,:,1,jpphd) * xthetapd 
     1384#elif defined key_fabm 
     1385               ! Diatom surface phytoplankton carbon from ERSEM 
     1386               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_p1c) 
     1387#else 
     1388               CALL ctl_stop( ' Trying to run slphydia observation operator', & 
     1389                  &           ' but no biogeochemical model appears to have been defined' ) 
     1390#endif 
     1391               llog10 = .TRUE. 
     1392 
     1393            CASE('slphynon') 
     1394#if defined key_hadocc 
     1395               CALL ctl_stop( ' Trying to run slphynon observation operator', & 
     1396                  &           ' but HadOCC does not explicitly simulate non-diatoms' ) 
     1397#elif defined key_medusa 
     1398               ! Non-diatom surface phytoplankton nitrogen from MEDUSA multiplied by C:N ratio 
     1399               zsurfvar(:,:) = trn(:,:,1,jpphn) * xthetapn 
     1400#elif defined key_fabm 
     1401               ! Add all non-diatom surface phytoplankton carbon groups from ERSEM 
     1402               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_p2c) + & 
     1403                  &            trn(:,:,1,jp_fabm_m1+jp_fabm_p3c) + trn(:,:,1,jp_fabm_m1+jp_fabm_p4c) 
     1404#else 
     1405               CALL ctl_stop( ' Trying to run slphynon observation operator', & 
     1406                  &           ' but no biogeochemical model appears to have been defined' ) 
     1407#endif 
     1408               llog10 = .TRUE. 
     1409 
     1410            CASE('sspm') 
     1411#if defined key_spm 
     1412               zsurfvar(:,:) = 0.0 
     1413               DO jn = 1, jp_spm 
     1414                  zsurfvar(:,:) = zsurfvar(:,:) + trn(:,:,1,jn)   ! sum SPM sizes 
     1415               END DO 
     1416#else 
     1417               CALL ctl_stop( ' Trying to run sspm observation operator', & 
     1418                  &           ' but no spm model appears to have been defined' ) 
     1419#endif 
     1420 
     1421            CASE('sfco2') 
     1422#if defined key_hadocc 
     1423               zsurfvar(:,:) = HADOCC_FCO2(:,:)    ! fCO2 from HadOCC 
     1424               IF ( ( MINVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) .AND. & 
     1425                  & ( MAXVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) ) THEN 
     1426                  zsurfvar(:,:) = obfillflt 
     1427                  zsurfmask(:,:) = 0 
     1428                  CALL ctl_warn( ' HadOCC fCO2 values masked out for observation operator', & 
     1429                     &           ' as HADOCC_FCO2(:,:) == HADOCC_FILL_FLT' ) 
     1430               ENDIF 
     1431#elif defined key_medusa && defined key_roam 
     1432               zsurfvar(:,:) = f2_fco2w(:,:) 
     1433#elif defined key_fabm 
     1434               ! First, get pCO2 from FABM 
     1435               pco2_3d(:,:,:) = fabm_get_interior_diagnostic_data(model, jp_fabm_o3pc) 
     1436               zsurfvar(:,:) = pco2_3d(:,:,1) 
     1437               ! Now, convert pCO2 to fCO2, based on SST in K. This follows the standard methodology of: 
     1438               ! Pierrot et al. (2009), Recommendations for autonomous underway pCO2 measuring systems 
     1439               ! and data reduction routines, Deep-Sea Research II, 56: 512-522. 
     1440               ! and 
     1441               ! Weiss (1974), Carbon dioxide in water and seawater: the solubility of a non-ideal gas, 
     1442               ! Marine Chemistry, 2: 203-215. 
     1443               ! In the implementation below, atmospheric pressure has been assumed to be 1 atm and so 
     1444               ! not explicitly included - atmospheric pressure is not necessarily available so this is 
     1445               ! the best assumption. 
     1446               ! Further, the (1-xCO2)^2 term has been neglected. This is common practice 
     1447               ! (see e.g. Zeebe and Wolf-Gladrow (2001), CO2 in Seawater: Equilibrium, Kinetics, Isotopes) 
     1448               ! because xCO2 in atm is ~0, and so this term will only affect the result to the 3rd decimal 
     1449               ! place for typical values, and xCO2 would need to be approximated from pCO2 anyway. 
     1450               zsurfvar(:,:) = zsurfvar(:,:) * EXP((-1636.75                                                          + & 
     1451                  &            12.0408      * (tsn(:,:,1,jp_tem)+rt0)                                                 - & 
     1452                  &            0.0327957    * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0)                         + & 
     1453                  &            0.0000316528 * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0) + & 
     1454                  &            2.0 * (57.7 - 0.118 * (tsn(:,:,1,jp_tem)+rt0)))                                        / & 
     1455                  &            (82.0578 * (tsn(:,:,1,jp_tem)+rt0))) 
     1456#else 
     1457               CALL ctl_stop( ' Trying to run sfco2 observation operator', & 
     1458                  &           ' but no biogeochemical model appears to have been defined' ) 
     1459#endif 
     1460 
     1461            CASE('spco2') 
     1462#if defined key_hadocc 
     1463               zsurfvar(:,:) = HADOCC_PCO2(:,:)    ! pCO2 from HadOCC 
     1464               IF ( ( MINVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) .AND. & 
     1465                  & ( MAXVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) ) THEN 
     1466                  zsurfvar(:,:) = obfillflt 
     1467                  zsurfmask(:,:) = 0 
     1468                  CALL ctl_warn( ' HadOCC pCO2 values masked out for observation operator', & 
     1469                     &           ' as HADOCC_PCO2(:,:) == HADOCC_FILL_FLT' ) 
     1470               ENDIF 
     1471#elif defined key_medusa && defined key_roam 
     1472               zsurfvar(:,:) = f2_pco2w(:,:) 
     1473#elif defined key_fabm 
     1474               pco2_3d(:,:,:) = fabm_get_interior_diagnostic_data(model, jp_fabm_o3pc) 
     1475               zsurfvar(:,:) = pco2_3d(:,:,1) 
     1476#else 
     1477               CALL ctl_stop( ' Trying to run spco2 observation operator', & 
     1478                  &           ' but no biogeochemical model appears to have been defined' ) 
     1479#endif 
     1480 
     1481            CASE DEFAULT 
     1482 
     1483               CALL ctl_stop( 'Unknown surface observation type '//TRIM(cobstypessurf(jtype))//' in dia_obs' ) 
     1484 
    6851485            END SELECT 
     1486             
     1487            IF ( llog10 ) THEN 
     1488               ! Take the log10 where we can, otherwise exclude 
     1489               tiny = 1.0e-20 
     1490               WHERE(zsurfvar(:,:) > tiny .AND. zsurfvar(:,:) /= obfillflt ) 
     1491                  zsurfvar(:,:)  = LOG10(zsurfvar(:,:)) 
     1492               ELSEWHERE 
     1493                  zsurfvar(:,:)  = obfillflt 
     1494                  zsurfmask(:,:) = 0 
     1495               END WHERE 
     1496            ENDIF 
    6861497 
    6871498            CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj,       & 
    688                &               nit000, idaystp, zsurfvar, tmask(:,:,1), & 
    689                &               nn_2dint, llnightav ) 
     1499               &               nit000, idaystp, zsurfvar, zsurfmask,    & 
     1500               &               n2dintsurf(jtype), llnightav(jtype),     & 
     1501               &               ravglamscl(jtype), ravgphiscl(jtype),     & 
     1502               &               lfpindegs(jtype) ) 
    6901503 
    6911504         END DO 
    6921505 
    693       ENDIF 
    694  
    695       CALL wrk_dealloc( jpi, jpj, jpk, zprofvar1 ) 
    696       CALL wrk_dealloc( jpi, jpj, jpk, zprofvar2 ) 
    697       CALL wrk_dealloc( jpi, jpj, jpk, zprofmask1 ) 
    698       CALL wrk_dealloc( jpi, jpj, jpk, zprofmask2 ) 
    699       CALL wrk_dealloc( jpi, jpj, zsurfvar ) 
    700       CALL wrk_dealloc( jpi, jpj, zglam1 ) 
    701       CALL wrk_dealloc( jpi, jpj, zglam2 ) 
    702       CALL wrk_dealloc( jpi, jpj, zgphi1 ) 
    703       CALL wrk_dealloc( jpi, jpj, zgphi2 ) 
    704 #if ! defined key_lim2 && ! defined key_lim3 
    705       CALL wrk_dealloc(jpi,jpj,frld) 
    706 #endif 
     1506         CALL wrk_dealloc( jpi, jpj, zsurfvar ) 
     1507         CALL wrk_dealloc( jpi, jpj, zsurfmask ) 
     1508#if defined key_fabm 
     1509         CALL wrk_dealloc( jpi, jpj, jpk, pco2_3d ) 
     1510#endif 
     1511 
     1512      ENDIF 
    7071513 
    7081514   END SUBROUTINE dia_obs 
    … …  
    7551561                  & ) 
    7561562 
    757                CALL obs_rotvel( profdataqc(jtype), nn_2dint, zu, zv ) 
     1563               CALL obs_rotvel( profdataqc(jtype), nn_2dint_default, zu, zv ) 
    7581564 
    7591565               DO jo = 1, profdataqc(jtype)%nprof 
    … …  
    8211627 
    8221628      IF ( nsurftypes > 0 ) & 
    823          &   DEALLOCATE( cobstypessurf, surfdata, surfdataqc, nvarssurf, nextrsurf ) 
     1629         &   DEALLOCATE( cobstypessurf, surfdata, surfdataqc, nvarssurf, nextrsurf, & 
     1630         &               n2dintsurf, ravglamscl, ravgphiscl, lfpindegs, llnightav ) 
    8241631 
    8251632   END SUBROUTINE dia_obs_dealloc 
    … …  
    9701777   END SUBROUTINE fin_date 
    9711778    
     1779   SUBROUTINE obs_settypefiles( ntypes, jpmaxnfiles, ifiles, cobstypes, cfiles ) 
     1780 
     1781      INTEGER, INTENT(IN) :: ntypes      ! Total number of obs types 
     1782      INTEGER, INTENT(IN) :: jpmaxnfiles ! Maximum number of files allowed for each type 
     1783      INTEGER, DIMENSION(ntypes), INTENT(OUT) :: & 
     1784         &                   ifiles      ! Out number of files for each type 
     1785      CHARACTER(len=8), DIMENSION(ntypes), INTENT(IN) :: & 
     1786         &                   cobstypes   ! List of obs types 
     1787      CHARACTER(len=128), DIMENSION(ntypes, jpmaxnfiles), INTENT(IN) :: & 
     1788         &                   cfiles      ! List of files for all types 
     1789 
     1790      !Local variables 
     1791      INTEGER :: jfile 
     1792      INTEGER :: jtype 
     1793 
     1794      DO jtype = 1, ntypes 
     1795 
     1796         ifiles(jtype) = 0 
     1797         DO jfile = 1, jpmaxnfiles 
     1798            IF ( trim(cfiles(jtype,jfile)) /= '' ) & 
     1799                      ifiles(jtype) = ifiles(jtype) + 1 
     1800         END DO 
     1801 
     1802         IF ( ifiles(jtype) == 0 ) THEN 
     1803              CALL ctl_stop( 'Logical for observation type '//TRIM(cobstypes(jtype))//   & 
     1804                 &           ' set to true but no files available to read' ) 
     1805         ENDIF 
     1806 
     1807         IF(lwp) THEN     
     1808            WRITE(numout,*) '             '//cobstypes(jtype)//' input observation file names:' 
     1809            DO jfile = 1, ifiles(jtype) 
     1810               WRITE(numout,*) '                '//TRIM(cfiles(jtype,jfile)) 
     1811            END DO 
     1812         ENDIF 
     1813 
     1814      END DO 
     1815 
     1816   END SUBROUTINE obs_settypefiles 
     1817 
     1818   SUBROUTINE obs_setinterpopts( ntypes, jtype, ctypein,             & 
     1819              &                  n2dint_default, n2dint_type,        & 
     1820              &                  ravglamscl_type, ravgphiscl_type,   & 
     1821              &                  lfp_indegs_type, lavnight_type,     & 
     1822              &                  n2dint, ravglamscl, ravgphiscl,     & 
     1823              &                  lfpindegs, lavnight ) 
     1824 
     1825      INTEGER, INTENT(IN)  :: ntypes             ! Total number of obs types 
     1826      INTEGER, INTENT(IN)  :: jtype              ! Index of the current type of obs 
     1827      INTEGER, INTENT(IN)  :: n2dint_default     ! Default option for interpolation type 
     1828      INTEGER, INTENT(IN)  :: n2dint_type        ! Option for interpolation type 
     1829      REAL(wp), INTENT(IN) :: & 
     1830         &                    ravglamscl_type, & !E/W diameter of obs footprint for this type 
     1831         &                    ravgphiscl_type    !N/S diameter of obs footprint for this type 
     1832      LOGICAL, INTENT(IN)  :: lfp_indegs_type    !T=> footprint in degrees, F=> in metres 
     1833      LOGICAL, INTENT(IN)  :: lavnight_type      !T=> obs represent night time average 
     1834      CHARACTER(len=8), INTENT(IN) :: ctypein  
     1835 
     1836      INTEGER, DIMENSION(ntypes), INTENT(INOUT) :: & 
     1837         &                    n2dint  
     1838      REAL(wp), DIMENSION(ntypes), INTENT(INOUT) :: & 
     1839         &                    ravglamscl, ravgphiscl 
     1840      LOGICAL, DIMENSION(ntypes), INTENT(INOUT) :: & 
     1841         &                    lfpindegs, lavnight 
     1842 
     1843      lavnight(jtype) = lavnight_type 
     1844 
     1845      IF ( (n2dint_type >= 0) .AND. (n2dint_type <= 6) ) THEN 
     1846         n2dint(jtype) = n2dint_type 
     1847      ELSE IF ( n2dint_type == -1 ) THEN 
     1848         n2dint(jtype) = n2dint_default 
     1849      ELSE 
     1850         CALL ctl_stop(' Choice of '//TRIM(ctypein)//' horizontal (2D) interpolation method', & 
     1851           &                    ' is not available') 
     1852      ENDIF 
     1853 
     1854      ! For averaging observation footprints set options for size of footprint  
     1855      IF ( (n2dint(jtype) > 4) .AND. (n2dint(jtype) <= 6) ) THEN 
     1856         IF ( ravglamscl_type > 0._wp ) THEN 
     1857            ravglamscl(jtype) = ravglamscl_type 
     1858         ELSE 
     1859            CALL ctl_stop( 'Incorrect value set for averaging footprint '// & 
     1860                           'scale (ravglamscl) for observation type '//TRIM(ctypein) )       
     1861         ENDIF 
     1862 
     1863         IF ( ravgphiscl_type > 0._wp ) THEN 
     1864            ravgphiscl(jtype) = ravgphiscl_type 
     1865         ELSE 
     1866            CALL ctl_stop( 'Incorrect value set for averaging footprint '// & 
     1867                           'scale (ravgphiscl) for observation type '//TRIM(ctypein) )       
     1868         ENDIF 
     1869 
     1870         lfpindegs(jtype) = lfp_indegs_type  
     1871 
     1872      ENDIF 
     1873 
     1874      ! Write out info  
     1875      IF(lwp) THEN 
     1876         IF ( n2dint(jtype) <= 4 ) THEN 
     1877            WRITE(numout,*) '             '//TRIM(ctypein)// & 
     1878               &            ' model counterparts will be interpolated horizontally' 
     1879         ELSE IF ( n2dint(jtype) <= 6 ) THEN 
     1880            WRITE(numout,*) '             '//TRIM(ctypein)// & 
     1881               &            ' model counterparts will be averaged horizontally' 
     1882            WRITE(numout,*) '             '//'    with E/W scale: ',ravglamscl(jtype) 
     1883            WRITE(numout,*) '             '//'    with N/S scale: ',ravgphiscl(jtype) 
     1884            IF ( lfpindegs(jtype) ) THEN 
     1885                WRITE(numout,*) '             '//'    (in degrees)' 
     1886            ELSE 
     1887                WRITE(numout,*) '             '//'    (in metres)' 
     1888            ENDIF 
     1889         ENDIF 
     1890      ENDIF 
     1891 
     1892   END SUBROUTINE obs_setinterpopts 
     1893 
    9721894END MODULE diaobs 

  • NEMO/branches/UKMO/r8395_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_mpp.F90

    r11350 r11361  
    9898      ! 
    9999      INTEGER :: ierr  
    100       INTEGER, DIMENSION(kno) ::   ivals 
    101       ! 
    102 INCLUDE 'mpif.h' 
    103       !!---------------------------------------------------------------------- 
     100      INTEGER, DIMENSION(:), ALLOCATABLE ::   ivals 
     101      ! 
     102INCLUDE 'mpif.h' 
     103      !!---------------------------------------------------------------------- 
     104 
     105      ALLOCATE( ivals(kno) ) 
    104106 
    105107      ! Call the MPI library to find the maximum across processors 
    … …  
    107109         &                mpi_max, mpi_comm_opa, ierr ) 
    108110      kvals(:) = ivals(:) 
     111 
     112      DEALLOCATE( ivals ) 
    109113#else 
    110114      ! no MPI: empty routine 
    … …  
    138142      ! 
    139143      INTEGER :: ji, isum 
    140       INTEGER, DIMENSION(kno) ::   iobsp 
    141       !! 
    142       !! 
    143  
    144       iobsp=kobsp 
     144      INTEGER, DIMENSION(:), ALLOCATABLE ::   iobsp 
     145      !! 
     146      !! 
     147 
     148      ALLOCATE( iobsp(kno) ) 
     149 
     150      iobsp(:)=kobsp(:) 
    145151 
    146152      WHERE( iobsp(:) == -1 ) 
    … …  
    148154      END WHERE 
    149155 
    150       iobsp=-1*iobsp 
     156      iobsp(:)=-1*iobsp(:) 
    151157 
    152158      CALL obs_mpp_max_integer( iobsp, kno ) 
    153159 
    154       kobsp=-1*iobsp 
     160      kobsp(:)=-1*iobsp(:) 
    155161 
    156162      isum=0 
    … …  
    168174      ENDIF 
    169175 
     176      DEALLOCATE( iobsp ) 
     177 
    170178#else 
    171179      ! no MPI: empty routine 

  • NEMO/branches/UKMO/r8395_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90

    r11350 r11361  
    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 
    60    SUBROUTINE obs_prof_opt( prodatqc, kt, kpi, kpj, kpk,          & 
    61       &                     kit000, kdaystp,                      & 
    62       &                     pvar1, pvar2, pgdept, pmask1, pmask2, & 
    63       &                     plam1, plam2, pphi1, pphi2,           & 
     59 
     60   SUBROUTINE obs_prof_opt( prodatqc, kt, kpi, kpj, kpk, & 
     61      &                     kit000, kdaystp, kvar,       & 
     62      &                     pvar, pgdept, pgdepw,        & 
     63      &                     pmask,                       &   
     64      &                     plam, pphi,                  & 
    6465      &                     k1dint, k2dint, kdailyavtypes ) 
    6566 
    … …  
    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 
    … …  
    130132      INTEGER, INTENT(IN) :: k2dint  ! Horizontal interpolation type (see header) 
    131133      INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day 
     134      INTEGER, INTENT(IN) :: kvar    ! Number of variable in prodatqc 
    132135      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 
    133          & pvar1,    &               ! Model field 1 
    134          & pvar2,    &               ! Model field 2 
    135          & pmask1,   &               ! Land-sea mask 1 
    136          & pmask2                    ! Land-sea mask 2 
     136         & pvar,   &                 ! Model field for variable 
     137         & pmask                     ! Land-sea mask for variable 
    137138      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 
    138          & plam1,    &               ! Model longitudes for variable 1 
    139          & plam2,    &               ! Model longitudes for variable 2 
    140          & pphi1,    &               ! Model latitudes for variable 1 
    141          & pphi2                     ! Model latitudes for variable 2 
    142       REAL(KIND=wp), INTENT(IN), DIMENSION(kpk) :: & 
    143          & pgdept                    ! Model array of depth levels 
     139         & plam,   &                 ! Model longitudes for variable 
     140         & pphi                      ! Model latitudes for variable 
     141      REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: &  
     142         & pgdept, &                 ! Model array of depth T levels  
     143         & pgdepw                    ! Model array of depth W levels  
    144144      INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 
    145145         & kdailyavtypes             ! Types for daily averages 
    … …  
    156156      INTEGER ::   iend 
    157157      INTEGER ::   iobs 
     158      INTEGER ::   iin, ijn, ikn, ik   ! looping indices over interpolation nodes  
     159      INTEGER ::   inum_obs 
    158160      INTEGER, DIMENSION(imaxavtypes) :: & 
    159161         & idailyavtypes 
    160162      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 
    161          & igrdi1, & 
    162          & igrdi2, & 
    163          & igrdj1, & 
    164          & igrdj2 
     163         & igrdi, & 
     164         & igrdj 
     165      INTEGER, ALLOCATABLE, DIMENSION(:) :: iv_indic 
     166 
    165167      REAL(KIND=wp) :: zlam 
    166168      REAL(KIND=wp) :: zphi 
    167169      REAL(KIND=wp) :: zdaystp 
    168170      REAL(KIND=wp), DIMENSION(kpk) :: & 
    169          & zobsmask1, & 
    170          & zobsmask2, & 
    171171         & zobsk,    & 
    172172         & zobs2k 
    173       REAL(KIND=wp), DIMENSION(2,2,kpk) :: & 
     173      REAL(KIND=wp), DIMENSION(2,2,1) :: & 
    174174         & zweig1, & 
    175          & zweig2 
     175         & zweig 
    176176      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: & 
    177          & zmask1, & 
    178          & zmask2, & 
    179          & zint1, & 
    180          & zint2, & 
    181          & zinm1, & 
    182          & zinm2 
     177         & zmask,  & 
     178         & zint,   & 
     179         & zinm,   & 
     180         & zgdept, &  
     181         & zgdepw 
    183182      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 
    184          & zglam1, & 
    185          & zglam2, & 
    186          & zgphi1, & 
    187          & zgphi2 
     183         & zglam,  & 
     184         & zgphi 
     185      REAL(KIND=wp), DIMENSION(1) :: zmsk 
     186      REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: interp_corner 
     187 
    188188      LOGICAL :: ld_dailyav 
    189189 
    … …  
    215215               DO jj = 1, jpj 
    216216                  DO ji = 1, jpi 
    217                      prodatqc%vdmean(ji,jj,jk,1) = 0.0 
    218                      prodatqc%vdmean(ji,jj,jk,2) = 0.0 
     217                     prodatqc%vdmean(ji,jj,jk,kvar) = 0.0 
    219218                  END DO 
    220219               END DO 
    … …  
    225224            DO jj = 1, jpj 
    226225               DO ji = 1, jpi 
    227                   ! Increment field 1 for computing daily mean 
    228                   prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) & 
    229                      &                        + pvar1(ji,jj,jk) 
    230                   ! Increment field 2 for computing daily mean 
    231                   prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 
    232                      &                        + pvar2(ji,jj,jk) 
     226                  ! Increment field for computing daily mean 
     227                  prodatqc%vdmean(ji,jj,jk,kvar) = prodatqc%vdmean(ji,jj,jk,kvar) & 
     228                     &                           + pvar(ji,jj,jk) 
    233229               END DO 
    234230            END DO 
    … …  
    243239               DO jj = 1, jpj 
    244240                  DO ji = 1, jpi 
    245                      prodatqc%vdmean(ji,jj,jk,1) = prodatqc%vdmean(ji,jj,jk,1) & 
    246                         &                        * zdaystp 
    247                      prodatqc%vdmean(ji,jj,jk,2) = prodatqc%vdmean(ji,jj,jk,2) & 
    248                         &                        * zdaystp 
     241                     prodatqc%vdmean(ji,jj,jk,kvar) = prodatqc%vdmean(ji,jj,jk,kvar) & 
     242                        &                           * zdaystp 
    249243                  END DO 
    250244               END DO 
    … …  
    256250      ! Get the data for interpolation 
    257251      ALLOCATE( & 
    258          & igrdi1(2,2,ipro),      & 
    259          & igrdi2(2,2,ipro),      & 
    260          & igrdj1(2,2,ipro),      & 
    261          & igrdj2(2,2,ipro),      & 
    262          & zglam1(2,2,ipro),      & 
    263          & zglam2(2,2,ipro),      & 
    264          & zgphi1(2,2,ipro),      & 
    265          & zgphi2(2,2,ipro),      & 
    266          & zmask1(2,2,kpk,ipro),  & 
    267          & zmask2(2,2,kpk,ipro),  & 
    268          & zint1(2,2,kpk,ipro),  & 
    269          & zint2(2,2,kpk,ipro)   & 
     252         & igrdi(2,2,ipro),       & 
     253         & igrdj(2,2,ipro),       & 
     254         & zglam(2,2,ipro),       & 
     255         & zgphi(2,2,ipro),       & 
     256         & zmask(2,2,kpk,ipro),   & 
     257         & zint(2,2,kpk,ipro),    & 
     258         & zgdept(2,2,kpk,ipro),  &  
     259         & zgdepw(2,2,kpk,ipro)   &  
    270260         & ) 
    271261 
    272262      DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro 
    273263         iobs = jobs - prodatqc%nprofup 
    274          igrdi1(1,1,iobs) = prodatqc%mi(jobs,1)-1 
    275          igrdj1(1,1,iobs) = prodatqc%mj(jobs,1)-1 
    276          igrdi1(1,2,iobs) = prodatqc%mi(jobs,1)-1 
    277          igrdj1(1,2,iobs) = prodatqc%mj(jobs,1) 
    278          igrdi1(2,1,iobs) = prodatqc%mi(jobs,1) 
    279          igrdj1(2,1,iobs) = prodatqc%mj(jobs,1)-1 
    280          igrdi1(2,2,iobs) = prodatqc%mi(jobs,1) 
    281          igrdj1(2,2,iobs) = prodatqc%mj(jobs,1) 
    282          igrdi2(1,1,iobs) = prodatqc%mi(jobs,2)-1 
    283          igrdj2(1,1,iobs) = prodatqc%mj(jobs,2)-1 
    284          igrdi2(1,2,iobs) = prodatqc%mi(jobs,2)-1 
    285          igrdj2(1,2,iobs) = prodatqc%mj(jobs,2) 
    286          igrdi2(2,1,iobs) = prodatqc%mi(jobs,2) 
    287          igrdj2(2,1,iobs) = prodatqc%mj(jobs,2)-1 
    288          igrdi2(2,2,iobs) = prodatqc%mi(jobs,2) 
    289          igrdj2(2,2,iobs) = prodatqc%mj(jobs,2) 
     264         igrdi(1,1,iobs) = prodatqc%mi(jobs,kvar)-1 
     265         igrdj(1,1,iobs) = prodatqc%mj(jobs,kvar)-1 
     266         igrdi(1,2,iobs) = prodatqc%mi(jobs,kvar)-1 
     267         igrdj(1,2,iobs) = prodatqc%mj(jobs,kvar) 
     268         igrdi(2,1,iobs) = prodatqc%mi(jobs,kvar) 
     269         igrdj(2,1,iobs) = prodatqc%mj(jobs,kvar)-1 
     270         igrdi(2,2,iobs) = prodatqc%mi(jobs,kvar) 
     271         igrdj(2,2,iobs) = prodatqc%mj(jobs,kvar) 
    290272      END DO 
    291273 
    292       CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi1, igrdj1, plam1, zglam1 ) 
    293       CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi1, igrdj1, pphi1, zgphi1 ) 
    294       CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pmask1, zmask1 ) 
    295       CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, pvar1,   zint1 ) 
    296        
    297       CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi2, igrdj2, plam2, zglam2 ) 
    298       CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi2, igrdj2, pphi2, zgphi2 ) 
    299       CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, pmask2, zmask2 ) 
    300       CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, pvar2,   zint2 ) 
     274      ! Initialise depth arrays 
     275      zgdept(:,:,:,:) = 0.0 
     276      zgdepw(:,:,:,:) = 0.0 
     277 
     278      CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi, igrdj, plam, zglam ) 
     279      CALL obs_int_comm_2d( 2, 2, ipro, kpi, kpj, igrdi, igrdj, pphi, zgphi ) 
     280      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pmask, zmask ) 
     281      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pvar,   zint ) 
     282 
     283      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pgdept, zgdept )  
     284      CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, pgdepw, zgdepw )  
    301285 
    302286      ! At the end of the day also get interpolated means 
    303287      IF ( ld_dailyav .AND. idayend == 0 ) THEN 
    304288 
    305          ALLOCATE( & 
    306             & zinm1(2,2,kpk,ipro),  & 
    307             & zinm2(2,2,kpk,ipro)   & 
    308             & ) 
    309  
    310          CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi1, igrdj1, & 
    311             &                  prodatqc%vdmean(:,:,:,1), zinm1 ) 
    312          CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi2, igrdj2, & 
    313             &                  prodatqc%vdmean(:,:,:,2), zinm2 ) 
     289         ALLOCATE( zinm(2,2,kpk,ipro) ) 
     290 
     291         CALL obs_int_comm_3d( 2, 2, ipro, kpi, kpj, kpk, igrdi, igrdj, & 
     292            &                  prodatqc%vdmean(:,:,:,kvar), zinm ) 
    314293 
    315294      ENDIF 
     295 
     296      ! Return if no observations to process  
     297      ! Has to be done after comm commands to ensure processors  
     298      ! stay in sync  
     299      IF ( ipro == 0 ) RETURN  
    316300 
    317301      DO jobs = prodatqc%nprofup + 1, prodatqc%nprofup + ipro 
    … …  
    339323         zphi = prodatqc%rphi(jobs) 
    340324 
    341          ! Horizontal weights and vertical mask 
    342  
    343          IF ( prodatqc%npvend(jobs,1) > 0 ) THEN 
    344  
    345             CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi,     & 
    346                &                   zglam1(:,:,iobs), zgphi1(:,:,iobs), & 
    347                &                   zmask1(:,:,:,iobs), zweig1, zobsmask1 ) 
    348  
    349          ENDIF 
    350  
    351          IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 
    352  
    353             CALL obs_int_h2d_init( kpk, kpk, k2dint, zlam, zphi,     & 
    354                &                   zglam2(:,:,iobs), zgphi2(:,:,iobs), & 
    355                &                   zmask2(:,:,:,iobs), zweig2, zobsmask2 ) 
    356   
    357          ENDIF 
    358  
    359          IF ( prodatqc%npvend(jobs,1) > 0 ) THEN 
     325         ! Horizontal weights  
     326         ! Masked values are calculated later.   
     327         IF ( prodatqc%npvend(jobs,kvar) > 0 ) THEN 
     328 
     329            CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,     & 
     330               &                   zglam(:,:,iobs), zgphi(:,:,iobs), & 
     331               &                   zmask(:,:,1,iobs), zweig1, zmsk ) 
     332 
     333         ENDIF 
     334 
     335         IF ( prodatqc%npvend(jobs,kvar) > 0 ) THEN 
    360336 
    361337            zobsk(:) = obfillflt 
    … …  
    365341               IF ( idayend == 0 )  THEN 
    366342                  ! 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  
    405          IF ( prodatqc%npvend(jobs,2) > 0 ) THEN 
    406  
    407             zobsk(:) = obfillflt 
    408  
    409             IF ( ANY (idailyavtypes(:) == prodatqc%ntyp(jobs)) ) THEN 
    410  
    411                IF ( idayend == 0 )  THEN 
    412  
    413                   ! 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 
    454  
    455       ! Deallocate the data for interpolation 
    456       DEALLOCATE( & 
    457          & igrdi1, & 
    458          & igrdi2, & 
    459          & igrdj1, & 
    460          & igrdj2, & 
    461          & zglam1, & 
    462          & zglam2, & 
    463          & zgphi1, & 
    464          & zgphi2, & 
    465          & zmask1, & 
    466          & zmask2, & 
    467          & zint1,  & 
    468          & zint2   & 
    469          & ) 
    470  
    471       ! At the end of the day also get interpolated means 
    472       IF ( ld_dailyav .AND. idayend == 0 ) THEN 
    473          DEALLOCATE( & 
    474             & zinm1,  & 
    475             & zinm2   & 
    476             & ) 
    477       ENDIF 
    478  
    479       prodatqc%nprofup = prodatqc%nprofup + ipro  
    480  
    481    END SUBROUTINE obs_prof_opt 
    482  
    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       
     343 
    831344                  ! vertically interpolate all 4 corners  
    832                   ista = prodatqc%npvsta(jobs,1)  
    833                   iend = prodatqc%npvend(jobs,1)  
     345                  ista = prodatqc%npvsta(jobs,kvar)  
     346                  iend = prodatqc%npvend(jobs,kvar)  
    834347                  inum_obs = iend - ista + 1  
    835                   ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs))  
    836                   DO iin=1,2   
     348                  ALLOCATE(interp_corner(2,2,inum_obs),iv_indic(inum_obs))  
     349 
     350                  DO iin=1,2  
    837351                     DO ijn=1,2  
    838                                      
    839                                      
     352 
    840353                        IF ( k1dint == 1 ) THEN  
    841354                           CALL obs_int_z1d_spl( kpk, &  
    842                               &    zintt(iin,ijn,:,iobs),&  
    843                               &    zobs2k, zgdept(iin,ijn,:,iobs), &  
    844                               &    zmask(iin,ijn,:,iobs))  
    845    
     355                              &     zinm(iin,ijn,:,iobs), &  
     356                              &     zobs2k, zgdept(iin,ijn,:,iobs), &  
     357                              &     zmask(iin,ijn,:,iobs))  
    846358                        ENDIF  
    847359        
    848360                        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           
     361                           &    zgdept(iin,ijn,:,iobs), &  
     362                           &    inum_obs, prodatqc%var(kvar)%vdep(ista:iend), &  
     363                           &    iv_indic)  
     364 
     365                        CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs, &  
     366                           &    prodatqc%var(kvar)%vdep(ista:iend), &  
     367                           &    zinm(iin,ijn,:,iobs), &  
     368                           &    zobs2k, interp_corner(iin,ijn,:), &  
     369                           &    zgdept(iin,ijn,:,iobs), &  
     370                           &    zmask(iin,ijn,:,iobs))  
     371        
    859372                     ENDDO  
    860373                  ENDDO  
     374 
     375               ENDIF !idayend 
     376 
     377            ELSE    
     378 
     379               ! Point data  
     380      
     381               ! vertically interpolate all 4 corners  
     382               ista = prodatqc%npvsta(jobs,kvar)  
     383               iend = prodatqc%npvend(jobs,kvar)  
     384               inum_obs = iend - ista + 1  
     385               ALLOCATE(interp_corner(2,2,inum_obs), iv_indic(inum_obs))  
     386               DO iin=1,2   
     387                  DO ijn=1,2  
     388                     
     389                     IF ( k1dint == 1 ) THEN  
     390                        CALL obs_int_z1d_spl( kpk, &  
     391                           &    zint(iin,ijn,:,iobs),&  
     392                           &    zobs2k, zgdept(iin,ijn,:,iobs), &  
     393                           &    zmask(iin,ijn,:,iobs))  
     394   
     395                     ENDIF  
     396        
     397                     CALL obs_level_search(kpk, &  
     398                         &        zgdept(iin,ijn,:,iobs),&  
     399                         &        inum_obs, prodatqc%var(kvar)%vdep(ista:iend), &  
     400                         &        iv_indic)  
     401 
     402                     CALL obs_int_z1d(kpk, iv_indic, k1dint, inum_obs,     &  
     403                         &          prodatqc%var(kvar)%vdep(ista:iend),     &  
     404                         &          zint(iin,ijn,:,iobs),            &  
     405                         &          zobs2k,interp_corner(iin,ijn,:), &  
     406                         &          zgdept(iin,ijn,:,iobs),         &  
     407                         &          zmask(iin,ijn,:,iobs) )       
     408          
     409                  ENDDO  
     410               ENDDO  
    861411              
    862                ENDIF  
    863         
    864                !-------------------------------------------------------------  
    865                ! Compute the horizontal interpolation for every profile level  
    866                !-------------------------------------------------------------  
     412            ENDIF  
     413 
     414            !-------------------------------------------------------------  
     415            ! Compute the horizontal interpolation for every profile level  
     416            !-------------------------------------------------------------  
    867417              
    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  
     418            DO ikn=1,inum_obs  
     419               iend=ista+ikn-1 
     420                   
     421               zweig(:,:,1) = 0._wp  
     422    
     423               ! This code forces the horizontal weights to be   
     424               ! zero IF the observation is below the bottom of the   
     425               ! corners of the interpolation nodes, Or if it is in   
     426               ! the mask. This is important for observations near   
     427               ! steep bathymetry  
     428               DO iin=1,2  
     429                  DO ijn=1,2  
    880430      
    881                         depth_loop1: DO ik=kpk,2,-1  
    882                            IF(zmask(iin,ijn,ik-1,iobs ) > 0.9 )THEN    
     431                     depth_loop: DO ik=kpk,2,-1  
     432                        IF(zmask(iin,ijn,ik-1,iobs ) > 0.9 )THEN    
    883433                             
    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)  
     434                           zweig(iin,ijn,1) = &   
     435                              & zweig1(iin,ijn,1) * &  
     436                              & MAX( SIGN(1._wp,(zgdepw(iin,ijn,ik,iobs) ) &  
     437                              &  - prodatqc%var(kvar)%vdep(iend)),0._wp)  
    888438                             
    889                               EXIT depth_loop1  
    890                            ENDIF  
    891                         ENDDO depth_loop1  
     439                           EXIT depth_loop 
     440 
     441                        ENDIF  
     442 
     443                     ENDDO depth_loop 
    892444      
    893                      ENDDO  
    894445                  ENDDO  
     446               ENDDO  
    895447    
    896                   CALL obs_int_h2d( 1, 1, l_zweig, interp_corner(:,:,ikn), &  
    897                   &          prodatqc%var(1)%vmod(iend:iend) )  
     448               CALL obs_int_h2d( 1, 1, zweig, interp_corner(:,:,ikn), &  
     449                  &              prodatqc%var(kvar)%vmod(iend:iend) )  
     450 
     451                  ! Set QC flag for any observations found below the bottom 
     452                  ! needed as the check here is more strict than that in obs_prep 
     453               IF (sum(zweig) == 0.0_wp) prodatqc%var(kvar)%nvqc(iend:iend)=4 
    898454  
    899                ENDDO  
     455            ENDDO  
    900456  
    901   
    902                DEALLOCATE(interp_corner,iv_indic)  
     457            DEALLOCATE(interp_corner,iv_indic)  
    903458           
    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 ) 
     459         ENDIF 
     460 
     461      ENDDO 
     462 
     463      ! Deallocate the data for interpolation 
     464      DEALLOCATE(  & 
     465         & igrdi,  & 
     466         & igrdj,  & 
     467         & zglam,  & 
     468         & zgphi,  & 
     469         & zmask,  & 
     470         & zint,   & 
     471         & zgdept, & 
     472         & zgdepw  & 
     473         & ) 
     474 
     475      ! At the end of the day also get interpolated means 
     476      IF ( ld_dailyav .AND. idayend == 0 ) THEN 
     477         DEALLOCATE( zinm ) 
     478      ENDIF 
     479 
     480      IF ( kvar == prodatqc%nvar ) THEN 
     481         prodatqc%nprofup = prodatqc%nprofup + ipro  
     482      ENDIF 
     483 
     484   END SUBROUTINE obs_prof_opt 
     485 
     486   SUBROUTINE obs_surf_opt( surfdataqc, kt, kpi, kpj,            & 
     487      &                     kit000, kdaystp, psurf, psurfmask,   & 
     488      &                     k2dint, ldnightav, plamscl, pphiscl, & 
     489      &                     lindegrees ) 
    1071490 
    1072491      !!----------------------------------------------------------------------- 
    … …  
    1090509      !!        - polynomial (quadrilateral grid)  (k2dint = 4) 
    1091510      !! 
     511      !!    Two horizontal averaging schemes are also available: 
     512      !!        - weighted radial footprint        (k2dint = 5) 
     513      !!        - weighted rectangular footprint   (k2dint = 6) 
     514      !! 
    1092515      !! 
    1093516      !! ** Action  : 
    … …  
    1096519      !!      ! 07-03 (A. Weaver) 
    1097520      !!      ! 15-02 (M. Martin) Combined routine for surface types 
     521      !!      ! 17-03 (M. Martin) Added horizontal averaging options 
    1098522      !!----------------------------------------------------------------------- 
    1099523 
    … …  
    1117541         & psurfmask                   ! Land-sea mask 
    1118542      LOGICAL, INTENT(IN) :: ldnightav ! Logical for averaging night-time data 
     543      REAL(KIND=wp), INTENT(IN) :: & 
     544         & plamscl, &                  ! Diameter in metres of obs footprint in E/W, N/S directions 
     545         & pphiscl                     ! This is the full width (rather than half-width) 
     546      LOGICAL, INTENT(IN) :: & 
     547         & lindegrees                  ! T=> plamscl and pphiscl are specified in degrees, F=> in metres 
    1119548 
    1120549      !! * Local declarations 
    … …  
    1125554      INTEGER :: isurf 
    1126555      INTEGER :: iobs 
     556      INTEGER :: imaxifp, imaxjfp 
     557      INTEGER :: imodi, imodj 
    1127558      INTEGER :: idayend 
    1128559      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: & 
    1129          & igrdi, & 
    1130          & igrdj 
     560         & igrdi,   & 
     561         & igrdj,   & 
     562         & igrdip1, & 
     563         & igrdjp1 
    1131564      INTEGER, DIMENSION(:,:), SAVE, ALLOCATABLE :: & 
    1132565         & icount_night,      & 
    … …  
    1136569      REAL(wp), DIMENSION(1) :: zext, zobsmask 
    1137570      REAL(wp) :: zdaystp 
    1138       REAL(wp), DIMENSION(2,2,1) :: & 
    1139          & zweig 
    1140571      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 
     572         & zweig,  & 
    1141573         & zmask,  & 
    1142574         & zsurf,  & 
    1143575         & zsurfm, & 
     576         & zsurftmp, & 
    1144577         & zglam,  & 
    1145          & zgphi 
     578         & zgphi,  & 
     579         & zglamf, & 
     580         & zgphif 
     581 
    1146582      REAL(wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: & 
    1147583         & zintmp,  & 
    … …  
    1155591      inrc = kt - kit000 + 2 
    1156592      isurf = surfdataqc%nsstp(inrc) 
     593 
     594      ! Work out the maximum footprint size for the  
     595      ! interpolation/averaging in model grid-points - has to be even. 
     596 
     597      CALL obs_max_fpsize( k2dint, plamscl, pphiscl, lindegrees, psurfmask, imaxifp, imaxjfp ) 
     598 
    1157599 
    1158600      IF ( ldnightav ) THEN 
    … …  
    1221663 
    1222664      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)  & 
     665         & zweig(imaxifp,imaxjfp,1),      & 
     666         & igrdi(imaxifp,imaxjfp,isurf), & 
     667         & igrdj(imaxifp,imaxjfp,isurf), & 
     668         & zglam(imaxifp,imaxjfp,isurf), & 
     669         & zgphi(imaxifp,imaxjfp,isurf), & 
     670         & zmask(imaxifp,imaxjfp,isurf), & 
     671         & zsurf(imaxifp,imaxjfp,isurf), & 
     672         & zsurftmp(imaxifp,imaxjfp,isurf),  & 
     673         & zglamf(imaxifp+1,imaxjfp+1,isurf), & 
     674         & zgphif(imaxifp+1,imaxjfp+1,isurf), & 
     675         & igrdip1(imaxifp+1,imaxjfp+1,isurf), & 
     676         & igrdjp1(imaxifp+1,imaxjfp+1,isurf) & 
    1229677         & ) 
    1230678 
    1231679      DO jobs = surfdataqc%nsurfup + 1, surfdataqc%nsurfup + isurf 
    1232680         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) 
     681         DO ji = 0, imaxifp 
     682            imodi = surfdataqc%mi(jobs) - int(imaxifp/2) + ji - 1 
     683             
     684            !Deal with wrap around in longitude 
     685            IF ( imodi < 1      ) imodi = imodi + jpiglo 
     686            IF ( imodi > jpiglo ) imodi = imodi - jpiglo 
     687             
     688            DO jj = 0, imaxjfp 
     689               imodj = surfdataqc%mj(jobs) - int(imaxjfp/2) + jj - 1 
     690               !If model values are out of the domain to the north/south then 
     691               !set them to be the edge of the domain 
     692               IF ( imodj < 1      ) imodj = 1 
     693               IF ( imodj > jpjglo ) imodj = jpjglo 
     694 
     695               igrdip1(ji+1,jj+1,iobs) = imodi 
     696               igrdjp1(ji+1,jj+1,iobs) = imodj 
     697                
     698               IF ( ji >= 1 .AND. jj >= 1 ) THEN 
     699                  igrdi(ji,jj,iobs) = imodi 
     700                  igrdj(ji,jj,iobs) = imodj 
     701               ENDIF 
     702                
     703            END DO 
     704         END DO 
    1241705      END DO 
    1242706 
    1243       CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, & 
     707      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 
    1244708         &                  igrdi, igrdj, glamt, zglam ) 
    1245       CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, & 
     709      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 
    1246710         &                  igrdi, igrdj, gphit, zgphi ) 
    1247       CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, & 
     711      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 
    1248712         &                  igrdi, igrdj, psurfmask, zmask ) 
    1249       CALL obs_int_comm_2d( 2, 2, isurf, kpi, kpj, & 
     713      CALL obs_int_comm_2d( imaxifp, imaxjfp, isurf, kpi, kpj, & 
    1250714         &                  igrdi, igrdj, psurf, zsurf ) 
     715      CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, & 
     716         &                  igrdip1, igrdjp1, glamf, zglamf ) 
     717      CALL obs_int_comm_2d( imaxifp+1, imaxjfp+1, isurf, kpi, kpj, & 
     718         &                  igrdip1, igrdjp1, gphif, zgphif ) 
    1251719 
    1252720      ! 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 
     721      IF ( idayend == 0 .AND. ldnightav ) THEN 
     722 
     723         ALLOCATE( & 
     724            & zsurfm(imaxifp,imaxjfp,isurf)  & 
     725            & ) 
     726 
     727         CALL obs_int_comm_2d( imaxifp,imaxjfp, isurf, kpi, kpj, igrdi, igrdj, & 
     728            &               surfdataqc%vdmean(:,:), zsurfm ) 
     729 
    1264730      ENDIF 
    1265731 
    … …  
    1290756         zphi = surfdataqc%rphi(jobs) 
    1291757 
    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 
    1298758         IF ( ldnightav .AND. idayend == 0 ) THEN 
    1299759            ! Night-time averaged data 
    1300             CALL obs_int_h2d( 1, 1, zweig, zsurfm(:,:,iobs), zext ) 
     760            zsurftmp(:,:,iobs) = zsurfm(:,:,iobs) 
    1301761         ELSE 
    1302             CALL obs_int_h2d( 1, 1, zweig, zsurf(:,:,iobs),  zext ) 
     762            zsurftmp(:,:,iobs) = zsurf(:,:,iobs) 
     763         ENDIF 
     764 
     765         IF ( k2dint <= 4 ) THEN 
     766 
     767            ! Get weights to interpolate the model value to the observation point 
     768            CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         & 
     769               &                   zglam(:,:,iobs), zgphi(:,:,iobs), & 
     770               &                   zmask(:,:,iobs), zweig, zobsmask ) 
     771 
     772            ! Interpolate the model value to the observation point  
     773            CALL obs_int_h2d( 1, 1, zweig, zsurftmp(:,:,iobs), zext ) 
     774 
     775         ELSE 
     776 
     777            ! Get weights to average the model SLA to the observation footprint 
     778            CALL obs_avg_h2d_init( 1, 1, imaxifp, imaxjfp, k2dint, zlam,  zphi, & 
     779               &                   zglam(:,:,iobs), zgphi(:,:,iobs), & 
     780               &                   zglamf(:,:,iobs), zgphif(:,:,iobs), & 
     781               &                   zmask(:,:,iobs), plamscl, pphiscl, & 
     782               &                   lindegrees, zweig, zobsmask ) 
     783 
     784            ! Average the model SST to the observation footprint 
     785            CALL obs_avg_h2d( 1, 1, imaxifp, imaxjfp, & 
     786               &              zweig, zsurftmp(:,:,iobs),  zext ) 
     787 
    1303788         ENDIF 
    1304789 
    … …  
    1310795            surfdataqc%rmod(jobs,1) = zext(1) 
    1311796         ENDIF 
     797          
     798         IF ( zext(1) == obfillflt ) THEN 
     799            ! If the observation value is a fill value, set QC flag to bad 
     800            surfdataqc%nqc(jobs) = 4 
     801         ENDIF 
    1312802 
    1313803      END DO 
    … …  
    1315805      ! Deallocate the data for interpolation 
    1316806      DEALLOCATE( & 
     807         & zweig, & 
    1317808         & igrdi, & 
    1318809         & igrdj, & 
    … …  
    1320811         & zgphi, & 
    1321812         & zmask, & 
    1322          & zsurf  & 
     813         & zsurf, & 
     814         & zsurftmp, & 
     815         & zglamf, & 
     816         & zgphif, & 
     817         & igrdip1,& 
     818         & igrdjp1 & 
    1323819         & ) 
    1324820 
    1325821      ! 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 
     822      IF ( idayend == 0 .AND. ldnightav ) THEN 
     823         DEALLOCATE( & 
     824            & zsurfm  & 
     825            & ) 
    1332826      ENDIF 
    1333827 

  • NEMO/branches/UKMO/r8395_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90

    r11350 r11361  
    1313   !!   obs_sor       : Sort the observation arrays 
    1414   !!--------------------------------------------------------------------- 
    15    USE par_kind, ONLY : wp ! Precision variables 
     15   !! * Modules used 
     16   USE par_kind, ONLY : & ! Precision variables 
     17      & wp    
    1618   USE in_out_manager     ! I/O manager 
    1719   USE obs_profiles_def   ! Definitions for storage arrays for profiles 
    … …  
    2224   USE obs_inter_sup      ! Interpolation support 
    2325   USE obs_oper           ! Observation operators 
    24    USE lib_mpp, ONLY :   ctl_warn, ctl_stop 
     26#if defined key_bdy 
     27   USE bdy_oce, ONLY : &        ! Boundary information 
     28      idx_bdy, nb_bdy 
     29#endif 
     30   USE lib_mpp, ONLY : & 
     31      & ctl_warn, ctl_stop 
    2532 
    2633   IMPLICIT NONE 
     34 
     35   !! * Routine accessibility 
    2736   PRIVATE 
    2837 
    29    PUBLIC   obs_pre_prof     ! First level check and screening of profile obs 
    30    PUBLIC   obs_pre_surf     ! First level check and screening of surface obs 
    31    PUBLIC   calc_month_len   ! Calculate the number of days in the months of a year 
     38   PUBLIC & 
     39      & obs_pre_prof, &    ! First level check and screening of profile obs 
     40      & obs_pre_surf, &    ! First level check and screening of surface obs 
     41      & calc_month_len     ! Calculate the number of days in the months of a year 
    3242 
    3343   !!---------------------------------------------------------------------- 
    … …  
    3747   !!---------------------------------------------------------------------- 
    3848 
    39  
    4049CONTAINS 
    4150 
    42    SUBROUTINE obs_pre_surf( surfdata, surfdataqc, ld_nea ) 
     51   SUBROUTINE obs_pre_surf( surfdata, surfdataqc, ld_nea, ld_bound_reject, & 
     52                            kqc_cutoff ) 
    4353      !!---------------------------------------------------------------------- 
    4454      !!                    ***  ROUTINE obs_pre_sla  *** 
    … …  
    5767      !!        !  2015-02  (M. Martin) Combined routine for surface types. 
    5868      !!---------------------------------------------------------------------- 
     69      !! * Modules used 
    5970      USE par_oce             ! Ocean parameters 
    60       USE dom_oce, ONLY       :   glamt, gphit, tmask, nproc   ! Geographical information 
     71      USE dom_oce, ONLY : &   ! Geographical information 
     72         & glamt,   & 
     73         & gphit,   & 
     74         & tmask,   & 
     75         & nproc 
    6176      !! * Arguments 
    6277      TYPE(obs_surf), INTENT(INOUT) :: surfdata    ! Full set of surface data 
    63       TYPE(obs_surf), INTENT(INOUT) :: surfdataqc   ! Subset of surface data not failing screening 
    64       LOGICAL, INTENT(IN) :: ld_nea         ! Switch for rejecting observation near land 
    65       ! 
     78      TYPE(obs_surf), INTENT(INOUT) :: surfdataqc  ! Subset of surface data not failing screening 
     79      LOGICAL, INTENT(IN) :: ld_nea                ! Switch for rejecting observation near land 
     80      LOGICAL, INTENT(IN) :: ld_bound_reject       ! Switch for rejecting obs near the boundary 
     81      INTEGER, INTENT(IN), OPTIONAL :: kqc_cutoff   ! cut off for QC value 
     82      !! * Local declarations 
     83      INTEGER :: iqc_cutoff = 255   ! cut off for QC value 
    6684      INTEGER :: iyea0        ! Initial date 
    6785      INTEGER :: imon0        !  - (year, month, day, hour, minute) 
    … …  
    7694      INTEGER :: inlasobs     !  - close to land 
    7795      INTEGER :: igrdobs      !  - fail the grid search 
     96      INTEGER :: ibdysobs     !  - close to open boundary 
    7897                              ! Global counters for observations that 
    7998      INTEGER :: iotdobsmpp     !  - outside time domain 
    … …  
    82101      INTEGER :: inlasobsmpp    !  - close to land 
    83102      INTEGER :: igrdobsmpp     !  - fail the grid search 
    84       LOGICAL, DIMENSION(:), ALLOCATABLE ::   llvalid            ! SLA data selection 
     103      INTEGER :: ibdysobsmpp  !  - close to open boundary 
     104      LOGICAL, DIMENSION(:), ALLOCATABLE :: &  
     105         & llvalid            ! SLA data selection 
    85106      INTEGER :: jobs         ! Obs. loop variable 
    86107      INTEGER :: jstp         ! Time loop variable 
    87108      INTEGER :: inrc         ! Time index variable 
    88       !!---------------------------------------------------------------------- 
    89  
    90       IF(lwp) WRITE(numout,*) 'obs_pre_surf : Preparing the surface observations...' 
    91       IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 
     109 
     110      IF(lwp) WRITE(numout,*)'obs_pre_surf : Preparing the surface observations...' 
     111      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 
    92112       
    93113      ! Initial date initialization (year, month, day, hour, minute) 
    … …  
    95115      imon0 = ( ndate0 - iyea0 * 10000 ) / 100 
    96116      iday0 =   ndate0 - iyea0 * 10000 - imon0 * 100 
    97       ihou0 = nn_time0 / 100 
    98       imin0 = ( nn_time0 - ihou0 * 100 ) 
     117      ihou0 = 0 
     118      imin0 = 0 
    99119 
    100120      icycle = no     ! Assimilation cycle 
    … …  
    107127      ilansobs = 0 
    108128      inlasobs = 0 
     129      ibdysobs = 0  
     130 
     131      ! Set QC cutoff to optional value if provided 
     132      IF ( PRESENT(kqc_cutoff) ) iqc_cutoff=kqc_cutoff 
    109133 
    110134      ! ----------------------------------------------------------------------- 
    … …  
    140164         &                 tmask(:,:,1), surfdata%nqc,  & 
    141165         &                 iosdsobs,     ilansobs,     & 
    142          &                 inlasobs,     ld_nea        ) 
     166         &                 inlasobs,     ld_nea,       & 
     167         &                 ibdysobs,     ld_bound_reject, & 
     168         &                 iqc_cutoff                     ) 
    143169 
    144170      CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp ) 
    145171      CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp ) 
    146172      CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp ) 
     173      CALL obs_mpp_sum_integer( ibdysobs, ibdysobsmpp ) 
    147174 
    148175      ! ----------------------------------------------------------------------- 
    … …  
    155182      ALLOCATE( llvalid(surfdata%nsurf) ) 
    156183       
    157       ! We want all data which has qc flags <= 10 
    158  
    159       llvalid(:)  = ( surfdata%nqc(:)  <= 10 ) 
     184      ! We want all data which has qc flags <= iqc_cutoff 
     185 
     186      llvalid(:)  = ( surfdata%nqc(:)  <= iqc_cutoff ) 
    160187 
    161188      ! The actual copying 
    … …  
    190217               &            inlasobsmpp 
    191218         ENDIF 
     219         WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data near open boundary (removed) = ', & 
     220            &            ibdysobsmpp   
    192221         WRITE(numout,*) ' '//surfdataqc%cvars(1)//' data accepted                             = ', & 
    193222            &            surfdataqc%nsurfmpp 
    … …  
    222251 
    223252 
    224    SUBROUTINE obs_pre_prof( profdata, prodatqc, ld_var1, ld_var2, & 
     253   SUBROUTINE obs_pre_prof( profdata, prodatqc, ld_var, & 
    225254      &                     kpi, kpj, kpk, & 
    226       &                     zmask1, pglam1, pgphi1, zmask2, pglam2, pgphi2,  & 
    227       &                     ld_nea, kdailyavtypes ) 
     255      &                     zmask, pglam, pgphi,  & 
     256      &                     ld_nea, ld_bound_reject, kdailyavtypes,  kqc_cutoff ) 
    228257 
    229258!!---------------------------------------------------------------------- 
    … …  
    241270      !! 
    242271      !!---------------------------------------------------------------------- 
    243       USE par_oce           ! Ocean parameters 
    244       USE dom_oce, ONLY :   gdept_1d, nproc   ! Geographical information 
     272      !! * Modules used 
     273      USE par_oce             ! Ocean parameters 
     274      USE dom_oce, ONLY : &   ! Geographical information 
     275         & gdept_1d,             & 
     276         & nproc 
    245277 
    246278      !! * Arguments 
    247279      TYPE(obs_prof), INTENT(INOUT) :: profdata   ! Full set of profile data 
    248280      TYPE(obs_prof), INTENT(INOUT) :: prodatqc   ! Subset of profile data not failing screening 
    249       LOGICAL, INTENT(IN) :: ld_var1              ! Observed variables switches 
    250       LOGICAL, INTENT(IN) :: ld_var2 
     281      LOGICAL, DIMENSION(profdata%nvar), INTENT(IN) :: & 
     282         & ld_var                                 ! Observed variables switches 
    251283      LOGICAL, INTENT(IN) :: ld_nea               ! Switch for rejecting observation near land 
     284      LOGICAL, INTENT(IN) :: ld_bound_reject      ! Switch for rejecting observations near the boundary 
    252285      INTEGER, INTENT(IN) :: kpi, kpj, kpk        ! Local domain sizes 
    253286      INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 
    254287         & kdailyavtypes                          ! Types for daily averages 
    255       REAL(wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 
    256          & zmask1, & 
    257          & zmask2 
    258       REAL(wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 
    259          & pglam1, & 
    260          & pglam2, & 
    261          & pgphi1, & 
    262          & pgphi2 
     288      REAL(wp), INTENT(IN), DIMENSION(kpi,kpj,kpk,profdata%nvar) :: & 
     289         & zmask 
     290      REAL(wp), INTENT(IN), DIMENSION(kpi,kpj,profdata%nvar) :: & 
     291         & pglam, & 
     292         & pgphi 
     293      INTEGER, INTENT(IN), OPTIONAL :: kqc_cutoff   ! cut off for QC value 
    263294 
    264295      !! * Local declarations 
     296      INTEGER :: iqc_cutoff = 255   ! cut off for QC value 
    265297      INTEGER :: iyea0        ! Initial date 
    266298      INTEGER :: imon0        !  - (year, month, day, hour, minute) 
    … …  
    269301      INTEGER :: imin0 
    270302      INTEGER :: icycle       ! Current assimilation cycle 
    271                               ! Counters for observations that are 
    272       INTEGER :: iotdobs      !  - outside time domain 
    273       INTEGER :: iosdv1obs    !  - outside space domain (variable 1) 
    274       INTEGER :: iosdv2obs    !  - outside space domain (variable 2) 
    275       INTEGER :: ilanv1obs    !  - within a model land cell (variable 1) 
    276       INTEGER :: ilanv2obs    !  - within a model land cell (variable 2) 
    277       INTEGER :: inlav1obs    !  - close to land (variable 1) 
    278       INTEGER :: inlav2obs    !  - close to land (variable 2) 
    279       INTEGER :: igrdobs      !  - fail the grid search 
    280       INTEGER :: iuvchku      !  - reject u if v rejected and vice versa 
    281       INTEGER :: iuvchkv      ! 
    282                               ! Global counters for observations that are 
    283       INTEGER :: iotdobsmpp   !  - outside time domain 
    284       INTEGER :: iosdv1obsmpp !  - outside space domain (variable 1) 
    285       INTEGER :: iosdv2obsmpp !  - outside space domain (variable 2) 
    286       INTEGER :: ilanv1obsmpp !  - within a model land cell (variable 1) 
    287       INTEGER :: ilanv2obsmpp !  - within a model land cell (variable 2) 
    288       INTEGER :: inlav1obsmpp !  - close to land (variable 1) 
    289       INTEGER :: inlav2obsmpp !  - close to land (variable 2) 
    290       INTEGER :: igrdobsmpp   !  - fail the grid search 
    291       INTEGER :: iuvchkumpp   !  - reject var1 if var2 rejected and vice versa 
    292       INTEGER :: iuvchkvmpp   ! 
     303                                                       ! Counters for observations that are 
     304      INTEGER                           :: iotdobs     !  - outside time domain 
     305      INTEGER, DIMENSION(profdata%nvar) :: iosdvobs    !  - outside space domain 
     306      INTEGER, DIMENSION(profdata%nvar) :: ilanvobs    !  - within a model land cell 
     307      INTEGER, DIMENSION(profdata%nvar) :: inlavobs    !  - close to land 
     308      INTEGER, DIMENSION(profdata%nvar) :: ibdyvobs    !  - boundary    
     309      INTEGER                           :: igrdobs     !  - fail the grid search 
     310      INTEGER                           :: iuvchku     !  - reject UVEL if VVEL rejected 
     311      INTEGER                           :: iuvchkv     !  - reject VVEL if UVEL rejected 
     312                                                       ! Global counters for observations that are 
     313      INTEGER                           :: iotdobsmpp  !  - outside time domain 
     314      INTEGER, DIMENSION(profdata%nvar) :: iosdvobsmpp !  - outside space domain 
     315      INTEGER, DIMENSION(profdata%nvar) :: ilanvobsmpp !  - within a model land cell 
     316      INTEGER, DIMENSION(profdata%nvar) :: inlavobsmpp !  - close to land 
     317      INTEGER, DIMENSION(profdata%nvar) :: ibdyvobsmpp !  - boundary 
     318      INTEGER :: igrdobsmpp                            !  - fail the grid search 
     319      INTEGER :: iuvchkumpp                            !  - reject UVEL if VVEL rejected 
     320      INTEGER :: iuvchkvmpp                            !  - reject VVEL if UVEL rejected 
    293321      TYPE(obs_prof_valid) ::  llvalid      ! Profile selection  
    294322      TYPE(obs_prof_valid), DIMENSION(profdata%nvar) :: & 
    295          & llvvalid           ! var1,var2 selection  
     323         & llvvalid           ! vars selection  
    296324      INTEGER :: jvar         ! Variable loop variable 
    297325      INTEGER :: jobs         ! Obs. loop variable 
    298326      INTEGER :: jstp         ! Time loop variable 
    299327      INTEGER :: inrc         ! Time index variable 
    300       !!---------------------------------------------------------------------- 
     328      CHARACTER(LEN=256) :: cout1  ! Diagnostic output line 
     329      CHARACTER(LEN=256) :: cout2  ! Diagnostic output line 
    301330 
    302331      IF(lwp) WRITE(numout,*)'obs_pre_prof: Preparing the profile data...' 
    … …  
    307336      imon0 = ( ndate0 - iyea0 * 10000 ) / 100 
    308337      iday0 =   ndate0 - iyea0 * 10000 - imon0 * 100 
    309       ihou0 = nn_time0 / 100 
    310       imin0 = ( nn_time0 - ihou0 * 100 ) 
     338      ihou0 = 0 
     339      imin0 = 0 
    311340 
    312341      icycle = no     ! Assimilation cycle 
    313342 
    314       ! Diagnotics counters for various failures. 
    315  
    316       iotdobs  = 0 
    317       igrdobs  = 0 
    318       iosdv1obs = 0 
    319       iosdv2obs = 0 
    320       ilanv1obs = 0 
    321       ilanv2obs = 0 
    322       inlav1obs = 0 
    323       inlav2obs = 0 
    324       iuvchku  = 0 
    325       iuvchkv = 0 
     343      ! Diagnostics counters for various failures. 
     344 
     345      iotdobs     = 0 
     346      igrdobs     = 0 
     347      iosdvobs(:) = 0 
     348      ilanvobs(:) = 0 
     349      inlavobs(:) = 0 
     350      ibdyvobs(:) = 0 
     351      iuvchku     = 0 
     352      iuvchkv     = 0 
     353 
     354 
     355      ! Set QC cutoff to optional value if provided 
     356      IF ( PRESENT(kqc_cutoff) ) iqc_cutoff=kqc_cutoff 
    326357 
    327358      ! ----------------------------------------------------------------------- 
    … …  
    335366            &              profdata%nday,    profdata%nhou, profdata%nmin, & 
    336367            &              profdata%ntyp,    profdata%nqc,  profdata%mstp, & 
    337             &              iotdobs, kdailyavtypes = kdailyavtypes ) 
     368            &              iotdobs, kdailyavtypes = kdailyavtypes,         & 
     369            &              kqc_cutoff = iqc_cutoff ) 
    338370      ELSE 
    339371         CALL obs_coo_tim_prof( icycle, & 
    … …  
    342374            &              profdata%nday,    profdata%nhou, profdata%nmin, & 
    343375            &              profdata%ntyp,    profdata%nqc,  profdata%mstp, & 
    344             &              iotdobs ) 
     376            &              iotdobs,          kqc_cutoff = iqc_cutoff ) 
    345377      ENDIF 
    346378 
    … …  
    351383      ! ----------------------------------------------------------------------- 
    352384 
    353       CALL obs_coo_grd( profdata%nprof,   profdata%mi(:,1), profdata%mj(:,1), & 
    354          &              profdata%nqc,     igrdobs                         ) 
    355       CALL obs_coo_grd( profdata%nprof,   profdata%mi(:,2), profdata%mj(:,2), & 
    356          &              profdata%nqc,     igrdobs                         ) 
     385      DO jvar = 1, profdata%nvar 
     386         CALL obs_coo_grd( profdata%nprof,   profdata%mi(:,jvar), profdata%mj(:,jvar), & 
     387            &              profdata%nqc,     igrdobs ) 
     388      END DO 
    357389 
    358390      CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp ) 
    359391 
    360392      ! ----------------------------------------------------------------------- 
    361       ! Reject all observations for profiles with nqc > 10 
    362       ! ----------------------------------------------------------------------- 
    363  
    364       CALL obs_pro_rej( profdata ) 
     393      ! Reject all observations for profiles with nqc > iqc_cutoff 
     394      ! ----------------------------------------------------------------------- 
     395 
     396      CALL obs_pro_rej( profdata, kqc_cutoff = iqc_cutoff ) 
    365397 
    366398      ! ----------------------------------------------------------------------- 
    … …  
    369401      ! ----------------------------------------------------------------------- 
    370402 
    371       ! Variable 1 
    372       CALL obs_coo_spc_3d( profdata%nprof,        profdata%nvprot(1),   & 
    373          &                 profdata%npvsta(:,1),  profdata%npvend(:,1), & 
    374          &                 jpi,                   jpj,                  & 
    375          &                 jpk,                                         & 
    376          &                 profdata%mi,           profdata%mj,          & 
    377          &                 profdata%var(1)%mvk,                         & 
    378          &                 profdata%rlam,         profdata%rphi,        & 
    379          &                 profdata%var(1)%vdep,                        & 
    380          &                 pglam1,                pgphi1,               & 
    381          &                 gdept_1d,              zmask1,               & 
    382          &                 profdata%nqc,          profdata%var(1)%nvqc, & 
    383          &                 iosdv1obs,              ilanv1obs,           & 
    384          &                 inlav1obs,              ld_nea                ) 
    385  
    386       CALL obs_mpp_sum_integer( iosdv1obs, iosdv1obsmpp ) 
    387       CALL obs_mpp_sum_integer( ilanv1obs, ilanv1obsmpp ) 
    388       CALL obs_mpp_sum_integer( inlav1obs, inlav1obsmpp ) 
    389  
    390       ! Variable 2 
    391       CALL obs_coo_spc_3d( profdata%nprof,        profdata%nvprot(2),   & 
    392          &                 profdata%npvsta(:,2),  profdata%npvend(:,2), & 
    393          &                 jpi,                   jpj,                  & 
    394          &                 jpk,                                         & 
    395          &                 profdata%mi,           profdata%mj,          &  
    396          &                 profdata%var(2)%mvk,                         & 
    397          &                 profdata%rlam,         profdata%rphi,        & 
    398          &                 profdata%var(2)%vdep,                        & 
    399          &                 pglam2,                pgphi2,               & 
    400          &                 gdept_1d,              zmask2,               & 
    401          &                 profdata%nqc,          profdata%var(2)%nvqc, & 
    402          &                 iosdv2obs,              ilanv2obs,           & 
    403          &                 inlav2obs,              ld_nea                ) 
    404  
    405       CALL obs_mpp_sum_integer( iosdv2obs, iosdv2obsmpp ) 
    406       CALL obs_mpp_sum_integer( ilanv2obs, ilanv2obsmpp ) 
    407       CALL obs_mpp_sum_integer( inlav2obs, inlav2obsmpp ) 
     403      DO jvar = 1, profdata%nvar 
     404         CALL obs_coo_spc_3d( profdata%nprof,          profdata%nvprot(jvar),   & 
     405            &                 profdata%npvsta(:,jvar), profdata%npvend(:,jvar), & 
     406            &                 jpi,                     jpj,                     & 
     407            &                 jpk,                                              & 
     408            &                 profdata%mi,             profdata%mj,             & 
     409            &                 profdata%var(jvar)%mvk,                           & 
     410            &                 profdata%rlam,           profdata%rphi,           & 
     411            &                 profdata%var(jvar)%vdep,                          & 
     412            &                 pglam(:,:,jvar),         pgphi(:,:,jvar),         & 
     413            &                 gdept_1d,                zmask(:,:,:,jvar),       & 
     414            &                 profdata%nqc,            profdata%var(jvar)%nvqc, & 
     415            &                 iosdvobs(jvar),          ilanvobs(jvar),          & 
     416            &                 inlavobs(jvar),          ld_nea,                  & 
     417            &                 ibdyvobs(jvar),          ld_bound_reject,         & 
     418            &                 iqc_cutoff       ) 
     419 
     420         CALL obs_mpp_sum_integer( iosdvobs(jvar), iosdvobsmpp(jvar) ) 
     421         CALL obs_mpp_sum_integer( ilanvobs(jvar), ilanvobsmpp(jvar) ) 
     422         CALL obs_mpp_sum_integer( inlavobs(jvar), inlavobsmpp(jvar) ) 
     423         CALL obs_mpp_sum_integer( ibdyvobs(jvar), ibdyvobsmpp(jvar) ) 
     424      END DO 
    408425 
    409426      ! ----------------------------------------------------------------------- 
    … …  
    412429 
    413430      IF ( TRIM(profdata%cvars(1)) == 'UVEL' ) THEN 
    414          CALL obs_uv_rej( profdata, iuvchku, iuvchkv ) 
     431         CALL obs_uv_rej( profdata, iuvchku, iuvchkv, iqc_cutoff ) 
    415432         CALL obs_mpp_sum_integer( iuvchku, iuvchkumpp ) 
    416433         CALL obs_mpp_sum_integer( iuvchkv, iuvchkvmpp ) 
    … …  
    429446      END DO 
    430447 
    431       ! We want all data which has qc flags = 0 
    432  
    433       llvalid%luse(:) = ( profdata%nqc(:)  <= 10 ) 
     448      ! We want all data which has qc flags <= iqc_cutoff 
     449 
     450      llvalid%luse(:) = ( profdata%nqc(:)  <= iqc_cutoff ) 
    434451      DO jvar = 1,profdata%nvar 
    435          llvvalid(jvar)%luse(:) = ( profdata%var(jvar)%nvqc(:) <= 10 ) 
     452         llvvalid(jvar)%luse(:) = ( profdata%var(jvar)%nvqc(:) <= iqc_cutoff ) 
    436453      END DO 
    437454 
    … …  
    456473       
    457474         WRITE(numout,*) 
    458          WRITE(numout,*) ' Profiles outside time domain                     = ', & 
     475         WRITE(numout,*) ' Profiles outside time domain                       = ', & 
    459476            &            iotdobsmpp 
    460          WRITE(numout,*) ' Remaining profiles that failed grid search       = ', & 
     477         WRITE(numout,*) ' Remaining profiles that failed grid search         = ', & 
    461478            &            igrdobsmpp 
    462          WRITE(numout,*) ' Remaining '//prodatqc%cvars(1)//' data outside space domain       = ', & 
    463             &            iosdv1obsmpp 
    464          WRITE(numout,*) ' Remaining '//prodatqc%cvars(1)//' data at land points             = ', & 
    465             &            ilanv1obsmpp 
    466          IF (ld_nea) THEN 
    467             WRITE(numout,*) ' Remaining '//prodatqc%cvars(1)//' data near land points (removed) = ',& 
    468                &            inlav1obsmpp 
    469          ELSE 
    470             WRITE(numout,*) ' Remaining '//prodatqc%cvars(1)//' data near land points (kept)    = ',& 
    471                &            inlav1obsmpp 
    472          ENDIF 
    473          IF ( TRIM(profdata%cvars(1)) == 'UVEL' ) THEN 
    474             WRITE(numout,*) ' U observation rejected since V rejected     = ', & 
    475                &            iuvchku 
    476          ENDIF 
    477          WRITE(numout,*) ' '//prodatqc%cvars(1)//' data accepted                             = ', & 
    478             &            prodatqc%nvprotmpp(1) 
    479          WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data outside space domain       = ', & 
    480             &            iosdv2obsmpp 
    481          WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data at land points             = ', & 
    482             &            ilanv2obsmpp 
    483          IF (ld_nea) THEN 
    484             WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data near land points (removed) = ',& 
    485                &            inlav2obsmpp 
    486          ELSE 
    487             WRITE(numout,*) ' Remaining '//prodatqc%cvars(2)//' data near land points (kept)    = ',& 
    488                &            inlav2obsmpp 
    489          ENDIF 
    490          IF ( TRIM(profdata%cvars(1)) == 'UVEL' ) THEN 
    491             WRITE(numout,*) ' V observation rejected since U rejected     = ', & 
    492                &            iuvchkv 
    493          ENDIF 
    494          WRITE(numout,*) ' '//prodatqc%cvars(2)//' data accepted                             = ', & 
    495             &            prodatqc%nvprotmpp(2) 
     479         DO jvar = 1, profdata%nvar 
     480            WRITE(numout,*) ' Remaining '//prodatqc%cvars(jvar)//' data outside space domain       = ', & 
     481               &            iosdvobsmpp(jvar) 
     482            WRITE(numout,*) ' Remaining '//prodatqc%cvars(jvar)//' data at land points             = ', & 
     483               &            ilanvobsmpp(jvar) 
     484            IF (ld_nea) THEN 
     485               WRITE(numout,*) ' Remaining '//prodatqc%cvars(jvar)//' data near land points (removed) = ',& 
     486                  &            inlavobsmpp(jvar) 
     487            ELSE 
     488               WRITE(numout,*) ' Remaining '//prodatqc%cvars(jvar)//' data near land points (kept)    = ',& 
     489                  &            inlavobsmpp(jvar) 
     490            ENDIF 
     491            IF ( TRIM(profdata%cvars(jvar)) == 'UVEL' ) THEN 
     492               WRITE(numout,*) ' U observation rejected since V rejected     = ', & 
     493                  &            iuvchku 
     494            ELSE IF ( TRIM(profdata%cvars(jvar)) == 'VVEL' ) THEN 
     495               WRITE(numout,*) ' V observation rejected since U rejected     = ', & 
     496                  &            iuvchkv 
     497            ENDIF 
     498            WRITE(numout,*) ' Remaining '//prodatqc%cvars(jvar)//' data near open boundary (removed) = ',& 
     499                  &            ibdyvobsmpp(jvar) 
     500            WRITE(numout,*) ' '//prodatqc%cvars(jvar)//' data accepted                             = ', & 
     501               &            prodatqc%nvprotmpp(jvar) 
     502         END DO 
    496503 
    497504         WRITE(numout,*) 
    498505         WRITE(numout,*) ' Number of observations per time step :' 
    499506         WRITE(numout,*) 
    500          WRITE(numout,'(10X,A,5X,A,5X,A,A)')'Time step','Profiles', & 
    501             &                               '     '//prodatqc%cvars(1)//'     ', & 
    502             &                               '     '//prodatqc%cvars(2)//'     ' 
    503          WRITE(numout,998) 
     507         WRITE(cout1,'(10X,A9,5X,A8)') 'Time step', 'Profiles' 
     508         WRITE(cout2,'(10X,A9,5X,A8)') '---------', '--------' 
     509         DO jvar = 1, prodatqc%nvar 
     510            WRITE(cout1,'(A,5X,A11)') TRIM(cout1), TRIM(prodatqc%cvars(jvar)) 
     511            WRITE(cout2,'(A,5X,A11)') TRIM(cout2), '-----------' 
     512         END DO 
     513         WRITE(numout,*) cout1 
     514         WRITE(numout,*) cout2 
    504515      ENDIF 
    505516       
    … …  
    528539         DO jstp = nit000 - 1, nitend 
    529540            inrc = jstp - nit000 + 2 
    530             WRITE(numout,999) jstp, prodatqc%npstpmpp(inrc), & 
    531                &                    prodatqc%nvstpmpp(inrc,1), & 
    532                &                    prodatqc%nvstpmpp(inrc,2) 
     541            WRITE(cout1,'(10X,I9,5X,I8)') jstp, prodatqc%npstpmpp(inrc) 
     542            DO jvar = 1, prodatqc%nvar 
     543               WRITE(cout1,'(A,5X,I11)') TRIM(cout1), prodatqc%nvstpmpp(inrc,jvar) 
     544            END DO 
     545            WRITE(numout,*) cout1 
    533546         END DO 
    534547      ENDIF 
    535  
    536 998   FORMAT(10X,'---------',5X,'--------',5X,'-----------',5X,'----------------') 
    537 999   FORMAT(10X,I9,5X,I8,5X,I11,5X,I8) 
    538548 
    539549   END SUBROUTINE obs_pre_prof 
    … …  
    644654            &        .AND. ( kobsmin(jobs) <= kmin0 ) ) ) THEN 
    645655            kobsstp(jobs) = -1 
    646             kobsqc(jobs)  = kobsqc(jobs) + 11 
     656            kobsqc(jobs)  = IBSET(kobsqc(jobs),13) 
    647657            kotdobs       = kotdobs + 1 
    648658            CYCLE 
    … …  
    695705         IF ( ( kobsstp(jobs) < ( nit000 - 1 ) ) & 
    696706            & .OR.( kobsstp(jobs) > nitend ) ) THEN 
    697             kobsqc(jobs) = kobsqc(jobs) + 12 
     707            kobsqc(jobs) = IBSET(kobsqc(jobs),13) 
    698708            kotdobs = kotdobs + 1 
    699709            CYCLE 
    … …  
    739749      &                    kobsno,                                        & 
    740750      &                    kobsyea, kobsmon, kobsday, kobshou, kobsmin,   & 
    741       &                    ktyp,    kobsqc,  kobsstp, kotdobs, kdailyavtypes ) 
     751      &                    ktyp,    kobsqc,  kobsstp, kotdobs, kdailyavtypes, & 
     752      &                    kqc_cutoff ) 
    742753      !!---------------------------------------------------------------------- 
    743754      !!                    ***  ROUTINE obs_coo_tim *** 
    … …  
    783794      INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 
    784795         & kdailyavtypes    ! Types for daily averages 
     796      INTEGER, OPTIONAL, INTENT(IN) :: kqc_cutoff           ! QC cutoff value 
     797 
    785798      !! * Local declarations 
    786799      INTEGER :: jobs 
     800      INTEGER :: iqc_cutoff=255 
    787801 
    788802      !----------------------------------------------------------------------- 
    … …  
    803817         DO jobs = 1, kobsno 
    804818             
    805             IF ( kobsqc(jobs) <= 10 ) THEN 
     819            IF ( kobsqc(jobs) <= iqc_cutoff ) THEN 
    806820                
    807821               IF ( ( kobsstp(jobs) == (nit000 - 1) ).AND.& 
    808822                  & ( ANY (kdailyavtypes(:) == ktyp(jobs)) ) ) THEN 
    809                   kobsqc(jobs) = kobsqc(jobs) + 14 
     823                  kobsqc(jobs) = IBSET(kobsqc(jobs),13) 
    810824                  kotdobs      = kotdobs + 1 
    811825                  CYCLE 
    … …  
    850864      DO jobs = 1, kobsno 
    851865         IF ( ( kobsi(jobs) <= 0 ) .AND. ( kobsj(jobs) <= 0 ) ) THEN 
    852             kobsqc(jobs) = kobsqc(jobs) + 18 
     866            kobsqc(jobs) = IBSET(kobsqc(jobs),12) 
    853867            kgrdobs = kgrdobs + 1 
    854868         ENDIF 
    … …  
    861875      &                       plam,   pphi,    pmask,            & 
    862876      &                       kobsqc, kosdobs, klanobs,          & 
    863       &                       knlaobs,ld_nea                     ) 
     877      &                       knlaobs,ld_nea,                    & 
     878      &                       kbdyobs,ld_bound_reject,           & 
     879      &                       kqc_cutoff                         ) 
    864880      !!---------------------------------------------------------------------- 
    865881      !!                    ***  ROUTINE obs_coo_spc_2d  *** 
    … …  
    894910      INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: & 
    895911         & 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 
     912      INTEGER, INTENT(INOUT) :: kosdobs          ! Observations outside space domain 
     913      INTEGER, INTENT(INOUT) :: klanobs          ! Observations within a model land cell 
     914      INTEGER, INTENT(INOUT) :: knlaobs          ! Observations near land 
     915      INTEGER, INTENT(INOUT) :: kbdyobs          ! Observations near boundary 
     916      LOGICAL, INTENT(IN)    :: ld_nea           ! Flag observations near land 
     917      LOGICAL, INTENT(IN)    :: ld_bound_reject  ! Flag observations near open boundary  
     918      INTEGER, INTENT(IN)    :: kqc_cutoff       ! Cutoff QC value 
     919 
    900920      !! * Local declarations 
    901921      REAL(KIND=wp), DIMENSION(2,2,kobsno) :: & 
    902922         & zgmsk              ! Grid mask 
     923#if defined key_bdy  
     924      REAL(KIND=wp), DIMENSION(2,2,kobsno) :: & 
     925         & zbmsk              ! Boundary mask 
     926      REAL(KIND=wp), DIMENSION(jpi,jpj) :: zbdymask 
     927#endif  
    903928      REAL(KIND=wp), DIMENSION(2,2,kobsno) :: & 
    904929         & zglam, &           ! Model longitude at grid points 
    … …  
    917942         ! For invalid points use 2,2 
    918943 
    919          IF ( kobsqc(jobs) >= 10 ) THEN 
     944         IF ( kobsqc(jobs) >= kqc_cutoff ) THEN 
    920945 
    921946            igrdi(1,1,jobs) = 1 
    … …  
    942967 
    943968      END DO 
     969 
     970#if defined key_bdy              
     971      ! Create a mask grid points in boundary rim 
     972      IF (ld_bound_reject) THEN 
     973         zbdymask(:,:) = 1.0_wp 
     974         DO ji = 1, nb_bdy 
     975            DO jj = 1, idx_bdy(ji)%nblen(1) 
     976               zbdymask(idx_bdy(ji)%nbi(jj,1),idx_bdy(ji)%nbj(jj,1)) = 0.0_wp 
     977            ENDDO 
     978         ENDDO 
     979  
     980         CALL obs_int_comm_2d( 2, 2, kobsno, kpi, kpj, igrdi, igrdj, zbdymask, zbmsk ) 
     981      ENDIF 
     982#endif        
    944983       
    945984      CALL obs_int_comm_2d( 2, 2, kobsno, kpi, kpj, igrdi, igrdj, pmask, zgmsk ) 
    … …  
    950989 
    951990         ! Skip bad observations 
    952          IF ( kobsqc(jobs) >= 10 ) CYCLE 
     991         IF ( kobsqc(jobs) >= kqc_cutoff ) CYCLE 
    953992 
    954993         ! Flag if the observation falls outside the model spatial domain 
    … …  
    957996            &  .OR. ( pobsphi(jobs) <  -90. ) & 
    958997            &  .OR. ( pobsphi(jobs) >   90. ) ) THEN 
    959             kobsqc(jobs) = kobsqc(jobs) + 11 
     998            kobsqc(jobs) = IBSET(kobsqc(jobs),11) 
    960999            kosdobs = kosdobs + 1 
    9611000            CYCLE 
    … …  
    9641003         ! Flag if the observation falls with a model land cell 
    9651004         IF ( SUM( zgmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN 
    966             kobsqc(jobs) = kobsqc(jobs)  + 12 
     1005            kobsqc(jobs) = IBSET(kobsqc(jobs),10) 
    9671006            klanobs = klanobs + 1 
    9681007            CYCLE 
    … …  
    9781017               IF ( ( ABS( zgphi(ji,jj,jobs) - pobsphi(jobs) ) < 1.0e-6_wp ) & 
    9791018                  & .AND. & 
    980                   & ( ABS( zglam(ji,jj,jobs) - pobslam(jobs) ) < 1.0e-6_wp ) & 
    981                   & ) THEN 
     1019                  & ( ABS( MOD( zglam(ji,jj,jobs) - pobslam(jobs),360.0) ) & 
     1020                  & < 1.0e-6_wp ) ) THEN 
    9821021                  lgridobs = .TRUE. 
    9831022                  iig = ji 
    … …  
    9861025            END DO 
    9871026         END DO 
    988    
    989          ! For observations on the grid reject them if their are at 
    990          ! a masked point 
    991           
     1027  
    9921028         IF (lgridobs) THEN 
    9931029            IF (zgmsk(iig,ijg,jobs) == 0.0_wp ) THEN 
    994                kobsqc(jobs) = kobsqc(jobs) + 12 
     1030               kobsqc(jobs) = IBSET(kobsqc(jobs),10) 
    9951031               klanobs = klanobs + 1 
    9961032               CYCLE 
    9971033            ENDIF 
    9981034         ENDIF 
    999                        
     1035 
     1036  
    10001037         ! Flag if the observation falls is close to land 
    10011038         IF ( MINVAL( zgmsk(1:2,1:2,jobs) ) == 0.0_wp) THEN 
    1002             IF (ld_nea) kobsqc(jobs) = kobsqc(jobs) + 14 
    10031039            knlaobs = knlaobs + 1 
    1004             CYCLE 
     1040            IF (ld_nea) THEN 
     1041               kobsqc(jobs) = IBSET(kobsqc(jobs),9) 
     1042               CYCLE 
     1043            ENDIF 
    10051044         ENDIF 
     1045 
     1046#if defined key_bdy 
     1047         ! Flag if the observation falls close to the boundary rim 
     1048         IF (ld_bound_reject) THEN 
     1049            IF ( MINVAL( zbmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN 
     1050               kobsqc(jobs) = IBSET(kobsqc(jobs),8) 
     1051               kbdyobs = kbdyobs + 1 
     1052               CYCLE 
     1053            ENDIF 
     1054            ! for observations on the grid... 
     1055            IF (lgridobs) THEN 
     1056               IF (zbmsk(iig,ijg,jobs) == 0.0_wp ) THEN 
     1057                  kobsqc(jobs) = IBSET(kobsqc(jobs),8) 
     1058                  kbdyobs = kbdyobs + 1 
     1059                  CYCLE 
     1060               ENDIF 
     1061            ENDIF 
     1062         ENDIF 
     1063#endif  
    10061064             
    10071065      END DO 
    … …  
    10151073      &                       plam,    pphi,    pdep,    pmask, & 
    10161074      &                       kpobsqc, kobsqc,  kosdobs,        & 
    1017       &                       klanobs, knlaobs, ld_nea          ) 
     1075      &                       klanobs, knlaobs, ld_nea,         & 
     1076      &                       kbdyobs, ld_bound_reject,         & 
     1077      &                       kqc_cutoff                        ) 
    10181078      !!---------------------------------------------------------------------- 
    10191079      !!                    ***  ROUTINE obs_coo_spc_3d  *** 
    … …  
    10401100         & gdepw_1d,      & 
    10411101         & gdepw_0,       &                        
    1042          & gdepw_n,       & 
     1102         & gdepw_n,       &  
     1103#if defined key_vvl 
    10431104         & gdept_n,       & 
     1105#endif 
    10441106         & ln_zco,        & 
    1045          & ln_zps              
     1107         & ln_zps,        & 
     1108         & ln_linssh  
    10461109 
    10471110      !! * Arguments 
    … …  
    10771140      INTEGER, INTENT(INOUT) :: klanobs     ! Observations within a model land cell 
    10781141      INTEGER, INTENT(INOUT) :: knlaobs     ! Observations near land 
     1142      INTEGER, INTENT(INOUT) :: kbdyobs     ! Observations near boundary 
    10791143      LOGICAL, INTENT(IN) :: ld_nea         ! Flag observations near land 
     1144      LOGICAL, INTENT(IN) :: ld_bound_reject  ! Flag observations near open boundary 
     1145      INTEGER, INTENT(IN) :: kqc_cutoff     ! Cutoff QC value 
     1146 
    10801147      !! * Local declarations 
    10811148      REAL(KIND=wp), DIMENSION(2,2,kpk,kprofno) :: & 
    10821149         & zgmsk              ! Grid mask 
     1150#if defined key_bdy  
     1151      REAL(KIND=wp), DIMENSION(2,2,kprofno) :: & 
     1152         & zbmsk              ! Boundary mask 
     1153      REAL(KIND=wp), DIMENSION(jpi,jpj) :: zbdymask 
     1154#endif  
    10831155      REAL(KIND=wp), DIMENSION(2,2,kpk,kprofno) :: & 
    1084          & zgdepw 
     1156         & zgdepw          
    10851157      REAL(KIND=wp), DIMENSION(2,2,kprofno) :: & 
    10861158         & zglam, &           ! Model longitude at grid points 
    … …  
    11001172         ! For invalid points use 2,2 
    11011173 
    1102          IF ( kpobsqc(jobs) >= 10 ) THEN 
     1174         IF ( kpobsqc(jobs) >= kqc_cutoff ) THEN 
    11031175             
    11041176            igrdi(1,1,jobs) = 1 
    … …  
    11251197          
    11261198      END DO 
     1199 
     1200#if defined key_bdy  
     1201      ! Create a mask grid points in boundary rim 
     1202      IF (ld_bound_reject) THEN            
     1203         zbdymask(:,:) = 1.0_wp 
     1204         DO ji = 1, nb_bdy 
     1205            DO jj = 1, idx_bdy(ji)%nblen(1) 
     1206               zbdymask(idx_bdy(ji)%nbi(jj,1),idx_bdy(ji)%nbj(jj,1)) = 0.0_wp 
     1207            ENDDO 
     1208         ENDDO 
     1209      ENDIF 
     1210  
     1211      CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, zbdymask, zbmsk ) 
     1212#endif  
    11271213       
    11281214      CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, pmask, zgmsk ) 
    11291215      CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, plam, zglam ) 
    11301216      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(:,:,:), & 
    1134         &                     zgdepw ) 
    1135       ENDIF 
     1217      ! Need to know the bathy depth for each observation for sco 
     1218      CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, gdepw_n(:,:,:), & 
     1219         &                  zgdepw ) 
    11361220 
    11371221      DO jobs = 1, kprofno 
    11381222 
    11391223         ! Skip bad profiles 
    1140          IF ( kpobsqc(jobs) >= 10 ) CYCLE 
     1224         IF ( kpobsqc(jobs) >= kqc_cutoff ) CYCLE 
    11411225 
    11421226         ! Check if this observation is on a grid point 
    … …  
    11491233               IF ( ( ABS( zgphi(ji,jj,jobs) - pobsphi(jobs) ) < 1.0e-6_wp ) & 
    11501234                  & .AND. & 
    1151                   & ( ABS( zglam(ji,jj,jobs) - pobslam(jobs) ) < 1.0e-6_wp ) & 
     1235                  & ( ABS( MOD( zglam(ji,jj,jobs) - pobslam(jobs),360.0) ) < 1.0e-6_wp ) & 
    11521236                  & ) THEN 
    11531237                  lgridobs = .TRUE. 
    … …  
    11581242         END DO 
    11591243 
    1160          ! Check if next to land 
    1161          IF (  ANY( zgmsk(1:2,1:2,1,jobs) == 0.0_wp ) ) THEN 
    1162             ll_next_to_land=.TRUE. 
    1163          ELSE 
    1164             ll_next_to_land=.FALSE. 
     1244         ! Check if next to land  
     1245         IF (  ANY( zgmsk(1:2,1:2,1,jobs) == 0.0_wp ) ) THEN  
     1246            ll_next_to_land=.TRUE.  
     1247         ELSE  
     1248            ll_next_to_land=.FALSE.  
    11651249         ENDIF  
    1166  
     1250          
    11671251         ! Reject observations 
    11681252 
    … …  
    11761260               &  .OR. ( pobsdep(jobsp) < 0.0          )       & 
    11771261               &  .OR. ( pobsdep(jobsp) > gdepw_1d(kpk)) ) THEN 
    1178                kobsqc(jobsp) = kobsqc(jobsp) + 11 
     1262               kobsqc(jobsp) = IBSET(kobsqc(jobsp),11) 
    11791263               kosdobs = kosdobs + 1 
    11801264               CYCLE 
    11811265            ENDIF 
    11821266 
    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 
     1267            ! To check if an observations falls within land there are two cases:  
     1268            ! 1: z-coordibnates, where the check uses the mask  
     1269            ! 2: terrain following (eg s-coordinates),   
     1270            !    where we use the depth of the bottom cell to mask observations  
    11871271              
    1188             IF( ln_zps .OR. ln_zco ) THEN !(CASE 1) 
    1189                 
    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 
     1272            IF( (ln_linssh) .AND. ( ln_zps .OR. ln_zco )  ) THEN !(CASE 1)  
     1273                 
     1274               ! Flag if the observation falls with a model land cell  
     1275               IF ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) &  
     1276                  &  == 0.0_wp ) THEN  
     1277                  kobsqc(jobsp) = IBSET(kobsqc(jobsp),10) 
     1278                  klanobs = klanobs + 1  
     1279                  CYCLE  
     1280               ENDIF  
     1281              
     1282               ! Flag if the observation is close to land  
     1283               IF ( MINVAL( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) == &  
     1284                  &  0.0_wp) THEN  
     1285                  knlaobs = knlaobs + 1  
     1286                  IF (ld_nea) THEN    
     1287                     kobsqc(jobsp) = IBSET(kobsqc(jobsp),10) 
     1288                  ENDIF   
     1289               ENDIF  
     1290              
     1291            ELSE ! Case 2  
     1292               ! Flag if the observation is deeper than the bathymetry  
     1293               ! Or if it is within the mask  
     1294               IF ( ANY( zgdepw(1:2,1:2,kpk,jobs) < pobsdep(jobsp) ) & 
     1295                  &     .OR. &  
     1296                  &  ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) & 
     1297                  &  == 0.0_wp) ) THEN 
     1298                  kobsqc(jobsp) = IBSET(kobsqc(jobsp),10) 
     1299                  klanobs = klanobs + 1  
     1300                  CYCLE  
     1301               ENDIF  
     1302                 
     1303               ! Flag if the observation is close to land  
     1304               IF ( ll_next_to_land ) THEN  
     1305                  knlaobs = knlaobs + 1  
     1306                  IF (ld_nea) THEN    
     1307                     kobsqc(jobsp) = IBSET(kobsqc(jobsp),10) 
     1308                  ENDIF   
     1309               ENDIF  
     1310              
     1311            ENDIF 
     1312 
     1313            ! For observations on the grid reject them if their are at 
     1314            ! a masked point 
     1315             
     1316            IF (lgridobs) THEN 
     1317               IF (zgmsk(iig,ijg,kobsk(jobsp)-1,jobs) == 0.0_wp ) THEN 
     1318                  kobsqc(jobsp) = IBSET(kobsqc(jobs),10) 
    11941319                  klanobs = klanobs + 1 
    11951320                  CYCLE 
    11961321               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 
    1227             ENDIF 
    1228              
    1229             ! For observations on the grid reject them if their are at 
    1230             ! a masked point 
    1231              
    1232             IF (lgridobs) THEN 
    1233                IF (zgmsk(iig,ijg,kobsk(jobsp)-1,jobs) == 0.0_wp ) THEN 
    1234                   kobsqc(jobsp) = kobsqc(jobsp) + 12 
    1235                   klanobs = klanobs + 1 
    1236                   CYCLE 
    1237                ENDIF 
    1238             ENDIF 
    1239              
    1240             ! Flag if the observation falls is close to land 
    1241             IF ( MINVAL( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) == & 
    1242                &  0.0_wp) THEN 
    1243                IF (ld_nea) kobsqc(jobsp) = kobsqc(jobsp) + 14 
    1244                knlaobs = knlaobs + 1 
    1245             ENDIF 
    1246  
     1322            ENDIF 
     1323             
    12471324            ! Set observation depth equal to that of the first model depth 
    12481325            IF ( pobsdep(jobsp) <= pdep(1) ) THEN 
    … …  
    12501327            ENDIF 
    12511328             
     1329#if defined key_bdy 
     1330            ! Flag if the observation falls close to the boundary rim 
     1331            IF (ld_bound_reject) THEN 
     1332               IF ( MINVAL( zbmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN 
     1333                  kobsqc(jobsp) = IBSET(kobsqc(jobs),8) 
     1334                  kbdyobs = kbdyobs + 1 
     1335                  CYCLE 
     1336               ENDIF 
     1337               ! for observations on the grid... 
     1338               IF (lgridobs) THEN 
     1339                  IF (zbmsk(iig,ijg,jobs) == 0.0_wp ) THEN 
     1340                     kobsqc(jobsp) = IBSET(kobsqc(jobs),8) 
     1341                     kbdyobs = kbdyobs + 1 
     1342                     CYCLE 
     1343                  ENDIF 
     1344               ENDIF 
     1345            ENDIF 
     1346#endif  
     1347             
    12521348         END DO 
    12531349      END DO 
    … …  
    12551351   END SUBROUTINE obs_coo_spc_3d 
    12561352 
    1257    SUBROUTINE obs_pro_rej( profdata ) 
     1353   SUBROUTINE obs_pro_rej( profdata, kqc_cutoff ) 
    12581354      !!---------------------------------------------------------------------- 
    12591355      !!                    ***  ROUTINE obs_pro_rej *** 
    … …  
    12731369      !! * Arguments 
    12741370      TYPE(obs_prof), INTENT(INOUT) :: profdata     ! Profile data 
     1371      INTEGER, INTENT(IN) :: kqc_cutoff             ! QC cutoff value 
     1372 
    12751373      !! * Local declarations 
    12761374      INTEGER :: jprof 
    … …  
    12821380      DO jprof = 1, profdata%nprof 
    12831381 
    1284          IF ( profdata%nqc(jprof) > 10 ) THEN 
     1382         IF ( profdata%nqc(jprof) > kqc_cutoff ) THEN 
    12851383             
    12861384            DO jvar = 1, profdata%nvar 
    … …  
    12901388                   
    12911389                  profdata%var(jvar)%nvqc(jobs) = & 
    1292                      & profdata%var(jvar)%nvqc(jobs) + 26 
     1390                     & IBSET(profdata%var(jvar)%nvqc(jobs),14) 
    12931391 
    12941392               END DO 
    … …  
    13021400   END SUBROUTINE obs_pro_rej 
    13031401 
    1304    SUBROUTINE obs_uv_rej( profdata, knumu, knumv ) 
     1402   SUBROUTINE obs_uv_rej( profdata, knumu, knumv, kqc_cutoff ) 
    13051403      !!---------------------------------------------------------------------- 
    13061404      !!                    ***  ROUTINE obs_uv_rej *** 
    … …  
    13221420      INTEGER, INTENT(INOUT) :: knumu             ! Number of u rejected 
    13231421      INTEGER, INTENT(INOUT) :: knumv             ! Number of v rejected 
     1422      INTEGER, INTENT(IN) :: kqc_cutoff           ! QC cutoff value 
     1423 
    13241424      !! * Local declarations 
    13251425      INTEGER :: jprof 
    … …  
    13411441         DO jobs = profdata%npvsta(jprof,1), profdata%npvend(jprof,1) 
    13421442             
    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 
     1443            IF ( ( profdata%var(1)%nvqc(jobs) >  kqc_cutoff ) .AND. & 
     1444               & ( profdata%var(2)%nvqc(jobs) <=  kqc_cutoff) ) THEN 
     1445               profdata%var(2)%nvqc(jobs) = IBSET(profdata%var(1)%nvqc(jobs),15) 
    13461446               knumv = knumv + 1 
    13471447            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 
     1448            IF ( ( profdata%var(2)%nvqc(jobs) >  kqc_cutoff ) .AND. & 
     1449               & ( profdata%var(1)%nvqc(jobs) <=  kqc_cutoff) ) THEN 
     1450               profdata%var(1)%nvqc(jobs) = IBSET(profdata%var(1)%nvqc(jobs),15) 
    13511451               knumu = knumu + 1 
    13521452            ENDIF 

  • NEMO/branches/UKMO/r8395_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_profiles_def.F90

    r11350 r11361  
    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 

  • NEMO/branches/UKMO/r8395_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_altbias.F90

    r11350 r11361  
    128128         ! Get the Alt bias data 
    129129          
    130          CALL iom_get( numaltbias, jpdom_data, 'altbias', z_altbias(:,:), 1 ) 
     130         CALL iom_get( numaltbias, jpdom_autoglo, 'altbias', z_altbias(:,:), 1 ) 
    131131          
    132132         ! Close the file 

  • NEMO/branches/UKMO/r8395_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_prof.F90

    r11350 r11361  
    4545   SUBROUTINE obs_rea_prof( profdata, knumfiles, cdfilenames, & 
    4646      &                     kvars, kextr, kstp, ddobsini, ddobsend, & 
    47       &                     ldvar1, ldvar2, ldignmis, ldsatt, & 
     47      &                     ldvar, ldignmis, ldsatt, & 
    4848      &                     ldmod, kdailyavtypes ) 
    4949      !!--------------------------------------------------------------------- 
    … …  
    7474      INTEGER, INTENT(IN) :: kextr      ! Number of extra fields for each var 
    7575      INTEGER, INTENT(IN) :: kstp       ! Ocean time-step index 
    76       LOGICAL, INTENT(IN) :: ldvar1     ! Observed variables switches 
    77       LOGICAL, INTENT(IN) :: ldvar2 
     76      LOGICAL, DIMENSION(kvars), INTENT(IN) :: ldvar     ! Observed variables switches 
    7877      LOGICAL, INTENT(IN) :: ldignmis   ! Ignore missing files 
    7978      LOGICAL, INTENT(IN) :: ldsatt     ! Compute salinity at all temperature points 
    … …  
    8786      CHARACTER(LEN=15), PARAMETER :: cpname='obs_rea_prof' 
    8887      CHARACTER(len=8) :: clrefdate 
    89       CHARACTER(len=6), DIMENSION(:), ALLOCATABLE :: clvars 
     88      CHARACTER(len=8), DIMENSION(:), ALLOCATABLE :: clvars 
    9089      INTEGER :: jvar 
    9190      INTEGER :: ji 
    … …  
    105104      INTEGER :: iprof 
    106105      INTEGER :: iproftot 
    107       INTEGER :: ivar1t0 
    108       INTEGER :: ivar2t0 
    109       INTEGER :: ivar1t 
    110       INTEGER :: ivar2t 
     106      INTEGER, DIMENSION(kvars) :: ivart0 
     107      INTEGER, DIMENSION(kvars) :: ivart 
    111108      INTEGER :: ip3dt 
    112109      INTEGER :: ios 
    113110      INTEGER :: ioserrcount 
    114       INTEGER :: ivar1tmpp 
    115       INTEGER :: ivar2tmpp 
     111      INTEGER, DIMENSION(kvars) :: ivartmpp 
    116112      INTEGER :: ip3dtmpp 
    117113      INTEGER :: itype 
    118114      INTEGER, DIMENSION(knumfiles) :: & 
    119115         & irefdate 
    120       INTEGER, DIMENSION(ntyp1770+1) :: & 
    121          & itypvar1,    & 
    122          & itypvar1mpp, & 
    123          & itypvar2,    & 
    124          & itypvar2mpp  
     116      INTEGER, DIMENSION(ntyp1770+1,kvars) :: & 
     117         & itypvar,    & 
     118         & itypvarmpp 
     119      INTEGER, DIMENSION(:,:), ALLOCATABLE :: & 
     120         & iobsi,    & 
     121         & iobsj,    & 
     122         & iproc 
    125123      INTEGER, DIMENSION(:), ALLOCATABLE :: & 
    126          & iobsi1,    & 
    127          & iobsj1,    & 
    128          & iproc1,    & 
    129          & iobsi2,    & 
    130          & iobsj2,    & 
    131          & iproc2,    & 
    132124         & iindx,    & 
    133125         & ifileidx, & 
    … …  
    147139      LOGICAL :: llvalprof 
    148140      LOGICAL :: lldavtimset 
     141      LOGICAL :: llcycle 
    149142      TYPE(obfbdata), POINTER, DIMENSION(:) :: & 
    150143         & inpfiles 
    … …  
    152145      ! Local initialization 
    153146      iprof = 0 
    154       ivar1t0 = 0 
    155       ivar2t0 = 0 
     147      ivart0(:) = 0 
    156148      ip3dt = 0 
    157149 
    … …  
    219211               &                ldgrid = .TRUE. ) 
    220212 
    221             IF ( inpfiles(jj)%nvar < 2 ) THEN 
     213            IF ( inpfiles(jj)%nvar /= kvars ) THEN 
    222214               CALL ctl_stop( 'Feedback format error: ', & 
    223                   &           ' less than 2 vars in profile file' ) 
     215                  &           ' unexpected number of vars in profile file' ) 
    224216            ENDIF 
    225217 
    … …  
    307299            inowin = 0 
    308300            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 
     301               IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE 
     302               llcycle = .TRUE. 
     303               DO jvar = 1, kvars 
     304                  IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN 
     305                     llcycle = .FALSE. 
     306                     EXIT 
     307                  ENDIF 
     308               END DO 
     309               IF ( llcycle ) CYCLE 
    312310               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. & 
    313311                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN 
    … …  
    317315            ALLOCATE( zlam(inowin)  ) 
    318316            ALLOCATE( zphi(inowin)  ) 
    319             ALLOCATE( iobsi1(inowin) ) 
    320             ALLOCATE( iobsj1(inowin) ) 
    321             ALLOCATE( iproc1(inowin) ) 
    322             ALLOCATE( iobsi2(inowin) ) 
    323             ALLOCATE( iobsj2(inowin) ) 
    324             ALLOCATE( iproc2(inowin) ) 
     317            ALLOCATE( iobsi(inowin,kvars) ) 
     318            ALLOCATE( iobsj(inowin,kvars) ) 
     319            ALLOCATE( iproc(inowin,kvars) ) 
    325320            inowin = 0 
    326321            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 
     322               IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE 
     323               llcycle = .TRUE. 
     324               DO jvar = 1, kvars 
     325                  IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN 
     326                     llcycle = .FALSE. 
     327                     EXIT 
     328                  ENDIF 
     329               END DO 
     330               IF ( llcycle ) CYCLE 
    330331               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. & 
    331332                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN 
    … …  
    336337            END DO 
    337338 
    338             IF ( TRIM( inpfiles(jj)%cname(1) ) == 'POTM' ) THEN 
    339                CALL obs_grid_search( inowin, zlam, zphi, iobsi1, iobsj1, & 
    340                   &                  iproc1, 'T' ) 
    341                iobsi2(:) = iobsi1(:) 
    342                iobsj2(:) = iobsj1(:) 
    343                iproc2(:) = iproc1(:) 
    344             ELSEIF ( TRIM( inpfiles(jj)%cname(1) ) == 'UVEL' ) THEN 
    345                CALL obs_grid_search( inowin, zlam, zphi, iobsi1, iobsj1, & 
    346                   &                  iproc1, 'U' ) 
    347                CALL obs_grid_search( inowin, zlam, zphi, iobsi2, iobsj2, & 
    348                   &                  iproc2, 'V' ) 
     339            ! Assume anything other than velocity is on T grid 
     340            IF ( TRIM( inpfiles(jj)%cname(1) ) == 'UVEL' ) THEN 
     341               CALL obs_grid_search( inowin, zlam, zphi, iobsi(:,1), iobsj(:,1), & 
     342                  &                  iproc(:,1), 'U' ) 
     343               CALL obs_grid_search( inowin, zlam, zphi, iobsi(:,2), iobsj(:,2), & 
     344                  &                  iproc(:,2), 'V' ) 
     345            ELSE 
     346               CALL obs_grid_search( inowin, zlam, zphi, iobsi(:,1), iobsj(:,1), & 
     347                  &                  iproc(:,1), 'T' ) 
     348               IF ( kvars > 1 ) THEN 
     349                  DO jvar = 2, kvars 
     350                     iobsi(:,jvar) = iobsi(:,1) 
     351                     iobsj(:,jvar) = iobsj(:,1) 
     352                     iproc(:,jvar) = iproc(:,1) 
     353                  END DO 
     354               ENDIF 
    349355            ENDIF 
    350356 
    351357            inowin = 0 
    352358            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 
     359               IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE 
     360               llcycle = .TRUE. 
     361               DO jvar = 1, kvars 
     362                  IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN 
     363                     llcycle = .FALSE. 
     364                     EXIT 
     365                  ENDIF 
     366               END DO 
     367               IF ( llcycle ) CYCLE 
    356368               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. & 
    357369                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN 
    358370                  inowin = inowin + 1 
    359                   inpfiles(jj)%iproc(ji,1) = iproc1(inowin) 
    360                   inpfiles(jj)%iobsi(ji,1) = iobsi1(inowin) 
    361                   inpfiles(jj)%iobsj(ji,1) = iobsj1(inowin) 
    362                   inpfiles(jj)%iproc(ji,2) = iproc2(inowin) 
    363                   inpfiles(jj)%iobsi(ji,2) = iobsi2(inowin) 
    364                   inpfiles(jj)%iobsj(ji,2) = iobsj2(inowin) 
    365                   IF ( inpfiles(jj)%iproc(ji,1) /= & 
    366                      & inpfiles(jj)%iproc(ji,2) ) THEN 
    367                      CALL ctl_stop( 'Error in obs_read_prof:', & 
    368                         & 'var1 and var2 observation on different processors') 
     371                  DO jvar = 1, kvars 
     372                     inpfiles(jj)%iproc(ji,jvar) = iproc(inowin,jvar) 
     373                     inpfiles(jj)%iobsi(ji,jvar) = iobsi(inowin,jvar) 
     374                     inpfiles(jj)%iobsj(ji,jvar) = iobsj(inowin,jvar) 
     375                  END DO 
     376                  IF ( kvars > 1 ) THEN 
     377                     DO jvar = 2, kvars 
     378                        IF ( inpfiles(jj)%iproc(ji,jvar) /= & 
     379                           & inpfiles(jj)%iproc(ji,1) ) THEN 
     380                           CALL ctl_stop( 'Error in obs_read_prof:', & 
     381                              & 'observation on different processors for different vars') 
     382                        ENDIF 
     383                     END DO 
    369384                  ENDIF 
    370385               ENDIF 
    371386            END DO 
    372             DEALLOCATE( zlam, zphi, iobsi1, iobsj1, iproc1, iobsi2, iobsj2, iproc2 ) 
     387            DEALLOCATE( zlam, zphi, iobsi, iobsj, iproc ) 
    373388 
    374389            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 
     390               IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE 
     391               llcycle = .TRUE. 
     392               DO jvar = 1, kvars 
     393                  IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN 
     394                     llcycle = .FALSE. 
     395                     EXIT 
     396                  ENDIF 
     397               END DO 
     398               IF ( llcycle ) CYCLE 
    378399               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. & 
    379400                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN 
    … …  
    384405                  ENDIF 
    385406                  llvalprof = .FALSE. 
    386                   IF ( ldvar1 ) THEN 
    387                      loop_t_count : DO ij = 1,inpfiles(jj)%nlev 
    388                         IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) & 
    389                            & CYCLE 
    390                         IF ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. & 
    391                            & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN 
    392                            ivar1t0 = ivar1t0 + 1 
    393                         ENDIF 
    394                      END DO loop_t_count 
    395                   ENDIF 
    396                   IF ( ldvar2 ) THEN 
    397                      loop_s_count : DO ij = 1,inpfiles(jj)%nlev 
    398                         IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) & 
    399                            & CYCLE 
    400                         IF ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. & 
    401                            & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN 
    402                            ivar2t0 = ivar2t0 + 1 
    403                         ENDIF 
    404                      END DO loop_s_count 
    405                   ENDIF 
    406                   loop_p_count : DO ij = 1,inpfiles(jj)%nlev 
     407                  DO jvar = 1, kvars 
     408                     IF ( ldvar(jvar) ) THEN 
     409                        DO ij = 1,inpfiles(jj)%nlev 
     410                           IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) & 
     411                              & CYCLE 
     412                           IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,jvar),2) .AND. & 
     413                              & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN 
     414                              ivart0(jvar) = ivart0(jvar) + 1 
     415                           ENDIF 
     416                        END DO 
     417                     ENDIF 
     418                  END DO 
     419                  DO ij = 1,inpfiles(jj)%nlev 
    407420                     IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) & 
    408421                        & 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. & 
    414                         &     ldvar2 ) ) THEN 
    415                         ip3dt = ip3dt + 1 
    416                         llvalprof = .TRUE. 
    417                      ENDIF 
    418                   END DO loop_p_count 
     422                     DO jvar = 1, kvars 
     423                        IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,jvar),2) .AND. & 
     424                           &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. & 
     425                           &    ldvar(jvar) ) ) THEN 
     426                           ip3dt = ip3dt + 1 
     427                           llvalprof = .TRUE. 
     428                           EXIT 
     429                        ENDIF 
     430                     END DO 
     431                  END DO 
    419432 
    420433                  IF ( llvalprof ) iprof = iprof + 1 
    … …  
    437450      DO jj = 1, inobf 
    438451         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 
     452            IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE 
     453            llcycle = .TRUE. 
     454            DO jvar = 1, kvars 
     455               IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN 
     456                  llcycle = .FALSE. 
     457                  EXIT 
     458               ENDIF 
     459            END DO 
     460            IF ( llcycle ) CYCLE 
    442461            IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. & 
    443462               & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN 
    … …  
    452471      DO jj = 1, inobf 
    453472         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 
     473            IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE 
     474            llcycle = .TRUE. 
     475            DO jvar = 1, kvars 
     476               IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN 
     477                  llcycle = .FALSE. 
     478                  EXIT 
     479               ENDIF 
     480            END DO 
     481            IF ( llcycle ) CYCLE 
    457482            IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. & 
    458483               & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN 
    … …  
    470495      iv3dt(:) = -1 
    471496      IF (ldsatt) THEN 
    472          iv3dt(1) = ip3dt 
    473          iv3dt(2) = ip3dt 
     497         iv3dt(:) = ip3dt 
    474498      ELSE 
    475          iv3dt(1) = ivar1t0 
    476          iv3dt(2) = ivar2t0 
     499         iv3dt(:) = ivart0(:) 
    477500      ENDIF 
    478501      CALL obs_prof_alloc( profdata, kvars, kextr, iprof, iv3dt, & 
    … …  
    487510 
    488511      ip3dt = 0 
    489       ivar1t = 0 
    490       ivar2t = 0 
    491       itypvar1   (:) = 0 
    492       itypvar1mpp(:) = 0 
    493  
    494       itypvar2   (:) = 0 
    495       itypvar2mpp(:) = 0 
     512      ivart(:) = 0 
     513      itypvar   (:,:) = 0 
     514      itypvarmpp(:,:) = 0 
    496515 
    497516      ioserrcount = 0 
    … …  
    501520         ji = iprofidx(iindx(jk)) 
    502521 
    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 
     522         IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE 
     523         llcycle = .TRUE. 
     524         DO jvar = 1, kvars 
     525            IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN 
     526               llcycle = .FALSE. 
     527               EXIT 
     528            ENDIF 
     529         END DO 
     530         IF ( llcycle ) CYCLE 
    506531 
    507532         IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND.  & 
    … …  
    518543            IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE 
    519544 
    520             IF ( ( inpfiles(jj)%ivqc(ji,1) > 2 ) .AND. & 
    521                & ( inpfiles(jj)%ivqc(ji,2) > 2 ) ) CYCLE 
     545            IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE 
     546            llcycle = .TRUE. 
     547            DO jvar = 1, kvars 
     548               IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN 
     549                  llcycle = .FALSE. 
     550                  EXIT 
     551               ENDIF 
     552            END DO 
     553            IF ( llcycle ) CYCLE 
    522554 
    523555            loop_prof : DO ij = 1, inpfiles(jj)%nlev 
    … …  
    526558                  & CYCLE 
    527559 
    528                IF ( ( inpfiles(jj)%ivlqc(ij,ji,1) <= 2 ) .AND. & 
    529                   & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN 
    530  
    531                   llvalprof = .TRUE.  
    532                   EXIT loop_prof 
    533  
    534                ENDIF 
    535  
    536                IF ( ( inpfiles(jj)%ivlqc(ij,ji,2) <= 2 ) .AND. & 
    537                   & ( inpfiles(jj)%idqc(ij,ji) <= 2 ) ) THEN 
    538  
    539                   llvalprof = .TRUE.  
    540                   EXIT loop_prof 
    541  
    542                ENDIF 
     560               DO jvar = 1, kvars 
     561                  IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,jvar),2) .AND. & 
     562                     & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN 
     563 
     564                     llvalprof = .TRUE.  
     565                     EXIT loop_prof 
     566 
     567                  ENDIF 
     568               END DO 
    543569 
    544570            END DO loop_prof 
    … …  
    572598 
    573599               ! Coordinate search parameters 
    574                profdata%mi  (iprof,1) = inpfiles(jj)%iobsi(ji,1) 
    575                profdata%mj  (iprof,1) = inpfiles(jj)%iobsj(ji,1) 
    576                profdata%mi  (iprof,2) = inpfiles(jj)%iobsi(ji,2) 
    577                profdata%mj  (iprof,2) = inpfiles(jj)%iobsj(ji,2) 
     600               DO jvar = 1, kvars 
     601                  profdata%mi  (iprof,jvar) = inpfiles(jj)%iobsi(ji,jvar) 
     602                  profdata%mj  (iprof,jvar) = inpfiles(jj)%iobsj(ji,jvar) 
     603               END DO 
    578604 
    579605               ! Profile WMO number 
    … …  
    615641                  IF (ldsatt) THEN 
    616642 
    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 
    623                         ip3dt = ip3dt + 1 
    624                      ELSE 
    625                         CYCLE 
     643                     DO jvar = 1, kvars 
     644                        IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,jvar),2) .AND. & 
     645                           &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. & 
     646                           &    ldvar(jvar) ) ) THEN 
     647                           ip3dt = ip3dt + 1 
     648                           EXIT 
     649                        ELSE IF ( jvar == kvars ) THEN 
     650                           CYCLE loop_p 
     651                        ENDIF 
     652                     END DO 
     653 
     654                  ENDIF 
     655 
     656                  DO jvar = 1, kvars 
     657                   
     658                     IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,jvar),2) .AND. & 
     659                       &   .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. & 
     660                       &    ldvar(jvar) ) .OR. ldsatt ) THEN 
     661 
     662                        IF (ldsatt) THEN 
     663 
     664                           ivart(jvar) = ip3dt 
     665 
     666                        ELSE 
     667 
     668                           ivart(jvar) = ivart(jvar) + 1 
     669 
     670                        ENDIF 
     671 
     672                        ! Depth of jvar observation 
     673                        profdata%var(jvar)%vdep(ivart(jvar)) = & 
     674                           &                inpfiles(jj)%pdep(ij,ji) 
     675 
     676                        ! Depth of jvar observation QC 
     677                        profdata%var(jvar)%idqc(ivart(jvar)) = & 
     678                           &                inpfiles(jj)%idqc(ij,ji) 
     679 
     680                        ! Depth of jvar observation QC flags 
     681                        profdata%var(jvar)%idqcf(:,ivart(jvar)) = & 
     682                           &                inpfiles(jj)%idqcf(:,ij,ji) 
     683 
     684                        ! Profile index 
     685                        profdata%var(jvar)%nvpidx(ivart(jvar)) = iprof 
     686 
     687                        ! Vertical index in original profile 
     688                        profdata%var(jvar)%nvlidx(ivart(jvar)) = ij 
     689 
     690                        ! Profile jvar value 
     691                        IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,jvar),2) .AND. & 
     692                           & .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN 
     693                           profdata%var(jvar)%vobs(ivart(jvar)) = & 
     694                              &                inpfiles(jj)%pob(ij,ji,jvar) 
     695                           IF ( ldmod ) THEN 
     696                              profdata%var(jvar)%vmod(ivart(jvar)) = & 
     697                                 &                inpfiles(jj)%padd(ij,ji,1,jvar) 
     698                           ENDIF 
     699                           ! Count number of profile var1 data as function of type 
     700                           itypvar( profdata%ntyp(iprof) + 1, jvar ) = & 
     701                              & itypvar( profdata%ntyp(iprof) + 1, jvar ) + 1 
     702                        ELSE 
     703                           profdata%var(jvar)%vobs(ivart(jvar)) = fbrmdi 
     704                        ENDIF 
     705 
     706                        ! Profile jvar qc 
     707                        profdata%var(jvar)%nvqc(ivart(jvar)) = & 
     708                           & inpfiles(jj)%ivlqc(ij,ji,jvar) 
     709 
     710                        ! Profile jvar qc flags 
     711                        profdata%var(jvar)%nvqcf(:,ivart(jvar)) = & 
     712                           & inpfiles(jj)%ivlqcf(:,ij,ji,jvar) 
     713 
     714                        ! Profile insitu T value 
     715                        IF ( TRIM( inpfiles(jj)%cname(jvar) ) == 'POTM' ) THEN 
     716                           profdata%var(jvar)%vext(ivart(jvar),1) = & 
     717                              &                inpfiles(jj)%pext(ij,ji,1) 
     718                        ENDIF