mosaic2/mosaic2.F90

   1 !***********************************************************************
   2 !*                   GNU Lesser General Public License
   3 !*
   4 !* This file is part of the GFDL Flexible Modeling System (FMS).
   5 !*
   6 !* FMS is free software: you can redistribute it and/or modify it under
   7 !* the terms of the GNU Lesser General Public License as published by
   8 !* the Free Software Foundation, either version 3 of the License, or (at
   9 !* your option) any later version.
  10 !*
  11 !* FMS is distributed in the hope that it will be useful, but WITHOUT
  12 !* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  13 !* FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  14 !* for more details.
  15 !*
  16 !* You should have received a copy of the GNU Lesser General Public
  17 !* License along with FMS.  If not, see <http://www.gnu.org/licenses/>.
  18 !***********************************************************************
  19 !> @defgroup mosaic2_mod mosaic2_mod
  20 !> @ingroup mosaic2
  21 !> @brief Implements some utility routines to read mosaic information.
  22 !!
  23 !> @author Zhi Liang
  24 !!
  25 !! Implements some utility routines to read mosaic information.
  26 !! The information includes number of tiles and contacts in the mosaic,
  27 !! mosaic grid resolution of each tile, mosaic contact information, mosaic exchange
  28 !! grid information. Each routine will call a C-version routine to get these information.
  29
  30 !> @addtogroup mosaic2_mod
  31 !> @{
  32 module mosaic2_mod
  33
  34 !use fms_mod,    only : write_version_number
  35
  36 use mpp_mod,    only : mpp_error, FATAL, mpp_pe, mpp_root_pe
  37 use mpp_domains_mod, only : domain2D, mpp_get_current_ntile, mpp_get_tile_id
  38 use constants_mod, only : PI, RADIUS
  39 use fms2_io_mod,   only : FmsNetcdfFile_t, open_file, close_file, get_dimension_size
  40 use fms2_io_mod,   only : read_data, variable_exists
  41 use platform_mod,  only : r4_kind, r8_kind, FMS_PATH_LEN
  42
  43 implicit none
  44 private
  45
  46 character(len=*), parameter :: &
  47      grid_dir  = 'INPUT/'      !> root directory for all grid files
  48
  49 integer, parameter :: &
  50      MAX_NAME = 256,  & !> max length of the variable names
  51      X_REFINE = 2,    & !> supergrid size/model grid size in x-direction
  52      Y_REFINE = 2       !> supergrid size/model grid size in y-direction
  53
  54 ! --- public interface
  55
  56
  57 public :: get_mosaic_ntiles
  58 public :: get_mosaic_ncontacts
  59 public :: get_mosaic_grid_sizes
  60 public :: get_mosaic_contact
  61 public :: get_mosaic_xgrid_size
  62 public :: get_mosaic_xgrid
  63 public :: get_mosaic_tile_grid
  64 public :: calc_mosaic_grid_area
  65 public :: calc_mosaic_grid_great_circle_area
  66 public :: is_inside_polygon
  67
  68
  69 interface get_mosaic_xgrid
  70    module procedure get_mosaic_xgrid_r4
  71    module procedure get_mosaic_xgrid_r8
  72 end interface get_mosaic_xgrid
  73
  74 interface calc_mosaic_grid_area
  75    module procedure calc_mosaic_grid_area_r4
  76    module procedure calc_mosaic_grid_area_r8
  77 end interface calc_mosaic_grid_area
  78
  79 interface calc_mosaic_grid_great_circle_area
  80    module procedure calc_mosaic_grid_great_circle_area_r4
  81    module procedure calc_mosaic_grid_great_circle_area_r8
  82 end interface calc_mosaic_grid_great_circle_area
  83
  84 interface is_inside_polygon
  85    module procedure is_inside_polygon_r4
  86    module procedure is_inside_polygon_r8
  87 end interface is_inside_polygon
  88
  89
  90 logical :: module_is_initialized = .true.
  91 ! Include variable "version" to be written to log file.
  92 #include<file_version.h>
  93
  94 contains
  95
  96 !#######################################################################
  97
  98 !> @brief Initialize the mosaic_mod.
  99 !> Initialization routine for the mosaic module. It writes the
 100 !! version information to the log file.
 101 !! <br>Example usage:
 102 !!     call mosaic_init ( )
 103 subroutine mosaic_init()
 104
 105   if (module_is_initialized) return
 106   module_is_initialized = .TRUE.
 107
 108 !--------- write version number and namelist ------------------
 109 !  call write_version_number("MOSAIC_MOD", version)
 110
 111 end subroutine mosaic_init
 112
 113 !#######################################################################
 114 !> @brief return exchange grid size of mosaic xgrid file.
 115 !!
 116 !> <br>Example usage:
 117 !!              nxgrid = get_mosaic_xgrid_size(xgrid_file)
 118   function get_mosaic_xgrid_size(fileobj)
 119     type(FmsNetcdfFile_t), intent(in) :: fileobj !> File that contains exchange grid information.
 120     integer                           :: get_mosaic_xgrid_size
 121
 122     call get_dimension_size(fileobj, "ncells", get_mosaic_xgrid_size)
 123
 124     return
 125
 126   end function get_mosaic_xgrid_size
 127   !###############################################################################
 128   !> Get number of tiles in the mosaic_file.
 129   !> @param fileobj mosaic file object
 130   !> @return Number of tiles in given file
 131   !> <br>Example usage:
 132   !!     ntiles = get_mosaic_ntiles( mosaic_file)
 133   !!
 134   function get_mosaic_ntiles(fileobj)
 135     type(FmsNetcdfFile_t), intent(in) :: fileobj
 136     integer                           :: get_mosaic_ntiles
 137
 138     call get_dimension_size(fileobj, "ntiles", get_mosaic_ntiles)
 139
 140     return
 141
 142   end function get_mosaic_ntiles
 143
 144   !###############################################################################
 145   !> Get number of contacts in the mosaic_file.
 146   !> @param fileobj mosaic file object
 147   !> @return number of contacts in a given file
 148   !> <br>Example usage:
 149   !!         ntiles = get_mosaic_ncontacts( mosaic_file)
 150   function get_mosaic_ncontacts(fileobj)
 151     type(FmsNetcdfFile_t), intent(in) :: fileobj
 152     integer                           :: get_mosaic_ncontacts
 153
 154     if(variable_exists(fileobj, "contacts") ) then
 155       call get_dimension_size(fileobj, "ncontact", get_mosaic_ncontacts)
 156     else
 157       get_mosaic_ncontacts = 0
 158     endif
 159
 160     return
 161
 162   end function get_mosaic_ncontacts
 163
 164
 165   !###############################################################################
 166   !> Get grid size of each tile from mosaic_file
 167   !> @param fileobj mosaic file object
 168   !> @param[inout] nx List of grid size in x-direction of each tile
 169   !> @param[inout] ny List of grid size in y-direction of each tile
 170   !> <br>Example usage:
 171   !!            call get_mosaic_grid_sizes(mosaic_file, nx, ny)
 172   subroutine get_mosaic_grid_sizes( fileobj, nx, ny)
 173     type(FmsNetcdfFile_t), intent(in)    :: fileobj
 174     integer, dimension(:), intent(inout) :: nx, ny
 175
 176     character(len=FMS_PATH_LEN) :: gridfile
 177     integer                 :: ntiles, n
 178     type(FmsNetcdfFile_t)   :: gridobj
 179
 180     ntiles = get_mosaic_ntiles(fileobj)
 181     if(ntiles .NE. size(nx(:)) .OR. ntiles .NE. size(ny(:)) ) then
 182       call mpp_error(FATAL, "get_mosaic_grid_sizes: size of nx/ny does not equal to ntiles")
 183     endif
 184
 185     do n = 1, ntiles
 186       call read_data(fileobj, 'gridfiles', gridfile, corner=n)
 187       gridfile = grid_dir//trim(gridfile)
 188
 189       if(.not. open_file(gridobj, gridfile, 'read')) then
 190          call mpp_error(FATAL, 'mosaic_mod(get_mosaic_grid_sizes):Error in opening file '//trim(gridfile))
 191       endif
 192       call get_dimension_size(gridobj, "nx", nx(n))
 193       call get_dimension_size(gridobj, "ny", ny(n))
 194       call close_file(gridobj)
 195       if(mod(nx(n),x_refine) .NE. 0) call mpp_error(FATAL, "get_mosaic_grid_sizes: nx is not divided by x_refine");
 196       if(mod(ny(n),y_refine) .NE. 0) call mpp_error(FATAL, "get_mosaic_grid_sizes: ny is not divided by y_refine");
 197       nx(n) = nx(n)/x_refine;
 198       ny(n) = ny(n)/y_refine;
 199     enddo
 200
 201     return
 202
 203   end subroutine get_mosaic_grid_sizes
 204
 205   !###############################################################################
 206   !> Get contact information from mosaic_file
 207   !> <br>Example usage:
 208   !!         call get_mosaic_contact(mosaic_file, tile1, tile2, istart1, iend1, jstart1, jend1,
 209   !!                                 istart2, iend2, jstart2, jend2)
 210   !> @param fileobj mosaic file object
 211   !> @param[inout] tile1 list of tile numbers in tile 1 of each contact
 212   !> @param[inout] tile2 list of tile numbers in tile 2 of each contact
 213   !> @param[inout] istart1 list of starting i-index in tile 1 of each contact
 214   !> @param[inout] iend1 list of ending i-index in tile 1 of each contact
 215   !> @param[inout] jstart1 list of starting j-index in tile 1 of each contact
 216   !> @param[inout] jend1 list of ending j-index in tile 1 of each contact
 217   !> @param[inout] istart2 list of starting i-index in tile 2 of each contact
 218   !> @param[inout] iend2 list of ending i-index in tile 2 of each contact
 219   !> @param[inout] jstart2 list of starting j-index in tile 2 of each contact
 220   !> @param[inout] jend2 list of ending j-index in tile 2 of each contact
 221   subroutine get_mosaic_contact( fileobj, tile1, tile2, istart1, iend1, jstart1, jend1, &
 222                                    istart2, iend2, jstart2, jend2)
 223     type(FmsNetcdfFile_t),    intent(in) :: fileobj
 224     integer, dimension(:), intent(inout) :: tile1, tile2
 225     integer, dimension(:), intent(inout) :: istart1, iend1, jstart1, jend1
 226     integer, dimension(:), intent(inout) :: istart2, iend2, jstart2, jend2
 227     character(len=MAX_NAME), allocatable :: gridtiles(:)
 228     character(len=MAX_NAME), allocatable :: contacts(:)
 229     character(len=MAX_NAME), allocatable :: contacts_index(:)
 230     character(len=MAX_NAME)              :: strlist(8)
 231     integer :: ntiles, n, m, ncontacts, nstr, ios
 232     integer :: i1_type, j1_type, i2_type, j2_type
 233     logical :: found
 234
 235     ntiles = get_mosaic_ntiles(fileobj)
 236     allocate(gridtiles(ntiles))
 237     if(mpp_pe()==mpp_root_pe()) then
 238       do n = 1, ntiles
 239         do m = 1,MAX_NAME
 240           gridtiles(n)(m:m) = " "
 241         enddo
 242       enddo
 243     endif
 244     call read_data(fileobj, 'gridtiles', gridtiles)
 245
 246     ncontacts = get_mosaic_ncontacts(fileobj)
 247     if(ncontacts>0) then
 248        allocate(contacts(ncontacts), contacts_index(ncontacts))
 249        if(mpp_pe()==mpp_root_pe()) then
 250          do n = 1, ncontacts
 251            do m = 1,MAX_NAME
 252              contacts(n)(m:m) = " "
 253              contacts_index(n)(m:m) = " "
 254            enddo
 255          enddo
 256        endif
 257        call read_data(fileobj, "contacts", contacts)
 258        call read_data(fileobj, "contact_index", contacts_index)
 259     endif
 260     do n = 1, ncontacts
 261       nstr = parse_string(contacts(n), ":", strlist)
 262       if(nstr .NE. 4) call mpp_error(FATAL, &
 263          "mosaic_mod(get_mosaic_contact): number of elements in contact seperated by :/:: should be 4")
 264       found = .false.
 265       do m = 1, ntiles
 266         if(trim(gridtiles(m)) == trim(strlist(2)) ) then !found the tile name
 267           found = .true.
 268           tile1(n) = m
 269           exit
 270         endif
 271       enddo
 272
 273       if(.not.found) call mpp_error(FATAL, &
 274          "mosaic_mod(get_mosaic_contact):the first tile name specified in contact is not found in tile list")
 275
 276       found = .false.
 277       do m = 1, ntiles
 278         if(trim(gridtiles(m)) == trim(strlist(4)) ) then !found the tile name
 279           found = .true.
 280           tile2(n) = m
 281           exit
 282         endif
 283       enddo
 284
 285       if(.not.found) call mpp_error(FATAL, &
 286          "mosaic_mod(get_mosaic_contact):the second tile name specified in contact is not found in tile list")
 287
 288       nstr = parse_string(contacts_index(n), ":,", strlist)
 289       if(nstr .NE. 8) then
 290         if(mpp_pe()==mpp_root_pe()) then
 291           print*, "nstr is ", nstr
 292           print*, "contacts is ", contacts_index(n)
 293           do m = 1, nstr
 294             print*, "strlist is ", trim(strlist(m))
 295           enddo
 296         endif
 297         call mpp_error(FATAL, &
 298                "mosaic_mod(get_mosaic_contact): number of elements in contact_index seperated by :/, should be 8")
 299       endif
 300       read(strlist(1), *, iostat=ios) istart1(n)
 301       if(ios .NE. 0) call mpp_error(FATAL, &
 302          "mosaic_mod(get_mosaic_contact): Error in reading istart1")
 303       read(strlist(2), *, iostat=ios) iend1(n)
 304       if(ios .NE. 0) call mpp_error(FATAL, &
 305          "mosaic_mod(get_mosaic_contact): Error in reading iend1")
 306       read(strlist(3), *, iostat=ios) jstart1(n)
 307       if(ios .NE. 0) call mpp_error(FATAL, &
 308          "mosaic_mod(get_mosaic_contact): Error in reading jstart1")
 309       read(strlist(4), *, iostat=ios) jend1(n)
 310       if(ios .NE. 0) call mpp_error(FATAL, &
 311          "mosaic_mod(get_mosaic_contact): Error in reading jend1")
 312       read(strlist(5), *, iostat=ios) istart2(n)
 313       if(ios .NE. 0) call mpp_error(FATAL, &
 314          "mosaic_mod(get_mosaic_contact): Error in reading istart2")
 315       read(strlist(6), *, iostat=ios) iend2(n)
 316       if(ios .NE. 0) call mpp_error(FATAL, &
 317          "mosaic_mod(get_mosaic_contact): Error in reading iend2")
 318       read(strlist(7), *, iostat=ios) jstart2(n)
 319       if(ios .NE. 0) call mpp_error(FATAL, &
 320          "mosaic_mod(get_mosaic_contact): Error in reading jstart2")
 321       read(strlist(8), *, iostat=ios) jend2(n)
 322       if(ios .NE. 0) call mpp_error(FATAL, &
 323          "mosaic_mod(get_mosaic_contact): Error in reading jend2")
 324
 325       i1_type = transfer_to_model_index(istart1(n), iend1(n), x_refine)
 326       j1_type = transfer_to_model_index(jstart1(n), jend1(n), y_refine)
 327       i2_type = transfer_to_model_index(istart2(n), iend2(n), x_refine)
 328       j2_type = transfer_to_model_index(jstart2(n), jend2(n), y_refine)
 329
 330       if( i1_type == 0 .AND. j1_type == 0 ) call mpp_error(FATAL, &
 331          "mosaic_mod(get_mosaic_contact): istart1==iend1 and jstart1==jend1")
 332       if( i2_type == 0 .AND. j2_type == 0 ) call mpp_error(FATAL, &
 333          "mosaic_mod(get_mosaic_contact): istart2==iend2 and jstart2==jend2")
 334       if( i1_type + j1_type .NE. i2_type + j2_type ) call mpp_error(FATAL, &
 335          "mosaic_mod(get_mosaic_contact): It is not a line or overlap contact")
 336
 337    enddo
 338
 339    deallocate(gridtiles)
 340    if(ncontacts>0) deallocate(contacts, contacts_index)
 341
 342   end subroutine get_mosaic_contact
 343
 344
 345 function transfer_to_model_index(istart, iend, refine_ratio)
 346    integer, intent(inout) :: istart, iend
 347    integer                :: refine_ratio
 348    integer                :: transfer_to_model_index
 349    integer                :: istart_in, iend_in
 350
 351    istart_in = istart
 352    iend_in = iend
 353
 354    if( istart_in == iend_in ) then
 355       transfer_to_model_index = 0
 356       istart = (istart_in + 1)/refine_ratio
 357       iend   = istart
 358    else
 359       transfer_to_model_index = 1
 360       if( iend_in > istart_in ) then
 361         istart = istart_in + 1
 362         iend   = iend_in
 363       else
 364         istart = istart_in
 365         iend   = iend_in + 1
 366       endif
 367       if( mod(istart, refine_ratio) .NE. 0 .OR. mod(iend,refine_ratio) .NE. 0) call mpp_error(FATAL, &
 368          "mosaic_mod(transfer_to_model_index): mismatch between refine_ratio and istart/iend")
 369       istart = istart/refine_ratio
 370       iend = iend/refine_ratio
 371
 372    endif
 373
 374    return
 375
 376 end function transfer_to_model_index
 377 !#####################################################################
 378 function parse_string(string, set, sval)
 379    character(len=*),  intent(in) :: string
 380    character(len=*),  intent(in) :: set
 381    character(len=*), intent(out) :: sval(:)
 382    integer                       :: parse_string
 383    integer :: nelem, length, first, last
 384
 385    nelem = size(sval(:))
 386    length = len_trim(string)
 387
 388    first = 1; last = 0
 389    parse_string = 0
 390
 391    do while(first .LE. length)
 392       parse_string = parse_string + 1
 393       if(parse_string>nelem) then
 394          call mpp_error(FATAL, "mosaic_mod(parse_string) : number of element is greater than size(sval(:))")
 395       endif
 396       last = first - 1 + scan(string(first:length), set)
 397       if(last == first-1 ) then  ! not found, end of string
 398          sval(parse_string) = string(first:length)
 399          exit
 400       else
 401          if(last <= first) then
 402             call mpp_error(FATAL, "mosaic_mod(parse_string) : last <= first")
 403          endif
 404          sval(parse_string) = string(first:(last-1))
 405          first = last + 1
 406          ! scan to make sure the next is not the character in the set
 407          do while (first == last+1)
 408             last = first - 1 + scan(string(first:length), set)
 409             if(last == first) then
 410                first = first+1
 411             else
 412                exit
 413             endif
 414          end do
 415       endif
 416    enddo
 417
 418    return
 419
 420 end function parse_string
 421
 422   !#############################################################################
 423   !> Gets the name of a mosaic tile grid file
 424   !> @param[out] grid_file name of grid file
 425   !> @param fileobj mosaic file object
 426   !> @param domain current domain
 427   !> @param tile_count optional count of tiles
 428 subroutine get_mosaic_tile_grid(grid_file, fileobj, domain, tile_count)
 429    character(len=*), intent(out)          :: grid_file
 430    type(FmsNetcdfFile_t), intent(in)      :: fileobj
 431    type(domain2D),   intent(in)           :: domain
 432    integer,          intent(in), optional :: tile_count
 433    integer                                :: tile, ntileMe
 434    integer, dimension(:), allocatable     :: tile_id
 435    character(len=256), allocatable        :: filelist(:)
 436    integer :: ntiles
 437
 438    ntiles = get_mosaic_ntiles(fileobj)
 439    allocate(filelist(ntiles))
 440    tile = 1
 441    if(present(tile_count)) tile = tile_count
 442    ntileMe = mpp_get_current_ntile(domain)
 443    allocate(tile_id(ntileMe))
 444    tile_id = mpp_get_tile_id(domain)
 445    call read_data(fileobj, "gridfiles", filelist)
 446    grid_file = 'INPUT/'//trim(filelist(tile_id(tile)))
 447    deallocate(tile_id, filelist)
 448
 449 end subroutine get_mosaic_tile_grid
 450
 451 #include "mosaic2_r4.fh"
 452 #include "mosaic2_r8.fh"
 453
 454 end module mosaic2_mod
 455 !> @}
 456 ! close documentation grouping