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

Context Navigation

source: branches/UKMO/dev_r5518_obs_oper_update_SlaAssim/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 11863

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

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

Download in other formats: