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) :: epsmax = 1.e-12 ! convergence limite value |
---|
20 | |
---|
21 | !!---------------------------------------------------------------------- |
---|
22 | !! OPA 9.0 ! LODYC-IPSL (2003) |
---|
23 | !!---------------------------------------------------------------------- |
---|
24 | |
---|
25 | !! $Id$ |
---|
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) :: kt ! time step |
---|
47 | |
---|
48 | ! |
---|
49 | !---Local variables |
---|
50 | INTEGER :: jiter, ji, jk, ipt ! dummy loop indices |
---|
51 | |
---|
52 | INTEGER :: itermax ! maximum number of Newton-Raphson iterations |
---|
53 | INTEGER :: itime ! number of time to perform Newton-Raphson algorithm |
---|
54 | LOGICAL :: lconv = .FALSE. ! flag for convergence |
---|
55 | REAL(wp) :: brems ! relaxation. parameter |
---|
56 | REAL(wp) :: zresm, zresm1, zhipor_min |
---|
57 | REAL(wp) :: zalk, zbor, zsil, zpo4, zdic |
---|
58 | REAL(wp) :: zh_old, zh_old2, zh_old3, zh_old4 |
---|
59 | REAL(wp) :: zuu, zvv, zduu, zdvv |
---|
60 | REAL(wp) :: zup, zvp, zdup, zdvp |
---|
61 | REAL(wp) :: zah_old, zah_olds |
---|
62 | REAL(wp) :: zh_new, zh_new2, zco3 |
---|
63 | |
---|
64 | !!---------------------------------------------------------------------- |
---|
65 | |
---|
66 | IF( kt == nitsed000 ) THEN |
---|
67 | WRITE(numsed,*) ' sed_co3 : carbonate ion and proton concentration calculation ' |
---|
68 | WRITE(numsed,*) ' ' |
---|
69 | ENDIF |
---|
70 | |
---|
71 | itermax = 30 |
---|
72 | brems = 1. |
---|
73 | itime = 0 |
---|
74 | |
---|
75 | |
---|
76 | DO jk = 1, jpksed |
---|
77 | 10001 CONTINUE |
---|
78 | IF( itime <= 2 ) THEN |
---|
79 | lconv = .FALSE. |
---|
80 | IF( itime > 0 ) THEN |
---|
81 | ! increase max number of iterations and relaxation parameter |
---|
82 | itermax = 200 |
---|
83 | !! brems = 0.3 |
---|
84 | IF( itime == 2 ) hipor(1:jpoce,jk) = 3.e-9 ! re-initilazation of [H] values |
---|
85 | ENDIF |
---|
86 | |
---|
87 | iflag: DO jiter = 1, itermax |
---|
88 | |
---|
89 | ! Store previous hi field. |
---|
90 | zresm = -1.e10 |
---|
91 | ipt = 1 |
---|
92 | DO ji = 1, jpoce |
---|
93 | ! dissociation constant are in mol/kg of solution |
---|
94 | ! convert pwcp concentration [mol/l] in mol/kg for solution |
---|
95 | zalk = pwcp(ji,jk,jwalk) / densSW(ji) |
---|
96 | zh_old = hipor(ji,jk) / densSW(ji) |
---|
97 | zh_old2 = zh_old * zh_old |
---|
98 | zh_old3 = zh_old2 * zh_old |
---|
99 | zh_old4 = zh_old3 * zh_old |
---|
100 | zbor = borats(ji) / densSW(ji) |
---|
101 | zsil = pwcp(ji,jk,jwsil) / densSW(ji) |
---|
102 | zpo4 = pwcp(ji,jk,jwpo4) / densSW(ji) |
---|
103 | zdic = pwcp(ji,jk,jwdic) / densSW(ji) |
---|
104 | ! intermediate calculation |
---|
105 | zuu = zdic * ( ak1s(ji) / zh_old + 2.* ak12s(ji) / zh_old2 ) |
---|
106 | zvv = 1. + ak1s(ji) / zh_old + ak12s(ji) / zh_old2 |
---|
107 | zduu = zdic * ( -ak1s(ji) / zh_old2 - 4. * ak12s(ji) / zh_old3 ) |
---|
108 | zdvv = -ak1s(ji) / zh_old2 - 2. * ak12s(ji) / zh_old3 |
---|
109 | zup = zpo4 * ( ak12ps(ji) / zh_old2 + 2. * ak123ps(ji) / zh_old3 - 1.) |
---|
110 | zvp = 1. + ak1ps(ji) / zh_old + ak12ps(ji) / zh_old2 + ak123ps(ji) / zh_old3 |
---|
111 | zdup = zpo4 * ( -2. * ak12ps(ji) / zh_old3 - 6. * ak123ps(ji) / zh_old4 ) |
---|
112 | zdvp = -ak1ps(ji) / zh_old2 - 2.* ak12ps(ji) / zh_old3 - 3. * ak123ps(ji) / zh_old4 |
---|
113 | |
---|
114 | zah_old = zuu / zvv + zbor / ( 1. + zh_old / akbs(ji) ) + & |
---|
115 | & akws(ji) / zh_old - zh_old + zsil / ( 1. + zh_old / aksis(ji) ) + & |
---|
116 | & zup / zvp |
---|
117 | |
---|
118 | zah_olds = ( ( zduu * zvv - zdvv * zuu ) / ( zvv * zvv ) ) - & |
---|
119 | & zbor / akbs(ji) * ( 1. + zh_old / akbs(ji) )**(-2) - & |
---|
120 | & akws(ji) / zh_old2 - 1. - & |
---|
121 | & zsil / aksis(ji) * ( 1. + zh_old / aksis(ji) )**(-2) + & |
---|
122 | & ( ( zdup * zvp - zdvp * zup ) / ( zvp * zvp ) ) |
---|
123 | ! |
---|
124 | zh_new = zh_old - brems * ( zah_old - zalk ) / zah_olds |
---|
125 | ! |
---|
126 | zresm1 = ABS( zh_new - zh_old ) |
---|
127 | IF( zresm1 > zresm ) THEN |
---|
128 | zresm = zresm1 |
---|
129 | ENDIF |
---|
130 | ! |
---|
131 | zh_new2 = zh_new * zh_new |
---|
132 | zco3 = ( ak12s(ji) * zdic ) / ( ak12s(ji) + ak1s(ji) * zh_new + zh_new2) |
---|
133 | ! again in mol/l |
---|
134 | hipor (ji,jk) = zh_new * densSW(ji) |
---|
135 | co3por(ji,jk) = zco3 * densSW(ji) |
---|
136 | |
---|
137 | ENDDO ! end loop ji |
---|
138 | |
---|
139 | ! convergence test |
---|
140 | IF( zresm <= epsmax ) THEN |
---|
141 | lconv = .TRUE. |
---|
142 | !minimum value of hipor |
---|
143 | zhipor_min = MINVAL( hipor(1:jpoce,jk ) ) |
---|
144 | EXIT iflag |
---|
145 | ENDIF |
---|
146 | |
---|
147 | ENDDO iflag |
---|
148 | |
---|
149 | IF( lconv ) THEN |
---|
150 | ! WRITE(numsed,*) ' convergence after iter =', jiter, ' iterations ; res =',zresm |
---|
151 | IF( zhipor_min < 0. ) THEN |
---|
152 | IF ( itime == 0 ) THEN |
---|
153 | ! WRITE(numsed,*) ' but hipor < 0 ; try one more time for jk = ', jk |
---|
154 | ! WRITE(numsed,*) ' with re-initialization of initial PH field ' |
---|
155 | itime = 2 |
---|
156 | GOTO 10001 |
---|
157 | ELSE |
---|
158 | ! WRITE(numsed,*) ' convergence after iter =', jiter, ' iterations ; res =',zresm |
---|
159 | ! WRITE(numsed,*) ' but hipor < 0, again for second time for jk = ', jk |
---|
160 | ! WRITE(numsed,*) ' We stop : STOP ' |
---|
161 | STOP |
---|
162 | ENDIF |
---|
163 | ELSE |
---|
164 | ! WRITE(numsed,*) ' successfull convergence for level jk = ',jk,& |
---|
165 | ! & ' after iter =', jiter, ' iterations ; res =',zresm |
---|
166 | ! WRITE(numsed,*) ' ' |
---|
167 | itime = 0 |
---|
168 | ENDIF |
---|
169 | ELSE |
---|
170 | itime = itime + 1 |
---|
171 | WRITE(numsed,*) ' No convergence for jk = ', jk, ' after ', itime, ' try' |
---|
172 | IF ( itime == 1 ) THEN |
---|
173 | WRITE(numsed,*) ' try one more time with more iterations and higher relax. value' |
---|
174 | GOTO 10001 |
---|
175 | ELSE IF ( itime == 2 ) THEN |
---|
176 | WRITE(numsed,*) ' try one more time for with more iterations, higher relax. value' |
---|
177 | WRITE(numsed,*) ' and with re-initialization of initial PH field ' |
---|
178 | ELSE |
---|
179 | WRITE(numsed,*) ' No more... we stop ' |
---|
180 | STOP |
---|
181 | ENDIF |
---|
182 | ENDIF |
---|
183 | ENDIF |
---|
184 | ENDDO |
---|
185 | |
---|
186 | END SUBROUTINE sed_co3 |
---|
187 | #else |
---|
188 | !!====================================================================== |
---|
189 | !! MODULE sedco3 : Dummy module |
---|
190 | !!====================================================================== |
---|
191 | !! $Id$ |
---|
192 | CONTAINS |
---|
193 | SUBROUTINE sed_co3( kt ) ! Empty routine |
---|
194 | INTEGER, INTENT(in) :: kt |
---|
195 | WRITE(*,*) 'sed_co3: You should not have seen this print! error?', kt |
---|
196 | END SUBROUTINE sed_co3 |
---|
197 | |
---|
198 | !!====================================================================== |
---|
199 | |
---|
200 | #endif |
---|
201 | |
---|
202 | END MODULE sedco3 |
---|