The TELEMAC-MASCARET system  trunk
calcs3d_o2.f
Go to the documentation of this file.
1 ! *********************
2  SUBROUTINE calcs3d_o2
3 ! *********************
4 !
5  & (npoin3,npoin2,nplan,wattemp,o2satu,demben,formk2,k1,k44,k22,
6  & photo,resp,tn,texp,timp,t31,t32,t21,hprop,zprop,un,vn)
7 !
8 !***********************************************************************
9 ! WAQTEL V8P1
10 !***********************************************************************
11 !
12 !brief COMPUTES SOURCE TERMS FOR O2 WAQ PROCESS FOR 3D CASE
13 !
14 !history R. ATA
15 !+ 21/02/2016
16 !+ V7P2
17 !+ CREATION
18 !
19 !history S.E. BOURBAN (HRW)
20 !+ 07/06/2017
21 !+ V7P3
22 !+ Indexing tracer (IND_*) to avoid conflicting naming convention
23 !+ between user defined tracers, water quality processes and
24 !+ ice processes. Introduction of the array RANK_*.
25 !
26 !history S.E. BOURBAN (HRW)
27 !+ 25/09/2017
28 !+ V7P3
29 !+ TEXP and TIMP are now additive to account for a variety of
30 !+ of sources / sinks on a given TRACER
31 !
32 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
33 !| DEMBEN |-->| BENTIC DEMAND
34 !| DT |-->| TIME STEP
35 !| FORMK2 |-->| FORMULA FOR COMPUTINK K2
36 !| HPROP |-->| PROPAGATION DEPTH
37 !| K1 |-->| CONST. OF DEGRADATION KINETIC OF ORGANIC LOAD
38 !| K22 |<--| CONST. OF REAERATION
39 !| K44 |-->| CONST. OF NITRIFICATION KINETIC
40 !| MASSOU |<--| MASS OF TRACER ADDED BY SOURCE TERM
41 !| NPOIN |-->| NUMBER OF NODES IN THE MESH
42 !| O2SATU |<--| O2 SATURATION DENSITY OF WATER (CS)
43 !| PHOTO |-->| PHOTOSYNTHESIS
44 !| WATTEMP |-->| TEMPERATURE OF WATER (DEG CELSIUS)
45 !| T1,T2 |-->| WORKING ARRAYS
46 !| TEXP |<--| EXPLICIT SOURCE TERM.
47 !| TIMP |<--| IMPLICIT SOURCE TERM.
48 !| TN |-->| TRACERS AT TIME N
49 !| UN,VN |-->| VELOCITY COMPONENTS AT TIME N
50 !| ZPROP |-->| PROPAGATION ELEVATION
51 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
52 !
53 
54  USE bief
57  & ind_t,ind_o2,ind_ol,ind_nh4
59  USE interface_waqtel, ex_calcs3d_o2 => calcs3d_o2
60 !
61  IMPLICIT NONE
62 !
63 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
64 !
65 !
66  INTEGER , INTENT(IN ) :: npoin2,npoin3,nplan
67  INTEGER , INTENT(IN ) :: formk2
68  DOUBLE PRECISION , INTENT(IN ) :: demben,wattemp
69  DOUBLE PRECISION , INTENT(IN ) :: photo,resp,k1,k44
70  DOUBLE PRECISION , INTENT(INOUT) :: o2satu,k22
71  TYPE(bief_obj) , INTENT(IN ) :: tn,hprop,un,vn,zprop
72  TYPE(bief_obj) , INTENT(INOUT) :: texp,timp,t31,t32,t21
73 !
74 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
75 !
76 ! LOCAL VARIABLES
77  INTEGER :: i,j
78  DOUBLE PRECISION, PARAMETER :: eps=1.d-3
79  DOUBLE PRECISION, PARAMETER :: corr1=1.065d0
80  DOUBLE PRECISION, PARAMETER :: corr2=1.0241d0
81  DOUBLE PRECISION :: bencorr,pmoinr,conv_k1,conv_k44
82  DOUBLE PRECISION :: power
83  INTRINSIC max
84 !
85  pmoinr = (photo-resp)*sectoday
86 ! CONVERT DAYS TO SECONDS
87  conv_k44= k44 *sectoday
88  conv_k1 = k1 *sectoday
89  bencorr = demben*sectoday
90  power = wattemp-20.d0
91  DO i=1,npoin3
92  IF( ind_t.GT.0 ) power=tn%ADR(ind_t)%P%R(i)-20.d0
93 ! CORR2T AND BENCOR STOCKED HERE IN T31,T32
94  t31%R(i)=corr2**power
95  t32%R(i)=bencorr*(corr1**power)
96  ENDDO
97 !
98 ! COMPUTE CS (O2SATU, STOCKED IN T21)
99 !
100  IF( ind_t.EQ.0 )THEN
101  CALL satur_o2(o2satu,formcs,wattemp,eps)
102  CALL os('X=C ',x=t21,c=o2satu )
103  ELSE
104  DO i=1,npoin2
105  j=(nplan-1)*npoin2+i
106  CALL satur_o2(t21%R(i),formcs,tn%ADR(ind_t)%P%R(j),eps)
107  ENDDO
108  ENDIF
109 !
110 ! COMPUTE K2
111 !
112  CALL reaer(formk2,k2,k22,npoin2,nplan,un,vn,hprop,eps)
113 ! CONVERT DAY TO SEC
114  CALL os('X=CX ',x=k2,c=sectoday)
115 !
116 ! COMPUTE RS (DONE IN DIFSOU)
117 !
118 !----------------------------------------------------------------------
119 ! LET'S NOW COMPUTE SOURCE TERMS
120 !
121 ! FIRST TRACER O2 (RANK IND_O2)
122 ! warning: here there are lots of choices that can be changed
123 ! by the user:
124 ! 1- reareration is considered for only water surface,
125 ! 2- bentic demand concerns only the bottom or all the water
126 ! column (here the second choice is considered)
127 ! these choices and many others, have to be decided by WAQ
128 ! experts
129 !
130 ! EXPLICIT PART
131 ! =============
132 !
133 ! SURFACE SOURCES
134  DO i=1,npoin2
135  j=(nplan-1)*npoin2+i
136  texp%ADR(ind_o2)%P%R(j) = texp%ADR(ind_o2)%P%R(j) +
137  & t31%R(j)*k2%R(i)*max((t21%R(i)-tn%ADR(ind_o2)%P%R(j)),0.d0) -
138 ! O2 CONCENTRATION CAN NOT BE GREATER THAN SATURATED
139  & t32%R(i)/max(eps,zprop%R(j))
140  ENDDO
141 !
142 ! ADD PHOTOSYNTHESIS - RESPIRATION
143 !
144  CALL os('X=X+C ',x=texp%ADR(ind_o2)%P,c=pmoinr)
145 !
146 ! THE REMAINING TERMS
147 !
148  CALL os('X=X+CY ',x=texp%ADR(ind_o2)%P,y=tn%ADR(ind_ol)%P,
149  & c=-conv_k1 )
150  CALL os('X=X+CY ',x=texp%ADR(ind_o2)%P,y=tn%ADR(ind_nh4)%P,
151  & c=-conv_k44)
152 ! IMPLICIT PART:
153 ! ==============
154 ! SOFAR, ALL TERMS ARE EXPLICIT FOR O2. CHANGE
155 ! IF THERE ARE DISAPPOINTING RESULTS
156 !
157 ! SECOND TRACER [L] ORGANIC LOAD
158 !
159 ! CALL OS('X=C ',X=TIMP%ADR(IND_OL)%P,C=CONV_K1)
160  CALL os('X=X+C ',x=timp%ADR(ind_ol)%P,c=conv_k1)
161 !
162 ! THIRD TRACER [NH4]
163 !
164 ! CALL OS('X=C ',X=TIMP%ADR(IND_NH4)%P,C=CONV_K44)
165  CALL os('X=X+C ',x=timp%ADR(ind_nh4)%P,c=conv_k44)
166 !
167 ! MASS BALANCE: MASS ADDED BY EXPLICIT TERMS (TO CHECK)
168 !
169 !
170 !-----------------------------------------------------------------------
171 !
172  RETURN
173  END
double precision, parameter sectoday
type(bief_obj), target k2
subroutine reaer(FORMK2, K2, K22, NPOIN2, NPLAN, UN, VN, H, EPS)
Definition: reaer.f:7
subroutine os(OP, X, Y, Z, C, IOPT, INFINI, ZERO)
Definition: os.f:7
subroutine satur_o2(SATO2, FORMCS, WATTEMP, EPS)
Definition: satur_o2.f:7
Definition: bief.f:3