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

Context Navigation

source: branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 10247

Last change on this file since 10247 was 10247, checked in by kingr, 6 years ago
Merged dev_r5518_obs_oper_update@10168
File size: 86.7 KB

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

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

Download in other formats: