The TELEMAC-MASCARET system  trunk
calcs2d_o2.f
Go to the documentation of this file.
1 ! *********************
2  SUBROUTINE calcs2d_o2
3 ! *********************
4 !
5  & (npoin,wattemp,o2satu,demben,formk2,k1,k44,k22,
6  & photo,resp,tn,texp,timp,t2,t3,t4,hprop,un,vn,debug)
7 !
8 !***********************************************************************
9 ! WAQTEL V8P1
10 !***********************************************************************
11 !
12 !brief COMPUTES SOURCE TERMS FOR O2 WAQ PROCESS
13 !
14 !history R. ATA
15 !+ 21/09/2014
16 !+ V7P0
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 !| DEBUG |-->| FOR DEBUGGING
34 !| DEMBEN |-->| BENTIC DEMAND
35 !| FORMK2 |-->| FORMULA FOR COMPUTINK K2
36 !| K1 |-->| CONST. OF DEGRADATION KINETIC OF ORGANIC LOAD
37 !| HPROP |-->| PROPAGATION DEPTH
38 !| K22 |<--| CONST. OF REAERATION
39 !| K44 |-->| CONST. OF NITRIFICATION KINETIC
40 !| NPOIN |-->| NUMBER OF NODES IN THE MESH
41 !| O2SATU |<--| O2 SATURATION DENSITY OF WATER (CS)
42 !| PHOTO |-->| PHOTOSYNTHESIS
43 !| WATTEMP |-->| TEMPERATURE OF WATER (DEG CELSIUS)
44 !| T2,T3,T4 |-->| WORKING ARRAY
45 !| TEXP |<--| EXPLICIT SOURCE TERM.
46 !| TIMP |<--| IMPLICIT SOURCE TERM.
47 !| TN |-->| TRACERS AT TIME N
48 !| UN,VN |-->| VELOCITY COMPONENTS AT TIME N
49 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
50 !
51 
52  USE bief
55  & ind_t,ind_o2,ind_ol,ind_nh4
57  USE interface_waqtel, ex_calcs2d_o2 => calcs2d_o2
58 !
59  IMPLICIT NONE
60 !
61 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
62 !
63 !
64  INTEGER , INTENT(IN ) :: formk2,npoin,debug
65  DOUBLE PRECISION , INTENT(IN ) :: demben,wattemp
66  DOUBLE PRECISION , INTENT(IN ) :: photo,resp,k1,k44
67  DOUBLE PRECISION , INTENT(INOUT) :: o2satu,k22
68  TYPE(bief_obj) , INTENT(IN ) :: tn,hprop,un,vn
69  TYPE(bief_obj) , INTENT(INOUT) :: texp,timp,t2,t3,t4
70 !
71 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
72 !
73 ! LOCAL VARIABLES
74  INTEGER :: i
75  DOUBLE PRECISION, PARAMETER :: eps=1.d-3
76  DOUBLE PRECISION, PARAMETER :: corr1=1.065d0
77  DOUBLE PRECISION, PARAMETER :: corr2=1.0241d0
78  DOUBLE PRECISION :: bencorr,pmoinr,conv_k1,conv_k44
79  DOUBLE PRECISION :: power
80  INTRINSIC max
81 !
82 !-----------------------------------------------------------------------
83 !
84 !
85 ! PRELIMINARY COMPUTATIONS
86 !
87  pmoinr = (photo-resp)*sectoday
88 ! CONVERT DAYS TO SECONDS
89  conv_k44= k44 *sectoday
90  conv_k1 = k1 *sectoday
91  bencorr = demben*sectoday
92  power = wattemp-20.d0
93  DO i=1,npoin
94  IF( ind_t.GT.0 ) power=tn%ADR(ind_t)%P%R(i)-20.d0
95 ! CORR2T AND BENCOR STOCKED HERE IN T3,T4
96  t3%R(i)=corr2**power
97  t4%R(i)=bencorr*(corr1**power)
98  ENDDO
99 !
100 ! COMPUTE CS (O2SATU, STOCKED IN T2)
101 !
102  IF( ind_t.EQ.0 )THEN
103  CALL satur_o2(o2satu,formcs,wattemp,eps)
104  CALL os('X=C ',x=t2,c=o2satu )
105  ELSE
106  DO i=1,npoin
107  CALL satur_o2(t2%R(i),formcs,tn%ADR(ind_t)%P%R(i),eps)
108  ENDDO
109  ENDIF
110 !
111 ! COMPUTE K2
112 !
113  CALL reaer(formk2,k2,k22,npoin,1,un,vn,hprop,eps)
114 ! CONVERT DAY TO SEC
115  CALL os('X=CX ',x=k2,c=sectoday)
116 !
117 ! COMPUTE RS (DONE IN DIFSOU)
118 !
119 !----------------------------------------------------------------------
120 ! LET'S NOW COMPUTE SOURCE TERMS
121 !
122 ! FIRST TRACER O2 (RANK IND_O2)
123 ! warning: here there are lots of choices that can be changed
124 ! by the user:
125 ! 1- reareration is considered for only water surface,
126 ! 2- bentic demand concerns only the bottom or all the water
127 ! column (here the second choice is considered)
128 ! these choices and many others, have to be decided by WAQ
129 ! experts
130 !
131 ! EXPLICIT PART
132 ! =============
133 !
134 ! SURFACE SOURCES
135  DO i=1,npoin
136  texp%ADR(ind_o2)%P%R(i) = texp%ADR(ind_o2)%P%R(i)
137  & + t3%R(i)*k2%R(i) * max((t2%R(i)-tn%ADR(ind_o2)%P%R(i)),0.d0)
138 ! O2 CONCENTRATION CAN NOT BE GREATER THAN SATURATED
139  & - t4%R(i)/max(eps,hprop%R(i))
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  IF(debug.GT.0)WRITE(lu,*)'IN O2, STEP1:',texp%ADR(ind_o2)%P%R(1)
153 
154 ! IMPLICIT PART:
155 ! ==============
156 ! SOFAR, ALL TERMS ARE EXPLICIT FOR O2. CHANGE
157 ! IF THERE ARE DISAPPOINTING RESULTS
158 !
159 ! SECOND TRACER [L] ORGANIC LOAD
160 !
161 ! /!\ contrary to 3D formulation
162  CALL os('X=X+CY ',x=timp%ADR(ind_ol)%P,y=hprop,c=-conv_k1)
163 ! CALL OS('X=X+C ',X=TIMP%ADR(IND_OL)%P,C=CONV_K1)
164  IF(debug.GT.0)WRITE(lu,*)'IN O2, STEP2:',timp%ADR(ind_ol)%P%R(1)
165 !
166 ! THIRD TRACER [NH4]
167 !
168 ! /!\ contrary to 3D formulation
169  CALL os('X=X+CY ',x=timp%ADR(ind_nh4)%P,y=hprop,c=-conv_k44)
170 ! CALL OS('X=X+C ',X=TIMP%ADR(IND_NH4)%P,C=CONV_K44)
171  IF(debug.GT.0)WRITE(lu,*)'IN O2, STEP3:',timp%ADR(ind_nh4)%P%R(1)
172 !
173 !-----------------------------------------------------------------------
174 !
175  RETURN
176  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