source: branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 11468

Last change on this file since 11468 was 11468, checked in by mattmartin, 2 years ago

Merged changes to allow writing of climatological information to feedback files.

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