Context Navigation

diaobs.F90 @ 7915

Last change on this file since 7915 was 7862, checked in by mattmartin, 7 years ago
Change standard output writing.
File size: 59.6 KB

Line
1	MODULE diaobs
2	!!======================================================================
3	!! * MODULE diaobs *
4	!! Observation diagnostics: Computation of the misfit between data and
5	!! their model equivalent
6	!!======================================================================
7
8	!!----------------------------------------------------------------------
9	!! dia_obs_init : Reading and prepare observations
10	!! dia_obs : Compute model equivalent to observations
11	!! dia_obs_wri : Write observational diagnostics
12	!! ini_date : Compute the initial date YYYYMMDD.HHMMSS
13	!! fin_date : Compute the final date YYYYMMDD.HHMMSS
14	!!----------------------------------------------------------------------
15	!! * Modules used
16	USE wrk_nemo ! Memory Allocation
17	USE par_kind ! Precision variables
18	USE in_out_manager ! I/O manager
19	USE par_oce
20	USE dom_oce ! Ocean space and time domain variables
21	USE obs_read_prof ! Reading and allocation of profile obs
22	USE obs_read_surf ! Reading and allocation of surface obs
23	USE obs_readmdt ! Reading and allocation of MDT for SLA.
24	USE obs_prep ! Preparation of obs. (grid search etc).
25	USE obs_oper ! Observation operators
26	USE obs_write ! Writing of observation related diagnostics
27	USE obs_grid ! Grid searching
28	USE obs_read_altbias ! Bias treatment for altimeter
29	USE obs_sstbias ! Bias correction routine for SST
30	USE obs_profiles_def ! Profile data definitions
31	USE obs_surf_def ! Surface data definitions
32	USE obs_types ! Definitions for observation types
33	USE mpp_map ! MPP mapping
34	USE lib_mpp ! For ctl_warn/stop
35
36	IMPLICIT NONE
37
38	!! * Routine accessibility
39	PRIVATE
40	PUBLIC dia_obs_init, & ! Initialize and read observations
41	& dia_obs, & ! Compute model equivalent to observations
42	& dia_obs_wri, & ! Write model equivalent to observations
43	& dia_obs_dealloc ! Deallocate dia_obs data
44
45	!! * Module variables
46	LOGICAL, PUBLIC :: &
47	& lk_diaobs = .TRUE. !: Include this for backwards compatibility at NEMO 3.6.
48	LOGICAL :: ln_diaobs !: Logical switch for the obs operator
49	LOGICAL :: ln_sstnight !: Logical switch for night mean SST obs
50	LOGICAL :: ln_sla_fp_indegs !: T=> SLA obs footprint size specified in degrees, F=> in metres
51	LOGICAL :: ln_sst_fp_indegs !: T=> SST obs footprint size specified in degrees, F=> in metres
52	LOGICAL :: ln_sss_fp_indegs !: T=> SSS obs footprint size specified in degrees, F=> in metres
53	LOGICAL :: ln_sic_fp_indegs !: T=> sea-ice obs footprint size specified in degrees, F=> in metres
54
55	REAL(wp) :: rn_sla_avglamscl !: E/W diameter of SLA observation footprint (metres)
56	REAL(wp) :: rn_sla_avgphiscl !: N/S diameter of SLA observation footprint (metres)
57	REAL(wp) :: rn_sst_avglamscl !: E/W diameter of SST observation footprint (metres)
58	REAL(wp) :: rn_sst_avgphiscl !: N/S diameter of SST observation footprint (metres)
59	REAL(wp) :: rn_sss_avglamscl !: E/W diameter of SSS observation footprint (metres)
60	REAL(wp) :: rn_sss_avgphiscl !: N/S diameter of SSS observation footprint (metres)
61	REAL(wp) :: rn_sic_avglamscl !: E/W diameter of sea-ice observation footprint (metres)
62	REAL(wp) :: rn_sic_avgphiscl !: N/S diameter of sea-ice observation footprint (metres)
63
64	INTEGER :: nn_1dint !: Vertical interpolation method
65	INTEGER :: nn_2dint !: Default horizontal interpolation method
66	INTEGER :: nn_2dint_sla !: SLA horizontal interpolation method
67	INTEGER :: nn_2dint_sst !: SST horizontal interpolation method
68	INTEGER :: nn_2dint_sss !: SSS horizontal interpolation method
69	INTEGER :: nn_2dint_sic !: Seaice horizontal interpolation method
70
71	INTEGER, DIMENSION(imaxavtypes) :: &
72	& nn_profdavtypes !: Profile data types representing a daily average
73	INTEGER :: nproftypes !: Number of profile obs types
74	INTEGER :: nsurftypes !: Number of surface obs types
75	INTEGER, DIMENSION(:), ALLOCATABLE :: &
76	& nvarsprof, & !: Number of profile variables
77	& nvarssurf !: Number of surface variables
78	INTEGER, DIMENSION(:), ALLOCATABLE :: &
79	& nextrprof, & !: Number of profile extra variables
80	& nextrsurf !: Number of surface extra variables
81	INTEGER, DIMENSION(:), ALLOCATABLE :: &
82	& n2dintsurf !: Interpolation option for surface variables
83	REAL(wp), DIMENSION(:), ALLOCATABLE :: &
84	& ravglamscl, & !: E/W diameter of averaging footprint for surface variables
85	& ravgphiscl !: N/S diameter of averaging footprint for surface variables
86	LOGICAL, DIMENSION(:), ALLOCATABLE :: &
87	& lfpindegs, & !: T=> surface obs footprint size specified in degrees, F=> in metres
88	& llnightav !: Logical for calculating night-time averages
89
90	TYPE(obs_surf), PUBLIC, POINTER, DIMENSION(:) :: &
91	& surfdata, & !: Initial surface data
92	& surfdataqc !: Surface data after quality control
93	TYPE(obs_prof), PUBLIC, POINTER, DIMENSION(:) :: &
94	& profdata, & !: Initial profile data
95	& profdataqc !: Profile data after quality control
96
97	CHARACTER(len=6), PUBLIC, DIMENSION(:), ALLOCATABLE :: &
98	& cobstypesprof, & !: Profile obs types
99	& cobstypessurf !: Surface obs types
100
101	!!----------------------------------------------------------------------
102	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
103	!! $Id$
104	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
105	!!----------------------------------------------------------------------
106
107	!! * Substitutions
108	# include "domzgr_substitute.h90"
109	CONTAINS
110
111	SUBROUTINE dia_obs_init
112	!!----------------------------------------------------------------------
113	!! * ROUTINE dia_obs_init *
114	!!
115	!! ** Purpose : Initialize and read observations
116	!!
117	!! ** Method : Read the namelist and call reading routines
118	!!
119	!! ** Action : Read the namelist and call reading routines
120	!!
121	!! History :
122	!! ! 06-03 (K. Mogensen) Original code
123	!! ! 06-05 (A. Weaver) Reformatted
124	!! ! 06-10 (A. Weaver) Cleaning and add controls
125	!! ! 07-03 (K. Mogensen) General handling of profiles
126	!! ! 14-08 (J.While) Incorporated SST bias correction
127	!! ! 15-02 (M. Martin) Simplification of namelist and code
128	!!----------------------------------------------------------------------
129
130	IMPLICIT NONE
131
132	!! * Local declarations
133	INTEGER, PARAMETER :: &
134	& jpmaxnfiles = 1000 ! Maximum number of files for each obs type
135	INTEGER, DIMENSION(:), ALLOCATABLE :: &
136	& ifilesprof, & ! Number of profile files
137	& ifilessurf ! Number of surface files
138	INTEGER :: ios ! Local integer output status for namelist read
139	INTEGER :: jtype ! Counter for obs types
140	INTEGER :: jvar ! Counter for variables
141	INTEGER :: jfile ! Counter for files
142	INTEGER :: jnumsstbias ! Number of SST bias files to read and apply
143
144	CHARACTER(len=128), DIMENSION(jpmaxnfiles) :: &
145	& cn_profbfiles, & ! T/S profile input filenames
146	& cn_sstfbfiles, & ! Sea surface temperature input filenames
147	& cn_slafbfiles, & ! Sea level anomaly input filenames
148	& cn_sicfbfiles, & ! Seaice concentration input filenames
149	& cn_velfbfiles, & ! Velocity profile input filenames
150	& cn_sssfbfiles, & ! Sea surface salinity input filenames
151	& cn_logchlfbfiles, & ! Log(Chl) input filenames
152	& cn_spmfbfiles, & ! Sediment input filenames
153	& cn_fco2fbfiles, & ! fco2 input filenames
154	& cn_pco2fbfiles, & ! pco2 input filenames
155	& cn_sstbiasfiles ! SST bias input filenames
156
157	CHARACTER(LEN=128) :: &
158	& cn_altbiasfile ! Altimeter bias input filename
159
160
161	LOGICAL :: ln_t3d ! Logical switch for temperature profiles
162	LOGICAL :: ln_s3d ! Logical switch for salinity profiles
163	LOGICAL :: ln_sla ! Logical switch for sea level anomalies
164	LOGICAL :: ln_sst ! Logical switch for sea surface temperature
165	LOGICAL :: ln_sic ! Logical switch for sea ice concentration
166	LOGICAL :: ln_sss ! Logical switch for sea surface salinity obs
167	LOGICAL :: ln_vel3d ! Logical switch for velocity (u,v) obs
168	LOGICAL :: ln_logchl ! Logical switch for log(Chl) obs
169	LOGICAL :: ln_spm ! Logical switch for sediment obs
170	LOGICAL :: ln_fco2 ! Logical switch for fco2 obs
171	LOGICAL :: ln_pco2 ! Logical switch for pco2 obs
172	LOGICAL :: ln_nea ! Logical switch to remove obs near land
173	LOGICAL :: ln_altbias ! Logical switch for altimeter bias
174	LOGICAL :: ln_sstbias ! Logical switch for bias correction of SST
175	LOGICAL :: ln_ignmis ! Logical switch for ignoring missing files
176	LOGICAL :: ln_s_at_t ! Logical switch to compute model S at T obs
177	LOGICAL :: ln_bound_reject ! Logical switch for rejecting obs near the boundary
178
179	REAL(dp) :: rn_dobsini ! Obs window start date YYYYMMDD.HHMMSS
180	REAL(dp) :: rn_dobsend ! Obs window end date YYYYMMDD.HHMMSS
181
182	CHARACTER(len=128), DIMENSION(:,:), ALLOCATABLE :: &
183	& clproffiles, & ! Profile filenames
184	& clsurffiles ! Surface filenames
185
186	LOGICAL :: llvar1 ! Logical for profile variable 1
187	LOGICAL :: llvar2 ! Logical for profile variable 1
188
189	REAL(wp), POINTER, DIMENSION(:,:) :: &
190	& zglam1, & ! Model longitudes for profile variable 1
191	& zglam2 ! Model longitudes for profile variable 2
192	REAL(wp), POINTER, DIMENSION(:,:) :: &
193	& zgphi1, & ! Model latitudes for profile variable 1
194	& zgphi2 ! Model latitudes for profile variable 2
195	REAL(wp), POINTER, DIMENSION(:,:,:) :: &
196	& zmask1, & ! Model land/sea mask associated with variable 1
197	& zmask2 ! Model land/sea mask associated with variable 2
198
199
200	NAMELIST/namobs/ln_diaobs, ln_t3d, ln_s3d, ln_sla, &
201	& ln_sst, ln_sic, ln_sss, ln_vel3d, &
202	& ln_logchl, ln_spm, ln_fco2, ln_pco2, &
203	& ln_altbias, ln_sstbias, ln_nea, &
204	& ln_grid_global, ln_grid_search_lookup, &
205	& ln_ignmis, ln_s_at_t, ln_bound_reject, &
206	& ln_sstnight, &
207	& ln_sla_fp_indegs, ln_sst_fp_indegs, &
208	& ln_sss_fp_indegs, ln_sic_fp_indegs, &
209	& cn_profbfiles, cn_slafbfiles, &
210	& cn_sstfbfiles, cn_sicfbfiles, &
211	& cn_velfbfiles, cn_sssfbfiles, &
212	& cn_logchlfbfiles, cn_spmfbfiles, &
213	& cn_fco2fbfiles, cn_pco2fbfiles, &
214	& cn_sstbiasfiles, cn_altbiasfile, &
215	& cn_gridsearchfile, rn_gridsearchres, &
216	& rn_dobsini, rn_dobsend, &
217	& rn_sla_avglamscl, rn_sla_avgphiscl, &
218	& rn_sst_avglamscl, rn_sst_avgphiscl, &
219	& rn_sss_avglamscl, rn_sss_avgphiscl, &
220	& rn_sic_avglamscl, rn_sic_avgphiscl, &
221	& nn_1dint, nn_2dint, &
222	& nn_2dint_sla, nn_2dint_sst, &
223	& nn_2dint_sss, nn_2dint_sic, &
224	& nn_msshc, rn_mdtcorr, rn_mdtcutoff, &
225	& nn_profdavtypes
226
227	CALL wrk_alloc( jpi, jpj, zglam1 )
228	CALL wrk_alloc( jpi, jpj, zglam2 )
229	CALL wrk_alloc( jpi, jpj, zgphi1 )
230	CALL wrk_alloc( jpi, jpj, zgphi2 )
231	CALL wrk_alloc( jpi, jpj, jpk, zmask1 )
232	CALL wrk_alloc( jpi, jpj, jpk, zmask2 )
233
234	!-----------------------------------------------------------------------
235	! Read namelist parameters
236	!-----------------------------------------------------------------------
237
238	! Some namelist arrays need initialising
239	cn_profbfiles(:) = ''
240	cn_slafbfiles(:) = ''
241	cn_sstfbfiles(:) = ''
242	cn_sicfbfiles(:) = ''
243	cn_velfbfiles(:) = ''
244	cn_sssfbfiles(:) = ''
245	cn_logchlfbfiles(:) = ''
246	cn_spmfbfiles(:) = ''
247	cn_fco2fbfiles(:) = ''
248	cn_pco2fbfiles(:) = ''
249	cn_sstbiasfiles(:) = ''
250	nn_profdavtypes(:) = -1
251
252	CALL ini_date( rn_dobsini )
253	CALL fin_date( rn_dobsend )
254
255	! Read namelist namobs : control observation diagnostics
256	REWIND( numnam_ref ) ! Namelist namobs in reference namelist
257	READ ( numnam_ref, namobs, IOSTAT = ios, ERR = 901)
258	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in reference namelist', lwp )
259
260	REWIND( numnam_cfg ) ! Namelist namobs in configuration namelist
261	READ ( numnam_cfg, namobs, IOSTAT = ios, ERR = 902 )
262	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in configuration namelist', lwp )
263	IF(lwm) WRITE ( numond, namobs )
264
265	lk_diaobs = .FALSE.
266	#if defined key_diaobs
267	IF ( ln_diaobs ) lk_diaobs = .TRUE.
268	#endif
269
270	IF ( .NOT. lk_diaobs ) THEN
271	IF(lwp) WRITE(numout,cform_war)
272	IF(lwp) WRITE(numout,*)' ln_diaobs is set to false or key_diaobs is not set, so not calling dia_obs'
273	RETURN
274	ENDIF
275
276	IF(lwp) THEN
277	WRITE(numout,*)
278	WRITE(numout,*) 'dia_obs_init : Observation diagnostic initialization'
279	WRITE(numout,*) '~~~~~~~~~~~~'
280	WRITE(numout,*) ' Namelist namobs : set observation diagnostic parameters'
281	WRITE(numout,*) ' Logical switch for T profile observations ln_t3d = ', ln_t3d
282	WRITE(numout,*) ' Logical switch for S profile observations ln_s3d = ', ln_s3d
283	WRITE(numout,*) ' Logical switch for SLA observations ln_sla = ', ln_sla
284	WRITE(numout,*) ' Logical switch for SST observations ln_sst = ', ln_sst
285	WRITE(numout,*) ' Logical switch for Sea Ice observations ln_sic = ', ln_sic
286	WRITE(numout,*) ' Logical switch for velocity observations ln_vel3d = ', ln_vel3d
287	WRITE(numout,*) ' Logical switch for SSS observations ln_sss = ', ln_sss
288	WRITE(numout,*) ' Logical switch for log(Chl) observations ln_logchl = ', ln_logchl
289	WRITE(numout,*) ' Logical switch for SPM observations ln_spm = ', ln_spm
290	WRITE(numout,*) ' Logical switch for FCO2 observations ln_fco2 = ', ln_fco2
291	WRITE(numout,*) ' Logical switch for PCO2 observations ln_pco2 = ', ln_pco2
292	WRITE(numout,*) ' Global distribution of observations ln_grid_global = ', ln_grid_global
293	WRITE(numout,*) ' Logical switch for obs grid search lookup ln_grid_search_lookup = ', ln_grid_search_lookup
294	IF (ln_grid_search_lookup) &
295	WRITE(numout,*) ' Grid search lookup file header cn_gridsearchfile = ', cn_gridsearchfile
296	WRITE(numout,*) ' Initial date in window YYYYMMDD.HHMMSS rn_dobsini = ', rn_dobsini
297	WRITE(numout,*) ' Final date in window YYYYMMDD.HHMMSS rn_dobsend = ', rn_dobsend
298	WRITE(numout,*) ' Type of vertical interpolation method nn_1dint = ', nn_1dint
299	WRITE(numout,*) ' Type of horizontal interpolation method nn_2dint = ', nn_2dint
300	WRITE(numout,*) ' Rejection of observations near land switch ln_nea = ', ln_nea
301	WRITE(numout,*) ' Rejection of obs near open bdys ln_bound_reject = ', ln_bound_reject
302	WRITE(numout,*) ' MSSH correction scheme nn_msshc = ', nn_msshc
303	WRITE(numout,*) ' MDT correction rn_mdtcorr = ', rn_mdtcorr
304	WRITE(numout,*) ' MDT cutoff for computed correction rn_mdtcutoff = ', rn_mdtcutoff
305	WRITE(numout,*) ' Logical switch for alt bias ln_altbias = ', ln_altbias
306	WRITE(numout,*) ' Logical switch for sst bias ln_sstbias = ', ln_sstbias
307	WRITE(numout,*) ' Logical switch for ignoring missing files ln_ignmis = ', ln_ignmis
308	WRITE(numout,*) ' Daily average types nn_profdavtypes = ', nn_profdavtypes
309	WRITE(numout,*) ' Logical switch for night-time SST obs ln_sstnight = ', ln_sstnight
310	ENDIF
311	!-----------------------------------------------------------------------
312	! Set up list of observation types to be used
313	! and the files associated with each type
314	!-----------------------------------------------------------------------
315
316	nproftypes = COUNT( (/ln_t3d .OR. ln_s3d, ln_vel3d /) )
317	nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic, ln_sss, &
318	& ln_logchl, ln_spm, ln_fco2, ln_pco2 /) )
319
320	IF ( nproftypes == 0 .AND. nsurftypes == 0 ) THEN
321	IF(lwp) WRITE(numout,cform_war)
322	IF(lwp) WRITE(numout,*) ' ln_diaobs is set to true, but all obs operator logical flags', &
323	& ' are set to .FALSE. so turning off calls to dia_obs'
324	nwarn = nwarn + 1
325	lk_diaobs = .FALSE.
326	RETURN
327	ENDIF
328
329	IF(lwp) WRITE(numout,*) ' Number of profile obs types: ',nproftypes
330	IF ( nproftypes > 0 ) THEN
331
332	ALLOCATE( cobstypesprof(nproftypes) )
333	ALLOCATE( ifilesprof(nproftypes) )
334	ALLOCATE( clproffiles(nproftypes,jpmaxnfiles) )
335
336	jtype = 0
337	IF (ln_t3d .OR. ln_s3d) THEN
338	jtype = jtype + 1
339	CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'prof ', &
340	& cn_profbfiles, ifilesprof, cobstypesprof, clproffiles )
341	ENDIF
342	IF (ln_vel3d) THEN
343	jtype = jtype + 1
344	CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'vel ', &
345	& cn_velfbfiles, ifilesprof, cobstypesprof, clproffiles )
346	ENDIF
347
348	ENDIF
349
350	IF(lwp) WRITE(numout,*)' Number of surface obs types: ',nsurftypes
351	IF ( nsurftypes > 0 ) THEN
352
353	ALLOCATE( cobstypessurf(nsurftypes) )
354	ALLOCATE( ifilessurf(nsurftypes) )
355	ALLOCATE( clsurffiles(nsurftypes, jpmaxnfiles) )
356	ALLOCATE(n2dintsurf(nsurftypes))
357	ALLOCATE(ravglamscl(nsurftypes))
358	ALLOCATE(ravgphiscl(nsurftypes))
359	ALLOCATE(lfpindegs(nsurftypes))
360	ALLOCATE(llnightav(nsurftypes))
361
362	jtype = 0
363	IF (ln_sla) THEN
364	jtype = jtype + 1
365	CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sla ', &
366	& cn_slafbfiles, ifilessurf, cobstypessurf, clsurffiles )
367	CALL obs_setinterpopts( nsurftypes, jtype, 'sla ', &
368	& nn_2dint, nn_2dint_sla, &
369	& rn_sla_avglamscl, rn_sla_avgphiscl, &
370	& ln_sla_fp_indegs, .FALSE., &
371	& n2dintsurf, ravglamscl, ravgphiscl, &
372	& lfpindegs, llnightav )
373	ENDIF
374	IF (ln_sst) THEN
375	jtype = jtype + 1
376	CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sst ', &
377	& cn_sstfbfiles, ifilessurf, cobstypessurf, clsurffiles )
378	CALL obs_setinterpopts( nsurftypes, jtype, 'sst ', &
379	& nn_2dint, nn_2dint_sst, &
380	& rn_sst_avglamscl, rn_sst_avgphiscl, &
381	& ln_sst_fp_indegs, ln_sstnight, &
382	& n2dintsurf, ravglamscl, ravgphiscl, &
383	& lfpindegs, llnightav )
384	ENDIF
385	#if defined key_lim2 \|\| defined key_lim3 \|\| defined key_cice
386	IF (ln_sic) THEN
387	jtype = jtype + 1
388	CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sic ', &
389	& cn_sicfbfiles, ifilessurf, cobstypessurf, clsurffiles )
390	CALL obs_setinterpopts( nsurftypes, jtype, 'sic ', &
391	& nn_2dint, nn_2dint_sic, &
392	& rn_sic_avglamscl, rn_sic_avgphiscl, &
393	& ln_sic_fp_indegs, .FALSE., &
394	& n2dintsurf, ravglamscl, ravgphiscl, &
395	& lfpindegs, llnightav )
396	ENDIF
397	#endif
398	IF (ln_sss) THEN
399	jtype = jtype + 1
400	CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sss ', &
401	& cn_sssfbfiles, ifilessurf, cobstypessurf, clsurffiles )
402	CALL obs_setinterpopts( nsurftypes, jtype, 'sss ', &
403	& nn_2dint, nn_2dint_sss, &
404	& rn_sss_avglamscl, rn_sss_avgphiscl, &
405	& ln_sss_fp_indegs, .FALSE., &
406	& n2dintsurf, ravglamscl, ravgphiscl, &
407	& lfpindegs, llnightav )
408	ENDIF
409
410	IF (ln_logchl) THEN
411	jtype = jtype + 1
412	CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'logchl', &
413	& cn_logchlfbfiles, ifilessurf, cobstypessurf, clsurffiles )
414	CALL obs_setinterpopts( nsurftypes, jtype, 'logchl', &
415	& nn_2dint, -1, 0., 0., .TRUE., .FALSE., &
416	& n2dintsurf, ravglamscl, ravgphiscl, &
417	& lfpindegs, llnightav )
418	ENDIF
419
420	IF (ln_spm) THEN
421	jtype = jtype + 1
422	CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'spm ', &
423	& cn_spmfbfiles, ifilessurf, cobstypessurf, clsurffiles )
424	CALL obs_setinterpopts( nsurftypes, jtype, 'spm ', &
425	& nn_2dint, -1, 0., 0., .TRUE., .FALSE., &
426	& n2dintsurf, ravglamscl, ravgphiscl, &
427	& lfpindegs, llnightav )
428	ENDIF
429
430	IF (ln_fco2) THEN
431	jtype = jtype + 1
432	CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'fco2 ', &
433	& cn_fco2fbfiles, ifilessurf, cobstypessurf, clsurffiles )
434	CALL obs_setinterpopts( nsurftypes, jtype, 'fco2 ', &
435	& nn_2dint, -1, 0., 0., .TRUE., .FALSE., &
436	& n2dintsurf, ravglamscl, ravgphiscl, &
437	& lfpindegs, llnightav )
438	ENDIF
439
440	IF (ln_pco2) THEN
441	jtype = jtype + 1
442	CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'pco2 ', &
443	& cn_pco2fbfiles, ifilessurf, cobstypessurf, clsurffiles )
444	CALL obs_setinterpopts( nsurftypes, jtype, 'pco2 ', &
445	& nn_2dint, -1, 0., 0., .TRUE., .FALSE., &
446	& n2dintsurf, ravglamscl, ravgphiscl, &
447	& lfpindegs, llnightav )
448	ENDIF
449
450	ENDIF
451
452	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
453
454
455	!-----------------------------------------------------------------------
456	! Obs operator parameter checking and initialisations
457	!-----------------------------------------------------------------------
458
459	IF ( ln_vel3d .AND. ( .NOT. ln_grid_global ) ) THEN
460	CALL ctl_stop( 'Velocity data only works with ln_grid_global=.true.' )
461	RETURN
462	ENDIF
463
464	IF ( ( nn_1dint < 0 ) .OR. ( nn_1dint > 1 ) ) THEN
465	CALL ctl_stop(' Choice of vertical (1D) interpolation method', &
466	& ' is not available')
467	ENDIF
468
469	IF ( ( nn_2dint < 0 ) .OR. ( nn_2dint > 6 ) ) THEN
470	CALL ctl_stop(' Choice of horizontal (2D) interpolation method', &
471	& ' is not available')
472	ENDIF
473
474	CALL obs_typ_init
475
476	CALL mppmap_init
477
478	CALL obs_grid_setup( )
479
480	!-----------------------------------------------------------------------
481	! Depending on switches read the various observation types
482	!-----------------------------------------------------------------------
483
484	IF ( nproftypes > 0 ) THEN
485
486	ALLOCATE(profdata(nproftypes))
487	ALLOCATE(profdataqc(nproftypes))
488	ALLOCATE(nvarsprof(nproftypes))
489	ALLOCATE(nextrprof(nproftypes))
490
491	DO jtype = 1, nproftypes
492
493	nvarsprof(jtype) = 2
494	IF ( TRIM(cobstypesprof(jtype)) == 'prof' ) THEN
495	nextrprof(jtype) = 1
496	llvar1 = ln_t3d
497	llvar2 = ln_s3d
498	zglam1 = glamt
499	zgphi1 = gphit
500	zmask1 = tmask
501	zglam2 = glamt
502	zgphi2 = gphit
503	zmask2 = tmask
504	ENDIF
505	IF ( TRIM(cobstypesprof(jtype)) == 'vel' ) THEN
506	nextrprof(jtype) = 2
507	llvar1 = ln_vel3d
508	llvar2 = ln_vel3d
509	zglam1 = glamu
510	zgphi1 = gphiu
511	zmask1 = umask
512	zglam2 = glamv
513	zgphi2 = gphiv
514	zmask2 = vmask
515	ENDIF
516
517	!Read in profile or profile obs types
518	CALL obs_rea_prof( profdata(jtype), ifilesprof(jtype), &
519	& clproffiles(jtype,1:ifilesprof(jtype)), &
520	& nvarsprof(jtype), nextrprof(jtype), nitend-nit000+2, &
521	& rn_dobsini, rn_dobsend, llvar1, llvar2, &
522	& ln_ignmis, ln_s_at_t, .FALSE., &
523	& kdailyavtypes = nn_profdavtypes )
524
525	DO jvar = 1, nvarsprof(jtype)
526	CALL obs_prof_staend( profdata(jtype), jvar )
527	END DO
528
529	CALL obs_pre_prof( profdata(jtype), profdataqc(jtype), &
530	& llvar1, llvar2, &
531	& jpi, jpj, jpk, &
532	& zmask1, zglam1, zgphi1, zmask2, zglam2, zgphi2, &
533	& ln_nea, ln_bound_reject, &
534	& kdailyavtypes = nn_profdavtypes )
535
536	END DO
537
538	DEALLOCATE( ifilesprof, clproffiles )
539
540	ENDIF
541
542	IF ( nsurftypes > 0 ) THEN
543
544	ALLOCATE(surfdata(nsurftypes))
545	ALLOCATE(surfdataqc(nsurftypes))
546	ALLOCATE(nvarssurf(nsurftypes))
547	ALLOCATE(nextrsurf(nsurftypes))
548
549	DO jtype = 1, nsurftypes
550
551	nvarssurf(jtype) = 1
552	nextrsurf(jtype) = 0
553	IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) nextrsurf(jtype) = 2
554
555	!Read in surface obs types
556	CALL obs_rea_surf( surfdata(jtype), ifilessurf(jtype), &
557	& clsurffiles(jtype,1:ifilessurf(jtype)), &
558	& nvarssurf(jtype), nextrsurf(jtype), nitend-nit000+2, &
559	& rn_dobsini, rn_dobsend, ln_ignmis, .FALSE., llnightav(jtype) )
560
561	CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea, ln_bound_reject )
562
563	IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN
564	CALL obs_rea_mdt( surfdataqc(jtype), n2dintsurf(jtype) )
565	IF ( ln_altbias ) &
566	& CALL obs_rea_altbias ( surfdataqc(jtype), n2dintsurf(jtype), cn_altbiasfile )
567	ENDIF
568
569	IF ( TRIM(cobstypessurf(jtype)) == 'sst' .AND. ln_sstbias ) THEN
570	jnumsstbias = 0
571	DO jfile = 1, jpmaxnfiles
572	IF ( TRIM(cn_sstbiasfiles(jfile)) /= '' ) &
573	& jnumsstbias = jnumsstbias + 1
574	END DO
575	IF ( jnumsstbias == 0 ) THEN
576	CALL ctl_stop("ln_sstbias set but no bias files to read in")
577	ENDIF
578
579	CALL obs_app_sstbias( surfdataqc(jtype), n2dintsurf(jtype), &
580	& jnumsstbias, cn_sstbiasfiles(1:jnumsstbias) )
581
582	ENDIF
583
584	END DO
585
586	DEALLOCATE( ifilessurf, clsurffiles )
587
588	ENDIF
589
590	CALL wrk_dealloc( jpi, jpj, zglam1 )
591	CALL wrk_dealloc( jpi, jpj, zglam2 )
592	CALL wrk_dealloc( jpi, jpj, zgphi1 )
593	CALL wrk_dealloc( jpi, jpj, zgphi2 )
594	CALL wrk_dealloc( jpi, jpj, jpk, zmask1 )
595	CALL wrk_dealloc( jpi, jpj, jpk, zmask2 )
596
597	END SUBROUTINE dia_obs_init
598
599	SUBROUTINE dia_obs( kstp )
600	!!----------------------------------------------------------------------
601	!! * ROUTINE dia_obs *
602	!!
603	!! ** Purpose : Call the observation operators on each time step
604	!!
605	!! ** Method : Call the observation operators on each time step to
606	!! compute the model equivalent of the following data:
607	!! - Profile data, currently T/S or U/V
608	!! - Surface data, currently SST, SLA or sea-ice concentration.
609	!!
610	!! ** Action :
611	!!
612	!! History :
613	!! ! 06-03 (K. Mogensen) Original code
614	!! ! 06-05 (K. Mogensen) Reformatted
615	!! ! 06-10 (A. Weaver) Cleaning
616	!! ! 07-03 (K. Mogensen) General handling of profiles
617	!! ! 07-04 (G. Smith) Generalized surface operators
618	!! ! 08-10 (M. Valdivieso) obs operator for velocity profiles
619	!! ! 15-08 (M. Martin) Combined surface/profile routines.
620	!!----------------------------------------------------------------------
621	!! * Modules used
622	USE phycst, ONLY : & ! Physical constants
623	& rday
624	USE oce, ONLY : & ! Ocean dynamics and tracers variables
625	& tsn, &
626	& un, &
627	& vn, &
628	& sshn
629	#if defined key_lim3
630	USE ice, ONLY : & ! LIM3 Ice model variables
631	& frld
632	#endif
633	#if defined key_lim2
634	USE ice_2, ONLY : & ! LIM2 Ice model variables
635	& frld
636	#endif
637	#if defined key_cice
638	USE sbc_oce, ONLY : fr_i ! ice fraction
639	#endif
640	#if defined key_hadocc
641	USE trc, ONLY : & ! HadOCC chlorophyll, fCO2 and pCO2
642	& HADOCC_CHL, &
643	& HADOCC_FCO2, &
644	& HADOCC_PCO2, &
645	& HADOCC_FILL_FLT
646	#elif defined key_medusa && defined key_foam_medusa
647	USE trc, ONLY : & ! MEDUSA chlorophyll, fCO2 and pCO2
648	& MEDUSA_CHL, &
649	& MEDUSA_FCO2, &
650	& MEDUSA_PCO2, &
651	& MEDUSA_FILL_FLT
652	#elif defined key_fabm
653	USE fabm
654	USE par_fabm
655	#endif
656	#if defined key_spm
657	USE par_spm, ONLY: & ! ERSEM/SPM sediments
658	& jp_spm
659	USE trc, ONLY : &
660	& trn
661	#endif
662
663	IMPLICIT NONE
664
665	!! * Arguments
666	INTEGER, INTENT(IN) :: kstp ! Current timestep
667	!! * Local declarations
668	INTEGER :: idaystp ! Number of timesteps per day
669	INTEGER :: jtype ! Data loop variable
670	INTEGER :: jvar ! Variable number
671	INTEGER :: ji, jj ! Loop counters
672	REAL(wp) :: tiny ! small number
673	REAL(wp), POINTER, DIMENSION(:,:,:) :: &
674	& zprofvar1, & ! Model values for 1st variable in a prof ob
675	& zprofvar2 ! Model values for 2nd variable in a prof ob
676	REAL(wp), POINTER, DIMENSION(:,:,:) :: &
677	& zprofmask1, & ! Mask associated with zprofvar1
678	& zprofmask2 ! Mask associated with zprofvar2
679	REAL(wp), POINTER, DIMENSION(:,:) :: &
680	& zsurfvar, & ! Model values equivalent to surface ob.
681	& zsurfmask ! Mask associated with surface variable
682	REAL(wp), POINTER, DIMENSION(:,:) :: &
683	& zglam1, & ! Model longitudes for prof variable 1
684	& zglam2, & ! Model longitudes for prof variable 2
685	& zgphi1, & ! Model latitudes for prof variable 1
686	& zgphi2 ! Model latitudes for prof variable 2
687
688
689	!Allocate local work arrays
690	CALL wrk_alloc( jpi, jpj, jpk, zprofvar1 )
691	CALL wrk_alloc( jpi, jpj, jpk, zprofvar2 )
692	CALL wrk_alloc( jpi, jpj, jpk, zprofmask1 )
693	CALL wrk_alloc( jpi, jpj, jpk, zprofmask2 )
694	CALL wrk_alloc( jpi, jpj, zsurfvar )
695	CALL wrk_alloc( jpi, jpj, zsurfmask )
696	CALL wrk_alloc( jpi, jpj, zglam1 )
697	CALL wrk_alloc( jpi, jpj, zglam2 )
698	CALL wrk_alloc( jpi, jpj, zgphi1 )
699	CALL wrk_alloc( jpi, jpj, zgphi2 )
700
701	IF(lwp) THEN
702	WRITE(numout,*)
703	WRITE(numout,*) 'dia_obs : Call the observation operators', kstp
704	WRITE(numout,*) '~~~~~~~'
705	CALL FLUSH(numout)
706	ENDIF
707
708	idaystp = NINT( rday / rdt )
709
710	!-----------------------------------------------------------------------
711	! Call the profile and surface observation operators
712	!-----------------------------------------------------------------------
713
714	IF ( nproftypes > 0 ) THEN
715
716	DO jtype = 1, nproftypes
717
718	SELECT CASE ( TRIM(cobstypesprof(jtype)) )
719	CASE('prof')
720	zprofvar1(:,:,:) = tsn(:,:,:,jp_tem)
721	zprofvar2(:,:,:) = tsn(:,:,:,jp_sal)
722	zprofmask1(:,:,:) = tmask(:,:,:)
723	zprofmask2(:,:,:) = tmask(:,:,:)
724	zglam1(:,:) = glamt(:,:)
725	zglam2(:,:) = glamt(:,:)
726	zgphi1(:,:) = gphit(:,:)
727	zgphi2(:,:) = gphit(:,:)
728	CASE('vel')
729	zprofvar1(:,:,:) = un(:,:,:)
730	zprofvar2(:,:,:) = vn(:,:,:)
731	zprofmask1(:,:,:) = umask(:,:,:)
732	zprofmask2(:,:,:) = vmask(:,:,:)
733	zglam1(:,:) = glamu(:,:)
734	zglam2(:,:) = glamv(:,:)
735	zgphi1(:,:) = gphiu(:,:)
736	zgphi2(:,:) = gphiv(:,:)
737	CASE DEFAULT
738	CALL ctl_stop( 'Unknown profile observation type '//TRIM(cobstypesprof(jtype))//' in dia_obs' )
739	END SELECT
740
741	CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk, &
742	& nit000, idaystp, &
743	& zprofvar1, zprofvar2, &
744	& fsdept(:,:,:), fsdepw(:,:,:), &
745	& zprofmask1, zprofmask2, &
746	& zglam1, zglam2, zgphi1, zgphi2, &
747	& nn_1dint, nn_2dint, &
748	& kdailyavtypes = nn_profdavtypes )
749
750	END DO
751
752	ENDIF
753
754	IF ( nsurftypes > 0 ) THEN
755
756	DO jtype = 1, nsurftypes
757
758	!Defaults which might be changed
759	zsurfmask(:,:) = tmask(:,:,1)
760
761	SELECT CASE ( TRIM(cobstypessurf(jtype)) )
762	CASE('sst')
763	zsurfvar(:,:) = tsn(:,:,1,jp_tem)
764	CASE('sla')
765	zsurfvar(:,:) = sshn(:,:)
766	CASE('sss')
767	zsurfvar(:,:) = tsn(:,:,1,jp_sal)
768	CASE('sic')
769	IF ( kstp == 0 ) THEN
770	IF ( lwp .AND. surfdataqc(jtype)%nsstpmpp(1) > 0 ) THEN
771	CALL ctl_warn( 'Sea-ice not initialised on zeroth '// &
772	& 'time-step but some obs are valid then.' )
773	WRITE(numout,*)surfdataqc(jtype)%nsstpmpp(1), &
774	& ' sea-ice obs will be missed'
775	ENDIF
776	surfdataqc(jtype)%nsurfup = surfdataqc(jtype)%nsurfup + &
777	& surfdataqc(jtype)%nsstp(1)
778	CYCLE
779	ELSE
780	#if defined key_cice
781	zsurfvar(:,:) = fr_i(:,:)
782	#elif defined key_lim2 \|\| defined key_lim3
783	zsurfvar(:,:) = 1._wp - frld(:,:)
784	#else
785	CALL ctl_stop( ' Trying to run sea-ice observation operator', &
786	& ' but no sea-ice model appears to have been defined' )
787	#endif
788	ENDIF
789
790	CASE('logchl')
791	#if defined key_hadocc
792	zsurfvar(:,:) = HADOCC_CHL(:,:,1) ! (not log) chlorophyll from HadOCC
793	#elif defined key_medusa && defined key_foam_medusa
794	zsurfvar(:,:) = MEDUSA_CHL(:,:,1) ! (not log) chlorophyll from HadOCC
795	#elif defined key_fabm
796	chl_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabmdia_chltot)
797	zsurfvar(:,:) = chl_3d(:,:,1)
798	#else
799	CALL ctl_stop( ' Trying to run logchl observation operator', &
800	& ' but no biogeochemical model appears to have been defined' )
801	#endif
802	zsurfmask(:,:) = tmask(:,:,1) ! create a special mask to exclude certain things
803	! Take the log10 where we can, otherwise exclude
804	tiny = 1.0e-20
805	WHERE(zsurfvar(:,:) > tiny .AND. zsurfvar(:,:) /= obfillflt )
806	zsurfvar(:,:) = LOG10(zsurfvar(:,:))
807	ELSEWHERE
808	zsurfvar(:,:) = obfillflt
809	zsurfmask(:,:) = 0
810	END WHERE
811	CASE('spm')
812	#if defined key_spm
813	zsurfvar(:,:) = 0.0
814	DO jn = 1, jp_spm
815	zsurfvar(:,:) = zsurfvar(:,:) + trn(:,:,1,jn) ! sum SPM sizes
816	END DO
817	#else
818	CALL ctl_stop( ' Trying to run spm observation operator', &
819	& ' but no spm model appears to have been defined' )
820	#endif
821	CASE('fco2')
822	#if defined key_hadocc
823	zsurfvar(:,:) = HADOCC_FCO2(:,:) ! fCO2 from HadOCC
824	IF ( ( MINVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) .AND. &
825	& ( MAXVAL( HADOCC_FCO2 ) == HADOCC_FILL_FLT ) ) THEN
826	zsurfvar(:,:) = obfillflt
827	zsurfmask(:,:) = 0
828	CALL ctl_warn( ' HadOCC fCO2 values masked out for observation operator', &
829	& ' on timestep ' // TRIM(STR(kstp)), &
830	& ' as HADOCC_FCO2(:,:) == HADOCC_FILL_FLT' )
831	ENDIF
832	#elif defined key_medusa && defined key_foam_medusa
833	zsurfmask(:,:) = MEDUSA_FCO2(:,:) ! fCO2 from MEDUSA
834	IF ( ( MINVAL( MEDUSA_FCO2 ) == MEDUSA_FILL_FLT ) .AND. &
835	& ( MAXVAL( MEDUSA_FCO2 ) == MEDUSA_FILL_FLT ) ) THEN
836	zsurfvar(:,:) = obfillflt
837	zsurfmask(:,:) = 0
838	CALL ctl_warn( ' MEDUSA fCO2 values masked out for observation operator', &
839	& ' on timestep ' // TRIM(STR(kstp)), &
840	& ' as MEDUSA_FCO2(:,:) == MEDUSA_FILL_FLT' )
841	ENDIF
842	#elif defined key_fabm
843	! First, get pCO2 from FABM
844	pco2_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabm_o3pc)
845	zsurfvar(:,:) = pco2_3d(:,:,1)
846	! Now, convert pCO2 to fCO2, based on SST in K. This follows the standard methodology of:
847	! Pierrot et al. (2009), Recommendations for autonomous underway pCO2 measuring systems
848	! and data reduction routines, Deep-Sea Research II, 56: 512-522.
849	! and
850	! Weiss (1974), Carbon dioxide in water and seawater: the solubility of a non-ideal gas,
851	! Marine Chemistry, 2: 203-215.
852	! In the implementation below, atmospheric pressure has been assumed to be 1 atm and so
853	! not explicitly included - atmospheric pressure is not necessarily available so this is
854	! the best assumption.
855	! Further, the (1-xCO2)^2 term has been neglected. This is common practice
856	! (see e.g. Zeebe and Wolf-Gladrow (2001), CO2 in Seawater: Equilibrium, Kinetics, Isotopes)
857	! because xCO2 in atm is ~0, and so this term will only affect the result to the 3rd decimal
858	! place for typical values, and xCO2 would need to be approximated from pCO2 anyway.
859	zsurfvar(:,:) = zsurfvar(:,:) * EXP((-1636.75 + &
860	& 12.0408 * (tsn(:,:,1,jp_tem)+rt0) - &
861	& 0.0327957 * (tsn(:,:,1,jp_tem)+rt0)*(tsn(:,:,1,jp_tem)+rt0) + &
862	& 0.0000316528 * (tsn(:,:,1,jp_tem)+rt0)(tsn(:,:,1,jp_tem)+rt0)(tsn(:,:,1,jp_tem)+rt0) + &
863	& 2.0 * (57.7 - 0.118 * (tsn(:,:,1,jp_tem)+rt0))) / &
864	& (82.0578 * (tsn(:,:,1,jp_tem)+rt0)))
865	#else
866	CALL ctl_stop( ' Trying to run fco2 observation operator', &
867	& ' but no biogeochemical model appears to have been defined' )
868	#endif
869	CASE('pco2')
870	#if defined key_hadocc
871	zsurfvar(:,:) = HADOCC_PCO2(:,:) ! pCO2 from HadOCC
872	IF ( ( MINVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) .AND. &
873	& ( MAXVAL( HADOCC_PCO2 ) == HADOCC_FILL_FLT ) ) THEN
874	zsurfvar(:,:) = obfillflt
875	zsurfmask(:,:) = 0
876	CALL ctl_warn( ' HadOCC pCO2 values masked out for observation operator', &
877	& ' on timestep ' // TRIM(STR(kstp)), &
878	& ' as HADOCC_PCO2(:,:) == HADOCC_FILL_FLT' )
879	ENDIF
880	#elif defined key_medusa && defined key_foam_medusa
881	zsurfvar(:,:) = MEDUSA_PCO2(:,:) ! pCO2 from MEDUSA
882	IF ( ( MINVAL( MEDUSA_PCO2 ) == MEDUSA_FILL_FLT ) .AND. &
883	& ( MAXVAL( MEDUSA_PCO2 ) == MEDUSA_FILL_FLT ) ) THEN
884	zsurfvar(:,:) = obfillflt
885	zsurfmask(:,:) = 0
886	CALL ctl_warn( ' MEDUSA pCO2 values masked out for observation operator', &
887	& ' on timestep ' // TRIM(STR(kstp)), &
888	& ' as MEDUSA_PCO2(:,:) == MEDUSA_FILL_FLT' )
889	ENDIF
890	#elif defined key_fabm
891	pco2_3d(:,:,:) = fabm_get_bulk_diagnostic_data(model, jp_fabm_o3pc)
892	zsurfvar(:,:) = pco2_3d(:,:,1)
893	#else
894	CALL ctl_stop( ' Trying to run pCO2 observation operator', &
895	& ' but no biogeochemical model appears to have been defined' )
896	#endif
897
898	CASE DEFAULT
899
900	CALL ctl_stop( 'Unknown surface observation type '//TRIM(cobstypessurf(jtype))//' in dia_obs' )
901
902	END SELECT
903
904	CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj, &
905	& nit000, idaystp, zsurfvar, zsurfmask, &
906	& n2dintsurf(jtype), llnightav(jtype), &
907	& ravglamscl(jtype), ravgphiscl(jtype), &
908	& lfpindegs(jtype) )
909
910	END DO
911
912	ENDIF
913
914	CALL wrk_dealloc( jpi, jpj, jpk, zprofvar1 )
915	CALL wrk_dealloc( jpi, jpj, jpk, zprofvar2 )
916	CALL wrk_dealloc( jpi, jpj, jpk, zprofmask1 )
917	CALL wrk_dealloc( jpi, jpj, jpk, zprofmask2 )
918	CALL wrk_dealloc( jpi, jpj, zsurfvar )
919	CALL wrk_dealloc( jpi, jpj, zsurfmask )
920	CALL wrk_dealloc( jpi, jpj, zglam1 )
921	CALL wrk_dealloc( jpi, jpj, zglam2 )
922	CALL wrk_dealloc( jpi, jpj, zgphi1 )
923	CALL wrk_dealloc( jpi, jpj, zgphi2 )
924
925	END SUBROUTINE dia_obs
926
927	SUBROUTINE dia_obs_wri
928	!!----------------------------------------------------------------------
929	!! * ROUTINE dia_obs_wri *
930	!!
931	!! ** Purpose : Call observation diagnostic output routines
932	!!
933	!! ** Method : Call observation diagnostic output routines
934	!!
935	!! ** Action :
936	!!
937	!! History :
938	!! ! 06-03 (K. Mogensen) Original code
939	!! ! 06-05 (K. Mogensen) Reformatted
940	!! ! 06-10 (A. Weaver) Cleaning
941	!! ! 07-03 (K. Mogensen) General handling of profiles
942	!! ! 08-09 (M. Valdivieso) Velocity component (U,V) profiles
943	!! ! 15-08 (M. Martin) Combined writing for prof and surf types
944	!!----------------------------------------------------------------------
945	!! * Modules used
946	USE obs_rot_vel ! Rotation of velocities
947
948	IMPLICIT NONE
949
950	!! * Local declarations
951	INTEGER :: jtype ! Data set loop variable
952	INTEGER :: jo, jvar, jk
953	REAL(wp), DIMENSION(:), ALLOCATABLE :: &
954	& zu, &
955	& zv
956
957	!-----------------------------------------------------------------------
958	! Depending on switches call various observation output routines
959	!-----------------------------------------------------------------------
960
961	IF ( nproftypes > 0 ) THEN
962
963	DO jtype = 1, nproftypes
964
965	IF ( TRIM(cobstypesprof(jtype)) == 'vel' ) THEN
966
967	! For velocity data, rotate the model velocities to N/S, E/W
968	! using the compressed data structure.
969	ALLOCATE( &
970	& zu(profdataqc(jtype)%nvprot(1)), &
971	& zv(profdataqc(jtype)%nvprot(2)) &
972	& )
973
974	CALL obs_rotvel( profdataqc(jtype), nn_2dint, zu, zv )
975
976	DO jo = 1, profdataqc(jtype)%nprof
977	DO jvar = 1, 2
978	DO jk = profdataqc(jtype)%npvsta(jo,jvar), profdataqc(jtype)%npvend(jo,jvar)
979
980	IF ( jvar == 1 ) THEN
981	profdataqc(jtype)%var(jvar)%vmod(jk) = zu(jk)
982	ELSE
983	profdataqc(jtype)%var(jvar)%vmod(jk) = zv(jk)
984	ENDIF
985
986	END DO
987	END DO
988	END DO
989
990	DEALLOCATE( zu )
991	DEALLOCATE( zv )
992
993	END IF
994
995	CALL obs_prof_decompress( profdataqc(jtype), &
996	& profdata(jtype), .TRUE., numout )
997
998	CALL obs_wri_prof( profdata(jtype) )
999
1000	END DO
1001
1002	ENDIF
1003
1004	IF ( nsurftypes > 0 ) THEN
1005
1006	DO jtype = 1, nsurftypes
1007
1008	CALL obs_surf_decompress( surfdataqc(jtype), &
1009	& surfdata(jtype), .TRUE., numout )
1010
1011	CALL obs_wri_surf( surfdata(jtype) )
1012
1013	END DO
1014
1015	ENDIF
1016
1017	END SUBROUTINE dia_obs_wri
1018
1019	SUBROUTINE dia_obs_dealloc
1020	IMPLICIT NONE
1021	!!----------------------------------------------------------------------
1022	!! * ROUTINE dia_obs_dealloc *
1023	!!
1024	!! ** Purpose : To deallocate data to enable the obs_oper online loop.
1025	!! Specifically: dia_obs_init --> dia_obs --> dia_obs_wri
1026	!!
1027	!! ** Method : Clean up various arrays left behind by the obs_oper.
1028	!!
1029	!! ** Action :
1030	!!
1031	!!----------------------------------------------------------------------
1032	! obs_grid deallocation
1033	CALL obs_grid_deallocate
1034
1035	! diaobs deallocation
1036	IF ( nproftypes > 0 ) &
1037	& DEALLOCATE( cobstypesprof, profdata, profdataqc, nvarsprof, nextrprof )
1038
1039	IF ( nsurftypes > 0 ) &
1040	& DEALLOCATE( cobstypessurf, surfdata, surfdataqc, nvarssurf, nextrsurf, &
1041	& n2dintsurf, ravglamscl, ravgphiscl, lfpindegs, llnightav )
1042
1043	END SUBROUTINE dia_obs_dealloc
1044
1045	SUBROUTINE ini_date( ddobsini )
1046	!!----------------------------------------------------------------------
1047	!! * ROUTINE ini_date *
1048	!!
1049	!! ** Purpose : Get initial date in double precision YYYYMMDD.HHMMSS format
1050	!!
1051	!! ** Method : Get initial date in double precision YYYYMMDD.HHMMSS format
1052	!!
1053	!! ** Action : Get initial date in double precision YYYYMMDD.HHMMSS format
1054	!!
1055	!! History :
1056	!! ! 06-03 (K. Mogensen) Original code
1057	!! ! 06-05 (K. Mogensen) Reformatted
1058	!! ! 06-10 (A. Weaver) Cleaning
1059	!! ! 06-10 (G. Smith) Calculates initial date the same as method for final date
1060	!! ! 10-05 (D. Lea) Update to month length calculation for NEMO vn3.2
1061	!!----------------------------------------------------------------------
1062	USE phycst, ONLY : & ! Physical constants
1063	& rday
1064	USE dom_oce, ONLY : & ! Ocean space and time domain variables
1065	& rdt
1066
1067	IMPLICIT NONE
1068
1069	!! * Arguments
1070	REAL(dp), INTENT(OUT) :: ddobsini ! Initial date in YYYYMMDD.HHMMSS
1071
1072	!! * Local declarations
1073	INTEGER :: iyea ! date - (year, month, day, hour, minute)
1074	INTEGER :: imon
1075	INTEGER :: iday
1076	INTEGER :: ihou
1077	INTEGER :: imin
1078	INTEGER :: imday ! Number of days in month.
1079	INTEGER, DIMENSION(12) :: &
1080	& imonth_len ! Length in days of the months of the current year
1081	REAL(wp) :: zdayfrc ! Fraction of day
1082
1083	!----------------------------------------------------------------------
1084	! Initial date initialization (year, month, day, hour, minute)
1085	! (This assumes that the initial date is for 00z))
1086	!----------------------------------------------------------------------
1087	iyea = ndate0 / 10000
1088	imon = ( ndate0 - iyea * 10000 ) / 100
1089	iday = ndate0 - iyea * 10000 - imon * 100
1090	ihou = 0
1091	imin = 0
1092
1093	!----------------------------------------------------------------------
1094	! Compute number of days + number of hours + min since initial time
1095	!----------------------------------------------------------------------
1096	iday = iday + ( nit000 -1 ) * rdt / rday
1097	zdayfrc = ( nit000 -1 ) * rdt / rday
1098	zdayfrc = zdayfrc - aint(zdayfrc)
1099	ihou = int( zdayfrc * 24 )
1100	imin = int( (zdayfrc * 24 - ihou) * 60 )
1101
1102	!-----------------------------------------------------------------------
1103	! Convert number of days (iday) into a real date
1104	!----------------------------------------------------------------------
1105
1106	CALL calc_month_len( iyea, imonth_len )
1107
1108	DO WHILE ( iday > imonth_len(imon) )
1109	iday = iday - imonth_len(imon)
1110	imon = imon + 1
1111	IF ( imon > 12 ) THEN
1112	imon = 1
1113	iyea = iyea + 1
1114	CALL calc_month_len( iyea, imonth_len ) ! update month lengths
1115	ENDIF
1116	END DO
1117
1118	!----------------------------------------------------------------------
1119	! Convert it into YYYYMMDD.HHMMSS format.
1120	!----------------------------------------------------------------------
1121	ddobsini = iyea * 10000_dp + imon * 100_dp + &
1122	& iday + ihou * 0.01_dp + imin * 0.0001_dp
1123
1124
1125	END SUBROUTINE ini_date
1126
1127	SUBROUTINE fin_date( ddobsfin )
1128	!!----------------------------------------------------------------------
1129	!! * ROUTINE fin_date *
1130	!!
1131	!! ** Purpose : Get final date in double precision YYYYMMDD.HHMMSS format
1132	!!
1133	!! ** Method : Get final date in double precision YYYYMMDD.HHMMSS format
1134	!!
1135	!! ** Action : Get final date in double precision YYYYMMDD.HHMMSS format
1136	!!
1137	!! History :
1138	!! ! 06-03 (K. Mogensen) Original code
1139	!! ! 06-05 (K. Mogensen) Reformatted
1140	!! ! 06-10 (A. Weaver) Cleaning
1141	!! ! 10-05 (D. Lea) Update to month length calculation for NEMO vn3.2
1142	!!----------------------------------------------------------------------
1143	USE phycst, ONLY : & ! Physical constants
1144	& rday
1145	USE dom_oce, ONLY : & ! Ocean space and time domain variables
1146	& rdt
1147
1148	IMPLICIT NONE
1149
1150	!! * Arguments
1151	REAL(dp), INTENT(OUT) :: ddobsfin ! Final date in YYYYMMDD.HHMMSS
1152
1153	!! * Local declarations
1154	INTEGER :: iyea ! date - (year, month, day, hour, minute)
1155	INTEGER :: imon
1156	INTEGER :: iday
1157	INTEGER :: ihou
1158	INTEGER :: imin
1159	INTEGER :: imday ! Number of days in month.
1160	INTEGER, DIMENSION(12) :: &
1161	& imonth_len ! Length in days of the months of the current year
1162	REAL(wp) :: zdayfrc ! Fraction of day
1163
1164	!-----------------------------------------------------------------------
1165	! Initial date initialization (year, month, day, hour, minute)
1166	! (This assumes that the initial date is for 00z)
1167	!-----------------------------------------------------------------------
1168	iyea = ndate0 / 10000
1169	imon = ( ndate0 - iyea * 10000 ) / 100
1170	iday = ndate0 - iyea * 10000 - imon * 100
1171	ihou = 0
1172	imin = 0
1173
1174	!-----------------------------------------------------------------------
1175	! Compute number of days + number of hours + min since initial time
1176	!-----------------------------------------------------------------------
1177	iday = iday + nitend * rdt / rday
1178	zdayfrc = nitend * rdt / rday
1179	zdayfrc = zdayfrc - AINT( zdayfrc )
1180	ihou = INT( zdayfrc * 24 )
1181	imin = INT( ( zdayfrc * 24 - ihou ) * 60 )
1182
1183	!-----------------------------------------------------------------------
1184	! Convert number of days (iday) into a real date
1185	!----------------------------------------------------------------------
1186
1187	CALL calc_month_len( iyea, imonth_len )
1188
1189	DO WHILE ( iday > imonth_len(imon) )
1190	iday = iday - imonth_len(imon)
1191	imon = imon + 1
1192	IF ( imon > 12 ) THEN
1193	imon = 1
1194	iyea = iyea + 1
1195	CALL calc_month_len( iyea, imonth_len ) ! update month lengths
1196	ENDIF
1197	END DO
1198
1199	!-----------------------------------------------------------------------
1200	! Convert it into YYYYMMDD.HHMMSS format
1201	!-----------------------------------------------------------------------
1202	ddobsfin = iyea * 10000_dp + imon * 100_dp + iday &
1203	& + ihou * 0.01_dp + imin * 0.0001_dp
1204
1205	END SUBROUTINE fin_date
1206
1207	SUBROUTINE obs_settypefiles( ntypes, jpmaxnfiles, jtype, ctypein, &
1208	& cfilestype, ifiles, cobstypes, cfiles )
1209
1210	INTEGER, INTENT(IN) :: ntypes ! Total number of obs types
1211	INTEGER, INTENT(IN) :: jpmaxnfiles ! Maximum number of files allowed for each type
1212	INTEGER, INTENT(IN) :: jtype ! Index of the current type of obs
1213	INTEGER, DIMENSION(ntypes), INTENT(INOUT) :: &
1214	& ifiles ! Out appended number of files for this type
1215
1216	CHARACTER(len=6), INTENT(IN) :: ctypein
1217	CHARACTER(len=128), DIMENSION(jpmaxnfiles), INTENT(IN) :: &
1218	& cfilestype ! In list of files for this obs type
1219	CHARACTER(len=6), DIMENSION(ntypes), INTENT(INOUT) :: &
1220	& cobstypes ! Out appended list of obs types
1221	CHARACTER(len=128), DIMENSION(ntypes, jpmaxnfiles), INTENT(INOUT) :: &
1222	& cfiles ! Out appended list of files for all types
1223
1224	!Local variables
1225	INTEGER :: jfile
1226
1227	cfiles(jtype,:) = cfilestype(:)
1228	cobstypes(jtype) = ctypein
1229	ifiles(jtype) = 0
1230	DO jfile = 1, jpmaxnfiles
1231	IF ( trim(cfiles(jtype,jfile)) /= '' ) &
1232	ifiles(jtype) = ifiles(jtype) + 1
1233	END DO
1234
1235	IF ( ifiles(jtype) == 0 ) THEN
1236	CALL ctl_stop( 'Logical for observation type '//TRIM(ctypein)// &
1237	& ' set to true but no files available to read' )
1238	ENDIF
1239
1240	IF(lwp) THEN
1241	WRITE(numout,*) ' '//cobstypes(jtype)//' input observation file names:'
1242	DO jfile = 1, ifiles(jtype)
1243	WRITE(numout,*) ' '//TRIM(cfiles(jtype,jfile))
1244	END DO
1245	ENDIF
1246
1247	END SUBROUTINE obs_settypefiles
1248
1249	SUBROUTINE obs_setinterpopts( ntypes, jtype, ctypein, &
1250	& n2dint_default, n2dint_type, &
1251	& ravglamscl_type, ravgphiscl_type, &
1252	& lfp_indegs_type, lavnight_type, &
1253	& n2dint, ravglamscl, ravgphiscl, &
1254	& lfpindegs, lavnight )
1255
1256	INTEGER, INTENT(IN) :: ntypes ! Total number of obs types
1257	INTEGER, INTENT(IN) :: jtype ! Index of the current type of obs
1258	INTEGER, INTENT(IN) :: n2dint_default ! Default option for interpolation type
1259	INTEGER, INTENT(IN) :: n2dint_type ! Option for interpolation type
1260	REAL(wp), INTENT(IN) :: &
1261	& ravglamscl_type, & !E/W diameter of obs footprint for this type
1262	& ravgphiscl_type !N/S diameter of obs footprint for this type
1263	LOGICAL, INTENT(IN) :: lfp_indegs_type !T=> footprint in degrees, F=> in metres
1264	LOGICAL, INTENT(IN) :: lavnight_type !T=> obs represent night time average
1265	CHARACTER(len=6), INTENT(IN) :: ctypein
1266
1267	INTEGER, DIMENSION(ntypes), INTENT(INOUT) :: &
1268	& n2dint
1269	REAL(wp), DIMENSION(ntypes), INTENT(INOUT) :: &
1270	& ravglamscl, ravgphiscl
1271	LOGICAL, DIMENSION(ntypes), INTENT(INOUT) :: &
1272	& lfpindegs, lavnight
1273
1274	lavnight(jtype) = lavnight_type
1275
1276	IF ( (n2dint_type >= 1) .AND. (n2dint_type <= 6) ) THEN
1277	n2dint(jtype) = n2dint_type
1278	ELSE
1279	n2dint(jtype) = n2dint_default
1280	ENDIF
1281
1282	! For averaging observation footprints set options for size of footprint
1283	IF ( (n2dint(jtype) > 4) .AND. (n2dint(jtype) <= 6) ) THEN
1284	IF ( ravglamscl_type > 0._wp ) THEN
1285	ravglamscl(jtype) = ravglamscl_type
1286	ELSE
1287	CALL ctl_stop( 'Incorrect value set for averaging footprint '// &
1288	'scale (ravglamscl) for observation type '//TRIM(ctypein) )
1289	ENDIF
1290
1291	IF ( ravgphiscl_type > 0._wp ) THEN
1292	ravgphiscl(jtype) = ravgphiscl_type
1293	ELSE
1294	CALL ctl_stop( 'Incorrect value set for averaging footprint '// &
1295	'scale (ravgphiscl) for observation type '//TRIM(ctypein) )
1296	ENDIF
1297
1298	lfpindegs(jtype) = lfp_indegs_type
1299
1300	ENDIF
1301
1302	! Write out info
1303	IF(lwp) THEN
1304	IF ( n2dint(jtype) <= 4 ) THEN
1305	WRITE(numout,*) ' '//TRIM(ctypein)// &
1306	& ' model counterparts will be interpolated horizontally'
1307	ELSE IF ( n2dint(jtype) <= 6 ) THEN
1308	WRITE(numout,*) ' '//TRIM(ctypein)// &
1309	& ' model counterparts will be averaged horizontally'
1310	WRITE(numout,*) ' '//' with E/W scale: ',ravglamscl(jtype)
1311	WRITE(numout,*) ' '//' with N/S scale: ',ravgphiscl(jtype)
1312	IF ( lfpindegs(jtype) ) THEN
1313	WRITE(numout,*) ' '//' (in degrees)'
1314	ELSE
1315	WRITE(numout,*) ' '//' (in metres)'
1316	ENDIF
1317	ENDIF
1318	ENDIF
1319
1320	END SUBROUTINE obs_setinterpopts
1321
1322	END MODULE diaobs

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: branches/UKMO/obs_oper_do_not_assim_update/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 7915

Download in other formats: