Context Navigation

obs_read_prof.F90

Last change on this file was 15799, checked in by dford, 2 years ago
More generic interface and structure for OBS code. See Met Office utils tickets 471 and 530.
File size: 39.6 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
[6140]	28	USE obs_fbm ! Feedback routines
[15799]	29	USE obs_group_def, ONLY : & ! Observation variable information
	30	& cobsname_uvel, &
	31	& cobsname_vvel, &
	32	& imaxavtypes
[2128]	33
	34	IMPLICIT NONE
	35
	36	!! * Routine accessibility
	37	PRIVATE
	38
[6140]	39	PUBLIC obs_rea_prof ! Read the profile observations
[2128]	40
[2287]	41	!!----------------------------------------------------------------------
[9598]	42	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[2287]	43	!! $Id$
[10068]	44	!! Software governed by the CeCILL license (see ./LICENSE)
[2287]	45	!!----------------------------------------------------------------------
	46
[2128]	47	CONTAINS
[6140]	48
	49	SUBROUTINE obs_rea_prof( profdata, knumfiles, cdfilenames, &
[15799]	50	& kvars, kadd, kextr, kstp, ddobsini, ddobsend, &
	51	& ldvar, ldignmis, ldallatall, &
	52	& ldmod, cdvars, kdailyavtypes )
[2128]	53	!!---------------------------------------------------------------------
	54	!!
[6140]	55	!! * ROUTINE obs_rea_prof *
[2128]	56	!!
	57	!! ** Purpose : Read from file the profile observations
	58	!!
[6140]	59	!! ** Method : Read feedback data in and transform to NEMO internal
	60	!! profile data structure
[2128]	61	!!
	62	!! ** Action :
	63	!!
	64	!! References :
	65	!!
	66	!! History :
	67	!! ! : 2009-09 (K. Mogensen) : New merged version of old routines
[6140]	68	!! ! : 2015-08 (M. Martin) : Merged profile and velocity routines
[2128]	69	!!----------------------------------------------------------------------
[6140]	70
[2128]	71	!! * Arguments
[6140]	72	TYPE(obs_prof), INTENT(OUT) :: &
	73	& profdata ! Profile data to be read
	74	INTEGER, INTENT(IN) :: knumfiles ! Number of files to read
[2128]	75	CHARACTER(LEN=128), INTENT(IN) :: &
[6140]	76	& cdfilenames(knumfiles) ! File names to read in
[2128]	77	INTEGER, INTENT(IN) :: kvars ! Number of variables in profdata
[15799]	78	INTEGER, INTENT(IN) :: kadd ! Number of additional fields
	79	! in addition to those in the input file(s)
	80	INTEGER, INTENT(IN) :: kextr ! Number of extra fields
	81	! in addition to those in the input file(s)
[6140]	82	INTEGER, INTENT(IN) :: kstp ! Ocean time-step index
[15799]	83	LOGICAL, DIMENSION(kvars), INTENT(IN) :: ldvar ! Observed variables switches
[6140]	84	LOGICAL, INTENT(IN) :: ldignmis ! Ignore missing files
[15799]	85	LOGICAL, INTENT(IN) :: ldallatall ! Compute salinity at all temperature points
[6140]	86	LOGICAL, INTENT(IN) :: ldmod ! Initialize model from input data
	87	REAL(dp), INTENT(IN) :: ddobsini ! Obs. ini time in YYYYMMDD.HHMMSS
	88	REAL(dp), INTENT(IN) :: ddobsend ! Obs. end time in YYYYMMDD.HHMMSS
[15799]	89	CHARACTER(len=8), DIMENSION(kvars), INTENT(IN) :: cdvars
[2128]	90	INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: &
[6140]	91	& kdailyavtypes ! Types of daily average observations
[2128]	92
	93	!! * Local declarations
[6140]	94	CHARACTER(LEN=15), PARAMETER :: cpname='obs_rea_prof'
	95	CHARACTER(len=8) :: clrefdate
[15799]	96	CHARACTER(len=ilenname), DIMENSION(:), ALLOCATABLE :: clvarsin
	97	CHARACTER(len=ilenlong), DIMENSION(:), ALLOCATABLE :: cllongin
	98	CHARACTER(len=ilenunit), DIMENSION(:), ALLOCATABLE :: clunitin
	99	CHARACTER(len=ilengrid), DIMENSION(:), ALLOCATABLE :: clgridin
	100	CHARACTER(len=ilenname), DIMENSION(:), ALLOCATABLE :: claddvarsin
	101	CHARACTER(len=ilenlong), DIMENSION(:,:), ALLOCATABLE :: claddlongin
	102	CHARACTER(len=ilenunit), DIMENSION(:,:), ALLOCATABLE :: claddunitin
	103	CHARACTER(len=ilenname), DIMENSION(:), ALLOCATABLE :: clextvarsin
	104	CHARACTER(len=ilenlong), DIMENSION(:), ALLOCATABLE :: clextlongin
	105	CHARACTER(len=ilenunit), DIMENSION(:), ALLOCATABLE :: clextunitin
[2128]	106	INTEGER :: ji
	107	INTEGER :: jj
	108	INTEGER :: jk
	109	INTEGER :: ij
[15799]	110	INTEGER :: jind
	111	INTEGER :: jext
	112	INTEGER :: jvar
	113	INTEGER :: jadd
	114	INTEGER :: jadd2
	115	INTEGER :: iadd
	116	INTEGER :: iaddin
	117	INTEGER :: iextr
[2128]	118	INTEGER :: iflag
	119	INTEGER :: inobf
	120	INTEGER :: i_file_id
	121	INTEGER :: inowin
	122	INTEGER :: iyea
	123	INTEGER :: imon
	124	INTEGER :: iday
	125	INTEGER :: ihou
	126	INTEGER :: imin
	127	INTEGER :: isec
