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