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

Context Navigation

source: branches/2012/dev_r3342_MERCATOR7_SST/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 3585

Last change on this file since 3585 was 3585, checked in by cbricaud, 11 years ago
commit modifications due to review for SST night operator
Property svn:keywords set to `Id`
File size: 59.7 KB

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

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

Download in other formats: