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

Context Navigation

source: branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 4897

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

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

Download in other formats: