Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

diaobs.F90

Last change on this file was 9013, checked in by dford, 6 years ago
Update to read in and write out logchl obs STD values. See internal Met Office NEMO ticket 736.
File size: 89.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	!! 'key_diaobs' : Switch on the observation diagnostic computation
10	!!----------------------------------------------------------------------
11	!! dia_obs_init : Reading and prepare observations
12	!! dia_obs : Compute model equivalent to observations
13	!! dia_obs_wri : Write observational diagnostics
14	!! ini_date : Compute the initial date YYYYMMDD.HHMMSS
15	!! fin_date : Compute the final date YYYYMMDD.HHMMSS
16	!!----------------------------------------------------------------------
17	!! * Modules used
18	USE wrk_nemo ! Memory Allocation
19	USE par_kind ! Precision variables
20	USE in_out_manager ! I/O manager
21	USE par_oce
22	USE dom_oce ! Ocean space and time domain variables
23	USE obs_const, ONLY: obfillflt ! Fill value
24	USE obs_fbm, ONLY: ln_cl4 ! Class 4 diagnostic switch
25	USE obs_read_prof ! Reading and allocation of observations (Coriolis)
26	USE obs_read_sla ! Reading and allocation of SLA observations
27	USE obs_read_sst ! Reading and allocation of SST observations
28	USE obs_sstbias ! Bias correction routine for SST
29	USE obs_readmdt ! Reading and allocation of MDT for SLA.
30	USE obs_read_seaice ! Reading and allocation of Sea Ice observations
31	USE obs_read_vel ! Reading and allocation of velocity component observations
32	USE obs_read_logchl ! Reading and allocation of logchl observations
33	USE obs_read_spm ! Reading and allocation of spm observations
34	USE obs_read_fco2 ! Reading and allocation of fco2 observations
35	USE obs_read_pco2 ! Reading and allocation of pco2 observations
36	USE obs_prep ! Preparation of obs. (grid search etc).
37	USE obs_oper ! Observation operators
38	USE obs_write ! Writing of observation related diagnostics
39	USE obs_grid ! Grid searching
40	USE obs_read_altbias ! Bias treatment for altimeter
41	USE obs_profiles_def ! Profile data definitions
42	USE obs_profiles ! Profile data storage
43	USE obs_surf_def ! Surface data definitions
44	USE obs_sla ! SLA data storage
45	USE obs_sst ! SST data storage
46	USE obs_seaice ! Sea Ice data storage
47	USE obs_logchl ! logchl data storage
48	USE obs_spm ! spm data storage
49	USE obs_fco2 ! fco2 data storage
50	USE obs_pco2 ! pco2 data storage
51	USE obs_types ! Definitions for observation types
52	USE mpp_map ! MPP mapping
53	USE lib_mpp ! For ctl_warn/stop
54
55	IMPLICIT NONE
56
57	!! * Routine accessibility
58	PRIVATE
59	PUBLIC dia_obs_init, & ! Initialize and read observations
60	& dia_obs, & ! Compute model equivalent to observations
61	& dia_obs_wri, & ! Write model equivalent to observations
62	& dia_obs_dealloc ! Deallocate dia_obs data
63
64	!! * Shared Module variables
65	LOGICAL, PUBLIC, PARAMETER :: &
66	#if defined key_diaobs
67	& lk_diaobs = .TRUE. !: Logical switch for observation diangostics
68	#else
69	& lk_diaobs = .FALSE. !: Logical switch for observation diangostics
70	#endif
71
72	!! * Module variables
73	LOGICAL, PUBLIC :: ln_t3d !: Logical switch for temperature profiles
74	LOGICAL, PUBLIC :: ln_s3d !: Logical switch for salinity profiles
75	LOGICAL, PUBLIC :: ln_ena !: Logical switch for the ENACT data set
76	LOGICAL, PUBLIC :: ln_cor !: Logical switch for the Coriolis data set
77	LOGICAL, PUBLIC :: ln_profb !: Logical switch for profile feedback datafiles
78	LOGICAL, PUBLIC :: ln_sla !: Logical switch for sea level anomalies
79	LOGICAL, PUBLIC :: ln_sladt !: Logical switch for SLA from AVISO files
80	LOGICAL, PUBLIC :: ln_slafb !: Logical switch for SLA from feedback files
81	LOGICAL, PUBLIC :: ln_sst !: Logical switch for sea surface temperature
82	LOGICAL, PUBLIC :: ln_reysst !: Logical switch for Reynolds sea surface temperature
83	LOGICAL, PUBLIC :: ln_ghrsst !: Logical switch for GHRSST data
84	LOGICAL, PUBLIC :: ln_sstfb !: Logical switch for SST from feedback files
85	LOGICAL, PUBLIC :: ln_seaice !: Logical switch for sea ice concentration
86	LOGICAL, PUBLIC :: ln_vel3d !: Logical switch for velocity component (u,v) observations
87	LOGICAL, PUBLIC :: ln_velavcur !: Logical switch for raw daily averaged netCDF current meter vel. data
88	LOGICAL, PUBLIC :: ln_velhrcur !: Logical switch for raw high freq netCDF current meter vel. data
89	LOGICAL, PUBLIC :: ln_velavadcp !: Logical switch for raw daily averaged netCDF ADCP vel. data
90	LOGICAL, PUBLIC :: ln_velhradcp !: Logical switch for raw high freq netCDF ADCP vel. data
91	LOGICAL, PUBLIC :: ln_velfb !: Logical switch for velocities from feedback files
92	LOGICAL, PUBLIC :: ln_logchl !: Logical switch for log10(chlorophyll)
93	LOGICAL, PUBLIC :: ln_logchlfb !: Logical switch for logchl from feedback files
94	LOGICAL, PUBLIC :: ln_spm !: Logical switch for spm
95	LOGICAL, PUBLIC :: ln_spmfb !: Logical switch for spm from feedback files
96	LOGICAL, PUBLIC :: ln_fco2 !: Logical switch for fco2
97	LOGICAL, PUBLIC :: ln_fco2fb !: Logical switch for fco2 from feedback files
98	LOGICAL, PUBLIC :: ln_pco2 !: Logical switch for pco2
99	LOGICAL, PUBLIC :: ln_pco2fb !: Logical switch for pco2 from feedback files
100	LOGICAL, PUBLIC :: ln_ssh !: Logical switch for sea surface height
101	LOGICAL, PUBLIC :: ln_sss !: Logical switch for sea surface salinity
102	LOGICAL, PUBLIC :: ln_sstnight !: Logical switch for night mean SST observations
103	LOGICAL, PUBLIC :: ln_nea !: Remove observations near land
104	LOGICAL, PUBLIC :: ln_altbias !: Logical switch for altimeter bias
105	LOGICAL, PUBLIC :: ln_ignmis !: Logical switch for ignoring missing files
106	LOGICAL, PUBLIC :: ln_s_at_t !: Logical switch to compute model S at T observations
107	LOGICAL, PUBLIC :: ln_sstbias !: Logical switch for bias corection of SST
108
109	REAL(KIND=dp), PUBLIC :: dobsini !: Observation window start date YYYYMMDD.HHMMSS
110	REAL(KIND=dp), PUBLIC :: dobsend !: Observation window end date YYYYMMDD.HHMMSS
111
112	INTEGER, PUBLIC :: n1dint !: Vertical interpolation method
113	INTEGER, PUBLIC :: n2dint !: Horizontal interpolation method
114
115	INTEGER, DIMENSION(imaxavtypes) :: &
116	& endailyavtypes !: ENACT data types which are daily average
117
118	INTEGER, PARAMETER :: MaxNumFiles = 1000
119	LOGICAL, DIMENSION(MaxNumFiles) :: &
120	& ln_profb_ena, & !: Is the feedback files from ENACT data ?
121	! !: If so use endailyavtypes
122	& ln_profb_enatim !: Change tim for 820 enact data set.
123
124	INTEGER, DIMENSION(MaxNumFiles), PUBLIC :: sstbias_type !SST bias type
125
126	LOGICAL, DIMENSION(MaxNumFiles) :: &
127	& ln_velfb_av !: Is the velocity feedback files daily average?
128	LOGICAL, DIMENSION(:), ALLOCATABLE :: &
129	& ld_enact !: Profile data is ENACT so use endailyavtypes
130	LOGICAL, DIMENSION(:), ALLOCATABLE :: &
131	& ld_velav !: Velocity data is daily averaged
132	LOGICAL, DIMENSION(:), ALLOCATABLE :: &
133	& ld_sstnight !: SST observation corresponds to night mean
134
135	!!----------------------------------------------------------------------
136	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
137	!! $Id$
138	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
139	!!----------------------------------------------------------------------
140
141	!! * Substitutions
142	# include "domzgr_substitute.h90"
143	CONTAINS
144
145	SUBROUTINE dia_obs_init
146	!!----------------------------------------------------------------------
147	!! * ROUTINE dia_obs_init *
148	!!
149	!! ** Purpose : Initialize and read observations
150	!!
151	!! ** Method : Read the namelist and call reading routines
152	!!
153	!! ** Action : Read the namelist and call reading routines
154	!!
155	!! History :
156	!! ! 06-03 (K. Mogensen) Original code
157	!! ! 06-05 (A. Weaver) Reformatted
158	!! ! 06-10 (A. Weaver) Cleaning and add controls
159	!! ! 07-03 (K. Mogensen) General handling of profiles
160	!! ! 14-08 (J.While) Incorporated SST bias correction
161	!!----------------------------------------------------------------------
162
163	IMPLICIT NONE
164
165	!! * Local declarations
166	CHARACTER(len=128) :: enactfiles(MaxNumFiles)
167	CHARACTER(len=128) :: coriofiles(MaxNumFiles)
168	CHARACTER(len=128) :: profbfiles(MaxNumFiles)
169	CHARACTER(len=128) :: sstfiles(MaxNumFiles)
170	CHARACTER(len=128) :: sstfbfiles(MaxNumFiles)
171	CHARACTER(len=128) :: sstbias_files(MaxNumFiles)
172	CHARACTER(len=128) :: slafilesact(MaxNumFiles)
173	CHARACTER(len=128) :: slafilespas(MaxNumFiles)
174	CHARACTER(len=128) :: slafbfiles(MaxNumFiles)
175	CHARACTER(len=128) :: seaicefiles(MaxNumFiles)
176	CHARACTER(len=128) :: velcurfiles(MaxNumFiles)
177	CHARACTER(len=128) :: veladcpfiles(MaxNumFiles)
178	CHARACTER(len=128) :: velavcurfiles(MaxNumFiles)
179	CHARACTER(len=128) :: velhrcurfiles(MaxNumFiles)
180	CHARACTER(len=128) :: velavadcpfiles(MaxNumFiles)
181	CHARACTER(len=128) :: velhradcpfiles(MaxNumFiles)
182	CHARACTER(len=128) :: velfbfiles(MaxNumFiles)
183	CHARACTER(len=128) :: logchlfiles(MaxNumFiles)
184	CHARACTER(len=128) :: logchlfbfiles(MaxNumFiles)
185	CHARACTER(len=128) :: spmfiles(MaxNumFiles)
186	CHARACTER(len=128) :: spmfbfiles(MaxNumFiles)
187	CHARACTER(len=128) :: fco2files(MaxNumFiles)
188	CHARACTER(len=128) :: fco2fbfiles(MaxNumFiles)
189	CHARACTER(len=128) :: pco2files(MaxNumFiles)
190	CHARACTER(len=128) :: pco2fbfiles(MaxNumFiles)
191	CHARACTER(LEN=128) :: reysstname
192	CHARACTER(LEN=12) :: reysstfmt
193	CHARACTER(LEN=128) :: bias_file
194	CHARACTER(LEN=20) :: datestr=" ", timestr=" "
195	NAMELIST/namobs/ln_ena, ln_cor, ln_profb, ln_t3d, ln_s3d, &
196	& ln_sla, ln_sladt, ln_slafb, &
197	& ln_ssh, ln_sst, ln_sstfb, ln_sss, ln_nea, &
198	& ln_bound_reject, &
199	& enactfiles, coriofiles, profbfiles, &
200	& slafilesact, slafilespas, slafbfiles, &
201	& sstfiles, sstfbfiles, &
202	& ln_seaice, seaicefiles, &
203	& dobsini, dobsend, n1dint, n2dint, &
204	& nmsshc, mdtcorr, mdtcutoff, &
205	& ln_reysst, ln_ghrsst, reysstname, reysstfmt, &
206	& ln_sstnight, &
207	& ln_grid_search_lookup, &
208	& grid_search_file, grid_search_res, &
209	& ln_grid_global, bias_file, ln_altbias, &
210	& endailyavtypes, ln_s_at_t, ln_profb_ena, &
211	& ln_vel3d, ln_velavcur, velavcurfiles, &
212	& ln_velhrcur, velhrcurfiles, &
213	& ln_velavadcp, velavadcpfiles, &
214	& ln_velhradcp, velhradcpfiles, &
215	& ln_velfb, velfbfiles, ln_velfb_av, &
216	& ln_logchl, ln_logchlfb, &
217	& logchlfiles, logchlfbfiles, &
218	& ln_spm, ln_spmfb, &
219	& spmfiles, spmfbfiles, &
220	& ln_fco2, ln_fco2fb, &
221	& fco2files, fco2fbfiles, &
222	& ln_pco2, ln_pco2fb, &
223	& pco2files, pco2fbfiles, &
224	& ln_profb_enatim, ln_ignmis, ln_cl4, &
225	& ln_sstbias, sstbias_files
226
227	INTEGER :: jprofset
228	INTEGER :: jveloset
229	INTEGER :: jvar
230	INTEGER :: jnumenact
231	INTEGER :: jnumcorio
232	INTEGER :: jnumprofb
233	INTEGER :: jnumslaact
234	INTEGER :: jnumslapas
235	INTEGER :: jnumslafb
236	INTEGER :: jnumsst
237	INTEGER :: jnumsstfb
238	INTEGER :: jnumsstbias
239	INTEGER :: jnumseaice
240	INTEGER :: jnumvelavcur
241	INTEGER :: jnumvelhrcur
242	INTEGER :: jnumvelavadcp
243	INTEGER :: jnumvelhradcp
244	INTEGER :: jnumvelfb
245	INTEGER :: jnumlogchl
246	INTEGER :: jnumlogchlfb
247	INTEGER :: jnumspm
248	INTEGER :: jnumspmfb
249	INTEGER :: jnumfco2
250	INTEGER :: jnumfco2fb
251	INTEGER :: jnumpco2
252	INTEGER :: jnumpco2fb
253	INTEGER :: ji
254	INTEGER :: jset
255	INTEGER :: ios ! Local integer output status for namelist read
256	LOGICAL :: lmask(MaxNumFiles), ll_u3d, ll_v3d
257
258	!-----------------------------------------------------------------------
259	! Read namelist parameters
260	!-----------------------------------------------------------------------
261
262	ln_logchl = .FALSE.
263	ln_logchlfb = .FALSE.
264	ln_spm = .FALSE.
265	ln_spmfb = .FALSE.
266	ln_fco2 = .FALSE.
267	ln_fco2fb = .FALSE.
268	ln_pco2 = .FALSE.
269	ln_pco2fb = .FALSE.
270
271	!Initalise all values in namelist arrays
272	enactfiles(:) = ''
273	coriofiles(:) = ''
274	profbfiles(:) = ''
275	slafilesact(:) = ''
276	slafilespas(:) = ''
277	slafbfiles(:) = ''
278	sstfiles(:) = ''
279	sstfbfiles(:) = ''
280	seaicefiles(:) = ''
281	velcurfiles(:) = ''
282	veladcpfiles(:) = ''
283	velavcurfiles(:) = ''
284	velhrcurfiles(:) = ''
285	velavadcpfiles(:) = ''
286	velhradcpfiles(:) = ''
287	velfbfiles(:) = ''
288	velcurfiles(:) = ''
289	veladcpfiles(:) = ''
290	logchlfiles(:) = ''
291	logchlfbfiles(:) = ''
292	spmfiles(:) = ''
293	spmfbfiles(:) = ''
294	fco2files(:) = ''
295	fco2fbfiles(:) = ''
296	pco2files(:) = ''
297	pco2fbfiles(:) = ''
298	sstbias_files(:) = ''
299	endailyavtypes(:) = -1
300	endailyavtypes(1) = 820
301	ln_profb_ena(:) = .FALSE.
302	ln_profb_enatim(:) = .TRUE.
303	ln_velfb_av(:) = .FALSE.
304	ln_ignmis = .FALSE.
305	ln_bound_reject = .TRUE.
306
307	! Read Namelist namobs : control observation diagnostics
308	REWIND( numnam_ref ) ! Namelist namobs in reference namelist : Diagnostic: control observation
309	READ ( numnam_ref, namobs, IOSTAT = ios, ERR = 901)
310	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in reference namelist', lwp )
311
312	CALL ini_date( dobsini )
313	CALL fin_date( dobsend )
314
315	REWIND( numnam_cfg ) ! Namelist namobs in configuration namelist : Diagnostic: control observation
316	READ ( numnam_cfg, namobs, IOSTAT = ios, ERR = 902 )
317	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in configuration namelist', lwp )
318	IF(lwm) WRITE ( numond, namobs )
319
320	! Count number of files for each type
321	IF (ln_ena) THEN
322	lmask(:) = .FALSE.
323	WHERE (enactfiles(:) /= '') lmask(:) = .TRUE.
324	jnumenact = COUNT(lmask)
325	ENDIF
326	IF (ln_cor) THEN
327	lmask(:) = .FALSE.
328	WHERE (coriofiles(:) /= '') lmask(:) = .TRUE.
329	jnumcorio = COUNT(lmask)
330	ENDIF
331	IF (ln_profb) THEN
332	lmask(:) = .FALSE.
333	WHERE (profbfiles(:) /= '') lmask(:) = .TRUE.
334	jnumprofb = COUNT(lmask)
335	ENDIF
336	IF (ln_sladt) THEN
337	lmask(:) = .FALSE.
338	WHERE (slafilesact(:) /= '') lmask(:) = .TRUE.
339	jnumslaact = COUNT(lmask)
340	lmask(:) = .FALSE.
341	WHERE (slafilespas(:) /= '') lmask(:) = .TRUE.
342	jnumslapas = COUNT(lmask)
343	ENDIF
344	IF (ln_slafb) THEN
345	lmask(:) = .FALSE.
346	WHERE (slafbfiles(:) /= '') lmask(:) = .TRUE.
347	jnumslafb = COUNT(lmask)
348	lmask(:) = .FALSE.
349	ENDIF
350	IF (ln_ghrsst) THEN
351	lmask(:) = .FALSE.
352	WHERE (sstfiles(:) /= '') lmask(:) = .TRUE.
353	jnumsst = COUNT(lmask)
354	ENDIF
355	IF (ln_sstfb) THEN
356	lmask(:) = .FALSE.
357	WHERE (sstfbfiles(:) /= '') lmask(:) = .TRUE.
358	jnumsstfb = COUNT(lmask)
359	lmask(:) = .FALSE.
360	ENDIF
361	IF (ln_sstbias) THEN
362	lmask(:) = .FALSE.
363	WHERE (sstbias_files(:) /= '') lmask(:) = .TRUE.
364	jnumsstbias = COUNT(lmask)
365	lmask(:) = .FALSE.
366	ENDIF
367	IF (ln_seaice) THEN
368	lmask(:) = .FALSE.
369	WHERE (seaicefiles(:) /= '') lmask(:) = .TRUE.
370	jnumseaice = COUNT(lmask)
371	ENDIF
372	IF (ln_velavcur) THEN
373	lmask(:) = .FALSE.
374	WHERE (velavcurfiles(:) /= '') lmask(:) = .TRUE.
375	jnumvelavcur = COUNT(lmask)
376	ENDIF
377	IF (ln_velhrcur) THEN
378	lmask(:) = .FALSE.
379	WHERE (velhrcurfiles(:) /= '') lmask(:) = .TRUE.
380	jnumvelhrcur = COUNT(lmask)
381	ENDIF
382	IF (ln_velavadcp) THEN
383	lmask(:) = .FALSE.
384	WHERE (velavadcpfiles(:) /= '') lmask(:) = .TRUE.
385	jnumvelavadcp = COUNT(lmask)
386	ENDIF
387	IF (ln_velhradcp) THEN
388	lmask(:) = .FALSE.
389	WHERE (velhradcpfiles(:) /= '') lmask(:) = .TRUE.
390	jnumvelhradcp = COUNT(lmask)
391	ENDIF
392	IF (ln_velfb) THEN
393	lmask(:) = .FALSE.
394	WHERE (velfbfiles(:) /= '') lmask(:) = .TRUE.
395	jnumvelfb = COUNT(lmask)
396	lmask(:) = .FALSE.
397	ENDIF
398	IF (ln_logchl) THEN
399	lmask(:) = .FALSE.
400	WHERE (logchlfiles(:) /= '') lmask(:) = .TRUE.
401	jnumlogchl = COUNT(lmask)
402	ENDIF
403	IF (ln_logchlfb) THEN
404	lmask(:) = .FALSE.
405	WHERE (logchlfbfiles(:) /= '') lmask(:) = .TRUE.
406	jnumlogchlfb = COUNT(lmask)
407	ENDIF
408	IF (ln_spm) THEN
409	lmask(:) = .FALSE.
410	WHERE (spmfiles(:) /= '') lmask(:) = .TRUE.
411	jnumspm = COUNT(lmask)
412	ENDIF
413	IF (ln_spmfb) THEN
414	lmask(:) = .FALSE.
415	WHERE (spmfbfiles(:) /= '') lmask(:) = .TRUE.
416	jnumspmfb = COUNT(lmask)
417	ENDIF
418	IF (ln_fco2) THEN
419	lmask(:) = .FALSE.
420	WHERE (fco2files(:) /= '') lmask(:) = .TRUE.
421	jnumfco2 = COUNT(lmask)
422	ENDIF
423	IF (ln_fco2fb) THEN
424	lmask(:) = .FALSE.
425	WHERE (fco2fbfiles(:) /= '') lmask(:) = .TRUE.
426	jnumfco2fb = COUNT(lmask)
427	ENDIF
428	IF (ln_pco2) THEN
429	lmask(:) = .FALSE.
430	WHERE (pco2files(:) /= '') lmask(:) = .TRUE.
431	jnumpco2 = COUNT(lmask)
432	ENDIF
433	IF (ln_pco2fb) THEN
434	lmask(:) = .FALSE.
435	WHERE (pco2fbfiles(:) /= '') lmask(:) = .TRUE.
436	jnumpco2fb = COUNT(lmask)
437	ENDIF
438
439	! Control print
440	IF(lwp) THEN
441	WRITE(numout,*)
442	WRITE(numout,*) 'dia_obs_init : Observation diagnostic initialization'
443	WRITE(numout,*) '~~~~~~~~~~~~'
444	WRITE(numout,*) ' Namelist namobs : set observation diagnostic parameters'
445	WRITE(numout,*) ' Logical switch for T profile observations ln_t3d = ', ln_t3d
446	WRITE(numout,*) ' Logical switch for S profile observations ln_s3d = ', ln_s3d
447	WRITE(numout,*) ' Logical switch for ENACT insitu data set ln_ena = ', ln_ena
448	WRITE(numout,*) ' Logical switch for Coriolis insitu data set ln_cor = ', ln_cor
449	WRITE(numout,*) ' Logical switch for feedback insitu data set ln_profb = ', ln_profb
450	WRITE(numout,*) ' Logical switch for SLA observations ln_sla = ', ln_sla
451	WRITE(numout,*) ' Logical switch for AVISO SLA data ln_sladt = ', ln_sladt
452	WRITE(numout,*) ' Logical switch for feedback SLA data ln_slafb = ', ln_slafb
453	WRITE(numout,*) ' Logical switch for SSH observations ln_ssh = ', ln_ssh
454	WRITE(numout,*) ' Logical switch for SST observations ln_sst = ', ln_sst
455	WRITE(numout,*) ' Logical switch for Reynolds observations ln_reysst = ', ln_reysst
456	WRITE(numout,*) ' Logical switch for GHRSST observations ln_ghrsst = ', ln_ghrsst
457	WRITE(numout,*) ' Logical switch for feedback SST data ln_sstfb = ', ln_sstfb
458	WRITE(numout,*) ' Logical switch for SST bias correction ln_sstbias = ', ln_sstbias
459	WRITE(numout,*) ' Logical switch for night-time SST obs ln_sstnight = ', ln_sstnight
460	WRITE(numout,*) ' Logical switch for SSS observations ln_sss = ', ln_sss
461	WRITE(numout,*) ' Logical switch for Sea Ice observations ln_seaice = ', ln_seaice
462	WRITE(numout,*) ' Logical switch for velocity observations ln_vel3d = ', ln_vel3d
463	WRITE(numout,*) ' Logical switch for velocity daily av. cur. ln_velavcur = ', ln_velavcur
464	WRITE(numout,*) ' Logical switch for velocity high freq. cur. ln_velhrcur = ', ln_velhrcur
465	WRITE(numout,*) ' Logical switch for velocity daily av. ADCP ln_velavadcp = ', ln_velavadcp
466	WRITE(numout,*) ' Logical switch for velocity high freq. ADCP ln_velhradcp = ', ln_velhradcp
467	WRITE(numout,*) ' Logical switch for feedback velocity data ln_velfb = ', ln_velfb
468	WRITE(numout,*) ' Logical switch for logchl observations ln_logchl = ', ln_logchl
469	WRITE(numout,*) ' Logical switch for feedback logchl data ln_logchlfb = ', ln_logchlfb
470	WRITE(numout,*) ' Logical switch for spm observations ln_spm = ', ln_spm
471	WRITE(numout,*) ' Logical switch for feedback spm data ln_spmfb = ', ln_spmfb
472	WRITE(numout,*) ' Logical switch for fco2 observations ln_fco2 = ', ln_fco2
473	WRITE(numout,*) ' Logical switch for pco2 observations ln_pco2 = ', ln_pco2
474	WRITE(numout,*) ' Logical switch for feedback pco2 data ln_pco2fb = ', ln_pco2fb
475	WRITE(numout,*) ' Logical switch for feedback fco2 data ln_fco2fb = ', ln_fco2fb
476	WRITE(numout,*) ' Global distribtion of observations ln_grid_global = ',ln_grid_global
477	WRITE(numout,*) &
478	' Logical switch for obs grid search w/lookup table ln_grid_search_lookup = ',ln_grid_search_lookup
479	IF (ln_grid_search_lookup) &
480	WRITE(numout,*) ' Grid search lookup file header grid_search_file = ', grid_search_file
481	IF (ln_ena) THEN
482	DO ji = 1, jnumenact
483	WRITE(numout,'(1X,2A)') ' ENACT input observation file name enactfiles = ', &
484	TRIM(enactfiles(ji))
485	END DO
486	ENDIF
487	IF (ln_cor) THEN
488	DO ji = 1, jnumcorio
489	WRITE(numout,'(1X,2A)') ' Coriolis input observation file name coriofiles = ', &
490	TRIM(coriofiles(ji))
491	END DO
492	ENDIF
493	IF (ln_profb) THEN
494	DO ji = 1, jnumprofb
495	IF (ln_profb_ena(ji)) THEN
496	WRITE(numout,'(1X,2A)') ' Enact feedback input observation file name profbfiles = ', &
497	TRIM(profbfiles(ji))
498	ELSE
499	WRITE(numout,'(1X,2A)') ' Feedback input observation file name profbfiles = ', &
500	TRIM(profbfiles(ji))
501	ENDIF
502	WRITE(numout,'(1X,2A)') ' Enact feedback input time setting switch ln_profb_enatim = ', ln_profb_enatim(ji)
503	END DO
504	ENDIF
505	IF (ln_sladt) THEN
506	DO ji = 1, jnumslaact
507	WRITE(numout,'(1X,2A)') ' Active SLA input observation file name slafilesact = ', &
508	TRIM(slafilesact(ji))
509	END DO
510	DO ji = 1, jnumslapas
511	WRITE(numout,'(1X,2A)') ' Passive SLA input observation file name slafilespas = ', &
512	TRIM(slafilespas(ji))
513	END DO
514	ENDIF
515	IF (ln_slafb) THEN
516	DO ji = 1, jnumslafb
517	WRITE(numout,'(1X,2A)') ' Feedback SLA input observation file name slafbfiles = ', &
518	TRIM(slafbfiles(ji))
519	END DO
520	ENDIF
521	IF (ln_ghrsst) THEN
522	DO ji = 1, jnumsst
523	WRITE(numout,'(1X,2A)') ' GHRSST input observation file name sstfiles = ', &
524	TRIM(sstfiles(ji))
525	END DO
526	ENDIF
527	IF (ln_sstfb) THEN
528	DO ji = 1, jnumsstfb
529	WRITE(numout,'(1X,2A)') ' Feedback SST input observation file name sstfbfiles = ', &
530	TRIM(sstfbfiles(ji))
531	END DO
532	ENDIF
533	IF (ln_seaice) THEN
534	DO ji = 1, jnumseaice
535	WRITE(numout,'(1X,2A)') ' Sea Ice input observation file name seaicefiles = ', &
536	TRIM(seaicefiles(ji))
537	END DO
538	ENDIF
539	IF (ln_velavcur) THEN
540	DO ji = 1, jnumvelavcur
541	WRITE(numout,'(1X,2A)') ' Vel. cur. daily av. input file name velavcurfiles = ', &
542	TRIM(velavcurfiles(ji))
543	END DO
544	ENDIF
545	IF (ln_velhrcur) THEN
546	DO ji = 1, jnumvelhrcur
547	WRITE(numout,'(1X,2A)') ' Vel. cur. high freq. input file name velhvcurfiles = ', &
548	TRIM(velhrcurfiles(ji))
549	END DO
550	ENDIF
551	IF (ln_velavadcp) THEN
552	DO ji = 1, jnumvelavadcp
553	WRITE(numout,'(1X,2A)') ' Vel. ADCP daily av. input file name velavadcpfiles = ', &
554	TRIM(velavadcpfiles(ji))
555	END DO
556	ENDIF
557	IF (ln_velhradcp) THEN
558	DO ji = 1, jnumvelhradcp
559	WRITE(numout,'(1X,2A)') ' Vel. ADCP high freq. input file name velhvadcpfiles = ', &
560	TRIM(velhradcpfiles(ji))
561	END DO
562	ENDIF
563	IF (ln_velfb) THEN
564	DO ji = 1, jnumvelfb
565	IF (ln_velfb_av(ji)) THEN
566	WRITE(numout,'(1X,2A)') ' Vel. feedback daily av. input file name velfbfiles = ', &
567	TRIM(velfbfiles(ji))
568	ELSE
569	WRITE(numout,'(1X,2A)') ' Vel. feedback input observation file name velfbfiles = ', &
570	TRIM(velfbfiles(ji))
571	ENDIF
572	END DO
573	ENDIF
574	IF (ln_logchl) THEN
575	DO ji = 1, jnumlogchl
576	WRITE(numout,'(1X,2A)') ' logchl input observation file name logchlfiles = ', &
577	TRIM(logchlfiles(ji))
578	END DO
579	ENDIF
580	IF (ln_logchlfb) THEN
581	DO ji = 1, jnumlogchlfb
582	WRITE(numout,'(1X,2A)') ' Feedback logchl input observation file name logchlfbfiles = ', &
583	TRIM(logchlfbfiles(ji))
584	END DO
585	ENDIF
586	IF (ln_spm) THEN
587	DO ji = 1, jnumspm
588	WRITE(numout,'(1X,2A)') ' spm input observation file name spmfiles = ', &
589	TRIM(spmfiles(ji))
590	END DO
591	ENDIF
592	IF (ln_spmfb) THEN
593	DO ji = 1, jnumspmfb
594	WRITE(numout,'(1X,2A)') ' Feedback spm input observation file name spmfbfiles = ', &
595	TRIM(spmfbfiles(ji))
596	END DO
597	ENDIF
598	IF (ln_fco2) THEN
599	DO ji = 1, jnumfco2
600	WRITE(numout,'(1X,2A)') ' fco2 input observation file name fco2files = ', &
601	TRIM(fco2files(ji))
602	END DO
603	ENDIF
604	IF (ln_fco2fb) THEN
605	DO ji = 1, jnumfco2fb
606	WRITE(numout,'(1X,2A)') ' Feedback fco2 input observation file name fco2fbfiles = ', &
607	TRIM(fco2fbfiles(ji))
608	END DO
609	ENDIF
610	IF (ln_pco2) THEN
611	DO ji = 1, jnumpco2
612	WRITE(numout,'(1X,2A)') ' pco2 input observation file name pco2files = ', &
613	TRIM(pco2files(ji))
614	END DO
615	ENDIF
616	IF (ln_pco2fb) THEN
617	DO ji = 1, jnumpco2fb
618	WRITE(numout,'(1X,2A)') ' Feedback pco2 input observation file name pco2fbfiles = ', &
619	TRIM(pco2fbfiles(ji))
620	END DO
621	ENDIF
622	WRITE(numout,*) ' Initial date in window YYYYMMDD.HHMMSS dobsini = ', dobsini
623	WRITE(numout,*) ' Final date in window YYYYMMDD.HHMMSS dobsend = ', dobsend
624	WRITE(numout,*) ' Type of vertical interpolation method n1dint = ', n1dint
625	WRITE(numout,*) ' Type of horizontal interpolation method n2dint = ', n2dint
626	WRITE(numout,*) ' Rejection of observations near land swithch ln_nea = ', ln_nea
627	WRITE(numout,*) ' Rejection of obs near open bdys ln_bound_reject = ', ln_bound_reject
628	WRITE(numout,*) ' MSSH correction scheme nmsshc = ', nmsshc
629	WRITE(numout,*) ' MDT correction mdtcorr = ', mdtcorr
630	WRITE(numout,*) ' MDT cutoff for computed correction mdtcutoff = ', mdtcutoff
631	WRITE(numout,*) ' Logical switch for alt bias ln_altbias = ', ln_altbias
632	WRITE(numout,*) ' Logical switch for ignoring missing files ln_ignmis = ', ln_ignmis
633	WRITE(numout,*) ' ENACT daily average types = ',endailyavtypes
634
635	ENDIF
636
637	IF ( ln_vel3d .AND. ( .NOT. ln_grid_global ) ) THEN
638	CALL ctl_stop( 'Velocity data only works with ln_grid_global=.true.' )
639	RETURN
640	ENDIF
641
642	IF ( ln_grid_global ) THEN
643	CALL ctl_warn( 'ln_grid_global=T may cause memory issues when used with a large number of processors' )
644	ENDIF
645
646	CALL obs_typ_init
647
648	IF ( ln_grid_global ) THEN
649	CALL mppmap_init
650	ENDIF
651
652	! Parameter control
653	#if defined key_diaobs
654	IF ( ( .NOT. ln_t3d ).AND.( .NOT. ln_s3d ).AND.( .NOT. ln_sla ).AND. &
655	& ( .NOT. ln_vel3d ).AND. &
656	& ( .NOT. ln_ssh ).AND.( .NOT. ln_sst ).AND.( .NOT. ln_sss ).AND. &
657	& ( .NOT. ln_seaice ).AND.( .NOT. ln_vel3d ).AND.( .NOT. ln_logchl ).AND. &
658	& ( .NOT. ln_spm ).AND.( .NOT. ln_fco2 ).AND.( .NOT. ln_pco2 ) ) THEN
659	IF(lwp) WRITE(numout,cform_war)
660	IF(lwp) WRITE(numout,*) ' key_diaobs is activated but logical flags', &
661	& ' ln_t3d, ln_s3d, ln_sla, ln_ssh, ln_sst, ln_sss, ln_seaice, ln_vel3d,', &
662	& ' ln_logchl, ln_spm, ln_fco2, ln_pco2 are all set to .FALSE.'
663	nwarn = nwarn + 1
664	ENDIF
665	#endif
666
667	CALL obs_grid_setup( )
668	IF ( ( n1dint < 0 ).OR.( n1dint > 1 ) ) THEN
669	CALL ctl_stop(' Choice of vertical (1D) interpolation method', &
670	& ' is not available')
671	ENDIF
672	IF ( ( n2dint < 0 ).OR.( n2dint > 4 ) ) THEN
673	CALL ctl_stop(' Choice of horizontal (2D) interpolation method', &
674	& ' is not available')
675	ENDIF
676
677	!-----------------------------------------------------------------------
678	! Depending on switches read the various observation types
679	!-----------------------------------------------------------------------
680	! - Temperature/salinity profiles
681
682	IF ( ln_t3d .OR. ln_s3d ) THEN
683
684	! Set the number of variables for profiles to 2 (T and S)
685	nprofvars = 2
686	! Set the number of extra variables for profiles to 1 (insitu temp).
687	nprofextr = 1
688
689	! Count how may insitu data sets we have and allocate data.
690	jprofset = 0
691	IF ( ln_ena ) jprofset = jprofset + 1
692	IF ( ln_cor ) jprofset = jprofset + 1
693	IF ( ln_profb ) jprofset = jprofset + jnumprofb
694	nprofsets = jprofset
695	IF ( nprofsets > 0 ) THEN
696	ALLOCATE(ld_enact(nprofsets))
697	ALLOCATE(profdata(nprofsets))
698	ALLOCATE(prodatqc(nprofsets))
699	ENDIF
700
701	jprofset = 0
702
703	! ENACT insitu data
704
705	IF ( ln_ena ) THEN
706
707	jprofset = jprofset + 1
708
709	ld_enact(jprofset) = .TRUE.
710
711	CALL obs_rea_pro_dri( 1, profdata(jprofset), &
712	& jnumenact, enactfiles(1:jnumenact), &
713	& nprofvars, nprofextr, &
714	& nitend-nit000+2, &
715	& dobsini, dobsend, ln_t3d, ln_s3d, &
716	& ln_ignmis, ln_s_at_t, .TRUE., .FALSE., &
717	& kdailyavtypes = endailyavtypes )
718
719	DO jvar = 1, 2
720
721	CALL obs_prof_staend( profdata(jprofset), jvar )
722
723	END DO
724
725	CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset), &
726	& ln_t3d, ln_s3d, ln_nea, &
727	& kdailyavtypes=endailyavtypes )
728
729	ENDIF
730
731	! Coriolis insitu data
732
733	IF ( ln_cor ) THEN
734
735	jprofset = jprofset + 1
736
737	ld_enact(jprofset) = .FALSE.
738
739	CALL obs_rea_pro_dri( 2, profdata(jprofset), &
740	& jnumcorio, coriofiles(1:jnumcorio), &
741	& nprofvars, nprofextr, &
742	& nitend-nit000+2, &
743	& dobsini, dobsend, ln_t3d, ln_s3d, &
744	& ln_ignmis, ln_s_at_t, .FALSE., .FALSE. )
745
746	DO jvar = 1, 2
747
748	CALL obs_prof_staend( profdata(jprofset), jvar )
749
750	END DO
751
752	CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset), &
753	& ln_t3d, ln_s3d, ln_nea )
754
755	ENDIF
756
757	! Feedback insitu data
758
759	IF ( ln_profb ) THEN
760
761	DO jset = 1, jnumprofb
762
763	jprofset = jprofset + 1
764	ld_enact (jprofset) = ln_profb_ena(jset)
765
766	CALL obs_rea_pro_dri( 0, profdata(jprofset), &
767	& 1, profbfiles(jset:jset), &
768	& nprofvars, nprofextr, &
769	& nitend-nit000+2, &
770	& dobsini, dobsend, ln_t3d, ln_s3d, &
771	& ln_ignmis, ln_s_at_t, &
772	& ld_enact(jprofset).AND.&
773	& ln_profb_enatim(jset), &
774	& .FALSE., kdailyavtypes = endailyavtypes )
775
776	DO jvar = 1, 2
777
778	CALL obs_prof_staend( profdata(jprofset), jvar )
779
780	END DO
781
782	IF ( ld_enact(jprofset) ) THEN
783	CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset), &
784	& ln_t3d, ln_s3d, ln_nea, &
785	& kdailyavtypes = endailyavtypes )
786	ELSE
787	CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset), &
788	& ln_t3d, ln_s3d, ln_nea )
789	ENDIF
790
791	END DO
792
793	ENDIF
794
795	ENDIF
796
797	! - Sea level anomalies
798	IF ( ln_sla ) THEN
799	! Set the number of variables for sla to 1
800	nslavars = 1
801
802	! Set the number of extra variables for sla to 2
803	nslaextr = 2
804
805	! Set the number of sla data sets to 2
806	nslasets = 0
807	IF ( ln_sladt ) THEN
808	nslasets = nslasets + 2
809	ENDIF
810	IF ( ln_slafb ) THEN
811	nslasets = nslasets + jnumslafb
812	ENDIF
813
814	ALLOCATE(sladata(nslasets))
815	ALLOCATE(sladatqc(nslasets))
816	sladata(:)%nsurf=0
817	sladatqc(:)%nsurf=0
818
819	nslasets = 0
820
821	! AVISO SLA data
822
823	IF ( ln_sladt ) THEN
824
825	! Active SLA observations
826
827	nslasets = nslasets + 1
828
829	CALL obs_rea_sla( 1, sladata(nslasets), jnumslaact, &
830	& slafilesact(1:jnumslaact), &
831	& nslavars, nslaextr, nitend-nit000+2, &
832	& dobsini, dobsend, ln_ignmis, .FALSE. )
833	CALL obs_pre_sla( sladata(nslasets), sladatqc(nslasets), &
834	& ln_sla, ln_nea )
835
836	! Passive SLA observations
837
838	nslasets = nslasets + 1
839
840	CALL obs_rea_sla( 1, sladata(nslasets), jnumslapas, &
841	& slafilespas(1:jnumslapas), &
842	& nslavars, nslaextr, nitend-nit000+2, &
843	& dobsini, dobsend, ln_ignmis, .FALSE. )
844
845	CALL obs_pre_sla( sladata(nslasets), sladatqc(nslasets), &
846	& ln_sla, ln_nea )
847
848	ENDIF
849
850	! Feedback SLA data
851
852	IF ( ln_slafb ) THEN
853
854	DO jset = 1, jnumslafb
855
856	nslasets = nslasets + 1
857
858	CALL obs_rea_sla( 0, sladata(nslasets), 1, &
859	& slafbfiles(jset:jset), &
860	& nslavars, nslaextr, nitend-nit000+2, &
861	& dobsini, dobsend, ln_ignmis, .FALSE. )
862	CALL obs_pre_sla( sladata(nslasets), sladatqc(nslasets), &
863	& ln_sla, ln_nea )
864
865	END DO
866
867	ENDIF
868
869	CALL obs_rea_mdt( nslasets, sladatqc, n2dint )
870
871	! read in altimeter bias
872
873	IF ( ln_altbias ) THEN
874	CALL obs_rea_altbias ( nslasets, sladatqc, n2dint, bias_file )
875	ENDIF
876
877	ENDIF
878
879	! - Sea surface height
880	IF ( ln_ssh ) THEN
881	IF(lwp) WRITE(numout,*) ' SSH currently not available'
882	ENDIF
883
884	! - Sea surface temperature
885	IF ( ln_sst ) THEN
886
887	! Set the number of variables for sst to 1
888	nsstvars = 1
889
890	! Set the number of extra variables for sst to 0
891	nsstextr = 0
892
893	nsstsets = 0
894
895	IF (ln_reysst) nsstsets = nsstsets + 1
896	IF (ln_ghrsst) nsstsets = nsstsets + 1
897	IF ( ln_sstfb ) THEN
898	nsstsets = nsstsets + jnumsstfb
899	ENDIF
900
901	ALLOCATE(sstdata(nsstsets))
902	ALLOCATE(sstdatqc(nsstsets))
903	ALLOCATE(ld_sstnight(nsstsets))
904	sstdata(:)%nsurf=0
905	sstdatqc(:)%nsurf=0
906	ld_sstnight(:)=.false.
907
908	nsstsets = 0
909
910	IF (ln_reysst) THEN
911
912	nsstsets = nsstsets + 1
913
914	ld_sstnight(nsstsets) = ln_sstnight
915
916	CALL obs_rea_sst_rey( reysstname, reysstfmt, sstdata(nsstsets), &
917	& nsstvars, nsstextr, &
918	& nitend-nit000+2, dobsini, dobsend )
919	CALL obs_pre_sst( sstdata(nsstsets), sstdatqc(nsstsets), ln_sst, &
920	& ln_nea )
921
922	ENDIF
923
924	IF (ln_ghrsst) THEN
925
926	nsstsets = nsstsets + 1
927
928	ld_sstnight(nsstsets) = ln_sstnight
929
930	CALL obs_rea_sst( 1, sstdata(nsstsets), jnumsst, &
931	& sstfiles(1:jnumsst), &
932	& nsstvars, nsstextr, nitend-nit000+2, &
933	& dobsini, dobsend, ln_ignmis, .FALSE. )
934	CALL obs_pre_sst( sstdata(nsstsets), sstdatqc(nsstsets), ln_sst, &
935	& ln_nea )
936
937	ENDIF
938
939	! Feedback SST data
940
941	IF ( ln_sstfb ) THEN
942
943	DO jset = 1, jnumsstfb
944
945	nsstsets = nsstsets + 1
946
947	ld_sstnight(nsstsets) = ln_sstnight
948
949	CALL obs_rea_sst( 0, sstdata(nsstsets), 1, &
950	& sstfbfiles(jset:jset), &
951	& nsstvars, nsstextr, nitend-nit000+2, &
952	& dobsini, dobsend, ln_ignmis, .FALSE. )
953	CALL obs_pre_sst( sstdata(nsstsets), sstdatqc(nsstsets), &
954	& ln_sst, ln_nea )
955
956	END DO
957
958	ENDIF
959
960	!Read in bias field and correct SST.
961	IF ( ln_sstbias ) THEN
962	IF ( jnumsstbias == 0 ) CALL ctl_stop("ln_sstbias set,"// &
963	" but no bias"// &
964	" files to read in")
965	CALL obs_app_sstbias( nsstsets, sstdatqc, n2dint, &
966	jnumsstbias, &
967	sstbias_files(1:jnumsstbias) )
968	ENDIF
969
970	ENDIF
971
972	! - Sea surface salinity
973	IF ( ln_sss ) THEN
974	IF(lwp) WRITE(numout,*) ' SSS currently not available'
975	ENDIF
976
977	! - Sea Ice Concentration
978
979	IF ( ln_seaice ) THEN
980
981	! Set the number of variables for seaice to 1
982	nseaicevars = 1
983
984	! Set the number of extra variables for seaice to 0
985	nseaiceextr = 0
986
987	! Set the number of data sets to 1
988	nseaicesets = 1
989
990	ALLOCATE(seaicedata(nseaicesets))
991	ALLOCATE(seaicedatqc(nseaicesets))
992	seaicedata(:)%nsurf=0
993	seaicedatqc(:)%nsurf=0
994
995	CALL obs_rea_seaice( 1, seaicedata(nseaicesets), jnumseaice, &
996	& seaicefiles(1:jnumseaice), &
997	& nseaicevars, nseaiceextr, nitend-nit000+2, &
998	& dobsini, dobsend, ln_ignmis, .FALSE. )
999
1000	CALL obs_pre_seaice( seaicedata(nseaicesets), seaicedatqc(nseaicesets), &
1001	& ln_seaice, ln_nea )
1002
1003	ENDIF
1004
1005	IF (ln_vel3d) THEN
1006
1007	! Set the number of variables for profiles to 2 (U and V)
1008	nvelovars = 2
1009
1010	! Set the number of extra variables for profiles to 2 to store
1011	! rotation parameters
1012	nveloextr = 2
1013
1014	jveloset = 0
1015
1016	IF ( ln_velavcur ) jveloset = jveloset + 1
1017	IF ( ln_velhrcur ) jveloset = jveloset + 1
1018	IF ( ln_velavadcp ) jveloset = jveloset + 1
1019	IF ( ln_velhradcp ) jveloset = jveloset + 1
1020	IF (ln_velfb) jveloset = jveloset + jnumvelfb
1021
1022	nvelosets = jveloset
1023	IF ( nvelosets > 0 ) THEN
1024	ALLOCATE( velodata(nvelosets) )
1025	ALLOCATE( veldatqc(nvelosets) )
1026	ALLOCATE( ld_velav(nvelosets) )
1027	ENDIF
1028
1029	jveloset = 0
1030
1031	! Daily averaged data
1032
1033	IF ( ln_velavcur ) THEN
1034
1035	jveloset = jveloset + 1
1036
1037	ld_velav(jveloset) = .TRUE.
1038
1039	CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelavcur, &
1040	& velavcurfiles(1:jnumvelavcur), &
1041	& nvelovars, nveloextr, &
1042	& nitend-nit000+2, &
1043	& dobsini, dobsend, ln_ignmis, &
1044	& ld_velav(jveloset), &
1045	& .FALSE. )
1046
1047	DO jvar = 1, 2
1048	CALL obs_prof_staend( velodata(jveloset), jvar )
1049	END DO
1050
1051	CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), &
1052	& ln_vel3d, ln_nea, ld_velav(jveloset) )
1053
1054	ENDIF
1055
1056	! High frequency data
1057
1058	IF ( ln_velhrcur ) THEN
1059
1060	jveloset = jveloset + 1
1061
1062	ld_velav(jveloset) = .FALSE.
1063
1064	CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelhrcur, &
1065	& velhrcurfiles(1:jnumvelhrcur), &
1066	& nvelovars, nveloextr, &
1067	& nitend-nit000+2, &
1068	& dobsini, dobsend, ln_ignmis, &
1069	& ld_velav(jveloset), &
1070	& .FALSE. )
1071
1072	DO jvar = 1, 2
1073	CALL obs_prof_staend( velodata(jveloset), jvar )
1074	END DO
1075
1076	CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), &
1077	& ln_vel3d, ln_nea, ld_velav(jveloset) )
1078
1079	ENDIF
1080
1081	! Daily averaged data
1082
1083	IF ( ln_velavadcp ) THEN
1084
1085	jveloset = jveloset + 1
1086
1087	ld_velav(jveloset) = .TRUE.
1088
1089	CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelavadcp, &
1090	& velavadcpfiles(1:jnumvelavadcp), &
1091	& nvelovars, nveloextr, &
1092	& nitend-nit000+2, &
1093	& dobsini, dobsend, ln_ignmis, &
1094	& ld_velav(jveloset), &
1095	& .FALSE. )
1096
1097	DO jvar = 1, 2
1098	CALL obs_prof_staend( velodata(jveloset), jvar )
1099	END DO
1100
1101	CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), &
1102	& ln_vel3d, ln_nea, ld_velav(jveloset) )
1103
1104	ENDIF
1105
1106	! High frequency data
1107
1108	IF ( ln_velhradcp ) THEN
1109
1110	jveloset = jveloset + 1
1111
1112	ld_velav(jveloset) = .FALSE.
1113
1114	CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelhradcp, &
1115	& velhradcpfiles(1:jnumvelhradcp), &
1116	& nvelovars, nveloextr, &
1117	& nitend-nit000+2, &
1118	& dobsini, dobsend, ln_ignmis, &
1119	& ld_velav(jveloset), &
1120	& .FALSE. )
1121
1122	DO jvar = 1, 2
1123	CALL obs_prof_staend( velodata(jveloset), jvar )
1124	END DO
1125
1126	CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), &
1127	& ln_vel3d, ln_nea, ld_velav(jveloset) )
1128
1129	ENDIF
1130
1131	IF ( ln_velfb ) THEN
1132
1133	DO jset = 1, jnumvelfb
1134
1135	jveloset = jveloset + 1
1136
1137	ld_velav(jveloset) = ln_velfb_av(jset)
1138
1139	CALL obs_rea_vel_dri( 0, velodata(jveloset), 1, &
1140	& velfbfiles(jset:jset), &
1141	& nvelovars, nveloextr, &
1142	& nitend-nit000+2, &
1143	& dobsini, dobsend, ln_ignmis, &
1144	& ld_velav(jveloset), &
1145	& .FALSE. )
1146
1147	DO jvar = 1, 2
1148	CALL obs_prof_staend( velodata(jveloset), jvar )
1149	END DO
1150
1151	CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), &
1152	& ln_vel3d, ln_nea, ld_velav(jveloset) )
1153
1154
1155	END DO
1156
1157	ENDIF
1158
1159	ENDIF
1160
1161	! - log10(chlorophyll)
1162
1163	IF ( ln_logchl ) THEN
1164
1165	! Set the number of variables for logchl to 1
1166	nlogchlvars = 1
1167
1168	! Set the number of extra variables for logchl to 1 (STD)
1169	nlogchlextr = 1
1170
1171	IF ( ln_logchlfb ) THEN
1172	nlogchlsets = jnumlogchlfb
1173	ELSE
1174	nlogchlsets = 1
1175	ENDIF
1176
1177	ALLOCATE(logchldata(nlogchlsets))
1178	ALLOCATE(logchldatqc(nlogchlsets))
1179	logchldata(:)%nsurf=0
1180	logchldatqc(:)%nsurf=0
1181
1182	nlogchlsets = 0
1183
1184	IF ( ln_logchlfb ) THEN ! Feedback file format
1185
1186	DO jset = 1, jnumlogchlfb
1187
1188	nlogchlsets = nlogchlsets + 1
1189
1190	CALL obs_rea_logchl( 0, logchldata(nlogchlsets), 1, &
1191	& logchlfbfiles(jset:jset), &
1192	& nlogchlvars, nlogchlextr, nitend-nit000+2, &
1193	& dobsini, dobsend, ln_ignmis, .FALSE. )
1194
1195	CALL obs_pre_logchl( logchldata(nlogchlsets), logchldatqc(nlogchlsets), &
1196	& ln_logchl, ln_nea )
1197
1198	ENDDO
1199
1200	ELSE ! Original file format
1201
1202	nlogchlsets = nlogchlsets + 1
1203
1204	CALL obs_rea_logchl( 1, logchldata(nlogchlsets), jnumlogchl, &
1205	& logchlfiles(1:jnumlogchl), &
1206	& nlogchlvars, nlogchlextr, nitend-nit000+2, &
1207	& dobsini, dobsend, ln_ignmis, .FALSE. )
1208
1209	CALL obs_pre_logchl( logchldata(nlogchlsets), logchldatqc(nlogchlsets), &
1210	& ln_logchl, ln_nea )
1211
1212	ENDIF
1213
1214	ENDIF
1215
1216	! - spm
1217
1218	IF ( ln_spm ) THEN
1219
1220	! Set the number of variables for spm to 1
1221	nspmvars = 1
1222
1223	! Set the number of extra variables for spm to 0
1224	nspmextr = 0
1225
1226	IF ( ln_spmfb ) THEN
1227	nspmsets = jnumspmfb
1228	ELSE
1229	nspmsets = 1
1230	ENDIF
1231
1232	ALLOCATE(spmdata(nspmsets))
1233	ALLOCATE(spmdatqc(nspmsets))
1234	spmdata(:)%nsurf=0
1235	spmdatqc(:)%nsurf=0
1236
1237	nspmsets = 0
1238
1239	IF ( ln_spmfb ) THEN ! Feedback file format
1240
1241	DO jset = 1, jnumspmfb
1242
1243	nspmsets = nspmsets + 1
1244
1245	CALL obs_rea_spm( 0, spmdata(nspmsets), 1, &
1246	& spmfbfiles(jset:jset), &
1247	& nspmvars, nspmextr, nitend-nit000+2, &
1248	& dobsini, dobsend, ln_ignmis, .FALSE. )
1249
1250	CALL obs_pre_spm( spmdata(nspmsets), spmdatqc(nspmsets), &
1251	& ln_spm, ln_nea )
1252
1253	ENDDO
1254
1255	ELSE ! Original file format
1256
1257	nspmsets = nspmsets + 1
1258
1259	CALL obs_rea_spm( 1, spmdata(nspmsets), jnumspm, &
1260	& spmfiles(1:jnumspm), &
1261	& nspmvars, nspmextr, nitend-nit000+2, &
1262	& dobsini, dobsend, ln_ignmis, .FALSE. )
1263
1264	CALL obs_pre_spm( spmdata(nspmsets), spmdatqc(nspmsets), &
1265	& ln_spm, ln_nea )
1266
1267	ENDIF
1268
1269	ENDIF
1270
1271	! - fco2
1272
1273	IF ( ln_fco2 ) THEN
1274
1275	! Set the number of variables for fco2 to 1
1276	nfco2vars = 1
1277
1278	! Set the number of extra variables for fco2 to 0
1279	nfco2extr = 0
1280
1281	IF ( ln_fco2fb ) THEN
1282	nfco2sets = jnumfco2fb
1283	ELSE
1284	nfco2sets = 1
1285	ENDIF
1286
1287	ALLOCATE(fco2data(nfco2sets))
1288	ALLOCATE(fco2datqc(nfco2sets))
1289	fco2data(:)%nsurf=0
1290	fco2datqc(:)%nsurf=0
1291
1292	nfco2sets = 0
1293
1294	IF ( ln_fco2fb ) THEN ! Feedback file format
1295
1296	DO jset = 1, jnumfco2fb
1297
1298	nfco2sets = nfco2sets + 1
1299
1300	CALL obs_rea_fco2( 0, fco2data(nfco2sets), 1, &
1301	& fco2fbfiles(jset:jset), &
1302	& nfco2vars, nfco2extr, nitend-nit000+2, &
1303	& dobsini, dobsend, ln_ignmis, .FALSE. )
1304
1305	CALL obs_pre_fco2( fco2data(nfco2sets), fco2datqc(nfco2sets), &
1306	& ln_fco2, ln_nea )
1307
1308	ENDDO
1309
1310	ELSE ! Original file format
1311
1312	nfco2sets = nfco2sets + 1
1313
1314	CALL obs_rea_fco2( 1, fco2data(nfco2sets), jnumfco2, &
1315	& fco2files(1:jnumfco2), &
1316	& nfco2vars, nfco2extr, nitend-nit000+2, &
1317	& dobsini, dobsend, ln_ignmis, .FALSE. )
1318
1319	CALL obs_pre_fco2( fco2data(nfco2sets), fco2datqc(nfco2sets), &
1320	& ln_fco2, ln_nea )
1321
1322	ENDIF
1323
1324	ENDIF
1325
1326	! - pco2
1327
1328	IF ( ln_pco2 ) THEN
1329
1330	! Set the number of variables for pco2 to 1
1331	npco2vars = 1
1332
1333	! Set the number of extra variables for pco2 to 0
1334	npco2extr = 0
1335
1336	IF ( ln_pco2fb ) THEN
1337	npco2sets = jnumpco2fb
1338	ELSE
1339	npco2sets = 1
1340	ENDIF
1341
1342	ALLOCATE(pco2data(npco2sets))
1343	ALLOCATE(pco2datqc(npco2sets))
1344	pco2data(:)%nsurf=0
1345	pco2datqc(:)%nsurf=0
1346
1347	npco2sets = 0
1348
1349	IF ( ln_pco2fb ) THEN ! Feedback file format
1350
1351	DO jset = 1, jnumpco2fb
1352
1353	npco2sets = npco2sets + 1
1354
1355	CALL obs_rea_pco2( 0, pco2data(npco2sets), 1, &
1356	& pco2fbfiles(jset:jset), &
1357	& npco2vars, npco2extr, nitend-nit000+2, &
1358	& dobsini, dobsend, ln_ignmis, .FALSE. )
1359
1360	CALL obs_pre_pco2( pco2data(npco2sets), pco2datqc(npco2sets), &
1361	& ln_pco2, ln_nea )
1362
1363	ENDDO
1364
1365	ELSE ! Original file format
1366
1367	npco2sets = npco2sets + 1
1368
1369	CALL obs_rea_pco2( 1, pco2data(npco2sets), jnumpco2, &
1370	& pco2files(1:jnumpco2), &
1371	& npco2vars, npco2extr, nitend-nit000+2, &
1372	& dobsini, dobsend, ln_ignmis, .FALSE. )
1373
1374	CALL obs_pre_pco2( pco2data(npco2sets), pco2datqc(npco2sets), &
1375	& ln_pco2, ln_nea )
1376
1377	ENDIF
1378
1379	ENDIF
1380
1381	END SUBROUTINE dia_obs_init
1382
1383	SUBROUTINE dia_obs( kstp )
1384	!!----------------------------------------------------------------------
1385	!! * ROUTINE dia_obs *
1386	!!
1387	!! ** Purpose : Call the observation operators on each time step
1388	!!
1389	!! ** Method : Call the observation operators on each time step to
1390	!! compute the model equivalent of the following date:
1391	!! - T profiles
1392	!! - S profiles
1393	!! - Sea surface height (referenced to a mean)
1394	!! - Sea surface temperature
1395	!! - Sea surface salinity
1396	!! - Velocity component (U,V) profiles
1397	!! - Sea surface log10(chlorophyll)
1398	!! - Sea surface spm
1399	!! - Sea surface fco2
1400	!! - Sea surface pco2
1401	!!
1402	!! ** Action :
1403	!!
1404	!! History :
1405	!! ! 06-03 (K. Mogensen) Original code
1406	!! ! 06-05 (K. Mogensen) Reformatted
1407	!! ! 06-10 (A. Weaver) Cleaning
1408	!! ! 07-03 (K. Mogensen) General handling of profiles
1409	!! ! 07-04 (G. Smith) Generalized surface operators
1410	!! ! 08-10 (M. Valdivieso) obs operator for velocity profiles
1411	!! ! 14-08 (J. While) observation operator for profiles in
1412	!! generalised vertical coordinates
1413	!!----------------------------------------------------------------------
1414	!! * Modules used
1415	USE dom_oce, ONLY : & ! Ocean space and time domain variables
1416	& rdt, &
1417	& gdept_1d, &
1418	#if defined key_vvl
1419	& gdept_n, &
1420	#else
1421	& gdept_1d, &
1422	#endif
1423	& tmask, umask, vmask
1424	USE phycst, ONLY : & ! Physical constants
1425	& rday, &
1426	& rt0
1427	USE oce, ONLY : & ! Ocean dynamics and tracers variables
1428	& tsn, &
1429	& un, vn, &
1430	& sshn
1431	#if defined key_lim3
1432	USE ice, ONLY : & ! LIM Ice model variables
1433	& frld
1434	#endif
1435	#if defined key_lim2
1436	USE ice_2, ONLY : & ! LIM Ice model variables
1437	& frld
1438	#endif
1439	#if defined key_hadocc
1440	USE trc, ONLY : & ! HadOCC chlorophyll, fCO2 and pCO2
1441	& HADOCC_CHL, &
1442	& HADOCC_FCO2, &
1443	& HADOCC_PCO2, &
1444	& HADOCC_FILL_FLT
1445	#elif defined key_medusa && defined key_foam_medusa
1446	USE trc, ONLY : & ! MEDUSA chlorophyll, fCO2 and pCO2
1447	& MEDUSA_CHL, &
1448	& MEDUSA_FCO2, &
1449	& MEDUSA_PCO2, &
1450	& MEDUSA_FILL_FLT
1451	#elif defined key_fabm
1452	USE fabm
1453	USE par_fabm
1454	#endif
1455	#if defined key_spm
1456	USE par_spm, ONLY: & ! ERSEM/SPM sediments
1457	& jp_spm
1458	USE trc, ONLY : &
1459	& trn
1460	#endif
1461	IMPLICIT NONE
1462
1463	!! * Arguments
1464	INTEGER, INTENT(IN) :: kstp ! Current timestep
1465	!! * Local declarations
1466	INTEGER :: idaystp ! Number of timesteps per day
1467	INTEGER :: jprofset ! Profile data set loop variable
1468	INTEGER :: jslaset ! SLA data set loop variable
1469	INTEGER :: jsstset ! SST data set loop variable
1470	INTEGER :: jseaiceset ! sea ice data set loop variable
1471	INTEGER :: jveloset ! velocity profile data loop variable
1472	INTEGER :: jlogchlset ! logchl data set loop variable
1473	INTEGER :: jspmset ! spm data set loop variable
1474	INTEGER :: jfco2set ! fco2 data set loop variable
1475	INTEGER :: jpco2set ! pco2 data set loop variable
1476	INTEGER :: jvar ! Variable number
1477	#if ! defined key_lim2 && ! defined key_lim3
1478	REAL(wp), POINTER, DIMENSION(:,:) :: frld
1479	#endif
1480	REAL(wp) :: tiny ! small number
1481	REAL(wp), DIMENSION(jpi,jpj) :: &
1482	logchl ! array for log chlorophyll
1483	REAL(wp), DIMENSION(jpi,jpj) :: &
1484	maskchl ! array for special chlorophyll mask
1485	REAL(wp), DIMENSION(jpi,jpj) :: &
1486	spm ! array for spm
1487	REAL(wp), DIMENSION(jpi,jpj) :: &
1488	fco2 ! array for fco2
1489	REAL(wp), DIMENSION(jpi,jpj) :: &
1490	maskfco2 ! array for special fco2 mask
1491	REAL(wp), DIMENSION(jpi,jpj) :: &
1492	pco2 ! array for pco2
1493	REAL(wp), DIMENSION(jpi,jpj) :: &
1494	maskpco2 ! array for special pco2 mask
1495	INTEGER :: jn ! loop index
1496	#if defined key_fabm
1497	REAL(wp), DIMENSION(jpi,jpj,jpk) :: logchl_3d
1498	REAL(wp), DIMENSION(jpi,jpj,jpk) :: pco2_3d
1499	#endif
1500	CHARACTER(LEN=20) :: datestr=" ",timestr=" "
1501
1502	#if ! defined key_lim2 && ! defined key_lim3
1503	CALL wrk_alloc(jpi,jpj,frld)
1504	#endif
1505
1506	IF(lwp) THEN
1507	WRITE(numout,*)
1508	WRITE(numout,*) 'dia_obs : Call the observation operators', kstp
1509	WRITE(numout,*) '~~~~~~~'
1510	ENDIF
1511
1512	idaystp = NINT( rday / rdt )
1513
1514	!-----------------------------------------------------------------------
1515	! No LIM => frld == 0.0_wp
1516	!-----------------------------------------------------------------------
1517	#if ! defined key_lim2 && ! defined key_lim3
1518	frld(:,:) = 0.0_wp
1519	#endif
1520	!-----------------------------------------------------------------------
1521	! Depending on switches call various observation operators
1522	!-----------------------------------------------------------------------
1523
1524	! - Temperature/salinity profiles
1525	IF ( ln_t3d .OR. ln_s3d ) THEN
1526	DO jprofset = 1, nprofsets
1527	IF( (.NOT. lk_vvl) .AND. (ln_zco .OR. ln_zps) ) THEN
1528	IF(lwp) THEN
1529	WRITE(numout,*) 'dia_obs : calling obs_pro_opt'
1530	ENDIF
1531	IF ( ld_enact(jprofset) ) THEN
1532	CALL obs_pro_opt( prodatqc(jprofset), &
1533	& kstp, jpi, jpj, jpk, nit000, idaystp, &
1534	& tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal), &
1535	& gdept_1d, tmask, n1dint, n2dint, &
1536	& kdailyavtypes = endailyavtypes )
1537	ELSE
1538	CALL obs_pro_opt( prodatqc(jprofset), &
1539	& kstp, jpi, jpj, jpk, nit000, idaystp, &
1540	& tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal), &
1541	& gdept_1d, tmask, n1dint, n2dint )
1542	ENDIF
1543	ELSE
1544	IF(lwp) THEN
1545	WRITE(numout,*) 'dia_obs : calling obs_pro_sco_opt'
1546	ENDIF
1547	IF ( ld_enact(jprofset) ) THEN
1548	CALL obs_pro_sco_opt( prodatqc(jprofset), &
1549	& kstp, jpi, jpj, jpk, nit000, idaystp, &
1550	& tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal), &
1551	& fsdept(:,:,:), fsdepw(:,:,:), &
1552	& tmask, n1dint, n2dint, &
1553	& kdailyavtypes = endailyavtypes )
1554	ELSE
1555	CALL obs_pro_sco_opt( prodatqc(jprofset), &
1556	& kstp, jpi, jpj, jpk, nit000, idaystp, &
1557	& tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal), &
1558	& fsdept(:,:,:), fsdepw(:,:,:), &
1559	& tmask, n1dint, n2dint )
1560	ENDIF
1561	ENDIF
1562	END DO
1563	ENDIF
1564
1565	! - Sea surface anomaly
1566	IF ( ln_sla ) THEN
1567	DO jslaset = 1, nslasets
1568	CALL obs_sla_opt( sladatqc(jslaset), &
1569	& kstp, jpi, jpj, nit000, sshn, &
1570	& tmask(:,:,1), n2dint )
1571	END DO
1572	ENDIF
1573
1574	! - Sea surface temperature
1575	IF ( ln_sst ) THEN
1576	DO jsstset = 1, nsstsets
1577	CALL obs_sst_opt( sstdatqc(jsstset), &
1578	& kstp, jpi, jpj, nit000, idaystp, &
1579	& tsn(:,:,1,jp_tem), tmask(:,:,1), &
1580	& n2dint, ld_sstnight(jsstset) )
1581	END DO
1582	ENDIF
1583
1584	! - Sea surface salinity
1585	IF ( ln_sss ) THEN
1586	IF(lwp) WRITE(numout,*) ' SSS currently not available'
1587	ENDIF
1588
1589	#if defined key_lim2 \|\| defined key_lim3
1590	IF ( ln_seaice ) THEN
1591	DO jseaiceset = 1, nseaicesets
1592	CALL obs_seaice_opt( seaicedatqc(jseaiceset), &
1593	& kstp, jpi, jpj, nit000, 1.-frld, &
1594	& tmask(:,:,1), n2dint )
1595	END DO
1596	ENDIF
1597	#endif
1598
1599	! - Velocity profiles
1600	IF ( ln_vel3d ) THEN
1601	DO jveloset = 1, nvelosets
1602	! zonal component of velocity
1603	CALL obs_vel_opt( veldatqc(jveloset), kstp, jpi, jpj, jpk, &
1604	& nit000, idaystp, un, vn, gdept_1d, umask, vmask, &
1605	n1dint, n2dint, ld_velav(jveloset) )
1606	END DO
1607	ENDIF
1608
1609	IF ( ln_logchl ) THEN
1610
1611	#if defined key_hadocc
1612	logchl(:,:) = HADOCC_CHL(:,:,1) ! (not log) chlorophyll from HadOCC
1613	#elif defined key_medusa && defined key_foam_medusa
1614	logchl(:,:) = MEDUSA_CHL(:,:,1) ! (not log) chlorophyll from HadOCC
1615	#elif defined key_fabm
1616	logchl_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabmdia_chltot)
1617	logchl(:,:) = logchl_3d(:,:,1)
1618	#else
1619	CALL ctl_stop( ' Trying to run logchl observation operator', &
1620	& ' but no biogeochemical model appears to have been defined' )
1621	#endif
1622
1623	maskchl(:,:) = tmask(:,:,1) ! create a special mask to exclude certain things
1624
1625	! Take the log10 where we can, otherwise exclude
1626	tiny = 1.0e-20
1627	WHERE(logchl(:,:) > tiny .AND. logchl(:,:) /= obfillflt )
1628	logchl(:,:) = LOG10(logchl(:,:))
1629	ELSEWHERE
1630	logchl(:,:) = obfillflt
1631	maskchl(:,:) = 0
1632	END WHERE
1633
1634	DO jlogchlset = 1, nlogchlsets
1635	CALL obs_logchl_opt( logchldatqc(jlogchlset), &
1636	& kstp, jpi, jpj, nit000, logchl(:,:), &
1637	& maskchl(:,:), n2dint )
1638	END DO
1639	ENDIF
1640
1641	IF ( ln_spm ) THEN
1642	#if defined key_spm
1643	spm(:,:) = 0.0
1644	DO jn = 1, jp_spm
1645	spm(:,:) = spm(:,:) + trn(:,:,1,jn) ! sum SPM sizes
1646	END DO
1647	#else
1648	CALL ctl_stop( ' Trying to run spm observation operator', &
1649	& ' but no spm model appears to have been defined' )
1650	#endif
1651
1652	DO jspmset = 1, nspmsets
1653	CALL obs_spm_opt( spmdatqc(jspmset), &
1654	& kstp, jpi, jpj, nit000, spm(:,:), &
1655	& tmask(:,:,1), n2dint )
1656	END DO
1657	ENDIF
1658
1659	IF ( ln_fco2 ) THEN
1660	maskfco2(:,:) = tmask(:,:,1) ! create a special mask to exclude certain things
1661	#if defined key_hadocc
1662	fco2(:,:) = HADOCC_FCO2(:,:) ! fCO2 from HadOCC
1663	IF ( ( MINVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ).AND.( MAXVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) ) THEN
1664	fco2(:,:) = obfillflt
1665	maskfco2(:,:) = 0
1666	CALL ctl_warn( ' HadOCC fCO2 values masked out for observation operator', &
1667	& ' on timestep ' // TRIM(STR(kstp)), &
1668	& ' as HADOCC_FCO2(:,:) == HADOCC_FILL_FLT' )
1669	ENDIF
1670	#elif defined key_medusa && defined key_foam_medusa
1671	fco2(:,:) = MEDUSA_FCO2(:,:) ! fCO2 from MEDUSA
1672	IF ( ( MINVAL( MEDUSA_FCO2 ) == MEDUSA_FILL_FLT ).AND.( MAXVAL( MEDUSA_FCO2 ) == MEDUSA_FILL_FLT ) ) THEN
1673	fco2(:,:) = obfillflt
1674	maskfco2(:,:) = 0
1675	CALL ctl_warn( ' MEDUSA fCO2 values masked out for observation operator', &
1676	& ' on timestep ' // TRIM(STR(kstp)), &
1677	& ' as MEDUSA_FCO2(:,:) == MEDUSA_FILL_FLT' )
1678	ENDIF
1679	#elif defined key_fabm
1680	! First, get pCO2 from FABM
1681	pco2_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabm_o3pc)
1682	pco2(:,:) = pco2_3d(:,:,1)
1683	! Now, convert pCO2 to fCO2, based on SST in K. This follows the standard methodology of:
1684	! Pierrot et al. (2009), Recommendations for autonomous underway pCO2 measuring systems
1685	! and data reduction routines, Deep-Sea Research II, 56: 512-522.
1686	! and
1687	! Weiss (1974), Carbon dioxide in water and seawater: the solubility of a non-ideal gas,
1688	! Marine Chemistry, 2: 203-215.
1689	! In the implementation below, atmospheric pressure has been assumed to be 1 atm and so
1690	! not explicitly included - atmospheric pressure is not necessarily available so this is
1691	! the best assumption.
1692	! Further, the (1-xCO2)^2 term has been neglected. This is common practice
1693	! (see e.g. Zeebe and Wolf-Gladrow (2001), CO2 in Seawater: Equilibrium, Kinetics, Isotopes)
1694	! because xCO2 in atm is ~0, and so this term will only affect the result to the 3rd decimal
1695	! place for typical values, and xCO2 would need to be approximated from pCO2 anyway.
1696	fco2(:,:) = pco2(:,:) * EXP((-1636.75 + &
1697	& 12.0408 * (tsn(:,:,1,jp_tem)+rt0) - &
1698	& 0.0327957 * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0) + &
1699	& 0.0000316528 * (tsn(:,:,1,jp_tem)+rt0)(tsn(:,:,1,jp_tem)+rt0)(tsn(:,:,1,jp_tem)+rt0) + &
1700	& 2.0 * (57.7 - 0.118 * (tsn(:,:,1,jp_tem)+rt0))) / &
1701	& (82.0578 * (tsn(:,:,1,jp_tem)+rt0)))
1702	#else
1703	CALL ctl_stop( ' Trying to run fco2 observation operator', &
1704	& ' but no biogeochemical model appears to have been defined' )
1705	#endif
1706
1707	DO jfco2set = 1, nfco2sets
1708	CALL obs_fco2_opt( fco2datqc(jfco2set), &
1709	& kstp, jpi, jpj, nit000, fco2(:,:), &
1710	& maskfco2(:,:), n2dint )
1711	END DO
1712	ENDIF
1713
1714	IF ( ln_pco2 ) THEN
1715	maskpco2(:,:) = tmask(:,:,1) ! create a special mask to exclude certain things
1716	#if defined key_hadocc
1717	pco2(:,:) = HADOCC_PCO2(:,:) ! pCO2 from HadOCC
1718	IF ( ( MINVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ).AND.( MAXVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) ) THEN
1719	pco2(:,:) = obfillflt
1720	maskpco2(:,:) = 0
1721	CALL ctl_warn( ' HadOCC pCO2 values masked out for observation operator', &
1722	& ' on timestep ' // TRIM(STR(kstp)), &
1723	& ' as HADOCC_PCO2(:,:) == HADOCC_FILL_FLT' )
1724	ENDIF
1725	#elif defined key_medusa && defined key_foam_medusa
1726	pco2(:,:) = MEDUSA_PCO2(:,:) ! pCO2 from MEDUSA
1727	IF ( ( MINVAL( MEDUSA_PCO2 ) == MEDUSA_FILL_FLT ).AND.( MAXVAL( MEDUSA_PCO2 ) == MEDUSA_FILL_FLT ) ) THEN
1728	pco2(:,:) = obfillflt
1729	maskpco2(:,:) = 0
1730	CALL ctl_warn( ' MEDUSA pCO2 values masked out for observation operator', &
1731	& ' on timestep ' // TRIM(STR(kstp)), &
1732	& ' as MEDUSA_PCO2(:,:) == MEDUSA_FILL_FLT' )
1733	ENDIF
1734	#elif defined key_fabm
1735	pco2_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabm_o3pc)
1736	pco2(:,:) = pco2_3d(:,:,1)
1737	#else
1738	CALL ctl_stop( ' Trying to run pCO2 observation operator', &
1739	& ' but no biogeochemical model appears to have been defined' )
1740	#endif
1741
1742	DO jpco2set = 1, npco2sets
1743	CALL obs_pco2_opt( pco2datqc(jpco2set), &
1744	& kstp, jpi, jpj, nit000, pco2(:,:), &
1745	& maskpco2(:,:), n2dint )
1746	END DO
1747	ENDIF
1748
1749	#if ! defined key_lim2 && ! defined key_lim3
1750	CALL wrk_dealloc(jpi,jpj,frld)
1751	#endif
1752
1753	END SUBROUTINE dia_obs
1754
1755	SUBROUTINE dia_obs_wri
1756	!!----------------------------------------------------------------------
1757	!! * ROUTINE dia_obs_wri *
1758	!!
1759	!! ** Purpose : Call observation diagnostic output routines
1760	!!
1761	!! ** Method : Call observation diagnostic output routines
1762	!!
1763	!! ** Action :
1764	!!
1765	!! History :
1766	!! ! 06-03 (K. Mogensen) Original code
1767	!! ! 06-05 (K. Mogensen) Reformatted
1768	!! ! 06-10 (A. Weaver) Cleaning
1769	!! ! 07-03 (K. Mogensen) General handling of profiles
1770	!! ! 08-09 (M. Valdivieso) Velocity component (U,V) profiles
1771	!!----------------------------------------------------------------------
1772	IMPLICIT NONE
1773
1774	!! * Local declarations
1775
1776	INTEGER :: jprofset ! Profile data set loop variable
1777	INTEGER :: jveloset ! Velocity data set loop variable
1778	INTEGER :: jslaset ! SLA data set loop variable
1779	INTEGER :: jsstset ! SST data set loop variable
1780	INTEGER :: jseaiceset ! Sea Ice data set loop variable
1781	INTEGER :: jlogchlset ! logchl data set loop variable
1782	INTEGER :: jspmset ! spm data set loop variable
1783	INTEGER :: jfco2set ! fco2 data set loop variable
1784	INTEGER :: jpco2set ! pco2 data set loop variable
1785	INTEGER :: jset
1786	INTEGER :: jfbini
1787	CHARACTER(LEN=20) :: datestr=" ",timestr=" "
1788	CHARACTER(LEN=20) :: cdtmp
1789	!-----------------------------------------------------------------------
1790	! Depending on switches call various observation output routines
1791	!-----------------------------------------------------------------------
1792
1793	! - Temperature/salinity profiles
1794
1795	IF( ln_t3d .OR. ln_s3d ) THEN
1796
1797	! Copy data from prodatqc to profdata structures
1798	DO jprofset = 1, nprofsets
1799
1800	CALL obs_prof_decompress( prodatqc(jprofset), &
1801	& profdata(jprofset), .TRUE., numout )
1802
1803	END DO
1804
1805	! Write the profiles.
1806
1807	jprofset = 0
1808
1809	! ENACT insitu data
1810
1811	IF ( ln_ena ) THEN
1812
1813	jprofset = jprofset + 1
1814
1815	CALL obs_wri_p3d( 'enact', profdata(jprofset) )
1816
1817	ENDIF
1818
1819	! Coriolis insitu data
1820
1821	IF ( ln_cor ) THEN
1822
1823	jprofset = jprofset + 1
1824
1825	CALL obs_wri_p3d( 'corio', profdata(jprofset) )
1826
1827	ENDIF
1828
1829	! Feedback insitu data
1830
1831	IF ( ln_profb ) THEN
1832
1833	jfbini = jprofset + 1
1834
1835	DO jprofset = jfbini, nprofsets
1836
1837	jset = jprofset - jfbini + 1
1838	WRITE(cdtmp,'(A,I2.2)')'profb_',jset
1839	CALL obs_wri_p3d( cdtmp, profdata(jprofset) )
1840
1841	END DO
1842
1843	ENDIF
1844
1845	ENDIF
1846
1847	! - Sea surface anomaly
1848	IF ( ln_sla ) THEN
1849
1850	! Copy data from sladatqc to sladata structures
1851	DO jslaset = 1, nslasets
1852
1853	CALL obs_surf_decompress( sladatqc(jslaset), &
1854	& sladata(jslaset), .TRUE., numout )
1855
1856	END DO
1857
1858	jslaset = 0
1859
1860	! Write the AVISO SLA data
1861
1862	IF ( ln_sladt ) THEN
1863
1864	jslaset = 1
1865	CALL obs_wri_sla( 'aviso_act', sladata(jslaset) )
1866	jslaset = 2
1867	CALL obs_wri_sla( 'aviso_pas', sladata(jslaset) )
1868
1869	ENDIF
1870
1871	IF ( ln_slafb ) THEN
1872
1873	jfbini = jslaset + 1
1874
1875	DO jslaset = jfbini, nslasets
1876
1877	jset = jslaset - jfbini + 1
1878	WRITE(cdtmp,'(A,I2.2)')'slafb_',jset
1879	CALL obs_wri_sla( cdtmp, sladata(jslaset) )
1880
1881	END DO
1882
1883	ENDIF
1884
1885	ENDIF
1886
1887	! - Sea surface temperature
1888	IF ( ln_sst ) THEN
1889
1890	! Copy data from sstdatqc to sstdata structures
1891	DO jsstset = 1, nsstsets
1892
1893	CALL obs_surf_decompress( sstdatqc(jsstset), &
1894	& sstdata(jsstset), .TRUE., numout )
1895
1896	END DO
1897
1898	jsstset = 0
1899
1900	! Write the AVISO SST data
1901
1902	IF ( ln_reysst ) THEN
1903
1904	jsstset = jsstset + 1
1905	CALL obs_wri_sst( 'reynolds', sstdata(jsstset) )
1906
1907	ENDIF
1908
1909	IF ( ln_ghrsst ) THEN
1910
1911	jsstset = jsstset + 1
1912	CALL obs_wri_sst( 'ghr', sstdata(jsstset) )
1913
1914	ENDIF
1915
1916	IF ( ln_sstfb ) THEN
1917
1918	jfbini = jsstset + 1
1919
1920	DO jsstset = jfbini, nsstsets
1921
1922	jset = jsstset - jfbini + 1
1923	WRITE(cdtmp,'(A,I2.2)')'sstfb_',jset
1924	CALL obs_wri_sst( cdtmp, sstdata(jsstset) )
1925
1926	END DO
1927
1928	ENDIF
1929
1930	ENDIF
1931
1932	! - Sea surface salinity
1933	IF ( ln_sss ) THEN
1934	IF(lwp) WRITE(numout,*) ' SSS currently not available'
1935	ENDIF
1936
1937	! - Sea Ice Concentration
1938	IF ( ln_seaice ) THEN
1939
1940	! Copy data from seaicedatqc to seaicedata structures
1941	DO jseaiceset = 1, nseaicesets
1942
1943	CALL obs_surf_decompress( seaicedatqc(jseaiceset), &
1944	& seaicedata(jseaiceset), .TRUE., numout )
1945
1946	END DO
1947
1948	! Write the Sea Ice data
1949	DO jseaiceset = 1, nseaicesets
1950
1951	WRITE(cdtmp,'(A,I2.2)')'seaicefb_',jseaiceset
1952	CALL obs_wri_seaice( cdtmp, seaicedata(jseaiceset) )
1953
1954	END DO
1955
1956	ENDIF
1957
1958	! Velocity data
1959	IF( ln_vel3d ) THEN
1960
1961	! Copy data from veldatqc to velodata structures
1962	DO jveloset = 1, nvelosets
1963
1964	CALL obs_prof_decompress( veldatqc(jveloset), &
1965	& velodata(jveloset), .TRUE., numout )
1966
1967	END DO
1968
1969	! Write the profiles.
1970
1971	jveloset = 0
1972
1973	! Daily averaged data
1974
1975	IF ( ln_velavcur ) THEN
1976
1977	jveloset = jveloset + 1
1978
1979	CALL obs_wri_vel( 'velavcurr', velodata(jveloset), n2dint )
1980
1981	ENDIF
1982
1983	! High frequency data
1984
1985	IF ( ln_velhrcur ) THEN
1986
1987	jveloset = jveloset + 1
1988
1989	CALL obs_wri_vel( 'velhrcurr', velodata(jveloset), n2dint )
1990
1991	ENDIF
1992
1993	! Daily averaged data
1994
1995	IF ( ln_velavadcp ) THEN
1996
1997	jveloset = jveloset + 1
1998
1999	CALL obs_wri_vel( 'velavadcp', velodata(jveloset), n2dint )
2000
2001	ENDIF
2002
2003	! High frequency data
2004
2005	IF ( ln_velhradcp ) THEN
2006
2007	jveloset = jveloset + 1
2008
2009	CALL obs_wri_vel( 'velhradcp', velodata(jveloset), n2dint )
2010
2011	ENDIF
2012
2013	! Feedback velocity data
2014
2015	IF ( ln_velfb ) THEN
2016
2017	jfbini = jveloset + 1
2018
2019	DO jveloset = jfbini, nvelosets
2020
2021	jset = jveloset - jfbini + 1
2022	WRITE(cdtmp,'(A,I2.2)')'velfb_',jset
2023	CALL obs_wri_vel( cdtmp, velodata(jveloset), n2dint )
2024
2025	END DO
2026
2027	ENDIF
2028
2029	ENDIF
2030
2031	! - log10(chlorophyll)
2032	IF ( ln_logchl ) THEN
2033
2034	! Copy data from logchldatqc to logchldata structures
2035	DO jlogchlset = 1, nlogchlsets
2036
2037	CALL obs_surf_decompress( logchldatqc(jlogchlset), &
2038	& logchldata(jlogchlset), .TRUE., numout )
2039
2040	END DO
2041
2042	! Mark as bad observations with no valid model counterpart due to activities in dia_obs
2043	! Seem to need to set to fill value rather than marking as bad to be effective, so do both
2044	DO jlogchlset = 1, nlogchlsets
2045	WHERE ( logchldata(jlogchlset)%rmod(:,1) == obfillflt )
2046	logchldata(jlogchlset)%nqc(:) = 1
2047	logchldata(jlogchlset)%robs(:,1) = obfillflt
2048	END WHERE
2049	END DO
2050
2051	! Write the logchl data
2052	! If padd is ever defined and added as an optional argument
2053	! to the obs_wri_logchl call, then the observed STD values
2054	! will need adding to padd in order to be written out
2055	DO jlogchlset = 1, nlogchlsets
2056
2057	WRITE(cdtmp,'(A,I2.2)')'logchlfb_',jlogchlset
2058	CALL obs_wri_logchl( cdtmp, logchldata(jlogchlset) )
2059
2060	END DO
2061
2062	ENDIF
2063
2064	! - spm
2065	IF ( ln_spm ) THEN
2066
2067	! Copy data from spmdatqc to spmdata structures
2068	DO jspmset = 1, nspmsets
2069
2070	CALL obs_surf_decompress( spmdatqc(jspmset), &
2071	& spmdata(jspmset), .TRUE., numout )
2072
2073	END DO
2074
2075	! Write the spm data
2076	DO jspmset = 1, nspmsets
2077
2078	WRITE(cdtmp,'(A,I2.2)')'spmfb_',jspmset
2079	CALL obs_wri_spm( cdtmp, spmdata(jspmset) )
2080
2081	END DO
2082
2083	ENDIF
2084
2085	! - fco2
2086	IF ( ln_fco2 ) THEN
2087
2088	! Copy data from fco2datqc to fco2data structures
2089	DO jfco2set = 1, nfco2sets
2090
2091	CALL obs_surf_decompress( fco2datqc(jfco2set), &
2092	& fco2data(jfco2set), .TRUE., numout )
2093
2094	END DO
2095
2096	! Mark as bad observations with no valid model counterpart due to fCO2 not being in the restart
2097	! Seem to need to set to fill value rather than marking as bad to be effective, so do both
2098	DO jfco2set = 1, nfco2sets
2099	WHERE ( fco2data(jfco2set)%rmod(:,1) == obfillflt )
2100	fco2data(jfco2set)%nqc(:) = 1
2101	fco2data(jfco2set)%robs(:,1) = obfillflt
2102	END WHERE
2103	END DO
2104
2105	! Write the fco2 data
2106	DO jfco2set = 1, nfco2sets
2107
2108	WRITE(cdtmp,'(A,I2.2)')'fco2fb_',jfco2set
2109	CALL obs_wri_fco2( cdtmp, fco2data(jfco2set) )
2110
2111	END DO
2112
2113	ENDIF
2114
2115	! - pco2
2116	IF ( ln_pco2 ) THEN
2117
2118	! Copy data from pco2datqc to pco2data structures
2119	DO jpco2set = 1, npco2sets
2120
2121	CALL obs_surf_decompress( pco2datqc(jpco2set), &
2122	& pco2data(jpco2set), .TRUE., numout )
2123
2124	END DO
2125
2126	! Mark as bad observations with no valid model counterpart due to pco2 not being in the restart
2127	! Seem to need to set to fill value rather than marking as bad to be effective, so do both
2128	DO jpco2set = 1, npco2sets
2129	WHERE ( pco2data(jpco2set)%rmod(:,1) == obfillflt )
2130	pco2data(jpco2set)%nqc(:) = 1
2131	pco2data(jpco2set)%robs(:,1) = obfillflt
2132	END WHERE
2133	END DO
2134
2135	! Write the pco2 data
2136	DO jpco2set = 1, npco2sets
2137
2138	WRITE(cdtmp,'(A,I2.2)')'pco2fb_',jpco2set
2139	CALL obs_wri_pco2( cdtmp, pco2data(jpco2set) )
2140
2141	END DO
2142
2143	ENDIF
2144
2145	END SUBROUTINE dia_obs_wri
2146
2147	SUBROUTINE dia_obs_dealloc
2148	IMPLICIT NONE
2149	!!----------------------------------------------------------------------
2150	!! * ROUTINE dia_obs_dealloc *
2151	!!
2152	!! ** Purpose : To deallocate data to enable the obs_oper online loop.
2153	!! Specifically: dia_obs_init --> dia_obs --> dia_obs_wri
2154	!!
2155	!! ** Method : Clean up various arrays left behind by the obs_oper.
2156	!!
2157	!! ** Action :
2158	!!
2159	!!----------------------------------------------------------------------
2160	!! obs_grid deallocation
2161	CALL obs_grid_deallocate
2162
2163	!! diaobs deallocation
2164	IF ( nprofsets > 0 ) THEN
2165	DEALLOCATE(ld_enact, &
2166	& profdata, &
2167	& prodatqc)
2168	END IF
2169	IF ( ln_sla ) THEN
2170	DEALLOCATE(sladata, &
2171	& sladatqc)
2172	END IF
2173	IF ( ln_seaice ) THEN
2174	DEALLOCATE(sladata, &
2175	& sladatqc)
2176	END IF
2177	IF ( ln_sst ) THEN
2178	DEALLOCATE(sstdata, &
2179	& sstdatqc)
2180	END IF
2181	IF ( ln_vel3d ) THEN
2182	DEALLOCATE(ld_velav, &
2183	& velodata, &
2184	& veldatqc)
2185	END IF
2186	END SUBROUTINE dia_obs_dealloc
2187
2188	SUBROUTINE ini_date( ddobsini )
2189	!!----------------------------------------------------------------------
2190	!! * ROUTINE ini_date *
2191	!!
2192	!! ** Purpose : Get initial data in double precision YYYYMMDD.HHMMSS format
2193	!!
2194	!! ** Method : Get initial data in double precision YYYYMMDD.HHMMSS format
2195	!!
2196	!! ** Action : Get initial data in double precision YYYYMMDD.HHMMSS format
2197	!!
2198	!! History :
2199	!! ! 06-03 (K. Mogensen) Original code
2200	!! ! 06-05 (K. Mogensen) Reformatted
2201	!! ! 06-10 (A. Weaver) Cleaning
2202	!! ! 06-10 (G. Smith) Calculates initial date the same as method for final date
2203	!! ! 10-05 (D. Lea) Update to month length calculation for NEMO vn3.2
2204	!!----------------------------------------------------------------------
2205	USE phycst, ONLY : & ! Physical constants
2206	& rday
2207	! USE daymod, ONLY : & ! Time variables
2208	! & nmonth_len
2209	USE dom_oce, ONLY : & ! Ocean space and time domain variables
2210	& rdt
2211
2212	IMPLICIT NONE
2213
2214	!! * Arguments
2215	REAL(KIND=dp), INTENT(OUT) :: ddobsini ! Initial date in YYYYMMDD.HHMMSS
2216
2217	!! * Local declarations
2218	INTEGER :: iyea ! date - (year, month, day, hour, minute)
2219	INTEGER :: imon
2220	INTEGER :: iday
2221	INTEGER :: ihou
2222	INTEGER :: imin
2223	INTEGER :: imday ! Number of days in month.
2224	REAL(KIND=wp) :: zdayfrc ! Fraction of day
2225
2226	INTEGER, DIMENSION(12) :: imonth_len !: length in days of the months of the current year
2227
2228	!!----------------------------------------------------------------------
2229	!! Initial date initialization (year, month, day, hour, minute)
2230	!! (This assumes that the initial date is for 00z))
2231	!!----------------------------------------------------------------------
2232	iyea = ndate0 / 10000
2233	imon = ( ndate0 - iyea * 10000 ) / 100
2234	iday = ndate0 - iyea * 10000 - imon * 100
2235	ihou = 0
2236	imin = 0
2237
2238	!!----------------------------------------------------------------------
2239	!! Compute number of days + number of hours + min since initial time
2240	!!----------------------------------------------------------------------
2241	iday = iday + ( nit000 -1 ) * rdt / rday
2242	zdayfrc = ( nit000 -1 ) * rdt / rday
2243	zdayfrc = zdayfrc - aint(zdayfrc)
2244	ihou = int( zdayfrc * 24 )
2245	imin = int( (zdayfrc * 24 - ihou) * 60 )
2246
2247	!!-----------------------------------------------------------------------
2248	!! Convert number of days (iday) into a real date
2249	!!----------------------------------------------------------------------
2250
2251	CALL calc_month_len( iyea, imonth_len )
2252
2253	DO WHILE ( iday > imonth_len(imon) )
2254	iday = iday - imonth_len(imon)
2255	imon = imon + 1
2256	IF ( imon > 12 ) THEN
2257	imon = 1
2258	iyea = iyea + 1
2259	CALL calc_month_len( iyea, imonth_len ) ! update month lengths
2260	ENDIF
2261	END DO
2262
2263	!!----------------------------------------------------------------------
2264	!! Convert it into YYYYMMDD.HHMMSS format.
2265	!!----------------------------------------------------------------------
2266	ddobsini = iyea * 10000_dp + imon * 100_dp + &
2267	& iday + ihou * 0.01_dp + imin * 0.0001_dp
2268
2269
2270	END SUBROUTINE ini_date
2271
2272	SUBROUTINE fin_date( ddobsfin )
2273	!!----------------------------------------------------------------------
2274	!! * ROUTINE fin_date *
2275	!!
2276	!! ** Purpose : Get final data in double precision YYYYMMDD.HHMMSS format
2277	!!
2278	!! ** Method : Get final data in double precision YYYYMMDD.HHMMSS format
2279	!!
2280	!! ** Action : Get final data in double precision YYYYMMDD.HHMMSS format
2281	!!
2282	!! History :
2283	!! ! 06-03 (K. Mogensen) Original code
2284	!! ! 06-05 (K. Mogensen) Reformatted
2285	!! ! 06-10 (A. Weaver) Cleaning
2286	!! ! 10-05 (D. Lea) Update to month length calculation for NEMO vn3.2
2287	!!----------------------------------------------------------------------
2288	USE phycst, ONLY : & ! Physical constants
2289	& rday
2290	! USE daymod, ONLY : & ! Time variables
2291	! & nmonth_len
2292	USE dom_oce, ONLY : & ! Ocean space and time domain variables
2293	& rdt
2294
2295	IMPLICIT NONE
2296
2297	!! * Arguments
2298	REAL(KIND=dp), INTENT(OUT) :: ddobsfin ! Final date in YYYYMMDD.HHMMSS
2299
2300	!! * Local declarations
2301	INTEGER :: iyea ! date - (year, month, day, hour, minute)
2302	INTEGER :: imon
2303	INTEGER :: iday
2304	INTEGER :: ihou
2305	INTEGER :: imin
2306	INTEGER :: imday ! Number of days in month.
2307	REAL(KIND=wp) :: zdayfrc ! Fraction of day
2308
2309	INTEGER, DIMENSION(12) :: imonth_len !: length in days of the months of the current year
2310
2311	!-----------------------------------------------------------------------
2312	! Initial date initialization (year, month, day, hour, minute)
2313	! (This assumes that the initial date is for 00z)
2314	!-----------------------------------------------------------------------
2315	iyea = ndate0 / 10000
2316	imon = ( ndate0 - iyea * 10000 ) / 100
2317	iday = ndate0 - iyea * 10000 - imon * 100
2318	ihou = 0
2319	imin = 0
2320
2321	!-----------------------------------------------------------------------
2322	! Compute number of days + number of hours + min since initial time
2323	!-----------------------------------------------------------------------
2324	iday = iday + nitend * rdt / rday
2325	zdayfrc = nitend * rdt / rday
2326	zdayfrc = zdayfrc - AINT( zdayfrc )
2327	ihou = INT( zdayfrc * 24 )
2328	imin = INT( ( zdayfrc * 24 - ihou ) * 60 )
2329
2330	!-----------------------------------------------------------------------
2331	! Convert number of days (iday) into a real date
2332	!----------------------------------------------------------------------
2333
2334	CALL calc_month_len( iyea, imonth_len )
2335
2336	DO WHILE ( iday > imonth_len(imon) )
2337	iday = iday - imonth_len(imon)
2338	imon = imon + 1
2339	IF ( imon > 12 ) THEN
2340	imon = 1
2341	iyea = iyea + 1
2342	CALL calc_month_len( iyea, imonth_len ) ! update month lengths
2343	ENDIF
2344	END DO
2345
2346	!-----------------------------------------------------------------------
2347	! Convert it into YYYYMMDD.HHMMSS format
2348	!-----------------------------------------------------------------------
2349	ddobsfin = iyea * 10000_dp + imon * 100_dp + iday &
2350	& + ihou * 0.01_dp + imin * 0.0001_dp
2351
2352	END SUBROUTINE fin_date
2353
2354	END MODULE diaobs

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90

Download in other formats: