- Timestamp:
- 2021-02-03T16:03:34+01:00 (3 years ago)
- Location:
- NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED
- Files:
-
- 1 added
- 2 deleted
- 11 edited
-
par_sed.F90 (modified) (1 diff)
-
sed.F90 (modified) (4 diffs)
-
sedadv.F90 (modified) (4 diffs)
-
sedarr.F90 (modified) (1 diff)
-
seddiff.F90 (modified) (2 diffs)
-
seddsr.F90 (modified) (21 diffs)
-
seddta.F90 (modified) (6 diffs)
-
sedini.F90 (modified) (7 diffs)
-
sedinorg.F90 (modified) (9 diffs)
-
sedmbc.F90 (deleted)
-
sedorg.F90 (deleted)
-
sedsol.F90 (added)
-
sedstp.F90 (modified) (3 diffs)
-
sedwri.F90 (modified) (5 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/par_sed.F90
r14086 r14385 60 60 INTEGER, PUBLIC, PARAMETER :: & 61 61 jptrased = jpsol + jpwat , & 62 jpdia3dsed = 2, &63 jpdia2dsed = 1262 jpdia3dsed = 4 , & 63 jpdia2dsed = 20 64 64 65 65 END MODULE par_sed -
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/sed.F90
r14086 r14385 70 70 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: fromsed !: 71 71 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: tosed !: 72 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: rloss !:73 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: tokbot74 !75 72 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: temp !: temperature 76 73 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: salt !: salinity 77 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: press !: pressure78 74 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: raintg !: total massic flux rained in each cell (sum of sol. comp.) 79 75 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: fecratio !: Fe/C ratio in falling particles to the sediments 80 76 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: dzdep !: total thickness of solid material rained [cm] in each cell 81 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: zkbot !: total thickness of solid material rained [cm] in each cell77 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: zkbot !: Total water column depth 82 78 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: wacc !: total thickness of solid material rained [cm] in each cell 83 79 ! … … 88 84 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: vols3d !: ??? 89 85 90 91 86 !! Chemistry 92 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: densSW 93 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: borats 87 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: densSW 88 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: borats 94 89 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: calcon2 95 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: akbs 96 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: ak1s 97 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: ak2s 98 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: akws 99 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: ak12s 100 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: ak1ps 101 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: ak2ps 102 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: ak3ps 103 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: ak12ps 90 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: akbs 91 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: ak1s 92 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: ak2s 93 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: akws 94 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: ak12s 95 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: ak1ps 96 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: ak2ps 97 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: ak3ps 98 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: ak12ps 104 99 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: ak123ps 105 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: aksis 106 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: aksps 100 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: aksis 101 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: aksps 107 102 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: sieqs 108 103 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: aks3s … … 127 122 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: db !: bioturbation ceofficient 128 123 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: irrig !: bioturbation ceofficient 129 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: rdtsed !: sediment model time-step130 124 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: sedligand 125 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: saturco3 131 126 REAL(wp) :: dens !: density of solid material 132 127 !! Inputs / Outputs … … 154 149 !!------------------------------------------------------------------- 155 150 ! 156 ALLOCATE( trc_data(jpi,jpj,jpdta) , & 157 & epkbot(jpi,jpj), gdepbot(jpi,jpj) , & 158 & dz(jpksed) , por(jpksed) , por1(jpksed) , & 159 & volw(jpksed), vols(jpksed), rdtsed(jpksed) , & 160 & trcsedi (jpi,jpj,jpksed,jptrased) , & 161 & flxsedi3d(jpi,jpj,jpksed,jpdia3dsed) , & 162 & flxsedi2d(jpi,jpj,jpdia2dsed) , & 163 & mol_wgt(jpsol), STAT=sed_alloc ) 151 ALLOCATE( trc_data(jpi,jpj,jpdta), trcsedi (jpi,jpj,jpksed,jptrased) , & 152 & epkbot(jpi,jpj), gdepbot(jpi,jpj), dz(jpksed), por(jpksed) , & 153 & por1(jpksed), volw(jpksed), vols(jpksed), mol_wgt(jpsol) , & 154 & flxsedi2d(jpi,jpj,jpdia2dsed), flxsedi3d(jpi,jpj,jpksed,jpdia3dsed), STAT=sed_alloc ) 164 155 165 156 IF( sed_alloc /= 0 ) CALL ctl_stop( 'STOP', 'sed_alloc: failed to allocate arrays' ) -
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/sedadv.F90
r10425 r14385 88 88 fromsed(:,:) = 0. 89 89 tosed (:,:) = 0. 90 rloss (:,:) = 0.91 90 ikwneg = 1 92 91 nztime = jpksed … … 254 253 255 254 ELSE IF( raintg(ji) < eps ) THEN ! rain = 0 256 !! Nadia rloss(:,:) = rainrm(:,:) bug ??????257 258 rloss(ji,1:jpsol) = rainrm(ji,1:jpsol)259 255 260 256 zfull (2) = zfilled(2) … … 300 296 IF( ikwneg == 2 ) THEN ! advection is reversed in the first sediment layer 301 297 302 zkwnup = rdtsed(ikwneg) * raintg(ji) / dz(ikwneg)298 zkwnup = dtsed * raintg(ji) / ( denssol * por1(ikwneg) * dz(ikwneg) ) 303 299 zkwnlo = ABS( zwb(ji,ikwneg) ) / dz(ikwneg) 304 300 zfull (ikwneg+1) = zfilled(ikwneg+1) - zkwnlo * dvolsm(ikwneg+1) … … 347 343 ! Heinze fromsed(ji,jsclay) = zkwnlo * 1. * denssol * por1(jpksed) / mol_wgt(jsclay) 348 344 fromsed(ji,jsclay) = zkwnlo * 1.* por1clay 349 ELSE ! 2 < ikwneg(ji) <= jpksedm1345 ELSE IF( ikwneg > 2 .AND. ikwneg < jpksed-1 ) THEN 350 346 351 347 zkwnup = ABS( zwb(ji,ikwneg-1) ) * por1(ikwneg-1) / ( dz(ikwneg) * por1(ikwneg) ) -
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/sedarr.F90
r14086 r14385 51 51 jid = MOD( tab_ind(jn) - 1, jpi ) + 1 52 52 jjd = ( tab_ind(jn) - 1 ) / jpi + 1 53 ! IF ( mig(jid) == 112 .and. mjg(jjd) == 25) write(0,*) 'plante ',jn,ndim1d 53 54 tab1d(jn) = tab2d(jid, jjd) 54 55 END DO -
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/seddiff.F90
r10225 r14385 43 43 INTEGER, INTENT(in) :: kt, knt ! number of iteration 44 44 ! --- local variables 45 INTEGER :: ji, jk, j s! dummy looop indices45 INTEGER :: ji, jk, jw ! dummy looop indices 46 46 47 47 REAL(wp), DIMENSION(jpoce,jpksed) :: zrearat1, zrearat2 ! reaction rate in pore water … … 60 60 ! Initializations 61 61 !---------------------- 62 zrearat1(:,:) 62 zrearat1(:,:) = 0. 63 63 zrearat2(:,:) = 0. 64 64 65 ! ---------------------------66 ! Solves PO4 diffusion67 ! ----------------------------65 ! -------------------------------------- 66 ! Solves solute diffusion in pore waters 67 ! -------------------------------------- 68 68 69 ! solves tridiagonal system 70 CALL sed_mat( jwpo4, jpoce, jpksed, zrearat1, zrearat2, pwcp(:,:,jwpo4), dtsed2 / 2.0 ) 71 72 !--------------------------- 73 ! Solves NH4 diffusion 74 !---------------------------- 75 76 ! solves tridiagonal system 77 CALL sed_mat( jwnh4, jpoce, jpksed, zrearat1, zrearat2, pwcp(:,:,jwnh4), dtsed2 / 2.0 ) 78 79 !--------------------------- 80 ! Solves Fe2+ diffusion 81 !---------------------------- 82 83 ! solves tridiagonal system 84 CALL sed_mat( jwfe2, jpoce, jpksed, zrearat1, zrearat2, pwcp(:,:,jwfe2), dtsed2 / 2.0 ) 85 86 !--------------------------- 87 ! Solves H2S diffusion 88 !---------------------------- 89 90 ! solves tridiagonal system 91 CALL sed_mat( jwh2s, jpoce, jpksed, zrearat1, zrearat2, pwcp(:,:,jwh2s), dtsed2 / 2.0 ) 92 93 !--------------------------- 94 ! Solves SO4 diffusion 95 !---------------------------- 96 97 ! solves tridiagonal system 98 CALL sed_mat( jwso4, jpoce, jpksed, zrearat1, zrearat2, pwcp(:,:,jwso4), dtsed2 / 2.0 ) 99 100 !--------------------------- 101 ! Solves O2 diffusion 102 !---------------------------- 103 104 ! solves tridiagonal system 105 CALL sed_mat( jwoxy, jpoce, jpksed, zrearat1, zrearat2, pwcp(:,:,jwoxy), dtsed2 / 2.0 ) 106 107 !--------------------------- 108 ! Solves NO3 diffusion 109 !---------------------------- 110 111 ! solves tridiagonal system 112 CALL sed_mat( jwno3, jpoce, jpksed, zrearat1, zrearat2, pwcp(:,:,jwno3), dtsed2 / 2.0 ) 69 ! Solves tridiagonal system 70 DO jw = 1, jpwat 71 CALL sed_mat( jw, jpoce, jpksed, zrearat1, zrearat2, pwcp(:,:,jw), dtsed2 / 2.0 ) 72 END DO 113 73 114 74 CALL sed_mat( jwdic, jpoce, jpksed, zrearat1, zrearat2, sedligand(:,:), dtsed2 / 2.0 ) -
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/seddsr.F90
r10362 r14385 55 55 REAL(wp), DIMENSION(jpoce,jpksed) :: zrearat1, zrearat2, zrearat3 ! reaction rate in pore water 56 56 REAL(wp), DIMENSION(jpoce,jpksed) :: zundsat ! undersaturation ; indice jpwatp1 is for calcite 57 REAL(wp), DIMENSION(jpoce,jpksed) :: z kpoc, zkpos, zkpor, zlimo2, zlimno3, zlimso4, zlimfeo ! undersaturation ; indice jpwatp1 is for calcite57 REAL(wp), DIMENSION(jpoce,jpksed) :: zlimo2, zlimno3, zlimso4, zlimfeo ! undersaturation ; indice jpwatp1 is for calcite 58 58 REAL(wp), DIMENSION(jpoce) :: zsumtot 59 59 REAL(wp) :: zsolid1, zsolid2, zsolid3, zvolw, zreasat 60 REAL(wp) :: zsatur, zsatur2, znusil, zkpoca, zkpocb, zkpocc 61 REAL(wp) :: zratio, zgamma, zbeta, zlimtmp, zundsat2 60 REAL(wp) :: zgamma, zbeta, zlimtmp 62 61 !! 63 62 !!---------------------------------------------------------------------- … … 75 74 !---------------------- 76 75 77 zrearat1(:,:) = 0. ; zundsat(:,:) = 0. ; zkpoc(:,:) = 0. 78 zlimo2 (:,:) = 0. ; zlimno3(:,:) = 0. ; zrearat2(:,:) = 0. 79 zlimso4(:,:) = 0. ; zkpor(:,:) = 0. ; zrearat3(:,:) = 0. 80 zkpos (:,:) = 0. 76 zrearat1(:,:) = 0. ; zundsat(:,:) = 0. 77 zlimo2 (:,:) = 0. ; zlimno3(:,:) = 0. ; zrearat2(:,:) = 0. 78 zlimso4(:,:) = 0. ; zrearat3(:,:) = 0. 81 79 zsumtot(:) = rtrn 82 80 … … 84 82 zvolc(:,:,:) = 0. 85 83 zadsnh4 = 1.0 / ( 1.0 + adsnh4 ) 86 87 ! Inhibition terms for the different redox equations88 ! --------------------------------------------------89 DO jk = 1, jpksed90 DO ji = 1, jpoce91 zkpoc(ji,jk) = reac_pocl92 zkpos(ji,jk) = reac_pocs93 zkpor(ji,jk) = reac_pocr94 END DO95 END DO96 84 97 85 ! Conversion of volume units … … 116 104 zrearat3(:,:) = 0. 117 105 118 zundsat(:,:) = pwcp(:,:,jwoxy)106 zundsat(:,:) = MAX( pwcp(:,:,jwoxy) - rtrn, 0. ) 119 107 120 108 DO jk = 2, jpksed 121 109 DO ji = 1, jpoce 122 zlimo2(ji,jk) = 1.0 / ( zundsat(ji,jk) + xksedo2 )123 110 zsolid1 = zvolc(ji,jk,jspoc) * solcp(ji,jk,jspoc) 124 111 zsolid2 = zvolc(ji,jk,jspos) * solcp(ji,jk,jspos) 125 112 zsolid3 = zvolc(ji,jk,jspor) * solcp(ji,jk,jspor) 126 zkpoca = zkpoc(ji,jk) * zlimo2(ji,jk) 127 zkpocb = zkpos(ji,jk) * zlimo2(ji,jk) 128 zkpocc = zkpor(ji,jk) * zlimo2(ji,jk) 129 zrearat1(ji,jk) = ( zkpoc(ji,jk) * dtsed2 * zsolid1 ) / & 130 & ( 1. + zkpoca * zundsat(ji,jk ) * dtsed2 ) 131 zrearat2(ji,jk) = ( zkpos(ji,jk) * dtsed2 * zsolid2 ) / & 132 & ( 1. + zkpocb * zundsat(ji,jk ) * dtsed2 ) 133 zrearat3(ji,jk) = ( zkpor(ji,jk) * dtsed2 * zsolid3 ) / & 134 & ( 1. + zkpocc * zundsat(ji,jk ) * dtsed2 ) 135 ENDDO 136 ENDDO 137 138 ! left hand side of coefficient matrix 139 ! DO jn = 1, 5 140 DO jk = 2, jpksed 141 DO ji = 1, jpoce 142 jflag1: DO jn = 1, 10 143 zsolid1 = zvolc(ji,jk,jspoc) * solcp(ji,jk,jspoc) 144 zsolid2 = zvolc(ji,jk,jspos) * solcp(ji,jk,jspos) 145 zsolid3 = zvolc(ji,jk,jspor) * solcp(ji,jk,jspor) 146 zbeta = xksedo2 - pwcp(ji,jk,jwoxy) + so2ut * ( zrearat1(ji,jk) & 147 & + zrearat2(ji,jk) + zrearat3(ji,jk) ) 148 zgamma = - xksedo2 * pwcp(ji,jk,jwoxy) 149 zundsat2 = zundsat(ji,jk) 150 zundsat(ji,jk) = ( - zbeta + SQRT( zbeta**2 - 4.0 * zgamma ) ) / 2.0 151 zlimo2(ji,jk) = 1.0 / ( zundsat(ji,jk) + xksedo2 ) 152 zkpoca = zkpoc(ji,jk) * zlimo2(ji,jk) 153 zkpocb = zkpos(ji,jk) * zlimo2(ji,jk) 154 zkpocc = zkpor(ji,jk) * zlimo2(ji,jk) 155 zrearat1(ji,jk) = ( zkpoc(ji,jk) * dtsed2 * zsolid1 ) / & 156 & ( 1. + zkpoca * zundsat(ji,jk ) * dtsed2 ) 157 zrearat2(ji,jk) = ( zkpos(ji,jk) * dtsed2 * zsolid2 ) / & 158 & ( 1. + zkpocb * zundsat(ji,jk ) * dtsed2 ) 159 zrearat3(ji,jk) = ( zkpor(ji,jk) * dtsed2 * zsolid3 ) / & 160 & ( 1. + zkpocc * zundsat(ji,jk ) * dtsed2 ) 161 IF ( ABS( (zundsat(ji,jk)-zundsat2)/(zundsat2+rtrn)) < 1E-8 ) THEN 162 EXIT jflag1 163 ENDIF 164 END DO jflag1 165 END DO 166 END DO 167 168 ! New solid concentration values (jk=2 to jksed) for each couple 169 DO jk = 2, jpksed 170 DO ji = 1, jpoce 171 zreasat = zrearat1(ji,jk) * zlimo2(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspoc) 172 solcp(ji,jk,jspoc) = solcp(ji,jk,jspoc) - zreasat 173 zreasat = zrearat2(ji,jk) * zlimo2(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspos) 174 solcp(ji,jk,jspos) = solcp(ji,jk,jspos) - zreasat 175 zreasat = zrearat3(ji,jk) * zlimo2(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspor) 176 solcp(ji,jk,jspor) = solcp(ji,jk,jspor) - zreasat 177 ENDDO 178 ENDDO 179 180 ! New pore water concentrations 181 DO jk = 2, jpksed 182 DO ji = 1, jpoce 183 ! Acid Silicic 184 pwcp(ji,jk,jwoxy) = zundsat(ji,jk) 185 zreasat = ( zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) ) * zlimo2(ji,jk) * zundsat(ji,jk) ! oxygen 113 zrearat1(ji,jk) = ( reac_pocl * dtsed2 * zsolid1 ) / ( 1. + reac_pocl * dtsed2 ) 114 zrearat2(ji,jk) = ( reac_pocs * dtsed2 * zsolid2 ) / ( 1. + reac_pocs * dtsed2 ) 115 zrearat3(ji,jk) = ( reac_pocr * dtsed2 * zsolid3 ) / ( 1. + reac_pocr * dtsed2 ) 116 solcp(ji,jk,jspoc) = solcp(ji,jk,jspoc) - zrearat1(ji,jk) / zvolc(ji,jk,jspoc) 117 solcp(ji,jk,jspos) = solcp(ji,jk,jspos) - zrearat2(ji,jk) / zvolc(ji,jk,jspos) 118 solcp(ji,jk,jspor) = solcp(ji,jk,jspor) - zrearat3(ji,jk) / zvolc(ji,jk,jspor) 119 zreasat = zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) 120 186 121 ! For DIC 187 122 pwcp(ji,jk,jwdic) = pwcp(ji,jk,jwdic) + zreasat 188 123 zsumtot(ji) = zsumtot(ji) + zreasat / dtsed2 * volw3d(ji,jk) * 1.e-3 * 86400. * 365. * 1E3 124 125 ! For alkalinity 126 pwcp(ji,jk,jwalk) = pwcp(ji,jk,jwalk) + zreasat * ( srno3 - 2.* spo4r ) 127 189 128 ! For Phosphate (in mol/l) 190 129 pwcp(ji,jk,jwpo4) = pwcp(ji,jk,jwpo4) + zreasat * spo4r 130 191 131 ! For iron (in mol/l) 192 132 pwcp(ji,jk,jwfe2) = pwcp(ji,jk,jwfe2) + fecratio(ji) * zreasat 193 ! For alkalinity 194 pwcp(ji,jk,jwalk) = pwcp(ji,jk,jwalk) + zreasat * ( srno3 * zadsnh4 - 2.* spo4r ) 133 195 134 ! Ammonium 196 135 pwcp(ji,jk,jwnh4) = pwcp(ji,jk,jwnh4) + zreasat * srno3 * zadsnh4 136 197 137 ! Ligands 198 sedligand(ji,jk) = sedligand(ji,jk) + ratligc * zreasat - reac_ligc * sedligand(ji,jk) 138 sedligand(ji,jk) = sedligand(ji,jk) + ratligc * zreasat 139 140 ENDDO 141 ENDDO 142 143 ! left hand side of coefficient matrix 144 DO jk = 2, jpksed 145 DO ji = 1, jpoce 146 zreasat = zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) 147 zbeta = xksedo2 - pwcp(ji,jk,jwoxy) + so2ut * zreasat 148 zgamma = - xksedo2 * pwcp(ji,jk,jwoxy) 149 zundsat(ji,jk) = ( - zbeta + SQRT( zbeta**2 - 4.0 * zgamma ) ) / 2.0 150 zlimo2(ji,jk) = zundsat(ji,jk) / ( zundsat(ji,jk) + xksedo2 ) 151 152 ! Oxygen 153 pwcp(ji,jk,jwoxy) = zundsat(ji,jk) + MIN( pwcp(ji,jk,jwoxy), rtrn ) 154 zreasat = zreasat * zlimo2(ji,jk) 155 156 ! Ligands 157 sedligand(ji,jk) = sedligand(ji,jk) - reac_ligc * sedligand(ji,jk) 199 158 ENDDO 200 159 ENDDO … … 205 164 !-------------------------------------------------------------------- 206 165 207 zrearat1(:,:) = 0. 208 zrearat2(:,:) = 0. 209 zrearat3(:,:) = 0. 210 211 zundsat(:,:) = pwcp(:,:,jwno3) 166 zundsat(:,:) = MAX(0., pwcp(:,:,jwno3) - rtrn) 212 167 213 168 DO jk = 2, jpksed 214 169 DO ji = 1, jpoce 215 zlimno3(ji,jk) = ( 1.0 - pwcp(ji,jk,jwoxy) * zlimo2(ji,jk) ) / ( zundsat(ji,jk) + xksedno3 ) 216 zsolid1 = zvolc(ji,jk,jspoc) * solcp(ji,jk,jspoc) 217 zsolid2 = zvolc(ji,jk,jspos) * solcp(ji,jk,jspos) 218 zsolid3 = zvolc(ji,jk,jspor) * solcp(ji,jk,jspor) 219 zkpoca = zkpoc(ji,jk) * zlimno3(ji,jk) 220 zkpocb = zkpos(ji,jk) * zlimno3(ji,jk) 221 zkpocc = zkpor(ji,jk) * zlimno3(ji,jk) 222 zrearat1(ji,jk) = ( zkpoc(ji,jk) * dtsed2 * zsolid1 ) / & 223 & ( 1. + zkpoca * zundsat(ji,jk ) * dtsed2 ) 224 zrearat2(ji,jk) = ( zkpos(ji,jk) * dtsed2 * zsolid2 ) / & 225 & ( 1. + zkpocb * zundsat(ji,jk ) * dtsed2 ) 226 zrearat3(ji,jk) = ( zkpor(ji,jk) * dtsed2 * zsolid3 ) / & 227 & ( 1. + zkpocc * zundsat(ji,jk ) * dtsed2 ) 228 END DO 229 END DO 230 231 ! DO jn = 1, 5 232 DO jk = 2, jpksed 233 DO ji = 1, jpoce 234 jflag2: DO jn = 1, 10 235 zlimtmp = ( 1.0 - pwcp(ji,jk,jwoxy) * zlimo2(ji,jk) ) 236 zsolid1 = zvolc(ji,jk,jspoc) * solcp(ji,jk,jspoc) 237 zsolid2 = zvolc(ji,jk,jspos) * solcp(ji,jk,jspos) 238 zsolid3 = zvolc(ji,jk,jspor) * solcp(ji,jk,jspor) 239 zbeta = xksedno3 - pwcp(ji,jk,jwno3) + srDnit * ( zrearat1(ji,jk) & 240 & + zrearat2(ji,jk) + zrearat3(ji,jk) ) * zlimtmp 241 zgamma = - xksedno3 * pwcp(ji,jk,jwno3) 242 zundsat2 = zundsat(ji,jk) 243 zundsat(ji,jk) = ( - zbeta + SQRT( zbeta**2 - 4.0 * zgamma ) ) / 2.0 244 zlimno3(ji,jk) = ( 1.0 - pwcp(ji,jk,jwoxy) * zlimo2(ji,jk) ) / ( zundsat(ji,jk) + xksedno3 ) 245 zkpoca = zkpoc(ji,jk) * zlimno3(ji,jk) 246 zkpocb = zkpos(ji,jk) * zlimno3(ji,jk) 247 zkpocc = zkpor(ji,jk) * zlimno3(ji,jk) 248 zrearat1(ji,jk) = ( zkpoc(ji,jk) * dtsed2 * zsolid1 ) / & 249 & ( 1. + zkpoca * zundsat(ji,jk ) * dtsed2 ) 250 zrearat2(ji,jk) = ( zkpos(ji,jk) * dtsed2 * zsolid2 ) / & 251 & ( 1. + zkpocb * zundsat(ji,jk ) * dtsed2 ) 252 zrearat3(ji,jk) = ( zkpor(ji,jk) * dtsed2 * zsolid3 ) / & 253 & ( 1. + zkpocc * zundsat(ji,jk ) * dtsed2 ) 254 IF ( ABS( (zundsat(ji,jk)-zundsat2)/(zundsat2+rtrn)) < 1E-8 ) THEN 255 EXIT jflag2 256 ENDIF 257 END DO jflag2 258 END DO 259 END DO 260 261 262 ! New solid concentration values (jk=2 to jksed) for each couple 263 DO jk = 2, jpksed 264 DO ji = 1, jpoce 265 zreasat = zrearat1(ji,jk) * zlimno3(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspoc) 266 solcp(ji,jk,jspoc) = solcp(ji,jk,jspoc) - zreasat 267 zreasat = zrearat2(ji,jk) * zlimno3(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspos) 268 solcp(ji,jk,jspos) = solcp(ji,jk,jspos) - zreasat 269 zreasat = zrearat3(ji,jk) * zlimno3(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspor) 270 solcp(ji,jk,jspor) = solcp(ji,jk,jspor) - zreasat 271 ENDDO 272 ENDDO 273 274 ! New dissolved concentrations 275 DO jk = 2, jpksed 276 DO ji = 1, jpoce 170 zlimtmp = ( 1.0 - zlimo2(ji,jk) ) 171 zreasat = zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) 172 zbeta = xksedno3 - pwcp(ji,jk,jwno3) + srDnit * zreasat * zlimtmp 173 zgamma = - xksedno3 * pwcp(ji,jk,jwno3) 174 zundsat(ji,jk) = ( - zbeta + SQRT( zbeta**2 - 4.0 * zgamma ) ) / 2.0 175 zlimno3(ji,jk) = zlimtmp * zundsat(ji,jk) / ( zundsat(ji,jk) + xksedno3 ) 176 277 177 ! For nitrates 278 pwcp(ji,jk,jwno3) = zundsat(ji,jk) 279 zreasat = ( zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) ) * zlimno3(ji,jk) * zundsat(ji,jk) 280 ! For DIC 281 pwcp(ji,jk,jwdic) = pwcp(ji,jk,jwdic) + zreasat 282 zsumtot(ji) = zsumtot(ji) + zreasat / dtsed2 * volw3d(ji,jk) * 1.e-3 * 86400. * 365. * 1E3 283 ! For Phosphate (in mol/l) 284 pwcp(ji,jk,jwpo4) = pwcp(ji,jk,jwpo4) + zreasat * spo4r 285 ! Ligands 286 sedligand(ji,jk) = sedligand(ji,jk) + ratligc * zreasat 287 ! For iron (in mol/l) 288 pwcp(ji,jk,jwfe2) = pwcp(ji,jk,jwfe2) + fecratio(ji) * zreasat 178 pwcp(ji,jk,jwno3) = zundsat(ji,jk) + MIN(pwcp(ji,jk,jwno3), rtrn) 179 zreasat = zreasat * zlimno3(ji,jk) 180 289 181 ! For alkalinity 290 pwcp(ji,jk,jwalk) = pwcp(ji,jk,jwalk) + zreasat * ( srDnit + srno3 * zadsnh4 - 2.* spo4r ) 291 ! Ammonium 292 pwcp(ji,jk,jwnh4) = pwcp(ji,jk,jwnh4) + zreasat * srno3 * zadsnh4 182 pwcp(ji,jk,jwalk) = pwcp(ji,jk,jwalk) + zreasat * srDnit 293 183 ENDDO 294 184 ENDDO … … 299 189 !-------------------------------------------------------------------- 300 190 301 zrearat1(:,:) = 0. 302 zrearat2(:,:) = 0. 303 zrearat3(:,:) = 0. 304 305 zundsat(:,:) = solcp(:,:,jsfeo) 191 zundsat(:,:) = MAX(0., solcp(:,:,jsfeo) - rtrn) 306 192 307 193 DO jk = 2, jpksed 308 194 DO ji = 1, jpoce 309 zlimfeo(ji,jk) = ( 1.0 - pwcp(ji,jk,jwoxy) * zlimo2(ji,jk) ) * ( 1.0 - pwcp(ji,jk,jwno3) & 310 & / ( pwcp(ji,jk,jwno3) + xksedno3 ) ) / ( zundsat(ji,jk) + xksedfeo ) 311 zsolid1 = zvolc(ji,jk,jspoc) * solcp(ji,jk,jspoc) 312 zsolid2 = zvolc(ji,jk,jspos) * solcp(ji,jk,jspos) 313 zsolid3 = zvolc(ji,jk,jspor) * solcp(ji,jk,jspor) 314 zkpoca = zkpoc(ji,jk) * zlimfeo(ji,jk) 315 zkpocb = zkpos(ji,jk) * zlimfeo(ji,jk) 316 zkpocc = zkpor(ji,jk) * zlimfeo(ji,jk) 317 zrearat1(ji,jk) = ( zkpoc(ji,jk) * dtsed2 * zsolid1 ) / & 318 & ( 1. + zkpoca * zundsat(ji,jk) * dtsed2 ) 319 zrearat2(ji,jk) = ( zkpos(ji,jk) * dtsed2 * zsolid2 ) / & 320 & ( 1. + zkpocb * zundsat(ji,jk) * dtsed2 ) 321 zrearat3(ji,jk) = ( zkpor(ji,jk) * dtsed2 * zsolid3 ) / & 322 & ( 1. + zkpocc * zundsat(ji,jk) * dtsed2 ) 323 END DO 324 END DO 325 326 ! DO jn = 1, 5 327 DO jk = 2, jpksed 328 DO ji = 1, jpoce 329 jflag3: DO jn = 1, 10 330 zlimtmp = ( 1.0 - pwcp(ji,jk,jwoxy) * zlimo2(ji,jk) ) * ( 1.0 - pwcp(ji,jk,jwno3) & 331 & / ( pwcp(ji,jk,jwno3) + xksedno3 ) ) 332 zsolid1 = zvolc(ji,jk,jspoc) * solcp(ji,jk,jspoc) 333 zsolid2 = zvolc(ji,jk,jspos) * solcp(ji,jk,jspos) 334 zsolid3 = zvolc(ji,jk,jspor) * solcp(ji,jk,jspor) 335 zreasat = ( zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) ) / zvolc(ji,jk,jsfeo) 336 zbeta = xksedfeo - solcp(ji,jk,jsfeo) + 4.0 * zreasat * zlimtmp 337 zgamma = -xksedfeo * solcp(ji,jk,jsfeo) 338 zundsat2 = zundsat(ji,jk) 339 zundsat(ji,jk) = ( - zbeta + SQRT( zbeta**2 - 4.0 * zgamma ) ) / 2.0 340 zlimfeo(ji,jk) = ( 1.0 - pwcp(ji,jk,jwoxy) * zlimo2(ji,jk) ) * ( 1.0 - pwcp(ji,jk,jwno3) & 341 & / ( pwcp(ji,jk,jwno3) + xksedno3 ) ) / ( zundsat(ji,jk) + xksedfeo ) 342 zkpoca = zkpoc(ji,jk) * zlimfeo(ji,jk) 343 zkpocb = zkpos(ji,jk) * zlimfeo(ji,jk) 344 zkpocc = zkpor(ji,jk) * zlimfeo(ji,jk) 345 zrearat1(ji,jk) = ( zkpoc(ji,jk) * dtsed2 * zsolid1 ) / & 346 & ( 1. + zkpoca * zundsat(ji,jk) * dtsed2 ) 347 zrearat2(ji,jk) = ( zkpos(ji,jk) * dtsed2 * zsolid2 ) / & 348 & ( 1. + zkpocb * zundsat(ji,jk) * dtsed2 ) 349 zrearat3(ji,jk) = ( zkpor(ji,jk) * dtsed2 * zsolid3 ) / & 350 & ( 1. + zkpocc * zundsat(ji,jk) * dtsed2 ) 351 IF ( ABS( (zundsat(ji,jk)-zundsat2)/( MAX(0.,zundsat2)+rtrn)) < 1E-8 ) THEN 352 EXIT jflag3 353 ENDIF 354 END DO jflag3 355 END DO 356 END DO 357 358 359 ! New solid concentration values (jk=2 to jksed) for each couple 360 DO jk = 2, jpksed 361 DO ji = 1, jpoce 362 zreasat = zrearat1(ji,jk) * zlimfeo(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspoc) 363 solcp(ji,jk,jspoc) = solcp(ji,jk,jspoc) - zreasat 364 zreasat = zrearat2(ji,jk) * zlimfeo(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspos) 365 solcp(ji,jk,jspos) = solcp(ji,jk,jspos) - zreasat 366 zreasat = zrearat3(ji,jk) * zlimfeo(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspor) 367 solcp(ji,jk,jspor) = solcp(ji,jk,jspor) - zreasat 368 END DO 369 END DO 370 371 ! New dissolved concentrations 372 DO jk = 2, jpksed 373 DO ji = 1, jpoce 374 zreasat = ( zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) ) * zlimfeo(ji,jk) * zundsat(ji,jk) 195 zlimtmp = 1.0 - zlimo2(ji,jk) - zlimno3(ji,jk) 196 zreasat = ( zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) ) / zvolc(ji,jk,jsfeo) 197 zbeta = xksedfeo - solcp(ji,jk,jsfeo) + 4.0 * zreasat * zlimtmp 198 zgamma = -xksedfeo * solcp(ji,jk,jsfeo) 199 zundsat(ji,jk) = ( - zbeta + SQRT( zbeta**2 - 4.0 * zgamma ) ) / 2.0 200 zlimfeo(ji,jk) = zlimtmp * zundsat(ji,jk) / ( zundsat(ji,jk) + xksedfeo ) 201 zreasat = ( zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) ) * zlimfeo(ji,jk) 202 375 203 ! For FEOH 376 solcp(ji,jk,jsfeo) = zundsat(ji,jk) 377 ! For DIC 378 pwcp(ji,jk,jwdic) = pwcp(ji,jk,jwdic) + zreasat 379 zsumtot(ji) = zsumtot(ji) + zreasat / dtsed2 * volw3d(ji,jk) * 1.e-3 * 86400. * 365. * 1E3 204 solcp(ji,jk,jsfeo) = zundsat(ji,jk) + MIN(solcp(ji,jk,jsfeo), rtrn) 205 380 206 ! For Phosphate (in mol/l) 381 pwcp(ji,jk,jwpo4) = pwcp(ji,jk,jwpo4) + zreasat * ( spo4r + 4.0 * redfep ) 382 ! Ligands 383 sedligand(ji,jk) = sedligand(ji,jk) + ratligc * zreasat 384 ! For iron (in mol/l) 385 pwcp(ji,jk,jwfe2) = pwcp(ji,jk,jwfe2) + fecratio(ji) * zreasat 207 pwcp(ji,jk,jwpo4) = pwcp(ji,jk,jwpo4) + zreasat * 4.0 * redfep 208 386 209 ! For alkalinity 387 pwcp(ji,jk,jwalk) = pwcp(ji,jk,jwalk) + zreasat * ( srno3 * zadsnh4 - 2.* spo4r ) +8.0 * zreasat388 ! Ammonium 389 pwcp(ji,jk,jwnh4) = pwcp(ji,jk,jwnh4) + zreasat * srno3 * zadsnh4210 pwcp(ji,jk,jwalk) = pwcp(ji,jk,jwalk) + 8.0 * zreasat 211 212 ! Iron 390 213 pwcp(ji,jk,jwfe2) = pwcp(ji,jk,jwfe2) + zreasat * 4.0 391 214 ENDDO 392 215 ENDDO 393 216 394 !-------------------------------------------------------------------- 395 ! Begining POC denitrification and NO3- diffusion 396 ! (indice n�5 for couple POC/NO3- ie solcp(:,:,jspoc)/pwcp(:,:,jwno3)) 397 !-------------------------------------------------------------------- 398 399 zrearat1(:,:) = 0. 400 zrearat2(:,:) = 0. 401 zrearat3(:,:) = 0. 402 403 zundsat(:,:) = pwcp(:,:,jwso4) 217 ! !-------------------------------------------------------------------- 218 ! ! Begining POC denitrification and NO3- diffusion 219 ! ! (indice n�5 for couple POC/NO3- ie solcp(:,:,jspoc)/pwcp(:,:,jwno3)) 220 ! !-------------------------------------------------------------------- 221 ! 222 zundsat(:,:) = MAX(0., pwcp(:,:,jwso4) - rtrn) 404 223 405 224 DO jk = 2, jpksed 406 225 DO ji = 1, jpoce 407 zlimso4(ji,jk) = ( 1.0 - pwcp(ji,jk,jwoxy) * zlimo2(ji,jk) ) * ( 1.0 - pwcp(ji,jk,jwno3) & 408 & / ( pwcp(ji,jk,jwno3) + xksedno3 ) ) * ( 1. - solcp(ji,jk,jsfeo) & 409 & / ( solcp(ji,jk,jsfeo) + xksedfeo ) ) / ( zundsat(ji,jk) + xksedso4 ) 410 zsolid1 = zvolc(ji,jk,jspoc) * solcp(ji,jk,jspoc) 411 zsolid2 = zvolc(ji,jk,jspos) * solcp(ji,jk,jspos) 412 zsolid3 = zvolc(ji,jk,jspor) * solcp(ji,jk,jspor) 413 zkpoca = zkpoc(ji,jk) * zlimso4(ji,jk) 414 zkpocb = zkpos(ji,jk) * zlimso4(ji,jk) 415 zkpocc = zkpor(ji,jk) * zlimso4(ji,jk) 416 zrearat1(ji,jk) = ( zkpoc(ji,jk) * dtsed2 * zsolid1 ) / & 417 & ( 1. + zkpoca * zundsat(ji,jk ) * dtsed2 ) 418 zrearat2(ji,jk) = ( zkpos(ji,jk) * dtsed2 * zsolid2 ) / & 419 & ( 1. + zkpocb * zundsat(ji,jk ) * dtsed2 ) 420 zrearat3(ji,jk) = ( zkpor(ji,jk) * dtsed2 * zsolid3 ) / & 421 & ( 1. + zkpocc * zundsat(ji,jk ) * dtsed2 ) 422 END DO 423 END DO 424 ! 425 ! DO jn = 1, 5 426 DO jk = 2, jpksed 427 DO ji = 1, jpoce 428 jflag4: DO jn = 1, 10 429 zlimtmp = ( 1.0 - pwcp(ji,jk,jwoxy) * zlimo2(ji,jk) ) * ( 1.0 - pwcp(ji,jk,jwno3) & 430 & / ( pwcp(ji,jk,jwno3) + xksedno3 ) ) * ( 1. - solcp(ji,jk,jsfeo) & 431 & / ( solcp(ji,jk,jsfeo) + xksedfeo ) ) 432 zsolid1 = zvolc(ji,jk,jspoc) * solcp(ji,jk,jspoc) 433 zsolid2 = zvolc(ji,jk,jspos) * solcp(ji,jk,jspos) 434 zsolid3 = zvolc(ji,jk,jspor) * solcp(ji,jk,jspor) 435 zreasat = ( zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) ) 436 zbeta = xksedso4 - pwcp(ji,jk,jwso4) + 0.5 * zreasat * zlimtmp 437 zgamma = - xksedso4 * pwcp(ji,jk,jwso4) 438 zundsat2 = zundsat(ji,jk) 439 zundsat(ji,jk) = ( - zbeta + SQRT( zbeta**2 - 4.0 * zgamma ) ) / 2.0 440 zlimso4(ji,jk) = ( 1.0 - pwcp(ji,jk,jwoxy) * zlimo2(ji,jk) ) * ( 1.0 - pwcp(ji,jk,jwno3) & 441 & / ( pwcp(ji,jk,jwno3) + xksedno3 ) ) * ( 1. - solcp(ji,jk,jsfeo) & 442 & / ( solcp(ji,jk,jsfeo) + xksedfeo ) ) / ( zundsat(ji,jk) + xksedso4 ) 443 zkpoca = zkpoc(ji,jk) * zlimso4(ji,jk) 444 zkpocb = zkpos(ji,jk) * zlimso4(ji,jk) 445 zkpocc = zkpor(ji,jk) * zlimso4(ji,jk) 446 zrearat1(ji,jk) = ( zkpoc(ji,jk) * dtsed2 * zsolid1 ) / & 447 & ( 1. + zkpoca * zundsat(ji,jk ) * dtsed2 ) 448 zrearat2(ji,jk) = ( zkpos(ji,jk) * dtsed2 * zsolid2 ) / & 449 & ( 1. + zkpocb * zundsat(ji,jk ) * dtsed2 ) 450 zrearat3(ji,jk) = ( zkpor(ji,jk) * dtsed2 * zsolid3 ) / & 451 & ( 1. + zkpocc * zundsat(ji,jk ) * dtsed2 ) 452 IF ( ABS( (zundsat(ji,jk)-zundsat2)/(zundsat2+rtrn)) < 1E-8 ) THEN 453 EXIT jflag4 454 ENDIF 455 END DO jflag4 456 END DO 457 END DO 458 459 ! New solid concentration values (jk=2 to jksed) for each couple 460 DO jk = 2, jpksed 461 DO ji = 1, jpoce 462 zreasat = zrearat1(ji,jk) * zlimso4(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspoc) 463 solcp(ji,jk,jspoc) = solcp(ji,jk,jspoc) - zreasat 464 zreasat = zrearat2(ji,jk) * zlimso4(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspos) 465 solcp(ji,jk,jspos) = solcp(ji,jk,jspos) - zreasat 466 zreasat = zrearat3(ji,jk) * zlimso4(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspor) 467 solcp(ji,jk,jspor) = solcp(ji,jk,jspor) - zreasat 468 ENDDO 469 ENDDO 470 ! 471 ! New dissolved concentrations 472 DO jk = 2, jpksed 473 DO ji = 1, jpoce 226 zlimtmp = 1.0 - zlimo2(ji,jk) - zlimno3(ji,jk) - zlimfeo(ji,jk) 227 zreasat = zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) 228 zbeta = xksedso4 - pwcp(ji,jk,jwso4) + 0.5 * zreasat * zlimtmp 229 zgamma = - xksedso4 * pwcp(ji,jk,jwso4) 230 zundsat(ji,jk) = ( - zbeta + SQRT( zbeta**2 - 4.0 * zgamma ) ) / 2.0 231 zlimso4(ji,jk) = zlimtmp * zundsat(ji,jk) / ( zundsat(ji,jk) + xksedso4 ) 232 474 233 ! For sulfur 475 pwcp(ji,jk,jwh2s) = pwcp(ji,jk,jwh2s) - ( zundsat(ji,jk) - pwcp(ji,jk,jwso4) ) 476 pwcp(ji,jk,jwso4) = zundsat(ji,jk) 477 zreasat = ( zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) ) * zlimso4(ji,jk) * zundsat(ji,jk) 478 ! For DIC 479 pwcp(ji,jk,jwdic) = pwcp(ji,jk,jwdic) + zreasat 480 zsumtot(ji) = zsumtot(ji) + zreasat / dtsed2 * volw3d(ji,jk) * 1.e-3 * 86400. * 365. * 1E3 481 ! For Phosphate (in mol/l) 482 pwcp(ji,jk,jwpo4) = pwcp(ji,jk,jwpo4) + zreasat * spo4r 483 ! Ligands 484 sedligand(ji,jk) = sedligand(ji,jk) + ratligc * zreasat 485 ! For iron (in mol/l) 486 pwcp(ji,jk,jwfe2) = pwcp(ji,jk,jwfe2) + fecratio(ji) * zreasat 234 pwcp(ji,jk,jwso4) = zundsat(ji,jk) + MIN(pwcp(ji,jk,jwso4), rtrn) 235 zreasat = zreasat * zlimso4(ji,jk) 236 pwcp(ji,jk,jwh2s) = pwcp(ji,jk,jwh2s) + 0.5 * zreasat 237 487 238 ! For alkalinity 488 pwcp(ji,jk,jwalk) = pwcp(ji,jk,jwalk) + zreasat * ( srno3 * zadsnh4 - 2.* spo4r ) + zreasat 489 ! Ammonium 490 pwcp(ji,jk,jwnh4) = pwcp(ji,jk,jwnh4) + zreasat * srno3 * zadsnh4 239 pwcp(ji,jk,jwalk) = pwcp(ji,jk,jwalk) + zreasat 491 240 ENDDO 492 241 ENDDO … … 496 245 ! Patankar scheme 497 246 ! ------------------------------------------------------------------------------ 498 499 247 call sed_dsr_redoxb 500 248 … … 564 312 !! Arguments 565 313 ! --- local variables 566 INTEGER :: ji, jk, jn ! dummy looop indices314 INTEGER :: ji, jk, jn, jw ! dummy looop indices 567 315 568 316 REAL, DIMENSION(6) :: zsedtrn, zsedtra … … 582 330 zsedtrn(5) = solcp(ji,jk,jsfeo) * zvolc(ji,jk,jsfeo) 583 331 zsedtrn(6) = solcp(ji,jk,jsfes) * zvolc(ji,jk,jsfes) 584 zsedtra(:) = zsedtrn(:)332 zsedtra(:) = MAX(0., zsedtrn(:) - rtrn ) 585 333 586 334 ! First pass of the scheme. At the end, it is 1st order … … 592 340 zdelta = pwcp(ji,jk,jwpo4) - redfep * zsedtra(4) 593 341 IF ( zalpha == 0. ) THEN 594 zsedtra(4) = zsedtra(4) / ( 1.0 + zsedtra(4) * reac_fe2 * dtsed2 / 2.0 ) 595 ELSE 596 zsedtra(4) = ( zsedtra(4) * zalpha ) / ( 0.25 * zsedtra(4) * ( exp( reac_fe2 * zalpha * dtsed2 / 2. ) - 1.0 ) & 342 zsedtra(4) = zsedtra(4) / ( 1.0 + 0.25 * zsedtra(4) * reac_fe2 * dtsed2 / 2.0 ) 343 ELSE 344 zsedtra(4) = ( zsedtra(4) * zalpha ) / ( 0.25 * zsedtra(4) & 345 & * ( exp( reac_fe2 * zalpha * dtsed2 / 2. ) - 1.0 ) & 597 346 & + zalpha * exp( reac_fe2 * zalpha * dtsed2 / 2. ) ) 598 347 ENDIF 599 zsedtra(1) = zalpha + 0.25 * zsedtra(4) 348 zsedtra(1) = zalpha + 0.25 * zsedtra(4) 600 349 zsedtra(5) = zbeta - zsedtra(4) 601 350 pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(4) … … 606 355 zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(2) 607 356 IF ( zalpha == 0. ) THEN 608 zsedtra(2) = zsedtra(2) / ( 1.0 + zsedtra(2) * reac_h2s * dtsed2 / 2.0 )357 zsedtra(2) = zsedtra(2) / ( 1.0 + 2.0 * zsedtra(2) * reac_h2s * dtsed2 / 2.0 ) 609 358 ELSE 610 359 zsedtra(2) = ( zsedtra(2) * zalpha ) / ( 2.0 * zsedtra(2) * ( exp( reac_h2s * zalpha * dtsed2 / 2. ) - 1.0 ) & … … 615 364 pwcp(ji,jk,jwso4) = zbeta - zsedtra(2) 616 365 ! NH4 + O2 617 zalpha = zsedtra(1) - 2.0 * zsedtra(3) 618 zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(3) 619 IF ( zalpha == 0. ) THEN 620 zsedtra(3) = zsedtra(3) / ( 1.0 + zsedtra(3) * reac_nh4 *zadsnh4 * dtsed2 / 2.0 )621 ELSE 622 zsedtra(3) = ( zsedtra(3) * zalpha ) / ( 2.0 * zsedtra(3) * ( exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 / 2. ) - 1.0 ) &623 & + zalpha * exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 /2. ) )624 ENDIF 625 zsedtra(1) = zalpha + 2.0 * zsedtra(3) 626 pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(3) 366 zalpha = zsedtra(1) - 2.0 * zsedtra(3) / zadsnh4 367 zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(3) /zadsnh4 368 IF ( zalpha == 0. ) THEN 369 zsedtra(3) = zsedtra(3) / ( 1.0 + 2.0 * zsedtra(3) * reac_nh4 / zadsnh4 * dtsed2 / 2.0 ) 370 ELSE 371 zsedtra(3) = ( zsedtra(3) * zalpha * zadsnh4 ) / ( 2.0 * zsedtra(3) * ( exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 / 2. ) - 1.0 ) & 372 & + zalpha * zadsnh4 * exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 /2. ) ) 373 ENDIF 374 zsedtra(1) = zalpha + 2.0 * zsedtra(3) / zadsnh4 375 pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(3) / zadsnh4 627 376 ! FeS - O2 628 377 zalpha = zsedtra(1) - 2.0 * zsedtra(6) … … 630 379 zgamma = pwcp(ji,jk,jwso4) + zsedtra(6) 631 380 IF ( zalpha == 0. ) THEN 632 zsedtra(6) = zsedtra(6) / ( 1.0 + zsedtra(6) * reac_feso * dtsed2 / 2.0 )381 zsedtra(6) = zsedtra(6) / ( 1.0 + 2.0 * zsedtra(6) * reac_feso * dtsed2 / 2.0 ) 633 382 ELSE 634 383 zsedtra(6) = ( zsedtra(6) * zalpha ) / ( 2.0 * zsedtra(6) * ( exp( reac_feso * zalpha * dtsed2 / 2. ) - 1.0 ) & … … 658 407 zepsi = pwcp(ji,jk,jwpo4) + redfep * zsedtra(5) 659 408 IF ( zalpha == 0. ) THEN 660 zsedtra(2) = zsedtra(2) / ( 1.0 + zsedtra(2) * reac_feh2s * dtsed2 )409 zsedtra(2) = zsedtra(2) / ( 1.0 + 2.0 * zsedtra(2) * reac_feh2s * dtsed2 ) 661 410 ELSE 662 411 zsedtra(2) = ( zsedtra(2) * zalpha ) / ( 2.0 * zsedtra(2) * ( exp( reac_feh2s * zalpha * dtsed2 ) - 1.0 ) & … … 668 417 pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(4) 669 418 pwcp(ji,jk,jwpo4) = zepsi - redfep * zsedtra(5) 670 419 ! Fe - H2S 671 420 zalpha = zsedtra(2) - zsedtra(4) 672 421 zbeta = zsedtra(4) + zsedtra(6) … … 686 435 zgamma = pwcp(ji,jk,jwso4) + zsedtra(6) 687 436 IF (zalpha == 0.) THEN 688 zsedtra(6) = zsedtra(6) / ( 1.0 + zsedtra(6) * reac_feso * dtsed2 / 2. )437 zsedtra(6) = zsedtra(6) / ( 1.0 + 2.0 * zsedtra(6) * reac_feso * dtsed2 / 2. ) 689 438 ELSE 690 439 zsedtra(6) = ( zsedtra(6) * zalpha ) / ( 2.0 * zsedtra(6) * ( exp( reac_feso * zalpha * dtsed2 / 2. ) - 1.0 ) & … … 694 443 zsedtra(4) = zbeta - zsedtra(6) 695 444 pwcp(ji,jk,jwso4) = zgamma - zsedtra(6) 445 696 446 ! NH4 + O2 697 zalpha = zsedtra(1) - 2.0 * zsedtra(3) 698 zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(3) 699 IF ( zalpha == 0.) THEN700 zsedtra(3) = zsedtra(3) / ( 1.0 + zsedtra(3) * reac_nh4 * zadsnh4 * dtsed2 / 2.0)701 ELSE 702 zsedtra(3) = ( zsedtra(3) * zalpha ) / ( 2.0 * zsedtra(3) * ( exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 / 2. ) - 1.0 ) &703 & + zalpha * exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 /2. ) )704 ENDIF 705 zsedtra(1) = zalpha + 2.0 * zsedtra(3) 706 pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(3) 447 zalpha = zsedtra(1) - 2.0 * zsedtra(3) / zadsnh4 448 zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(3) / zadsnh4 449 IF ( zalpha == 0. ) THEN 450 zsedtra(3) = zsedtra(3) / ( 1.0 + 2.0 * zsedtra(3) * reac_nh4 / zadsnh4 * dtsed2 / 2.0 ) 451 ELSE 452 zsedtra(3) = ( zsedtra(3) * zalpha * zadsnh4 ) / ( 2.0 * zsedtra(3) * ( exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 / 2. ) - 1.0 ) & 453 & + zalpha * zadsnh4 * exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 /2. ) ) 454 ENDIF 455 zsedtra(1) = zalpha + 2.0 * zsedtra(3) / zadsnh4 456 pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(3) / zadsnh4 707 457 ! H2S + O2 708 458 zalpha = zsedtra(1) - 2.0 * zsedtra(2) … … 710 460 zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(2) 711 461 IF ( zalpha == 0. ) THEN 712 zsedtra(2) = zsedtra(2) / ( 1.0 + zsedtra(2) * reac_h2s * dtsed2 / 2.0 )462 zsedtra(2) = zsedtra(2) / ( 1.0 + 2.0 * zsedtra(2) * reac_h2s * dtsed2 / 2.0 ) 713 463 ELSE 714 464 zsedtra(2) = ( zsedtra(2) * zalpha ) / ( 2.0 * zsedtra(2) * ( exp( reac_h2s * zalpha * dtsed2 / 2. ) - 1.0 ) & … … 718 468 pwcp(ji,jk,jwso4) = zbeta - zsedtra(2) 719 469 pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(2) 470 720 471 ! Fe + O2 721 472 zalpha = zsedtra(1) - 0.25 * zsedtra(4) … … 724 475 zdelta = pwcp(ji,jk,jwpo4) - redfep * zsedtra(4) 725 476 IF ( zalpha == 0. ) THEN 726 zsedtra(4) = zsedtra(4) / ( 1.0 + zsedtra(4) * reac_fe2 * dtsed2 / 2.0 ) 727 ELSE 728 zsedtra(4) = ( zsedtra(4) * zalpha ) / ( 0.25 * zsedtra(4) * ( exp( reac_fe2 * zalpha * dtsed2 / 2. ) - 1.0 ) & 477 zsedtra(4) = zsedtra(4) / ( 1.0 + 0.25 * zsedtra(4) * reac_fe2 * dtsed2 / 2.0 ) 478 ELSE 479 zsedtra(4) = ( zsedtra(4) * zalpha ) / ( 0.25 * zsedtra(4) & 480 & * ( exp( reac_fe2 * zalpha * dtsed2 / 2. ) - 1.0 ) & 729 481 & + zalpha * exp( reac_fe2 * zalpha * dtsed2 / 2. ) ) 730 482 ENDIF … … 733 485 pwcp(ji,jk,jwpo4) = zdelta + redfep * zsedtra(4) 734 486 pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(4) 487 488 ! Update the concentrations after the secondary reactions 489 zsedtra(:) = zsedtra(:) + MIN( zsedtrn(:), rtrn ) 735 490 pwcp(ji,jk,jwoxy) = zsedtra(1) 736 491 pwcp(ji,jk,jwh2s) = zsedtra(2) -
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/seddta.F90
r13295 r14385 53 53 54 54 REAL(wp), DIMENSION(jpoce) :: zdtap, zdtag 55 REAL(wp), DIMENSION(jpi,jpj) :: zwsbio4, zwsbio3 55 REAL(wp), DIMENSION(jpi,jpj) :: zwsbio4, zwsbio3, zddust 56 56 REAL(wp) :: zf0, zf1, zf2, zkapp, zratio, zdep 57 57 … … 80 80 ! conv2 = 1.0e+3 / ( 1.0e+4 * rsecday * 30. ) 81 81 conv2 = 1.0e+3 / 1.0e+4 82 rdtsed(2:jpksed) = dtsed / ( denssol * por1(2:jpksed) )83 82 ENDIF 84 83 … … 86 85 zdtap(:) = 0. 87 86 zdtag(:) = 0. 87 zddust(:,:) = 0.0 88 88 89 89 ! reading variables … … 156 156 ! zf1 = MIN(0.98, MAX(0., zf1 ) ) 157 157 zf1 = 0.48 158 zf0 = 1.0 - 0.02 - zf1159 158 zf2 = 0.02 159 zf0 = 1.0 - zf2 - zf1 160 160 rainrm_dta(ji,jspoc) = ( zdtap(ji) + zdtag(ji) ) * 1e-4 * zf0 161 161 rainrm_dta(ji,jspos) = ( zdtap(ji) + zdtag(ji) ) * 1e-4 * zf1 162 162 rainrm_dta(ji,jspor) = ( zdtap(ji) + zdtag(ji) ) * 1e-4 * zf2 163 163 END DO 164 164 165 ! Sinking fluxes for Calcite in mol.m-2.s-1 ; conversion in mol.cm-2.s-1 165 166 CALL pack_arr ( jpoce, rainrm_dta(1:jpoce,jscal), trc_data(1:jpi,1:jpj,14), iarroce(1:jpoce) ) … … 175 176 rainrm_dta(1:jpoce,jsclay) = rainrm_dta(1:jpoce,jsclay) * conv2 / mol_wgt(jsclay) & 176 177 & + wacc(1:jpoce) * por1(2) * denssol / mol_wgt(jsclay) / ( rsecday * 365.0 ) 177 rainrm_dta(1:jpoce,jsclay) = rainrm_dta(1:jpoce,jsclay) * 0.965 178 rainrm_dta(1:jpoce,jsfeo) = rainrm_dta(1:jpoce,jsclay) * mol_wgt(jsclay) / mol_wgt(jsfeo) * 0.035 / 0.965 178 rainrm_dta(1:jpoce,jsfeo) = rainrm_dta(1:jpoce,jsclay) * mol_wgt(jsclay) / mol_wgt(jsfeo) * 0.035 * 0.5 * 0.333 179 rainrm_dta(1:jpoce,jsclay) = rainrm_dta(1:jpoce,jsclay) * (1.0 - 0.035 * 0.5 * 0.333 ) 180 CALL unpack_arr ( jpoce, zddust(1:jpi,1:jpj), iarroce(1:jpoce), wacc(1:jpoce) ) 181 zddust(:,:) = dust(:,:) + zddust(:,:) / ( rsecday * 365.0 ) * por1(2) * denssol / conv2 182 179 183 ! rainrm_dta(1:jpoce,jsclay) = 1.0E-4 * conv2 / mol_wgt(jsclay) 180 184 … … 209 213 210 214 ! computation of dzdep = total thickness of solid material rained [cm] in each cell 211 dzdep(1:jpoce) = raintg(1:jpoce) * rdtsed(2)215 dzdep(1:jpoce) = raintg(1:jpoce) * dtsed / ( denssol * por1(2) ) 212 216 213 217 IF( lk_iomput ) THEN 214 IF( iom_use("sflxclay" ) ) CALL iom_put( "sflxclay", dust(:,:) * conv2 * 1E4 )215 IF( iom_use("sflxcal" ) ) CALL iom_put( "sflxcal", trc_data(:,:,1 3))216 IF( iom_use("sflxbsi" ) ) CALL iom_put( "sflxbsi", trc_data(:,:,1 0))217 IF( iom_use("sflxpoc" ) ) CALL iom_put( "sflxpoc", trc_data(:,:,11) + trc_data(:,:,12))218 IF( iom_use("sflxclay" ) ) CALL iom_put( "sflxclay", zddust(:,:) * 1E3 / 1.E4 ) 219 IF( iom_use("sflxcal" ) ) CALL iom_put( "sflxcal", trc_data(:,:,14) / 1.E4 ) 220 IF( iom_use("sflxbsi" ) ) CALL iom_put( "sflxbsi", trc_data(:,:,11) / 1.E4 ) 221 IF( iom_use("sflxpoc" ) ) CALL iom_put( "sflxpoc", ( trc_data(:,:,12) + trc_data(:,:,13) ) / 1.E4 ) 218 222 ENDIF 219 223 -
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/sedini.F90
r14086 r14385 43 43 44 44 REAL(wp), PUBLIC :: & 45 redO2 = 1 72. , & !: Redfield coef for Oxygen45 redO2 = 133. , & !: Redfield coef for Oxygen 46 46 redNo3 = 16. , & !: Redfield coef for Nitrate 47 47 redPo4 = 1. , & !: Redfield coef for Phosphate … … 73 73 74 74 REAL(wp), DIMENSION(jpwat), PUBLIC :: diff1 75 DATA diff1/4.59E-6, 1.104E-5, 4.81E-6 , 9.78E-6, 3.58E-6, 4. 01E-6, 9.8E-6, 9.73E-6, 5.0E-6, 3.31E-6 /75 DATA diff1/4.59E-6, 1.104E-5, 4.81E-6 , 9.78E-6, 3.58E-6, 4.81E-6, 9.8E-6, 9.73E-6, 5.0E-6, 3.31E-6 / 76 76 77 77 REAL(wp), DIMENSION(jpwat), PUBLIC :: diff2 78 DATA diff2/1.74E-7, 4.47E-7, 2.51E-7, 3.89E-7, 1.77E-7, 2.5 E-7, 3.89E-7, 3.06E-7, 2.5E-7, 1.5E-7 /78 DATA diff2/1.74E-7, 4.47E-7, 2.51E-7, 3.89E-7, 1.77E-7, 2.51E-7, 3.89E-7, 3.06E-7, 2.5E-7, 1.5E-7 / 79 79 80 80 … … 162 162 ALLOCATE( diff(jpoce,jpksed,jpwat ) ) ; ALLOCATE( irrig(jpoce, jpksed) ) 163 163 ALLOCATE( wacc(jpoce) ) ; ALLOCATE( fecratio(jpoce) ) 164 ALLOCATE( press(jpoce) ) ; ALLOCATE( densSW(jpoce) )165 164 ALLOCATE( hipor(jpoce,jpksed) ) ; ALLOCATE( co3por(jpoce,jpksed) ) 166 165 ALLOCATE( dz3d(jpoce,jpksed) ) ; ALLOCATE( volw3d(jpoce,jpksed) ) ; ALLOCATE( vols3d(jpoce,jpksed) ) 167 ALLOCATE( sedligand(jpoce, jpksed) ) 166 ALLOCATE( sedligand(jpoce, jpksed) ) ; ALLOCATE( saturco3(jpoce,jpksed) ) ; ALLOCATE( densSW(jpoce) ) 168 167 169 168 ! Initialization of global variables 170 pwcp (:,:,:) = 0. ; pwcp0 (:,:,:) = 0. ; pwcp_dta (:,:) = 0.171 solcp (:,:,:) = 0. ; solcp0 (:,:,:) = 0. ; rainrm_dta(:,:) = 0.172 rainrm(:,: ) = 0. ; rainrg (:,: ) = 0. ; raintg (: ) = 0.173 dzdep (: ) = 0. ; iarroce(: ) = 0 ; dzkbot (: ) = 0.174 temp (: ) = 0. ; salt (: ) = 0. ; zkbot (: ) = 0.175 press (: ) = 0. ; densSW (: ) = 0. ; db (:,:)= 0.176 hipor (:,: ) = 0. ; co3por (:,: ) = 0. ; irrig (:,:) = 0.177 dz3d (:,: ) = 0. ; volw3d (:,: ) = 0. ; vols3d (:,:) = 0.169 pwcp (:,:,:) = 0. ; pwcp0 (:,:,:) = 0. ; pwcp_dta (:,:) = 0. 170 solcp (:,:,:) = 0. ; solcp0 (:,:,:) = 0. ; rainrm_dta(:,:) = 0. 171 rainrm(:,: ) = 0. ; rainrg (:,: ) = 0. ; raintg (: ) = 0. 172 dzdep (: ) = 0. ; iarroce(: ) = 0 ; dzkbot (: ) = 0. 173 temp (: ) = 0. ; salt (: ) = 0. ; zkbot (: ) = 0. 174 db (:,: ) = 0. ; densSW (: ) = 0. 175 hipor (:,: ) = 0. ; co3por (:,: ) = 0. ; irrig (:,:) = 0. 176 dz3d (:,: ) = 0. ; volw3d (:,: ) = 0. ; vols3d (:,:) = 0. 178 177 fecratio(:) = 1E-5 179 178 sedligand(:,:) = 0.6E-9 180 179 181 180 ! Chemical variables 182 ALLOCATE( akbs (jpoce) ) ; ALLOCATE( ak1s (jpoce) ) ; ALLOCATE( ak2s (jpoce) ) ; ALLOCATE( akws (jpoce) ) 183 ALLOCATE( ak1ps (jpoce) ) ; ALLOCATE( ak2ps (jpoce) ) ; ALLOCATE( ak3ps (jpoce) ) ; ALLOCATE( aksis (jpoce) ) 184 ALLOCATE( aksps (jpoce) ) ; ALLOCATE( ak12s (jpoce) ) ; ALLOCATE( ak12ps(jpoce) ) ; ALLOCATE( ak123ps(jpoce) ) 185 ALLOCATE( borats(jpoce) ) ; ALLOCATE( calcon2(jpoce) ) ; ALLOCATE( sieqs (jpoce) ) 181 ALLOCATE( akbs (jpoce) ) ; ALLOCATE( ak1s (jpoce) ) ; ALLOCATE( ak2s (jpoce) ) ; ALLOCATE( akws (jpoce) ) 182 ALLOCATE( ak1ps (jpoce) ) ; ALLOCATE( ak2ps (jpoce) ) ; ALLOCATE( ak3ps (jpoce) ) ; ALLOCATE( aksis (jpoce) ) 183 ALLOCATE( aksps (jpoce) ) ; ALLOCATE( ak12s (jpoce) ) ; ALLOCATE( ak12ps(jpoce) ) ; ALLOCATE( ak123ps(jpoce) ) 184 ALLOCATE( borats(jpoce) ) ; ALLOCATE( calcon2(jpoce) ) ; ALLOCATE( sieqs (jpoce) ) 186 185 ALLOCATE( aks3s(jpoce) ) ; ALLOCATE( akf3s(jpoce) ) ; ALLOCATE( sulfats(jpoce) ) 187 186 ALLOCATE( fluorids(jpoce) ) … … 194 193 195 194 ! Mass balance calculation 196 ALLOCATE( fromsed(jpoce, jpsol) ) ; ALLOCATE( tosed(jpoce, jpsol) ) ; ALLOCATE( rloss(jpoce, jpsol) ) 197 ALLOCATE( tokbot (jpoce, jpwat) ) 198 199 fromsed(:,:) = 0. ; tosed(:,:) = 0. ; rloss(:,:) = 0. ; tokbot(:,:) = 0. 195 ALLOCATE( fromsed(jpoce, jpsol) ) ; ALLOCATE( tosed(jpoce, jpsol) ) 196 197 fromsed(:,:) = 0. ; tosed(:,:) = 0. 200 198 201 199 ! Initialization of sediment geometry … … 218 216 ! -------------------------------- 219 217 IF (ln_sediment_offline) THEN 220 CALL trc_dta_ini(jptra)221 218 CALL trc_dmp_sed_ini 222 219 ENDIF … … 370 367 denssol = 2.6 371 368 372 ! Initialization of diffusion coefficient as function of porosity [cm**2/s]373 !--------------------------------------------------------------------374 ! DO jn = 1, jpsol375 ! DO jk = 1, jpksed376 ! DO ji = 1, jpoce377 ! diff(ji,jk,jn) = dcoef / ( 1.0 - 2.0 * log(por(jk)) )378 ! END DO379 ! END DO380 ! END DO381 382 369 ! Accumulation rate from Burwicz et al. (2011). This is used to 383 370 ! compute the flux of clays and minerals … … 386 373 ztmp1 = 0.117 / ( 1.0 + ( zkbot(ji) / 200.)**3 ) 387 374 ztmp2 = 0.006 / ( 1.0 + ( zkbot(ji) / 4000.)**10 ) 388 wacc(ji) = ztmp1 + ztmp2 375 wacc(ji) = (ztmp1 + ztmp2)/10.0 376 wacc(ji) = ztmp2 / 10.0 389 377 END DO 390 378 -
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/sedinorg.F90
r12839 r14385 23 23 CONTAINS 24 24 25 SUBROUTINE sed_inorg( kt )25 SUBROUTINE sed_inorg( kt, knt ) 26 26 !!---------------------------------------------------------------------- 27 27 !! *** ROUTINE sed_inorg *** … … 46 46 !!---------------------------------------------------------------------- 47 47 !! Arguments 48 INTEGER, INTENT(in) :: kt ! number of iteration48 INTEGER, INTENT(in) :: kt, knt ! number of iteration 49 49 ! --- local variables 50 50 INTEGER :: ji, jk, js, jw ! dummy looop indices … … 54 54 REAL(wp), DIMENSION(jpoce,jpksed,jpsol) :: zvolc ! temp. variables 55 55 REAL(wp), DIMENSION(jpoce) :: zsieq 56 REAL(wp) :: zsolid1, zvolw, zreasat 56 REAL(wp) :: zsolid1, zvolw, zreasat, zbeta, zgamma 57 57 REAL(wp) :: zsatur, zsatur2, znusil, zsolcpcl, zsolcpsi 58 58 !! … … 61 61 IF( ln_timing ) CALL timing_start('sed_inorg') 62 62 ! 63 IF( kt == nitsed000 ) THEN63 IF( kt == nitsed000 .AND. knt == 1 ) THEN 64 64 IF (lwp) THEN 65 65 WRITE(numsed,*) ' sed_inorg : Dissolution reaction ' … … 73 73 74 74 zrearat1(:,:) = 0. ; zundsat(:,:) = 0. 75 zrearat2(:,:) = 0. ; zrearat2(:,:) = 0.75 zrearat2(:,:) = 0. 76 76 zco3eq(:) = rtrn 77 77 zvolc(:,:,:) = 0. … … 86 86 zsolcpsi = 0.0 87 87 DO jk = 1, jpksed 88 zsolcpsi = zsolcpsi + solcp(ji,jk,jsopal) * dz(jk)89 zsolcpcl = zsolcpcl + solcp(ji,jk,jsclay) * dz(jk)88 zsolcpsi = zsolcpsi + solcp(ji,jk,jsopal) * vols3d(ji,jk) 89 zsolcpcl = zsolcpcl + solcp(ji,jk,jsclay) * vols3d(ji,jk) 90 90 END DO 91 91 zsolcpsi = MAX( zsolcpsi, rtrn ) … … 135 135 DO jk = 2, jpksed 136 136 DO ji = 1, jpoce 137 zsolid1 = zvolc(ji,jk,jsopal) * solcp(ji,jk,jsopal) 138 zsatur = MAX(0., zundsat(ji,jk) / zsieq(ji) ) 139 zsatur2 = (1.0 + temp(ji) / 400.0 )**37 140 znusil = ( 0.225 * ( 1.0 + temp(ji) / 15.) + 0.775 * zsatur2 * zsatur**2.25 ) / zsieq(ji) 141 zrearat1(ji,jk) = ( reac_sil * znusil * dtsed * zsolid1 ) / & 142 & ( 1. + reac_sil * znusil * dtsed * zundsat(ji,jk) ) 143 ENDDO 144 ENDDO 145 146 CALL sed_mat( jwsil, jpoce, jpksed, zrearat1, zrearat2, zundsat, dtsed ) 147 148 ! New solid concentration values (jk=2 to jksed) for each couple 149 DO jk = 2, jpksed 150 DO ji = 1, jpoce 151 zreasat = zrearat1(ji,jk) * zundsat(ji,jk) / ( zvolc(ji,jk,jsopal) ) 152 solcp(ji,jk,jsopal) = solcp(ji,jk,jsopal) - zreasat 153 ENDDO 154 ENDDO 155 156 ! New pore water concentrations 157 DO jk = 1, jpksed 158 DO ji = 1, jpoce 159 ! Acid Silicic 160 pwcp(ji,jk,jwsil) = zsieq(ji) - zundsat(ji,jk) 137 IF ( zundsat(ji,jk) > 0. ) THEN 138 zsolid1 = zvolc(ji,jk,jsopal) * solcp(ji,jk,jsopal) 139 zsatur = zundsat(ji,jk) / zsieq(ji) 140 zsatur2 = (1.0 + temp(ji) / 400.0 )**37 141 znusil = ( 0.225 * ( 1.0 + temp(ji) / 15.) + 0.775 * zsatur2 * zsatur**8.25 ) / zsieq(ji) 142 zbeta = 1.0 - reac_sil * znusil * dtsed2 * zundsat(ji,jk) + reac_sil * znusil * dtsed2 * zsolid1 143 zgamma = - reac_sil * znusil * dtsed2 * zundsat(ji,jk) 144 zundsat(ji,jk) = ( - zbeta + SQRT( zbeta**2 - 4.0 * zgamma ) ) / ( 2.0 * reac_sil * znusil * dtsed2 ) 145 solcp(ji,jk,jsopal ) = zsolid1 / ( 1. + reac_sil * znusil * dtsed2 * zundsat(ji,jk) ) / zvolc(ji,jk,jsopal) 146 pwcp(ji,jk,jwsil) = zsieq(ji) - zundsat(ji,jk) 147 ENDIF 161 148 ENDDO 162 149 ENDDO … … 175 162 zco3eq(ji) = aksps(ji) * densSW(ji) * densSW(ji) / ( calcon2(ji) + rtrn ) 176 163 zco3eq(ji) = MAX( rtrn, zco3eq(ji) ) 177 zundsat(ji,jk) = MAX(0., zco3eq(ji) - co3por(ji,jk) ) 164 zundsat(ji,jk) = ( zco3eq(ji) - co3por(ji,jk) ) / zco3eq(ji) 165 saturco3(ji,jk) = zundsat(ji,jk) 178 166 ENDDO 179 167 ENDDO … … 181 169 DO jk = 2, jpksed 182 170 DO ji = 1, jpoce 183 zsolid1 = zvolc(ji,jk,jscal) * solcp(ji,jk,jscal) 184 zrearat1(ji,jk) = ( reac_cal * dtsed * zsolid1 / zco3eq(ji) ) / & 185 & ( 1. + reac_cal * dtsed * zundsat(ji,jk) / zco3eq(ji) ) 171 IF ( zundsat(ji,jk) > 0. ) THEN 172 zsolid1 = zvolc(ji,jk,jscal) * solcp(ji,jk,jscal) 173 zbeta = 1.0 - reac_cal * dtsed2 * zundsat(ji,jk) + reac_cal * dtsed2 * zsolid1 174 zgamma = -reac_cal * dtsed2 * zundsat(ji,jk) 175 zundsat(ji,jk) = ( -zbeta + SQRT( zbeta**2 - 4.0 * zgamma ) ) / ( 2.0 * reac_cal * dtsed2 ) 176 saturco3(ji,jk) = zundsat(ji,jk) 177 zreasat = reac_cal * dtsed2 * zundsat(ji,jk) * zsolid1 / ( 1. + reac_cal * dtsed2 * zundsat(ji,jk) ) 178 solcp(ji,jk,jscal) = solcp(ji,jk,jscal) - zreasat / zvolc(ji,jk,jscal) 179 ! For DIC 180 pwcp(ji,jk,jwdic) = pwcp(ji,jk,jwdic) + zreasat 181 ! For alkalinity 182 pwcp(ji,jk,jwalk) = pwcp(ji,jk,jwalk) + 2.0 * zreasat 183 ENDIF 186 184 END DO 187 185 END DO 188 189 ! solves tridiagonal system190 CALL sed_mat( jwdic, jpoce, jpksed, zrearat1, zrearat2, zundsat, dtsed )191 192 ! New solid concentration values (jk=2 to jksed) for cacO3193 DO jk = 2, jpksed194 DO ji = 1, jpoce195 zreasat = zrearat1(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jscal)196 solcp(ji,jk,jscal) = solcp(ji,jk,jscal) - zreasat197 ENDDO198 ENDDO199 200 ! New dissolved concentrations201 DO jk = 1, jpksed202 DO ji = 1, jpoce203 zreasat = zrearat1(ji,jk) * zundsat(ji,jk)204 ! For DIC205 pwcp(ji,jk,jwdic) = pwcp(ji,jk,jwdic) + zreasat206 ! For alkalinity207 pwcp(ji,jk,jwalk) = pwcp(ji,jk,jwalk) + 2.0 * zreasat208 ENDDO209 ENDDO210 211 !-------------------------------------------------212 ! Beginning DIC, Alkalinity213 !-------------------------------------------------214 215 DO jk = 1, jpksed216 DO ji = 1, jpoce217 zundsat(ji,jk) = pwcp(ji,jk,jwdic)218 zrearat1(ji,jk) = 0.219 ENDDO220 ENDDO221 222 ! solves tridiagonal system223 CALL sed_mat( jwdic, jpoce, jpksed, zrearat1, zrearat2, zundsat, dtsed )224 225 ! New dissolved concentrations226 DO jk = 1, jpksed227 DO ji = 1, jpoce228 pwcp(ji,jk,jwdic) = zundsat(ji,jk)229 ENDDO230 ENDDO231 232 !-------------------------------------------------233 ! Beginning DIC, Alkalinity234 !-------------------------------------------------235 236 DO jk = 1, jpksed237 DO ji = 1, jpoce238 zundsat(ji,jk) = pwcp(ji,jk,jwalk)239 zrearat1(ji,jk) = 0.240 ENDDO241 ENDDO242 !243 ! ! solves tridiagonal system244 CALL sed_mat( jwalk, jpoce, jpksed, zrearat1, zrearat2, zundsat, dtsed )245 !246 ! ! New dissolved concentrations247 DO jk = 1, jpksed248 DO ji = 1, jpoce249 pwcp(ji,jk,jwalk) = zundsat(ji,jk)250 ENDDO251 ENDDO252 253 !----------------------------------254 ! Back to initial geometry255 !-----------------------------256 257 !---------------------------------------------------------------------258 ! 1/ Compensation for ajustement of the bottom water concentrations259 ! (see note n° 1 about *por(2))260 !--------------------------------------------------------------------261 DO jw = 1, jpwat262 DO ji = 1, jpoce263 pwcp(ji,1,jw) = pwcp(ji,1,jw) + &264 & pwcp(ji,2,jw) * dzdep(ji) * por(2) / dzkbot(ji)265 END DO266 ENDDO267 268 !-----------------------------------------------------------------------269 ! 2/ Det of new rainrg taking account of the new weight fraction obtained270 ! in dz3d(2) after diffusion/reaction (react/diffu are also in dzdep!)271 ! This new rain (rgntg rm) will be used in advection/burial routine272 !------------------------------------------------------------------------273 DO js = 1, jpsol274 DO ji = 1, jpoce275 rainrg(ji,js) = raintg(ji) * solcp(ji,2,js)276 rainrm(ji,js) = rainrg(ji,js) / mol_wgt(js)277 END DO278 ENDDO279 280 ! New raintg281 raintg(:) = 0.282 DO js = 1, jpsol283 DO ji = 1, jpoce284 raintg(ji) = raintg(ji) + rainrg(ji,js)285 END DO286 ENDDO287 288 !--------------------------------289 ! 3/ back to initial geometry290 !--------------------------------291 DO ji = 1, jpoce292 dz3d (ji,2) = dz(2)293 volw3d(ji,2) = dz3d(ji,2) * por(2)294 vols3d(ji,2) = dz3d(ji,2) * por1(2)295 ENDDO296 186 297 187 IF( ln_timing ) CALL timing_stop('sed_inorg') -
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/sedstp.F90
r13970 r14385 8 8 USE sedchem ! chemical constant 9 9 USE sedco3 ! carbonate in sediment pore water 10 USE sedorg ! Organic reactions and diffusion 11 USE sedinorg ! Inorganic dissolution 10 USE sedsol ! Organic reactions and diffusion 12 11 USE sedbtb ! bioturbation 13 12 USE sedadv ! vertical advection 14 USE sedmbc ! mass balance calculation15 13 USE sedsfc ! sediment surface data 16 14 USE sedrst ! restart … … 46 44 INTEGER, INTENT(in) :: kt ! number of iteration 47 45 INTEGER, INTENT(in) :: Kbb, Kmm, Krhs ! time level indices 48 INTEGER :: ji,jk,js,jn,jw49 46 !!---------------------------------------------------------------------- 50 IF( ln_timing ) CALL timing_start('sed_stp')47 IF( ln_timing ) CALL timing_start('sed_stp') 51 48 ! 52 49 CALL sed_rst_opn ( kt ) ! Open tracer restart file … … 66 63 ENDIF 67 64 68 CALL sed_btb( kt ) ! 1st pass of bioturbation at t+1/2 69 CALL sed_org( kt ) ! Organic related reactions and diffusion 70 CALL sed_inorg( kt ) ! Dissolution reaction 71 CALL sed_btb( kt ) ! 2nd pass of bioturbation at t+1 72 tokbot(:,:) = 0.0 73 DO jw = 1, jpwat 74 DO ji = 1, jpoce 75 tokbot(ji,jw) = pwcp(ji,1,jw) * 1.e-3 * dzkbot(ji) 76 END DO 77 ENDDO 65 CALL sed_btb( kt ) ! 1st pass of bioturbation at t+1/2 66 CALL sed_sol( kt ) ! Solute diffusion and reactions 67 CALL sed_btb( kt ) ! 2nd pass of bioturbation at t+1 78 68 CALL sed_adv( kt ) ! advection 79 69 CALL sed_co3( kt ) ! pH actualization for saving 80 ! This routine is commented out since it does not work at all81 CALL sed_mbc( kt ) ! cumulation for mass balance calculation82 83 70 IF (ln_sed_2way) CALL sed_sfc( kt, Kbb ) ! Give back new bottom wat chem to tracer model 84 71 ENDIF -
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/sedwri.F90
r14086 r14385 54 54 IF (lwp) WRITE(numsed,*) ' ' 55 55 56 ALLOCATE( zdta(jpoce,jpksed) ) ; ALLOCATE( zflx(jpoce,jp watp1) )56 ALLOCATE( zdta(jpoce,jpksed) ) ; ALLOCATE( zflx(jpoce,jptrased+1) ) 57 57 58 58 ! Initialize variables … … 88 88 CALL unpack_arr( jpoce, flxsedi3d(1:jpi,1:jpj,1:jpksed,2) , iarroce(1:jpoce), & 89 89 & co3por(1:jpoce,1:jpksed) ) 90 91 CALL unpack_arr( jpoce, flxsedi3d(1:jpi,1:jpj,1:jpksed,3) , iarroce(1:jpoce), & 92 & sedligand(1:jpoce,1:jpksed) ) 93 94 CALL unpack_arr( jpoce, flxsedi3d(1:jpi,1:jpj,1:jpksed,4) , iarroce(1:jpoce), & 95 & saturco3(1:jpoce,1:jpksed) ) 96 90 97 91 98 ! flxsedi3d = 0. … … 95 102 DO ji = 1, jpoce 96 103 zflx(ji,jw) = ( pwcp(ji,1,jw) - pwcp_dta(ji,jw) ) & 97 & * 1.e3 / 1.e2 * dzkbot(ji) / rDt_trc 104 & * 1.e3 * ( 1.e-2 * dzkbot(ji) ) / 1.E4 / rDt_trc 105 ENDDO 106 ENDDO 107 108 ! Calculation of fluxes g/cm2/s 109 DO js = 1, jpsol 110 zrate = 1.0 / rDt_trc 111 DO ji = 1, jpoce 112 zflx(ji,jpwat+js) = zflx(ji,jpwat+js) + ( tosed(ji,js) - fromsed(ji,js) ) * zrate 98 113 ENDDO 99 114 ENDDO … … 101 116 ! Calculation of accumulation rate per dt 102 117 DO js = 1, jpsol 103 zrate = 1.0 / ( denssol * por1(jpksed) ) /rDt_trc118 zrate = 1.0 / rDt_trc 104 119 DO ji = 1, jpoce 105 zflx(ji,jp watp1) = zflx(ji,jpwatp1) + ( tosed(ji,js) - fromsed(ji,js) ) * zrate120 zflx(ji,jptrased+1) = zflx(ji,jptrased+1) + ( tosed(ji,js) - fromsed(ji,js) ) * zrate 106 121 ENDDO 107 122 ENDDO … … 110 125 CALL unpack_arr( jpoce, flxsedi2d(1:jpi,1:jpj,jn), iarroce(1:jpoce), zflx(1:jpoce,jn) ) 111 126 END DO 127 112 128 zflx(:,1) = dzdep(:) / dtsed 113 129 CALL unpack_arr( jpoce, flxsedi2d(1:jpi,1:jpj,jpdia2dsed), iarroce(1:jpoce), zflx(1:jpoce,1) ) 114 130 115 116 117 118 119 120 131 ! Start writing data 132 ! --------------------- 133 DO jn = 1, jptrased 134 cltra = sedtrcd(jn) ! short title for 3D diagnostic 135 CALL iom_put( cltra, trcsedi(:,:,:,jn) ) 136 END DO 121 137 122 123 124 125 138 DO jn = 1, jpdia3dsed 139 cltra = seddia3d(jn) ! short title for 3D diagnostic 140 CALL iom_put( cltra, flxsedi3d(:,:,:,jn) ) 141 END DO 126 142 127 128 129 130 143 DO jn = 1, jpdia2dsed 144 cltra = seddia2d(jn) ! short title for 2D diagnostic 145 CALL iom_put( cltra, flxsedi2d(:,:,jn) ) 146 END DO 131 147 132 148
Note: See TracChangeset
for help on using the changeset viewer.
