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