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

Context Navigation

source: branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 6406

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

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

Download in other formats: