The TELEMAC-MASCARET system  trunk
cvsp_integrate_volume_gaia.f
Go to the documentation of this file.
1 ! ****************************************************
2  DOUBLE PRECISION FUNCTION cvsp_integrate_volume_gaia
3 ! ****************************************************
4 !
5  &(j,i,z_high,z_low,a)
6 !
7 !***********************************************************************
8 ! GAIA V8P1 16/05/2017
9 !***********************************************************************
10 !
13 !
17 !
22 !
27 !
32 !
39 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
45 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
46 !
47  USE interface_gaia,
48  & ex_cvsp_integrate_volume_gaia => cvsp_integrate_volume_gaia
50  USE cvsp_outputfiles_gaia, ONLY: cp
51 !
53  IMPLICIT NONE
54 !
55 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
56 !
57  INTEGER, INTENT(IN) :: J
58  INTEGER, INTENT(IN) :: I
59  DOUBLE PRECISION, INTENT(IN) :: Z_HIGH
60  DOUBLE PRECISION, INTENT(IN) :: Z_LOW
61  DOUBLE PRECISION, INTENT(OUT) :: A(nsicla)
62 !
63 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
64 !
65  DOUBLE PRECISION TEMP,TEMPOLD, DHIG, DLOW, AT, SUMUP, FLOW, FHIG
66  DOUBLE PRECISION TEMP2, TEMP2MAX, SUMUP2,TEMP3, TEMP3MAX, SUMUP3
67  DOUBLE PRECISION CORRECT, CHSUM
68  INTEGER L_CNT, F_CNT, REVCNT, HELPER, LASTCASE, JG, K
69  INTEGER(KIND=K8) :: MYCASE
70  LOGICAL RET
71 !
72 !-----------------------------------------------------------------------
73 !
74  jg = j
75  IF (ncsize.GT.1) jg = mesh%KNOLG%I(j)
76 !
77  at = dt*lt/percou
78  sumup = 0.d0
79  sumup2 = 0.d0
80  sumup3 = 0.d0
81  mycase = 0
82 !
83  DO f_cnt = 1, nsicla
84  a(f_cnt) = 0.d0
85  END DO
87  IF ((z_high-z_low).LT.zero) THEN
88  RETURN
89  ENDIF
90 !-----------------------------------------------------------------------
91 ! DOING ALL FRACTIONS
92 !-----------------------------------------------------------------------
93  DO f_cnt = 1, nsicla
94  temp = 0.d0
95  temp2 = 0.d0
96  temp2max = 0.d0
97  temp3 = 0.d0
98  temp3max = 0.d0
99  chsum = 0.d0
100 !-----------------------------------------------------------------------
101 ! GOING THROUGH ALL SECTIONS
102 !-----------------------------------------------------------------------
103  DO l_cnt = 0,(pro_max(j)-2)
104  tempold = temp
105  revcnt = pro_max(j)-l_cnt
106 !-----------------------------------------------------------------------
107 ! DEPTH COORDINATES OF THE SECTION TO CHECK
108 !-----------------------------------------------------------------------
109  dhig = pro_d(j,revcnt,f_cnt)
110  dlow = pro_d(j,revcnt-1,f_cnt)
111 !-----------------------------------------------------------------------
112 ! SORTING PROFILE SECTION TOTALLY INSIDE (CASE ZDDZ)
113 !-----------------------------------------------------------------------
114  IF ( (dhig <= z_high ) .AND.
115  & (dlow >= z_low ) ) THEN
116  flow = pro_f(j,revcnt-1,f_cnt)
117  fhig = pro_f(j,revcnt,f_cnt)
118  mycase = mycase + 1
119  lastcase = 1
120  temp2 = 0.5d0*(fhig+flow)*(dhig-dlow) + temp2
121  temp2max = z_high - dlow
122  DO helper = 1, nsicla
123  chsum = pro_f(j,revcnt, helper) + chsum
124  ENDDO
125  chsum = 1.d0 - chsum
126 !-----------------------------------------------------------------------
127 ! SORTING PROFILE SECTION PARTIALLY LOWER (CASE ZDZD)
128 !-----------------------------------------------------------------------
129  ELSEIF ((dhig <= z_high) .AND.
130  & (dhig > z_low ) .AND.
131  & (dlow < z_low ) ) THEN
132 
133  fhig = pro_f(j,revcnt,f_cnt)
134  flow = pro_f(j,revcnt-1,f_cnt)
135  flow = - ((fhig-flow)/(dhig-dlow))*(dhig-z_low) + fhig
136 !-----------------------------------------------------------------------
137 ! CUT THE SECTION
138 !-----------------------------------------------------------------------
139  dlow = z_low
140  mycase = mycase + 1000
141  lastcase = 2
142  temp3 = 0.5d0*(fhig+flow)*(dhig-dlow) + temp3
143  temp3max = dhig - dlow
144 !-----------------------------------------------------------------------
145 ! SORTING PROFILE SECTION PARTIALLY HIGHER (CASE DZDZ)
146 !-----------------------------------------------------------------------
147 !
148  ELSEIF ((dhig > z_high) .AND.
149  & (dlow >= z_low) .AND.
150  & (dlow < z_high) ) THEN
151  flow = pro_f(j,revcnt-1,f_cnt)
152  fhig = pro_f(j,revcnt,f_cnt)
153  fhig = ((fhig-flow)/(dhig-dlow))*(z_high-dlow) + flow
154 !-----------------------------------------------------------------------
155 !INSERT
156 !-----------------------------------------------------------------------
157 !CUT THE SECTION
158  dhig = z_high
159  mycase = mycase + 1000000
160  lastcase = 3
161 !-----------------------------------------------------------------------
162 ! LAYER TOTALLY INSIDE ONE SECTION (CASE DZZD)
163 !-----------------------------------------------------------------------
164  ELSEIF ((dhig >= z_high) .AND.
165  & (dlow <= z_low) ) THEN
166  fhig =
167  & - (pro_f(j,revcnt,f_cnt)-pro_f(j,revcnt-1,f_cnt)) /
168  & (dhig-dlow) *
169  & (dhig-z_high)
170  & + pro_f(j,revcnt,f_cnt)
171  flow =
172  & - (pro_f(j,revcnt,f_cnt)-pro_f(j,revcnt-1,f_cnt)) /
173  & (dhig-dlow) *
174  & (dhig-z_low)
175  & + pro_f(j,revcnt,f_cnt)
176 !-----------------------------------------------------------------------
177 ! CUT THE SECTION
178 !-----------------------------------------------------------------------
179  dhig = z_high
180  dlow = z_low
181  lastcase = 8
182  mycase = mycase + 100000000
183 !-----------------------------------------------------------------------
184 ! SECTION WITH 0 STRENGTH
185 !-----------------------------------------------------------------------
186  ELSEIF (dhig == dlow) THEN
187  flow = 0.d0
188  fhig = 0.d0
189  mycase = mycase + 1000000000
190  lastcase = 4
191 !A TRUE BUG!
192  ELSEIF (z_low.GE.z_high) THEN
193  flow = 0.d0
194  fhig = 0.d0
195  IF(cp)WRITE(lu,*)'Z_LOW>=Z_HIGH',dhig,dlow,z_high,z_low
196  CALL cvsp_p_gaia('./','ZLOHI',jg)
197  mycase = mycase + 1000000000
198  lastcase = 5
199 !A TRUE BUG!
200  ELSEIF (dhig < dlow) THEN
201  flow = 0.d0
202  fhig = 0.d0
203 
204  IF(cp)WRITE(lu,*)'DHIG<=DLOW',j,dhig,dlow,z_high,z_low
205  CALL cvsp_p_gaia('./','DLOHI',jg)
206  mycase = mycase + 1000000000
207  lastcase = 6
208 !A SECTION THAT IS NOT INVOLVED / NOT A BUG!!
209  ELSE
210  flow = 0.d0
211  fhig = 0.d0
212  mycase = mycase + 1000000000
213  lastcase = 7
214  ENDIF
215 !-----------------------------------------------------------------------
216 ! TRAPEZOID FORMULA
217 !-----------------------------------------------------------------------
218  temp = 0.5d0*(fhig+flow)*(dhig-dlow) + temp
219 !DEBUG
220  IF (0.5d0*(fhig+flow)*(dhig-dlow) < 0.d0) THEN
221  WRITE(lu,fmt='(A,1X,2(I11,1X),11(G20.10,1X),1X,I11)')
222  & 'INTEGRATE_VOL_ER_TMP:<0:'
223  & ,jg, i, at, fhig,flow,dhig,dlow, dhig-dlow, revcnt,
224  & pro_f(j,revcnt-1,f_cnt),pro_f(j,revcnt,f_cnt),
225  & pro_d(j,revcnt-1,f_cnt),pro_d(j,revcnt,f_cnt),
226  & lastcase
227  CALL cvsp_p_gaia('./','IVKT',jg)
228  CALL plante(1)
229  ENDIF
230 !
231  ENDDO !SECTION
232 !-----------------------------------------------------------------------
233 ! ADDING UP FRACTIONS FOR DEBUGGING PURPOSES
234 !-----------------------------------------------------------------------
235  sumup = temp + sumup
236  sumup2 = temp2 + sumup2
237  sumup3 = temp3 + sumup3
238  a(f_cnt) = temp
239  ENDDO !FRACTION
241 !
242 ! DEBUG
243  IF (cvsp_integrate_volume_gaia < 0.d0) THEN
244  CALL cvsp_p_gaia('./','IVK0',jg)
245  WRITE(*,fmt='(A,2(I11),14(G20.10))')'INTEGRATE_VOL_ER:<0:'
246  & ,jg, i, at, cvsp_integrate_volume_gaia,
247  & mycase ,z_high,z_low,
248  & sumup,sumup2,sumup3,chsum,a(1),a(2),a(3),a(4),a(5)
249  CALL plante(1)
250  ENDIF
251 !-----------------------------------------------------------------------
252 ! CHECKSUM OVER ALL FRACTIONS AND LAYERS
253 !-----------------------------------------------------------------------
254  sumup = sumup - abs(z_high-z_low)
255  sumup2 = sumup2 - abs(temp2max)
256  sumup3 = sumup3 - abs(temp3max)
257  chsum = chsum / 5.d0
258 !-----------------------------------------------------------------------
259 ! BEWARE: 1.D-5 IS UP TO 10G OF A 1000KG BUT HIGHER ACCURACY
260 ! LEADS AGAIN TO FLOATING POINT TRUNCATION ERRORS ...
261 !-----------------------------------------------------------------------
262  IF (abs(sumup).GT.1.d-5) THEN
263  IF(cp) CALL cvsp_p_gaia('./','IV_E',jg)
264  IF(cp)WRITE(lu,*) 'INTEGRATE VOLUME ACCURRACY!!!',sumup,jg
265  ENDIF
266 !
267 !-----------------------------------------------------------------------
268 !
269  RETURN
270  END FUNCTION cvsp_integrate_volume_gaia
double precision, target dt
Time step It may be different from the one in TELEMAC because of the morphological factor...
subroutine cvsp_p_gaia(PATH_PRE, FILE_PRE, JG)
Definition: cvsp_p_gaia.f:7
logical cp
Logical for debug printouts.
double precision zero
Parameter used for clipping variables or testing values against zero.
integer, target nsicla
Number of sediment classes of bed material (less than NISCLM)
double precision, dimension(:,:,:), allocatable, target pro_f
Vertical sorting profile: fraction for each layer, class, point.
double precision function cvsp_integrate_volume_gaia(J, I, Z_HIGH, Z_LOW, A)
integer, target lt
Numero du pas de temps.
double precision, dimension(:,:,:), allocatable, target pro_d
Vertical sorting profile: depth for each layer, class, point.
integer percou
COUPLING PERIOD USED IN CVSM TO CALCULATE THE TIME, SHOULD COME FROM TELEMAC ONE DAY.
type(bief_mesh), target mesh
Mesh structure.
integer, dimension(:), allocatable pro_max
Maximum layer number in a vertical sorting profile for each point.