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 @ 6042

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

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

Download in other formats: