source: lmdz_wrf/trunk/tools/trajectories_overlap.f90 @ 2343

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