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_surf.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_read_surf.F90 @ 14644

Last change on this file since 14644 was 14644, checked in by sparonuz, 3 years ago
Merge trunk -r14642:HEAD
Property svn:keywords set to `Id`
File size: 19.2 KB

Rev	Line
[5659]	1	MODULE obs_read_surf
	2	!!======================================================================
	3	!! * MODULE obs_read_surf *
	4	!! Observation diagnostics: Read the surface data from feedback files
	5	!!======================================================================
	6
	7	!!----------------------------------------------------------------------
	8	!! obs_rea_surf : Driver for reading surface data from feedback files
	9	!!----------------------------------------------------------------------
	10
[5682]	11	!! * Modules used
[5659]	12	USE par_kind ! Precision variables
	13	USE in_out_manager ! I/O manager
	14	USE dom_oce ! Ocean space and time domain variables
	15	USE obs_mpp ! MPP support routines for observation diagnostics
	16	USE julian ! Julian date routines
	17	USE obs_utils ! Observation operator utility functions
	18	USE obs_grid ! Grid search
	19	USE obs_sort ! Sorting observation arrays
	20	USE obs_surf_def ! Surface observation definitions
	21	USE obs_types ! Observation type definitions
[5682]	22	USE obs_fbm ! Feedback routines
[5659]	23	USE netcdf ! NetCDF library
	24
	25	IMPLICIT NONE
	26
	27	!! * Routine accessibility
	28	PRIVATE
	29
	30	PUBLIC obs_rea_surf ! Read the surface observations from the point data
[5682]	31
[5659]	32	!!----------------------------------------------------------------------
[9598]	33	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[10069]	34	!! $Id$
[10068]	35	!! Software governed by the CeCILL license (see ./LICENSE)
[5659]	36	!!----------------------------------------------------------------------
	37
	38	CONTAINS
	39
[5682]	40	SUBROUTINE obs_rea_surf( surfdata, knumfiles, cdfilenames, &
[5659]	41	& kvars, kextr, kstp, ddobsini, ddobsend, &
[14056]	42	& ldignmis, ldmod, ldnightav, cdvars )
[5659]	43	!!---------------------------------------------------------------------
	44	!!
	45	!! * ROUTINE obs_rea_surf *
	46	!!
	47	!! ** Purpose : Read from file the surface data
	48	!!
	49	!! ** Method : Read in the data from feedback format files and
	50	!! put into the NEMO internal surface data structure
	51	!!
	52	!! ** Action :
	53	!!
	54	!!
	55	!! History :
	56	!! ! : 2009-01 (K. Mogensen) Initial version based on old versions
	57	!! ! : 2015-02 (M. Martin) Unify the different surface data type reading.
	58	!!----------------------------------------------------------------------
	59	!! * Modules used
	60
	61	!! * Arguments
[5682]	62	TYPE(obs_surf), INTENT(INOUT) :: &
	63	& surfdata ! Surface data to be read
	64	INTEGER, INTENT(IN) :: knumfiles ! Number of corio format files to read
	65	CHARACTER(LEN=128), INTENT(IN) :: &
	66	& cdfilenames(knumfiles) ! File names to read in
[5659]	67	INTEGER, INTENT(IN) :: kvars ! Number of variables in surfdata
[5682]	68	INTEGER, INTENT(IN) :: kextr ! Number of extra fields for each var
[5659]	69	INTEGER, INTENT(IN) :: kstp ! Ocean time-step index
	70	LOGICAL, INTENT(IN) :: ldignmis ! Ignore missing files
	71	LOGICAL, INTENT(IN) :: ldmod ! Initialize model from input data
[5682]	72	LOGICAL, INTENT(IN) :: ldnightav ! Observations represent a night-time average
	73	REAL(dp), INTENT(IN) :: ddobsini ! Obs. ini time in YYYYMMDD.HHMMSS
	74	REAL(dp), INTENT(IN) :: ddobsend ! Obs. end time in YYYYMMDD.HHMMSS
[14056]	75	CHARACTER(len=8), DIMENSION(kvars), INTENT(IN) :: cdvars
[5682]	76
[5659]	77	!! * Local declarations
	78	CHARACTER(LEN=11), PARAMETER :: cpname='obs_rea_surf'
[5682]	79	CHARACTER(len=8) :: clrefdate
[14056]	80	CHARACTER(len=8), DIMENSION(:), ALLOCATABLE :: clvarsin
[5659]	81	INTEGER :: ji
	82	INTEGER :: jj
	83	INTEGER :: jk
	84	INTEGER :: iflag
	85	INTEGER :: inobf
	86	INTEGER :: i_file_id
	87	INTEGER :: inowin
	88	INTEGER :: iyea
	89	INTEGER :: imon
	90	INTEGER :: iday
	91	INTEGER :: ihou
	92	INTEGER :: imin
	93	INTEGER :: isec
[5682]	94	INTEGER :: itype
	95	INTEGER :: iobsmpp
	96	INTEGER :: iobs
	97	INTEGER :: iobstot
	98	INTEGER :: ios
	99	INTEGER :: ioserrcount
	100	INTEGER, PARAMETER :: jpsurfmaxtype = 1024
[5659]	101	INTEGER, DIMENSION(knumfiles) :: irefdate
[5682]	102	INTEGER, DIMENSION(jpsurfmaxtype+1) :: &
[5659]	103	& ityp, &
	104	& itypmpp
	105	INTEGER, DIMENSION(:), ALLOCATABLE :: &
	106	& iobsi, &
	107	& iobsj, &
	108	& iproc, &
	109	& iindx, &
	110	& ifileidx, &
	111	& isurfidx
	112	REAL(wp), DIMENSION(:), ALLOCATABLE :: &
	113	& zphi, &
	114	& zlam
[13226]	115	REAL(dp), DIMENSION(:), ALLOCATABLE :: &
[5659]	116	& zdat
[13226]	117	REAL(dp), DIMENSION(knumfiles) :: &
[5682]	118	& djulini, &
	119	& djulend
[5659]	120	LOGICAL :: llvalprof
	121	TYPE(obfbdata), POINTER, DIMENSION(:) :: &
	122	& inpfiles
[5682]	123
[5659]	124	! Local initialization
	125	iobs = 0
	126
	127	!-----------------------------------------------------------------------
	128	! Count the number of files needed and allocate the obfbdata type
	129	!-----------------------------------------------------------------------
[5682]	130
[5659]	131	inobf = knumfiles
[5682]	132
[5659]	133	ALLOCATE( inpfiles(inobf) )
	134
	135	surf_files : DO jj = 1, inobf
	136
	137	!---------------------------------------------------------------------
	138	! Prints
	139	!---------------------------------------------------------------------
	140	IF(lwp) THEN
	141	WRITE(numout,*)
	142	WRITE(numout,*) ' obs_rea_surf : Reading from file = ', &
[5682]	143	& TRIM( TRIM( cdfilenames(jj) ) )
[5659]	144	WRITE(numout,*) ' ~~~~~~~~~~~'
	145	WRITE(numout,*)
	146	ENDIF
	147
	148	!---------------------------------------------------------------------
	149	! Initialization: Open file and get dimensions only
	150	!---------------------------------------------------------------------
[5682]	151
	152	iflag = nf90_open( TRIM( TRIM( cdfilenames(jj) ) ), nf90_nowrite, &
[5659]	153	& i_file_id )
[5682]	154
[5659]	155	IF ( iflag /= nf90_noerr ) THEN
	156
	157	IF ( ldignmis ) THEN
	158	inpfiles(jj)%nobs = 0
[5682]	159	CALL ctl_warn( 'File ' // TRIM( TRIM( cdfilenames(jj) ) ) // &
[5659]	160	& ' not found' )
	161	ELSE
[5682]	162	CALL ctl_stop( 'File ' // TRIM( TRIM( cdfilenames(jj) ) ) // &
[5659]	163	& ' not found' )
	164	ENDIF
	165
	166	ELSE
[5682]	167
[5659]	168	!------------------------------------------------------------------
[5682]	169	! Close the file since it is opened in read_obfbdata
[5659]	170	!------------------------------------------------------------------
[5682]	171
[5659]	172	iflag = nf90_close( i_file_id )
	173
	174	!------------------------------------------------------------------
[9023]	175	! Read the surface file into inpfiles
[5659]	176	!------------------------------------------------------------------
[5682]	177	CALL init_obfbdata( inpfiles(jj) )
	178	CALL read_obfbdata( TRIM( cdfilenames(jj) ), inpfiles(jj), &
[5659]	179	& ldgrid = .TRUE. )
[5682]	180
[14056]	181	IF ( inpfiles(jj)%nvar /= kvars ) THEN
	182	CALL ctl_stop( 'Feedback format error: ', &
	183	& ' unexpected number of vars in feedback file' )
	184	ENDIF
	185
[5659]	186	IF ( ldmod .AND. ( inpfiles(jj)%nadd == 0 ) ) THEN
	187	CALL ctl_stop( 'Model not in input data' )
	188	RETURN
	189	ENDIF
	190
[5682]	191	IF ( jj == 1 ) THEN
[14056]	192	ALLOCATE( clvarsin( inpfiles(jj)%nvar ) )
[5682]	193	DO ji = 1, inpfiles(jj)%nvar
[14056]	194	clvarsin(ji) = inpfiles(jj)%cname(ji)
	195	IF ( clvarsin(ji) /= cdvars(ji) ) THEN
	196	CALL ctl_stop( 'Feedback file variables do not match', &
	197	& ' expected variable names for this type' )
	198	ENDIF
[5682]	199	END DO
	200	ELSE
	201	DO ji = 1, inpfiles(jj)%nvar
[14056]	202	IF ( inpfiles(jj)%cname(ji) /= clvarsin(ji) ) THEN
[5682]	203	CALL ctl_stop( 'Feedback file variables not consistent', &
	204	& ' with previous files for this type' )
	205	ENDIF
	206	END DO
	207	ENDIF
	208
[5704]	209	IF (lwp) WRITE(numout,*)'Observation file contains ',inpfiles(jj)%nobs,' observations'
	210
[5659]	211	!------------------------------------------------------------------
	212	! Change longitude (-180,180)
	213	!------------------------------------------------------------------
	214
[5682]	215	DO ji = 1, inpfiles(jj)%nobs
[5659]	216
	217	IF ( inpfiles(jj)%plam(ji) < -180. ) &
	218	& inpfiles(jj)%plam(ji) = inpfiles(jj)%plam(ji) + 360.
	219
	220	IF ( inpfiles(jj)%plam(ji) > 180. ) &
	221	& inpfiles(jj)%plam(ji) = inpfiles(jj)%plam(ji) - 360.
	222
	223	END DO
	224
	225	!------------------------------------------------------------------
	226	! Calculate the date (change eventually)
	227	!------------------------------------------------------------------
[5682]	228	clrefdate=inpfiles(jj)%cdjuldref(1:8)
	229	READ(clrefdate,'(I8)') irefdate(jj)
	230
[5659]	231	CALL ddatetoymdhms( ddobsini, iyea, imon, iday, ihou, imin, isec )
	232	CALL greg2jul( isec, imin, ihou, iday, imon, iyea, djulini(jj), &
	233	& krefdate = irefdate(jj) )
	234	CALL ddatetoymdhms( ddobsend, iyea, imon, iday, ihou, imin, isec )
	235	CALL greg2jul( isec, imin, ihou, iday, imon, iyea, djulend(jj), &
	236	& krefdate = irefdate(jj) )
[5682]	237
	238	IF ( ldnightav ) THEN
	239
	240	IF ( lwp ) THEN
	241	WRITE(numout,*)'Resetting time of night-time averaged observations', &
	242	& ' to the end of the day'
	243	ENDIF
	244
	245	DO ji = 1, inpfiles(jj)%nobs
	246	! for night-time averaged data force the time
	247	! to be the last time-step of the day, but still within the day.
	248	IF ( inpfiles(jj)%ptim(ji) >= 0. ) THEN
	249	inpfiles(jj)%ptim(ji) = &
	250	& INT(inpfiles(jj)%ptim(ji)) + 0.9999
	251	ELSE
	252	inpfiles(jj)%ptim(ji) = &
	253	& INT(inpfiles(jj)%ptim(ji)) - 0.0001
	254	ENDIF
	255	END DO
	256	ENDIF
	257
[5659]	258	IF ( inpfiles(jj)%nobs > 0 ) THEN
	259	inpfiles(jj)%iproc = -1
	260	inpfiles(jj)%iobsi = -1
	261	inpfiles(jj)%iobsj = -1
	262	ENDIF
	263	inowin = 0
	264	DO ji = 1, inpfiles(jj)%nobs
	265	IF ( ( inpfiles(jj)%ptim(ji) > djulini(jj) ) .AND. &
	266	& ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN
	267	inowin = inowin + 1
	268	ENDIF
	269	END DO
	270	ALLOCATE( zlam(inowin) )
	271	ALLOCATE( zphi(inowin) )
	272	ALLOCATE( iobsi(inowin) )
	273	ALLOCATE( iobsj(inowin) )
	274	ALLOCATE( iproc(inowin) )
	275	inowin = 0
	276	DO ji = 1, inpfiles(jj)%nobs
	277	IF ( ( inpfiles(jj)%ptim(ji) > djulini(jj) ) .AND. &
	278	& ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN
	279	inowin = inowin + 1
	280	zlam(inowin) = inpfiles(jj)%plam(ji)
	281	zphi(inowin) = inpfiles(jj)%pphi(ji)
	282	ENDIF
	283	END DO
	284
	285	CALL obs_grid_search( inowin, zlam, zphi, iobsi, iobsj, iproc, 'T' )
	286
	287	inowin = 0
	288	DO ji = 1, inpfiles(jj)%nobs
	289	IF ( ( inpfiles(jj)%ptim(ji) > djulini(jj) ) .AND. &
	290	& ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN
	291	inowin = inowin + 1
	292	inpfiles(jj)%iproc(ji,1) = iproc(inowin)
	293	inpfiles(jj)%iobsi(ji,1) = iobsi(inowin)
	294	inpfiles(jj)%iobsj(ji,1) = iobsj(inowin)
	295	ENDIF
	296	END DO
	297	DEALLOCATE( zlam, zphi, iobsi, iobsj, iproc )
	298
	299	DO ji = 1, inpfiles(jj)%nobs
	300	IF ( ( inpfiles(jj)%ptim(ji) > djulini(jj) ) .AND. &
	301	& ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN
[14644]	302	IF ( narea == 1 ) THEN
	303	IF ( inpfiles(jj)%iproc(ji,1) > narea-1 ) CYCLE
[5659]	304	ELSE
[14644]	305	IF ( inpfiles(jj)%iproc(ji,1) /= narea-1 ) CYCLE
[5659]	306	ENDIF
	307	llvalprof = .FALSE.
[9023]	308	IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(1,ji,1),2) ) THEN
[5659]	309	iobs = iobs + 1
	310	ENDIF
	311	ENDIF
	312	END DO
	313
	314	ENDIF
	315
	316	END DO surf_files
	317
	318	!-----------------------------------------------------------------------
	319	! Get the time ordered indices of the input data
	320	!-----------------------------------------------------------------------
	321
	322	!---------------------------------------------------------------------
	323	! Loop over input data files to count total number of profiles
	324	!---------------------------------------------------------------------
	325	iobstot = 0
	326	DO jj = 1, inobf
	327	DO ji = 1, inpfiles(jj)%nobs
	328	IF ( ( inpfiles(jj)%ptim(ji) > djulini(jj) ) .AND. &
	329	& ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN
	330	iobstot = iobstot + 1
	331	ENDIF
	332	END DO
	333	END DO
	334
	335	ALLOCATE( iindx(iobstot), ifileidx(iobstot), &
	336	& isurfidx(iobstot), zdat(iobstot) )
	337	jk = 0
	338	DO jj = 1, inobf
	339	DO ji = 1, inpfiles(jj)%nobs
	340	IF ( ( inpfiles(jj)%ptim(ji) > djulini(jj) ) .AND. &
	341	& ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN
	342	jk = jk + 1
	343	ifileidx(jk) = jj
	344	isurfidx(jk) = ji
	345	zdat(jk) = inpfiles(jj)%ptim(ji)
	346	ENDIF
	347	END DO
	348	END DO
	349	CALL sort_dp_indx( iobstot, &
	350	& zdat, &
	351	& iindx )
[5682]	352
[5659]	353	CALL obs_surf_alloc( surfdata, iobs, kvars, kextr, kstp, jpi, jpj )
[5682]	354
	355	! Read obs/positions, QC, all variable and assign to surfdata
	356
[5659]	357	iobs = 0
	358
[14056]	359	surfdata%cvars(:) = clvarsin(:)
[5682]	360
[5659]	361	ityp (:) = 0
	362	itypmpp(:) = 0
[5682]	363
[5659]	364	ioserrcount = 0
[5682]	365
[5659]	366	DO jk = 1, iobstot
[5682]	367
[5659]	368	jj = ifileidx(iindx(jk))
	369	ji = isurfidx(iindx(jk))
	370	IF ( ( inpfiles(jj)%ptim(ji) > djulini(jj) ) .AND. &
	371	& ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN
[5682]	372
[14644]	373	IF ( narea == 1 ) THEN
	374	IF ( inpfiles(jj)%iproc(ji,1) > narea-1 ) CYCLE
[5659]	375	ELSE
[14644]	376	IF ( inpfiles(jj)%iproc(ji,1) /= narea-1 ) CYCLE
[5659]	377	ENDIF
[5682]	378
[5659]	379	! Set observation information
[5682]	380
[9023]	381	IF ( .NOT. BTEST(inpfiles(jj)%ivlqc(1,ji,1),2) ) THEN
[5659]	382
	383	iobs = iobs + 1
	384
	385	CALL jul2greg( isec, &
	386	& imin, &
	387	& ihou, &
	388	& iday, &
	389	& imon, &
	390	& iyea, &
	391	& inpfiles(jj)%ptim(ji), &
	392	& irefdate(jj) )
	393
	394
	395	! Surface time coordinates
	396	surfdata%nyea(iobs) = iyea
	397	surfdata%nmon(iobs) = imon
	398	surfdata%nday(iobs) = iday
	399	surfdata%nhou(iobs) = ihou
	400	surfdata%nmin(iobs) = imin
[5682]	401
[5659]	402	! Surface space coordinates
	403	surfdata%rlam(iobs) = inpfiles(jj)%plam(ji)
	404	surfdata%rphi(iobs) = inpfiles(jj)%pphi(ji)
	405
	406	! Coordinate search parameters
	407	surfdata%mi (iobs) = inpfiles(jj)%iobsi(ji,1)
	408	surfdata%mj (iobs) = inpfiles(jj)%iobsj(ji,1)
[5682]	409
[5704]	410	! WMO number
	411	surfdata%cwmo(iobs) = inpfiles(jj)%cdwmo(ji)
	412
[5659]	413	! Instrument type
	414	READ( inpfiles(jj)%cdtyp(ji), '(I4)', IOSTAT = ios, ERR = 901 ) itype
	415	901 IF ( ios /= 0 ) THEN
[5682]	416	IF (ioserrcount == 0) THEN
	417	CALL ctl_warn ( 'Problem converting an instrument type ', &
	418	& 'to integer. Setting type to zero' )
	419	ENDIF
[5659]	420	ioserrcount = ioserrcount + 1
	421	itype = 0
	422	ENDIF
	423	surfdata%ntyp(iobs) = itype
[5682]	424	IF ( itype < jpsurfmaxtype + 1 ) THEN
[5659]	425	ityp(itype+1) = ityp(itype+1) + 1
	426	ELSE
[5682]	427	IF(lwp)WRITE(numout,*)'WARNING:Increase jpsurfmaxtype in ',&
[5659]	428	& cpname
	429	ENDIF
	430
	431	! Bookkeeping data to match observations
	432	surfdata%nsidx(iobs) = iobs
	433	surfdata%nsfil(iobs) = iindx(jk)
	434
	435	! QC flags
	436	surfdata%nqc(iobs) = inpfiles(jj)%ivqc(ji,1)
	437
	438	! Observed value
	439	surfdata%robs(iobs,1) = inpfiles(jj)%pob(1,ji,1)
	440
	441
	442	! Model and MDT is set to fbrmdi unless read from file
	443	IF ( ldmod ) THEN
	444	surfdata%rmod(iobs,1) = inpfiles(jj)%padd(1,ji,1,1)
[5682]	445	IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) THEN
	446	surfdata%rext(iobs,1) = inpfiles(jj)%padd(1,ji,2,1)
	447	surfdata%rext(iobs,2) = inpfiles(jj)%pext(1,ji,1)
	448	ENDIF
	449	ELSE
[5659]	450	surfdata%rmod(iobs,1) = fbrmdi
[5682]	451	IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) surfdata%rext(iobs,:) = fbrmdi
[5659]	452	ENDIF
	453	ENDIF
	454	ENDIF
	455
	456	END DO
	457
	458	!-----------------------------------------------------------------------
	459	! Sum up over processors
	460	!-----------------------------------------------------------------------
[5682]	461
[5659]	462	CALL obs_mpp_sum_integer( iobs, iobsmpp )
[5682]	463	CALL obs_mpp_sum_integers( ityp, itypmpp, jpsurfmaxtype + 1 )
	464
[5659]	465	!-----------------------------------------------------------------------
	466	! Output number of observations.
	467	!-----------------------------------------------------------------------
	468	IF (lwp) THEN
	469
	470	WRITE(numout,*)
[5682]	471	WRITE(numout,'(1X,A)')TRIM( surfdata%cvars(1) )//' data'
[5659]	472	WRITE(numout,'(1X,A)')'--------------'
	473	DO jj = 1,8
	474	IF ( itypmpp(jj) > 0 ) THEN
	475	WRITE(numout,'(1X,A4,I4,A3,I10)')'Type ', jj,' = ',itypmpp(jj)
	476	ENDIF
	477	END DO
[5682]	478	WRITE(numout,'(1X,A)') &
	479	& '---------------------------------------------------------------'
	480	WRITE(numout,'(1X,A,I8)') &
	481	& 'Total data for variable '//TRIM( surfdata%cvars(1) )// &
	482	& ' = ', iobsmpp
	483	WRITE(numout,'(1X,A)') &
	484	& '---------------------------------------------------------------'
[5659]	485	WRITE(numout,*)
	486
	487	ENDIF
	488
	489	!-----------------------------------------------------------------------
	490	! Deallocate temporary data
	491	!-----------------------------------------------------------------------
[14056]	492	DEALLOCATE( ifileidx, isurfidx, zdat, clvarsin )
[5659]	493
	494	!-----------------------------------------------------------------------
	495	! Deallocate input data
	496	!-----------------------------------------------------------------------
	497	DO jj = 1, inobf
	498	IF ( inpfiles(jj)%lalloc ) THEN
	499	CALL dealloc_obfbdata( inpfiles(jj) )
	500	ENDIF
	501	END DO
	502	DEALLOCATE( inpfiles )
	503
	504	END SUBROUTINE obs_rea_surf
	505
	506	END MODULE obs_read_surf

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

Download in other formats: