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

Context Navigation

source: branches/UKMO/dev_r5518_obs_oper_update_utils366_fabmv1/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 13240

Last change on this file since 13240 was 13240, checked in by dford, 4 years ago
Update FABM syntax for v1.0, and fix bug where ERSEM pH was being taken from the wrong place.
File size: 99.5 KB

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

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

Download in other formats: