ecearth3/sources/nemo-3.6/NEMO/TOP_SRC/PISCES/SED/sedco3.F90

MODULE sedco3
#if defined key_sed
   !!======================================================================
   !!              ***  MODULE  sedco3  ***
   !!    Sediment : carbonate in sediment pore water
   !!=====================================================================
   !! * Modules used
   USE sed     ! sediment global variable


   IMPLICIT NONE
   PRIVATE

   !! *  Routine accessibility
   PUBLIC sed_co3     


   !! * Module variables
   REAL(wp) :: epsmax   =  1.e-12      ! convergence limite value

   !!----------------------------------------------------------------------
   !!   OPA 9.0   !   LODYC-IPSL   (2003)
   !!----------------------------------------------------------------------

   !! $Id: sedco3.F90 2355 2015-05-20 07:11:50Z ufla $
CONTAINS


   SUBROUTINE sed_co3( kt )
      !!----------------------------------------------------------------------
      !!                   ***  ROUTINE sed_co3  ***
      !!
      !! ** Purpose :  carbonate ion and proton concentration 
      !!               in sediment pore water
      !!
      !! ** Methode :  - solving nonlinear equation for [H+] with given alkalinity
      !!               and total co2 
      !!               - one dimensional newton-raphson algorithm for [H+])
      !!
      !!   History :
      !!        !  98-08 (E. Maier-Reimer, Christoph Heinze )  Original code
      !!        !  04-10 (N. Emprin, M. Gehlen ) coupled with PISCES
      !!        !  06-04 (C. Ethe)  Re-organization
      !!----------------------------------------------------------------------
      !! * Arguments
      INTEGER, INTENT(in)  :: kt   ! time step

      !
      !---Local variables
      INTEGER  :: jiter, ji, jk, ipt  ! dummy loop indices

      INTEGER  :: itermax             ! maximum number of Newton-Raphson iterations
      INTEGER  :: itime               ! number of time to perform Newton-Raphson algorithm
      LOGICAL  :: lconv = .FALSE.     ! flag for convergence
      REAL(wp) :: brems               !  relaxation. parameter
      REAL(wp) :: zresm, zresm1, zhipor_min 
      REAL(wp) :: zalk, zbor, zsil, zpo4, zdic
      REAL(wp) :: zh_old, zh_old2, zh_old3, zh_old4
      REAL(wp) :: zuu, zvv, zduu, zdvv 
      REAL(wp) :: zup, zvp, zdup, zdvp
      REAL(wp) :: zah_old, zah_olds
      REAL(wp) :: zh_new, zh_new2, zco3

     !!----------------------------------------------------------------------

      IF( kt == nitsed000 ) THEN
         WRITE(numsed,*) ' sed_co3 : carbonate ion and proton concentration calculation  '
         WRITE(numsed,*) ' '
      ENDIF

      itermax     = 30
      brems       = 1.
      itime       = 0


      DO jk = 1, jpksed
10001    CONTINUE
         IF( itime <= 2 ) THEN
            lconv  = .FALSE.
            IF( itime > 0 ) THEN  
               ! increase max number of iterations and relaxation parameter
               itermax = 200
!!               brems   = 0.3
               IF( itime == 2 ) hipor(1:jpoce,jk) = 3.e-9 ! re-initilazation of [H] values
            ENDIF

iflag:      DO jiter = 1, itermax

                ! Store previous hi field.   
               zresm = -1.e10
               ipt = 1
               DO ji = 1, jpoce
                  ! dissociation constant are in mol/kg of solution
                  ! convert pwcp concentration [mol/l] in mol/kg for solution
                  zalk    = pwcp(ji,jk,jwalk) / densSW(ji)
                  zh_old  = hipor(ji,jk) / densSW(ji)
                  zh_old2 = zh_old  * zh_old
                  zh_old3 = zh_old2 * zh_old
                  zh_old4 = zh_old3 * zh_old
                  zbor    = borats(ji) / densSW(ji)
                  zsil    = pwcp(ji,jk,jwsil) / densSW(ji)
                  zpo4    = pwcp(ji,jk,jwpo4) / densSW(ji)
                  zdic    = pwcp(ji,jk,jwdic) / densSW(ji)               
                  ! intermediate calculation 
                  zuu     = zdic * ( ak1s(ji) / zh_old + 2.* ak12s(ji) / zh_old2 )
                  zvv     = 1. + ak1s(ji) / zh_old + ak12s(ji) / zh_old2
                  zduu    = zdic * ( -ak1s(ji) / zh_old2 - 4. * ak12s(ji) / zh_old3 )
                  zdvv    = -ak1s(ji) / zh_old2 - 2. * ak12s(ji) / zh_old3
                  zup     = zpo4 * ( ak12ps(ji) / zh_old2 + 2. * ak123ps(ji) / zh_old3 - 1.)
                  zvp     = 1. + ak1ps(ji) / zh_old + ak12ps(ji) / zh_old2 + ak123ps(ji) / zh_old3
                  zdup    = zpo4 * ( -2. * ak12ps(ji) / zh_old3 - 6. * ak123ps(ji) / zh_old4 )
                  zdvp    = -ak1ps(ji) / zh_old2 - 2.* ak12ps(ji) / zh_old3 - 3. * ak123ps(ji) / zh_old4
                  
                  zah_old  = zuu / zvv + zbor / ( 1. + zh_old / akbs(ji) ) + &
                     &      akws(ji) / zh_old - zh_old + zsil / ( 1. + zh_old / aksis(ji) ) + &
                     &      zup / zvp
                  
                  zah_olds = ( ( zduu * zvv - zdvv * zuu ) / ( zvv * zvv ) )      - &
                     &        zbor / akbs(ji) * ( 1. + zh_old / akbs(ji) )**(-2) - &
                     &        akws(ji) / zh_old2 - 1. -                            &
                     &        zsil / aksis(ji) * ( 1. + zh_old / aksis(ji) )**(-2) + &
                     &       ( ( zdup * zvp - zdvp * zup ) / ( zvp * zvp ) )
                  !
                  zh_new = zh_old - brems * ( zah_old - zalk ) / zah_olds
                  ! 
                  zresm1 = ABS( zh_new - zh_old )
                  IF( zresm1 > zresm ) THEN 
                     zresm = zresm1  
                  ENDIF
                  !
                  zh_new2  = zh_new * zh_new
                  zco3   = ( ak12s(ji) * zdic ) / ( ak12s(ji) + ak1s(ji) * zh_new + zh_new2)
                  ! again in mol/l
                  hipor (ji,jk) = zh_new * densSW(ji)
                  co3por(ji,jk) = zco3   * densSW(ji)
                  
               ENDDO  ! end loop ji
               
               ! convergence test
               IF( zresm <= epsmax ) THEN
                  lconv = .TRUE.
                  !minimum value of hipor
                  zhipor_min = MINVAL( hipor(1:jpoce,jk ) )
                  EXIT iflag
               ENDIF

            ENDDO iflag

            IF( lconv ) THEN
!               WRITE(numsed,*) ' convergence after iter =', jiter, ' iterations ;  res =',zresm  
               IF( zhipor_min < 0. ) THEN 
                  IF ( itime == 0 ) THEN
!                     WRITE(numsed,*) '    but hipor < 0 ; try one more time for jk = ', jk 
!                     WRITE(numsed,*) '    with re-initialization of initial PH field '       
                     itime = 2
                     GOTO 10001
                  ELSE
!                     WRITE(numsed,*) ' convergence after iter =', jiter, ' iterations ;  res =',zresm 
!                     WRITE(numsed,*) '    but hipor < 0, again for second time for jk = ', jk 
!                     WRITE(numsed,*) ' We stop : STOP '
                     STOP
                  ENDIF
               ELSE
!                  WRITE(numsed,*) ' successfull convergence for level jk = ',jk,&
!                     &               '  after iter =', jiter, ' iterations ;  res =',zresm  
!                  WRITE(numsed,*) ' '
                  itime = 0
               ENDIF
            ELSE
               itime = itime + 1
               WRITE(numsed,*) ' No convergence for jk = ', jk, ' after ', itime, '  try'            
               IF ( itime == 1 ) THEN
                  WRITE(numsed,*) ' try one more time with more iterations and higher relax. value'
                  GOTO 10001
               ELSE IF ( itime == 2 ) THEN
                  WRITE(numsed,*) ' try one more time for with more iterations, higher relax. value'               
                  WRITE(numsed,*) ' and with re-initialization of initial PH field ' 
               ELSE       
                  WRITE(numsed,*) ' No more... we stop '
                  STOP
               ENDIF
            ENDIF
         ENDIF
     ENDDO

   END SUBROUTINE sed_co3
#else
   !!======================================================================
   !! MODULE sedco3  :   Dummy module 
   !!======================================================================
   !! $Id: sedco3.F90 2355 2015-05-20 07:11:50Z ufla $
CONTAINS
   SUBROUTINE sed_co3( kt )         ! Empty routine
      INTEGER, INTENT(in) :: kt
      WRITE(*,*) 'sed_co3: You should not have seen this print! error?', kt
   END SUBROUTINE sed_co3

   !!======================================================================

#endif

END MODULE sedco3