Context Navigation

← Previous Revision
Next Revision →
Normal
Revision Log

obs_prep.F90

Last change on this file was 11361, checked in by jcastill, 5 years ago
First version of the new observation branch - it compiles, but has not been tested
File size: 58.7 KB

Rev	Line
[2128]	1	MODULE obs_prep
	2	!!=====================================================================
	3	!! * MODULE obs_prep *
	4	!! Observation diagnostics: Prepare observation arrays: screening,
	5	!! sorting, coordinate search
	6	!!=====================================================================
	7
	8	!!---------------------------------------------------------------------
[6140]	9	!! obs_pre_prof : First level check and screening of profile observations
	10	!! obs_pre_surf : First level check and screening of surface observations
	11	!! obs_scr : Basic screening of the observations
	12	!! obs_coo_tim : Compute number of time steps to the observation time
	13	!! obs_sor : Sort the observation arrays
[2128]	14	!!---------------------------------------------------------------------
[11361]	15	!! * Modules used
	16	USE par_kind, ONLY : & ! Precision variables
	17	& wp
[2128]	18	USE in_out_manager ! I/O manager
	19	USE obs_profiles_def ! Definitions for storage arrays for profiles
	20	USE obs_surf_def ! Definitions for storage arrays for surface data
	21	USE obs_mpp, ONLY : & ! MPP support routines for observation diagnostics
	22	& obs_mpp_sum_integer, &
	23	& obs_mpp_sum_integers
	24	USE obs_inter_sup ! Interpolation support
	25	USE obs_oper ! Observation operators
[11361]	26	#if defined key_bdy
	27	USE bdy_oce, ONLY : & ! Boundary information
	28	idx_bdy, nb_bdy
	29	#endif
	30	USE lib_mpp, ONLY : &
	31	& ctl_warn, ctl_stop
[2128]	32
	33	IMPLICIT NONE
[11361]	34
	35	!! * Routine accessibility
[2128]	36	PRIVATE
	37
[11361]	38	PUBLIC &
	39	& obs_pre_prof, & ! First level check and screening of profile obs
	40	& obs_pre_surf, & ! First level check and screening of surface obs
	41	& calc_month_len ! Calculate the number of days in the months of a year
[2128]	42
[2287]	43	!!----------------------------------------------------------------------
	44	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
	45	!! $Id$
	46	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
	47	!!----------------------------------------------------------------------
	48
[2128]	49	CONTAINS
	50
[11361]	51	SUBROUTINE obs_pre_surf( surfdata, surfdataqc, ld_nea, ld_bound_reject, &
	52	kqc_cutoff )
[2128]	53	!!----------------------------------------------------------------------
	54	!! * ROUTINE obs_pre_sla *
	55	!!
[6140]	56	!! ** Purpose : First level check and screening of surface observations
[2128]	57	!!
[6140]	58	!! ** Method : First level check and screening of surface observations
[2128]	59	!!
	60	!! ** Action :
	61	!!
	62	!! References :
	63	!!
	64	!! History :
	65	!! ! 2007-03 (A. Weaver, K. Mogensen) Original
	66	!! ! 2007-06 (K. Mogensen et al) Reject obs. near land.
[6140]	67	!! ! 2015-02 (M. Martin) Combined routine for surface types.
[2128]	68	!!----------------------------------------------------------------------
[11361]	69	!! * Modules used
[2128]	70	USE par_oce ! Ocean parameters
[11361]	71	USE dom_oce, ONLY : & ! Geographical information
	72	& glamt, &
	73	& gphit, &
	74	& tmask, &
	75	& nproc
[2128]	76	!! * Arguments
[6140]	77	TYPE(obs_surf), INTENT(INOUT) :: surfdata ! Full set of surface data
[11361]	78	TYPE(obs_surf), INTENT(INOUT) :: surfdataqc ! Subset of surface data not failing screening
	79	LOGICAL, INTENT(IN) :: ld_nea ! Switch for rejecting observation near land
	80	LOGICAL, INTENT(IN) :: ld_bound_reject ! Switch for rejecting obs near the boundary
	81	INTEGER, INTENT(IN), OPTIONAL :: kqc_cutoff ! cut off for QC value
	82	!! * Local declarations
	83	INTEGER :: iqc_cutoff = 255 ! cut off for QC value
[2128]	84	INTEGER :: iyea0 ! Initial date
	85	INTEGER :: imon0 ! - (year, month, day, hour, minute)
	86	INTEGER :: iday0
	87	INTEGER :: ihou0
	88	INTEGER :: imin0
	89	INTEGER :: icycle ! Current assimilation cycle
	90	! Counters for observations that
	91	INTEGER :: iotdobs ! - outside time domain
	92	INTEGER :: iosdsobs ! - outside space domain
	93	INTEGER :: ilansobs ! - within a model land cell
	94	INTEGER :: inlasobs ! - close to land
	95	INTEGER :: igrdobs ! - fail the grid search
[11361]	96	INTEGER :: ibdysobs ! - close to open boundary
[2128]	97	! Global counters for observations that
	98	INTEGER :: iotdobsmpp ! - outside time domain
	99	INTEGER :: iosdsobsmpp ! - outside space domain
	100	INTEGER :: ilansobsmpp ! - within a model land cell
	101	INTEGER :: inlasobsmpp ! - close to land
	102	INTEGER :: igrdobsmpp ! - fail the grid search
[11361]	103	INTEGER :: ibdysobsmpp ! - close to open boundary
	104	LOGICAL, DIMENSION(:), ALLOCATABLE :: &
	105	& llvalid ! SLA data selection
[2128]	106	INTEGER :: jobs ! Obs. loop variable
	107	INTEGER :: jstp ! Time loop variable
	108	INTEGER :: inrc ! Time index variable
	109
[11361]	110	IF(lwp) WRITE(numout,*)'obs_pre_surf : Preparing the surface observations...'
	111	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
[6140]	112
[2128]	113	! Initial date initialization (year, month, day, hour, minute)
	114	iyea0 = ndate0 / 10000
	115	imon0 = ( ndate0 - iyea0 * 10000 ) / 100
	116	iday0 = ndate0 - iyea0 * 10000 - imon0 * 100
[11361]	117	ihou0 = 0
	118	imin0 = 0
[2128]	119
	120	icycle = no ! Assimilation cycle
	121
	122	! Diagnotics counters for various failures.
	123
	124	iotdobs = 0
	125	igrdobs = 0
	126	iosdsobs = 0
	127	ilansobs = 0
	128	inlasobs = 0
[11361]	129	ibdysobs = 0
[2128]	130
[11361]	131	! Set QC cutoff to optional value if provided
	132	IF ( PRESENT(kqc_cutoff) ) iqc_cutoff=kqc_cutoff
	133
[2128]	134	! -----------------------------------------------------------------------
[6140]	135	! Find time coordinate for surface data
[2128]	136	! -----------------------------------------------------------------------
	137
	138	CALL obs_coo_tim( icycle, &
	139	& iyea0, imon0, iday0, ihou0, imin0, &
[6140]	140	& surfdata%nsurf, surfdata%nyea, surfdata%nmon, &
	141	& surfdata%nday, surfdata%nhou, surfdata%nmin, &
	142	& surfdata%nqc, surfdata%mstp, iotdobs )
[2128]	143
	144	CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp )
	145
	146	! -----------------------------------------------------------------------
[6140]	147	! Check for surface data failing the grid search
[2128]	148	! -----------------------------------------------------------------------
	149
[6140]	150	CALL obs_coo_grd( surfdata%nsurf, surfdata%mi, surfdata%mj, &
	151	& surfdata%nqc, igrdobs )
[2128]	152
	153	CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp )
	154
	155	! -----------------------------------------------------------------------
	156	! Check for land points.
	157	! -----------------------------------------------------------------------
	158
[6140]	159	CALL obs_coo_spc_2d( surfdata%nsurf, &
[2128]	160	& jpi, jpj, &
[6140]	161	& surfdata%mi, surfdata%mj, &
	162	& surfdata%rlam, surfdata%rphi, &
[2128]	163	& glamt, gphit, &
[6140]	164	& tmask(:,:,1), surfdata%nqc, &
[2128]	165	& iosdsobs, ilansobs, &
[11361]	166	& inlasobs, ld_nea, &
	167	& ibdysobs, ld_bound_reject, &
	168	& iqc_cutoff )
[2128]	169
	170	CALL obs_mpp_sum_integer( iosdsobs, iosdsobsmpp )
	171	CALL obs_mpp_sum_integer( ilansobs, ilansobsmpp )
	172	CALL obs_mpp_sum_integer( inlasobs, inlasobsmpp )
[11361]	173	CALL obs_mpp_sum_integer( ibdysobs, ibdysobsmpp )
[2128]	174
	175	! -----------------------------------------------------------------------
[6140]	176	! Copy useful data from the surfdata data structure to
	177	! the surfdataqc data structure
[2128]	178	! -----------------------------------------------------------------------
	179
	180	! Allocate the selection arrays
	181
[6140]	182	ALLOCATE( llvalid(surfdata%nsurf) )
[2128]	183
[11361]	184	! We want all data which has qc flags <= iqc_cutoff
[2128]	185
[11361]	186	llvalid(:) = ( surfdata%nqc(:) <= iqc_cutoff )
[2128]	187
	188	! The actual copying
	189
[6140]	190	CALL obs_surf_compress( surfdata, surfdataqc, .TRUE., numout, &
[2128]	191	& lvalid=llvalid )
	192
	193	! Dellocate the selection arrays
	194	DEALLOCATE( llvalid )
	195
	196	! -----------------------------------------------------------------------
	197	! Print information about what observations are left after qc
	198	! -----------------------------------------------------------------------
	199
	200	! Update the total observation counter array
	201
	202	IF(lwp) THEN
	203	WRITE(numout,*)
[6140]	204	WRITE(numout,*) ' '//surfdataqc%cvars(1)//' data outside time domain = ', &
[2128]	205	& iotdobsmpp
[6140]	206	WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data that failed grid search = ', &
[2128]	207	& igrdobsmpp
[6140]	208	WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data outside space domain = ', &
[2128]	209	& iosdsobsmpp
[6140]	210	WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data at land points = ', &
[2128]	211	& ilansobsmpp
	212	IF (ld_nea) THEN
[6140]	213	WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data near land points (removed) = ', &
[2128]	214	& inlasobsmpp
	215	ELSE
[6140]	216	WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data near land points (kept) = ', &
[2128]	217	& inlasobsmpp
	218	ENDIF
[11361]	219	WRITE(numout,*) ' Remaining '//surfdataqc%cvars(1)//' data near open boundary (removed) = ', &
	220	& ibdysobsmpp
[6140]	221	WRITE(numout,*) ' '//surfdataqc%cvars(1)//' data accepted = ', &
	222	& surfdataqc%nsurfmpp
[2128]	223
	224	WRITE(numout,*)
	225	WRITE(numout,*) ' Number of observations per time step :'
	226	WRITE(numout,*)
[6140]	227	WRITE(numout,'(10X,A,10X,A)')'Time step',surfdataqc%cvars(1)
	228	WRITE(numout,'(10X,A,5X,A)')'---------','-----------------'
	229	CALL FLUSH(numout)
[2128]	230	ENDIF
	231
[6140]	232	DO jobs = 1, surfdataqc%nsurf
	233	inrc = surfdataqc%mstp(jobs) + 2 - nit000
	234	surfdataqc%nsstp(inrc) = surfdataqc%nsstp(inrc) + 1
[2128]	235	END DO
	236
[6140]	237	CALL obs_mpp_sum_integers( surfdataqc%nsstp, surfdataqc%nsstpmpp, &
[2128]	238	& nitend - nit000 + 2 )
	239
	240	IF ( lwp ) THEN
	241	DO jstp = nit000 - 1, nitend
	242	inrc = jstp - nit000 + 2
[6140]	243	WRITE(numout,1999) jstp, surfdataqc%nsstpmpp(inrc)
	244	CALL FLUSH(numout)
[2128]	245	END DO
	246	ENDIF
	247
	248	1999 FORMAT(10X,I9,5X,I17)
	249
[6140]	250	END SUBROUTINE obs_pre_surf
[2128]	251
	252
[11361]	253	SUBROUTINE obs_pre_prof( profdata, prodatqc, ld_var, &
[6140]	254	& kpi, kpj, kpk, &
[11361]	255	& zmask, pglam, pgphi, &
	256	& ld_nea, ld_bound_reject, kdailyavtypes, kqc_cutoff )
[2128]	257
[6140]	258	!!----------------------------------------------------------------------
	259	!! * ROUTINE obs_pre_prof *
[2128]	260	!!
[6140]	261	!! ** Purpose : First level check and screening of profiles
[2128]	262	!!
[6140]	263	!! ** Method : First level check and screening of profiles
[2128]	264	!!
	265	!! History :
	266	!! ! 2007-06 (K. Mogensen) original : T and S profile data
	267	!! ! 2008-09 (M. Valdivieso) : TAO velocity data
[6140]	268	!! ! 2009-01 (K. Mogensen) : New feedback stricture
	269	!! ! 2015-02 (M. Martin) : Combined profile routine.
[2128]	270	!!
	271	!!----------------------------------------------------------------------
[11361]	272	!! * Modules used
	273	USE par_oce ! Ocean parameters
	274	USE dom_oce, ONLY : & ! Geographical information
	275	& gdept_1d, &
	276	& nproc
[6140]	277
[2128]	278	!! * Arguments
	279	TYPE(obs_prof), INTENT(INOUT) :: profdata ! Full set of profile data
	280	TYPE(obs_prof), INTENT(INOUT) :: prodatqc ! Subset of profile data not failing screening
[11361]	281	LOGICAL, DIMENSION(profdata%nvar), INTENT(IN) :: &
	282	& ld_var ! Observed variables switches
[6140]	283	LOGICAL, INTENT(IN) :: ld_nea ! Switch for rejecting observation near land
[11361]	284	LOGICAL, INTENT(IN) :: ld_bound_reject ! Switch for rejecting observations near the boundary
[6140]	285	INTEGER, INTENT(IN) :: kpi, kpj, kpk ! Local domain sizes
	286	INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &
	287	& kdailyavtypes ! Types for daily averages
[11361]	288	REAL(wp), INTENT(IN), DIMENSION(kpi,kpj,kpk,profdata%nvar) :: &
	289	& zmask
	290	REAL(wp), INTENT(IN), DIMENSION(kpi,kpj,profdata%nvar) :: &
	291	& pglam, &
	292	& pgphi
	293	INTEGER, INTENT(IN), OPTIONAL :: kqc_cutoff ! cut off for QC value
[6140]	294
[2128]	295	!! * Local declarations
[11361]	296	INTEGER :: iqc_cutoff = 255 ! cut off for QC value
[2128]	297	INTEGER :: iyea0 ! Initial date
	298	INTEGER :: imon0 ! - (year, month, day, hour, minute)
	299	INTEGER :: iday0
	300	INTEGER :: ihou0
	301	INTEGER :: imin0
	302	INTEGER :: icycle ! Current assimilation cycle
[11361]	303	! Counters for observations that are
	304	INTEGER :: iotdobs ! - outside time domain
	305	INTEGER, DIMENSION(profdata%nvar) :: iosdvobs ! - outside space domain
	306	INTEGER, DIMENSION(profdata%nvar) :: ilanvobs ! - within a model land cell
	307	INTEGER, DIMENSION(profdata%nvar) :: inlavobs ! - close to land
	308	INTEGER, DIMENSION(profdata%nvar) :: ibdyvobs ! - boundary
	309	INTEGER :: igrdobs ! - fail the grid search
	310	INTEGER :: iuvchku ! - reject UVEL if VVEL rejected
	311	INTEGER :: iuvchkv ! - reject VVEL if UVEL rejected
	312	! Global counters for observations that are
	313	INTEGER :: iotdobsmpp ! - outside time domain
	314	INTEGER, DIMENSION(profdata%nvar) :: iosdvobsmpp ! - outside space domain
	315	INTEGER, DIMENSION(profdata%nvar) :: ilanvobsmpp ! - within a model land cell
	316	INTEGER, DIMENSION(profdata%nvar) :: inlavobsmpp ! - close to land
	317	INTEGER, DIMENSION(profdata%nvar) :: ibdyvobsmpp ! - boundary
	318	INTEGER :: igrdobsmpp ! - fail the grid search
	319	INTEGER :: iuvchkumpp ! - reject UVEL if VVEL rejected
	320	INTEGER :: iuvchkvmpp ! - reject VVEL if UVEL rejected
[2128]	321	TYPE(obs_prof_valid) :: llvalid ! Profile selection
	322	TYPE(obs_prof_valid), DIMENSION(profdata%nvar) :: &
[11361]	323	& llvvalid ! vars selection
[2128]	324	INTEGER :: jvar ! Variable loop variable
	325	INTEGER :: jobs ! Obs. loop variable
	326	INTEGER :: jstp ! Time loop variable
	327	INTEGER :: inrc ! Time index variable
[11361]	328	CHARACTER(LEN=256) :: cout1 ! Diagnostic output line
	329	CHARACTER(LEN=256) :: cout2 ! Diagnostic output line
[2128]	330
[6140]	331	IF(lwp) WRITE(numout,*)'obs_pre_prof: Preparing the profile data...'
	332	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
[2128]	333
	334	! Initial date initialization (year, month, day, hour, minute)
	335	iyea0 = ndate0 / 10000
	336	imon0 = ( ndate0 - iyea0 * 10000 ) / 100
	337	iday0 = ndate0 - iyea0 * 10000 - imon0 * 100
[11361]	338	ihou0 = 0
	339	imin0 = 0
[2128]	340
	341	icycle = no ! Assimilation cycle
	342
[11361]	343	! Diagnostics counters for various failures.
[2128]	344
[11361]	345	iotdobs = 0
	346	igrdobs = 0
	347	iosdvobs(:) = 0
	348	ilanvobs(:) = 0
	349	inlavobs(:) = 0
	350	ibdyvobs(:) = 0
	351	iuvchku = 0
	352	iuvchkv = 0
[2128]	353
[11361]	354
	355	! Set QC cutoff to optional value if provided
	356	IF ( PRESENT(kqc_cutoff) ) iqc_cutoff=kqc_cutoff
	357
[2128]	358	! -----------------------------------------------------------------------
	359	! Find time coordinate for profiles
	360	! -----------------------------------------------------------------------
	361
[6140]	362	IF ( PRESENT(kdailyavtypes) ) THEN
	363	CALL obs_coo_tim_prof( icycle, &
	364	& iyea0, imon0, iday0, ihou0, imin0, &
	365	& profdata%nprof, profdata%nyea, profdata%nmon, &
	366	& profdata%nday, profdata%nhou, profdata%nmin, &
	367	& profdata%ntyp, profdata%nqc, profdata%mstp, &
[11361]	368	& iotdobs, kdailyavtypes = kdailyavtypes, &
	369	& kqc_cutoff = iqc_cutoff )
[6140]	370	ELSE
	371	CALL obs_coo_tim_prof( icycle, &
	372	& iyea0, imon0, iday0, ihou0, imin0, &
	373	& profdata%nprof, profdata%nyea, profdata%nmon, &
	374	& profdata%nday, profdata%nhou, profdata%nmin, &
	375	& profdata%ntyp, profdata%nqc, profdata%mstp, &
[11361]	376	& iotdobs, kqc_cutoff = iqc_cutoff )
[6140]	377	ENDIF
	378
[2128]	379	CALL obs_mpp_sum_integer( iotdobs, iotdobsmpp )
	380
	381	! -----------------------------------------------------------------------
	382	! Check for profiles failing the grid search
	383	! -----------------------------------------------------------------------
	384
[11361]	385	DO jvar = 1, profdata%nvar
	386	CALL obs_coo_grd( profdata%nprof, profdata%mi(:,jvar), profdata%mj(:,jvar), &
	387	& profdata%nqc, igrdobs )
	388	END DO
[2128]	389
	390	CALL obs_mpp_sum_integer( igrdobs, igrdobsmpp )
	391
	392	! -----------------------------------------------------------------------
[11361]	393	! Reject all observations for profiles with nqc > iqc_cutoff
[2128]	394	! -----------------------------------------------------------------------
	395
[11361]	396	CALL obs_pro_rej( profdata, kqc_cutoff = iqc_cutoff )
[2128]	397
	398	! -----------------------------------------------------------------------
	399	! Check for land points. This includes points below the model
	400	! bathymetry so this is done for every point in the profile
	401	! -----------------------------------------------------------------------
	402
[11361]	403	DO jvar = 1, profdata%nvar
	404	CALL obs_coo_spc_3d( profdata%nprof, profdata%nvprot(jvar), &
	405	& profdata%npvsta(:,jvar), profdata%npvend(:,jvar), &
	406	& jpi, jpj, &
	407	& jpk, &
	408	& profdata%mi, profdata%mj, &
	409	& profdata%var(jvar)%mvk, &
	410	& profdata%rlam, profdata%rphi, &
	411	& profdata%var(jvar)%vdep, &
	412	& pglam(:,:,jvar), pgphi(:,:,jvar), &
	413	& gdept_1d, zmask(:,:,:,jvar), &
	414	& profdata%nqc, profdata%var(jvar)%nvqc, &
	415	& iosdvobs(jvar), ilanvobs(jvar), &
	416	& inlavobs(jvar), ld_nea, &
	417	& ibdyvobs(jvar), ld_bound_reject, &
	418	& iqc_cutoff )
[2128]	419
[11361]	420	CALL obs_mpp_sum_integer( iosdvobs(jvar), iosdvobsmpp(jvar) )
	421	CALL obs_mpp_sum_integer( ilanvobs(jvar), ilanvobsmpp(jvar) )
	422	CALL obs_mpp_sum_integer( inlavobs(jvar), inlavobsmpp(jvar) )
	423	CALL obs_mpp_sum_integer( ibdyvobs(jvar), ibdyvobsmpp(jvar) )
	424	END DO
[2128]	425
	426	! -----------------------------------------------------------------------
	427	! Reject u if v is rejected and vice versa
	428	! -----------------------------------------------------------------------
	429
[6140]	430	IF ( TRIM(profdata%cvars(1)) == 'UVEL' ) THEN
[11361]	431	CALL obs_uv_rej( profdata, iuvchku, iuvchkv, iqc_cutoff )
[6140]	432	CALL obs_mpp_sum_integer( iuvchku, iuvchkumpp )
	433	CALL obs_mpp_sum_integer( iuvchkv, iuvchkvmpp )
	434	ENDIF
[2128]	435
	436	! -----------------------------------------------------------------------
	437	! Copy useful data from the profdata data structure to
	438	! the prodatqc data structure
	439	! -----------------------------------------------------------------------
	440
	441	! Allocate the selection arrays
	442
	443	ALLOCATE( llvalid%luse(profdata%nprof) )
	444	DO jvar = 1,profdata%nvar
	445	ALLOCATE( llvvalid(jvar)%luse(profdata%nvprot(jvar)) )
	446	END DO
	447
[11361]	448	! We want all data which has qc flags <= iqc_cutoff
[2128]	449
[11361]	450	llvalid%luse(:) = ( profdata%nqc(:) <= iqc_cutoff )
[2128]	451	DO jvar = 1,profdata%nvar
[11361]	452	llvvalid(jvar)%luse(:) = ( profdata%var(jvar)%nvqc(:) <= iqc_cutoff )
[2128]	453	END DO
	454
	455	! The actual copying
	456
	457	CALL obs_prof_compress( profdata, prodatqc, .TRUE., numout, &
	458	& lvalid=llvalid, lvvalid=llvvalid )
	459
	460	! Dellocate the selection arrays
	461	DEALLOCATE( llvalid%luse )
	462	DO jvar = 1,profdata%nvar
	463	DEALLOCATE( llvvalid(jvar)%luse )
	464	END DO
	465
	466	! -----------------------------------------------------------------------
	467	! Print information about what observations are left after qc
	468	! -----------------------------------------------------------------------
	469
	470	! Update the total observation counter array
	471
	472	IF(lwp) THEN
[6140]	473
[2128]	474	WRITE(numout,*)
[11361]	475	WRITE(numout,*) ' Profiles outside time domain = ', &
[2128]	476	& iotdobsmpp
[11361]	477	WRITE(numout,*) ' Remaining profiles that failed grid search = ', &
[2128]	478	& igrdobsmpp
[11361]	479	DO jvar = 1, profdata%nvar
	480	WRITE(numout,*) ' Remaining '//prodatqc%cvars(jvar)//' data outside space domain = ', &
	481	& iosdvobsmpp(jvar)
	482	WRITE(numout,*) ' Remaining '//prodatqc%cvars(jvar)//' data at land points = ', &
	483	& ilanvobsmpp(jvar)
	484	IF (ld_nea) THEN
	485	WRITE(numout,*) ' Remaining '//prodatqc%cvars(jvar)//' data near land points (removed) = ',&
	486	& inlavobsmpp(jvar)
	487	ELSE
	488	WRITE(numout,*) ' Remaining '//prodatqc%cvars(jvar)//' data near land points (kept) = ',&
	489	& inlavobsmpp(jvar)
	490	ENDIF
	491	IF ( TRIM(profdata%cvars(jvar)) == 'UVEL' ) THEN
	492	WRITE(numout,*) ' U observation rejected since V rejected = ', &
	493	& iuvchku
	494	ELSE IF ( TRIM(profdata%cvars(jvar)) == 'VVEL' ) THEN
	495	WRITE(numout,*) ' V observation rejected since U rejected = ', &
	496	& iuvchkv
	497	ENDIF
	498	WRITE(numout,*) ' Remaining '//prodatqc%cvars(jvar)//' data near open boundary (removed) = ',&
	499	& ibdyvobsmpp(jvar)
	500	WRITE(numout,*) ' '//prodatqc%cvars(jvar)//' data accepted = ', &
	501	& prodatqc%nvprotmpp(jvar)
	502	END DO
[2128]	503
	504	WRITE(numout,*)
	505	WRITE(numout,*) ' Number of observations per time step :'
	506	WRITE(numout,*)
[11361]	507	WRITE(cout1,'(10X,A9,5X,A8)') 'Time step', 'Profiles'
	508	WRITE(cout2,'(10X,A9,5X,A8)') '---------', '--------'
	509	DO jvar = 1, prodatqc%nvar
	510	WRITE(cout1,'(A,5X,A11)') TRIM(cout1), TRIM(prodatqc%cvars(jvar))
	511	WRITE(cout2,'(A,5X,A11)') TRIM(cout2), '-----------'
	512	END DO
	513	WRITE(numout,*) cout1
	514	WRITE(numout,*) cout2
[2128]	515	ENDIF
	516
	517	DO jobs = 1, prodatqc%nprof
	518	inrc = prodatqc%mstp(jobs) + 2 - nit000
	519	prodatqc%npstp(inrc) = prodatqc%npstp(inrc) + 1
	520	DO jvar = 1, prodatqc%nvar
	521	IF ( prodatqc%npvend(jobs,jvar) > 0 ) THEN
	522	prodatqc%nvstp(inrc,jvar) = prodatqc%nvstp(inrc,jvar) + &
	523	& ( prodatqc%npvend(jobs,jvar) - &
	524	& prodatqc%npvsta(jobs,jvar) + 1 )
	525	ENDIF
	526	END DO
	527	END DO
	528
	529
	530	CALL obs_mpp_sum_integers( prodatqc%npstp, prodatqc%npstpmpp, &
	531	& nitend - nit000 + 2 )
	532	DO jvar = 1, prodatqc%nvar
	533	CALL obs_mpp_sum_integers( prodatqc%nvstp(:,jvar), &
	534	& prodatqc%nvstpmpp(:,jvar), &
	535	& nitend - nit000 + 2 )
	536	END DO
	537
	538	IF ( lwp ) THEN
	539	DO jstp = nit000 - 1, nitend
	540	inrc = jstp - nit000 + 2
[11361]	541	WRITE(cout1,'(10X,I9,5X,I8)') jstp, prodatqc%npstpmpp(inrc)
	542	DO jvar = 1, prodatqc%nvar
	543	WRITE(cout1,'(A,5X,I11)') TRIM(cout1), prodatqc%nvstpmpp(inrc,jvar)
	544	END DO
	545	WRITE(numout,*) cout1
[2128]	546	END DO
	547	ENDIF
	548
[6140]	549	END SUBROUTINE obs_pre_prof
[2128]	550
	551	SUBROUTINE obs_coo_tim( kcycle, &
	552	& kyea0, kmon0, kday0, khou0, kmin0, &
	553	& kobsno, &
	554	& kobsyea, kobsmon, kobsday, kobshou, kobsmin, &
	555	& kobsqc, kobsstp, kotdobs )
	556	!!----------------------------------------------------------------------
	557	!! * ROUTINE obs_coo_tim *
	558	!!
	559	!! ** Purpose : Compute the number of time steps to the observation time.
	560	!!
	561	!! ** Method : For time coordinates ( yea_obs, mon_obs, day_obs,
	562	!! hou_obs, min_obs ), this routine locates the time step
	563	!! that is closest to this time.
	564	!!
	565	!! ** Action :
	566	!!
	567	!! References :
	568	!!
	569	!! History :
	570	!! ! 1997-07 (A. Weaver) Original
	571	!! ! 2006-08 (A. Weaver) NEMOVAR migration
	572	!! ! 2006-10 (A. Weaver) Cleanup
	573	!! ! 2007-01 (K. Mogensen) Rewritten with loop
	574	!! ! 2010-05 (D. Lea) Fix in leap year calculation for NEMO vn3.2
	575	!!----------------------------------------------------------------------
	576	!! * Modules used
	577	USE dom_oce, ONLY : & ! Geographical information
	578	& rdt
	579	USE phycst, ONLY : & ! Physical constants
	580	& rday, &
	581	& rmmss, &
	582	& rhhmm
	583	!! * Arguments
	584	INTEGER, INTENT(IN) :: kcycle ! Current cycle
	585	INTEGER, INTENT(IN) :: kyea0 ! Initial date coordinates
	586	INTEGER, INTENT(IN) :: kmon0
	587	INTEGER, INTENT(IN) :: kday0
	588	INTEGER, INTENT(IN) :: khou0
	589	INTEGER, INTENT(IN) :: kmin0
	590	INTEGER, INTENT(IN) :: kobsno ! Number of observations
	591	INTEGER, INTENT(INOUT) :: kotdobs ! Number of observations failing time check
	592	INTEGER, DIMENSION(kobsno), INTENT(IN ) :: &
	593	& kobsyea, & ! Observation time coordinates
	594	& kobsmon, &
	595	& kobsday, &
	596	& kobshou, &
	597	& kobsmin
	598	INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: &
	599	& kobsqc ! Quality control flag
	600	INTEGER, DIMENSION(kobsno), INTENT(OUT) :: &
	601	& kobsstp ! Number of time steps up to the
	602	! observation time
	603
	604	!! * Local declarations
	605	INTEGER :: jyea
	606	INTEGER :: jmon
	607	INTEGER :: jday
	608	INTEGER :: jobs
	609	INTEGER :: iyeastr
	610	INTEGER :: iyeaend
	611	INTEGER :: imonstr
	612	INTEGER :: imonend
	613	INTEGER :: idaystr
	614	INTEGER :: idayend
	615	INTEGER :: iskip
	616	INTEGER :: idaystp
	617	REAL(KIND=wp) :: zminstp
	618	REAL(KIND=wp) :: zhoustp
	619	REAL(KIND=wp) :: zobsstp
	620	INTEGER, DIMENSION(12) :: imonth_len !: length in days of the months of the current year
	621
	622	!-----------------------------------------------------------------------
	623	! Initialization
	624	!-----------------------------------------------------------------------
	625
	626	! Intialize the number of time steps per day
	627	idaystp = NINT( rday / rdt )
	628
	629	!---------------------------------------------------------------------
	630	! Locate the model time coordinates for interpolation
	631	!---------------------------------------------------------------------
	632
	633	DO jobs = 1, kobsno
	634
	635	! Initialize the time step counter
	636	kobsstp(jobs) = nit000 - 1
	637
	638	! Flag if observation date is less than the initial date
	639
	640	IF ( ( kobsyea(jobs) < kyea0 ) &
	641	& .OR. ( ( kobsyea(jobs) == kyea0 ) &
	642	& .AND. ( kobsmon(jobs) < kmon0 ) ) &
	643	& .OR. ( ( kobsyea(jobs) == kyea0 ) &
	644	& .AND. ( kobsmon(jobs) == kmon0 ) &
	645	& .AND. ( kobsday(jobs) < kday0 ) ) &
	646	& .OR. ( ( kobsyea(jobs) == kyea0 ) &
	647	& .AND. ( kobsmon(jobs) == kmon0 ) &
	648	& .AND. ( kobsday(jobs) == kday0 ) &
	649	& .AND. ( kobshou(jobs) < khou0 ) ) &
	650	& .OR. ( ( kobsyea(jobs) == kyea0 ) &
	651	& .AND. ( kobsmon(jobs) == kmon0 ) &
	652	& .AND. ( kobsday(jobs) == kday0 ) &
	653	& .AND. ( kobshou(jobs) == khou0 ) &
	654	& .AND. ( kobsmin(jobs) <= kmin0 ) ) ) THEN
	655	kobsstp(jobs) = -1
[11361]	656	kobsqc(jobs) = IBSET(kobsqc(jobs),13)
[2128]	657	kotdobs = kotdobs + 1
	658	CYCLE
	659	ENDIF
	660
	661	! Compute the number of time steps to the observation day
	662	iyeastr = kyea0
	663	iyeaend = kobsyea(jobs)
	664
	665	!---------------------------------------------------------------------
	666	! Year loop
	667	!---------------------------------------------------------------------
	668	DO jyea = iyeastr, iyeaend
	669
	670	CALL calc_month_len( jyea, imonth_len )
	671
	672	imonstr = 1
	673	IF ( jyea == kyea0 ) imonstr = kmon0
	674	imonend = 12
	675	IF ( jyea == kobsyea(jobs) ) imonend = kobsmon(jobs)
	676
	677	! Month loop
	678	DO jmon = imonstr, imonend
	679
	680	idaystr = 1
	681	IF ( ( jmon == kmon0 ) &
	682	& .AND. ( jyea == kyea0 ) ) idaystr = kday0
	683	idayend = imonth_len(jmon)
	684	IF ( ( jmon == kobsmon(jobs) ) &
	685	& .AND. ( jyea == kobsyea(jobs) ) ) idayend = kobsday(jobs) - 1
	686
	687	! Day loop
	688	DO jday = idaystr, idayend
	689	kobsstp(jobs) = kobsstp(jobs) + idaystp
	690	END DO
	691
	692	END DO
	693
	694	END DO
	695
	696	! Add in the number of time steps to the observation minute
	697	zminstp = rmmss / rdt
	698	zhoustp = rhhmm * zminstp
	699
	700	zobsstp = REAL( kobsmin(jobs) - kmin0, KIND=wp ) * zminstp &
	701	& + REAL( kobshou(jobs) - khou0, KIND=wp ) * zhoustp
	702	kobsstp(jobs) = kobsstp(jobs) + NINT( zobsstp )
	703
	704	! Flag if observation step outside the time window
	705	IF ( ( kobsstp(jobs) < ( nit000 - 1 ) ) &
	706	& .OR.( kobsstp(jobs) > nitend ) ) THEN
