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

Context Navigation

source: branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 5682

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

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

Download in other formats: