1 | include(dom.inc)
   2 | 
   3 |       SUBROUTINE FSCK_CASE(all_k_abs, w, celldata, nelmts, all_WVNB,    &
   4 |      &                     all_DWVNB, ngbands,                          &
   5 |      &                     nallbandes, nkabs, data_ref, WSGG_W,         &
   6 |      &                     ieltd, ieltf, WSGG_Wb, Tf, nbfaces)
   7 | 
   8 |       USE mod_inout
   9 |       USE mod_pmm
  10 | #ifdef USEPALM
  11 |       USE palmlib
  12 | #endif
  13 | 
  14 |       IMPLICIT NONE
  15 | 
  16 |       include 'dom_constants.h'
  17 | 
  18 | !     IN
  19 |       DOM_INT                          :: nallbandes
  20 |       DOM_INT                          :: nbfaces, ibnd
  21 |       DOM_INT                          :: nkabs, ieltd, ieltf
  22 |       DOM_INT                          :: nelmts, ngbands
  23 |       DOM_REAL,DIMENSION(nallbandes)   :: all_WVNB,all_DWVNB
  24 |       DOM_REAL,DIMENSION(8,nelmts)     :: celldata
  25 |       DOM_REAL,DIMENSION(8)            :: data_ref
  26 |       DOM_REAL, DIMENSION(nbfaces)     :: Tf
  27 | 
  28 | !     OUT
  29 |       DOM_REAL,DIMENSION(nkabs)        :: w
  30 |       DOM_REAL,DIMENSION(nkabs,nelmts) :: all_k_abs
  31 |       DOM_REAL,DIMENSION(nkabs,nelmts) :: WSGG_W
  32 |       DOM_REAL,DIMENSION(nkabs,nbfaces):: WSGG_Wb
  33 | 
  34 | !     LOCAL
  35 |       DOM_REAL,DIMENSION(nkabs,nelmts) :: local_all_k_abs, local_WSGG_W
  36 |       DOM_REAL,DIMENSION(ngbands,2)    :: gasdata
  37 |       DOM_REAL                         :: FSK_SOOT, K_SOOT 
  38 |       DOM_REAL,DIMENSION(nkabs)        :: k_absmel
  39 |       DOM_REAL,DIMENSION(ngbands)      :: F, F_ref
  40 |       DOM_REAL                         :: blae, planck 
  41 |       DOM_REAL                         :: WVNB,DWVNB,WVNB_SI,DWVNB_SI
  42 |       DOM_INT                          :: ielt, i, ierr, j
  43 |       DOM_INT                          :: i_bande, n_cl
  44 |       DOM_INT                          :: ICO,ICO2,IH2O
  45 |       LOGICAL                          :: LICO,LICO2,LIH2O
  46 | 
  47 |       DOM_REAL,DIMENSION(ngbands,2)    :: gasdata_ref
  48 |       DOM_REAL,DIMENSION(nkabs)        :: k_ref, gdek, fdek, fdek_ref
  49 |       DOM_REAL                         :: gdek_i, fdek_i, Tp
  50 | 
  51 | !     --------------------------------------------!
  52 | !     STEP 1 : Calculation of the reference state !
  53 | !              k_ref and fdek_ref                 !
  54 | !     --------------------------------------------!
  55 |   
  56 |       F           = 0.
  57 |       gasdata_ref = 0.
  58 |       fdek_ref    = 0.
  59 | 
  60 |       Tp = s_Tp
  61 | 
  62 |       DO i_bande=1,ngbands
  63 |         WVNB    = all_WVNB(i_bande)
  64 |         DWVNB   = all_DWVNB(i_bande)
  65 |         WVNB_SI = 100.*WVNB
  66 |         DWVNB_SI= 100.*DWVNB
  67 | 
  68 |         IF (WVNB.le.9300.0) THEN
  69 |           CALL FINDI(LICO,LICO2,LIH2O,ICO,ICO2,IH2O,WVNB,DWVNB)
  70 | 
  71 |           CALL KBARANDPHI(data_ref(:),LICO,LICO2,LIH2O,ICO,             &
  72 |      &                    ICO2,IH2O,                                    &
  73 |      &                    gasdata_ref(i_bande,1:2), nallbandes)
  74 |         ENDIF
  75 | 
  76 |         IF (s_Tp.eq.0) Tp = data_ref(1)
  77 |         F(i_bande) = blae(WVNB_SI, Tp)*DWVNB_SI 
  78 | 
  79 |       ENDDO
  80 | 
  81 |       F     = F / SUM (F)
  82 |       F_ref = F
  83 | 
  84 |       CALL k_distributeur(F, ngbands, gasdata_ref, nkabs, k_ref,w)
  85 | 
  86 |       DO j = 1, nkabs
  87 |         DO i_bande=1,ngbands
  88 | 
  89 |           CALL pdf(gasdata_ref(i_bande,2), gasdata_ref(i_bande,1),      &
  90 |      &             k_ref(j),fdek_i)
  91 | 
  92 |           fdek_ref(j) = fdek_ref(j) + F(i_bande)*fdek_i
  93 | 
  94 |         ENDDO
  95 |       ENDDO
  96 | 
  97 |       IF (pmm_rank.eq.0) THEN
  98 |         PRINT*, "      DATA_REF : T    = ",data_ref(1)
  99 |         PRINT*, "                 P    = ",data_ref(2)
 100 |         PRINT*, "                 XH2O = ",data_ref(3)
 101 |         PRINT*, "                 XCO2 = ", data_ref(4)
 102 |         PRINT*, "                 XCO  = ", data_ref(5)
 103 |         PRINT*
 104 | !       PRINT*, "K_ref    :", k_ref
 105 | !       PRINT*, "fdek_ref :", fdek_ref
 106 |       ENDIF
 107 | 
 108 | !     ----------------------------------------------------------!
 109 | !     Calculation of the absorption coefficient over the domain !
 110 | !     ----------------------------------------------------------!
 111 | 
 112 |       local_all_k_abs = 0.
 113 |       all_k_abs       = 0.
 114 |       local_WSGG_W    = 0.
 115 |       WSGG_W          = 0.
 116 | 
 117 | !     ---------------------------------------------!
 118 | !     Non-homogeneous system or homogeneous system !
 119 | !     ---------------------------------------------!
 120 |       IF (homosyst.eq.'NO') THEN
 121 |         n_cl=ieltf
 122 |       ELSEIF (homosyst.eq.'YES') THEN
 123 |         n_cl=ieltd
 124 |       ENDIF
 125 | 
 126 |       DO ielt=ieltd,n_cl
 127 | 
 128 |         k_absmel = 0.
 129 |         FSK_SOOT = 0.
 130 |         gasdata  = 0.
 131 |         F        = 0.
 132 | 
 133 | !       ------------------------------------!
 134 | !       Setting spectral data for each band !
 135 | !       ------------------------------------!
 136 | 
 137 |         DO i_bande=1,ngbands
 138 | 
 139 |           WVNB    = all_WVNB(i_bande)
 140 |           DWVNB   = all_DWVNB(i_bande)
 141 |           WVNB_SI = 100.*WVNB
 142 |           DWVNB_SI= 100.*DWVNB
 143 | 
 144 | !         ----------------------------------------!
 145 | !         Spectral index lecture for each species !
 146 | !         ----------------------------------------!
 147 |           IF (WVNB.le.9300.0) THEN
 148 | 
 149 |             CALL FINDI(LICO,LICO2,LIH2O,ICO,ICO2,IH2O,WVNB,DWVNB)
 150 | 
 151 |             CALL KBARANDPHI(celldata(:,ielt),LICO,LICO2,LIH2O,ICO,      &
 152 |      &                      ICO2,IH2O,                                  &
 153 |      &                      gasdata(i_bande,1:2), nallbandes)
 154 | 
 155 |           ENDIF
 156 | 
 157 | !         -------------------------------------!
 158 | !         Band groupment multiplicative factor !
 159 | !         -------------------------------------!
 160 | 
 161 |           IF (s_Tp.eq.0) Tp = celldata(1,ielt)
 162 |           F(i_bande) = blae(WVNB_SI,Tp)*DWVNB_SI 
 163 | 
 164 |           K_SOOT     = 5.5*celldata(8,ielt)*WVNB_SI
 165 |           FSK_SOOT   = FSK_SOOT + F(i_bande)*K_SOOT
 166 |  
 167 |         ENDDO
 168 | 
 169 |         FSK_SOOT = FSK_SOOT / SUM(F)
 170 |         F        = F        / SUM(F)
 171 | 
 172 | !       --------------------------------------!
 173 | !       STEP 2 : Calculate g(T,phi_ref,k_ref) !
 174 | !                and f(T,phi_ref,k_ref)       !
 175 | !       --------------------------------------!
 176 | 
 177 |         gdek = 0.
 178 |         fdek = 0.
 179 | 
 180 |         DO j = 1, nkabs
 181 |           DO i_bande=1,ngbands
 182 | 
 183 |             CALL cdss(gasdata_ref(i_bande,2), gasdata_ref(i_bande,1),   &
 184 |      &                k_ref(j), gdek_i)
 185 | 
 186 |             CALL pdf(gasdata_ref(i_bande,2), gasdata_ref(i_bande:,1),   &
 187 |      &               k_ref(j), fdek_i)
 188 | 
 189 |             gdek(j) = gdek(j) + F_ref(i_bande)*gdek_i
 190 |             fdek(j) = fdek(j) + F(i_bande)*fdek_i
 191 | 
 192 |           ENDDO
 193 |         ENDDO
 194 | 
 195 | !         --------------------------------------------------------------!
 196 | !         STEP 3 : g(T_ref,phi_ref,k) = g(T_ref,phi,k*) => Correlated k !
 197 | !                  Invert g in the RHS to k* (k_absmel)                 !
 198 | !         --------------------------------------------------------------!
 199 | 
 200 |         DO j = 1, nkabs
 201 | 
 202 |           CALL COFG (F_ref, ngbands, gasdata(:,2), gasdata(:,1),gdek(j),&
 203 |      &               gdek(nkabs), k_absmel(j))
 204 | 
 205 |         ENDDO
 206 | 
 207 |         local_all_k_abs(:,ielt) = k_absmel + FSK_SOOT
 208 | 
 209 | !       ----------------------------------------------------------------!
 210 | !       STEP 4 : Calculation of the weight function                     !
 211 | !                a(T,T_ref,g_ref)= fdek(T,phi_ref)/fdek_(T_ref,phi_ref) !
 212 | !       ----------------------------------------------------------------!
 213 | 
 214 |         DO i = 1, nkabs
 215 |           IF(fdek_ref(i).ne.0) THEN
 216 |             local_WSGG_W(i,ielt) = fdek(i)/fdek_ref(i)
 217 |           ELSE
 218 |             local_WSGG_W(i,ielt) = 1.
 219 |           ENDIF
 220 |         ENDDO
 221 | 
 222 |         IF ((MOD(ielt,100).eq.0).and.(pmm_rank.eq.0)) THEN
 223 |           print*, "      >> Node:",ielt-ieltd,"/",n_cl-ieltd
 224 |         ENDIF
 225 | 
 226 |       ENDDO
 227 | 
 228 |       IF (homosyst.eq.'YES') THEN
 229 |         DO i=ieltd+1,ieltf
 230 |           local_all_k_abs(:,i) = local_all_k_abs(:,ieltd)
 231 |           local_WSGG_W(:,ielt) = local_WSGG_W(:,ieltd)
 232 |         ENDDO
 233 |       ENDIF
 234 | 
 235 | !     ------------------------------------------------------------------!
 236 | !     STEP 5 : Calculation of the weight function at walls              !
 237 | !              a(Tw,T_ref,g_ref)= fdek(Tw,phi_ref)/fdek_(T_ref,phi_ref) !
 238 | !     ------------------------------------------------------------------!
 239 | 
 240 |       WSGG_Wb = 0.
 241 | 
 242 |       DO ibnd=1,nbfaces
 243 | 
 244 |         fdek = 0.
 245 | 
 246 |         DO i_bande=1,ngbands
 247 |           DWVNB   = all_DWVNB(i_bande)
 248 |           DWVNB_SI= 100.*DWVNB
 249 | 
 250 |           IF(s_TP.eq.0) Tp = Tf(ibnd)
 251 |           F(i_bande) = blae(WVNB_SI,Tp)*DWVNB_SI 
 252 | 
 253 |         ENDDO
 254 | 
 255 |         F = F / SUM (F)
 256 | 
 257 |         DO j = 1, nkabs
 258 |           DO i_bande=1,ngbands
 259 | 
 260 |             CALL pdf(gasdata_ref(i_bande,2), gasdata_ref(i_bande:,1),   &
 261 |      &               k_ref(j),fdek_i)
 262 | 
 263 |             fdek(j) = fdek(j) + F(i_bande)*fdek_i
 264 | 
 265 |           ENDDO
 266 | 
 267 |           IF(fdek_ref(j).ne.0) THEN
 268 |             WSGG_Wb(j,ibnd) = fdek(j)/fdek_ref(j)
 269 |           ELSE
 270 |             WSGG_Wb(j,ibnd) = 1.
 271 |           ENDIF          
 272 | 
 273 |         ENDDO
 274 |       ENDDO
 275 | 
 276 | !     ------------------------------------------------!
 277 | !     sending all_k_abs....change to a pmm subroutine !
 278 | !     ------------------------------------------------!
 279 | 
 280 |        CALL MPI_ALLREDUCE(local_all_k_abs, all_k_abs, nkabs*nelmts,     &
 281 |      &                    MPI_DOUBLE_PRECISION, MPI_SUM, COMM_PARA ,    &
 282 |      &                    ierr)
 283 | 
 284 |        CALL MPI_ALLREDUCE(local_WSGG_W, WSGG_W, nkabs*nelmts,           &
 285 |      &                    MPI_DOUBLE_PRECISION, MPI_SUM, COMM_PARA ,    &
 286 |      &                    ierr)
 287 | 
 288 |       IF (homosyst.eq.'YES') THEN
 289 |         DO j = 1, nkabs
 290 |           all_k_abs(j,:) = k_ref(j) 
 291 |         ENDDO
 292 |         WSGG_W = 1.d0
 293 |       ENDIF
 294 | 
 295 |       END SUBROUTINE FSCK_CASE