Blame | Last modification | View Log | Download | RSS feed
program contrack! ----------------------------------------------------------------------!! Date: 1st January 2011!!! ------->> Purpose: Spatial and temporal contour tracking programm!! This programm is a completely re-written new version of the tracking! programm used for the studies in e.g. Schwierz et al. (2005),! Croci-Maspoli et al. (2007a,b), Croci-Maspoli et al. (2009).!! ------->> What is necessary to compile contrack?! - contrack has been compiled with gfortran! (included in the GCC family, http://gcc.gnu.org/ ). Feel free! to test other fortran compilers.! - contrack requires the netcdf-4.0.1 library (also compiled with gfortran).! - contrack comes with a simple make-file (you have to adjust this file! before compiling).!! ------->> What netCDF-file is necessary as input?!! - a single CF-netCDF-file (http://cf-pcmdi.llnl.gov/) including all required! time steps.! - a longitude, latitude and time dimension.! - contrack does not calculate a climatological anomaly field. This needs to! be calculated before running contrack.!! ------->> How does the output file look like?!! - the output CF-netCDF file has the same dimensions as the input netCDF file! - the output has a binary format (0 and 1), whereas 1 captures regions of! the blocks.!! ------->> Known issues or features not implemented yet in contrack 0.8:!! 1) The calculation and saving of common statistical blocking features! (e.g. blocking size, blocking length, blocking path, ...) is not! implemented yet. My intention is to save these information directly! on the netCDF file.! 2) Only two-dimensional netCDF-files can be used as input file (longitude,! latitude and time). The code can not handle a third dimension (e.g.! vertical) on the same netCDF-file yet.! 3) The asymmetry of the blocking features is not implemented yet.! 3) ... more to come!!! ------->> Version control! Version 0.1-0.8: May 2010 - fortran 90 compatible! - CF-netCDF compatible (instead of IVE-netCDF)! - new tracking algorithm! Oct 2011 - started with handling also different vertical! levels but not finished yet.! Jan 2011 - added these explanations!! If you have any comments please send them to mischa.croci-maspoli@alumni.ethz.ch.!! This programm code comes with no guarantee and is currently still! under construction. Be aware that no extensive testing and comparison to the! old code has been performed yet.!! Copyright: mischa.croci-maspoli@alumni.ethz.ch!! Modifications SP: omit first and last n time steps in output,! where n=persistence;! comment out allocation etc. related to vertical dimension!! ----------------------------------------------------------------------USE netcdfimplicit noneinteger, external :: iargcinteger :: i,j,t,k,ii,jj,tt,counterinteger :: persistence, verbose, calcmeranom, continteger :: titleLength, commentLength, sourceLength, institLengthinteger :: errmsg, outarrinteger :: nrcontoursinteger :: ttbegin, ttendreal :: fraction, areacon, areaover, overlap, thresholdcharacter*(100) infilename,outfilename,argcharacter*(100) outvarname,outstandvarname, outvarunitcharacter*(100) invarname, inlat, inlon, intime, invercharacter*(100) contrackversioncharacter*(5) gorlcharacter*(1000) cfcomment, cftitle, cfsource, cfinstitcharacter*(1000) cfoutcomment, cfouttitle, cfoutsource, cfoutinstitinteger :: ncid, status, nDims, nVars, nGlobalAtts, unlimDimIDinteger :: label, nrsmooth, nrg, nrgcinteger :: LatDimID,LonDimID, TimeDimID,timeVarID,lonVarID,latVarID,vorpotVarIDinteger :: VerDimIDinteger :: varLatID,varLonID,varPVID,varTimeIDinteger :: nLats,nLong,nTimes,nVertinteger, dimension(:), allocatable :: times, latitude, longitude, longiinteger, dimension(:), allocatable :: verticalinteger, dimension(:), allocatable :: iiv, jjv, ttvinteger, dimension(:), allocatable :: time_arrayreal, dimension(:,:,:), allocatable :: inarr, arr, arr1, arr2real, dimension(:,:), allocatable :: latmean! --------------------------------------------------------------! input handling! --------------------------------------------------------------contrackversion = 'version 0.8'status = 0print*, '-------------------------------------------------------'print*, 'contrack ', trim(contrackversion)print*, '-------------------------------------------------------'999 continueif ( (iargc().ne.1).or.(status.gt.0) ) thenprint*, '-------------------------------------------------------'print*, 'contrack ', trim(contrackversion)print*, '-------------------------------------------------------'print*, 'Usage: contrack inputfile'print*, '-------------------------------------------------------'print*, 'The inputfile needs the following structure:'print*, ''print*, 'infilename -> name of the input CF-netCDF file'print*, 'invarname -> name of the netCDF variable name'print*, 'inlat -> name of the netCDF latitude variable'print*, 'inlon -> name of the netCDF longitude variable'print*, 'intime -> name of the netCDF time variable'print*, 'threshold -> threshold value to detect contours'print*, 'gorl -> find contours that are greater or 'print*, ' lower threshold value [ge,le,gt,lt]'print*, 'persistence -> temporal persistence (in time steps)'print*, ' of the contour life time'print*, 'overlap -> overlapping fraction of two contours'print*, ' between two time steps [0-1]'print*, 'nrsmooth -> temporal smoothing of the infilename'print*, ' (in time steps)'print*, 'outarr -> values of the contours in the output array'print*, ' 0: contours are numbered by each contour'print*, ' 1: contours are labeled 0 and 1'print*, 'outfilename -> name of the output CF-netCDF file'print*, 'outvarname -> name of the output netCDF variable name'print*, 'outstvarname -> name of the output netCDF standard'print*, ' variable name'print*, 'outvarunit -> name of the outuput netCDF unit'print*, 'verbose -> [0 or 1]'print*, 'calcmeranom -> flag whether to calc meriodonal anomalies'print*, ' or not [0 or 1]'print*, '-------------------------------------------------------'print*, 'Please send any comments / bugs to:'print*, 'mischa.croci-maspoli@alumni.ethz.ch'print*, '-------------------------------------------------------'print*, ''print*, '-->> contrack not executed, an error occured!'print*, ''call exit(1)endif! Read inputfile! --------------call getarg(1,arg)infilename=trim(arg)! Open inputfile and read variables! ---------------------------------open(2,file=infilename,status='old',iostat=status)if ( status.gt.0 ) thenprint*, ''print*, '-> Problem in contrack: Your ASCII input file "', trim(infilename), &'" can not be found'print*, '-> For further information use ./contrack only'print*, ''stopendifread(2,*,iostat=status) infilenameif ( status.gt.0 ) goto 999read(2,*,iostat=status) invarnameif ( status.gt.0 ) goto 999read(2,*,iostat=status) inlatif ( status.gt.0 ) goto 999read(2,*,iostat=status) inlonif ( status.gt.0 ) goto 999read(2,*,iostat=status) intimeif ( status.gt.0 ) goto 999read(2,*,iostat=status) thresholdif ( status.gt.0 ) goto 999read(2,*,iostat=status) gorlif ( status.gt.0 ) goto 999read(2,*,iostat=status) persistenceif ( status.gt.0 ) goto 999read(2,*,iostat=status) overlapif ( status.gt.0 ) goto 999read(2,*,iostat=status) nrsmoothif ( status.gt.0 ) goto 999read(2,*,iostat=status) outarrif ( status.gt.0 ) goto 999read(2,*,iostat=status) outfilenameif ( status.gt.0 ) goto 999read(2,*,iostat=status) outvarnameif ( status.gt.0 ) goto 999read(2,*,iostat=status) outstandvarnameif ( status.gt.0 ) goto 999read(2,*,iostat=status) outvarunitif ( status.gt.0 ) goto 999read(2,*,iostat=status) verboseif ( status.gt.0 ) goto 999read(2,*,iostat=status) calcmeranomif ( status.gt.0 ) goto 999close(2)! --------------------------------------------------------------! read CF-netCDF file! --------------------------------------------------------------allocate(time_array(8))call date_and_time (values=time_array)write(*,100), "-> Read ", trim(infilename), time_array(5), ':', time_array(6), ':', time_array(7)100 format(A,A,I5,A,I2,A,I2)! Initially set error to indicate no errors.status = 0! open netCDF filestatus = nf90_open(infilename, nf90_nowrite, ncid)! get dimensionsstatus = nf90_inq_dimid(ncid, inlat, LatDimID)status = nf90_inquire_dimension(ncid, LatDimID, len = nLats)status = nf90_inq_dimid(ncid, inlon, LonDimID)status = nf90_inquire_dimension(ncid, LonDimID, len = nLong)!status = nf90_inq_dimid(ncid, inver, VerDimID)!status = nf90_inquire_dimension(ncid, VerDimID, len = nVert)status = nf90_inq_dimid(ncid, intime, TimeDimID)status = nf90_inquire_dimension(ncid, TimeDimID, len = nTimes)! allocate arraysallocate(times(nTimes))allocate(latitude(nLats))allocate(longitude(nLong))!allocate(vertical(nVert))allocate(inarr(nLong,nLats,nTimes))allocate(arr(nLong,nLats,nTimes))allocate(arr1(nLong,nLats,nTimes))allocate(arr2(nLong,nLats,nTimes))allocate(latmean(nLats,nTimes))allocate(iiv(nLats*nLong*nTimes))allocate(jjv(nLats*nLong*nTimes))allocate(ttv(nLats*nLong*nTimes))! get variablesstatus = nf90_inq_varid(ncid, intime, timeVarID)status = nf90_get_var(ncid, timeVarID, times)status = nf90_inq_varid(ncid, inlon, lonVarID)status = nf90_get_var(ncid, lonVarID, longitude)status = nf90_inq_varid(ncid, inlat, latVarID)status = nf90_get_var(ncid, latVarID, latitude)status = nf90_inq_varid(ncid, invarname, vorpotVarID)status = nf90_get_var(ncid, vorpotVarID, inarr)status = nf90_inquire_attribute(ncid, nf90_global, "comment", len = commentLength)status = nf90_get_att(ncid, NF90_GLOBAL, 'comment', cfcomment)status = nf90_inquire_attribute(ncid, nf90_global, "title", len = titleLength)status = nf90_get_att(ncid, NF90_GLOBAL, 'title', cftitle)status = nf90_inquire_attribute(ncid, nf90_global, "source", len = sourceLength)status = nf90_get_att(ncid, NF90_GLOBAL, 'source', cfsource)status = nf90_inquire_attribute(ncid, nf90_global, "institution", len = institLength)status = nf90_get_att(ncid, NF90_GLOBAL, 'institution', cfinstit)status = nf90_close(ncid)! --------------------------------------------------------------! Calculate meridional anomaly! --------------------------------------------------------------if ( calcmeranom.eq.1 ) thenprint*, '-->> calculate meridional mean anomaly'! Calculate latitudinal mean! --------------------------------------------------------------do t = 1,nTimesdo j = 1,nLatsdo i = 1,nLonglatmean(j,t) = latmean(j,t) + inarr(i,j,t)enddoenddoenddo! Subtract latitudinal mean! --------------------------------------------------------------do t = 1,nTimesdo j = 1,nLatsdo i = 1,nLongarr(i,j,t) = inarr(i,j,t) - latmean(j,t)/nLongenddoenddoenddo! Calculate running means (centered at the first file)! --------------------------------------------------------------do t = 1,nTimesdo j = 1,nLatsdo i = 1,nLongdo k = 0,nrsmootharr1(i,j,t) = arr1(i,j,t) + arr(i,j,t+k)enddoenddoenddoenddoarr = arr1/nrsmoothelseprint*, '-->> no anomaly calculated'arr(:,:,:) = inarr(:,:,:)endif! Define closed contours by a threshold value! --------------------------------------------------------------print*, '-->> define contours with threshold value'do t = 1,nTimesdo j = 1,nLatsdo i = 1,nLongif ( trim(gorl).eq.'ge' ) thenif ( arr(i,j,t).ge.threshold) thenarr(i,j,t) = 1elsearr(i,j,t) = 0endifendifif ( trim(gorl).eq.'le' ) thenif ( arr(i,j,t).le.threshold) thenarr(i,j,t) = 1elsearr(i,j,t) = 0endifendifif ( trim(gorl).eq.'gt' ) thenif ( arr(i,j,t).gt.threshold) thenarr(i,j,t) = 1elsearr(i,j,t) = 0endifendifif ( trim(gorl).eq.'lt' ) thenif ( arr(i,j,t).lt.threshold) thenarr(i,j,t) = 1elsearr(i,j,t) = 0endifendifenddoenddoenddo! Identify individual contour areas on a 2D-space (only x-y, no time)! -------------------------------------------------------------------print*, '-->> identify individual contours'arr2(:,:,:) = 0 ! set arr2 to zerodo t = 1,nTimeslabel = 1 ! begin labeling with 1do j = 1,nLatsdo i = 1,nLongif ( (arr(i,j,t).eq.1).and.(arr2(i,j,t).lt.1) ) thenarr2(i,j,t) = label ! first identified grid point in countournrg = 1 ! number of grid pointsnrgc = 1 ! number of grid points counteriiv(nrgc) = i ! x-dir coordinates within contourjjv(nrgc) = j ! y-dir coordinates within contourdo while ( nrgc.le.nrg )do ii=-1,1do jj=-1,1if ( (arr(iiv(nrgc)+ii,jjv(nrgc)+jj,t).eq.1).and. &(arr2(iiv(nrgc)+ii,jjv(nrgc)+jj,t).ne.label) ) thenarr2(iiv(nrgc)+ii,jjv(nrgc)+jj,t) = labelnrg = nrg + 1iiv(nrg) = iiv(nrgc)+iijjv(nrg) = jjv(nrgc)+jjendifenddoenddonrgc = nrgc + 1enddolabel = label + 1endifenddoenddoenddo! Define temporal overlapping! --------------------------------------------------------------print*, '-->> overlapping'do t = 1,nTimesdo ii = 1,int(maxval(arr2(:,:,t))) ! loop over the number of individual contoursareacon = 0.areaover = 0.fraction = 0.do i = 1,nLongdo j = 1,nLatsif ( ( arr2(i,j,t).eq.ii) ) thenareacon = areacon+(111*(111*cos(2*3.14159/360*j)))endifif ( ( arr2(i,j,t).eq.ii).and.(arr2(i,j,t+1).ge.1) ) thenareaover = areaover+(111*(111*cos(2*3.14159/360*j)))endifenddoenddofraction = 1 / areacon * areaover!print*, ii, int(maxval(arr2(:,:,t))), fractiondo i = 1,nLongdo j = 1,nLatsif ( (fraction.lt.overlap).and.(arr2(i,j,t).eq.ii) ) thenarr2(i,j,t) = 0.endifenddoenddoenddoenddo! Identify individual contour areas (3D including time)! --------------------------------------------------------------print*, '-->> identify 3D contours'arr(:,:,:) = arr2(:,:,:)do t = 1,nTimesdo j = 1,nLatsdo i = 1,nLongif ( arr(i,j,t).ge.1 ) thenarr(i,j,t) = 1endifenddoenddoenddoarr2(:,:,:) = 0 ! set arr2 to zerolabel = 1 ! begin labeling with 1do t = 1,nTimesdo j = 1,nLatsdo i = 1,nLongif ( (arr(i,j,t).eq.1).and.(arr2(i,j,t).lt.1) ) thenarr2(i,j,t) = label ! first identified grid point in countournrg = 1 ! number of grid pointsnrgc = 1 ! number of grid points counteriiv(nrgc) = i ! x-dir coordinates within contourjjv(nrgc) = j ! y-dir coordinates within contourttv(nrgc) = t ! t-dir coordinates within contourttbegin = t ! genesis time! print*, i,j,t,labeldo while ( nrgc.le.nrg )do ii=-1,1do jj=-1,1do tt=-1,1if ( (arr(iiv(nrgc)+ii,jjv(nrgc)+jj,ttv(nrgc)+tt).eq.1).and. &(arr2(iiv(nrgc)+ii,jjv(nrgc)+jj,ttv(nrgc)+tt).ne.label) ) thenarr2(iiv(nrgc)+ii,jjv(nrgc)+jj,ttv(nrgc)+tt) = labelnrg = nrg + 1iiv(nrg) = iiv(nrgc)+iijjv(nrg) = jjv(nrgc)+jjttv(nrg) = ttv(nrgc)+tt! get time of lysisif (ttv(nrg).gt.ttend) thenttend = ttend + 1endifendifenddoenddoenddo! print*, iiv(nrg),jjv(nrg),ttv(nrg), label! print*, nrg, nrgc, ttv(nrgc)! if (nrgc.gt.20) stopnrgc = nrgc + 1enddoprint*, label, ttbegin, ttend,ttend-ttbeginlabel = label + 1endifenddoenddoenddoif ( verbose.eq.1 ) print*, ' -> number of individual contours: ', label-1! Temporal persistence! --------------------------------------------------------------print*, '-->> temporal persistence'do ii = 1,labelcounter = 0do t = 1,nTimesdo j = 1,nLatsdo i = 1,nLongif ( arr2(i,j,t).eq.ii) thencounter = counter + 1goto 234endifenddoenddo234 continueenddodo t = 1,nTimesdo j = 1,nLatsdo i = 1,nLongif ( (counter.lt.persistence ).and.(arr2(i,j,t).eq.ii) ) thenarr2(i,j,t) = 0.endifenddoenddoenddoenddo! get contour information! --------------------------------------------------------------print*, '-->> get contour information'! number of individual contours! length of individual contours! mean area size of individual contours! genesis/lysis date! genesis/lysis area at specific datenrcontours = 0do t = 1,nTimesdo j = 1,nLatsdo i = 1,nLongif ( arr2(i,j,t).gt.nrcontours ) thennrcontours = nrcontours + 1endifenddoenddoenddoprint*, ' -> number of tracked contours: ', nrcontours! change output array to a 0/1 array! --------------------------------------------------------------if ( outarr.eq.1 ) thendo t = 1,nTimesdo j = 1,nLatsdo i = 1,nLongif ( arr2(i,j,t).gt.0 ) thenarr2(i,j,t) = 1.endifenddoenddoenddoendifarr = arr211 continue! --------------------------------------------------------------! write CF-netCDF file! --------------------------------------------------------------call date_and_time (values=time_array)write(*,100), "-> Write ", trim(outfilename), time_array(5), ':', time_array(6), ':', time_array(7)! definitions! --------------------------------------------------------------cfoutinstit = cfinstitcfouttitle = 'contrack '//trim(contrackversion)cfoutsource = cfsourcewrite(cfoutcomment, 101) 'contrack is based upon the following input file attributes: ', &trim(cfcomment), &' -->> Contrack specifications:: ',&'contours identified that are -', trim(gorl), '- threshold value, ', &'threshold value for contour:', threshold, &', overlapping fraction:', overlap, &', persistence time steps:', persistence, &', anomalies calculated by subtracting meridional mean [1] or ', &'from the input file directly [0]: ', calcmeranom101 format(A,A,A,A,A,A,A,F8.3,A,F7.3,A,I5,A,A,I3)! start writing! ----------------------------------------------------------------! Initially set error to indicate no errors.status = 0! create the netCDF filestatus = nf90_create(trim(outfilename), NF90_CLOBBER, ncID)! IF (ierr.NE.0) GOTO 920! Dimensions -----------------------------------------------------status=nf90_def_dim(ncID,'longitude',nLong,LonDimID)status=nf90_def_dim(ncID,'latitude',nLats,LatDimID)status=nf90_def_dim(ncID,'time',nf90_unlimited, TimeDimID)! Variables ------------------------------------------------------status = nf90_def_var(ncID,'longitude',NF90_FLOAT,(/ LonDimID /),varLonID)status = nf90_put_att(ncID, varLonID, "standard_name", "longitude")status = nf90_put_att(ncID, varLonID, "units", "degree_east")status = nf90_def_var(ncID,'latitude',NF90_FLOAT,(/ LatDimID /),varLatID)status = nf90_put_att(ncID, varLatID, "standard_name", "latitude")status = nf90_put_att(ncID, varLatID, "units", "degree_north")status = nf90_def_var(ncID,'time',NF90_FLOAT, (/ TimeDimID /), varTimeID)status = nf90_put_att(ncID, varTimeID, "axis", "T")status = nf90_put_att(ncID, varTimeID, "calendar", "standard")status = nf90_put_att(ncID, varTimeID, "long_name", "time")status = nf90_put_att(ncID, varTimeID, "units", "hours since &1950-01-01 00:OO UTC")! Global Attributes -----------------------------------------------status = nf90_put_att(ncID, NF90_GLOBAL, 'Conventions', 'CF-1.0')status = nf90_put_att(ncID, NF90_GLOBAL, 'title', cfouttitle)status = nf90_put_att(ncID, NF90_GLOBAL, 'source', cfoutsource)status = nf90_put_att(ncID, NF90_GLOBAL, 'institution', cfoutinstit)status = nf90_put_att(ncID, NF90_GLOBAL, 'comment', cfoutcomment)! Specific variable -----------------------------------------------status = nf90_def_var(ncID,trim(outvarname),NF90_FLOAT,&(/ LonDimID, LatDimID, varTimeID /),varPVID)!status = nf90_def_var(ncID,trim(outvarname),NF90_FLOAT,&! (/ LatDimID, varTimeID /),varPVID)status = nf90_put_att(ncID, varPVID, "standard_name",trim(outstandvarname))status = nf90_put_att(ncID, varPVID, "units", trim(outvarunit))status = nf90_put_att(ncID, varPVID, '_FillValue', -999.99)! END variable definitions.status = nf90_enddef(ncID)IF (status.GT.0) THENprint*, 'An error occurred while attempting to ', &'finish definition mode.'!GOTO 920ENDIF! write DATA! -------------------------------------------------------------------status = nf90_put_var(ncID,varLonID,longitude)status = nf90_put_var(ncID,varLatID,latitude)!status = nf90_put_var(ncID,varTimeID, times)status = nf90_put_var(ncID,varTimeID, ×((persistence+1):(nTimes-persistence)))! Write block!status = nf90_put_var(ncID,varPVID, arr(:,:,:))status = nf90_put_var(ncID,varPVID, &arr(:,:,(persistence+1):(nTimes-persistence)))status = nf90_close(ncID)IF (status.NE.0) THENWRITE(0,*) trim(nf90_strerror(status))WRITE(0,*) 'An error occurred while attempting to ', &'close the netcdf file.'WRITE(0,*) 'in clscdf_CF'ENDIFclose(ncID)END