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

Context Navigation

source: branches/UKMO/dev_r5518_obs_oper_update_sit/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 11953

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

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

Download in other formats: