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.11_SETTE_OBSASM/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

Context Navigation

source: branches/2014/dev_r4650_UKMO14.11_SETTE_OBSASM/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 4781

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

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

Download in other formats: