Context Navigation

diaobs.F90 @ 11101

Last change on this file since 11101 was 11101, checked in by frrh, 5 years ago

Merge changes from Met Office GMED ticket 450 to reduce unnecessary
text output from NEMO.
This output, which is typically not switchable, is rarely of interest
in normal (non-debugging) runs and simply redunantley consumes extra
file space.
Further, the presence of this text output has been shown to
significantly degrade performance of models which are run during
Met Office HPC RAID (disk) checks.
The new code introduces switches which are configurable via the
changes made in the associated Met Office MOCI ticket 399.

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

Download in other formats: