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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

diaobs.F90 in branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

Context Navigation

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

Last change on this file since 2382 was 2358, checked in by rblod, 14 years ago
Changes to be able to compile v3_3_beta with key_agrif,see ticket #753 ; just compilation fixes, I was to scared to try to run AGRIF
Property svn:keywords set to `Id`
File size: 61.8 KB

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

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

Download in other formats: