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

Context Navigation

source: trunk/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 3294

Last change on this file since 3294 was 3294, checked in by rblod, 12 years ago
Merge of 3.4beta into the trunk
Property svn:keywords set to `Id`
File size: 58.3 KB

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

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

Download in other formats: