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 NEMO/branches/UKMO/r8395_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

source: NEMO/branches/UKMO/r8395_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 11361

Last change on this file since 11361 was 11361, checked in by jcastill, 5 years ago

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