[11361]	707	kobsqc(jobs) = IBSET(kobsqc(jobs),13)
[2128]	708	kotdobs = kotdobs + 1
	709	CYCLE
	710	ENDIF
	711
	712	END DO
	713
	714	END SUBROUTINE obs_coo_tim
	715
	716	SUBROUTINE calc_month_len( iyear, imonth_len )
	717	!!----------------------------------------------------------------------
	718	!! * ROUTINE calc_month_len *
	719	!!
	720	!! ** Purpose : Compute the number of days in a months given a year.
	721	!!
	722	!! ** Method :
	723	!!
	724	!! ** Action :
	725	!!
	726	!! History :
	727	!! ! 10-05 (D. Lea) New routine based on day_init
	728	!!----------------------------------------------------------------------
	729
	730	INTEGER, DIMENSION(12) :: imonth_len !: length in days of the months of the current year
	731	INTEGER :: iyear !: year
	732
	733	! length of the month of the current year (from nleapy, read in namelist)
	734	IF ( nleapy < 2 ) THEN
	735	imonth_len(:) = (/ 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 /)
	736	IF ( nleapy == 1 ) THEN ! we are using calendar with leap years
	737	IF ( MOD(iyear, 4) == 0 .AND. ( MOD(iyear, 400) == 0 .OR. MOD(iyear, 100) /= 0 ) ) THEN
	738	imonth_len(2) = 29
	739	ENDIF
	740	ENDIF
	741	ELSE
	742	imonth_len(:) = nleapy ! all months with nleapy days per year
	743	ENDIF
	744
	745	END SUBROUTINE
	746
	747	SUBROUTINE obs_coo_tim_prof( kcycle, &
	748	& kyea0, kmon0, kday0, khou0, kmin0, &
	749	& kobsno, &
	750	& kobsyea, kobsmon, kobsday, kobshou, kobsmin, &
[11361]	751	& ktyp, kobsqc, kobsstp, kotdobs, kdailyavtypes, &
	752	& kqc_cutoff )
[2128]	753	!!----------------------------------------------------------------------
	754	!! * ROUTINE obs_coo_tim *
	755	!!
	756	!! ** Purpose : Compute the number of time steps to the observation time.
	757	!!
	758	!! ** Method : For time coordinates ( yea_obs, mon_obs, day_obs,
	759	!! hou_obs, min_obs ), this routine locates the time step
	760	!! that is closest to this time.
	761	!!
	762	!! ** Action :
	763	!!
	764	!! References :
	765	!!
	766	!! History :
	767	!! ! 1997-07 (A. Weaver) Original
	768	!! ! 2006-08 (A. Weaver) NEMOVAR migration
	769	!! ! 2006-10 (A. Weaver) Cleanup
	770	!! ! 2007-01 (K. Mogensen) Rewritten with loop
	771	!!----------------------------------------------------------------------
	772	!! * Modules used
	773	!! * Arguments
	774	INTEGER, INTENT(IN) :: kcycle ! Current cycle
	775	INTEGER, INTENT(IN) :: kyea0 ! Initial date coordinates
	776	INTEGER, INTENT(IN) :: kmon0
	777	INTEGER, INTENT(IN) :: kday0
	778	INTEGER, INTENT(IN) :: khou0
	779	INTEGER, INTENT(IN) :: kmin0
	780	INTEGER, INTENT(IN) :: kobsno ! Number of observations
	781	INTEGER, INTENT(INOUT) :: kotdobs ! Number of observations failing time check
	782	INTEGER, DIMENSION(kobsno), INTENT(IN ) :: &
	783	& kobsyea, & ! Observation time coordinates
	784	& kobsmon, &
	785	& kobsday, &
	786	& kobshou, &
	787	& kobsmin, &
	788	& ktyp ! Observation type.
	789	INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: &
	790	& kobsqc ! Quality control flag
	791	INTEGER, DIMENSION(kobsno), INTENT(OUT) :: &
	792	& kobsstp ! Number of time steps up to the
	793	! observation time
	794	INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &
	795	& kdailyavtypes ! Types for daily averages
[11361]	796	INTEGER, OPTIONAL, INTENT(IN) :: kqc_cutoff ! QC cutoff value
	797
[2128]	798	!! * Local declarations
	799	INTEGER :: jobs
[11361]	800	INTEGER :: iqc_cutoff=255
[2128]	801
	802	!-----------------------------------------------------------------------
	803	! Call standard obs_coo_tim
	804	!-----------------------------------------------------------------------
	805
	806	CALL obs_coo_tim( kcycle, &
	807	& kyea0, kmon0, kday0, khou0, kmin0, &
	808	& kobsno, &
	809	& kobsyea, kobsmon, kobsday, kobshou, kobsmin, &
	810	& kobsqc, kobsstp, kotdobs )
	811
	812	!------------------------------------------------------------------------
	813	! Always reject daily averaged data (e.g. MRB data (820)) at initial time
	814	!------------------------------------------------------------------------
	815
	816	IF ( PRESENT(kdailyavtypes) ) THEN
	817	DO jobs = 1, kobsno
	818
[11361]	819	IF ( kobsqc(jobs) <= iqc_cutoff ) THEN
[2128]	820
	821	IF ( ( kobsstp(jobs) == (nit000 - 1) ).AND.&
	822	& ( ANY (kdailyavtypes(:) == ktyp(jobs)) ) ) THEN
[11361]	823	kobsqc(jobs) = IBSET(kobsqc(jobs),13)
[2128]	824	kotdobs = kotdobs + 1
	825	CYCLE
	826	ENDIF
	827
	828	ENDIF
	829	END DO
	830	ENDIF
	831
	832
	833	END SUBROUTINE obs_coo_tim_prof
	834
	835	SUBROUTINE obs_coo_grd( kobsno, kobsi, kobsj, kobsqc, kgrdobs )
	836	!!----------------------------------------------------------------------
	837	!! * ROUTINE obs_coo_grd *
	838	!!
	839	!! ** Purpose : Verify that the grid search has not failed
	840	!!
	841	!! ** Method : The previously computed i,j indeces are checked
	842	!!
	843	!! ** Action :
	844	!!
	845	!! References :
	846	!!
	847	!! History :
	848	!! ! 2007-01 (K. Mogensen) Original
	849	!!----------------------------------------------------------------------
	850	!! * Arguments
	851	INTEGER, INTENT(IN) :: kobsno ! Number of observations
	852	INTEGER, DIMENSION(kobsno), INTENT(IN ) :: &
	853	& kobsi, & ! i,j indeces previously computed
	854	& kobsj
	855	INTEGER, INTENT(INOUT) :: kgrdobs ! Number of observations failing the check
	856	INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: &
	857	& kobsqc ! Quality control flag
	858
	859	!! * Local declarations
	860	INTEGER :: jobs ! Loop variable
	861
	862	! Flag if the grid search failed
	863
	864	DO jobs = 1, kobsno
	865	IF ( ( kobsi(jobs) <= 0 ) .AND. ( kobsj(jobs) <= 0 ) ) THEN
[11361]	866	kobsqc(jobs) = IBSET(kobsqc(jobs),12)
[2128]	867	kgrdobs = kgrdobs + 1
	868	ENDIF
	869	END DO
	870
	871	END SUBROUTINE obs_coo_grd
	872
	873	SUBROUTINE obs_coo_spc_2d( kobsno, kpi, kpj, &
	874	& kobsi, kobsj, pobslam, pobsphi, &
	875	& plam, pphi, pmask, &
	876	& kobsqc, kosdobs, klanobs, &
[11361]	877	& knlaobs,ld_nea, &
	878	& kbdyobs,ld_bound_reject, &
	879	& kqc_cutoff )
[2128]	880	!!----------------------------------------------------------------------
	881	!! * ROUTINE obs_coo_spc_2d *
	882	!!
	883	!! ** Purpose : Check for points outside the domain and land points
	884	!!
	885	!! ** Method : Remove the observations that are outside the model space
	886	!! and time domain or located within model land cells.
	887	!!
	888	!! ** Action :
	889	!!
	890	!! History :
	891	!! ! 2007-03 (A. Weaver, K. Mogensen) Original
	892	!! ! 2007-06 (K. Mogensen et al) Reject obs. near land.
	893	!!----------------------------------------------------------------------
	894	!! * Modules used
	895
	896	!! * Arguments
	897	INTEGER, INTENT(IN) :: kobsno ! Total number of observations
	898	INTEGER, INTENT(IN) :: kpi ! Number of grid points in (i,j)
	899	INTEGER, INTENT(IN) :: kpj
	900	INTEGER, DIMENSION(kobsno), INTENT(IN) :: &
	901	& kobsi, & ! Observation (i,j) coordinates
	902	& kobsj
	903	REAL(KIND=wp), DIMENSION(kobsno), INTENT(IN) :: &
	904	& pobslam, & ! Observation (lon,lat) coordinates
	905	& pobsphi
	906	REAL(KIND=wp), DIMENSION(kpi,kpj), INTENT(IN) :: &
	907	& plam, pphi ! Model (lon,lat) coordinates
	908	REAL(KIND=wp), DIMENSION(kpi,kpj), INTENT(IN) :: &
	909	& pmask ! Land mask array
	910	INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: &
	911	& kobsqc ! Observation quality control
[11361]	912	INTEGER, INTENT(INOUT) :: kosdobs ! Observations outside space domain
	913	INTEGER, INTENT(INOUT) :: klanobs ! Observations within a model land cell
	914	INTEGER, INTENT(INOUT) :: knlaobs ! Observations near land
	915	INTEGER, INTENT(INOUT) :: kbdyobs ! Observations near boundary
	916	LOGICAL, INTENT(IN) :: ld_nea ! Flag observations near land
	917	LOGICAL, INTENT(IN) :: ld_bound_reject ! Flag observations near open boundary
	918	INTEGER, INTENT(IN) :: kqc_cutoff ! Cutoff QC value
	919
[2128]	920	!! * Local declarations
	921	REAL(KIND=wp), DIMENSION(2,2,kobsno) :: &
	922	& zgmsk ! Grid mask
[11361]	923	#if defined key_bdy
[2128]	924	REAL(KIND=wp), DIMENSION(2,2,kobsno) :: &
[11361]	925	& zbmsk ! Boundary mask
	926	REAL(KIND=wp), DIMENSION(jpi,jpj) :: zbdymask
	927	#endif
	928	REAL(KIND=wp), DIMENSION(2,2,kobsno) :: &
[2128]	929	& zglam, & ! Model longitude at grid points
	930	& zgphi ! Model latitude at grid points
	931	INTEGER, DIMENSION(2,2,kobsno) :: &
	932	& igrdi, & ! Grid i,j
	933	& igrdj
	934	LOGICAL :: lgridobs ! Is observation on a model grid point.
	935	INTEGER :: iig, ijg ! i,j of observation on model grid point.
	936	INTEGER :: jobs, ji, jj
	937
	938	! Get grid point indices
	939
	940	DO jobs = 1, kobsno
	941
	942	! For invalid points use 2,2
	943
[11361]	944	IF ( kobsqc(jobs) >= kqc_cutoff ) THEN
[2128]	945
	946	igrdi(1,1,jobs) = 1
	947	igrdj(1,1,jobs) = 1
	948	igrdi(1,2,jobs) = 1
	949	igrdj(1,2,jobs) = 2
	950	igrdi(2,1,jobs) = 2
	951	igrdj(2,1,jobs) = 1
	952	igrdi(2,2,jobs) = 2
	953	igrdj(2,2,jobs) = 2
	954
	955	ELSE
	956
	957	igrdi(1,1,jobs) = kobsi(jobs)-1
	958	igrdj(1,1,jobs) = kobsj(jobs)-1
	959	igrdi(1,2,jobs) = kobsi(jobs)-1
	960	igrdj(1,2,jobs) = kobsj(jobs)
	961	igrdi(2,1,jobs) = kobsi(jobs)
	962	igrdj(2,1,jobs) = kobsj(jobs)-1
	963	igrdi(2,2,jobs) = kobsi(jobs)
	964	igrdj(2,2,jobs) = kobsj(jobs)
	965
	966	ENDIF
	967
	968	END DO
[11361]	969
	970	#if defined key_bdy
	971	! Create a mask grid points in boundary rim
	972	IF (ld_bound_reject) THEN
	973	zbdymask(:,:) = 1.0_wp
	974	DO ji = 1, nb_bdy
	975	DO jj = 1, idx_bdy(ji)%nblen(1)
	976	zbdymask(idx_bdy(ji)%nbi(jj,1),idx_bdy(ji)%nbj(jj,1)) = 0.0_wp
	977	ENDDO
	978	ENDDO
	979
	980	CALL obs_int_comm_2d( 2, 2, kobsno, kpi, kpj, igrdi, igrdj, zbdymask, zbmsk )
	981	ENDIF
	982	#endif
[2128]	983
[6140]	984	CALL obs_int_comm_2d( 2, 2, kobsno, kpi, kpj, igrdi, igrdj, pmask, zgmsk )
	985	CALL obs_int_comm_2d( 2, 2, kobsno, kpi, kpj, igrdi, igrdj, plam, zglam )
	986	CALL obs_int_comm_2d( 2, 2, kobsno, kpi, kpj, igrdi, igrdj, pphi, zgphi )
[2128]	987
	988	DO jobs = 1, kobsno
	989
	990	! Skip bad observations
[11361]	991	IF ( kobsqc(jobs) >= kqc_cutoff ) CYCLE
[2128]	992
	993	! Flag if the observation falls outside the model spatial domain
	994	IF ( ( pobslam(jobs) < -180. ) &
	995	& .OR. ( pobslam(jobs) > 180. ) &
	996	& .OR. ( pobsphi(jobs) < -90. ) &
	997	& .OR. ( pobsphi(jobs) > 90. ) ) THEN
[11361]	998	kobsqc(jobs) = IBSET(kobsqc(jobs),11)
[2128]	999	kosdobs = kosdobs + 1
	1000	CYCLE
	1001	ENDIF
	1002
	1003	! Flag if the observation falls with a model land cell
	1004	IF ( SUM( zgmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN
[11361]	1005	kobsqc(jobs) = IBSET(kobsqc(jobs),10)
[2128]	1006	klanobs = klanobs + 1
	1007	CYCLE
	1008	ENDIF
	1009
	1010	! Check if this observation is on a grid point
	1011
	1012	lgridobs = .FALSE.
	1013	iig = -1
	1014	ijg = -1
	1015	DO jj = 1, 2
	1016	DO ji = 1, 2
	1017	IF ( ( ABS( zgphi(ji,jj,jobs) - pobsphi(jobs) ) < 1.0e-6_wp ) &
	1018	& .AND. &
[11361]	1019	& ( ABS( MOD( zglam(ji,jj,jobs) - pobslam(jobs),360.0) ) &
	1020	& < 1.0e-6_wp ) ) THEN
[2128]	1021	lgridobs = .TRUE.
	1022	iig = ji
	1023	ijg = jj
	1024	ENDIF
	1025	END DO
	1026	END DO
[11361]	1027
[2128]	1028	IF (lgridobs) THEN
	1029	IF (zgmsk(iig,ijg,jobs) == 0.0_wp ) THEN
[11361]	1030	kobsqc(jobs) = IBSET(kobsqc(jobs),10)
[2128]	1031	klanobs = klanobs + 1
	1032	CYCLE
	1033	ENDIF
	1034	ENDIF
[11361]	1035
	1036
[2128]	1037	! Flag if the observation falls is close to land
	1038	IF ( MINVAL( zgmsk(1:2,1:2,jobs) ) == 0.0_wp) THEN
	1039	knlaobs = knlaobs + 1
[11361]	1040	IF (ld_nea) THEN
	1041	kobsqc(jobs) = IBSET(kobsqc(jobs),9)
	1042	CYCLE
	1043	ENDIF
[2128]	1044	ENDIF
[11361]	1045
	1046	#if defined key_bdy
	1047	! Flag if the observation falls close to the boundary rim
	1048	IF (ld_bound_reject) THEN
	1049	IF ( MINVAL( zbmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN
	1050	kobsqc(jobs) = IBSET(kobsqc(jobs),8)
	1051	kbdyobs = kbdyobs + 1
	1052	CYCLE
	1053	ENDIF
	1054	! for observations on the grid...
	1055	IF (lgridobs) THEN
	1056	IF (zbmsk(iig,ijg,jobs) == 0.0_wp ) THEN
	1057	kobsqc(jobs) = IBSET(kobsqc(jobs),8)
	1058	kbdyobs = kbdyobs + 1
	1059	CYCLE
	1060	ENDIF
	1061	ENDIF
	1062	ENDIF
	1063	#endif
[2128]	1064
	1065	END DO
	1066
	1067	END SUBROUTINE obs_coo_spc_2d
	1068
	1069	SUBROUTINE obs_coo_spc_3d( kprofno, kobsno, kpstart, kpend, &
	1070	& kpi, kpj, kpk, &
	1071	& kobsi, kobsj, kobsk, &
	1072	& pobslam, pobsphi, pobsdep, &
	1073	& plam, pphi, pdep, pmask, &
	1074	& kpobsqc, kobsqc, kosdobs, &
[11361]	1075	& klanobs, knlaobs, ld_nea, &
	1076	& kbdyobs, ld_bound_reject, &
	1077	& kqc_cutoff )
[2128]	1078	!!----------------------------------------------------------------------
	1079	!! * ROUTINE obs_coo_spc_3d *
	1080	!!
	1081	!! ** Purpose : Check for points outside the domain and land points
	1082	!! Reset depth of observation above highest model level
	1083	!! to the value of highest model level
	1084	!!
	1085	!! ** Method : Remove the observations that are outside the model space
	1086	!! and time domain or located within model land cells.
	1087	!!
	1088	!! NB. T and S profile observations lying between the ocean
	1089	!! surface and the depth of the first model T point are
	1090	!! assigned a depth equal to that of the first model T pt.
	1091	!!
	1092	!! ** Action :
	1093	!!
	1094	!! History :
	1095	!! ! 2007-01 (K. Mogensen) Rewrite of parts of obs_scr
	1096	!! ! 2007-06 (K. Mogensen et al) Reject obs. near land.
	1097	!!----------------------------------------------------------------------
	1098	!! * Modules used
	1099	USE dom_oce, ONLY : & ! Geographical information
[6140]	1100	& gdepw_1d, &
	1101	& gdepw_0, &
[11361]	1102	& gdepw_n, &
	1103	#if defined key_vvl
[6140]	1104	& gdept_n, &
[11361]	1105	#endif
[6140]	1106	& ln_zco, &
[11361]	1107	& ln_zps, &
	1108	& ln_linssh
[2128]	1109
	1110	!! * Arguments
	1111	INTEGER, INTENT(IN) :: kprofno ! Number of profiles
	1112	INTEGER, INTENT(IN) :: kobsno ! Total number of observations
	1113	INTEGER, INTENT(IN) :: kpi ! Number of grid points in (i,j,k)
	1114	INTEGER, INTENT(IN) :: kpj
	1115	INTEGER, INTENT(IN) :: kpk
	1116	INTEGER, DIMENSION(kprofno), INTENT(IN) :: &
	1117	& kpstart, & ! Start of individual profiles
	1118	& kpend ! End of individual profiles
	1119	INTEGER, DIMENSION(kprofno), INTENT(IN) :: &
	1120	& kobsi, & ! Observation (i,j) coordinates
	1121	& kobsj
	1122	INTEGER, DIMENSION(kobsno), INTENT(IN) :: &
	1123	& kobsk ! Observation k coordinate
	1124	REAL(KIND=wp), DIMENSION(kprofno), INTENT(IN) :: &
	1125	& pobslam, & ! Observation (lon,lat) coordinates
	1126	& pobsphi
	1127	REAL(KIND=wp), DIMENSION(kobsno), INTENT(INOUT) :: &
	1128	& pobsdep ! Observation depths
	1129	REAL(KIND=wp), DIMENSION(kpi,kpj), INTENT(IN) :: &
	1130	& plam, pphi ! Model (lon,lat) coordinates
	1131	REAL(KIND=wp), DIMENSION(kpk), INTENT(IN) :: &
	1132	& pdep ! Model depth coordinates
	1133	REAL(KIND=wp), DIMENSION(kpi,kpj,kpk), INTENT(IN) :: &
	1134	& pmask ! Land mask array
	1135	INTEGER, DIMENSION(kprofno), INTENT(INOUT) :: &
	1136	& kpobsqc ! Profile quality control
	1137	INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: &
	1138	& kobsqc ! Observation quality control
	1139	INTEGER, INTENT(INOUT) :: kosdobs ! Observations outside space domain
	1140	INTEGER, INTENT(INOUT) :: klanobs ! Observations within a model land cell
	1141	INTEGER, INTENT(INOUT) :: knlaobs ! Observations near land
[11361]	1142	INTEGER, INTENT(INOUT) :: kbdyobs ! Observations near boundary
[2128]	1143	LOGICAL, INTENT(IN) :: ld_nea ! Flag observations near land
[11361]	1144	LOGICAL, INTENT(IN) :: ld_bound_reject ! Flag observations near open boundary
	1145	INTEGER, INTENT(IN) :: kqc_cutoff ! Cutoff QC value
	1146
[2128]	1147	!! * Local declarations
	1148	REAL(KIND=wp), DIMENSION(2,2,kpk,kprofno) :: &
	1149	& zgmsk ! Grid mask
[11361]	1150	#if defined key_bdy
	1151	REAL(KIND=wp), DIMENSION(2,2,kprofno) :: &
	1152	& zbmsk ! Boundary mask
	1153	REAL(KIND=wp), DIMENSION(jpi,jpj) :: zbdymask
	1154	#endif
[6140]	1155	REAL(KIND=wp), DIMENSION(2,2,kpk,kprofno) :: &
[11361]	1156	& zgdepw
[2128]	1157	REAL(KIND=wp), DIMENSION(2,2,kprofno) :: &
	1158	& zglam, & ! Model longitude at grid points
	1159	& zgphi ! Model latitude at grid points
	1160	INTEGER, DIMENSION(2,2,kprofno) :: &
	1161	& igrdi, & ! Grid i,j
	1162	& igrdj
	1163	LOGICAL :: lgridobs ! Is observation on a model grid point.
[6140]	1164	LOGICAL :: ll_next_to_land ! Is a profile next to land
[2128]	1165	INTEGER :: iig, ijg ! i,j of observation on model grid point.
	1166	INTEGER :: jobs, jobsp, jk, ji, jj
	1167
	1168	! Get grid point indices
	1169
	1170	DO jobs = 1, kprofno
	1171
	1172	! For invalid points use 2,2
	1173
[11361]	1174	IF ( kpobsqc(jobs) >= kqc_cutoff ) THEN
[2128]	1175
	1176	igrdi(1,1,jobs) = 1
	1177	igrdj(1,1,jobs) = 1
	1178	igrdi(1,2,jobs) = 1
	1179	igrdj(1,2,jobs) = 2
	1180	igrdi(2,1,jobs) = 2
	1181	igrdj(2,1,jobs) = 1
	1182	igrdi(2,2,jobs) = 2
	1183	igrdj(2,2,jobs) = 2
	1184
	1185	ELSE
	1186
	1187	igrdi(1,1,jobs) = kobsi(jobs)-1
	1188	igrdj(1,1,jobs) = kobsj(jobs)-1
	1189	igrdi(1,2,jobs) = kobsi(jobs)-1
	1190	igrdj(1,2,jobs) = kobsj(jobs)
	1191	igrdi(2,1,jobs) = kobsi(jobs)
	1192	igrdj(2,1,jobs) = kobsj(jobs)-1
	1193	igrdi(2,2,jobs) = kobsi(jobs)
	1194	igrdj(2,2,jobs) = kobsj(jobs)
	1195
	1196	ENDIF
	1197
	1198	END DO
[11361]	1199
	1200	#if defined key_bdy
	1201	! Create a mask grid points in boundary rim
	1202	IF (ld_bound_reject) THEN
	1203	zbdymask(:,:) = 1.0_wp
	1204	DO ji = 1, nb_bdy
	1205	DO jj = 1, idx_bdy(ji)%nblen(1)
	1206	zbdymask(idx_bdy(ji)%nbi(jj,1),idx_bdy(ji)%nbj(jj,1)) = 0.0_wp
	1207	ENDDO
	1208	ENDDO
	1209	ENDIF
	1210
	1211	CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, zbdymask, zbmsk )
	1212	#endif
[2128]	1213
[6140]	1214	CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, pmask, zgmsk )
	1215	CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, plam, zglam )
	1216	CALL obs_int_comm_2d( 2, 2, kprofno, kpi, kpj, igrdi, igrdj, pphi, zgphi )
[11361]	1217	! Need to know the bathy depth for each observation for sco
	1218	CALL obs_int_comm_3d( 2, 2, kprofno, kpi, kpj, kpk, igrdi, igrdj, gdepw_n(:,:,:), &
	1219	& zgdepw )
[2128]	1220
	1221	DO jobs = 1, kprofno
	1222
	1223	! Skip bad profiles
[11361]	1224	IF ( kpobsqc(jobs) >= kqc_cutoff ) CYCLE
[2128]	1225
	1226	! Check if this observation is on a grid point
	1227
	1228	lgridobs = .FALSE.
	1229	iig = -1
	1230	ijg = -1
	1231	DO jj = 1, 2
	1232	DO ji = 1, 2
	1233	IF ( ( ABS( zgphi(ji,jj,jobs) - pobsphi(jobs) ) < 1.0e-6_wp ) &
	1234	& .AND. &
[11361]	1235	& ( ABS( MOD( zglam(ji,jj,jobs) - pobslam(jobs),360.0) ) < 1.0e-6_wp ) &
[2128]	1236	& ) THEN
	1237	lgridobs = .TRUE.
	1238	iig = ji
	1239	ijg = jj
	1240	ENDIF
	1241	END DO
	1242	END DO
	1243
[11361]	1244	! Check if next to land
	1245	IF ( ANY( zgmsk(1:2,1:2,1,jobs) == 0.0_wp ) ) THEN
	1246	ll_next_to_land=.TRUE.
	1247	ELSE
	1248	ll_next_to_land=.FALSE.
[6140]	1249	ENDIF
[11361]	1250
[2128]	1251	! Reject observations
	1252
	1253	DO jobsp = kpstart(jobs), kpend(jobs)
	1254
	1255	! Flag if the observation falls outside the model spatial domain
	1256	IF ( ( pobslam(jobs) < -180. ) &
	1257	& .OR. ( pobslam(jobs) > 180. ) &
	1258	& .OR. ( pobsphi(jobs) < -90. ) &
	1259	& .OR. ( pobsphi(jobs) > 90. ) &
	1260	& .OR. ( pobsdep(jobsp) < 0.0 ) &
[4292]	1261	& .OR. ( pobsdep(jobsp) > gdepw_1d(kpk)) ) THEN
[11361]	1262	kobsqc(jobsp) = IBSET(kobsqc(jobsp),11)
[2128]	1263	kosdobs = kosdobs + 1
	1264	CYCLE
	1265	ENDIF
	1266
[11361]	1267	! To check if an observations falls within land there are two cases:
	1268	! 1: z-coordibnates, where the check uses the mask
	1269	! 2: terrain following (eg s-coordinates),
	1270	! where we use the depth of the bottom cell to mask observations
[6140]	1271
[11361]	1272	IF( (ln_linssh) .AND. ( ln_zps .OR. ln_zco ) ) THEN !(CASE 1)
	1273
	1274	! Flag if the observation falls with a model land cell
	1275	IF ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) &
	1276	& == 0.0_wp ) THEN
	1277	kobsqc(jobsp) = IBSET(kobsqc(jobsp),10)
	1278	klanobs = klanobs + 1
	1279	CYCLE
	1280	ENDIF
[6140]	1281
[11361]	1282	! Flag if the observation is close to land
	1283	IF ( MINVAL( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) == &
	1284	& 0.0_wp) THEN
	1285	knlaobs = knlaobs + 1
[6140]	1286	IF (ld_nea) THEN
[11361]	1287	kobsqc(jobsp) = IBSET(kobsqc(jobsp),10)
	1288	ENDIF
	1289	ENDIF
[6140]	1290
[11361]	1291	ELSE ! Case 2
	1292	! Flag if the observation is deeper than the bathymetry
	1293	! Or if it is within the mask
	1294	IF ( ANY( zgdepw(1:2,1:2,kpk,jobs) < pobsdep(jobsp) ) &
	1295	& .OR. &
[6140]	1296	& ( SUM( zgmsk(1:2,1:2,kobsk(jobsp)-1:kobsk(jobsp),jobs) ) &
	1297	& == 0.0_wp) ) THEN
[11361]	1298	kobsqc(jobsp) = IBSET(kobsqc(jobsp),10)
	1299	klanobs = klanobs + 1
	1300	CYCLE
	1301	ENDIF
	1302
	1303	! Flag if the observation is close to land
	1304	IF ( ll_next_to_land ) THEN
	1305	knlaobs = knlaobs + 1
[6140]	1306	IF (ld_nea) THEN
[11361]	1307	kobsqc(jobsp) = IBSET(kobsqc(jobsp),10)
	1308	ENDIF
	1309	ENDIF
	1310
[2128]	1311	ENDIF
[11361]	1312
[2128]	1313	! For observations on the grid reject them if their are at
	1314	! a masked point
	1315
	1316	IF (lgridobs) THEN
	1317	IF (zgmsk(iig,ijg,kobsk(jobsp)-1,jobs) == 0.0_wp ) THEN
[11361]	1318	kobsqc(jobsp) = IBSET(kobsqc(jobs),10)
[2128]	1319	klanobs = klanobs + 1
	1320	CYCLE
	1321	ENDIF
	1322	ENDIF
	1323
	1324	! Set observation depth equal to that of the first model depth
	1325	IF ( pobsdep(jobsp) <= pdep(1) ) THEN
	1326	pobsdep(jobsp) = pdep(1)
	1327	ENDIF
	1328
[11361]	1329	#if defined key_bdy
	1330	! Flag if the observation falls close to the boundary rim
	1331	IF (ld_bound_reject) THEN
	1332	IF ( MINVAL( zbmsk(1:2,1:2,jobs) ) == 0.0_wp ) THEN
	1333	kobsqc(jobsp) = IBSET(kobsqc(jobs),8)
	1334	kbdyobs = kbdyobs + 1
	1335	CYCLE
	1336	ENDIF
	1337	! for observations on the grid...
	1338	IF (lgridobs) THEN
	1339	IF (zbmsk(iig,ijg,jobs) == 0.0_wp ) THEN
	1340	kobsqc(jobsp) = IBSET(kobsqc(jobs),8)
	1341	kbdyobs = kbdyobs + 1
	1342	CYCLE
	1343	ENDIF
	1344	ENDIF
	1345	ENDIF
	1346	#endif
	1347
[2128]	1348	END DO
	1349	END DO
	1350
	1351	END SUBROUTINE obs_coo_spc_3d
	1352
[11361]	1353	SUBROUTINE obs_pro_rej( profdata, kqc_cutoff )
[2128]	1354	!!----------------------------------------------------------------------
	1355	!! * ROUTINE obs_pro_rej *
	1356	!!
	1357	!! ** Purpose : Reject all data within a rejected profile
	1358	!!
	1359	!! ** Method :
	1360	!!
	1361	!! ** Action :
	1362	!!
	1363	!! References :
	1364	!!
	1365	!! History :
	1366	!! ! 2007-10 (K. Mogensen) Original code
	1367	!!----------------------------------------------------------------------
	1368	!! * Modules used
	1369	!! * Arguments
	1370	TYPE(obs_prof), INTENT(INOUT) :: profdata ! Profile data
[11361]	1371	INTEGER, INTENT(IN) :: kqc_cutoff ! QC cutoff value
	1372
[2128]	1373	!! * Local declarations
	1374	INTEGER :: jprof
	1375	INTEGER :: jvar
	1376	INTEGER :: jobs
	1377
	1378	! Loop over profiles
	1379
	1380	DO jprof = 1, profdata%nprof
	1381
[11361]	1382	IF ( profdata%nqc(jprof) > kqc_cutoff ) THEN
[2128]	1383
	1384	DO jvar = 1, profdata%nvar
	1385
	1386	DO jobs = profdata%npvsta(jprof,jvar), &
	1387	& profdata%npvend(jprof,jvar)
	1388
	1389	profdata%var(jvar)%nvqc(jobs) = &
[11361]	1390	& IBSET(profdata%var(jvar)%nvqc(jobs),14)
[2128]	1391
	1392	END DO
	1393
	1394	END DO
	1395
	1396	ENDIF
	1397
	1398	END DO
	1399
	1400	END SUBROUTINE obs_pro_rej
	1401
[11361]	1402	SUBROUTINE obs_uv_rej( profdata, knumu, knumv, kqc_cutoff )
[2128]	1403	!!----------------------------------------------------------------------
	1404	!! * ROUTINE obs_uv_rej *
	1405	!!
	1406	!! ** Purpose : Reject u if v is rejected and vice versa
	1407	!!
	1408	!! ** Method :
	1409	!!
	1410	!! ** Action :
	1411	!!
	1412	!! References :
	1413	!!
	1414	!! History :
	1415	!! ! 2009-2 (K. Mogensen) Original code
	1416	!!----------------------------------------------------------------------
	1417	!! * Modules used
	1418	!! * Arguments
	1419	TYPE(obs_prof), INTENT(INOUT) :: profdata ! Profile data
	1420	INTEGER, INTENT(INOUT) :: knumu ! Number of u rejected
	1421	INTEGER, INTENT(INOUT) :: knumv ! Number of v rejected
[11361]	1422	INTEGER, INTENT(IN) :: kqc_cutoff ! QC cutoff value
	1423
[2128]	1424	!! * Local declarations
	1425	INTEGER :: jprof
	1426	INTEGER :: jvar
	1427	INTEGER :: jobs
	1428
	1429	! Loop over profiles
	1430
	1431	DO jprof = 1, profdata%nprof
	1432
	1433	IF ( ( profdata%npvsta(jprof,1) /= profdata%npvsta(jprof,2) ) .OR. &
	1434	& ( profdata%npvend(jprof,1) /= profdata%npvend(jprof,2) ) ) THEN
	1435
	1436	CALL ctl_stop('U,V profiles inconsistent in obs_uv_rej')
	1437	RETURN
	1438
	1439	ENDIF
	1440
	1441	DO jobs = profdata%npvsta(jprof,1), profdata%npvend(jprof,1)
	1442
[11361]	1443	IF ( ( profdata%var(1)%nvqc(jobs) > kqc_cutoff ) .AND. &
	1444	& ( profdata%var(2)%nvqc(jobs) <= kqc_cutoff) ) THEN
	1445	profdata%var(2)%nvqc(jobs) = IBSET(profdata%var(1)%nvqc(jobs),15)
[2128]	1446	knumv = knumv + 1
	1447	ENDIF
[11361]	1448	IF ( ( profdata%var(2)%nvqc(jobs) > kqc_cutoff ) .AND. &
	1449	& ( profdata%var(1)%nvqc(jobs) <= kqc_cutoff) ) THEN
	1450	profdata%var(1)%nvqc(jobs) = IBSET(profdata%var(1)%nvqc(jobs),15)
[2128]	1451	knumu = knumu + 1
	1452	ENDIF
	1453
	1454	END DO
	1455
	1456	END DO
	1457
	1458	END SUBROUTINE obs_uv_rej
	1459
	1460	END MODULE obs_prep

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

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

Context Navigation

source: NEMO/branches/UKMO/r8395_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90

Download in other formats: