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_UKMO13_CICE_changes_take2/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

Context Navigation

source: branches/2014/dev_r4650_UKMO13_CICE_changes_take2/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 4921

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

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

Download in other formats: