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

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

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/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 
Note: See TracChangeset for help on using the changeset viewer.