Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

diaobs.F90 @ 11274

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

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

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

Context Navigation

source: branches/UKMO/dev_r5518_GO6_package_STARTHOUR/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 11274

Download in other formats: