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

Context Navigation

source: branches/2013/dev_MERCATOR_UKMO_2013/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 9088

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

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

Download in other formats: