[Dart-dev] [5695] DART/branches/development/location: 2D version of the annulus code.

nancy at ucar.edu nancy at ucar.edu
Thu Apr 12 15:29:59 MDT 2012


Revision: 5695
Author:   nancy
Date:     2012-04-12 15:29:58 -0600 (Thu, 12 Apr 2012)
Log Message:
-----------
2D version of the annulus code.  supports limits for the
azimuth angle of -90/90 (latitude-like), 0-360 (longitude-
like with wrap around 360/0), or unlimited.  radius limits
are set with a min/max value in a namelist and are assumed
to be in meters.

Modified Paths:
--------------
    DART/branches/development/location/twod_annulus/README
    DART/branches/development/location/twod_annulus/location_mod.f90
    DART/branches/development/location/twod_annulus/test/input.nml
    DART/branches/development/location/twod_annulus/test/path_names_location_test
    DART/branches/development/location/twod_annulus/test/test.in

Added Paths:
-----------
    DART/branches/development/location/twod_annulus/

-------------- next part --------------
Modified: DART/branches/development/location/twod_annulus/README
===================================================================
--- DART/branches/development/location/annulus/README	2012-04-11 17:28:19 UTC (rev 5692)
+++ DART/branches/development/location/twod_annulus/README	2012-04-12 21:29:58 UTC (rev 5695)
@@ -4,26 +4,15 @@
 #
 # DART $Id$
 
-Last modified 29 Jan, 2010 by nancy collins
+A location consists of the pair (azimuth, radius).
 
-made consistent with other location modules, including the
-revised get_close() routines.   this code didn't have any
-min/max limits on the radius, and was still using lon/lat
-in a few subroutines.  i put in a stub of a namelist for
-the radius limits, and changed the remaining lon/lats
-to azm and rad.
+This is copied from the cylindrical version of the
+annulus location module which has an azimuthal angle,
+a radius, and a height.
 
-Last modified June 28, 2004 by Jim Hansen.
-
-The location_mod.f90 for the annulus is a modified version
-of ~/DART/location/threed_sphere/location_mod.f90.  The most
-important thing to keep in mind in the annulus version is that
-longitude actually mean azimuthal angle, and latitude actually
-means radius.
-
 The location_mod.f90 requires the parameters inner_rad (inner
 radius of the annulus) and outer_rad (outer radius of the annulus)
 from the input.nml.
 
-Longitude (azimuthal angle) is in degrees, latitude (radius) is
-in meters.
+Azimuthal angle is in degrees, radius is in meters.
+

Modified: DART/branches/development/location/twod_annulus/location_mod.f90
===================================================================
--- DART/branches/development/location/annulus/location_mod.f90	2012-04-11 17:28:19 UTC (rev 5692)
+++ DART/branches/development/location/twod_annulus/location_mod.f90	2012-04-12 21:29:58 UTC (rev 5695)
@@ -10,12 +10,12 @@
 ! $Revision$
 ! $Date$
 
-! Implements location interfaces for a three dimensional annulus
-! with a vertical coordinate based on the models native set of
-! discrete levels. The internal representation of the location is 
-! currently implemented as radians from 0 to 2 PI for the azimuthal
-! direction (longitude-like).  The radial distance is latitude-like,
-! and the vertical coordinate is zero at the bottom of the annulus.
+! Implements location interfaces for a two dimensional annulus
+! discrete levels. The internal representation of the azimuth is 
+! implemented as radians and is namelist selectable to be valid 
+! from -PI/2 to PI/2 (like latitudes) or from 0 to 2 PI (like longitudes)
+! The radial distance is in meters and the inner and outer boundary
+! distances are namelist settable.
 !
 
 use      types_mod, only : r8, PI, RAD2DEG, DEG2RAD, MISSING_R8, MISSING_I
@@ -36,8 +36,7 @@
           operator(==), operator(/=), get_dist, get_close_obs_destroy, &
           nc_write_location_atts, nc_get_location_varids, nc_write_location, &
           vert_is_height, vert_is_pressure, vert_is_undef, vert_is_level, &
-          vert_is_surface, has_vertical_localization, VERTISSURFACE, &
-          VERTISLEVEL, VERTISHEIGHT
+          vert_is_surface, has_vertical_localization
 
 ! version controlled file description for error handling, do not edit
 character(len=128), parameter :: &
@@ -47,17 +46,17 @@
 
 type location_type
    private
-   real(r8) :: azm, rad, vloc
-   integer  :: which_vert
-   ! which_vert determines if the location is by level or by height
-   ! -1 ==> obs is on surface
-   ! 1 ===> obs is by level
-   ! 3 ===> obs is by height
+   real(r8) :: azm, rad
+   integer  :: which_azm
+   ! which_azm determines the valid boundaries of the azimuth
+   ! -1 ==> no limits
+   !  0 ==> -PI/2 to PI/2 (like latitudes)
+   !  1 ==> 0 to 2PI (like longitudes)
 end type location_type
 
-integer, parameter :: VERTISSURFACE  = -1 ! surface
-integer, parameter :: VERTISLEVEL    =  1 ! by level
-integer, parameter :: VERTISHEIGHT   =  3 ! by height
+integer, parameter :: AZMISUNBOUND  = -1 ! unbounded
+integer, parameter :: AZMISLAT      =  0 ! like lats
+integer, parameter :: AZMISLON      =  1 ! like lons
 
 type get_close_type
    private
@@ -69,14 +68,12 @@
 logical :: ran_seq_init = .false.
 logical, save :: module_initialized = .false.
 
-integer,              parameter :: LocationDims = 3
-character(len = 129), parameter :: LocationName = "loc_annulus"
+integer,              parameter :: LocationDims = 2
+character(len = 129), parameter :: LocationName = "loc_2d_annulus"
 character(len = 129), parameter :: LocationLName = &
-                                   "Annulus location: azimuthal angle, radius, and height"
+                                   "2D Annulus location: azimuthal angle, radius"
 
-! really just a placeholder.  there was a comment that this code
-! needs a namelist with a min & max limit on the radius, but 
-! the code no longer has one.  
+! limits on the radius (sets inner and outer radius values)
 real(r8) :: min_radius = 0.0_r8
 real(r8) :: max_radius = 100000.0_r8
 
@@ -121,7 +118,7 @@
 if(do_nml_term()) write(     *     , nml=location_nml)
 
 ! copy code from threed sphere module for handing the
-! options in the vertical?  e.g. distances?
+! distances?
 
 end subroutine initialize_module
 
@@ -129,9 +126,7 @@
 
 function get_dist(loc1, loc2, kind1, kind2)
 
-! Compute distance between 2 locations.  Right now the distance only
-! depends on the horizontal.  A namelist option might need to be added
-! that supports computing a true 3d distance.
+! Compute distance between 2 locations.
 
 type(location_type), intent(in) :: loc1, loc2
 integer, optional,   intent(in) :: kind1, kind2
@@ -141,8 +136,6 @@
 
 if ( .not. module_initialized ) call initialize_module
 
-! FIXME: this does not take into account any vertical separation
-! convert from cylindrical to cartesian coordinates
 x1 = loc1%rad * cos(loc1%azm)
 y1 = loc1%rad * sin(loc1%azm)
 x2 = loc2%rad * cos(loc2%azm)
@@ -168,10 +161,9 @@
 
 loc_eq = .false.
 
-if ( loc1%which_vert /= loc2%which_vert ) return
+if ( loc1%which_azm /= loc2%which_azm ) return
 if ( abs(loc1%azm  - loc2%azm ) > epsilon(loc1%azm ) ) return
 if ( abs(loc1%rad  - loc2%rad ) > epsilon(loc1%rad ) ) return
-if ( abs(loc1%vloc - loc2%vloc) > epsilon(loc1%vloc) ) return
 
 loc_eq = .true.
 
@@ -198,33 +190,32 @@
 function get_location(loc)
  
 ! Given a location type (where the azimuthal angle is in radians), this 
-! routine return the azimuthal angle in degrees, the radius, and the vert 
+! routine return the azimuthal angle in degrees, and the radius.
 
 type(location_type), intent(in) :: loc
-real(r8), dimension(3) :: get_location
+real(r8), dimension(2) :: get_location
 
 if ( .not. module_initialized ) call initialize_module
 
 get_location(1) = loc%azm * RAD2DEG                 
 get_location(2) = loc%rad
-get_location(3) = loc%vloc     
 
 end function get_location
 
 !----------------------------------------------------------------------------
 
-function set_location_single(azm, rad, vert_loc, which_vert)
+function set_location_single(azm, rad, which_azm)
  
-! Puts the given azimuthal angle (in degrees), radius, and vertical 
+! Puts the given azimuthal angle (in degrees), and radius
 ! location into a location datatype.
 
 real(r8), intent(in) :: azm, rad
-real(r8), intent(in) :: vert_loc
-integer,  intent(in) :: which_vert
+integer,  intent(in) :: which_azm
 type (location_type) :: set_location_single
 
 if ( .not. module_initialized ) call initialize_module
 
+! FIXME: test range based on which_azm
 if(azm < 0.0_r8 .or. azm > 360.0_r8) then
    write(errstring,*)'azimuthal angle (',azm,') is not within range [0,360]'
    call error_handler(E_ERR, 'set_location', errstring, source, revision, revdate)
@@ -233,16 +224,8 @@
 set_location_single%azm = azm * DEG2RAD
 set_location_single%rad = rad 
 
-if(which_vert /= VERTISSURFACE .and. &
-   which_vert /= VERTISLEVEL   .and. &
-   which_vert /= VERTISHEIGHT  ) then
-   write(errstring,*)'which_vert (',which_vert,') must be -1, 1 or 3'
-   call error_handler(E_ERR,'set_location', errstring, source, revision, revdate)
-endif
+set_location_single%which_azm = which_azm
 
-set_location_single%which_vert = which_vert
-set_location_single%vloc = vert_loc
-
 end function set_location_single
 
 !----------------------------------------------------------------------------
@@ -250,19 +233,19 @@
 function set_location_array(list)
  
 ! Location semi-independent interface routine
-! given 4 float numbers, call the underlying set_location routine
+! given 3 float numbers, call the underlying set_location routine
 
 real(r8), intent(in) :: list(:)
 type (location_type) :: set_location_array
 
 if ( .not. module_initialized ) call initialize_module
 
-if (size(list) < 4) then
-   write(errstring,*)'requires 4 input values'
+if (size(list) < 3) then
+   write(errstring,*)'requires 3 input values'
    call error_handler(E_ERR, 'set_location', errstring, source, revision, revdate)
 endif
 
-set_location_array = set_location_single(list(1), list(2), list(3), nint(list(4)))
+set_location_array = set_location_single(list(1), list(2), nint(list(3)))
 
 end function set_location_array
 
@@ -278,8 +261,7 @@
 
 set_location_missing%azm        = MISSING_R8
 set_location_missing%rad        = MISSING_R8
-set_location_missing%vloc       = MISSING_R8
-set_location_missing%which_vert = MISSING_I
+set_location_missing%which_azm  = MISSING_I
 
 end function set_location_missing
 
@@ -300,22 +282,20 @@
 ! module for warnings about compiler bugs before you change
 ! this code.
 
-query_location = loc%which_vert
+query_location = loc%which_azm
 
 if (.not. present(attr)) return
 
 select case(attr)
-   case ('which_vert','WHICH_VERT')
-      query_location = loc%which_vert
+   case ('which_azm','WHICH_AZM')
+      query_location = loc%which_azm
    case ('rad','RAD','radius','RADIUS')
       query_location = loc%rad
    case ('azm','AZM','azimuth','AZIMUTH')
       query_location = loc%azm
-   case ('vloc','VLOC')
-      query_location = loc%vloc
    case default
        call error_handler(E_ERR, 'query_location:', &
-          'Only "azm","rad","vloc","which_vert" are legal attributes to request from location', &
+          'Only "azm","rad","which_azm" are legal attributes to request from location', &
           source, revision, revdate)
 end select
 
@@ -340,8 +320,7 @@
 logical             :: writebuf
 character(len = 128) :: string1
 
-! 10 format(1x,3(f22.14,1x),i4)   ! old
-10 format(1X,F21.16,2(1X,G25.16),1X,I2)
+10 format(1X,F21.16,1X,G25.16,1X,I2)
 
 if ( .not. module_initialized ) call initialize_module
 
@@ -351,12 +330,11 @@
 ! output file; test for ascii or binary, write what's asked, and return
 if (.not. writebuf) then
    if (ascii_file_format(fform)) then
-      ! Write out pressure or level along with integer tag
-      ! we know azm is between 0, 360, and which_vert is a single digit.
-      write(locfile, '(''loc3a'')' ) 
-      write(locfile, 10) loc%azm, loc%rad, loc%vloc, loc%which_vert
+      ! azm is between -90, 360, and which_vert is a single digit.
+      write(locfile, '(''loc2a'')' ) 
+      write(locfile, 10) loc%azm, loc%rad, loc%which_azm
    else
-      write(locfile) loc%azm, loc%rad, loc%vloc, loc%which_vert
+      write(locfile) loc%azm, loc%rad, loc%which_azm
    endif
    return
 endif
@@ -374,31 +352,16 @@
 ! hectopascals instead of pascals for pressure, etc.
 
 ! this must be the sum of the longest of the formats below.
-charlength = 85
+charlength = 48
 
 if (len(charstring) < charlength) then
    write(errstring, *) 'charstring buffer must be at least ', charlength, ' chars long'
    call error_handler(E_ERR, 'write_location', errstring, source, revision, revdate)
 endif
 
-write(string1, '(A,F12.8,1X,G15.8,A)') 'Azm(deg)/Radius(m): ',  &
-   loc%azm*RAD2DEG, loc%rad, ' Depth:'
+write(charstring, '(A,F12.8,1X,G15.8)') 'Azm(deg)/Radius(m): ',  &
+   loc%azm*RAD2DEG, loc%rad
 
-! i am attempting to make these line up so if you have a list of mixed
-! vertical units, they all take the same number of columns.  thus the extra
-! white space around some of the labels below.
-select case  (loc%which_vert)
-   case (VERTISSURFACE)
-      write(charstring, '(A,1X,G15.6,A)') trim(string1), loc%vloc, ' surface (hPa)'
-   case (VERTISLEVEL)
-      write(charstring, '(A,1X,F6.0,A)')  trim(string1), loc%vloc, '          level'
-   case (VERTISHEIGHT)
-      write(charstring, '(A,1X,G15.6,A)') trim(string1), loc%vloc / 1000.0_r8, ' km'
-   case default
-      write(errstring, *) 'unrecognized key for vertical type: ', loc%which_vert
-      call error_handler(E_ERR, 'write_location', errstring, source, revision, revdate)
-end select
-
 end subroutine write_location
 
 !----------------------------------------------------------------------------
@@ -419,16 +382,16 @@
 if (ascii_file_format(fform)) then
    read(locfile, '(a5)' ) header
 
-   if(header /= 'loc3a') then
-      write(errstring,*)'Expected location header "loc3a" in input file, got ', header
+   if(header /= 'loc2a') then
+      write(errstring,*)'Expected location header "loc2a" in input file, got ', header
       call error_handler(E_ERR, 'read_location', errstring, source, revision, revdate)
    endif
    ! Now read the location data value
    read(locfile, *) read_location%azm, read_location%rad, &
-                    read_location%vloc, read_location%which_vert
+                    read_location%which_azm
 else
    read(locfile) read_location%azm, read_location%rad, &
-                 read_location%vloc, read_location%which_vert
+                 read_location%which_azm
 endif
 
 end function read_location
@@ -443,7 +406,8 @@
 type(location_type), intent(out) :: location
 logical, intent(in), optional    :: set_to_default
 
-real(r8) :: azm, rad, minazm, maxazm, minrad, maxrad
+real(r8) :: azm, rad, minazm, maxazm, minrad, maxrad, minv, maxv
+integer :: r
 
 if ( .not. module_initialized ) call initialize_module
 
@@ -451,47 +415,63 @@
 if(present(set_to_default)) then
    if(set_to_default) then
       location%azm        = 0.0_r8
-      location%rad        = 0.0_r8
-      location%vloc       = 0.0_r8
-      location%which_vert = 0    ! zero is an invalid vert type
+      location%rad        = min_radius
+      location%which_azm = -1   ! unlimited angle
       return
    endif
 endif
 
-write(*, *)'Vertical coordinate options'
-write(*, *)'-1 -> surface, 1 -> model level, 3 -> depth'
+write(*, *)'Azimuth coordinate limit options'
+write(*, *)'-1 -> none, 0 -> -90, 90, 1 -> 0, 360'
 
-100   read(*, *) location%which_vert
-if(location%which_vert == VERTISLEVEL ) then
-   write(*, *) 'Vertical coordinate model level'
-   read(*, *) location%vloc
-else if(location%which_vert == VERTISHEIGHT ) then
-   write(*, *) 'Vertical coordinate depth (in negative m)'
-   read(*, *) location%vloc
-   do while (location%vloc > 0)
-      write(*, *) 'Depth must be negative (zero at top of fluid), please try again'
-      read(*, *) location%vloc
-   end do
-else if(location%which_vert == VERTISSURFACE ) then
-   write(*, *) 'Vertical coordinate surface pressure (in hPa)'
-   read(*, *) location%vloc
-   location%vloc = 100.0 * location%vloc
+100   read(*, *) location%which_azm
+if(location%which_azm == AZMISUNBOUND ) then
+   minazm = -HUGE(r8)
+   maxazm =  HUGE(r8)
+else if(location%which_azm == AZMISLAT ) then
+   minazm = -PI/2
+   maxazm =  PI/2
+else if(location%which_azm == AZMISLON ) then
+   minazm = 0.0_r8
+   maxazm = 2.0_r8*PI
 else
-   write(*, *) 'Wrong choice of which_vert try again between -1, 1, and 3'
+   write(*, *) 'Wrong choice of which_azm try again, valid: -1, 0, and 1'
    go to 100
 end if
 
-write(*, *) 'Input azimuthal angle: value 0 to 360.0 or a negative number '
-write(*, *) 'for uniformly distributed random location in the horizontal.'
-read(*, *) azm
+r = 1
+do while (r > 0)
 
-do while(azm > 360.0_r8)
-   write(*, *) 'Input value greater than 360.0 is illegal, please try again'
-   read(*, *) azm
-end do
+   write(*, *) 'Input 0 to specify a value for the location, or'
+   write(*, *) '-1 for a uniformly distributed random location'
+   read(*, *) r
 
-if(azm < 0.0_r8) then
+   if (r > 0) write(*, *) 'Please input 0 or -1 for selection'
+enddo
 
+if (r == 0) then
+101 continue
+   if (location%which_azm == AZMISUNBOUND) then
+      write(*, *) 'Input value for aziumth in degrees '
+      read(*,*) azm
+   else if (location%which_azm == AZMISLAT) then
+      write(*, *) 'Input value for aziumth in degrees (-90 to 90) '
+      read(*,*) azm
+      if (azm < -90.0_r8 .or. azm > 90.0_r8) then
+         write(*,*) 'Illegal value; must be between -90 and 90'
+         goto 101
+      endif
+   else if (location%which_azm == AZMISLON) then
+      write(*, *) 'Input value for aziumth in degrees (0 to 360)'
+      read(*,*) azm
+      if (azm < 0.0_r8 .or. azm > 360.0_r8) then
+         write(*,*) 'Illegal value; must be between 0 and 360'
+         goto 101
+      endif
+   endif
+   location%azm = azm * DEG2RAD
+
+else
    ! Need to make sure random sequence is initialized
 
    if(.not. ran_seq_init) then
@@ -499,31 +479,66 @@
       ran_seq_init = .TRUE.
    endif
 
-   minazm = -1.0
-   do while (minazm < 0.0 .or. minazm > 360.0) 
-      write(*, *) 'Input minimum azimuthal angle (0 to 360.0)'
-      read(*, *) minazm
-      if (minazm < 0.0 .or. minazm > 360.0) then
-         write(*, *) 'Angle must be between 0 to 360.0'
-      endif
-   enddo
-   minazm = minazm * DEG2RAD
+102 continue
+   write(*, *) 'Input minimum azimuth value in degrees '
+   read(*, *) minv
+   if (location%which_azm == AZMISLAT) then
+      if (minv < -90.0_r8) then
+         write(*,*) 'Illegal value; minimum must be >= -90'
+         goto 102
+      endif   
+   else if (location%which_azm == AZMISLON) then
+      if (minv < 0.0_r8) then
+         write(*,*) 'Illegal value; minimum must be >= 0'
+         goto 102
+      endif   
+   endif
 
-   maxazm = -1.0
-   do while (maxazm < 0.0 .or. maxazm > 360.0) 
-      write(*, *) 'Input maximum azimuthal angle (0 to 360.0)'
-      read(*, *) maxazm
-      if (maxazm < 0.0 .or. maxazm > 360.0) then
-         write(*, *) 'Angle must be between 0 to 360.0'
+103 continue
+   write(*, *) 'Input maximum azimuth value in degrees '
+   read(*, *) maxv
+   if (location%which_azm == AZMISLAT) then
+      if (maxv > 90.0_r8) then
+         write(*,*) 'Illegal value; maximum must be <= 90'
+         goto 103
+      endif   
+   else if (location%which_azm == AZMISLON) then
+      if (maxv > 360.0_r8) then
+         write(*,*) 'Illegal value; maximum must be <= 360'
+         goto 103
+      endif   
+   endif
+
+   minv = minv * DEG2RAD
+   maxv = maxv * DEG2RAD
+
+   ! Azimuth is random from minazm to maxazm, handle wrap around 360.0
+   if (location%which_azm == AZMISLON) then
+      if (minv > maxv) maxv = maxv + 2.0_r8 * PI
+   endif
+
+   location%azm = random_uniform(ran_seq) * (maxv-minv) + minv
+
+   if (location%which_azm == AZMISLON) then
+      if (location%azm > 2.0_r8 * PI) location%azm = location%azm - 2.0_r8 * PI
+   endif
+
+endif
+
+if (r == 0) then
+   rad = -1.0
+   do while (rad < min_radius .or. rad > max_radius) 
+      write(*, *) 'Input radius '
+      read(*, *) rad
+      if (rad < min_radius .or. rad > max_radius) then
+         write(*, *) 'Radius must be between ', min_radius, ' and ', max_radius
       endif
    enddo
-   maxazm = maxazm * DEG2RAD
 
-   ! Azimuth is random from minazm to maxazm, handle wrap around 360.0
-   if (minazm > maxazm) maxazm = maxazm + 2.0_r8 * PI
-   location%azm = random_uniform(ran_seq) * (maxazm-minazm) + minazm
-   if (location%azm > 2.0_r8 * PI) location%azm = location%azm - 2.0_r8 * PI
+   location%rad = rad
 
+else
+
    minrad = -1.0
    do while (minrad < min_radius) 
       write(*, *) 'Input minimum radius '
@@ -532,7 +547,7 @@
          write(*, *) 'Radius must be larger or equal to ', min_radius
       endif
    enddo
-
+   
    maxrad = -1.0
    do while (maxrad > max_radius .or. maxrad <= minrad) 
       write(*, *) 'Input maximum radius '
@@ -545,25 +560,10 @@
    ! Radius must be area weighted to obtain proper random realizations
    location%rad = sqrt(random_uniform(ran_seq)) * (maxrad-minrad) + minrad
 
-   write(*, *) 'random location is ', location%azm / DEG2RAD, &
+   write(*, *) 'random location is ', location%azm * RAD2DEG, &
                                       location%rad 
+endif
 
-else
-
-   rad = -1.0
-   do while (rad < min_radius .or. rad > max_radius) 
-      write(*, *) 'Input radius '
-      read(*, *) rad
-      if (rad < min_radius .or. rad > max_radius) then
-         write(*, *) 'Radius must be between ', min_radius, ' and ', max_radius
-      endif
-   enddo
-
-   location%rad = rad
-   location%azm = azm*DEG2RAD
-
-end if
-
 end subroutine interactive_location
 
 !----------------------------------------------------------------------------
@@ -604,25 +604,12 @@
 call nc_check(nf90_put_att(ncFileID, VarID, 'long_name', &
         trim(LocationLName)), 'nc_write_location_atts', 'location:long_name')
 call nc_check(nf90_put_att(ncFileID, VarID, 'storage_order',     &
-        'Azimuth Radius Vertical'), 'nc_write_location_atts', 'location:storage_order')
+        'Azimuth Radius'), 'nc_write_location_atts', 'location:storage_order')
 call nc_check(nf90_put_att(ncFileID, VarID, 'units',     &
-        'degrees meters which_vert'), 'nc_write_location_atts', 'location:units')
+        'degrees meters'), 'nc_write_location_atts', 'location:units')
 
 ! Define the ancillary vertical array and attributes
 
-call nc_check(nf90_def_var(ncid=ncFileID, name='which_vert', xtype=nf90_int, &
-          dimids=(/ ObsNumDimID /), varid=VarID), &
-            'nc_write_location_atts', 'which_vert:def_var')
-
-call nc_check(nf90_put_att(ncFileID, VarID, 'long_name', 'vertical coordinate system code'), &
-           'nc_write_location_atts', 'which_vert:long_name')
-call nc_check(nf90_put_att(ncFileID, VarID, 'VERTISSURFACE', VERTISSURFACE), &
-           'nc_write_location_atts', 'which_vert:VERTISSURFACE')
-call nc_check(nf90_put_att(ncFileID, VarID, 'VERTISLEVEL', VERTISLEVEL), &
-           'nc_write_location_atts', 'which_vert:VERTISLEVEL')
-call nc_check(nf90_put_att(ncFileID, VarID, 'VERTISHEIGHT', VERTISHEIGHT), &
-           'nc_write_location_atts', 'which_vert:VERTISHEIGHT')
-
 ! If we made it to here without error-ing out ... we're good.
 
 ierr = 0
@@ -652,8 +639,7 @@
 call nc_check(nf90_inq_varid(ncFileID, 'location', varid=LocationVarID), &
           'nc_get_location_varids', 'inq_varid:location '//trim(fname))
 
-call nc_check(nf90_inq_varid(ncFileID, 'which_vert', varid=WhichVertVarID), &
-          'nc_get_location_varids', 'inq_varid:which_vert '//trim(fname))
+WhichVertVarID = -99
 
 end subroutine nc_get_location_varids
 
@@ -684,11 +670,6 @@
           start=(/ 1, obsindex /), count=(/ LocationDims, 1 /) ), &
             'nc_write_location', 'put_var:location')
 
-intval = loc%which_vert
-call nc_check(nf90_put_var(ncFileID, WhichVertVarID, intval, &
-          start=(/ obsindex /), count=(/ 1 /) ), &
-            'nc_write_location','put_var:vert' )
-
 end subroutine nc_write_location
 
 !----------------------------------------------------------------------------
@@ -779,9 +760,9 @@
 
 if ( .not. module_initialized ) call initialize_module
 
-if ((minl%which_vert /= maxl%which_vert) .or. &
-    (minl%which_vert /= loc%which_vert)) then
-   write(errstring,*)'which_vert (',loc%which_vert,') must be same in all args'
+if ((minl%which_azm /= maxl%which_azm) .or. &
+    (minl%which_azm /= loc%which_azm)) then
+   write(errstring,*)'which_azm (',loc%which_azm,') must be same in all args'
    call error_handler(E_ERR, 'is_location_in_region', errstring, source, revision, revdate)
 endif
 
@@ -789,10 +770,13 @@
 ! set to success only at the bottom after all tests have passed.
 is_location_in_region = .false.
 
-! use the code in the utils module that knows how to wrap longitude/radians.
-if (.not. is_longitude_between(loc%azm, minl%azm, maxl%azm, doradians=.true.)) return
+if (loc%which_azm /= AZMISLON) then
+   if ((loc%azm < minl%azm) .or. (loc%azm > maxl%azm)) return
+else
+   ! use the code in the utils module that knows how to wrap longitude/radians.
+   if (.not. is_longitude_between(loc%azm, minl%azm, maxl%azm, doradians=.true.)) return
+endif
 if ((loc%rad < minl%rad) .or. (loc%rad > maxl%rad)) return
-if ((loc%vloc < minl%vloc) .or. (loc%vloc > maxl%vloc)) return
  
 is_location_in_region = .true.
 
@@ -815,18 +799,14 @@
 
 function vert_is_surface(loc)
  
-! Given a location, return true if vertical coordinate is surface, else false.
+! Stub, always returns false.
 
 logical                          :: vert_is_surface
 type(location_type), intent(in)  :: loc
 
 if ( .not. module_initialized ) call initialize_module
 
-if(loc%which_vert == VERTISSURFACE ) then
-   vert_is_surface = .true.
-else
-   vert_is_surface = .false.
-endif
+vert_is_surface = .false.
 
 end function vert_is_surface
 
@@ -834,7 +814,7 @@
 
 function vert_is_pressure(loc)
  
-! Always returns false, as vertical coordinate is never pressure for the annulus.
+! Stub, always returns false.
 
 logical                          :: vert_is_pressure
 type(location_type), intent(in)  :: loc
@@ -847,18 +827,14 @@
 
 function vert_is_height(loc)
  
-! Given a location, return true if vertical coordinate is height, else false.
+! Stub, always returns false.
 
 logical                          :: vert_is_height
 type(location_type), intent(in)  :: loc
 
 if ( .not. module_initialized ) call initialize_module
 
-if(loc%which_vert == VERTISHEIGHT ) then
-   vert_is_height = .true.
-else
-   vert_is_height = .false.
-endif
+vert_is_height = .false.
 
 end function vert_is_height
 
@@ -866,18 +842,14 @@
 
 function vert_is_level(loc)
  
-! Given a location, return true if vertical coordinate is level, else false.
+! Stub, always returns false.
 
 logical                          :: vert_is_level
 type(location_type), intent(in)  :: loc
 
 if ( .not. module_initialized ) call initialize_module
 
-if(loc%which_vert == VERTISLEVEL ) then
-   vert_is_level = .true.
-else
-   vert_is_level = .false.
-endif
+vert_is_level = .false.
 
 end function vert_is_level
 
@@ -898,7 +870,7 @@
 
 
 !----------------------------------------------------------------------------
-! end of location/annulus/location_mod.f90
+! end of location/twod_annulus/location_mod.f90
 !----------------------------------------------------------------------------
 
 end module location_mod

Modified: DART/branches/development/location/twod_annulus/test/input.nml
===================================================================
--- DART/branches/development/location/annulus/test/input.nml	2012-04-11 17:28:19 UTC (rev 5692)
+++ DART/branches/development/location/twod_annulus/test/input.nml	2012-04-12 21:29:58 UTC (rev 5695)
@@ -1,4 +1,6 @@
 &location_nml
+ min_radius =    100.0
+ max_radius = 100000.0
  /
 
 &utilities_nml

Modified: DART/branches/development/location/twod_annulus/test/path_names_location_test
===================================================================
--- DART/branches/development/location/annulus/test/path_names_location_test	2012-04-11 17:28:19 UTC (rev 5692)
+++ DART/branches/development/location/twod_annulus/test/path_names_location_test	2012-04-12 21:29:58 UTC (rev 5695)
@@ -3,5 +3,5 @@
 random_seq/random_seq_mod.f90
 time_manager/time_manager_mod.f90
 utilities/utilities_mod.f90
-location/annulus/location_mod.f90
+location/twod_annulus/location_mod.f90
 location/location_test.f90

Modified: DART/branches/development/location/twod_annulus/test/test.in
===================================================================
--- DART/branches/development/location/annulus/test/test.in	2012-04-11 17:28:19 UTC (rev 5692)
+++ DART/branches/development/location/twod_annulus/test/test.in	2012-04-12 21:29:58 UTC (rev 5695)
@@ -1,96 +1,66 @@
-3
--300
-240
-100
-3
--7000
+0
+0
+45
+1000
+0
 -1
-240
-40
+0
+90
 1000
-12000
-3
--30
+10000 
+0
+0
+30
+10000
+0
+0
+33
+1001
+0
+0
+50
+3000
+0
+0
+-89
+1000
+0
 -1
-140
-240
-10
-20
-3
--5000
+0
+90
+100
+10000
+0
+0
+45
+3456
+0
 -1
-340
-359
+0
+90
+100
 10000
-20000
-3
--2000
+0
 -1
-140
-250
-1
-20
-3
--200
+0
+90
+100
+10000
+0
 -1
-15
-25
+0
+90
+1000
 10000
-100000
-3
--8000
+0
 -1
-140
-250
-10
-200
-3
--100
--1
-130
-250
-10
-200
-3
--9000
--1
-240
-250
-51
-60
--1
-998
--1
-50
-60
+0
+90
+100
 10000
-10200
-1
-10
+0
 -1
-40
-50
-600
-900
-3
--10
--1
-240
-250
-10
-50
-1
-15
--1
-20
-25
-10
-20
-3
--100
--1
-40
-250
-10
-200
-
+0
+90
+100 
+10000


More information about the Dart-dev mailing list