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

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

source: branches/UKMO/dev_r5518_obs_oper_update_addclim/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 11455

Last change on this file since 11455 was 11455, checked in by mattmartin, 5 years ago

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