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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

diaobs.F90 in branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

Context Navigation

source: branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 6043

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

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

Download in other formats: