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

Context Navigation

source: branches/dev_1784_OBS/NEMO/OPA_SRC/OBS/diaobs.F90 @ 2001

Last change on this file since 2001 was 2001, checked in by djlea, 14 years ago
Adding observation operator code
File size: 61.0 KB

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

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

Download in other formats: