psmile_bbcells_gauss2_dble.F90

Go to the documentation of this file.
00001 !-----------------------------------------------------------------------
00002 ! Copyright 2006-2010, NEC Europe Ltd., London, UK.
00003 ! All rights reserved. Use is subject to OASIS4 license terms.
00004 !-----------------------------------------------------------------------
00005 !BOP
00006 !
00007 ! !ROUTINE: PSMILe_Bbcells_gauss2_dble
00008 !
00009 ! !INTERFACE:
00010 
00011       subroutine psmile_bbcells_gauss2_dble (array_x, array_y,     &
00012                            shape, range, nbr_lats, points_per_lat, &
00013                            corners_y, corner_shape, nbr_corners,   &
00014                            chmin_x, chmax_x, midp_x, &
00015                            chmin_y, chmax_y, midp_y, &
00016                            nbrs, levdim, ierror )
00017 !
00018 ! !USES:
00019 !
00020       use PRISM_constants
00021 !
00022       use PSMILe, dummy_interface => PSMILe_Bbcells_gauss2_dble
00023 
00024       implicit none
00025 !
00026 ! !INPUT PARAMETERS:
00027 !
00028       Integer, Intent (In)            :: shape (2)
00029 
00030 !     Dimension of (method/grid) coordinate "array_*"
00031 
00032       Double Precision, Intent (In)   :: array_x(shape(1):shape(2))
00033 
00034 !     Coordinates of the longitude method
00035 
00036       Double Precision, Intent (In)   :: array_y(shape(1):shape(2))
00037 
00038 !     Coordinates of the latitude method
00039 
00040       Integer, Intent (In)            :: range (2)
00041 
00042 !     Definintion of the subgrid for which the (bounding) cells should
00043 !     be computed.
00044 
00045       Integer, Intent (In)            :: corner_shape (2)
00046 
00047 !     Dimension of the array points_per_lat 
00048 
00049       Integer,           Intent (In)  :: nbr_lats
00050 
00051 !     Container for the number of points for each latitude band
00052 
00053       Integer,           Intent (In)  :: points_per_lat(nbr_lats,2)
00054 
00055 !     Dimension of the array "corners" containing the corner points
00056 !     of the control volume.
00057 
00058       Integer,           Intent (In)  :: nbr_corners
00059 
00060 !     Number of corners per control volume
00061 !     Must be 2.
00062 
00063       Double Precision,  Intent (In)  :: corners_y (corner_shape(1): 
00064                                                     corner_shape(2), 
00065                                                     nbr_corners)
00066 
00067 !     Coordinates of corner points of control volume
00068 
00069       Integer,           Intent (In)  :: levdim(1)
00070 
00071 !     Dimension of chmin, chmax (nur zur kontrolle)
00072 
00073 ! !OUTPUT PARAMETERS:
00074 !
00075       Double Precision,  Intent (Out) :: chmin_x (range(1)-1:range(2))
00076 !
00077 !     Minimum value of bounding box of coordinate.
00078 !
00079       Double Precision,  Intent (Out) :: chmax_x (range(1)-1:range(2))
00080 !
00081 !     Maximum value of bounding box of coordinate.
00082 !
00083       Double Precision,  Intent (Out) :: midp_x  (range(1)-1:range(2))
00084 !
00085 !     Maxi of bounding box of coordinate.
00086       Double Precision,  Intent (Out) :: chmin_y (range(1)-1:range(2))
00087 !
00088 !     Minimum value of bounding box of coordinate.
00089 !
00090       Double Precision,  Intent (Out) :: chmax_y (range(1)-1:range(2))
00091 !
00092 !     Maximum value of bounding box of coordinate.
00093 !
00094       Double Precision,  Intent (Out) :: midp_y  (range(1)-1:range(2))
00095 !
00096 !     Mid point of bounding box of coordinate.
00097 !
00098       Integer,  Intent (Out)          :: nbrs    (range(2)-range(1)+1,4)
00099 !
00100 !     List of neighbour points
00101 !
00102       Integer,           Intent (Out) :: ierror
00103 
00104 !     Returns the error code of PSMILe_Bbcells_gauss2_dble
00105 !             ierror = 0 : No error
00106 !             ierror > 0 : Severe error
00107 !
00108 ! !LOCAL VARIABLES
00109 !
00110       Integer                         :: i, ii, ibeg, iend, jlat
00111       Integer                         :: nbr, isubs, isube, ilow, ihigh, last
00112 !
00113       Double Precision                :: dist, dist_low, dist_high
00114       Double Precision                :: chmin_cy
00115       Double Precision                :: chmax_cx, chmax_cy
00116 !
00117 !
00118 ! !DESCRIPTION:
00119 !
00120 ! Subroutine "PSMILe_Bbcells_gauss2_dble" computes the bounding box for each
00121 ! horizontal cell of an gaussian reduced grid. Probably the worst code I've ever
00122 ! written with lots of options for improvements. rr
00123 !
00124 ! chmin (i) returns the minimal value of bounding box of
00125 !           method cell (i:i+1)
00126 ! chmax (i) returns the maximal value of bounding box of
00127 !           method cell (i:i+1)
00128 !
00129 ! The first box
00130 ! chmin (range(1)-1) returns the minimal value of the bounding box from
00131 !                    the virtual cell (array(range(1)):bounding box of
00132 !                                                      corner control volume)
00133 ! The last box
00134 ! chmin (range(2))   returns the minimal value of the bounding box from
00135 !                    the virtual cell (array(range(2)):bounding box of
00136 !                                                      corner control volume)
00137 !
00138 ! !REVISION HISTORY:
00139 !
00140 !   Date      Programmer   Description
00141 ! ----------  ----------   -----------
00142 ! 05.07.10    R. Redler    created
00143 !
00144 !EOP
00145 !----------------------------------------------------------------------
00146 !
00147 !  $Id: psmile_bbcells_gauss2_dble.F90 2325 2010-04-21 15:00:07Z valcke $
00148 !  $Author: valcke $
00149 !
00150    Character(len=len_cvs_string), save :: mycvs = 
00151        '$Id: psmile_bbcells_gauss2_dble.F90 2325 2010-04-21 15:00:07Z valcke $'
00152 !
00153 !----------------------------------------------------------------------
00154 !
00155 #ifdef VERBOSE
00156       print 9990, trim(ch_id)
00157 
00158       call psmile_flushstd
00159 #endif /* VERBOSE */
00160 !
00161 #ifdef PRISM_ASSERTION
00162       if (range(1) < corner_shape(1) .or. range(2) > corner_shape(2)) then
00163          print *, trim(ch_id), "PSMILe_bbcells_gauss2_dble: range, corner", &
00164                   range, corner_shape
00165          call psmile_assert (__FILE__, __LINE__, &
00166               "levdim_corner must be >= range(2)-range(1)")
00167       endif
00168 !
00169       if (levdim(1) < range(2)-range(1)+2) then
00170          print *, trim(ch_id), "PSMILe_bbcells_gauss2_dble: levdim", levdim
00171          call psmile_assert (__FILE__, __LINE__, &
00172               "levdim(1) < range(2)-range(1)+2")
00173       endif
00174 !
00175       if (nbr_corners /= 2) then
00176          print *, trim(ch_id), "PSMILe_bbcells_gauss2_dble: nbr_corner", &
00177                   nbr_corners
00178          call psmile_assert (__FILE__, __LINE__, &
00179               "nbr_corners != 2")
00180       endif
00181 #endif
00182 !
00183 !  Initialization
00184 !
00185       ierror = 0
00186       nbrs   = 0
00187 !
00188       jlat = 1
00189       ibeg = 1
00190       iend = points_per_lat(jlat,1)
00191 !
00192 !  Get bounding box of the first and last corner cell
00193 !
00194 !  For southern- and northernmost lines take corners
00195 !
00196       chmin_cy = min(corners_y(range(1),1),corners_y(range(1),2))
00197       chmax_cy = max(corners_y(range(2),1),corners_y(range(2),2))
00198 !
00199 !  For process boundaries values from the neighbor process have to be taken
00200 !
00201 !  ... to be implemented ...
00202 !
00203 !  Latitudes of bounding boxes 
00204 !
00205       chmin_x(range(1)-1) = PSMILE_DUNDEF
00206       chmin_y(range(1)-1) = chmin_cy
00207       chmax_x(range(1)-1) = PSMILE_DUNDEF
00208       chmax_y(range(1)-1) = array_y(ibeg)
00209       midp_x(range(1)-1)  = PSMILE_DUNDEF
00210       midp_y(range(1)-1)  = PSMILE_DUNDEF
00211 
00212       do i = range(1), range(2)
00213 
00214          if ( jlat < nbr_lats ) then
00215             chmin_y(i) = array_y(ibeg)
00216             chmax_y(i) = array_y(iend+1)
00217          else
00218             chmin_y(i) = array_y(iend)
00219             chmax_y(i) = chmax_cy
00220          endif
00221 
00222          if ( i == iend .and. i < range(2) ) then
00223             jlat = jlat + 1
00224             ibeg = iend + 1
00225             iend = iend + points_per_lat(jlat,1)
00226          endif
00227 
00228       enddo
00229 
00230       jlat = 1
00231       ibeg = 1
00232       iend = points_per_lat(jlat,1)
00233 
00234       do i = range(1), range(2)
00235 !
00236 !        For cyclic boundaries take values from the other end of the x-axis
00237 !
00238          chmax_cx = array_x(ibeg) + dble(360.0)
00239 !
00240 !  For process boundaries values from the neighbor process have to be taken
00241 !
00242 !  ... to be implemented ...
00243 !
00244          chmin_x(i) = array_x(i)
00245          nbrs(i,1) = i
00246          nbrs(i,4) = i
00247 
00248          if ( i < iend ) then
00249             chmax_x(i) = array_x(i+1)
00250             nbrs(i,2) = i+1
00251             nbrs(i,3) = i+1
00252          else
00253             chmax_x(i) = chmax_cx
00254             nbrs(i,2) = -ibeg
00255             nbrs(i,3) = -ibeg
00256          endif
00257          
00258          if ( jlat < nbr_lats ) then
00259             isubs = ibeg       !needs to be optimised
00260             isube = iend       !needs to be optimised
00261          else
00262             isubs = max(i  ,ibeg)
00263             isube = min(i+1,iend)
00264          endif
00265 
00266          nbr = isube-isubs+1
00267 
00268          ilow = isubs 
00269 
00270          if ( nbr == 1 ) then
00271 
00272             dist_low  = abs(array_x(isubs)-array_x(i))
00273             dist_high = abs(array_x(isube)-array_x(i))
00274             ihigh     = isubs 
00275 
00276          else
00277 
00278             dist_low  = huge(chmax_cx)
00279             dist_high = huge(chmax_cx)
00280             ihigh     = isube
00281 
00282          endif
00283 
00284          do ii = isubs, isube
00285 
00286             dist = abs(array_x(i)-array_x(ii))
00287 
00288             if ( dist < dist_low ) then
00289                dist_low = dist
00290                ilow     = ii
00291             endif
00292 
00293          enddo
00294 
00295          do ii = isubs, isube
00296 
00297             if ( i < iend ) then
00298                dist = abs(array_x(i+1)-array_x(ii))
00299             else
00300                dist = abs(array_x(i  )-array_x(ii))
00301             endif
00302 
00303             if ( dist < dist_high ) then
00304                ihigh     = ii
00305                dist_high = dist
00306             endif
00307 
00308          enddo
00309 
00310          if ( ilow == ihigh ) then
00311 
00312             if ( dist_low <= dist_high ) then
00313                ihigh = ihigh+1
00314             else
00315                ilow  = ilow -1
00316             endif
00317 
00318             if ( ilow < ibeg ) then
00319                ilow  = ibeg
00320                ihigh = ibeg + 1
00321             endif
00322 
00323             if ( ihigh > iend ) then
00324                ihigh = -iend
00325                ilow  = ilow-1
00326             endif
00327 
00328          endif
00329 
00330          nbrs(i,3) = ihigh
00331          nbrs(i,4) = ilow
00332 
00333          if ( i == iend .and. i < range(2) ) then
00334             jlat = jlat + 1
00335             ibeg = iend + 1
00336             iend = iend + points_per_lat(jlat,1)
00337          endif
00338 
00339       enddo ! i-loop
00340 
00341       jlat = 1
00342       ibeg = 1
00343       iend = points_per_lat(jlat,1)
00344 
00345       do i = range(1), range(2)
00346          if ( nbrs(i,2) < 0 .or. nbrs(i,3) < 0 ) then
00347             chmax_x(i) = array_x(ibeg) + dble(360.0)
00348             chmin_x(i) = array_x(nbrs(i,1))
00349          else
00350             chmin_x(i) = min(array_x(nbrs(i,1)),array_x(nbrs(i,4)))
00351             chmax_x(i) = max(array_x(nbrs(i,2)),array_x(nbrs(i,3)))
00352         endif
00353 
00354          if ( i == iend .and. i < range(2) ) then
00355             jlat = jlat + 1
00356             ibeg = iend + 1
00357             iend = iend + points_per_lat(jlat,1)
00358          endif
00359 
00360       enddo
00361 
00362 
00363 !
00364 !  Now determine a approximation of the mid point of the cell
00365 !
00366       do i = range(1), range(2)
00367          midp_y(i) = dble(0.5) * ( chmax_y(i) + chmin_y(i) )
00368       enddo
00369 
00370       last = sum(points_per_lat(1:nbr_lats-1,1))
00371 
00372       jlat = 1
00373       ibeg = 1
00374       iend = points_per_lat(jlat,1)
00375 
00376       do i = range(1), range(2)
00377 
00378          if ( i < last ) then
00379 
00380             if ( nbrs(i,2) < 0 .and. nbrs(i,3) > 0 ) then
00381 
00382                midp_x(i) = dble(0.25) * ( array_x(nbrs(i,1))                &
00383                                         + array_x(-nbrs(i,2)) + dble(360.0) &
00384                                         + array_x(nbrs(i,3))                &
00385                                         + array_x(nbrs(i,4)) )
00386 
00387             else if ( nbrs(i,2) > 0 .and. nbrs(i,3) < 0 ) then
00388 
00389                midp_x(i) = dble(0.25) * ( array_x(nbrs(i,1))                &
00390                                         + array_x(nbrs(i,2))                &
00391                                         + array_x(-nbrs(i,3)) + dble(360.0) &
00392                                         + array_x(nbrs(i,4)) )
00393 
00394             else if (  nbrs(i,2) < 0 .and. nbrs(i,3) < 0 ) then
00395 
00396                midp_x(i) = dble(0.5) * ( array_x(nbrs(i,1)) &
00397                                        + array_x(-nbrs(i,2)) + dble(360.0) )
00398             else
00399 
00400                midp_x(i) = dble(0.25) * ( array_x(nbrs(i,1)) &
00401                                         + array_x(nbrs(i,2)) &
00402                                         + array_x(nbrs(i,3)) &
00403                                         + array_x(nbrs(i,4)) )
00404             endif
00405 
00406          else
00407 
00408             if ( nbrs(i,2) < 0 ) then
00409                midp_x(i) = dble(0.5) * ( array_x(nbrs(i,1)) &
00410                                        + array_x(ibeg) + dble(360.0) )
00411             else
00412                midp_x(i) = dble(0.5) * ( array_x(nbrs(i,1)) &
00413                                        + array_x(nbrs(i,2)) )
00414             endif        
00415 
00416          endif
00417 
00418          if ( i == iend .and. i < range(2) ) then
00419             jlat = jlat + 1
00420             ibeg = iend + 1
00421             iend = iend + points_per_lat(jlat,1)
00422          endif
00423       enddo
00424 
00425 #undef __DEBUGX
00426 #ifdef __DEBUGX
00427 
00428       print *, ' bb_cells_gauss2 min, midp, max: '
00429       ii = 0
00430       do jlat = 1, nbr_lats
00431          do i = 1, points_per_lat(jlat,1)
00432             ii = ii + 1
00433             print '(i6,i4,x,6f9.4)', i, jlat, &
00434                  real(chmin_x(ii)), real(midp_x(ii)), real(chmax_x(ii)), 
00435                  real(chmin_y(ii)), real(midp_y(ii)), real(chmax_y(ii))
00436          enddo
00437       enddo
00438 
00439       print *, ' bb_cells_gauss2 x arrays: ', nbr_lats
00440       ii = 0
00441       do jlat = 1, nbr_lats
00442          do i = 1, points_per_lat(jlat,1)
00443             ii = ii + 1
00444             if (min(nbrs(ii,1),nbrs(ii,2),nbrs(ii,3),nbrs(ii,4)) > 0 ) &
00445                  print *, i, jlat, array_x(nbrs(ii,1)), &
00446                                    array_x(nbrs(ii,2)), &
00447                                    array_x(nbrs(ii,3)), &
00448                                    array_x(nbrs(ii,4))
00449          enddo
00450       enddo
00451 
00452       print *, ' bb_cells_gauss2 nbrs: ', nbr_lats
00453       ii = 0
00454       do jlat = 1, nbr_lats
00455          do i = 1, points_per_lat(jlat,1)
00456             ii = ii + 1
00457             print '(7i10)', i, jlat, ii, nbrs(ii,1),nbrs(ii,2),nbrs(ii,3),nbrs(ii,4)
00458          enddo
00459       enddo
00460 
00461 #endif
00462 
00463 !
00464 !===> All done
00465 !
00466 #ifdef VERBOSE
00467       print 9980, trim(ch_id), ierror
00468 
00469       call psmile_flushstd
00470 #endif /* VERBOSE */
00471 !
00472 !  Formats:
00473 !
00474 9990 format (1x, a, ': psmile_bbcells_gauss2_dble')
00475 9980 format (1x, a, ': psmile_bbcells_gauss2_dble: eof ierror =', i3)
00476 !
00477       end subroutine PSMILe_Bbcells_gauss2_dble