The TELEMAC-MASCARET system  trunk
cvsp_make_actlay_gaia.f
Go to the documentation of this file.
1 ! ********************************
2  SUBROUTINE cvsp_make_actlay_gaia
3 ! ********************************
4 !
5 !
6 !***********************************************************************
7 ! GAIA V8P1 16/05/2017
8 !***********************************************************************
9 !
12 !
13 !
17 !
23 !
28 !
33 !
38 !
43 !
44 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
45 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
46 
47  USE bief_def
48  USE bief
50  USE cvsp_outputfiles_gaia, ONLY: cp
51  USE interface_parallel, ONLY : p_dsum
53 
55  IMPLICIT NONE
56 !
57 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
58 !
59  INTEGER I,J,K,JG, II
60  DOUBLE PRECISION TEMP, Z_HIGH, Z_LOW, AT
61  DOUBLE PRECISION ASUM, NEW_ALT
62 !
63 !-----------------------------------------------------------------------
64 !
65  DOUBLE PRECISION SUMES, SUMAV, ZSPACE
66 !
67 !-----------------------------------------------------------------------
68 ! TIME
69 !-----------------------------------------------------------------------
70  at = dt*lt/percou
72 !-----------------------------------------------------------------------
73 ! CHECKING MASS AND CORRECTING
74 !-----------------------------------------------------------------------
76 !-----------------------------------------------------------------------
77 ! NODES LOOP
78 !-----------------------------------------------------------------------
79  DO j=1,npoin
80  jg = j
81  IF (ncsize > 1) jg = mesh%KNOLG%I(j)
82 !-----------------------------------------------------------------------
83 ! CALCULATING ACTIVE LAYER THICKNESS
84 !-----------------------------------------------------------------------
85  zspace = zf%R(j)-zr%R(j)
86  nlayer%I(j) = nomblay
87  new_alt = min( zspace, cvsp_alt_gaia(j,alt_model))
88 !
89 !-----------------------------------------------------------------------
90 ! CALCULATING ES PER LAYER
91 !-----------------------------------------------------------------------
92  es(j,1) = new_alt
93  sumes = es(j,1)
94  DO k=2,nomblay-1
95  IF ((zspace-sumes).GE.elay0) THEN
96  es(j,k) = elay0
97  ELSE
98  es(j,k) = zspace-sumes
99  ENDIF
100  sumes = sumes + es(j,k)
101  END DO
102  es(j,nomblay) = zspace-sumes
103  sumes = sumes + es(j,nomblay)
104 
105 !
106  IF(abs(zspace-sumes).GT.zero) THEN
107  WRITE(lu,*)'ERR in MAKE_ACTLAY',j,(es(j,k),k=1,nomblay),
108  & sumes,zspace,zf%R(j),zr%R(j)
109  stop
110  ENDIF
111 !
112 !-----------------------------------------------------------------------
113 ! CALCULATING AVAIL PER LAYER
114 !-----------------------------------------------------------------------
115  sumes = 0.d0
116  DO k=1,nomblay
117  sumes = sumes + es(j,k)
118  IF (es(j,k).GT.0.d0) THEN
119  !DEPTH COORDINATES OF NEW ES: Z_HIGH & Z_LOW
120  IF (k == 1) THEN
121  z_high = pro_d(j,pro_max(j),1)
122  ELSE
123  z_high = -sumes + pro_d(j,pro_max(j),1) + es(j,k)
124  ENDIF
125  z_low = z_high - es(j,k)
126  !FETCH IT
127  temp = cvsp_integrate_volume_gaia(j,1,z_high,z_low,t1%R)
128  !ASSIGN IT
129  asum = 0.d0
130  DO i=1,nsicla
131  avail(j,k,i) = t1%R(i) / (z_high-z_low)
132  asum = avail(j,k,i) + asum
133  ENDDO
134 !
135 ! DEBUG
136  IF (abs(asum-1.d0).GT.1.0d-7.AND.asum.GT.0.d0.AND.cp) THEN
137  WRITE(lu,*)'ERR in MAKE_ACTLAY',es(j,1),es(j,2),
138  & es(j,3),sumes,zspace,zf%R(j),zr%R(j)
139  WRITE(lu,*)'ASUM in MAKE_ACTLAY',j,k,asum, z_high,
140  & z_low
141  DO i=1,nsicla
142  WRITE(lu,*)' AVAIL,T1',i, avail(j,k,i), t1%R(i)
143  ENDDO
144  ENDIF
145  ELSE !DUMMY VALUES FOR EMPTY LAYERS!!!
146  DO i=1,nsicla
147  avail(j,k,i) = 0.d0
148  ENDDO
149  ENDIF
150  ENDDO
151  elay%R(j) = es(j,1)
152  estrat%R(j) = es(j,2)
153  ENDDO ! J
154 !-----------------------------------------------------------------------
155 ! CHECKS, EXPENSIVE! One might try to remove this one day!
156 !-----------------------------------------------------------------------
157  DO j=1,npoin
158  zspace = zf%R(j)-zr%R(j) ! Available Space
159  sumes = 0.d0 ! Sum of layer thicknesses
160  jg = j
161  IF (ncsize > 1) jg = mesh%KNOLG%I(j)
162 ! LOOP ALL LAYERS
163  DO k=1,nlayer%I(j) ! = NOMBLAY-NUMBER OF EMPTY LAYERS
164  IF (es(j,k).LT.0.d0) THEN
165  IF (es(j,k).GT.0.d0 - zero) THEN
166  es(j,k) = zero
167  ELSE
168  WRITE(lu,*) 'Lethal_ERROR ES',j,k,i,avail(j,k,i)
169  CALL cvsp_p_gaia('./','ERR_ES_VSP_', jg)
170  CALL layers_p_gaia('./','ERR_ES_LAY_', jg)
171  CALL plante(1)
172  ENDIF
173  END IF
174  IF (k.EQ.1) elay%R(j) = es(j,k)
175  IF (k.EQ.2) estrat%R(j) = es(j,k)
176  sumes = sumes+es(j,k)
177 ! LOOP ALL FRACTIONS
178  sumav = 0.d0
179  DO i=1,nsicla
180  !STOP DUE TO WRONG FRACTION CALCULATION
181  IF((avail(j,k,i) .GT. 1.d0 + zero) .OR.
182  & (avail(j,k,i) .LT. 0.d0 - zero)) THEN
183  WRITE(lu,*) 'Lethal_ERROR AVAIL',j,k,i,avail(j,k,i)
184  CALL cvsp_p_gaia('./','ERR_AVAIL_VSP_', jg)
185  CALL layers_p_gaia('./','ERR_AVAIL_LAY_', jg)
186  CALL plante(1)
187  ENDIF
188 ! ONLY TRUNCATION ERROR
189  IF (avail(j,k,i) .LT. 0.d0) avail(j,k,i) = 0.0d0
190  IF (avail(j,k,i) .GT. 1.d0) avail(j,k,i) = 1.0d0
191  sumav=avail(j,k,i)+sumav
192  ENDDO !I = Fractions
193 !
194 ! ELIMINTATION OF TRUNCATION ERRORS
195  IF(abs(sumav-1.d0).GT.zero.AND.sumav.GT.0.d0) THEN
196  IF(cp) WRITE(lu,*) 'ERR SUMAV:J,K,SUMAV<>1:',j,k,sumav
197  DO ii=1,nsicla
198  avail(j,k,ii) = avail(j,k,ii) / sumav
199  ENDDO
200  ENDIF
201  ENDDO !K = LAYERS
202 
203 ! FINAL DEPTH CHECK
204  IF( abs(sumes-zspace).GT.zero) THEN
205  WRITE(lu,*) 'ERROR SUM ES',j,sumes,zspace,nlayer%I(j)
206  CALL cvsp_p_gaia('./','ERR_SUM_ES_VSP_', jg)
207  CALL layers_p_gaia('./','ERR_SUM_ES_LAY_', jg)
208  CALL plante(1)
209  ENDIF
210  ENDDO ! J = NODES
211 !
212 ! COMPUTES THE TOTAL VOLUME OF SEDIMENTS IN THE DOMAIN
213  DO i = 1, nsicla
214  voltot(i) = 0.d0
215  ENDDO
216 !
217  DO i=1,nsicla
218  DO j=1,npoin
219  DO k=1,pro_max(j)-1
220  voltot(i) = voltot(i) + (pro_f(j,k,i)+pro_f(j,k+1,i))/2.d0
221  & *(pro_d(j,k+1,i)-pro_d(j,k,i))*volu2d%R(j)
222  ENDDO
223  ENDDO
224  ENDDO
225 !
226  IF(ncsize.GT.1) THEN
227  DO i=1,nsicla
228  voltot(i) = p_dsum(voltot(i))
229  ENDDO
230  ENDIF
231 !
232 !-----------------------------------------------------------------------
233 !
234  RETURN
235  END SUBROUTINE cvsp_make_actlay_gaia
subroutine layers_p_gaia(PATH_PRE, FILE_PRE, JG)
Definition: layers_p_gaia.f:7
double precision elay0
Wanted active layer thickness; ELAYO is a target value for ELAY, but ELAY may be lower than ELAY0 if ...
subroutine cvsp_check_anything_gaia
type(bief_obj), target nlayer
Number of layers for each point.
integer, target nomblay
Number of bed load model layers = NUMSTRAT+1 to take the active layer into account.
type(bief_obj), target zr
Non erodable (rigid) bottom elevation.
type(bief_obj), target volu2d
Integral of bases.
double precision, target dt
Time step It may be different from the one in TELEMAC because of the morphological factor...
double precision, dimension(:,:), allocatable, target es
Layer thicknesses as double precision.
subroutine cvsp_p_gaia(PATH_PRE, FILE_PRE, JG)
Definition: cvsp_p_gaia.f:7
subroutine cvsp_check_mass_bilan_gaia
integer ncsize
Definition: bief_def.f:49
logical cp
Logical for debug printouts.
type(bief_obj), target elay
Active layer thickness.
double precision, dimension(nsiclm) voltot
Total volume of sediment of each class.
double precision zero
Parameter used for clipping variables or testing values against zero.
type(bief_obj), target estrat
Active stratum thickness.
integer alt_model
CHOOSE A MODEL FOR ESTIMATION OF A DYNAMIC ACTIVE LAYER THICKNESS.
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.
subroutine cvsp_make_actlay_gaia
double precision function p_dsum(MYPART)
Definition: p_dsum.F:7
type(bief_obj), pointer t1
Aliases for work vectors in tb.
double precision function cvsp_integrate_volume_gaia(J, I, Z_HIGH, Z_LOW, A)
integer, target lt
Numero du pas de temps.
double precision function cvsp_alt_gaia(J, FORMULA)
Definition: cvsp_alt_gaia.f:7
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_obj), target zf
Bottom elevation.
type(bief_mesh), target mesh
Mesh structure.
integer, dimension(:), allocatable pro_max
Maximum layer number in a vertical sorting profile for each point.
double precision, dimension(:,:,:), allocatable, target avail
Sediment fraction for each layer, class, point.
integer, pointer npoin
Number of 2d points in the mesh.
Definition: bief.f:3