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

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

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

Context Navigation

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

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

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

Download in other formats: