Context Navigation

diaobs.F90 @ 2576

Last change on this file since 2576 was 2474, checked in by djlea, 13 years ago
Update to OBS and ASM documentation. Removal of cpp key options to OBS code. Also moved the diaobs call to after the timestep code in step.
Property svn:keywords set to `Id`
File size: 58.1 KB

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

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

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

Context Navigation

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

Download in other formats: