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

Last change on this file since 12610 was 12610, checked in by dcarneir, 19 months ago

Inclusion of sea ice thickness in OBS branch

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