Context Navigation

diaobs.F90 @ 8738

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

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

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

Context Navigation

source: branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 8738

Download in other formats: