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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

diaobs.F90 in branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

Context Navigation

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 9023

Last change on this file since 9023 was 9023, checked in by timgraham, 7 years ago
Merged METO_MERCATOR branch and resolved all conflicts in OPA_SRC
Property svn:keywords set to `Id`
File size: 48.2 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	!! calc_date : Compute the date of timestep in YYYYMMDD.HHMMSS format
13	!! ini_date : Compute the initial date YYYYMMDD.HHMMSS
14	!! fin_date : Compute the final date YYYYMMDD.HHMMSS
15	!!----------------------------------------------------------------------
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_sstbias ! Bias correction routine for SST
24	USE obs_readmdt ! Reading and allocation of MDT for SLA.
25	USE obs_prep ! Preparation of obs. (grid search etc).
26	USE obs_oper ! Observation operators
27	USE obs_write ! Writing of observation related diagnostics
28	USE obs_grid ! Grid searching
29	USE obs_read_altbias ! Bias treatment for altimeter
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	PRIVATE
38
39	PUBLIC dia_obs_init ! Initialize and read observations
40	PUBLIC dia_obs ! Compute model equivalent to observations
41	PUBLIC dia_obs_wri ! Write model equivalent to observations
42	PUBLIC dia_obs_dealloc ! Deallocate dia_obs data
43	PUBLIC calc_date ! Compute the date of a timestep
44
45	!! * Module variables
46	LOGICAL, PUBLIC :: ln_diaobs !: Logical switch for the obs operator
47	LOGICAL :: ln_sstnight !: Logical switch for night mean SST obs
48	LOGICAL :: ln_sla_fp_indegs !: T=> SLA obs footprint size specified in degrees, F=> in metres
49	LOGICAL :: ln_sst_fp_indegs !: T=> SST obs footprint size specified in degrees, F=> in metres
50	LOGICAL :: ln_sss_fp_indegs !: T=> SSS obs footprint size specified in degrees, F=> in metres
51	LOGICAL :: ln_sic_fp_indegs !: T=> sea-ice obs footprint size specified in degrees, F=> in metres
52
53	REAL(wp) :: rn_sla_avglamscl !: E/W diameter of SLA observation footprint (metres)
54	REAL(wp) :: rn_sla_avgphiscl !: N/S diameter of SLA observation footprint (metres)
55	REAL(wp) :: rn_sst_avglamscl !: E/W diameter of SST observation footprint (metres)
56	REAL(wp) :: rn_sst_avgphiscl !: N/S diameter of SST observation footprint (metres)
57	REAL(wp) :: rn_sss_avglamscl !: E/W diameter of SSS observation footprint (metres)
58	REAL(wp) :: rn_sss_avgphiscl !: N/S diameter of SSS observation footprint (metres)
59	REAL(wp) :: rn_sic_avglamscl !: E/W diameter of sea-ice observation footprint (metres)
60	REAL(wp) :: rn_sic_avgphiscl !: N/S diameter of sea-ice observation footprint (metres)
61
62	INTEGER :: nn_1dint !: Vertical interpolation method
63	INTEGER :: nn_2dint !: Default horizontal interpolation method
64	INTEGER :: nn_2dint_sla !: SLA horizontal interpolation method
65	INTEGER :: nn_2dint_sst !: SST horizontal interpolation method
66	INTEGER :: nn_2dint_sss !: SSS horizontal interpolation method
67	INTEGER :: nn_2dint_sic !: Seaice horizontal interpolation method
68	INTEGER, DIMENSION(imaxavtypes) :: nn_profdavtypes !: Profile data types representing a daily average
69	INTEGER :: nproftypes !: Number of profile obs types
70	INTEGER :: nsurftypes !: Number of surface obs types
71	INTEGER, DIMENSION(:), ALLOCATABLE :: &
72	& nvarsprof, & !: Number of profile variables
73	& nvarssurf !: Number of surface variables
74	INTEGER, DIMENSION(:), ALLOCATABLE :: &
75	& nextrprof, & !: Number of profile extra variables
76	& nextrsurf !: Number of surface extra variables
77	INTEGER, DIMENSION(:), ALLOCATABLE :: &
78	& n2dintsurf !: Interpolation option for surface variables
79	REAL(wp), DIMENSION(:), ALLOCATABLE :: &
80	& zavglamscl, & !: E/W diameter of averaging footprint for surface variables
81	& zavgphiscl !: N/S diameter of averaging footprint for surface variables
82	LOGICAL, DIMENSION(:), ALLOCATABLE :: &
83	& lfpindegs, & !: T=> surface obs footprint size specified in degrees, F=> in metres
84	& llnightav !: Logical for calculating night-time averages
85
86	TYPE(obs_surf), PUBLIC, POINTER, DIMENSION(:) :: &
87	& surfdata, & !: Initial surface data
88	& surfdataqc !: Surface data after quality control
89	TYPE(obs_prof), PUBLIC, POINTER, DIMENSION(:) :: &
90	& profdata, & !: Initial profile data
91	& profdataqc !: Profile data after quality control
92
93	CHARACTER(len=6), PUBLIC, DIMENSION(:), ALLOCATABLE :: cobstypesprof, cobstypessurf !: Profile & surface obs types
94
95	!!----------------------------------------------------------------------
96	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
97	!! $Id$
98	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
99	!!----------------------------------------------------------------------
100	CONTAINS
101
102	SUBROUTINE dia_obs_init
103	!!----------------------------------------------------------------------
104	!! * ROUTINE dia_obs_init *
105	!!
106	!! ** Purpose : Initialize and read observations
107	!!
108	!! ** Method : Read the namelist and call reading routines
109	!!
110	!! ** Action : Read the namelist and call reading routines
111	!!
112	!! History :
113	!! ! 06-03 (K. Mogensen) Original code
114	!! ! 06-05 (A. Weaver) Reformatted
115	!! ! 06-10 (A. Weaver) Cleaning and add controls
116	!! ! 07-03 (K. Mogensen) General handling of profiles
117	!! ! 14-08 (J.While) Incorporated SST bias correction
118	!! ! 15-02 (M. Martin) Simplification of namelist and code
119	!!----------------------------------------------------------------------
120	INTEGER, PARAMETER :: jpmaxnfiles = 1000 ! Maximum number of files for each obs type
121	INTEGER, DIMENSION(:), ALLOCATABLE :: ifilesprof, ifilessurf ! Number of profile & surface files
122	INTEGER :: ios ! Local integer output status for namelist read
123	INTEGER :: jtype ! Counter for obs types
124	INTEGER :: jvar ! Counter for variables
125	INTEGER :: jfile ! Counter for files
126
127	CHARACTER(len=128), DIMENSION(jpmaxnfiles) :: &
128	& cn_profbfiles, & ! T/S profile input filenames
129	& cn_sstfbfiles, & ! Sea surface temperature input filenames
130	& cn_sssfbfiles, & ! Sea surface salinity input filenames
131	& cn_slafbfiles, & ! Sea level anomaly input filenames
132	& cn_sicfbfiles, & ! Seaice concentration input filenames
133	& cn_velfbfiles, & ! Velocity profile input filenames
134	& cn_sstbiasfiles ! SST bias input filenames
135	CHARACTER(LEN=128) :: &
136	& cn_altbiasfile ! Altimeter bias input filename
137	CHARACTER(len=128), DIMENSION(:,:), ALLOCATABLE :: &
138	& clproffiles, & ! Profile filenames
139	& clsurffiles ! Surface filenames
140
141	LOGICAL :: ln_t3d ! Logical switch for temperature profiles
142	LOGICAL :: ln_s3d ! Logical switch for salinity profiles
143	LOGICAL :: ln_sla ! Logical switch for sea level anomalies
144	LOGICAL :: ln_sst ! Logical switch for sea surface temperature
145	LOGICAL :: ln_sss ! Logical switch for sea surface salinity
146	LOGICAL :: ln_sic ! Logical switch for sea ice concentration
147	LOGICAL :: ln_vel3d ! Logical switch for velocity (u,v) obs
148	LOGICAL :: ln_nea ! Logical switch to remove obs near land
149	LOGICAL :: ln_altbias ! Logical switch for altimeter bias
150	LOGICAL :: ln_sstbias ! Logical switch for bias corection of SST
151	LOGICAL :: ln_ignmis ! Logical switch for ignoring missing files
152	LOGICAL :: ln_s_at_t ! Logical switch to compute model S at T obs
153	LOGICAL :: ln_bound_reject ! Logical to remove obs near boundaries in LAMs.
154	LOGICAL :: llvar1 ! Logical for profile variable 1
155	LOGICAL :: llvar2 ! Logical for profile variable 1
156	LOGICAL, DIMENSION(jpmaxnfiles) :: lmask ! Used for finding number of sstbias files
157
158	REAL(dp) :: rn_dobsini ! Obs window start date YYYYMMDD.HHMMSS
159	REAL(dp) :: rn_dobsend ! Obs window end date YYYYMMDD.HHMMSS
160	REAL(wp), POINTER, DIMENSION(:,:) :: &
161	& zglam1, & ! Model longitudes for profile variable 1
162	& zglam2 ! Model longitudes for profile variable 2
163	REAL(wp), POINTER, DIMENSION(:,:) :: &
164	& zgphi1, & ! Model latitudes for profile variable 1
165	& zgphi2 ! Model latitudes for profile variable 2
166	REAL(wp), POINTER, DIMENSION(:,:,:) :: &
167	& zmask1, & ! Model land/sea mask associated with variable 1
168	& zmask2 ! Model land/sea mask associated with variable 2
169
170	NAMELIST/namobs/ln_diaobs, ln_t3d, ln_s3d, ln_sla, &
171	& ln_sst, ln_sic, ln_sss, ln_vel3d, &
172	& ln_altbias, ln_sstbias, ln_nea, &
173	& ln_grid_global, ln_grid_search_lookup, &
174	& ln_ignmis, ln_s_at_t, ln_bound_reject, &
175	& ln_sstnight, &
176	& ln_sla_fp_indegs, ln_sst_fp_indegs, &
177	& ln_sss_fp_indegs, ln_sic_fp_indegs, &
178	& cn_profbfiles, cn_slafbfiles, &
179	& cn_sstfbfiles, cn_sicfbfiles, &
180	& cn_velfbfiles, cn_sssfbfiles, &
181	& cn_sstbiasfiles, cn_altbiasfile, &
182	& cn_gridsearchfile, rn_gridsearchres, &
183	& rn_dobsini, rn_dobsend, &
184	& rn_sla_avglamscl, rn_sla_avgphiscl, &
185	& rn_sst_avglamscl, rn_sst_avgphiscl, &
186	& rn_sss_avglamscl, rn_sss_avgphiscl, &
187	& rn_sic_avglamscl, rn_sic_avgphiscl, &
188	& nn_1dint, nn_2dint, &
189	& nn_2dint_sla, nn_2dint_sst, &
190	& nn_2dint_sss, nn_2dint_sic, &
191	& nn_msshc, rn_mdtcorr, rn_mdtcutoff, &
192	& nn_profdavtypes
193
194	INTEGER :: jnumsstbias
195	CALL wrk_alloc( jpi, jpj, zglam1 )
196	CALL wrk_alloc( jpi, jpj, zglam2 )
197	CALL wrk_alloc( jpi, jpj, zgphi1 )
198	CALL wrk_alloc( jpi, jpj, zgphi2 )
199	CALL wrk_alloc( jpi, jpj, jpk, zmask1 )
200	CALL wrk_alloc( jpi, jpj, jpk, zmask2 )
201
202	!-----------------------------------------------------------------------
203	! Read namelist parameters
204	!-----------------------------------------------------------------------
205	! Some namelist arrays need initialising
206	cn_profbfiles(:) = ''
207	cn_slafbfiles(:) = ''
208	cn_sstfbfiles(:) = ''
209	cn_sicfbfiles(:) = ''
210	cn_velfbfiles(:) = ''
211	cn_sssfbfiles(:) = ''
212	cn_sstbiasfiles(:) = ''
213	nn_profdavtypes(:) = -1
214
215	CALL ini_date( rn_dobsini )
216	CALL fin_date( rn_dobsend )
217
218	! Read namelist namobs : control observation diagnostics
219	REWIND( numnam_ref ) ! Namelist namobs in reference namelist
220	READ ( numnam_ref, namobs, IOSTAT = ios, ERR = 901)
221	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in reference namelist', lwp )
222
223	REWIND( numnam_cfg ) ! Namelist namobs in configuration namelist
224	READ ( numnam_cfg, namobs, IOSTAT = ios, ERR = 902 )
225	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in configuration namelist', lwp )
226	IF(lwm) WRITE ( numond, namobs )
227
228	IF ( .NOT. ln_diaobs ) THEN
229	IF(lwp) WRITE(numout,cform_war)
230	IF(lwp) WRITE(numout,*)' ln_diaobs is set to false so not calling dia_obs'
231	RETURN
232	ENDIF
233
234	IF(lwp) THEN
235	WRITE(numout,*)
236	WRITE(numout,*) 'dia_obs_init : Observation diagnostic initialization'
237	WRITE(numout,*) '~~~~~~~~~~~~'
238	WRITE(numout,*) ' Namelist namobs : set observation diagnostic parameters'
239	WRITE(numout,*) ' Logical switch for T profile observations ln_t3d = ', ln_t3d
240	WRITE(numout,*) ' Logical switch for S profile observations ln_s3d = ', ln_s3d
241	WRITE(numout,*) ' Logical switch for SLA observations ln_sla = ', ln_sla
242	WRITE(numout,*) ' Logical switch for SST observations ln_sst = ', ln_sst
243	WRITE(numout,*) ' Logical switch for Sea Ice observations ln_sic = ', ln_sic
244	WRITE(numout,*) ' Logical switch for velocity observations ln_vel3d = ', ln_vel3d
245	WRITE(numout,*) ' Logical switch for SSS observations ln_sss = ', ln_sss
246	WRITE(numout,*) ' Global distribution of observations ln_grid_global = ', ln_grid_global
247	WRITE(numout,*) ' Logical switch for obs grid search lookup ln_grid_search_lookup = ', ln_grid_search_lookup
248	IF (ln_grid_search_lookup) &
249	WRITE(numout,*) ' Grid search lookup file header cn_gridsearchfile = ', cn_gridsearchfile
250	WRITE(numout,*) ' Initial date in window YYYYMMDD.HHMMSS rn_dobsini = ', rn_dobsini
251	WRITE(numout,*) ' Final date in window YYYYMMDD.HHMMSS rn_dobsend = ', rn_dobsend
252	WRITE(numout,*) ' Type of vertical interpolation method nn_1dint = ', nn_1dint
253	WRITE(numout,*) ' Type of horizontal interpolation method nn_2dint = ', nn_2dint
254	WRITE(numout,*) ' Rejection of observations near land switch ln_nea = ', ln_nea
255	WRITE(numout,*) ' Rejection of obs near open bdys ln_bound_reject = ', ln_bound_reject
256	WRITE(numout,*) ' MSSH correction scheme nn_msshc = ', nn_msshc
257	WRITE(numout,*) ' MDT correction rn_mdtcorr = ', rn_mdtcorr
258	WRITE(numout,*) ' MDT cutoff for computed correction rn_mdtcutoff = ', rn_mdtcutoff
259	WRITE(numout,*) ' Logical switch for alt bias ln_altbias = ', ln_altbias
260	WRITE(numout,*) ' Logical switch for sst bias ln_sstbias = ', ln_sstbias
261	WRITE(numout,*) ' Logical switch for ignoring missing files ln_ignmis = ', ln_ignmis
262	WRITE(numout,*) ' Daily average types nn_profdavtypes = ', nn_profdavtypes
263	WRITE(numout,*) ' Logical switch for night-time SST obs ln_sstnight = ', ln_sstnight
264	ENDIF
265	!-----------------------------------------------------------------------
266	! Set up list of observation types to be used
267	! and the files associated with each type
268	!-----------------------------------------------------------------------
269
270	nproftypes = COUNT( (/ln_t3d .OR. ln_s3d, ln_vel3d /) )
271	nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic, ln_sss /) )
272
273	IF (ln_sstbias) THEN
274	lmask(:) = .FALSE.
275	WHERE (cn_sstbiasfiles(:) /= '') lmask(:) = .TRUE.
276	jnumsstbias = COUNT(lmask)
277	lmask(:) = .FALSE.
278	ENDIF
279
280	IF ( nproftypes == 0 .AND. nsurftypes == 0 ) THEN
281	IF(lwp) WRITE(numout,cform_war)
282	IF(lwp) WRITE(numout,*) ' ln_diaobs is set to true, but all obs operator logical flags', &
283	& ' ln_t3d, ln_s3d, ln_sla, ln_sst, ln_sic, ln_vel3d', &
284	& ' are set to .FALSE. so turning off calls to dia_obs'
285	nwarn = nwarn + 1
286	ln_diaobs = .FALSE.
287	RETURN
288	ENDIF
289
290	IF ( nproftypes > 0 ) THEN
291
292	ALLOCATE( cobstypesprof(nproftypes) )
293	ALLOCATE( ifilesprof(nproftypes) )
294	ALLOCATE( clproffiles(nproftypes,jpmaxnfiles) )
295
296	jtype = 0
297	IF (ln_t3d .OR. ln_s3d) THEN
298	jtype = jtype + 1
299	CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'prof ', &
300	& cn_profbfiles, ifilesprof, cobstypesprof, clproffiles )
301	ENDIF
302	IF (ln_vel3d) THEN
303	jtype = jtype + 1
304	CALL obs_settypefiles( nproftypes, jpmaxnfiles, jtype, 'vel ', &
305	& cn_velfbfiles, ifilesprof, cobstypesprof, clproffiles )
306	ENDIF
307
308	ENDIF
309
310	IF ( nsurftypes > 0 ) THEN
311
312	ALLOCATE( cobstypessurf(nsurftypes) )
313	ALLOCATE( ifilessurf(nsurftypes) )
314	ALLOCATE( clsurffiles(nsurftypes, jpmaxnfiles) )
315	ALLOCATE(n2dintsurf(nsurftypes))
316	ALLOCATE(zavglamscl(nsurftypes))
317	ALLOCATE(zavgphiscl(nsurftypes))
318	ALLOCATE(lfpindegs(nsurftypes))
319	ALLOCATE(llnightav(nsurftypes))
320
321	jtype = 0
322	IF (ln_sla) THEN
323	jtype = jtype + 1
324	CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sla ', &
325	& cn_slafbfiles, ifilessurf, cobstypessurf, clsurffiles )
326	CALL obs_setinterpopts( nsurftypes, jtype, 'sla ', &
327	& nn_2dint, nn_2dint_sla, &
328	& rn_sla_avglamscl, rn_sla_avgphiscl, &
329	& ln_sla_fp_indegs, .FALSE., &
330	& n2dintsurf, zavglamscl, zavgphiscl, &
331	& lfpindegs, llnightav )
332	ENDIF
333	IF (ln_sst) THEN
334	jtype = jtype + 1
335	CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sst ', &
336	& cn_sstfbfiles, ifilessurf, cobstypessurf, clsurffiles )
337	CALL obs_setinterpopts( nsurftypes, jtype, 'sst ', &
338	& nn_2dint, nn_2dint_sst, &
339	& rn_sst_avglamscl, rn_sst_avgphiscl, &
340	& ln_sst_fp_indegs, ln_sstnight, &
341	& n2dintsurf, zavglamscl, zavgphiscl, &
342	& lfpindegs, llnightav )
343	ENDIF
344	#if defined key_lim3 \|\| defined key_cice
345	IF (ln_sic) THEN
346	jtype = jtype + 1
347	CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sic ', &
348	& cn_sicfbfiles, ifilessurf, cobstypessurf, clsurffiles )
349	CALL obs_setinterpopts( nsurftypes, jtype, 'sic ', &
350	& nn_2dint, nn_2dint_sic, &
351	& rn_sic_avglamscl, rn_sic_avgphiscl, &
352	& ln_sic_fp_indegs, .FALSE., &
353	& n2dintsurf, zavglamscl, zavgphiscl, &
354	& lfpindegs, llnightav )
355	ENDIF
356	#endif
357	IF (ln_sss) THEN
358	jtype = jtype + 1
359	CALL obs_settypefiles( nsurftypes, jpmaxnfiles, jtype, 'sss ', &
360	& cn_sssfbfiles, ifilessurf, cobstypessurf, clsurffiles )
361	CALL obs_setinterpopts( nsurftypes, jtype, 'sss ', &
362	& nn_2dint, nn_2dint_sss, &
363	& rn_sss_avglamscl, rn_sss_avgphiscl, &
364	& ln_sss_fp_indegs, .FALSE., &
365	& n2dintsurf, zavglamscl, zavgphiscl, &
366	& lfpindegs, llnightav )
367	ENDIF
368
369	ENDIF
370
371
372
373	!-----------------------------------------------------------------------
374	! Obs operator parameter checking and initialisations
375	!-----------------------------------------------------------------------
376
377	IF ( ln_vel3d .AND. ( .NOT. ln_grid_global ) ) THEN
378	CALL ctl_stop( 'Velocity data only works with ln_grid_global=.true.' )
379	RETURN
380	ENDIF
381
382	IF ( ln_grid_global ) THEN
383	CALL ctl_warn( 'ln_grid_global=T may cause memory issues when used with a large number of processors' )
384	ENDIF
385
386	IF ( ( nn_1dint < 0 ) .OR. ( nn_1dint > 1 ) ) THEN
387	CALL ctl_stop(' Choice of vertical (1D) interpolation method', &
388	& ' is not available')
389	ENDIF
390
391	IF ( ( nn_2dint < 0 ) .OR. ( nn_2dint > 6 ) ) THEN
392	CALL ctl_stop(' Choice of horizontal (2D) interpolation method', &
393	& ' is not available')
394	ENDIF
395
396	CALL obs_typ_init
397	IF(ln_grid_global) THEN
398	CALL mppmap_init
399	ENDIF
400
401	CALL obs_grid_setup( )
402
403	!-----------------------------------------------------------------------
404	! Depending on switches read the various observation types
405	!-----------------------------------------------------------------------
406
407	IF ( nproftypes > 0 ) THEN
408
409	ALLOCATE(profdata(nproftypes))
410	ALLOCATE(profdataqc(nproftypes))
411	ALLOCATE(nvarsprof(nproftypes))
412	ALLOCATE(nextrprof(nproftypes))
413
414	DO jtype = 1, nproftypes
415
416	nvarsprof(jtype) = 2
417	IF ( TRIM(cobstypesprof(jtype)) == 'prof' ) THEN
418	nextrprof(jtype) = 1
419	llvar1 = ln_t3d
420	llvar2 = ln_s3d
421	zglam1 = glamt
422	zgphi1 = gphit
423	zmask1 = tmask
424	zglam2 = glamt
425	zgphi2 = gphit
426	zmask2 = tmask
427	ENDIF
428	IF ( TRIM(cobstypesprof(jtype)) == 'vel' ) THEN
429	nextrprof(jtype) = 2
430	llvar1 = ln_vel3d
431	llvar2 = ln_vel3d
432	zglam1 = glamu
433	zgphi1 = gphiu
434	zmask1 = umask
435	zglam2 = glamv
436	zgphi2 = gphiv
437	zmask2 = vmask
438	ENDIF
439
440	!Read in profile or profile obs types
441	CALL obs_rea_prof( profdata(jtype), ifilesprof(jtype), &
442	& clproffiles(jtype,1:ifilesprof(jtype)), &
443	& nvarsprof(jtype), nextrprof(jtype), nitend-nit000+2, &
444	& rn_dobsini, rn_dobsend, llvar1, llvar2, &
445	& ln_ignmis, ln_s_at_t, .FALSE., &
446	& kdailyavtypes = nn_profdavtypes )
447
448	DO jvar = 1, nvarsprof(jtype)
449	CALL obs_prof_staend( profdata(jtype), jvar )
450	END DO
451
452	CALL obs_pre_prof( profdata(jtype), profdataqc(jtype), &
453	& llvar1, llvar2, &
454	& jpi, jpj, jpk, &
455	& zmask1, zglam1, zgphi1, zmask2, zglam2, zgphi2, &
456	& ln_nea, ln_bound_reject, &
457	& kdailyavtypes = nn_profdavtypes )
458
459	END DO
460
461	DEALLOCATE( ifilesprof, clproffiles )
462
463	ENDIF
464
465	IF ( nsurftypes > 0 ) THEN
466
467	ALLOCATE(surfdata(nsurftypes))
468	ALLOCATE(surfdataqc(nsurftypes))
469	ALLOCATE(nvarssurf(nsurftypes))
470	ALLOCATE(nextrsurf(nsurftypes))
471
472	DO jtype = 1, nsurftypes
473
474	nvarssurf(jtype) = 1
475	nextrsurf(jtype) = 0
476	llnightav(jtype) = .FALSE.
477	IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) nextrsurf(jtype) = 2
478	IF ( TRIM(cobstypessurf(jtype)) == 'sst' ) llnightav(jtype) = ln_sstnight
479
480	!Read in surface obs types
481	CALL obs_rea_surf( surfdata(jtype), ifilessurf(jtype), &
482	& clsurffiles(jtype,1:ifilessurf(jtype)), &
483	& nvarssurf(jtype), nextrsurf(jtype), nitend-nit000+2, &
484	& rn_dobsini, rn_dobsend, ln_ignmis, .FALSE., llnightav(jtype) )
485
486	CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea, ln_bound_reject )
487
488	IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN
489	CALL obs_rea_mdt( surfdataqc(jtype), n2dintsurf(jtype) )
490	IF ( ln_altbias ) &
491	& CALL obs_rea_altbias ( surfdataqc(jtype), n2dintsurf(jtype), cn_altbiasfile )
492	ENDIF
493
494	IF ( TRIM(cobstypessurf(jtype)) == 'sst' .AND. ln_sstbias ) THEN
495	jnumsstbias = 0
496	DO jfile = 1, jpmaxnfiles
497	IF ( TRIM(cn_sstbiasfiles(jfile)) /= '' ) &
498	& jnumsstbias = jnumsstbias + 1
499	END DO
500	IF ( jnumsstbias == 0 ) THEN
501	CALL ctl_stop("ln_sstbias set but no bias files to read in")
502	ENDIF
503
504	CALL obs_app_sstbias( surfdataqc(jtype), n2dintsurf(jtype), &
505	& jnumsstbias, cn_sstbiasfiles(1:jnumsstbias) )
506
507	ENDIF
508	END DO
509
510	DEALLOCATE( ifilessurf, clsurffiles )
511
512	ENDIF
513
514	CALL wrk_dealloc( jpi, jpj, zglam1 )
515	CALL wrk_dealloc( jpi, jpj, zglam2 )
516	CALL wrk_dealloc( jpi, jpj, zgphi1 )
517	CALL wrk_dealloc( jpi, jpj, zgphi2 )
518	CALL wrk_dealloc( jpi, jpj, jpk, zmask1 )
519	CALL wrk_dealloc( jpi, jpj, jpk, zmask2 )
520
521	END SUBROUTINE dia_obs_init
522
523
524	SUBROUTINE dia_obs( kstp )
525	!!----------------------------------------------------------------------
526	!! * ROUTINE dia_obs *
527	!!
528	!! ** Purpose : Call the observation operators on each time step
529	!!
530	!! ** Method : Call the observation operators on each time step to
531	!! compute the model equivalent of the following data:
532	!! - Profile data, currently T/S or U/V
533	!! - Surface data, currently SST, SLA or sea-ice concentration.
534	!!
535	!! ** Action :
536	!!
537	!! History :
538	!! ! 06-03 (K. Mogensen) Original code
539	!! ! 06-05 (K. Mogensen) Reformatted
540	!! ! 06-10 (A. Weaver) Cleaning
541	!! ! 07-03 (K. Mogensen) General handling of profiles
542	!! ! 07-04 (G. Smith) Generalized surface operators
543	!! ! 08-10 (M. Valdivieso) obs operator for velocity profiles
544	!! ! 14-08 (J. While) observation operator for profiles in
545	!! generalised vertical coordinates
546	!! ! 15-08 (M. Martin) Combined surface/profile routines.
547	!!----------------------------------------------------------------------
548	USE dom_oce, ONLY : gdept_n, gdept_1d ! Ocean space and time domain variables
549	USE phycst , ONLY : rday ! Physical constants
550	USE oce , ONLY : tsn, un, vn, sshn ! Ocean dynamics and tracers variables
551	USE phycst , ONLY : rday ! Physical constants
552	#if defined key_lim3
553	USE ice , ONLY : at_i ! LIM3 Ice model variables
554	#endif
555	#if defined key_cice
556	USE sbc_oce, ONLY : fr_i ! ice fraction
557	#endif
558
559	IMPLICIT NONE
560
561	!! * Arguments
562	INTEGER, INTENT(IN) :: kstp ! Current timestep
563	!! * Local declarations
564	INTEGER :: idaystp ! Number of timesteps per day
565	INTEGER :: jtype ! Data loop variable
566	INTEGER :: jvar ! Variable number
567	INTEGER :: ji, jj ! Loop counters
568	REAL(wp), POINTER, DIMENSION(:,:,:) :: &
569	& zprofvar1, & ! Model values for 1st variable in a prof ob
570	& zprofvar2 ! Model values for 2nd variable in a prof ob
571	REAL(wp), POINTER, DIMENSION(:,:,:) :: &
572	& zprofmask1, & ! Mask associated with zprofvar1
573	& zprofmask2 ! Mask associated with zprofvar2
574	REAL(wp), POINTER, DIMENSION(:,:) :: &
575	& zsurfvar, & ! Model values equivalent to surface ob.
576	& zsurfmask ! Mask associated with surface variable
577	REAL(wp), POINTER, DIMENSION(:,:) :: &
578	& zglam1, & ! Model longitudes for prof variable 1
579	& zglam2, & ! Model longitudes for prof variable 2
580	& zgphi1, & ! Model latitudes for prof variable 1
581	& zgphi2 ! Model latitudes for prof variable 2
582
583	!Allocate local work arrays
584	CALL wrk_alloc( jpi, jpj, jpk, zprofvar1 )
585	CALL wrk_alloc( jpi, jpj, jpk, zprofvar2 )
586	CALL wrk_alloc( jpi, jpj, jpk, zprofmask1 )
587	CALL wrk_alloc( jpi, jpj, jpk, zprofmask2 )
588	CALL wrk_alloc( jpi, jpj, zsurfvar )
589	CALL wrk_alloc( jpi, jpj, zsurfmask )
590	CALL wrk_alloc( jpi, jpj, zglam1 )
591	CALL wrk_alloc( jpi, jpj, zglam2 )
592	CALL wrk_alloc( jpi, jpj, zgphi1 )
593	CALL wrk_alloc( jpi, jpj, zgphi2 )
594	!-----------------------------------------------------------------------
595
596	IF(lwp) THEN
597	WRITE(numout,*)
598	WRITE(numout,*) 'dia_obs : Call the observation operators', kstp
599	WRITE(numout,*) '~~~~~~~'
600	ENDIF
601
602	idaystp = NINT( rday / rdt )
603
604	!-----------------------------------------------------------------------
605	! Call the profile and surface observation operators
606	!-----------------------------------------------------------------------
607
608	IF ( nproftypes > 0 ) THEN
609
610	DO jtype = 1, nproftypes
611
612	SELECT CASE ( TRIM(cobstypesprof(jtype)) )
613	CASE('prof')
614	zprofvar1(:,:,:) = tsn(:,:,:,jp_tem)
615	zprofvar2(:,:,:) = tsn(:,:,:,jp_sal)
616	zprofmask1(:,:,:) = tmask(:,:,:)
617	zprofmask2(:,:,:) = tmask(:,:,:)
618	zglam1(:,:) = glamt(:,:)
619	zglam2(:,:) = glamt(:,:)
620	zgphi1(:,:) = gphit(:,:)
621	zgphi2(:,:) = gphit(:,:)
622	CASE('vel')
623	zprofvar1(:,:,:) = un(:,:,:)
624	zprofvar2(:,:,:) = vn(:,:,:)
625	zprofmask1(:,:,:) = umask(:,:,:)
626	zprofmask2(:,:,:) = vmask(:,:,:)
627	zglam1(:,:) = glamu(:,:)
628	zglam2(:,:) = glamv(:,:)
629	zgphi1(:,:) = gphiu(:,:)
630	zgphi2(:,:) = gphiv(:,:)
631	CASE DEFAULT
632	CALL ctl_stop( 'Unknown profile observation type '//TRIM(cobstypesprof(jtype))//' in dia_obs' )
633	END SELECT
634
635	CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk, &
636	& nit000, idaystp, &
637	& zprofvar1, zprofvar2, &
638	& gdept_n(:,:,:), gdepw_n(:,:,:), &
639	& zprofmask1, zprofmask2, &
640	& zglam1, zglam2, zgphi1, zgphi2, &
641	& nn_1dint, nn_2dint, &
642	& kdailyavtypes = nn_profdavtypes )
643
644	END DO
645
646	ENDIF
647
648	IF ( nsurftypes > 0 ) THEN
649
650	DO jtype = 1, nsurftypes
651
652	!Defaults which might be changed
653	zsurfmask(:,:) = tmask(:,:,1)
654
655	SELECT CASE ( TRIM(cobstypessurf(jtype)) )
656	CASE('sst')
657	zsurfvar(:,:) = tsn(:,:,1,jp_tem)
658	CASE('sla')
659	zsurfvar(:,:) = sshn(:,:)
660	CASE('sss')
661	zsurfvar(:,:) = tsn(:,:,1,jp_sal)
662	CASE('sic')
663	IF ( kstp == 0 ) THEN
664	IF ( lwp .AND. surfdataqc(jtype)%nsstpmpp(1) > 0 ) THEN
665	CALL ctl_warn( 'Sea-ice not initialised on zeroth '// &
666	& 'time-step but some obs are valid then.' )
667	WRITE(numout,*)surfdataqc(jtype)%nsstpmpp(1), &
668	& ' sea-ice obs will be missed'
669	ENDIF
670	surfdataqc(jtype)%nsurfup = surfdataqc(jtype)%nsurfup + &
671	& surfdataqc(jtype)%nsstp(1)
672	CYCLE
673	ELSE
674	#if defined key_cice
675	zsurfvar(:,:) = fr_i(:,:)
676	#elif defined key_lim2 \|\| defined key_lim3
677	zsurfvar(:,:) = 1._wp - frld(:,:)
678	#else
679	CALL ctl_stop( ' Trying to run sea-ice observation operator', &
680	& ' but no sea-ice model appears to have been defined' )
681	#endif
682	ENDIF
683
684	END SELECT
685
686	CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj, &
687	& nit000, idaystp, zsurfvar, zsurfmask, &
688	& n2dintsurf(jtype), llnightav(jtype), &
689	& zavglamscl(jtype), zavgphiscl(jtype), &
690	& lfpindegs(jtype) )
691
692	END DO
693
694	ENDIF
695
696	CALL wrk_dealloc( jpi, jpj, jpk, zprofvar1 )
697	CALL wrk_dealloc( jpi, jpj, jpk, zprofvar2 )
698	CALL wrk_dealloc( jpi, jpj, jpk, zprofmask1 )
699	CALL wrk_dealloc( jpi, jpj, jpk, zprofmask2 )
700	CALL wrk_dealloc( jpi, jpj, zsurfvar )
701	CALL wrk_dealloc( jpi, jpj, zsurfmask )
702	CALL wrk_dealloc( jpi, jpj, zglam1 )
703	CALL wrk_dealloc( jpi, jpj, zglam2 )
704	CALL wrk_dealloc( jpi, jpj, zgphi1 )
705	CALL wrk_dealloc( jpi, jpj, zgphi2 )
706
707	END SUBROUTINE dia_obs
708
709	SUBROUTINE dia_obs_wri
710	!!----------------------------------------------------------------------
711	!! * ROUTINE dia_obs_wri *
712	!!
713	!! ** Purpose : Call observation diagnostic output routines
714	!!
715	!! ** Method : Call observation diagnostic output routines
716	!!
717	!! ** Action :
718	!!
719	!! History :
720	!! ! 06-03 (K. Mogensen) Original code
721	!! ! 06-05 (K. Mogensen) Reformatted
722	!! ! 06-10 (A. Weaver) Cleaning
723	!! ! 07-03 (K. Mogensen) General handling of profiles
724	!! ! 08-09 (M. Valdivieso) Velocity component (U,V) profiles
725	!! ! 15-08 (M. Martin) Combined writing for prof and surf types
726	!!----------------------------------------------------------------------
727	!! * Modules used
728	USE obs_rot_vel ! Rotation of velocities
729
730	IMPLICIT NONE
731
732	!! * Local declarations
733	INTEGER :: jtype ! Data set loop variable
734	INTEGER :: jo, jvar, jk
735	REAL(wp), DIMENSION(:), ALLOCATABLE :: &
736	& zu, &
737	& zv
738
739	!-----------------------------------------------------------------------
740	! Depending on switches call various observation output routines
741	!-----------------------------------------------------------------------
742
743	IF ( nproftypes > 0 ) THEN
744
745	DO jtype = 1, nproftypes
746
747	IF ( TRIM(cobstypesprof(jtype)) == 'vel' ) THEN
748
749	! For velocity data, rotate the model velocities to N/S, E/W
750	! using the compressed data structure.
751	ALLOCATE( &
752	& zu(profdataqc(jtype)%nvprot(1)), &
753	& zv(profdataqc(jtype)%nvprot(2)) &
754	& )
755
756	CALL obs_rotvel( profdataqc(jtype), nn_2dint, zu, zv )
757
758	DO jo = 1, profdataqc(jtype)%nprof
759	DO jvar = 1, 2
760	DO jk = profdataqc(jtype)%npvsta(jo,jvar), profdataqc(jtype)%npvend(jo,jvar)
761
762	IF ( jvar == 1 ) THEN
763	profdataqc(jtype)%var(jvar)%vmod(jk) = zu(jk)
764	ELSE
765	profdataqc(jtype)%var(jvar)%vmod(jk) = zv(jk)
766	ENDIF
767
768	END DO
769	END DO
770	END DO
771
772	DEALLOCATE( zu )
773	DEALLOCATE( zv )
774
775	END IF
776
777	CALL obs_prof_decompress( profdataqc(jtype), &
778	& profdata(jtype), .TRUE., numout )
779
780	CALL obs_wri_prof( profdata(jtype) )
781
782	END DO
783
784	ENDIF
785
786	IF ( nsurftypes > 0 ) THEN
787
788	DO jtype = 1, nsurftypes
789
790	CALL obs_surf_decompress( surfdataqc(jtype), &
791	& surfdata(jtype), .TRUE., numout )
792
793	CALL obs_wri_surf( surfdata(jtype) )
794
795	END DO
796
797	ENDIF
798
799	END SUBROUTINE dia_obs_wri
800
801	SUBROUTINE dia_obs_dealloc
802	IMPLICIT NONE
803	!!----------------------------------------------------------------------
804	!! * ROUTINE dia_obs_dealloc *
805	!!
806	!! ** Purpose : To deallocate data to enable the obs_oper online loop.
807	!! Specifically: dia_obs_init --> dia_obs --> dia_obs_wri
808	!!
809	!! ** Method : Clean up various arrays left behind by the obs_oper.
810	!!
811	!! ** Action :
812	!!
813	!!----------------------------------------------------------------------
814	! obs_grid deallocation
815	CALL obs_grid_deallocate
816
817	! diaobs deallocation
818	IF ( nproftypes > 0 ) &
819	& DEALLOCATE( cobstypesprof, profdata, profdataqc, nvarsprof, nextrprof )
820
821	IF ( nsurftypes > 0 ) &
822	& DEALLOCATE( cobstypessurf, surfdata, surfdataqc, nvarssurf, nextrsurf, &
823	& n2dintsurf, zavglamscl, zavgphiscl, lfpindegs, llnightav )
824
825	END SUBROUTINE dia_obs_dealloc
826
827	SUBROUTINE calc_date( kstp, ddobs )
828	!!----------------------------------------------------------------------
829	!! * ROUTINE calc_date *
830	!!
831	!! ** Purpose : Get date in double precision YYYYMMDD.HHMMSS format
832	!!
833	!! ** Method : Get date in double precision YYYYMMDD.HHMMSS format
834	!!
835	!! ** Action : Get date in double precision YYYYMMDD.HHMMSS format
836	!!
837	!! ** Action : Get initial date in double precision YYYYMMDD.HHMMSS format
838	!!
839	!! History :
840	!! ! 06-03 (K. Mogensen) Original code
841	!! ! 06-05 (K. Mogensen) Reformatted
842	!! ! 06-10 (A. Weaver) Cleaning
843	!! ! 06-10 (G. Smith) Calculates initial date the same as method for final date
844	!! ! 10-05 (D. Lea) Update to month length calculation for NEMO vn3.2
845	!! ! 2014-09 (D. Lea) New generic routine now deals with arbitrary initial time of day
846	!!----------------------------------------------------------------------
847	USE phycst, ONLY : & ! Physical constants
848	& rday
849	USE dom_oce, ONLY : & ! Ocean space and time domain variables
850	& rdt
851
852	IMPLICIT NONE
853
854	!! * Arguments
855	REAL(KIND=dp), INTENT(OUT) :: ddobs ! Date in YYYYMMDD.HHMMSS
856	INTEGER :: kstp
857
858	!! * Local declarations
859	INTEGER :: iyea ! date - (year, month, day, hour, minute)
860	INTEGER :: imon
861	INTEGER :: iday
862	INTEGER :: ihou
863	INTEGER :: imin
864	INTEGER :: imday ! Number of days in month.
865	REAL(wp) :: zdayfrc ! Fraction of day
866
867	INTEGER, DIMENSION(12) :: imonth_len !: length in days of the months of the current year
868
869	!!----------------------------------------------------------------------
870	!! Initial date initialization (year, month, day, hour, minute)
871	!!----------------------------------------------------------------------
872	iyea = ndate0 / 10000
873	imon = ( ndate0 - iyea * 10000 ) / 100
874	iday = ndate0 - iyea * 10000 - imon * 100
875	ihou = nn_time0 / 100
876	imin = ( nn_time0 - ihou * 100 )
877
878	!!----------------------------------------------------------------------
879	!! Compute number of days + number of hours + min since initial time
880	!!----------------------------------------------------------------------
881	zdayfrc = kstp * rdt / rday
882	zdayfrc = zdayfrc - aint(zdayfrc)
883	imin = imin + int( zdayfrc * 24 * 60 )
884	DO WHILE (imin >= 60)
885	imin=imin-60
886	ihou=ihou+1
887	END DO
888	DO WHILE (ihou >= 24)
889	ihou=ihou-24
890	iday=iday+1
891	END DO
892	iday = iday + kstp * rdt / rday
893
894	!-----------------------------------------------------------------------
895	! Convert number of days (iday) into a real date
896	!----------------------------------------------------------------------
897
898	CALL calc_month_len( iyea, imonth_len )
899
900	DO WHILE ( iday > imonth_len(imon) )
901	iday = iday - imonth_len(imon)
902	imon = imon + 1
903	IF ( imon > 12 ) THEN
904	imon = 1
905	iyea = iyea + 1
906	CALL calc_month_len( iyea, imonth_len ) ! update month lengths
907	ENDIF
908	END DO
909
910	!----------------------------------------------------------------------
911	! Convert it into YYYYMMDD.HHMMSS format.
912	!----------------------------------------------------------------------
913	ddobs = iyea * 10000_dp + imon * 100_dp + &
914	& iday + ihou * 0.01_dp + imin * 0.0001_dp
915
916	END SUBROUTINE calc_date
917
918	SUBROUTINE ini_date( ddobsini )
919	!!----------------------------------------------------------------------
920	!! * ROUTINE ini_date *
921	!!
922	!! ** Purpose : Get initial date in double precision YYYYMMDD.HHMMSS format
923	!!
924	!! ** Method :
925	!!
926	!! ** Action :
927	!!
928	!! History :
929	!! ! 06-03 (K. Mogensen) Original code
930	!! ! 06-05 (K. Mogensen) Reformatted
931	!! ! 06-10 (A. Weaver) Cleaning
932	!! ! 10-05 (D. Lea) Update to month length calculation for NEMO vn3.2
933	!! ! 2014-09 (D. Lea) Change to call generic routine calc_date
934	!!----------------------------------------------------------------------
935
936	IMPLICIT NONE
937
938	!! * Arguments
939	REAL(KIND=dp), INTENT(OUT) :: ddobsini ! Initial date in YYYYMMDD.HHMMSS
940
941	CALL calc_date( nit000 - 1, ddobsini )
942
943	END SUBROUTINE ini_date
944
945	SUBROUTINE fin_date( ddobsfin )
946	!!----------------------------------------------------------------------
947	!! * ROUTINE fin_date *
948	!!
949	!! ** Purpose : Get final date in double precision YYYYMMDD.HHMMSS format
950	!!
951	!! ** Method :
952	!!
953	!! ** Action :
954	!!
955	!! History :
956	!! ! 06-03 (K. Mogensen) Original code
957	!! ! 06-05 (K. Mogensen) Reformatted
958	!! ! 06-10 (A. Weaver) Cleaning
959	!! ! 10-05 (D. Lea) Update to month length calculation for NEMO vn3.2
960	!! ! 2014-09 (D. Lea) Change to call generic routine calc_date
961	!!----------------------------------------------------------------------
962
963	IMPLICIT NONE
964
965	!! * Arguments
966	REAL(dp), INTENT(OUT) :: ddobsfin ! Final date in YYYYMMDD.HHMMSS
967
968	CALL calc_date( nitend, ddobsfin )
969
970	END SUBROUTINE fin_date
971
972	SUBROUTINE obs_settypefiles( ntypes, jpmaxnfiles, jtype, ctypein, &
973	& cfilestype, ifiles, cobstypes, cfiles )
974
975	INTEGER, INTENT(IN) :: ntypes ! Total number of obs types
976	INTEGER, INTENT(IN) :: jpmaxnfiles ! Maximum number of files allowed for each type
977	INTEGER, INTENT(IN) :: jtype ! Index of the current type of obs
978	INTEGER, DIMENSION(ntypes), INTENT(INOUT) :: &
979	& ifiles ! Out appended number of files for this type
980
981	CHARACTER(len=6), INTENT(IN) :: ctypein
982	CHARACTER(len=128), DIMENSION(jpmaxnfiles), INTENT(IN) :: &
983	& cfilestype ! In list of files for this obs type
984	CHARACTER(len=6), DIMENSION(ntypes), INTENT(INOUT) :: &
985	& cobstypes ! Out appended list of obs types
986	CHARACTER(len=128), DIMENSION(ntypes, jpmaxnfiles), INTENT(INOUT) :: &
987	& cfiles ! Out appended list of files for all types
988
989	!Local variables
990	INTEGER :: jfile
991
992	cfiles(jtype,:) = cfilestype(:)
993	cobstypes(jtype) = ctypein
994	ifiles(jtype) = 0
995	DO jfile = 1, jpmaxnfiles
996	IF ( trim(cfiles(jtype,jfile)) /= '' ) &
997	ifiles(jtype) = ifiles(jtype) + 1
998	END DO
999
1000	IF ( ifiles(jtype) == 0 ) THEN
1001	CALL ctl_stop( 'Logical for observation type '//TRIM(ctypein)// &
1002	& ' set to true but no files available to read' )
1003	ENDIF
1004
1005	IF(lwp) THEN
1006	WRITE(numout,*) ' '//cobstypes(jtype)//' input observation file names:'
1007	DO jfile = 1, ifiles(jtype)
1008	WRITE(numout,*) ' '//TRIM(cfiles(jtype,jfile))
1009	END DO
1010	ENDIF
1011
1012	END SUBROUTINE obs_settypefiles
1013
1014	SUBROUTINE obs_setinterpopts( ntypes, jtype, ctypein, &
1015	& n2dint_default, n2dint_type, &
1016	& zavglamscl_type, zavgphiscl_type, &
1017	& lfp_indegs_type, lavnight_type, &
1018	& n2dint, zavglamscl, zavgphiscl, &
1019	& lfpindegs, lavnight )
1020
1021	INTEGER, INTENT(IN) :: ntypes ! Total number of obs types
1022	INTEGER, INTENT(IN) :: jtype ! Index of the current type of obs
1023	INTEGER, INTENT(IN) :: n2dint_default ! Default option for interpolation type
1024	INTEGER, INTENT(IN) :: n2dint_type ! Option for interpolation type
1025	REAL(wp), INTENT(IN) :: &
1026	& zavglamscl_type, & !E/W diameter of obs footprint for this type
1027	& zavgphiscl_type !N/S diameter of obs footprint for this type
1028	LOGICAL, INTENT(IN) :: lfp_indegs_type !T=> footprint in degrees, F=> in metres
1029	LOGICAL, INTENT(IN) :: lavnight_type !T=> obs represent night time average
1030	CHARACTER(len=6), INTENT(IN) :: ctypein
1031
1032	INTEGER, DIMENSION(ntypes), INTENT(INOUT) :: &
1033	& n2dint
1034	REAL(wp), DIMENSION(ntypes), INTENT(INOUT) :: &
1035	& zavglamscl, zavgphiscl
1036	LOGICAL, DIMENSION(ntypes), INTENT(INOUT) :: &
1037	& lfpindegs, lavnight
1038
1039	lavnight(jtype) = lavnight_type
1040
1041	IF ( (n2dint_type >= 1) .AND. (n2dint_type <= 6) ) THEN
1042	n2dint(jtype) = n2dint_type
1043	ELSE
1044	n2dint(jtype) = n2dint_default
1045	ENDIF
1046
1047	! For averaging observation footprints set options for size of footprint
1048	IF ( (n2dint(jtype) > 4) .AND. (n2dint(jtype) <= 6) ) THEN
1049	IF ( zavglamscl_type > 0._wp ) THEN
1050	zavglamscl(jtype) = zavglamscl_type
1051	ELSE
1052	CALL ctl_stop( 'Incorrect value set for averaging footprint '// &
1053	'scale (zavglamscl) for observation type '//TRIM(ctypein) )
1054	ENDIF
1055
1056	IF ( zavgphiscl_type > 0._wp ) THEN
1057	zavgphiscl(jtype) = zavgphiscl_type
1058	ELSE
1059	CALL ctl_stop( 'Incorrect value set for averaging footprint '// &
1060	'scale (zavgphiscl) for observation type '//TRIM(ctypein) )
1061	ENDIF
1062
1063	lfpindegs(jtype) = lfp_indegs_type
1064
1065	ENDIF
1066
1067	! Write out info
1068	IF(lwp) THEN
1069	IF ( n2dint(jtype) <= 4 ) THEN
1070	WRITE(numout,*) ' '//TRIM(ctypein)// &
1071	& ' model counterparts will be interpolated horizontally'
1072	ELSE IF ( n2dint(jtype) <= 6 ) THEN
1073	WRITE(numout,*) ' '//TRIM(ctypein)// &
1074	& ' model counterparts will be averaged horizontally'
1075	WRITE(numout,*) ' '//' with E/W scale: ',zavglamscl(jtype)
1076	WRITE(numout,*) ' '//' with N/S scale: ',zavgphiscl(jtype)
1077	IF ( lfpindegs(jtype) ) THEN
1078	WRITE(numout,*) ' '//' (in degrees)'
1079	ELSE
1080	WRITE(numout,*) ' '//' (in metres)'
1081	ENDIF
1082	ENDIF
1083	ENDIF
1084
1085	END SUBROUTINE obs_setinterpopts
1086
1087	END MODULE diaobs

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

Download in other formats: