prismtrs_apply_grads.F90

Go to the documentation of this file.
00001 !------------------------------------------------------------------------
00002 ! Copyright 2006-2010, NEC Europe Ltd., London, UK.E
00003 ! All rights reserved. Use is subject to OASIS4 license terms.
00004 !-----------------------------------------------------------------------
00005 !BOP
00006 !
00007 ! !ROUTINE: PRISMTrs_Apply_grads
00008 !
00009 ! !INTERFACE
00010 
00011 subroutine prismtrs_apply_grads(  il_src_size,         &
00012                                   dda_trans_out,       &
00013                                   ida_src_mask,        &
00014                                   il_tgt_size,         &
00015                                   dda_trans_in,        &
00016                                   ida_tgt_mask,        &
00017                                   id_nb_neighbors,     &
00018                                   ida_neighbors,       &
00019                                   dda_weights,         &
00020                                   nbr_fields,          &
00021                                   id_err)
00022 !
00023 ! !USES:
00024 !
00025   USE PRISMDrv, dummy_interface => PRISMTrs_Apply_grads
00026 
00027   IMPLICIT NONE
00028 
00029 !
00030 ! !PARAMETERS:
00031 !
00032 ! size of the transient in
00033   INTEGER, INTENT (IN)                              :: il_src_size
00034 
00035 ! Nbr of bundle components
00036   INTEGER, INTENT (IN)                              :: nbr_fields
00037 
00038 ! transient in field
00039   DOUBLE PRECISION, DIMENSION(il_src_size*nbr_fields), INTENT(IN) :: dda_trans_out
00040  
00041 ! src mask
00042   INTEGER, DIMENSION(il_src_size), INTENT(IN)       :: ida_src_mask
00043 
00044 ! size of the transient out
00045   INTEGER, INTENT (IN)                              :: il_tgt_size
00046 
00047 ! tgt mask
00048   INTEGER, DIMENSION(il_tgt_size), INTENT(IN)       :: ida_tgt_mask
00049   
00050 
00051 ! number of neighbors involved in the interpolation
00052   INTEGER, INTENT (IN)                              :: id_nb_neighbors
00053  
00054 ! Neighbors and associated weights
00055   INTEGER, DIMENSION(il_tgt_size,id_nb_neighbors), INTENT (IN) :: 
00056      ida_neighbors
00057   DOUBLE PRECISION, DIMENSION(il_tgt_size,id_nb_neighbors), INTENT(IN) :: dda_weights
00058 
00059 !
00060 ! ! RETURN VALUE
00061 !
00062 ! transient out field
00063   DOUBLE PRECISION, DIMENSION(il_tgt_size*nbr_fields), INTENT (Out) :: dda_trans_in
00064 
00065   INTEGER, INTENT (Out)               :: id_err   ! error value
00066 
00067   DOUBLE PRECISION :: ctrl
00068 
00069   DOUBLE PRECISION, DIMENSION(16)   :: data    ! only for id_nb_neighbors==16
00070   INTEGER, DIMENSION(16)            :: tmpmask ! only for id_nb_neighbors==16
00071   DOUBLE PRECISION, DIMENSION(4)    :: grad_i  ! only for id_nb_neighbors==16
00072   DOUBLE PRECISION, DIMENSION(4)    :: grad_j  ! only for id_nb_neighbors==16
00073   DOUBLE PRECISION, DIMENSION(4)    :: grad_ij ! only for id_nb_neighbors==16
00074 
00075 ! !DESCRIPTION
00076 ! Subroutine "PRISMTrs_Interp_Apply_grads" performs a weighted gradient 
00077 ! computation for bicubic interpolation to get the field on the target
00078 ! grid.
00079 !
00080 ! !REVISED HISTORY
00081 !   Date      Programmer   Description
00082 ! ----------  ----------   -----------
00083 ! 06/05/2005  R. Redler    Creation
00084 ! 06/03/2008  J. Charles   Modifications added for the use of bundle fields
00085 !
00086 !
00087 ! TODO: clean implementation at non-cyclic edges
00088 !
00089 ! EOP
00090 !----------------------------------------------------------------------
00091 ! $Id: prismtrs_apply_grads.F90 2685 2010-10-28 14:05:10Z coquart $
00092 ! $Author: coquart $
00093 !----------------------------------------------------------------------
00094 !
00095 ! Local declarations
00096 !
00097   CHARACTER(LEN=len_cvs_string), SAVE  :: mycvs = 
00098      '$Id: prismtrs_apply_grads.F90 2685 2010-10-28 14:05:10Z coquart $'
00099 
00100   INTEGER :: ib, ib_bis, i
00101   INTEGER :: il_count_neighbors
00102 !
00103 ! ---------------------------------------------------------------------
00104 !
00105 #ifdef VERBOSE
00106   PRINT *, '| | | | | | | Enter PRISMTrs_Apply_grads'
00107   call psmile_flushstd
00108 #endif
00109 
00110   id_err = 0
00111 
00112   IF ( id_nb_neighbors /= 16 ) THEN
00113      PRINT *, &
00114            '| | | | | | | | application of bicubic grads restricted to 16 points.'
00115      PRINT *, &
00116            '| | | | | | | | but number of points is ', id_nb_neighbors
00117 
00118      call psmile_abort
00119 
00120   ENDIF
00121 !
00122 ! 1. For each non masked target point Apply the weights and gradients
00123 !    on the source field to get the target field
00124 !
00125   dda_trans_in(:) = 0.
00126   
00127   DO i = 1, nbr_fields
00128 
00129      DO ib = 1, il_tgt_size
00130 
00131         ctrl = SUM(dda_weights(ib,:))
00132 
00133 #ifdef DEBUGX
00134         IF ( (ib == 30) .or. (ib == 1667) ) THEN
00135             PRINT*, '++++++++++++++++++++++++++++++++++++++++++++++++++'
00136             PRINT*, ' SUM(weights(ib,1:16)) for EPIOT number ib :',ib
00137             PRINT*, '++++++++++++++++++++++++++++++++++++++++++++++++++'
00138             PRINT*, ctrl
00139         ENDIF
00140 #endif
00141 
00142         IF (ida_tgt_mask(ib) == 1) THEN
00143 
00144         IF ( ctrl < 0.0 ) THEN
00145 
00146 ! We have to perform a nearest neighbor weight since weights
00147 ! are marked as negative. Note that for this reason the partial
00148 ! contributions are subtracted instead of added.
00149 
00150            il_count_neighbors = 0
00151 
00152            DO ib_bis = 1, id_nb_neighbors
00153               IF (ida_neighbors(ib,ib_bis) /= 0) THEN
00154                  dda_trans_in(ib+(i-1)*il_tgt_size) = dda_trans_in(ib+(i-1)*il_tgt_size) -  &
00155                       dda_weights(ib,ib_bis) *          &
00156                       dda_trans_out(ida_neighbors(ib,ib_bis)+(i-1)*il_src_size)
00157               ELSE
00158                  il_count_neighbors = il_count_neighbors + 1
00159               END IF
00160 
00161               IF (il_count_neighbors == id_nb_neighbors) &
00162                    dda_trans_in(ib+(i-1)*il_tgt_size) = PSMILe_dundef
00163 
00164            END DO
00165 
00166         ELSE IF ( ctrl > 0 ) THEN
00167 
00168 ! Apply weights and gradients
00169 !
00170 !           13   14   15   16
00171 !
00172 !     ^      9   10   11   12
00173 !     |
00174 !     |      5    6    7    8
00175 !
00176 !     j      1    2    3    4
00177 !
00178 !       i ---->
00179 
00180 ! We go counterclockwise over the points
00181 ! The "inner" weights refer to indices 6, 7, 11 and 10
00182 ! i,j reference is on index 6, ==> (im1,jp1) = (1,9)
00183 ! Then we do the gradients for 6, 7, 11 and 10  
00184 
00185            do ib_bis = 1, 16
00186              data(ib_bis)    = dda_trans_out(ida_neighbors(ib,ib_bis)+(i-1)*il_src_size)
00187              tmpmask(ib_bis) = ida_src_mask (ida_neighbors(ib,ib_bis))
00188            enddo
00189 
00190            call gradient  ( tmpmask, data, grad_i, grad_j, grad_ij  )
00191 
00192            dda_trans_in(ib+(i-1)*il_tgt_size) = dda_weights(ib, 1) * dda_trans_out(ida_neighbors(ib, 6)+(i-1)*il_src_size)   &
00193                             + dda_weights(ib, 5) * dda_trans_out(ida_neighbors(ib, 7)+(i-1)*il_src_size)   &
00194                             + dda_weights(ib, 9) * dda_trans_out(ida_neighbors(ib,11)+(i-1)*il_src_size)   &
00195                             + dda_weights(ib,13) * dda_trans_out(ida_neighbors(ib,10)+(i-1)*il_src_size)   &
00196 ! gradient i       (ip1,j) - (im1,j)
00197                             + dda_weights(ib, 2) * grad_i(1) &
00198                             + dda_weights(ib, 6) * grad_i(2) &
00199                             + dda_weights(ib,10) * grad_i(3) &
00200                             + dda_weights(ib,14) * grad_i(4) &
00201 ! gradient j       (i,jp1) - (i,jm1)
00202                             + dda_weights(ib, 3) * grad_j(1) &
00203                             + dda_weights(ib, 7) * grad_j(2) &
00204                             + dda_weights(ib,11) * grad_j(3) &
00205                             + dda_weights(ib,15) * grad_j(4) &
00206 ! gradient ij      ((ip1,jp1) - (im1,jp1)) - ((ip1,jm1) - (im1,jm1))
00207                             + dda_weights(ib, 4) * grad_ij(1) &
00208                             + dda_weights(ib, 8) * grad_ij(2) &
00209                             + dda_weights(ib,12) * grad_ij(3) &
00210                             + dda_weights(ib,16) * grad_ij(4)
00211 
00212 
00213 #ifdef DEBUGX
00214         IF ( (ib == 30) .or. (ib == 1667) ) THEN
00215             PRINT* 
00216             PRINT*, '++++++++++++++++++++++++++++++++++++++++++++'
00217             PRINT*, 'weights(ib,1:16) for EPIOT number ib :',ib
00218             PRINT*, '++++++++++++++++++++++++++++++++++++++++++++'
00219             PRINT*, dda_weights(ib,1),dda_weights(ib,2),dda_weights(ib,3),dda_weights(ib,4)
00220             PRINT*, dda_weights(ib,5),dda_weights(ib,6),dda_weights(ib,7),dda_weights(ib,8)
00221             PRINT*, dda_weights(ib,9),dda_weights(ib,10),dda_weights(ib,11),dda_weights(ib,12)
00222             PRINT*, dda_weights(ib,13),dda_weights(ib,14),dda_weights(ib,15),dda_weights(ib,16)
00223             PRINT*, '++++++++++++++++++++++++++++++++++++++++++++'
00224             PRINT*, 'gradi,gradj,gradij for EPIOT number ib :',ib
00225             PRINT*, '++++++++++++++++++++++++++++++++++++++++++++'
00226             PRINT*, grad_i(1),grad_i(2),grad_i(3),grad_i(4)
00227             PRINT*, grad_j(1),grad_j(2),grad_j(3),grad_j(4)
00228             PRINT*, grad_ij(1),grad_ij(2),grad_ij(3),grad_ij(4)
00229             PRINT*
00230         ENDIF
00231 #endif
00232 
00233         ELSE
00234 
00235            dda_trans_in(ib+(i-1)*il_tgt_size) = PSMILe_dundef
00236            
00237         ENDIF
00238 
00239      ELSE
00240 !
00241 ! 2. Else give the mask value to the point
00242 !
00243         dda_trans_in(ib+(i-1)*il_tgt_size) = PSMILe_dundef
00244 
00245      END IF
00246 
00247   END DO
00248 
00249   END DO
00250 !
00251 #ifdef VERBOSE
00252   PRINT *, '| | | | | | | Quit PRISMTrs_Apply_grads'
00253   call psmile_flushstd
00254 #endif
00255 
00256 END SUBROUTINE PRISMTrs_Apply_grads
00257 
00258 
00259 
00260 subroutine gradient (mask, data, grad_i, grad_j, grad_ij )
00261 
00262    implicit none
00263 
00264    double precision :: grad_i  (1:4)
00265    double precision :: grad_j  (1:4)
00266    double precision :: grad_ij (1:4)
00267    double precision :: data   (1:16)
00268 
00269    double precision :: weight, weight_i, weight_j 
00270 
00271    double precision :: grad_jp1  ! (ip1,jp1) - (im1,jp1)
00272    double precision :: grad_jm1  ! (ip1,jm1) - (im1,jm1)
00273 
00274    integer :: index (1:4)
00275    integer :: mask(16)
00276 
00277    integer :: i, ip1, im1, jp1, jm1
00278 
00279    integer :: index_jp1 (1:4)
00280    integer :: index_jm1 (1:4)
00281 
00282 
00283    data index / 6, 7, 11, 10 / 
00284 
00285    data index_jp1 / 10, 11, 15, 14 /
00286    data index_jm1 /  2,  3,  7,  6 / 
00287 !
00288 !  gradient_i
00289 !
00290 
00291 !!$                            + dda_weights(ib, 2) * (dda_trans_out(ida_neighbors(ib, 7))  &
00292 !!$                                                  - dda_trans_out(ida_neighbors(ib, 5))) &
00293 !!$                                                 * 0.5_8                                 &
00294 !!$                            + dda_weights(ib, 6) * (dda_trans_out(ida_neighbors(ib, 8))  &
00295 !!$                                                  - dda_trans_out(ida_neighbors(ib, 6))) &
00296 !!$                                                 * 0.5_8                                 &
00297 !!$                            + dda_weights(ib,10) * (dda_trans_out(ida_neighbors(ib,12))  &
00298 !!$                                                  - dda_trans_out(ida_neighbors(ib,10))) &
00299 !!$                                                 * 0.5_8                                 &
00300 !!$                            + dda_weights(ib,14) * (dda_trans_out(ida_neighbors(ib,11))  &
00301 !!$                                                  - dda_trans_out(ida_neighbors(ib, 9))) &
00302 !!$                                                 * 0.5_8                                 &
00303 
00304 
00305    do i = 1, 4
00306 
00307      ip1 = index(i)+1
00308      im1 = index(i)-1
00309 
00310      weight = 0.5d0
00311 
00312      if ( mask(ip1) == 1 .or. mask(im1) == 1 ) then
00313         if ( mask(index(i)+1) /= 1 ) then
00314            ip1 = index(i)
00315            weight = 1.0d0
00316         else if ( mask(index(i)-1) /= 1 ) then
00317            im1 = index(i)
00318            weight = 1.0d0
00319         endif
00320         grad_i(i) = weight * ( data(ip1) - data(im1) )
00321      else
00322         grad_i(i) = 0.0d0
00323      endif
00324 
00325    enddo
00326 
00327 !
00328 !  gradient_j
00329 !
00330 
00331 !!$                            + dda_weights(ib, 3) * (dda_trans_out(ida_neighbors(ib,10))  &
00332 !!$                                                  - dda_trans_out(ida_neighbors(ib, 2))) &
00333 !!$                                                 * 0.5d0                                 &
00334 !!$                            + dda_weights(ib, 7) * (dda_trans_out(ida_neighbors(ib,11))  &
00335 !!$                                                  - dda_trans_out(ida_neighbors(ib, 3))) &
00336 !!$                                                 * 0.5d0                                 &
00337 !!$                            + dda_weights(ib,11) * (dda_trans_out(ida_neighbors(ib,15))  &
00338 !!$                                                  - dda_trans_out(ida_neighbors(ib, 7))) &
00339 !!$                                                 * 0.5d0                                 &
00340 !!$                            + dda_weights(ib,15) * (dda_trans_out(ida_neighbors(ib,14))  &
00341 !!$                                                  - dda_trans_out(ida_neighbors(ib, 6))) &
00342 !!$                                                 * 0.5d0                                 &
00343 
00344 
00345    do i = 1, 4
00346 
00347      jp1 = index(i)+4
00348      jm1 = index(i)-4
00349 
00350      weight = 0.5d0
00351 
00352      if ( mask(jp1) == 1 .or. mask(jm1) == 1 ) then
00353         if ( mask(index(i)+4) /= 1 ) then
00354            jp1 = index(i)
00355            weight = 1.0d0
00356         else if ( mask(index(i)-4) /= 1 ) then
00357            jm1 = index(i)
00358            weight = 1.0d0
00359         endif
00360         grad_j(i) = weight * ( data(jp1) - data(jm1) )
00361      else
00362         grad_j(i) = 0.0d0
00363      endif
00364 
00365    enddo
00366 
00367 !
00368 !  gradient_ij
00369 !
00370 
00371 !!$                            + dda_weights(ib, 4) * (dda_trans_out(ida_neighbors(ib,11))  &
00372 !!$                                                  - dda_trans_out(ida_neighbors(ib, 9))  &
00373 !!$                                                  - dda_trans_out(ida_neighbors(ib, 3))  &
00374 !!$                                                  + dda_trans_out(ida_neighbors(ib, 1))) &
00375 !!$                                                 * 0.25d0                                &
00376 !!$                            + dda_weights(ib, 8) * (dda_trans_out(ida_neighbors(ib,12))  &
00377 !!$                                                  - dda_trans_out(ida_neighbors(ib,10))  &
00378 !!$                                                  - dda_trans_out(ida_neighbors(ib, 4))  &
00379 !!$                                                  + dda_trans_out(ida_neighbors(ib, 2))) &
00380 !!$                                                 * 0.25d0                                &
00381 !!$                            + dda_weights(ib,12) * (dda_trans_out(ida_neighbors(ib,16))  &
00382 !!$                                                  - dda_trans_out(ida_neighbors(ib,14))  &
00383 !!$                                                  - dda_trans_out(ida_neighbors(ib, 8))  &
00384 !!$                                                  + dda_trans_out(ida_neighbors(ib, 6))) &
00385 !!$                                                 * 0.25d0                                &
00386 !!$                            + dda_weights(ib,16) * (dda_trans_out(ida_neighbors(ib,15))  &
00387 !!$                                                  - dda_trans_out(ida_neighbors(ib,13))  &
00388 !!$                                                  - dda_trans_out(ida_neighbors(ib, 7))  &
00389 !!$                                                  + dda_trans_out(ida_neighbors(ib, 5))) &
00390 !!$                                                 * 0.25d0
00391 
00392 
00393    weight_j = 0.5d0
00394    grad_ij  = 0.0d0
00395 
00396    do i = 1, 4
00397 
00398      ip1 = index_jp1(i)+1
00399      im1 = index_jp1(i)-1
00400 
00401      weight_i = 0.5d0
00402      grad_jp1 = 0.0d0
00403 
00404      if ( mask(index_jp1(i)+1) == 1 .or. mask(index_jp1(i)-1) == 1 ) then
00405         if ( mask(index_jp1(i)+1) /= 1 .and. mask(index_jp1(i)) == 1 ) then
00406            ip1 = index_jp1(i)
00407            weight_i = 1.0d0
00408         else if ( mask(index_jp1(i)-1) /= 1 .and. mask(index_jp1(i)) == 1 ) then
00409            im1 = index_jp1(i)
00410            weight_i = 1.0d0
00411         else
00412            weight_i = 0.0d0
00413         endif
00414         grad_jp1 = weight_i * ( data(ip1) - data(im1) )
00415      endif
00416 
00417      ip1 = index_jm1(i)+1
00418      im1 = index_jm1(i)-1
00419 
00420      weight_i = 0.5d0
00421      grad_jm1 = 0.0d0
00422 
00423      if ( mask(index_jm1(i)+1) == 1 .or. mask(index_jm1(i)-1) == 1 ) then
00424         if ( mask(index_jm1(i)+1) /= 1 .and. mask(index_jm1(i)) == 1 ) then
00425            ip1 = index_jm1(i)
00426            weight_i = 1.0d0
00427         else if ( mask(index_jm1(i)-1) /= 1 .and. mask(index_jm1(i)) == 1 ) then
00428            im1 = index_jm1(i)
00429            weight_i = 1.0d0
00430         else
00431            weight_i = 0.0d0
00432         endif
00433         grad_jm1 = weight_i * ( data(ip1) - data(im1) )
00434      endif
00435 
00436      if ( grad_jp1 /= 0.0d0 .and. grad_jm1 /= 0.0d0 ) then
00437        grad_ij(i) = weight_j * ( grad_jp1 - grad_jm1 )
00438      endif
00439 
00440    enddo
00441 
00442 end subroutine gradient