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.

obs_prep.F90 in NEMO/trunk/src/OCE/OBS – NEMO

Context Navigation

source: NEMO/trunk/src/OCE/OBS/obs_prep.F90 @ 14587

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

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

Download in other formats: