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

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

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

Context Navigation

source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 7 years ago
All wrk_alloc removed
Property svn:keywords set to `Id`
File size: 39.2 KB

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

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

Download in other formats: