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_icethick/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

source: branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 10276

Last change on this file since 10276 was 10276, checked in by emmafiedler, 4 years ago

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