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/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

Context Navigation

source: branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 5106

Last change on this file since 5106 was 5106, checked in by jwhile, 9 years ago
Minor fixes due to review
Property svn:keywords set to `Id`
File size: 61.4 KB

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

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

Download in other formats: