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

Context Navigation

source: trunk/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 4624

Last change on this file since 4624 was 4624, checked in by acc, 10 years ago
#1305. Fix slow start-up problems on some systems by introducing and using lwm logical to restrict output of merged namelists to the first (or only) processor. lwm is true only on the first processor regardless of ln_ctl. Small changes to all flavours of nemogcm.F90 are also required to write namctl and namcfg after the call to mynode which now opens output.namelist.dyn and writes nammpp.
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: