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/branches/2020/dev_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/src/OCE/OBS – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/src/OCE/OBS/obs_prep.F90 @ 14652

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

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

Download in other formats: