1 | include(dom.inc)
2 |
3 | SUBROUTINE DERI_KAPPA(ivar, delta, order, dkptot)
4 | ! ================================================================!
5 | ! !
6 | ! deri_kappa.F : Calculate kappa derivative !
7 | ! !
8 | ! author : D. Poitou !
9 | ! !
10 | ! ================================================================!
11 |
12 | USE mod_slave
13 | USE mod_inout
14 |
15 | IMPLICIT NONE
16 |
17 | include 'dom_constants.h'
18 |
19 | ! IN
20 | DOM_INT :: ivar, order
21 | DOM_REAL :: delta
22 |
23 | ! LOCAL
24 | DOM_INT :: iquad, idelta, i
25 | DOM_INT :: ngbands, ios
26 | DOM_REAL :: DWVNB_SI, planck
27 | DOM_REAL,DIMENSION(order,is_nkabs,is_nnodes) :: all_kabs, allWSGG
28 | DOM_REAL,DIMENSION(is_nkabs) :: wkabs
29 | DOM_REAL,DIMENSION(is_nnodes) :: Lb, dkptot, dkabs
30 | DOM_REAL,DIMENSION(8,is_nnodes) :: data
31 | DOM_REAL,DIMENSION(order) :: delta_data
32 | DOM_REAL,DIMENSION(is_ncells) :: dkabs_cell, Lb_cell
33 | DOM_REAL, ALLOCATABLE, DIMENSION(:,:) :: data_cell
34 |
35 | all_kabs = 0.
36 | data = s_celldata
37 | DO i=1, is_nnodes
38 | Lb(i) = planck(s_celldata(1,i))
39 | ENDDO
40 |
41 | IF(order.eq. 2) THEN
42 | delta_data(1) = -delta
43 | delta_data(2) = delta
44 | ENDIF
45 |
46 | DO idelta = 1, order
47 |
48 | data(ivar, :) = s_celldata(ivar,:) + delta_data(idelta)
49 |
50 | DWVNB_SI = 1.0
51 | ngbands = is_nallbandes
52 |
53 | ! ---------------------!
54 | ! WSGG case treatement !
55 | ! ---------------------!
56 | IF (mediumtype.eq.'WSGG') THEN
57 |
58 | wkabs = 1.0
59 |
60 | CALL WSGG_CASE(all_kabs(idelta,:,:),s_WSGG_W,Lb,data, &
61 | & is_nnodes, s_alpha,s_kwsgg,s_all_WVNB, &
62 | & s_all_DWVNB, is_ngg,is_nallbandes, 1, &
63 | & is_nnodes)
64 |
65 | ! --------------------!
66 | ! FSK CASE TREATEMENT !
67 | ! --------------------!
68 | ELSEIF (mediumtype.eq.'FSK') THEN
69 |
70 | s_WSGG_W = 1.0
71 |
72 | CALL FSK_CASE(all_kabs(idelta,:,:), wkabs, data,is_nnodes, &
73 | & s_all_WVNB, s_all_DWVNB, s_KCO, s_DCO, s_KC, &
74 | & s_DC, s_KH, s_DH, ngbands, is_nallbandes, &
75 | & is_nkabs, is_lbcd, is_lbcf)
76 |
77 | ! ---------------------!
78 | ! FSCK CASE TREATEMENT !
79 | ! ---------------------!
80 | ELSEIF (mediumtype.eq.'FSCK') THEN
81 |
82 | IF(.not.ALLOCATED(data_cell)) &
83 | & ALLOCATE(data_cell(8,is_ncells))
84 | CALL GATHER(s_celldata, data_cell, 8)
85 |
86 | ! -------------------------------!
87 | ! Choice of the reference state !
88 | ! 0 -> from input_fsck.dat !
89 | ! 1 -> average state !
90 | ! 2 -> Maximal temperature !
91 | ! 3 -> Maximal absorbing species !
92 | ! -------------------------------!
93 |
94 | IF (i_refstate.eq.0) THEN
95 | OPEN(UNIT=1,FILE='input_fsck.dat',FORM='FORMATTED', &
96 | & STATUS='old',IOSTAT=ios)
97 | DO i=1,8
98 | READ(1,*,IOSTAT=ios) s_dataref(i)
99 | ENDDO
100 | IF (ios.ne.0) THEN
101 | PRINT*, "Error reading input_fsck.dat"
102 | STOP
103 | ENDIF
104 |
105 | ELSEIF (i_refstate.eq.1) THEN
106 | CALL VOLAVERAGE(data_cell,s_dataref,8,1)
107 | DEALLOCATE(data_cell)
108 | ELSEIF (i_refstate.eq.2) THEN
109 | i = MAXVAL(MAXLOC(s_celldata(1,:)))
110 | s_dataref(:) = s_celldata(:,i)
111 | ELSEIF (i_refstate.eq.3) THEN
112 | i = MAXVAL(MAXLOC(s_celldata(3:5,:)))
113 | s_dataref(:) = s_celldata(:,i)
114 | ENDIF
115 |
116 | IF (.not.ALLOCATED(s_WSGG_Wb)) THEN
117 | ALLOCATE(s_WSGG_Wb(is_nkabs,is_nbfaces))
118 | ENDIF
119 |
120 | CALL FSCK_CASE(all_kabs(idelta,:,:), wkabs, data,is_nnodes, &
121 | & s_all_WVNB, s_all_DWVNB, s_KCO, s_DCO, s_KC, &
122 | & s_DC, s_KH, s_DH, ngbands, is_nallbandes, &
123 | & is_nkabs, s_dataref, s_WSGG_W, &
124 | & is_lbcd, is_lbcf, s_WSGG_Wb, ts_boundary%Tf, &
125 | & is_nbfaces)
126 |
127 | ! ---------------------------------------!
128 | ! Tabulated FSCK and FSK CASE TREATEMENT !
129 | ! ---------------------------------------!
130 |
131 | ELSEIF ((mediumtype.eq.'TFSCK').or.(mediumtype.eq.'TFSK')) THEN
132 |
133 | IF(mediumtype.eq.'TFSCK') THEN
134 |
135 | IF (.not.ALLOCATED(s_WSGG_Wb)) THEN
136 | ALLOCATE(s_WSGG_Wb(is_nkabs,is_nbfaces))
137 | ENDIF
138 |
139 | CALL TFSCK_CASE(data ,is_nnodes, &
140 | & s_all_WVNB, s_all_DWVNB, s_KCO, s_DCO, s_KC, &
141 | & s_DC, s_KH, s_DH, ngbands, is_nallbandes, &
142 | & is_nkabs, s_dataref, s_WSGG_W, &
143 | & 1, is_nnodes, s_WSGG_Wb, ts_boundary%Tf, &
144 | & is_nbfaces)
145 |
146 | ELSE
147 |
148 | s_WSGG_W = 1.0
149 |
150 | ENDIF
151 |
152 | CALL TFSK_CASE(data, all_kabs(idelta,:,:), wkabs, s_tabkabs,&
153 | & s_DYH, s_DYC, s_DYCO, is_nYH, is_nYC, is_nYCO,&
154 | & s_DT, is_nT, is_nkabs, &
155 | & is_nnodes, is_nallbandes, s_all_WVNB, &
156 | & s_all_DWVNB, 1, is_nnodes)
157 |
158 | ELSE
159 | PRINT*," Unknown spectral type to derive kappa",mediumtype
160 | ENDIF
161 |
162 | allWSGG(idelta,:,:) = s_WSGG_W(:,:)
163 |
164 | ENDDO
165 |
166 | dkptot = 0
167 |
168 | call gather(Lb, Lb_cell, 1)
169 | call scatter(Lb_cell, Lb,1)
170 |
171 | DO iquad = 1, is_nkabs
172 |
173 | IF(delta.eq.0.or.order.lt.2) THEN
174 | dkabs(:) = all_kabs(1,iquad,:)*allWSGG(1,iquad,:)
175 | ELSEIF (order.eq.2) THEN
176 | dkabs(:) = (all_kabs(2,iquad,:) &
177 | & - all_kabs(1,iquad,:))/delta
178 | ENDIF
179 |
180 | ! With Gather/scatter
181 | call gather(dkabs*(allWSGG(1,iquad,:)+allWSGG(2,iquad,:))/2, &
182 | & dkabs_cell,1)
183 | call scatter(dkabs_cell, dkabs, 1)
184 |
185 | dkptot(:) = dkptot(:) + Lb(:)*4*pi*dkabs(:)*wkabs(iquad) &
186 | & *DWVNB_SI
187 |
188 | ! Without Gather/Scatter
189 | ! dkptot(:) = dkptot(:) + Lb(:)*4*pi*dkabs(:)*wkabs(iquad) &
190 | ! & *DWVNB_SI*s_WSGG_W(iquad,:)
191 |
192 | ENDDO
193 |
194 | dkptot(:) = dkptot(:) / (4*pi*sigma*s_celldata(1,:)**4)
195 |
196 | ENDSUBROUTINE DERI_KAPPA
deri_kappa.F could be called by: