source: branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 11546

Last change on this file since 11546 was 11546, checked in by dford, 23 months ago

Add an observation operator for Kd490, and add some error checking so that the observation variable name read in from a feedback file matches a defined name for each observation type. Previously it could read in an incorrect variable without failure. See https://code.metoffice.gov.uk/trac/utils/ticket/247.

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.