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 branches/UKMO/dev_r5518_nemo2cice_prints/NEMOGCM/NEMO/TOP_SRC/PISCES/SED – NEMO

source: branches/UKMO/dev_r5518_nemo2cice_prints/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/sedco3.F90 @ 9816

Last change on this file since 9816 was 9816, checked in by dancopsey, 6 years ago

Strip out SVN keywords.

File size: 8.3 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   !   LODYC-IPSL   (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
7710001    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
87iflag:      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$
192CONTAINS
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
202END MODULE sedco3
Note: See TracBrowser for help on using the repository browser.