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

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

source: branches/UKMO/dev_r5518_obs_oper_update_bgc3d/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 9191

Last change on this file since 9191 was 9191, checked in by dford, 6 years ago

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