Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

diaobs.F90 @ 13384

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

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

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

Context Navigation

source: branches/UKMO/dev_r5518_obs_oper_update_vel_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 13384

Download in other formats: