1 | MODULE sedco3 |
---|
2 | #if defined key_sed |
---|
3 | !!====================================================================== |
---|
4 | !! *** MODULE sedco3 *** |
---|
5 | !! Sediment : carbonate in sediment pore water |
---|
6 | !!===================================================================== |
---|
7 | !! * Modules used |
---|
8 | USE sed ! sediment global variable |
---|
9 | |
---|
10 | |
---|
11 | IMPLICIT NONE |
---|
12 | PRIVATE |
---|
13 | |
---|
14 | !! * Routine accessibility |
---|
15 | PUBLIC sed_co3 |
---|
16 | |
---|
17 | |
---|
18 | !! * Module variables |
---|
19 | REAL(wp) :: & !! |
---|
20 | epsmax = 1.e-12 ! convergence limite value |
---|
21 | |
---|
22 | !!---------------------------------------------------------------------- |
---|
23 | !! OPA 9.0 ! LODYC-IPSL (2003) |
---|
24 | !!---------------------------------------------------------------------- |
---|
25 | |
---|
26 | CONTAINS |
---|
27 | |
---|
28 | |
---|
29 | SUBROUTINE sed_co3( kt ) |
---|
30 | !!---------------------------------------------------------------------- |
---|
31 | !! *** ROUTINE sed_co3 *** |
---|
32 | !! |
---|
33 | !! ** Purpose : carbonate ion and proton concentration |
---|
34 | !! in sediment pore water |
---|
35 | !! |
---|
36 | !! ** Methode : - solving nonlinear equation for [H+] with given alkalinity |
---|
37 | !! and total co2 |
---|
38 | !! - one dimensional newton-raphson algorithm for [H+]) |
---|
39 | !! |
---|
40 | !! History : |
---|
41 | !! ! 98-08 (E. Maier-Reimer, Christoph Heinze ) Original code |
---|
42 | !! ! 04-10 (N. Emprin, M. Gehlen ) coupled with PISCES |
---|
43 | !! ! 06-04 (C. Ethe) Re-organization |
---|
44 | !!---------------------------------------------------------------------- |
---|
45 | !! * Arguments |
---|
46 | INTEGER, INTENT(in) :: & |
---|
47 | kt ! time step |
---|
48 | |
---|
49 | ! |
---|
50 | !---Local variables |
---|
51 | INTEGER :: & |
---|
52 | jiter, ji, jk, ipt ! dummy loop indices |
---|
53 | |
---|
54 | INTEGER :: & |
---|
55 | itermax , & ! maximum number of Newton-Raphson iterations |
---|
56 | itime ! number of time to perform Newton-Raphson algorithm |
---|
57 | |
---|
58 | LOGICAL :: & |
---|
59 | lconv = .FALSE. ! flag for convergence |
---|
60 | |
---|
61 | REAL(wp) :: & |
---|
62 | brems ! relaxation. parameter |
---|
63 | |
---|
64 | REAL(wp) :: & ! temporay variables |
---|
65 | zresm, zresm1, zhipor_min , & |
---|
66 | zalk, zbor, zsil, zpo4, zdic, & |
---|
67 | zh_old, zh_old2, zh_old3, zh_old4, & |
---|
68 | zuu, zvv, zduu, zdvv , & |
---|
69 | zup, zvp, zdup, zdvp , & |
---|
70 | zah_old, zah_olds, & |
---|
71 | zh_new, zh_new2, & |
---|
72 | zco3 |
---|
73 | |
---|
74 | |
---|
75 | !!---------------------------------------------------------------------- |
---|
76 | |
---|
77 | IF( kt == nitsed000 ) THEN |
---|
78 | WRITE(numsed,*) ' sed_co3 : carbonate ion and proton concentration calculation ' |
---|
79 | WRITE(numsed,*) ' ' |
---|
80 | ENDIF |
---|
81 | |
---|
82 | itermax = 30 |
---|
83 | brems = 1. |
---|
84 | itime = 0 |
---|
85 | |
---|
86 | |
---|
87 | DO jk = 1, jpksed |
---|
88 | 10001 CONTINUE |
---|
89 | IF( itime <= 2 ) THEN |
---|
90 | lconv = .FALSE. |
---|
91 | IF( itime > 0 ) THEN |
---|
92 | ! increase max number of iterations and relaxation parameter |
---|
93 | itermax = 200 |
---|
94 | !! brems = 0.3 |
---|
95 | IF( itime == 2 ) hipor(1:jpoce,jk) = 3.e-9 ! re-initilazation of [H] values |
---|
96 | ENDIF |
---|
97 | |
---|
98 | iflag: DO jiter = 1, itermax |
---|
99 | |
---|
100 | ! Store previous hi field. |
---|
101 | zresm = -1.e10 |
---|
102 | ipt = 1 |
---|
103 | DO ji = 1, jpoce |
---|
104 | ! dissociation constant are in mol/kg of solution |
---|
105 | ! convert pwcp concentration [mol/l] in mol/kg for solution |
---|
106 | zalk = pwcp(ji,jk,jwalk) / densSW(ji) |
---|
107 | zh_old = hipor(ji,jk) / densSW(ji) |
---|
108 | zh_old2 = zh_old * zh_old |
---|
109 | zh_old3 = zh_old2 * zh_old |
---|
110 | zh_old4 = zh_old3 * zh_old |
---|
111 | zbor = borats(ji) / densSW(ji) |
---|
112 | zsil = pwcp(ji,jk,jwsil) / densSW(ji) |
---|
113 | zpo4 = pwcp(ji,jk,jwpo4) / densSW(ji) |
---|
114 | zdic = pwcp(ji,jk,jwdic) / densSW(ji) |
---|
115 | ! intermediate calculation |
---|
116 | zuu = zdic * ( ak1s(ji) / zh_old + 2.* ak12s(ji) / zh_old2 ) |
---|
117 | zvv = 1. + ak1s(ji) / zh_old + ak12s(ji) / zh_old2 |
---|
118 | zduu = zdic * ( -ak1s(ji) / zh_old2 - 4. * ak12s(ji) / zh_old3 ) |
---|
119 | zdvv = -ak1s(ji) / zh_old2 - 2. * ak12s(ji) / zh_old3 |
---|
120 | zup = zpo4 * ( ak12ps(ji) / zh_old2 + 2. * ak123ps(ji) / zh_old3 - 1.) |
---|
121 | zvp = 1. + ak1ps(ji) / zh_old + ak12ps(ji) / zh_old2 + ak123ps(ji) / zh_old3 |
---|
122 | zdup = zpo4 * ( -2. * ak12ps(ji) / zh_old3 - 6. * ak123ps(ji) / zh_old4 ) |
---|
123 | zdvp = -ak1ps(ji) / zh_old2 - 2.* ak12ps(ji) / zh_old3 - 3. * ak123ps(ji) / zh_old4 |
---|
124 | |
---|
125 | zah_old = zuu / zvv + zbor / ( 1. + zh_old / akbs(ji) ) + & |
---|
126 | & akws(ji) / zh_old - zh_old + zsil / ( 1. + zh_old / aksis(ji) ) + & |
---|
127 | & zup / zvp |
---|
128 | |
---|
129 | zah_olds = ( ( zduu * zvv - zdvv * zuu ) / ( zvv * zvv ) ) - & |
---|
130 | & zbor / akbs(ji) * ( 1. + zh_old / akbs(ji) )**(-2) - & |
---|
131 | & akws(ji) / zh_old2 - 1. - & |
---|
132 | & zsil / aksis(ji) * ( 1. + zh_old / aksis(ji) )**(-2) + & |
---|
133 | & ( ( zdup * zvp - zdvp * zup ) / ( zvp * zvp ) ) |
---|
134 | ! |
---|
135 | zh_new = zh_old - brems * ( zah_old - zalk ) / zah_olds |
---|
136 | ! |
---|
137 | zresm1 = ABS( zh_new - zh_old ) |
---|
138 | IF( zresm1 > zresm ) THEN |
---|
139 | zresm = zresm1 |
---|
140 | ENDIF |
---|
141 | ! |
---|
142 | zh_new2 = zh_new * zh_new |
---|
143 | zco3 = ( ak12s(ji) * zdic ) / ( ak12s(ji) + ak1s(ji) * zh_new + zh_new2) |
---|
144 | ! again in mol/l |
---|
145 | hipor (ji,jk) = zh_new * densSW(ji) |
---|
146 | co3por(ji,jk) = zco3 * densSW(ji) |
---|
147 | |
---|
148 | ENDDO ! end loop ji |
---|
149 | |
---|
150 | ! convergence test |
---|
151 | IF( zresm <= epsmax ) THEN |
---|
152 | lconv = .TRUE. |
---|
153 | !minimum value of hipor |
---|
154 | zhipor_min = MINVAL( hipor(1:jpoce,jk ) ) |
---|
155 | EXIT iflag |
---|
156 | ENDIF |
---|
157 | |
---|
158 | ENDDO iflag |
---|
159 | |
---|
160 | IF( lconv ) THEN |
---|
161 | ! WRITE(numsed,*) ' convergence after iter =', jiter, ' iterations ; res =',zresm |
---|
162 | IF( zhipor_min < 0. ) THEN |
---|
163 | IF ( itime == 0 ) THEN |
---|
164 | ! WRITE(numsed,*) ' but hipor < 0 ; try one more time for jk = ', jk |
---|
165 | ! WRITE(numsed,*) ' with re-initialization of initial PH field ' |
---|
166 | itime = 2 |
---|
167 | GOTO 10001 |
---|
168 | ELSE |
---|
169 | ! WRITE(numsed,*) ' convergence after iter =', jiter, ' iterations ; res =',zresm |
---|
170 | ! WRITE(numsed,*) ' but hipor < 0, again for second time for jk = ', jk |
---|
171 | ! WRITE(numsed,*) ' We stop : STOP ' |
---|
172 | STOP |
---|
173 | ENDIF |
---|
174 | ELSE |
---|
175 | ! WRITE(numsed,*) ' successfull convergence for level jk = ',jk,& |
---|
176 | ! & ' after iter =', jiter, ' iterations ; res =',zresm |
---|
177 | ! WRITE(numsed,*) ' ' |
---|
178 | itime = 0 |
---|
179 | ENDIF |
---|
180 | ELSE |
---|
181 | itime = itime + 1 |
---|
182 | WRITE(numsed,*) ' No convergence for jk = ', jk, ' after ', itime, ' try' |
---|
183 | IF ( itime == 1 ) THEN |
---|
184 | WRITE(numsed,*) ' try one more time with more iterations and higher relax. value' |
---|
185 | GOTO 10001 |
---|
186 | ELSE IF ( itime == 2 ) THEN |
---|
187 | WRITE(numsed,*) ' try one more time for with more iterations, higher relax. value' |
---|
188 | WRITE(numsed,*) ' and with re-initialization of initial PH field ' |
---|
189 | ELSE |
---|
190 | WRITE(numsed,*) ' No more... we stop ' |
---|
191 | STOP |
---|
192 | ENDIF |
---|
193 | ENDIF |
---|
194 | ENDIF |
---|
195 | ENDDO |
---|
196 | |
---|
197 | END SUBROUTINE sed_co3 |
---|
198 | #else |
---|
199 | !!====================================================================== |
---|
200 | !! MODULE sedco3 : Dummy module |
---|
201 | !!====================================================================== |
---|
202 | CONTAINS |
---|
203 | SUBROUTINE sed_co3( kt ) ! Empty routine |
---|
204 | INTEGER, INTENT(in) :: kt |
---|
205 | WRITE(*,*) 'sed_co3: You should not have seen this print! error?', kt |
---|
206 | END SUBROUTINE sed_co3 |
---|
207 | |
---|
208 | !!====================================================================== |
---|
209 | |
---|
210 | #endif |
---|
211 | |
---|
212 | END MODULE sedco3 |
---|