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

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

source: branches/UKMO/dev_r5518_obs_oper_update_utils366_fabmv1_v2/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 13487

Last change on this file since 13487 was 13487, checked in by dford, 4 years ago

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