Changeset 1659 in lmdz_wrf

Timestamp:

Sep 27, 2017, 6:22:33 PM (8 years ago)

Author:

lfita

Message:

First compiling version of the code

File:

• trunk/tools/trajectories_overlap.f90 (modified) (11 diffs)

trunk/tools/trajectories_overlap.f90

@@ -24 +24 @@
 
 ! Local
-  INTEGER                                                :: ncid, rcode
+  INTEGER                                                :: funit, ncunit, rcode
+  INTEGER                                                :: ierr, ios
   INTEGER                                                :: Nvdims
   INTEGER, DIMENSION(1)                                  :: dims1D
@@ -30 +31 @@
   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               :: polys, polywctrs, finaltracks
+  REAL(r_k), DIMENSION(:,:,:), ALLOCATABLE               :: polywctrs, finaltracks
   REAL(r_k), DIMENSION(:,:,:,:), ALLOCATABLE             :: polytracks
   LOGICAL                                                :: gotvar
@@ -39 +41 @@
     Nmaxtrack
 
-  NAMELIST /files/ polyfilename, poly2Dn, poluyNn, polywctrn, polywarean, varnresx, varnresy, timen
+  NAMELIST /files/ polyfilename, poly2Dn, polyNn, polywctrn, polywarean, varnresx, varnresy, timen
   NAMELIST /config/ minarea, projectn, debug
 
@@ -62 +64 @@
 
   ! Reading namelist
-  funit = Fileunit()
+  funit = freeunit()
   OPEN(funit, file='trajectories_overlap.namelist', status='old', form='formatted', iostat=ios)
   msg = "Problems reading nameslit file: 'trajectories_overlap.namelist'" 
@@ -72 +74 @@
 
   ! Opening polygon file
-  rcode = nf90_open(TRIM(filename), NF90_NOWRITE, ncid)
+  rcode = nf90_open(TRIM(polyfilename), NF90_NOWRITE, ncunit)
   IF (rcode /= NF90_NOERR) CALL handle_err(rcode)
 
@@ -79 +81 @@
 
   ! polygons
-  gotvar = isin_ncunit(ncid, poly2Dn)
-  IF (.NOT.gotvar) THEN
-    msg = "File '" // TRIM() // "' does not have 2D polygons variable '" // TRIM(poly2Dn) // "'"
-    CALL ErrMsg(msg, fname, -1)
-  END IF
-  Nvdims = get_varNdims_ncunit(ncid, poly2Dn)
+  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
@@ -96 +98 @@
   msg = "Problems allocating 'polys'"
   CALL ErrMsg(msg, fname, ierr)
-  CALL get_varRK3D_ncunit(ncunit, dims3D(1), dims3D(2), dims3D(3), poly2Dn, polys)
+  CALL get_varI3D_ncunit(ncunit, dims3D(1), dims3D(2), dims3D(3), poly2Dn, polys)
+  IF (debug) PRINT *,"  Got variable 'polys':", UBOUND(polys)
 
   dimx = dims3D(1)
@@ -103 +106 @@
 
   ! N polygons
-  gotvar = isin_ncunit(ncid, polyNn)
-  IF (.NOT.gotvar) THEN
-    msg = "File '" // TRIM() // "' does not have number of polygons variable '" // TRIM(polyNn) // "'"
-    CALL ErrMsg(msg, fname, -1)
-  END IF
-  Nvdims = get_varNdims_ncunit(ncid, polyNn)
-  IF (Nvdims /= 2) THEN
+  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 2 dimensions !!"
-    CALL ErrMsg(msg, fname, -1)
-  END IF
-  dims2D = get_var2dims_ncunit(ncunit, polyNn)
+      " 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(dims2D(1), dims2D(2)), STAT=ierr)
+  ALLOCATE(polyN(dims1D(1)), STAT=ierr)
   msg = "Problems allocating 'polyN'"
   CALL ErrMsg(msg, fname, ierr)
-  CALL get_varI2D_ncunit(ncunit, dims1D(1), dims1D(2), polyNn, polyN)
+  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(ncid, polywctrn)
-  IF (.NOT.gotvar) THEN
-    msg = "File '" // TRIM() // "' does not have weighted-centerd polygon variable '" //              &
+  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(ncid, polywctrn)
+  Nvdims = get_varNdims_ncunit(ncunit, polywctrn)
   IF (Nvdims /= 3) THEN
     WRITE(num3S,'(I3)')Nvdims
@@ -141 +146 @@
   msg = "Problems allocating 'polywctrs'"
   CALL ErrMsg(msg, fname, ierr)
-  CALL get_varRK3D_ncunit(ncunit, dims3D(1), dims3D(2), dims3D(3), poly2Dn, polywctrs)
+  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(ncid, polywarean)
-  IF (.NOT.gotvar) THEN
-    msg = "File '" // TRIM() // "' does not have variable '" // TRIM(polywarean) // "'"
-    CALL ErrMsg(msg, fname, -1)
-  END IF
-  Nvdims = get_varNdims_ncunit(ncid, polywarean)
-  IF (Nvdims /= 3) THEN
+  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 //                &
@@ -163 +169 @@
   msg = "Problems allocating 'polyas'"
   CALL ErrMsg(msg, fname, ierr)
-  CALL get_varRK2D_ncunit(ncunit, dims2D(1), dims2D(2), poly2Dn, polyas)
+  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(ncid, varnresx)
-  IF (.NOT.gotvar) THEN
-    msg = "File '" // TRIM() // "' does not have variable '" // TRIM(varnresx) // "'"
-    CALL ErrMsg(msg, fname, -1)
-  END IF
-  CALL get_varRK0dims_ncunit(ncunit, varnresx, resx)
-  gotvar = isin_ncunit(ncid, varnresy)
-  IF (.NOT.gotvar) THEN
-    msg = "File '" // TRIM() // "' does not have variable '" // TRIM(varnresy) // "'"
-    CALL ErrMsg(msg, fname, -1)
-  END IF
-  CALL get_varRK0dims_ncunit(ncunit, varnresy, resy)
-
-  rcode = nf90_close(ncid)
+  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)
 
@@ -205 +215 @@
   ! Computing trajectories
   CALL poly_overlap_tracks_area(debug, dimx, dimy, dimt, minarea, polyN, polys, polywctrs, polyas,    &
-    Nmaxpoly, Nmaxtracks, polytracks, finaltracks)
+    Nmaxpoly, Nmaxtrack, polytracks, finaltracks)
 
   ! Writting trajectory files
 
   DEALLOCATE(polys)
-  DEALLOCATE(polyNs)
+  DEALLOCATE(polyN)
   DEALLOCATE(polywctrs)
   DEALLOCATE(polyas)