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 NEMO/trunk/src/OCE/OBS – NEMO

Context Navigation

source: NEMO/trunk/src/OCE/OBS/diaobs.F90 @ 14996

Last change on this file since 14996 was 14982, checked in by hadcv, 3 years ago
#2665: Various fixes for code enabled with key_qco/key_linssh
Property svn:keywords set to `Id`
File size: 52.6 KB

Rev	Line
[2128]	1	MODULE diaobs
	2	!!======================================================================
	3	!! * MODULE diaobs *
	4	!! Observation diagnostics: Computation of the misfit between data and
	5	!! their model equivalent
	6	!!======================================================================
[9168]	7	!! History : 1.0 ! 2006-03 (K. Mogensen) Original code
	8	!! - ! 2006-05 (K. Mogensen, A. Weaver) Reformatted
	9	!! - ! 2006-10 (A. Weaver) Cleaning and add controls
	10	!! - ! 2007-03 (K. Mogensen) General handling of profiles
	11	!! - ! 2007-04 (G. Smith) Generalized surface operators
	12	!! 2.0 ! 2008-10 (M. Valdivieso) obs operator for velocity profiles
	13	!! 3.4 ! 2014-08 (J. While) observation operator for profiles in all vertical coordinates
	14	!! - ! Incorporated SST bias correction
	15	!! 3.6 ! 2015-02 (M. Martin) Simplification of namelist and code
	16	!! - ! 2015-08 (M. Martin) Combined surface/profile routines.
	17	!! 4.0 ! 2017-11 (G. Madec) style only
	18	!!----------------------------------------------------------------------
[2128]	19
	20	!!----------------------------------------------------------------------
[9168]	21	!! dia_obs_init : Reading and prepare observations
	22	!! dia_obs : Compute model equivalent to observations
	23	!! dia_obs_wri : Write observational diagnostics
	24	!! calc_date : Compute the date of timestep in YYYYMMDD.HHMMSS format
	25	!! ini_date : Compute the initial date YYYYMMDD.HHMMSS
	26	!! fin_date : Compute the final date YYYYMMDD.HHMMSS
[2128]	27	!!----------------------------------------------------------------------
[9168]	28	USE par_kind ! Precision variables
	29	USE in_out_manager ! I/O manager
	30	USE par_oce ! ocean parameter
	31	USE dom_oce ! Ocean space and time domain variables
	32	USE sbc_oce ! Sea-ice fraction
	33	!
	34	USE obs_read_prof ! Reading and allocation of profile obs
	35	USE obs_read_surf ! Reading and allocation of surface obs
	36	USE obs_sstbias ! Bias correction routine for SST
	37	USE obs_readmdt ! Reading and allocation of MDT for SLA.
	38	USE obs_prep ! Preparation of obs. (grid search etc).
	39	USE obs_oper ! Observation operators
	40	USE obs_write ! Writing of observation related diagnostics
	41	USE obs_grid ! Grid searching
	42	USE obs_read_altbias ! Bias treatment for altimeter
	43	USE obs_profiles_def ! Profile data definitions
	44	USE obs_surf_def ! Surface data definitions
	45	USE obs_types ! Definitions for observation types
	46	!
	47	USE mpp_map ! MPP mapping
	48	USE lib_mpp ! For ctl_warn/stop
[2128]	49
	50	IMPLICIT NONE
	51	PRIVATE
	52
[9019]	53	PUBLIC dia_obs_init ! Initialize and read observations
	54	PUBLIC dia_obs ! Compute model equivalent to observations
	55	PUBLIC dia_obs_wri ! Write model equivalent to observations
	56	PUBLIC dia_obs_dealloc ! Deallocate dia_obs data
	57	PUBLIC calc_date ! Compute the date of a timestep
	58
[14056]	59	LOGICAL, PUBLIC :: ln_diaobs !: Logical switch for the obs operator
	60	LOGICAL :: ln_sstnight ! Logical switch for night mean SST obs
	61	LOGICAL :: ln_default_fp_indegs ! T=> Default obs footprint size specified in degrees, F=> in metres
	62	LOGICAL :: ln_sla_fp_indegs ! T=> SLA obs footprint size specified in degrees, F=> in metres
	63	LOGICAL :: ln_sst_fp_indegs ! T=> SST obs footprint size specified in degrees, F=> in metres
	64	LOGICAL :: ln_sss_fp_indegs ! T=> SSS obs footprint size specified in degrees, F=> in metres
	65	LOGICAL :: ln_sic_fp_indegs ! T=> sea-ice obs footprint size specified in degrees, F=> in metres
[9023]	66
[14056]	67	REAL(wp) :: rn_default_avglamscl ! E/W diameter of SLA observation footprint (metres)
	68	REAL(wp) :: rn_default_avgphiscl ! N/S diameter of SLA observation footprint (metre
	69	REAL(wp) :: rn_sla_avglamscl ! E/W diameter of SLA observation footprint (metres)
	70	REAL(wp) :: rn_sla_avgphiscl ! N/S diameter of SLA observation footprint (metres)
	71	REAL(wp) :: rn_sst_avglamscl ! E/W diameter of SST observation footprint (metres)
	72	REAL(wp) :: rn_sst_avgphiscl ! N/S diameter of SST observation footprint (metres)
	73	REAL(wp) :: rn_sss_avglamscl ! E/W diameter of SSS observation footprint (metres)
	74	REAL(wp) :: rn_sss_avgphiscl ! N/S diameter of SSS observation footprint (metres)
	75	REAL(wp) :: rn_sic_avglamscl ! E/W diameter of sea-ice observation footprint (metres)
	76	REAL(wp) :: rn_sic_avgphiscl ! N/S diameter of sea-ice observation footprint (metres)
[9023]	77
[14056]	78	INTEGER :: nn_1dint ! Vertical interpolation method
	79	INTEGER :: nn_2dint_default ! Default horizontal interpolation method
	80	INTEGER :: nn_2dint_sla ! SLA horizontal interpolation method
	81	INTEGER :: nn_2dint_sst ! SST horizontal interpolation method
	82	INTEGER :: nn_2dint_sss ! SSS horizontal interpolation method
	83	INTEGER :: nn_2dint_sic ! Seaice horizontal interpolation method
[9168]	84	INTEGER, DIMENSION(imaxavtypes) :: nn_profdavtypes ! Profile data types representing a daily average
	85	INTEGER :: nproftypes ! Number of profile obs types
	86	INTEGER :: nsurftypes ! Number of surface obs types
	87	INTEGER , DIMENSION(:), ALLOCATABLE :: nvarsprof, nvarssurf ! Number of profile & surface variables
	88	INTEGER , DIMENSION(:), ALLOCATABLE :: nextrprof, nextrsurf ! Number of profile & surface extra variables
	89	INTEGER , DIMENSION(:), ALLOCATABLE :: n2dintsurf ! Interpolation option for surface variables
	90	REAL(wp), DIMENSION(:), ALLOCATABLE :: zavglamscl, zavgphiscl ! E/W & N/S diameter of averaging footprint for surface variables
	91	LOGICAL , DIMENSION(:), ALLOCATABLE :: lfpindegs ! T=> surface obs footprint size specified in degrees, F=> in metres
	92	LOGICAL , DIMENSION(:), ALLOCATABLE :: llnightav ! Logical for calculating night-time averages
[2128]	93
[9168]	94	TYPE(obs_surf), PUBLIC, POINTER, DIMENSION(:) :: surfdata !: Initial surface data
	95	TYPE(obs_surf), PUBLIC, POINTER, DIMENSION(:) :: surfdataqc !: Surface data after quality control
	96	TYPE(obs_prof), PUBLIC, POINTER, DIMENSION(:) :: profdata !: Initial profile data
	97	TYPE(obs_prof), PUBLIC, POINTER, DIMENSION(:) :: profdataqc !: Profile data after quality control
[9023]	98
[14056]	99	CHARACTER(len=8), PUBLIC, DIMENSION(:), ALLOCATABLE :: cobstypesprof, cobstypessurf !: Profile & surface obs types
[2128]	100
[14982]	101	# include "domzgr_substitute.h90"
[2287]	102	!!----------------------------------------------------------------------
[9598]	103	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[2287]	104	!! $Id$
[10068]	105	!! Software governed by the CeCILL license (see ./LICENSE)
[2287]	106	!!----------------------------------------------------------------------
[2128]	107	CONTAINS
	108
[12377]	109	SUBROUTINE dia_obs_init( Kmm )
[2128]	110	!!----------------------------------------------------------------------
	111	!! * ROUTINE dia_obs_init *
	112	!!
	113	!! ** Purpose : Initialize and read observations
	114	!!
	115	!! ** Method : Read the namelist and call reading routines
	116	!!
	117	!! ** Action : Read the namelist and call reading routines
	118	!!
	119	!!----------------------------------------------------------------------
[12377]	120	INTEGER, INTENT(in) :: Kmm ! ocean time level indices
	121	INTEGER, PARAMETER :: jpmaxnfiles = 1000 ! Maximum number of files for each obs type
[9019]	122	INTEGER, DIMENSION(:), ALLOCATABLE :: ifilesprof, ifilessurf ! Number of profile & surface files
[6140]	123	INTEGER :: ios ! Local integer output status for namelist read
	124	INTEGER :: jtype ! Counter for obs types
	125	INTEGER :: jvar ! Counter for variables
	126	INTEGER :: jfile ! Counter for files
[14056]	127	INTEGER :: jnumsstbias ! Number of SST bias files to read and apply
	128	INTEGER :: n2dint_type ! Local version of nn_2dint*
[9168]	129	!
[6140]	130	CHARACTER(len=128), DIMENSION(jpmaxnfiles) :: &
	131	& cn_profbfiles, & ! T/S profile input filenames
	132	& cn_sstfbfiles, & ! Sea surface temperature input filenames
[9023]	133	& cn_sssfbfiles, & ! Sea surface salinity input filenames
[6140]	134	& cn_slafbfiles, & ! Sea level anomaly input filenames
	135	& cn_sicfbfiles, & ! Seaice concentration input filenames
	136	& cn_velfbfiles, & ! Velocity profile input filenames
[14056]	137	& cn_sstbiasfiles ! SST bias input filenames
[6140]	138	CHARACTER(LEN=128) :: &
	139	& cn_altbiasfile ! Altimeter bias input filename
	140	CHARACTER(len=128), DIMENSION(:,:), ALLOCATABLE :: &
	141	& clproffiles, & ! Profile filenames
	142	& clsurffiles ! Surface filenames
[14056]	143	CHARACTER(len=8), DIMENSION(:), ALLOCATABLE :: &
	144	& clvars ! Expected variable names
[9168]	145	!
[6140]	146	LOGICAL :: ln_t3d ! Logical switch for temperature profiles
	147	LOGICAL :: ln_s3d ! Logical switch for salinity profiles
	148	LOGICAL :: ln_sla ! Logical switch for sea level anomalies
	149	LOGICAL :: ln_sst ! Logical switch for sea surface temperature
[9023]	150	LOGICAL :: ln_sss ! Logical switch for sea surface salinity
[6140]	151	LOGICAL :: ln_sic ! Logical switch for sea ice concentration
	152	LOGICAL :: ln_vel3d ! Logical switch for velocity (u,v) obs
	153	LOGICAL :: ln_nea ! Logical switch to remove obs near land
	154	LOGICAL :: ln_altbias ! Logical switch for altimeter bias
[9019]	155	LOGICAL :: ln_sstbias ! Logical switch for bias corection of SST
[6140]	156	LOGICAL :: ln_ignmis ! Logical switch for ignoring missing files
	157	LOGICAL :: ln_s_at_t ! Logical switch to compute model S at T obs
[9023]	158	LOGICAL :: ln_bound_reject ! Logical to remove obs near boundaries in LAMs.
[14056]	159	LOGICAL :: ltype_fp_indegs ! Local version of ln_*_fp_indegs
	160	LOGICAL :: ltype_night ! Local version of ln_sstnight (false for other variables)
	161	LOGICAL, DIMENSION(:), ALLOCATABLE :: llvar ! Logical for profile variable read
[6140]	162	LOGICAL, DIMENSION(jpmaxnfiles) :: lmask ! Used for finding number of sstbias files
[9168]	163	!
[14056]	164	REAL(dp) :: rn_dobsini ! Obs window start date YYYYMMDD.HHMMSS
	165	REAL(dp) :: rn_dobsend ! Obs window end date YYYYMMDD.HHMMSS
	166	REAL(wp) :: ztype_avglamscl ! Local version of rn_*_avglamscl
	167	REAL(wp) :: ztype_avgphiscl ! Local version of rn_*_avgphiscl
	168	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zglam ! Model longitudes for profile variables
	169	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zgphi ! Model latitudes for profile variables
	170	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: zmask ! Model land/sea mask associated with variables
[9168]	171	!!
[6140]	172	NAMELIST/namobs/ln_diaobs, ln_t3d, ln_s3d, ln_sla, &
[9023]	173	& ln_sst, ln_sic, ln_sss, ln_vel3d, &
	174	& ln_altbias, ln_sstbias, ln_nea, &
	175	& ln_grid_global, ln_grid_search_lookup, &
	176	& ln_ignmis, ln_s_at_t, ln_bound_reject, &
[14056]	177	& ln_sstnight, ln_default_fp_indegs, &
[9023]	178	& ln_sla_fp_indegs, ln_sst_fp_indegs, &
	179	& ln_sss_fp_indegs, ln_sic_fp_indegs, &
[6140]	180	& cn_profbfiles, cn_slafbfiles, &
	181	& cn_sstfbfiles, cn_sicfbfiles, &
[9023]	182	& cn_velfbfiles, cn_sssfbfiles, &
	183	& cn_sstbiasfiles, cn_altbiasfile, &
[6140]	184	& cn_gridsearchfile, rn_gridsearchres, &
[9023]	185	& rn_dobsini, rn_dobsend, &
[14056]	186	& rn_default_avglamscl, rn_default_avgphiscl, &
[9023]	187	& rn_sla_avglamscl, rn_sla_avgphiscl, &
	188	& rn_sst_avglamscl, rn_sst_avgphiscl, &
	189	& rn_sss_avglamscl, rn_sss_avgphiscl, &
	190	& rn_sic_avglamscl, rn_sic_avgphiscl, &
[14056]	191	& nn_1dint, nn_2dint_default, &
[9023]	192	& nn_2dint_sla, nn_2dint_sst, &
	193	& nn_2dint_sss, nn_2dint_sic, &
[6140]	194	& nn_msshc, rn_mdtcorr, rn_mdtcutoff, &
[9023]	195	& nn_profdavtypes
[9168]	196	!-----------------------------------------------------------------------
[2128]	197
	198	!-----------------------------------------------------------------------
	199	! Read namelist parameters
	200	!-----------------------------------------------------------------------
[6140]	201	! Some namelist arrays need initialising
[9168]	202	cn_profbfiles (:) = ''
	203	cn_slafbfiles (:) = ''
	204	cn_sstfbfiles (:) = ''
	205	cn_sicfbfiles (:) = ''
	206	cn_velfbfiles (:) = ''
	207	cn_sssfbfiles (:) = ''
[9023]	208	cn_sstbiasfiles(:) = ''
[6140]	209	nn_profdavtypes(:) = -1
[2128]	210
[6140]	211	CALL ini_date( rn_dobsini )
	212	CALL fin_date( rn_dobsend )
	213
	214	! Read namelist namobs : control observation diagnostics
[4147]	215	READ ( numnam_ref, namobs, IOSTAT = ios, ERR = 901)
[11536]	216	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in reference namelist' )
[4147]	217	READ ( numnam_cfg, namobs, IOSTAT = ios, ERR = 902 )
[11536]	218	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namobs in configuration namelist' )
[4624]	219	IF(lwm) WRITE ( numond, namobs )
[4147]	220
[9168]	221	IF( .NOT.ln_diaobs ) THEN
	222	IF(lwp) WRITE(numout,*)
	223	IF(lwp) WRITE(numout,*) 'dia_obs_init : NO Observation diagnostic used'
	224	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
[2128]	225	RETURN
	226	ENDIF
[9023]	227
	228	IF(lwp) THEN
	229	WRITE(numout,*)
	230	WRITE(numout,*) 'dia_obs_init : Observation diagnostic initialization'
	231	WRITE(numout,*) '~~~~~~~~~~~~'
[9168]	232	WRITE(numout,*) ' Namelist namobs : set observation diagnostic parameters'
	233	WRITE(numout,*) ' Logical switch for T profile observations ln_t3d = ', ln_t3d
	234	WRITE(numout,*) ' Logical switch for S profile observations ln_s3d = ', ln_s3d
	235	WRITE(numout,*) ' Logical switch for SLA observations ln_sla = ', ln_sla
	236	WRITE(numout,*) ' Logical switch for SST observations ln_sst = ', ln_sst
	237	WRITE(numout,*) ' Logical switch for Sea Ice observations ln_sic = ', ln_sic
	238	WRITE(numout,*) ' Logical switch for velocity observations ln_vel3d = ', ln_vel3d
	239	WRITE(numout,*) ' Logical switch for SSS observations ln_sss = ', ln_sss
	240	WRITE(numout,*) ' Global distribution of observations ln_grid_global = ', ln_grid_global
	241	WRITE(numout,*) ' Logical switch for obs grid search lookup ln_grid_search_lookup = ', ln_grid_search_lookup
[9023]	242	IF (ln_grid_search_lookup) &
[9168]	243	WRITE(numout,*) ' Grid search lookup file header cn_gridsearchfile = ', cn_gridsearchfile
	244	WRITE(numout,*) ' Initial date in window YYYYMMDD.HHMMSS rn_dobsini = ', rn_dobsini
	245	WRITE(numout,*) ' Final date in window YYYYMMDD.HHMMSS rn_dobsend = ', rn_dobsend
	246	WRITE(numout,*) ' Type of vertical interpolation method nn_1dint = ', nn_1dint
[14056]	247	WRITE(numout,*) ' Default horizontal interpolation method nn_2dint_default = ', nn_2dint_default
	248	WRITE(numout,*) ' Type of horizontal interpolation method for SLA nn_2dint_sla = ', nn_2dint_sla
	249	WRITE(numout,*) ' Type of horizontal interpolation method for SST nn_2dint_sst = ', nn_2dint_sst
	250	WRITE(numout,*) ' Type of horizontal interpolation method for SSS nn_2dint_sss = ', nn_2dint_sss
	251	WRITE(numout,*) ' Type of horizontal interpolation method for SIC nn_2dint_sic = ', nn_2dint_sic
	252	WRITE(numout,*) ' Default E/W diameter of obs footprint rn_default_avglamscl = ', rn_default_avglamscl
	253	WRITE(numout,*) ' Default N/S diameter of obs footprint rn_default_avgphiscl = ', rn_default_avgphiscl
	254	WRITE(numout,*) ' Default obs footprint in deg [T] or m [F] ln_default_fp_indegs = ', ln_default_fp_indegs
	255	WRITE(numout,*) ' SLA E/W diameter of obs footprint rn_sla_avglamscl = ', rn_sla_avglamscl
	256	WRITE(numout,*) ' SLA N/S diameter of obs footprint rn_sla_avgphiscl = ', rn_sla_avgphiscl
	257	WRITE(numout,*) ' SLA obs footprint in deg [T] or m [F] ln_sla_fp_indegs = ', ln_sla_fp_indegs
	258	WRITE(numout,*) ' SST E/W diameter of obs footprint rn_sst_avglamscl = ', rn_sst_avglamscl
	259	WRITE(numout,*) ' SST N/S diameter of obs footprint rn_sst_avgphiscl = ', rn_sst_avgphiscl
	260	WRITE(numout,*) ' SST obs footprint in deg [T] or m [F] ln_sst_fp_indegs = ', ln_sst_fp_indegs
	261	WRITE(numout,*) ' SIC E/W diameter of obs footprint rn_sic_avglamscl = ', rn_sic_avglamscl
	262	WRITE(numout,*) ' SIC N/S diameter of obs footprint rn_sic_avgphiscl = ', rn_sic_avgphiscl
	263	WRITE(numout,*) ' SIC obs footprint in deg [T] or m [F] ln_sic_fp_indegs = ', ln_sic_fp_indegs
[9168]	264	WRITE(numout,*) ' Rejection of observations near land switch ln_nea = ', ln_nea
	265	WRITE(numout,*) ' Rejection of obs near open bdys ln_bound_reject = ', ln_bound_reject
	266	WRITE(numout,*) ' MSSH correction scheme nn_msshc = ', nn_msshc
	267	WRITE(numout,*) ' MDT correction rn_mdtcorr = ', rn_mdtcorr
	268	WRITE(numout,*) ' MDT cutoff for computed correction rn_mdtcutoff = ', rn_mdtcutoff
	269	WRITE(numout,*) ' Logical switch for alt bias ln_altbias = ', ln_altbias
	270	WRITE(numout,*) ' Logical switch for sst bias ln_sstbias = ', ln_sstbias
	271	WRITE(numout,*) ' Logical switch for ignoring missing files ln_ignmis = ', ln_ignmis
	272	WRITE(numout,*) ' Daily average types nn_profdavtypes = ', nn_profdavtypes
	273	WRITE(numout,*) ' Logical switch for night-time SST obs ln_sstnight = ', ln_sstnight
[9023]	274	ENDIF
[2128]	275	!-----------------------------------------------------------------------
[6140]	276	! Set up list of observation types to be used
	277	! and the files associated with each type
[2128]	278	!-----------------------------------------------------------------------
	279
[6140]	280	nproftypes = COUNT( (/ln_t3d .OR. ln_s3d, ln_vel3d /) )
[9023]	281	nsurftypes = COUNT( (/ln_sla, ln_sst, ln_sic, ln_sss /) )
[2128]	282
[9168]	283	IF( ln_sstbias ) THEN
[6140]	284	lmask(:) = .FALSE.
[9168]	285	WHERE( cn_sstbiasfiles(:) /= '' ) lmask(:) = .TRUE.
[6140]	286	jnumsstbias = COUNT(lmask)
	287	lmask(:) = .FALSE.
	288	ENDIF
[2128]	289
[9168]	290	IF( nproftypes == 0 .AND. nsurftypes == 0 ) THEN
	291	CALL ctl_warn( 'dia_obs_init: ln_diaobs is set to true, but all obs operator logical flags', &
	292	& ' (ln_t3d, ln_s3d, ln_sla, ln_sst, ln_sic, ln_vel3d)', &
	293	& ' are set to .FALSE. so turning off calls to dia_obs' )
[6140]	294	ln_diaobs = .FALSE.
	295	RETURN
	296	ENDIF
[2128]	297
[9168]	298	IF( nproftypes > 0 ) THEN
	299	!
	300	ALLOCATE( cobstypesprof(nproftypes) )
	301	ALLOCATE( ifilesprof (nproftypes) )
	302	ALLOCATE( clproffiles (nproftypes,jpmaxnfiles) )
	303	!
[6140]	304	jtype = 0
[9168]	305	IF( ln_t3d .OR. ln_s3d ) THEN
[6140]	306	jtype = jtype + 1
[14056]	307	cobstypesprof(jtype) = 'prof'
	308	clproffiles(jtype,:) = cn_profbfiles
[2128]	309	ENDIF
[9168]	310	IF( ln_vel3d ) THEN
[6140]	311	jtype = jtype + 1
[14056]	312	cobstypesprof(jtype) = 'vel'
	313	clproffiles(jtype,:) = cn_velfbfiles
[6140]	314	ENDIF
[9168]	315	!
[14056]	316	CALL obs_settypefiles( nproftypes, jpmaxnfiles, ifilesprof, cobstypesprof, clproffiles )
	317	!
[6140]	318	ENDIF
[2128]	319
[9168]	320	IF( nsurftypes > 0 ) THEN
	321	!
	322	ALLOCATE( cobstypessurf(nsurftypes) )
	323	ALLOCATE( ifilessurf (nsurftypes) )
	324	ALLOCATE( clsurffiles (nsurftypes,jpmaxnfiles) )
	325	ALLOCATE( n2dintsurf (nsurftypes) )
	326	ALLOCATE( zavglamscl (nsurftypes) )
	327	ALLOCATE( zavgphiscl (nsurftypes) )
	328	ALLOCATE( lfpindegs (nsurftypes) )
	329	ALLOCATE( llnightav (nsurftypes) )
	330	!
[6140]	331	jtype = 0
[9168]	332	IF( ln_sla ) THEN
[6140]	333	jtype = jtype + 1
[14056]	334	cobstypessurf(jtype) = 'sla'
	335	clsurffiles(jtype,:) = cn_slafbfiles
[2128]	336	ENDIF
[9168]	337	IF( ln_sst ) THEN
[6140]	338	jtype = jtype + 1
[14056]	339	cobstypessurf(jtype) = 'sst'
	340	clsurffiles(jtype,:) = cn_sstfbfiles
[2128]	341	ENDIF
[9570]	342	#if defined key_si3 \|\| defined key_cice
[9168]	343	IF( ln_sic ) THEN
[6140]	344	jtype = jtype + 1
[14056]	345	cobstypessurf(jtype) = 'sic'
	346	clsurffiles(jtype,:) = cn_sicfbfiles
[6140]	347	ENDIF
	348	#endif
[9168]	349	IF( ln_sss ) THEN
[9023]	350	jtype = jtype + 1
[14056]	351	cobstypessurf(jtype) = 'sss'
	352	clsurffiles(jtype,:) = cn_sssfbfiles
[9023]	353	ENDIF
[9168]	354	!
[14056]	355	CALL obs_settypefiles( nsurftypes, jpmaxnfiles, ifilessurf, cobstypessurf, clsurffiles )
	356
	357	DO jtype = 1, nsurftypes
	358
	359	IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN
	360	IF ( nn_2dint_sla == -1 ) THEN
	361	n2dint_type = nn_2dint_default
	362	ELSE
	363	n2dint_type = nn_2dint_sla
	364	ENDIF
	365	ztype_avglamscl = rn_sla_avglamscl
	366	ztype_avgphiscl = rn_sla_avgphiscl
	367	ltype_fp_indegs = ln_sla_fp_indegs
	368	ltype_night = .FALSE.
	369	ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sst' ) THEN
	370	IF ( nn_2dint_sst == -1 ) THEN
	371	n2dint_type = nn_2dint_default
	372	ELSE
	373	n2dint_type = nn_2dint_sst
	374	ENDIF
	375	ztype_avglamscl = rn_sst_avglamscl
	376	ztype_avgphiscl = rn_sst_avgphiscl
	377	ltype_fp_indegs = ln_sst_fp_indegs
	378	ltype_night = ln_sstnight
	379	ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sic' ) THEN
	380	IF ( nn_2dint_sic == -1 ) THEN
	381	n2dint_type = nn_2dint_default
	382	ELSE
	383	n2dint_type = nn_2dint_sic
	384	ENDIF
	385	ztype_avglamscl = rn_sic_avglamscl
	386	ztype_avgphiscl = rn_sic_avgphiscl
	387	ltype_fp_indegs = ln_sic_fp_indegs
	388	ltype_night = .FALSE.
	389	ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sss' ) THEN
	390	IF ( nn_2dint_sss == -1 ) THEN
	391	n2dint_type = nn_2dint_default
	392	ELSE
	393	n2dint_type = nn_2dint_sss
	394	ENDIF
	395	ztype_avglamscl = rn_sss_avglamscl
	396	ztype_avgphiscl = rn_sss_avgphiscl
	397	ltype_fp_indegs = ln_sss_fp_indegs
	398	ltype_night = .FALSE.
	399	ELSE
	400	n2dint_type = nn_2dint_default
	401	ztype_avglamscl = rn_default_avglamscl
	402	ztype_avgphiscl = rn_default_avgphiscl
	403	ltype_fp_indegs = ln_default_fp_indegs
	404	ltype_night = .FALSE.
	405	ENDIF
	406
	407	CALL obs_setinterpopts( nsurftypes, jtype, TRIM(cobstypessurf(jtype)), &
	408	& nn_2dint_default, n2dint_type, &
	409	& ztype_avglamscl, ztype_avgphiscl, &
	410	& ltype_fp_indegs, ltype_night, &
	411	& n2dintsurf, zavglamscl, zavgphiscl, &
	412	& lfpindegs, llnightav )
	413
	414	END DO
	415	!
[2128]	416	ENDIF
	417
	418
[6140]	419	!-----------------------------------------------------------------------
	420	! Obs operator parameter checking and initialisations
	421	!-----------------------------------------------------------------------
[9168]	422	!
	423	IF( ln_vel3d .AND. .NOT.ln_grid_global ) THEN
[6140]	424	CALL ctl_stop( 'Velocity data only works with ln_grid_global=.true.' )
	425	RETURN
	426	ENDIF
[9168]	427	!
	428	IF( ln_grid_global ) THEN
	429	CALL ctl_warn( 'dia_obs_init: ln_grid_global=T may cause memory issues when used with a large number of processors' )
[2128]	430	ENDIF
[9168]	431	!
	432	IF( nn_1dint < 0 .OR. nn_1dint > 1 ) THEN
	433	CALL ctl_stop('dia_obs_init: Choice of vertical (1D) interpolation method is not available')
[2128]	434	ENDIF
[9168]	435	!
[14056]	436	IF( nn_2dint_default < 0 .OR. nn_2dint_default > 6 ) THEN
	437	CALL ctl_stop('dia_obs_init: Choice of default horizontal (2D) interpolation method is not available')
[6140]	438	ENDIF
[9168]	439	!
[6140]	440	CALL obs_typ_init
[9168]	441	IF( ln_grid_global ) CALL mppmap_init
	442	!
[6140]	443	CALL obs_grid_setup( )
[2128]	444
[6140]	445	!-----------------------------------------------------------------------
	446	! Depending on switches read the various observation types
	447	!-----------------------------------------------------------------------
[9168]	448	!
	449	IF( nproftypes > 0 ) THEN
	450	!
	451	ALLOCATE( profdata (nproftypes) , nvarsprof (nproftypes) )
	452	ALLOCATE( profdataqc(nproftypes) , nextrprof (nproftypes) )
	453	!
[6140]	454	DO jtype = 1, nproftypes
[9168]	455	!
[6140]	456	IF ( TRIM(cobstypesprof(jtype)) == 'prof' ) THEN
[14056]	457	nvarsprof(jtype) = 2
	458	nextrprof(jtype) = 1
	459	ALLOCATE( llvar (nvarsprof(jtype)) )
	460	ALLOCATE( clvars(nvarsprof(jtype)) )
	461	ALLOCATE( zglam(jpi, jpj, nvarsprof(jtype)) )
	462	ALLOCATE( zgphi(jpi, jpj, nvarsprof(jtype)) )
	463	ALLOCATE( zmask(jpi, jpj, jpk, nvarsprof(jtype)) )
	464	llvar(1) = ln_t3d
	465	llvar(2) = ln_s3d
	466	clvars(1) = 'POTM'
	467	clvars(2) = 'PSAL'
	468	zglam(:,:,1) = glamt(:,:)
	469	zglam(:,:,2) = glamt(:,:)
	470	zgphi(:,:,1) = gphit(:,:)
	471	zgphi(:,:,2) = gphit(:,:)
	472	zmask(:,:,:,1) = tmask(:,:,:)
	473	zmask(:,:,:,2) = tmask(:,:,:)
	474	ELSE IF ( TRIM(cobstypesprof(jtype)) == 'vel' ) THEN
	475	nvarsprof(jtype) = 2
[6140]	476	nextrprof(jtype) = 2
[14056]	477	ALLOCATE( llvar (nvarsprof(jtype)) )
	478	ALLOCATE( clvars(nvarsprof(jtype)) )
	479	ALLOCATE( zglam(jpi, jpj, nvarsprof(jtype)) )
	480	ALLOCATE( zgphi(jpi, jpj, nvarsprof(jtype)) )
	481	ALLOCATE( zmask(jpi, jpj, jpk, nvarsprof(jtype)) )
	482	llvar(1) = ln_vel3d
	483	llvar(2) = ln_vel3d
	484	clvars(1) = 'UVEL'
	485	clvars(2) = 'VVEL'
	486	zglam(:,:,1) = glamu(:,:)
	487	zglam(:,:,2) = glamv(:,:)
	488	zgphi(:,:,1) = gphiu(:,:)
	489	zgphi(:,:,2) = gphiv(:,:)
	490	zmask(:,:,:,1) = umask(:,:,:)
	491	zmask(:,:,:,2) = vmask(:,:,:)
	492	ELSE
	493	nvarsprof(jtype) = 1
	494	nextrprof(jtype) = 0
	495	ALLOCATE( llvar (nvarsprof(jtype)) )
	496	ALLOCATE( clvars(nvarsprof(jtype)) )
	497	ALLOCATE( zglam(jpi, jpj, nvarsprof(jtype)) )
	498	ALLOCATE( zgphi(jpi, jpj, nvarsprof(jtype)) )
	499	ALLOCATE( zmask(jpi, jpj, jpk, nvarsprof(jtype)) )
	500	llvar(1) = .TRUE.
	501	zglam(:,:,1) = glamt(:,:)
	502	zgphi(:,:,1) = gphit(:,:)
	503	zmask(:,:,:,1) = tmask(:,:,:)
[6140]	504	ENDIF
[9168]	505	!
	506	! Read in profile or profile obs types
[6140]	507	CALL obs_rea_prof( profdata(jtype), ifilesprof(jtype), &
	508	& clproffiles(jtype,1:ifilesprof(jtype)), &
	509	& nvarsprof(jtype), nextrprof(jtype), nitend-nit000+2, &
[14056]	510	& rn_dobsini, rn_dobsend, llvar, &
	511	& ln_ignmis, ln_s_at_t, .FALSE., clvars, &
[6140]	512	& kdailyavtypes = nn_profdavtypes )
[9168]	513	!
[6140]	514	DO jvar = 1, nvarsprof(jtype)
	515	CALL obs_prof_staend( profdata(jtype), jvar )
	516	END DO
[9168]	517	!
[6140]	518	CALL obs_pre_prof( profdata(jtype), profdataqc(jtype), &
[14056]	519	& llvar, &
[6140]	520	& jpi, jpj, jpk, &
[14056]	521	& zmask, zglam, zgphi, &
[12377]	522	& ln_nea, ln_bound_reject, Kmm, &
[9023]	523	& kdailyavtypes = nn_profdavtypes )
[14056]	524	!
	525	DEALLOCATE( llvar, clvars, zglam, zgphi, zmask )
	526	!
[6140]	527	END DO
[9168]	528	!
[6140]	529	DEALLOCATE( ifilesprof, clproffiles )
[9168]	530	!
[2128]	531	ENDIF
[9168]	532	!
	533	IF( nsurftypes > 0 ) THEN
	534	!
	535	ALLOCATE( surfdata (nsurftypes) , nvarssurf(nsurftypes) )
	536	ALLOCATE( surfdataqc(nsurftypes) , nextrsurf(nsurftypes) )
	537	!
[6140]	538	DO jtype = 1, nsurftypes
[9168]	539	!
[6140]	540	nvarssurf(jtype) = 1
	541	nextrsurf(jtype) = 0
[9023]	542	llnightav(jtype) = .FALSE.
[9168]	543	IF( TRIM(cobstypessurf(jtype)) == 'sla' ) nextrsurf(jtype) = 2
	544	IF( TRIM(cobstypessurf(jtype)) == 'sst' ) llnightav(jtype) = ln_sstnight
	545	!
[14056]	546	ALLOCATE( clvars( nvarssurf(jtype) ) )
	547	IF ( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN
	548	clvars(1) = 'SLA'
	549	ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sst' ) THEN
	550	clvars(1) = 'SST'
	551	ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sic' ) THEN
	552	clvars(1) = 'ICECONC'
	553	ELSE IF ( TRIM(cobstypessurf(jtype)) == 'sss' ) THEN
	554	clvars(1) = 'SSS'
	555	ENDIF
	556	!
[9168]	557	! Read in surface obs types
[6140]	558	CALL obs_rea_surf( surfdata(jtype), ifilessurf(jtype), &
	559	& clsurffiles(jtype,1:ifilessurf(jtype)), &
	560	& nvarssurf(jtype), nextrsurf(jtype), nitend-nit000+2, &
[14056]	561	& rn_dobsini, rn_dobsend, ln_ignmis, .FALSE., llnightav(jtype), &
	562	& clvars )
[9168]	563	!
[9023]	564	CALL obs_pre_surf( surfdata(jtype), surfdataqc(jtype), ln_nea, ln_bound_reject )
[9168]	565	!
	566	IF( TRIM(cobstypessurf(jtype)) == 'sla' ) THEN
[12377]	567	CALL obs_rea_mdt( surfdataqc(jtype), n2dintsurf(jtype), Kmm )
[9168]	568	IF( ln_altbias ) &
[9023]	569	& CALL obs_rea_altbias ( surfdataqc(jtype), n2dintsurf(jtype), cn_altbiasfile )
[6140]	570	ENDIF
[9168]	571	!
	572	IF( TRIM(cobstypessurf(jtype)) == 'sst' .AND. ln_sstbias ) THEN
[9023]	573	jnumsstbias = 0
	574	DO jfile = 1, jpmaxnfiles
[9168]	575	IF( TRIM(cn_sstbiasfiles(jfile)) /= '' ) jnumsstbias = jnumsstbias + 1
[9023]	576	END DO
[9168]	577	IF( jnumsstbias == 0 ) CALL ctl_stop("ln_sstbias set but no bias files to read in")
	578	!
	579	CALL obs_app_sstbias( surfdataqc(jtype), n2dintsurf(jtype) , &
	580	& jnumsstbias , cn_sstbiasfiles(1:jnumsstbias) )
[6140]	581	ENDIF
[14056]	582	!
	583	DEALLOCATE( clvars )
[6140]	584	END DO
[9168]	585	!
[6140]	586	DEALLOCATE( ifilessurf, clsurffiles )
[9168]	587	!
[6140]	588	ENDIF
[9168]	589	!
[2128]	590	END SUBROUTINE dia_obs_init
	591
[9019]	592
[12377]	593	SUBROUTINE dia_obs( kstp, Kmm )
[2128]	594	!!----------------------------------------------------------------------
	595	!! * ROUTINE dia_obs *
	596	!!
	597	!! ** Purpose : Call the observation operators on each time step
	598	!!
	599	!! ** Method : Call the observation operators on each time step to
[6140]	600	!! compute the model equivalent of the following data:
	601	!! - Profile data, currently T/S or U/V
	602	!! - Surface data, currently SST, SLA or sea-ice concentration.
[2128]	603	!!
[6140]	604	!! ** Action :
[2128]	605	!!----------------------------------------------------------------------
[12377]	606	USE dom_oce, ONLY : gdept, gdept_1d ! Ocean space domain variables (Kmm time-level only)
[9019]	607	USE phycst , ONLY : rday ! Physical constants
[12377]	608	USE oce , ONLY : ts, uu, vv, ssh ! Ocean dynamics and tracers variables (Kmm time-level only)
[9019]	609	USE phycst , ONLY : rday ! Physical constants
[9570]	610	#if defined key_si3
[9656]	611	USE ice , ONLY : at_i ! SI3 Ice model variables
[2128]	612	#endif
[9023]	613	#if defined key_cice
	614	USE sbc_oce, ONLY : fr_i ! ice fraction
	615	#endif
	616
[2128]	617	IMPLICIT NONE
	618
	619	!! * Arguments
[6140]	620	INTEGER, INTENT(IN) :: kstp ! Current timestep
[12377]	621	INTEGER, INTENT(in) :: Kmm ! ocean time level indices
[2128]	622	!! * Local declarations
[6140]	623	INTEGER :: idaystp ! Number of timesteps per day
	624	INTEGER :: jtype ! Data loop variable
	625	INTEGER :: jvar ! Variable number
[14982]	626	INTEGER :: ji, jj, jk ! Loop counters
[14056]	627	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
	628	& zprofvar ! Model values for variables in a prof ob
	629	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
	630	& zprofmask ! Mask associated with zprofvar
[9125]	631	REAL(wp), DIMENSION(jpi,jpj) :: &
[9023]	632	& zsurfvar, & ! Model values equivalent to surface ob.
	633	& zsurfmask ! Mask associated with surface variable
[14056]	634	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
	635	& zglam, & ! Model longitudes for prof variables
	636	& zgphi ! Model latitudes for prof variables
[14982]	637	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zdept, zdepw
[6140]	638
[9019]	639	!-----------------------------------------------------------------------
[2715]	640
[2128]	641	IF(lwp) THEN
	642	WRITE(numout,*)
	643	WRITE(numout,*) 'dia_obs : Call the observation operators', kstp
	644	WRITE(numout,*) '~~~~~~~'
	645	ENDIF
	646
[12489]	647	idaystp = NINT( rday / rn_Dt )
[2128]	648
	649	!-----------------------------------------------------------------------
[6140]	650	! Call the profile and surface observation operators
[2128]	651	!-----------------------------------------------------------------------
	652
[6140]	653	IF ( nproftypes > 0 ) THEN
	654
[14982]	655	ALLOCATE( zdept(jpi,jpj,jpk), zdepw(jpi,jpj,jpk) )
	656	DO jk = 1, jpk
	657	zdept(:,:,jk) = gdept(:,:,jk,Kmm)
	658	zdepw(:,:,jk) = gdepw(:,:,jk,Kmm)
	659	END DO
	660
[6140]	661	DO jtype = 1, nproftypes
	662
[14056]	663	! Allocate local work arrays
	664	ALLOCATE( zprofvar (jpi, jpj, jpk, profdataqc(jtype)%nvar) )
	665	ALLOCATE( zprofmask(jpi, jpj, jpk, profdataqc(jtype)%nvar) )
	666	ALLOCATE( zglam (jpi, jpj, profdataqc(jtype)%nvar) )
	667	ALLOCATE( zgphi (jpi, jpj, profdataqc(jtype)%nvar) )
	668
	669	! Defaults which might change
	670	DO jvar = 1, profdataqc(jtype)%nvar
	671	zprofmask(:,:,:,jvar) = tmask(:,:,:)
	672	zglam(:,:,jvar) = glamt(:,:)
	673	zgphi(:,:,jvar) = gphit(:,:)
	674	END DO
	675
[6140]	676	SELECT CASE ( TRIM(cobstypesprof(jtype)) )
	677	CASE('prof')
[14056]	678	zprofvar(:,:,:,1) = ts(:,:,:,jp_tem,Kmm)
	679	zprofvar(:,:,:,2) = ts(:,:,:,jp_sal,Kmm)
[6140]	680	CASE('vel')
[14056]	681	zprofvar(:,:,:,1) = uu(:,:,:,Kmm)
	682	zprofvar(:,:,:,2) = vv(:,:,:,Kmm)
	683	zprofmask(:,:,:,1) = umask(:,:,:)
	684	zprofmask(:,:,:,2) = vmask(:,:,:)
	685	zglam(:,:,1) = glamu(:,:)
	686	zglam(:,:,2) = glamv(:,:)
	687	zgphi(:,:,1) = gphiu(:,:)
	688	zgphi(:,:,2) = gphiv(:,:)
[9023]	689	CASE DEFAULT
	690	CALL ctl_stop( 'Unknown profile observation type '//TRIM(cobstypesprof(jtype))//' in dia_obs' )
[6140]	691	END SELECT
	692
[14056]	693	DO jvar = 1, profdataqc(jtype)%nvar
	694	CALL obs_prof_opt( profdataqc(jtype), kstp, jpi, jpj, jpk, &
	695	& nit000, idaystp, jvar, &
	696	& zprofvar(:,:,:,jvar), &
[14982]	697	& zdept(:,:,:), zdepw(:,:,:), &
[14056]	698	& zprofmask(:,:,:,jvar), &
	699	& zglam(:,:,jvar), zgphi(:,:,jvar), &
	700	& nn_1dint, nn_2dint_default, &
	701	& kdailyavtypes = nn_profdavtypes )
	702	END DO
	703
	704	DEALLOCATE( zprofvar, zprofmask, zglam, zgphi )
[6140]	705
[2128]	706	END DO
	707
[14982]	708	DEALLOCATE( zdept, zdepw )
	709
[2128]	710	ENDIF
	711
[6140]	712	IF ( nsurftypes > 0 ) THEN
[2128]	713
[6140]	714	DO jtype = 1, nsurftypes
[2128]	715
[9023]	716	!Defaults which might be changed
	717	zsurfmask(:,:) = tmask(:,:,1)
	718
[6140]	719	SELECT CASE ( TRIM(cobstypessurf(jtype)) )
	720	CASE('sst')
[12377]	721	zsurfvar(:,:) = ts(:,:,1,jp_tem,Kmm)
[6140]	722	CASE('sla')
[12377]	723	zsurfvar(:,:) = ssh(:,:,Kmm)
[9023]	724	CASE('sss')
[12377]	725	zsurfvar(:,:) = ts(:,:,1,jp_sal,Kmm)
[6140]	726	CASE('sic')
	727	IF ( kstp == 0 ) THEN
	728	IF ( lwp .AND. surfdataqc(jtype)%nsstpmpp(1) > 0 ) THEN
	729	CALL ctl_warn( 'Sea-ice not initialised on zeroth '// &
	730	& 'time-step but some obs are valid then.' )
	731	WRITE(numout,*)surfdataqc(jtype)%nsstpmpp(1), &
	732	& ' sea-ice obs will be missed'
	733	ENDIF
	734	surfdataqc(jtype)%nsurfup = surfdataqc(jtype)%nsurfup + &
	735	& surfdataqc(jtype)%nsstp(1)
	736	CYCLE
	737	ELSE
[9570]	738	#if defined key_cice \|\| defined key_si3
[9023]	739	zsurfvar(:,:) = fr_i(:,:)
[9354]	740	#else
[9040]	741	CALL ctl_stop( ' Trying to run sea-ice observation operator', &
	742	& ' but no sea-ice model appears to have been defined' )
[9023]	743	#endif
[6140]	744	ENDIF
	745
	746	END SELECT
[2128]	747
[6140]	748	CALL obs_surf_opt( surfdataqc(jtype), kstp, jpi, jpj, &
[9023]	749	& nit000, idaystp, zsurfvar, zsurfmask, &
	750	& n2dintsurf(jtype), llnightav(jtype), &
	751	& zavglamscl(jtype), zavgphiscl(jtype), &
	752	& lfpindegs(jtype) )
[6140]	753
[2128]	754	END DO
[6140]	755
[2128]	756	ENDIF
	757
	758	END SUBROUTINE dia_obs
[6140]	759
	760	SUBROUTINE dia_obs_wri
[2128]	761	!!----------------------------------------------------------------------
	762	!! * ROUTINE dia_obs_wri *
	763	!!
	764	!! ** Purpose : Call observation diagnostic output routines
	765	!!
	766	!! ** Method : Call observation diagnostic output routines
	767	!!
[6140]	768	!! ** Action :
[2128]	769	!!
	770	!! History :
	771	!! ! 06-03 (K. Mogensen) Original code
	772	!! ! 06-05 (K. Mogensen) Reformatted
	773	!! ! 06-10 (A. Weaver) Cleaning
	774	!! ! 07-03 (K. Mogensen) General handling of profiles
	775	!! ! 08-09 (M. Valdivieso) Velocity component (U,V) profiles
[6140]	776	!! ! 15-08 (M. Martin) Combined writing for prof and surf types
[2128]	777	!!----------------------------------------------------------------------
[6140]	778	!! * Modules used
	779	USE obs_rot_vel ! Rotation of velocities
	780
[2128]	781	IMPLICIT NONE
	782
	783	!! * Local declarations
[6140]	784	INTEGER :: jtype ! Data set loop variable
	785	INTEGER :: jo, jvar, jk
	786	REAL(wp), DIMENSION(:), ALLOCATABLE :: &
	787	& zu, &
	788	& zv
[2128]	789
	790	!-----------------------------------------------------------------------
	791	! Depending on switches call various observation output routines
	792	!-----------------------------------------------------------------------
	793
[6140]	794	IF ( nproftypes > 0 ) THEN
[2128]	795
[6140]	796	DO jtype = 1, nproftypes
[2128]	797
[6140]	798	IF ( TRIM(cobstypesprof(jtype)) == 'vel' ) THEN
[2128]	799
[6140]	800	! For velocity data, rotate the model velocities to N/S, E/W
	801	! using the compressed data structure.
	802	ALLOCATE( &
	803	& zu(profdataqc(jtype)%nvprot(1)), &
	804	& zv(profdataqc(jtype)%nvprot(2)) &
	805	& )
[2128]	806
[14056]	807	CALL obs_rotvel( profdataqc(jtype), nn_2dint_default, zu, zv )
[2128]	808
[6140]	809	DO jo = 1, profdataqc(jtype)%nprof
	810	DO jvar = 1, 2
	811	DO jk = profdataqc(jtype)%npvsta(jo,jvar), profdataqc(jtype)%npvend(jo,jvar)
[2128]	812
[6140]	813	IF ( jvar == 1 ) THEN
	814	profdataqc(jtype)%var(jvar)%vmod(jk) = zu(jk)
	815	ELSE
	816	profdataqc(jtype)%var(jvar)%vmod(jk) = zv(jk)
	817	ENDIF
[2128]	818
[6140]	819	END DO
	820	END DO
	821	END DO
[2128]	822
[6140]	823	DEALLOCATE( zu )
	824	DEALLOCATE( zv )
[2128]	825
[6140]	826	END IF
[2128]	827
[6140]	828	CALL obs_prof_decompress( profdataqc(jtype), &
	829	& profdata(jtype), .TRUE., numout )
[2128]	830
[6140]	831	CALL obs_wri_prof( profdata(jtype) )
[2128]	832
	833	END DO
	834
	835	ENDIF
	836
[6140]	837	IF ( nsurftypes > 0 ) THEN
[2128]	838
[6140]	839	DO jtype = 1, nsurftypes
[2128]	840
[6140]	841	CALL obs_surf_decompress( surfdataqc(jtype), &
	842	& surfdata(jtype), .TRUE., numout )
[2128]	843
[6140]	844	CALL obs_wri_surf( surfdata(jtype) )
[2128]	845
	846	END DO
	847
	848	ENDIF
	849
	850	END SUBROUTINE dia_obs_wri
	851
[4245]	852	SUBROUTINE dia_obs_dealloc
	853	IMPLICIT NONE
	854	!!----------------------------------------------------------------------
	855	!! * ROUTINE dia_obs_dealloc *
	856	!!
	857	!! ** Purpose : To deallocate data to enable the obs_oper online loop.
	858	!! Specifically: dia_obs_init --> dia_obs --> dia_obs_wri
	859	!!
	860	!! ** Method : Clean up various arrays left behind by the obs_oper.
	861	!!
	862	!! ** Action :
	863	!!
	864	!!----------------------------------------------------------------------
[6140]	865	! obs_grid deallocation
[4245]	866	CALL obs_grid_deallocate
	867
[6140]	868	! diaobs deallocation
	869	IF ( nproftypes > 0 ) &
	870	& DEALLOCATE( cobstypesprof, profdata, profdataqc, nvarsprof, nextrprof )
	871
	872	IF ( nsurftypes > 0 ) &
[9023]	873	& DEALLOCATE( cobstypessurf, surfdata, surfdataqc, nvarssurf, nextrsurf, &
	874	& n2dintsurf, zavglamscl, zavgphiscl, lfpindegs, llnightav )
[6140]	875
[4245]	876	END SUBROUTINE dia_obs_dealloc
	877
[6140]	878	SUBROUTINE calc_date( kstp, ddobs )
[2128]	879	!!----------------------------------------------------------------------
[6140]	880	!! * ROUTINE calc_date *
[2128]	881	!!
[6140]	882	!! ** Purpose : Get date in double precision YYYYMMDD.HHMMSS format
[2128]	883	!!
[6140]	884	!! ** Method : Get date in double precision YYYYMMDD.HHMMSS format
[2128]	885	!!
[6140]	886	!! ** Action : Get date in double precision YYYYMMDD.HHMMSS format
[2128]	887	!!
[6140]	888	!! ** Action : Get initial date in double precision YYYYMMDD.HHMMSS format
	889	!!
[2128]	890	!! History :
	891	!! ! 06-03 (K. Mogensen) Original code
	892	!! ! 06-05 (K. Mogensen) Reformatted
	893	!! ! 06-10 (A. Weaver) Cleaning
	894	!! ! 06-10 (G. Smith) Calculates initial date the same as method for final date
	895	!! ! 10-05 (D. Lea) Update to month length calculation for NEMO vn3.2
[6140]	896	!! ! 2014-09 (D. Lea) New generic routine now deals with arbitrary initial time of day
[2128]	897	!!----------------------------------------------------------------------
	898	USE phycst, ONLY : & ! Physical constants
	899	& rday
	900	USE dom_oce, ONLY : & ! Ocean space and time domain variables
[12489]	901	& rn_Dt
[2128]	902
	903	IMPLICIT NONE
	904
	905	!! * Arguments
[6140]	906	REAL(KIND=dp), INTENT(OUT) :: ddobs ! Date in YYYYMMDD.HHMMSS
	907	INTEGER :: kstp
[2128]	908
	909	!! * Local declarations
	910	INTEGER :: iyea ! date - (year, month, day, hour, minute)
	911	INTEGER :: imon
	912	INTEGER :: iday
	913	INTEGER :: ihou
	914	INTEGER :: imin
[6140]	915	INTEGER :: imday ! Number of days in month.
	916	REAL(wp) :: zdayfrc ! Fraction of day
[2128]	917
	918	INTEGER, DIMENSION(12) :: imonth_len !: length in days of the months of the current year
	919
	920	!!----------------------------------------------------------------------
	921	!! Initial date initialization (year, month, day, hour, minute)
	922	!!----------------------------------------------------------------------
	923	iyea = ndate0 / 10000
	924	imon = ( ndate0 - iyea * 10000 ) / 100
	925	iday = ndate0 - iyea * 10000 - imon * 100
[6140]	926	ihou = nn_time0 / 100
	927	imin = ( nn_time0 - ihou * 100 )
[2128]	928
	929	!!----------------------------------------------------------------------
	930	!! Compute number of days + number of hours + min since initial time
	931	!!----------------------------------------------------------------------
[12489]	932	zdayfrc = kstp * rn_Dt / rday
[2128]	933	zdayfrc = zdayfrc - aint(zdayfrc)
[6140]	934	imin = imin + int( zdayfrc * 24 * 60 )
	935	DO WHILE (imin >= 60)
	936	imin=imin-60
	937	ihou=ihou+1
	938	END DO
	939	DO WHILE (ihou >= 24)
	940	ihou=ihou-24
	941	iday=iday+1
	942	END DO
[12489]	943	iday = iday + kstp * rn_Dt / rday
[2128]	944
[6140]	945	!-----------------------------------------------------------------------
	946	! Convert number of days (iday) into a real date
	947	!----------------------------------------------------------------------
[2128]	948
	949	CALL calc_month_len( iyea, imonth_len )
[6140]	950
[2128]	951	DO WHILE ( iday > imonth_len(imon) )
	952	iday = iday - imonth_len(imon)
	953	imon = imon + 1
	954	IF ( imon > 12 ) THEN
	955	imon = 1
	956	iyea = iyea + 1
	957	CALL calc_month_len( iyea, imonth_len ) ! update month lengths
	958	ENDIF
	959	END DO
	960
[6140]	961	!----------------------------------------------------------------------
	962	! Convert it into YYYYMMDD.HHMMSS format.
	963	!----------------------------------------------------------------------
	964	ddobs = iyea * 10000_dp + imon * 100_dp + &
	965	& iday + ihou * 0.01_dp + imin * 0.0001_dp
	966
	967	END SUBROUTINE calc_date
	968
	969	SUBROUTINE ini_date( ddobsini )
[2128]	970	!!----------------------------------------------------------------------
[6140]	971	!! * ROUTINE ini_date *
	972	!!
	973	!! ** Purpose : Get initial date in double precision YYYYMMDD.HHMMSS format
	974	!!
	975	!! ** Method :
	976	!!
	977	!! ** Action :
	978	!!
	979	!! History :
	980	!! ! 06-03 (K. Mogensen) Original code
	981	!! ! 06-05 (K. Mogensen) Reformatted
	982	!! ! 06-10 (A. Weaver) Cleaning
	983	!! ! 10-05 (D. Lea) Update to month length calculation for NEMO vn3.2
	984	!! ! 2014-09 (D. Lea) Change to call generic routine calc_date
[2128]	985	!!----------------------------------------------------------------------
	986
[6140]	987	IMPLICIT NONE
[2128]	988
[6140]	989	!! * Arguments
	990	REAL(KIND=dp), INTENT(OUT) :: ddobsini ! Initial date in YYYYMMDD.HHMMSS
	991
	992	CALL calc_date( nit000 - 1, ddobsini )
	993
[2128]	994	END SUBROUTINE ini_date
	995
	996	SUBROUTINE fin_date( ddobsfin )
	997	!!----------------------------------------------------------------------
	998	!! * ROUTINE fin_date *
	999	!!
[6140]	1000	!! ** Purpose : Get final date in double precision YYYYMMDD.HHMMSS format
[2128]	1001	!!
[6140]	1002	!! ** Method :
[2128]	1003	!!
[6140]	1004	!! ** Action :
[2128]	1005	!!
	1006	!! History :
	1007	!! ! 06-03 (K. Mogensen) Original code
	1008	!! ! 06-05 (K. Mogensen) Reformatted
	1009	!! ! 06-10 (A. Weaver) Cleaning
	1010	!! ! 10-05 (D. Lea) Update to month length calculation for NEMO vn3.2
[6140]	1011	!! ! 2014-09 (D. Lea) Change to call generic routine calc_date
[2128]	1012	!!----------------------------------------------------------------------
	1013
	1014	IMPLICIT NONE
	1015
	1016	!! * Arguments
[6140]	1017	REAL(dp), INTENT(OUT) :: ddobsfin ! Final date in YYYYMMDD.HHMMSS
[2128]	1018
[6140]	1019	CALL calc_date( nitend, ddobsfin )
[2128]	1020
[6140]	1021	END SUBROUTINE fin_date
	1022
[14056]	1023	SUBROUTINE obs_settypefiles( ntypes, jpmaxnfiles, ifiles, cobstypes, cfiles )
[9023]	1024
[14056]	1025	INTEGER, INTENT(IN) :: ntypes ! Total number of obs types
	1026	INTEGER, INTENT(IN) :: jpmaxnfiles ! Maximum number of files allowed for each type
	1027	INTEGER, DIMENSION(ntypes), INTENT(OUT) :: &
	1028	& ifiles ! Out number of files for each type
	1029	CHARACTER(len=8), DIMENSION(ntypes), INTENT(IN) :: &
	1030	& cobstypes ! List of obs types
	1031	CHARACTER(len=128), DIMENSION(ntypes, jpmaxnfiles), INTENT(IN) :: &
	1032	& cfiles ! List of files for all types
[9023]	1033
[14056]	1034	!Local variables
	1035	INTEGER :: jfile
	1036	INTEGER :: jtype
[9023]	1037
[14056]	1038	DO jtype = 1, ntypes
[9023]	1039
[14056]	1040	ifiles(jtype) = 0
	1041	DO jfile = 1, jpmaxnfiles
	1042	IF ( trim(cfiles(jtype,jfile)) /= '' ) &
	1043	ifiles(jtype) = ifiles(jtype) + 1
	1044	END DO
[9023]	1045
[14056]	1046	IF ( ifiles(jtype) == 0 ) THEN
	1047	CALL ctl_stop( 'Logical for observation type '//TRIM(cobstypes(jtype))// &
	1048	& ' set to true but no files available to read' )
	1049	ENDIF
[9023]	1050
[14056]	1051	IF(lwp) THEN
	1052	WRITE(numout,*) ' '//cobstypes(jtype)//' input observation file names:'
	1053	DO jfile = 1, ifiles(jtype)
	1054	WRITE(numout,*) ' '//TRIM(cfiles(jtype,jfile))
	1055	END DO
	1056	ENDIF
[9023]	1057
[14056]	1058	END DO
[9023]	1059
[14056]	1060	END SUBROUTINE obs_settypefiles
[9023]	1061
[14056]	1062	SUBROUTINE obs_setinterpopts( ntypes, jtype, ctypein, &
	1063	& n2dint_default, n2dint_type, &
	1064	& ravglamscl_type, ravgphiscl_type, &
	1065	& lfp_indegs_type, lavnight_type, &
	1066	& n2dint, ravglamscl, ravgphiscl, &
	1067	& lfpindegs, lavnight )
[9023]	1068
[14056]	1069	INTEGER, INTENT(IN) :: ntypes ! Total number of obs types
	1070	INTEGER, INTENT(IN) :: jtype ! Index of the current type of obs
	1071	INTEGER, INTENT(IN) :: n2dint_default ! Default option for interpolation type
	1072	INTEGER, INTENT(IN) :: n2dint_type ! Option for interpolation type
	1073	REAL(wp), INTENT(IN) :: &
	1074	& ravglamscl_type, & !E/W diameter of obs footprint for this type
	1075	& ravgphiscl_type !N/S diameter of obs footprint for this type
	1076	LOGICAL, INTENT(IN) :: lfp_indegs_type !T=> footprint in degrees, F=> in metres
	1077	LOGICAL, INTENT(IN) :: lavnight_type !T=> obs represent night time average
	1078	CHARACTER(len=8), INTENT(IN) :: ctypein
[9023]	1079
[14056]	1080	INTEGER, DIMENSION(ntypes), INTENT(INOUT) :: &
	1081	& n2dint
	1082	REAL(wp), DIMENSION(ntypes), INTENT(INOUT) :: &
	1083	& ravglamscl, ravgphiscl
	1084	LOGICAL, DIMENSION(ntypes), INTENT(INOUT) :: &
	1085	& lfpindegs, lavnight
[9023]	1086
[14056]	1087	lavnight(jtype) = lavnight_type
[9023]	1088
[14056]	1089	IF ( (n2dint_type >= 0) .AND. (n2dint_type <= 6) ) THEN
	1090	n2dint(jtype) = n2dint_type
	1091	ELSE IF ( n2dint_type == -1 ) THEN
	1092	n2dint(jtype) = n2dint_default
	1093	ELSE
	1094	CALL ctl_stop(' Choice of '//TRIM(ctypein)//' horizontal (2D) interpolation method', &
	1095	& ' is not available')
	1096	ENDIF
[9023]	1097
[14056]	1098	! For averaging observation footprints set options for size of footprint
	1099	IF ( (n2dint(jtype) > 4) .AND. (n2dint(jtype) <= 6) ) THEN
	1100	IF ( ravglamscl_type > 0._wp ) THEN
	1101	ravglamscl(jtype) = ravglamscl_type
	1102	ELSE
	1103	CALL ctl_stop( 'Incorrect value set for averaging footprint '// &
	1104	'scale (ravglamscl) for observation type '//TRIM(ctypein) )
	1105	ENDIF
[9023]	1106
[14056]	1107	IF ( ravgphiscl_type > 0._wp ) THEN
	1108	ravgphiscl(jtype) = ravgphiscl_type
	1109	ELSE
	1110	CALL ctl_stop( 'Incorrect value set for averaging footprint '// &
	1111	'scale (ravgphiscl) for observation type '//TRIM(ctypein) )
	1112	ENDIF
[9023]	1113
[14056]	1114	lfpindegs(jtype) = lfp_indegs_type
[9023]	1115
[14056]	1116	ENDIF
[9023]	1117
[14056]	1118	! Write out info
	1119	IF(lwp) THEN
	1120	IF ( n2dint(jtype) <= 4 ) THEN
	1121	WRITE(numout,*) ' '//TRIM(ctypein)// &
	1122	& ' model counterparts will be interpolated horizontally'
	1123	ELSE IF ( n2dint(jtype) <= 6 ) THEN
	1124	WRITE(numout,*) ' '//TRIM(ctypein)// &
	1125	& ' model counterparts will be averaged horizontally'
	1126	WRITE(numout,*) ' '//' with E/W scale: ',ravglamscl(jtype)
	1127	WRITE(numout,*) ' '//' with N/S scale: ',ravgphiscl(jtype)
	1128	IF ( lfpindegs(jtype) ) THEN
	1129	WRITE(numout,*) ' '//' (in degrees)'
	1130	ELSE
	1131	WRITE(numout,*) ' '//' (in metres)'
	1132	ENDIF
	1133	ENDIF
	1134	ENDIF
[9023]	1135
[14056]	1136	END SUBROUTINE obs_setinterpopts
	1137
[2128]	1138	END MODULE diaobs

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

Download in other formats: