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

Context Navigation

source: branches/UKMO/dev_r5518_obs_oper_update_sit_SMOS/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 12820

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

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

Download in other formats: