New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
sedco3.F90 in tags/nemo_v3_2/nemo_v3_2/NEMO/TOP_SRC/SED – NEMO

source: tags/nemo_v3_2/nemo_v3_2/NEMO/TOP_SRC/SED/sedco3.F90 @ 1878

Last change on this file since 1878 was 1878, checked in by flavoni, 14 years ago

initial test for nemogcm
File size: 8.4 KB
Line 
1MODULE 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
26CONTAINS
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
8810001    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
98iflag:      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   !!======================================================================
202CONTAINS
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
212END MODULE sedco3
Note: See TracBrowser for help on using the repository browser.