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

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

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

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

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