[6140]	128	INTEGER :: iprof
	129	INTEGER :: iproftot
[15799]	130	INTEGER, DIMENSION(kvars) :: ivart0
	131	INTEGER, DIMENSION(kvars) :: ivart
[6140]	132	INTEGER :: ip3dt
	133	INTEGER :: ios
	134	INTEGER :: ioserrcount
[15799]	135	INTEGER, DIMENSION(kvars) :: ivartmpp
[6140]	136	INTEGER :: ip3dtmpp
	137	INTEGER :: itype
[2128]	138	INTEGER, DIMENSION(knumfiles) :: &
	139	& irefdate
[15799]	140	INTEGER, DIMENSION(ntyp1770+1,kvars) :: &
	141	& itypvar, &
	142	& itypvarmpp
	143	INTEGER, DIMENSION(:,:), ALLOCATABLE :: &
	144	& iobsi, &
	145	& iobsj, &
	146	& iproc
[2128]	147	INTEGER, DIMENSION(:), ALLOCATABLE :: &
	148	& iindx, &
	149	& ifileidx, &
	150	& iprofidx
	151	INTEGER, DIMENSION(imaxavtypes) :: &
	152	& idailyavtypes
[6140]	153	INTEGER, DIMENSION(kvars) :: &
	154	& iv3dt
[2128]	155	REAL(wp), DIMENSION(:), ALLOCATABLE :: &
	156	& zphi, &
	157	& zlam
[6140]	158	REAL(wp), DIMENSION(:), ALLOCATABLE :: &
[2128]	159	& zdat
[6140]	160	REAL(wp), DIMENSION(knumfiles) :: &
	161	& djulini, &
	162	& djulend
[2128]	163	LOGICAL :: llvalprof
[6140]	164	LOGICAL :: lldavtimset
[15799]	165	LOGICAL :: llcycle
	166	LOGICAL :: llpotm
[2128]	167	TYPE(obfbdata), POINTER, DIMENSION(:) :: &
	168	& inpfiles
[6140]	169
[2128]	170	! Local initialization
	171	iprof = 0
[15799]	172	ivart0(:) = 0
[2128]	173	ip3dt = 0
	174
	175	! Daily average types
[6140]	176	lldavtimset = .FALSE.
[2128]	177	IF ( PRESENT(kdailyavtypes) ) THEN
	178	idailyavtypes(:) = kdailyavtypes(:)
[6140]	179	IF ( ANY (idailyavtypes(:) /= -1) ) lldavtimset = .TRUE.
[2128]	180	ELSE
	181	idailyavtypes(:) = -1
	182	ENDIF
	183
	184	!-----------------------------------------------------------------------
[6140]	185	! Count the number of files needed and allocate the obfbdata type
[2128]	186	!-----------------------------------------------------------------------
	187
	188	inobf = knumfiles
[6140]	189
[2128]	190	ALLOCATE( inpfiles(inobf) )
	191
[15799]	192	iadd = 0
	193	iextr = 0
	194
[2128]	195	prof_files : DO jj = 1, inobf
[6140]	196
[2128]	197	!---------------------------------------------------------------------
	198	! Prints
	199	!---------------------------------------------------------------------
	200	IF(lwp) THEN
	201	WRITE(numout,*)
	202	WRITE(numout,*) ' obs_rea_pro_dri : Reading from file = ', &
[6140]	203	& TRIM( TRIM( cdfilenames(jj) ) )
[2128]	204	WRITE(numout,*) ' ~~~~~~~~~~~~~~~'
	205	WRITE(numout,*)
	206	ENDIF
	207
	208	!---------------------------------------------------------------------
	209	! Initialization: Open file and get dimensions only
	210	!---------------------------------------------------------------------
[6140]	211
	212	iflag = nf90_open( TRIM( cdfilenames(jj) ), nf90_nowrite, &
[2128]	213	& i_file_id )
[6140]	214
[2128]	215	IF ( iflag /= nf90_noerr ) THEN
	216
	217	IF ( ldignmis ) THEN
	218	inpfiles(jj)%nobs = 0
[6140]	219	CALL ctl_warn( 'File ' // TRIM( cdfilenames(jj) ) // &
[2128]	220	& ' not found' )
	221	ELSE
[6140]	222	CALL ctl_stop( 'File ' // TRIM( cdfilenames(jj) ) // &
[2128]	223	& ' not found' )
	224	ENDIF
	225
	226	ELSE
[6140]	227
[2128]	228	!------------------------------------------------------------------
[6140]	229	! Close the file since it is opened in read_obfbdata
[2128]	230	!------------------------------------------------------------------
[6140]	231
[2128]	232	iflag = nf90_close( i_file_id )
	233
	234	!------------------------------------------------------------------
	235	! Read the profile file into inpfiles
	236	!------------------------------------------------------------------
[6140]	237	CALL init_obfbdata( inpfiles(jj) )
	238	CALL read_obfbdata( TRIM( cdfilenames(jj) ), inpfiles(jj), &
	239	& ldgrid = .TRUE. )
	240
[15799]	241	IF ( inpfiles(jj)%nvar /= kvars ) THEN
[6140]	242	CALL ctl_stop( 'Feedback format error: ', &
[15799]	243	& ' unexpected number of vars in feedback file', &
	244	& TRIM(cdfilenames(jj)) )
[6140]	245	ENDIF
	246
	247	IF ( ldmod .AND. ( inpfiles(jj)%nadd == 0 ) ) THEN
[15799]	248	CALL ctl_stop( 'Model not in input data in', &
	249	& TRIM(cdfilenames(jj)) )
	250	RETURN
[6140]	251	ENDIF
	252
[15799]	253	IF ( (iextr > 0) .AND. (inpfiles(jj)%next /= iextr) ) THEN
	254	CALL ctl_stop( 'Number of extra variables not consistent', &
	255	& ' with previous files for this type in', &
	256	& TRIM(cdfilenames(jj)) )
	257	ELSE
	258	iextr = inpfiles(jj)%next
	259	ENDIF
	260
	261	! Ignore model counterpart
	262	iaddin = inpfiles(jj)%nadd
	263	DO ji = 1, iaddin
	264	IF ( TRIM(inpfiles(jj)%caddname(ji)) == 'Hx' ) THEN
	265	iaddin = iaddin - 1
	266	EXIT
	267	ENDIF
	268	END DO
	269	IF ( ldmod .AND. ( inpfiles(jj)%nadd == iaddin ) ) THEN
	270	CALL ctl_stop( 'Model not in input data', &
	271	& TRIM(cdfilenames(jj)) )
	272	ENDIF
	273
	274	IF ( (iadd > 0) .AND. (iaddin /= iadd) ) THEN
	275	CALL ctl_stop( 'Number of additional variables not consistent', &
	276	& ' with previous files for this type in', &
	277	& TRIM(cdfilenames(jj)) )
	278	ELSE
	279	iadd = iaddin
	280	ENDIF
	281
[6140]	282	IF ( jj == 1 ) THEN
[15799]	283	ALLOCATE( clvarsin( inpfiles(jj)%nvar ) )
	284	ALLOCATE( cllongin( inpfiles(jj)%nvar ) )
	285	ALLOCATE( clunitin( inpfiles(jj)%nvar ) )
	286	ALLOCATE( clgridin( inpfiles(jj)%nvar ) )
[6140]	287	DO ji = 1, inpfiles(jj)%nvar
[15799]	288	clvarsin(ji) = inpfiles(jj)%cname(ji)
	289	cllongin(ji) = inpfiles(jj)%coblong(ji)
	290	clunitin(ji) = inpfiles(jj)%cobunit(ji)
	291	clgridin(ji) = inpfiles(jj)%cgrid(ji)
	292	IF ( clvarsin(ji) /= cdvars(ji) ) THEN
	293	CALL ctl_stop( 'Feedback file variables do not match', &
	294	& ' expected variable names for this type' )
	295	ENDIF
[6140]	296	END DO
[15799]	297	IF ( iadd > 0 ) THEN
	298	ALLOCATE( claddvarsin( iadd ) )
	299	ALLOCATE( claddlongin( iadd, inpfiles(jj)%nvar ) )
	300	ALLOCATE( claddunitin( iadd, inpfiles(jj)%nvar ) )
	301	jadd = 0
	302	DO ji = 1, inpfiles(jj)%nadd
	303	IF ( TRIM(inpfiles(jj)%caddname(ji)) /= 'Hx' ) THEN
	304	jadd = jadd + 1
	305	claddvarsin(jadd) = inpfiles(jj)%caddname(ji)
	306	DO jk = 1, inpfiles(jj)%nvar
	307	claddlongin(jadd,jk) = inpfiles(jj)%caddlong(ji,jk)
	308	claddunitin(jadd,jk) = inpfiles(jj)%caddunit(ji,jk)
	309	END DO
	310	ENDIF
	311	END DO
	312	ENDIF
	313	IF ( iextr > 0 ) THEN
	314	ALLOCATE( clextvarsin( iextr ) )
	315	ALLOCATE( clextlongin( iextr ) )
	316	ALLOCATE( clextunitin( iextr ) )
	317	DO ji = 1, iextr
	318	clextvarsin(ji) = inpfiles(jj)%cextname(ji)
	319	clextlongin(ji) = inpfiles(jj)%cextlong(ji)
	320	clextunitin(ji) = inpfiles(jj)%cextunit(ji)
	321	END DO
	322	ENDIF
[2128]	323	ELSE
[6140]	324	DO ji = 1, inpfiles(jj)%nvar
[15799]	325	IF ( inpfiles(jj)%cname(ji) /= clvarsin(ji) ) THEN
[6140]	326	CALL ctl_stop( 'Feedback file variables not consistent', &
[15799]	327	& ' with previous files for this type in', &
	328	& TRIM(cdfilenames(jj)) )
[6140]	329	ENDIF
	330	END DO
[15799]	331	IF ( iadd > 0 ) THEN
	332	jadd = 0
	333	DO ji = 1, inpfiles(jj)%nadd
	334	IF ( TRIM(inpfiles(jj)%caddname(ji)) /= 'Hx' ) THEN
	335	jadd = jadd + 1
	336	IF ( inpfiles(jj)%caddname(ji) /= claddvarsin(jadd) ) THEN
	337	CALL ctl_stop( 'Feedback file additional variables not consistent', &
	338	& ' with previous files for this type in', &
	339	& TRIM(cdfilenames(jj)) )
	340	ENDIF
	341	ENDIF
	342	END DO
	343	ENDIF
	344	IF ( iextr > 0 ) THEN
	345	DO ji = 1, iextr
	346	IF ( inpfiles(jj)%cextname(ji) /= clextvarsin(ji) ) THEN
	347	CALL ctl_stop( 'Feedback file extra variables not consistent', &
	348	& ' with previous files for this type in', &
	349	& TRIM(cdfilenames(jj)) )
	350	ENDIF
	351	END DO
	352	ENDIF
[2128]	353	ENDIF
[6140]	354
[2128]	355	!------------------------------------------------------------------
	356	! Change longitude (-180,180)
	357	!------------------------------------------------------------------
	358
	359	DO ji = 1, inpfiles(jj)%nobs
	360
	361	IF ( inpfiles(jj)%plam(ji) < -180. ) &
	362	& inpfiles(jj)%plam(ji) = inpfiles(jj)%plam(ji) + 360.
	363
	364	IF ( inpfiles(jj)%plam(ji) > 180. ) &
	365	& inpfiles(jj)%plam(ji) = inpfiles(jj)%plam(ji) - 360.
	366
	367	END DO
	368
	369	!------------------------------------------------------------------
	370	! Calculate the date (change eventually)
	371	!------------------------------------------------------------------
[6140]	372	clrefdate=inpfiles(jj)%cdjuldref(1:8)
	373	READ(clrefdate,'(I8)') irefdate(jj)
	374
[2128]	375	CALL ddatetoymdhms( ddobsini, iyea, imon, iday, ihou, imin, isec )
	376	CALL greg2jul( isec, imin, ihou, iday, imon, iyea, djulini(jj), &
	377	& krefdate = irefdate(jj) )
	378	CALL ddatetoymdhms( ddobsend, iyea, imon, iday, ihou, imin, isec )
	379	CALL greg2jul( isec, imin, ihou, iday, imon, iyea, djulend(jj), &
	380	& krefdate = irefdate(jj) )
	381
[4990]	382	ioserrcount=0
[6140]	383	IF ( lldavtimset ) THEN
	384
	385	IF ( ANY ( idailyavtypes(:) /= -1 ) .AND. lwp) THEN
	386	WRITE(numout,*)' Resetting time of daily averaged', &
	387	& ' observations to the end of the day'
	388	ENDIF
	389
[2128]	390	DO ji = 1, inpfiles(jj)%nobs
[4990]	391	READ( inpfiles(jj)%cdtyp(ji), '(I4)', IOSTAT = ios, ERR = 900 ) itype
	392	900 IF ( ios /= 0 ) THEN
[6140]	393	! Set type to zero if there is a problem in the string conversion
	394	itype = 0
[4990]	395	ENDIF
[6140]	396
	397	IF ( ANY ( idailyavtypes(:) == itype ) ) THEN
	398	! for daily averaged data force the time
	399	! to be the last time-step of the day, but still within the day.
	400	IF ( inpfiles(jj)%ptim(ji) >= 0. ) THEN
	401	inpfiles(jj)%ptim(ji) = &
	402	& INT(inpfiles(jj)%ptim(ji)) + 0.9999
	403	ELSE
	404	inpfiles(jj)%ptim(ji) = &
	405	& INT(inpfiles(jj)%ptim(ji)) - 0.0001
	406	ENDIF
[2128]	407	ENDIF
[6140]	408
[2128]	409	END DO
[6140]	410
[2128]	411	ENDIF
[6140]	412
[2128]	413	IF ( inpfiles(jj)%nobs > 0 ) THEN
[6140]	414	inpfiles(jj)%iproc(:,:) = -1
	415	inpfiles(jj)%iobsi(:,:) = -1
	416	inpfiles(jj)%iobsj(:,:) = -1
[2128]	417	ENDIF
	418	inowin = 0
	419	DO ji = 1, inpfiles(jj)%nobs
[9023]	420	IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE
[15799]	421	llcycle = .TRUE.
	422	DO jvar = 1, kvars
	423	IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN
	424	llcycle = .FALSE.
	425	EXIT
	426	ENDIF
	427	END DO
	428	IF ( llcycle ) CYCLE
[2128]	429	IF ( ( inpfiles(jj)%ptim(ji) > djulini(jj) ) .AND. &
	430	& ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN
	431	inowin = inowin + 1
	432	ENDIF
	433	END DO
	434	ALLOCATE( zlam(inowin) )
	435	ALLOCATE( zphi(inowin) )
[15799]	436	ALLOCATE( iobsi(inowin,kvars) )
	437	ALLOCATE( iobsj(inowin,kvars) )
	438	ALLOCATE( iproc(inowin,kvars) )
[2128]	439	inowin = 0
	440	DO ji = 1, inpfiles(jj)%nobs
[9023]	441	IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE
[15799]	442	llcycle = .TRUE.
	443	DO jvar = 1, kvars
	444	IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN
	445	llcycle = .FALSE.
	446	EXIT
	447	ENDIF
	448	END DO
	449	IF ( llcycle ) CYCLE
[2128]	450	IF ( ( inpfiles(jj)%ptim(ji) > djulini(jj) ) .AND. &
	451	& ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN
	452	inowin = inowin + 1
	453	zlam(inowin) = inpfiles(jj)%plam(ji)
	454	zphi(inowin) = inpfiles(jj)%pphi(ji)
	455	ENDIF
	456	END DO
	457
[15799]	458	! Do grid search
	459	! Assume anything other than velocity is on T grid
	460	! Save resource by not repeating for the same grid
	461	jind = 0
	462	DO jvar = 1, kvars
	463	IF ( TRIM(inpfiles(jj)%cname(jvar)) == cobsname_uvel ) THEN
	464	CALL obs_grid_search( inowin, zlam, zphi, iobsi(:,jvar), iobsj(:,jvar), &
	465	& iproc(:,jvar), 'U' )
	466	ELSE IF ( TRIM(inpfiles(jj)%cname(jvar)) == cobsname_vvel ) THEN
	467	CALL obs_grid_search( inowin, zlam, zphi, iobsi(:,jvar), iobsj(:,jvar), &
	468	& iproc(:,jvar), 'V' )
	469	ELSE
	470	IF ( jind > 0 ) THEN
	471	iobsi(:,jvar) = iobsi(:,jind)
	472	iobsj(:,jvar) = iobsj(:,jind)
	473	iproc(:,jvar) = iproc(:,jind)
	474	ELSE
	475	jind = jvar
	476	CALL obs_grid_search( inowin, zlam, zphi, iobsi(:,jvar), iobsj(:,jvar), &
	477	& iproc(:,jvar), 'T' )
	478	ENDIF
	479	ENDIF
	480	END DO
[2128]	481
	482	inowin = 0
	483	DO ji = 1, inpfiles(jj)%nobs
[9023]	484	IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE
[15799]	485	llcycle = .TRUE.
	486	DO jvar = 1, kvars
	487	IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN
	488	llcycle = .FALSE.
	489	EXIT
	490	ENDIF
	491	END DO
	492	IF ( llcycle ) CYCLE
[2128]	493	IF ( ( inpfiles(jj)%ptim(ji) > djulini(jj) ) .AND. &
	494	& ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN
	495	inowin = inowin + 1
[15799]	496	DO jvar = 1, kvars
	497	inpfiles(jj)%iproc(ji,jvar) = iproc(inowin,jvar)
	498	inpfiles(jj)%iobsi(ji,jvar) = iobsi(inowin,jvar)
	499	inpfiles(jj)%iobsj(ji,jvar) = iobsj(inowin,jvar)
	500	END DO
	501	IF ( kvars > 1 ) THEN
	502	DO jvar = 2, kvars
	503	IF ( inpfiles(jj)%iproc(ji,jvar) /= &
	504	& inpfiles(jj)%iproc(ji,1) ) THEN
	505	CALL ctl_stop( 'Error in obs_read_prof:', &
	506	& 'observation on different processors for different vars')
	507	ENDIF
	508	END DO
[6140]	509	ENDIF
[2128]	510	ENDIF
	511	END DO
[15799]	512	DEALLOCATE( zlam, zphi, iobsi, iobsj, iproc )
[2128]	513
	514	DO ji = 1, inpfiles(jj)%nobs
[9023]	515	IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE
[15799]	516	llcycle = .TRUE.
	517	DO jvar = 1, kvars
	518	IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN
	519	llcycle = .FALSE.
	520	EXIT
	521	ENDIF
	522	END DO
	523	IF ( llcycle ) CYCLE
[2128]	524	IF ( ( inpfiles(jj)%ptim(ji) > djulini(jj) ) .AND. &
	525	& ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN
	526	IF ( nproc == 0 ) THEN
	527	IF ( inpfiles(jj)%iproc(ji,1) > nproc ) CYCLE
	528	ELSE
	529	IF ( inpfiles(jj)%iproc(ji,1) /= nproc ) CYCLE
	530	ENDIF
	531	llvalprof = .FALSE.
[15799]	532	DO jvar = 1, kvars
	533	IF ( ldvar(jvar) ) THEN
	534	DO ij = 1,inpfiles(jj)%nlev
	535	IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) &
	536	& CYCLE
	537	IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,jvar),2) .AND. &
	538	& .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN
	539	ivart0(jvar) = ivart0(jvar) + 1
	540	ENDIF
	541	END DO
	542	ENDIF
	543	END DO
	544	DO ij = 1,inpfiles(jj)%nlev
[2128]	545	IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) &
	546	& CYCLE
[15799]	547	DO jvar = 1, kvars
	548	IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,jvar),2) .AND. &
	549	& .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. &
	550	& ldvar(jvar) ) ) THEN
	551	ip3dt = ip3dt + 1
	552	llvalprof = .TRUE.
	553	EXIT
	554	ENDIF
	555	END DO
	556	END DO
[2128]	557
	558	IF ( llvalprof ) iprof = iprof + 1
	559
	560	ENDIF
	561	END DO
	562
	563	ENDIF
	564
	565	END DO prof_files
[6140]	566
[2128]	567	!-----------------------------------------------------------------------
	568	! Get the time ordered indices of the input data
	569	!-----------------------------------------------------------------------
	570
	571	!---------------------------------------------------------------------
	572	! Loop over input data files to count total number of profiles
	573	!---------------------------------------------------------------------
	574	iproftot = 0
	575	DO jj = 1, inobf
	576	DO ji = 1, inpfiles(jj)%nobs
[9023]	577	IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE
[15799]	578	llcycle = .TRUE.
	579	DO jvar = 1, kvars
	580	IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN
	581	llcycle = .FALSE.
	582	EXIT
	583	ENDIF
	584	END DO
	585	IF ( llcycle ) CYCLE
[2128]	586	IF ( ( inpfiles(jj)%ptim(ji) > djulini(jj) ) .AND. &
	587	& ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN
	588	iproftot = iproftot + 1
	589	ENDIF
	590	END DO
	591	END DO
	592
	593	ALLOCATE( iindx(iproftot), ifileidx(iproftot), &
	594	& iprofidx(iproftot), zdat(iproftot) )
	595	jk = 0
	596	DO jj = 1, inobf
	597	DO ji = 1, inpfiles(jj)%nobs
[9023]	598	IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE
[15799]	599	llcycle = .TRUE.
	600	DO jvar = 1, kvars
	601	IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN
	602	llcycle = .FALSE.
	603	EXIT
	604	ENDIF
	605	END DO
	606	IF ( llcycle ) CYCLE
[2128]	607	IF ( ( inpfiles(jj)%ptim(ji) > djulini(jj) ) .AND. &
	608	& ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN
	609	jk = jk + 1
	610	ifileidx(jk) = jj
	611	iprofidx(jk) = ji
	612	zdat(jk) = inpfiles(jj)%ptim(ji)
	613	ENDIF
	614	END DO
	615	END DO
	616	CALL sort_dp_indx( iproftot, &
	617	& zdat, &
	618	& iindx )
[6140]	619
[2128]	620	iv3dt(:) = -1
[15799]	621	IF (ldallatall) THEN
	622	iv3dt(:) = ip3dt
[2128]	623	ELSE
[15799]	624	iv3dt(:) = ivart0(:)
[2128]	625	ENDIF
[15799]	626	CALL obs_prof_alloc( profdata, kvars, kadd+iadd, kextr+iextr, iprof, iv3dt, &
	627	& ip3dt, kstp, jpi, jpj, jpk )
[6140]	628
[2128]	629	! * Read obs/positions, QC, all variable and assign to profdata
	630
	631	profdata%nprof = 0
	632	profdata%nvprot(:) = 0
[15799]	633	profdata%cvars(:) = clvarsin(:)
	634	profdata%clong(:) = cllongin(:)
	635	profdata%cunit(:) = clunitin(:)
	636	profdata%cgrid(:) = clgridin(:)
	637	IF ( iadd > 0 ) THEN
	638	profdata%caddvars(kadd+1:) = claddvarsin(:)
	639	profdata%caddlong(kadd+1:,:) = claddlongin(:,:)
	640	profdata%caddunit(kadd+1:,:) = claddunitin(:,:)
	641	ENDIF
	642	IF ( iextr > 0 ) THEN
	643	profdata%cextvars(kextr+1:) = clextvarsin(:)
	644	profdata%cextlong(kextr+1:) = clextlongin(:)
	645	profdata%cextunit(kextr+1:) = clextunitin(:)
	646	ENDIF
[2128]	647	iprof = 0
	648
	649	ip3dt = 0
[15799]	650	ivart(:) = 0
	651	itypvar (:,:) = 0
	652	itypvarmpp(:,:) = 0
[6140]	653
	654	ioserrcount = 0
[2128]	655	DO jk = 1, iproftot
[6140]	656
[2128]	657	jj = ifileidx(iindx(jk))
	658	ji = iprofidx(iindx(jk))
	659
[15799]	660	IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE
	661	llcycle = .TRUE.
	662	DO jvar = 1, kvars
	663	IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN
	664	llcycle = .FALSE.
	665	EXIT
	666	ENDIF
	667	END DO
	668	IF ( llcycle ) CYCLE
[2128]	669
	670	IF ( ( inpfiles(jj)%ptim(ji) > djulini(jj) ) .AND. &
	671	& ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN
[6140]	672
[2128]	673	IF ( nproc == 0 ) THEN
	674	IF ( inpfiles(jj)%iproc(ji,1) > nproc ) CYCLE
	675	ELSE
	676	IF ( inpfiles(jj)%iproc(ji,1) /= nproc ) CYCLE
	677	ENDIF
[6140]	678
[2128]	679	llvalprof = .FALSE.
	680
	681	IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE
	682
[9023]	683	IF ( BTEST(inpfiles(jj)%ioqc(ji),2 ) ) CYCLE
[15799]	684	llcycle = .TRUE.
	685	DO jvar = 1, kvars
	686	IF ( .NOT. ( BTEST(inpfiles(jj)%ivqc(ji,jvar),2) ) ) THEN
	687	llcycle = .FALSE.
	688	EXIT
	689	ENDIF
	690	END DO
	691	IF ( llcycle ) CYCLE
[2128]	692
	693	loop_prof : DO ij = 1, inpfiles(jj)%nlev
[6140]	694
[2128]	695	IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) &
	696	& CYCLE
[6140]	697
[15799]	698	DO jvar = 1, kvars
	699	IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,jvar),2) .AND. &
	700	& .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN
[6140]	701
[15799]	702	llvalprof = .TRUE.
	703	EXIT loop_prof
[6140]	704
[15799]	705	ENDIF
	706	END DO
[6140]	707
[2128]	708	END DO loop_prof
[6140]	709
[2128]	710	! Set profile information
[6140]	711
[2128]	712	IF ( llvalprof ) THEN
[6140]	713
[2128]	714	iprof = iprof + 1
	715
	716	CALL jul2greg( isec, &
	717	& imin, &
	718	& ihou, &
	719	& iday, &
	720	& imon, &
	721	& iyea, &
	722	& inpfiles(jj)%ptim(ji), &
	723	& irefdate(jj) )
	724
	725
	726	! Profile time coordinates
	727	profdata%nyea(iprof) = iyea
	728	profdata%nmon(iprof) = imon
	729	profdata%nday(iprof) = iday
	730	profdata%nhou(iprof) = ihou
	731	profdata%nmin(iprof) = imin
[6140]	732
[2128]	733	! Profile space coordinates
	734	profdata%rlam(iprof) = inpfiles(jj)%plam(ji)
	735	profdata%rphi(iprof) = inpfiles(jj)%pphi(ji)
	736
	737	! Coordinate search parameters
[15799]	738	DO jvar = 1, kvars
	739	profdata%mi (iprof,jvar) = inpfiles(jj)%iobsi(ji,jvar)
	740	profdata%mj (iprof,jvar) = inpfiles(jj)%iobsj(ji,jvar)
	741	END DO
[6140]	742
[2128]	743	! Profile WMO number
	744	profdata%cwmo(iprof) = inpfiles(jj)%cdwmo(ji)
[6140]	745
[2128]	746	! Instrument type
[4990]	747	READ( inpfiles(jj)%cdtyp(ji), '(I4)', IOSTAT = ios, ERR = 901 ) itype
	748	901 IF ( ios /= 0 ) THEN
	749	IF (ioserrcount == 0) CALL ctl_warn ( 'Problem converting an instrument type to integer. Setting type to zero' )
	750	ioserrcount = ioserrcount + 1
	751	itype = 0
	752	ENDIF
[6140]	753
[2128]	754	profdata%ntyp(iprof) = itype
[6140]	755
[2128]	756	! QC stuff
	757
	758	profdata%nqc(iprof) = inpfiles(jj)%ioqc(ji)
	759	profdata%nqcf(:,iprof) = inpfiles(jj)%ioqcf(:,ji)
	760	profdata%ipqc(iprof) = inpfiles(jj)%ipqc(ji)
	761	profdata%ipqcf(:,iprof) = inpfiles(jj)%ipqcf(:,ji)
	762	profdata%itqc(iprof) = inpfiles(jj)%itqc(ji)
	763	profdata%itqcf(:,iprof) = inpfiles(jj)%itqcf(:,ji)
	764	profdata%ivqc(iprof,:) = inpfiles(jj)%ivqc(ji,:)
	765	profdata%ivqcf(:,iprof,:) = inpfiles(jj)%ivqcf(:,ji,:)
	766
	767	! Bookkeeping data to match profiles
	768	profdata%npidx(iprof) = iprof
	769	profdata%npfil(iprof) = iindx(jk)
	770
	771	! Observation QC flag (whole profile)
	772	profdata%nqc(iprof) = 0 !TODO
	773
[6140]	774	loop_p : DO ij = 1, inpfiles(jj)%nlev
	775
[2128]	776	IF ( inpfiles(jj)%pdep(ij,ji) >= 6000. ) &
	777	& CYCLE
	778
[15799]	779	IF ( ldallatall .OR. (iextr > 0) ) THEN
[2128]	780
[15799]	781	DO jvar = 1, kvars
	782	IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,jvar),2) .AND. &
	783	& .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. &
	784	& ldvar(jvar) ) ) THEN
	785	ip3dt = ip3dt + 1
	786	EXIT
	787	ELSE IF ( jvar == kvars ) THEN
	788	CYCLE loop_p
	789	ENDIF
	790	END DO
[6140]	791
[2128]	792	ENDIF
	793
[15799]	794	DO jvar = 1, kvars
	795
	796	IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,jvar),2) .AND. &
	797	& .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. &
	798	& ldvar(jvar) ) .OR. ldallatall ) THEN
[6140]	799
[15799]	800	IF (ldallatall) THEN
[2128]	801
[15799]	802	ivart(jvar) = ip3dt
[2128]	803
[15799]	804	ELSE
[2128]	805
[15799]	806	ivart(jvar) = ivart(jvar) + 1
[6140]	807
[15799]	808	ENDIF
[2128]	809
[15799]	810	! Depth of jvar observation
	811	profdata%var(jvar)%vdep(ivart(jvar)) = &
	812	& inpfiles(jj)%pdep(ij,ji)
[6140]	813
[15799]	814	! Depth of jvar observation QC
	815	profdata%var(jvar)%idqc(ivart(jvar)) = &
	816	& inpfiles(jj)%idqc(ij,ji)
[6140]	817
[15799]	818	! Depth of jvar observation QC flags
	819	profdata%var(jvar)%idqcf(:,ivart(jvar)) = &
	820	& inpfiles(jj)%idqcf(:,ij,ji)
[6140]	821
[15799]	822	! Profile index
	823	profdata%var(jvar)%nvpidx(ivart(jvar)) = iprof
[6140]	824
[15799]	825	! Vertical index in original profile
	826	profdata%var(jvar)%nvlidx(ivart(jvar)) = ij
[2128]	827
[15799]	828	! Profile jvar value
	829	IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,jvar),2) .AND. &
	830	& .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) ) THEN
	831	profdata%var(jvar)%vobs(ivart(jvar)) = &
	832	& inpfiles(jj)%pob(ij,ji,jvar)
	833	! Count number of profile var1 data as function of type
	834	itypvar( profdata%ntyp(iprof) + 1, jvar ) = &
	835	& itypvar( profdata%ntyp(iprof) + 1, jvar ) + 1
	836	ELSE
	837	profdata%var(jvar)%vobs(ivart(jvar)) = fbrmdi
[2128]	838	ENDIF
	839
[15799]	840	! Profile jvar qc
	841	profdata%var(jvar)%nvqc(ivart(jvar)) = &
	842	& inpfiles(jj)%ivlqc(ij,ji,jvar)
[2128]	843
[15799]	844	! Profile jvar qc flags
	845	profdata%var(jvar)%nvqcf(:,ivart(jvar)) = &
	846	& inpfiles(jj)%ivlqcf(:,ij,ji,jvar)
[2128]	847
[15799]	848	! Additional variables
	849	IF ( iadd > 0 ) THEN
	850	jadd2 = 0
	851	DO jadd = 1, inpfiles(jj)%nadd
	852	IF ( TRIM(inpfiles(jj)%caddname(jadd)) == 'Hx' ) THEN
	853	IF ( ldmod ) THEN
	854	profdata%var(jvar)%vmod(ivart(jvar)) = &
	855	& inpfiles(jj)%padd(ij,ji,jadd,jvar)
	856	ENDIF
	857	ELSE
	858	jadd2 = jadd2 + 1
	859	profdata%var(jvar)%vadd(ivart(jvar),kadd+jadd2) = &
	860	& inpfiles(jj)%padd(ij,ji,jadd,jvar)
	861	ENDIF
	862	END DO
	863	ENDIF
[6140]	864
[2128]	865	ENDIF
[15799]	866
	867	END DO
	868
	869	! Extra variables
	870	! Special consideration for if the extra variable is called TEMP
	871	! and there's a regular variable called POTM. These are in situ
	872	! and potential temperature respectively, and need the same QC checks
	873	IF ( iextr > 0 ) THEN
	874	profdata%vext%nepidx(ip3dt) = iprof
	875	profdata%vext%nelidx(ip3dt) = ij
	876	DO jext = 1, iextr
	877	IF ( TRIM(inpfiles(jj)%cextname(jext)) == 'TEMP' ) THEN
	878	llpotm = .false.
	879	DO jvar = 1, kvars
	880	IF ( TRIM(inpfiles(jj)%cname(jvar)) == 'POTM' ) THEN
	881	IF ( ( .NOT. BTEST(inpfiles(jj)%ivlqc(ij,ji,jvar),2) .AND. &
	882	& .NOT. BTEST(inpfiles(jj)%idqc(ij,ji),2) .AND. &
	883	& ldvar(jvar) ) ) THEN
	884	profdata%vext%eobs(ip3dt,kextr+jext) = inpfiles(jj)%pext(ij,ji,jext)
	885	ELSE
	886	profdata%vext%eobs(ip3dt,kextr+jext) = fbrmdi
	887	ENDIF
	888	llpotm = .true.
	889	EXIT
	890	ENDIF
	891	END DO
	892	IF ( .NOT. llpotm ) THEN
	893	profdata%vext%eobs(ip3dt,kextr+jext) = inpfiles(jj)%pext(ij,ji,jext)
	894	ENDIF
	895	ELSE
	896	profdata%vext%eobs(ip3dt,kextr+jext) = inpfiles(jj)%pext(ij,ji,jext)
[2128]	897	ENDIF
[15799]	898	END DO
[2128]	899	ENDIF
[6140]	900
[2128]	901	END DO loop_p
	902
	903	ENDIF
	904
	905	ENDIF
	906
	907	END DO
	908
	909	!-----------------------------------------------------------------------
	910	! Sum up over processors
	911	!-----------------------------------------------------------------------
[6140]	912
[15799]	913	DO jvar = 1, kvars
	914	CALL obs_mpp_sum_integer ( ivart0(jvar), ivartmpp(jvar) )
	915	END DO
[6140]	916	CALL obs_mpp_sum_integer ( ip3dt, ip3dtmpp )
	917
[15799]	918	DO jvar = 1, kvars
	919	CALL obs_mpp_sum_integers( itypvar(:,jvar), itypvarmpp(:,jvar), ntyp1770 + 1 )
	920	END DO
[6140]	921
[2128]	922	!-----------------------------------------------------------------------
	923	! Output number of observations.
	924	!-----------------------------------------------------------------------
	925	IF(lwp) THEN
	926	WRITE(numout,*)
[6140]	927	WRITE(numout,'(A)') ' Profile data'
[2128]	928	WRITE(numout,'(1X,A)') '------------'
	929	WRITE(numout,*)
[15799]	930	DO jvar = 1, kvars
	931	WRITE(numout,'(1X,A)') 'Profile data, '//TRIM( profdata%cvars(jvar) )
	932	WRITE(numout,'(1X,A)') '------------------------'
	933	DO ji = 0, ntyp1770
	934	IF ( itypvarmpp(ji+1,jvar) > 0 ) THEN
	935	WRITE(numout,'(1X,A3,1X,A48,A3,I8)') ctypshort(ji), &
	936	& cwmonam1770(ji)(1:52),' = ', &
	937	& itypvarmpp(ji+1,jvar)
	938	ENDIF
	939	END DO
	940	WRITE(numout,'(1X,A)') &
	941	& '---------------------------------------------------------------'
	942	WRITE(numout,'(1X,A55,I8)') &
	943	& 'Total profile data for variable '//TRIM( profdata%cvars(jvar) )// &
	944	& ' = ', ivartmpp(jvar)
	945	WRITE(numout,'(1X,A)') &
	946	& '---------------------------------------------------------------'
	947	WRITE(numout,*)
[2128]	948	END DO
	949	ENDIF
[6140]	950
[15799]	951	IF (ldallatall) THEN
	952	profdata%nvprot(:) = ip3dt
	953	profdata%nvprotmpp(:) = ip3dtmpp
[2128]	954	ELSE
[15799]	955	DO jvar = 1, kvars
	956	profdata%nvprot(jvar) = ivart(jvar)
	957	profdata%nvprotmpp(jvar) = ivartmpp(jvar)
	958	END DO
[2128]	959	ENDIF
	960	profdata%nprof = iprof
	961
	962	!-----------------------------------------------------------------------
	963	! Model level search
	964	!-----------------------------------------------------------------------
[15799]	965	DO jvar = 1, kvars
	966	IF ( ldvar(jvar) ) THEN
	967	CALL obs_level_search( jpk, gdept_1d, &
	968	& profdata%nvprot(jvar), profdata%var(jvar)%vdep, &
	969	& profdata%var(jvar)%mvk )
	970	ENDIF
	971	END DO
[6140]	972
[2128]	973	!-----------------------------------------------------------------------
	974	! Set model equivalent to 99999
	975	!-----------------------------------------------------------------------
	976	IF ( .NOT. ldmod ) THEN
	977	DO jvar = 1, kvars
	978	profdata%var(jvar)%vmod(:) = fbrmdi
	979	END DO
	980	ENDIF
	981	!-----------------------------------------------------------------------
	982	! Deallocate temporary data
	983	!-----------------------------------------------------------------------
[15799]	984	DEALLOCATE( ifileidx, iprofidx, zdat, &
	985	& clvarsin, cllongin, clunitin, clgridin )
	986	IF ( iadd > 0 ) THEN
	987	DEALLOCATE( claddvarsin, claddlongin, claddunitin)
	988	ENDIF
	989	IF ( iextr > 0 ) THEN
	990	DEALLOCATE( clextvarsin, clextlongin, clextunitin )
	991	ENDIF
[2128]	992
	993	!-----------------------------------------------------------------------
	994	! Deallocate input data
	995	!-----------------------------------------------------------------------
	996	DO jj = 1, inobf
	997	CALL dealloc_obfbdata( inpfiles(jj) )
	998	END DO
	999	DEALLOCATE( inpfiles )
	1000
[6140]	1001	END SUBROUTINE obs_rea_prof
[2128]	1002
	1003	END MODULE obs_read_prof

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

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

Context Navigation

source: NEMO/branches/UKMO/NEMO_4.0.4_FOAM_package/src/OCE/OBS/obs_read_prof.F90

Download in other formats: