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

Context Navigation

source: branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 4147

Last change on this file since 4147 was 4147, checked in by cetlod, 10 years ago
merge in dev_LOCEAN_2013, the 1st development branch dev_r3853_CNRS9_Confsetting, from its starting point ( r3853 ) on the trunk: see ticket #1169
Property svn:keywords set to `Id`
File size: 57.9 KB

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

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

Download in other formats: