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

Context Navigation

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 2636

Last change on this file since 2636 was 2633, checked in by trackstand2, 13 years ago
Renamed wrk_use => wrk_in_use and wrk_release => wrk_not_released
Property svn:keywords set to `Id`
File size: 58.5 KB

Line
1	MODULE diaobs
2	!!======================================================================
3	!! * MODULE diaobs *
4	!! Observation diagnostics: Computation of the misfit between data and
5	!! their model equivalent
6	!!======================================================================
7
8	!!----------------------------------------------------------------------
9	!! 'key_diaobs' : Switch on the observation diagnostic computation
10	!!----------------------------------------------------------------------
11	!! dia_obs_init : Reading and prepare observations
12	!! dia_obs : Compute model equivalent to observations
13	!! dia_obs_wri : Write observational diagnostics
14	!! ini_date : Compute the initial date YYYYMMDD.HHMMSS
15	!! fin_date : Compute the final date YYYYMMDD.HHMMSS
16	!!----------------------------------------------------------------------
17	!! * Modules used
18	USE par_kind ! Precision variables
19	USE in_out_manager ! I/O manager
20	USE par_oce
21	USE dom_oce ! Ocean space and time domain variables
22	USE obs_read_prof ! Reading and allocation of observations (Coriolis)
23	USE obs_read_sla ! Reading and allocation of SLA observations
24	USE obs_read_sst ! Reading and allocation of SST observations
25	USE obs_readmdt ! Reading and allocation of MDT for SLA.
26	USE obs_read_seaice ! Reading and allocation of Sea Ice observations
27	USE obs_read_vel ! Reading and allocation of velocity component observations
28	USE obs_prep ! Preparation of obs. (grid search etc).
29	USE obs_oper ! Observation operators
30	USE obs_write ! Writing of observation related diagnostics
31	USE obs_grid ! Grid searching
32	USE obs_read_altbias ! Bias treatment for altimeter
33	USE obs_profiles_def ! Profile data definitions
34	USE obs_profiles ! Profile data storage
35	USE obs_surf_def ! Surface data definitions
36	USE obs_sla ! SLA data storage
37	USE obs_sst ! SST data storage
38	USE obs_seaice ! Sea Ice data storage
39	USE obs_types ! Definitions for observation types
40	USE mpp_map ! MPP mapping
41
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	USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
1029	#if ! defined key_ice_lim
1030	USE wrk_nemo, ONLY: frld => wrk_2d_1
1031	#endif
1032	IMPLICIT NONE
1033
1034	!! * Arguments
1035	INTEGER, INTENT(IN) :: kstp ! Current timestep
1036	!! * Local declarations
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 ! defined key_ice_lim
1047	IF(wrk_in_use(2, 1))THEN
1048	CALL ctl_stop('dia_obs : requested workspace array unavailable.')
1049	RETURN
1050	END IF
1051	#endif
1052
1053	IF(lwp) THEN
1054	WRITE(numout,*)
1055	WRITE(numout,*) 'dia_obs : Call the observation operators', kstp
1056	WRITE(numout,*) '~~~~~~~'
1057	ENDIF
1058
1059	idaystp = NINT( rday / rdt )
1060
1061	!-----------------------------------------------------------------------
1062	! No LIM => frld == 0.0_wp
1063	!-----------------------------------------------------------------------
1064	#if ! defined key_ice_lim
1065	frld(:,:) = 0.0_wp
1066	#endif
1067	!-----------------------------------------------------------------------
1068	! Depending on switches call various observation operators
1069	!-----------------------------------------------------------------------
1070
1071	! - Temperature/salinity profiles
1072	IF ( ln_t3d .OR. ln_s3d ) THEN
1073	DO jprofset = 1, nprofsets
1074	IF ( ld_enact(jprofset) ) THEN
1075	CALL obs_pro_opt( prodatqc(jprofset), &
1076	& kstp, jpi, jpj, jpk, nit000, idaystp, tn, sn,&
1077	& gdept_0, tmask, n1dint, n2dint, &
1078	& kdailyavtypes = endailyavtypes )
1079	ELSE
1080	CALL obs_pro_opt( prodatqc(jprofset), &
1081	& kstp, jpi, jpj, jpk, nit000, idaystp, tn, sn,&
1082	& gdept_0, tmask, n1dint, n2dint )
1083	ENDIF
1084	END DO
1085	ENDIF
1086
1087	! - Sea surface anomaly
1088	IF ( ln_sla ) THEN
1089	DO jslaset = 1, nslasets
1090	CALL obs_sla_opt( sladatqc(jslaset), &
1091	& kstp, jpi, jpj, nit000, sshn, &
1092	& tmask(:,:,1), n2dint )
1093	END DO
1094	ENDIF
1095
1096	! - Sea surface temperature
1097	IF ( ln_sst ) THEN
1098	DO jsstset = 1, nsstsets
1099	CALL obs_sst_opt( sstdatqc(jsstset), &
1100	& kstp, jpi, jpj, nit000, tn(:,:,1), &
1101	& tmask(:,:,1), n2dint )
1102	END DO
1103	ENDIF
1104
1105	! - Sea surface salinity
1106	IF ( ln_sss ) THEN
1107	IF(lwp) WRITE(numout,*) ' SSS currently not available'
1108	ENDIF
1109
1110	#if defined key_ice_lim
1111	IF ( ln_seaice ) THEN
1112	DO jseaiceset = 1, nseaicesets
1113	CALL obs_seaice_opt( seaicedatqc(jseaiceset), &
1114	& kstp, jpi, jpj, nit000, 1.-frld, &
1115	& tmask(:,:,1), n2dint )
1116	END DO
1117	ENDIF
1118	#endif
1119
1120	! - Velocity profiles
1121	IF ( ln_vel3d ) THEN
1122	DO jveloset = 1, nvelosets
1123	! zonal component of velocity
1124	CALL obs_vel_opt( veldatqc(jveloset), kstp, jpi, jpj, jpk, &
1125	& nit000, idaystp, un, vn, gdept_0, umask, vmask, &
1126	n1dint, n2dint, ld_velav(jveloset) )
1127	END DO
1128	ENDIF
1129
1130	#if ! defined key_ice_lim
1131	IF(wrk_not_released(2, 1))THEN
1132	CALL ctl_stop('dia_obs : failed to release workspace array.')
1133	END IF
1134	#endif
1135
1136	END SUBROUTINE dia_obs
1137
1138	SUBROUTINE dia_obs_wri
1139	!!----------------------------------------------------------------------
1140	!! * ROUTINE dia_obs_wri *
1141	!!
1142	!! ** Purpose : Call observation diagnostic output routines
1143	!!
1144	!! ** Method : Call observation diagnostic output routines
1145	!!
1146	!! ** Action :
1147	!!
1148	!! History :
1149	!! ! 06-03 (K. Mogensen) Original code
1150	!! ! 06-05 (K. Mogensen) Reformatted
1151	!! ! 06-10 (A. Weaver) Cleaning
1152	!! ! 07-03 (K. Mogensen) General handling of profiles
1153	!! ! 08-09 (M. Valdivieso) Velocity component (U,V) profiles
1154	!!----------------------------------------------------------------------
1155	IMPLICIT NONE
1156
1157	!! * Local declarations
1158
1159	INTEGER :: jprofset ! Profile data set loop variable
1160	INTEGER :: jveloset ! Velocity data set loop variable
1161	INTEGER :: jslaset ! SLA data set loop variable
1162	INTEGER :: jsstset ! SST data set loop variable
1163	INTEGER :: jseaiceset ! Sea Ice data set loop variable
1164	INTEGER :: jset
1165	INTEGER :: jfbini
1166	CHARACTER(LEN=20) :: datestr=" ",timestr=" "
1167	CHARACTER(LEN=10) :: cdtmp
1168	!-----------------------------------------------------------------------
1169	! Depending on switches call various observation output routines
1170	!-----------------------------------------------------------------------
1171
1172	! - Temperature/salinity profiles
1173
1174	IF( ln_t3d .OR. ln_s3d ) THEN
1175
1176	! Copy data from prodatqc to profdata structures
1177	DO jprofset = 1, nprofsets
1178
1179	CALL obs_prof_decompress( prodatqc(jprofset), &
1180	& profdata(jprofset), .TRUE., numout )
1181
1182	END DO
1183
1184	! Write the profiles.
1185
1186	jprofset = 0
1187
1188	! ENACT insitu data
1189
1190	IF ( ln_ena ) THEN
1191
1192	jprofset = jprofset + 1
1193
1194	CALL obs_wri_p3d( 'enact', profdata(jprofset) )
1195
1196	ENDIF
1197
1198	! Coriolis insitu data
1199
1200	IF ( ln_cor ) THEN
1201
1202	jprofset = jprofset + 1
1203
1204	CALL obs_wri_p3d( 'corio', profdata(jprofset) )
1205
1206	ENDIF
1207
1208	! Feedback insitu data
1209
1210	IF ( ln_profb ) THEN
1211
1212	jfbini = jprofset + 1
1213
1214	DO jprofset = jfbini, nprofsets
1215
1216	jset = jprofset - jfbini + 1
1217	WRITE(cdtmp,'(A,I2.2)')'profb_',jset
1218	CALL obs_wri_p3d( cdtmp, profdata(jprofset) )
1219
1220	END DO
1221
1222	ENDIF
1223
1224	ENDIF
1225
1226	! - Sea surface anomaly
1227	IF ( ln_sla ) THEN
1228
1229	! Copy data from sladatqc to sladata structures
1230	DO jslaset = 1, nslasets
1231
1232	CALL obs_surf_decompress( sladatqc(jslaset), &
1233	& sladata(jslaset), .TRUE., numout )
1234
1235	END DO
1236
1237	jslaset = 0
1238
1239	! Write the AVISO SLA data
1240
1241	IF ( ln_sladt ) THEN
1242
1243	jslaset = 1
1244	CALL obs_wri_sla( 'aviso_act', sladata(jslaset) )
1245	jslaset = 2
1246	CALL obs_wri_sla( 'aviso_pas', sladata(jslaset) )
1247
1248	ENDIF
1249
1250	IF ( ln_slafb ) THEN
1251
1252	jfbini = jslaset + 1
1253
1254	DO jslaset = jfbini, nslasets
1255
1256	jset = jslaset - jfbini + 1
1257	WRITE(cdtmp,'(A,I2.2)')'slafb_',jset
1258	CALL obs_wri_sla( cdtmp, sladata(jslaset) )
1259
1260	END DO
1261
1262	ENDIF
1263
1264	ENDIF
1265
1266	! - Sea surface temperature
1267	IF ( ln_sst ) THEN
1268
1269	! Copy data from sstdatqc to sstdata structures
1270	DO jsstset = 1, nsstsets
1271
1272	CALL obs_surf_decompress( sstdatqc(jsstset), &
1273	& sstdata(jsstset), .TRUE., numout )
1274
1275	END DO
1276
1277	jsstset = 0
1278
1279	! Write the AVISO SST data
1280
1281	IF ( ln_reysst ) THEN
1282
1283	jsstset = jsstset + 1
1284	CALL obs_wri_sst( 'reynolds', sstdata(jsstset) )
1285
1286	ENDIF
1287
1288	IF ( ln_ghrsst ) THEN
1289
1290	jsstset = jsstset + 1
1291	CALL obs_wri_sst( 'ghr', sstdata(jsstset) )
1292
1293	ENDIF
1294
1295	IF ( ln_sstfb ) THEN
1296
1297	jfbini = jsstset + 1
1298
1299	DO jsstset = jfbini, nsstsets
1300
1301	jset = jsstset - jfbini + 1
1302	WRITE(cdtmp,'(A,I2.2)')'sstfb_',jset
1303	CALL obs_wri_sst( cdtmp, sstdata(jsstset) )
1304
1305	END DO
1306
1307	ENDIF
1308
1309	ENDIF
1310
1311	! - Sea surface salinity
1312	IF ( ln_sss ) THEN
1313	IF(lwp) WRITE(numout,*) ' SSS currently not available'
1314	ENDIF
1315
1316	! - Sea Ice Concentration
1317	IF ( ln_seaice ) THEN
1318
1319	! Copy data from seaicedatqc to seaicedata structures
1320	DO jseaiceset = 1, nseaicesets
1321
1322	CALL obs_surf_decompress( seaicedatqc(jseaiceset), &
1323	& seaicedata(jseaiceset), .TRUE., numout )
1324
1325	END DO
1326
1327	! Write the Sea Ice data
1328	DO jseaiceset = 1, nseaicesets
1329
1330	WRITE(cdtmp,'(A,I2.2)')'seaicefb_',jseaiceset
1331	CALL obs_wri_seaice( cdtmp, seaicedata(jseaiceset) )
1332
1333	END DO
1334
1335	ENDIF
1336
1337	! Velocity data
1338	IF( ln_vel3d ) THEN
1339
1340	! Copy data from veldatqc to velodata structures
1341	DO jveloset = 1, nvelosets
1342
1343	CALL obs_prof_decompress( veldatqc(jveloset), &
1344	& velodata(jveloset), .TRUE., numout )
1345
1346	END DO
1347
1348	! Write the profiles.
1349
1350	jveloset = 0
1351
1352	! Daily averaged data
1353
1354	IF ( ln_velavcur ) THEN
1355
1356	jveloset = jveloset + 1
1357
1358	CALL obs_wri_vel( 'velavcurr', velodata(jveloset), n2dint )
1359
1360	ENDIF
1361
1362	! High frequency data
1363
1364	IF ( ln_velhrcur ) THEN
1365
1366	jveloset = jveloset + 1
1367
1368	CALL obs_wri_vel( 'velhrcurr', velodata(jveloset), n2dint )
1369
1370	ENDIF
1371
1372	! Daily averaged data
1373
1374	IF ( ln_velavadcp ) THEN
1375
1376	jveloset = jveloset + 1
1377
1378	CALL obs_wri_vel( 'velavadcp', velodata(jveloset), n2dint )
1379
1380	ENDIF
1381
1382	! High frequency data
1383
1384	IF ( ln_velhradcp ) THEN
1385
1386	jveloset = jveloset + 1
1387
1388	CALL obs_wri_vel( 'velhradcp', velodata(jveloset), n2dint )
1389
1390	ENDIF
1391
1392	! Feedback velocity data
1393
1394	IF ( ln_velfb ) THEN
1395
1396	jfbini = jveloset + 1
1397
1398	DO jveloset = jfbini, nvelosets
1399
1400	jset = jveloset - jfbini + 1
1401	WRITE(cdtmp,'(A,I2.2)')'velfb_',jset
1402	CALL obs_wri_vel( cdtmp, velodata(jveloset), n2dint )
1403
1404	END DO
1405
1406	ENDIF
1407
1408	ENDIF
1409
1410	END SUBROUTINE dia_obs_wri
1411
1412	SUBROUTINE ini_date( ddobsini )
1413	!!----------------------------------------------------------------------
1414	!! * ROUTINE ini_date *
1415	!!
1416	!! ** Purpose : Get initial data in double precision YYYYMMDD.HHMMSS format
1417	!!
1418	!! ** Method : Get initial data in double precision YYYYMMDD.HHMMSS format
1419	!!
1420	!! ** Action : Get initial data in double precision YYYYMMDD.HHMMSS format
1421	!!
1422	!! History :
1423	!! ! 06-03 (K. Mogensen) Original code
1424	!! ! 06-05 (K. Mogensen) Reformatted
1425	!! ! 06-10 (A. Weaver) Cleaning
1426	!! ! 06-10 (G. Smith) Calculates initial date the same as method for final date
1427	!! ! 10-05 (D. Lea) Update to month length calculation for NEMO vn3.2
1428	!!----------------------------------------------------------------------
1429	USE phycst, ONLY : & ! Physical constants
1430	& rday
1431	! USE daymod, ONLY : & ! Time variables
1432	! & nmonth_len
1433	USE dom_oce, ONLY : & ! Ocean space and time domain variables
1434	& rdt
1435
1436	IMPLICIT NONE
1437
1438	!! * Arguments
1439	REAL(KIND=dp), INTENT(OUT) :: ddobsini ! Initial date in YYYYMMDD.HHMMSS
1440
1441	!! * Local declarations
1442	INTEGER :: iyea ! date - (year, month, day, hour, minute)
1443	INTEGER :: imon
1444	INTEGER :: iday
1445	INTEGER :: ihou
1446	INTEGER :: imin
1447	INTEGER :: imday ! Number of days in month.
1448	REAL(KIND=wp) :: zdayfrc ! Fraction of day
1449
1450	INTEGER, DIMENSION(12) :: imonth_len !: length in days of the months of the current year
1451
1452	!!----------------------------------------------------------------------
1453	!! Initial date initialization (year, month, day, hour, minute)
1454	!! (This assumes that the initial date is for 00z))
1455	!!----------------------------------------------------------------------
1456	iyea = ndate0 / 10000
1457	imon = ( ndate0 - iyea * 10000 ) / 100
1458	iday = ndate0 - iyea * 10000 - imon * 100
1459	ihou = 0
1460	imin = 0
1461
1462	!!----------------------------------------------------------------------
1463	!! Compute number of days + number of hours + min since initial time
1464	!!----------------------------------------------------------------------
1465	iday = iday + ( nit000 -1 ) * rdt / rday
1466	zdayfrc = ( nit000 -1 ) * rdt / rday
1467	zdayfrc = zdayfrc - aint(zdayfrc)
1468	ihou = int( zdayfrc * 24 )
1469	imin = int( (zdayfrc * 24 - ihou) * 60 )
1470
1471	!!-----------------------------------------------------------------------
1472	!! Convert number of days (iday) into a real date
1473	!!----------------------------------------------------------------------
1474
1475	CALL calc_month_len( iyea, imonth_len )
1476
1477	DO WHILE ( iday > imonth_len(imon) )
1478	iday = iday - imonth_len(imon)
1479	imon = imon + 1
1480	IF ( imon > 12 ) THEN
1481	imon = 1
1482	iyea = iyea + 1
1483	CALL calc_month_len( iyea, imonth_len ) ! update month lengths
1484	ENDIF
1485	END DO
1486
1487	!!----------------------------------------------------------------------
1488	!! Convert it into YYYYMMDD.HHMMSS format.
1489	!!----------------------------------------------------------------------
1490	ddobsini = iyea * 10000_dp + imon * 100_dp + &
1491	& iday + ihou * 0.01_dp + imin * 0.0001_dp
1492
1493
1494	END SUBROUTINE ini_date
1495
1496	SUBROUTINE fin_date( ddobsfin )
1497	!!----------------------------------------------------------------------
1498	!! * ROUTINE fin_date *
1499	!!
1500	!! ** Purpose : Get final data in double precision YYYYMMDD.HHMMSS format
1501	!!
1502	!! ** Method : Get final data in double precision YYYYMMDD.HHMMSS format
1503	!!
1504	!! ** Action : Get final data in double precision YYYYMMDD.HHMMSS format
1505	!!
1506	!! History :
1507	!! ! 06-03 (K. Mogensen) Original code
1508	!! ! 06-05 (K. Mogensen) Reformatted
1509	!! ! 06-10 (A. Weaver) Cleaning
1510	!! ! 10-05 (D. Lea) Update to month length calculation for NEMO vn3.2
1511	!!----------------------------------------------------------------------
1512	USE phycst, ONLY : & ! Physical constants
1513	& rday
1514	! USE daymod, ONLY : & ! Time variables
1515	! & nmonth_len
1516	USE dom_oce, ONLY : & ! Ocean space and time domain variables
1517	& rdt
1518
1519	IMPLICIT NONE
1520
1521	!! * Arguments
1522	REAL(KIND=dp), INTENT(OUT) :: ddobsfin ! Final date in YYYYMMDD.HHMMSS
1523
1524	!! * Local declarations
1525	INTEGER :: iyea ! date - (year, month, day, hour, minute)
1526	INTEGER :: imon
1527	INTEGER :: iday
1528	INTEGER :: ihou
1529	INTEGER :: imin
1530	INTEGER :: imday ! Number of days in month.
1531	REAL(KIND=wp) :: zdayfrc ! Fraction of day
1532
1533	INTEGER, DIMENSION(12) :: imonth_len !: length in days of the months of the current year
1534
1535	!-----------------------------------------------------------------------
1536	! Initial date initialization (year, month, day, hour, minute)
1537	! (This assumes that the initial date is for 00z)
1538	!-----------------------------------------------------------------------
1539	iyea = ndate0 / 10000
1540	imon = ( ndate0 - iyea * 10000 ) / 100
1541	iday = ndate0 - iyea * 10000 - imon * 100
1542	ihou = 0
1543	imin = 0
1544
1545	!-----------------------------------------------------------------------
1546	! Compute number of days + number of hours + min since initial time
1547	!-----------------------------------------------------------------------
1548	iday = iday + nitend * rdt / rday
1549	zdayfrc = nitend * rdt / rday
1550	zdayfrc = zdayfrc - AINT( zdayfrc )
1551	ihou = INT( zdayfrc * 24 )
1552	imin = INT( ( zdayfrc * 24 - ihou ) * 60 )
1553
1554	!-----------------------------------------------------------------------
1555	! Convert number of days (iday) into a real date
1556	!----------------------------------------------------------------------
1557
1558	CALL calc_month_len( iyea, imonth_len )
1559
1560	DO WHILE ( iday > imonth_len(imon) )
1561	iday = iday - imonth_len(imon)
1562	imon = imon + 1
1563	IF ( imon > 12 ) THEN
1564	imon = 1
1565	iyea = iyea + 1
1566	CALL calc_month_len( iyea, imonth_len ) ! update month lengths
1567	ENDIF
1568	END DO
1569
1570	!-----------------------------------------------------------------------
1571	! Convert it into YYYYMMDD.HHMMSS format
1572	!-----------------------------------------------------------------------
1573	ddobsfin = iyea * 10000_dp + imon * 100_dp + iday &
1574	& + ihou * 0.01_dp + imin * 0.0001_dp
1575
1576	END SUBROUTINE fin_date
1577
1578	END MODULE diaobs

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

Download in other formats: