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

Last change on this file since 13504 was 13504, checked in by dford, 5 months ago

Update for getting kd490 from spectral optical model.

File size: 105.7 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_Kd_band3 diagnostic variable if using spectral optical model
1685               ! light_xEPS diagnostic variable if using standard ERSEM light model
1686               IF ( jp_fabm_kd490 /= -1 ) THEN
1687                  fabm_3d(:,:,:) = model%get_interior_diagnostic_data(jp_fabm_kd490)
1688               ELSEIF ( jp_fabm_xeps /= -1 ) THEN
1689                  fabm_3d(:,:,:) = model%get_interior_diagnostic_data(jp_fabm_xeps)
1690               ELSE
1691                  CALL ctl_stop( ' Trying to run skd490 observation operator', &
1692                     &           ' but cannot access Kd490 from ERSEM' )
1693               ENDIF
1694               zsurfvar(:,:,1) = fabm_3d(:,:,1)
1695#else
1696               CALL ctl_stop( ' Trying to run skd490 observation operator', &
1697                  &           ' but no biogeochemical model appears to have been defined' )
1698#endif
1699
1700            CASE('sfco2')
1701#if defined key_hadocc
1702               zsurfvar(:,:,1) = HADOCC_FCO2(:,:)    ! fCO2 from HadOCC
1703               IF ( ( MINVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) .AND. &
1704                  & ( MAXVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) ) THEN
1705                  zsurfvar(:,:) = obfillflt
1706                  zsurfmask(:,:,1) = 0
1707                  CALL ctl_warn( ' HadOCC fCO2 values masked out for observation operator', &
1708                     &           ' as HADOCC_FCO2(:,:) == HADOCC_FILL_FLT' )
1709               ENDIF
1710#elif defined key_medusa && defined key_roam
1711               zsurfvar(:,:,1) = f2_fco2w(:,:)
1712#elif defined key_fabm
1713               ! First, get pCO2 from FABM
1714               fabm_3d(:,:,:) = model%get_interior_diagnostic_data(jp_fabm_o3pc)
1715               zsurfvar(:,:,1) = fabm_3d(:,:,1)
1716               ! Now, convert pCO2 to fCO2, based on SST in K. This follows the standard methodology of:
1717               ! Pierrot et al. (2009), Recommendations for autonomous underway pCO2 measuring systems
1718               ! and data reduction routines, Deep-Sea Research II, 56: 512-522.
1719               ! and
1720               ! Weiss (1974), Carbon dioxide in water and seawater: the solubility of a non-ideal gas,
1721               ! Marine Chemistry, 2: 203-215.
1722               ! In the implementation below, atmospheric pressure has been assumed to be 1 atm and so
1723               ! not explicitly included - atmospheric pressure is not necessarily available so this is
1724               ! the best assumption.
1725               ! Further, the (1-xCO2)^2 term has been neglected. This is common practice
1726               ! (see e.g. Zeebe and Wolf-Gladrow (2001), CO2 in Seawater: Equilibrium, Kinetics, Isotopes)
1727               ! because xCO2 in atm is ~0, and so this term will only affect the result to the 3rd decimal
1728               ! place for typical values, and xCO2 would need to be approximated from pCO2 anyway.
1729               zsurfvar(:,:,1) = zsurfvar(:,:,1) * EXP((-1636.75                                                          + &
1730                  &              12.0408      * (tsn(:,:,1,jp_tem)+rt0)                                                 - &
1731                  &              0.0327957    * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0)                         + &
1732                  &              0.0000316528 * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0) + &
1733                  &              2.0 * (57.7 - 0.118 * (tsn(:,:,1,jp_tem)+rt0)))                                        / &
1734                  &              (82.0578 * (tsn(:,:,1,jp_tem)+rt0)))
1735#else
1736               CALL ctl_stop( ' Trying to run sfco2 observation operator', &
1737                  &           ' but no biogeochemical model appears to have been defined' )
1738#endif
1739
1740            CASE('spco2')
1741#if defined key_hadocc
1742               zsurfvar(:,:,1) = HADOCC_PCO2(:,:)    ! pCO2 from HadOCC
1743               IF ( ( MINVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) .AND. &
1744                  & ( MAXVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) ) THEN
1745                  zsurfvar(:,:,1) = obfillflt
1746                  zsurfmask(:,:,1) = 0
1747                  CALL ctl_warn( ' HadOCC pCO2 values masked out for observation operator', &
1748                     &           ' as HADOCC_PCO2(:,:) == HADOCC_FILL_FLT' )
1749               ENDIF
1750#elif defined key_medusa && defined key_roam
1751               zsurfvar(:,:,1) = f2_pco2w(:,:)
1752#elif defined key_fabm
1753               fabm_3d(:,:,:) = model%get_interior_diagnostic_data(jp_fabm_o3pc)
1754               zsurfvar(:,:,1) = fabm_3d(:,:,1)
1755#else
1756               CALL ctl_stop( ' Trying to run spco2 observation operator', &
1757                  &           ' but no biogeochemical model appears to have been defined' )
1758#endif
1759
1760            CASE DEFAULT
1761
1762               CALL ctl_stop( 'Unknown surface observation type '//TRIM(cobstypessurf(jtype))//' in dia_obs' )
1763
1764            END SELECT
1765           
1766            IF ( llog10 ) THEN
1767               ! Take the log10 where we can, otherwise exclude
1768               tiny = 1.0e-20
1769               WHERE(zsurfvar(:,:,1) > tiny .AND. zsurfvar(:,:,1) /= obfillflt )
1770                  zsurfvar(:,:,1)  = LOG10(zsurfvar(:,:,1))
1771               ELSEWHERE
1772                  zsurfvar(:,:,1)  = obfillflt
1773                  zsurfmask(:,:,1) = 0
1774               END WHERE
1775            ENDIF
1776
1777            DO jvar = 1, surfdataqc(jtype)%nvar
1778
1779               IF ( TRIM(cobstypessurf(jtype)) == 'sla' .AND.                 &
1780                     &  ln_time_mean_sla_bkg ) THEN
1781                  !Number of time-steps in meaning period
1782                  imeanstp = NINT( ( MeanPeriodHours * 60. * 60. ) / rdt )
1783                  CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj,       &
1784                     &               nit000, idaystp, jvar,                   &
1785                     &               zsurfvar(:,:,jvar),                      &
1786                     &               zsurfclim(:,:,jvar),                     &
1787                     &               zsurfmask(:,:,jvar),                     &
1788                     &               zglam(:,:,jvar), zgphi(:,:,jvar),        &                     
1789                     &               n2dintsurf(jtype), llnightav(jtype),     &
1790                     &               ravglamscl(jtype), ravgphiscl(jtype),    &
1791                     &               lfpindegs(jtype), kmeanstp = imeanstp )
1792
1793               ELSE
1794                  CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj,       &
1795                     &               nit000, idaystp, jvar,                   &
1796                     &               zsurfvar(:,:,jvar),                      &
1797                     &               zsurfclim(:,:,jvar),                     &
1798                     &               zsurfmask(:,:,jvar),                     &
1799                     &               zglam(:,:,jvar), zgphi(:,:,jvar),        &                                         
1800                     &               n2dintsurf(jtype), llnightav(jtype),     &
1801                     &               ravglamscl(jtype), ravgphiscl(jtype),    &
1802                     &               lfpindegs(jtype) )
1803               ENDIF
1804
1805            END DO
1806           
1807            ! Change label of data from FBD ("freeboard") to SIT ("Sea Ice
1808            ! Thickness")
1809            IF ( TRIM(surfdataqc(jtype)%cvars(1)) == 'FBD' ) THEN
1810                 surfdataqc(jtype)%cvars(1) = 'SIT'
1811            ENDIF
1812
1813            CALL wrk_dealloc( jpi, jpj, surfdataqc(jtype)%nvar, zsurfvar  )
1814            CALL wrk_dealloc( jpi, jpj, surfdataqc(jtype)%nvar, zsurfclim )                 
1815            CALL wrk_dealloc( jpi, jpj, surfdataqc(jtype)%nvar, zsurfmask )
1816            CALL wrk_dealloc( jpi, jpj, surfdataqc(jtype)%nvar, zglam     )
1817            CALL wrk_dealloc( jpi, jpj, surfdataqc(jtype)%nvar, zgphi     )         
1818
1819         END DO
1820#if defined key_fabm
1821         CALL wrk_dealloc( jpi, jpj, jpk, fabm_3d )
1822#endif
1823
1824      ENDIF
1825
1826   END SUBROUTINE dia_obs
1827
1828   SUBROUTINE dia_obs_wri
1829      !!----------------------------------------------------------------------
1830      !!                    ***  ROUTINE dia_obs_wri  ***
1831      !!         
1832      !! ** Purpose : Call observation diagnostic output routines
1833      !!
1834      !! ** Method  : Call observation diagnostic output routines
1835      !!
1836      !! ** Action  :
1837      !!
1838      !! History :
1839      !!        !  06-03  (K. Mogensen) Original code
1840      !!        !  06-05  (K. Mogensen) Reformatted
1841      !!        !  06-10  (A. Weaver) Cleaning
1842      !!        !  07-03  (K. Mogensen) General handling of profiles
1843      !!        !  08-09  (M. Valdivieso) Velocity component (U,V) profiles
1844      !!        !  15-08  (M. Martin) Combined writing for prof and surf types
1845      !!----------------------------------------------------------------------
1846      !! * Modules used
1847      USE obs_rot_vel          ! Rotation of velocities
1848
1849      IMPLICIT NONE
1850
1851      !! * Local declarations
1852      INTEGER :: jtype                    ! Data set loop variable
1853      INTEGER :: jo, jvar, jk
1854      REAL(wp), DIMENSION(:), ALLOCATABLE :: &
1855         & zu, &
1856         & zv
1857
1858      !-----------------------------------------------------------------------
1859      ! Depending on switches call various observation output routines
1860      !-----------------------------------------------------------------------
1861
1862      IF ( nproftypes > 0 ) THEN
1863
1864         DO jtype = 1, nproftypes
1865
1866            IF ( TRIM(cobstypesprof(jtype)) == 'vel' ) THEN
1867
1868               ! For velocity data, rotate the model velocities to N/S, E/W
1869               ! using the compressed data structure.
1870               ALLOCATE( &
1871                  & zu(profdataqc(jtype)%nvprot(1)), &
1872                  & zv(profdataqc(jtype)%nvprot(2))  &
1873                  & )
1874
1875               CALL obs_rotvel_pro( profdataqc(jtype), nn_2dint_default, zu, zv )
1876
1877               DO jo = 1, profdataqc(jtype)%nprof
1878                  DO jvar = 1, 2
1879                     DO jk = profdataqc(jtype)%npvsta(jo,jvar), profdataqc(jtype)%npvend(jo,jvar)
1880
1881                        IF ( jvar == 1 ) THEN
1882                           profdataqc(jtype)%var(jvar)%vmod(jk) = zu(jk)
1883                        ELSE
1884                           profdataqc(jtype)%var(jvar)%vmod(jk) = zv(jk)
1885                        ENDIF
1886
1887                     END DO
1888                  END DO
1889               END DO
1890
1891               DEALLOCATE( zu )
1892               DEALLOCATE( zv )
1893
1894            END IF
1895
1896            CALL obs_prof_decompress( profdataqc(jtype), &
1897               &                      profdata(jtype), .TRUE., numout )
1898
1899            CALL obs_wri_prof( profdata(jtype) )
1900
1901         END DO
1902
1903      ENDIF
1904
1905      IF ( nsurftypes > 0 ) THEN
1906
1907         DO jtype = 1, nsurftypes
1908
1909            IF ( TRIM(cobstypessurf(jtype)) == 'vel' ) THEN
1910
1911               ! For velocity data, rotate the model velocities to N/S, E/W
1912               ! using the compressed data structure.
1913               ALLOCATE( &
1914                  & zu(surfdataqc(jtype)%nsurf), &
1915                  & zv(surfdataqc(jtype)%nsurf)  &
1916                  & )
1917
1918               CALL obs_rotvel_surf( surfdataqc(jtype), nn_2dint_default, zu, zv )
1919
1920               DO jo = 1, surfdataqc(jtype)%nsurf
1921                  surfdataqc(jtype)%rmod(jo,1) = zu(jo)
1922                  surfdataqc(jtype)%rmod(jo,2) = zv(jo)
1923               END DO
1924
1925               DEALLOCATE( zu )
1926               DEALLOCATE( zv )
1927
1928            END IF
1929
1930
1931            CALL obs_surf_decompress( surfdataqc(jtype), &
1932               &                      surfdata(jtype), .TRUE., numout )
1933
1934            CALL obs_wri_surf( surfdata(jtype) )
1935
1936         END DO
1937
1938      ENDIF
1939
1940   END SUBROUTINE dia_obs_wri
1941
1942   SUBROUTINE dia_obs_dealloc
1943      IMPLICIT NONE
1944      !!----------------------------------------------------------------------
1945      !!                    *** ROUTINE dia_obs_dealloc ***
1946      !!
1947      !!  ** Purpose : To deallocate data to enable the obs_oper online loop.
1948      !!               Specifically: dia_obs_init --> dia_obs --> dia_obs_wri
1949      !!
1950      !!  ** Method : Clean up various arrays left behind by the obs_oper.
1951      !!
1952      !!  ** Action :
1953      !!
1954      !!----------------------------------------------------------------------
1955      ! obs_grid deallocation
1956      CALL obs_grid_deallocate
1957
1958      ! diaobs deallocation
1959      IF ( nproftypes > 0 ) &
1960         &   DEALLOCATE( cobstypesprof, profdata, profdataqc, nvarsprof, nextrprof )
1961
1962      IF ( nsurftypes > 0 ) &
1963         &   DEALLOCATE( cobstypessurf, surfdata, surfdataqc, nvarssurf, nextrsurf, &
1964         &               n2dintsurf, ravglamscl, ravgphiscl, lfpindegs, llnightav )
1965
1966   END SUBROUTINE dia_obs_dealloc
1967
1968   SUBROUTINE ini_date( ddobsini )
1969      !!----------------------------------------------------------------------
1970      !!                    ***  ROUTINE ini_date  ***
1971      !!
1972      !! ** Purpose : Get initial date in double precision YYYYMMDD.HHMMSS format
1973      !!
1974      !! ** Method  : Get initial date in double precision YYYYMMDD.HHMMSS format
1975      !!
1976      !! ** Action  : Get initial date in double precision YYYYMMDD.HHMMSS format
1977      !!
1978      !! History :
1979      !!        !  06-03  (K. Mogensen)  Original code
1980      !!        !  06-05  (K. Mogensen)  Reformatted
1981      !!        !  06-10  (A. Weaver) Cleaning
1982      !!        !  06-10  (G. Smith) Calculates initial date the same as method for final date
1983      !!        !  10-05  (D. Lea) Update to month length calculation for NEMO vn3.2
1984      !!----------------------------------------------------------------------
1985      USE phycst, ONLY : &            ! Physical constants
1986         & rday
1987      USE dom_oce, ONLY : &           ! Ocean space and time domain variables
1988         & rdt
1989
1990      IMPLICIT NONE
1991
1992      !! * Arguments
1993      REAL(dp), INTENT(OUT) :: ddobsini  ! Initial date in YYYYMMDD.HHMMSS
1994
1995      !! * Local declarations
1996      INTEGER :: iyea        ! date - (year, month, day, hour, minute)
1997      INTEGER :: imon
1998      INTEGER :: iday
1999      INTEGER :: ihou
2000      INTEGER :: imin
2001      INTEGER :: imday       ! Number of days in month.
2002      INTEGER, DIMENSION(12) :: &
2003         &       imonth_len  ! Length in days of the months of the current year
2004      REAL(wp) :: zdayfrc    ! Fraction of day
2005
2006      !----------------------------------------------------------------------
2007      ! Initial date initialization (year, month, day, hour, minute)
2008      ! (This assumes that the initial date is for 00z))
2009      !----------------------------------------------------------------------
2010      iyea =   ndate0 / 10000
2011      imon = ( ndate0 - iyea * 10000 ) / 100
2012      iday =   ndate0 - iyea * 10000 - imon * 100
2013      ihou = 0
2014      imin = 0
2015
2016      !----------------------------------------------------------------------
2017      ! Compute number of days + number of hours + min since initial time
2018      !----------------------------------------------------------------------
2019      iday = iday + ( nit000 -1 ) * rdt / rday
2020      zdayfrc = ( nit000 -1 ) * rdt / rday
2021      zdayfrc = zdayfrc - aint(zdayfrc)
2022      ihou = int( zdayfrc * 24 )
2023      imin = int( (zdayfrc * 24 - ihou) * 60 )
2024
2025      !-----------------------------------------------------------------------
2026      ! Convert number of days (iday) into a real date
2027      !----------------------------------------------------------------------
2028
2029      CALL calc_month_len( iyea, imonth_len )
2030
2031      DO WHILE ( iday > imonth_len(imon) )
2032         iday = iday - imonth_len(imon)
2033         imon = imon + 1 
2034         IF ( imon > 12 ) THEN
2035            imon = 1
2036            iyea = iyea + 1
2037            CALL calc_month_len( iyea, imonth_len )  ! update month lengths
2038         ENDIF
2039      END DO
2040
2041      !----------------------------------------------------------------------
2042      ! Convert it into YYYYMMDD.HHMMSS format.
2043      !----------------------------------------------------------------------
2044      ddobsini = iyea * 10000_dp + imon * 100_dp + &
2045         &       iday + ihou * 0.01_dp + imin * 0.0001_dp
2046
2047
2048   END SUBROUTINE ini_date
2049
2050   SUBROUTINE fin_date( ddobsfin )
2051      !!----------------------------------------------------------------------
2052      !!                    ***  ROUTINE fin_date  ***
2053      !!
2054      !! ** Purpose : Get final date in double precision YYYYMMDD.HHMMSS format
2055      !!
2056      !! ** Method  : Get final date in double precision YYYYMMDD.HHMMSS format
2057      !!
2058      !! ** Action  : Get final date in double precision YYYYMMDD.HHMMSS format
2059      !!
2060      !! History :
2061      !!        !  06-03  (K. Mogensen)  Original code
2062      !!        !  06-05  (K. Mogensen)  Reformatted
2063      !!        !  06-10  (A. Weaver) Cleaning
2064      !!        !  10-05  (D. Lea) Update to month length calculation for NEMO vn3.2
2065      !!----------------------------------------------------------------------
2066      USE phycst, ONLY : &            ! Physical constants
2067         & rday
2068      USE dom_oce, ONLY : &           ! Ocean space and time domain variables
2069         & rdt
2070
2071      IMPLICIT NONE
2072
2073      !! * Arguments
2074      REAL(dp), INTENT(OUT) :: ddobsfin ! Final date in YYYYMMDD.HHMMSS
2075
2076      !! * Local declarations
2077      INTEGER :: iyea        ! date - (year, month, day, hour, minute)
2078      INTEGER :: imon
2079      INTEGER :: iday
2080      INTEGER :: ihou
2081      INTEGER :: imin
2082      INTEGER :: imday       ! Number of days in month.
2083      INTEGER, DIMENSION(12) :: &
2084         &       imonth_len  ! Length in days of the months of the current year
2085      REAL(wp) :: zdayfrc    ! Fraction of day
2086
2087      !-----------------------------------------------------------------------
2088      ! Initial date initialization (year, month, day, hour, minute)
2089      ! (This assumes that the initial date is for 00z)
2090      !-----------------------------------------------------------------------
2091      iyea =   ndate0 / 10000
2092      imon = ( ndate0 - iyea * 10000 ) / 100
2093      iday =   ndate0 - iyea * 10000 - imon * 100
2094      ihou = 0
2095      imin = 0
2096
2097      !-----------------------------------------------------------------------
2098      ! Compute number of days + number of hours + min since initial time
2099      !-----------------------------------------------------------------------
2100      iday    = iday +  nitend  * rdt / rday
2101      zdayfrc =  nitend  * rdt / rday
2102      zdayfrc = zdayfrc - AINT( zdayfrc )
2103      ihou    = INT( zdayfrc * 24 )
2104      imin    = INT( ( zdayfrc * 24 - ihou ) * 60 )
2105
2106      !-----------------------------------------------------------------------
2107      ! Convert number of days (iday) into a real date
2108      !----------------------------------------------------------------------
2109
2110      CALL calc_month_len( iyea, imonth_len )
2111
2112      DO WHILE ( iday > imonth_len(imon) )
2113         iday = iday - imonth_len(imon)
2114         imon = imon + 1 
2115         IF ( imon > 12 ) THEN
2116            imon = 1
2117            iyea = iyea + 1
2118            CALL calc_month_len( iyea, imonth_len )  ! update month lengths
2119         ENDIF
2120      END DO
2121
2122      !-----------------------------------------------------------------------
2123      ! Convert it into YYYYMMDD.HHMMSS format
2124      !-----------------------------------------------------------------------
2125      ddobsfin = iyea * 10000_dp + imon * 100_dp    + iday &
2126         &     + ihou * 0.01_dp  + imin * 0.0001_dp
2127
2128    END SUBROUTINE fin_date
2129
2130    SUBROUTINE obs_settypefiles( ntypes, jpmaxnfiles, ifiles, cobstypes, cfiles )
2131
2132       INTEGER, INTENT(IN) :: ntypes      ! Total number of obs types
2133       INTEGER, INTENT(IN) :: jpmaxnfiles ! Maximum number of files allowed for each type
2134       INTEGER, DIMENSION(ntypes), INTENT(OUT) :: &
2135          &                   ifiles      ! Out number of files for each type
2136       CHARACTER(len=8), DIMENSION(ntypes), INTENT(IN) :: &
2137          &                   cobstypes   ! List of obs types
2138       CHARACTER(len=128), DIMENSION(ntypes, jpmaxnfiles), INTENT(IN) :: &
2139          &                   cfiles      ! List of files for all types
2140
2141       !Local variables
2142       INTEGER :: jfile
2143       INTEGER :: jtype
2144
2145       DO jtype = 1, ntypes
2146
2147          ifiles(jtype) = 0
2148          DO jfile = 1, jpmaxnfiles
2149             IF ( trim(cfiles(jtype,jfile)) /= '' ) &
2150                       ifiles(jtype) = ifiles(jtype) + 1
2151          END DO
2152
2153          IF ( ifiles(jtype) == 0 ) THEN
2154               CALL ctl_stop( 'Logical for observation type '//TRIM(cobstypes(jtype))//   &
2155                  &           ' set to true but no files available to read' )
2156          ENDIF
2157
2158          IF(lwp) THEN   
2159             WRITE(numout,*) '             '//cobstypes(jtype)//' input observation file names:'
2160             DO jfile = 1, ifiles(jtype)
2161                WRITE(numout,*) '                '//TRIM(cfiles(jtype,jfile))
2162             END DO
2163          ENDIF
2164
2165       END DO
2166
2167    END SUBROUTINE obs_settypefiles
2168
2169    SUBROUTINE obs_setinterpopts( ntypes, jtype, ctypein,             &
2170               &                  n2dint_default, n2dint_type,        &
2171               &                  ravglamscl_type, ravgphiscl_type,   &
2172               &                  lfp_indegs_type, lavnight_type,     &
2173               &                  n2dint, ravglamscl, ravgphiscl,     &
2174               &                  lfpindegs, lavnight )
2175
2176       INTEGER, INTENT(IN)  :: ntypes             ! Total number of obs types
2177       INTEGER, INTENT(IN)  :: jtype              ! Index of the current type of obs
2178       INTEGER, INTENT(IN)  :: n2dint_default     ! Default option for interpolation type
2179       INTEGER, INTENT(IN)  :: n2dint_type        ! Option for interpolation type
2180       REAL(wp), INTENT(IN) :: &
2181          &                    ravglamscl_type, & !E/W diameter of obs footprint for this type
2182          &                    ravgphiscl_type    !N/S diameter of obs footprint for this type
2183       LOGICAL, INTENT(IN)  :: lfp_indegs_type    !T=> footprint in degrees, F=> in metres
2184       LOGICAL, INTENT(IN)  :: lavnight_type      !T=> obs represent night time average
2185       CHARACTER(len=8), INTENT(IN) :: ctypein 
2186
2187       INTEGER, DIMENSION(ntypes), INTENT(INOUT) :: &
2188          &                    n2dint 
2189       REAL(wp), DIMENSION(ntypes), INTENT(INOUT) :: &
2190          &                    ravglamscl, ravgphiscl
2191       LOGICAL, DIMENSION(ntypes), INTENT(INOUT) :: &
2192          &                    lfpindegs, lavnight
2193
2194       lavnight(jtype) = lavnight_type
2195
2196       IF ( (n2dint_type >= 0) .AND. (n2dint_type <= 6) ) THEN
2197          n2dint(jtype) = n2dint_type
2198       ELSE IF ( n2dint_type == -1 ) THEN
2199          n2dint(jtype) = n2dint_default
2200       ELSE
2201          CALL ctl_stop(' Choice of '//TRIM(ctypein)//' horizontal (2D) interpolation method', &
2202            &                    ' is not available')
2203       ENDIF
2204
2205       ! For averaging observation footprints set options for size of footprint
2206       IF ( (n2dint(jtype) > 4) .AND. (n2dint(jtype) <= 6) ) THEN
2207          IF ( ravglamscl_type > 0._wp ) THEN
2208             ravglamscl(jtype) = ravglamscl_type
2209          ELSE
2210             CALL ctl_stop( 'Incorrect value set for averaging footprint '// &
2211                            'scale (ravglamscl) for observation type '//TRIM(ctypein) )     
2212          ENDIF
2213
2214          IF ( ravgphiscl_type > 0._wp ) THEN
2215             ravgphiscl(jtype) = ravgphiscl_type
2216          ELSE
2217             CALL ctl_stop( 'Incorrect value set for averaging footprint '// &
2218                            'scale (ravgphiscl) for observation type '//TRIM(ctypein) )     
2219          ENDIF
2220
2221          lfpindegs(jtype) = lfp_indegs_type 
2222
2223       ENDIF
2224
2225       ! Write out info
2226       IF(lwp) THEN
2227          IF ( n2dint(jtype) <= 4 ) THEN
2228             WRITE(numout,*) '             '//TRIM(ctypein)// &
2229                &            ' model counterparts will be interpolated horizontally'
2230          ELSE IF ( n2dint(jtype) <= 6 ) THEN
2231             WRITE(numout,*) '             '//TRIM(ctypein)// &
2232                &            ' model counterparts will be averaged horizontally'
2233             WRITE(numout,*) '             '//'    with E/W scale: ',ravglamscl(jtype)
2234             WRITE(numout,*) '             '//'    with N/S scale: ',ravgphiscl(jtype)
2235             IF ( lfpindegs(jtype) ) THEN
2236                 WRITE(numout,*) '             '//'    (in degrees)'
2237             ELSE
2238                 WRITE(numout,*) '             '//'    (in metres)'
2239             ENDIF
2240          ENDIF
2241       ENDIF
2242
2243    END SUBROUTINE obs_setinterpopts
2244
2245END MODULE diaobs
Note: See TracBrowser for help on using the repository browser.