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

Context Navigation

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 2287

Last change on this file since 2287 was 2287, checked in by smasson, 14 years ago
update licence of all NEMO files...
Property svn:keywords set to `Id`
File size: 61.4 KB

Line
1	MODULE diaobs
2	!!======================================================================
3	!! * MODULE diaobs *
4	!! Observation diagnostics: Computation of the misfit between data and
5	!! their model equivalent
6	!!======================================================================
7
8	!!----------------------------------------------------------------------
9	!! '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 par_kind ! Precision variables
19	USE in_out_manager ! I/O manager
20	USE par_oce
21	USE dom_oce ! Ocean space and time domain variables
22	USE obs_read_prof ! Reading and allocation of observations (Coriolis)
23	USE obs_read_sla ! Reading and allocation of SLA observations
24	USE obs_read_sst ! Reading and allocation of SST observations
25	USE obs_readmdt ! Reading and allocation of MDT for SLA.
26	USE obs_read_seaice ! Reading and allocation of Sea Ice observations
27	USE obs_read_vel ! Reading and allocation of velocity component observations
28	USE obs_prep ! Preparation of obs. (grid search etc).
29	USE obs_oper ! Observation operators
30	USE obs_write ! Writing of observation related diagnostics
31	USE obs_grid ! Grid searching
32	USE obs_read_altbias ! Bias treatment for altimeter
33	USE obs_profiles_def ! Profile data definitions
34	USE obs_profiles ! Profile data storage
35	USE obs_surf_def ! Surface data definitions
36	USE obs_sla ! SLA data storage
37	USE obs_sst ! SST data storage
38	USE obs_seaice ! Sea Ice data storage
39	USE obs_types ! Definitions for observation types
40	USE mpp_map ! MPP mapping
41
42	IMPLICIT NONE
43
44	!! * Routine accessibility
45	PRIVATE
46	PUBLIC dia_obs_init, & ! Initialize and read observations
47	& dia_obs, & ! Compute model equivalent to observations
48	& dia_obs_wri ! Write model equivalent to observations
49
50	!! * Shared Module variables
51	LOGICAL, PUBLIC, PARAMETER :: &
52	#if defined key_diaobs
53	& lk_diaobs = .TRUE. !: Logical switch for observation diangostics
54	#else
55	& lk_diaobs = .FALSE. !: Logical switch for observation diangostics
56	#endif
57
58	!! * Module variables
59	LOGICAL, PUBLIC :: ln_t3d !: Logical switch for temperature profiles
60	LOGICAL, PUBLIC :: ln_s3d !: Logical switch for salinity profiles
61	LOGICAL, PUBLIC :: ln_ena !: Logical switch for the ENACT data set
62	LOGICAL, PUBLIC :: ln_cor !: Logical switch for the Coriolis data set
63	LOGICAL, PUBLIC :: ln_profb !: Logical switch for profile feedback datafiles
64	LOGICAL, PUBLIC :: ln_sla !: Logical switch for sea level anomalies
65	LOGICAL, PUBLIC :: ln_sladt !: Logical switch for SLA from AVISO files
66	LOGICAL, PUBLIC :: ln_slafb !: Logical switch for SLA from feedback files
67	LOGICAL, PUBLIC :: ln_sst !: Logical switch for sea surface temperature
68	LOGICAL, PUBLIC :: ln_reysst !: Logical switch for Reynolds sea surface temperature
69	LOGICAL, PUBLIC :: ln_ghrsst !: Logical switch for GHRSST data
70	LOGICAL, PUBLIC :: ln_sstfb !: Logical switch for SST from feedback files
71	LOGICAL, PUBLIC :: ln_seaice !: Logical switch for sea ice concentration
72	LOGICAL, PUBLIC :: ln_vel3d !: Logical switch for velocity component (u,v) observations
73	LOGICAL, PUBLIC :: ln_velavcur !: Logical switch for raw daily averaged netCDF current meter vel. data
74	LOGICAL, PUBLIC :: ln_velhrcur !: Logical switch for raw high freq netCDF current meter vel. data
75	LOGICAL, PUBLIC :: ln_velavadcp !: Logical switch for raw daily averaged netCDF ADCP vel. data
76	LOGICAL, PUBLIC :: ln_velhradcp !: Logical switch for raw high freq netCDF ADCP vel. data
77	LOGICAL, PUBLIC :: ln_velfb !: Logical switch for velocities from feedback files
78	LOGICAL, PUBLIC :: ln_ssh !: Logical switch for sea surface height
79	LOGICAL, PUBLIC :: ln_sss !: Logical switch for sea surface salinity
80	LOGICAL, PUBLIC :: ln_nea !: Remove observations near land
81	LOGICAL, PUBLIC :: ln_altbias !: Logical switch for altimeter bias
82	LOGICAL, PUBLIC :: ln_ignmis !: Logical switch for ignoring missing files
83	LOGICAL, PUBLIC :: ln_s_at_t !: Logical switch to compute model S at T observations
84
85	REAL(KIND=dp), PUBLIC :: dobsini !: Observation window start date YYYYMMDD.HHMMSS
86	REAL(KIND=dp), PUBLIC :: dobsend !: Observation window end date YYYYMMDD.HHMMSS
87
88	INTEGER, PUBLIC :: n1dint !: Vertical interpolation method
89	INTEGER, PUBLIC :: n2dint !: Horizontal interpolation method
90
91	INTEGER, DIMENSION(imaxavtypes) :: &
92	& endailyavtypes !: ENACT data types which are daily average
93
94	INTEGER, PARAMETER :: MaxNumFiles = 1000
95	LOGICAL, DIMENSION(MaxNumFiles) :: &
96	& ln_profb_ena, & !: Is the feedback files from ENACT data ?
97	! !: If so use endailyavtypes
98	& ln_profb_enatim !: Change tim for 820 enact data set.
99
100	LOGICAL, DIMENSION(MaxNumFiles) :: &
101	& ln_velfb_av !: Is the velocity feedback files daily average?
102	LOGICAL, DIMENSION(:), ALLOCATABLE :: &
103	& ld_enact !: Profile data is ENACT so use endailyavtypes
104	LOGICAL, DIMENSION(:), ALLOCATABLE :: &
105	& ld_velav !: Velocity data is daily averaged
106
107	!!----------------------------------------------------------------------
108	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
109	!! $Id$
110	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
111	!!----------------------------------------------------------------------
112
113	CONTAINS
114
115	SUBROUTINE dia_obs_init
116	!!----------------------------------------------------------------------
117	!! * ROUTINE dia_obs_init *
118	!!
119	!! ** Purpose : Initialize and read observations
120	!!
121	!! ** Method : Read the namelist and call reading routines
122	!!
123	!! ** Action : Read the namelist and call reading routines
124	!!
125	!! History :
126	!! ! 06-03 (K. Mogensen) Original code
127	!! ! 06-05 (A. Weaver) Reformatted
128	!! ! 06-10 (A. Weaver) Cleaning and add controls
129	!! ! 07-03 (K. Mogensen) General handling of profiles
130	!!----------------------------------------------------------------------
131
132	IMPLICIT NONE
133
134	!! * Local declarations
135	CHARACTER(len=128) :: enactfiles(MaxNumFiles)
136	CHARACTER(len=128) :: coriofiles(MaxNumFiles)
137	CHARACTER(len=128) :: profbfiles(MaxNumFiles)
138	CHARACTER(len=128) :: sstfiles(MaxNumFiles)
139	CHARACTER(len=128) :: sstfbfiles(MaxNumFiles)
140	CHARACTER(len=128) :: slafilesact(MaxNumFiles)
141	CHARACTER(len=128) :: slafilespas(MaxNumFiles)
142	CHARACTER(len=128) :: slafbfiles(MaxNumFiles)
143	CHARACTER(len=128) :: seaicefiles(MaxNumFiles)
144	CHARACTER(len=128) :: velcurfiles(MaxNumFiles)
145	CHARACTER(len=128) :: veladcpfiles(MaxNumFiles)
146	CHARACTER(len=128) :: velavcurfiles(MaxNumFiles)
147	CHARACTER(len=128) :: velhrcurfiles(MaxNumFiles)
148	CHARACTER(len=128) :: velavadcpfiles(MaxNumFiles)
149	CHARACTER(len=128) :: velhradcpfiles(MaxNumFiles)
150	CHARACTER(len=128) :: velfbfiles(MaxNumFiles)
151	CHARACTER(LEN=128) :: reysstname
152	CHARACTER(LEN=12) :: reysstfmt
153	CHARACTER(LEN=128) :: bias_file
154	CHARACTER(LEN=20) :: datestr=" ", timestr=" "
155	NAMELIST/namobs/ln_ena, ln_cor, ln_profb, ln_t3d, ln_s3d, &
156	& ln_sla, ln_sladt, ln_slafb, &
157	& ln_ssh, ln_sst, ln_sstfb, ln_sss, ln_nea, &
158	& enactfiles, coriofiles, profbfiles, &
159	& slafilesact, slafilespas, slafbfiles, &
160	& sstfiles, sstfbfiles, &
161	& ln_seaice, seaicefiles, &
162	& dobsini, dobsend, n1dint, n2dint, &
163	& nmsshc, mdtcorr, mdtcutoff, &
164	& ln_reysst, ln_ghrsst, reysstname, reysstfmt, &
165	& ln_grid_search_lookup, &
166	& grid_search_file, grid_search_res, &
167	& ln_grid_global, bias_file, ln_altbias, &
168	& endailyavtypes, ln_s_at_t, ln_profb_ena, &
169	& ln_vel3d, ln_velavcur, velavcurfiles, &
170	& ln_velhrcur, velhrcurfiles, &
171	& ln_velavadcp, velavadcpfiles, &
172	& ln_velhradcp, velhradcpfiles, &
173	& ln_velfb, velfbfiles, ln_velfb_av, &
174	& ln_profb_enatim, ln_ignmis
175
176	INTEGER :: jprofset
177	INTEGER :: jveloset
178	INTEGER :: jvar
179	INTEGER :: jnumenact
180	INTEGER :: jnumcorio
181	INTEGER :: jnumprofb
182	INTEGER :: jnumslaact
183	INTEGER :: jnumslapas
184	INTEGER :: jnumslafb
185	INTEGER :: jnumsst
186	INTEGER :: jnumsstfb
187	INTEGER :: jnumseaice
188	INTEGER :: jnumvelavcur
189	INTEGER :: jnumvelhrcur
190	INTEGER :: jnumvelavadcp
191	INTEGER :: jnumvelhradcp
192	INTEGER :: jnumvelfb
193	INTEGER :: ji
194	INTEGER :: jset
195	LOGICAL :: lmask(MaxNumFiles), ll_u3d, ll_v3d
196
197	!-----------------------------------------------------------------------
198	! Read namelist parameters
199	!-----------------------------------------------------------------------
200	! default values already set except:
201
202	ln_t3d = .FALSE.
203	ln_s3d = .FALSE.
204	ln_vel3d = .FALSE.
205	ln_sla = .FALSE.
206	ln_altbias = .FALSE.
207	ln_ssh = .FALSE.
208	ln_sst = .FALSE.
209	ln_seaice = .FALSE.
210	ln_reysst = .FALSE.
211	ln_ghrsst = .FALSE.
212	ln_sss = .FALSE.
213	ln_profb = .FALSE.
214	ln_ena = .TRUE.
215	ln_cor = .FALSE.
216	ln_sladt = .TRUE.
217	ln_slafb = .FALSE.
218	ln_sstfb = .FALSE.
219	ln_velavcur = .FALSE.
220	ln_velhrcur = .FALSE.
221	ln_velavadcp = .FALSE.
222	ln_velhradcp = .FALSE.
223	ln_velfb = .FALSE.
224	ln_nea = .FALSE.
225	ln_grid_search_lookup = .FALSE.
226	ln_grid_global = .FALSE.
227	ln_s_at_t = .TRUE.
228	grid_search_file = 'xypos'
229	bias_file='bias.nc'
230	enactfiles(:) = ''
231	coriofiles(:) = ''
232	profbfiles(:) = ''
233	slafilesact(:) = ''
234	slafilespas(:) = ''
235	slafbfiles(:) = ''
236	sstfiles(:) = ''
237	sstfbfiles(:) = ''
238	seaicefiles(:) = ''
239	velcurfiles(:) = ''
240	veladcpfiles(:) = ''
241	velavcurfiles(:) = ''
242	velhrcurfiles(:) = ''
243	velavadcpfiles(:) = ''
244	velhradcpfiles(:) = ''
245	velfbfiles(:) = ''
246	reysstname = 'sst_yYYYYmMM.nc'
247	reysstfmt = 'monthly'
248	endailyavtypes(:) = -1
249	endailyavtypes(1) = 820
250	ln_profb_ena(:) = .FALSE.
251	ln_profb_enatim(:) = .TRUE.
252	ln_velfb_av(:) = .FALSE.
253	ln_ignmis = .FALSE.
254	CALL ini_date( dobsini )
255	CALL fin_date( dobsend )
256	n1dint = 1
257	n2dint = 3
258
259	! Read Namelist namobs : control observation diagnostics
260	REWIND( numnam )
261	READ ( numnam, namobs )
262
263	! Count number of files for each type
264	IF (ln_ena) THEN
265	lmask(:) = .FALSE.
266	WHERE (enactfiles(:) /= '') lmask(:) = .TRUE.
267	jnumenact = COUNT(lmask)
268	ENDIF
269	IF (ln_cor) THEN
270	lmask(:) = .FALSE.
271	WHERE (coriofiles(:) /= '') lmask(:) = .TRUE.
272	jnumcorio = COUNT(lmask)
273	ENDIF
274	IF (ln_profb) THEN
275	lmask(:) = .FALSE.
276	WHERE (profbfiles(:) /= '') lmask(:) = .TRUE.
277	jnumprofb = COUNT(lmask)
278	ENDIF
279	IF (ln_sladt) THEN
280	lmask(:) = .FALSE.
281	WHERE (slafilesact(:) /= '') lmask(:) = .TRUE.
282	jnumslaact = COUNT(lmask)
283	lmask(:) = .FALSE.
284	WHERE (slafilespas(:) /= '') lmask(:) = .TRUE.
285	jnumslapas = COUNT(lmask)
286	ENDIF
287	IF (ln_slafb) THEN
288	lmask(:) = .FALSE.
289	WHERE (slafbfiles(:) /= '') lmask(:) = .TRUE.
290	jnumslafb = COUNT(lmask)
291	lmask(:) = .FALSE.
292	ENDIF
293	IF (ln_ghrsst) THEN
294	lmask(:) = .FALSE.
295	WHERE (sstfiles(:) /= '') lmask(:) = .TRUE.
296	jnumsst = COUNT(lmask)
297	ENDIF
298	IF (ln_sstfb) THEN
299	lmask(:) = .FALSE.
300	WHERE (sstfbfiles(:) /= '') lmask(:) = .TRUE.
301	jnumsstfb = COUNT(lmask)
302	lmask(:) = .FALSE.
303	ENDIF
304	IF (ln_seaice) THEN
305	lmask(:) = .FALSE.
306	WHERE (seaicefiles(:) /= '') lmask(:) = .TRUE.
307	jnumseaice = COUNT(lmask)
308	ENDIF
309	IF (ln_velavcur) THEN
310	lmask(:) = .FALSE.
311	WHERE (velavcurfiles(:) /= '') lmask(:) = .TRUE.
312	jnumvelavcur = COUNT(lmask)
313	ENDIF
314	IF (ln_velhrcur) THEN
315	lmask(:) = .FALSE.
316	WHERE (velhrcurfiles(:) /= '') lmask(:) = .TRUE.
317	jnumvelhrcur = COUNT(lmask)
318	ENDIF
319	IF (ln_velavadcp) THEN
320	lmask(:) = .FALSE.
321	WHERE (velavadcpfiles(:) /= '') lmask(:) = .TRUE.
322	jnumvelavadcp = COUNT(lmask)
323	ENDIF
324	IF (ln_velhradcp) THEN
325	lmask(:) = .FALSE.
326	WHERE (velhradcpfiles(:) /= '') lmask(:) = .TRUE.
327	jnumvelhradcp = COUNT(lmask)
328	ENDIF
329	IF (ln_velfb) THEN
330	lmask(:) = .FALSE.
331	WHERE (velfbfiles(:) /= '') lmask(:) = .TRUE.
332	jnumvelfb = COUNT(lmask)
333	lmask(:) = .FALSE.
334	ENDIF
335
336	! Control print
337	IF(lwp) THEN
338	WRITE(numout,*)
339	WRITE(numout,*) 'dia_obs_init : Observation diagnostic initialization'
340	WRITE(numout,*) '~~~~~~~~~~~~'
341	#if defined key_datetime_out
342	CALL DATE_AND_TIME(datestr,timestr)
343	WRITE(numout,*) 'dia_obs_init date_and_time ',datestr,' ',timestr
344	#endif
345	WRITE(numout,*) ' Namelist namobs : set observation diagnostic parameters'
346	WRITE(numout,*) ' Logical switch for T profile observations ln_t3d = ', ln_t3d
347	WRITE(numout,*) ' Logical switch for S profile observations ln_s3d = ', ln_s3d
348	WRITE(numout,*) ' Logical switch for ENACT insitu data set ln_ena = ', ln_ena
349	WRITE(numout,*) ' Logical switch for Coriolis insitu data set ln_cor = ', ln_cor
350	WRITE(numout,*) ' Logical switch for feedback insitu data set ln_profb = ', ln_profb
351	WRITE(numout,*) ' Logical switch for SLA observations ln_sla = ', ln_sla
352	WRITE(numout,*) ' Logical switch for AVISO SLA data ln_sladt = ', ln_sladt
353	WRITE(numout,*) ' Logical switch for feedback SLA data ln_slafb = ', ln_slafb
354	WRITE(numout,*) ' Logical switch for SSH observations ln_ssh = ', ln_ssh
355	WRITE(numout,*) ' Logical switch for SST observations ln_sst = ', ln_sst
356	WRITE(numout,*) ' Logical switch for Reynolds observations ln_reysst = ', ln_reysst
357	WRITE(numout,*) ' Logical switch for GHRSST observations ln_ghrsst = ', ln_ghrsst
358	WRITE(numout,*) ' Logical switch for feedback SST data ln_sstfb = ', ln_sstfb
359	WRITE(numout,*) ' Logical switch for SSS observations ln_sss = ', ln_sss
360	WRITE(numout,*) ' Logical switch for Sea Ice observations ln_seaice = ', ln_seaice
361	WRITE(numout,*) ' Logical switch for velocity observations ln_vel3d = ', ln_vel3d
362	WRITE(numout,*) ' Logical switch for velocity daily av. cur. ln_velavcur = ', ln_velavcur
363	WRITE(numout,*) ' Logical switch for velocity high freq. cur. ln_velhrcur = ', ln_velhrcur
364	WRITE(numout,*) ' Logical switch for velocity daily av. ADCP ln_velavadcp = ', ln_velavadcp
365	WRITE(numout,*) ' Logical switch for velocity high freq. ADCP ln_velhradcp = ', ln_velhradcp
366	WRITE(numout,*) ' Logical switch for feedback velocity data ln_velfb = ', ln_velfb
367	WRITE(numout,*) ' Global distribtion of observations ln_grid_global = ',ln_grid_global
368	WRITE(numout,*) &
369	' Logical switch for obs grid search w/lookup table ln_grid_search_lookup = ',ln_grid_search_lookup
370	IF (ln_grid_search_lookup) &
371	WRITE(numout,*) ' Grid search lookup file header grid_search_file = ', grid_search_file
372	IF (ln_ena) THEN
373	DO ji = 1, jnumenact
374	WRITE(numout,'(1X,2A)') ' ENACT input observation file name enactfiles = ', &
375	TRIM(enactfiles(ji))
376	END DO
377	ENDIF
378	IF (ln_cor) THEN
379	DO ji = 1, jnumcorio
380	WRITE(numout,'(1X,2A)') ' Coriolis input observation file name coriofiles = ', &
381	TRIM(coriofiles(ji))
382	END DO
383	ENDIF
384	IF (ln_profb) THEN
385	DO ji = 1, jnumprofb
386	IF (ln_profb_ena(ji)) THEN
387	WRITE(numout,'(1X,2A)') ' Enact feedback input observation file name profbfiles = ', &
388	TRIM(profbfiles(ji))
389	ELSE
390	WRITE(numout,'(1X,2A)') ' Feedback input observation file name profbfiles = ', &
391	TRIM(profbfiles(ji))
392	ENDIF
393	WRITE(numout,'(1X,2A)') ' Enact feedback input time setting switch ln_profb_enatim = ', ln_profb_enatim(ji)
394	END DO
395	ENDIF
396	IF (ln_sladt) THEN
397	DO ji = 1, jnumslaact
398	WRITE(numout,'(1X,2A)') ' Active SLA input observation file name slafilesact = ', &
399	TRIM(slafilesact(ji))
400	END DO
401	DO ji = 1, jnumslapas
402	WRITE(numout,'(1X,2A)') ' Passive SLA input observation file name slafilespas = ', &
403	TRIM(slafilespas(ji))
404	END DO
405	ENDIF
406	IF (ln_slafb) THEN
407	DO ji = 1, jnumslafb
408	WRITE(numout,'(1X,2A)') ' Feedback SLA input observation file name slafbfiles = ', &
409	TRIM(slafbfiles(ji))
410	END DO
411	ENDIF
412	IF (ln_ghrsst) THEN
413	DO ji = 1, jnumsst
414	WRITE(numout,'(1X,2A)') ' GHRSST input observation file name sstfiles = ', &
415	TRIM(sstfiles(ji))
416	END DO
417	ENDIF
418	IF (ln_sstfb) THEN
419	DO ji = 1, jnumsstfb
420	WRITE(numout,'(1X,2A)') ' Feedback SST input observation file name sstfbfiles = ', &
421	TRIM(sstfbfiles(ji))
422	END DO
423	ENDIF
424	IF (ln_seaice) THEN
425	DO ji = 1, jnumseaice
426	WRITE(numout,'(1X,2A)') ' Sea Ice input observation file name seaicefiles = ', &
427	TRIM(seaicefiles(ji))
428	END DO
429	ENDIF
430	IF (ln_velavcur) THEN
431	DO ji = 1, jnumvelavcur
432	WRITE(numout,'(1X,2A)') ' Vel. cur. daily av. input file name velavcurfiles = ', &
433	TRIM(velavcurfiles(ji))
434	END DO
435	ENDIF
436	IF (ln_velhrcur) THEN
437	DO ji = 1, jnumvelhrcur
438	WRITE(numout,'(1X,2A)') ' Vel. cur. high freq. input file name velhvcurfiles = ', &
439	TRIM(velhrcurfiles(ji))
440	END DO
441	ENDIF
442	IF (ln_velavadcp) THEN
443	DO ji = 1, jnumvelavadcp
444	WRITE(numout,'(1X,2A)') ' Vel. ADCP daily av. input file name velavadcpfiles = ', &
445	TRIM(velavadcpfiles(ji))
446	END DO
447	ENDIF
448	IF (ln_velhradcp) THEN
449	DO ji = 1, jnumvelhradcp
450	WRITE(numout,'(1X,2A)') ' Vel. ADCP high freq. input file name velhvadcpfiles = ', &
451	TRIM(velhradcpfiles(ji))
452	END DO
453	ENDIF
454	IF (ln_velfb) THEN
455	DO ji = 1, jnumvelfb
456	IF (ln_velfb_av(ji)) THEN
457	WRITE(numout,'(1X,2A)') ' Vel. feedback daily av. input file name velfbfiles = ', &
458	TRIM(velfbfiles(ji))
459	ELSE
460	WRITE(numout,'(1X,2A)') ' Vel. feedback input observation file name velfbfiles = ', &
461	TRIM(velfbfiles(ji))
462	ENDIF
463	END DO
464	ENDIF
465	WRITE(numout,*) ' Initial date in window YYYYMMDD.HHMMSS dobsini = ', dobsini
466	WRITE(numout,*) ' Final date in window YYYYMMDD.HHMMSS dobsend = ', dobsend
467	WRITE(numout,*) ' Type of vertical interpolation method n1dint = ', n1dint
468	WRITE(numout,*) ' Type of horizontal interpolation method n2dint = ', n2dint
469	WRITE(numout,*) ' Rejection of observations near land swithch ln_nea = ', ln_nea
470	WRITE(numout,*) ' MSSH correction scheme nmsshc = ', nmsshc
471	WRITE(numout,*) ' MDT correction mdtcorr = ', mdtcorr
472	WRITE(numout,*) ' MDT cutoff for computed correction mdtcutoff = ', mdtcutoff
473	WRITE(numout,*) ' Logical switch for alt bias ln_altbias = ', ln_altbias
474	WRITE(numout,*) ' Logical switch for ignoring missing files ln_ignmis = ', ln_ignmis
475	WRITE(numout,*) ' ENACT daily average types = ',endailyavtypes
476
477	ENDIF
478
479	IF ( ln_vel3d .AND. ( .NOT. ln_grid_global ) ) THEN
480	CALL ctl_stop( 'Velocity data only works with ln_grid_global=.true.' )
481	RETURN
482	ENDIF
483
484	CALL obs_typ_init
485
486	CALL mppmap_init
487
488	! Parameter control
489	#if defined key_diaobs
490	IF ( ( .NOT. ln_t3d ).AND.( .NOT. ln_s3d ).AND.( .NOT. ln_sla ).AND. &
491	& ( .NOT. ln_vel3d ).AND. &
492	& ( .NOT. ln_ssh ).AND.( .NOT. ln_sst ).AND.( .NOT. ln_sss ).AND. &
493	& ( .NOT. ln_seaice ).AND.( .NOT. ln_vel3d ) ) THEN
494	IF(lwp) WRITE(numout,cform_war)
495	IF(lwp) WRITE(numout,*) ' key_diaobs is activated but logical flags', &
496	& ' ln_t3d, ln_s3d, ln_sla, ln_ssh, ln_sst, ln_sss, ln_seaice, ln_vel3d are all set to .FALSE.'
497	nwarn = nwarn + 1
498	ENDIF
499	#endif
500
501	CALL obs_grid_setup( )
502	IF ( ( n1dint < 0 ).OR.( n1dint > 1 ) ) THEN
503	IF(lwp) WRITE(numout,cform_err)
504	IF(lwp) WRITE(numout,*) ' Choice of vertical (1D) interpolation method', &
505	& ' is not available'
506	nstop = nstop + 1
507	ENDIF
508	IF ( ( n2dint < 0 ).OR.( n2dint > 4 ) ) THEN
509	IF(lwp) WRITE(numout,cform_err)
510	IF(lwp) WRITE(numout,*) ' Choice of horizontal (2D) interpolation method', &
511	& ' is not available'
512	nstop = nstop + 1
513	ENDIF
514
515	!-----------------------------------------------------------------------
516	! Depending on switches read the various observation types
517	!-----------------------------------------------------------------------
518	! - Temperature/salinity profiles
519
520	IF ( ln_t3d .OR. ln_s3d ) THEN
521
522	! Set the number of variables for profiles to 2 (T and S)
523	nprofvars = 2
524	! Set the number of extra variables for profiles to 1 (insitu temp).
525	nprofextr = 1
526
527	! Count how may insitu data sets we have and allocate data.
528	jprofset = 0
529	IF ( ln_ena ) jprofset = jprofset + 1
530	IF ( ln_cor ) jprofset = jprofset + 1
531	IF ( ln_profb ) jprofset = jprofset + jnumprofb
532	nprofsets = jprofset
533	IF ( nprofsets > 0 ) THEN
534	ALLOCATE(ld_enact(nprofsets))
535	ALLOCATE(profdata(nprofsets))
536	ALLOCATE(prodatqc(nprofsets))
537	ENDIF
538
539	jprofset = 0
540
541	! ENACT insitu data
542
543	IF ( ln_ena ) THEN
544
545	#if defined key_datetime_out
546	IF (lwp) THEN
547	CALL DATE_AND_TIME(datestr,timestr)
548	WRITE(numout,*) 'Read enact date_and_time ',datestr,' ',timestr
549	ENDIF
550	#endif
551	jprofset = jprofset + 1
552
553	ld_enact(jprofset) = .TRUE.
554
555	CALL obs_rea_pro_dri( 1, profdata(jprofset), &
556	& jnumenact, enactfiles(1:jnumenact), &
557	& nprofvars, nprofextr, &
558	& nitend-nit000+2, &
559	& dobsini, dobsend, ln_t3d, ln_s3d, &
560	& ln_ignmis, ln_s_at_t, .TRUE., .FALSE., &
561	& kdailyavtypes = endailyavtypes )
562
563	DO jvar = 1, 2
564
565	CALL obs_prof_staend( profdata(jprofset), jvar )
566
567	END DO
568
569	CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset), &
570	& ln_t3d, ln_s3d, ln_nea, &
571	& kdailyavtypes=endailyavtypes )
572
573	ENDIF
574
575	! Coriolis insitu data
576
577	IF ( ln_cor ) THEN
578
579	jprofset = jprofset + 1
580
581	ld_enact(jprofset) = .FALSE.
582
583	CALL obs_rea_pro_dri( 2, profdata(jprofset), &
584	& jnumcorio, coriofiles(1:jnumcorio), &
585	& nprofvars, nprofextr, &
586	& nitend-nit000+2, &
587	& dobsini, dobsend, ln_t3d, ln_s3d, &
588	& ln_ignmis, ln_s_at_t, .FALSE., .FALSE. )
589
590	DO jvar = 1, 2
591
592	CALL obs_prof_staend( profdata(jprofset), jvar )
593
594	END DO
595
596	CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset), &
597	& ln_t3d, ln_s3d, ln_nea )
598
599	ENDIF
600
601	! Feedback insitu data
602
603	IF ( ln_profb ) THEN
604
605	DO jset = 1, jnumprofb
606
607	jprofset = jprofset + 1
608	ld_enact (jprofset) = ln_profb_ena(jset)
609
610	CALL obs_rea_pro_dri( 0, profdata(jprofset), &
611	& 1, profbfiles(jset:jset), &
612	& nprofvars, nprofextr, &
613	& nitend-nit000+2, &
614	& dobsini, dobsend, ln_t3d, ln_s3d, &
615	& ln_ignmis, ln_s_at_t, &
616	& ld_enact(jprofset).AND.&
617	& ln_profb_enatim(jset), &
618	& .FALSE., kdailyavtypes = endailyavtypes )
619
620	DO jvar = 1, 2
621
622	CALL obs_prof_staend( profdata(jprofset), jvar )
623
624	END DO
625
626	IF ( ld_enact(jprofset) ) THEN
627	CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset), &
628	& ln_t3d, ln_s3d, ln_nea, &
629	& kdailyavtypes = endailyavtypes )
630	ELSE
631	CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset), &
632	& ln_t3d, ln_s3d, ln_nea )
633	ENDIF
634
635	END DO
636
637	ENDIF
638
639	ENDIF
640
641	! - Sea level anomalies
642	IF ( ln_sla ) THEN
643	#if defined key_datetime_out
644	IF (lwp) THEN
645	CALL DATE_AND_TIME(datestr,timestr)
646	WRITE(numout,*) 'Read SLA date_and_time ',datestr,' ',timestr
647	ENDIF
648	#endif
649	! Set the number of variables for sla to 1
650	nslavars = 1
651
652	! Set the number of extra variables for sla to 2
653	nslaextr = 2
654
655	! Set the number of sla data sets to 2
656	nslasets = 0
657	IF ( ln_sladt ) THEN
658	nslasets = nslasets + 2
659	ENDIF
660	IF ( ln_slafb ) THEN
661	nslasets = nslasets + jnumslafb
662	ENDIF
663
664	ALLOCATE(sladata(nslasets))
665	ALLOCATE(sladatqc(nslasets))
666	sladata(:)%nsurf=0
667	sladatqc(:)%nsurf=0
668
669	nslasets = 0
670
671	! AVISO SLA data
672
673	IF ( ln_sladt ) THEN
674
675	! Active SLA observations
676
677	nslasets = nslasets + 1
678
679	CALL obs_rea_sla( 1, sladata(nslasets), jnumslaact, &
680	& slafilesact(1:jnumslaact), &
681	& nslavars, nslaextr, nitend-nit000+2, &
682	& dobsini, dobsend, ln_ignmis, .FALSE. )
683	#if defined key_datetime_out
684	IF (lwp) THEN
685	CALL DATE_AND_TIME(datestr,timestr)
686	WRITE(numout,*) 'obs_pre_sla date_and_time ',datestr,' ',timestr
687	ENDIF
688	#endif
689	CALL obs_pre_sla( sladata(nslasets), sladatqc(nslasets), &
690	& ln_sla, ln_nea )
691
692	! Passive SLA observations
693
694	nslasets = nslasets + 1
695
696	#if defined key_datetime_out
697	IF (lwp) THEN
698	CALL DATE_AND_TIME(datestr,timestr)
699	WRITE(numout,*) 'obs_pre_sla date_and_time ',datestr,' ',timestr
700	ENDIF
701	#endif
702	CALL obs_rea_sla( 1, sladata(nslasets), jnumslapas, &
703	& slafilespas(1:jnumslapas), &
704	& nslavars, nslaextr, nitend-nit000+2, &
705	& dobsini, dobsend, ln_ignmis, .FALSE. )
706
707	#if defined key_datetime_out
708	IF (lwp) THEN
709	CALL DATE_AND_TIME(datestr,timestr)
710	WRITE(numout,*) 'obs_pre_sla date_and_time ',datestr,' ',timestr
711	ENDIF
712	#endif
713	CALL obs_pre_sla( sladata(nslasets), sladatqc(nslasets), &
714	& ln_sla, ln_nea )
715
716	ENDIF
717
718	! Feedback SLA data
719
720	IF ( ln_slafb ) THEN
721
722	DO jset = 1, jnumslafb
723
724	nslasets = nslasets + 1
725
726	CALL obs_rea_sla( 0, sladata(nslasets), 1, &
727	& slafbfiles(jset:jset), &
728	& nslavars, nslaextr, nitend-nit000+2, &
729	& dobsini, dobsend, ln_ignmis, .FALSE. )
730	#if defined key_datetime_out
731	IF (lwp) THEN
732	CALL DATE_AND_TIME(datestr,timestr)
733	WRITE(numout,*) 'obs_pre_sla date_and_time ',datestr,' ',timestr
734	ENDIF
735	#endif
736	CALL obs_pre_sla( sladata(nslasets), sladatqc(nslasets), &
737	& ln_sla, ln_nea )
738
739	END DO
740
741	ENDIF
742
743	#if defined key_datetime_out
744	IF (lwp) THEN
745	CALL DATE_AND_TIME(datestr,timestr)
746	WRITE(numout,*) 'obs_rea_mdt date_and_time ',datestr,' ',timestr
747	ENDIF
748	#endif
749	CALL obs_rea_mdt( nslasets, sladatqc, n2dint )
750
751	! read in altimeter bias
752
753	IF ( ln_altbias ) THEN
754	CALL obs_rea_altbias ( nslasets, sladatqc, n2dint, bias_file )
755	ENDIF
756
757	ENDIF
758
759	! - Sea surface height
760	IF ( ln_ssh ) THEN
761	IF(lwp) WRITE(numout,*) ' SSH currently not available'
762	ENDIF
763
764	! - Sea surface temperature
765	IF ( ln_sst ) THEN
766	#if defined key_datetime_out
767	IF (lwp) THEN
768	CALL DATE_AND_TIME(datestr,timestr)
769	WRITE(numout,*) 'Read SST date_and_time ',datestr,' ',timestr
770	ENDIF
771	#endif
772
773	! Set the number of variables for sst to 1
774	nsstvars = 1
775
776	! Set the number of extra variables for sst to 0
777	nsstextr = 0
778
779	nsstsets = 0
780
781	IF (ln_reysst) nsstsets = nsstsets + 1
782	IF (ln_ghrsst) nsstsets = nsstsets + 1
783	IF ( ln_sstfb ) THEN
784	nsstsets = nsstsets + jnumsstfb
785	ENDIF
786
787	ALLOCATE(sstdata(nsstsets))
788	ALLOCATE(sstdatqc(nsstsets))
789	sstdata(:)%nsurf=0
790	sstdatqc(:)%nsurf=0
791
792	nsstsets = 0
793
794	IF (ln_reysst) THEN
795
796	nsstsets = nsstsets + 1
797
798	CALL obs_rea_sst_rey( reysstname, reysstfmt, sstdata(nsstsets), &
799	& nsstvars, nsstextr, &
800	& nitend-nit000+2, dobsini, dobsend )
801	CALL obs_pre_sst( sstdata(nsstsets), sstdatqc(nsstsets), ln_sst, &
802	& ln_nea )
803
804	ENDIF
805
806	IF (ln_ghrsst) THEN
807
808	nsstsets = nsstsets + 1
809
810	CALL obs_rea_sst( 1, sstdata(nsstsets), jnumsst, &
811	& sstfiles(1:jnumsst), &
812	& nsstvars, nsstextr, nitend-nit000+2, &
813	& dobsini, dobsend, ln_ignmis, .FALSE. )
814	#if defined key_datetime_out
815	IF (lwp) THEN
816	CALL DATE_AND_TIME(datestr,timestr)
817	WRITE(numout,*) 'obs_pre_sst date_and_time ',datestr,' ',timestr
818	ENDIF
819	#endif
820	CALL obs_pre_sst( sstdata(nsstsets), sstdatqc(nsstsets), ln_sst, &
821	& ln_nea )
822
823	ENDIF
824
825	! Feedback SST data
826
827	IF ( ln_sstfb ) THEN
828
829	DO jset = 1, jnumsstfb
830
831	nsstsets = nsstsets + 1
832
833	CALL obs_rea_sst( 0, sstdata(nsstsets), 1, &
834	& sstfbfiles(jset:jset), &
835	& nsstvars, nsstextr, nitend-nit000+2, &
836	& dobsini, dobsend, ln_ignmis, .FALSE. )
837	#if defined key_datetime_out
838	IF (lwp) THEN
839	CALL DATE_AND_TIME(datestr,timestr)
840	WRITE(numout,*) 'obs_pre_sst date_and_time ',datestr,' ',timestr
841	ENDIF
842	#endif
843	CALL obs_pre_sst( sstdata(nsstsets), sstdatqc(nsstsets), &
844	& ln_sst, ln_nea )
845
846	END DO
847
848	ENDIF
849
850	ENDIF
851
852	! - Sea surface salinity
853	IF ( ln_sss ) THEN
854	IF(lwp) WRITE(numout,*) ' SSS currently not available'
855	ENDIF
856
857	! - Sea Ice Concentration
858
859	IF ( ln_seaice ) THEN
860
861	#if defined key_datetime_out
862	IF (lwp) THEN
863	CALL DATE_AND_TIME(datestr,timestr)
864	WRITE(numout,*) 'Read sea ice data date_and_time ',datestr,' ',timestr
865	ENDIF
866	#endif
867
868	! Set the number of variables for seaice to 1
869	nseaicevars = 1
870
871	! Set the number of extra variables for seaice to 0
872	nseaiceextr = 0
873
874	! Set the number of data sets to 1
875	nseaicesets = 1
876
877	ALLOCATE(seaicedata(nseaicesets))
878	ALLOCATE(seaicedatqc(nseaicesets))
879	seaicedata(:)%nsurf=0
880	seaicedatqc(:)%nsurf=0
881
882	CALL obs_rea_seaice( 1, seaicedata(nseaicesets), jnumseaice, &
883	& seaicefiles(1:jnumseaice), &
884	& nseaicevars, nseaiceextr, nitend-nit000+2, &
885	& dobsini, dobsend, ln_ignmis, .FALSE. )
886
887	#if defined key_datetime_out
888	IF (lwp) THEN
889	CALL DATE_AND_TIME(datestr,timestr)
890	WRITE(numout,*) 'obs_pre_seaice date_and_time ',datestr,' ',timestr
891	ENDIF
892	#endif
893
894	CALL obs_pre_seaice( seaicedata(nseaicesets), seaicedatqc(nseaicesets), &
895	& ln_seaice, ln_nea )
896
897	ENDIF
898
899	IF (ln_vel3d) THEN
900
901	! Set the number of variables for profiles to 2 (U and V)
902	nvelovars = 2
903
904	! Set the number of extra variables for profiles to 2 to store
905	! rotation parameters
906	nveloextr = 2
907
908	jveloset = 0
909
910	IF ( ln_velavcur ) jveloset = jveloset + 1
911	IF ( ln_velhrcur ) jveloset = jveloset + 1
912	IF ( ln_velavadcp ) jveloset = jveloset + 1
913	IF ( ln_velhradcp ) jveloset = jveloset + 1
914	IF (ln_velfb) jveloset = jveloset + jnumvelfb
915
916	nvelosets = jveloset
917	IF ( nvelosets > 0 ) THEN
918	ALLOCATE( velodata(nvelosets) )
919	ALLOCATE( veldatqc(nvelosets) )
920	ALLOCATE( ld_velav(nvelosets) )
921	ENDIF
922
923	jveloset = 0
924
925	! Daily averaged data
926
927	IF ( ln_velavcur ) THEN
928
929	jveloset = jveloset + 1
930
931	ld_velav(jveloset) = .TRUE.
932
933	CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelavcur, &
934	& velavcurfiles(1:jnumvelavcur), &
935	& nvelovars, nveloextr, &
936	& nitend-nit000+2, &
937	& dobsini, dobsend, ln_ignmis, &
938	& ld_velav(jveloset), &
939	& .FALSE. )
940
941	DO jvar = 1, 2
942	CALL obs_prof_staend( velodata(jveloset), jvar )
943	END DO
944
945	CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), &
946	& ln_vel3d, ln_nea, ld_velav(jveloset) )
947
948	ENDIF
949
950	! High frequency data
951
952	IF ( ln_velhrcur ) THEN
953
954	jveloset = jveloset + 1
955
956	ld_velav(jveloset) = .FALSE.
957
958	CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelhrcur, &
959	& velhrcurfiles(1:jnumvelhrcur), &
960	& nvelovars, nveloextr, &
961	& nitend-nit000+2, &
962	& dobsini, dobsend, ln_ignmis, &
963	& ld_velav(jveloset), &
964	& .FALSE. )
965
966	DO jvar = 1, 2
967	CALL obs_prof_staend( velodata(jveloset), jvar )
968	END DO
969
970	CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), &
971	& ln_vel3d, ln_nea, ld_velav(jveloset) )
972
973	ENDIF
974
975	! Daily averaged data
976
977	IF ( ln_velavadcp ) THEN
978
979	jveloset = jveloset + 1
980
981	ld_velav(jveloset) = .TRUE.
982
983	CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelavadcp, &
984	& velavadcpfiles(1:jnumvelavadcp), &
985	& nvelovars, nveloextr, &
986	& nitend-nit000+2, &
987	& dobsini, dobsend, ln_ignmis, &
988	& ld_velav(jveloset), &
989	& .FALSE. )
990
991	DO jvar = 1, 2
992	CALL obs_prof_staend( velodata(jveloset), jvar )
993	END DO
994
995	CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), &
996	& ln_vel3d, ln_nea, ld_velav(jveloset) )
997
998	ENDIF
999
1000	! High frequency data
1001
1002	IF ( ln_velhradcp ) THEN
1003
1004	jveloset = jveloset + 1
1005
1006	ld_velav(jveloset) = .FALSE.
1007
1008	CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelhradcp, &
1009	& velhradcpfiles(1:jnumvelhradcp), &
1010	& nvelovars, nveloextr, &
1011	& nitend-nit000+2, &
1012	& dobsini, dobsend, ln_ignmis, &
1013	& ld_velav(jveloset), &
1014	& .FALSE. )
1015
1016	DO jvar = 1, 2
1017	CALL obs_prof_staend( velodata(jveloset), jvar )
1018	END DO
1019
1020	CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), &
1021	& ln_vel3d, ln_nea, ld_velav(jveloset) )
1022
1023	ENDIF
1024
1025	IF ( ln_velfb ) THEN
1026
1027	DO jset = 1, jnumvelfb
1028
1029	jveloset = jveloset + 1
1030
1031	ld_velav(jveloset) = ln_velfb_av(jset)
1032
1033	CALL obs_rea_vel_dri( 0, velodata(jveloset), 1, &
1034	& velfbfiles(jset:jset), &
1035	& nvelovars, nveloextr, &
1036	& nitend-nit000+2, &
1037	& dobsini, dobsend, ln_ignmis, &
1038	& ld_velav(jveloset), &
1039	& .FALSE. )
1040
1041	DO jvar = 1, 2
1042	CALL obs_prof_staend( velodata(jveloset), jvar )
1043	END DO
1044
1045	CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), &
1046	& ln_vel3d, ln_nea, ld_velav(jveloset) )
1047
1048
1049	END DO
1050
1051	ENDIF
1052
1053	ENDIF
1054
1055	#if defined key_datetime_out
1056	IF (lwp) THEN
1057	CALL DATE_AND_TIME(datestr,timestr)
1058	WRITE(numout,*) 'end dia_obs_init date_and_time ',datestr,' ',timestr
1059	ENDIF
1060	#endif
1061
1062	END SUBROUTINE dia_obs_init
1063
1064	SUBROUTINE dia_obs( kstp )
1065	!!----------------------------------------------------------------------
1066	!! * ROUTINE dia_obs *
1067	!!
1068	!! ** Purpose : Call the observation operators on each time step
1069	!!
1070	!! ** Method : Call the observation operators on each time step to
1071	!! compute the model equivalent of the following date:
1072	!! - T profiles
1073	!! - S profiles
1074	!! - Sea surface height (referenced to a mean)
1075	!! - Sea surface temperature
1076	!! - Sea surface salinity
1077	!! - Velocity component (U,V) profiles
1078	!!
1079	!! ** Action :
1080	!!
1081	!! History :
1082	!! ! 06-03 (K. Mogensen) Original code
1083	!! ! 06-05 (K. Mogensen) Reformatted
1084	!! ! 06-10 (A. Weaver) Cleaning
1085	!! ! 07-03 (K. Mogensen) General handling of profiles
1086	!! ! 07-04 (G. Smith) Generalized surface operators
1087	!! ! 08-10 (M. Valdivieso) obs operator for velocity profiles
1088	!!----------------------------------------------------------------------
1089	!! * Modules used
1090	USE dom_oce, ONLY : & ! Ocean space and time domain variables
1091	& rdt, &
1092	& gdept_0, &
1093	& tmask, umask, vmask
1094	USE phycst, ONLY : & ! Physical constants
1095	& rday
1096	USE oce, ONLY : & ! Ocean dynamics and tracers variables
1097	& tn, &
1098	& sn, &
1099	& un, vn, &
1100	& sshn
1101	#if defined key_ice_lim
1102	USE ice, ONLY : & ! LIM Ice model variables
1103	& frld
1104	#endif
1105
1106	IMPLICIT NONE
1107
1108	!! * Arguments
1109	INTEGER, INTENT(IN) :: kstp ! Current timestep
1110	!! * Local declarations
1111	#if ! defined key_ice_lim
1112	REAL(wp), DIMENSION(jpi,jpj) :: frld
1113	#endif
1114	INTEGER :: idaystp ! Number of timesteps per day
1115	INTEGER :: jprofset ! Profile data set loop variable
1116	INTEGER :: jslaset ! SLA data set loop variable
1117	INTEGER :: jsstset ! SST data set loop variable
1118	INTEGER :: jseaiceset ! sea ice data set loop variable
1119	INTEGER :: jveloset ! velocity profile data loop variable
1120	INTEGER :: jvar ! Variable number
1121	CHARACTER(LEN=20) :: datestr=" ",timestr=" "
1122
1123	IF(lwp) THEN
1124	WRITE(numout,*)
1125	WRITE(numout,*) 'dia_obs : Call the observation operators', kstp
1126	WRITE(numout,*) '~~~~~~~'
1127	#if defined key_datetime_out
1128	CALL DATE_AND_TIME(datestr,timestr)
1129	WRITE(numout,*) 'dia_obs date_and_time ',datestr,' ',timestr
1130	#endif
1131	ENDIF
1132
1133	idaystp = NINT( rday / rdt )
1134
1135	!-----------------------------------------------------------------------
1136	! No LIM => frld == 0.0_wp
1137	!-----------------------------------------------------------------------
1138	#if ! defined key_ice_lim
1139	frld(:,:) = 0.0_wp
1140	#endif
1141	!-----------------------------------------------------------------------
1142	! Depending on switches call various observation operators
1143	!-----------------------------------------------------------------------
1144
1145	! - Temperature/salinity profiles
1146	IF ( ln_t3d .OR. ln_s3d ) THEN
1147	DO jprofset = 1, nprofsets
1148	IF ( ld_enact(jprofset) ) THEN
1149	CALL obs_pro_opt( prodatqc(jprofset), &
1150	& kstp, jpi, jpj, jpk, nit000, idaystp, tn, sn,&
1151	& gdept_0, tmask, n1dint, n2dint, &
1152	& kdailyavtypes = endailyavtypes )
1153	ELSE
1154	CALL obs_pro_opt( prodatqc(jprofset), &
1155	& kstp, jpi, jpj, jpk, nit000, idaystp, tn, sn,&
1156	& gdept_0, tmask, n1dint, n2dint )
1157	ENDIF
1158	END DO
1159	ENDIF
1160
1161	! - Sea surface anomaly
1162	IF ( ln_sla ) THEN
1163	DO jslaset = 1, nslasets
1164	CALL obs_sla_opt( sladatqc(jslaset), &
1165	& kstp, jpi, jpj, nit000, sshn, &
1166	& tmask(:,:,1), n2dint )
1167	END DO
1168	ENDIF
1169
1170	! - Sea surface temperature
1171	IF ( ln_sst ) THEN
1172	DO jsstset = 1, nsstsets
1173	CALL obs_sst_opt( sstdatqc(jsstset), &
1174	& kstp, jpi, jpj, nit000, tn(:,:,1), &
1175	& tmask(:,:,1), n2dint )
1176	END DO
1177	ENDIF
1178
1179	! - Sea surface salinity
1180	IF ( ln_sss ) THEN
1181	IF(lwp) WRITE(numout,*) ' SSS currently not available'
1182	ENDIF
1183
1184	#if defined key_ice_lim
1185	IF ( ln_seaice ) THEN
1186	DO jseaiceset = 1, nseaicesets
1187	CALL obs_seaice_opt( seaicedatqc(jseaiceset), &
1188	& kstp, jpi, jpj, nit000, 1.-frld, &
1189	& tmask(:,:,1), n2dint )
1190	END DO
1191	ENDIF
1192	#endif
1193
1194	! - Velocity profiles
1195	IF ( ln_vel3d ) THEN
1196	DO jveloset = 1, nvelosets
1197	! zonal component of velocity
1198	CALL obs_vel_opt( veldatqc(jveloset), kstp, jpi, jpj, jpk, &
1199	& nit000, idaystp, un, vn, gdept_0, umask, vmask, &
1200	n1dint, n2dint, ld_velav(jveloset) )
1201	END DO
1202	ENDIF
1203
1204	#if defined key_datetime_out
1205	IF(lwp) THEN
1206	CALL DATE_AND_TIME(datestr,timestr)
1207	WRITE(numout,*) 'dia_obs finished date_and_time ',datestr,' ',timestr
1208	ENDIF
1209	#endif
1210
1211	END SUBROUTINE dia_obs
1212
1213	SUBROUTINE dia_obs_wri
1214	!!----------------------------------------------------------------------
1215	!! * ROUTINE dia_obs_wri *
1216	!!
1217	!! ** Purpose : Call observation diagnostic output routines
1218	!!
1219	!! ** Method : Call observation diagnostic output routines
1220	!!
1221	!! ** Action :
1222	!!
1223	!! History :
1224	!! ! 06-03 (K. Mogensen) Original code
1225	!! ! 06-05 (K. Mogensen) Reformatted
1226	!! ! 06-10 (A. Weaver) Cleaning
1227	!! ! 07-03 (K. Mogensen) General handling of profiles
1228	!! ! 08-09 (M. Valdivieso) Velocity component (U,V) profiles
1229	!!----------------------------------------------------------------------
1230	IMPLICIT NONE
1231
1232	!! * Local declarations
1233
1234	INTEGER :: jprofset ! Profile data set loop variable
1235	INTEGER :: jveloset ! Velocity data set loop variable
1236	INTEGER :: jslaset ! SLA data set loop variable
1237	INTEGER :: jsstset ! SST data set loop variable
1238	INTEGER :: jseaiceset ! Sea Ice data set loop variable
1239	INTEGER :: jset
1240	INTEGER :: jfbini
1241	CHARACTER(LEN=20) :: datestr=" ",timestr=" "
1242	CHARACTER(LEN=10) :: cdtmp
1243	!-----------------------------------------------------------------------
1244	! Depending on switches call various observation output routines
1245	!-----------------------------------------------------------------------
1246
1247	#if defined key_datetime_out
1248	IF (lwp) THEN
1249	CALL DATE_AND_TIME(datestr,timestr)
1250	WRITE(numout,*) 'dia_obs_wri date_and_time ',datestr,' ',timestr
1251	ENDIF
1252	#endif
1253
1254	! - Temperature/salinity profiles
1255
1256	IF( ln_t3d .OR. ln_s3d ) THEN
1257
1258	! Copy data from prodatqc to profdata structures
1259	DO jprofset = 1, nprofsets
1260
1261	CALL obs_prof_decompress( prodatqc(jprofset), &
1262	& profdata(jprofset), .TRUE., numout )
1263
1264	END DO
1265
1266	! Write the profiles.
1267
1268	jprofset = 0
1269
1270
1271	! ENACT insitu data
1272
1273	IF ( ln_ena ) THEN
1274
1275	jprofset = jprofset + 1
1276
1277	CALL obs_wri_p3d( 'enact', profdata(jprofset) )
1278
1279	ENDIF
1280
1281	! Coriolis insitu data
1282
1283	IF ( ln_cor ) THEN
1284
1285	jprofset = jprofset + 1
1286
1287	CALL obs_wri_p3d( 'corio', profdata(jprofset) )
1288
1289	ENDIF
1290
1291	! Feedback insitu data
1292
1293	IF ( ln_profb ) THEN
1294
1295	jfbini = jprofset + 1
1296
1297	DO jprofset = jfbini, nprofsets
1298
1299	jset = jprofset - jfbini + 1
1300	WRITE(cdtmp,'(A,I2.2)')'profb_',jset
1301	CALL obs_wri_p3d( cdtmp, profdata(jprofset) )
1302
1303	END DO
1304
1305	ENDIF
1306
1307	ENDIF
1308
1309	! - Sea surface anomaly
1310	IF ( ln_sla ) THEN
1311
1312	! Copy data from sladatqc to sladata structures
1313	DO jslaset = 1, nslasets
1314
1315	CALL obs_surf_decompress( sladatqc(jslaset), &
1316	& sladata(jslaset), .TRUE., numout )
1317
1318	END DO
1319
1320	jslaset = 0
1321
1322	! Write the AVISO SLA data
1323
1324	IF ( ln_sladt ) THEN
1325
1326	jslaset = 1
1327	CALL obs_wri_sla( 'aviso_act', sladata(jslaset) )
1328	jslaset = 2
1329	CALL obs_wri_sla( 'aviso_pas', sladata(jslaset) )
1330
1331	ENDIF
1332
1333	IF ( ln_slafb ) THEN
1334
1335	jfbini = jslaset + 1
1336
1337	DO jslaset = jfbini, nslasets
1338
1339	jset = jslaset - jfbini + 1
1340	WRITE(cdtmp,'(A,I2.2)')'slafb_',jset
1341	CALL obs_wri_sla( cdtmp, sladata(jslaset) )
1342
1343	END DO
1344
1345	ENDIF
1346
1347	ENDIF
1348
1349	! - Sea surface temperature
1350	IF ( ln_sst ) THEN
1351
1352	! Copy data from sstdatqc to sstdata structures
1353	DO jsstset = 1, nsstsets
1354
1355	CALL obs_surf_decompress( sstdatqc(jsstset), &
1356	& sstdata(jsstset), .TRUE., numout )
1357
1358	END DO
1359
1360	jsstset = 0
1361
1362	! Write the AVISO SST data
1363
1364	IF ( ln_reysst ) THEN
1365
1366	jsstset = jsstset + 1
1367	CALL obs_wri_sst( 'reynolds', sstdata(jsstset) )
1368
1369	ENDIF
1370
1371	IF ( ln_ghrsst ) THEN
1372
1373	jsstset = jsstset + 1
1374	CALL obs_wri_sst( 'ghr', sstdata(jsstset) )
1375
1376	ENDIF
1377
1378	IF ( ln_sstfb ) THEN
1379
1380	jfbini = jsstset + 1
1381
1382	DO jsstset = jfbini, nsstsets
1383
1384	jset = jsstset - jfbini + 1
1385	WRITE(cdtmp,'(A,I2.2)')'sstfb_',jset
1386	CALL obs_wri_sst( cdtmp, sstdata(jsstset) )
1387
1388	END DO
1389
1390	ENDIF
1391
1392	ENDIF
1393
1394	! - Sea surface salinity
1395	IF ( ln_sss ) THEN
1396	IF(lwp) WRITE(numout,*) ' SSS currently not available'
1397	ENDIF
1398
1399	! - Sea Ice Concentration
1400	IF ( ln_seaice ) THEN
1401
1402	! Copy data from seaicedatqc to seaicedata structures
1403	DO jseaiceset = 1, nseaicesets
1404
1405	CALL obs_surf_decompress( seaicedatqc(jseaiceset), &
1406	& seaicedata(jseaiceset), .TRUE., numout )
1407
1408	END DO
1409
1410	! Write the Sea Ice data
1411	DO jseaiceset = 1, nseaicesets
1412
1413	WRITE(cdtmp,'(A,I2.2)')'seaicefb_',jseaiceset
1414	CALL obs_wri_seaice( cdtmp, seaicedata(jseaiceset) )
1415
1416	END DO
1417
1418	ENDIF
1419
1420	! Velocity data
1421	IF( ln_vel3d ) THEN
1422
1423	! Copy data from veldatqc to velodata structures
1424	DO jveloset = 1, nvelosets
1425
1426	CALL obs_prof_decompress( veldatqc(jveloset), &
1427	& velodata(jveloset), .TRUE., numout )
1428
1429	END DO
1430
1431	! Write the profiles.
1432
1433	jveloset = 0
1434
1435	! Daily averaged data
1436
1437	IF ( ln_velavcur ) THEN
1438
1439	jveloset = jveloset + 1
1440
1441	CALL obs_wri_vel( 'velavcurr', velodata(jveloset), n2dint )
1442
1443	ENDIF
1444
1445	! High frequency data
1446
1447	IF ( ln_velhrcur ) THEN
1448
1449	jveloset = jveloset + 1
1450
1451	CALL obs_wri_vel( 'velhrcurr', velodata(jveloset), n2dint )
1452
1453	ENDIF
1454
1455	! Daily averaged data
1456
1457	IF ( ln_velavadcp ) THEN
1458
1459	jveloset = jveloset + 1
1460
1461	CALL obs_wri_vel( 'velavadcp', velodata(jveloset), n2dint )
1462
1463	ENDIF
1464
1465	! High frequency data
1466
1467	IF ( ln_velhradcp ) THEN
1468
1469	jveloset = jveloset + 1
1470
1471	CALL obs_wri_vel( 'velhradcp', velodata(jveloset), n2dint )
1472
1473	ENDIF
1474
1475	! Feedback velocity data
1476
1477	IF ( ln_velfb ) THEN
1478
1479	jfbini = jveloset + 1
1480
1481	DO jveloset = jfbini, nvelosets
1482
1483	jset = jveloset - jfbini + 1
1484	WRITE(cdtmp,'(A,I2.2)')'velfb_',jset
1485	CALL obs_wri_vel( cdtmp, velodata(jveloset), n2dint )
1486
1487	END DO
1488
1489	ENDIF
1490
1491	ENDIF
1492
1493	#if defined key_datetime_out
1494	IF (lwp) THEN
1495	CALL DATE_AND_TIME(datestr,timestr)
1496	WRITE(numout,*) 'dia_obs_wri end date_and_time ',datestr,' ',timestr
1497	ENDIF
1498	#endif
1499
1500	END SUBROUTINE dia_obs_wri
1501
1502	SUBROUTINE ini_date( ddobsini )
1503	!!----------------------------------------------------------------------
1504	!! * ROUTINE ini_date *
1505	!!
1506	!! ** Purpose : Get initial data in double precision YYYYMMDD.HHMMSS format
1507	!!
1508	!! ** Method : Get initial data in double precision YYYYMMDD.HHMMSS format
1509	!!
1510	!! ** Action : Get initial data in double precision YYYYMMDD.HHMMSS format
1511	!!
1512	!! History :
1513	!! ! 06-03 (K. Mogensen) Original code
1514	!! ! 06-05 (K. Mogensen) Reformatted
1515	!! ! 06-10 (A. Weaver) Cleaning
1516	!! ! 06-10 (G. Smith) Calculates initial date the same as method for final date
1517	!! ! 10-05 (D. Lea) Update to month length calculation for NEMO vn3.2
1518	!!----------------------------------------------------------------------
1519	USE phycst, ONLY : & ! Physical constants
1520	& rday
1521	! USE daymod, ONLY : & ! Time variables
1522	! & nmonth_len
1523	USE dom_oce, ONLY : & ! Ocean space and time domain variables
1524	& rdt
1525
1526	IMPLICIT NONE
1527
1528	!! * Arguments
1529	REAL(KIND=dp), INTENT(OUT) :: ddobsini ! Initial date in YYYYMMDD.HHMMSS
1530
1531	!! * Local declarations
1532	INTEGER :: iyea ! date - (year, month, day, hour, minute)
1533	INTEGER :: imon
1534	INTEGER :: iday
1535	INTEGER :: ihou
1536	INTEGER :: imin
1537	INTEGER :: imday ! Number of days in month.
1538	REAL(KIND=wp) :: zdayfrc ! Fraction of day
1539
1540	INTEGER, DIMENSION(12) :: imonth_len !: length in days of the months of the current year
1541
1542	!!----------------------------------------------------------------------
1543	!! Initial date initialization (year, month, day, hour, minute)
1544	!! (This assumes that the initial date is for 00z))
1545	!!----------------------------------------------------------------------
1546	iyea = ndate0 / 10000
1547	imon = ( ndate0 - iyea * 10000 ) / 100
1548	iday = ndate0 - iyea * 10000 - imon * 100
1549	ihou = 0
1550	imin = 0
1551
1552	!!----------------------------------------------------------------------
1553	!! Compute number of days + number of hours + min since initial time
1554	!!----------------------------------------------------------------------
1555	iday = iday + ( nit000 -1 ) * rdt / rday
1556	zdayfrc = ( nit000 -1 ) * rdt / rday
1557	zdayfrc = zdayfrc - aint(zdayfrc)
1558	ihou = int( zdayfrc * 24 )
1559	imin = int( (zdayfrc * 24 - ihou) * 60 )
1560
1561	!!-----------------------------------------------------------------------
1562	!! Convert number of days (iday) into a real date
1563	!!----------------------------------------------------------------------
1564
1565	CALL calc_month_len( iyea, imonth_len )
1566
1567	DO WHILE ( iday > imonth_len(imon) )
1568	iday = iday - imonth_len(imon)
1569	imon = imon + 1
1570	IF ( imon > 12 ) THEN
1571	imon = 1
1572	iyea = iyea + 1
1573	CALL calc_month_len( iyea, imonth_len ) ! update month lengths
1574	ENDIF
1575	END DO
1576
1577	!!----------------------------------------------------------------------
1578	!! Convert it into YYYYMMDD.HHMMSS format.
1579	!!----------------------------------------------------------------------
1580	ddobsini = iyea * 10000_dp + imon * 100_dp + &
1581	& iday + ihou * 0.01_dp + imin * 0.0001_dp
1582
1583
1584	END SUBROUTINE ini_date
1585
1586	SUBROUTINE fin_date( ddobsfin )
1587	!!----------------------------------------------------------------------
1588	!! * ROUTINE fin_date *
1589	!!
1590	!! ** Purpose : Get final data in double precision YYYYMMDD.HHMMSS format
1591	!!
1592	!! ** Method : Get final data in double precision YYYYMMDD.HHMMSS format
1593	!!
1594	!! ** Action : Get final data in double precision YYYYMMDD.HHMMSS format
1595	!!
1596	!! History :
1597	!! ! 06-03 (K. Mogensen) Original code
1598	!! ! 06-05 (K. Mogensen) Reformatted
1599	!! ! 06-10 (A. Weaver) Cleaning
1600	!! ! 10-05 (D. Lea) Update to month length calculation for NEMO vn3.2
1601	!!----------------------------------------------------------------------
1602	USE phycst, ONLY : & ! Physical constants
1603	& rday
1604	! USE daymod, ONLY : & ! Time variables
1605	! & nmonth_len
1606	USE dom_oce, ONLY : & ! Ocean space and time domain variables
1607	& rdt
1608
1609	IMPLICIT NONE
1610
1611	!! * Arguments
1612	REAL(KIND=dp), INTENT(OUT) :: ddobsfin ! Final date in YYYYMMDD.HHMMSS
1613
1614	!! * Local declarations
1615	INTEGER :: iyea ! date - (year, month, day, hour, minute)
1616	INTEGER :: imon
1617	INTEGER :: iday
1618	INTEGER :: ihou
1619	INTEGER :: imin
1620	INTEGER :: imday ! Number of days in month.
1621	REAL(KIND=wp) :: zdayfrc ! Fraction of day
1622
1623	INTEGER, DIMENSION(12) :: imonth_len !: length in days of the months of the current year
1624
1625	!-----------------------------------------------------------------------
1626	! Initial date initialization (year, month, day, hour, minute)
1627	! (This assumes that the initial date is for 00z)
1628	!-----------------------------------------------------------------------
1629	iyea = ndate0 / 10000
1630	imon = ( ndate0 - iyea * 10000 ) / 100
1631	iday = ndate0 - iyea * 10000 - imon * 100
1632	ihou = 0
1633	imin = 0
1634
1635	!-----------------------------------------------------------------------
1636	! Compute number of days + number of hours + min since initial time
1637	!-----------------------------------------------------------------------
1638	iday = iday + nitend * rdt / rday
1639	zdayfrc = nitend * rdt / rday
1640	zdayfrc = zdayfrc - AINT( zdayfrc )
1641	ihou = INT( zdayfrc * 24 )
1642	imin = INT( ( zdayfrc * 24 - ihou ) * 60 )
1643
1644	!-----------------------------------------------------------------------
1645	! Convert number of days (iday) into a real date
1646	!----------------------------------------------------------------------
1647
1648	CALL calc_month_len( iyea, imonth_len )
1649
1650	DO WHILE ( iday > imonth_len(imon) )
1651	iday = iday - imonth_len(imon)
1652	imon = imon + 1
1653	IF ( imon > 12 ) THEN
1654	imon = 1
1655	iyea = iyea + 1
1656	CALL calc_month_len( iyea, imonth_len ) ! update month lengths
1657	ENDIF
1658	END DO
1659
1660	!-----------------------------------------------------------------------
1661	! Convert it into YYYYMMDD.HHMMSS format
1662	!-----------------------------------------------------------------------
1663	ddobsfin = iyea * 10000_dp + imon * 100_dp + iday &
1664	& + ihou * 0.01_dp + imin * 0.0001_dp
1665
1666	END SUBROUTINE fin_date
1667
1668	END MODULE diaobs

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

Download in other formats: