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 @ 2281

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

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

Download in other formats: