module mdl_grd_def implicit none private public :: model_grid_definition real, parameter :: ZERO = 0. real, parameter :: ZERO_8 = 0._8 real, parameter :: ONEHALF = .5 real, parameter :: ONEHALF_8 = .5_8 real, parameter :: ONE = 1. real, parameter :: ONE_8 = 1._8 real, parameter :: FOUR = 4. real, parameter :: FOUR_8 = 4._8 real, parameter :: TEN_8 = 10._8 real, parameter :: NINETY = 90. real, parameter :: ONE80 = 180. real, parameter :: ONE80_8 = 180._8 real, parameter :: THREE60 = 360. real, parameter :: THREE60_8 = 360._8 real, parameter :: ROUNDOFF = 10.*epsilon(ROUNDOFF) real :: PI, R2D character(len=16) :: nPolygons_varnm, maxVtx_varnm character(len=16) :: CenterLon_varnm, CenterLat_varnm, Area_varnm character(len=16) :: VertexLon_varnm, VertexLat_varnm include 'netcdf.inc' CONTAINS subroutine model_grid_definition( ncid, mdlType, nPolygons, maxVtx, & CenterLon, CenterLat, Area, & VertexLon, VertexLat ) use netcdf_utils !----------------------------------------------------------------------- ! dummy arguments !----------------------------------------------------------------------- integer, intent(in) :: ncid ! model netcdf file id integer, intent(out) :: nPolygons ! number of polygons in grid integer, intent(out) :: maxVtx ! max number of vertices per polygon real, allocatable, intent(out) :: CenterLon(:) real, allocatable, intent(out) :: CenterLat(:) real, allocatable, intent(out) :: Area(:) real, allocatable, intent(out) :: VertexLon(:,:) real, allocatable, intent(out) :: VertexLat(:,:) character(len=*), intent(in) :: mdlType !----------------------------------------------------------------------- ! local variables !----------------------------------------------------------------------- integer :: astat integer :: dimid, varid integer :: index, nVertex integer :: polygon, vertex integer, allocatable :: vtxOnCell(:,:) real, allocatable :: LonVtx(:), LatVtx(:) character(len=64) :: mess character(len=16) :: varnm select case( trim(mdlType) ) case( 'CAMSE' ) nPolygons_varnm = 'grid_size' maxVtx_varnm = 'grid_corners' CenterLon_varnm = 'grid_center_lon' CenterLat_varnm = 'grid_center_lat' Area_varnm = 'grid_area' VertexLon_varnm = 'grid_corner_lon' VertexLat_varnm = 'grid_corner_lat' case( 'MPAS' ) nPolygons_varnm = 'nCells' maxVtx_varnm = 'maxEdges' CenterLon_varnm = 'lonCell' CenterLat_varnm = 'latCell' Area_varnm = 'areaCell' VertexLon_varnm = 'lonVertex' VertexLat_varnm = 'latVertex' PI = FOUR*atan(ONE) R2D = ONE80/PI case default write(*,*) 'model_grid_definition: model type ' // trim(mdlType) // ' is not CAMSE or MPAS' stop 'ArgError' end select !----------------------------------------------------------------------- ! get number polygons !----------------------------------------------------------------------- varnm = trim(nPolygons_varnm) mess = 'model_grid_definition: Failed to get ' // trim(varnm) // ' dimension id' call handle_ncerr( nf_inq_dimid( ncid, trim(varnm), dimid ), mess ) mess = 'model_grid_definition: Failed to get ' // trim(varnm) // ' dimension' call handle_ncerr( nf_inq_dimlen( ncid, dimid, nPolygons ), mess ) !----------------------------------------------------------------------- ! get max number vertices/polygon !----------------------------------------------------------------------- varnm = trim(maxVtx_varnm) mess = 'model_grid_definition: Failed to get ' // trim(varnm) // ' dimension id' call handle_ncerr( nf_inq_dimid( ncid, trim(varnm), dimid ), mess ) mess = 'model_grid_definition: Failed to get ' // trim(varnm) // ' dimension' call handle_ncerr( nf_inq_dimlen( ncid, dimid, maxVtx ), mess ) !----------------------------------------------------------------------- ! allocate Center lon,lat, and area arrays !----------------------------------------------------------------------- allocate( CenterLon(nPolygons),CenterLat(nPolygons),Area(nPolygons),stat=astat ) if( astat /= 0 ) then write(*,*) 'model_grid_definition: Failed to allocate Center{Lon,Lat},Area arrays; error = ',astat stop 'AllocErr' endif !----------------------------------------------------------------------- ! get polygon center longitudes !----------------------------------------------------------------------- varnm = trim(CenterLon_varnm) mess = 'model_grid_definition: Failed to get ' // trim(varnm) // ' variable id' call handle_ncerr( nf_inq_varid( ncid, trim(varnm), varid ), mess ) mess = 'model_grid_definition: Failed to read ' // trim(varnm) call handle_ncerr( nf_get_var_real( ncid, varid, CenterLon ), mess ) !----------------------------------------------------------------------- ! get polygon center latitudes !----------------------------------------------------------------------- varnm = trim(CenterLat_varnm) mess = 'model_grid_definition: Failed to get ' // trim(varnm) // ' variable id' call handle_ncerr( nf_inq_varid( ncid, trim(varnm), varid ), mess ) mess = 'model_grid_definition: Failed to read ' // trim(varnm) call handle_ncerr( nf_get_var_real( ncid, varid, CenterLat ), mess ) !----------------------------------------------------------------------- ! get polygon areas (square radians) !----------------------------------------------------------------------- varnm = trim(Area_varnm) mess = 'model_grid_definition: Failed to get ' // trim(varnm) // ' variable id' call handle_ncerr( nf_inq_varid( ncid, trim(varnm), varid ), mess ) mess = 'model_grid_definition: Failed to read ' // trim(varnm) call handle_ncerr( nf_get_var_real( ncid, varid, Area ), mess ) CenterLon(:) = R2D*CenterLon(:) CenterLat(:) = R2D*CenterLat(:) !----------------------------------------------------------------------- ! allocate Vertex lon,lat arrays !----------------------------------------------------------------------- allocate( VertexLon(maxVtx,nPolygons),VertexLat(maxVtx,nPolygons),stat=astat ) if( astat /= 0 ) then write(*,*) 'model_grid_definition: Failed to allocate VertexLon,Lat},Area arrays; error = ',astat stop 'AllocErr' endif select case( trim(mdlType) ) case( 'CAMSE' ) !----------------------------------------------------------------------- ! get polygon vertex longitudes !----------------------------------------------------------------------- varnm = trim(VertexLon_varnm) mess = 'model_grid_definition: Failed to get ' // trim(varnm) // ' variable id' call handle_ncerr( nf_inq_varid( ncid, trim(varnm), varid ), mess ) mess = 'model_grid_definition: Failed to read ' // trim(varnm) call handle_ncerr( nf_get_var_real( ncid, varid, VertexLon ), mess ) !----------------------------------------------------------------------- ! get polygon vertex latitudes !----------------------------------------------------------------------- varnm = trim(VertexLat_varnm) mess = 'model_grid_definition: Failed to get ' // trim(varnm) // ' variable id' call handle_ncerr( nf_inq_varid( ncid, trim(varnm), varid ), mess ) mess = 'model_grid_definition: Failed to read ' // trim(varnm) call handle_ncerr( nf_get_var_real( ncid, varid, VertexLat ), mess ) case( 'MPAS' ) !----------------------------------------------------------------------- ! get number vertices !----------------------------------------------------------------------- varnm = 'nVertices' mess = 'model_grid_definition: Failed to get ' // trim(varnm) // ' dimension id' call handle_ncerr( nf_inq_dimid( ncid, trim(varnm), dimid ), mess ) mess = 'model_grid_definition: Failed to get ' // trim(varnm) // ' dimension' call handle_ncerr( nf_inq_dimlen( ncid, dimid, nVertex ), mess ) allocate( vtxOnCell(maxVtx,nPolygons),LonVtx(nVertex),LatVtx(nVertex),stat=astat ) if( astat /= 0 ) then write(*,*) 'model_grid_definition: Failed to allocate vtxOnCell,LonVtx,LatVtx arrays; error = ',astat stop 'AllocErr' endif !----------------------------------------------------------------------- ! get MPAS vertex mapping array !----------------------------------------------------------------------- varnm = 'verticesOnCell' mess = 'model_grid_definition: Failed to get ' // trim(varnm) // ' variable id' call handle_ncerr( nf_inq_varid( ncid, trim(varnm), varid ), mess ) mess = 'model_grid_definition: Failed to read ' // trim(varnm) call handle_ncerr( nf_get_var_int( ncid, varid, vtxOnCell ), mess ) !----------------------------------------------------------------------- ! get MPAS vertex longitudes !----------------------------------------------------------------------- varnm = 'lonVertex' mess = 'model_grid_definition: Failed to get ' // trim(varnm) // ' variable id' call handle_ncerr( nf_inq_varid( ncid, trim(varnm), varid ), mess ) mess = 'model_grid_definition: Failed to read ' // trim(varnm) call handle_ncerr( nf_get_var_real( ncid, varid, LonVtx ), mess ) !----------------------------------------------------------------------- ! get MPAS vertex latitudes !----------------------------------------------------------------------- varnm = 'latVertex' mess = 'model_grid_definition: Failed to get ' // trim(varnm) // ' variable id' call handle_ncerr( nf_inq_varid( ncid, trim(varnm), varid ), mess ) mess = 'model_grid_definition: Failed to read ' // trim(varnm) call handle_ncerr( nf_get_var_real( ncid, varid, LatVtx ), mess ) do polygon = 1,nPolygons do vertex = 1,maxVtx index = vtxOnCell(vertex,polygon) if( index /= 0 ) then VertexLon(vertex,polygon) = LonVtx(index) VertexLat(vertex,polygon) = LatVtx(index) else VertexLon(vertex:maxVtx,polygon) = VertexLon(vertex-1,polygon) VertexLat(vertex:maxVtx,polygon) = VertexLat(vertex-1,polygon) exit endif enddo VertexLon(:,polygon) = R2D*VertexLon(:,polygon) where( VertexLon(:,polygon) < ZERO ) VertexLon(:,polygon) = VertexLon(:,polygon) + THREE60 endwhere VertexLat(:,polygon) = R2D*VertexLat(:,polygon) enddo deallocate( vtxOnCell,LonVtx,LatVtx ) end select end subroutine model_grid_definition end module mdl_grd_def