1 | include(dom.inc)
   2 | 
   3 |       SUBROUTINE BAND_INTEG_SNB(wkabs, kabs, kscat, Lb, Lo, DWVNB_SI,   &
   4 |      &                      Gtot, Lbtot,  Htot, Q_rtot, Q_ptot, ifreq,  &
   5 |      &                      ntot_iter)
   6 | 
   7 |         USE mod_slave
   8 |         USE mod_inout
   9 |         USE mod_pmm
  10 | 
  11 |         IMPLICIT NONE
  12 | 
  13 |         include 'dom_constants.h'
  14 |         include 'pmm_constants.h'
  15 | 
  16 | !       IN
  17 |         DOM_INT                               :: ifreq
  18 |         DOM_REAL                              :: DWVNB_SI
  19 |         DOM_REAL,DIMENSION(is_nkabs)          :: wkabs
  20 |         DOM_REAL,DIMENSION(is_nkabs,is_ncells):: kabs
  21 |         DOM_REAL,DIMENSION(is_ncells)         :: Lb, kscat
  22 |         DOM_REAL,DIMENSION(is_nbfaces)        :: Lo
  23 | 
  24 | !       OUT
  25 |         DOM_REAL,DIMENSION(3,is_ncells)          :: Q_rtot
  26 |         DOM_REAL,DIMENSION(is_ncells)            :: Gtot, Lbtot
  27 |         DOM_REAL,DIMENSION(is_nbfaces)           :: Htot
  28 |         DOM_REAL,DIMENSION(3,is_nprobes)         :: Q_ptot
  29 | 
  30 | !       LOCAL
  31 |         DOM_INT                                  :: iquad, i, ntot_iter
  32 |         DOM_INT                                  :: n_iter, maxiter
  33 |         DOM_REAL                                 :: error,error1,error2
  34 |         DOM_REAL,DIMENSION(is_ncells)            :: G, Gi
  35 |         DOM_REAL,DIMENSION(is_nbfaces)           :: H, Hi
  36 |         DOM_REAL,DIMENSION(3,is_ncells)          :: Qr
  37 |         DOM_REAL,DIMENSION(3,is_nprobes)         :: Qp
  38 |         DOM_REAL,ALLOCATABLE, DIMENSION(:,:,:)   :: Livirt
  39 | !       DOM_INT                                  :: mloc(1) , mloc2(3)
  40 | 
  41 |         IF((i_dom_npart.gt.1).and.(.not.ALLOCATED(Livirt))) THEN
  42 |           ALLOCATE(Livirt(is_ntot_vfaces,is_ntotdir, 1))
  43 |         ENDIF
  44 | 
  45 | !       print*, " (",pmm_rank,") Doing band integration"
  46 | 
  47 | !       --------------------------!
  48 | !       Integrate the narrow band !
  49 | !       --------------------------!
  50 | 
  51 |         ntot_iter = 0
  52 | 
  53 |         DO iquad=1,is_nkabs
  54 | 
  55 | 
  56 |          IF( (pmm_rank.eq.0)                                            &
  57 |      &     .and.( ( (SUM(ts_boundary%epsil)/is_nbfaces.ne.1)            &
  58 |      &            .or.(i_dom_npart.ne.1) ) ))                          &
  59 |      &   print*, "       - Nq=",iquad,"/",is_nkabs
  60 | 
  61 | !         -----------------------------------------------------!
  62 | !         Initialisation values before reflection subiteration !
  63 | !         -----------------------------------------------------!
  64 | 
  65 |           G  = 0.
  66 |           H  = 0.
  67 | 
  68 |           error   = 1000.
  69 |           n_iter  = 0
  70 |           maxiter = 50
  71 | 
  72 |           DO WHILE ((error.gt.critconv).and.(n_iter.lt.maxiter))
  73 | 
  74 | !           -------------------------------------------!
  75 | !           Update values for reflection sub iteration !
  76 | !           -------------------------------------------!
  77 | 
  78 |             Gi = G
  79 |             Hi = H
  80 |             IF(i_dom_npart.gt.1) Livirt(:,:,1)=s_Lvirt(:,:,iquad,ifreq)
  81 | 
  82 | !           ---------------!
  83 | !           SPATIAL SCHEME !
  84 | !           ---------------!
  85 | 
  86 |             CALL SPATIAL_SCHEME(I_SCHEME, Lb, Gi, G, kabs(iquad,:),     &
  87 |      &                          kscat, Qr, H, alpha, Qp, Lo, ifreq,     &
  88 |      &                          iquad, i_dom_nthread, i_dom_npart)
  89 | 
  90 | !           --------------------------------!
  91 | !           Convergence test for reflection !
  92 | !           --------------------------------!
  93 | 
  94 |             IF(i_dom_npart.eq.1) THEN
  95 | 
  96 | !              mloc  = MAXLOC(abs(H-Hi))
  97 | !              IF (abs(Hi(mloc(1))).ne.0.) THEN
  98 | !                error1 = (abs(H(mloc(1))-Hi(mloc(1))))/abs(Hi(mloc(1)))
  99 |               error2 = 0.
 100 |               IF(SUM(Hi).ne.0) THEN
 101 |                 error1 = abs(SUM(Hi)-SUM(H))/SUM(Hi)
 102 |               ELSEIF(SUM(ts_boundary%epsil)/is_nbfaces.eq.1) THEN
 103 |                 error1 = 0.
 104 |               ELSE
 105 |                 error1 = 1000.
 106 |               ENDIF
 107 | 
 108 |             ELSE
 109 | !              mloc2  = MAXLOC(abs(s_Lvirt(:,:,:,ifreq)-Livirt))
 110 | !              IF (abs(Livirt(mloc2(1),mloc2(2),mloc2(3))).ne.0.) THEN
 111 | !                error2 = (abs(s_Lvirt(mloc2(1),mloc2(2),mloc2(3),1)-    &
 112 | !     &                      Livirt(mloc2(1),mloc2(2),mloc2(3)))   )     &
 113 | !     &                 /abs(Livirt(mloc2(1),mloc2(2),mloc2(3)))
 114 |               IF (SUM(Livirt).ne.0.) THEN
 115 |                error2 = abs(SUM(s_Lvirt(:,:,iquad,ifreq))-SUM(Livirt))   &
 116 |      &                  /SUM(Livirt)
 117 |               ELSE
 118 |                 error2 = 1000.
 119 |               ENDIF
 120 |           ENDIF
 121 | 
 122 |           error=MAX(error1,error2)
 123 |           IF(SUM(kabs(iquad,:)).eq.0) error = 0
 124 |           n_iter = n_iter + 1
 125 | 
 126 |            IF ((pmm_rank.eq.0).and.(error.ge.1e-15))                    &
 127 |      &    WRITE(*,'(A19XI2XA7XF12.10)')"   *** n_iter=",n_iter,         &
 128 |      &    " - err=", error
 129 |           IF((i_inter.eq.1).and.(pmm_rank.eq.0)) print*
 130 | 
 131 | !           IF (pmm_rank.eq.0) print*, "mloc=",mloc(1),"- H=",          &
 132 | !    &                         H(mloc(1)),"- Hi=",Hi(mloc(1))
 133 | 
 134 |           ENDDO
 135 | 
 136 | !         -------------------------------!
 137 | !         Sum over band quadrature point !
 138 | !         -------------------------------!
 139 | 
 140 |           Gtot(:)  = Gtot(:)  + G(:)*kabs(iquad,:)*                   &
 141 |      &                 wkabs(iquad)*DWVNB_SI
 142 |           Lbtot(:) = Lbtot(:) + Lb(:)*4*pi*kabs(iquad,:)*             &
 143 |      &                 wkabs(iquad)*DWVNB_SI
 144 |           Htot(:)  = Htot(:)  + H(:)*wkabs(iquad)*DWVNB_SI
 145 |           DO i=1,3
 146 |             Q_rtot(i,:) = Q_rtot(i,:)+Qr(i,:)*wkabs(iquad)*DWVNB_SI
 147 |             Q_ptot(i,:) = Q_ptot(i,:)+Qp(i,:)*wkabs(iquad)*DWVNB_SI
 148 |           ENDDO
 149 | 
 150 |         ntot_iter = ntot_iter + n_iter
 151 | 
 152 |         ENDDO
 153 | 
 154 |       END SUBROUTINE BAND_INTEG_SNB