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

Context Navigation

source: branches/UKMO/dev_r5518_obs_oper_update_bgc3d/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 9192

Last change on this file since 9192 was 9192, checked in by dford, 6 years ago
Reduce duplication of log transforms.
File size: 83.6 KB

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

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

Download in other formats: