metgrid/src/target_mesh.F

   1 module target_mesh
   2
   3     type target_mesh_type
   4         logical :: valid = .false.
   5         integer :: irank = 0
   6         integer :: nLat = 0
   7         integer :: nLon = 0
   8         real :: startLat = 0.0
   9         real :: endLat = 0.0
  10         real :: startLon = 0.0
  11         real :: endLon = 0.0
  12         real, dimension(:,:), pointer :: lats => null()
  13         real, dimension(:,:), pointer :: lons => null()
  14     end type target_mesh_type
  15
  16
  17     contains
  18
  19
  20     integer function target_mesh_setup(mesh, lat2d, lon2d) result(stat)
  21
  22         implicit none
  23
  24         type (target_mesh_type), intent(out) :: mesh
  25         real, dimension(:,:), target, optional :: lat2d
  26         real, dimension(:,:), target, optional :: lon2d
  27
  28         integer :: i, j
  29         integer :: iostatus
  30         integer :: eqIdx
  31         real :: delta
  32         logical :: exists
  33         character (len=64) :: spec
  34         real, parameter :: pi_const = 2.0 * asin(1.0)
  35
  36         stat = 0
  37
  38         !
  39         ! If 2-d arrays of latitude and longitude are provided, we can just
  40         ! point to those arrays rather than generate lat/lon values based on
  41         ! a specified target domain
  42         !
  43         if (present(lat2d) .and. present(lon2d)) then
  44
  45             mesh % irank = 1
  46             mesh % nLat = size(lat2d,2)
  47             mesh % nLon = size(lon2d,1)
  48             mesh % lats => lat2d
  49             mesh % lons => lon2d
  50             mesh % valid = .true.
  51
  52             return
  53         end if
  54
  55
  56         !
  57         ! Try to parse nLat, nLon from target_domain file
  58         !
  59         inquire(file='target_domain', exist=exists)
  60         if (exists) then
  61             write(0,*) ' '
  62             write(0,*) 'Reading target domain specification from file ''target_domain'''
  63
  64             mesh % startLat = -90.0
  65             mesh % endLat   =  90.0
  66             mesh % startLon = -180.0
  67             mesh % endLon   =  180.0
  68             mesh % nLat = 360
  69             mesh % nLon = 720
  70
  71             open(22, file='target_domain', form='formatted')
  72             read(22,fmt='(a)',iostat=iostatus) spec
  73             j = 1
  74             do while (iostatus >= 0)
  75                 call despace(spec)
  76                 eqIdx = index(spec, '=')
  77                 if (eqIdx /= 0) then
  78                     if (spec(1:eqIdx-1) == 'nlat') then
  79                         read(spec(eqIdx+1:len_trim(spec)),fmt=*) mesh % nLat
  80                         write(0,*) 'Setting nlat = ', mesh % nLat
  81                     else if (spec(1:eqIdx-1) == 'nlon') then
  82                         read(spec(eqIdx+1:len_trim(spec)),fmt=*) mesh % nLon
  83                         write(0,*) 'Setting nlon = ', mesh % nLon
  84                     else if (spec(1:eqIdx-1) == 'startlat') then
  85                         read(spec(eqIdx+1:len_trim(spec)),fmt=*) mesh % startLat
  86                         write(0,*) 'Setting startlat = ', mesh % startLat
  87                     else if (spec(1:eqIdx-1) == 'endlat') then
  88                         read(spec(eqIdx+1:len_trim(spec)),fmt=*) mesh % endLat
  89                         write(0,*) 'Setting endlat = ', mesh % endLat
  90                     else if (spec(1:eqIdx-1) == 'startlon') then
  91                         read(spec(eqIdx+1:len_trim(spec)),fmt=*) mesh % startLon
  92                         write(0,*) 'Setting startlon = ', mesh % startLon
  93                     else if (spec(1:eqIdx-1) == 'endlon') then
  94                         read(spec(eqIdx+1:len_trim(spec)),fmt=*) mesh % endLon
  95                         write(0,*) 'Setting endlon = ', mesh % endLon
  96                     else
  97                         write(0,*) 'Unrecognized keyword on line ', j, ' of file ''target_domain'': '//spec(1:eqIdx-1)
  98                         stat = 1
  99                         close(22)
 100                         return
 101                     end if
 102                 else
 103                     write(0,*) 'Syntax error on line ', j, ' of file ''target_domain'': ''='' not found'
 104                     stat = 1
 105                     close(22)
 106                     return
 107                 end if
 108                 read(22,fmt='(a)',iostat=iostatus) spec
 109                 j = j + 1
 110             end do
 111             close(22)
 112         else
 113             write(0,*) ' '
 114             write(0,*) 'Target domain specification file ''target_domain'' not found.'
 115             write(0,*) 'Default 0.5-degree global target domain will be used.'
 116             write(0,*) ' '
 117
 118             mesh % startLat = -90.0
 119             mesh % endLat   =  90.0
 120             mesh % startLon = -180.0
 121             mesh % endLon   =  180.0
 122             mesh % nLat = 360
 123             mesh % nLon = 720
 124         end if
 125
 126
 127         allocate(mesh % lats(1, mesh % nLat))
 128         allocate(mesh % lons(mesh % nLon, 1))
 129
 130         delta = (mesh % endLat - mesh % startLat) / real(mesh % nLat)
 131         do i=0,mesh % nLat-1
 132            mesh % lats(1,i+1) = mesh % startLat + 0.5 * delta + real(i) * delta
 133            mesh % lats(1,i+1) = mesh % lats(1,i+1) * pi_const / 180.0
 134         end do
 135
 136         delta = (mesh % endLon - mesh % startLon) / real(mesh % nLon)
 137         do i=0,mesh % nLon-1
 138            mesh % lons(i+1,1) = mesh % startLon + 0.5 * delta + real(i) * delta
 139            mesh % lons(i+1,1) = mesh % lons(i+1,1) * pi_const / 180.0
 140         end do
 141
 142         mesh % valid = .true.
 143
 144     end function target_mesh_setup
 145
 146
 147     integer function target_mesh_free(mesh) result(stat)
 148
 149         implicit none
 150
 151         type (target_mesh_type), intent(inout) :: mesh
 152
 153
 154         stat = 0
 155
 156         mesh % valid = .false.
 157         mesh % nLat = 0
 158         mesh % nLon = 0
 159         mesh % startLat = 0.0
 160         mesh % endLat = 0.0
 161         mesh % startLon = 0.0
 162         mesh % endLon = 0.0
 163
 164         !
 165         ! When irank == 0, we allocated the lats and lons arrays
 166         ! internally and should therefore deallocate them
 167         !
 168         if (mesh % irank == 0) then
 169            if (associated(mesh % lats)) deallocate(mesh % lats)
 170            if (associated(mesh % lons)) deallocate(mesh % lons)
 171         end if
 172
 173     end function target_mesh_free
 174
 175
 176    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 177    ! Name: despace
 178    !
 179    ! Purpose: Remove all space and tab characters from a string, thus
 180    !          compressing the string to the left.
 181    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 182    subroutine despace(string)
 183
 184        implicit none
 185
 186        ! Arguments
 187        character (len=*), intent(inout) :: string
 188
 189        ! Local variables
 190        integer :: i, j, length, iquoted
 191
 192        length = len(string)
 193
 194        iquoted = 0
 195        j = 1
 196        do i=1,length
 197            ! Check for a quote mark
 198            if (string(i:i) == '"' .or. string(i:i) == '''') iquoted = mod(iquoted+1,2)
 199
 200            ! Check for non-space, non-tab character, or if we are inside quoted
 201            ! text
 202            if ((string(i:i) /= ' ' .and. string(i:i) /= achar(9)) .or. iquoted == 1) then
 203                string(j:j) = string(i:i)
 204                j = j + 1
 205            end if
 206        end do
 207
 208        do i=j,length
 209            string(i:i) = ' '
 210        end do
 211
 212    end subroutine despace
 213
 214 end module target_mesh