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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

diaobs.F90 in branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

Context Navigation

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

Last change on this file since 12140 was 12140, checked in by kingr, 4 years ago
Added option to remove tides from SLA bkg by taking average over 24h50m.
File size: 95.9 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: