source: branches/UKMO/dev_r5518_obs_oper_update_SlaAssim/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 11863

Last change on this file since 11863 was 11863, checked in by kingr, 2 years ago

Merged recent changes to dev_r5518_obs_oper_update

File size: 94.6 KB
Line 
1MODULE diaobs
2   !!======================================================================
3   !!                       ***  MODULE diaobs  ***
4   !! Observation diagnostics: Computation of the misfit between data and
5   !!                          their model equivalent
6   !!======================================================================
7
8   !!----------------------------------------------------------------------
9   !!   dia_obs_init : Reading and prepare observations
10   !!   dia_obs      : Compute model equivalent to observations
11   !!   dia_obs_wri  : Write observational diagnostics
12   !!   ini_date     : Compute the initial date YYYYMMDD.HHMMSS
13   !!   fin_date     : Compute the final date YYYYMMDD.HHMMSS
14   !!----------------------------------------------------------------------
15   !! * Modules used
16   USE wrk_nemo                 ! Memory Allocation
17   USE par_kind                 ! Precision variables
18   USE in_out_manager           ! I/O manager
19   USE par_oce
20   USE dom_oce                  ! Ocean space and time domain variables
21   USE obs_read_prof            ! Reading and allocation of profile obs
22   USE obs_read_surf            ! Reading and allocation of surface obs
23   USE obs_readmdt              ! Reading and allocation of MDT for SLA.
24   USE obs_prep                 ! Preparation of obs. (grid search etc).
25   USE obs_oper                 ! Observation operators
26   USE obs_write                ! Writing of observation related diagnostics
27   USE obs_grid                 ! Grid searching
28   USE obs_read_altbias         ! Bias treatment for altimeter
29   USE obs_sstbias              ! Bias correction routine for SST
30   USE obs_profiles_def         ! Profile data definitions
31   USE obs_surf_def             ! Surface data definitions
32   USE obs_types                ! Definitions for observation types
33   USE mpp_map                  ! MPP mapping
34   USE lib_mpp                  ! For ctl_warn/stop
35   USE tradmp                   ! For climatological temperature & salinity
36
37   IMPLICIT NONE
38
39   !! * Routine accessibility
40   PRIVATE
41   PUBLIC dia_obs_init, &  ! Initialize and read observations
42      &   dia_obs,      &  ! Compute model equivalent to observations
43      &   dia_obs_wri,  &  ! Write model equivalent to observations
44      &   dia_obs_dealloc  ! Deallocate dia_obs data
45
46   !! * Module variables
47   LOGICAL, PUBLIC :: &
48      &       lk_diaobs = .TRUE.   !: Include this for backwards compatibility at NEMO 3.6.
49   LOGICAL :: ln_diaobs            !: Logical switch for the obs operator
50   LOGICAL :: ln_sstnight          !: Logical switch for night mean SST obs
51   LOGICAL :: ln_default_fp_indegs !: T=> Default obs footprint size specified in degrees, F=> in metres
52   LOGICAL :: ln_sla_fp_indegs     !: T=>     SLA obs footprint size specified in degrees, F=> in metres
53   LOGICAL :: ln_sst_fp_indegs     !: T=>     SST obs footprint size specified in degrees, F=> in metres
54   LOGICAL :: ln_sss_fp_indegs     !: T=>     SSS obs footprint size specified in degrees, F=> in metres
55   LOGICAL :: ln_sic_fp_indegs     !: T=> sea-ice obs footprint size specified in degrees, F=> in metres
56   LOGICAL :: ln_output_clim       !: Logical switch for interpolating and writing T/S climatology
57
58   REAL(wp) :: rn_default_avglamscl !: Default E/W diameter of observation footprint
59   REAL(wp) :: rn_default_avgphiscl !: Default N/S diameter of observation footprint
60   REAL(wp) :: rn_sla_avglamscl     !: E/W diameter of SLA observation footprint
61   REAL(wp) :: rn_sla_avgphiscl     !: N/S diameter of SLA observation footprint
62   REAL(wp) :: rn_sst_avglamscl     !: E/W diameter of SST observation footprint
63   REAL(wp) :: rn_sst_avgphiscl     !: N/S diameter of SST observation footprint
64   REAL(wp) :: rn_sss_avglamscl     !: E/W diameter of SSS observation footprint
65   REAL(wp) :: rn_sss_avgphiscl     !: N/S diameter of SSS observation footprint
66   REAL(wp) :: rn_sic_avglamscl     !: E/W diameter of sea-ice observation footprint
67   REAL(wp) :: rn_sic_avgphiscl     !: N/S diameter of sea-ice observation footprint
68
69   INTEGER :: nn_1dint         !: Vertical interpolation method
70   INTEGER :: nn_2dint_default !: Default horizontal interpolation method
71   INTEGER :: nn_2dint_sla     !: SLA horizontal interpolation method (-1 = default)
72   INTEGER :: nn_2dint_sst     !: SST horizontal interpolation method (-1 = default)
73   INTEGER :: nn_2dint_sss     !: SSS horizontal interpolation method (-1 = default)
74   INTEGER :: nn_2dint_sic     !: Seaice horizontal interpolation method (-1 = default)
75 
76   INTEGER, DIMENSION(imaxavtypes) :: &
77      & nn_profdavtypes      !: Profile data types representing a daily average
78   INTEGER :: nproftypes     !: Number of profile obs types
79   INTEGER :: nsurftypes     !: Number of surface obs types
80   INTEGER, DIMENSION(:), ALLOCATABLE :: &
81      & nvarsprof, &         !: Number of profile variables
82      & nvarssurf            !: Number of surface variables
83   INTEGER, DIMENSION(:), ALLOCATABLE :: &
84      & nextrprof, &         !: Number of profile extra variables
85      & nextrsurf            !: Number of surface extra variables
86   INTEGER, DIMENSION(:), ALLOCATABLE :: &
87      & n2dintsurf           !: Interpolation option for surface variables
88   REAL(wp), DIMENSION(:), ALLOCATABLE :: &
89      & ravglamscl, &        !: E/W diameter of averaging footprint for surface variables
90      & ravgphiscl           !: N/S diameter of averaging footprint for surface variables
91   LOGICAL, DIMENSION(:), ALLOCATABLE :: &
92      & lfpindegs, &         !: T=> surface obs footprint size specified in degrees, F=> in metres
93      & llnightav            !: Logical for calculating night-time averages
94
95   TYPE(obs_surf), PUBLIC, POINTER, DIMENSION(:) :: &
96      & surfdata, &          !: Initial surface data
97      & surfdataqc           !: Surface data after quality control
98   TYPE(obs_prof), PUBLIC, POINTER, DIMENSION(:) :: &
99      & profdata, &          !: Initial profile data
100      & profdataqc           !: Profile data after quality control
101
102   CHARACTER(len=8), PUBLIC, DIMENSION(:), ALLOCATABLE :: &
103      & cobstypesprof, &     !: Profile obs types
104      & cobstypessurf        !: Surface obs types
105
106   !!----------------------------------------------------------------------
107   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
108   !! $Id$
109   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
110   !!----------------------------------------------------------------------
111
112   !! * Substitutions
113#  include "domzgr_substitute.h90"
114CONTAINS
115
116   SUBROUTINE dia_obs_init
117      !!----------------------------------------------------------------------
118      !!                    ***  ROUTINE dia_obs_init  ***
119      !!         
120      !! ** Purpose : Initialize and read observations
121      !!
122      !! ** Method  : Read the namelist and call reading routines
123      !!
124      !! ** Action  : Read the namelist and call reading routines
125      !!
126      !! History :
127      !!        !  06-03  (K. Mogensen) Original code
128      !!        !  06-05  (A. Weaver) Reformatted
129      !!        !  06-10  (A. Weaver) Cleaning and add controls
130      !!        !  07-03  (K. Mogensen) General handling of profiles
131      !!        !  14-08  (J.While) Incorporated SST bias correction
132      !!        !  15-02  (M. Martin) Simplification of namelist and code
133      !!----------------------------------------------------------------------
134
135      IMPLICIT NONE
136
137      !! * Local declarations
138      INTEGER, PARAMETER :: &
139         & jpmaxnfiles = 1000    ! Maximum number of files for each obs type
140      INTEGER, DIMENSION(:), ALLOCATABLE :: &
141         & ifilesprof, &         ! Number of profile files
142         & ifilessurf            ! Number of surface files
143      INTEGER :: ios             ! Local integer output status for namelist read
144      INTEGER :: jtype           ! Counter for obs types
145      INTEGER :: jvar            ! Counter for variables
146      INTEGER :: jfile           ! Counter for files
147      INTEGER :: jnumsstbias     ! Number of SST bias files to read and apply
148      INTEGER :: n2dint_type     ! Local version of nn_2dint*
149
150      CHARACTER(len=128), DIMENSION(jpmaxnfiles) :: &
151         & cn_profbfiles,      & ! T/S profile input filenames
152         & cn_sstfbfiles,      & ! Sea surface temperature input filenames
153         & cn_slafbfiles,      & ! Sea level anomaly input filenames
154         & cn_sicfbfiles,      & ! Seaice concentration input filenames
155         & cn_velfbfiles,      & ! Velocity profile input filenames
156         & cn_sssfbfiles,      & ! Sea surface salinity input filenames
157         & cn_slchltotfbfiles, & ! Surface total              log10(chlorophyll) input filenames
158         & cn_slchldiafbfiles, & ! Surface diatom             log10(chlorophyll) input filenames
159         & cn_slchlnonfbfiles, & ! Surface non-diatom         log10(chlorophyll) input filenames
160         & cn_slchldinfbfiles, & ! Surface dinoflagellate     log10(chlorophyll) input filenames
161         & cn_slchlmicfbfiles, & ! Surface microphytoplankton log10(chlorophyll) input filenames
162         & cn_slchlnanfbfiles, & ! Surface nanophytoplankton  log10(chlorophyll) input filenames
163         & cn_slchlpicfbfiles, & ! Surface picophytoplankton  log10(chlorophyll) input filenames
164         & cn_schltotfbfiles,  & ! Surface total              chlorophyll        input filenames
165         & cn_slphytotfbfiles, & ! Surface total      log10(phytoplankton carbon) input filenames
166         & cn_slphydiafbfiles, & ! Surface diatom     log10(phytoplankton carbon) input filenames
167         & cn_slphynonfbfiles, & ! Surface non-diatom log10(phytoplankton carbon) input filenames
168         & cn_sspmfbfiles,     & ! Surface suspended particulate matter input filenames
169         & cn_skd490fbfiles,   & ! Surface Kd490 input filenames
170         & cn_sfco2fbfiles,    & ! Surface fugacity         of carbon dioxide input filenames
171         & cn_spco2fbfiles,    & ! Surface partial pressure of carbon dioxide input filenames
172         & cn_plchltotfbfiles, & ! Profile total log10(chlorophyll) input filenames
173         & cn_pchltotfbfiles,  & ! Profile total chlorophyll input filenames
174         & cn_pno3fbfiles,     & ! Profile nitrate input filenames
175         & cn_psi4fbfiles,     & ! Profile silicate input filenames
176         & cn_ppo4fbfiles,     & ! Profile phosphate input filenames
177         & cn_pdicfbfiles,     & ! Profile dissolved inorganic carbon input filenames
178         & cn_palkfbfiles,     & ! Profile alkalinity input filenames
179         & cn_pphfbfiles,      & ! Profile pH input filenames
180         & cn_po2fbfiles,      & ! Profile dissolved oxygen input filenames
181         & cn_sstbiasfiles       ! SST bias input filenames
182
183      CHARACTER(LEN=128) :: &
184         & cn_altbiasfile        ! Altimeter bias input filename
185
186
187      LOGICAL :: ln_t3d          ! Logical switch for temperature profiles
188      LOGICAL :: ln_s3d          ! Logical switch for salinity profiles
189      LOGICAL :: ln_sla          ! Logical switch for sea level anomalies
190      LOGICAL :: ln_sst          ! Logical switch for sea surface temperature
191      LOGICAL :: ln_sic          ! Logical switch for sea ice concentration
192      LOGICAL :: ln_sss          ! Logical switch for sea surface salinity obs
193      LOGICAL :: ln_vel3d        ! Logical switch for velocity (u,v) obs
194      LOGICAL :: ln_slchltot     ! Logical switch for surface total              log10(chlorophyll) obs
195      LOGICAL :: ln_slchldia     ! Logical switch for surface diatom             log10(chlorophyll) obs
196      LOGICAL :: ln_slchlnon     ! Logical switch for surface non-diatom         log10(chlorophyll) obs
197      LOGICAL :: ln_slchldin     ! Logical switch for surface dinoflagellate     log10(chlorophyll) obs
198      LOGICAL :: ln_slchlmic     ! Logical switch for surface microphytoplankton log10(chlorophyll) obs
199      LOGICAL :: ln_slchlnan     ! Logical switch for surface nanophytoplankton  log10(chlorophyll) obs
200      LOGICAL :: ln_slchlpic     ! Logical switch for surface picophytoplankton  log10(chlorophyll) obs
201      LOGICAL :: ln_schltot      ! Logical switch for surface total              chlorophyll        obs
202      LOGICAL :: ln_slphytot     ! Logical switch for surface total      log10(phytoplankton carbon) obs
203      LOGICAL :: ln_slphydia     ! Logical switch for surface diatom     log10(phytoplankton carbon) obs
204      LOGICAL :: ln_slphynon     ! Logical switch for surface non-diatom log10(phytoplankton carbon) obs
205      LOGICAL :: ln_sspm         ! Logical switch for surface suspended particulate matter obs
206      LOGICAL :: ln_skd490       ! Logical switch for surface Kd490
207      LOGICAL :: ln_sfco2        ! Logical switch for surface fugacity         of carbon dioxide obs
208      LOGICAL :: ln_spco2        ! Logical switch for surface partial pressure of carbon dioxide obs
209      LOGICAL :: ln_plchltot     ! Logical switch for profile total log10(chlorophyll) obs
210      LOGICAL :: ln_pchltot      ! Logical switch for profile total chlorophyll obs
211      LOGICAL :: ln_pno3         ! Logical switch for profile nitrate obs
212      LOGICAL :: ln_psi4         ! Logical switch for profile silicate obs
213      LOGICAL :: ln_ppo4         ! Logical switch for profile phosphate obs
214      LOGICAL :: ln_pdic         ! Logical switch for profile dissolved inorganic carbon obs
215      LOGICAL :: ln_palk         ! Logical switch for profile alkalinity obs
216      LOGICAL :: ln_pph          ! Logical switch for profile pH obs
217      LOGICAL :: ln_po2          ! Logical switch for profile dissolved oxygen obs
218      LOGICAL :: ln_nea          ! Logical switch to remove obs near land
219      LOGICAL :: ln_altbias      ! Logical switch for altimeter bias
220      LOGICAL :: ln_sstbias      ! Logical switch for bias correction of SST
221      LOGICAL :: ln_ignmis       ! Logical switch for ignoring missing files
222      LOGICAL :: ln_s_at_t       ! Logical switch to compute model S at T obs
223      LOGICAL :: ln_bound_reject ! Logical switch for rejecting obs near the boundary
224
225      REAL(dp) :: rn_dobsini     ! Obs window start date YYYYMMDD.HHMMSS
226      REAL(dp) :: rn_dobsend     ! Obs window end date   YYYYMMDD.HHMMSS
227
228      REAL(wp) :: ztype_avglamscl ! Local version of rn_*_avglamscl
229      REAL(wp) :: ztype_avgphiscl ! Local version of rn_*_avgphiscl
230
231      CHARACTER(len=128), DIMENSION(:,:), ALLOCATABLE :: &
232         & clproffiles, &        ! Profile filenames
233         & clsurffiles           ! Surface filenames
234      CHARACTER(len=8), DIMENSION(:), ALLOCATABLE :: clvars  ! Expected variable names
235
236      LOGICAL, DIMENSION(:), ALLOCATABLE :: llvar   ! Logical for profile variable read
237      LOGICAL :: ltype_fp_indegs ! Local version of ln_*_fp_indegs
238      LOGICAL :: ltype_night     ! Local version of ln_sstnight (false for other variables)
239      LOGICAL :: ltype_clim      ! Local version of ln_output_clim
240
241      REAL(wp), POINTER, DIMENSION(:,:,:) :: &
242         & zglam                 ! Model longitudes for profile variables
243      REAL(wp), POINTER, DIMENSION(:,:,:) :: &
244         & zgphi                 ! Model latitudes for profile variables
245      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: &
246         & zmask                 ! Model land/sea mask associated with variables
247
248
249      NAMELIST/namobs/ln_diaobs, ln_t3d, ln_s3d, ln_sla,              &
250         &            ln_sst, ln_sic, ln_sss, ln_vel3d,               &
251         &            ln_slchltot, ln_slchldia, ln_slchlnon,          &
252         &            ln_slchldin, ln_slchlmic, ln_slchlnan,          &
253         &            ln_slchlpic, ln_schltot,                        &
254         &            ln_slphytot, ln_slphydia, ln_slphynon,          &
255         &            ln_sspm,     ln_sfco2,    ln_spco2,             &
256         &            ln_skd490,                                      &
257         &            ln_plchltot, ln_pchltot,  ln_pno3,              &
258         &            ln_psi4,     ln_ppo4,     ln_pdic,              &
259         &            ln_palk,     ln_pph,      ln_po2,               &
260         &            ln_altbias, ln_sstbias, ln_nea,                 &
261         &            ln_grid_global, ln_grid_search_lookup,          &
262         &            ln_ignmis, ln_s_at_t, ln_bound_reject,          &
263         &            ln_sstnight,  ln_output_clim,                   &
264         &            ln_default_fp_indegs,                           &
265         &            ln_sla_fp_indegs, ln_sst_fp_indegs,             &
266         &            ln_sss_fp_indegs, ln_sic_fp_indegs,             &
267         &            cn_profbfiles, cn_slafbfiles,                   &
268         &            cn_sstfbfiles, cn_sicfbfiles,                   &
269         &            cn_velfbfiles, cn_sssfbfiles,                   &
270         &            cn_slchltotfbfiles, cn_slchldiafbfiles,         &
271         &            cn_slchlnonfbfiles, cn_slchldinfbfiles,         &
272         &            cn_slchlmicfbfiles, cn_slchlnanfbfiles,         &
273         &            cn_slchlpicfbfiles, cn_schltotfbfiles,          &
274         &            cn_slphytotfbfiles, cn_slphydiafbfiles,         &
275         &            cn_slphynonfbfiles, cn_sspmfbfiles,             &
276         &            cn_skd490fbfiles,                               &
277         &            cn_sfco2fbfiles, cn_spco2fbfiles,               &
278         &            cn_plchltotfbfiles, cn_pchltotfbfiles,          &
279         &            cn_pno3fbfiles, cn_psi4fbfiles, cn_ppo4fbfiles, &
280         &            cn_pdicfbfiles, cn_palkfbfiles, cn_pphfbfiles,  &
281         &            cn_po2fbfiles,                                  &
282         &            cn_sstbiasfiles, cn_altbiasfile,                &
283         &            cn_gridsearchfile, rn_gridsearchres,            &
284         &            rn_dobsini, rn_dobsend,                         &
285         &            rn_default_avglamscl, rn_default_avgphiscl,     &
286         &            rn_sla_avglamscl, rn_sla_avgphiscl,             &
287         &            rn_sst_avglamscl, rn_sst_avgphiscl,             &
288         &            rn_sss_avglamscl, rn_sss_avgphiscl,             &
289         &            rn_sic_avglamscl, rn_sic_avgphiscl,             &
290         &            nn_1dint, nn_2dint_default,                     &
291         &            nn_2dint_sla, nn_2dint_sst,                     &
292         &            nn_2dint_sss, nn_2dint_sic,                     &
293         &            nn_msshc, rn_mdtcorr, rn_mdtcutoff,             &
294         &            nn_profdavtypes
295
296      !-----------------------------------------------------------------------
297      ! Read namelist parameters
298      !-----------------------------------------------------------------------
299
300      ! Some namelist arrays need initialising
301      cn_profbfiles(:)      = ''
302      cn_slafbfiles(:)      = ''
303      cn_sstfbfiles(:)      = ''
304      cn_sicfbfiles(:)      = ''
305      cn_velfbfiles(:)      = ''
306      cn_sssfbfiles(:)      = ''
307      cn_slchltotfbfiles(:) = ''
308      cn_slchldiafbfiles(:) = ''
309      cn_slchlnonfbfiles(:) = ''
310      cn_slchldinfbfiles(:) = ''
311      cn_slchlmicfbfiles(:) = ''
312      cn_slchlnanfbfiles(:) = ''
313      cn_slchlpicfbfiles(:) = ''
314      cn_schltotfbfiles(:)  = ''
315      cn_slphytotfbfiles(:) = ''
316      cn_slphydiafbfiles(:) = ''
317      cn_slphynonfbfiles(:) = ''
318      cn_sspmfbfiles(:)     = ''
319      cn_skd490fbfiles(:)   = ''
320      cn_sfco2fbfiles(:)    = ''
321      cn_spco2fbfiles(:)    = ''
322      cn_plchltotfbfiles(:) = ''
323      cn_pchltotfbfiles(:)  = ''
324      cn_pno3fbfiles(:)     = ''
325      cn_psi4fbfiles(:)     = ''
326      cn_ppo4fbfiles(:)     = ''
327      cn_pdicfbfiles(:)     = ''
328      cn_palkfbfiles(:)     = ''
329      cn_pphfbfiles(:)      = ''
330      cn_po2fbfiles(:)      = ''
331      cn_sstbiasfiles(:)    = ''
332      nn_profdavtypes(:)    = -1
333
334      CALL ini_date( rn_dobsini )
335      CALL fin_date( rn_dobsend )
336
337      ! Read namelist namobs : control observation diagnostics
338      REWIND( numnam_ref )   ! Namelist namobs in reference namelist
339      READ  ( numnam_ref, namobs, IOSTAT = ios, ERR = 901)
340901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in reference namelist', lwp )
341
342      REWIND( numnam_cfg )   ! Namelist namobs in configuration namelist
343      READ  ( numnam_cfg, namobs, IOSTAT = ios, ERR = 902 )
344902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in configuration namelist', lwp )
345      IF(lwm) WRITE ( numond, namobs )
346
347      lk_diaobs = .FALSE.
348#if defined key_diaobs
349      IF ( ln_diaobs ) lk_diaobs = .TRUE.
350#endif
351
352      IF ( .NOT. lk_diaobs ) THEN
353         IF(lwp) WRITE(numout,cform_war)
354         IF(lwp) WRITE(numout,*)' ln_diaobs is set to false or key_diaobs is not set, so not calling dia_obs'
355         RETURN
356      ENDIF
357
358      IF(lwp) THEN
359         WRITE(numout,*)
360         WRITE(numout,*) 'dia_obs_init : Observation diagnostic initialization'
361         WRITE(numout,*) '~~~~~~~~~~~~'
362         WRITE(numout,*) '          Namelist namobs : set observation diagnostic parameters' 
363         WRITE(numout,*) '             Logical switch for T profile observations                ln_t3d = ', ln_t3d
364         WRITE(numout,*) '             Logical switch for S profile observations                ln_s3d = ', ln_s3d
365         WRITE(numout,*) '             Logical switch for SLA observations                      ln_sla = ', ln_sla
366         WRITE(numout,*) '             Logical switch for SST observations                      ln_sst = ', ln_sst
367         WRITE(numout,*) '             Logical switch for Sea Ice observations                  ln_sic = ', ln_sic
368         WRITE(numout,*) '             Logical switch for velocity observations               ln_vel3d = ', ln_vel3d
369         WRITE(numout,*) '             Logical switch for SSS observations                      ln_sss = ', ln_sss
370         WRITE(numout,*) '             Logical switch for surface total logchl obs         ln_slchltot = ', ln_slchltot
371         WRITE(numout,*) '             Logical switch for surface diatom logchl obs        ln_slchldia = ', ln_slchldia
372         WRITE(numout,*) '             Logical switch for surface non-diatom logchl obs    ln_slchlnon = ', ln_slchlnon
373         WRITE(numout,*) '             Logical switch for surface dino logchl obs          ln_slchldin = ', ln_slchldin
374         WRITE(numout,*) '             Logical switch for surface micro logchl obs         ln_slchlmic = ', ln_slchlmic
375         WRITE(numout,*) '             Logical switch for surface nano logchl obs          ln_slchlnan = ', ln_slchlnan
376         WRITE(numout,*) '             Logical switch for surface pico logchl obs          ln_slchlpic = ', ln_slchlpic
377         WRITE(numout,*) '             Logical switch for surface total chl obs             ln_schltot = ', ln_schltot
378         WRITE(numout,*) '             Logical switch for surface total log(phyC) obs      ln_slphytot = ', ln_slphytot
379         WRITE(numout,*) '             Logical switch for surface diatom log(phyC) obs     ln_slphydia = ', ln_slphydia
380         WRITE(numout,*) '             Logical switch for surface non-diatom log(phyC) obs ln_slphynon = ', ln_slphynon
381         WRITE(numout,*) '             Logical switch for surface SPM observations             ln_sspm = ', ln_sspm
382         WRITE(numout,*) '             Logical switch for surface Kd490 observations         ln_skd490 = ', ln_skd490
383         WRITE(numout,*) '             Logical switch for surface fCO2 observations           ln_sfco2 = ', ln_sfco2
384         WRITE(numout,*) '             Logical switch for surface pCO2 observations           ln_spco2 = ', ln_spco2
385         WRITE(numout,*) '             Logical switch for profile total logchl obs         ln_plchltot = ', ln_plchltot
386         WRITE(numout,*) '             Logical switch for profile total chl obs             ln_pchltot = ', ln_pchltot
387         WRITE(numout,*) '             Logical switch for profile nitrate obs                  ln_pno3 = ', ln_pno3
388         WRITE(numout,*) '             Logical switch for profile silicate obs                 ln_psi4 = ', ln_psi4
389         WRITE(numout,*) '             Logical switch for profile phosphate obs                ln_ppo4 = ', ln_ppo4
390         WRITE(numout,*) '             Logical switch for profile DIC obs                      ln_pdic = ', ln_pdic
391         WRITE(numout,*) '             Logical switch for profile alkalinity obs               ln_palk = ', ln_palk
392         WRITE(numout,*) '             Logical switch for profile pH obs                        ln_pph = ', ln_pph
393         WRITE(numout,*) '             Logical switch for profile oxygen obs                    ln_po2 = ', ln_po2
394         WRITE(numout,*) '             Global distribution of observations              ln_grid_global = ', ln_grid_global
395         WRITE(numout,*) '             Logical switch for obs grid search lookup ln_grid_search_lookup = ', ln_grid_search_lookup
396         IF (ln_grid_search_lookup) &
397            WRITE(numout,*) '             Grid search lookup file header                cn_gridsearchfile = ', cn_gridsearchfile
398         WRITE(numout,*) '             Initial date in window YYYYMMDD.HHMMSS               rn_dobsini = ', rn_dobsini
399         WRITE(numout,*) '             Final date in window YYYYMMDD.HHMMSS                 rn_dobsend = ', rn_dobsend
400         WRITE(numout,*) '             Type of vertical interpolation method                  nn_1dint = ', nn_1dint
401         WRITE(numout,*) '             Default horizontal interpolation method        nn_2dint_default = ', nn_2dint_default
402         WRITE(numout,*) '             Type of horizontal interpolation method for SLA    nn_2dint_sla = ', nn_2dint_sla
403         WRITE(numout,*) '             Type of horizontal interpolation method for SST    nn_2dint_sst = ', nn_2dint_sst
404         WRITE(numout,*) '             Type of horizontal interpolation method for SSS    nn_2dint_sss = ', nn_2dint_sss
405         WRITE(numout,*) '             Type of horizontal interpolation method for SIC    nn_2dint_sic = ', nn_2dint_sic
406         WRITE(numout,*) '             Default E/W diameter of obs footprint      rn_default_avglamscl = ', rn_default_avglamscl
407         WRITE(numout,*) '             Default N/S diameter of obs footprint      rn_default_avgphiscl = ', rn_default_avgphiscl
408         WRITE(numout,*) '             Default obs footprint in deg [T] or m [F]  ln_default_fp_indegs = ', ln_default_fp_indegs
409         WRITE(numout,*) '             SLA E/W diameter of obs footprint              rn_sla_avglamscl = ', rn_sla_avglamscl
410         WRITE(numout,*) '             SLA N/S diameter of obs footprint              rn_sla_avgphiscl = ', rn_sla_avgphiscl
411         WRITE(numout,*) '             SLA obs footprint in deg [T] or m [F]          ln_sla_fp_indegs = ', ln_sla_fp_indegs
412         WRITE(numout,*) '             SST E/W diameter of obs footprint              rn_sst_avglamscl = ', rn_sst_avglamscl
413         WRITE(numout,*) '             SST N/S diameter of obs footprint              rn_sst_avgphiscl = ', rn_sst_avgphiscl
414         WRITE(numout,*) '             SST obs footprint in deg [T] or m [F]          ln_sst_fp_indegs = ', ln_sst_fp_indegs
415         WRITE(numout,*) '             SIC E/W diameter of obs footprint              rn_sic_avglamscl = ', rn_sic_avglamscl
416         WRITE(numout,*) '             SIC N/S diameter of obs footprint              rn_sic_avgphiscl = ', rn_sic_avgphiscl
417         WRITE(numout,*) '             SIC obs footprint in deg [T] or m [F]          ln_sic_fp_indegs = ', ln_sic_fp_indegs
418         WRITE(numout,*) '             Rejection of observations near land switch               ln_nea = ', ln_nea
419         WRITE(numout,*) '             Rejection of obs near open bdys                 ln_bound_reject = ', ln_bound_reject
420         WRITE(numout,*) '             MSSH correction scheme                                 nn_msshc = ', nn_msshc
421         WRITE(numout,*) '             MDT  correction                                      rn_mdtcorr = ', rn_mdtcorr
422         WRITE(numout,*) '             MDT cutoff for computed correction                 rn_mdtcutoff = ', rn_mdtcutoff
423         WRITE(numout,*) '             Logical switch for alt bias                          ln_altbias = ', ln_altbias
424         WRITE(numout,*) '             Logical switch for sst bias                          ln_sstbias = ', ln_sstbias
425         WRITE(numout,*) '             Logical switch for ignoring missing files             ln_ignmis = ', ln_ignmis
426         WRITE(numout,*) '             Daily average types                             nn_profdavtypes = ', nn_profdavtypes
427         WRITE(numout,*) '             Logical switch for night-time SST obs               ln_sstnight = ', ln_sstnight
428         WRITE(numout,*) '             Logical switch for writing climat. at obs points ln_output_clim = ', ln_output_clim
429      ENDIF
430      !-----------------------------------------------------------------------
431      ! Set up list of observation types to be used
432      ! and the files associated with each type
433      !-----------------------------------------------------------------------
434
435      nproftypes = COUNT( (/ln_t3d .OR. ln_s3d, ln_vel3d, ln_plchltot,          &
436         &                  ln_pchltot,  ln_pno3,     ln_psi4,     ln_ppo4,     &
437         &                  ln_pdic,     ln_palk,     ln_pph,      ln_po2 /) )
438      nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic, ln_sss,                     &
439         &                  ln_slchltot, ln_slchldia, ln_slchlnon, ln_slchldin, &
440         &                  ln_slchlmic, ln_slchlnan, ln_slchlpic, ln_schltot,  &
441         &                  ln_slphytot, ln_slphydia, ln_slphynon, ln_sspm,     &
442         &                  ln_skd490,   ln_sfco2,    ln_spco2 /) )
443
444      IF ( nproftypes == 0 .AND. nsurftypes == 0 ) THEN
445         IF(lwp) WRITE(numout,cform_war)
446         IF(lwp) WRITE(numout,*) ' ln_diaobs is set to true, but all obs operator logical flags', &
447            &                    ' are set to .FALSE. so turning off calls to dia_obs'
448         nwarn = nwarn + 1
449         lk_diaobs = .FALSE.
450         RETURN
451      ENDIF
452
453      IF ( ln_output_clim .AND. ( .NOT. ln_tradmp ) ) THEN
454         IF(lwp) WRITE(numout,cform_war)
455         IF(lwp) WRITE(numout,*) ' ln_output_clim is true, but ln_tradmp is false', &
456            &                    ' so climatological T/S not available and will not be output'
457         nwarn = nwarn + 1
458         ln_output_clim = .FALSE.
459      ENDIF
460     
461
462      IF(lwp) WRITE(numout,*) '          Number of profile obs types: ',nproftypes
463      IF ( nproftypes > 0 ) THEN
464
465         ALLOCATE( cobstypesprof(nproftypes) )
466         ALLOCATE( ifilesprof(nproftypes) )
467         ALLOCATE( clproffiles(nproftypes,jpmaxnfiles) )
468
469         jtype = 0
470         IF (ln_t3d .OR. ln_s3d) THEN
471            jtype = jtype + 1
472            cobstypesprof(jtype) = 'prof'
473            clproffiles(jtype,:) = cn_profbfiles
474         ENDIF
475         IF (ln_vel3d) THEN
476            jtype = jtype + 1
477            cobstypesprof(jtype) =  'vel'
478            clproffiles(jtype,:) = cn_velfbfiles
479         ENDIF
480         IF (ln_plchltot) THEN
481            jtype = jtype + 1
482            cobstypesprof(jtype) = 'plchltot'
483            clproffiles(jtype,:) = cn_plchltotfbfiles
484         ENDIF
485         IF (ln_pchltot) THEN
486            jtype = jtype + 1
487            cobstypesprof(jtype) = 'pchltot'
488            clproffiles(jtype,:) = cn_pchltotfbfiles
489         ENDIF
490         IF (ln_pno3) THEN
491            jtype = jtype + 1
492            cobstypesprof(jtype) = 'pno3'
493            clproffiles(jtype,:) = cn_pno3fbfiles
494         ENDIF
495         IF (ln_psi4) THEN
496            jtype = jtype + 1
497            cobstypesprof(jtype) = 'psi4'
498            clproffiles(jtype,:) = cn_psi4fbfiles
499         ENDIF
500         IF (ln_ppo4) THEN
501            jtype = jtype + 1
502            cobstypesprof(jtype) = 'ppo4'
503            clproffiles(jtype,:) = cn_ppo4fbfiles
504         ENDIF
505         IF (ln_pdic) THEN
506            jtype = jtype + 1
507            cobstypesprof(jtype) = 'pdic'
508            clproffiles(jtype,:) = cn_pdicfbfiles
509         ENDIF
510         IF (ln_palk) THEN
511            jtype = jtype + 1
512            cobstypesprof(jtype) = 'palk'
513            clproffiles(jtype,:) = cn_palkfbfiles
514         ENDIF
515         IF (ln_pph) THEN
516            jtype = jtype + 1
517            cobstypesprof(jtype) = 'pph'
518            clproffiles(jtype,:) = cn_pphfbfiles
519         ENDIF
520         IF (ln_po2) THEN
521            jtype = jtype + 1
522            cobstypesprof(jtype) = 'po2'
523            clproffiles(jtype,:) = cn_po2fbfiles
524         ENDIF
525
526         CALL obs_settypefiles( nproftypes, jpmaxnfiles, ifilesprof, cobstypesprof, clproffiles )
527
528      ENDIF
529
530      IF(lwp) WRITE(numout,*)'          Number of surface obs types: ',nsurftypes
531      IF ( nsurftypes > 0 ) THEN
532
533         ALLOCATE( cobstypessurf(nsurftypes) )
534         ALLOCATE( ifilessurf(nsurftypes) )
535         ALLOCATE( clsurffiles(nsurftypes, jpmaxnfiles) )
536         ALLOCATE(n2dintsurf(nsurftypes))
537         ALLOCATE(ravglamscl(nsurftypes))
538         ALLOCATE(ravgphiscl(nsurftypes))
539         ALLOCATE(lfpindegs(nsurftypes))
540         ALLOCATE(llnightav(nsurftypes))
541
542         jtype = 0
543         IF (ln_sla) THEN
544            jtype = jtype + 1
545            cobstypessurf(jtype) = 'sla'
546            clsurffiles(jtype,:) = cn_slafbfiles
547         ENDIF
548         IF (ln_sst) THEN
549            jtype = jtype + 1
550            cobstypessurf(jtype) = 'sst'
551            clsurffiles(jtype,:) = cn_sstfbfiles
552         ENDIF
553         IF (ln_sic) THEN
554            jtype = jtype + 1
555            cobstypessurf(jtype) = 'sic'
556            clsurffiles(jtype,:) = cn_sicfbfiles
557         ENDIF
558         IF (ln_sss) THEN
559            jtype = jtype + 1
560            cobstypessurf(jtype) = 'sss'
561            clsurffiles(jtype,:) = cn_sssfbfiles
562         ENDIF
563         IF (ln_slchltot) THEN
564            jtype = jtype + 1
565            cobstypessurf(jtype) = 'slchltot'
566            clsurffiles(jtype,:) = cn_slchltotfbfiles
567         ENDIF
568         IF (ln_slchldia) THEN
569            jtype = jtype + 1
570            cobstypessurf(jtype) = 'slchldia'
571            clsurffiles(jtype,:) = cn_slchldiafbfiles
572         ENDIF
573         IF (ln_slchlnon) THEN
574            jtype = jtype + 1
575            cobstypessurf(jtype) = 'slchlnon'
576            clsurffiles(jtype,:) = cn_slchlnonfbfiles
577         ENDIF
578         IF (ln_slchldin) THEN
579            jtype = jtype + 1
580            cobstypessurf(jtype) = 'slchldin'
581            clsurffiles(jtype,:) = cn_slchldinfbfiles
582         ENDIF
583         IF (ln_slchlmic) THEN
584            jtype = jtype + 1
585            cobstypessurf(jtype) = 'slchlmic'
586            clsurffiles(jtype,:) = cn_slchlmicfbfiles
587         ENDIF
588         IF (ln_slchlnan) THEN
589            jtype = jtype + 1
590            cobstypessurf(jtype) = 'slchlnan'
591            clsurffiles(jtype,:) = cn_slchlnanfbfiles
592         ENDIF
593         IF (ln_slchlpic) THEN
594            jtype = jtype + 1
595            cobstypessurf(jtype) = 'slchlpic'
596            clsurffiles(jtype,:) = cn_slchlpicfbfiles
597         ENDIF
598         IF (ln_schltot) THEN
599            jtype = jtype + 1
600            cobstypessurf(jtype) = 'schltot'
601            clsurffiles(jtype,:) = cn_schltotfbfiles
602         ENDIF
603         IF (ln_slphytot) THEN
604            jtype = jtype + 1
605            cobstypessurf(jtype) = 'slphytot'
606            clsurffiles(jtype,:) = cn_slphytotfbfiles
607         ENDIF
608         IF (ln_slphydia) THEN
609            jtype = jtype + 1
610            cobstypessurf(jtype) = 'slphydia'
611            clsurffiles(jtype,:) = cn_slphydiafbfiles
612         ENDIF
613         IF (ln_slphynon) THEN
614            jtype = jtype + 1
615            cobstypessurf(jtype) = 'slphynon'
616            clsurffiles(jtype,:) = cn_slphynonfbfiles
617         ENDIF
618         IF (ln_sspm) THEN
619            jtype = jtype + 1
620            cobstypessurf(jtype) = 'sspm'
621            clsurffiles(jtype,:) = cn_sspmfbfiles
622         ENDIF
623         IF (ln_skd490) THEN
624            jtype = jtype + 1
625            cobstypessurf(jtype) = 'skd490'
626            clsurffiles(jtype,:) = cn_skd490fbfiles
627         ENDIF
628         IF (ln_sfco2) THEN
629            jtype = jtype + 1
630            cobstypessurf(jtype) = 'sfco2'
631            clsurffiles(jtype,:) = cn_sfco2fbfiles
632         ENDIF
633         IF (ln_spco2) THEN
634            jtype = jtype + 1
635            cobstypessurf(jtype) = 'spco2'
636            clsurffiles(jtype,:) = cn_spco2fbfiles
637         ENDIF
638
639         CALL obs_settypefiles( nsurftypes, jpmaxnfiles, ifilessurf, cobstypessurf, clsurffiles )
640
641         DO jtype = 1, nsurftypes
642
643            IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN
644               IF ( nn_2dint_sla == -1 ) THEN
645                  n2dint_type  = nn_2dint_default
646               ELSE
647                  n2dint_type  = nn_2dint_sla
648               ENDIF
649               ztype_avglamscl = rn_sla_avglamscl
650               ztype_avgphiscl = rn_sla_avgphiscl
651               ltype_fp_indegs = ln_sla_fp_indegs
652               ltype_night     = .FALSE.
653            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sst' ) THEN
654               IF ( nn_2dint_sst == -1 ) THEN
655                  n2dint_type  = nn_2dint_default
656               ELSE
657                  n2dint_type  = nn_2dint_sst
658               ENDIF
659               ztype_avglamscl = rn_sst_avglamscl
660               ztype_avgphiscl = rn_sst_avgphiscl
661               ltype_fp_indegs = ln_sst_fp_indegs
662               ltype_night     = ln_sstnight
663            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sic' ) THEN
664               IF ( nn_2dint_sic == -1 ) THEN
665                  n2dint_type  = nn_2dint_default
666               ELSE
667                  n2dint_type  = nn_2dint_sic
668               ENDIF
669               ztype_avglamscl = rn_sic_avglamscl
670               ztype_avgphiscl = rn_sic_avgphiscl
671               ltype_fp_indegs = ln_sic_fp_indegs
672               ltype_night     = .FALSE.
673            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sss' ) THEN
674               IF ( nn_2dint_sss == -1 ) THEN
675                  n2dint_type  = nn_2dint_default
676               ELSE
677                  n2dint_type  = nn_2dint_sss
678               ENDIF
679               ztype_avglamscl = rn_sss_avglamscl
680               ztype_avgphiscl = rn_sss_avgphiscl
681               ltype_fp_indegs = ln_sss_fp_indegs
682               ltype_night     = .FALSE.
683            ELSE
684               n2dint_type     = nn_2dint_default
685               ztype_avglamscl = rn_default_avglamscl
686               ztype_avgphiscl = rn_default_avgphiscl
687               ltype_fp_indegs = ln_default_fp_indegs
688               ltype_night     = .FALSE.
689            ENDIF
690           
691            CALL obs_setinterpopts( nsurftypes, jtype, TRIM(cobstypessurf(jtype)), &
692               &                    nn_2dint_default, n2dint_type,                 &
693               &                    ztype_avglamscl, ztype_avgphiscl,              &
694               &                    ltype_fp_indegs, ltype_night,                  &
695               &                    n2dintsurf, ravglamscl, ravgphiscl,            &
696               &                    lfpindegs, llnightav )
697
698         END DO
699
700      ENDIF
701
702      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
703
704
705      !-----------------------------------------------------------------------
706      ! Obs operator parameter checking and initialisations
707      !-----------------------------------------------------------------------
708
709      IF ( ln_vel3d .AND. ( .NOT. ln_grid_global ) ) THEN
710         CALL ctl_stop( 'Velocity data only works with ln_grid_global=.true.' )
711         RETURN
712      ENDIF
713
714      IF ( ( nn_1dint < 0 ) .OR. ( nn_1dint > 1 ) ) THEN
715         CALL ctl_stop(' Choice of vertical (1D) interpolation method', &
716            &                    ' is not available')
717      ENDIF
718
719      IF ( ( nn_2dint_default < 0 ) .OR. ( nn_2dint_default > 6 ) ) THEN
720         CALL ctl_stop(' Choice of default horizontal (2D) interpolation method', &
721            &                    ' is not available')
722      ENDIF
723
724      CALL obs_typ_init
725
726      CALL mppmap_init
727
728      CALL obs_grid_setup( )
729
730      !-----------------------------------------------------------------------
731      ! Depending on switches read the various observation types
732      !-----------------------------------------------------------------------
733
734      IF ( nproftypes > 0 ) THEN
735
736         ALLOCATE(profdata(nproftypes))
737         ALLOCATE(profdataqc(nproftypes))
738         ALLOCATE(nvarsprof(nproftypes))
739         ALLOCATE(nextrprof(nproftypes))
740
741         DO jtype = 1, nproftypes
742           
743            ltype_clim = .FALSE. 
744           
745            IF ( TRIM(cobstypesprof(jtype)) == 'prof' ) THEN
746               nvarsprof(jtype) = 2
747               nextrprof(jtype) = 1
748               IF ( ln_output_clim ) ltype_clim = .TRUE.             
749               ALLOCATE(llvar(nvarsprof(jtype)))
750               ALLOCATE(clvars(nvarsprof(jtype)))
751               CALL wrk_alloc( jpi, jpj,      nvarsprof(jtype), zglam )
752               CALL wrk_alloc( jpi, jpj,      nvarsprof(jtype), zgphi )
753               CALL wrk_alloc( jpi, jpj, jpk, nvarsprof(jtype), zmask )
754               llvar(1)       = ln_t3d
755               llvar(2)       = ln_s3d
756               clvars(1)      = 'POTM'
757               clvars(2)      = 'PSAL'
758               zglam(:,:,1)   = glamt(:,:)
759               zglam(:,:,2)   = glamt(:,:)
760               zgphi(:,:,1)   = gphit(:,:)
761               zgphi(:,:,2)   = gphit(:,:)
762               zmask(:,:,:,1) = tmask(:,:,:)
763               zmask(:,:,:,2) = tmask(:,:,:)
764            ELSE IF ( TRIM(cobstypesprof(jtype)) == 'vel' )  THEN
765               nvarsprof(jtype) = 2
766               nextrprof(jtype) = 2
767               ALLOCATE(llvar(nvarsprof(jtype)))
768               ALLOCATE(clvars(nvarsprof(jtype)))
769               CALL wrk_alloc( jpi, jpj,      nvarsprof(jtype), zglam )
770               CALL wrk_alloc( jpi, jpj,      nvarsprof(jtype), zgphi )
771               CALL wrk_alloc( jpi, jpj, jpk, nvarsprof(jtype), zmask )
772               llvar(1)       = ln_vel3d
773               llvar(2)       = ln_vel3d
774               clvars(1)      = 'UVEL'
775               clvars(2)      = 'VVEL'
776               zglam(:,:,1)   = glamu(:,:)
777               zglam(:,:,2)   = glamv(:,:)
778               zgphi(:,:,1)   = gphiu(:,:)
779               zgphi(:,:,2)   = gphiv(:,:)
780               zmask(:,:,:,1) = umask(:,:,:)
781               zmask(:,:,:,2) = vmask(:,:,:)
782            ELSE
783               nvarsprof(jtype) = 1
784               nextrprof(jtype) = 0
785               ALLOCATE(llvar(nvarsprof(jtype)))
786               ALLOCATE(clvars(nvarsprof(jtype)))
787               CALL wrk_alloc( jpi, jpj,      nvarsprof(jtype), zglam )
788               CALL wrk_alloc( jpi, jpj,      nvarsprof(jtype), zgphi )
789               CALL wrk_alloc( jpi, jpj, jpk, nvarsprof(jtype), zmask )
790               llvar(1)       = .TRUE.
791               zglam(:,:,1)   = glamt(:,:)
792               zgphi(:,:,1)   = gphit(:,:)
793               zmask(:,:,:,1) = tmask(:,:,:)
794               IF ( TRIM(cobstypesprof(jtype)) == 'plchltot' )  THEN
795                  clvars(1) = 'PLCHLTOT'
796               ELSE IF ( TRIM(cobstypesprof(jtype)) == 'pchltot' )  THEN
797                  clvars(1) = 'PCHLTOT'
798               ELSE IF ( TRIM(cobstypesprof(jtype)) == 'pno3' )  THEN
799                  clvars(1) = 'PNO3'
800               ELSE IF ( TRIM(cobstypesprof(jtype)) == 'psi4' )  THEN
801                  clvars(1) = 'PSI4'
802               ELSE IF ( TRIM(cobstypesprof(jtype)) == 'ppo4' )  THEN
803                  clvars(1) = 'PPO4'
804               ELSE IF ( TRIM(cobstypesprof(jtype)) == 'pdic' )  THEN
805                  clvars(1) = 'PDIC'
806               ELSE IF ( TRIM(cobstypesprof(jtype)) == 'palk' )  THEN
807                  clvars(1) = 'PALK'
808               ELSE IF ( TRIM(cobstypesprof(jtype)) == 'pph' )  THEN
809                  clvars(1) = 'PPH'
810               ELSE IF ( TRIM(cobstypesprof(jtype)) == 'po2' )  THEN
811                  clvars(1) = 'PO2'
812               ENDIF
813            ENDIF
814
815            !Read in profile or profile obs types
816            CALL obs_rea_prof( profdata(jtype), ifilesprof(jtype),       &
817               &               clproffiles(jtype,1:ifilesprof(jtype)), &
818               &               nvarsprof(jtype), nextrprof(jtype), nitend-nit000+2, &
819               &               rn_dobsini, rn_dobsend, llvar, &
820               &               ln_ignmis, ln_s_at_t, .FALSE., ltype_clim, clvars, &
821               &               kdailyavtypes = nn_profdavtypes )
822
823            DO jvar = 1, nvarsprof(jtype)
824               CALL obs_prof_staend( profdata(jtype), jvar )
825            END DO
826
827            CALL obs_pre_prof( profdata(jtype), profdataqc(jtype), &
828               &               llvar, &
829               &               jpi, jpj, jpk, &
830               &               zmask, zglam, zgphi,  &
831               &               ln_nea, ln_bound_reject, &
832               &               kdailyavtypes = nn_profdavtypes )
833           
834            DEALLOCATE( llvar, clvars )
835            CALL wrk_dealloc( jpi, jpj,      nvarsprof(jtype), zglam )
836            CALL wrk_dealloc( jpi, jpj,      nvarsprof(jtype), zgphi )
837            CALL wrk_dealloc( jpi, jpj, jpk, nvarsprof(jtype), zmask )
838
839         END DO
840
841         DEALLOCATE( ifilesprof, clproffiles )
842
843      ENDIF
844
845      IF ( nsurftypes > 0 ) THEN
846
847         ALLOCATE(surfdata(nsurftypes))
848         ALLOCATE(surfdataqc(nsurftypes))
849         ALLOCATE(nvarssurf(nsurftypes))
850         ALLOCATE(nextrsurf(nsurftypes))
851
852         DO jtype = 1, nsurftypes
853
854            ltype_clim = .FALSE.
855            IF ( ln_output_clim .AND. &
856               & ( ( TRIM(cobstypessurf(jtype)) == 'sst' ) .OR. &
857               &   ( TRIM(cobstypessurf(jtype)) == 'sss' ) ) ) &
858               & ltype_clim = .TRUE.
859
860            IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN
861               nvarssurf(jtype) = 1
862               nextrsurf(jtype) = 2
863            ELSE
864               nvarssurf(jtype) = 1
865               nextrsurf(jtype) = 0
866            ENDIF
867           
868            ALLOCATE( clvars( nvarssurf(jtype) ) )
869
870            IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN
871               clvars(1) = 'SLA'
872            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sst' ) THEN
873               clvars(1) = 'SST'
874            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sic' ) THEN
875               clvars(1) = 'ICECONC'
876            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sss' ) THEN
877               clvars(1) = 'SSS'
878            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'slchltot' ) THEN
879               clvars(1) = 'SLCHLTOT'
880            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'slchldia' ) THEN
881               clvars(1) = 'SLCHLDIA'
882            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'slchlnon' ) THEN
883               clvars(1) = 'SLCHLNON'
884            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'slchldin' ) THEN
885               clvars(1) = 'SLCHLDIN'
886            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'slchlmic' ) THEN
887               clvars(1) = 'SLCHLMIC'
888            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'slchlnan' ) THEN
889               clvars(1) = 'SLCHLNAN'
890            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'slchlpic' ) THEN
891               clvars(1) = 'SLCHLPIC'
892            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'schltot' ) THEN
893               clvars(1) = 'SCHLTOT'
894            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'slphytot' ) THEN
895               clvars(1) = 'SLPHYTOT'
896            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'slphydia' ) THEN
897               clvars(1) = 'SLPHYDIA'
898            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'slphynon' ) THEN
899               clvars(1) = 'SLPHYNON'
900            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sspm' ) THEN
901               clvars(1) = 'SSPM'
902            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'skd490' ) THEN
903               clvars(1) = 'SKD490'
904            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sfco2' ) THEN
905               clvars(1) = 'SFCO2'
906            ELSE IF ( TRIM(cobstypessurf(jtype)) == 'spco2' ) THEN
907               clvars(1) = 'SPCO2'
908            ENDIF
909
910            !Read in surface obs types
911            CALL obs_rea_surf( surfdata(jtype), ifilessurf(jtype), &
912               &               clsurffiles(jtype,1:ifilessurf(jtype)), &
913               &               nvarssurf(jtype), nextrsurf(jtype), nitend-nit000+2, &
914               &               rn_dobsini, rn_dobsend, ln_ignmis, .FALSE., &
915               &               llnightav(jtype), ltype_clim, clvars )
916
917            CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea, ln_bound_reject )
918
919            IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN
920               CALL obs_rea_mdt( surfdataqc(jtype), n2dintsurf(jtype) )
921               IF ( ln_altbias ) &
922                  & CALL obs_rea_altbias ( surfdataqc(jtype), n2dintsurf(jtype), cn_altbiasfile )
923            ENDIF
924
925            IF ( TRIM(cobstypessurf(jtype)) == 'sst' .AND. ln_sstbias ) THEN
926               jnumsstbias = 0
927               DO jfile = 1, jpmaxnfiles
928                  IF ( TRIM(cn_sstbiasfiles(jfile)) /= '' ) &
929                     &  jnumsstbias = jnumsstbias + 1
930               END DO
931               IF ( jnumsstbias == 0 ) THEN
932                  CALL ctl_stop("ln_sstbias set but no bias files to read in")   
933               ENDIF
934
935               CALL obs_app_sstbias( surfdataqc(jtype), n2dintsurf(jtype), & 
936                  &                  jnumsstbias, cn_sstbiasfiles(1:jnumsstbias) ) 
937
938            ENDIF
939           
940            DEALLOCATE( clvars )
941
942         END DO
943
944         DEALLOCATE( ifilessurf, clsurffiles )
945
946      ENDIF
947
948   END SUBROUTINE dia_obs_init
949
950   SUBROUTINE dia_obs( kstp )
951      !!----------------------------------------------------------------------
952      !!                    ***  ROUTINE dia_obs  ***
953      !!         
954      !! ** Purpose : Call the observation operators on each time step
955      !!
956      !! ** Method  : Call the observation operators on each time step to
957      !!              compute the model equivalent of the following data:
958      !!               - Profile data, currently T/S or U/V
959      !!               - Surface data, currently SST, SLA or sea-ice concentration.
960      !!
961      !! ** Action  :
962      !!
963      !! History :
964      !!        !  06-03  (K. Mogensen) Original code
965      !!        !  06-05  (K. Mogensen) Reformatted
966      !!        !  06-10  (A. Weaver) Cleaning
967      !!        !  07-03  (K. Mogensen) General handling of profiles
968      !!        !  07-04  (G. Smith) Generalized surface operators
969      !!        !  08-10  (M. Valdivieso) obs operator for velocity profiles
970      !!        !  15-08  (M. Martin) Combined surface/profile routines.
971      !!----------------------------------------------------------------------
972      !! * Modules used
973      USE phycst, ONLY : &         ! Physical constants
974#if defined key_fabm
975         & rt0,          &
976#endif
977         & rday
978      USE oce, ONLY : &            ! Ocean dynamics and tracers variables
979         & tsn,       &
980         & un,        &
981         & vn,        &
982         & sshn
983#if defined  key_lim3
984      USE ice, ONLY : &            ! LIM3 Ice model variables
985         & frld
986#endif
987#if defined key_lim2
988      USE ice_2, ONLY : &          ! LIM2 Ice model variables
989         & frld
990#endif
991#if defined key_cice
992      USE sbc_oce, ONLY : fr_i     ! ice fraction
993#endif
994#if defined key_top
995      USE trc, ONLY :  &           ! Biogeochemical state variables
996         & trn
997#endif
998#if defined key_hadocc
999      USE par_hadocc               ! HadOCC parameters
1000      USE trc, ONLY :  &
1001         & HADOCC_CHL, &
1002         & HADOCC_FCO2, &
1003         & HADOCC_PCO2, &
1004         & HADOCC_FILL_FLT
1005      USE had_bgc_const, ONLY: c2n_p
1006#elif defined key_medusa
1007      USE par_medusa               ! MEDUSA parameters
1008      USE sms_medusa, ONLY: &
1009         & xthetapn, &
1010         & xthetapd
1011#if defined key_roam
1012      USE sms_medusa, ONLY: &
1013         & f2_pco2w, &
1014         & f2_fco2w, &
1015         & f3_pH
1016#endif
1017#elif defined key_fabm
1018      USE par_fabm                 ! FABM parameters
1019      USE fabm, ONLY: &
1020         & fabm_get_interior_diagnostic_data
1021#endif
1022#if defined key_spm
1023      USE par_spm, ONLY: &         ! Sediment parameters
1024         & jp_spm
1025#endif
1026
1027      IMPLICIT NONE
1028
1029      !! * Arguments
1030      INTEGER, INTENT(IN) :: kstp  ! Current timestep
1031      !! * Local declarations
1032      INTEGER :: idaystp           ! Number of timesteps per day
1033      INTEGER :: jtype             ! Data loop variable
1034      INTEGER :: jvar              ! Variable number
1035      INTEGER :: ji, jj, jk        ! Loop counters
1036      REAL(wp) :: tiny             ! small number
1037      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: &
1038         & zprofvar, &             ! Model values for variables in a prof ob
1039         & zprofclim               ! Climatology values for variables in a prof ob
1040      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: &
1041         & zprofmask               ! Mask associated with zprofvar
1042      REAL(wp), POINTER, DIMENSION(:,:) :: &
1043         & zsurfvar, &             ! Model values equivalent to surface ob.
1044         & zsurfclim, &            ! Climatology values for variables in a surface ob.
1045         & zsurfmask               ! Mask associated with surface variable
1046      REAL(wp), POINTER, DIMENSION(:,:,:) :: &
1047         & zglam,    &             ! Model longitudes for prof variables
1048         & zgphi                   ! Model latitudes for prof variables
1049      LOGICAL :: llog10            ! Perform log10 transform of variable
1050#if defined key_fabm
1051      REAL(wp), POINTER, DIMENSION(:,:,:) :: &
1052         & fabm_3d                 ! 3D variable from FABM
1053#endif
1054     
1055      IF(lwp) THEN
1056         WRITE(numout,*)
1057         WRITE(numout,*) 'dia_obs : Call the observation operators', kstp
1058         WRITE(numout,*) '~~~~~~~'
1059         CALL FLUSH(numout)
1060      ENDIF
1061
1062      idaystp = NINT( rday / rdt )
1063
1064      !-----------------------------------------------------------------------
1065      ! Call the profile and surface observation operators
1066      !-----------------------------------------------------------------------
1067
1068      IF ( nproftypes > 0 ) THEN
1069
1070         DO jtype = 1, nproftypes
1071
1072            ! Allocate local work arrays
1073            CALL wrk_alloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofvar  )
1074            CALL wrk_alloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofmask )
1075            CALL wrk_alloc( jpi, jpj,      profdataqc(jtype)%nvar, zglam     )
1076            CALL wrk_alloc( jpi, jpj,      profdataqc(jtype)%nvar, zgphi     )
1077            CALL wrk_alloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofclim )   
1078                             
1079            ! Defaults which might change
1080            DO jvar = 1, profdataqc(jtype)%nvar
1081               zprofmask(:,:,:,jvar) = tmask(:,:,:)
1082               zglam(:,:,jvar)       = glamt(:,:)
1083               zgphi(:,:,jvar)       = gphit(:,:)
1084               zprofclim(:,:,:,jvar) = 0._wp
1085            END DO
1086
1087            SELECT CASE ( TRIM(cobstypesprof(jtype)) )
1088
1089            CASE('prof')
1090               zprofvar(:,:,:,1) = tsn(:,:,:,jp_tem)
1091               zprofvar(:,:,:,2) = tsn(:,:,:,jp_sal)
1092               IF ( ln_output_clim ) THEN         
1093                  zprofclim(:,:,:,1) = tclim(:,:,:)
1094                  zprofclim(:,:,:,2) = sclim(:,:,:)
1095               ENDIF
1096               
1097            CASE('vel')
1098               zprofvar(:,:,:,1) = un(:,:,:)
1099               zprofvar(:,:,:,2) = vn(:,:,:)
1100               zprofmask(:,:,:,1) = umask(:,:,:)
1101               zprofmask(:,:,:,2) = vmask(:,:,:)
1102               zglam(:,:,1) = glamu(:,:)
1103               zglam(:,:,2) = glamv(:,:)
1104               zgphi(:,:,1) = gphiu(:,:)
1105               zgphi(:,:,2) = gphiv(:,:)
1106
1107            CASE('plchltot')
1108#if defined key_hadocc
1109               ! Chlorophyll from HadOCC
1110               zprofvar(:,:,:,1) = HADOCC_CHL(:,:,:)
1111#elif defined key_medusa
1112               ! Add non-diatom and diatom chlorophyll from MEDUSA
1113               zprofvar(:,:,:,1) = trn(:,:,:,jpchn) + trn(:,:,:,jpchd)
1114#elif defined key_fabm
1115               ! Add all chlorophyll groups from ERSEM
1116               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl2) + &
1117                  &                trn(:,:,:,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl4)
1118#else
1119               CALL ctl_stop( ' Trying to run plchltot observation operator', &
1120                  &           ' but no biogeochemical model appears to have been defined' )
1121#endif
1122               ! Take the log10 where we can, otherwise exclude
1123               tiny = 1.0e-20
1124               WHERE(zprofvar(:,:,:,:) > tiny .AND. zprofvar(:,:,:,:) /= obfillflt )
1125                  zprofvar(:,:,:,:)  = LOG10(zprofvar(:,:,:,:))
1126               ELSEWHERE
1127                  zprofvar(:,:,:,:)  = obfillflt
1128                  zprofmask(:,:,:,:) = 0
1129               END WHERE
1130               ! Mask out model below any excluded values,
1131               ! to avoid interpolation issues
1132               DO jvar = 1, profdataqc(jtype)%nvar
1133                 DO jj = 1, jpj
1134                    DO ji = 1, jpi
1135                       depth_loop: DO jk = 1, jpk
1136                          IF ( zprofmask(ji,jj,jk,jvar) == 0 ) THEN
1137                             zprofmask(ji,jj,jk:jpk,jvar) = 0
1138                             EXIT depth_loop
1139                          ENDIF
1140                       END DO depth_loop
1141                    END DO
1142                 END DO
1143              END DO
1144
1145            CASE('pchltot')
1146#if defined key_hadocc
1147               ! Chlorophyll from HadOCC
1148               zprofvar(:,:,:,1) = HADOCC_CHL(:,:,:)
1149#elif defined key_medusa
1150               ! Add non-diatom and diatom chlorophyll from MEDUSA
1151               zprofvar(:,:,:,1) = trn(:,:,:,jpchn) + trn(:,:,:,jpchd)
1152#elif defined key_fabm
1153               ! Add all chlorophyll groups from ERSEM
1154               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl2) + &
1155                  &                trn(:,:,:,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,:,jp_fabm_m1+jp_fabm_chl4)
1156#else
1157               CALL ctl_stop( ' Trying to run pchltot observation operator', &
1158                  &           ' but no biogeochemical model appears to have been defined' )
1159#endif
1160
1161            CASE('pno3')
1162#if defined key_hadocc
1163               ! Dissolved inorganic nitrogen from HadOCC
1164               zprofvar(:,:,:,1) = trn(:,:,:,jp_had_nut)
1165#elif defined key_medusa
1166               ! Dissolved inorganic nitrogen from MEDUSA
1167               zprofvar(:,:,:,1) = trn(:,:,:,jpdin)
1168#elif defined key_fabm
1169               ! Nitrate from ERSEM
1170               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_n3n)
1171#else
1172               CALL ctl_stop( ' Trying to run pno3 observation operator', &
1173                  &           ' but no biogeochemical model appears to have been defined' )
1174#endif
1175
1176            CASE('psi4')
1177#if defined key_hadocc
1178               CALL ctl_stop( ' Trying to run psi4 observation operator', &
1179                  &           ' but HadOCC does not simulate silicate' )
1180#elif defined key_medusa
1181               ! Silicate from MEDUSA
1182               zprofvar(:,:,:,1) = trn(:,:,:,jpsil)
1183#elif defined key_fabm
1184               ! Silicate from ERSEM
1185               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_n5s)
1186#else
1187               CALL ctl_stop( ' Trying to run psi4 observation operator', &
1188                  &           ' but no biogeochemical model appears to have been defined' )
1189#endif
1190
1191            CASE('ppo4')
1192#if defined key_hadocc
1193               CALL ctl_stop( ' Trying to run ppo4 observation operator', &
1194                  &           ' but HadOCC does not simulate phosphate' )
1195#elif defined key_medusa
1196               CALL ctl_stop( ' Trying to run ppo4 observation operator', &
1197                  &           ' but MEDUSA does not simulate phosphate' )
1198#elif defined key_fabm
1199               ! Phosphate from ERSEM
1200               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_n1p)
1201#else
1202               CALL ctl_stop( ' Trying to run ppo4 observation operator', &
1203                  &           ' but no biogeochemical model appears to have been defined' )
1204#endif
1205
1206            CASE('pdic')
1207#if defined key_hadocc
1208               ! Dissolved inorganic carbon from HadOCC
1209               zprofvar(:,:,:,1) = trn(:,:,:,jp_had_dic)
1210#elif defined key_medusa
1211               ! Dissolved inorganic carbon from MEDUSA
1212               zprofvar(:,:,:,1) = trn(:,:,:,jpdic)
1213#elif defined key_fabm
1214               ! Dissolved inorganic carbon from ERSEM
1215               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_o3c)
1216#else
1217               CALL ctl_stop( ' Trying to run pdic observation operator', &
1218                  &           ' but no biogeochemical model appears to have been defined' )
1219#endif
1220
1221            CASE('palk')
1222#if defined key_hadocc
1223               ! Alkalinity from HadOCC
1224               zprofvar(:,:,:,1) = trn(:,:,:,jp_had_alk)
1225#elif defined key_medusa
1226               ! Alkalinity from MEDUSA
1227               zprofvar(:,:,:,1) = trn(:,:,:,jpalk)
1228#elif defined key_fabm
1229               ! Alkalinity from ERSEM
1230               zprofvar(:,:,:,1) = fabm_get_interior_diagnostic_data(model, jp_fabm_o3ta)
1231#else
1232               CALL ctl_stop( ' Trying to run palk observation operator', &
1233                  &           ' but no biogeochemical model appears to have been defined' )
1234#endif
1235
1236            CASE('pph')
1237#if defined key_hadocc
1238               CALL ctl_stop( ' Trying to run pph observation operator', &
1239                  &           ' but HadOCC has no pH diagnostic defined' )
1240#elif defined key_medusa && defined key_roam
1241               ! pH from MEDUSA
1242               zprofvar(:,:,:,1) = f3_pH(:,:,:)
1243#elif defined key_fabm
1244               ! pH from ERSEM
1245               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_o3ph)
1246#else
1247               CALL ctl_stop( ' Trying to run pph observation operator', &
1248                  &           ' but no biogeochemical model appears to have been defined' )
1249#endif
1250
1251            CASE('po2')
1252#if defined key_hadocc
1253               CALL ctl_stop( ' Trying to run po2 observation operator', &
1254                  &           ' but HadOCC does not simulate oxygen' )
1255#elif defined key_medusa
1256               ! Oxygen from MEDUSA
1257               zprofvar(:,:,:,1) = trn(:,:,:,jpoxy)
1258#elif defined key_fabm
1259               ! Oxygen from ERSEM
1260               zprofvar(:,:,:,1) = trn(:,:,:,jp_fabm_m1+jp_fabm_o2o)
1261#else
1262               CALL ctl_stop( ' Trying to run po2 observation operator', &
1263                  &           ' but no biogeochemical model appears to have been defined' )
1264#endif
1265
1266            CASE DEFAULT
1267               CALL ctl_stop( 'Unknown profile observation type '//TRIM(cobstypesprof(jtype))//' in dia_obs' )
1268
1269            END SELECT
1270
1271            DO jvar = 1, profdataqc(jtype)%nvar
1272               CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk,  &
1273                  &               nit000, idaystp, jvar,                   &
1274                  &               zprofvar(:,:,:,jvar),                    &
1275                  &               zprofclim(:,:,:,jvar),                   &
1276                  &               fsdept(:,:,:), fsdepw(:,:,:),            & 
1277                  &               zprofmask(:,:,:,jvar),                   &
1278                  &               zglam(:,:,jvar), zgphi(:,:,jvar),        &
1279                  &               nn_1dint, nn_2dint_default,              &
1280                  &               kdailyavtypes = nn_profdavtypes )
1281            END DO
1282
1283            CALL wrk_dealloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofvar  )
1284            CALL wrk_dealloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofmask )
1285            CALL wrk_dealloc( jpi, jpj,      profdataqc(jtype)%nvar, zglam     )
1286            CALL wrk_dealloc( jpi, jpj,      profdataqc(jtype)%nvar, zgphi     )
1287            CALL wrk_dealloc( jpi, jpj, jpk, profdataqc(jtype)%nvar, zprofclim  )           
1288
1289         END DO
1290
1291      ENDIF
1292
1293      IF ( nsurftypes > 0 ) THEN
1294
1295         !Allocate local work arrays
1296         CALL wrk_alloc( jpi, jpj, zsurfvar )
1297         CALL wrk_alloc( jpi, jpj, zsurfclim )         
1298         CALL wrk_alloc( jpi, jpj, zsurfmask )
1299#if defined key_fabm
1300         CALL wrk_alloc( jpi, jpj, jpk, fabm_3d )
1301#endif
1302
1303         DO jtype = 1, nsurftypes
1304
1305            !Defaults which might be changed
1306            zsurfmask(:,:) = tmask(:,:,1)
1307            zsurfclim(:,:) = 0._wp         
1308            llog10 = .FALSE.
1309
1310            SELECT CASE ( TRIM(cobstypessurf(jtype)) )
1311            CASE('sst')
1312               zsurfvar(:,:) = tsn(:,:,1,jp_tem)
1313               IF ( ln_output_clim ) zsurfclim(:,:) = tclim(:,:,1)
1314            CASE('sla')
1315               zsurfvar(:,:) = sshn(:,:)
1316            CASE('sss')
1317               zsurfvar(:,:) = tsn(:,:,1,jp_sal)
1318               IF ( ln_output_clim ) zsurfclim(:,:) = sclim(:,:,1)             
1319            CASE('sic')
1320               IF ( kstp == 0 ) THEN
1321                  IF ( lwp .AND. surfdataqc(jtype)%nsstpmpp(1) > 0 ) THEN
1322                     CALL ctl_warn( 'Sea-ice not initialised on zeroth '// &
1323                        &           'time-step but some obs are valid then.' )
1324                     WRITE(numout,*)surfdataqc(jtype)%nsstpmpp(1), &
1325                        &           ' sea-ice obs will be missed'
1326                  ENDIF
1327                  surfdataqc(jtype)%nsurfup = surfdataqc(jtype)%nsurfup + &
1328                     &                        surfdataqc(jtype)%nsstp(1)
1329                  CYCLE
1330               ELSE
1331#if defined key_cice
1332                  zsurfvar(:,:) = fr_i(:,:)
1333#elif defined key_lim2 || defined key_lim3
1334                  zsurfvar(:,:) = 1._wp - frld(:,:)
1335#else
1336               CALL ctl_stop( ' Trying to run sea-ice observation operator', &
1337                  &           ' but no sea-ice model appears to have been defined' )
1338#endif
1339               ENDIF
1340
1341            CASE('slchltot')
1342#if defined key_hadocc
1343               ! Surface chlorophyll from HadOCC
1344               zsurfvar(:,:) = HADOCC_CHL(:,:,1)
1345#elif defined key_medusa
1346               ! Add non-diatom and diatom surface chlorophyll from MEDUSA
1347               zsurfvar(:,:) = trn(:,:,1,jpchn) + trn(:,:,1,jpchd)
1348#elif defined key_fabm
1349               ! Add all surface chlorophyll groups from ERSEM
1350               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) + &
1351                  &            trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4)
1352#else
1353               CALL ctl_stop( ' Trying to run slchltot observation operator', &
1354                  &           ' but no biogeochemical model appears to have been defined' )
1355#endif
1356               llog10 = .TRUE.
1357
1358            CASE('slchldia')
1359#if defined key_hadocc
1360               CALL ctl_stop( ' Trying to run slchldia observation operator', &
1361                  &           ' but HadOCC does not explicitly simulate diatoms' )
1362#elif defined key_medusa
1363               ! Diatom surface chlorophyll from MEDUSA
1364               zsurfvar(:,:) = trn(:,:,1,jpchd)
1365#elif defined key_fabm
1366               ! Diatom surface chlorophyll from ERSEM
1367               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1)
1368#else
1369               CALL ctl_stop( ' Trying to run slchldia observation operator', &
1370                  &           ' but no biogeochemical model appears to have been defined' )
1371#endif
1372               llog10 = .TRUE.
1373
1374            CASE('slchlnon')
1375#if defined key_hadocc
1376               CALL ctl_stop( ' Trying to run slchlnon observation operator', &
1377                  &           ' but HadOCC does not explicitly simulate non-diatoms' )
1378#elif defined key_medusa
1379               ! Non-diatom surface chlorophyll from MEDUSA
1380               zsurfvar(:,:) = trn(:,:,1,jpchn)
1381#elif defined key_fabm
1382               ! Add all non-diatom surface chlorophyll groups from ERSEM
1383               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) + &
1384                  &            trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4)
1385#else
1386               CALL ctl_stop( ' Trying to run slchlnon observation operator', &
1387                  &           ' but no biogeochemical model appears to have been defined' )
1388#endif
1389               llog10 = .TRUE.
1390
1391            CASE('slchldin')
1392#if defined key_hadocc
1393               CALL ctl_stop( ' Trying to run slchldin observation operator', &
1394                  &           ' but HadOCC does not explicitly simulate dinoflagellates' )
1395#elif defined key_medusa
1396               CALL ctl_stop( ' Trying to run slchldin observation operator', &
1397                  &           ' but MEDUSA does not explicitly simulate dinoflagellates' )
1398#elif defined key_fabm
1399               ! Dinoflagellate surface chlorophyll from ERSEM
1400               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl4)
1401#else
1402               CALL ctl_stop( ' Trying to run slchldin observation operator', &
1403                  &           ' but no biogeochemical model appears to have been defined' )
1404#endif
1405               llog10 = .TRUE.
1406
1407            CASE('slchlmic')
1408#if defined key_hadocc
1409               CALL ctl_stop( ' Trying to run slchlmic observation operator', &
1410                  &           ' but HadOCC does not explicitly simulate microphytoplankton' )
1411#elif defined key_medusa
1412               CALL ctl_stop( ' Trying to run slchlmic observation operator', &
1413                  &           ' but MEDUSA does not explicitly simulate microphytoplankton' )
1414#elif defined key_fabm
1415               ! Add diatom and dinoflagellate surface chlorophyll from ERSEM
1416               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4)
1417#else
1418               CALL ctl_stop( ' Trying to run slchlmic observation operator', &
1419                  &           ' but no biogeochemical model appears to have been defined' )
1420#endif
1421               llog10 = .TRUE.
1422
1423            CASE('slchlnan')
1424#if defined key_hadocc
1425               CALL ctl_stop( ' Trying to run slchlnan observation operator', &
1426                  &           ' but HadOCC does not explicitly simulate nanophytoplankton' )
1427#elif defined key_medusa
1428               CALL ctl_stop( ' Trying to run slchlnan observation operator', &
1429                  &           ' but MEDUSA does not explicitly simulate nanophytoplankton' )
1430#elif defined key_fabm
1431               ! Nanophytoplankton surface chlorophyll from ERSEM
1432               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl2)
1433#else
1434               CALL ctl_stop( ' Trying to run slchlnan observation operator', &
1435                  &           ' but no biogeochemical model appears to have been defined' )
1436#endif
1437               llog10 = .TRUE.
1438
1439            CASE('slchlpic')
1440#if defined key_hadocc
1441               CALL ctl_stop( ' Trying to run slchlpic observation operator', &
1442                  &           ' but HadOCC does not explicitly simulate picophytoplankton' )
1443#elif defined key_medusa
1444               CALL ctl_stop( ' Trying to run slchlpic observation operator', &
1445                  &           ' but MEDUSA does not explicitly simulate picophytoplankton' )
1446#elif defined key_fabm
1447               ! Picophytoplankton surface chlorophyll from ERSEM
1448               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl3)
1449#else
1450               CALL ctl_stop( ' Trying to run slchlpic observation operator', &
1451                  &           ' but no biogeochemical model appears to have been defined' )
1452#endif
1453               llog10 = .TRUE.
1454
1455            CASE('schltot')
1456#if defined key_hadocc
1457               ! Surface chlorophyll from HadOCC
1458               zsurfvar(:,:) = HADOCC_CHL(:,:,1)
1459#elif defined key_medusa
1460               ! Add non-diatom and diatom surface chlorophyll from MEDUSA
1461               zsurfvar(:,:) = trn(:,:,1,jpchn) + trn(:,:,1,jpchd)
1462#elif defined key_fabm
1463               ! Add all surface chlorophyll groups from ERSEM
1464               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_chl1) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl2) + &
1465                  &            trn(:,:,1,jp_fabm_m1+jp_fabm_chl3) + trn(:,:,1,jp_fabm_m1+jp_fabm_chl4)
1466#else
1467               CALL ctl_stop( ' Trying to run schltot observation operator', &
1468                  &           ' but no biogeochemical model appears to have been defined' )
1469#endif
1470
1471            CASE('slphytot')
1472#if defined key_hadocc
1473               ! Surface phytoplankton nitrogen from HadOCC multiplied by C:N ratio
1474               zsurfvar(:,:) = trn(:,:,1,jp_had_phy) * c2n_p
1475#elif defined key_medusa
1476               ! Add non-diatom and diatom surface phytoplankton nitrogen from MEDUSA
1477               ! multiplied by C:N ratio for each
1478               zsurfvar(:,:) = (trn(:,:,1,jpphn) * xthetapn) + (trn(:,:,1,jpphd) * xthetapd)
1479#elif defined key_fabm
1480               ! Add all surface phytoplankton carbon groups from ERSEM
1481               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_p1c) + trn(:,:,1,jp_fabm_m1+jp_fabm_p2c) + &
1482                  &            trn(:,:,1,jp_fabm_m1+jp_fabm_p3c) + trn(:,:,1,jp_fabm_m1+jp_fabm_p4c)
1483#else
1484               CALL ctl_stop( ' Trying to run slphytot observation operator', &
1485                  &           ' but no biogeochemical model appears to have been defined' )
1486#endif
1487               llog10 = .TRUE.
1488
1489            CASE('slphydia')
1490#if defined key_hadocc
1491               CALL ctl_stop( ' Trying to run slphydia observation operator', &
1492                  &           ' but HadOCC does not explicitly simulate diatoms' )
1493#elif defined key_medusa
1494               ! Diatom surface phytoplankton nitrogen from MEDUSA multiplied by C:N ratio
1495               zsurfvar(:,:) = trn(:,:,1,jpphd) * xthetapd
1496#elif defined key_fabm
1497               ! Diatom surface phytoplankton carbon from ERSEM
1498               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_p1c)
1499#else
1500               CALL ctl_stop( ' Trying to run slphydia observation operator', &
1501                  &           ' but no biogeochemical model appears to have been defined' )
1502#endif
1503               llog10 = .TRUE.
1504
1505            CASE('slphynon')
1506#if defined key_hadocc
1507               CALL ctl_stop( ' Trying to run slphynon observation operator', &
1508                  &           ' but HadOCC does not explicitly simulate non-diatoms' )
1509#elif defined key_medusa
1510               ! Non-diatom surface phytoplankton nitrogen from MEDUSA multiplied by C:N ratio
1511               zsurfvar(:,:) = trn(:,:,1,jpphn) * xthetapn
1512#elif defined key_fabm
1513               ! Add all non-diatom surface phytoplankton carbon groups from ERSEM
1514               zsurfvar(:,:) = trn(:,:,1,jp_fabm_m1+jp_fabm_p2c) + &
1515                  &            trn(:,:,1,jp_fabm_m1+jp_fabm_p3c) + trn(:,:,1,jp_fabm_m1+jp_fabm_p4c)
1516#else
1517               CALL ctl_stop( ' Trying to run slphynon observation operator', &
1518                  &           ' but no biogeochemical model appears to have been defined' )
1519#endif
1520               llog10 = .TRUE.
1521
1522            CASE('sspm')
1523#if defined key_spm
1524               zsurfvar(:,:) = 0.0
1525               DO jn = 1, jp_spm
1526                  zsurfvar(:,:) = zsurfvar(:,:) + trn(:,:,1,jn)   ! sum SPM sizes
1527               END DO
1528#else
1529               CALL ctl_stop( ' Trying to run sspm observation operator', &
1530                  &           ' but no spm model appears to have been defined' )
1531#endif
1532
1533            CASE('skd490')
1534#if defined key_hadocc
1535               CALL ctl_stop( ' Trying to run skd490 observation operator', &
1536                  &           ' but HadOCC does not explicitly simulate Kd490' )
1537#elif defined key_medusa
1538               CALL ctl_stop( ' Trying to run skd490 observation operator', &
1539                  &           ' but MEDUSA does not explicitly simulate Kd490' )
1540#elif defined key_fabm
1541               ! light_xEPS diagnostic variable
1542               fabm_3d(:,:,:) = fabm_get_interior_diagnostic_data(model, jp_fabm_xeps)
1543               zsurfvar(:,:) = fabm_3d(:,:,1)
1544#else
1545               CALL ctl_stop( ' Trying to run skd490 observation operator', &
1546                  &           ' but no biogeochemical model appears to have been defined' )
1547#endif
1548
1549            CASE('sfco2')
1550#if defined key_hadocc
1551               zsurfvar(:,:) = HADOCC_FCO2(:,:)    ! fCO2 from HadOCC
1552               IF ( ( MINVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) .AND. &
1553                  & ( MAXVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) ) THEN
1554                  zsurfvar(:,:) = obfillflt
1555                  zsurfmask(:,:) = 0
1556                  CALL ctl_warn( ' HadOCC fCO2 values masked out for observation operator', &
1557                     &           ' as HADOCC_FCO2(:,:) == HADOCC_FILL_FLT' )
1558               ENDIF
1559#elif defined key_medusa && defined key_roam
1560               zsurfvar(:,:) = f2_fco2w(:,:)
1561#elif defined key_fabm
1562               ! First, get pCO2 from FABM
1563               fabm_3d(:,:,:) = fabm_get_interior_diagnostic_data(model, jp_fabm_o3pc)
1564               zsurfvar(:,:) = fabm_3d(:,:,1)
1565               ! Now, convert pCO2 to fCO2, based on SST in K. This follows the standard methodology of:
1566               ! Pierrot et al. (2009), Recommendations for autonomous underway pCO2 measuring systems
1567               ! and data reduction routines, Deep-Sea Research II, 56: 512-522.
1568               ! and
1569               ! Weiss (1974), Carbon dioxide in water and seawater: the solubility of a non-ideal gas,
1570               ! Marine Chemistry, 2: 203-215.
1571               ! In the implementation below, atmospheric pressure has been assumed to be 1 atm and so
1572               ! not explicitly included - atmospheric pressure is not necessarily available so this is
1573               ! the best assumption.
1574               ! Further, the (1-xCO2)^2 term has been neglected. This is common practice
1575               ! (see e.g. Zeebe and Wolf-Gladrow (2001), CO2 in Seawater: Equilibrium, Kinetics, Isotopes)
1576               ! because xCO2 in atm is ~0, and so this term will only affect the result to the 3rd decimal
1577               ! place for typical values, and xCO2 would need to be approximated from pCO2 anyway.
1578               zsurfvar(:,:) = zsurfvar(:,:) * EXP((-1636.75                                                          + &
1579                  &            12.0408      * (tsn(:,:,1,jp_tem)+rt0)                                                 - &
1580                  &            0.0327957    * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0)                         + &
1581                  &            0.0000316528 * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0) + &
1582                  &            2.0 * (57.7 - 0.118 * (tsn(:,:,1,jp_tem)+rt0)))                                        / &
1583                  &            (82.0578 * (tsn(:,:,1,jp_tem)+rt0)))
1584#else
1585               CALL ctl_stop( ' Trying to run sfco2 observation operator', &
1586                  &           ' but no biogeochemical model appears to have been defined' )
1587#endif
1588
1589            CASE('spco2')
1590#if defined key_hadocc
1591               zsurfvar(:,:) = HADOCC_PCO2(:,:)    ! pCO2 from HadOCC
1592               IF ( ( MINVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) .AND. &
1593                  & ( MAXVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) ) THEN
1594                  zsurfvar(:,:) = obfillflt
1595                  zsurfmask(:,:) = 0
1596                  CALL ctl_warn( ' HadOCC pCO2 values masked out for observation operator', &
1597                     &           ' as HADOCC_PCO2(:,:) == HADOCC_FILL_FLT' )
1598               ENDIF
1599#elif defined key_medusa && defined key_roam
1600               zsurfvar(:,:) = f2_pco2w(:,:)
1601#elif defined key_fabm
1602               fabm_3d(:,:,:) = fabm_get_interior_diagnostic_data(model, jp_fabm_o3pc)
1603               zsurfvar(:,:) = fabm_3d(:,:,1)
1604#else
1605               CALL ctl_stop( ' Trying to run spco2 observation operator', &
1606                  &           ' but no biogeochemical model appears to have been defined' )
1607#endif
1608
1609            CASE DEFAULT
1610
1611               CALL ctl_stop( 'Unknown surface observation type '//TRIM(cobstypessurf(jtype))//' in dia_obs' )
1612
1613            END SELECT
1614           
1615            IF ( llog10 ) THEN
1616               ! Take the log10 where we can, otherwise exclude
1617               tiny = 1.0e-20
1618               WHERE(zsurfvar(:,:) > tiny .AND. zsurfvar(:,:) /= obfillflt )
1619                  zsurfvar(:,:)  = LOG10(zsurfvar(:,:))
1620               ELSEWHERE
1621                  zsurfvar(:,:)  = obfillflt
1622                  zsurfmask(:,:) = 0
1623               END WHERE
1624            ENDIF
1625
1626            CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj,       &
1627               &               nit000, idaystp, zsurfvar,               &
1628               &               zsurfclim, zsurfmask,                    &
1629               &               n2dintsurf(jtype), llnightav(jtype),     &
1630               &               ravglamscl(jtype), ravgphiscl(jtype),    &
1631               &               lfpindegs(jtype) )
1632
1633         END DO
1634
1635         CALL wrk_dealloc( jpi, jpj, zsurfvar )
1636         CALL wrk_dealloc( jpi, jpj, zsurfmask )
1637#if defined key_fabm
1638         CALL wrk_dealloc( jpi, jpj, jpk, fabm_3d )
1639#endif
1640
1641      ENDIF
1642
1643   END SUBROUTINE dia_obs
1644
1645   SUBROUTINE dia_obs_wri
1646      !!----------------------------------------------------------------------
1647      !!                    ***  ROUTINE dia_obs_wri  ***
1648      !!         
1649      !! ** Purpose : Call observation diagnostic output routines
1650      !!
1651      !! ** Method  : Call observation diagnostic output routines
1652      !!
1653      !! ** Action  :
1654      !!
1655      !! History :
1656      !!        !  06-03  (K. Mogensen) Original code
1657      !!        !  06-05  (K. Mogensen) Reformatted
1658      !!        !  06-10  (A. Weaver) Cleaning
1659      !!        !  07-03  (K. Mogensen) General handling of profiles
1660      !!        !  08-09  (M. Valdivieso) Velocity component (U,V) profiles
1661      !!        !  15-08  (M. Martin) Combined writing for prof and surf types
1662      !!----------------------------------------------------------------------
1663      !! * Modules used
1664      USE obs_rot_vel          ! Rotation of velocities
1665
1666      IMPLICIT NONE
1667
1668      !! * Local declarations
1669      INTEGER :: jtype                    ! Data set loop variable
1670      INTEGER :: jo, jvar, jk
1671      REAL(wp), DIMENSION(:), ALLOCATABLE :: &
1672         & zu, &
1673         & zv
1674
1675      !-----------------------------------------------------------------------
1676      ! Depending on switches call various observation output routines
1677      !-----------------------------------------------------------------------
1678
1679      IF ( nproftypes > 0 ) THEN
1680
1681         DO jtype = 1, nproftypes
1682
1683            IF ( TRIM(cobstypesprof(jtype)) == 'vel' ) THEN
1684
1685               ! For velocity data, rotate the model velocities to N/S, E/W
1686               ! using the compressed data structure.
1687               ALLOCATE( &
1688                  & zu(profdataqc(jtype)%nvprot(1)), &
1689                  & zv(profdataqc(jtype)%nvprot(2))  &
1690                  & )
1691
1692               CALL obs_rotvel( profdataqc(jtype), nn_2dint_default, zu, zv )
1693
1694               DO jo = 1, profdataqc(jtype)%nprof
1695                  DO jvar = 1, 2
1696                     DO jk = profdataqc(jtype)%npvsta(jo,jvar), profdataqc(jtype)%npvend(jo,jvar)
1697
1698                        IF ( jvar == 1 ) THEN
1699                           profdataqc(jtype)%var(jvar)%vmod(jk) = zu(jk)
1700                        ELSE
1701                           profdataqc(jtype)%var(jvar)%vmod(jk) = zv(jk)
1702                        ENDIF
1703
1704                     END DO
1705                  END DO
1706               END DO
1707
1708               DEALLOCATE( zu )
1709               DEALLOCATE( zv )
1710
1711            END IF
1712
1713            CALL obs_prof_decompress( profdataqc(jtype), &
1714               &                      profdata(jtype), .TRUE., numout )
1715
1716            CALL obs_wri_prof( profdata(jtype) )
1717
1718         END DO
1719
1720      ENDIF
1721
1722      IF ( nsurftypes > 0 ) THEN
1723
1724         DO jtype = 1, nsurftypes
1725
1726            CALL obs_surf_decompress( surfdataqc(jtype), &
1727               &                      surfdata(jtype), .TRUE., numout )
1728
1729            CALL obs_wri_surf( surfdata(jtype) )
1730
1731         END DO
1732
1733      ENDIF
1734
1735   END SUBROUTINE dia_obs_wri
1736
1737   SUBROUTINE dia_obs_dealloc
1738      IMPLICIT NONE
1739      !!----------------------------------------------------------------------
1740      !!                    *** ROUTINE dia_obs_dealloc ***
1741      !!
1742      !!  ** Purpose : To deallocate data to enable the obs_oper online loop.
1743      !!               Specifically: dia_obs_init --> dia_obs --> dia_obs_wri
1744      !!
1745      !!  ** Method : Clean up various arrays left behind by the obs_oper.
1746      !!
1747      !!  ** Action :
1748      !!
1749      !!----------------------------------------------------------------------
1750      ! obs_grid deallocation
1751      CALL obs_grid_deallocate
1752
1753      ! diaobs deallocation
1754      IF ( nproftypes > 0 ) &
1755         &   DEALLOCATE( cobstypesprof, profdata, profdataqc, nvarsprof, nextrprof )
1756
1757      IF ( nsurftypes > 0 ) &
1758         &   DEALLOCATE( cobstypessurf, surfdata, surfdataqc, nvarssurf, nextrsurf, &
1759         &               n2dintsurf, ravglamscl, ravgphiscl, lfpindegs, llnightav )
1760
1761   END SUBROUTINE dia_obs_dealloc
1762
1763   SUBROUTINE ini_date( ddobsini )
1764      !!----------------------------------------------------------------------
1765      !!                    ***  ROUTINE ini_date  ***
1766      !!
1767      !! ** Purpose : Get initial date in double precision YYYYMMDD.HHMMSS format
1768      !!
1769      !! ** Method  : Get initial date in double precision YYYYMMDD.HHMMSS format
1770      !!
1771      !! ** Action  : Get initial date in double precision YYYYMMDD.HHMMSS format
1772      !!
1773      !! History :
1774      !!        !  06-03  (K. Mogensen)  Original code
1775      !!        !  06-05  (K. Mogensen)  Reformatted
1776      !!        !  06-10  (A. Weaver) Cleaning
1777      !!        !  06-10  (G. Smith) Calculates initial date the same as method for final date
1778      !!        !  10-05  (D. Lea) Update to month length calculation for NEMO vn3.2
1779      !!----------------------------------------------------------------------
1780      USE phycst, ONLY : &            ! Physical constants
1781         & rday
1782      USE dom_oce, ONLY : &           ! Ocean space and time domain variables
1783         & rdt
1784
1785      IMPLICIT NONE
1786
1787      !! * Arguments
1788      REAL(dp), INTENT(OUT) :: ddobsini  ! Initial date in YYYYMMDD.HHMMSS
1789
1790      !! * Local declarations
1791      INTEGER :: iyea        ! date - (year, month, day, hour, minute)
1792      INTEGER :: imon
1793      INTEGER :: iday
1794      INTEGER :: ihou
1795      INTEGER :: imin
1796      INTEGER :: imday       ! Number of days in month.
1797      INTEGER, DIMENSION(12) :: &
1798         &       imonth_len  ! Length in days of the months of the current year
1799      REAL(wp) :: zdayfrc    ! Fraction of day
1800
1801      !----------------------------------------------------------------------
1802      ! Initial date initialization (year, month, day, hour, minute)
1803      ! (This assumes that the initial date is for 00z))
1804      !----------------------------------------------------------------------
1805      iyea =   ndate0 / 10000
1806      imon = ( ndate0 - iyea * 10000 ) / 100
1807      iday =   ndate0 - iyea * 10000 - imon * 100
1808      ihou = 0
1809      imin = 0
1810
1811      !----------------------------------------------------------------------
1812      ! Compute number of days + number of hours + min since initial time
1813      !----------------------------------------------------------------------
1814      iday = iday + ( nit000 -1 ) * rdt / rday
1815      zdayfrc = ( nit000 -1 ) * rdt / rday
1816      zdayfrc = zdayfrc - aint(zdayfrc)
1817      ihou = int( zdayfrc * 24 )
1818      imin = int( (zdayfrc * 24 - ihou) * 60 )
1819
1820      !-----------------------------------------------------------------------
1821      ! Convert number of days (iday) into a real date
1822      !----------------------------------------------------------------------
1823
1824      CALL calc_month_len( iyea, imonth_len )
1825
1826      DO WHILE ( iday > imonth_len(imon) )
1827         iday = iday - imonth_len(imon)
1828         imon = imon + 1 
1829         IF ( imon > 12 ) THEN
1830            imon = 1
1831            iyea = iyea + 1
1832            CALL calc_month_len( iyea, imonth_len )  ! update month lengths
1833         ENDIF
1834      END DO
1835
1836      !----------------------------------------------------------------------
1837      ! Convert it into YYYYMMDD.HHMMSS format.
1838      !----------------------------------------------------------------------
1839      ddobsini = iyea * 10000_dp + imon * 100_dp + &
1840         &       iday + ihou * 0.01_dp + imin * 0.0001_dp
1841
1842
1843   END SUBROUTINE ini_date
1844
1845   SUBROUTINE fin_date( ddobsfin )
1846      !!----------------------------------------------------------------------
1847      !!                    ***  ROUTINE fin_date  ***
1848      !!
1849      !! ** Purpose : Get final date in double precision YYYYMMDD.HHMMSS format
1850      !!
1851      !! ** Method  : Get final date in double precision YYYYMMDD.HHMMSS format
1852      !!
1853      !! ** Action  : Get final date in double precision YYYYMMDD.HHMMSS format
1854      !!
1855      !! History :
1856      !!        !  06-03  (K. Mogensen)  Original code
1857      !!        !  06-05  (K. Mogensen)  Reformatted
1858      !!        !  06-10  (A. Weaver) Cleaning
1859      !!        !  10-05  (D. Lea) Update to month length calculation for NEMO vn3.2
1860      !!----------------------------------------------------------------------
1861      USE phycst, ONLY : &            ! Physical constants
1862         & rday
1863      USE dom_oce, ONLY : &           ! Ocean space and time domain variables
1864         & rdt
1865
1866      IMPLICIT NONE
1867
1868      !! * Arguments
1869      REAL(dp), INTENT(OUT) :: ddobsfin ! Final date in YYYYMMDD.HHMMSS
1870
1871      !! * Local declarations
1872      INTEGER :: iyea        ! date - (year, month, day, hour, minute)
1873      INTEGER :: imon
1874      INTEGER :: iday
1875      INTEGER :: ihou
1876      INTEGER :: imin
1877      INTEGER :: imday       ! Number of days in month.
1878      INTEGER, DIMENSION(12) :: &
1879         &       imonth_len  ! Length in days of the months of the current year
1880      REAL(wp) :: zdayfrc    ! Fraction of day
1881
1882      !-----------------------------------------------------------------------
1883      ! Initial date initialization (year, month, day, hour, minute)
1884      ! (This assumes that the initial date is for 00z)
1885      !-----------------------------------------------------------------------
1886      iyea =   ndate0 / 10000
1887      imon = ( ndate0 - iyea * 10000 ) / 100
1888      iday =   ndate0 - iyea * 10000 - imon * 100
1889      ihou = 0
1890      imin = 0
1891
1892      !-----------------------------------------------------------------------
1893      ! Compute number of days + number of hours + min since initial time
1894      !-----------------------------------------------------------------------
1895      iday    = iday +  nitend  * rdt / rday
1896      zdayfrc =  nitend  * rdt / rday
1897      zdayfrc = zdayfrc - AINT( zdayfrc )
1898      ihou    = INT( zdayfrc * 24 )
1899      imin    = INT( ( zdayfrc * 24 - ihou ) * 60 )
1900
1901      !-----------------------------------------------------------------------
1902      ! Convert number of days (iday) into a real date
1903      !----------------------------------------------------------------------
1904
1905      CALL calc_month_len( iyea, imonth_len )
1906
1907      DO WHILE ( iday > imonth_len(imon) )
1908         iday = iday - imonth_len(imon)
1909         imon = imon + 1 
1910         IF ( imon > 12 ) THEN
1911            imon = 1
1912            iyea = iyea + 1
1913            CALL calc_month_len( iyea, imonth_len )  ! update month lengths
1914         ENDIF
1915      END DO
1916
1917      !-----------------------------------------------------------------------
1918      ! Convert it into YYYYMMDD.HHMMSS format
1919      !-----------------------------------------------------------------------
1920      ddobsfin = iyea * 10000_dp + imon * 100_dp    + iday &
1921         &     + ihou * 0.01_dp  + imin * 0.0001_dp
1922
1923    END SUBROUTINE fin_date
1924
1925    SUBROUTINE obs_settypefiles( ntypes, jpmaxnfiles, ifiles, cobstypes, cfiles )
1926
1927       INTEGER, INTENT(IN) :: ntypes      ! Total number of obs types
1928       INTEGER, INTENT(IN) :: jpmaxnfiles ! Maximum number of files allowed for each type
1929       INTEGER, DIMENSION(ntypes), INTENT(OUT) :: &
1930          &                   ifiles      ! Out number of files for each type
1931       CHARACTER(len=8), DIMENSION(ntypes), INTENT(IN) :: &
1932          &                   cobstypes   ! List of obs types
1933       CHARACTER(len=128), DIMENSION(ntypes, jpmaxnfiles), INTENT(IN) :: &
1934          &                   cfiles      ! List of files for all types
1935
1936       !Local variables
1937       INTEGER :: jfile
1938       INTEGER :: jtype
1939
1940       DO jtype = 1, ntypes
1941
1942          ifiles(jtype) = 0
1943          DO jfile = 1, jpmaxnfiles
1944             IF ( trim(cfiles(jtype,jfile)) /= '' ) &
1945                       ifiles(jtype) = ifiles(jtype) + 1
1946          END DO
1947
1948          IF ( ifiles(jtype) == 0 ) THEN
1949               CALL ctl_stop( 'Logical for observation type '//TRIM(cobstypes(jtype))//   &
1950                  &           ' set to true but no files available to read' )
1951          ENDIF
1952
1953          IF(lwp) THEN   
1954             WRITE(numout,*) '             '//cobstypes(jtype)//' input observation file names:'
1955             DO jfile = 1, ifiles(jtype)
1956                WRITE(numout,*) '                '//TRIM(cfiles(jtype,jfile))
1957             END DO
1958          ENDIF
1959
1960       END DO
1961
1962    END SUBROUTINE obs_settypefiles
1963
1964    SUBROUTINE obs_setinterpopts( ntypes, jtype, ctypein,             &
1965               &                  n2dint_default, n2dint_type,        &
1966               &                  ravglamscl_type, ravgphiscl_type,   &
1967               &                  lfp_indegs_type, lavnight_type,     &
1968               &                  n2dint, ravglamscl, ravgphiscl,     &
1969               &                  lfpindegs, lavnight )
1970
1971       INTEGER, INTENT(IN)  :: ntypes             ! Total number of obs types
1972       INTEGER, INTENT(IN)  :: jtype              ! Index of the current type of obs
1973       INTEGER, INTENT(IN)  :: n2dint_default     ! Default option for interpolation type
1974       INTEGER, INTENT(IN)  :: n2dint_type        ! Option for interpolation type
1975       REAL(wp), INTENT(IN) :: &
1976          &                    ravglamscl_type, & !E/W diameter of obs footprint for this type
1977          &                    ravgphiscl_type    !N/S diameter of obs footprint for this type
1978       LOGICAL, INTENT(IN)  :: lfp_indegs_type    !T=> footprint in degrees, F=> in metres
1979       LOGICAL, INTENT(IN)  :: lavnight_type      !T=> obs represent night time average
1980       CHARACTER(len=8), INTENT(IN) :: ctypein 
1981
1982       INTEGER, DIMENSION(ntypes), INTENT(INOUT) :: &
1983          &                    n2dint 
1984       REAL(wp), DIMENSION(ntypes), INTENT(INOUT) :: &
1985          &                    ravglamscl, ravgphiscl
1986       LOGICAL, DIMENSION(ntypes), INTENT(INOUT) :: &
1987          &                    lfpindegs, lavnight
1988
1989       lavnight(jtype) = lavnight_type
1990
1991       IF ( (n2dint_type >= 0) .AND. (n2dint_type <= 6) ) THEN
1992          n2dint(jtype) = n2dint_type
1993       ELSE IF ( n2dint_type == -1 ) THEN
1994          n2dint(jtype) = n2dint_default
1995       ELSE
1996          CALL ctl_stop(' Choice of '//TRIM(ctypein)//' horizontal (2D) interpolation method', &
1997            &                    ' is not available')
1998       ENDIF
1999
2000       ! For averaging observation footprints set options for size of footprint
2001       IF ( (n2dint(jtype) > 4) .AND. (n2dint(jtype) <= 6) ) THEN
2002          IF ( ravglamscl_type > 0._wp ) THEN
2003             ravglamscl(jtype) = ravglamscl_type
2004          ELSE
2005             CALL ctl_stop( 'Incorrect value set for averaging footprint '// &
2006                            'scale (ravglamscl) for observation type '//TRIM(ctypein) )     
2007          ENDIF
2008
2009          IF ( ravgphiscl_type > 0._wp ) THEN
2010             ravgphiscl(jtype) = ravgphiscl_type
2011          ELSE
2012             CALL ctl_stop( 'Incorrect value set for averaging footprint '// &
2013                            'scale (ravgphiscl) for observation type '//TRIM(ctypein) )     
2014          ENDIF
2015
2016          lfpindegs(jtype) = lfp_indegs_type 
2017
2018       ENDIF
2019
2020       ! Write out info
2021       IF(lwp) THEN
2022          IF ( n2dint(jtype) <= 4 ) THEN
2023             WRITE(numout,*) '             '//TRIM(ctypein)// &
2024                &            ' model counterparts will be interpolated horizontally'
2025          ELSE IF ( n2dint(jtype) <= 6 ) THEN
2026             WRITE(numout,*) '             '//TRIM(ctypein)// &
2027                &            ' model counterparts will be averaged horizontally'
2028             WRITE(numout,*) '             '//'    with E/W scale: ',ravglamscl(jtype)
2029             WRITE(numout,*) '             '//'    with N/S scale: ',ravgphiscl(jtype)
2030             IF ( lfpindegs(jtype) ) THEN
2031                 WRITE(numout,*) '             '//'    (in degrees)'
2032             ELSE
2033                 WRITE(numout,*) '             '//'    (in metres)'
2034             ENDIF
2035          ENDIF
2036       ENDIF
2037
2038    END SUBROUTINE obs_setinterpopts
2039
2040END MODULE diaobs
Note: See TracBrowser for help on using the repository browser.