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

Context Navigation

source: branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 5998

Last change on this file since 5998 was 5998, checked in by timgraham, 8 years ago
Merge in changes from OBS simplification branch (branches/2015/dev_r5072_UKMO2_OBS_simplification)
Property svn:keywords set to `Id`
File size: 42.5 KB

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

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

Download in other formats: