metgrid/src/mpas_mesh.F

   1 module mpas_mesh
   2
   3     use scan_input
   4
   5     type mpas_mesh_type
   6         logical :: valid = .false.
   7         integer :: nCells = 0
   8         integer :: nVertices = 0
   9         integer :: nEdges = 0
  10         integer :: maxEdges = 0
  11         integer, dimension(:), pointer :: landmask => null()
  12         integer, dimension(:), pointer :: nEdgesOnCell => null()
  13         integer, dimension(:,:), pointer :: cellsOnCell => null()
  14         integer, dimension(:,:), pointer :: verticesOnCell => null()
  15         integer, dimension(:,:), pointer :: cellsOnVertex => null()
  16         integer, dimension(:,:), pointer :: edgesOnCell => null()
  17         integer, dimension(:,:), pointer :: cellsOnEdge => null()
  18         real, dimension(:), pointer :: latCell => null()
  19         real, dimension(:), pointer :: lonCell => null()
  20         real, dimension(:), pointer :: latVertex => null()
  21         real, dimension(:), pointer :: lonVertex => null()
  22         real, dimension(:), pointer :: latEdge => null()
  23         real, dimension(:), pointer :: lonEdge => null()
  24     end type mpas_mesh_type
  25
  26
  27     contains
  28
  29
  30     integer function mpas_mesh_setup(mesh_filename, mesh) result(stat)
  31
  32         implicit none
  33
  34         character (len=*), intent(in) :: mesh_filename
  35         type (mpas_mesh_type), intent(out) :: mesh
  36
  37         type (input_handle_type) :: handle
  38         type (input_field_type) :: field
  39
  40         stat = 0
  41
  42         if (scan_input_open(mesh_filename, handle) /= 0) then
  43             stat = 1
  44             return
  45         end if
  46
  47         !
  48         ! nEdgesOnCell
  49         !
  50         if (scan_input_for_field(handle, 'nEdgesOnCell', field) /= 0) then
  51             stat = 1
  52             return
  53         end if
  54
  55         mesh % nCells = field % dimlens(1)
  56         mesh % nVertices = 2 * (mesh % nCells - 2)
  57         mesh % nEdges = 3 * (mesh % nCells - 2)
  58
  59         stat = scan_input_read_field(field)
  60         allocate(mesh % nEdgesOnCell(mesh % nCells))
  61         mesh % nEdgesOnCell(:) = field % array1i(:)
  62         stat = scan_input_free_field(field)
  63
  64         !
  65         ! landmask
  66         !
  67         if (scan_input_for_field(handle, 'landmask', field) == 0) then
  68             stat = scan_input_read_field(field)
  69             allocate(mesh % landmask(mesh % nCells))
  70             mesh % landmask(:) = field % array1i(:)
  71             stat = scan_input_free_field(field)
  72         else     ! no landmask available
  73             allocate(mesh % landmask(mesh % nCells))
  74             mesh % landmask(:) = 1
  75         end if
  76
  77         !
  78         ! cellsOnCell
  79         !
  80         if (scan_input_for_field(handle, 'cellsOnCell', field) /= 0) then
  81             stat = 1
  82             return
  83         end if
  84
  85         mesh % maxEdges = field % dimlens(1)
  86
  87         stat = scan_input_read_field(field)
  88         allocate(mesh % cellsOnCell(mesh % maxEdges, mesh % nCells))
  89         mesh % cellsOnCell(:,:) = field % array2i(:,:)
  90         stat = scan_input_free_field(field)
  91
  92         !
  93         ! verticesOnCell
  94         !
  95         if (scan_input_for_field(handle, 'verticesOnCell', field) /= 0) then
  96             stat = 1
  97             return
  98         end if
  99
 100         stat = scan_input_read_field(field)
 101         allocate(mesh % verticesOnCell(mesh % maxEdges, mesh % nCells))
 102         mesh % verticesOnCell(:,:) = field % array2i(:,:)
 103         stat = scan_input_free_field(field)
 104
 105         !
 106         ! cellsOnVertex
 107         !
 108         if (scan_input_for_field(handle, 'cellsOnVertex', field) /= 0) then
 109             stat = 1
 110             return
 111         end if
 112
 113         stat = scan_input_read_field(field)
 114         allocate(mesh % cellsOnVertex(3, mesh % nVertices))
 115         mesh % cellsOnVertex(:,:) = field % array2i(:,:)
 116         stat = scan_input_free_field(field)
 117
 118         !
 119         ! edgesOnCell
 120         !
 121         if (scan_input_for_field(handle, 'edgesOnCell', field) /= 0) then
 122             stat = 1
 123             return
 124         end if
 125
 126         stat = scan_input_read_field(field)
 127         allocate(mesh % edgesOnCell(mesh % maxEdges, mesh % nCells))
 128         mesh % edgesOnCell(:,:) = field % array2i(:,:)
 129         stat = scan_input_free_field(field)
 130
 131         !
 132         ! cellsOnEdge
 133         !
 134         if (scan_input_for_field(handle, 'cellsOnEdge', field) /= 0) then
 135             stat = 1
 136             return
 137         end if
 138
 139         stat = scan_input_read_field(field)
 140         allocate(mesh % cellsOnEdge(2, mesh % nEdges))
 141         mesh % cellsOnEdge(:,:) = field % array2i(:,:)
 142         stat = scan_input_free_field(field)
 143
 144         !
 145         ! latCell
 146         !
 147         if (scan_input_for_field(handle, 'latCell', field) /= 0) then
 148             stat = 1
 149             return
 150         end if
 151
 152         stat = scan_input_read_field(field)
 153         allocate(mesh % latCell(mesh % nCells))
 154         if (field % xtype == FIELD_TYPE_REAL) then
 155             mesh % latCell(:) = field % array1r(:)
 156         else if (field % xtype == FIELD_TYPE_DOUBLE) then
 157             mesh % latCell(:) = real(field % array1d(:))
 158         end if
 159         stat = scan_input_free_field(field)
 160
 161         !
 162         ! lonCell
 163         !
 164         if (scan_input_for_field(handle, 'lonCell', field) /= 0) then
 165             stat = 1
 166             return
 167         end if
 168
 169         stat = scan_input_read_field(field)
 170         allocate(mesh % lonCell(mesh % nCells))
 171         if (field % xtype == FIELD_TYPE_REAL) then
 172             mesh % lonCell(:) = field % array1r(:)
 173         else if (field % xtype == FIELD_TYPE_DOUBLE) then
 174             mesh % lonCell(:) = real(field % array1d(:))
 175         end if
 176         stat = scan_input_free_field(field)
 177
 178         !
 179         ! latVertex
 180         !
 181         if (scan_input_for_field(handle, 'latVertex', field) /= 0) then
 182             stat = 1
 183             return
 184         end if
 185
 186         stat = scan_input_read_field(field)
 187         allocate(mesh % latVertex(mesh % nVertices))
 188         if (field % xtype == FIELD_TYPE_REAL) then
 189             mesh % latVertex(:) = field % array1r(:)
 190         else if (field % xtype == FIELD_TYPE_DOUBLE) then
 191             mesh % latVertex(:) = real(field % array1d(:))
 192         end if
 193         stat = scan_input_free_field(field)
 194
 195         !
 196         ! lonVertex
 197         !
 198         if (scan_input_for_field(handle, 'lonVertex', field) /= 0) then
 199             stat = 1
 200             return
 201         end if
 202
 203         stat = scan_input_read_field(field)
 204         allocate(mesh % lonVertex(mesh % nVertices))
 205         if (field % xtype == FIELD_TYPE_REAL) then
 206             mesh % lonVertex(:) = field % array1r(:)
 207         else if (field % xtype == FIELD_TYPE_DOUBLE) then
 208             mesh % lonVertex(:) = real(field % array1d(:))
 209         end if
 210         stat = scan_input_free_field(field)
 211
 212         !
 213         ! latEdge
 214         !
 215         if (scan_input_for_field(handle, 'latEdge', field) /= 0) then
 216             stat = 1
 217             return
 218         end if
 219
 220         stat = scan_input_read_field(field)
 221         allocate(mesh % latEdge(mesh % nEdges))
 222         if (field % xtype == FIELD_TYPE_REAL) then
 223             mesh % latEdge(:) = field % array1r(:)
 224         else if (field % xtype == FIELD_TYPE_DOUBLE) then
 225             mesh % latEdge(:) = real(field % array1d(:))
 226         end if
 227         stat = scan_input_free_field(field)
 228
 229         !
 230         ! lonEdge
 231         !
 232         if (scan_input_for_field(handle, 'lonEdge', field) /= 0) then
 233             stat = 1
 234             return
 235         end if
 236
 237         stat = scan_input_read_field(field)
 238         allocate(mesh % lonEdge(mesh % nEdges))
 239         if (field % xtype == FIELD_TYPE_REAL) then
 240             mesh % lonEdge(:) = field % array1r(:)
 241         else if (field % xtype == FIELD_TYPE_DOUBLE) then
 242             mesh % lonEdge(:) = real(field % array1d(:))
 243         end if
 244         stat = scan_input_free_field(field)
 245
 246         stat = scan_input_close(handle)
 247
 248         mesh % valid = .true.
 249
 250     end function mpas_mesh_setup
 251
 252
 253     integer function mpas_mesh_free(mesh) result(stat)
 254
 255         implicit none
 256
 257         type (mpas_mesh_type), intent(inout) :: mesh
 258
 259
 260         stat = 0
 261
 262         mesh % valid = .false.
 263         mesh % nCells = 0
 264         mesh % nVertices = 0
 265         mesh % nEdges = 0
 266         mesh % maxEdges = 0
 267
 268         if (associated(mesh % landmask)) then
 269             deallocate(mesh % landmask)
 270         end if
 271         if (associated(mesh % nEdgesOnCell)) then
 272             deallocate(mesh % nEdgesOnCell)
 273         end if
 274         if (associated(mesh % cellsOnCell)) then
 275             deallocate(mesh % cellsOnCell)
 276         end if
 277         if (associated(mesh % verticesOnCell)) then
 278             deallocate(mesh % verticesOnCell)
 279         end if
 280         if (associated(mesh % cellsOnVertex)) then
 281             deallocate(mesh % cellsOnVertex)
 282         end if
 283         if (associated(mesh % edgesOnCell)) then
 284             deallocate(mesh % edgesOnCell)
 285         end if
 286         if (associated(mesh % cellsOnEdge)) then
 287             deallocate(mesh % cellsOnEdge)
 288         end if
 289         if (associated(mesh % latCell)) then
 290             deallocate(mesh % latCell)
 291         end if
 292         if (associated(mesh % lonCell)) then
 293             deallocate(mesh % lonCell)
 294         end if
 295         if (associated(mesh % latVertex)) then
 296             deallocate(mesh % latVertex)
 297         end if
 298         if (associated(mesh % lonVertex)) then
 299             deallocate(mesh % lonVertex)
 300         end if
 301         if (associated(mesh % latEdge)) then
 302             deallocate(mesh % latEdge)
 303         end if
 304         if (associated(mesh % lonEdge)) then
 305             deallocate(mesh % lonEdge)
 306         end if
 307
 308     end function mpas_mesh_free
 309
 310 end module mpas_mesh