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

Context Navigation

source: branches/UKMO/dev_r5518_obs_oper_update_bgc3d/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 9205

Last change on this file since 9205 was 9205, checked in by dford, 6 years ago
Change obs_pre_prof to loop over an arbitrary number of variables.
File size: 84.0 KB

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

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

Download in other formats: