PROGRAM trajectories_overlap ! Porgram to compute trajectories from a series of consecutive maps of polygons ! L. Fita, CIMA. September 2017 ! ! To be used with a netcdf file which must contain: ! - Series of 2D maps of polygons (numbered by integers) ! - Areal weighted centers of the polygons ! - Weighted areas of the polygons ! USE module_definitions USE module_generic USE module_NCgeneric USE module_scientific USE netcdf IMPLICIT NONE ! namelist INTEGER :: missing_value REAL(r_k) :: minarea CHARACTER(len=50) :: polyfilename CHARACTER(len=50) :: poly2Dn, polyNn, polywctrn, polywarean, & varnresx, varnresy, timen, projectn, waycomp, waymulti LOGICAL :: debug ! Local INTEGER :: funit, ncunit, rcode INTEGER :: ierr, ios INTEGER :: Nvdims INTEGER, DIMENSION(1) :: dims1D INTEGER, DIMENSION(2) :: dims2D INTEGER, DIMENSION(3) :: dims3D INTEGER, DIMENSION(:), ALLOCATABLE :: polyN INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: polys REAL(r_k) :: resx, resy REAL(r_k), DIMENSION(:,:), ALLOCATABLE :: polyas REAL(r_k), DIMENSION(:,:,:), ALLOCATABLE :: polywctrs LOGICAL :: gotvar CHARACTER(len=3) :: num3S INTEGER :: dimx, dimy, dimt, Ncoord, Nmaxpoly, & Nmaxtrack NAMELIST /files/ polyfilename, poly2Dn, polyNn, polywctrn, polywarean, varnresx, varnresy, timen NAMELIST /config/ missing_value, minarea, projectn, waycomp, waymulti, debug !!!!!!! ! polyfilename: name of the file with the polygon information ! poly2Dn: name of the variable with the maps of 2D polygons ! polyNn: name of the variable with the number of 2D polygons per time-step ! polywctrn: name of the varibale with the coordinates of the weighted centers of all the polygons ! polywarean: name of the variable with the weighted areas of the polygons ! varnres[x/y]: name of the variable with the resolution along x an y axis (same units as polywarean) ! lonn: name of the variable with the longitudes ! latn: name of the variable with the latitudes ! timen: name of the variable with the times ! missing_value: value to determine the missing value ! minarea: minimal area of the polygons to determine which polygons to use ! projectn: name of the projection of the data in the file. Available ones: ! 'ctsarea': Constant Area ! 'lon/lat': for regular longitude-latitude ! 'nadir-sat,[lonNADIR],[latNADIR]': for satellite data with the resolution given at nadir (lonNADIR, latNADIR) ! waycomp: how to copmute when intermediate ASCII files are used ! 'scratch': everything from the beginning ! 'continue': skipt that parts which already have the ascii file written ! waymulti: methodology to follow when multiple polygons are given for the same track ! 'mean': get coordinates from the areal-weighted mean of the centers of the given polygons and their areas ! 'largest': get the coorindates of the largest polygon ! 'closest': get the coordinates of the closest polygon fname = 'trajectories_overlap' ! Reading namelist funit = freeunit() OPEN(funit, file='trajectories_overlap.namelist', status='old', form='formatted', iostat=ios) msg = "Problems reading nameslit file: 'trajectories_overlap.namelist'" CALL ErrMsg(msg, 'main', ios) READ(funit, files) READ(funit, config) CLOSE(funit) ! Opening polygon file rcode = nf90_open(TRIM(polyfilename), NF90_NOWRITE, ncunit) IF (rcode /= NF90_NOERR) CALL handle_err(rcode) ! Getting variables !! ! polygons gotvar = isin_ncunit(ncunit, poly2Dn) IF (.NOT.gotvar) THEN msg= "File '"//TRIM(polyfilename)//"' does not have 2D polygons variable '" // TRIM(poly2Dn) // "'" CALL ErrMsg(msg, fname, -1) END IF Nvdims = get_varNdims_ncunit(ncunit, poly2Dn) IF (Nvdims /= 3) THEN WRITE(num3S,'(I3)')Nvdims msg = "Variable with 2D polygons '" // TRIM(poly2Dn) // "' has: " // num3S // & " it must have 3 dimensions !!" CALL ErrMsg(msg, fname, -1) END IF dims3D = get_var3dims_ncunit(ncunit, poly2Dn) IF (ALLOCATED(polys)) DEALLOCATE(polys) ALLOCATE(polys(dims3D(1), dims3D(2), dims3D(3)), STAT=ierr) msg = "Problems allocating 'polys'" CALL ErrMsg(msg, fname, ierr) CALL get_varI3D_ncunit(ncunit, dims3D(1), dims3D(2), dims3D(3), poly2Dn, polys) IF (debug) PRINT *," Got variable 'polys':", UBOUND(polys) ! Getting rid of the mising values WHERE (polys == missing_value) polys = 0 END WHERE dimx = dims3D(1) dimy = dims3D(2) dimt = dims3D(3) ! N polygons gotvar = isin_ncunit(ncunit, polyNn) IF (.NOT.gotvar) THEN msg="File '"//TRIM(polyfilename)//"' does not have number of polygons variable '" // TRIM(polyNn) & // "'" CALL ErrMsg(msg, fname, -1) END IF Nvdims = get_varNdims_ncunit(ncunit, polyNn) IF (Nvdims /= 1) THEN WRITE(num3S,'(I3)')Nvdims msg = "Variable with N polygons '" // TRIM(polyNn) // "' has: " // num3S // & " it must have 1 dimensions !!" CALL ErrMsg(msg, fname, -1) END IF dims1D = get_var1dims_ncunit(ncunit, polyNn) IF (ALLOCATED(polyN)) DEALLOCATE(polyN) ALLOCATE(polyN(dims1D(1)), STAT=ierr) msg = "Problems allocating 'polyN'" CALL ErrMsg(msg, fname, ierr) CALL get_varI1D_ncunit(ncunit, dims1D(1), polyNn, polyN) IF (debug) PRINT *," Got variable 'polyN':", UBOUND(polyN) ! weighted-center of the polygons gotvar = isin_ncunit(ncunit, polywctrn) IF (.NOT.gotvar) THEN msg = "File '" // TRIM(polyfilename) // "' does not have weighted-centerd polygon variable '" // & TRIM(polywctrn) // "'" CALL ErrMsg(msg, fname, -1) END IF Nvdims = get_varNdims_ncunit(ncunit, polywctrn) IF (Nvdims /= 3) THEN WRITE(num3S,'(I3)')Nvdims msg = "Variable with weighted-centerd polygons '" // TRIM(polywctrn) // "' has: " // num3S // & " it must have 3 dimensions !!" CALL ErrMsg(msg, fname, -1) END IF dims3D = get_var3dims_ncunit(ncunit, polywctrn) IF (ALLOCATED(polywctrs)) DEALLOCATE(polywctrs) ALLOCATE(polywctrs(dims3D(1), dims3D(2), dims3D(3)), STAT=ierr) msg = "Problems allocating 'polywctrs'" CALL ErrMsg(msg, fname, ierr) CALL get_varRK3D_ncunit(ncunit, dims3D(1), dims3D(2), dims3D(3), polywctrn, polywctrs) IF (debug) PRINT *," Got variable 'polyectrs':", UBOUND(polywctrs) Nmaxpoly = dims3D(2) Ncoord = dims3D(3) ! area of polygons gotvar = isin_ncunit(ncunit, polywarean) IF (.NOT.gotvar) THEN msg = "File '" // TRIM(polyfilename) // "' does not have variable '" // TRIM(polywarean) // "'" CALL ErrMsg(msg, fname, -1) END IF Nvdims = get_varNdims_ncunit(ncunit, polywarean) IF (Nvdims /= 2) THEN WRITE(num3S,'(I3)')Nvdims msg = "Variable with area polygons '" // TRIM(polywarean) // "' has: " // num3S // & " it must have 2 dimensions !!" CALL ErrMsg(msg, fname, -1) END IF dims2D = get_var2dims_ncunit(ncunit, polywarean) IF (ALLOCATED(polyas)) DEALLOCATE(polyas) ALLOCATE(polyas(dims2D(1), dims2D(2)), STAT=ierr) msg = "Problems allocating 'polyas'" CALL ErrMsg(msg, fname, ierr) CALL get_varRK2D_ncunit(ncunit, dims2D(1), dims2D(2), polywarean, polyas) IF (debug) PRINT *," Got variable 'polyas':", UBOUND(polyas) ! Resolution along x/y axis gotvar = isin_ncunit(ncunit, varnresx) IF (.NOT.gotvar) THEN msg = "File '" // TRIM(polyfilename) // "' does not have variable '" // TRIM(varnresx) // "'" CALL ErrMsg(msg, fname, -1) END IF CALL get_varRK0D_ncunit(ncunit, varnresx, resx) IF (debug) PRINT *," Got variable 'resx':", resx gotvar = isin_ncunit(ncunit, varnresy) IF (.NOT.gotvar) THEN msg = "File '" // TRIM(polyfilename) // "' does not have variable '" // TRIM(varnresy) // "'" CALL ErrMsg(msg, fname, -1) END IF CALL get_varRK0D_ncunit(ncunit, varnresy, resy) IF (debug) PRINT *," Got variable 'resy':", resy rcode = nf90_close(ncunit) IF (rcode /= NF90_NOERR) CALL handle_err(rcode) ! Preparing outcome (trajectories) Nmaxtrack = Nmaxpoly*dimt/100 IF (debug) THEN PRINT *,'Computing trajectories _______' PRINT *,' 2D space:', dimx, ' x', dimy PRINT *,' number time-steps:', dimt PRINT *,' maximum number of polygons:', Nmaxpoly PRINT *,' maximum number of trajectories:', Nmaxtrack END IF ! Computing trajectories CALL poly_overlap_tracks_area_ascii(debug, waycomp, dimx, dimy, dimt, minarea, polyN, polys, & polywctrs, polyas, Nmaxpoly, Nmaxtrack, waymulti) ! Writting trajectory files DEALLOCATE(polys) DEALLOCATE(polyN) DEALLOCATE(polywctrs) DEALLOCATE(polyas) END PROGRAM trajectories_overlap