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_read_prof.F90 in branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

Context Navigation

source: branches/UKMO/dev_r5518_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_prof.F90 @ 11468

Last change on this file since 11468 was 11468, checked in by mattmartin, 5 years ago
Merged changes to allow writing of climatological information to feedback files.
File size: 31.1 KB

Rev	Line
[2128]	1	MODULE obs_read_prof
	2	!!======================================================================
	3	!! * MODULE obs_read_prof *
	4	!! Observation diagnostics: Read the T and S profile observations
	5	!!======================================================================
	6
	7	!!----------------------------------------------------------------------
	8	!! obs_rea_pro_dri : Driver for reading profile obs
	9	!!----------------------------------------------------------------------
	10
	11	!! * Modules used
	12	USE par_kind ! Precision variables
	13	USE par_oce ! Ocean parameters
	14	USE in_out_manager ! I/O manager
	15	USE dom_oce ! Ocean space and time domain variables
	16	USE obs_mpp ! MPP support routines for observation diagnostics
	17	USE julian ! Julian date routines
	18	USE obs_utils ! Observation operator utility functions
	19	USE obs_prep ! Prepare observation arrays
	20	USE obs_grid ! Grid search
	21	USE obs_sort ! Sorting observation arrays
	22	USE obs_profiles_def ! Profile definitions
	23	USE obs_conv ! Various conversion routines
	24	USE obs_types ! Observation type definitions
	25	USE netcdf ! NetCDF library
	26	USE obs_oper ! Observation operators
[2715]	27	USE lib_mpp ! For ctl_warn/stop
[7992]	28	USE obs_fbm ! Feedback routines
[2128]	29
	30	IMPLICIT NONE
	31
	32	!! * Routine accessibility
	33	PRIVATE
	34
[7992]	35	PUBLIC obs_rea_prof ! Read the profile observations
[2128]	36
[2287]	37	!!----------------------------------------------------------------------
	38	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
	39	!! $Id$
	40	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
	41	!!----------------------------------------------------------------------
	42
[2128]	43	CONTAINS
[7992]	44
	45	SUBROUTINE obs_rea_prof( profdata, knumfiles, cdfilenames, &
	46	& kvars, kextr, kstp, ddobsini, ddobsend, &
[9306]	47	& ldvar, ldignmis, ldsatt, &
[11468]	48	& ldmod, ldclim, kdailyavtypes )
[2128]	49	!!---------------------------------------------------------------------
	50	!!
[7992]	51	!! * ROUTINE obs_rea_prof *
[2128]	52	!!
	53	!! ** Purpose : Read from file the profile observations
	54	!!
[7992]	55	!! ** Method : Read feedback data in and transform to NEMO internal
	56	!! profile data structure
[2128]	57	!!
	58	!! ** Action :
	59	!!
	60	!! References :
	61	!!
	62	!! History :
	63	!! ! : 2009-09 (K. Mogensen) : New merged version of old routines
[7992]	64	!! ! : 2015-08 (M. Martin) : Merged profile and velocity routines
[2128]	65	!!----------------------------------------------------------------------
[7992]	66
[2128]	67	!! * Arguments
[7992]	68	TYPE(obs_prof), INTENT(OUT) :: &
	69	& profdata ! Profile data to be read
	70	INTEGER, INTENT(IN) :: knumfiles ! Number of files to read
[2128]	71	CHARACTER(LEN=128), INTENT(IN) :: &
[7992]	72	& cdfilenames(knumfiles) ! File names to read in
[2128]	73	INTEGER, INTENT(IN) :: kvars ! Number of variables in profdata
[7992]	74	INTEGER, INTENT(IN) :: kextr ! Number of extra fields for each var
	75	INTEGER, INTENT(IN) :: kstp ! Ocean time-step index
[9306]	76	LOGICAL, DIMENSION(kvars), INTENT(IN) :: ldvar ! Observed variables switches
[7992]	77	LOGICAL, INTENT(IN) :: ldignmis ! Ignore missing files
	78	LOGICAL, INTENT(IN) :: ldsatt ! Compute salinity at all temperature points
	79	LOGICAL, INTENT(IN) :: ldmod ! Initialize model from input data
[11468]	80	LOGICAL, INTENT(IN) :: ldclim ! Set flag to show climatology will be output
[7992]	81	REAL(dp), INTENT(IN) :: ddobsini ! Obs. ini time in YYYYMMDD.HHMMSS
	82	REAL(dp), INTENT(IN) :: ddobsend ! Obs. end time in YYYYMMDD.HHMMSS
[2128]	83	INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &
[7992]	84	& kdailyavtypes ! Types of daily average observations
[2128]	85
	86	!! * Local declarations
[7992]	87	CHARACTER(LEN=15), PARAMETER :: cpname='obs_rea_prof'
	88	CHARACTER(len=8) :: clrefdate
	89	CHARACTER(len=8), DIMENSION(:), ALLOCATABLE :: clvars
[2128]	90	INTEGER :: jvar
	91	INTEGER :: ji
	92	INTEGER :: jj
	93	INTEGER :: jk
	94	INTEGER :: ij
	95	INTEGER :: iflag
	96	INTEGER :: inobf
	97	INTEGER :: i_file_id
	98	INTEGER :: inowin
	99	INTEGER :: iyea
	100	INTEGER :: imon
	101	INTEGER :: iday
	102	INTEGER :: ihou
	103	INTEGER :: imin
	104	INTEGER :: isec
[7992]	105	INTEGER :: iprof
	106	INTEGER :: iproftot
[9306]	107	INTEGER, DIMENSION(kvars) :: ivart0
	108	INTEGER, DIMENSION(kvars) :: ivart
[7992]	109	INTEGER :: ip3dt
	110	INTEGER :: ios
	111	INTEGER :: ioserrcount
[9306]	112	INTEGER, DIMENSION(kvars) :: ivartmpp
[7992]	113	INTEGER :: ip3dtmpp
	114	INTEGER :: itype
[2128]	115	INTEGER, DIMENSION(knumfiles) :: &
	116	& irefdate
[9306]	117	INTEGER, DIMENSION(ntyp1770+1,kvars) :: &
	118	& itypvar, &
	119	& itypvarmpp
	120	INTEGER, DIMENSION(:,:), ALLOCATABLE :: &
	121	& iobsi, &
	122	& iobsj, &
	123	& iproc
[2128]	124	INTEGER, DIMENSION(:), ALLOCATABLE :: &
	125	& iindx, &
	126	& ifileidx, &
	127	& iprofidx
	128	INTEGER, DIMENSION(imaxavtypes) :: &
	129	& idailyavtypes
[7992]	130	INTEGER, DIMENSION(kvars) :: &
	131	& iv3dt
[2128]	132	REAL(wp), DIMENSION(:), ALLOCATABLE :: &
	133	& zphi, &
	134	& zlam
[7992]	135	REAL(wp), DIMENSION(:), ALLOCATABLE :: &
[2128]	136	& zdat
[7992]	137	REAL(wp), DIMENSION(knumfiles) :: &
	138	& djulini, &
	139	& djulend
[2128]	140	LOGICAL :: llvalprof
[7992]	141	LOGICAL :: lldavtimset
[9306]	142	LOGICAL :: llcycle
[2128]	143	TYPE(obfbdata), POINTER, DIMENSION(:) :: &
	144	& inpfiles
[7992]	145
[2128]	146	! Local initialization
	147	iprof = 0
[9306]	148	ivart0(:) = 0
[2128]	149	ip3dt = 0
	150
	151	! Daily average types
[7992]	152	lldavtimset = .FALSE.
[2128]	153	IF ( PRESENT(kdailyavtypes) ) THEN
	154	idailyavtypes(:) = kdailyavtypes(:)
[7992]	155	IF ( ANY (idailyavtypes(:) /= -1) ) lldavtimset = .TRUE.
[2128]	156	ELSE
	157	idailyavtypes(:) = -1
	158	ENDIF
	159
	160	!-----------------------------------------------------------------------
[7992]	161	! Count the number of files needed and allocate the obfbdata type
[2128]	162	!-----------------------------------------------------------------------
	163
	164	inobf = knumfiles
[7992]	165
[2128]	166	ALLOCATE( inpfiles(inobf) )
	167
	168	prof_files : DO jj = 1, inobf
[7992]	169
[2128]	170	!---------------------------------------------------------------------
	171	! Prints
	172	!---------------------------------------------------------------------
	173	IF(lwp) THEN
	174	WRITE(numout,*)
	175	WRITE(numout,*) ' obs_rea_pro_dri : Reading from file = ', &
[7992]	176	& TRIM( TRIM( cdfilenames(jj) ) )
[2128]	177	WRITE(numout,*) ' ~~~~~~~~~~~~~~~'
	178	WRITE(numout,*)
	179	ENDIF
	180
	181	!---------------------------------------------------------------------
	182	! Initialization: Open file and get dimensions only
	183	!---------------------------------------------------------------------
[7992]	184
	185	iflag = nf90_open( TRIM( cdfilenames(jj) ), nf90_nowrite, &
[2128]	186	& i_file_id )
[7992]	187
[2128]	188	IF ( iflag /= nf90_noerr ) THEN
	189
	190	IF ( ldignmis ) THEN
	191	inpfiles(jj)%nobs = 0
[7992]	192	CALL ctl_warn( 'File ' // TRIM( cdfilenames(jj) ) // &
[2128]	193	& ' not found' )
	194	ELSE
[7992]	195	CALL ctl_stop( 'File ' // TRIM( cdfilenames(jj) ) // &
[2128]	196	& ' not found' )
	197	ENDIF
	198
	199	ELSE
[7992]	200
[2128]	201	!------------------------------------------------------------------
[7992]	202	! Close the file since it is opened in read_obfbdata
[2128]	203	!------------------------------------------------------------------
[7992]	204
[2128]	205	iflag = nf90_close( i_file_id )
	206
	207	!------------------------------------------------------------------
	208	! Read the profile file into inpfiles
	209	!------------------------------------------------------------------
[7992]	210	CALL init_obfbdata( inpfiles(jj) )
	211	CALL read_obfbdata( TRIM( cdfilenames(jj) ), inpfiles(jj), &
	212	& ldgrid = .TRUE. )
	213
[9306]	214	IF ( inpfiles(jj)%nvar /= kvars ) THEN
[7992]	215	CALL ctl_stop( 'Feedback format error: ', &
[9306]	216	& ' unexpected number of vars in profile file' )
[7992]	217	ENDIF
	218
	219	IF ( ldmod .AND. ( inpfiles(jj)%nadd == 0 ) ) THEN
	220	CALL ctl_stop( 'Model not in input data' )
	221	ENDIF
	222
	223	IF ( jj == 1 ) THEN
	224	ALLOCATE( clvars( inpfiles(jj)%nvar ) )
	225	DO ji = 1, inpfiles(jj)%nvar
	226	clvars(ji) = inpfiles(jj)%cname(ji)
	227	END DO
[2128]	228	ELSE
[7992]	229	DO ji = 1, inpfiles(jj)%nvar
	230	IF ( inpfiles(jj)%cname(ji) /= clvars(ji) ) THEN
	231	CALL ctl_stop( 'Feedback file variables not consistent', &
	232	& ' with previous files for this type' )
	233	ENDIF
	234	END DO
[2128]	235	ENDIF
[7992]	236
[2128]	237	!------------------------------------------------------------------
	238	! Change longitude (-180,180)
	239	!------------------------------------------------------------------
	240
	241	DO ji = 1, inpfiles(jj)%nobs
	242
	243	IF ( inpfiles(jj)%plam(ji) < -180. ) &
	244	& inpfiles(jj)%plam(ji) = inpfiles(jj)%plam(ji) + 360.
	245
	246	IF ( inpfiles(jj)%plam(ji) > 180. ) &
	247	& inpfiles(jj)%plam(ji) = inpfiles(jj)%plam(ji) - 360.
	248
	249	END DO
	250
	251	!------------------------------------------------------------------
	252	! Calculate the date (change eventually)
	253	!------------------------------------------------------------------
[7992]	254	clrefdate=inpfiles(jj)%cdjuldref(1:8)
	255	READ(clrefdate,'(I8)') irefdate(jj)
	256
[2128]	257	CALL ddatetoymdhms( ddobsini, iyea, imon, iday, ihou, imin, isec )
	258	CALL greg2jul( isec, imin, ihou, iday, imon, iyea, djulini(jj), &
	259	& krefdate = irefdate(jj) )
	260	CALL ddatetoymdhms( ddobsend, iyea, imon, iday, ihou, imin, isec )
	261	CALL greg2jul( isec, imin, ihou, iday, imon, iyea, djulend(jj), &
	262	& krefdate = irefdate(jj) )
	263
[4990]	264	ioserrcount=0
[7992]	265	IF ( lldavtimset ) THEN
	266
	267	IF ( ANY ( idailyavtypes(:) /= -1 ) .AND. lwp) THEN
	268	WRITE(numout,*)' Resetting time of daily averaged', &
	269	& ' observations to the end of the day'
	270	ENDIF
	271
[2128]	272	DO ji = 1, inpfiles(jj)%nobs
[4990]	273	READ( inpfiles(jj)%cdtyp(ji), '(I4)', IOSTAT = ios, ERR = 900 ) itype
	274	900 IF ( ios /= 0 ) THEN
[7992]	275	! Set type to zero if there is a problem in the string conversion
	276	itype = 0
[4990]	277	ENDIF
[7992]	278
	279	IF ( ANY ( idailyavtypes(:) == itype ) ) THEN
	280	! for daily averaged data force the time
	281	! to be the last time-step of the day, but still within the day.
	282	IF ( inpfiles(jj)%ptim(ji) >= 0. ) THEN
	283	inpfiles(jj)%ptim(ji) = &
	284	& INT(inpfiles(jj)%ptim(ji)) + 0.9999
	285	ELSE
	286	inpfiles(jj)%ptim(ji) = &
	287	& INT(inpfiles(jj)%ptim(ji)) - 0.0001
	288	ENDIF
[2128]	289	ENDIF
[7992]	290
[2128]	291	END DO
[7992]	292
[2128]	293	ENDIF
[7992]	294
[2128]	295	IF ( inpfiles(jj)%nobs > 0 ) THEN
[7992]	296	inpfiles(jj)%iproc(:,:) = -1
	297	inpfiles(jj)%iobsi(:,:) = -1
	298	inpfiles(jj)%iobsj(:,:) = -1
[2128]	299	ENDIF
	300	inowin = 0
	301	DO ji = 1, inpfiles(jj)%nobs
[7992]	302	IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE
[9306]	303	llcycle = .TRUE.
	304	DO jvar = 1, kvars
	305	IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN
	306	llcycle = .FALSE.
	307	EXIT
	308	ENDIF
	309	END DO
	310	IF ( llcycle ) CYCLE
[2128]	311	IF ( ( inpfiles(jj)%ptim(ji) > djulini(jj) ) .AND. &
	312	& ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN
	313	inowin = inowin + 1
	314	ENDIF
	315	END DO
	316	ALLOCATE( zlam(inowin) )
	317	ALLOCATE( zphi(inowin) )
[9306]	318	ALLOCATE( iobsi(inowin,kvars) )
	319	ALLOCATE( iobsj(inowin,kvars) )
	320	ALLOCATE( iproc(inowin,kvars) )
[2128]	321	inowin = 0
	322	DO ji = 1, inpfiles(jj)%nobs
[7992]	323	IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE
[9306]	324	llcycle = .TRUE.
	325	DO jvar = 1, kvars
	326	IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN
	327	llcycle = .FALSE.
	328	EXIT
	329	ENDIF
	330	END DO
	331	IF ( llcycle ) CYCLE
[2128]	332	IF ( ( inpfiles(jj)%ptim(ji) > djulini(jj) ) .AND. &
	333	& ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN
	334	inowin = inowin + 1
	335	zlam(inowin) = inpfiles(jj)%plam(ji)
	336	zphi(inowin) = inpfiles(jj)%pphi(ji)
	337	ENDIF
	338	END DO
	339
[9306]	340	! Assume anything other than velocity is on T grid
	341	IF ( TRIM( inpfiles(jj)%cname(1) ) == 'UVEL' ) THEN
	342	CALL obs_grid_search( inowin, zlam, zphi, iobsi(:,1), iobsj(:,1), &
	343	& iproc(:,1), 'U' )
	344	CALL obs_grid_search( inowin, zlam, zphi, iobsi(:,2), iobsj(:,2), &
	345	& iproc(:,2), 'V' )
	346	ELSE
	347	CALL obs_grid_search( inowin, zlam, zphi, iobsi(:,1), iobsj(:,1), &
	348	& iproc(:,1), 'T' )
	349	IF ( kvars > 1 ) THEN
	350	DO jvar = 2, kvars
	351	iobsi(:,jvar) = iobsi(:,1)
	352	iobsj(:,jvar) = iobsj(:,1)
	353	iproc(:,jvar) = iproc(:,1)
	354	END DO
	355	ENDIF
[7992]	356	ENDIF
[2128]	357
	358	inowin = 0
	359	DO ji = 1, inpfiles(jj)%nobs
[7992]	360	IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE
[9306]	361	llcycle = .TRUE.
	362	DO jvar = 1, kvars
	363	IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN
	364	llcycle = .FALSE.
	365	EXIT
	366	ENDIF
	367	END DO
	368	IF ( llcycle ) CYCLE
[2128]	369	IF ( ( inpfiles(jj)%ptim(ji) > djulini(jj) ) .AND. &
	370	& ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN
	371	inowin = inowin + 1
[9306]	372	DO jvar = 1, kvars
	373	inpfiles(jj)%iproc(ji,jvar) = iproc(inowin,jvar)
	374	inpfiles(jj)%iobsi(ji,jvar) = iobsi(inowin,jvar)
	375	inpfiles(jj)%iobsj(ji,jvar) = iobsj(inowin,jvar)
	376	END DO
	377	IF ( kvars > 1 ) THEN
	378	DO jvar = 2, kvars
	379	IF ( inpfiles(jj)%iproc(ji,jvar) /= &
	380	& inpfiles(jj)%iproc(ji,1) ) THEN
	381	CALL ctl_stop( 'Error in obs_read_prof:', &
	382	& 'observation on different processors for different vars')
	383	ENDIF
	384	END DO
[7992]	385	ENDIF
[2128]	386	ENDIF
	387	END DO
[9306]	388	DEALLOCATE( zlam, zphi, iobsi, iobsj, iproc )
[2128]	389
	390	DO ji = 1, inpfiles(jj)%nobs
[7992]	391	IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE
[9306]	392	llcycle = .TRUE.
	393	DO jvar = 1, kvars
	394	IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN
	395	llcycle = .FALSE.
	396	EXIT
	397	ENDIF
	398	END DO
	399	IF ( llcycle ) CYCLE
[2128]	400	IF ( ( inpfiles(jj)%ptim(ji) > djulini(jj) ) .AND. &
	401	& ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN
	402	IF ( nproc == 0 ) THEN
	403	IF ( inpfiles(jj)%iproc(ji,1) > nproc ) CYCLE
	404	ELSE
	405	IF ( inpfiles(jj)%iproc(ji,1) /= nproc ) CYCLE
	406	ENDIF
	407	llvalprof = .FALSE.
[9306]	408	DO jvar = 1, kvars
	409	IF ( ldvar(jvar) ) THEN
	410	DO ij = 1,inpfiles(jj)%nlev
	411	IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) &
	412	& CYCLE
	413	IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,jvar),2) .AND. &
	414	& .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN
	415	ivart0(jvar) = ivart0(jvar) + 1
	416	ENDIF
	417	END DO
	418	ENDIF
	419	END DO
	420	DO ij = 1,inpfiles(jj)%nlev
[2128]	421	IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) &
	422	& CYCLE
[9306]	423	DO jvar = 1, kvars
	424	IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,jvar),2) .AND. &
	425	& .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. &
	426	& ldvar(jvar) ) ) THEN
	427	ip3dt = ip3dt + 1
	428	llvalprof = .TRUE.
	429	EXIT
	430	ENDIF
	431	END DO
	432	END DO
[2128]	433
	434	IF ( llvalprof ) iprof = iprof + 1
	435
	436	ENDIF
	437	END DO
	438
	439	ENDIF
	440
	441	END DO prof_files
[7992]	442
[2128]	443	!-----------------------------------------------------------------------
	444	! Get the time ordered indices of the input data
	445	!-----------------------------------------------------------------------
	446
	447	!---------------------------------------------------------------------
	448	! Loop over input data files to count total number of profiles
	449	!---------------------------------------------------------------------
	450	iproftot = 0
	451	DO jj = 1, inobf
	452	DO ji = 1, inpfiles(jj)%nobs
[7992]	453	IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE
[9306]	454	llcycle = .TRUE.
	455	DO jvar = 1, kvars
	456	IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN
	457	llcycle = .FALSE.
	458	EXIT
	459	ENDIF
	460	END DO
	461	IF ( llcycle ) CYCLE
[2128]	462	IF ( ( inpfiles(jj)%ptim(ji) > djulini(jj) ) .AND. &
	463	& ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN
	464	iproftot = iproftot + 1
	465	ENDIF
	466	END DO
	467	END DO
	468
	469	ALLOCATE( iindx(iproftot), ifileidx(iproftot), &
	470	& iprofidx(iproftot), zdat(iproftot) )
	471	jk = 0
	472	DO jj = 1, inobf
	473	DO ji = 1, inpfiles(jj)%nobs
[7992]	474	IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE
[9306]	475	llcycle = .TRUE.
	476	DO jvar = 1, kvars
	477	IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN
	478	llcycle = .FALSE.
	479	EXIT
	480	ENDIF
	481	END DO
	482	IF ( llcycle ) CYCLE
[2128]	483	IF ( ( inpfiles(jj)%ptim(ji) > djulini(jj) ) .AND. &
	484	& ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN
	485	jk = jk + 1
	486	ifileidx(jk) = jj
	487	iprofidx(jk) = ji
	488	zdat(jk) = inpfiles(jj)%ptim(ji)
	489	ENDIF
	490	END DO
	491	END DO
	492	CALL sort_dp_indx( iproftot, &
	493	& zdat, &
	494	& iindx )
[7992]	495
[2128]	496	iv3dt(:) = -1
	497	IF (ldsatt) THEN
[9306]	498	iv3dt(:) = ip3dt
[2128]	499	ELSE
[9306]	500	iv3dt(:) = ivart0(:)
[2128]	501	ENDIF
	502	CALL obs_prof_alloc( profdata, kvars, kextr, iprof, iv3dt, &
[11468]	503	& kstp, jpi, jpj, jpk, ldclim )
[7992]	504
[2128]	505	! * Read obs/positions, QC, all variable and assign to profdata
	506
	507	profdata%nprof = 0
	508	profdata%nvprot(:) = 0
[7992]	509	profdata%cvars(:) = clvars(:)
[2128]	510	iprof = 0
	511
	512	ip3dt = 0
[9306]	513	ivart(:) = 0
	514	itypvar (:,:) = 0
	515	itypvarmpp(:,:) = 0
[7992]	516
	517	ioserrcount = 0
[2128]	518	DO jk = 1, iproftot
[7992]	519
[2128]	520	jj = ifileidx(iindx(jk))
	521	ji = iprofidx(iindx(jk))
	522
[9306]	523	IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE
	524	llcycle = .TRUE.
	525	DO jvar = 1, kvars
	526	IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN
	527	llcycle = .FALSE.
	528	EXIT
	529	ENDIF
	530	END DO
	531	IF ( llcycle ) CYCLE
[2128]	532
	533	IF ( ( inpfiles(jj)%ptim(ji) > djulini(jj) ) .AND. &
	534	& ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN
[7992]	535
[2128]	536	IF ( nproc == 0 ) THEN
	537	IF ( inpfiles(jj)%iproc(ji,1) > nproc ) CYCLE
	538	ELSE
	539	IF ( inpfiles(jj)%iproc(ji,1) /= nproc ) CYCLE
	540	ENDIF
[7992]	541
[2128]	542	llvalprof = .FALSE.
	543
	544	IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE
	545
[7992]	546	IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE
[9306]	547	llcycle = .TRUE.
	548	DO jvar = 1, kvars
	549	IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN
	550	llcycle = .FALSE.
	551	EXIT
	552	ENDIF
	553	END DO
	554	IF ( llcycle ) CYCLE
[2128]	555
	556	loop_prof : DO ij = 1, inpfiles(jj)%nlev
[7992]	557
[2128]	558	IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) &
	559	& CYCLE
[7992]	560
[9306]	561	DO jvar = 1, kvars
	562	IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,jvar),2) .AND. &
	563	& .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN
[7992]	564
[9306]	565	llvalprof = .TRUE.
	566	EXIT loop_prof
[7992]	567
[9306]	568	ENDIF
	569	END DO
[7992]	570
[2128]	571	END DO loop_prof
[7992]	572
[2128]	573	! Set profile information
[7992]	574
[2128]	575	IF ( llvalprof ) THEN
[7992]	576
[2128]	577	iprof = iprof + 1
	578
	579	CALL jul2greg( isec, &
	580	& imin, &
	581	& ihou, &
	582	& iday, &
	583	& imon, &
	584	& iyea, &
	585	& inpfiles(jj)%ptim(ji), &
	586	& irefdate(jj) )
	587
	588
	589	! Profile time coordinates
	590	profdata%nyea(iprof) = iyea
	591	profdata%nmon(iprof) = imon
	592	profdata%nday(iprof) = iday
	593	profdata%nhou(iprof) = ihou
	594	profdata%nmin(iprof) = imin
[7992]	595
[2128]	596	! Profile space coordinates
	597	profdata%rlam(iprof) = inpfiles(jj)%plam(ji)
	598	profdata%rphi(iprof) = inpfiles(jj)%pphi(ji)
	599
	600	! Coordinate search parameters
[9306]	601	DO jvar = 1, kvars
	602	profdata%mi (iprof,jvar) = inpfiles(jj)%iobsi(ji,jvar)
	603	profdata%mj (iprof,jvar) = inpfiles(jj)%iobsj(ji,jvar)
	604	END DO
[7992]	605
[2128]	606	! Profile WMO number
	607	profdata%cwmo(iprof) = inpfiles(jj)%cdwmo(ji)
[7992]	608
[2128]	609	! Instrument type
[4990]	610	READ( inpfiles(jj)%cdtyp(ji), '(I4)', IOSTAT = ios, ERR = 901 ) itype
	611	901 IF ( ios /= 0 ) THEN
	612	IF (ioserrcount == 0) CALL ctl_warn ( 'Problem converting an instrument type to integer. Setting type to zero' )
	613	ioserrcount = ioserrcount + 1
	614	itype = 0
	615	ENDIF
[7992]	616
[2128]	617	profdata%ntyp(iprof) = itype
[7992]	618
[2128]	619	! QC stuff
	620
	621	profdata%nqc(iprof) = inpfiles(jj)%ioqc(ji)
	622	profdata%nqcf(:,iprof) = inpfiles(jj)%ioqcf(:,ji)
	623	profdata%ipqc(iprof) = inpfiles(jj)%ipqc(ji)
	624	profdata%ipqcf(:,iprof) = inpfiles(jj)%ipqcf(:,ji)
	625	profdata%itqc(iprof) = inpfiles(jj)%itqc(ji)
	626	profdata%itqcf(:,iprof) = inpfiles(jj)%itqcf(:,ji)
	627	profdata%ivqc(iprof,:) = inpfiles(jj)%ivqc(ji,:)
	628	profdata%ivqcf(:,iprof,:) = inpfiles(jj)%ivqcf(:,ji,:)
	629
	630	! Bookkeeping data to match profiles
	631	profdata%npidx(iprof) = iprof
	632	profdata%npfil(iprof) = iindx(jk)
	633
	634	! Observation QC flag (whole profile)
	635	profdata%nqc(iprof) = 0 !TODO
	636
[7992]	637	loop_p : DO ij = 1, inpfiles(jj)%nlev
	638
[2128]	639	IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) &
	640	& CYCLE
	641
	642	IF (ldsatt) THEN
	643
[9306]	644	DO jvar = 1, kvars
	645	IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,jvar),2) .AND. &
	646	& .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. &
	647	& ldvar(jvar) ) ) THEN
	648	ip3dt = ip3dt + 1
	649	EXIT
	650	ELSE IF ( jvar == kvars ) THEN
	651	CYCLE loop_p
	652	ENDIF
	653	END DO
[7992]	654
[2128]	655	ENDIF
	656
[9306]	657	DO jvar = 1, kvars
	658
	659	IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,jvar),2) .AND. &
	660	& .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. &
	661	& ldvar(jvar) ) .OR. ldsatt ) THEN
[7992]	662
[9306]	663	IF (ldsatt) THEN
[2128]	664
[9306]	665	ivart(jvar) = ip3dt
[2128]	666
[9306]	667	ELSE
[2128]	668
[9306]	669	ivart(jvar) = ivart(jvar) + 1
[7992]	670
[9306]	671	ENDIF
[2128]	672
[9306]	673	! Depth of jvar observation
	674	profdata%var(jvar)%vdep(ivart(jvar)) = &
	675	& inpfiles(jj)%pdep(ij,ji)
[7992]	676
[9306]	677	! Depth of jvar observation QC
	678	profdata%var(jvar)%idqc(ivart(jvar)) = &
	679	& inpfiles(jj)%idqc(ij,ji)
[7992]	680
[9306]	681	! Depth of jvar observation QC flags
	682	profdata%var(jvar)%idqcf(:,ivart(jvar)) = &
	683	& inpfiles(jj)%idqcf(:,ij,ji)
[7992]	684
[9306]	685	! Profile index
	686	profdata%var(jvar)%nvpidx(ivart(jvar)) = iprof
[7992]	687
[9306]	688	! Vertical index in original profile
	689	profdata%var(jvar)%nvlidx(ivart(jvar)) = ij
[2128]	690
[9306]	691	! Profile jvar value
	692	IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,jvar),2) .AND. &
	693	& .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN
	694	profdata%var(jvar)%vobs(ivart(jvar)) = &
	695	& inpfiles(jj)%pob(ij,ji,jvar)
	696	IF ( ldmod ) THEN
	697	profdata%var(jvar)%vmod(ivart(jvar)) = &
	698	& inpfiles(jj)%padd(ij,ji,1,jvar)
	699	ENDIF
[11468]	700	IF ( profdata%lclim ) THEN
	701	profdata%var(jvar)%vclm(ivart(jvar)) = fbrmdi
	702	ENDIF
[9306]	703	! Count number of profile var1 data as function of type
	704	itypvar( profdata%ntyp(iprof) + 1, jvar ) = &
	705	& itypvar( profdata%ntyp(iprof) + 1, jvar ) + 1
	706	ELSE
	707	profdata%var(jvar)%vobs(ivart(jvar)) = fbrmdi
[2128]	708	ENDIF
	709
[9306]	710	! Profile jvar qc
	711	profdata%var(jvar)%nvqc(ivart(jvar)) = &
	712	& inpfiles(jj)%ivlqc(ij,ji,jvar)
[2128]	713
[9306]	714	! Profile jvar qc flags
	715	profdata%var(jvar)%nvqcf(:,ivart(jvar)) = &
	716	& inpfiles(jj)%ivlqcf(:,ij,ji,jvar)
[2128]	717
[9306]	718	! Profile insitu T value
	719	IF ( TRIM( inpfiles(jj)%cname(jvar) ) == 'POTM' ) THEN
	720	profdata%var(jvar)%vext(ivart(jvar),1) = &
	721	& inpfiles(jj)%pext(ij,ji,1)
	722	ENDIF
[7992]	723
[2128]	724	ENDIF
[9306]	725
	726	END DO
[2128]	727
	728	END DO loop_p
	729
	730	ENDIF
	731
	732	ENDIF
	733
	734	END DO
	735
	736	!-----------------------------------------------------------------------
	737	! Sum up over processors
	738	!-----------------------------------------------------------------------
[7992]	739
[9306]	740	DO jvar = 1, kvars
	741	CALL obs_mpp_sum_integer ( ivart0(jvar), ivartmpp(jvar) )
	742	END DO
[7992]	743	CALL obs_mpp_sum_integer ( ip3dt, ip3dtmpp )
	744
[9306]	745	DO jvar = 1, kvars
	746	CALL obs_mpp_sum_integers( itypvar(:,jvar), itypvarmpp(:,jvar), ntyp1770 + 1 )
	747	END DO
[7992]	748
[2128]	749	!-----------------------------------------------------------------------
	750	! Output number of observations.
	751	!-----------------------------------------------------------------------
	752	IF(lwp) THEN
	753	WRITE(numout,*)
[7992]	754	WRITE(numout,'(A)') ' Profile data'
[2128]	755	WRITE(numout,'(1X,A)') '------------'
	756	WRITE(numout,*)
[9306]	757	DO jvar = 1, kvars
	758	WRITE(numout,'(1X,A)') 'Profile data, '//TRIM( profdata%cvars(jvar) )
	759	WRITE(numout,'(1X,A)') '------------------------'
	760	DO ji = 0, ntyp1770
	761	IF ( itypvarmpp(ji+1,jvar) > 0 ) THEN
	762	WRITE(numout,'(1X,A3,1X,A48,A3,I8)') ctypshort(ji), &
	763	& cwmonam1770(ji)(1:52),' = ', &
	764	& itypvarmpp(ji+1,jvar)
	765	ENDIF
	766	END DO
	767	WRITE(numout,'(1X,A)') &
	768	& '---------------------------------------------------------------'
	769	WRITE(numout,'(1X,A55,I8)') &
	770	& 'Total profile data for variable '//TRIM( profdata%cvars(jvar) )// &
	771	& ' = ', ivartmpp(jvar)
	772	WRITE(numout,'(1X,A)') &
	773	& '---------------------------------------------------------------'
	774	WRITE(numout,*)
[2128]	775	END DO
	776	ENDIF
[7992]	777
[2128]	778	IF (ldsatt) THEN
[9306]	779	profdata%nvprot(:) = ip3dt
	780	profdata%nvprotmpp(:) = ip3dtmpp
[2128]	781	ELSE
[9306]	782	DO jvar = 1, kvars
	783	profdata%nvprot(jvar) = ivart(jvar)
	784	profdata%nvprotmpp(jvar) = ivartmpp(jvar)
	785	END DO
[2128]	786	ENDIF
	787	profdata%nprof = iprof
	788
	789	!-----------------------------------------------------------------------
	790	! Model level search
	791	!-----------------------------------------------------------------------
[9306]	792	DO jvar = 1, kvars
	793	IF ( ldvar(jvar) ) THEN
	794	CALL obs_level_search( jpk, gdept_1d, &
	795	& profdata%nvprot(jvar), profdata%var(jvar)%vdep, &
	796	& profdata%var(jvar)%mvk )
	797	ENDIF
	798	END DO
[7992]	799
[2128]	800	!-----------------------------------------------------------------------
	801	! Set model equivalent to 99999
	802	!-----------------------------------------------------------------------
	803	IF ( .NOT. ldmod ) THEN
	804	DO jvar = 1, kvars
	805	profdata%var(jvar)%vmod(:) = fbrmdi
	806	END DO
	807	ENDIF
	808	!-----------------------------------------------------------------------
	809	! Deallocate temporary data
	810	!-----------------------------------------------------------------------
[7992]	811	DEALLOCATE( ifileidx, iprofidx, zdat, clvars )
[2128]	812
	813	!-----------------------------------------------------------------------
	814	! Deallocate input data
	815	!-----------------------------------------------------------------------
	816	DO jj = 1, inobf
	817	CALL dealloc_obfbdata( inpfiles(jj) )
	818	END DO
	819	DEALLOCATE( inpfiles )
	820
[7992]	821	END SUBROUTINE obs_rea_prof
[2128]	822
	823	END MODULE obs_read_prof

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

Download in other formats: