source: NEMO/trunk/src/TOP/PISCES/SED/sedco3.F90 @ 9598

Last change on this file since 9598 was 9598, checked in by nicolasmartin, 3 years ago

Reorganisation plan for NEMO repository: changes to make compilation succeed with new structure
Juste one issue left with AGRIF_NORDIC with AGRIF preprocessing
Standardisation of routines header with version 4.0 and year 2018
Fix for some broken symlinks

  • Property svn:keywords set to Id
File size: 8.2 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) :: epsmax   =  1.e-12      ! convergence limite value
20
21   !!----------------------------------------------------------------------
22   !!   OPA 9.0   !   NEMO Consortium   (2003)
23   !!----------------------------------------------------------------------
24
25   !! $Id$
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)  :: 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         DO WHILE( itime <= 2 )
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
86iflag:      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,*) ' '
165                  itime = 3
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
180         ENDDO ! End of WHILE LOOP
181     ENDDO
182
183   END SUBROUTINE sed_co3
184#else
185   !!======================================================================
186   !! MODULE sedco3  :   Dummy module
187   !!======================================================================
188   !! $Id$
189CONTAINS
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
199END MODULE sedco3
Note: See TracBrowser for help on using the repository browser.