Changeset 10345 for NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/SED/sedco3.F90
- Timestamp:
- 2018-11-21T11:25:53+01:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/SED/sedco3.F90
r9598 r10345 1 1 MODULE sedco3 2 #if defined key_sed3 2 !!====================================================================== 4 3 !! *** MODULE sedco3 *** … … 7 6 !! * Modules used 8 7 USE sed ! sediment global variable 8 USE sedchem 9 USE lib_mpp ! distribued memory computing library 9 10 10 11 … … 15 16 PUBLIC sed_co3 16 17 17 18 !! * Module variables19 REAL(wp) :: epsmax = 1.e-12 ! convergence limite value20 21 18 !!---------------------------------------------------------------------- 22 !! OPA 9.0 ! NEMO Consortium(2003)19 !! OPA 9.0 ! LODYC-IPSL (2003) 23 20 !!---------------------------------------------------------------------- 24 21 … … 45 42 !! * Arguments 46 43 INTEGER, INTENT(in) :: kt ! time step 47 48 44 ! 49 45 !---Local variables 50 INTEGER :: ji ter, ji, jk, ipt! dummy loop indices46 INTEGER :: ji, jk ! dummy loop indices 51 47 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 48 REAL(wp), DIMENSION(jpoce,jpksed) :: zhinit, zhi 64 49 !!---------------------------------------------------------------------- 65 50 51 IF( ln_timing ) CALL timing_start('sed_co3') 52 66 53 IF( kt == nitsed000 ) THEN 67 WRITE(numsed,*) ' sed_co3 : carbonate ion and proton concentration calculation '68 WRITE(numsed,*) ' '54 IF (lwp) WRITE(numsed,*) ' sed_co3 : carbonate ion and proton concentration calculation ' 55 IF (lwp) WRITE(numsed,*) ' ' 69 56 ENDIF 70 57 71 itermax = 3072 brems = 1.73 itime = 058 DO jk = 1, jpksed 59 zhinit(:,jk) = hipor(:,jk) / densSW(:) 60 END DO 74 61 62 ! ------------------------------------------- 63 ! COMPUTE [CO3--] and [H+] CONCENTRATIONS 64 ! ------------------------------------------- 65 66 CALL solve_at_general_sed(zhinit, zhi) 75 67 76 68 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 69 DO ji = 1, jpoce 70 co3por(ji,jk) = pwcp(ji,jk,jwdic) * ak1s(ji) * ak2s(ji) / (zhi(ji,jk)**2 & 71 & + ak1s(ji) * zhi(ji,jk) + ak1s(ji) * ak2s(ji) + rtrn ) 72 hipor(ji,jk) = zhi(ji,jk) * densSW(ji) 73 END DO 74 END DO 85 75 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,*) ' ' 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 76 IF( ln_timing ) CALL timing_stop('sed_co3') 182 77 183 78 END SUBROUTINE sed_co3 184 #else185 !!======================================================================186 !! MODULE sedco3 : Dummy module187 !!======================================================================188 !! $Id$189 CONTAINS190 SUBROUTINE sed_co3( kt ) ! Empty routine191 INTEGER, INTENT(in) :: kt192 WRITE(*,*) 'sed_co3: You should not have seen this print! error?', kt193 END SUBROUTINE sed_co3194 195 !!======================================================================196 197 #endif198 79 199 80 END MODULE sedco3
Note: See TracChangeset
for help on using the changeset viewer.