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

Context Navigation

source: trunk/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 2715

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

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

Download in other formats: