psmile_set_points_3d_double.F90

Go to the documentation of this file.
00001 !
00002 !-----------------------------------------------------------------------
00003 ! Copyright 2006-2010, CERFACS, Toulouse, France.
00004 ! Copyright 2006-2010, SGI Germany, Munich, Germany.
00005 ! Copyright 2006-2010, NEC Europe Ltd., London, UK.
00006 ! All rights reserved. Use is subject to OASIS4 license terms.
00007 !-----------------------------------------------------------------------
00008 !BOP
00009 !
00010 ! !ROUTINE: psmile_set_points_3d_double
00011 !
00012 ! !INTERFACE:
00013 
00014       subroutine psmile_set_points_3d_double ( method_id, point_name,  &
00015                 grid_id, points_actual_shape,                         &
00016                 points_1st_array, points_2nd_array, points_3rd_array, &
00017                 new_points, ierror )
00018 !
00019 ! !USES:
00020 !
00021       use PRISM
00022 !
00023       use PSMILe, dummy => psmile_set_points_3d_double
00024 !
00025       implicit none
00026 !
00027 ! !INPUT PARAMETERS:
00028 !
00029       Integer, Intent (In)                :: grid_id
00030 
00031 !     Specifies handle to the grid information created by routine
00032 !     PRISM_def_grid.
00033 
00034       Character(len=*), Intent (In)       :: point_name
00035 
00036 !     Specifies the name of the set of points specified, unique on each
00037 !     process within each component.
00038 
00039       Double Precision, Intent (In)       :: points_1st_array (*)
00040 
00041 !     Array defining the longitude of the grid points
00042 
00043       Double Precision, Intent (In)       :: points_2nd_array (*)
00044 
00045 !     Array defining the latitude of the grid points
00046 
00047       Double Precision, Intent (In)       :: points_3rd_array (*)
00048 !
00049       Integer, Intent (In)                :: points_actual_shape (1:2, *)
00050 
00051 !     Array defining the depth of the grid points
00052 
00053 !     Specifies the shape for the "n_dim"-dimensional
00054 !     (including halo/overlap regions) with
00055 
00056 !     points_actual_shape (1,     1) = Lowest  first   dimension of points_1st_array, ...
00057 !     points_actual_shape (2,     1) = Highest first   dimension of points_1st_array, ...
00058 !     points_actual_shape (1,     2) = Lowest  second  dimension of points_1st_array, ...
00059 !     points_actual_shape (2,     2) = Highest second  dimension of points_1st_array, ...
00060 !        ...
00061 !     points_actual_shape (1,     i) = Lowest  i-th    dimension of points_1st_array, ...
00062 !     points_actual_shape (2,     i) = Highest i-th    dimension of points_1st_array, ...
00063 !     points_actual_shape (1, n_dim) = Lowest  n_dim-th dimension of points_1st_array, ...
00064 !     points_actual_shape (2, n_dim) = Highest n_dim-th dimension of points_1st_array, ...
00065 
00066 !     where n_dim depends on the grid_type. n_dim is
00067 !     set to an appropriate value in PRISM_def_grid.
00068 
00069       Logical, Intent (In)                :: new_points
00070 
00071 !     Logical to indicate whether new points are specified (.true.) or old points
00072 !     shall be overwritten
00073 !
00074 ! !INPUT/OUTPUT PARAMETERS:
00075 !
00076       Integer, intent (InOut)             :: method_id
00077 
00078 !     handle to the defined coordinate points
00079 !
00080 ! !OUTPUT PARAMETERS:
00081 !
00082       Integer, intent (Out)               :: ierror
00083 
00084 !     Returns the error code of psmile_set_points_3d_double;
00085 !             ierror = 0 : No error
00086 !             ierror > 0 : Severe error
00087 !
00088 ! !LOCAL VARIABLES
00089 !
00090       Type (Coords_Block), pointer :: coords_pointer
00091 
00092       Integer(kind=int64)          :: length, len_alloc(ndim_3d)
00093       Integer                      :: i, n_dim
00094       Integer                      :: overlap (2, ndim_3d)
00095       Logical                      :: overlap_provided
00096 !
00097       Integer, parameter           :: nerrp = 3
00098       Integer                      :: ierrp (nerrp)
00099 !
00100 !
00101 ! !DESCRIPTION:
00102 !
00103 ! Subroutine "psmile_set_points_3d_double" defines the localisation of the
00104 ! center 3 coordinate arrays of a grid covering a 3D physical space
00105 ! for the grid which has grid id "grid_id".
00106 !
00107 ! !REVISION HISTORY:
00108 !
00109 !   Date      Programmer   Description
00110 ! ----------  ----------   -----------
00111 ! 01.12.03    H. Ritzdorf  created
00112 !
00113 !EOP
00114 !----------------------------------------------------------------------
00115 !
00116 ! $Id: psmile_set_points_3d_double.F90 2773 2010-11-25 14:48:32Z hanke $
00117 ! $Author: hanke $
00118 !
00119   Character(len=len_cvs_string), save :: mycvs = 
00120       '$Id: psmile_set_points_3d_double.F90 2773 2010-11-25 14:48:32Z hanke $'
00121 !
00122 !----------------------------------------------------------------------
00123 
00124 #ifdef VERBOSE
00125       print *, trim(ch_id), ': psmile_set_points_3d_double: grid_id', grid_id
00126 
00127       call psmile_flushstd
00128 #endif /* VERBOSE */
00129 
00130 !-----------------------------------------------------------------------
00131 !  Initialization
00132 !-----------------------------------------------------------------------
00133 
00134       ierror = 0
00135       length = 1       ! Total number of elements to be allocated
00136       len_alloc(:) = 1 ! Number of elements per direction
00137 !
00138       overlap(:, :) = 0 ! Extent of additional overlap
00139       overlap_provided = .false. ! Has application overlap area provided
00140 
00141       Nullify ( coords_pointer )
00142 
00143 !-----------------------------------------------------------------------
00144 !  Control input arguments
00145 !-----------------------------------------------------------------------
00146 
00147       if (grid_id < 1 .or. &
00148           grid_id > Number_of_Grids_allocated ) then
00149          ierrp (1) = grid_id
00150          ierrp (2) = Number_of_Grids_allocated
00151 
00152          ierror = PRISM_Error_Arg
00153 
00154          call psmile_error ( ierror, 'grid_id', &
00155                              ierrp, 2, __FILE__, __LINE__ )
00156          return
00157       endif
00158 !
00159       if (Grids(grid_id)%status /= PSMILe_status_defined) then
00160          ierrp (1) = grid_id
00161 
00162          ierror = PRISM_Error_Arg
00163 
00164          call psmile_error ( PRISM_Error_Arg, 'grid_id (not active)', &
00165                              ierrp, 1, __FILE__, __LINE__ )
00166          return
00167       endif
00168 
00169       if ( Grids(grid_id)%grid_type == PRISM_Gridless ) then
00170          print *, trim(ch_id), ': psmile_set_points_3d_double: No need to '
00171          print *, trim(ch_id), ':     call psmile_set_points for PRISM_Gridless.'
00172          return
00173       endif
00174 
00175 !-----------------------------------------------------------------------
00176 !  Store data on center coordinates
00177 !-----------------------------------------------------------------------
00178 
00179       if (Grids(grid_id)%grid_structure == PSMILe_Grid_Block .or. &
00180           Grids(grid_id)%grid_structure == PSMILe_Grid_Unstruct) then
00181 
00182          if ( new_points ) then
00183 
00184 !-----------------------------------------------------------------------
00185 !  get a new method Id
00186 !-----------------------------------------------------------------------
00187 
00188             call psmile_get_method_handle (grid_id, method_id, ierror)
00189 
00190             if (ierror > 0) then
00191                ierrp (1) = method_id
00192                ierror = PRISM_Error_Arg
00193                call psmile_error ( PRISM_Error_Arg, 'Failed to get a new method id', &
00194                                    ierrp, 1, __FILE__, __LINE__ )
00195                return
00196             endif
00197 
00198             ! Rene: Do we really have to allocate this pointer as it is not an array?
00199             ! HR: Yes, we don't know the grid structure type
00200 
00201             Allocate(Methods(method_id)%coords_pointer, STAT = ierror)
00202             if (ierror > 0) then
00203                ierrp (1) = ierror
00204                ierror = PRISM_Error_Alloc
00205                call psmile_error ( ierror, 'coords_pointer', &
00206                     ierrp, 1, __FILE__, __LINE__ )
00207                return
00208             endif
00209 
00210             coords_pointer => Methods(method_id)%coords_pointer
00211 
00212             coords_pointer%coords_datatype  = MPI_DATATYPE_NULL
00213 
00214                do i = 1, ndim_3d
00215                Nullify ( coords_pointer%coords_real(i)%vector )
00216                Nullify ( coords_pointer%coords_dble(i)%vector )
00217 #if defined ( PRISM_QUAD_TYPE )
00218                Nullify ( coords_pointer%coords_quad(i)%vector )
00219 #endif
00220                end do
00221 
00222 ! HR: Shouldn't this contained in Type Coords_Block;
00223 !     Otherwise move nullify into PSMILe_Get_method_handle
00224 
00225             Nullify ( Methods(method_id)%halo_pointer )
00226 
00227          else
00228 
00229 !           This method does already exist
00230 
00231             coords_pointer => Methods(method_id)%coords_pointer
00232 
00233             if ( Methods(method_id)%method_type /= PSMILe_PointMethod ) then
00234 
00235                ierror = PRISM_Error_Parameter
00236                ierrp (1) = ierror
00237                ierrp (2) = method_id
00238                call psmile_error ( ierror, 'not a method_id for points', &
00239                                    ierrp, 2, __FILE__, __LINE__ )
00240                return
00241             endif
00242 
00243             if ( coords_pointer%coords_datatype /= MPI_DATATYPE_NULL) then
00244 
00245                 print *, 'Warning: going to replace coordinates for method Id', method_id
00246 
00247             endif
00248 
00249          endif ! if (new_points)
00250 
00251 !-----------------------------------------------------------------------
00252 ! Get n_dim which was defined in PRISM_def_grid
00253 !-----------------------------------------------------------------------
00254 
00255          n_dim = Grids(grid_id)%n_dim
00256 
00257          if ( n_dim > ndim_3d ) then
00258             ierrp (1) = n_dim
00259 
00260             ierror = PRISM_Error_Internal
00261             call psmile_error ( ierror, 'unsupported dimension Grids(grid_id)%n_dim', &
00262                                 ierrp, 1, __FILE__, __LINE__ )
00263             return
00264           endif
00265 
00266 !-----------------------------------------------------------------------
00267 !  Free old arrays if already allocated
00268 !-----------------------------------------------------------------------
00269 
00270       do i = 1, ndim_3d
00271          if (Associated(coords_pointer%coords_dble(i)%vector)) then
00272 
00273             Deallocate (coords_pointer%coords_dble(i)%vector, STAT = ierror)
00274 
00275             if (ierror > 0) then
00276                ierrp (1) = ierror
00277                ierror = PRISM_Error_Dealloc
00278                call psmile_error ( ierror, 'coords_dble(i)%vector', &
00279                     ierrp, 1, __FILE__, __LINE__ )
00280                return
00281             endif
00282 
00283          endif
00284       end do
00285 
00286 !-----------------------------------------------------------------------
00287 !  Store shapes
00288 !-----------------------------------------------------------------------
00289 
00290       select case ( Grids(grid_id)%grid_type )
00291 
00292       case ( PRISM_Gaussreduced_regvrt )
00293 
00294          coords_pointer%coords_shape(1:2,1) = points_actual_shape (1:2,1)
00295          coords_pointer%coords_shape(1:2,2) = 1
00296          coords_pointer%coords_shape(1:2,3) = points_actual_shape (1:2,2)
00297 
00298          coords_pointer%actual_shape = coords_pointer%coords_shape
00299 
00300       case ( PRISM_Unstructlonlat_regvrt )
00301 
00302          coords_pointer%coords_shape(1:2,1) = points_actual_shape (1:2,1)
00303          coords_pointer%coords_shape(1:2,2) = 1
00304          coords_pointer%coords_shape(1:2,3) = points_actual_shape (1:2,2)
00305 
00306          coords_pointer%actual_shape = coords_pointer%coords_shape
00307 
00308       case ( PRISM_Irrlonlatvrt, PRISM_Reglonlatvrt )
00309 
00310          coords_pointer%coords_shape(1:2,1:n_dim) = &
00311                 points_actual_shape (1:2,1:n_dim)
00312 
00313          coords_pointer%actual_shape = coords_pointer%coords_shape
00314 
00315       case ( PRISM_Irrlonlat_regvrt )
00316 !        Extend (lon, lat) region by overlap region if necessary
00317 !        This extended region is used to store the boundaries of
00318 !        (virtual) method cell
00319 
00320          overlap (1, 1:ndim_2d) = min (0,                                  &
00321                              (Grids(grid_id)%grid_shape(1, 1:ndim_2d)-1) - &
00322                                     points_actual_shape(1, 1:ndim_2d)) 
00323          overlap (2, 1:ndim_2d) = max (0,                                  &
00324                              (Grids(grid_id)%grid_shape(2, 1:ndim_2d)+1) - &
00325                                     points_actual_shape(2, 1:ndim_2d))
00326 
00327          coords_pointer%coords_shape (1:2,1:n_dim) = &
00328                  points_actual_shape (1:2,1:n_dim) + overlap (1:2, 1:n_dim)
00329 
00330          overlap_provided = minval (overlap (1, 1:ndim_2d)) == 0 .and. &
00331                             maxval (overlap (2, 1:ndim_2d)) == 0
00332 
00333          coords_pointer%actual_shape (1:2,1:n_dim) =  points_actual_shape(1:2,1:n_dim)
00334 
00335       case DEFAULT
00336 
00337          ierrp (1) = Grids(grid_id)%grid_type
00338 
00339          ierror = PRISM_Error_Internal
00340 
00341          call psmile_error ( ierror, 'unsupported grid generation type', &
00342                               ierrp, 1, __FILE__, __LINE__ )
00343 
00344          return
00345 
00346       end select
00347 
00348 !-----------------------------------------------------------------------
00349 ! Determine length of arrays
00350 !-----------------------------------------------------------------------
00351 
00352          do i = 1, n_dim
00353             length = length * ( coords_pointer%coords_shape(2,i) - &
00354                                 coords_pointer%coords_shape(1,i) + 1 )
00355          enddo
00356 
00357       Methods(method_id)%size = length
00358 
00359 !-----------------------------------------------------------------------
00360 !  Determine lengths of center coordinates
00361 !-----------------------------------------------------------------------
00362 
00363       select case ( Grids(grid_id)%grid_type )
00364 
00365       case ( PRISM_Unstructlonlatvrt )
00366 
00367          len_alloc (:) = length
00368 
00369       case ( PRISM_Unstructlonlat_regvrt )
00370 
00371          len_alloc (3)   = coords_pointer%coords_shape(2,2) - &
00372                            coords_pointer%coords_shape(1,2) + 1
00373          len_alloc (1:2) = length / len_alloc (3)
00374 
00375       case ( PRISM_Unstructlonlat_sigmavrt )
00376 
00377          len_alloc (3)   = coords_pointer%coords_shape(2,2) - &
00378                            coords_pointer%coords_shape(1,2) + 1
00379          len_alloc (1:2) = length
00380 
00381       case ( PRISM_Irrlonlatvrt )
00382 
00383          len_alloc (:) = length
00384 
00385       case ( PRISM_Irrlonlat_sigmavrt )
00386 
00387          len_alloc (3) = coords_pointer%coords_shape(2,3) - &
00388                          coords_pointer%coords_shape(1,3) + 1
00389          len_alloc (1:2) = length
00390 
00391       case ( PRISM_Irrlonlat_regvrt )
00392 
00393          len_alloc (3) = coords_pointer%coords_shape(2,3) - &
00394                          coords_pointer%coords_shape(1,3) + 1
00395          len_alloc (1:2) = length / len_alloc (3)
00396 
00397       case ( PRISM_Reglonlatvrt )
00398 
00399          do i = 1, 3
00400          len_alloc (i) = coords_pointer%coords_shape(2,i) - &
00401                          coords_pointer%coords_shape(1,i) + 1
00402          enddo
00403 
00404       case ( PRISM_Gaussreduced_regvrt )
00405 
00406          len_alloc (3)   = coords_pointer%coords_shape(2,3) - &
00407                            coords_pointer%coords_shape(1,3) + 1
00408          len_alloc (1:2) = length / len_alloc (3)
00409 
00410       case DEFAULT
00411 
00412          ierror = PRISM_Error_Internal
00413 
00414          ierrp (1) = method_id
00415          ierrp (2) = grid_id
00416          ierrp (3) = Grids(grid_id)%grid_type
00417 
00418          call psmile_error ( ierror, 'unknown grid generation type', &
00419                              ierrp, 3, __FILE__, __LINE__ )
00420          return
00421 
00422       end select
00423 
00424 !-----------------------------------------------------------------------
00425 ! Check range of points_actual_shape.
00426 !-----------------------------------------------------------------------
00427 
00428          do i = 1, n_dim
00429             if (coords_pointer%coords_shape(1,i) > Grids(grid_id)%grid_shape(1,i) .or. &
00430                 coords_pointer%coords_shape(2,i) < Grids(grid_id)%grid_shape(2,i)) exit
00431          enddo
00432 
00433          if ( i <= n_dim ) then
00434 
00435              ierror = PRISM_Error_Arglist
00436 
00437              call psmile_error ( PRISM_Error_Arglist, &
00438                                  'points_actual_shape too small', &
00439                                  points_actual_shape, n_dim*2, &
00440                                  __FILE__, __LINE__ )
00441              return
00442          endif
00443 
00444       else
00445 
00446 !-----------------------------------------------------------------------
00447 !  ... 3-dimensional center coordinates for an invalid method type
00448 !-----------------------------------------------------------------------
00449 
00450          ierrp (1) = method_id
00451          ierror = PRISM_Error_Arg
00452 
00453          call psmile_error ( PRISM_Error_Arg, &
00454                              'grid structure not supported.', &
00455                              ierrp, 1, __FILE__, __LINE__ )
00456          return
00457 
00458       endif
00459 
00460 !-----------------------------------------------------------------------
00461 !  Allocate vectors
00462 !-----------------------------------------------------------------------
00463 
00464       do i = 1, ndim_3d
00465          Allocate (coords_pointer%coords_dble(i)%vector(1:len_alloc(i)), &
00466                    STAT = ierror)
00467 
00468          if ( ierror > 0 ) then
00469             ierrp (1) = ierror
00470             ierrp (2) = len_alloc (i)
00471 
00472             ierror = PRISM_Error_Alloc
00473             call psmile_error ( ierror, 'coords_dble(i)%vector', &
00474                  ierrp, 2, __FILE__, __LINE__ )
00475             return
00476          endif
00477       end do
00478 
00479 !-----------------------------------------------------------------------
00480 !     Copy data to internal Coordinates array
00481 !-----------------------------------------------------------------------
00482 !
00483       if ( Grids(grid_id)%grid_type == PRISM_Irrlonlatvrt .or. &
00484            Grids(grid_id)%grid_type == PRISM_Unstructlonlatvrt ) then
00485 !
00486 !   Irregular grid type: 3 full 3d Arrays
00487 !
00488          coords_pointer%coords_dble(1)%vector = points_1st_array (1:len_alloc(1))
00489          coords_pointer%coords_dble(2)%vector = points_2nd_array (1:len_alloc(1))
00490          coords_pointer%coords_dble(3)%vector = points_3rd_array (1:len_alloc(1))
00491 
00492       else if ( Grids(grid_id)%grid_type == PRISM_Unstructlonlat_sigmavrt .or. &
00493                 Grids(grid_id)%grid_type == PRISM_Unstructlonlat_regvrt ) then
00494 !
00495 !   Unstructured grid types with regular vertical axis:
00496 !   2 full 3d Arrays, 1d Array
00497 !
00498          coords_pointer%coords_dble(1)%vector = points_1st_array (1:len_alloc(1))
00499          coords_pointer%coords_dble(2)%vector = points_2nd_array (1:len_alloc(1))
00500          coords_pointer%coords_dble(3)%vector = points_3rd_array (1:len_alloc(3))
00501 
00502       else if ( Grids(grid_id)%grid_type == PRISM_Irrlonlat_regvrt .or. &
00503                 Grids(grid_id)%grid_type == PRISM_Irrlonlat_sigmavrt ) then
00504 !
00505 !   Regular grid type: Regular in Z-direction, 2 full 2d Arrays
00506 !
00507          if (overlap_provided) then
00508             coords_pointer%coords_dble(1)%vector = &
00509                points_1st_array (1:len_alloc(1))
00510             coords_pointer%coords_dble(2)%vector = &
00511                points_2nd_array (1:len_alloc(2))
00512          else
00513             call psmile_reshape_2d_dble (points_1st_array,       &
00514                points_actual_shape (1:2, 1:ndim_2d),             &
00515                points_actual_shape (1:2, 1:ndim_2d),             &
00516                coords_pointer%coords_dble(1)%vector,             &
00517                coords_pointer%coords_shape(1:2,1:ndim_2d), ierror)
00518 
00519             call psmile_reshape_2d_dble (points_2nd_array,       &
00520                points_actual_shape (1:2, 1:ndim_2d),             &
00521                points_actual_shape (1:2, 1:ndim_2d),             &
00522                coords_pointer%coords_dble(2)%vector,             &
00523                coords_pointer%coords_shape(1:2,1:ndim_2d), ierror)
00524          endif
00525 
00526          coords_pointer%coords_dble(3)%vector = points_3rd_array (1:len_alloc(3))
00527 
00528       else if ( Grids(grid_id)%grid_type == PRISM_Reglonlatvrt ) then
00529 !
00530 !   Regular grid type: Regular in all directions, 3 1d Arrays
00531 !
00532          coords_pointer%coords_dble(1)%vector = points_1st_array (1:len_alloc(1))
00533          coords_pointer%coords_dble(2)%vector = points_2nd_array (1:len_alloc(2))
00534          coords_pointer%coords_dble(3)%vector = points_3rd_array (1:len_alloc(3))
00535 
00536       else if ( Grids(grid_id)%grid_type == PRISM_Gaussreduced_sigmavrt .or. &
00537                 Grids(grid_id)%grid_type == PRISM_Gaussreduced_regvrt ) then
00538 !
00539 !   Gaussian reduced grid types with regular vertical axis:
00540 !   2 2d Arrays, 1d Array
00541 !
00542          coords_pointer%coords_dble(1)%vector = points_1st_array (1:len_alloc(1))
00543          coords_pointer%coords_dble(2)%vector = points_2nd_array (1:len_alloc(1))
00544          coords_pointer%coords_dble(3)%vector = points_3rd_array (1:len_alloc(3))
00545 
00546       else
00547 
00548          ierrp (1) = Grids(grid_id)%grid_type
00549 
00550          ierror = PRISM_Error_Internal
00551 
00552          call psmile_error ( ierror, 'unsupported grid generation type', &
00553                              ierrp, 1, __FILE__, __LINE__ )
00554          return
00555 
00556       endif
00557 !
00558 !-----------------------------------------------------------------------
00559 !     Update status and other info's
00560 !-----------------------------------------------------------------------
00561 !
00562 !     Mark status as changed
00563 !     ? Special status for ``unimportant data'' required?
00564 !
00565       Methods(method_id)%status          = PSMILe_status_defined
00566 !
00567       Methods(method_id)%used_for_coupling = .false.
00568 !
00569       Methods(method_id)%grid_id         = grid_id
00570 !
00571       Methods(method_id)%method_type     = PSMILe_PointMethod
00572 !
00573       Methods(method_id)%point_name      = point_name
00574 !
00575       coords_pointer%coords_datatype     = MPI_DOUBLE_PRECISION
00576 !     coords_pointer%overlap_provided    = overlap_provided
00577 
00578 #ifdef VERBOSE
00579       print *, trim(ch_id), ': psmile_set_points_3d_double: eof ierror =', ierror
00580 
00581       call psmile_flushstd
00582 #endif /* VERBOSE */
00583 
00584       end subroutine psmile_set_points_3d_double