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

Context Navigation

source: branches/2011/dev_r2787_LOCEAN3_TRA_TRP/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 9260

Last change on this file since 9260 was 2789, checked in by cetlod, 13 years ago
Implementation of the merge of TRA/TRP : first guess, see ticket #842
Property svn:keywords set to `Id`
File size: 58.3 KB

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

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

Download in other formats: