psmile_apply_user_data.F90

Go to the documentation of this file.
00001 !-----------------------------------------------------------------------
00002 ! Copyright 2010, DKRZ, Hamburg, Germany.
00003 ! All rights reserved. Use is subject to OASIS4 license terms.
00004 !-----------------------------------------------------------------------
00005 !
00006 ! !DESCRIPTION:
00007 !
00008 ! The routine psmile_apply_user_data preprocesses the data collected
00009 ! by the psmile_store_data_* routines and generates the psmile internal
00010 ! data structures.
00011 !
00012 !
00013 ! !REVISION HISTORY:
00014 !
00015 !   Date      Programmer   Description
00016 ! ----------  ----------   -----------
00017 ! 27.10.10    M. Hanke     created
00018 !
00019 !----------------------------------------------------------------------
00020 !
00021 !  $Id: psmile_apply_user_data.F90 2801 2010-12-06 13:27:54Z hanke $
00022 !  $Author: hanke $
00023 !
00024 !----------------------------------------------------------------------
00025 
00026 subroutine psmile_apply_user_data(ierror)
00027 
00028    use psmile_user_data, only : grid_data, point_data,      &
00029                                 mask_data, var_data,        &
00030                                 grid_id_map, point_id_map,  &
00031                                 mask_id_map, var_id_map,    &
00032                                 get_grid_type, get_comp_id, &
00033                                 generate_partition_data,    &
00034                                 partition_data_available,   &
00035                                 reducedgrid_map_available
00036    use psmile_multimap, only  : add_pair, get_num_values, &
00037                                 get_values, init_multimap
00038    use psmile_common, only    : ch_id, ndim_1d, ndim_2d, ndim_3d
00039    use psmile_grid, only      : PRISM_Gaussreduced_regvrt,   &
00040                                 PRISM_Gaussreduced_sigmavrt, &
00041                                 PRISM_Gridless,              &
00042                                 get_size_of_shape
00043    use prism_constants, only  : PRISM_Error_Grid, prism_undefined
00044    use psmile, only           : Grids,                        &
00045                                 psmile_def_grid,              &
00046                                 psmile_def_partition,         &
00047                                 psmile_reducedgrid_map,       &
00048                                 psmile_set_corners_3d_double, &
00049                                 psmile_set_corners_3d_real,   &
00050                                 psmile_set_mask,              &
00051                                 psmile_set_points_gridless,   &
00052                                 psmile_set_points_3d_double,  &
00053                                 psmile_set_points_3d_real,    &
00054                                 psmile_def_var
00055 
00056    implicit none
00057 
00058    integer, intent (out) :: ierror
00059 
00060    integer               :: i, j
00061    integer               :: grid_type, grid_id
00062    integer               :: psmile_id
00063    integer               :: size_of_shape(2)
00064    integer, allocatable  :: block_start_index(:), block_extent(:), grid_extent(:)
00065    integer, allocatable  :: block_valid_shape(:, :)
00066    integer, allocatable  :: grid_ids(:) ! contains all grid_id''s of grid (grid can have
00067                                         ! multiple grid_ids due to the splitting into blocks)
00068    integer, allocatable  :: point_ids(:), mask_ids(:)
00069    integer               :: num_blocks ! contains the number of blocks of a grid
00070 
00071    integer               :: block_actual_shape (2, ndim_3d)
00072    integer               :: grid_actual_shape (2, ndim_3d)
00073 
00074 #ifdef VERBOSE
00075    print 9990, trim(ch_id)
00076    call psmile_flushstd
00077 #endif /* VERBOSE */
00078 
00079    ierror = 0
00080 
00081    if (associated (grid_data)) then
00082 
00083       ! for all grids defined by the user
00084       do i = 1, size(grid_data)
00085 
00086          grid_type = get_grid_type (i)
00087 
00088          ! unfortunately, some grid types require special care...:(
00089          ! The partition information of these grid types has a special meaning,
00090          ! therefore we cannot split the grid into blocks.
00091          if (grid_type == PRISM_Gaussreduced_regvrt .or.&
00092              grid_type == PRISM_Gaussreduced_sigmavrt .or. &
00093              grid_type == PRISM_Gridless) then
00094 
00095             call psmile_def_grid (psmile_id, &
00096                                   grid_data(i)%grid_name, &
00097                                   grid_data(i)%comp_id, &
00098                                   grid_data(i)%grid_valid_shape, &
00099                                   grid_data(i)%grid_type, &
00100                                   ierror)
00101 
00102             call add_pair(grid_id_map, i, psmile_id)
00103 
00104             if (.not. partition_data_available (i)) then
00105 
00106                ierror = PRISM_Error_Grid
00107 
00108                call psmile_error (PRISM_Error_Grid, 'partition information is missing', &
00109                                  (/psmile_id, get_comp_id(i), grid_type/), 3,           &
00110                                  __FILE__, __LINE__ )
00111             else
00112 
00113                call psmile_def_partition(psmile_id, size (grid_data(i)%partition_array, 1), &
00114                                          grid_data(i)%partition_array(:, :),                &
00115                                          grid_data(i)%extent_array(:, :),                   &
00116                                          ierror)
00117             endif
00118 
00119             ! in addition to psmile_def_partition, reduced gauss grid require a call
00120             ! to psmile_reducedgrid_map and set_corners
00121             if (grid_type == PRISM_Gaussreduced_regvrt .or. &
00122                grid_type == PRISM_Gaussreduced_sigmavrt) then
00123 
00124                if (.not. reducedgrid_map_available(i)) then
00125 
00126                   ierror = PRISM_Error_Grid
00127 
00128                   call psmile_error (PRISM_Error_Grid, 'prism_reducedgrid_map was not called', &
00129                                        (/psmile_id, get_comp_id(i), grid_type/), 3,            &
00130                                        __FILE__, __LINE__ )
00131                else
00132 
00133                   call psmile_reducedgrid_map (psmile_id,                                 &
00134                                                size (grid_data(i)%nbr_points_per_lat, 1), &
00135                                                grid_data(i)%nbr_points_per_lat, ierror)
00136                endif
00137 
00138                ! double or real data?
00139                if (associated (grid_data(i)%corners%st_array_dble)) then
00140 
00141                   ! the data provided to set_corners also contains all blocks the current
00142                   ! grid. I currently do not know how to extract only the required data.
00143                   ! This can be optimised later.
00144                   call psmile_set_corners_3d_double (psmile_id,                          &
00145                                                      grid_data(i)%nbr_corners,           &
00146                                                      grid_data(i)%corners%actual_shape,  &
00147                                                      grid_data(i)%corners%st_array_dble, &
00148                                                      grid_data(i)%corners%nd_array_dble, &
00149                                                      grid_data(i)%corners%rd_array_dble, &
00150                                                      ierror)
00151                else ! we are processing real data
00152 
00153                   call psmile_set_corners_3d_real (psmile_id,                          &
00154                                                    grid_data(i)%nbr_corners,           &
00155                                                    grid_data(i)%corners%actual_shape,  &
00156                                                    grid_data(i)%corners%st_array_real, &
00157                                                    grid_data(i)%corners%nd_array_real, &
00158                                                    grid_data(i)%corners%rd_array_real, &
00159                                                    ierror)
00160 
00161                endif ! (associated (grid_data(i)%corners%st_array_dble))
00162             endif ! (grid_type == PRISM_Gaussreduced_regvrt .or. grid_type == PRISM_Gaussreduced_sigmavrt)
00163          else ! grid_type is not gridless or gaussreduced
00164 
00165             size_of_shape = get_size_of_shape(grid_data(i)%grid_type)
00166 
00167             allocate (block_start_index(size_of_shape(2)), &
00168                       block_extent(size_of_shape(2)),      &
00169                       grid_extent(size_of_shape(2)),       &
00170                       block_valid_shape(size_of_shape(1),  &
00171                       size_of_shape(2)))
00172 
00173             block_start_index = grid_data(i)%grid_valid_shape(1,:)
00174 
00175             grid_extent = grid_data(i)%grid_valid_shape(2,:) - &
00176                           grid_data(i)%grid_valid_shape(1,:) + 1
00177 
00178             if (.not. partition_data_available (i)) then
00179                ! if the user did not provide any partition information, we use the
00180                ! grid_valid_shape to generate it
00181                call generate_partition_data(i)
00182             endif
00183 
00184             ! for all blocks in current grid
00185             do j = 1, size(grid_data(i)%partition_array, 1)
00186 
00187                block_valid_shape(1,:) = block_start_index
00188                block_valid_shape(2,:) = block_start_index + &
00189                                         grid_data(i)%extent_array(j, :) - 1
00190 
00191                block_extent = grid_data(i)%extent_array(j,:)
00192 
00193                call psmile_def_grid (psmile_id,              &
00194                                      grid_data(i)%grid_name, &
00195                                      grid_data(i)%comp_id,   &
00196                                      block_valid_shape,      &
00197                                      grid_data(i)%grid_type, &
00198                                      ierror)
00199 
00200                call add_pair(grid_id_map, i, psmile_id)
00201 
00202                !def partition (can be omitted, because there is only one block anyway)
00203                call psmile_def_partition(psmile_id, 1,                     &
00204                                        grid_data(i)%partition_array(j, :), &
00205                                        grid_data(i)%extent_array(j, :),    &
00206                                        ierror)
00207 
00208                ! double or real data?
00209                if (associated (grid_data(i)%corners%st_array_dble)) then
00210 
00211                   ! the data provided to set_corners also contains all blocks the current
00212                   ! grid. I currently do not know how to extract only the required data.
00213                   ! This can be optimised later.
00214                   call psmile_set_corners_3d_double (psmile_id,                          &
00215                                                      grid_data(i)%nbr_corners,           &
00216                                                      grid_data(i)%corners%actual_shape,  &
00217                                                      grid_data(i)%corners%st_array_dble, &
00218                                                      grid_data(i)%corners%nd_array_dble, &
00219                                                      grid_data(i)%corners%rd_array_dble, &
00220                                                      ierror)
00221                else ! we are processing real data
00222 
00223                   call psmile_set_corners_3d_real (psmile_id,                          &
00224                                                    grid_data(i)%nbr_corners,           &
00225                                                    grid_data(i)%corners%actual_shape,  &
00226                                                    grid_data(i)%corners%st_array_real, &
00227                                                    grid_data(i)%corners%nd_array_real, &
00228                                                    grid_data(i)%corners%rd_array_real, &
00229                                                    ierror)
00230                endif
00231 
00232                where (get_varying_dimension(block_extent, grid_extent, size_of_shape(2)))
00233 
00234                   block_start_index = block_start_index + block_extent
00235                end where
00236 
00237             enddo ! j = 1, size(grid_data(i)%partition_array, 1)
00238 
00239             deallocate (block_start_index, block_extent, block_valid_shape, grid_extent)
00240 
00241          endif ! grid_type
00242       enddo ! i = 1, size(grid_data)
00243 
00244       allocate (grid_ids(get_max_num_blocks_per_grid ()),  &
00245                 point_ids(get_max_num_blocks_per_grid ()), &
00246                 mask_ids(get_max_num_blocks_per_grid ()))
00247 
00248    endif ! associated (grid_data)
00249 
00250    if (associated (mask_data)) then
00251 
00252       ! for all mask
00253       do i = 1, size (mask_data)
00254 
00255          if (associated (mask_data(i)%mask)) then
00256 
00257             num_blocks = get_num_values (grid_id_map, mask_data(i)%grid_id)
00258             ! get the grid_ids of the associated grid
00259             grid_ids(1:num_blocks) = get_values(grid_id_map, mask_data(i)%grid_id, &
00260                                                 num_blocks)
00261 
00262             size_of_shape = get_size_of_shape (grid_data(mask_data(i)%grid_id)%grid_type)
00263 
00264             ! get the mask actual shape
00265             grid_actual_shape = 1
00266             grid_actual_shape(:,1:size_of_shape(2)) = mask_data(i)%mask_actual_shape
00267 
00268             ! for all blocks of the respective grid
00269             do j = 1, num_blocks
00270 
00271                ! get the actual shape of the current block
00272                block_actual_shape = 1
00273                block_actual_shape(:,1:size_of_shape(2)) =                                      &
00274                      get_block_actual_shape (mask_data(i)%mask_actual_shape,                   &
00275                                              Grids(grid_ids(j))%grid_shape,                    &
00276                                              grid_data(mask_data(i)%grid_id)%grid_valid_shape, &
00277                                              size_of_shape(2))
00278 
00279                call psmile_set_mask (psmile_id, grid_ids (j), block_actual_shape,                 &
00280                      get_sub_block_3d (mask_data(i)%mask, grid_actual_shape, block_actual_shape), &
00281                      .true., ierror)
00282 
00283                call add_pair (mask_id_map, i, psmile_id)
00284 
00285             enddo ! j
00286          endif ! (associated (mask_data(i)%mask))
00287       enddo ! i
00288    endif ! associated (mask_data)
00289 
00290    if (associated (point_data)) then
00291 
00292       ! for all point sets
00293       do i = 1, size (point_data)
00294 
00295          num_blocks = get_num_values (grid_id_map, point_data(i)%grid_id)
00296          ! get the grid_ids of the associated grid
00297          grid_ids(1:num_blocks) = get_values(grid_id_map, point_data(i)%grid_id, &
00298                                              num_blocks)
00299 
00300          ! for all blocks of the respective grid
00301          do j = 1, num_blocks
00302 
00303             ! for gridless grids we call psmile_set_points_gridless instead of psmile_set_points
00304             if (get_grid_type(point_data(i)%grid_id) == PRISM_Gridless) then
00305 
00306                call psmile_set_points_gridless(psmile_id, point_data(i)%name, &
00307                                              grid_ids(j), .true., ierror)
00308             else
00309 
00310                ! double or real data?
00311                if (associated (point_data(i)%points%st_array_dble)) then
00312 
00313                   call psmile_set_points_3d_double (psmile_id, point_data(i)%name, grid_ids(j), &
00314                                                     point_data(i)%points%actual_shape,          &
00315                                                     point_data(i)%points%st_array_dble,         &
00316                                                     point_data(i)%points%nd_array_dble,         &
00317                                                     point_data(i)%points%rd_array_dble, .true., ierror)
00318                else
00319 
00320                   call psmile_set_points_3d_real (psmile_id, point_data(i)%name, grid_ids(j), &
00321                                                   point_data(i)%points%actual_shape,          &
00322                                                   point_data(i)%points%st_array_real,         &
00323                                                   point_data(i)%points%nd_array_real,         &
00324                                                   point_data(i)%points%rd_array_real, .true., ierror)
00325                endif
00326             endif ! (get_grid_type(point_data(i)%grid_id) == PRISM_Gridless)
00327 
00328             call add_pair (point_id_map, i, psmile_id)
00329 
00330          enddo ! j = 1, num_blocks
00331       enddo ! i = 1, size (point_data)
00332    endif ! (associated (point_data))
00333 
00334    if (associated (var_data)) then
00335 
00336       ! for all vars
00337       do i = 1, size (var_data)
00338 
00339          grid_id = point_data(var_data(i)%point_id)%grid_id
00340          num_blocks = get_num_values (grid_id_map, grid_id)
00341 
00342          ! get the grid_ids of the associated grid
00343          grid_ids(1:num_blocks) = get_values(grid_id_map, grid_id, num_blocks)
00344          point_ids(1:num_blocks) = get_values(point_id_map, var_data(i)%point_id, &
00345                                               num_blocks)
00346 
00347          if (var_data(i)%mask_id /= prism_undefined) then
00348 
00349             mask_ids(1:num_blocks) = get_values(mask_id_map, var_data(i)%mask_id, &
00350                                                 num_blocks)
00351 
00352          else
00353             mask_ids(1:num_blocks) = prism_undefined
00354          endif
00355 
00356          ! for all blocks of the respective grid
00357          do j = 1, num_blocks
00358 
00359             call psmile_def_var (psmile_id, var_data(i)%var_name,              &
00360                                  grid_ids(j),                                  &
00361                                  point_ids(j),                                 &
00362                                  mask_ids(j),                                  &
00363                                  (/size (var_data(i)%var_actual_shape, 2),0/), &
00364                                  var_data(i)%var_actual_shape,                 &
00365                                  var_data(i)%var_type, ierror)
00366 
00367             call add_pair (var_id_map, i, psmile_id)
00368 
00369          enddo ! j = 1, num_blocks
00370       enddo ! i = 1, size (var_data)
00371    endif ! associated (var_data)
00372 
00373    if (allocated (grid_ids)) deallocate (grid_ids)
00374 
00375 #ifdef VERBOSE
00376    print 9980, trim(ch_id), ierror
00377    call psmile_flushstd
00378 #endif /* VERBOSE */
00379 
00380 9990 format (1x, a, ': psmile_apply_user_data: ')
00381 9980 format (1x, a, ': psmile_apply_user_data: eof ierror =', i5)
00382 
00383 contains
00384 
00385    ! ===========================================================================
00386 
00387    function get_varying_dimension(block_extent, grid_extent, ndim)
00388 
00389       integer, intent(in) :: ndim
00390       integer, intent(in) :: block_extent(ndim), grid_extent(ndim)
00391       logical             :: get_varying_dimension(ndim)
00392 
00393       get_varying_dimension = block_extent /= grid_extent
00394 
00395 #ifdef PRISM_ASSERTION
00396       if (count (get_varying_dimension) > 1) then
00397          print *, "block extent", block_extent
00398          print *, "grid extent", grid_extent
00399          call psmile_assert (__FILE__, __LINE__, &
00400                              "more than one varying dimension")
00401       endif
00402 #endif
00403    end function get_varying_dimension
00404 
00405    ! ===========================================================================
00406 
00407    function get_max_num_blocks_per_grid ()
00408 
00409       integer :: get_max_num_blocks_per_grid
00410       integer :: i
00411 
00412       get_max_num_blocks_per_grid = 1
00413 
00414       do i = 1, size (grid_data)
00415 
00416          get_max_num_blocks_per_grid = max (get_max_num_blocks_per_grid, &
00417                                             get_num_values(grid_id_map, i))
00418       enddo
00419 
00420    end function get_max_num_blocks_per_grid
00421 
00422    ! ===========================================================================
00423 
00424    function get_block_actual_shape (grid_actual_shape, block_valid_shape, grid_valid_shape, ndim)
00425 
00426       integer, intent (in) :: ndim
00427       integer, intent (in) :: grid_actual_shape (2, ndim), 
00428                               block_valid_shape (2, ndim), 
00429                               grid_valid_shape (2, ndim)
00430 
00431       integer :: get_block_actual_shape (2, ndim)
00432 
00433       where (block_valid_shape(2,:) == grid_valid_shape(2,:))
00434          get_block_actual_shape(2,:) = max (block_valid_shape(2,:), grid_actual_shape(2,:))
00435       elsewhere
00436          get_block_actual_shape(2,:) = block_valid_shape(2,:)
00437       end where
00438 
00439       where (block_valid_shape(1,:) == grid_valid_shape(1,:))
00440          get_block_actual_shape(1,:) = min (block_valid_shape(1,:), grid_actual_shape(1,:))
00441       elsewhere
00442          get_block_actual_shape(1,:) = block_valid_shape(1,:)
00443       end where
00444 
00445    end function
00446 
00447    ! ===========================================================================
00448 
00449    function get_sub_block_3d (src, src_shape, sub_block_shape)
00450 
00451       integer, intent (in) :: src_shape (2, ndim_3d), 
00452                               sub_block_shape (2, ndim_3d)
00453       logical, intent (in) :: src (src_shape(1,1):src_shape(2,1), 
00454                                    src_shape(1,2):src_shape(2,2), 
00455                                    src_shape(1,3):src_shape(2,3))
00456 
00457       logical :: get_sub_block_3d (sub_block_shape(1,1):sub_block_shape(2,1), 
00458                                    sub_block_shape(1,2):sub_block_shape(2,2), 
00459                                    sub_block_shape(1,3):sub_block_shape(2,3))
00460 
00461       get_sub_block_3d = src(sub_block_shape(1,1):sub_block_shape(2,1), &
00462                              sub_block_shape(1,2):sub_block_shape(2,2), &
00463                              sub_block_shape(1,3):sub_block_shape(2,3))
00464 
00465    end function get_sub_block_3d
00466 
00467 end subroutine psmile_apply_user_data