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

Context Navigation

source: branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 10120

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

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

Download in other formats: