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

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

source: branches/UKMO/dev_r5518_obs_oper_update_sit/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 12364

Last change on this file since 12364 was 12364, checked in by dcarneir, 4 years ago

Modifications in diaobs.F90 to include SIT flags

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