1 | include(dom.inc)
   2 | 
   3 |       SUBROUTINE FSCK_CASE(all_k_abs, w, celldata, nelmts, all_WVNB,    &
   4 |      &                     all_DWVNB, KCO, DCO, KC, DC, KH, DH, ngbands,&
   5 |      &                     nallbandes, nkabs, data_ref, WSGG_W,homosyst,&
   6 |      &                     s_Tp)
   7 | 
   8 |       IMPLICIT NONE
   9 | 
  10 | !     IN
  11 |       DOM_REAL, PARAMETER              :: pi = 3.141592654
  12 |       CHARACTER*3                      :: homosyst 
  13 |       DOM_INT                          :: nallbandes
  14 |       DOM_INT                          :: nkabs
  15 |       DOM_INT                          :: nelmts, ngbands
  16 |       DOM_REAL,DIMENSION(14,nallbandes):: KCO,DCO,KC,DC,KH,DH
  17 |       DOM_REAL,DIMENSION(nallbandes)   :: all_WVNB,all_DWVNB
  18 |       DOM_REAL,DIMENSION(nelmts,8)     :: celldata
  19 |       DOM_REAL,DIMENSION(8)            :: data_ref
  20 | 
  21 | 
  22 | !     OUT
  23 |       DOM_REAL,DIMENSION(nkabs)        :: w
  24 |       DOM_REAL,DIMENSION(nelmts,nkabs) :: all_k_abs
  25 |       DOM_REAL,DIMENSION(nelmts,nkabs) :: WSGG_W
  26 | 
  27 | !     LOCAL
  28 |       DOM_REAL,DIMENSION(2,ngbands)    :: gasdata
  29 |       DOM_REAL                         :: FSK_SOOT, K_SOOT 
  30 |       DOM_REAL,DIMENSION(nkabs)        :: k_absmel
  31 |       DOM_REAL,DIMENSION(ngbands)      :: F, F_ref
  32 |       DOM_REAL                         :: blae, planck 
  33 |       DOM_REAL                         :: WVNB,DWVNB,WVNB_SI,DWVNB_SI
  34 |       DOM_INT                          :: ielt, i, ierr, j
  35 |       DOM_INT                          :: i_bande, n_cl
  36 |       DOM_INT                          :: ICO,ICO2,IH2O
  37 |       LOGICAL                          :: LICO,LICO2,LIH2O
  38 | 
  39 |       DOM_REAL,DIMENSION(2,ngbands)    :: gasdata_ref
  40 |       DOM_REAL,DIMENSION(nkabs)        :: k_ref, gdek, fdek, fdek_ref
  41 |       DOM_REAL                         :: gdek_i, fdek_i, s_Tp, Tp
  42 | 
  43 | !     --------------------------------------------!
  44 | !     STEP 1 : Calculation of the reference state !
  45 | !              k_ref and fdek_ref                 !
  46 | !     --------------------------------------------!
  47 |   
  48 |       F           = 0.
  49 |       gasdata_ref = 0.
  50 |       fdek_ref    = 0.
  51 | 
  52 |       Tp = s_Tp
  53 | 
  54 |       DO i_bande=1,ngbands
  55 |         WVNB    = all_WVNB(i_bande)
  56 |         DWVNB   = all_DWVNB(i_bande)
  57 |         WVNB_SI = 100.*WVNB
  58 |         DWVNB_SI= 100.*DWVNB
  59 | 
  60 |         IF (WVNB.le.9300.0) THEN
  61 |           CALL FINDI(LICO,LICO2,LIH2O,ICO,ICO2,IH2O,WVNB,DWVNB)
  62 | 
  63 |           CALL KBARANDPHI(data_ref(:),LICO,LICO2,LIH2O,ICO,             &
  64 |      &                    ICO2,IH2O,KCO,KC,KH,DCO,DC,DH,                &
  65 |      &                    gasdata_ref(1:2,i_bande), nallbandes)
  66 | 
  67 |         ENDIF
  68 | 
  69 |         IF (s_Tp.eq.0) Tp = data_ref(1)
  70 |         F(i_bande) = blae(WVNB_SI, Tp)*DWVNB_SI 
  71 | 
  72 |       ENDDO
  73 | 
  74 |       F     = F / SUM (F)
  75 |       F_ref = F
  76 | 
  77 |       CALL k_distributeur(F, ngbands, gasdata_ref, nkabs, k_ref,w)
  78 | 
  79 |       DO j = 1, nkabs
  80 |         DO i_bande=1,ngbands
  81 | 
  82 |           CALL pdf(gasdata_ref(2,i_bande), gasdata_ref(1,i_bande),      &
  83 |      &             k_ref(j),fdek_i)
  84 | 
  85 |           fdek_ref(j) = fdek_ref(j) + F(i_bande)*fdek_i
  86 | 
  87 |         ENDDO
  88 |       ENDDO
  89 | 
  90 | !     PRINT*, "      DATA_REF : T    = ",data_ref(1)
  91 | !     PRINT*, "                 P    = ",data_ref(2)
  92 | !     PRINT*, "                 XH2O = ",data_ref(3)
  93 | !     PRINT*, "                 XCO2 = ", data_ref(4)
  94 | !     PRINT*, "                 XCO  = ", data_ref(5)
  95 | !     PRINT*
  96 | !     PRINT*, "K_ref    :", k_ref
  97 | !     PRINT*, "fdek_ref :", fdek_ref
  98 | 
  99 | !     ----------------------------------------------------------!
 100 | !     Calculation of the absorption coefficient over the domain !
 101 | !     ----------------------------------------------------------!
 102 | 
 103 |       all_k_abs       = 0.
 104 | 
 105 | !     ---------------------------------------------!
 106 | !     Non-homogeneous system or homogeneous system !
 107 | !     ---------------------------------------------!
 108 |       IF (homosyst.eq.'NO') THEN
 109 |         n_cl=nelmts
 110 |       ELSEIF (homosyst.eq.'YES') THEN
 111 |         n_cl=1
 112 |       ENDIF
 113 | 
 114 |       DO ielt=1,n_cl
 115 | 
 116 |         k_absmel = 0.
 117 |         FSK_SOOT = 0.
 118 |         gasdata  = 0.
 119 |         F        = 0.
 120 | 
 121 | !       ------------------------------------!
 122 | !       Setting spectral data for each band !
 123 | !       ------------------------------------!
 124 | 
 125 |         DO i_bande=1,ngbands
 126 | 
 127 |           WVNB    = all_WVNB(i_bande)
 128 |           DWVNB   = all_DWVNB(i_bande)
 129 |           WVNB_SI = 100.*WVNB
 130 |           DWVNB_SI= 100.*DWVNB
 131 | 
 132 | !         ----------------------------------------!
 133 | !         Spectral index lecture for each species !
 134 | !         ----------------------------------------!
 135 |           IF (WVNB.le.9300.0) THEN
 136 | 
 137 |             CALL FINDI(LICO,LICO2,LIH2O,ICO,ICO2,IH2O,WVNB,DWVNB)
 138 | 
 139 |             CALL KBARANDPHI(celldata(ielt,:),LICO,LICO2,LIH2O,ICO,      &
 140 |      &                      ICO2,IH2O,KCO,KC,KH,DCO,DC,DH,              &
 141 |      &                      gasdata(1:2,i_bande), nallbandes)
 142 | 
 143 |           ENDIF
 144 | 
 145 | !         -------------------------------------!
 146 | !         Band groupment multiplicative factor !
 147 | !         -------------------------------------!
 148 | 
 149 |           IF (s_Tp.eq.0) Tp = celldata(ielt,1)
 150 |           F(i_bande) = blae(WVNB_SI,Tp)*DWVNB_SI
 151 | 
 152 |           K_SOOT     = 5.5*celldata(ielt,8)*WVNB_SI
 153 |           FSK_SOOT   = FSK_SOOT + F(i_bande)*K_SOOT
 154 |  
 155 |         ENDDO
 156 | 
 157 |         FSK_SOOT = FSK_SOOT / SUM(F)
 158 |         F        = F        / SUM(F)
 159 | 
 160 | !       --------------------------------------!
 161 | !       STEP 2 : Calculate g(T,phi_ref,k_ref) !
 162 | !                and f(T,phi_ref,k_ref)       !
 163 | !       --------------------------------------!
 164 | 
 165 |         gdek = 0.
 166 |         fdek = 0.
 167 | 
 168 |         DO j = 1, nkabs
 169 |           DO i_bande=1,ngbands
 170 | 
 171 |             CALL cdss(gasdata_ref(2,i_bande), gasdata_ref(1,i_bande),   &
 172 |      &                k_ref(j), gdek_i)
 173 | 
 174 |             CALL pdf(gasdata_ref(2,i_bande), gasdata_ref(1,i_bande),    &
 175 |      &               k_ref(j), fdek_i)
 176 | 
 177 |             gdek(j) = gdek(j) + F_ref(i_bande)*gdek_i
 178 |             fdek(j) = fdek(j) + F(i_bande)*fdek_i
 179 | 
 180 |           ENDDO
 181 |         ENDDO
 182 | 
 183 | !         --------------------------------------------------------------!
 184 | !         STEP 3 : g(T_ref,phi_ref,k) = g(T_ref,phi,k*) => Correlated k !
 185 | !                  Invert g in the RHS to k* (k_absmel)                 !
 186 | !         --------------------------------------------------------------!
 187 | 
 188 |         DO j = 1, nkabs
 189 | 
 190 |           CALL COFG (F_ref, ngbands, gasdata(2,:), gasdata(1,:),gdek(j),&
 191 |      &               gdek(nkabs), k_absmel(j))
 192 | 
 193 |         ENDDO
 194 | 
 195 |         all_k_abs(ielt,:) = k_absmel + FSK_SOOT
 196 | 
 197 | !       ----------------------------------------------------------------!
 198 | !       STEP 4 : Calculation of the weight function                     !
 199 | !                a(T,T_ref,g_ref)= fdek(T,phi_ref)/fdek_(T_ref,phi_ref) !
 200 | !       ----------------------------------------------------------------!
 201 | 
 202 |         DO i = 1, nkabs
 203 |           IF(fdek_ref(i).ne.0) THEN
 204 |             WSGG_W(ielt,i) = fdek(i)/fdek_ref(i)
 205 |           ELSE
 206 |             WSGG_W(ielt,i) = 1.
 207 |           ENDIF
 208 |         ENDDO
 209 | 
 210 |         IF (MOD(ielt,100).eq.0) print*, ">> Node:",ielt,"/",n_cl
 211 | 
 212 |       ENDDO
 213 | 
 214 |       IF (homosyst.eq.'YES') THEN
 215 |         DO i=2,nelmts
 216 |           all_k_abs(i,:) = all_k_abs(1,:)
 217 |         ENDDO
 218 |       ENDIF
 219 | 
 220 |       END SUBROUTINE FSCK_CASE