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

Changeset 11202


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

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

Location:
branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC
Files:
4 added
19 deleted
17 edited

Legend:

Unmodified
Added
Removed

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

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

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

    r4990 r11202  
    66   !!====================================================================== 
    77 
    8    !!---------------------------------------------------------------------- 
    9    !!   'key_diaobs' : Switch on the observation diagnostic computation 
    108   !!---------------------------------------------------------------------- 
    119   !!   dia_obs_init : Reading and prepare observations 
    … …  
    1513   !!   fin_date     : Compute the final date YYYYMMDD.HHMMSS 
    1614   !!---------------------------------------------------------------------- 
    17    !! * Modules used    
     15   !! * Modules used 
    1816   USE wrk_nemo                 ! Memory Allocation 
    1917   USE par_kind                 ! Precision variables 
    … …  
    2119   USE par_oce 
    2220   USE dom_oce                  ! Ocean space and time domain variables 
    23    USE obs_fbm, ONLY: ln_cl4    ! Class 4 diagnostic switch 
    24    USE obs_read_prof            ! Reading and allocation of observations (Coriolis) 
    25    USE obs_read_sla             ! Reading and allocation of SLA observations   
    26    USE obs_read_sst             ! Reading and allocation of SST observations   
     21   USE obs_read_prof            ! Reading and allocation of profile obs 
     22   USE obs_read_surf            ! Reading and allocation of surface obs 
    2723   USE obs_readmdt              ! Reading and allocation of MDT for SLA. 
    28    USE obs_read_seaice          ! Reading and allocation of Sea Ice observations   
    29    USE obs_read_vel             ! Reading and allocation of velocity component observations 
    3024   USE obs_prep                 ! Preparation of obs. (grid search etc). 
    3125   USE obs_oper                 ! Observation operators 
    … …  
    3327   USE obs_grid                 ! Grid searching 
    3428   USE obs_read_altbias         ! Bias treatment for altimeter 
     29   USE obs_sstbias              ! Bias correction routine for SST 
    3530   USE obs_profiles_def         ! Profile data definitions 
    36    USE obs_profiles             ! Profile data storage 
    3731   USE obs_surf_def             ! Surface data definitions 
    38    USE obs_sla                  ! SLA data storage 
    39    USE obs_sst                  ! SST data storage 
    40    USE obs_seaice               ! Sea Ice data storage 
    4132   USE obs_types                ! Definitions for observation types 
    4233   USE mpp_map                  ! MPP mapping 
    … …  
    5243      &   dia_obs_dealloc  ! Deallocate dia_obs data 
    5344 
    54    !! * Shared Module variables 
    55    LOGICAL, PUBLIC, PARAMETER :: & 
    56 #if defined key_diaobs 
    57       & lk_diaobs = .TRUE.   !: Logical switch for observation diangostics 
    58 #else 
    59       & lk_diaobs = .FALSE.  !: Logical switch for observation diangostics 
    60 #endif 
    61  
    6245   !! * Module variables 
    63    LOGICAL, PUBLIC :: ln_t3d         !: Logical switch for temperature profiles 
    64    LOGICAL, PUBLIC :: ln_s3d         !: Logical switch for salinity profiles 
    65    LOGICAL, PUBLIC :: ln_ena         !: Logical switch for the ENACT data set 
    66    LOGICAL, PUBLIC :: ln_cor         !: Logical switch for the Coriolis data set 
    67    LOGICAL, PUBLIC :: ln_profb       !: Logical switch for profile feedback datafiles 
    68    LOGICAL, PUBLIC :: ln_sla         !: Logical switch for sea level anomalies  
    69    LOGICAL, PUBLIC :: ln_sladt       !: Logical switch for SLA from AVISO files 
    70    LOGICAL, PUBLIC :: ln_slafb       !: Logical switch for SLA from feedback files 
    71    LOGICAL, PUBLIC :: ln_sst         !: Logical switch for sea surface temperature 
    72    LOGICAL, PUBLIC :: ln_reysst      !: Logical switch for Reynolds sea surface temperature 
    73    LOGICAL, PUBLIC :: ln_ghrsst      !: Logical switch for GHRSST data 
    74    LOGICAL, PUBLIC :: ln_sstfb       !: Logical switch for SST from feedback files 
    75    LOGICAL, PUBLIC :: ln_seaice      !: Logical switch for sea ice concentration 
    76    LOGICAL, PUBLIC :: ln_vel3d       !: Logical switch for velocity component (u,v) observations 
    77    LOGICAL, PUBLIC :: ln_velavcur    !: Logical switch for raw daily averaged netCDF current meter vel. data  
    78    LOGICAL, PUBLIC :: ln_velhrcur    !: Logical switch for raw high freq netCDF current meter vel. data  
    79    LOGICAL, PUBLIC :: ln_velavadcp   !: Logical switch for raw daily averaged netCDF ADCP vel. data  
    80    LOGICAL, PUBLIC :: ln_velhradcp   !: Logical switch for raw high freq netCDF ADCP vel. data  
    81    LOGICAL, PUBLIC :: ln_velfb       !: Logical switch for velocities from feedback files 
    82    LOGICAL, PUBLIC :: ln_ssh         !: Logical switch for sea surface height 
    83    LOGICAL, PUBLIC :: ln_sss         !: Logical switch for sea surface salinity 
    84    LOGICAL, PUBLIC :: ln_sstnight    !: Logical switch for night mean SST observations 
    85    LOGICAL, PUBLIC :: ln_nea         !: Remove observations near land 
    86    LOGICAL, PUBLIC :: ln_altbias     !: Logical switch for altimeter bias   
    87    LOGICAL, PUBLIC :: ln_ignmis      !: Logical switch for ignoring missing files 
    88    LOGICAL, PUBLIC :: ln_s_at_t      !: Logical switch to compute model S at T observations 
    89  
    90    REAL(KIND=dp), PUBLIC :: dobsini   !: Observation window start date YYYYMMDD.HHMMSS 
    91    REAL(KIND=dp), PUBLIC :: dobsend   !: Observation window end date YYYYMMDD.HHMMSS 
    92    
    93    INTEGER, PUBLIC :: n1dint       !: Vertical interpolation method 
    94    INTEGER, PUBLIC :: n2dint       !: Horizontal interpolation method  
    95  
     46   LOGICAL, PUBLIC :: & 
     47      &       lk_diaobs = .TRUE.   !: Include this for backwards compatibility at NEMO 3.6. 
     48   LOGICAL :: ln_diaobs            !: Logical switch for the obs operator 
     49   LOGICAL :: ln_sstnight          !: Logical switch for night mean SST obs 
     50   LOGICAL :: ln_default_fp_indegs !: T=> Default obs footprint size specified in degrees, F=> in metres 
     51   LOGICAL :: ln_sla_fp_indegs     !: T=>     SLA obs footprint size specified in degrees, F=> in metres 
     52   LOGICAL :: ln_sst_fp_indegs     !: T=>     SST obs footprint size specified in degrees, F=> in metres 
     53   LOGICAL :: ln_sss_fp_indegs     !: T=>     SSS obs footprint size specified in degrees, F=> in metres 
     54   LOGICAL :: ln_sic_fp_indegs     !: T=> sea-ice obs footprint size specified in degrees, F=> in metres 
     55 
     56   REAL(wp) :: rn_default_avglamscl !: Default E/W diameter of observation footprint 
     57   REAL(wp) :: rn_default_avgphiscl !: Default N/S diameter of observation footprint 
     58   REAL(wp) :: rn_sla_avglamscl     !: E/W diameter of SLA observation footprint 
     59   REAL(wp) :: rn_sla_avgphiscl     !: N/S diameter of SLA observation footprint 
     60   REAL(wp) :: rn_sst_avglamscl     !: E/W diameter of SST observation footprint 
     61   REAL(wp) :: rn_sst_avgphiscl     !: N/S diameter of SST observation footprint 
     62   REAL(wp) :: rn_sss_avglamscl     !: E/W diameter of SSS observation footprint 
     63   REAL(wp) :: rn_sss_avgphiscl     !: N/S diameter of SSS observation footprint 
     64   REAL(wp) :: rn_sic_avglamscl     !: E/W diameter of sea-ice observation footprint 
     65   REAL(wp) :: rn_sic_avgphiscl     !: N/S diameter of sea-ice observation footprint 
     66 
     67   INTEGER :: nn_1dint         !: Vertical interpolation method 
     68   INTEGER :: nn_2dint_default !: Default horizontal interpolation method 
     69   INTEGER :: nn_2dint_sla     !: SLA horizontal interpolation method (-1 = default) 
     70   INTEGER :: nn_2dint_sst     !: SST horizontal interpolation method (-1 = default) 
     71   INTEGER :: nn_2dint_sss     !: SSS horizontal interpolation method (-1 = default) 
     72   INTEGER :: nn_2dint_sic     !: Seaice horizontal interpolation method (-1 = default) 
     73  
    9674   INTEGER, DIMENSION(imaxavtypes) :: & 
    97       & endailyavtypes !: ENACT data types which are daily average 
    98  
    99    INTEGER, PARAMETER :: MaxNumFiles = 1000 
    100    LOGICAL, DIMENSION(MaxNumFiles) :: & 
    101       & ln_profb_ena, & !: Is the feedback files from ENACT data ? 
    102    !                    !: If so use endailyavtypes 
    103       & ln_profb_enatim !: Change tim for 820 enact data set. 
    104     
    105    LOGICAL, DIMENSION(MaxNumFiles) :: & 
    106       & ln_velfb_av   !: Is the velocity feedback files daily average? 
     75      & nn_profdavtypes      !: Profile data types representing a daily average 
     76   INTEGER :: nproftypes     !: Number of profile obs types 
     77   INTEGER :: nsurftypes     !: Number of surface obs types 
     78   INTEGER, DIMENSION(:), ALLOCATABLE :: & 
     79      & nvarsprof, &         !: Number of profile variables 
     80      & nvarssurf            !: Number of surface variables 
     81   INTEGER, DIMENSION(:), ALLOCATABLE :: & 
     82      & nextrprof, &         !: Number of profile extra variables 
     83      & nextrsurf            !: Number of surface extra variables 
     84   INTEGER, DIMENSION(:), ALLOCATABLE :: & 
     85      & n2dintsurf           !: Interpolation option for surface variables 
     86   REAL(wp), DIMENSION(:), ALLOCATABLE :: & 
     87      & ravglamscl, &        !: E/W diameter of averaging footprint for surface variables 
     88      & ravgphiscl           !: N/S diameter of averaging footprint for surface variables 
    10789   LOGICAL, DIMENSION(:), ALLOCATABLE :: & 
    108       & ld_enact     !: Profile data is ENACT so use endailyavtypes 
    109    LOGICAL, DIMENSION(:), ALLOCATABLE :: & 
    110       & ld_velav     !: Velocity data is daily averaged 
    111    LOGICAL, DIMENSION(:), ALLOCATABLE :: & 
    112       & ld_sstnight  !: SST observation corresponds to night mean 
     90      & lfpindegs, &         !: T=> surface obs footprint size specified in degrees, F=> in metres 
     91      & llnightav            !: Logical for calculating night-time averages 
     92 
     93   TYPE(obs_surf), PUBLIC, POINTER, DIMENSION(:) :: & 
     94      & surfdata, &          !: Initial surface data 
     95      & surfdataqc           !: Surface data after quality control 
     96   TYPE(obs_prof), PUBLIC, POINTER, DIMENSION(:) :: & 
     97      & profdata, &          !: Initial profile data 
     98      & profdataqc           !: Profile data after quality control 
     99 
     100   CHARACTER(len=8), PUBLIC, DIMENSION(:), ALLOCATABLE :: & 
     101      & cobstypesprof, &     !: Profile obs types 
     102      & cobstypessurf        !: Surface obs types 
    113103 
    114104   !!---------------------------------------------------------------------- 
    … …  
    118108   !!---------------------------------------------------------------------- 
    119109 
     110   !! * Substitutions  
     111#  include "domzgr_substitute.h90" 
    120112CONTAINS 
    121113 
    … …  
    135127      !!        !  06-10  (A. Weaver) Cleaning and add controls 
    136128      !!        !  07-03  (K. Mogensen) General handling of profiles 
     129      !!        !  14-08  (J.While) Incorporated SST bias correction 
     130      !!        !  15-02  (M. Martin) Simplification of namelist and code 
    137131      !!---------------------------------------------------------------------- 
    138132 
    … …  
    140134 
    141135      !! * Local declarations 
    142       CHARACTER(len=128) :: enactfiles(MaxNumFiles) 
    143       CHARACTER(len=128) :: coriofiles(MaxNumFiles) 
    144       CHARACTER(len=128) :: profbfiles(MaxNumFiles) 
    145       CHARACTER(len=128) :: sstfiles(MaxNumFiles)       
    146       CHARACTER(len=128) :: sstfbfiles(MaxNumFiles)  
    147       CHARACTER(len=128) :: slafilesact(MaxNumFiles)       
    148       CHARACTER(len=128) :: slafilespas(MaxNumFiles)       
    149       CHARACTER(len=128) :: slafbfiles(MaxNumFiles) 
    150       CHARACTER(len=128) :: seaicefiles(MaxNumFiles)            
    151       CHARACTER(len=128) :: velcurfiles(MaxNumFiles)   
    152       CHARACTER(len=128) :: veladcpfiles(MaxNumFiles)     
    153       CHARACTER(len=128) :: velavcurfiles(MaxNumFiles) 
    154       CHARACTER(len=128) :: velhrcurfiles(MaxNumFiles) 
    155       CHARACTER(len=128) :: velavadcpfiles(MaxNumFiles) 
    156       CHARACTER(len=128) :: velhradcpfiles(MaxNumFiles) 
    157       CHARACTER(len=128) :: velfbfiles(MaxNumFiles) 
    158       CHARACTER(LEN=128) :: reysstname 
    159       CHARACTER(LEN=12)  :: reysstfmt 
    160       CHARACTER(LEN=128) :: bias_file 
    161       CHARACTER(LEN=20)  :: datestr=" ", timestr=" " 
    162       NAMELIST/namobs/ln_ena, ln_cor, ln_profb, ln_t3d, ln_s3d,       & 
    163          &            ln_sla, ln_sladt, ln_slafb,                     & 
    164          &            ln_ssh, ln_sst, ln_sstfb, ln_sss, ln_nea,       & 
    165          &            enactfiles, coriofiles, profbfiles,             & 
    166          &            slafilesact, slafilespas, slafbfiles,           & 
    167          &            sstfiles, sstfbfiles,                           & 
    168          &            ln_seaice, seaicefiles,                         & 
    169          &            dobsini, dobsend, n1dint, n2dint,               & 
    170          &            nmsshc, mdtcorr, mdtcutoff,                     & 
    171          &            ln_reysst, ln_ghrsst, reysstname, reysstfmt,    & 
    172          &            ln_sstnight,                                    & 
    173          &            ln_grid_search_lookup,                          & 
    174          &            grid_search_file, grid_search_res,              & 
    175          &            ln_grid_global, bias_file, ln_altbias,          & 
    176          &            endailyavtypes, ln_s_at_t, ln_profb_ena,        & 
    177          &            ln_vel3d, ln_velavcur, velavcurfiles,           & 
    178          &            ln_velhrcur, velhrcurfiles,                     & 
    179          &            ln_velavadcp, velavadcpfiles,                   & 
    180          &            ln_velhradcp, velhradcpfiles,                   & 
    181          &            ln_velfb, velfbfiles, ln_velfb_av,              & 
    182          &            ln_profb_enatim, ln_ignmis, ln_cl4 
    183  
    184       INTEGER :: jprofset 
    185       INTEGER :: jveloset 
    186       INTEGER :: jvar 
    187       INTEGER :: jnumenact 
    188       INTEGER :: jnumcorio 
    189       INTEGER :: jnumprofb 
    190       INTEGER :: jnumslaact 
    191       INTEGER :: jnumslapas 
    192       INTEGER :: jnumslafb 
    193       INTEGER :: jnumsst 
    194       INTEGER :: jnumsstfb 
    195       INTEGER :: jnumseaice 
    196       INTEGER :: jnumvelavcur 
    197       INTEGER :: jnumvelhrcur   
    198       INTEGER :: jnumvelavadcp 
    199       INTEGER :: jnumvelhradcp    
    200       INTEGER :: jnumvelfb 
    201       INTEGER :: ji 
    202       INTEGER :: jset 
    203       INTEGER :: ios                 ! Local integer output status for namelist read 
    204       LOGICAL :: lmask(MaxNumFiles), ll_u3d, ll_v3d 
     136      INTEGER, PARAMETER :: & 
     137         & jpmaxnfiles = 1000    ! Maximum number of files for each obs type 
     138      INTEGER, DIMENSION(:), ALLOCATABLE :: & 
     139         & ifilesprof, &         ! Number of profile files 
     140         & ifilessurf            ! Number of surface files 
     141      INTEGER :: ios             ! Local integer output status for namelist read 
     142      INTEGER :: jtype           ! Counter for obs types 
     143      INTEGER :: jvar            ! Counter for variables 
     144      INTEGER :: jfile           ! Counter for files 
     145      INTEGER :: jnumsstbias     ! Number of SST bias files to read and apply 
     146      INTEGER :: n2dint_type     ! Local version of nn_2dint* 
     147 
     148      CHARACTER(len=128), DIMENSION(jpmaxnfiles) :: & 
     149         & cn_profbfiles,      & ! T/S profile input filenames 
     150         & cn_sstfbfiles,      & ! Sea surface temperature input filenames 
     151         & cn_slafbfiles,      & ! Sea level anomaly input filenames 
     152         & cn_sicfbfiles,      & ! Seaice concentration input filenames 
     153         & cn_velfbfiles,      & ! Velocity profile input filenames 
     154         & cn_sssfbfiles,      & ! Sea surface salinity input filenames 
     155         & cn_slchltotfbfiles, & ! Surface total              log10(chlorophyll) input filenames 
     156         & cn_slchldiafbfiles, & ! Surface diatom             log10(chlorophyll) input filenames 
     157         & cn_slchlnonfbfiles, & ! Surface non-diatom         log10(chlorophyll) input filenames 
     158         & cn_slchldinfbfiles, & ! Surface dinoflagellate     log10(chlorophyll) input filenames 
     159         & cn_slchlmicfbfiles, & ! Surface microphytoplankton log10(chlorophyll) input filenames 
     160         & cn_slchlnanfbfiles, & ! Surface nanophytoplankton  log10(chlorophyll) input filenames 
     161         & cn_slchlpicfbfiles, & ! Surface picophytoplankton  log10(chlorophyll) input filenames 
     162         & cn_schltotfbfiles,  & ! Surface total              chlorophyll        input filenames 
     163         & cn_slphytotfbfiles, & ! Surface total      log10(phytoplankton carbon) input filenames 
     164         & cn_slphydiafbfiles, & ! Surface diatom     log10(phytoplankton carbon) input filenames 
     165         & cn_slphynonfbfiles, & ! Surface non-diatom log10(phytoplankton carbon) input filenames 
     166         & cn_sspmfbfiles,     & ! Surface suspended particulate matter input filenames 
     167         & cn_sfco2fbfiles,    & ! Surface fugacity         of carbon dioxide input filenames 
     168         & cn_spco2fbfiles,    & ! Surface partial pressure of carbon dioxide input filenames 
     169         & cn_plchltotfbfiles, & ! Profile total log10(chlorophyll) input filenames 
     170         & cn_pchltotfbfiles,  & ! Profile total chlorophyll input filenames 
     171         & cn_pno3fbfiles,     & ! Profile nitrate input filenames 
     172         & cn_psi4fbfiles,     & ! Profile silicate input filenames 
     173         & cn_ppo4fbfiles,     & ! Profile phosphate input filenames 
     174         & cn_pdicfbfiles,     & ! Profile dissolved inorganic carbon input filenames 
     175         & cn_palkfbfiles,     & ! Profile alkalinity input filenames 
     176         & cn_pphfbfiles,      & ! Profile pH input filenames 
     177         & cn_po2fbfiles,      & ! Profile dissolved oxygen input filenames 
     178         & cn_sstbiasfiles       ! SST bias input filenames 
     179 
     180      CHARACTER(LEN=128) :: & 
     181         & cn_altbiasfile        ! Altimeter bias input filename 
     182 
     183 
     184      LOGICAL :: ln_t3d          ! Logical switch for temperature profiles 
     185      LOGICAL :: ln_s3d          ! Logical switch for salinity profiles 
     186      LOGICAL :: ln_sla          ! Logical switch for sea level anomalies  
     187      LOGICAL :: ln_sst          ! Logical switch for sea surface temperature 
     188      LOGICAL :: ln_sic          ! Logical switch for sea ice concentration 
     189      LOGICAL :: ln_sss          ! Logical switch for sea surface salinity obs 
     190      LOGICAL :: ln_vel3d        ! Logical switch for velocity (u,v) obs 
     191      LOGICAL :: ln_slchltot     ! Logical switch for surface total              log10(chlorophyll) obs 
     192      LOGICAL :: ln_slchldia     ! Logical switch for surface diatom             log10(chlorophyll) obs 
     193      LOGICAL :: ln_slchlnon     ! Logical switch for surface non-diatom         log10(chlorophyll) obs 
     194      LOGICAL :: ln_slchldin     ! Logical switch for surface dinoflagellate     log10(chlorophyll) obs 
     195      LOGICAL :: ln_slchlmic     ! Logical switch for surface microphytoplankton log10(chlorophyll) obs 
     196      LOGICAL :: ln_slchlnan     ! Logical switch for surface nanophytoplankton  log10(chlorophyll) obs 
     197      LOGICAL :: ln_slchlpic     ! Logical switch for surface picophytoplankton  log10(chlorophyll) obs 
     198      LOGICAL :: ln_schltot      ! Logical switch for surface total              chlorophyll        obs 
     199      LOGICAL :: ln_slphytot     ! Logical switch for surface total      log10(phytoplankton carbon) obs 
     200      LOGICAL :: ln_slphydia     ! Logical switch for surface diatom     log10(phytoplankton carbon) obs 
     201      LOGICAL :: ln_slphynon     ! Logical switch for surface non-diatom log10(phytoplankton carbon) obs 
     202      LOGICAL :: ln_sspm         ! Logical switch for surface suspended particulate matter obs 
     203      LOGICAL :: ln_sfco2        ! Logical switch for surface fugacity         of carbon dioxide obs 
     204      LOGICAL :: ln_spco2        ! Logical switch for surface partial pressure of carbon dioxide obs 
     205      LOGICAL :: ln_plchltot     ! Logical switch for profile total log10(chlorophyll) obs 
     206      LOGICAL :: ln_pchltot      ! Logical switch for profile total chlorophyll obs 
     207      LOGICAL :: ln_pno3         ! Logical switch for profile nitrate obs 
     208      LOGICAL :: ln_psi4         ! Logical switch for profile silicate obs 
     209      LOGICAL :: ln_ppo4         ! Logical switch for profile phosphate obs 
     210      LOGICAL :: ln_pdic         ! Logical switch for profile dissolved inorganic carbon obs 
     211      LOGICAL :: ln_palk         ! Logical switch for profile alkalinity obs 
     212      LOGICAL :: ln_pph          ! Logical switch for profile pH obs 
     213      LOGICAL :: ln_po2          ! Logical switch for profile dissolved oxygen obs 
     214      LOGICAL :: ln_nea          ! Logical switch to remove obs near land 
     215      LOGICAL :: ln_altbias      ! Logical switch for altimeter bias 
     216      LOGICAL :: ln_sstbias      ! Logical switch for bias correction of SST 
     217      LOGICAL :: ln_ignmis       ! Logical switch for ignoring missing files 
     218      LOGICAL :: ln_s_at_t       ! Logical switch to compute model S at T obs 
     219      LOGICAL :: ln_bound_reject ! Logical switch for rejecting obs near the boundary 
     220 
     221      REAL(dp) :: rn_dobsini     ! Obs window start date YYYYMMDD.HHMMSS 
     222      REAL(dp) :: rn_dobsend     ! Obs window end date   YYYYMMDD.HHMMSS 
     223 
     224      REAL(wp) :: ztype_avglamscl ! Local version of rn_*_avglamscl 
     225      REAL(wp) :: ztype_avgphiscl ! Local version of rn_*_avgphiscl 
     226 
     227      CHARACTER(len=128), DIMENSION(:,:), ALLOCATABLE :: & 
     228         & clproffiles, &        ! Profile filenames 
     229         & clsurffiles           ! Surface filenames 
     230 
     231      LOGICAL, DIMENSION(:), ALLOCATABLE :: llvar   ! Logical for profile variable read 
     232      LOGICAL :: ltype_fp_indegs ! Local version of ln_*_fp_indegs 
     233      LOGICAL :: ltype_night     ! Local version of ln_sstnight (false for other variables) 
     234 
     235      REAL(wp), POINTER, DIMENSION(:,:,:) :: & 
     236         & zglam                 ! Model longitudes for profile variables 
     237      REAL(wp), POINTER, DIMENSION(:,:,:) :: & 
     238         & zgphi                 ! Model latitudes for profile variables 
     239      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: & 
     240         & zmask                 ! Model land/sea mask associated with variables 
     241 
     242 
     243      NAMELIST/namobs/ln_diaobs, ln_t3d, ln_s3d, ln_sla,              & 
     244         &            ln_sst, ln_sic, ln_sss, ln_vel3d,               & 
     245         &            ln_slchltot, ln_slchldia, ln_slchlnon,          & 
     246         &            ln_slchldin, ln_slchlmic, ln_slchlnan,          & 
     247         &            ln_slchlpic, ln_schltot,                        & 
     248         &            ln_slphytot, ln_slphydia, ln_slphynon,          & 
     249         &            ln_sspm,     ln_sfco2,    ln_spco2,             & 
     250         &            ln_plchltot, ln_pchltot,  ln_pno3,              & 
     251         &            ln_psi4,     ln_ppo4,     ln_pdic,              & 
     252         &            ln_palk,     ln_pph,      ln_po2,               & 
     253         &            ln_altbias, ln_sstbias, ln_nea,                 & 
     254         &            ln_grid_global, ln_grid_search_lookup,          & 
     255         &            ln_ignmis, ln_s_at_t, ln_bound_reject,          & 
     256         &            ln_sstnight, ln_default_fp_indegs,              & 
     257         &            ln_sla_fp_indegs, ln_sst_fp_indegs,             & 
     258         &            ln_sss_fp_indegs, ln_sic_fp_indegs,             & 
     259         &            cn_profbfiles, cn_slafbfiles,                   & 
     260         &            cn_sstfbfiles, cn_sicfbfiles,                   & 
     261         &            cn_velfbfiles, cn_sssfbfiles,                   & 
     262         &            cn_slchltotfbfiles, cn_slchldiafbfiles,         & 
     263         &            cn_slchlnonfbfiles, cn_slchldinfbfiles,         & 
     264         &            cn_slchlmicfbfiles, cn_slchlnanfbfiles,         & 
     265         &            cn_slchlpicfbfiles, cn_schltotfbfiles,          & 
     266         &            cn_slphytotfbfiles, cn_slphydiafbfiles,         & 
     267         &            cn_slphynonfbfiles, cn_sspmfbfiles,             & 
     268         &            cn_sfco2fbfiles, cn_spco2fbfiles,               & 
     269         &            cn_plchltotfbfiles, cn_pchltotfbfiles,          & 
     270         &            cn_pno3fbfiles, cn_psi4fbfiles, cn_ppo4fbfiles, & 
     271         &            cn_pdicfbfiles, cn_palkfbfiles, cn_pphfbfiles,  & 
     272         &            cn_po2fbfiles,                                  & 
     273         &            cn_sstbiasfiles, cn_altbiasfile,                & 
     274         &            cn_gridsearchfile, rn_gridsearchres,            & 
     275         &            rn_dobsini, rn_dobsend,                         & 
     276         &            rn_default_avglamscl, rn_default_avgphiscl,     & 
     277         &            rn_sla_avglamscl, rn_sla_avgphiscl,             & 
     278         &            rn_sst_avglamscl, rn_sst_avgphiscl,             & 
     279         &            rn_sss_avglamscl, rn_sss_avgphiscl,             & 
     280         &            rn_sic_avglamscl, rn_sic_avgphiscl,             & 
     281         &            nn_1dint, nn_2dint_default,                     & 
     282         &            nn_2dint_sla, nn_2dint_sst,                     & 
     283         &            nn_2dint_sss, nn_2dint_sic,                     & 
     284         &            nn_msshc, rn_mdtcorr, rn_mdtcutoff,             & 
     285         &            nn_profdavtypes 
    205286 
    206287      !----------------------------------------------------------------------- 
    … …  
    208289      !----------------------------------------------------------------------- 
    209290 
    210       enactfiles(:) = '' 
    211       coriofiles(:) = '' 
    212       profbfiles(:) = '' 
    213       slafilesact(:) = '' 
    214       slafilespas(:) = '' 
    215       slafbfiles(:) = '' 
    216       sstfiles(:)   = '' 
    217       sstfbfiles(:) = '' 
    218       seaicefiles(:) = '' 
    219       velcurfiles(:) = '' 
    220       veladcpfiles(:) = '' 
    221       velavcurfiles(:) = '' 
    222       velhrcurfiles(:) = '' 
    223       velavadcpfiles(:) = '' 
    224       velhradcpfiles(:) = '' 
    225       velfbfiles(:) = '' 
    226       velcurfiles(:) = '' 
    227       veladcpfiles(:) = '' 
    228       endailyavtypes(:) = -1 
    229       endailyavtypes(1) = 820 
    230       ln_profb_ena(:) = .FALSE. 
    231       ln_profb_enatim(:) = .TRUE. 
    232       ln_velfb_av(:) = .FALSE. 
    233       ln_ignmis = .FALSE. 
    234        
    235       CALL ini_date( dobsini ) 
    236       CALL fin_date( dobsend ) 
    237   
    238       ! Read Namelist namobs : control observation diagnostics 
    239       REWIND( numnam_ref )              ! Namelist namobs in reference namelist : Diagnostic: control observation 
     291      ! Some namelist arrays need initialising 
     292      cn_profbfiles(:)      = '' 
     293      cn_slafbfiles(:)      = '' 
     294      cn_sstfbfiles(:)      = '' 
     295      cn_sicfbfiles(:)      = '' 
     296      cn_velfbfiles(:)      = '' 
     297      cn_sssfbfiles(:)      = '' 
     298      cn_slchltotfbfiles(:) = '' 
     299      cn_slchldiafbfiles(:) = '' 
     300      cn_slchlnonfbfiles(:) = '' 
     301      cn_slchldinfbfiles(:) = '' 
     302      cn_slchlmicfbfiles(:) = '' 
     303      cn_slchlnanfbfiles(:) = '' 
     304      cn_slchlpicfbfiles(:) = '' 
     305      cn_schltotfbfiles(:)  = '' 
     306      cn_slphytotfbfiles(:) = '' 
     307      cn_slphydiafbfiles(:) = '' 
     308      cn_slphynonfbfiles(:) = '' 
     309      cn_sspmfbfiles(:)     = '' 
     310      cn_sfco2fbfiles(:)    = '' 
     311      cn_spco2fbfiles(:)    = '' 
     312      cn_plchltotfbfiles(:) = '' 
     313      cn_pchltotfbfiles(:)  = '' 
     314      cn_pno3fbfiles(:)     = '' 
     315      cn_psi4fbfiles(:)     = '' 
     316      cn_ppo4fbfiles(:)     = '' 
     317      cn_pdicfbfiles(:)     = '' 
     318      cn_palkfbfiles(:)     = '' 
     319      cn_pphfbfiles(:)      = '' 
     320      cn_po2fbfiles(:)      = '' 
     321      cn_sstbiasfiles(:)    = '' 
     322      nn_profdavtypes(:)    = -1 
     323 
     324      CALL ini_date( rn_dobsini ) 
     325      CALL fin_date( rn_dobsend ) 
     326 
     327      ! Read namelist namobs : control observation diagnostics 
     328      REWIND( numnam_ref )   ! Namelist namobs in reference namelist 
    240329      READ  ( numnam_ref, namobs, IOSTAT = ios, ERR = 901) 
    241330901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in reference namelist', lwp ) 
    242331 
    243       REWIND( numnam_cfg )              ! Namelist namobs in configuration namelist : Diagnostic: control observation 
     332      REWIND( numnam_cfg )   ! Namelist namobs in configuration namelist 
    244333      READ  ( numnam_cfg, namobs, IOSTAT = ios, ERR = 902 ) 
    245334902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in configuration namelist', lwp ) 
    246335      IF(lwm) WRITE ( numond, namobs ) 
    247336 
    248       ! Count number of files for each type 
    249       IF (ln_ena) THEN 
    250          lmask(:) = .FALSE. 
    251          WHERE (enactfiles(:) /= '') lmask(:) = .TRUE. 
    252          jnumenact = COUNT(lmask) 
     337      lk_diaobs = .FALSE. 
     338#if defined key_diaobs 
     339      IF ( ln_diaobs ) lk_diaobs = .TRUE. 
     340#endif 
     341 
     342      IF ( .NOT. lk_diaobs ) THEN 
     343         IF(lwp) WRITE(numout,cform_war) 
     344         IF(lwp) WRITE(numout,*)' ln_diaobs is set to false or key_diaobs is not set, so not calling dia_obs' 
     345         RETURN 
    253346      ENDIF 
    254       IF (ln_cor) THEN 
    255          lmask(:) = .FALSE. 
    256          WHERE (coriofiles(:) /= '') lmask(:) = .TRUE. 
    257          jnumcorio = COUNT(lmask) 
    258       ENDIF 
    259       IF (ln_profb) THEN 
    260          lmask(:) = .FALSE. 
    261          WHERE (profbfiles(:) /= '') lmask(:) = .TRUE. 
    262          jnumprofb = COUNT(lmask) 
    263       ENDIF 
    264       IF (ln_sladt) THEN 
    265          lmask(:) = .FALSE. 
    266          WHERE (slafilesact(:) /= '') lmask(:) = .TRUE. 
    267          jnumslaact = COUNT(lmask) 
    268          lmask(:) = .FALSE. 
    269          WHERE (slafilespas(:) /= '') lmask(:) = .TRUE. 
    270          jnumslapas = COUNT(lmask) 
    271       ENDIF 
    272       IF (ln_slafb) THEN 
    273          lmask(:) = .FALSE. 
    274          WHERE (slafbfiles(:) /= '') lmask(:) = .TRUE. 
    275          jnumslafb = COUNT(lmask) 
    276          lmask(:) = .FALSE. 
    277       ENDIF 
    278       IF (ln_ghrsst) THEN 
    279          lmask(:) = .FALSE. 
    280          WHERE (sstfiles(:) /= '') lmask(:) = .TRUE. 
    281          jnumsst = COUNT(lmask) 
    282       ENDIF       
    283       IF (ln_sstfb) THEN 
    284          lmask(:) = .FALSE. 
    285          WHERE (sstfbfiles(:) /= '') lmask(:) = .TRUE. 
    286          jnumsstfb = COUNT(lmask) 
    287          lmask(:) = .FALSE. 
    288       ENDIF 
    289       IF (ln_seaice) THEN 
    290          lmask(:) = .FALSE. 
    291          WHERE (seaicefiles(:) /= '') lmask(:) = .TRUE. 
    292          jnumseaice = COUNT(lmask) 
    293       ENDIF 
    294       IF (ln_velavcur) THEN 
    295          lmask(:) = .FALSE. 
    296          WHERE (velavcurfiles(:) /= '') lmask(:) = .TRUE. 
    297          jnumvelavcur = COUNT(lmask) 
    298       ENDIF 
    299       IF (ln_velhrcur) THEN 
    300          lmask(:) = .FALSE. 
    301          WHERE (velhrcurfiles(:) /= '') lmask(:) = .TRUE. 
    302          jnumvelhrcur = COUNT(lmask) 
    303       ENDIF 
    304       IF (ln_velavadcp) THEN 
    305          lmask(:) = .FALSE. 
    306          WHERE (velavadcpfiles(:) /= '') lmask(:) = .TRUE. 
    307          jnumvelavadcp = COUNT(lmask) 
    308       ENDIF 
    309       IF (ln_velhradcp) THEN 
    310          lmask(:) = .FALSE. 
    311          WHERE (velhradcpfiles(:) /= '') lmask(:) = .TRUE. 
    312          jnumvelhradcp = COUNT(lmask) 
    313       ENDIF 
    314       IF (ln_velfb) THEN 
    315          lmask(:) = .FALSE. 
    316          WHERE (velfbfiles(:) /= '') lmask(:) = .TRUE. 
    317          jnumvelfb = COUNT(lmask) 
    318          lmask(:) = .FALSE. 
    319       ENDIF 
    320        
    321       ! Control print 
     347 
    322348      IF(lwp) THEN 
    323349         WRITE(numout,*) 
    … …  
    325351         WRITE(numout,*) '~~~~~~~~~~~~' 
    326352         WRITE(numout,*) '          Namelist namobs : set observation diagnostic parameters'  
    327          WRITE(numout,*) '             Logical switch for T profile observations          ln_t3d = ', ln_t3d 
    328          WRITE(numout,*) '             Logical switch for S profile observations          ln_s3d = ', ln_s3d 
    329          WRITE(numout,*) '             Logical switch for ENACT insitu data set           ln_ena = ', ln_ena 
    330          WRITE(numout,*) '             Logical switch for Coriolis insitu data set        ln_cor = ', ln_cor 
    331          WRITE(numout,*) '             Logical switch for feedback insitu data set      ln_profb = ', ln_profb 
    332          WRITE(numout,*) '             Logical switch for SLA observations                ln_sla = ', ln_sla 
    333          WRITE(numout,*) '             Logical switch for AVISO SLA data                ln_sladt = ', ln_sladt 
    334          WRITE(numout,*) '             Logical switch for feedback SLA data             ln_slafb = ', ln_slafb 
    335          WRITE(numout,*) '             Logical switch for SSH observations                ln_ssh = ', ln_ssh 
    336          WRITE(numout,*) '             Logical switch for SST observations                ln_sst = ', ln_sst 
    337          WRITE(numout,*) '             Logical switch for Reynolds observations        ln_reysst = ', ln_reysst     
    338          WRITE(numout,*) '             Logical switch for GHRSST observations          ln_ghrsst = ', ln_ghrsst 
    339          WRITE(numout,*) '             Logical switch for feedback SST data             ln_sstfb = ', ln_sstfb 
    340          WRITE(numout,*) '             Logical switch for night-time SST obs         ln_sstnight = ', ln_sstnight 
    341          WRITE(numout,*) '             Logical switch for SSS observations                ln_sss = ', ln_sss 
    342          WRITE(numout,*) '             Logical switch for Sea Ice observations         ln_seaice = ', ln_seaice 
    343          WRITE(numout,*) '             Logical switch for velocity observations         ln_vel3d = ', ln_vel3d 
    344          WRITE(numout,*) '             Logical switch for velocity daily av. cur.    ln_velavcur = ', ln_velavcur 
    345          WRITE(numout,*) '             Logical switch for velocity high freq. cur.   ln_velhrcur = ', ln_velhrcur 
    346          WRITE(numout,*) '             Logical switch for velocity daily av. ADCP   ln_velavadcp = ', ln_velavadcp 
    347          WRITE(numout,*) '             Logical switch for velocity high freq. ADCP  ln_velhradcp = ', ln_velhradcp 
    348          WRITE(numout,*) '             Logical switch for feedback velocity data        ln_velfb = ', ln_velfb 
    349          WRITE(numout,*) '             Global distribtion of observations         ln_grid_global = ',ln_grid_global 
    350          WRITE(numout,*) & 
    351    '             Logical switch for obs grid search w/lookup table  ln_grid_search_lookup = ',ln_grid_search_lookup 
     353         WRITE(numout,*) '             Logical switch for T profile observations                ln_t3d = ', ln_t3d 
     354         WRITE(numout,*) '             Logical switch for S profile observations                ln_s3d = ', ln_s3d 
     355         WRITE(numout,*) '             Logical switch for SLA observations                      ln_sla = ', ln_sla 
     356         WRITE(numout,*) '             Logical switch for SST observations                      ln_sst = ', ln_sst 
     357         WRITE(numout,*) '             Logical switch for Sea Ice observations                  ln_sic = ', ln_sic 
     358         WRITE(numout,*) '             Logical switch for velocity observations               ln_vel3d = ', ln_vel3d 
     359         WRITE(numout,*) '             Logical switch for SSS observations                      ln_sss = ', ln_sss 
     360         WRITE(numout,*) '             Logical switch for surface total logchl obs         ln_slchltot = ', ln_slchltot 
     361         WRITE(numout,*) '             Logical switch for surface diatom logchl obs        ln_slchldia = ', ln_slchldia 
     362         WRITE(numout,*) '             Logical switch for surface non-diatom logchl obs    ln_slchlnon = ', ln_slchlnon 
     363         WRITE(numout,*) '             Logical switch for surface dino logchl obs          ln_slchldin = ', ln_slchldin 
     364         WRITE(numout,*) '             Logical switch for surface micro logchl obs         ln_slchlmic = ', ln_slchlmic 
     365         WRITE(numout,*) '             Logical switch for surface nano logchl obs          ln_slchlnan = ', ln_slchlnan 
     366         WRITE(numout,*) '             Logical switch for surface pico logchl obs          ln_slchlpic = ', ln_slchlpic 
     367         WRITE(numout,*) '             Logical switch for surface total chl obs             ln_schltot = ', ln_schltot 
     368         WRITE(numout,*) '             Logical switch for surface total log(phyC) obs      ln_slphytot = ', ln_slphytot 
     369         WRITE(numout,*) '             Logical switch for surface diatom log(phyC) obs     ln_slphydia = ', ln_slphydia 
     370         WRITE(numout,*) '             Logical switch for surface non-diatom log(phyC) obs ln_slphynon = ', ln_slphynon 
     371         WRITE(numout,*) '             Logical switch for surface SPM observations             ln_sspm = ', ln_sspm 
     372         WRITE(numout,*) '             Logical switch for surface fCO2 observations           ln_sfco2 = ', ln_sfco2 
     373         WRITE(numout,*) '             Logical switch for surface pCO2 observations           ln_spco2 = ', ln_spco2 
     374         WRITE(numout,*) '             Logical switch for profile total logchl obs         ln_plchltot = ', ln_plchltot 
     375         WRITE(numout,*) '             Logical switch for profile total chl obs             ln_pchltot = ', ln_pchltot 
     376         WRITE(numout,*) '             Logical switch for profile nitrate obs                  ln_pno3 = ', ln_pno3 
     377         WRITE(numout,*) '             Logical switch for profile silicate obs                 ln_psi4 = ', ln_psi4 
     378         WRITE(numout,*) '             Logical switch for profile phosphate obs                ln_ppo4 = ', ln_ppo4 
     379         WRITE(numout,*) '             Logical switch for profile DIC obs                      ln_pdic = ', ln_pdic 
     380         WRITE(numout,*) '             Logical switch for profile alkalinity obs               ln_palk = ', ln_palk 
     381         WRITE(numout,*) '             Logical switch for profile pH obs                        ln_pph = ', ln_pph 
     382         WRITE(numout,*) '             Logical switch for profile oxygen obs                    ln_po2 = ', ln_po2 
     383         WRITE(numout,*) '             Global distribution of observations              ln_grid_global = ', ln_grid_global 
     384         WRITE(numout,*) '             Logical switch for obs grid search lookup ln_grid_search_lookup = ', ln_grid_search_lookup 
    352385         IF (ln_grid_search_lookup) & 
    353             WRITE(numout,*) '             Grid search lookup file header       grid_search_file = ', grid_search_file 
    354          IF (ln_ena) THEN 
    355             DO ji = 1, jnumenact 
    356                WRITE(numout,'(1X,2A)') '             ENACT input observation file name          enactfiles = ', & 
    357                   TRIM(enactfiles(ji)) 
    358             END DO 
    359          ENDIF 
    360          IF (ln_cor) THEN 
    361             DO ji = 1, jnumcorio 
    362                WRITE(numout,'(1X,2A)') '             Coriolis input observation file name       coriofiles = ', & 
    363                   TRIM(coriofiles(ji)) 
    364             END DO 
    365          ENDIF 
    366          IF (ln_profb) THEN 
    367             DO ji = 1, jnumprofb 
    368                IF (ln_profb_ena(ji)) THEN 
    369                   WRITE(numout,'(1X,2A)') '       Enact feedback input observation file name       profbfiles = ', & 
    370                      TRIM(profbfiles(ji)) 
     386            WRITE(numout,*) '             Grid search lookup file header                cn_gridsearchfile = ', cn_gridsearchfile 
     387         WRITE(numout,*) '             Initial date in window YYYYMMDD.HHMMSS               rn_dobsini = ', rn_dobsini 
     388         WRITE(numout,*) '             Final date in window YYYYMMDD.HHMMSS                 rn_dobsend = ', rn_dobsend 
     389         WRITE(numout,*) '             Type of vertical interpolation method                  nn_1dint = ', nn_1dint 
     390         WRITE(numout,*) '             Default horizontal interpolation method        nn_2dint_default = ', nn_2dint_default 
     391         WRITE(numout,*) '             Type of horizontal interpolation method for SLA    nn_2dint_sla = ', nn_2dint_sla 
     392         WRITE(numout,*) '             Type of horizontal interpolation method for SST    nn_2dint_sst = ', nn_2dint_sst 
     393         WRITE(numout,*) '             Type of horizontal interpolation method for SSS    nn_2dint_sss = ', nn_2dint_sss 
     394         WRITE(numout,*) '             Type of horizontal interpolation method for SIC    nn_2dint_sic = ', nn_2dint_sic 
     395         WRITE(numout,*) '             Default E/W diameter of obs footprint      rn_default_avglamscl = ', rn_default_avglamscl 
     396         WRITE(numout,*) '             Default N/S diameter of obs footprint      rn_default_avgphiscl = ', rn_default_avgphiscl 
     397         WRITE(numout,*) '             Default obs footprint in deg [T] or m [F]  ln_default_fp_indegs = ', ln_default_fp_indegs 
     398         WRITE(numout,*) '             SLA E/W diameter of obs footprint              rn_sla_avglamscl = ', rn_sla_avglamscl 
     399         WRITE(numout,*) '             SLA N/S diameter of obs footprint              rn_sla_avgphiscl = ', rn_sla_avgphiscl 
     400         WRITE(numout,*) '             SLA obs footprint in deg [T] or m [F]          ln_sla_fp_indegs = ', ln_sla_fp_indegs 
     401         WRITE(numout,*) '             SST E/W diameter of obs footprint              rn_sst_avglamscl = ', rn_sst_avglamscl 
     402         WRITE(numout,*) '             SST N/S diameter of obs footprint              rn_sst_avgphiscl = ', rn_sst_avgphiscl 
     403         WRITE(numout,*) '             SST obs footprint in deg [T] or m [F]          ln_sst_fp_indegs = ', ln_sst_fp_indegs 
     404         WRITE(numout,*) '             SIC E/W diameter of obs footprint              rn_sic_avglamscl = ', rn_sic_avglamscl 
     405         WRITE(numout,*) '             SIC N/S diameter of obs footprint              rn_sic_avgphiscl = ', rn_sic_avgphiscl 
     406         WRITE(numout,*) '             SIC obs footprint in deg [T] or m [F]          ln_sic_fp_indegs = ', ln_sic_fp_indegs 
     407         WRITE(numout,*) '             Rejection of observations near land switch               ln_nea = ', ln_nea 
     408         WRITE(numout,*) '             Rejection of obs near open bdys                 ln_bound_reject = ', ln_bound_reject 
     409         WRITE(numout,*) '             MSSH correction scheme                                 nn_msshc = ', nn_msshc 
     410         WRITE(numout,*) '             MDT  correction                                      rn_mdtcorr = ', rn_mdtcorr 
     411         WRITE(numout,*) '             MDT cutoff for computed correction                 rn_mdtcutoff = ', rn_mdtcutoff 
     412         WRITE(numout,*) '             Logical switch for alt bias                          ln_altbias = ', ln_altbias 
     413         WRITE(numout,*) '             Logical switch for sst bias                          ln_sstbias = ', ln_sstbias 
     414         WRITE(numout,*) '             Logical switch for ignoring missing files             ln_ignmis = ', ln_ignmis 
     415         WRITE(numout,*) '             Daily average types                             nn_profdavtypes = ', nn_profdavtypes 
     416         WRITE(numout,*) '             Logical switch for night-time SST obs               ln_sstnight = ', ln_sstnight 
     417      ENDIF 
     418      !----------------------------------------------------------------------- 
     419      ! Set up list of observation types to be used 
     420      ! and the files associated with each type 
     421      !----------------------------------------------------------------------- 
     422 
     423      nproftypes = COUNT( (/ln_t3d .OR. ln_s3d, ln_vel3d, ln_plchltot,          & 
     424         &                  ln_pchltot,  ln_pno3,     ln_psi4,     ln_ppo4,     & 
     425         &                  ln_pdic,     ln_palk,     ln_pph,      ln_po2 /) ) 
     426      nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic, ln_sss,                     & 
     427         &                  ln_slchltot, ln_slchldia, ln_slchlnon, ln_slchldin, & 
     428         &                  ln_slchlmic, ln_slchlnan, ln_slchlpic, ln_schltot,  & 
     429         &                  ln_slphytot, ln_slphydia, ln_slphynon, ln_sspm,     & 
     430         &                  ln_sfco2,    ln_spco2 /) ) 
     431 
     432      IF ( nproftypes == 0 .AND. nsurftypes == 0 ) THEN 
     433         IF(lwp) WRITE(numout,cform_war) 
     434         IF(lwp) WRITE(numout,*) ' ln_diaobs is set to true, but all obs operator logical flags', & 
     435            &                    ' are set to .FALSE. so turning off calls to dia_obs' 
     436         nwarn = nwarn + 1 
     437         lk_diaobs = .FALSE. 
     438         RETURN 
     439      ENDIF 
     440 
     441      IF(lwp) WRITE(numout,*) '          Number of profile obs types: ',nproftypes 
     442      IF ( nproftypes > 0 ) THEN 
     443 
     444         ALLOCATE( cobstypesprof(nproftypes) ) 
     445         ALLOCATE( ifilesprof(nproftypes) ) 
     446         ALLOCATE( clproffiles(nproftypes,jpmaxnfiles) ) 
     447 
     448         jtype = 0 
     449         IF (ln_t3d .OR. ln_s3d) THEN 
     450            jtype = jtype + 1 
     451            cobstypesprof(jtype) = 'prof' 
     452            clproffiles(jtype,:) = cn_profbfiles 
     453         ENDIF 
     454         IF (ln_vel3d) THEN 
     455            jtype = jtype + 1 
     456            cobstypesprof(jtype) =  'vel' 
     457            clproffiles(jtype,:) = cn_velfbfiles 
     458         ENDIF 
     459         IF (ln_plchltot) THEN 
     460            jtype = jtype + 1 
     461            cobstypesprof(jtype) = 'plchltot' 
     462            clproffiles(jtype,:) = cn_plchltotfbfiles 
     463         ENDIF 
     464         IF (ln_pchltot) THEN 
     465            jtype = jtype + 1 
     466            cobstypesprof(jtype) = 'pchltot' 
     467            clproffiles(jtype,:) = cn_pchltotfbfiles 
     468         ENDIF 
     469         IF (ln_pno3) THEN 
     470            jtype = jtype + 1 
     471            cobstypesprof(jtype) = 'pno3' 
     472            clproffiles(jtype,:) = cn_pno3fbfiles 
     473         ENDIF 
     474         IF (ln_psi4) THEN 
     475            jtype = jtype + 1 
     476            cobstypesprof(jtype) = 'psi4' 
     477            clproffiles(jtype,:) = cn_psi4fbfiles 
     478         ENDIF 
     479         IF (ln_ppo4) THEN 
     480            jtype = jtype + 1 
     481            cobstypesprof(jtype) = 'ppo4' 
     482            clproffiles(jtype,:) = cn_ppo4fbfiles 
     483         ENDIF 
     484         IF (ln_pdic) THEN 
     485            jtype = jtype + 1 
     486            cobstypesprof(jtype) = 'pdic' 
     487            clproffiles(jtype,:) = cn_pdicfbfiles 
     488         ENDIF 
     489         IF (ln_palk) THEN 
     490            jtype = jtype + 1 
     491            cobstypesprof(jtype) = 'palk' 
     492            clproffiles(jtype,:) = cn_palkfbfiles 
     493         ENDIF 
     494         IF (ln_pph) THEN 
     495            jtype = jtype + 1 
     496            cobstypesprof(jtype) = 'pph' 
     497            clproffiles(jtype,:) = cn_pphfbfiles 
     498         ENDIF 
     499         IF (ln_po2) THEN 
     500            jtype = jtype + 1 
     501            cobstypesprof(jtype) = 'po2' 
     502            clproffiles(jtype,:) = cn_po2fbfiles 
     503         ENDIF 
     504 
     505         CALL obs_settypefiles( nproftypes, jpmaxnfiles, ifilesprof, cobstypesprof, clproffiles ) 
     506 
     507      ENDIF 
     508 
     509      IF(lwp) WRITE(numout,*)'          Number of surface obs types: ',nsurftypes 
     510      IF ( nsurftypes > 0 ) THEN 
     511 
     512         ALLOCATE( cobstypessurf(nsurftypes) ) 
     513         ALLOCATE( ifilessurf(nsurftypes) ) 
     514         ALLOCATE( clsurffiles(nsurftypes, jpmaxnfiles) ) 
     515         ALLOCATE(n2dintsurf(nsurftypes)) 
     516         ALLOCATE(ravglamscl(nsurftypes)) 
     517         ALLOCATE(ravgphiscl(nsurftypes)) 
     518         ALLOCATE(lfpindegs(nsurftypes)) 
     519         ALLOCATE(llnightav(nsurftypes)) 
     520 
     521         jtype = 0 
     522         IF (ln_sla) THEN 
     523            jtype = jtype + 1 
     524            cobstypessurf(jtype) = 'sla' 
     525            clsurffiles(jtype,:) = cn_slafbfiles 
     526         ENDIF 
     527         IF (ln_sst) THEN 
     528            jtype = jtype + 1 
     529            cobstypessurf(jtype) = 'sst' 
     530            clsurffiles(jtype,:) = cn_sstfbfiles 
     531         ENDIF 
     532         IF (ln_sic) THEN 
     533            jtype = jtype + 1 
     534            cobstypessurf(jtype) = 'sic' 
     535            clsurffiles(jtype,:) = cn_sicfbfiles 
     536         ENDIF 
     537         IF (ln_sss) THEN 
     538            jtype = jtype + 1 
     539            cobstypessurf(jtype) = 'sss' 
     540            clsurffiles(jtype,:) = cn_sssfbfiles 
     541         ENDIF 
     542         IF (ln_slchltot) THEN 
     543            jtype = jtype + 1 
     544            cobstypessurf(jtype) = 'slchltot' 
     545            clsurffiles(jtype,:) = cn_slchltotfbfiles 
     546         ENDIF 
     547         IF (ln_slchldia) THEN 
     548            jtype = jtype + 1 
     549            cobstypessurf(jtype) = 'slchldia' 
     550            clsurffiles(jtype,:) = cn_slchldiafbfiles 
     551         ENDIF 
     552         IF (ln_slchlnon) THEN 
     553            jtype = jtype + 1 
     554            cobstypessurf(jtype) = 'slchlnon' 
     555            clsurffiles(jtype,:) = cn_slchlnonfbfiles 
     556         ENDIF 
     557         IF (ln_slchldin) THEN 
     558            jtype = jtype + 1 
     559            cobstypessurf(jtype) = 'slchldin' 
     560            clsurffiles(jtype,:) = cn_slchldinfbfiles 
     561         ENDIF 
     562         IF (ln_slchlmic) THEN 
     563            jtype = jtype + 1 
     564            cobstypessurf(jtype) = 'slchlmic' 
     565            clsurffiles(jtype,:) = cn_slchlmicfbfiles 
     566         ENDIF 
     567         IF (ln_slchlnan) THEN 
     568            jtype = jtype + 1 
     569            cobstypessurf(jtype) = 'slchlnan' 
     570            clsurffiles(jtype,:) = cn_slchlnanfbfiles 
     571         ENDIF 
     572         IF (ln_slchlpic) THEN 
     573            jtype = jtype + 1 
     574            cobstypessurf(jtype) = 'slchlpic' 
     575            clsurffiles(jtype,:) = cn_slchlpicfbfiles 
     576         ENDIF 
     577         IF (ln_schltot) THEN 
     578            jtype = jtype + 1 
     579            cobstypessurf(jtype) = 'schltot' 
     580            clsurffiles(jtype,:) = cn_schltotfbfiles 
     581         ENDIF 
     582         IF (ln_slphytot) THEN 
     583            jtype = jtype + 1 
     584            cobstypessurf(jtype) = 'slphytot' 
     585            clsurffiles(jtype,:) = cn_slphytotfbfiles 
     586         ENDIF 
     587         IF (ln_slphydia) THEN 
     588            jtype = jtype + 1 
     589            cobstypessurf(jtype) = 'slphydia' 
     590            clsurffiles(jtype,:) = cn_slphydiafbfiles 
     591         ENDIF 
     592         IF (ln_slphynon) THEN 
     593            jtype = jtype + 1 
     594            cobstypessurf(jtype) = 'slphynon' 
     595            clsurffiles(jtype,:) = cn_slphynonfbfiles 
     596         ENDIF 
     597         IF (ln_sspm) THEN 
     598            jtype = jtype + 1 
     599            cobstypessurf(jtype) = 'sspm' 
     600            clsurffiles(jtype,:) = cn_sspmfbfiles 
     601         ENDIF 
     602         IF (ln_sfco2) THEN 
     603            jtype = jtype + 1 
     604            cobstypessurf(jtype) = 'sfco2' 
     605            clsurffiles(jtype,:) = cn_sfco2fbfiles 
     606         ENDIF 
     607         IF (ln_spco2) THEN 
     608            jtype = jtype + 1 
     609            cobstypessurf(jtype) = 'spco2' 
     610            clsurffiles(jtype,:) = cn_spco2fbfiles 
     611         ENDIF 
     612 
     613         CALL obs_settypefiles( nsurftypes, jpmaxnfiles, ifilessurf, cobstypessurf, clsurffiles ) 
     614 
     615         DO jtype = 1, nsurftypes 
     616 
     617            IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN 
     618               IF ( nn_2dint_sla == -1 ) THEN 
     619                  n2dint_type  = nn_2dint_default 
    371620               ELSE 
    372                   WRITE(numout,'(1X,2A)') '             Feedback input observation file name       profbfiles = ', & 
    373                      TRIM(profbfiles(ji)) 
     621                  n2dint_type  = nn_2dint_sla 
    374622               ENDIF 
    375                WRITE(numout,'(1X,2A)') '       Enact feedback input time setting switch    ln_profb_enatim = ', ln_profb_enatim(ji) 
    376             END DO 
    377          ENDIF 
    378          IF (ln_sladt) THEN 
    379             DO ji = 1, jnumslaact 
    380                WRITE(numout,'(1X,2A)') '             Active SLA input observation file name    slafilesact = ', & 
    381                   TRIM(slafilesact(ji)) 
    382             END DO 
    383             DO ji = 1, jnumslapas 
    384                WRITE(numout,'(1X,2A)') '             Passive SLA input observation file name   slafilespas = ', & 
    385                   TRIM(slafilespas(ji)) 
    386             END DO 
    387          ENDIF 
    388          IF (ln_slafb) THEN 
    389             DO ji = 1, jnumslafb 
    390                WRITE(numout,'(1X,2A)') '             Feedback SLA input observation file name   slafbfiles = ', & 
    391                   TRIM(slafbfiles(ji)) 
    392             END DO 
    393          ENDIF 
    394          IF (ln_ghrsst) THEN 
    395             DO ji = 1, jnumsst 
    396                WRITE(numout,'(1X,2A)') '             GHRSST input observation file name           sstfiles = ', & 
    397                   TRIM(sstfiles(ji)) 
    398             END DO 
    399          ENDIF 
    400          IF (ln_sstfb) THEN 
    401             DO ji = 1, jnumsstfb 
    402                WRITE(numout,'(1X,2A)') '             Feedback SST input observation file name   sstfbfiles = ', & 
    403                   TRIM(sstfbfiles(ji)) 
    404             END DO 
    405          ENDIF 
    406          IF (ln_seaice) THEN 
    407             DO ji = 1, jnumseaice 
    408                WRITE(numout,'(1X,2A)') '             Sea Ice input observation file name       seaicefiles = ', & 
    409                   TRIM(seaicefiles(ji)) 
    410             END DO 
    411          ENDIF 
    412          IF (ln_velavcur) THEN 
    413             DO ji = 1, jnumvelavcur 
    414                WRITE(numout,'(1X,2A)') '             Vel. cur. daily av. input file name     velavcurfiles = ', & 
    415                   TRIM(velavcurfiles(ji)) 
    416             END DO 
    417          ENDIF 
    418          IF (ln_velhrcur) THEN 
    419             DO ji = 1, jnumvelhrcur 
    420                WRITE(numout,'(1X,2A)') '             Vel. cur. high freq. input file name    velhvcurfiles = ', & 
    421                   TRIM(velhrcurfiles(ji)) 
    422             END DO 
    423          ENDIF 
    424          IF (ln_velavadcp) THEN 
    425             DO ji = 1, jnumvelavadcp 
    426                WRITE(numout,'(1X,2A)') '             Vel. ADCP daily av. input file name    velavadcpfiles = ', & 
    427                   TRIM(velavadcpfiles(ji)) 
    428             END DO 
    429          ENDIF 
    430          IF (ln_velhradcp) THEN 
    431             DO ji = 1, jnumvelhradcp 
    432                WRITE(numout,'(1X,2A)') '             Vel. ADCP high freq. input file name   velhvadcpfiles = ', & 
    433                   TRIM(velhradcpfiles(ji)) 
    434             END DO 
    435          ENDIF 
    436          IF (ln_velfb) THEN 
    437             DO ji = 1, jnumvelfb 
    438                IF (ln_velfb_av(ji)) THEN 
    439                   WRITE(numout,'(1X,2A)') '             Vel. feedback daily av. input file name    velfbfiles = ', & 
    440                      TRIM(velfbfiles(ji)) 
     623               ztype_avglamscl = rn_sla_avglamscl 
     624               ztype_avgphiscl = rn_sla_avgphiscl 
     625               ltype_fp_indegs = ln_sla_fp_indegs 
     626               ltype_night     = .FALSE. 
     627            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sst' ) THEN 
     628               IF ( nn_2dint_sst == -1 ) THEN 
     629                  n2dint_type  = nn_2dint_default 
    441630               ELSE 
    442                   WRITE(numout,'(1X,2A)') '             Vel. feedback input observation file name  velfbfiles = ', & 
    443                      TRIM(velfbfiles(ji)) 
     631                  n2dint_type  = nn_2dint_sst 
    444632               ENDIF 
    445             END DO 
    446          ENDIF 
    447          WRITE(numout,*) '             Initial date in window YYYYMMDD.HHMMSS        dobsini = ', dobsini 
    448          WRITE(numout,*) '             Final date in window YYYYMMDD.HHMMSS          dobsend = ', dobsend 
    449          WRITE(numout,*) '             Type of vertical interpolation method          n1dint = ', n1dint 
    450          WRITE(numout,*) '             Type of horizontal interpolation method        n2dint = ', n2dint 
    451          WRITE(numout,*) '             Rejection of observations near land swithch    ln_nea = ', ln_nea 
    452          WRITE(numout,*) '             MSSH correction scheme                         nmsshc = ', nmsshc 
    453          WRITE(numout,*) '             MDT  correction                               mdtcorr = ', mdtcorr 
    454          WRITE(numout,*) '             MDT cutoff for computed correction          mdtcutoff = ', mdtcutoff 
    455          WRITE(numout,*) '             Logical switch for alt bias                ln_altbias = ', ln_altbias 
    456          WRITE(numout,*) '             Logical switch for ignoring missing files   ln_ignmis = ', ln_ignmis 
    457          WRITE(numout,*) '             ENACT daily average types                             = ',endailyavtypes 
     633               ztype_avglamscl = rn_sst_avglamscl 
     634               ztype_avgphiscl = rn_sst_avgphiscl 
     635               ltype_fp_indegs = ln_sst_fp_indegs 
     636               ltype_night     = ln_sstnight 
     637            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sic' ) THEN 
     638               IF ( nn_2dint_sic == -1 ) THEN 
     639                  n2dint_type  = nn_2dint_default 
     640               ELSE 
     641                  n2dint_type  = nn_2dint_sic 
     642               ENDIF 
     643               ztype_avglamscl = rn_sic_avglamscl 
     644               ztype_avgphiscl = rn_sic_avgphiscl 
     645               ltype_fp_indegs = ln_sic_fp_indegs 
     646               ltype_night     = .FALSE. 
     647            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sss' ) THEN 
     648               IF ( nn_2dint_sss == -1 ) THEN 
     649                  n2dint_type  = nn_2dint_default 
     650               ELSE 
     651                  n2dint_type  = nn_2dint_sss 
     652               ENDIF 
     653               ztype_avglamscl = rn_sss_avglamscl 
     654               ztype_avgphiscl = rn_sss_avgphiscl 
     655               ltype_fp_indegs = ln_sss_fp_indegs 
     656               ltype_night     = .FALSE. 
     657            ELSE 
     658               n2dint_type     = nn_2dint_default 
     659               ztype_avglamscl = rn_default_avglamscl 
     660               ztype_avgphiscl = rn_default_avgphiscl 
     661               ltype_fp_indegs = ln_default_fp_indegs 
     662               ltype_night     = .FALSE. 
     663            ENDIF 
     664             
     665            CALL obs_setinterpopts( nsurftypes, jtype, TRIM(cobstypessurf(jtype)), & 
     666               &                    nn_2dint_default, n2dint_type,                 & 
     667               &                    ztype_avglamscl, ztype_avgphiscl,              & 
     668               &                    ltype_fp_indegs, ltype_night,                  & 
     669               &                    n2dintsurf, ravglamscl, ravgphiscl,            & 
     670               &                    lfpindegs, llnightav ) 
     671 
     672         END DO 
    458673 
    459674      ENDIF 
    460        
     675 
     676      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 
     677 
     678 
     679      !----------------------------------------------------------------------- 
     680      ! Obs operator parameter checking and initialisations 
     681      !----------------------------------------------------------------------- 
     682 
    461683      IF ( ln_vel3d .AND. ( .NOT. ln_grid_global ) ) THEN 
    462684         CALL ctl_stop( 'Velocity data only works with ln_grid_global=.true.' ) 
    … …  
    464686      ENDIF 
    465687 
    466       CALL obs_typ_init 
    467        
    468       CALL mppmap_init 
    469        
    470       ! Parameter control 
    471 #if defined key_diaobs 
    472       IF ( ( .NOT. ln_t3d ).AND.( .NOT. ln_s3d ).AND.( .NOT. ln_sla ).AND. & 
    473          & ( .NOT. ln_vel3d ).AND.                                         & 
    474          & ( .NOT. ln_ssh ).AND.( .NOT. ln_sst ).AND.( .NOT. ln_sss ).AND. & 
    475          & ( .NOT. ln_seaice ).AND.( .NOT. ln_vel3d ) ) THEN 
    476          IF(lwp) WRITE(numout,cform_war) 
    477          IF(lwp) WRITE(numout,*) ' key_diaobs is activated but logical flags', & 
    478             &                    ' ln_t3d, ln_s3d, ln_sla, ln_ssh, ln_sst, ln_sss, ln_seaice, ln_vel3d are all set to .FALSE.' 
    479          nwarn = nwarn + 1 
    480       ENDIF 
    481 #endif 
    482  
    483       CALL obs_grid_setup( ) 
    484       IF ( ( n1dint < 0 ).OR.( n1dint > 1 ) ) THEN 
     688      IF ( ( nn_1dint < 0 ) .OR. ( nn_1dint > 1 ) ) THEN 
    485689         CALL ctl_stop(' Choice of vertical (1D) interpolation method', & 
    486690            &                    ' is not available') 
    487691      ENDIF 
    488       IF ( ( n2dint < 0 ).OR.( n2dint > 4 ) ) THEN 
    489          CALL ctl_stop(' Choice of horizontal (2D) interpolation method', & 
     692 
     693      IF ( ( nn_2dint_default < 0 ) .OR. ( nn_2dint_default > 6 ) ) THEN 
     694         CALL ctl_stop(' Choice of default horizontal (2D) interpolation method', & 
    490695            &                    ' is not available') 
    491696      ENDIF 
    492697 
     698      CALL obs_typ_init 
     699 
     700      CALL mppmap_init 
     701 
     702      CALL obs_grid_setup( ) 
     703 
    493704      !----------------------------------------------------------------------- 
    494705      ! Depending on switches read the various observation types 
    495706      !----------------------------------------------------------------------- 
    496       !  - Temperature/salinity profiles 
    497  
    498       IF ( ln_t3d .OR. ln_s3d ) THEN 
    499  
    500          ! Set the number of variables for profiles to 2 (T and S) 
    501          nprofvars = 2 
    502          ! Set the number of extra variables for profiles to 1 (insitu temp). 
    503          nprofextr = 1 
    504  
    505          ! Count how may insitu data sets we have and allocate data. 
    506          jprofset = 0 
    507          IF ( ln_ena ) jprofset = jprofset + 1 
    508          IF ( ln_cor ) jprofset = jprofset + 1 
    509          IF ( ln_profb ) jprofset = jprofset + jnumprofb 
    510          nprofsets = jprofset 
    511          IF ( nprofsets > 0 ) THEN 
    512             ALLOCATE(ld_enact(nprofsets)) 
    513             ALLOCATE(profdata(nprofsets)) 
    514             ALLOCATE(prodatqc(nprofsets)) 
    515          ENDIF 
    516  
    517          jprofset = 0 
    518            
    519          ! ENACT insitu data 
    520  
    521          IF ( ln_ena ) THEN 
    522  
    523             jprofset = jprofset + 1 
     707 
     708      IF ( nproftypes > 0 ) THEN 
     709 
     710         ALLOCATE(profdata(nproftypes)) 
     711         ALLOCATE(profdataqc(nproftypes)) 
     712         ALLOCATE(nvarsprof(nproftypes)) 
     713         ALLOCATE(nextrprof(nproftypes)) 
     714 
     715         DO jtype = 1, nproftypes 
     716 
     717            IF ( TRIM(cobstypesprof(jtype)) == 'prof' ) THEN 
     718               nvarsprof(jtype) = 2 
     719               nextrprof(jtype) = 1 
     720               ALLOCATE(llvar(nvarsprof(jtype))) 
     721               CALL wrk_alloc( jpi, jpj,      nvarsprof(jtype), zglam ) 
     722               CALL wrk_alloc( jpi, jpj,      nvarsprof(jtype), zgphi ) 
     723               CALL wrk_alloc( jpi, jpj, jpk, nvarsprof(jtype), zmask ) 
     724               llvar(1)       = ln_t3d 
     725               llvar(2)       = ln_s3d 
     726               zglam(:,:,1)   = glamt(:,:) 
     727               zglam(:,:,2)   = glamt(:,:) 
     728               zgphi(:,:,1)   = gphit(:,:) 
     729               zgphi(:,:,2)   = gphit(:,:) 
     730               zmask(:,:,:,1) = tmask(:,:,:) 
     731               zmask(:,:,:,2) = tmask(:,:,:) 
     732            ELSE IF ( TRIM(cobstypesprof(jtype)) == 'vel' )  THEN 
     733               nvarsprof(jtype) = 2 
     734               nextrprof(jtype) = 2 
     735               ALLOCATE(llvar(nvarsprof(jtype))) 
     736               CALL wrk_alloc( jpi, jpj,      nvarsprof(jtype), zglam ) 
     737               CALL wrk_alloc( jpi, jpj,      nvarsprof(jtype), zgphi ) 
     738               CALL wrk_alloc( jpi, jpj, jpk, nvarsprof(jtype), zmask ) 
     739               llvar(1)       = ln_vel3d 
     740               llvar(2)       = ln_vel3d 
     741               zglam(:,:,1)   = glamu(:,:) 
     742               zglam(:,:,2)   = glamv(:,:) 
     743               zgphi(:,:,1)   = gphiu(:,:) 
     744               zgphi(:,:,2)   = gphiv(:,:) 
     745               zmask(:,:,:,1) = umask(:,:,:) 
     746               zmask(:,:,:,2) = vmask(:,:,:) 
     747            ELSE 
     748               nvarsprof(jtype) = 1 
     749               nextrprof(jtype) = 0 
     750               ALLOCATE(llvar(nvarsprof(jtype))) 
     751               CALL wrk_alloc( jpi, jpj,      nvarsprof(jtype), zglam ) 
     752               CALL wrk_alloc( jpi, jpj,      nvarsprof(jtype), zgphi ) 
     753               CALL wrk_alloc( jpi, jpj, jpk, nvarsprof(jtype), zmask ) 
     754               llvar(1)       = .TRUE. 
     755               zglam(:,:,1)   = glamt(:,:) 
     756               zgphi(:,:,1)   = gphit(:,:) 
     757               zmask(:,:,:,1) = tmask(:,:,:) 
     758            ENDIF 
     759 
     760            !Read in profile or profile obs types 
     761            CALL obs_rea_prof( profdata(jtype), ifilesprof(jtype),       & 
     762               &               clproffiles(jtype,1:ifilesprof(jtype)), & 
     763               &               nvarsprof(jtype), nextrprof(jtype), nitend-nit000+2, & 
     764               &               rn_dobsini, rn_dobsend, llvar, & 
     765               &               ln_ignmis, ln_s_at_t, .FALSE., & 
     766               &               kdailyavtypes = nn_profdavtypes ) 
     767 
     768            DO jvar = 1, nvarsprof(jtype) 
     769               CALL obs_prof_staend( profdata(jtype), jvar ) 
     770            END DO 
     771 
     772            CALL obs_pre_prof( profdata(jtype), profdataqc(jtype), & 
     773               &               llvar, & 
     774               &               jpi, jpj, jpk, & 
     775               &               zmask, zglam, zgphi,  & 
     776               &               ln_nea, ln_bound_reject, & 
     777               &               kdailyavtypes = nn_profdavtypes ) 
    524778             
    525             ld_enact(jprofset) = .TRUE. 
    526  
    527             CALL obs_rea_pro_dri( 1, profdata(jprofset),          & 
    528                &                  jnumenact, enactfiles(1:jnumenact), & 
    529                &                  nprofvars, nprofextr,        & 
    530                &                  nitend-nit000+2,             & 
    531                &                  dobsini, dobsend, ln_t3d, ln_s3d, & 
    532                &                  ln_ignmis, ln_s_at_t, .TRUE., .FALSE., & 
    533                &                  kdailyavtypes = endailyavtypes ) 
    534  
    535             DO jvar = 1, 2 
    536  
    537                CALL obs_prof_staend( profdata(jprofset), jvar ) 
    538  
    539             END DO 
    540  
    541             CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset),   & 
    542                &              ln_t3d, ln_s3d, ln_nea, & 
    543                &              kdailyavtypes=endailyavtypes ) 
    544              
    545          ENDIF 
    546  
    547          ! Coriolis insitu data 
    548  
    549          IF ( ln_cor ) THEN 
    550             
    551             jprofset = jprofset + 1 
    552  
    553             ld_enact(jprofset) = .FALSE. 
    554  
    555             CALL obs_rea_pro_dri( 2, profdata(jprofset),          & 
    556                &                  jnumcorio, coriofiles(1:jnumcorio), & 
    557                &                  nprofvars, nprofextr,        & 
    558                &                  nitend-nit000+2,             & 
    559                &                  dobsini, dobsend, ln_t3d, ln_s3d, & 
    560                &                  ln_ignmis, ln_s_at_t, .FALSE., .FALSE. ) 
    561  
    562             DO jvar = 1, 2 
    563  
    564                CALL obs_prof_staend( profdata(jprofset), jvar ) 
    565  
    566             END DO 
    567  
    568             CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset),   & 
    569                  &            ln_t3d, ln_s3d, ln_nea ) 
    570              
    571          ENDIF 
    572   
    573          ! Feedback insitu data 
    574  
    575          IF ( ln_profb ) THEN 
    576             
    577             DO jset = 1, jnumprofb 
    578                 
    579                jprofset = jprofset + 1 
    580                ld_enact (jprofset) = ln_profb_ena(jset) 
    581  
    582                CALL obs_rea_pro_dri( 0, profdata(jprofset),          & 
    583                   &                  1, profbfiles(jset:jset), & 
    584                   &                  nprofvars, nprofextr,        & 
    585                   &                  nitend-nit000+2,             & 
    586                   &                  dobsini, dobsend, ln_t3d, ln_s3d, & 
    587                   &                  ln_ignmis, ln_s_at_t, & 
    588                   &                  ld_enact(jprofset).AND.& 
    589                   &                  ln_profb_enatim(jset), & 
    590                   &                  .FALSE., kdailyavtypes = endailyavtypes ) 
    591                 
    592                DO jvar = 1, 2 
    593                    
    594                   CALL obs_prof_staend( profdata(jprofset), jvar ) 
    595                    
     779            DEALLOCATE( llvar ) 
     780            CALL wrk_dealloc( jpi, jpj,      nvarsprof(jtype), zglam ) 
     781            CALL wrk_dealloc( jpi, jpj,      nvarsprof(jtype), zgphi ) 
     782            CALL wrk_dealloc( jpi, jpj, jpk, nvarsprof(jtype), zmask ) 
     783 
     784         END DO 
     785 
     786         DEALLOCATE( ifilesprof, clproffiles ) 
     787 
     788      ENDIF 
     789 
     790      IF ( nsurftypes > 0 ) THEN 
     791 
     792         ALLOCATE(surfdata(nsurftypes)) 
     793         ALLOCATE(surfdataqc(nsurftypes)) 
     794         ALLOCATE(nvarssurf(nsurftypes)) 
     795         ALLOCATE(nextrsurf(nsurftypes)) 
     796 
     797         DO jtype = 1, nsurftypes 
     798 
     799            nvarssurf(jtype) = 1 
     800            nextrsurf(jtype) = 0 
     801            IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) nextrsurf(jtype) = 2 
     802 
     803            !Read in surface obs types 
     804            CALL obs_rea_surf( surfdata(jtype), ifilessurf(jtype), & 
     805               &               clsurffiles(jtype,1:ifilessurf(jtype)), & 
     806               &               nvarssurf(jtype), nextrsurf(jtype), nitend-nit000+2, & 
     807               &               rn_dobsini, rn_dobsend, ln_ignmis, .FALSE., llnightav(jtype) ) 
     808 
     809            CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea, ln_bound_reject ) 
     810 
     811            IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN 
     812               CALL obs_rea_mdt( surfdataqc(jtype), n2dintsurf(jtype) ) 
     813               IF ( ln_altbias ) & 
     814                  & CALL obs_rea_altbias ( surfdataqc(jtype), n2dintsurf(jtype), cn_altbiasfile ) 
     815            ENDIF 
     816 
     817            IF ( TRIM(cobstypessurf(jtype)) == 'sst' .AND. ln_sstbias ) THEN 
     818               jnumsstbias = 0 
     819               DO jfile = 1, jpmaxnfiles 
     820                  IF ( TRIM(cn_sstbiasfiles(jfile)) /= '' ) & 
     821                     &  jnumsstbias = jnumsstbias + 1 
    596822               END DO 
    597                 
    598                IF ( ld_enact(jprofset) ) THEN 
    599                   CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset),   & 
    600                      &              ln_t3d, ln_s3d, ln_nea, & 
    601                      &              kdailyavtypes = endailyavtypes ) 
    602                ELSE 
    603                   CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset),   & 
    604                      &              ln_t3d, ln_s3d, ln_nea ) 
     823               IF ( jnumsstbias == 0 ) THEN 
     824                  CALL ctl_stop("ln_sstbias set but no bias files to read in")     
    605825               ENDIF 
    606                 
    607             END DO 
    608  
    609          ENDIF 
     826 
     827               CALL obs_app_sstbias( surfdataqc(jtype), n2dintsurf(jtype), &  
     828                  &                  jnumsstbias, cn_sstbiasfiles(1:jnumsstbias) )  
     829 
     830            ENDIF 
     831 
     832         END DO 
     833 
     834         DEALLOCATE( ifilessurf, clsurffiles ) 
    610835 
    611836      ENDIF 
    612837 
    613       !  - Sea level anomalies 
    614       IF ( ln_sla ) THEN 
    615         ! Set the number of variables for sla to 1 
    616          nslavars = 1 
    617  
    618          ! Set the number of extra variables for sla to 2 
    619          nslaextr = 2 
    620           
    621          ! Set the number of sla data sets to 2 
    622          nslasets = 0 
    623          IF ( ln_sladt ) THEN 
    624             nslasets = nslasets + 2 
    625          ENDIF 
    626          IF ( ln_slafb ) THEN 
    627             nslasets = nslasets + jnumslafb 
    628          ENDIF 
    629           
    630          ALLOCATE(sladata(nslasets)) 
    631          ALLOCATE(sladatqc(nslasets)) 
    632          sladata(:)%nsurf=0 
    633          sladatqc(:)%nsurf=0 
    634  
    635          nslasets = 0 
    636  
    637          ! AVISO SLA data 
    638  
    639          IF ( ln_sladt ) THEN 
    640  
    641             ! Active SLA observations 
    642              
    643             nslasets = nslasets + 1 
    644              
    645             CALL obs_rea_sla( 1, sladata(nslasets), jnumslaact, & 
    646                &              slafilesact(1:jnumslaact), & 
    647                &              nslavars, nslaextr, nitend-nit000+2, & 
    648                &              dobsini, dobsend, ln_ignmis, .FALSE. ) 
    649             CALL obs_pre_sla( sladata(nslasets), sladatqc(nslasets), & 
    650                &              ln_sla, ln_nea ) 
    651              
    652             ! Passive SLA observations 
    653              
    654             nslasets = nslasets + 1 
    655              
    656             CALL obs_rea_sla( 1, sladata(nslasets), jnumslapas, & 
    657                &              slafilespas(1:jnumslapas), & 
    658                &              nslavars, nslaextr, nitend-nit000+2, & 
    659                &              dobsini, dobsend, ln_ignmis, .FALSE. ) 
    660              
    661             CALL obs_pre_sla( sladata(nslasets), sladatqc(nslasets), & 
    662                &              ln_sla, ln_nea ) 
    663  
    664          ENDIF 
    665           
    666          ! Feedback SLA data 
    667  
    668          IF ( ln_slafb ) THEN 
    669  
    670             DO jset = 1, jnumslafb 
    671              
    672                nslasets = nslasets + 1 
    673              
    674                CALL obs_rea_sla( 0, sladata(nslasets), 1, & 
    675                   &              slafbfiles(jset:jset), & 
    676                   &              nslavars, nslaextr, nitend-nit000+2, & 
    677                   &              dobsini, dobsend, ln_ignmis, .FALSE. ) 
    678                CALL obs_pre_sla( sladata(nslasets), sladatqc(nslasets), & 
    679                   &              ln_sla, ln_nea ) 
    680  
    681             END DO                
    682  
    683          ENDIF 
    684           
    685          CALL obs_rea_mdt( nslasets, sladatqc, n2dint ) 
    686              
    687          ! read in altimeter bias 
    688           
    689          IF ( ln_altbias ) THEN      
    690             CALL obs_rea_altbias ( nslasets, sladatqc, n2dint, bias_file ) 
    691          ENDIF 
    692       
    693       ENDIF 
    694  
    695       !  - Sea surface height 
    696       IF ( ln_ssh ) THEN 
    697          IF(lwp) WRITE(numout,*) ' SSH currently not available' 
    698       ENDIF 
    699  
    700       !  - Sea surface temperature 
    701       IF ( ln_sst ) THEN 
    702  
    703          ! Set the number of variables for sst to 1 
    704          nsstvars = 1 
    705  
    706          ! Set the number of extra variables for sst to 0 
    707          nsstextr = 0 
    708  
    709          nsstsets = 0 
    710  
    711          IF (ln_reysst) nsstsets = nsstsets + 1 
    712          IF (ln_ghrsst) nsstsets = nsstsets + 1 
    713          IF ( ln_sstfb ) THEN 
    714             nsstsets = nsstsets + jnumsstfb 
    715          ENDIF 
    716  
    717          ALLOCATE(sstdata(nsstsets)) 
    718          ALLOCATE(sstdatqc(nsstsets)) 
    719          ALLOCATE(ld_sstnight(nsstsets)) 
    720          sstdata(:)%nsurf=0 
    721          sstdatqc(:)%nsurf=0     
    722          ld_sstnight(:)=.false. 
    723  
    724          nsstsets = 0 
    725  
    726          IF (ln_reysst) THEN 
    727  
    728             nsstsets = nsstsets + 1 
    729  
    730             ld_sstnight(nsstsets) = ln_sstnight 
    731  
    732             CALL obs_rea_sst_rey( reysstname, reysstfmt, sstdata(nsstsets), & 
    733                &                  nsstvars, nsstextr, & 
    734                &                  nitend-nit000+2, dobsini, dobsend ) 
    735             CALL obs_pre_sst( sstdata(nsstsets), sstdatqc(nsstsets), ln_sst, & 
    736                &              ln_nea ) 
    737  
    738         ENDIF 
    739          
    740         IF (ln_ghrsst) THEN 
    741          
    742             nsstsets = nsstsets + 1 
    743  
    744             ld_sstnight(nsstsets) = ln_sstnight 
    745            
    746             CALL obs_rea_sst( 1, sstdata(nsstsets), jnumsst, & 
    747                &              sstfiles(1:jnumsst), & 
    748                &              nsstvars, nsstextr, nitend-nit000+2, & 
    749                &              dobsini, dobsend, ln_ignmis, .FALSE. ) 
    750             CALL obs_pre_sst( sstdata(nsstsets), sstdatqc(nsstsets), ln_sst, & 
    751                &              ln_nea ) 
    752  
    753         ENDIF 
    754                 
    755          ! Feedback SST data 
    756  
    757          IF ( ln_sstfb ) THEN 
    758  
    759             DO jset = 1, jnumsstfb 
    760              
    761                nsstsets = nsstsets + 1 
    762  
    763                ld_sstnight(nsstsets) = ln_sstnight 
    764              
    765                CALL obs_rea_sst( 0, sstdata(nsstsets), 1, & 
    766                   &              sstfbfiles(jset:jset), & 
    767                   &              nsstvars, nsstextr, nitend-nit000+2, & 
    768                   &              dobsini, dobsend, ln_ignmis, .FALSE. ) 
    769                CALL obs_pre_sst( sstdata(nsstsets), sstdatqc(nsstsets), & 
    770                   &              ln_sst, ln_nea ) 
    771  
    772             END DO                
    773  
    774          ENDIF 
    775  
    776       ENDIF 
    777  
    778       !  - Sea surface salinity 
    779       IF ( ln_sss ) THEN 
    780          IF(lwp) WRITE(numout,*) ' SSS currently not available' 
    781       ENDIF 
    782  
    783       !  - Sea Ice Concentration 
    784        
    785       IF ( ln_seaice ) THEN 
    786  
    787          ! Set the number of variables for seaice to 1 
    788          nseaicevars = 1 
    789  
    790          ! Set the number of extra variables for seaice to 0 
    791          nseaiceextr = 0 
    792           
    793          ! Set the number of data sets to 1 
    794          nseaicesets = 1 
    795  
    796          ALLOCATE(seaicedata(nseaicesets)) 
    797          ALLOCATE(seaicedatqc(nseaicesets)) 
    798          seaicedata(:)%nsurf=0 
    799          seaicedatqc(:)%nsurf=0 
    800  
    801          CALL obs_rea_seaice( 1, seaicedata(nseaicesets), jnumseaice, & 
    802             &                 seaicefiles(1:jnumseaice), & 
    803             &                 nseaicevars, nseaiceextr, nitend-nit000+2, & 
    804             &                 dobsini, dobsend, ln_ignmis, .FALSE. ) 
    805  
    806          CALL obs_pre_seaice( seaicedata(nseaicesets), seaicedatqc(nseaicesets), & 
    807             &                 ln_seaice, ln_nea ) 
    808   
    809       ENDIF 
    810  
    811       IF (ln_vel3d) THEN 
    812  
    813          ! Set the number of variables for profiles to 2 (U and V) 
    814          nvelovars = 2 
    815  
    816          ! Set the number of extra variables for profiles to 2 to store  
    817          ! rotation parameters 
    818          nveloextr = 2 
    819  
    820          jveloset = 0 
    821           
    822          IF ( ln_velavcur ) jveloset = jveloset + 1 
    823          IF ( ln_velhrcur ) jveloset = jveloset + 1 
    824          IF ( ln_velavadcp ) jveloset = jveloset + 1 
    825          IF ( ln_velhradcp ) jveloset = jveloset + 1 
    826          IF (ln_velfb) jveloset = jveloset + jnumvelfb 
    827  
    828          nvelosets = jveloset 
    829          IF ( nvelosets > 0 ) THEN 
    830             ALLOCATE( velodata(nvelosets) ) 
    831             ALLOCATE( veldatqc(nvelosets) ) 
    832             ALLOCATE( ld_velav(nvelosets) ) 
    833          ENDIF 
    834           
    835          jveloset = 0 
    836           
    837          ! Daily averaged data 
    838  
    839          IF ( ln_velavcur ) THEN 
    840              
    841             jveloset = jveloset + 1 
    842              
    843             ld_velav(jveloset) = .TRUE. 
    844              
    845             CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelavcur, & 
    846                &                  velavcurfiles(1:jnumvelavcur), & 
    847                &                  nvelovars, nveloextr, & 
    848                &                  nitend-nit000+2,              & 
    849                &                  dobsini, dobsend, ln_ignmis, & 
    850                &                  ld_velav(jveloset), & 
    851                &                  .FALSE. ) 
    852              
    853             DO jvar = 1, 2 
    854                CALL obs_prof_staend( velodata(jveloset), jvar ) 
    855             END DO 
    856              
    857             CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), & 
    858                &              ln_vel3d, ln_nea, ld_velav(jveloset) ) 
    859              
    860          ENDIF 
    861  
    862          ! High frequency data 
    863  
    864          IF ( ln_velhrcur ) THEN 
    865              
    866             jveloset = jveloset + 1 
    867              
    868             ld_velav(jveloset) = .FALSE. 
    869                 
    870             CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelhrcur, & 
    871                &                  velhrcurfiles(1:jnumvelhrcur), & 
    872                &                  nvelovars, nveloextr, & 
    873                &                  nitend-nit000+2,              & 
    874                &                  dobsini, dobsend, ln_ignmis, & 
    875                &                  ld_velav(jveloset), & 
    876                &                  .FALSE. ) 
    877              
    878             DO jvar = 1, 2 
    879                CALL obs_prof_staend( velodata(jveloset), jvar ) 
    880             END DO 
    881              
    882             CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), & 
    883                &              ln_vel3d, ln_nea, ld_velav(jveloset) ) 
    884              
    885          ENDIF 
    886  
    887          ! Daily averaged data 
    888  
    889          IF ( ln_velavadcp ) THEN 
    890              
    891             jveloset = jveloset + 1 
    892              
    893             ld_velav(jveloset) = .TRUE. 
    894              
    895             CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelavadcp, & 
    896                &                  velavadcpfiles(1:jnumvelavadcp), & 
    897                &                  nvelovars, nveloextr, & 
    898                &                  nitend-nit000+2,              & 
    899                &                  dobsini, dobsend, ln_ignmis, & 
    900                &                  ld_velav(jveloset), & 
    901                &                  .FALSE. ) 
    902              
    903             DO jvar = 1, 2 
    904                CALL obs_prof_staend( velodata(jveloset), jvar ) 
    905             END DO 
    906              
    907             CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), & 
    908                &              ln_vel3d, ln_nea, ld_velav(jveloset) ) 
    909              
    910          ENDIF 
    911  
    912          ! High frequency data 
    913  
    914          IF ( ln_velhradcp ) THEN 
    915              
    916             jveloset = jveloset + 1 
    917              
    918             ld_velav(jveloset) = .FALSE. 
    919                 
    920             CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelhradcp, & 
    921                &                  velhradcpfiles(1:jnumvelhradcp), & 
    922                &                  nvelovars, nveloextr, & 
    923                &                  nitend-nit000+2,              & 
    924                &                  dobsini, dobsend, ln_ignmis, & 
    925                &                  ld_velav(jveloset), & 
    926                &                  .FALSE. ) 
    927              
    928             DO jvar = 1, 2 
    929                CALL obs_prof_staend( velodata(jveloset), jvar ) 
    930             END DO 
    931              
    932             CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), & 
    933                &              ln_vel3d, ln_nea, ld_velav(jveloset) ) 
    934              
    935          ENDIF 
    936  
    937          IF ( ln_velfb ) THEN 
    938  
    939             DO jset = 1, jnumvelfb 
    940              
    941                jveloset = jveloset + 1 
    942  
    943                ld_velav(jveloset) = ln_velfb_av(jset) 
    944                 
    945                CALL obs_rea_vel_dri( 0, velodata(jveloset), 1, & 
    946                   &                  velfbfiles(jset:jset), & 
    947                   &                  nvelovars, nveloextr, & 
    948                   &                  nitend-nit000+2,              & 
    949                   &                  dobsini, dobsend, ln_ignmis, & 
    950                   &                  ld_velav(jveloset), & 
    951                   &                  .FALSE. ) 
    952                 
    953                DO jvar = 1, 2 
    954                   CALL obs_prof_staend( velodata(jveloset), jvar ) 
    955                END DO 
    956                 
    957                CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), & 
    958                   &              ln_vel3d, ln_nea, ld_velav(jveloset) ) 
    959  
    960  
    961             END DO 
    962              
    963          ENDIF 
    964  
    965       ENDIF 
    966       
    967838   END SUBROUTINE dia_obs_init 
    968839 
    … …  
    974845      !! 
    975846      !! ** Method  : Call the observation operators on each time step to 
    976       !!              compute the model equivalent of the following date: 
    977       !!               - T profiles 
    978       !!               - S profiles 
    979       !!               - Sea surface height (referenced to a mean) 
    980       !!               - Sea surface temperature 
    981       !!               - Sea surface salinity 
    982       !!               - Velocity component (U,V) profiles 
    983       !! 
    984       !! ** Action  :  
     847      !!              compute the model equivalent of the following data: 
     848      !!               - Profile data, currently T/S or U/V 
     849      !!               - Surface data, currently SST, SLA or sea-ice concentration. 
     850      !! 
     851      !! ** Action  : 
    985852      !! 
    986853      !! History : 
    … …  
    991858      !!        !  07-04  (G. Smith) Generalized surface operators 
    992859      !!        !  08-10  (M. Valdivieso) obs operator for velocity profiles 
     860      !!        !  15-08  (M. Martin) Combined surface/profile routines. 
    993861      !!---------------------------------------------------------------------- 
    994862      !! * Modules used 
    995       USE dom_oce, ONLY : &             ! Ocean space and time domain variables 
    996          & rdt,           &                        
    997          & gdept_1d,       &              
    998          & tmask, umask, vmask                             
    999       USE phycst, ONLY : &              ! Physical constants 
    1000          & rday                          
    1001       USE oce, ONLY : &                 ! Ocean dynamics and tracers variables 
    1002          & tsn,  &              
    1003          & un, vn,  & 
     863      USE phycst, ONLY : &         ! Physical constants 
     864#if defined key_fabm 
     865         & rt0,          & 
     866#endif 
     867         & rday 
     868      USE oce, ONLY : &            ! Ocean dynamics and tracers variables 
     869         & tsn,       & 
     870         & un,        & 
     871         & vn,        & 
    1004872         & sshn 
    1005873#if defined  key_lim3 
    1006       USE ice, ONLY : &                     ! LIM Ice model variables 
     874      USE ice, ONLY : &            ! LIM3 Ice model variables 
    1007875         & frld 
    1008876#endif 
    1009877#if defined key_lim2 
    1010       USE ice_2, ONLY : &                     ! LIM Ice model variables 
     878      USE ice_2, ONLY : &          ! LIM2 Ice model variables 
    1011879         & frld 
    1012880#endif 
     881#if defined key_cice 
     882      USE sbc_oce, ONLY : fr_i     ! ice fraction 
     883#endif 
     884#if defined key_top 
     885      USE trc, ONLY :  &           ! Biogeochemical state variables 
     886         & trn 
     887#endif 
     888#if defined key_hadocc 
     889      USE par_hadocc               ! HadOCC parameters 
     890      USE trc, ONLY :  & 
     891         & HADOCC_CHL, & 
     892         & HADOCC_FCO2, & 
     893         & HADOCC_PCO2, & 
     894         & HADOCC_FILL_FLT 
     895      USE had_bgc_const, ONLY: c2n_p 
     896#elif defined key_medusa 
     897      USE par_medusa               ! MEDUSA parameters 
     898      USE sms_medusa, ONLY: & 
     899         & xthetapn, & 
     900         & xthetapd 
     901#if defined key_roam 
     902      USE sms_medusa, ONLY: & 
     903         & f2_pco2w, & 
     904         & f2_fco2w, & 
     905         & f3_pH 
     906#endif 
     907#elif defined key_fabm 
     908      USE par_fabm                 ! FABM parameters 
     909      USE fabm, ONLY: & 
     910         & fabm_get_interior_diagnostic_data 
     911#endif 
     912#if defined key_spm 
     913      USE par_spm, ONLY: &         ! Sediment parameters 
     914         & jp_spm 
     915#endif 
     916 
    1013917      IMPLICIT NONE 
    1014918 
    1015919      !! * Arguments 
    1016       INTEGER, INTENT(IN) :: kstp                         ! Current timestep 
     920      INTEGER, INTENT(IN) :: kstp  ! Current timestep 
    1017921      !! * Local declarations 
    1018       INTEGER :: idaystp                ! Number of timesteps per day 
    1019       INTEGER :: jprofset               ! Profile data set loop variable 
    1020       INTEGER :: jslaset                ! SLA data set loop variable 
    1021       INTEGER :: jsstset                ! SST data set loop variable 
    1022       INTEGER :: jseaiceset             ! sea ice data set loop variable 
    1023       INTEGER :: jveloset               ! velocity profile data loop variable 
    1024       INTEGER :: jvar                   ! Variable number     
    1025 #if ! defined key_lim2 && ! defined key_lim3 
    1026       REAL(wp), POINTER, DIMENSION(:,:) :: frld    
    1027 #endif 
    1028       CHARACTER(LEN=20) :: datestr=" ",timestr=" " 
    1029   
    1030 #if ! defined key_lim2 && ! defined key_lim3 
    1031       CALL wrk_alloc(jpi,jpj,frld)  
     922      INTEGER :: idaystp           ! Number of timesteps per day 
     923      INTEGER :: jtype             ! Data loop variable 
     924      INTEGER :: jvar              ! Variable number 
     925      INTEGER :: ji, jj, jk        ! Loop counters 
     926      REAL(wp) :: tiny             ! small number 
     927      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: & 
     928         & zprofvar                ! Model values for variables in a prof ob 
     929      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: & 
     930         & zprofmask               ! Mask associated with zprofvar 
     931      REAL(wp), POINTER, DIMENSION(:,:) :: & 
     932         & zsurfvar, &             ! Model values equivalent to surface ob. 
     933         & zsurfmask               ! Mask associated with surface variable 
     934      REAL(wp), POINTER, DIMENSION(:,:,:) :: & 
     935         & zglam,    &             ! Model longitudes for prof variables 
     936         & zgphi                   ! Model latitudes for prof variables 
     937      LOGICAL :: llog10            ! Perform log10 transform of variable 
     938#if defined key_fabm 
     939      REAL(wp), POINTER, DIMENSION(:,:,:) :: & 
     940         & pco2_3d                 ! 3D pCO2 from FABM 
    1032941#endif 
    1033942 
    … …  
    1036945         WRITE(numout,*) 'dia_obs : Call the observation operators', kstp 
    1037946         WRITE(numout,*) '~~~~~~~' 
     947         CALL FLUSH(numout) 
    1038948      ENDIF 
    1039949 
    … …  
    1041951 
    1042952      !----------------------------------------------------------------------- 
    1043       ! No LIM => frld == 0.0_wp 
    1044       !----------------------------------------------------------------------- 
    1045 #if ! defined key_lim2 && ! defined key_lim3 
    1046       frld(:,:) = 0.0_wp 
    1047 #endif 
    1048       !----------------------------------------------------------------------- 
    1049       ! Depending on switches call various observation operators 
    1050       !----------------------------------------------------------------------- 
    1051  
    1052       !  - Temperature/salinity profiles 
    1053       IF ( ln_t3d .OR. ln_s3d ) THEN 
    1054          DO jprofset = 1, nprofsets 
    1055             IF ( ld_enact(jprofset) ) THEN 
    1056                CALL obs_pro_opt( prodatqc(jprofset),                     & 
    1057                   &              kstp, jpi, jpj, jpk, nit000, idaystp,   & 
    1058                   &              tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal),   & 
    1059                   &              gdept_1d, tmask, n1dint, n2dint,        & 
    1060                   &              kdailyavtypes = endailyavtypes ) 
    1061             ELSE 
    1062                CALL obs_pro_opt( prodatqc(jprofset),                     & 
    1063                   &              kstp, jpi, jpj, jpk, nit000, idaystp,   & 
    1064                   &              tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal),   & 
    1065                   &              gdept_1d, tmask, n1dint, n2dint              ) 
     953      ! Call the profile and surface observation operators 
     954      !----------------------------------------------------------------------- 
     955 
     956      IF ( nproftypes > 0 ) THEN 
     957 
     958         DO jtype = 1, nproftypes 
     959 
     960            ! Allocate local work arrays 
     961            CALL wrk_alloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofvar  ) 
     962            CALL wrk_alloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofmask ) 
     963            CALL wrk_alloc( jpi, jpj,      profdataqc(jtype)%nvar, zglam     ) 
     964            CALL wrk_alloc( jpi, jpj,      profdataqc(jtype)%nvar, zgphi     ) 
     965             
     966            ! Defaults which might change 
     967            DO jvar = 1, profdataqc(jtype)%nvar 
     968               zprofmask(:,:,:,jvar) = tmask(:,:,:) 
     969               zglam(:,:,jvar)       = glamt(:,:) 
     970               zgphi(:,:,jvar)       = gphit(:,:) 
     971            END DO 
     972 
     973            SELECT CASE ( TRIM(cobstypesprof(jtype)) ) 
     974 
     975            CASE('prof') 
     976               zprofvar(:,:,:,1) = tsn(:,:,:,jp_tem) 
     977               zprofvar(:,:,:,2) = tsn(:,:,:,jp_sal) 
     978 
     979            CASE('vel') 
     980               zprofvar(:,:,:,1) = un(:,:,:) 
     981               zprofvar(:,:,:,2) = vn(:,:,:) 
     982               zprofmask(:,:,:,1) = umask(:,:,:) 
     983               zprofmask(:,:,:,2) = vmask(:,:,:) 
     984               zglam(:,:,1) = glamu(:,:) 
     985               zglam(:,:,2) = glamv(:,:) 
     986               zgphi(:,:,1) = gphiu(:,:) 
     987               zgphi(:,:,2) = gphiv(:,:) 
     988 
     989            CASE('plchltot') 
     990#if defined key_hadocc 
     991               ! Chlorophyll from HadOCC 
     992               zprofvar(:,:,:,1) = HADOCC_CHL(:,:,:) 
     993#elif defined key_medusa 
     994               ! Add non-diatom and diatom chlorophyll from MEDUSA 
     995               zprofvar(:,:,:,1) = trn(:,:,:,jpchn) + trn(:,:,:,jpchd) 
     996#elif defined key_fabm 
     997               ! Add all chlorophyll groups from ERSEM 
     998               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl2) + & 
     999                  &                trn(:,:,:,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl4) 
     1000#else 
     1001               CALL ctl_stop( ' Trying to run plchltot observation operator', & 
     1002                  &           ' but no biogeochemical model appears to have been defined' ) 
     1003#endif 
     1004               ! Take the log10 where we can, otherwise exclude 
     1005               tiny = 1.0e-20 
     1006               WHERE(zprofvar(:,:,:,:) > tiny .AND. zprofvar(:,:,:,:) /= obfillflt ) 
     1007                  zprofvar(:,:,:,:)  = LOG10(zprofvar(:,:,:,:)) 
     1008               ELSEWHERE 
     1009                  zprofvar(:,:,:,:)  = obfillflt 
     1010                  zprofmask(:,:,:,:) = 0 
     1011               END WHERE 
     1012               ! Mask out model below any excluded values, 
     1013               ! to avoid interpolation issues 
     1014               DO jvar = 1, profdataqc(jtype)%nvar 
     1015                 DO jj = 1, jpj 
     1016                    DO ji = 1, jpi 
     1017                       depth_loop: DO jk = 1, jpk 
     1018                          IF ( zprofmask(ji,jj,jk,jvar) == 0 ) THEN 
     1019                             zprofmask(ji,jj,jk:jpk,jvar) = 0 
     1020                             EXIT depth_loop 
     1021                          ENDIF 
     1022                       END DO depth_loop 
     1023                    END DO 
     1024                 END DO 
     1025              END DO 
     1026 
     1027            CASE('pchltot') 
     1028#if defined key_hadocc 
     1029               ! Chlorophyll from HadOCC 
     1030               zprofvar(:,:,:,1) = HADOCC_CHL(:,:,:) 
     1031#elif defined key_medusa 
     1032               ! Add non-diatom and diatom chlorophyll from MEDUSA 
     1033               zprofvar(:,:,:,1) = trn(:,:,:,jpchn) + trn(:,:,:,jpchd) 
     1034#elif defined key_fabm 
     1035               ! Add all chlorophyll groups from ERSEM 
     1036               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl2) + & 
     1037                  &                trn(:,:,:,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl4) 
     1038#else 
     1039               CALL ctl_stop( ' Trying to run pchltot observation operator', & 
     1040                  &           ' but no biogeochemical model appears to have been defined' ) 
     1041#endif 
     1042 
     1043            CASE('pno3') 
     1044#if defined key_hadocc 
     1045               ! Dissolved inorganic nitrogen from HadOCC 
     1046               zprofvar(:,:,:,1) = trn(:,:,:,jp_had_nut) 
     1047#elif defined key_medusa 
     1048               ! Dissolved inorganic nitrogen from MEDUSA 
     1049               zprofvar(:,:,:,1) = trn(:,:,:,jpdin) 
     1050#elif defined key_fabm 
     1051               ! Nitrate from ERSEM 
     1052               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_n3n) 
     1053#else 
     1054               CALL ctl_stop( ' Trying to run pno3 observation operator', & 
     1055                  &           ' but no biogeochemical model appears to have been defined' ) 
     1056#endif 
     1057 
     1058            CASE('psi4') 
     1059#if defined key_hadocc 
     1060               CALL ctl_stop( ' Trying to run psi4 observation operator', & 
     1061                  &           ' but HadOCC does not simulate silicate' ) 
     1062#elif defined key_medusa 
     1063               ! Silicate from MEDUSA 
     1064               zprofvar(:,:,:,1) = trn(:,:,:,jpsil) 
     1065#elif defined key_fabm 
     1066               ! Silicate from ERSEM 
     1067               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_n5s) 
     1068#else 
     1069               CALL ctl_stop( ' Trying to run psi4 observation operator', & 
     1070                  &           ' but no biogeochemical model appears to have been defined' ) 
     1071#endif 
     1072 
     1073            CASE('ppo4') 
     1074#if defined key_hadocc 
     1075               CALL ctl_stop( ' Trying to run ppo4 observation operator', & 
     1076                  &           ' but HadOCC does not simulate phosphate' ) 
     1077#elif defined key_medusa 
     1078               CALL ctl_stop( ' Trying to run ppo4 observation operator', & 
     1079                  &           ' but MEDUSA does not simulate phosphate' ) 
     1080#elif defined key_fabm 
     1081               ! Phosphate from ERSEM 
     1082               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_n1p) 
     1083#else 
     1084               CALL ctl_stop( ' Trying to run ppo4 observation operator', & 
     1085                  &           ' but no biogeochemical model appears to have been defined' ) 
     1086#endif 
     1087 
     1088            CASE('pdic') 
     1089#if defined key_hadocc 
     1090               ! Dissolved inorganic carbon from HadOCC 
     1091               zprofvar(:,:,:,1) = trn(:,:,:,jp_had_dic) 
     1092#elif defined key_medusa 
     1093               ! Dissolved inorganic carbon from MEDUSA 
     1094               zprofvar(:,:,:,1) = trn(:,:,:,jpdic) 
     1095#elif defined key_fabm 
     1096               ! Dissolved inorganic carbon from ERSEM 
     1097               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_o3c) 
     1098#else 
     1099               CALL ctl_stop( ' Trying to run pdic observation operator', & 
     1100                  &           ' but no biogeochemical model appears to have been defined' ) 
     1101#endif 
     1102 
     1103            CASE('palk') 
     1104#if defined key_hadocc 
     1105               ! Alkalinity from HadOCC 
     1106               zprofvar(:,:,:,1) = trn(:,:,:,jp_had_alk) 
     1107#elif defined key_medusa 
     1108               ! Alkalinity from MEDUSA 
     1109               zprofvar(:,:,:,1) = trn(:,:,:,jpalk) 
     1110#elif defined key_fabm 
     1111               ! Alkalinity from ERSEM 
     1112               zprofvar(:,:,:,1) = fabm_get_interior_diagnostic_data(model, jp_fabm_o3a) 
     1113#else 
     1114               CALL ctl_stop( ' Trying to run palk observation operator', & 
     1115                  &           ' but no biogeochemical model appears to have been defined' ) 
     1116#endif 
     1117 
     1118            CASE('pph') 
     1119#if defined key_hadocc 
     1120               CALL ctl_stop( ' Trying to run pph observation operator', & 
     1121                  &           ' but HadOCC has no pH diagnostic defined' ) 
     1122#elif defined key_medusa && defined key_roam 
     1123               ! pH from MEDUSA 
     1124               zprofvar(:,:,:,1) = f3_pH(:,:,:) 
     1125#elif defined key_fabm 
     1126               ! pH from ERSEM 
     1127               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_o3ph) 
     1128#else 
     1129               CALL ctl_stop( ' Trying to run pph observation operator', & 
     1130                  &           ' but no biogeochemical model appears to have been defined' ) 
     1131#endif 
     1132 
     1133            CASE('po2') 
     1134#if defined key_hadocc 
     1135               CALL ctl_stop( ' Trying to run po2 observation operator', & 
     1136                  &           ' but HadOCC does not simulate oxygen' ) 
     1137#elif defined key_medusa 
     1138               ! Oxygen from MEDUSA 
     1139               zprofvar(:,:,:,1) = trn(:,:,:,jpoxy) 
     1140#elif defined key_fabm 
     1141               ! Oxygen from ERSEM 
     1142               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_o2o) 
     1143#else 
     1144               CALL ctl_stop( ' Trying to run po2 observation operator', & 
     1145                  &           ' but no biogeochemical model appears to have been defined' ) 
     1146#endif 
     1147 
     1148            CASE DEFAULT 
     1149               CALL ctl_stop( 'Unknown profile observation type '//TRIM(cobstypesprof(jtype))//' in dia_obs' ) 
     1150 
     1151            END SELECT 
     1152 
     1153            DO jvar = 1, profdataqc(jtype)%nvar 
     1154               CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk,  & 
     1155                  &               nit000, idaystp, jvar,                   & 
     1156                  &               zprofvar(:,:,:,jvar),                    & 
     1157                  &               fsdept(:,:,:), fsdepw(:,:,:),            &  
     1158                  &               zprofmask(:,:,:,jvar),                   & 
     1159                  &               zglam(:,:,jvar), zgphi(:,:,jvar),        & 
     1160                  &               nn_1dint, nn_2dint_default,              & 
     1161                  &               kdailyavtypes = nn_profdavtypes ) 
     1162            END DO 
     1163 
     1164            CALL wrk_dealloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofvar  ) 
     1165            CALL wrk_dealloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofmask ) 
     1166            CALL wrk_dealloc( jpi, jpj,      profdataqc(jtype)%nvar, zglam     ) 
     1167            CALL wrk_dealloc( jpi, jpj,      profdataqc(jtype)%nvar, zgphi     ) 
     1168 
     1169         END DO 
     1170 
     1171      ENDIF 
     1172 
     1173      IF ( nsurftypes > 0 ) THEN 
     1174 
     1175         !Allocate local work arrays 
     1176         CALL wrk_alloc( jpi, jpj, zsurfvar ) 
     1177         CALL wrk_alloc( jpi, jpj, zsurfmask ) 
     1178#if defined key_fabm 
     1179         CALL wrk_alloc( jpi, jpj, jpk, pco2_3d ) 
     1180#endif 
     1181 
     1182         DO jtype = 1, nsurftypes 
     1183 
     1184            !Defaults which might be changed 
     1185            zsurfmask(:,:) = tmask(:,:,1) 
     1186            llog10 = .FALSE. 
     1187 
     1188            SELECT CASE ( TRIM(cobstypessurf(jtype)) ) 
     1189            CASE('sst') 
     1190               zsurfvar(:,:) = tsn(:,:,1,jp_tem) 
     1191            CASE('sla') 
     1192               zsurfvar(:,:) = sshn(:,:) 
     1193            CASE('sss') 
     1194               zsurfvar(:,:) = tsn(:,:,1,jp_sal) 
     1195            CASE('sic') 
     1196               IF ( kstp == 0 ) THEN 
     1197                  IF ( lwp .AND. surfdataqc(jtype)%nsstpmpp(1) > 0 ) THEN 
     1198                     CALL ctl_warn( 'Sea-ice not initialised on zeroth '// & 
     1199                        &           'time-step but some obs are valid then.' ) 
     1200                     WRITE(numout,*)surfdataqc(jtype)%nsstpmpp(1), & 
     1201                        &           ' sea-ice obs will be missed' 
     1202                  ENDIF 
     1203                  surfdataqc(jtype)%nsurfup = surfdataqc(jtype)%nsurfup + & 
     1204                     &                        surfdataqc(jtype)%nsstp(1) 
     1205                  CYCLE 
     1206               ELSE 
     1207#if defined key_cice 
     1208                  zsurfvar(:,:) = fr_i(:,:) 
     1209#elif defined key_lim2 || defined key_lim3 
     1210                  zsurfvar(:,:) = 1._wp - frld(:,:) 
     1211#else 
     1212               CALL ctl_stop( ' Trying to run sea-ice observation operator', & 
     1213                  &           ' but no sea-ice model appears to have been defined' ) 
     1214#endif 
     1215               ENDIF 
     1216 
     1217            CASE('slchltot') 
     1218#if defined key_hadocc 
     1219               ! Surface chlorophyll from HadOCC 
     1220               zsurfvar(:,:) = HADOCC_CHL(:,:,1) 
     1221#elif defined key_medusa 
     1222               ! Add non-diatom and diatom surface chlorophyll from MEDUSA 
     1223               zsurfvar(:,:) = trn(:,:,1,jpchn) + trn(:,:,1,jpchd) 
     1224#elif defined key_fabm 
     1225               ! Add all surface chlorophyll groups from ERSEM 
     1226               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) + & 
     1227                  &            trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4) 
     1228#else 
     1229               CALL ctl_stop( ' Trying to run slchltot observation operator', & 
     1230                  &           ' but no biogeochemical model appears to have been defined' ) 
     1231#endif 
     1232               llog10 = .TRUE. 
     1233 
     1234            CASE('slchldia') 
     1235#if defined key_hadocc 
     1236               CALL ctl_stop( ' Trying to run slchldia observation operator', & 
     1237                  &           ' but HadOCC does not explicitly simulate diatoms' ) 
     1238#elif defined key_medusa 
     1239               ! Diatom surface chlorophyll from MEDUSA 
     1240               zsurfvar(:,:) = trn(:,:,1,jpchd) 
     1241#elif defined key_fabm 
     1242               ! Diatom surface chlorophyll from ERSEM 
     1243               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) 
     1244#else 
     1245               CALL ctl_stop( ' Trying to run slchldia observation operator', & 
     1246                  &           ' but no biogeochemical model appears to have been defined' ) 
     1247#endif 
     1248               llog10 = .TRUE. 
     1249 
     1250            CASE('slchlnon') 
     1251#if defined key_hadocc 
     1252               CALL ctl_stop( ' Trying to run slchlnon observation operator', & 
     1253                  &           ' but HadOCC does not explicitly simulate non-diatoms' ) 
     1254#elif defined key_medusa 
     1255               ! Non-diatom surface chlorophyll from MEDUSA 
     1256               zsurfvar(:,:) = trn(:,:,1,jpchn) 
     1257#elif defined key_fabm 
     1258               ! Add all non-diatom surface chlorophyll groups from ERSEM 
     1259               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) + & 
     1260                  &            trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4) 
     1261#else 
     1262               CALL ctl_stop( ' Trying to run slchlnon observation operator', & 
     1263                  &           ' but no biogeochemical model appears to have been defined' ) 
     1264#endif 
     1265               llog10 = .TRUE. 
     1266 
     1267            CASE('slchldin') 
     1268#if defined key_hadocc 
     1269               CALL ctl_stop( ' Trying to run slchldin observation operator', & 
     1270                  &           ' but HadOCC does not explicitly simulate dinoflagellates' ) 
     1271#elif defined key_medusa 
     1272               CALL ctl_stop( ' Trying to run slchldin observation operator', & 
     1273                  &           ' but MEDUSA does not explicitly simulate dinoflagellates' ) 
     1274#elif defined key_fabm 
     1275               ! Dinoflagellate surface chlorophyll from ERSEM 
     1276               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl4) 
     1277#else 
     1278               CALL ctl_stop( ' Trying to run slchldin observation operator', & 
     1279                  &           ' but no biogeochemical model appears to have been defined' ) 
     1280#endif 
     1281               llog10 = .TRUE. 
     1282 
     1283            CASE('slchlmic') 
     1284#if defined key_hadocc 
     1285               CALL ctl_stop( ' Trying to run slchlmic observation operator', & 
     1286                  &           ' but HadOCC does not explicitly simulate microphytoplankton' ) 
     1287#elif defined key_medusa 
     1288               CALL ctl_stop( ' Trying to run slchlmic observation operator', & 
     1289                  &           ' but MEDUSA does not explicitly simulate microphytoplankton' ) 
     1290#elif defined key_fabm 
     1291               ! Add diatom and dinoflagellate surface chlorophyll from ERSEM 
     1292               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4) 
     1293#else 
     1294               CALL ctl_stop( ' Trying to run slchlmic observation operator', & 
     1295                  &           ' but no biogeochemical model appears to have been defined' ) 
     1296#endif 
     1297               llog10 = .TRUE. 
     1298 
     1299            CASE('slchlnan') 
     1300#if defined key_hadocc 
     1301               CALL ctl_stop( ' Trying to run slchlnan observation operator', & 
     1302                  &           ' but HadOCC does not explicitly simulate nanophytoplankton' ) 
     1303#elif defined key_medusa 
     1304               CALL ctl_stop( ' Trying to run slchlnan observation operator', & 
     1305                  &           ' but MEDUSA does not explicitly simulate nanophytoplankton' ) 
     1306#elif defined key_fabm 
     1307               ! Nanophytoplankton surface chlorophyll from ERSEM 
     1308               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) 
     1309#else 
     1310               CALL ctl_stop( ' Trying to run slchlnan observation operator', & 
     1311                  &           ' but no biogeochemical model appears to have been defined' ) 
     1312#endif 
     1313               llog10 = .TRUE. 
     1314 
     1315            CASE('slchlpic') 
     1316#if defined key_hadocc 
     1317               CALL ctl_stop( ' Trying to run slchlpic observation operator', & 
     1318                  &           ' but HadOCC does not explicitly simulate picophytoplankton' ) 
     1319#elif defined key_medusa 
     1320               CALL ctl_stop( ' Trying to run slchlpic observation operator', & 
     1321                  &           ' but MEDUSA does not explicitly simulate picophytoplankton' ) 
     1322#elif defined key_fabm 
     1323               ! Picophytoplankton surface chlorophyll from ERSEM 
     1324               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) 
     1325#else 
     1326               CALL ctl_stop( ' Trying to run slchlpic observation operator', & 
     1327                  &           ' but no biogeochemical model appears to have been defined' ) 
     1328#endif 
     1329               llog10 = .TRUE. 
     1330 
     1331            CASE('schltot') 
     1332#if defined key_hadocc 
     1333               ! Surface chlorophyll from HadOCC 
     1334               zsurfvar(:,:) = HADOCC_CHL(:,:,1) 
     1335#elif defined key_medusa 
     1336               ! Add non-diatom and diatom surface chlorophyll from MEDUSA 
     1337               zsurfvar(:,:) = trn(:,:,1,jpchn) + trn(:,:,1,jpchd) 
     1338#elif defined key_fabm 
     1339               ! Add all surface chlorophyll groups from ERSEM 
     1340               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) + & 
     1341                  &            trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4) 
     1342#else 
     1343               CALL ctl_stop( ' Trying to run schltot observation operator', & 
     1344                  &           ' but no biogeochemical model appears to have been defined' ) 
     1345#endif 
     1346 
     1347            CASE('slphytot') 
     1348#if defined key_hadocc 
     1349               ! Surface phytoplankton nitrogen from HadOCC multiplied by C:N ratio 
     1350               zsurfvar(:,:) = trn(:,:,1,jp_had_phy) * c2n_p 
     1351#elif defined key_medusa 
     1352               ! Add non-diatom and diatom surface phytoplankton nitrogen from MEDUSA 
     1353               ! multiplied by C:N ratio for each 
     1354               zsurfvar(:,:) = (trn(:,:,1,jpphn) * xthetapn) + (trn(:,:,1,jpphd) * xthetapd) 
     1355#elif defined key_fabm 
     1356               ! Add all surface phytoplankton carbon groups from ERSEM 
     1357               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_p1c) + trn(:,:,1,jp_fabm_m1+jp_fabm_p2c) + & 
     1358                  &            trn(:,:,1,jp_fabm_m1+jp_fabm_p3c) + trn(:,:,1,jp_fabm_m1+jp_fabm_p4c) 
     1359#else 
     1360               CALL ctl_stop( ' Trying to run slphytot observation operator', & 
     1361                  &           ' but no biogeochemical model appears to have been defined' ) 
     1362#endif 
     1363               llog10 = .TRUE. 
     1364 
     1365            CASE('slphydia') 
     1366#if defined key_hadocc 
     1367               CALL ctl_stop( ' Trying to run slphydia observation operator', & 
     1368                  &           ' but HadOCC does not explicitly simulate diatoms' ) 
     1369#elif defined key_medusa 
     1370               ! Diatom surface phytoplankton nitrogen from MEDUSA multiplied by C:N ratio 
     1371               zsurfvar(:,:) = trn(:,:,1,jpphd) * xthetapd 
     1372#elif defined key_fabm 
     1373               ! Diatom surface phytoplankton carbon from ERSEM 
     1374               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_p1c) 
     1375#else 
     1376               CALL ctl_stop( ' Trying to run slphydia observation operator', & 
     1377                  &           ' but no biogeochemical model appears to have been defined' ) 
     1378#endif 
     1379               llog10 = .TRUE. 
     1380 
     1381            CASE('slphynon') 
     1382#if defined key_hadocc 
     1383               CALL ctl_stop( ' Trying to run slphynon observation operator', & 
     1384                  &           ' but HadOCC does not explicitly simulate non-diatoms' ) 
     1385#elif defined key_medusa 
     1386               ! Non-diatom surface phytoplankton nitrogen from MEDUSA multiplied by C:N ratio 
     1387               zsurfvar(:,:) = trn(:,:,1,jpphn) * xthetapn 
     1388#elif defined key_fabm 
     1389               ! Add all non-diatom surface phytoplankton carbon groups from ERSEM 
     1390               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_p2c) + & 
     1391                  &            trn(:,:,1,jp_fabm_m1+jp_fabm_p3c) + trn(:,:,1,jp_fabm_m1+jp_fabm_p4c) 
     1392#else 
     1393               CALL ctl_stop( ' Trying to run slphynon observation operator', & 
     1394                  &           ' but no biogeochemical model appears to have been defined' ) 
     1395#endif 
     1396               llog10 = .TRUE. 
     1397 
     1398            CASE('sspm') 
     1399#if defined key_spm 
     1400               zsurfvar(:,:) = 0.0 
     1401               DO jn = 1, jp_spm 
     1402                  zsurfvar(:,:) = zsurfvar(:,:) + trn(:,:,1,jn)   ! sum SPM sizes 
     1403               END DO 
     1404#else 
     1405               CALL ctl_stop( ' Trying to run sspm observation operator', & 
     1406                  &           ' but no spm model appears to have been defined' ) 
     1407#endif 
     1408 
     1409            CASE('sfco2') 
     1410#if defined key_hadocc 
     1411               zsurfvar(:,:) = HADOCC_FCO2(:,:)    ! fCO2 from HadOCC 
     1412               IF ( ( MINVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) .AND. & 
     1413                  & ( MAXVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) ) THEN 
     1414                  zsurfvar(:,:) = obfillflt 
     1415                  zsurfmask(:,:) = 0 
     1416                  CALL ctl_warn( ' HadOCC fCO2 values masked out for observation operator', & 
     1417                     &           ' as HADOCC_FCO2(:,:) == HADOCC_FILL_FLT' ) 
     1418               ENDIF 
     1419#elif defined key_medusa && defined key_roam 
     1420               zsurfvar(:,:) = f2_fco2w(:,:) 
     1421#elif defined key_fabm 
     1422               ! First, get pCO2 from FABM 
     1423               pco2_3d(:,:,:) = fabm_get_interior_diagnostic_data(model, jp_fabm_o3pc) 
     1424               zsurfvar(:,:) = pco2_3d(:,:,1) 
     1425               ! Now, convert pCO2 to fCO2, based on SST in K. This follows the standard methodology of: 
     1426               ! Pierrot et al. (2009), Recommendations for autonomous underway pCO2 measuring systems 
     1427               ! and data reduction routines, Deep-Sea Research II, 56: 512-522. 
     1428               ! and 
     1429               ! Weiss (1974), Carbon dioxide in water and seawater: the solubility of a non-ideal gas, 
     1430               ! Marine Chemistry, 2: 203-215. 
     1431               ! In the implementation below, atmospheric pressure has been assumed to be 1 atm and so 
     1432               ! not explicitly included - atmospheric pressure is not necessarily available so this is 
     1433               ! the best assumption. 
     1434               ! Further, the (1-xCO2)^2 term has been neglected. This is common practice 
     1435               ! (see e.g. Zeebe and Wolf-Gladrow (2001), CO2 in Seawater: Equilibrium, Kinetics, Isotopes) 
     1436               ! because xCO2 in atm is ~0, and so this term will only affect the result to the 3rd decimal 
     1437               ! place for typical values, and xCO2 would need to be approximated from pCO2 anyway. 
     1438               zsurfvar(:,:) = zsurfvar(:,:) * EXP((-1636.75                                                          + & 
     1439                  &            12.0408      * (tsn(:,:,1,jp_tem)+rt0)                                                 - & 
     1440                  &            0.0327957    * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0)                         + & 
     1441                  &            0.0000316528 * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0) + & 
     1442                  &            2.0 * (57.7 - 0.118 * (tsn(:,:,1,jp_tem)+rt0)))                                        / & 
     1443                  &            (82.0578 * (tsn(:,:,1,jp_tem)+rt0))) 
     1444#else 
     1445               CALL ctl_stop( ' Trying to run sfco2 observation operator', & 
     1446                  &           ' but no biogeochemical model appears to have been defined' ) 
     1447#endif 
     1448 
     1449            CASE('spco2') 
     1450#if defined key_hadocc 
     1451               zsurfvar(:,:) = HADOCC_PCO2(:,:)    ! pCO2 from HadOCC 
     1452               IF ( ( MINVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) .AND. & 
     1453                  & ( MAXVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) ) THEN 
     1454                  zsurfvar(:,:) = obfillflt 
     1455                  zsurfmask(:,:) = 0 
     1456                  CALL ctl_warn( ' HadOCC pCO2 values masked out for observation operator', & 
     1457                     &           ' as HADOCC_PCO2(:,:) == HADOCC_FILL_FLT' ) 
     1458               ENDIF 
     1459#elif defined key_medusa && defined key_roam 
     1460               zsurfvar(:,:) = f2_pco2w(:,:) 
     1461#elif defined key_fabm 
     1462               pco2_3d(:,:,:) = fabm_get_interior_diagnostic_data(model, jp_fabm_o3pc) 
     1463               zsurfvar(:,:) = pco2_3d(:,:,1) 
     1464#else 
     1465               CALL ctl_stop( ' Trying to run spco2 observation operator', & 
     1466                  &           ' but no biogeochemical model appears to have been defined' ) 
     1467#endif 
     1468 
     1469            CASE DEFAULT 
     1470 
     1471               CALL ctl_stop( 'Unknown surface observation type '//TRIM(cobstypessurf(jtype))//' in dia_obs' ) 
     1472 
     1473            END SELECT 
     1474             
     1475            IF ( llog10 ) THEN 
     1476               ! Take the log10 where we can, otherwise exclude 
     1477               tiny = 1.0e-20 
     1478               WHERE(zsurfvar(:,:) > tiny .AND. zsurfvar(:,:) /= obfillflt ) 
     1479                  zsurfvar(:,:)  = LOG10(zsurfvar(:,:)) 
     1480               ELSEWHERE 
     1481                  zsurfvar(:,:)  = obfillflt 
     1482                  zsurfmask(:,:) = 0 
     1483               END WHERE 
    10661484            ENDIF 
     1485 
     1486            CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj,       & 
     1487               &               nit000, idaystp, zsurfvar, zsurfmask,    & 
     1488               &               n2dintsurf(jtype), llnightav(jtype),     & 
     1489               &               ravglamscl(jtype), ravgphiscl(jtype),     & 
     1490               &               lfpindegs(jtype) ) 
     1491 
    10671492         END DO 
     1493 
     1494         CALL wrk_dealloc( jpi, jpj, zsurfvar ) 
     1495         CALL wrk_dealloc( jpi, jpj, zsurfmask ) 
     1496#if defined key_fabm 
     1497         CALL wrk_dealloc( jpi, jpj, jpk, pco2_3d ) 
     1498#endif 
     1499 
    10681500      ENDIF 
    10691501 
    1070       !  - Sea surface anomaly 
    1071       IF ( ln_sla ) THEN 
    1072          DO jslaset = 1, nslasets 
    1073             CALL obs_sla_opt( sladatqc(jslaset),            & 
    1074                &              kstp, jpi, jpj, nit000, sshn, & 
    1075                &              tmask(:,:,1), n2dint ) 
    1076          END DO          
    1077       ENDIF 
    1078  
    1079       !  - Sea surface temperature 
    1080       IF ( ln_sst ) THEN 
    1081          DO jsstset = 1, nsstsets 
    1082             CALL obs_sst_opt( sstdatqc(jsstset),                & 
    1083                &              kstp, jpi, jpj, nit000, idaystp,  & 
    1084                &              tsn(:,:,1,jp_tem), tmask(:,:,1),  & 
    1085                &              n2dint, ld_sstnight(jsstset) ) 
    1086          END DO 
    1087       ENDIF 
    1088  
    1089       !  - Sea surface salinity 
    1090       IF ( ln_sss ) THEN 
    1091          IF(lwp) WRITE(numout,*) ' SSS currently not available' 
    1092       ENDIF 
    1093  
    1094 #if defined key_lim2 || defined key_lim3 
    1095       IF ( ln_seaice ) THEN 
    1096          DO jseaiceset = 1, nseaicesets 
    1097             CALL obs_seaice_opt( seaicedatqc(jseaiceset),      & 
    1098                &              kstp, jpi, jpj, nit000, 1.-frld, & 
    1099                &              tmask(:,:,1), n2dint ) 
    1100          END DO 
    1101       ENDIF       
    1102 #endif 
    1103  
    1104       !  - Velocity profiles 
    1105       IF ( ln_vel3d ) THEN 
    1106          DO jveloset = 1, nvelosets 
    1107            ! zonal component of velocity 
    1108            CALL obs_vel_opt( veldatqc(jveloset), kstp, jpi, jpj, jpk, & 
    1109               &              nit000, idaystp, un, vn, gdept_1d, umask, vmask, & 
    1110                              n1dint, n2dint, ld_velav(jveloset) ) 
    1111          END DO 
    1112       ENDIF 
    1113  
    1114 #if ! defined key_lim2 && ! defined key_lim3 
    1115       CALL wrk_dealloc(jpi,jpj,frld)  
    1116 #endif 
    1117  
    11181502   END SUBROUTINE dia_obs 
    1119    
    1120    SUBROUTINE dia_obs_wri  
     1503 
     1504   SUBROUTINE dia_obs_wri 
    11211505      !!---------------------------------------------------------------------- 
    11221506      !!                    ***  ROUTINE dia_obs_wri  *** 
    … …  
    11261510      !! ** Method  : Call observation diagnostic output routines 
    11271511      !! 
    1128       !! ** Action  :  
     1512      !! ** Action  : 
    11291513      !! 
    11301514      !! History : 
    … …  
    11341518      !!        !  07-03  (K. Mogensen) General handling of profiles 
    11351519      !!        !  08-09  (M. Valdivieso) Velocity component (U,V) profiles 
     1520      !!        !  15-08  (M. Martin) Combined writing for prof and surf types 
    11361521      !!---------------------------------------------------------------------- 
     1522      !! * Modules used 
     1523      USE obs_rot_vel          ! Rotation of velocities 
     1524 
    11371525      IMPLICIT NONE 
    11381526 
    11391527      !! * Local declarations 
    1140  
    1141       INTEGER :: jprofset                 ! Profile data set loop variable 
    1142       INTEGER :: jveloset                 ! Velocity data set loop variable 
    1143       INTEGER :: jslaset                  ! SLA data set loop variable 
    1144       INTEGER :: jsstset                  ! SST data set loop variable 
    1145       INTEGER :: jseaiceset               ! Sea Ice data set loop variable 
    1146       INTEGER :: jset 
    1147       INTEGER :: jfbini 
    1148       CHARACTER(LEN=20) :: datestr=" ",timestr=" " 
    1149       CHARACTER(LEN=10) :: cdtmp 
     1528      INTEGER :: jtype                    ! Data set loop variable 
     1529      INTEGER :: jo, jvar, jk 
     1530      REAL(wp), DIMENSION(:), ALLOCATABLE :: & 
     1531         & zu, & 
     1532         & zv 
     1533 
    11501534      !----------------------------------------------------------------------- 
    11511535      ! Depending on switches call various observation output routines 
    11521536      !----------------------------------------------------------------------- 
    11531537 
    1154       !  - Temperature/salinity profiles 
    1155  
    1156       IF( ln_t3d .OR. ln_s3d ) THEN 
    1157  
    1158          ! Copy data from prodatqc to profdata structures 
    1159          DO jprofset = 1, nprofsets 
    1160  
    1161             CALL obs_prof_decompress( prodatqc(jprofset), & 
    1162                  &                    profdata(jprofset), .TRUE., numout ) 
     1538      IF ( nproftypes > 0 ) THEN 
     1539 
     1540         DO jtype = 1, nproftypes 
     1541 
     1542            IF ( TRIM(cobstypesprof(jtype)) == 'vel' ) THEN 
     1543 
     1544               ! For velocity data, rotate the model velocities to N/S, E/W 
     1545               ! using the compressed data structure. 
     1546               ALLOCATE( & 
     1547                  & zu(profdataqc(jtype)%nvprot(1)), & 
     1548                  & zv(profdataqc(jtype)%nvprot(2))  & 
     1549                  & ) 
     1550 
     1551               CALL obs_rotvel( profdataqc(jtype), nn_2dint_default, zu, zv ) 
     1552 
     1553               DO jo = 1, profdataqc(jtype)%nprof 
     1554                  DO jvar = 1, 2 
     1555                     DO jk = profdataqc(jtype)%npvsta(jo,jvar), profdataqc(jtype)%npvend(jo,jvar) 
     1556 
     1557                        IF ( jvar == 1 ) THEN 
     1558                           profdataqc(jtype)%var(jvar)%vmod(jk) = zu(jk) 
     1559                        ELSE 
     1560                           profdataqc(jtype)%var(jvar)%vmod(jk) = zv(jk) 
     1561                        ENDIF 
     1562 
     1563                     END DO 
     1564                  END DO 
     1565               END DO 
     1566 
     1567               DEALLOCATE( zu ) 
     1568               DEALLOCATE( zv ) 
     1569 
     1570            END IF 
     1571 
     1572            CALL obs_prof_decompress( profdataqc(jtype), & 
     1573               &                      profdata(jtype), .TRUE., numout ) 
     1574 
     1575            CALL obs_wri_prof( profdata(jtype) ) 
    11631576 
    11641577         END DO 
    11651578 
    1166          ! Write the profiles. 
    1167  
    1168          jprofset = 0 
    1169  
    1170          ! ENACT insitu data 
    1171  
    1172          IF ( ln_ena ) THEN 
    1173             
    1174             jprofset = jprofset + 1 
    1175  
    1176             CALL obs_wri_p3d( 'enact', profdata(jprofset) ) 
    1177  
    1178          ENDIF 
    1179  
    1180          ! Coriolis insitu data 
    1181  
    1182          IF ( ln_cor ) THEN 
    1183              
    1184             jprofset = jprofset + 1 
    1185  
    1186             CALL obs_wri_p3d( 'corio', profdata(jprofset) ) 
    1187              
    1188          ENDIF 
    1189           
    1190          ! Feedback insitu data 
    1191  
    1192          IF ( ln_profb ) THEN 
    1193  
    1194             jfbini = jprofset + 1 
    1195  
    1196             DO jprofset = jfbini, nprofsets 
    1197                 
    1198                jset = jprofset - jfbini + 1 
    1199                WRITE(cdtmp,'(A,I2.2)')'profb_',jset 
    1200                CALL obs_wri_p3d( cdtmp, profdata(jprofset) ) 
    1201  
    1202             END DO 
    1203  
    1204          ENDIF 
    1205  
    12061579      ENDIF 
    12071580 
    1208       !  - Sea surface anomaly 
    1209       IF ( ln_sla ) THEN 
    1210  
    1211          ! Copy data from sladatqc to sladata structures 
    1212          DO jslaset = 1, nslasets 
    1213  
    1214               CALL obs_surf_decompress( sladatqc(jslaset), & 
    1215                  &                    sladata(jslaset), .TRUE., numout ) 
     1581      IF ( nsurftypes > 0 ) THEN 
     1582 
     1583         DO jtype = 1, nsurftypes 
     1584 
     1585            CALL obs_surf_decompress( surfdataqc(jtype), & 
     1586               &                      surfdata(jtype), .TRUE., numout ) 
     1587 
     1588            CALL obs_wri_surf( surfdata(jtype) ) 
    12161589 
    12171590         END DO 
    12181591 
    1219          jslaset = 0  
    1220  
    1221          ! Write the AVISO SLA data 
    1222  
    1223          IF ( ln_sladt ) THEN 
    1224              
    1225             jslaset = 1 
    1226             CALL obs_wri_sla( 'aviso_act', sladata(jslaset) ) 
    1227             jslaset = 2 
    1228             CALL obs_wri_sla( 'aviso_pas', sladata(jslaset) ) 
    1229  
    1230          ENDIF 
    1231  
    1232          IF ( ln_slafb ) THEN 
    1233              
    1234             jfbini = jslaset + 1 
    1235  
    1236             DO jslaset = jfbini, nslasets 
    1237                 
    1238                jset = jslaset - jfbini + 1 
    1239                WRITE(cdtmp,'(A,I2.2)')'slafb_',jset 
    1240                CALL obs_wri_sla( cdtmp, sladata(jslaset) ) 
    1241  
    1242             END DO 
    1243  
    1244          ENDIF 
    1245  
    1246       ENDIF 
    1247  
    1248       !  - Sea surface temperature 
    1249       IF ( ln_sst ) THEN 
    1250  
    1251          ! Copy data from sstdatqc to sstdata structures 
    1252          DO jsstset = 1, nsstsets 
    1253       
    1254               CALL obs_surf_decompress( sstdatqc(jsstset), & 
    1255                  &                    sstdata(jsstset), .TRUE., numout ) 
    1256  
    1257          END DO 
    1258  
    1259          jsstset = 0  
    1260  
    1261          ! Write the AVISO SST data 
    1262  
    1263          IF ( ln_reysst ) THEN 
    1264              
    1265             jsstset = jsstset + 1 
    1266             CALL obs_wri_sst( 'reynolds', sstdata(jsstset) ) 
    1267  
    1268          ENDIF 
    1269  
    1270          IF ( ln_ghrsst ) THEN 
    1271              
    1272             jsstset = jsstset + 1 
    1273             CALL obs_wri_sst( 'ghr', sstdata(jsstset) ) 
    1274  
    1275          ENDIF 
    1276  
    1277          IF ( ln_sstfb ) THEN 
    1278              
    1279             jfbini = jsstset + 1 
    1280  
    1281             DO jsstset = jfbini, nsstsets 
    1282                 
    1283                jset = jsstset - jfbini + 1 
    1284                WRITE(cdtmp,'(A,I2.2)')'sstfb_',jset 
    1285                CALL obs_wri_sst( cdtmp, sstdata(jsstset) ) 
    1286  
    1287             END DO 
    1288  
    1289          ENDIF 
    1290  
    1291       ENDIF 
    1292  
    1293       !  - Sea surface salinity 
    1294       IF ( ln_sss ) THEN 
    1295          IF(lwp) WRITE(numout,*) ' SSS currently not available' 
    1296       ENDIF 
    1297  
    1298       !  - Sea Ice Concentration 
    1299       IF ( ln_seaice ) THEN 
    1300  
    1301          ! Copy data from seaicedatqc to seaicedata structures 
    1302          DO jseaiceset = 1, nseaicesets 
    1303  
    1304               CALL obs_surf_decompress( seaicedatqc(jseaiceset), & 
    1305                  &                    seaicedata(jseaiceset), .TRUE., numout ) 
    1306  
    1307          END DO 
    1308  
    1309          ! Write the Sea Ice data 
    1310          DO jseaiceset = 1, nseaicesets 
    1311        
    1312             WRITE(cdtmp,'(A,I2.2)')'seaicefb_',jseaiceset 
    1313             CALL obs_wri_seaice( cdtmp, seaicedata(jseaiceset) ) 
    1314  
    1315          END DO 
    1316  
    1317       ENDIF 
    1318        
    1319       ! Velocity data 
    1320       IF( ln_vel3d ) THEN 
    1321  
    1322          ! Copy data from veldatqc to velodata structures 
    1323          DO jveloset = 1, nvelosets 
    1324  
    1325             CALL obs_prof_decompress( veldatqc(jveloset), & 
    1326                  &                    velodata(jveloset), .TRUE., numout ) 
    1327  
    1328          END DO 
    1329  
    1330          ! Write the profiles. 
    1331  
    1332          jveloset = 0 
    1333  
    1334          ! Daily averaged data 
    1335  
    1336          IF ( ln_velavcur ) THEN 
    1337              
    1338             jveloset = jveloset + 1 
    1339  
    1340             CALL obs_wri_vel( 'velavcurr', velodata(jveloset), n2dint ) 
    1341  
    1342          ENDIF 
    1343  
    1344          ! High frequency data 
    1345  
    1346          IF ( ln_velhrcur ) THEN 
    1347              
    1348             jveloset = jveloset + 1 
    1349  
    1350             CALL obs_wri_vel( 'velhrcurr', velodata(jveloset), n2dint ) 
    1351  
    1352          ENDIF 
    1353  
    1354          ! Daily averaged data 
    1355  
    1356          IF ( ln_velavadcp ) THEN 
    1357              
    1358             jveloset = jveloset + 1 
    1359  
    1360             CALL obs_wri_vel( 'velavadcp', velodata(jveloset), n2dint ) 
    1361  
    1362          ENDIF 
    1363  
    1364          ! High frequency data 
    1365  
    1366          IF ( ln_velhradcp ) THEN 
    1367              
    1368             jveloset = jveloset + 1 
    1369              
    1370             CALL obs_wri_vel( 'velhradcp', velodata(jveloset), n2dint ) 
    1371                 
    1372          ENDIF 
    1373  
    1374          ! Feedback velocity data 
    1375  
    1376          IF ( ln_velfb ) THEN 
    1377  
    1378             jfbini = jveloset + 1 
    1379  
    1380             DO jveloset = jfbini, nvelosets 
    1381                 
    1382                jset = jveloset - jfbini + 1 
    1383                WRITE(cdtmp,'(A,I2.2)')'velfb_',jset 
    1384                CALL obs_wri_vel( cdtmp, velodata(jveloset), n2dint ) 
    1385  
    1386             END DO 
    1387  
    1388          ENDIF 
    1389           
    13901592      ENDIF 
    13911593 
    … …  
    14051607      !! 
    14061608      !!---------------------------------------------------------------------- 
    1407       !! obs_grid deallocation 
     1609      ! obs_grid deallocation 
    14081610      CALL obs_grid_deallocate 
    14091611 
    1410       !! diaobs deallocation 
    1411       IF ( nprofsets > 0 ) THEN 
    1412           DEALLOCATE(ld_enact, & 
    1413                   &  profdata, & 
    1414                   &  prodatqc) 
    1415       END IF 
    1416       IF ( ln_sla ) THEN 
    1417           DEALLOCATE(sladata, & 
    1418                   &  sladatqc) 
    1419       END IF 
    1420       IF ( ln_seaice ) THEN 
    1421           DEALLOCATE(sladata, & 
    1422                   &  sladatqc) 
    1423       END IF 
    1424       IF ( ln_sst ) THEN 
    1425           DEALLOCATE(sstdata, & 
    1426                   &  sstdatqc) 
    1427       END IF 
    1428       IF ( ln_vel3d ) THEN 
    1429           DEALLOCATE(ld_velav, & 
    1430                   &  velodata, & 
    1431                   &  veldatqc) 
    1432       END IF 
     1612      ! diaobs deallocation 
     1613      IF ( nproftypes > 0 ) & 
     1614         &   DEALLOCATE( cobstypesprof, profdata, profdataqc, nvarsprof, nextrprof ) 
     1615 
     1616      IF ( nsurftypes > 0 ) & 
     1617         &   DEALLOCATE( cobstypessurf, surfdata, surfdataqc, nvarssurf, nextrsurf, & 
     1618         &               n2dintsurf, ravglamscl, ravgphiscl, lfpindegs, llnightav ) 
     1619 
    14331620   END SUBROUTINE dia_obs_dealloc 
    14341621 
    … …  
    14361623      !!---------------------------------------------------------------------- 
    14371624      !!                    ***  ROUTINE ini_date  *** 
    1438       !!           
    1439       !! ** Purpose : Get initial data in double precision YYYYMMDD.HHMMSS format 
    1440       !! 
    1441       !! ** Method  : Get initial data in double precision YYYYMMDD.HHMMSS format 
    1442       !! 
    1443       !! ** Action  : Get initial data in double precision YYYYMMDD.HHMMSS format 
     1625      !! 
     1626      !! ** Purpose : Get initial date in double precision YYYYMMDD.HHMMSS format 
     1627      !! 
     1628      !! ** Method  : Get initial date in double precision YYYYMMDD.HHMMSS format 
     1629      !! 
     1630      !! ** Action  : Get initial date in double precision YYYYMMDD.HHMMSS format 
    14441631      !! 
    14451632      !! History : 
    … …  
    14521639      USE phycst, ONLY : &            ! Physical constants 
    14531640         & rday 
    1454 !      USE daymod, ONLY : &            ! Time variables 
    1455 !         & nmonth_len            
    14561641      USE dom_oce, ONLY : &           ! Ocean space and time domain variables 
    14571642         & rdt 
    … …  
    14601645 
    14611646      !! * Arguments 
    1462       REAL(KIND=dp), INTENT(OUT) :: ddobsini                         ! Initial date in YYYYMMDD.HHMMSS 
     1647      REAL(dp), INTENT(OUT) :: ddobsini  ! Initial date in YYYYMMDD.HHMMSS 
    14631648 
    14641649      !! * Local declarations 
    … …  
    14681653      INTEGER :: ihou 
    14691654      INTEGER :: imin 
    1470       INTEGER :: imday         ! Number of days in month. 
    1471       REAL(KIND=wp) :: zdayfrc ! Fraction of day 
    1472  
    1473       INTEGER, DIMENSION(12) ::   imonth_len    !: length in days of the months of the current year 
    1474  
    1475       !!---------------------------------------------------------------------- 
    1476       !! Initial date initialization (year, month, day, hour, minute) 
    1477       !! (This assumes that the initial date is for 00z)) 
    1478       !!---------------------------------------------------------------------- 
     1655      INTEGER :: imday       ! Number of days in month. 
     1656      INTEGER, DIMENSION(12) :: & 
     1657         &       imonth_len  ! Length in days of the months of the current year 
     1658      REAL(wp) :: zdayfrc    ! Fraction of day 
     1659 
     1660      !---------------------------------------------------------------------- 
     1661      ! Initial date initialization (year, month, day, hour, minute) 
     1662      ! (This assumes that the initial date is for 00z)) 
     1663      !---------------------------------------------------------------------- 
    14791664      iyea =   ndate0 / 10000 
    14801665      imon = ( ndate0 - iyea * 10000 ) / 100 
    … …  
    14831668      imin = 0 
    14841669 
    1485       !!---------------------------------------------------------------------- 
    1486       !! Compute number of days + number of hours + min since initial time 
    1487       !!---------------------------------------------------------------------- 
     1670      !---------------------------------------------------------------------- 
     1671      ! Compute number of days + number of hours + min since initial time 
     1672      !---------------------------------------------------------------------- 
    14881673      iday = iday + ( nit000 -1 ) * rdt / rday 
    14891674      zdayfrc = ( nit000 -1 ) * rdt / rday 
    … …  
    14921677      imin = int( (zdayfrc * 24 - ihou) * 60 ) 
    14931678 
    1494       !!----------------------------------------------------------------------- 
    1495       !! Convert number of days (iday) into a real date 
    1496       !!---------------------------------------------------------------------- 
     1679      !----------------------------------------------------------------------- 
     1680      ! Convert number of days (iday) into a real date 
     1681      !---------------------------------------------------------------------- 
    14971682 
    14981683      CALL calc_month_len( iyea, imonth_len ) 
    1499        
     1684 
    15001685      DO WHILE ( iday > imonth_len(imon) ) 
    15011686         iday = iday - imonth_len(imon) 
    … …  
    15081693      END DO 
    15091694 
    1510       !!---------------------------------------------------------------------- 
    1511       !! Convert it into YYYYMMDD.HHMMSS format. 
    1512       !!---------------------------------------------------------------------- 
     1695      !---------------------------------------------------------------------- 
     1696      ! Convert it into YYYYMMDD.HHMMSS format. 
     1697      !---------------------------------------------------------------------- 
    15131698      ddobsini = iyea * 10000_dp + imon * 100_dp + & 
    15141699         &       iday + ihou * 0.01_dp + imin * 0.0001_dp 
    … …  
    15201705      !!---------------------------------------------------------------------- 
    15211706      !!                    ***  ROUTINE fin_date  *** 
    1522       !!           
    1523       !! ** Purpose : Get final data in double precision YYYYMMDD.HHMMSS format 
    1524       !! 
    1525       !! ** Method  : Get final data in double precision YYYYMMDD.HHMMSS format 
    1526       !! 
    1527       !! ** Action  : Get final data in double precision YYYYMMDD.HHMMSS format 
     1707      !! 
     1708      !! ** Purpose : Get final date in double precision YYYYMMDD.HHMMSS format 
     1709      !! 
     1710      !! ** Method  : Get final date in double precision YYYYMMDD.HHMMSS format 
     1711      !! 
     1712      !! ** Action  : Get final date in double precision YYYYMMDD.HHMMSS format 
    15281713      !! 
    15291714      !! History : 
    … …  
    15351720      USE phycst, ONLY : &            ! Physical constants 
    15361721         & rday 
    1537 !      USE daymod, ONLY : &            ! Time variables 
    1538 !         & nmonth_len                 
    15391722      USE dom_oce, ONLY : &           ! Ocean space and time domain variables 
    15401723         & rdt 
    … …  
    15431726 
    15441727      !! * Arguments 
    1545       REAL(KIND=dp), INTENT(OUT) :: ddobsfin                  ! Final date in YYYYMMDD.HHMMSS 
     1728      REAL(dp), INTENT(OUT) :: ddobsfin ! Final date in YYYYMMDD.HHMMSS 
    15461729 
    15471730      !! * Local declarations 
    … …  
    15511734      INTEGER :: ihou 
    15521735      INTEGER :: imin 
    1553       INTEGER :: imday         ! Number of days in month. 
    1554       REAL(KIND=wp) :: zdayfrc       ! Fraction of day 
    1555           
    1556       INTEGER, DIMENSION(12) ::   imonth_len    !: length in days of the months of the current year 
    1557              
     1736      INTEGER :: imday       ! Number of days in month. 
     1737      INTEGER, DIMENSION(12) :: & 
     1738         &       imonth_len  ! Length in days of the months of the current year 
     1739      REAL(wp) :: zdayfrc    ! Fraction of day 
     1740 
    15581741      !----------------------------------------------------------------------- 
    15591742      ! Initial date initialization (year, month, day, hour, minute) 
    … …  
    15651748      ihou = 0 
    15661749      imin = 0 
    1567        
     1750 
    15681751      !----------------------------------------------------------------------- 
    15691752      ! Compute number of days + number of hours + min since initial time 
    … …  
    15801763 
    15811764      CALL calc_month_len( iyea, imonth_len ) 
    1582        
     1765 
    15831766      DO WHILE ( iday > imonth_len(imon) ) 
    15841767         iday = iday - imonth_len(imon) 
    … …  
    15981781 
    15991782    END SUBROUTINE fin_date 
    1600      
     1783 
     1784    SUBROUTINE obs_settypefiles( ntypes, jpmaxnfiles, ifiles, cobstypes, cfiles ) 
     1785 
     1786       INTEGER, INTENT(IN) :: ntypes      ! Total number of obs types 
     1787       INTEGER, INTENT(IN) :: jpmaxnfiles ! Maximum number of files allowed for each type 
     1788       INTEGER, DIMENSION(ntypes), INTENT(OUT) :: & 
     1789          &                   ifiles      ! Out number of files for each type 
     1790       CHARACTER(len=8), DIMENSION(ntypes), INTENT(IN) :: & 
     1791          &                   cobstypes   ! List of obs types 
     1792       CHARACTER(len=128), DIMENSION(ntypes, jpmaxnfiles), INTENT(IN) :: & 
     1793          &                   cfiles      ! List of files for all types 
     1794 
     1795       !Local variables 
     1796       INTEGER :: jfile 
     1797       INTEGER :: jtype 
     1798 
     1799       DO jtype = 1, ntypes 
     1800 
     1801          ifiles(jtype) = 0 
     1802          DO jfile = 1, jpmaxnfiles 
     1803             IF ( trim(cfiles(jtype,jfile)) /= '' ) & 
     1804                       ifiles(jtype) = ifiles(jtype) + 1 
     1805          END DO 
     1806 
     1807          IF ( ifiles(jtype) == 0 ) THEN 
     1808               CALL ctl_stop( 'Logical for observation type '//TRIM(cobstypes(jtype))//   & 
     1809                  &           ' set to true but no files available to read' ) 
     1810          ENDIF 
     1811 
     1812          IF(lwp) THEN     
     1813             WRITE(numout,*) '             '//cobstypes(jtype)//' input observation file names:' 
     1814             DO jfile = 1, ifiles(jtype) 
     1815                WRITE(numout,*) '                '//TRIM(cfiles(jtype,jfile)) 
     1816             END DO 
     1817          ENDIF 
     1818 
     1819       END DO 
     1820 
     1821    END SUBROUTINE obs_settypefiles 
     1822 
     1823    SUBROUTINE obs_setinterpopts( ntypes, jtype, ctypein,             & 
     1824               &                  n2dint_default, n2dint_type,        & 
     1825               &                  ravglamscl_type, ravgphiscl_type,   & 
     1826               &                  lfp_indegs_type, lavnight_type,     & 
     1827               &                  n2dint, ravglamscl, ravgphiscl,     & 
     1828               &                  lfpindegs, lavnight ) 
     1829 
     1830       INTEGER, INTENT(IN)  :: ntypes             ! Total number of obs types 
     1831       INTEGER, INTENT(IN)  :: jtype              ! Index of the current type of obs 
     1832       INTEGER, INTENT(IN)  :: n2dint_default     ! Default option for interpolation type 
     1833       INTEGER, INTENT(IN)  :: n2dint_type        ! Option for interpolation type 
     1834       REAL(wp), INTENT(IN) :: & 
     1835          &                    ravglamscl_type, & !E/W diameter of obs footprint for this type 
     1836          &                    ravgphiscl_type    !N/S diameter of obs footprint for this type 
     1837       LOGICAL, INTENT(IN)  :: lfp_indegs_type    !T=> footprint in degrees, F=> in metres 
     1838       LOGICAL, INTENT(IN)  :: lavnight_type      !T=> obs represent night time average 
     1839       CHARACTER(len=8), INTENT(IN) :: ctypein  
     1840 
     1841       INTEGER, DIMENSION(ntypes), INTENT(INOUT) :: & 
     1842          &                    n2dint  
     1843       REAL(wp), DIMENSION(ntypes), INTENT(INOUT) :: & 
     1844          &                    ravglamscl, ravgphiscl 
     1845       LOGICAL, DIMENSION(ntypes), INTENT(INOUT) :: & 
     1846          &                    lfpindegs, lavnight 
     1847 
     1848       lavnight(jtype) = lavnight_type 
     1849 
     1850       IF ( (n2dint_type >= 0) .AND. (n2dint_type <= 6) ) THEN 
     1851          n2dint(jtype) = n2dint_type 
     1852       ELSE IF ( n2dint_type == -1 ) THEN 
     1853          n2dint(jtype) = n2dint_default 
     1854       ELSE 
     1855          CALL ctl_stop(' Choice of '//TRIM(ctypein)//' horizontal (2D) interpolation method', & 
     1856            &                    ' is not available') 
     1857       ENDIF 
     1858 
     1859       ! For averaging observation footprints set options for size of footprint  
     1860       IF ( (n2dint(jtype) > 4) .AND. (n2dint(jtype) <= 6) ) THEN 
     1861          IF ( ravglamscl_type > 0._wp ) THEN 
     1862             ravglamscl(jtype) = ravglamscl_type 
     1863          ELSE 
     1864             CALL ctl_stop( 'Incorrect value set for averaging footprint '// & 
     1865                            'scale (ravglamscl) for observation type '//TRIM(ctypein) )       
     1866          ENDIF 
     1867 
     1868          IF ( ravgphiscl_type > 0._wp ) THEN 
     1869             ravgphiscl(jtype) = ravgphiscl_type 
     1870          ELSE 
     1871             CALL ctl_stop( 'Incorrect value set for averaging footprint '// & 
     1872                            'scale (ravgphiscl) for observation type '//TRIM(ctypein) )       
     1873          ENDIF 
     1874 
     1875          lfpindegs(jtype) = lfp_indegs_type  
     1876 
     1877       ENDIF 
     1878 
     1879       ! Write out info  
     1880       IF(lwp) THEN 
     1881          IF ( n2dint(jtype) <= 4 ) THEN 
     1882             WRITE(numout,*) '             '//TRIM(ctypein)// & 
     1883                &            ' model counterparts will be interpolated horizontally' 
     1884          ELSE IF ( n2dint(jtype) <= 6 ) THEN 
     1885             WRITE(numout,*) '             '//TRIM(ctypein)// & 
     1886                &            ' model counterparts will be averaged horizontally' 
     1887             WRITE(numout,*) '             '//'    with E/W scale: ',ravglamscl(jtype) 
     1888             WRITE(numout,*) '             '//'    with N/S scale: ',ravgphiscl(jtype) 
     1889             IF ( lfpindegs(jtype) ) THEN 
     1890                 WRITE(numout,*) '             '//'    (in degrees)' 
     1891             ELSE 
     1892                 WRITE(numout,*) '             '//'    (in metres)' 
     1893             ENDIF 
     1894          ENDIF 
     1895       ENDIF 
     1896 
     1897    END SUBROUTINE obs_setinterpopts 
     1898 
    16011899END MODULE diaobs 

  • branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_grd_bruteforce.h90

    r2358 r11202  
    325325         CALL obs_mpp_max_integer( kobsj, kobs ) 
    326326      ELSE 
    327          CALL obs_mpp_find_obs_proc( kproc, kobsi, kobsj, kobs ) 
     327         CALL obs_mpp_find_obs_proc( kproc,kobs ) 
    328328      ENDIF 
    329329 

  • branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_grid.F90

    r4990 r11202  
    5252 
    5353   !! Default values 
    54    REAL, PUBLIC :: grid_search_res = 0.5    ! Resolution of grid 
     54   REAL, PUBLIC :: rn_gridsearchres = 0.5   ! Resolution of grid 
    5555   INTEGER, PRIVATE :: gsearch_nlons_def    ! Num of longitudes 
    5656   INTEGER, PRIVATE :: gsearch_nlats_def    ! Num of latitudes 
    … …  
    8383   LOGICAL, PUBLIC :: ln_grid_global         ! Use global distribution of observations 
    8484   CHARACTER(LEN=44), PUBLIC :: & 
    85       & grid_search_file    ! file name head for grid search lookup  
     85      & cn_gridsearchfile    ! file name head for grid search lookup  
    8686 
    8787   !!---------------------------------------------------------------------- 
    … …  
    613613         CALL obs_mpp_max_integer( kobsj, kobs ) 
    614614      ELSE 
    615          CALL obs_mpp_find_obs_proc( kproc, kobsi, kobsj, kobs ) 
     615         CALL obs_mpp_find_obs_proc( kproc, kobs ) 
    616616      ENDIF 
    617617 
    … …  
    690690          
    691691         IF(lwp) WRITE(numout,*) 
    692          IF(lwp) WRITE(numout,*)'Grid search resolution : ', grid_search_res 
    693           
    694          gsearch_nlons_def  = NINT( 360.0_wp / grid_search_res )  
    695          gsearch_nlats_def  = NINT( 180.0_wp / grid_search_res ) 
    696          gsearch_lonmin_def = -180.0_wp + 0.5_wp * grid_search_res 
    697          gsearch_latmin_def =  -90.0_wp + 0.5_wp * grid_search_res 
    698          gsearch_dlon_def   = grid_search_res 
    699          gsearch_dlat_def   = grid_search_res 
     692         IF(lwp) WRITE(numout,*)'Grid search resolution : ', rn_gridsearchres 
     693          
     694         gsearch_nlons_def  = NINT( 360.0_wp / rn_gridsearchres )  
     695         gsearch_nlats_def  = NINT( 180.0_wp / rn_gridsearchres ) 
     696         gsearch_lonmin_def = -180.0_wp + 0.5_wp * rn_gridsearchres 
     697         gsearch_latmin_def =  -90.0_wp + 0.5_wp * rn_gridsearchres 
     698         gsearch_dlon_def   = rn_gridsearchres 
     699         gsearch_dlat_def   = rn_gridsearchres 
    700700          
    701701         IF (lwp) THEN 
    … …  
    710710         IF ( ln_grid_global ) THEN 
    711711            WRITE(cfname, FMT="(A,'_',A)") & 
    712                &          TRIM(grid_search_file), 'global.nc' 
     712               &          TRIM(cn_gridsearchfile), 'global.nc' 
    713713         ELSE 
    714714            WRITE(cfname, FMT="(A,'_',I4.4,'of',I4.4,'by',I4.4,'.nc')") & 
    715                &          TRIM(grid_search_file), nproc, jpni, jpnj 
     715               &          TRIM(cn_gridsearchfile), nproc, jpni, jpnj 
    716716         ENDIF 
    717717 

  • branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_inter_sup.F90

    r3294 r11202  
    3535CONTAINS 
    3636 
    37    SUBROUTINE obs_int_comm_3d( kptsi, kptsj, kobs, kpk, kgrdi, kgrdj, & 
     37   SUBROUTINE obs_int_comm_3d( kptsi, kptsj, kobs, kpi, kpj, kpk, kgrdi, kgrdj, & 
    3838      &                        pval, pgval, kproc ) 
    3939      !!---------------------------------------------------------------------- 
    … …  
    5757      INTEGER, INTENT(IN) :: kptsj     ! Number of j horizontal points per stencil 
    5858      INTEGER, INTENT(IN) :: kobs      ! Local number of observations 
     59      INTEGER, INTENT(IN) :: kpi       ! Number of points in i direction 
     60      INTEGER, INTENT(IN) :: kpj       ! Number of points in j direction 
    5961      INTEGER, INTENT(IN) :: kpk       ! Number of levels 
    6062      INTEGER, DIMENSION(kptsi,kptsj,kobs), INTENT(IN) :: & 
    … …  
    6365      INTEGER, OPTIONAL, DIMENSION(kptsi,kptsj,kobs), INTENT(IN) :: & 
    6466         & kproc            ! Precomputed processor for each i,j,iobs points 
    65       REAL(KIND=wp), DIMENSION(jpi,jpj,kpk), INTENT(IN) ::& 
     67      REAL(KIND=wp), DIMENSION(kpi,kpj,kpk), INTENT(IN) ::& 
    6668         & pval             ! Local 3D array to extract data from 
    6769      REAL(KIND=wp), DIMENSION(kptsi,kptsj,kpk,kobs), INTENT(OUT) ::& 
    … …  
    7375         IF (PRESENT(kproc)) THEN 
    7476 
    75             CALL obs_int_comm_3d_global( kptsi, kptsj, kobs, kpk, kgrdi, & 
     77            CALL obs_int_comm_3d_global( kptsi, kptsj, kobs, kpi, kpj, kpk, kgrdi, & 
    7678               &                         kgrdj, pval, pgval, kproc=kproc ) 
    7779 
    7880         ELSE 
    7981 
    80             CALL obs_int_comm_3d_global( kptsi, kptsj, kobs, kpk, kgrdi, & 
     82            CALL obs_int_comm_3d_global( kptsi, kptsj, kobs, kpi, kpj, kpk, kgrdi, & 
    8183               &                         kgrdj, pval, pgval ) 
    8284 
    … …  
    8587      ELSE 
    8688 
    87          CALL obs_int_comm_3d_local( kptsi, kptsj, kobs, kpk, kgrdi, kgrdj, & 
     89         CALL obs_int_comm_3d_local( kptsi, kptsj, kobs, kpi, kpj, kpk, kgrdi, kgrdj, & 
    8890            &                        pval, pgval ) 
    8991 
    … …  
    9294   END SUBROUTINE obs_int_comm_3d 
    9395 
    94    SUBROUTINE obs_int_comm_2d( kptsi, kptsj, kobs, kgrdi, kgrdj, pval, pgval, & 
     96   SUBROUTINE obs_int_comm_2d( kptsi, kptsj, kobs, kpi, kpj, kgrdi, kgrdj, pval, pgval, & 
    9597      &                        kproc ) 
    9698      !!---------------------------------------------------------------------- 
    … …  
    111113      INTEGER, INTENT(IN) :: kptsj        ! Number of j horizontal points per stencil 
    112114      INTEGER, INTENT(IN) :: kobs          ! Local number of observations 
     115      INTEGER, INTENT(IN) :: kpi          ! Number of model grid points in i direction 
     116      INTEGER, INTENT(IN) :: kpj          ! Number of model grid points in j direction 
    113117      INTEGER, DIMENSION(kptsi,kptsj,kobs), INTENT(IN) :: & 
    114118         & kgrdi, &         ! i,j indicies for each stencil 
    … …  
    116120      INTEGER, OPTIONAL, DIMENSION(kptsi,kptsj,kobs), INTENT(IN) :: & 
    117121         & kproc            ! Precomputed processor for each i,j,iobs points 
    118       REAL(KIND=wp), DIMENSION(jpi,jpj), INTENT(IN) ::& 
     122      REAL(KIND=wp), DIMENSION(kpi,kpj), INTENT(IN) ::& 
    119123         & pval             ! Local 3D array to extra data from 
    120124      REAL(KIND=wp), DIMENSION(kptsi,kptsj,kobs), INTENT(OUT) ::& 
    … …  
    136140      IF (PRESENT(kproc)) THEN 
    137141 
    138          CALL obs_int_comm_3d( kptsi, kptsj, kobs, 1, kgrdi, kgrdj, zval, & 
     142         CALL obs_int_comm_3d( kptsi, kptsj, kobs, kpi, kpj, 1, kgrdi, kgrdj, zval, & 
    139143            &                  zgval, kproc=kproc ) 
    140144      ELSE 
    141145 
    142          CALL obs_int_comm_3d( kptsi, kptsj, kobs, 1, kgrdi, kgrdj, zval, & 
     146         CALL obs_int_comm_3d( kptsi, kptsj, kobs, kpi, kpj, 1, kgrdi, kgrdj, zval, & 
    143147            &                  zgval ) 
    144148 
    … …  
    154158   END SUBROUTINE obs_int_comm_2d 
    155159 
    156    SUBROUTINE obs_int_comm_3d_global( kptsi, kptsj, kobs, kpk, kgrdi, kgrdj, & 
     160   SUBROUTINE obs_int_comm_3d_global( kptsi, kptsj, kobs, kpi, kpj, kpk, kgrdi, kgrdj, & 
    157161      &                               pval, pgval, kproc ) 
    158162      !!---------------------------------------------------------------------- 
    … …  
    174178      INTEGER, INTENT(IN) :: kptsj     ! Number of j horizontal points per stencil 
    175179      INTEGER, INTENT(IN) :: kobs      ! Local number of observations 
     180      INTEGER, INTENT(IN) :: kpi       ! Number of model points in i direction 
     181      INTEGER, INTENT(IN) :: kpj       ! Number of model points in j direction 
    176182      INTEGER, INTENT(IN) :: kpk       ! Number of levels 
    177183      INTEGER, DIMENSION(kptsi,kptsj,kobs), INTENT(IN) :: & 
    … …  
    180186      INTEGER, OPTIONAL, DIMENSION(kptsi,kptsj,kobs), INTENT(IN) :: & 
    181187         & kproc            ! Precomputed processor for each i,j,iobs points 
    182       REAL(KIND=wp), DIMENSION(jpi,jpj,kpk), INTENT(IN) ::& 
     188      REAL(KIND=wp), DIMENSION(kpi,kpj,kpk), INTENT(IN) ::& 
    183189         & pval             ! Local 3D array to extract data from 
    184190      REAL(KIND=wp), DIMENSION(kptsi,kptsj,kpk,kobs), INTENT(OUT) ::& 
    … …  
    207213 
    208214      ! Check valid points 
    209        
     215 
    210216      IF ( ( MAXVAL(kgrdi) > jpiglo ) .OR. ( MINVAL(kgrdi) < 1 ) .OR. & 
    211217         & ( MAXVAL(kgrdj) > jpjglo ) .OR. ( MINVAL(kgrdj) < 1 ) ) THEN 
    212           
     218 
    213219         CALL ctl_stop( 'Error in obs_int_comm_3d_global', & 
    214220            &           'Point outside global domain' ) 
    215           
     221 
    216222      ENDIF 
    217223 
    … …  
    323329   END SUBROUTINE obs_int_comm_3d_global 
    324330    
    325    SUBROUTINE obs_int_comm_3d_local( kptsi, kptsj, kobs, kpk, kgrdi, kgrdj, & 
     331   SUBROUTINE obs_int_comm_3d_local( kptsi, kptsj, kobs, kpi, kpj, kpk, kgrdi, kgrdj, & 
    326332      &                              pval, pgval ) 
    327333      !!---------------------------------------------------------------------- 
    … …  
    343349      INTEGER, INTENT(IN) :: kptsj        ! Number of j horizontal points per stencil 
    344350      INTEGER, INTENT(IN) :: kobs         ! Local number of observations 
     351      INTEGER, INTENT(IN) :: kpi          ! Number of model points in i direction 
     352      INTEGER, INTENT(IN) :: kpj          ! Number of model points in j direction 
    345353      INTEGER, INTENT(IN) :: kpk          ! Number of levels 
    346354      INTEGER, DIMENSION(kptsi,kptsj,kobs), INTENT(IN) :: & 
    347355         & kgrdi, &         ! i,j indicies for each stencil 
    348356         & kgrdj 
    349       REAL(KIND=wp), DIMENSION(jpi,jpj,kpk), INTENT(IN) ::& 
     357      REAL(KIND=wp), DIMENSION(kpi,kpj,kpk), INTENT(IN) ::& 
    350358         & pval             ! Local 3D array to extract data from 
    351359      REAL(KIND=wp), DIMENSION(kptsi,kptsj,kpk,kobs), INTENT(OUT) ::& 

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

    r2513 r11202  
    77   !!             -   ! 2006-05  (K. Mogensen)  Reformatted 
    88   !!             -   ! 2008-01  (K. Mogensen)  add mpp_global_max 
     9   !!            3.6  ! 2015-01  (J. Waters) obs_mpp_find_obs_proc  
     10   !!                            rewritten to avoid global arrays 
    911   !!---------------------------------------------------------------------- 
    1012#  define mpivar mpi_double_precision 
    … …  
    1214   !! obs_mpp_bcast_integer : Broadcast an integer array from a processor to all processors 
    1315   !! obs_mpp_max_integer   : Find maximum on all processors of each value in an integer on all processors 
    14    !! obs_mpp_find_obs_proc : Find processors which should hold the observations 
     16   !! obs_mpp_find_obs_proc : Find processors which should hold the observations, avoiding global arrays 
    1517   !! obs_mpp_sum_integers  : Sum an integer array from all processors 
    1618   !! obs_mpp_sum_integer   : Sum an integer from all processors 
    … …  
    9698      ! 
    9799      INTEGER :: ierr  
    98       INTEGER, DIMENSION(kno) ::   ivals 
    99       ! 
    100 INCLUDE 'mpif.h' 
    101       !!---------------------------------------------------------------------- 
     100      INTEGER, DIMENSION(:), ALLOCATABLE ::   ivals 
     101      ! 
     102INCLUDE 'mpif.h' 
     103      !!---------------------------------------------------------------------- 
     104 
     105      ALLOCATE( ivals(kno) ) 
    102106 
    103107      ! Call the MPI library to find the maximum across processors 
    … …  
    105109         &                mpi_max, mpi_comm_opa, ierr ) 
    106110      kvals(:) = ivals(:) 
     111 
     112      DEALLOCATE( ivals ) 
    107113#else 
    108114      ! no MPI: empty routine 
    … …  
    111117 
    112118 
    113    SUBROUTINE obs_mpp_find_obs_proc( kobsp, kobsi, kobsj, kno ) 
    114       !!---------------------------------------------------------------------- 
    115       !!               ***  ROUTINE obs_mpp_find_obs_proc *** 
    116       !!           
    117       !! ** Purpose : From the array kobsp containing the results of the grid 
     119   SUBROUTINE obs_mpp_find_obs_proc( kobsp,kno ) 
     120      !!---------------------------------------------------------------------- 
     121      !!               ***  ROUTINE obs_mpp_find_obs_proc  *** 
     122      !!          
     123      !! ** Purpose : From the array kobsp containing the results of the 
    118124      !!              grid search on each processor the processor return a 
    119125      !!              decision of which processors should hold the observation. 
    120126      !! 
    121       !! ** Method  : A temporary 2D array holding all the decisions is 
    122       !!              constructed using mpi_allgather on each processor. 
    123       !!              If more than one processor has found the observation 
    124       !!              with the observation in the inner domain gets it 
    125       !! 
    126       !! ** Action  : This does only work for MPI.  
     127      !! ** Method  : Synchronize the processor number for each obs using 
     128      !!              obs_mpp_max_integer. If an observation exists on two  
     129      !!              processors it will be allocated to the lower numbered 
     130      !!              processor. 
     131      !! 
     132      !! ** Action  : This does only work for MPI. 
    127133      !!              It does not work for SHMEM. 
    128134      !! 
    … …  
    130136      !!---------------------------------------------------------------------- 
    131137      INTEGER                , INTENT(in   ) ::   kno 
    132       INTEGER, DIMENSION(kno), INTENT(in   ) ::   kobsi, kobsj 
    133138      INTEGER, DIMENSION(kno), INTENT(inout) ::   kobsp 
    134139      ! 
    135140#if defined key_mpp_mpi 
    136141      ! 
    137       INTEGER :: ji 
    138       INTEGER :: jj 
    139       INTEGER :: size 
    140       INTEGER :: ierr 
    141       INTEGER :: iobsip 
    142       INTEGER :: iobsjp 
    143       INTEGER :: num_sus_obs 
    144       INTEGER, DIMENSION(kno) ::   iobsig, iobsjg 
    145       INTEGER, ALLOCATABLE, DIMENSION(:,:) ::   iobsp, iobsi, iobsj 
    146       !! 
    147 INCLUDE 'mpif.h' 
    148       !!---------------------------------------------------------------------- 
    149  
    150       !----------------------------------------------------------------------- 
    151       ! Call the MPI library to find the maximum accross processors 
    152       !----------------------------------------------------------------------- 
    153       CALL mpi_comm_size( mpi_comm_opa, size, ierr ) 
    154       !----------------------------------------------------------------------- 
    155       ! Convert local grids points to global grid points 
    156       !----------------------------------------------------------------------- 
     142      ! 
     143      INTEGER :: ji, isum 
     144      INTEGER, DIMENSION(:), ALLOCATABLE ::   iobsp 
     145      !! 
     146      !! 
     147 
     148      ALLOCATE( iobsp(kno) ) 
     149 
     150      iobsp(:)=kobsp(:) 
     151 
     152      WHERE( iobsp(:) == -1 ) 
     153         iobsp(:) = 9999999 
     154      END WHERE 
     155 
     156      iobsp(:)=-1*iobsp(:) 
     157 
     158      CALL obs_mpp_max_integer( iobsp, kno ) 
     159 
     160      kobsp(:)=-1*iobsp(:) 
     161 
     162      isum=0 
    157163      DO ji = 1, kno 
    158          IF ( ( kobsi(ji) >= 1 ) .AND. ( kobsi(ji) <= jpi ) .AND. & 
    159             & ( kobsj(ji) >= 1 ) .AND. ( kobsj(ji) <= jpj ) ) THEN 
    160             iobsig(ji) = mig( kobsi(ji) ) 
    161             iobsjg(ji) = mjg( kobsj(ji) ) 
    162          ELSE 
    163             iobsig(ji) = -1 
    164             iobsjg(ji) = -1 
     164         IF ( kobsp(ji) == 9999999 ) THEN 
     165            isum=isum+1 
     166            kobsp(ji)=-1 
    165167         ENDIF 
    166       END DO 
    167       !----------------------------------------------------------------------- 
    168       ! Get the decisions from all processors 
    169       !----------------------------------------------------------------------- 
    170       ALLOCATE( iobsp(kno,size) ) 
    171       ALLOCATE( iobsi(kno,size) ) 
    172       ALLOCATE( iobsj(kno,size) ) 
    173       CALL mpi_allgather( kobsp, kno, mpi_integer, & 
    174          &                iobsp, kno, mpi_integer, & 
    175          &                mpi_comm_opa, ierr ) 
    176       CALL mpi_allgather( iobsig, kno, mpi_integer, & 
    177          &                iobsi, kno, mpi_integer, & 
    178          &                mpi_comm_opa, ierr ) 
    179       CALL mpi_allgather( iobsjg, kno, mpi_integer, & 
    180          &                iobsj, kno, mpi_integer, & 
    181          &                mpi_comm_opa, ierr ) 
    182  
    183       !----------------------------------------------------------------------- 
    184       ! Find the processor with observations from the lowest processor  
    185       ! number among processors holding the observation. 
    186       !----------------------------------------------------------------------- 
    187       kobsp(:) = -1 
    188       num_sus_obs = 0 
    189       DO ji = 1, kno 
    190          DO jj = 1, size 
    191             IF ( ( kobsp(ji) == -1 ) .AND. ( iobsp(ji,jj) /= -1 ) ) THEN 
    192                kobsp(ji) = iobsp(ji,jj) 
    193                iobsip = iobsi(ji,jj) 
    194                iobsjp = iobsj(ji,jj) 
    195             ENDIF 
    196             IF ( ( kobsp(ji) /= -1 ) .AND. ( iobsp(ji,jj) /= -1 ) ) THEN 
    197                IF ( ( iobsip /= iobsi(ji,jj) ) .OR. & 
    198                   & ( iobsjp /= iobsj(ji,jj) ) ) THEN 
    199                   IF ( ( kobsp(ji) < 1000000 ) .AND. & 
    200                      & ( iobsp(ji,jj) < 1000000 ) ) THEN 
    201                      num_sus_obs=num_sus_obs+1 
    202                   ENDIF 
    203                ENDIF 
    204                IF ( mppmap(iobsip,iobsjp) /= ( kobsp(ji)+1 ) ) THEN 
    205                   IF ( ( iobsi(ji,jj) /= -1 ) .AND. & 
    206                      & ( iobsj(ji,jj) /= -1 ) ) THEN 
    207                      IF ((mppmap(iobsi(ji,jj),iobsj(ji,jj)) == (iobsp(ji,jj)+1))& 
    208                         & .OR. ( iobsp(ji,jj) < kobsp(ji) ) ) THEN 
    209                         kobsp(ji) = iobsp(ji,jj) 
    210                         iobsip = iobsi(ji,jj) 
    211                         iobsjp = iobsj(ji,jj) 
    212                      ENDIF 
    213                   ENDIF 
    214                ENDIF 
    215             ENDIF 
    216          END DO 
    217       END DO 
    218       IF (lwp) WRITE(numout,*) 'Number of suspicious observations: ',num_sus_obs 
    219  
    220       DEALLOCATE( iobsj ) 
    221       DEALLOCATE( iobsi ) 
     168      ENDDO 
     169 
     170 
     171      IF ( isum > 0 ) THEN 
     172         IF (lwp) WRITE(numout,*) isum, ' observations failed the grid search.' 
     173         IF (lwp) WRITE(numout,*)'If ln_grid_search_lookup=.TRUE., try reducing grid_search_res' 
     174      ENDIF 
     175 
    222176      DEALLOCATE( iobsp ) 
     177 
    223178#else 
    224179      ! no MPI: empty routine 
    225 #endif 
    226       ! 
     180#endif      
     181       
    227182   END SUBROUTINE obs_mpp_find_obs_proc 
    228183 

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

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