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

Context Navigation

source: branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 5992

Last change on this file since 5992 was 5992, checked in by timgraham, 8 years ago
Merged branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION into branch
Property svn:keywords set to `Id`
File size: 62.6 KB

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

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

Download in other formats: