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 @ 11269

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

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

Download in other formats: