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/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

Context Navigation

source: branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 4768

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

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

Download in other formats: