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

Last change on this file since 10181 was 10181, checked in by emmafiedler, 3 years ago

Working version of ice thickness assimilation updates

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