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

Context Navigation

source: branches/UKMO/dev_r5518_obs_oper_update_kd490/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 11253

Last change on this file since 11253 was 11253, checked in by dford, 5 years ago
Add Kd490 obs operator.
File size: 88.9 KB

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

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

Download in other formats: