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

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

Ignore:
Timestamp:
2017-05-02T13:21:57+02:00 (7 years ago)
Author:
jwhile
Message:

This version of the code seems to work correctly after some bug fixes

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r7960 r7992  
    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_sla_fp_indegs    !: T=> SLA obs footprint size specified in degrees, F=> in metres 
     51   LOGICAL :: ln_sst_fp_indegs    !: T=> SST obs footprint size specified in degrees, F=> in metres 
     52   LOGICAL :: ln_sss_fp_indegs    !: T=> SSS obs footprint size specified in degrees, F=> in metres 
     53   LOGICAL :: ln_sic_fp_indegs    !: T=> sea-ice obs footprint size specified in degrees, F=> in metres 
     54 
     55   REAL(wp) :: rn_sla_avglamscl !: E/W diameter of SLA observation footprint (metres) 
     56   REAL(wp) :: rn_sla_avgphiscl !: N/S diameter of SLA observation footprint (metres) 
     57   REAL(wp) :: rn_sst_avglamscl !: E/W diameter of SST observation footprint (metres) 
     58   REAL(wp) :: rn_sst_avgphiscl !: N/S diameter of SST observation footprint (metres) 
     59   REAL(wp) :: rn_sss_avglamscl !: E/W diameter of SSS observation footprint (metres) 
     60   REAL(wp) :: rn_sss_avgphiscl !: N/S diameter of SSS observation footprint (metres) 
     61   REAL(wp) :: rn_sic_avglamscl !: E/W diameter of sea-ice observation footprint (metres) 
     62   REAL(wp) :: rn_sic_avgphiscl !: N/S diameter of sea-ice observation footprint (metres) 
     63 
     64   INTEGER :: nn_1dint       !: Vertical interpolation method 
     65   INTEGER :: nn_2dint       !: Default horizontal interpolation method 
     66   INTEGER :: nn_2dint_sla   !: SLA horizontal interpolation method  
     67   INTEGER :: nn_2dint_sst   !: SST horizontal interpolation method  
     68   INTEGER :: nn_2dint_sss   !: SSS horizontal interpolation method  
     69   INTEGER :: nn_2dint_sic   !: Seaice horizontal interpolation method  
     70  
    9671   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? 
     72      & nn_profdavtypes      !: Profile data types representing a daily average 
     73   INTEGER :: nproftypes     !: Number of profile obs types 
     74   INTEGER :: nsurftypes     !: Number of surface obs types 
     75   INTEGER, DIMENSION(:), ALLOCATABLE :: & 
     76      & nvarsprof, &         !: Number of profile variables 
     77      & nvarssurf            !: Number of surface variables 
     78   INTEGER, DIMENSION(:), ALLOCATABLE :: & 
     79      & nextrprof, &         !: Number of profile extra variables 
     80      & nextrsurf            !: Number of surface extra variables 
     81   INTEGER, DIMENSION(:), ALLOCATABLE :: & 
     82      & n2dintsurf           !: Interpolation option for surface variables 
     83   REAL(wp), DIMENSION(:), ALLOCATABLE :: & 
     84      & ravglamscl, &        !: E/W diameter of averaging footprint for surface variables 
     85      & ravgphiscl           !: N/S diameter of averaging footprint for surface variables 
    10786   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 
     87      & lfpindegs, &         !: T=> surface obs footprint size specified in degrees, F=> in metres 
     88      & llnightav            !: Logical for calculating night-time averages 
     89 
     90   TYPE(obs_surf), PUBLIC, POINTER, DIMENSION(:) :: & 
     91      & surfdata, &          !: Initial surface data 
     92      & surfdataqc           !: Surface data after quality control 
     93   TYPE(obs_prof), PUBLIC, POINTER, DIMENSION(:) :: & 
     94      & profdata, &          !: Initial profile data 
     95      & profdataqc           !: Profile data after quality control 
     96 
     97   CHARACTER(len=6), PUBLIC, DIMENSION(:), ALLOCATABLE :: & 
     98      & cobstypesprof, &     !: Profile obs types 
     99      & cobstypessurf        !: Surface obs types 
    113100 
    114101   !!---------------------------------------------------------------------- 
    … …  
    118105   !!---------------------------------------------------------------------- 
    119106 
     107   !! * Substitutions  
     108#  include "domzgr_substitute.h90" 
    120109CONTAINS 
    121110 
    … …  
    135124      !!        !  06-10  (A. Weaver) Cleaning and add controls 
    136125      !!        !  07-03  (K. Mogensen) General handling of profiles 
     126      !!        !  14-08  (J.While) Incorporated SST bias correction 
     127      !!        !  15-02  (M. Martin) Simplification of namelist and code 
    137128      !!---------------------------------------------------------------------- 
    138129 
    … …  
    140131 
    141132      !! * 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,    & 
     133      INTEGER, PARAMETER :: & 
     134         & jpmaxnfiles = 1000    ! Maximum number of files for each obs type 
     135      INTEGER, DIMENSION(:), ALLOCATABLE :: & 
     136         & ifilesprof, &         ! Number of profile files 
     137         & ifilessurf            ! Number of surface files 
     138      INTEGER :: ios             ! Local integer output status for namelist read 
     139      INTEGER :: jtype           ! Counter for obs types 
     140      INTEGER :: jvar            ! Counter for variables 
     141      INTEGER :: jfile           ! Counter for files 
     142      INTEGER :: jnumsstbias     ! Number of SST bias files to read and apply 
     143 
     144      CHARACTER(len=128), DIMENSION(jpmaxnfiles) :: & 
     145         & cn_profbfiles,    &   ! T/S profile input filenames 
     146         & cn_sstfbfiles,    &   ! Sea surface temperature input filenames 
     147         & cn_slafbfiles,    &   ! Sea level anomaly input filenames 
     148         & cn_sicfbfiles,    &   ! Seaice concentration input filenames 
     149         & cn_velfbfiles,    &   ! Velocity profile input filenames 
     150         & cn_sssfbfiles,    &   ! Sea surface salinity input filenames 
     151         & cn_logchlfbfiles, &   ! Log(Chl) input filenames 
     152         & cn_spmfbfiles,    &   ! Sediment input filenames 
     153         & cn_fco2fbfiles,   &   ! fco2 input filenames 
     154         & cn_pco2fbfiles,   &   ! pco2 input filenames 
     155         & cn_sstbiasfiles       ! SST bias input filenames 
     156 
     157      CHARACTER(LEN=128) :: & 
     158         & cn_altbiasfile        ! Altimeter bias input filename 
     159 
     160 
     161      LOGICAL :: ln_t3d          ! Logical switch for temperature profiles 
     162      LOGICAL :: ln_s3d          ! Logical switch for salinity profiles 
     163      LOGICAL :: ln_sla          ! Logical switch for sea level anomalies  
     164      LOGICAL :: ln_sst          ! Logical switch for sea surface temperature 
     165      LOGICAL :: ln_sic          ! Logical switch for sea ice concentration 
     166      LOGICAL :: ln_sss          ! Logical switch for sea surface salinity obs 
     167      LOGICAL :: ln_vel3d        ! Logical switch for velocity (u,v) obs 
     168      LOGICAL :: ln_logchl       ! Logical switch for log(Chl) obs 
     169      LOGICAL :: ln_spm          ! Logical switch for sediment obs 
     170      LOGICAL :: ln_fco2         ! Logical switch for fco2 obs 
     171      LOGICAL :: ln_pco2         ! Logical switch for pco2 obs 
     172      LOGICAL :: ln_nea          ! Logical switch to remove obs near land 
     173      LOGICAL :: ln_altbias      ! Logical switch for altimeter bias 
     174      LOGICAL :: ln_sstbias      ! Logical switch for bias correction of SST 
     175      LOGICAL :: ln_ignmis       ! Logical switch for ignoring missing files 
     176      LOGICAL :: ln_s_at_t       ! Logical switch to compute model S at T obs 
     177      LOGICAL :: ln_bound_reject ! Logical switch for rejecting obs near the boundary 
     178 
     179      REAL(dp) :: rn_dobsini     ! Obs window start date YYYYMMDD.HHMMSS 
     180      REAL(dp) :: rn_dobsend     ! Obs window end date   YYYYMMDD.HHMMSS 
     181 
     182      CHARACTER(len=128), DIMENSION(:,:), ALLOCATABLE :: & 
     183         & clproffiles, &        ! Profile filenames 
     184         & clsurffiles           ! Surface filenames 
     185 
     186      LOGICAL :: llvar1          ! Logical for profile variable 1 
     187      LOGICAL :: llvar2          ! Logical for profile variable 1 
     188 
     189      REAL(wp), POINTER, DIMENSION(:,:) :: & 
     190         & zglam1, &             ! Model longitudes for profile variable 1 
     191         & zglam2                ! Model longitudes for profile variable 2 
     192      REAL(wp), POINTER, DIMENSION(:,:) :: & 
     193         & zgphi1, &             ! Model latitudes for profile variable 1 
     194         & zgphi2                ! Model latitudes for profile variable 2 
     195      REAL(wp), POINTER, DIMENSION(:,:,:) :: & 
     196         & zmask1, &             ! Model land/sea mask associated with variable 1 
     197         & zmask2                ! Model land/sea mask associated with variable 2 
     198 
     199 
     200      NAMELIST/namobs/ln_diaobs, ln_t3d, ln_s3d, ln_sla,              & 
     201         &            ln_sst, ln_sic, ln_sss, ln_vel3d,               & 
     202         &            ln_logchl, ln_spm, ln_fco2, ln_pco2,            & 
     203         &            ln_altbias, ln_sstbias, ln_nea,                 & 
     204         &            ln_grid_global, ln_grid_search_lookup,          & 
     205         &            ln_ignmis, ln_s_at_t, ln_bound_reject,          & 
    172206         &            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 
     207         &            ln_sla_fp_indegs, ln_sst_fp_indegs,             & 
     208         &            ln_sss_fp_indegs, ln_sic_fp_indegs,             & 
     209         &            cn_profbfiles, cn_slafbfiles,                   & 
     210         &            cn_sstfbfiles, cn_sicfbfiles,                   & 
     211         &            cn_velfbfiles, cn_sssfbfiles,                   & 
     212         &            cn_logchlfbfiles, cn_spmfbfiles,                & 
     213         &            cn_fco2fbfiles, cn_pco2fbfiles,                 & 
     214         &            cn_sstbiasfiles, cn_altbiasfile,                & 
     215         &            cn_gridsearchfile, rn_gridsearchres,            & 
     216         &            rn_dobsini, rn_dobsend,                         & 
     217         &            rn_sla_avglamscl, rn_sla_avgphiscl,             & 
     218         &            rn_sst_avglamscl, rn_sst_avgphiscl,             & 
     219         &            rn_sss_avglamscl, rn_sss_avgphiscl,             & 
     220         &            rn_sic_avglamscl, rn_sic_avgphiscl,             & 
     221         &            nn_1dint, nn_2dint,                             & 
     222         &            nn_2dint_sla, nn_2dint_sst,                     & 
     223         &            nn_2dint_sss, nn_2dint_sic,                     & 
     224         &            nn_msshc, rn_mdtcorr, rn_mdtcutoff,             & 
     225         &            nn_profdavtypes 
     226 
     227      CALL wrk_alloc( jpi, jpj, zglam1 ) 
     228      CALL wrk_alloc( jpi, jpj, zglam2 ) 
     229      CALL wrk_alloc( jpi, jpj, zgphi1 ) 
     230      CALL wrk_alloc( jpi, jpj, zgphi2 ) 
     231      CALL wrk_alloc( jpi, jpj, jpk, zmask1 ) 
     232      CALL wrk_alloc( jpi, jpj, jpk, zmask2 ) 
    205233 
    206234      !----------------------------------------------------------------------- 
    … …  
    208236      !----------------------------------------------------------------------- 
    209237 
    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 
     238      ! Some namelist arrays need initialising 
     239      cn_profbfiles(:)    = '' 
     240      cn_slafbfiles(:)    = '' 
     241      cn_sstfbfiles(:)    = '' 
     242      cn_sicfbfiles(:)    = '' 
     243      cn_velfbfiles(:)    = '' 
     244      cn_sssfbfiles(:)    = '' 
     245      cn_logchlfbfiles(:) = '' 
     246      cn_spmfbfiles(:)    = '' 
     247      cn_fco2fbfiles(:)   = '' 
     248      cn_pco2fbfiles(:)   = '' 
     249      cn_sstbiasfiles(:)  = '' 
     250      nn_profdavtypes(:)  = -1 
     251 
     252      CALL ini_date( rn_dobsini ) 
     253      CALL fin_date( rn_dobsend ) 
     254 
     255      ! Read namelist namobs : control observation diagnostics 
     256      REWIND( numnam_ref )   ! Namelist namobs in reference namelist 
    240257      READ  ( numnam_ref, namobs, IOSTAT = ios, ERR = 901) 
    241258901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in reference namelist', lwp ) 
    242259 
    243       REWIND( numnam_cfg )              ! Namelist namobs in configuration namelist : Diagnostic: control observation 
     260      REWIND( numnam_cfg )   ! Namelist namobs in configuration namelist 
    244261      READ  ( numnam_cfg, namobs, IOSTAT = ios, ERR = 902 ) 
    245262902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in configuration namelist', lwp ) 
    246263      IF(lwm) WRITE ( numond, namobs ) 
    247264 
    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) 
    253       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 
     265      lk_diaobs = .FALSE. 
     266#if defined key_diaobs 
     267      IF ( ln_diaobs ) lk_diaobs = .TRUE. 
     268#endif 
     269 
     270      IF ( .NOT. lk_diaobs ) THEN 
     271         IF(lwp) WRITE(numout,cform_war) 
     272         IF(lwp) WRITE(numout,*)' ln_diaobs is set to false or key_diaobs is not set, so not calling dia_obs' 
     273         RETURN 
     274      ENDIF 
     275 
    322276      IF(lwp) THEN 
    323277         WRITE(numout,*) 
    … …  
    325279         WRITE(numout,*) '~~~~~~~~~~~~' 
    326280         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 
     281         WRITE(numout,*) '             Logical switch for T profile observations                ln_t3d = ', ln_t3d 
     282         WRITE(numout,*) '             Logical switch for S profile observations                ln_s3d = ', ln_s3d 
     283         WRITE(numout,*) '             Logical switch for SLA observations                      ln_sla = ', ln_sla 
     284         WRITE(numout,*) '             Logical switch for SST observations                      ln_sst = ', ln_sst 
     285         WRITE(numout,*) '             Logical switch for Sea Ice observations                  ln_sic = ', ln_sic 
     286         WRITE(numout,*) '             Logical switch for velocity observations               ln_vel3d = ', ln_vel3d 
     287         WRITE(numout,*) '             Logical switch for SSS observations                      ln_sss = ', ln_sss 
     288         WRITE(numout,*) '             Logical switch for log(Chl) observations              ln_logchl = ', ln_logchl 
     289         WRITE(numout,*) '             Logical switch for SPM observations                      ln_spm = ', ln_spm 
     290         WRITE(numout,*) '             Logical switch for FCO2 observations                    ln_fco2 = ', ln_fco2 
     291         WRITE(numout,*) '             Logical switch for PCO2 observations                    ln_pco2 = ', ln_pco2 
     292         WRITE(numout,*) '             Global distribution of observations              ln_grid_global = ', ln_grid_global 
     293         WRITE(numout,*) '             Logical switch for obs grid search lookup ln_grid_search_lookup = ', ln_grid_search_lookup 
    352294         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 
     295            WRITE(numout,*) '             Grid search lookup file header                cn_gridsearchfile = ', cn_gridsearchfile 
     296         WRITE(numout,*) '             Initial date in window YYYYMMDD.HHMMSS               rn_dobsini = ', rn_dobsini 
     297         WRITE(numout,*) '             Final date in window YYYYMMDD.HHMMSS                 rn_dobsend = ', rn_dobsend 
     298         WRITE(numout,*) '             Type of vertical interpolation method                  nn_1dint = ', nn_1dint 
     299         WRITE(numout,*) '             Type of horizontal interpolation method                nn_2dint = ', nn_2dint 
     300         WRITE(numout,*) '             Rejection of observations near land switch               ln_nea = ', ln_nea 
     301         WRITE(numout,*) '             Rejection of obs near open bdys                 ln_bound_reject = ', ln_bound_reject 
     302         WRITE(numout,*) '             MSSH correction scheme                                 nn_msshc = ', nn_msshc 
     303         WRITE(numout,*) '             MDT  correction                                      rn_mdtcorr = ', rn_mdtcorr 
     304         WRITE(numout,*) '             MDT cutoff for computed correction                 rn_mdtcutoff = ', rn_mdtcutoff 
     305         WRITE(numout,*) '             Logical switch for alt bias                          ln_altbias = ', ln_altbias 
     306         WRITE(numout,*) '             Logical switch for sst bias                          ln_sstbias = ', ln_sstbias 
     307         WRITE(numout,*) '             Logical switch for ignoring missing files             ln_ignmis = ', ln_ignmis 
     308         WRITE(numout,*) '             Daily average types                             nn_profdavtypes = ', nn_profdavtypes 
     309         WRITE(numout,*) '             Logical switch for night-time SST obs               ln_sstnight = ', ln_sstnight 
     310      ENDIF 
     311      !----------------------------------------------------------------------- 
     312      ! Set up list of observation types to be used 
     313      ! and the files associated with each type 
     314      !----------------------------------------------------------------------- 
     315 
     316      nproftypes = COUNT( (/ln_t3d .OR. ln_s3d, ln_vel3d /) ) 
     317      nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic, ln_sss, & 
     318         &                  ln_logchl, ln_spm, ln_fco2, ln_pco2 /) ) 
     319 
     320      IF ( nproftypes == 0 .AND. nsurftypes == 0 ) THEN 
     321         IF(lwp) WRITE(numout,cform_war) 
     322         IF(lwp) WRITE(numout,*) ' ln_diaobs is set to true, but all obs operator logical flags', & 
     323            &                    ' are set to .FALSE. so turning off calls to dia_obs' 
     324         nwarn = nwarn + 1 
     325         lk_diaobs = .FALSE. 
     326         RETURN 
     327      ENDIF 
     328 
     329      IF(lwp) WRITE(numout,*) '          Number of profile obs types: ',nproftypes 
     330      IF ( nproftypes > 0 ) THEN 
     331 
     332         ALLOCATE( cobstypesprof(nproftypes) ) 
     333         ALLOCATE( ifilesprof(nproftypes) ) 
     334         ALLOCATE( clproffiles(nproftypes,jpmaxnfiles) ) 
     335 
     336         jtype = 0 
     337         IF (ln_t3d .OR. ln_s3d) THEN 
     338            jtype = jtype + 1 
     339            CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'prof  ', & 
     340               &                   cn_profbfiles, ifilesprof, cobstypesprof, clproffiles ) 
    359341         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 
     342         IF (ln_vel3d) THEN 
     343            jtype = jtype + 1 
     344            CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'vel   ', & 
     345               &                   cn_velfbfiles, ifilesprof, cobstypesprof, clproffiles ) 
    365346         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)) 
    371                ELSE 
    372                   WRITE(numout,'(1X,2A)') '             Feedback input observation file name       profbfiles = ', & 
    373                      TRIM(profbfiles(ji)) 
    374                ENDIF 
    375                WRITE(numout,'(1X,2A)') '       Enact feedback input time setting switch    ln_profb_enatim = ', ln_profb_enatim(ji) 
    376             END DO 
     347 
     348      ENDIF 
     349 
     350      IF(lwp) WRITE(numout,*)'          Number of surface obs types: ',nsurftypes 
     351      IF ( nsurftypes > 0 ) THEN 
     352 
     353         ALLOCATE( cobstypessurf(nsurftypes) ) 
     354         ALLOCATE( ifilessurf(nsurftypes) ) 
     355         ALLOCATE( clsurffiles(nsurftypes, jpmaxnfiles) ) 
     356         ALLOCATE(n2dintsurf(nsurftypes)) 
     357         ALLOCATE(ravglamscl(nsurftypes)) 
     358         ALLOCATE(ravgphiscl(nsurftypes)) 
     359         ALLOCATE(lfpindegs(nsurftypes)) 
     360         ALLOCATE(llnightav(nsurftypes)) 
     361 
     362         jtype = 0 
     363         IF (ln_sla) THEN 
     364            jtype = jtype + 1 
     365            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sla   ', & 
     366               &                   cn_slafbfiles, ifilessurf, cobstypessurf, clsurffiles ) 
     367            CALL obs_setinterpopts( nsurftypes, jtype, 'sla   ',      & 
     368               &                  nn_2dint, nn_2dint_sla,             & 
     369               &                  rn_sla_avglamscl, rn_sla_avgphiscl, & 
     370               &                  ln_sla_fp_indegs, .FALSE.,          & 
     371               &                  n2dintsurf, ravglamscl, ravgphiscl, & 
     372               &                  lfpindegs, llnightav ) 
    377373         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 
     374         IF (ln_sst) THEN 
     375            jtype = jtype + 1 
     376            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sst   ', & 
     377               &                   cn_sstfbfiles, ifilessurf, cobstypessurf, clsurffiles ) 
     378            CALL obs_setinterpopts( nsurftypes, jtype, 'sst   ',      & 
     379               &                  nn_2dint, nn_2dint_sst,             & 
     380               &                  rn_sst_avglamscl, rn_sst_avgphiscl, & 
     381               &                  ln_sst_fp_indegs, ln_sstnight,      & 
     382               &                  n2dintsurf, ravglamscl, ravgphiscl, & 
     383               &                  lfpindegs, llnightav ) 
    387384         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 
     385#if defined key_lim2 || defined key_lim3 || defined key_cice 
     386         IF (ln_sic) THEN 
     387            jtype = jtype + 1 
     388            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sic   ', & 
     389               &                   cn_sicfbfiles, ifilessurf, cobstypessurf, clsurffiles ) 
     390            CALL obs_setinterpopts( nsurftypes, jtype, 'sic   ',      & 
     391               &                  nn_2dint, nn_2dint_sic,             & 
     392               &                  rn_sic_avglamscl, rn_sic_avgphiscl, & 
     393               &                  ln_sic_fp_indegs, .FALSE.,          & 
     394               &                  n2dintsurf, ravglamscl, ravgphiscl, & 
     395               &                  lfpindegs, llnightav ) 
    393396         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 
     397#endif 
     398         IF (ln_sss) THEN 
     399            jtype = jtype + 1 
     400            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sss   ', & 
     401               &                   cn_sssfbfiles, ifilessurf, cobstypessurf, clsurffiles ) 
     402            CALL obs_setinterpopts( nsurftypes, jtype, 'sss   ',      & 
     403               &                  nn_2dint, nn_2dint_sss,             & 
     404               &                  rn_sss_avglamscl, rn_sss_avgphiscl, & 
     405               &                  ln_sss_fp_indegs, .FALSE.,          & 
     406               &                  n2dintsurf, ravglamscl, ravgphiscl, & 
     407               &                  lfpindegs, llnightav ) 
    399408         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 
     409 
     410         IF (ln_logchl) THEN 
     411            jtype = jtype + 1 
     412            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'logchl', & 
     413               &                   cn_logchlfbfiles, ifilessurf, cobstypessurf, clsurffiles ) 
     414            CALL obs_setinterpopts( nsurftypes, jtype, 'logchl',         & 
     415               &                    nn_2dint, -1, 0., 0., .TRUE., .FALSE., & 
     416               &                    n2dintsurf, ravglamscl, ravgphiscl,    & 
     417               &                    lfpindegs, llnightav ) 
    405418         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 
     419 
     420         IF (ln_spm) THEN 
     421            jtype = jtype + 1 
     422            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'spm   ', & 
     423               &                   cn_spmfbfiles, ifilessurf, cobstypessurf, clsurffiles ) 
     424            CALL obs_setinterpopts( nsurftypes, jtype, 'spm   ',         & 
     425               &                    nn_2dint, -1, 0., 0., .TRUE., .FALSE., & 
     426               &                    n2dintsurf, ravglamscl, ravgphiscl,    & 
     427               &                    lfpindegs, llnightav ) 
    411428         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 
     429 
     430         IF (ln_fco2) THEN 
     431            jtype = jtype + 1 
     432            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'fco2  ', & 
     433               &                   cn_fco2fbfiles, ifilessurf, cobstypessurf, clsurffiles ) 
     434            CALL obs_setinterpopts( nsurftypes, jtype, 'fco2  ',         & 
     435               &                    nn_2dint, -1, 0., 0., .TRUE., .FALSE., & 
     436               &                    n2dintsurf, ravglamscl, ravgphiscl,    & 
     437               &                    lfpindegs, llnightav ) 
    417438         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 
     439 
     440         IF (ln_pco2) THEN 
     441            jtype = jtype + 1 
     442            CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'pco2  ', & 
     443               &                   cn_pco2fbfiles, ifilessurf, cobstypessurf, clsurffiles ) 
     444            CALL obs_setinterpopts( nsurftypes, jtype, 'pco2  ',         & 
     445               &                    nn_2dint, -1, 0., 0., .TRUE., .FALSE., & 
     446               &                    n2dintsurf, ravglamscl, ravgphiscl,    & 
     447               &                    lfpindegs, llnightav ) 
    423448         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)) 
    441                ELSE 
    442                   WRITE(numout,'(1X,2A)') '             Vel. feedback input observation file name  velfbfiles = ', & 
    443                      TRIM(velfbfiles(ji)) 
    444                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 
    458  
    459       ENDIF 
    460        
     449 
     450      ENDIF 
     451 
     452      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 
     453 
     454 
     455      !----------------------------------------------------------------------- 
     456      ! Obs operator parameter checking and initialisations 
     457      !----------------------------------------------------------------------- 
     458 
    461459      IF ( ln_vel3d .AND. ( .NOT. ln_grid_global ) ) THEN 
    462460         CALL ctl_stop( 'Velocity data only works with ln_grid_global=.true.' ) 
    … …  
    464462      ENDIF 
    465463 
    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 
     464      IF ( ( nn_1dint < 0 ) .OR. ( nn_1dint > 1 ) ) THEN 
    485465         CALL ctl_stop(' Choice of vertical (1D) interpolation method', & 
    486466            &                    ' is not available') 
    487467      ENDIF 
    488       IF ( ( n2dint < 0 ).OR.( n2dint > 4 ) ) THEN 
     468 
     469      IF ( ( nn_2dint < 0 ) .OR. ( nn_2dint > 6 ) ) THEN 
    489470         CALL ctl_stop(' Choice of horizontal (2D) interpolation method', & 
    490471            &                    ' is not available') 
    491472      ENDIF 
    492473 
     474      CALL obs_typ_init 
     475 
     476      CALL mppmap_init 
     477 
     478      CALL obs_grid_setup( ) 
     479 
    493480      !----------------------------------------------------------------------- 
    494481      ! Depending on switches read the various observation types 
    495482      !----------------------------------------------------------------------- 
    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 
    524              
    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  
     483 
     484      IF ( nproftypes > 0 ) THEN 
     485 
     486         ALLOCATE(profdata(nproftypes)) 
     487         ALLOCATE(profdataqc(nproftypes)) 
     488         ALLOCATE(nvarsprof(nproftypes)) 
     489         ALLOCATE(nextrprof(nproftypes)) 
     490 
     491         DO jtype = 1, nproftypes 
     492 
     493            nvarsprof(jtype) = 2 
     494            IF ( TRIM(cobstypesprof(jtype)) == 'prof' ) THEN 
     495               nextrprof(jtype) = 1 
     496               llvar1 = ln_t3d 
     497               llvar2 = ln_s3d 
     498               zglam1 = glamt 
     499               zgphi1 = gphit 
     500               zmask1 = tmask 
     501               zglam2 = glamt 
     502               zgphi2 = gphit 
     503               zmask2 = tmask 
     504            ENDIF 
     505            IF ( TRIM(cobstypesprof(jtype)) == 'vel' )  THEN 
     506               nextrprof(jtype) = 2 
     507               llvar1 = ln_vel3d 
     508               llvar2 = ln_vel3d 
     509               zglam1 = glamu 
     510               zgphi1 = gphiu 
     511               zmask1 = umask 
     512               zglam2 = glamv 
     513               zgphi2 = gphiv 
     514               zmask2 = vmask 
     515            ENDIF 
     516 
     517            !Read in profile or profile obs types 
     518            CALL obs_rea_prof( profdata(jtype), ifilesprof(jtype),       & 
     519               &               clproffiles(jtype,1:ifilesprof(jtype)), & 
     520               &               nvarsprof(jtype), nextrprof(jtype), nitend-nit000+2, & 
     521               &               rn_dobsini, rn_dobsend, llvar1, llvar2, & 
     522               &               ln_ignmis, ln_s_at_t, .FALSE., & 
     523               &               kdailyavtypes = nn_profdavtypes ) 
     524 
     525            DO jvar = 1, nvarsprof(jtype) 
     526               CALL obs_prof_staend( profdata(jtype), jvar ) 
    539527            END DO 
    540528 
    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                    
     529            CALL obs_pre_prof( profdata(jtype), profdataqc(jtype), & 
     530               &               llvar1, llvar2, & 
     531               &               jpi, jpj, jpk, & 
     532               &               zmask1, zglam1, zgphi1, zmask2, zglam2, zgphi2,  & 
     533               &               ln_nea, ln_bound_reject, & 
     534               &               kdailyavtypes = nn_profdavtypes ) 
     535 
     536         END DO 
     537 
     538         DEALLOCATE( ifilesprof, clproffiles ) 
     539 
     540      ENDIF 
     541 
     542      IF ( nsurftypes > 0 ) THEN 
     543 
     544         ALLOCATE(surfdata(nsurftypes)) 
     545         ALLOCATE(surfdataqc(nsurftypes)) 
     546         ALLOCATE(nvarssurf(nsurftypes)) 
     547         ALLOCATE(nextrsurf(nsurftypes)) 
     548 
     549         DO jtype = 1, nsurftypes 
     550 
     551            nvarssurf(jtype) = 1 
     552            nextrsurf(jtype) = 0 
     553            IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) nextrsurf(jtype) = 2 
     554 
     555            !Read in surface obs types 
     556            CALL obs_rea_surf( surfdata(jtype), ifilessurf(jtype), & 
     557               &               clsurffiles(jtype,1:ifilessurf(jtype)), & 
     558               &               nvarssurf(jtype), nextrsurf(jtype), nitend-nit000+2, & 
     559               &               rn_dobsini, rn_dobsend, ln_ignmis, .FALSE., llnightav(jtype) ) 
     560 
     561            CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea, ln_bound_reject ) 
     562 
     563            IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN 
     564               CALL obs_rea_mdt( surfdataqc(jtype), n2dintsurf(jtype) ) 
     565               IF ( ln_altbias ) & 
     566                  & CALL obs_rea_altbias ( surfdataqc(jtype), n2dintsurf(jtype), cn_altbiasfile ) 
     567            ENDIF 
     568 
     569            IF ( TRIM(cobstypessurf(jtype)) == 'sst' .AND. ln_sstbias ) THEN 
     570               jnumsstbias = 0 
     571               DO jfile = 1, jpmaxnfiles 
     572                  IF ( TRIM(cn_sstbiasfiles(jfile)) /= '' ) & 
     573                     &  jnumsstbias = jnumsstbias + 1 
    596574               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 ) 
     575               IF ( jnumsstbias == 0 ) THEN 
     576                  CALL ctl_stop("ln_sstbias set but no bias files to read in")     
    605577               ENDIF 
    606                 
    607             END DO 
    608  
    609          ENDIF 
    610  
    611       ENDIF 
    612  
    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       
     578 
     579               CALL obs_app_sstbias( surfdataqc(jtype), n2dintsurf(jtype), &  
     580                  &                  jnumsstbias, cn_sstbiasfiles(1:jnumsstbias) )  
     581 
     582            ENDIF 
     583 
     584         END DO 
     585 
     586         DEALLOCATE( ifilessurf, clsurffiles ) 
     587 
     588      ENDIF 
     589 
     590      CALL wrk_dealloc( jpi, jpj, zglam1 ) 
     591      CALL wrk_dealloc( jpi, jpj, zglam2 ) 
     592      CALL wrk_dealloc( jpi, jpj, zgphi1 ) 
     593      CALL wrk_dealloc( jpi, jpj, zgphi2 ) 
     594      CALL wrk_dealloc( jpi, jpj, jpk, zmask1 ) 
     595      CALL wrk_dealloc( jpi, jpj, jpk, zmask2 ) 
     596 
    967597   END SUBROUTINE dia_obs_init 
    968598 
    … …  
    974604      !! 
    975605      !! ** 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  :  
     606      !!              compute the model equivalent of the following data: 
     607      !!               - Profile data, currently T/S or U/V 
     608      !!               - Surface data, currently SST, SLA or sea-ice concentration. 
     609      !! 
     610      !! ** Action  : 
    985611      !! 
    986612      !! History : 
    … …  
    991617      !!        !  07-04  (G. Smith) Generalized surface operators 
    992618      !!        !  08-10  (M. Valdivieso) obs operator for velocity profiles 
     619      !!        !  15-08  (M. Martin) Combined surface/profile routines. 
    993620      !!---------------------------------------------------------------------- 
    994621      !! * 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,  & 
     622      USE phycst, ONLY : &         ! Physical constants 
     623         & rday 
     624      USE oce, ONLY : &            ! Ocean dynamics and tracers variables 
     625         & tsn,       & 
     626         & un,        & 
     627         & vn,        & 
    1004628         & sshn 
    1005629#if defined  key_lim3 
    1006       USE ice, ONLY : &                     ! LIM Ice model variables 
     630      USE ice, ONLY : &            ! LIM3 Ice model variables 
    1007631         & frld 
    1008632#endif 
    1009633#if defined key_lim2 
    1010       USE ice_2, ONLY : &                     ! LIM Ice model variables 
     634      USE ice_2, ONLY : &          ! LIM2 Ice model variables 
    1011635         & frld 
    1012636#endif 
     637#if defined key_cice 
     638      USE sbc_oce, ONLY : fr_i     ! ice fraction 
     639#endif 
     640#if defined key_hadocc 
     641      USE trc, ONLY :  &           ! HadOCC chlorophyll, fCO2 and pCO2 
     642         & HADOCC_CHL, & 
     643         & HADOCC_FCO2, & 
     644         & HADOCC_PCO2, & 
     645         & HADOCC_FILL_FLT 
     646#elif defined key_medusa && defined key_foam_medusa 
     647      USE trc, ONLY :  &           ! MEDUSA chlorophyll, fCO2 and pCO2 
     648         & MEDUSA_CHL, & 
     649         & MEDUSA_FCO2, & 
     650         & MEDUSA_PCO2, & 
     651         & MEDUSA_FILL_FLT 
     652#elif defined key_fabm 
     653      USE fabm 
     654      USE par_fabm 
     655#endif 
     656#if defined key_spm 
     657      USE par_spm, ONLY: &         ! ERSEM/SPM sediments 
     658         & jp_spm 
     659      USE trc, ONLY :  & 
     660         & trn 
     661#endif 
     662 
    1013663      IMPLICIT NONE 
    1014664 
    1015665      !! * Arguments 
    1016       INTEGER, INTENT(IN) :: kstp                         ! Current timestep 
     666      INTEGER, INTENT(IN) :: kstp  ! Current timestep 
    1017667      !! * 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)  
    1032 #endif 
     668      INTEGER :: idaystp           ! Number of timesteps per day 
     669      INTEGER :: jtype             ! Data loop variable 
     670      INTEGER :: jvar              ! Variable number 
     671      INTEGER :: ji, jj            ! Loop counters 
     672      REAL(wp) :: tiny                  ! small number 
     673      REAL(wp), POINTER, DIMENSION(:,:,:) :: & 
     674         & zprofvar1, &            ! Model values for 1st variable in a prof ob 
     675         & zprofvar2               ! Model values for 2nd variable in a prof ob 
     676      REAL(wp), POINTER, DIMENSION(:,:,:) :: & 
     677         & zprofmask1, &           ! Mask associated with zprofvar1 
     678         & zprofmask2              ! Mask associated with zprofvar2 
     679      REAL(wp), POINTER, DIMENSION(:,:) :: & 
     680         & zsurfvar, &             ! Model values equivalent to surface ob. 
     681         & zsurfmask               ! Mask associated with surface variable 
     682      REAL(wp), POINTER, DIMENSION(:,:) :: & 
     683         & zglam1,    &            ! Model longitudes for prof variable 1 
     684         & zglam2,    &            ! Model longitudes for prof variable 2 
     685         & zgphi1,    &            ! Model latitudes for prof variable 1 
     686         & zgphi2                  ! Model latitudes for prof variable 2 
     687 
     688 
     689      !Allocate local work arrays 
     690      CALL wrk_alloc( jpi, jpj, jpk, zprofvar1 ) 
     691      CALL wrk_alloc( jpi, jpj, jpk, zprofvar2 ) 
     692      CALL wrk_alloc( jpi, jpj, jpk, zprofmask1 ) 
     693      CALL wrk_alloc( jpi, jpj, jpk, zprofmask2 ) 
     694      CALL wrk_alloc( jpi, jpj, zsurfvar ) 
     695      CALL wrk_alloc( jpi, jpj, zsurfmask ) 
     696      CALL wrk_alloc( jpi, jpj, zglam1 ) 
     697      CALL wrk_alloc( jpi, jpj, zglam2 ) 
     698      CALL wrk_alloc( jpi, jpj, zgphi1 ) 
     699      CALL wrk_alloc( jpi, jpj, zgphi2 ) 
    1033700 
    1034701      IF(lwp) THEN 
    … …  
    1036703         WRITE(numout,*) 'dia_obs : Call the observation operators', kstp 
    1037704         WRITE(numout,*) '~~~~~~~' 
     705         CALL FLUSH(numout) 
    1038706      ENDIF 
    1039707 
    … …  
    1041709 
    1042710      !----------------------------------------------------------------------- 
    1043       ! No LIM => frld == 0.0_wp 
    1044       !----------------------------------------------------------------------- 
    1045 #if ! defined key_lim2 && ! defined key_lim3 
    1046       frld(:,:) = 0.0_wp 
     711      ! Call the profile and surface observation operators 
     712      !----------------------------------------------------------------------- 
     713 
     714      IF ( nproftypes > 0 ) THEN 
     715 
     716         DO jtype = 1, nproftypes 
     717 
     718            SELECT CASE ( TRIM(cobstypesprof(jtype)) ) 
     719            CASE('prof') 
     720               zprofvar1(:,:,:) = tsn(:,:,:,jp_tem) 
     721               zprofvar2(:,:,:) = tsn(:,:,:,jp_sal) 
     722               zprofmask1(:,:,:) = tmask(:,:,:) 
     723               zprofmask2(:,:,:) = tmask(:,:,:) 
     724               zglam1(:,:) = glamt(:,:) 
     725               zglam2(:,:) = glamt(:,:) 
     726               zgphi1(:,:) = gphit(:,:) 
     727               zgphi2(:,:) = gphit(:,:) 
     728            CASE('vel') 
     729               zprofvar1(:,:,:) = un(:,:,:) 
     730               zprofvar2(:,:,:) = vn(:,:,:) 
     731               zprofmask1(:,:,:) = umask(:,:,:) 
     732               zprofmask2(:,:,:) = vmask(:,:,:) 
     733               zglam1(:,:) = glamu(:,:) 
     734               zglam2(:,:) = glamv(:,:) 
     735               zgphi1(:,:) = gphiu(:,:) 
     736               zgphi2(:,:) = gphiv(:,:) 
     737            CASE DEFAULT 
     738               CALL ctl_stop( 'Unknown profile observation type '//TRIM(cobstypesprof(jtype))//' in dia_obs' ) 
     739            END SELECT 
     740 
     741            CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk,  & 
     742               &               nit000, idaystp,                         & 
     743               &               zprofvar1, zprofvar2,                    & 
     744               &               fsdept(:,:,:), fsdepw(:,:,:),            &  
     745               &               zprofmask1, zprofmask2,                  & 
     746               &               zglam1, zglam2, zgphi1, zgphi2,          & 
     747               &               nn_1dint, nn_2dint,                      & 
     748               &               kdailyavtypes = nn_profdavtypes ) 
     749 
     750         END DO 
     751 
     752      ENDIF 
     753 
     754      IF ( nsurftypes > 0 ) THEN 
     755 
     756         DO jtype = 1, nsurftypes 
     757 
     758            !Defaults which might be changed 
     759            zsurfmask(:,:) = tmask(:,:,1) 
     760 
     761            SELECT CASE ( TRIM(cobstypessurf(jtype)) ) 
     762            CASE('sst') 
     763               zsurfvar(:,:) = tsn(:,:,1,jp_tem) 
     764            CASE('sla') 
     765               zsurfvar(:,:) = sshn(:,:) 
     766            CASE('sss') 
     767               zsurfvar(:,:) = tsn(:,:,1,jp_sal) 
     768            CASE('sic') 
     769               IF ( kstp == 0 ) THEN 
     770                  IF ( lwp .AND. surfdataqc(jtype)%nsstpmpp(1) > 0 ) THEN 
     771                     CALL ctl_warn( 'Sea-ice not initialised on zeroth '// & 
     772                        &           'time-step but some obs are valid then.' ) 
     773                     WRITE(numout,*)surfdataqc(jtype)%nsstpmpp(1), & 
     774                        &           ' sea-ice obs will be missed' 
     775                  ENDIF 
     776                  surfdataqc(jtype)%nsurfup = surfdataqc(jtype)%nsurfup + & 
     777                     &                        surfdataqc(jtype)%nsstp(1) 
     778                  CYCLE 
     779               ELSE 
     780#if defined key_cice 
     781                  zsurfvar(:,:) = fr_i(:,:) 
     782#elif defined key_lim2 || defined key_lim3 
     783                  zsurfvar(:,:) = 1._wp - frld(:,:) 
     784#else 
     785               CALL ctl_stop( ' Trying to run sea-ice observation operator', & 
     786                  &           ' but no sea-ice model appears to have been defined' ) 
    1047787#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              ) 
    1066             ENDIF 
     788               ENDIF 
     789 
     790            CASE('logchl') 
     791#if defined key_hadocc 
     792               zsurfvar(:,:) = HADOCC_CHL(:,:,1)    ! (not log) chlorophyll from HadOCC 
     793#elif defined key_medusa && defined key_foam_medusa 
     794               zsurfvar(:,:) = MEDUSA_CHL(:,:,1)    ! (not log) chlorophyll from HadOCC 
     795#elif defined key_fabm 
     796               chl_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabmdia_chltot) 
     797               zsurfvar(:,:) = chl_3d(:,:,1) 
     798#else 
     799               CALL ctl_stop( ' Trying to run logchl observation operator', & 
     800                  &           ' but no biogeochemical model appears to have been defined' ) 
     801#endif 
     802               zsurfmask(:,:) = tmask(:,:,1)         ! create a special mask to exclude certain things 
     803               ! Take the log10 where we can, otherwise exclude 
     804               tiny = 1.0e-20 
     805               WHERE(zsurfvar(:,:) > tiny .AND. zsurfvar(:,:) /= obfillflt ) 
     806                  zsurfvar(:,:)  = LOG10(zsurfvar(:,:)) 
     807               ELSEWHERE 
     808                  zsurfvar(:,:)  = obfillflt 
     809                  zsurfmask(:,:) = 0 
     810               END WHERE 
     811            CASE('spm') 
     812#if defined key_spm 
     813               zsurfvar(:,:) = 0.0 
     814               DO jn = 1, jp_spm 
     815                  zsurfvar(:,:) = zsurfvar(:,:) + trn(:,:,1,jn)   ! sum SPM sizes 
     816               END DO 
     817#else 
     818               CALL ctl_stop( ' Trying to run spm observation operator', & 
     819                  &           ' but no spm model appears to have been defined' ) 
     820#endif 
     821            CASE('fco2') 
     822#if defined key_hadocc 
     823               zsurfvar(:,:) = HADOCC_FCO2(:,:)    ! fCO2 from HadOCC 
     824               IF ( ( MINVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) .AND. & 
     825                  & ( MAXVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) ) THEN 
     826                  zsurfvar(:,:) = obfillflt 
     827                  zsurfmask(:,:) = 0 
     828                  CALL ctl_warn( ' HadOCC fCO2 values masked out for observation operator', & 
     829                     &           ' on timestep ' // TRIM(STR(kstp)),                              & 
     830                     &           ' as HADOCC_FCO2(:,:) == HADOCC_FILL_FLT' ) 
     831               ENDIF 
     832#elif defined key_medusa && defined key_foam_medusa 
     833               zsurfmask(:,:) = MEDUSA_FCO2(:,:)    ! fCO2 from MEDUSA 
     834               IF ( ( MINVAL( MEDUSA_FCO2 ) == MEDUSA_FILL_FLT ) .AND. & 
     835                  & ( MAXVAL( MEDUSA_FCO2 ) == MEDUSA_FILL_FLT ) ) THEN 
     836                  zsurfvar(:,:) = obfillflt 
     837                  zsurfmask(:,:) = 0 
     838                  CALL ctl_warn( ' MEDUSA fCO2 values masked out for observation operator', & 
     839                     &           ' on timestep ' // TRIM(STR(kstp)),                              & 
     840                     &           ' as MEDUSA_FCO2(:,:) == MEDUSA_FILL_FLT' ) 
     841               ENDIF 
     842#elif defined key_fabm 
     843               ! First, get pCO2 from FABM 
     844               pco2_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabm_o3pc) 
     845               zsurfvar(:,:) = pco2_3d(:,:,1) 
     846               ! Now, convert pCO2 to fCO2, based on SST in K. This follows the standard methodology of: 
     847               ! Pierrot et al. (2009), Recommendations for autonomous underway pCO2 measuring systems 
     848               ! and data reduction routines, Deep-Sea Research II, 56: 512-522. 
     849               ! and 
     850               ! Weiss (1974), Carbon dioxide in water and seawater: the solubility of a non-ideal gas, 
     851               ! Marine Chemistry, 2: 203-215. 
     852               ! In the implementation below, atmospheric pressure has been assumed to be 1 atm and so 
     853               ! not explicitly included - atmospheric pressure is not necessarily available so this is 
     854               ! the best assumption. 
     855               ! Further, the (1-xCO2)^2 term has been neglected. This is common practice 
     856               ! (see e.g. Zeebe and Wolf-Gladrow (2001), CO2 in Seawater: Equilibrium, Kinetics, Isotopes) 
     857               ! because xCO2 in atm is ~0, and so this term will only affect the result to the 3rd decimal 
     858               ! place for typical values, and xCO2 would need to be approximated from pCO2 anyway. 
     859               zsurfvar(:,:) = zsurfvar(:,:) * EXP((-1636.75                                                          + & 
     860                  &            12.0408      * (tsn(:,:,1,jp_tem)+rt0)                                                 - & 
     861                  &            0.0327957    * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0)                         + & 
     862                  &            0.0000316528 * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0) + & 
     863                  &            2.0 * (57.7 - 0.118 * (tsn(:,:,1,jp_tem)+rt0)))                                        / & 
     864                  &            (82.0578 * (tsn(:,:,1,jp_tem)+rt0))) 
     865#else 
     866               CALL ctl_stop( ' Trying to run fco2 observation operator', & 
     867                  &           ' but no biogeochemical model appears to have been defined' ) 
     868#endif 
     869            CASE('pco2') 
     870#if defined key_hadocc 
     871               zsurfvar(:,:) = HADOCC_PCO2(:,:)    ! pCO2 from HadOCC 
     872               IF ( ( MINVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) .AND. & 
     873                  & ( MAXVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) ) THEN 
     874                  zsurfvar(:,:) = obfillflt 
     875                  zsurfmask(:,:) = 0 
     876                  CALL ctl_warn( ' HadOCC pCO2 values masked out for observation operator', & 
     877                     &           ' on timestep ' // TRIM(STR(kstp)),                              & 
     878                     &           ' as HADOCC_PCO2(:,:) == HADOCC_FILL_FLT' ) 
     879               ENDIF 
     880#elif defined key_medusa && defined key_foam_medusa 
     881               zsurfvar(:,:) = MEDUSA_PCO2(:,:)    ! pCO2 from MEDUSA 
     882               IF ( ( MINVAL( MEDUSA_PCO2 ) == MEDUSA_FILL_FLT ) .AND. & 
     883                  & ( MAXVAL( MEDUSA_PCO2 ) == MEDUSA_FILL_FLT ) ) THEN 
     884                  zsurfvar(:,:) = obfillflt 
     885                  zsurfmask(:,:) = 0 
     886                  CALL ctl_warn( ' MEDUSA pCO2 values masked out for observation operator', & 
     887                     &           ' on timestep ' // TRIM(STR(kstp)),                              & 
     888                     &           ' as MEDUSA_PCO2(:,:) == MEDUSA_FILL_FLT' ) 
     889               ENDIF 
     890#elif defined key_fabm 
     891               pco2_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabm_o3pc) 
     892               zsurfvar(:,:) = pco2_3d(:,:,1) 
     893#else 
     894               CALL ctl_stop( ' Trying to run pCO2 observation operator', & 
     895                  &           ' but no biogeochemical model appears to have been defined' ) 
     896#endif 
     897 
     898            CASE DEFAULT 
     899 
     900               CALL ctl_stop( 'Unknown surface observation type '//TRIM(cobstypessurf(jtype))//' in dia_obs' ) 
     901 
     902            END SELECT 
     903 
     904            CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj,       & 
     905               &               nit000, idaystp, zsurfvar, zsurfmask,    & 
     906               &               n2dintsurf(jtype), llnightav(jtype),     & 
     907               &               ravglamscl(jtype), ravgphiscl(jtype),     & 
     908               &               lfpindegs(jtype) ) 
     909 
    1067910         END DO 
    1068       ENDIF 
    1069  
    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 
     911 
     912      ENDIF 
     913 
     914      CALL wrk_dealloc( jpi, jpj, jpk, zprofvar1 ) 
     915      CALL wrk_dealloc( jpi, jpj, jpk, zprofvar2 ) 
     916      CALL wrk_dealloc( jpi, jpj, jpk, zprofmask1 ) 
     917      CALL wrk_dealloc( jpi, jpj, jpk, zprofmask2 ) 
     918      CALL wrk_dealloc( jpi, jpj, zsurfvar ) 
     919      CALL wrk_dealloc( jpi, jpj, zsurfmask ) 
     920      CALL wrk_dealloc( jpi, jpj, zglam1 ) 
     921      CALL wrk_dealloc( jpi, jpj, zglam2 ) 
     922      CALL wrk_dealloc( jpi, jpj, zgphi1 ) 
     923      CALL wrk_dealloc( jpi, jpj, zgphi2 ) 
    1117924 
    1118925   END SUBROUTINE dia_obs 
    1119    
    1120    SUBROUTINE dia_obs_wri  
     926 
     927   SUBROUTINE dia_obs_wri 
    1121928      !!---------------------------------------------------------------------- 
    1122929      !!                    ***  ROUTINE dia_obs_wri  *** 
    … …  
    1126933      !! ** Method  : Call observation diagnostic output routines 
    1127934      !! 
    1128       !! ** Action  :  
     935      !! ** Action  : 
    1129936      !! 
    1130937      !! History : 
    … …  
    1134941      !!        !  07-03  (K. Mogensen) General handling of profiles 
    1135942      !!        !  08-09  (M. Valdivieso) Velocity component (U,V) profiles 
     943      !!        !  15-08  (M. Martin) Combined writing for prof and surf types 
    1136944      !!---------------------------------------------------------------------- 
     945      !! * Modules used 
     946      USE obs_rot_vel          ! Rotation of velocities 
     947 
    1137948      IMPLICIT NONE 
    1138949 
    1139950      !! * 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 
     951      INTEGER :: jtype                    ! Data set loop variable 
     952      INTEGER :: jo, jvar, jk 
     953      REAL(wp), DIMENSION(:), ALLOCATABLE :: & 
     954         & zu, & 
     955         & zv 
     956 
    1150957      !----------------------------------------------------------------------- 
    1151958      ! Depending on switches call various observation output routines 
    1152959      !----------------------------------------------------------------------- 
    1153960 
    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 ) 
     961      IF ( nproftypes > 0 ) THEN 
     962 
     963         DO jtype = 1, nproftypes 
     964 
     965            IF ( TRIM(cobstypesprof(jtype)) == 'vel' ) THEN 
     966 
     967               ! For velocity data, rotate the model velocities to N/S, E/W 
     968               ! using the compressed data structure. 
     969               ALLOCATE( & 
     970                  & zu(profdataqc(jtype)%nvprot(1)), & 
     971                  & zv(profdataqc(jtype)%nvprot(2))  & 
     972                  & ) 
     973 
     974               CALL obs_rotvel( profdataqc(jtype), nn_2dint, zu, zv ) 
     975 
     976               DO jo = 1, profdataqc(jtype)%nprof 
     977                  DO jvar = 1, 2 
     978                     DO jk = profdataqc(jtype)%npvsta(jo,jvar), profdataqc(jtype)%npvend(jo,jvar) 
     979 
     980                        IF ( jvar == 1 ) THEN 
     981                           profdataqc(jtype)%var(jvar)%vmod(jk) = zu(jk) 
     982                        ELSE 
     983                           profdataqc(jtype)%var(jvar)%vmod(jk) = zv(jk) 
     984                        ENDIF 
     985 
     986                     END DO 
     987                  END DO 
     988               END DO 
     989 
     990               DEALLOCATE( zu ) 
     991               DEALLOCATE( zv ) 
     992 
     993            END IF 
     994 
     995            CALL obs_prof_decompress( profdataqc(jtype), & 
     996               &                      profdata(jtype), .TRUE., numout ) 
     997 
     998            CALL obs_wri_prof( profdata(jtype) ) 
    1163999 
    11641000         END DO 
    11651001 
    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  
    1206       ENDIF 
    1207  
    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 ) 
     1002      ENDIF 
     1003 
     1004      IF ( nsurftypes > 0 ) THEN 
     1005 
     1006         DO jtype = 1, nsurftypes 
     1007 
     1008            CALL obs_surf_decompress( surfdataqc(jtype), & 
     1009               &                      surfdata(jtype), .TRUE., numout ) 
     1010 
     1011            CALL obs_wri_surf( surfdata(jtype) ) 
    12161012 
    12171013         END DO 
    12181014 
    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           
    13901015      ENDIF 
    13911016 
    … …  
    14051030      !! 
    14061031      !!---------------------------------------------------------------------- 
    1407       !! obs_grid deallocation 
     1032      ! obs_grid deallocation 
    14081033      CALL obs_grid_deallocate 
    14091034 
    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 
     1035      ! diaobs deallocation 
     1036      IF ( nproftypes > 0 ) & 
     1037         &   DEALLOCATE( cobstypesprof, profdata, profdataqc, nvarsprof, nextrprof ) 
     1038 
     1039      IF ( nsurftypes > 0 ) & 
     1040         &   DEALLOCATE( cobstypessurf, surfdata, surfdataqc, nvarssurf, nextrsurf, & 
     1041         &               n2dintsurf, ravglamscl, ravgphiscl, lfpindegs, llnightav ) 
     1042 
    14331043   END SUBROUTINE dia_obs_dealloc 
    14341044 
    … …  
    14361046      !!---------------------------------------------------------------------- 
    14371047      !!                    ***  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 
     1048      !! 
     1049      !! ** Purpose : Get initial date in double precision YYYYMMDD.HHMMSS format 
     1050      !! 
     1051      !! ** Method  : Get initial date in double precision YYYYMMDD.HHMMSS format 
     1052      !! 
     1053      !! ** Action  : Get initial date in double precision YYYYMMDD.HHMMSS format 
    14441054      !! 
    14451055      !! History : 
    … …  
    14521062      USE phycst, ONLY : &            ! Physical constants 
    14531063         & rday 
    1454 !      USE daymod, ONLY : &            ! Time variables 
    1455 !         & nmonth_len            
    14561064      USE dom_oce, ONLY : &           ! Ocean space and time domain variables 
    14571065         & rdt 
    … …  
    14601068 
    14611069      !! * Arguments 
    1462       REAL(KIND=dp), INTENT(OUT) :: ddobsini                         ! Initial date in YYYYMMDD.HHMMSS 
     1070      REAL(dp), INTENT(OUT) :: ddobsini  ! Initial date in YYYYMMDD.HHMMSS 
    14631071 
    14641072      !! * Local declarations 
    … …  
    14681076      INTEGER :: ihou 
    14691077      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       !!---------------------------------------------------------------------- 
     1078      INTEGER :: imday       ! Number of days in month. 
     1079      INTEGER, DIMENSION(12) :: & 
     1080         &       imonth_len  ! Length in days of the months of the current year 
     1081      REAL(wp) :: zdayfrc    ! Fraction of day 
     1082 
     1083      !---------------------------------------------------------------------- 
     1084      ! Initial date initialization (year, month, day, hour, minute) 
     1085      ! (This assumes that the initial date is for 00z)) 
     1086      !---------------------------------------------------------------------- 
    14791087      iyea =   ndate0 / 10000 
    14801088      imon = ( ndate0 - iyea * 10000 ) / 100 
    … …  
    14831091      imin = 0 
    14841092 
    1485       !!---------------------------------------------------------------------- 
    1486       !! Compute number of days + number of hours + min since initial time 
    1487       !!---------------------------------------------------------------------- 
     1093      !---------------------------------------------------------------------- 
     1094      ! Compute number of days + number of hours + min since initial time 
     1095      !---------------------------------------------------------------------- 
    14881096      iday = iday + ( nit000 -1 ) * rdt / rday 
    14891097      zdayfrc = ( nit000 -1 ) * rdt / rday 
    … …  
    14921100      imin = int( (zdayfrc * 24 - ihou) * 60 ) 
    14931101 
    1494       !!----------------------------------------------------------------------- 
    1495       !! Convert number of days (iday) into a real date 
    1496       !!---------------------------------------------------------------------- 
     1102      !----------------------------------------------------------------------- 
     1103      ! Convert number of days (iday) into a real date 
     1104      !---------------------------------------------------------------------- 
    14971105 
    14981106      CALL calc_month_len( iyea, imonth_len ) 
    1499        
     1107 
    15001108      DO WHILE ( iday > imonth_len(imon) ) 
    15011109         iday = iday - imonth_len(imon) 
    … …  
    15081116      END DO 
    15091117 
    1510       !!---------------------------------------------------------------------- 
    1511       !! Convert it into YYYYMMDD.HHMMSS format. 
    1512       !!---------------------------------------------------------------------- 
     1118      !---------------------------------------------------------------------- 
     1119      ! Convert it into YYYYMMDD.HHMMSS format. 
     1120      !---------------------------------------------------------------------- 
    15131121      ddobsini = iyea * 10000_dp + imon * 100_dp + & 
    15141122         &       iday + ihou * 0.01_dp + imin * 0.0001_dp 
    … …  
    15201128      !!---------------------------------------------------------------------- 
    15211129      !!                    ***  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 
     1130      !! 
     1131      !! ** Purpose : Get final date in double precision YYYYMMDD.HHMMSS format 
     1132      !! 
     1133      !! ** Method  : Get final date in double precision YYYYMMDD.HHMMSS format 
     1134      !! 
     1135      !! ** Action  : Get final date in double precision YYYYMMDD.HHMMSS format 
    15281136      !! 
    15291137      !! History : 
    … …  
    15351143      USE phycst, ONLY : &            ! Physical constants 
    15361144         & rday 
    1537 !      USE daymod, ONLY : &            ! Time variables 
    1538 !         & nmonth_len                 
    15391145      USE dom_oce, ONLY : &           ! Ocean space and time domain variables 
    15401146         & rdt 
    … …  
    15431149 
    15441150      !! * Arguments 
    1545       REAL(KIND=dp), INTENT(OUT) :: ddobsfin                  ! Final date in YYYYMMDD.HHMMSS 
     1151      REAL(dp), INTENT(OUT) :: ddobsfin ! Final date in YYYYMMDD.HHMMSS 
    15461152 
    15471153      !! * Local declarations 
    … …  
    15511157      INTEGER :: ihou 
    15521158      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              
     1159      INTEGER :: imday       ! Number of days in month. 
     1160      INTEGER, DIMENSION(12) :: & 
     1161         &       imonth_len  ! Length in days of the months of the current year 
     1162      REAL(wp) :: zdayfrc    ! Fraction of day 
     1163 
    15581164      !----------------------------------------------------------------------- 
    15591165      ! Initial date initialization (year, month, day, hour, minute) 
    … …  
    15651171      ihou = 0 
    15661172      imin = 0 
    1567        
     1173 
    15681174      !----------------------------------------------------------------------- 
    15691175      ! Compute number of days + number of hours + min since initial time 
    … …  
    15801186 
    15811187      CALL calc_month_len( iyea, imonth_len ) 
    1582        
     1188 
    15831189      DO WHILE ( iday > imonth_len(imon) ) 
    15841190         iday = iday - imonth_len(imon) 
    … …  
    15981204 
    15991205    END SUBROUTINE fin_date 
    1600      
     1206 
     1207    SUBROUTINE obs_settypefiles( ntypes, jpmaxnfiles, jtype, ctypein, & 
     1208       &                         cfilestype, ifiles, cobstypes, cfiles ) 
     1209 
     1210    INTEGER, INTENT(IN) :: ntypes      ! Total number of obs types 
     1211    INTEGER, INTENT(IN) :: jpmaxnfiles ! Maximum number of files allowed for each type 
     1212    INTEGER, INTENT(IN) :: jtype       ! Index of the current type of obs 
     1213    INTEGER, DIMENSION(ntypes), INTENT(INOUT) :: & 
     1214       &                   ifiles      ! Out appended number of files for this type 
     1215 
     1216    CHARACTER(len=6), INTENT(IN) :: ctypein  
     1217    CHARACTER(len=128), DIMENSION(jpmaxnfiles), INTENT(IN) :: & 
     1218       &                   cfilestype  ! In list of files for this obs type 
     1219    CHARACTER(len=6), DIMENSION(ntypes), INTENT(INOUT) :: & 
     1220       &                   cobstypes   ! Out appended list of obs types 
     1221    CHARACTER(len=128), DIMENSION(ntypes, jpmaxnfiles), INTENT(INOUT) :: & 
     1222       &                   cfiles      ! Out appended list of files for all types 
     1223 
     1224    !Local variables 
     1225    INTEGER :: jfile 
     1226 
     1227    cfiles(jtype,:) = cfilestype(:) 
     1228    cobstypes(jtype) = ctypein 
     1229    ifiles(jtype) = 0 
     1230    DO jfile = 1, jpmaxnfiles 
     1231       IF ( trim(cfiles(jtype,jfile)) /= '' ) & 
     1232                 ifiles(jtype) = ifiles(jtype) + 1 
     1233    END DO 
     1234 
     1235    IF ( ifiles(jtype) == 0 ) THEN 
     1236         CALL ctl_stop( 'Logical for observation type '//TRIM(ctypein)//   & 
     1237            &           ' set to true but no files available to read' ) 
     1238    ENDIF 
     1239 
     1240    IF(lwp) THEN     
     1241       WRITE(numout,*) '             '//cobstypes(jtype)//' input observation file names:' 
     1242       DO jfile = 1, ifiles(jtype) 
     1243          WRITE(numout,*) '                '//TRIM(cfiles(jtype,jfile)) 
     1244       END DO 
     1245    ENDIF 
     1246 
     1247    END SUBROUTINE obs_settypefiles 
     1248 
     1249    SUBROUTINE obs_setinterpopts( ntypes, jtype, ctypein,             & 
     1250               &                  n2dint_default, n2dint_type,        & 
     1251               &                  ravglamscl_type, ravgphiscl_type,   & 
     1252               &                  lfp_indegs_type, lavnight_type,     & 
     1253               &                  n2dint, ravglamscl, ravgphiscl,     & 
     1254               &                  lfpindegs, lavnight ) 
     1255 
     1256    INTEGER, INTENT(IN)  :: ntypes             ! Total number of obs types 
     1257    INTEGER, INTENT(IN)  :: jtype              ! Index of the current type of obs 
     1258    INTEGER, INTENT(IN)  :: n2dint_default     ! Default option for interpolation type 
     1259    INTEGER, INTENT(IN)  :: n2dint_type        ! Option for interpolation type 
     1260    REAL(wp), INTENT(IN) :: & 
     1261       &                    ravglamscl_type, & !E/W diameter of obs footprint for this type 
     1262       &                    ravgphiscl_type    !N/S diameter of obs footprint for this type 
     1263    LOGICAL, INTENT(IN)  :: lfp_indegs_type    !T=> footprint in degrees, F=> in metres 
     1264    LOGICAL, INTENT(IN)  :: lavnight_type      !T=> obs represent night time average 
     1265    CHARACTER(len=6), INTENT(IN) :: ctypein  
     1266 
     1267    INTEGER, DIMENSION(ntypes), INTENT(INOUT) :: & 
     1268       &                    n2dint  
     1269    REAL(wp), DIMENSION(ntypes), INTENT(INOUT) :: & 
     1270       &                    ravglamscl, ravgphiscl 
     1271    LOGICAL, DIMENSION(ntypes), INTENT(INOUT) :: & 
     1272       &                    lfpindegs, lavnight 
     1273 
     1274    lavnight(jtype) = lavnight_type 
     1275 
     1276    IF ( (n2dint_type >= 1) .AND. (n2dint_type <= 6) ) THEN 
     1277       n2dint(jtype) = n2dint_type 
     1278    ELSE 
     1279       n2dint(jtype) = n2dint_default 
     1280    ENDIF 
     1281 
     1282    ! For averaging observation footprints set options for size of footprint  
     1283    IF ( (n2dint(jtype) > 4) .AND. (n2dint(jtype) <= 6) ) THEN 
     1284       IF ( ravglamscl_type > 0._wp ) THEN 
     1285          ravglamscl(jtype) = ravglamscl_type 
     1286       ELSE 
     1287          CALL ctl_stop( 'Incorrect value set for averaging footprint '// & 
     1288                         'scale (ravglamscl) for observation type '//TRIM(ctypein) )       
     1289       ENDIF 
     1290 
     1291       IF ( ravgphiscl_type > 0._wp ) THEN 
     1292          ravgphiscl(jtype) = ravgphiscl_type 
     1293       ELSE 
     1294          CALL ctl_stop( 'Incorrect value set for averaging footprint '// & 
     1295                         'scale (ravgphiscl) for observation type '//TRIM(ctypein) )       
     1296       ENDIF 
     1297 
     1298       lfpindegs(jtype) = lfp_indegs_type  
     1299 
     1300    ENDIF 
     1301 
     1302    ! Write out info  
     1303    IF(lwp) THEN 
     1304       IF ( n2dint(jtype) <= 4 ) THEN 
     1305          WRITE(numout,*) '             '//TRIM(ctypein)// & 
     1306             &            ' model counterparts will be interpolated horizontally' 
     1307       ELSE IF ( n2dint(jtype) <= 6 ) THEN 
     1308          WRITE(numout,*) '             '//TRIM(ctypein)// & 
     1309             &            ' model counterparts will be averaged horizontally' 
     1310          WRITE(numout,*) '             '//'    with E/W scale: ',ravglamscl(jtype) 
     1311          WRITE(numout,*) '             '//'    with N/S scale: ',ravgphiscl(jtype) 
     1312          IF ( lfpindegs(jtype) ) THEN 
     1313              WRITE(numout,*) '             '//'    (in degrees)' 
     1314          ELSE 
     1315              WRITE(numout,*) '             '//'    (in metres)' 
     1316          ENDIF 
     1317       ENDIF 
     1318    ENDIF 
     1319 
     1320    END SUBROUTINE obs_setinterpopts 
     1321 
    16011322END MODULE diaobs 
Note: See TracChangeset for help on using the changeset viewer.