- Timestamp:
- 2021-01-21T00:33:21+01:00 (3 years ago)
- Location:
- NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA
- Files:
-
- 1 added
- 1 deleted
- 12 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/cfgs/SHARED/field_def_nemo-pisces.xml
r14276 r14323 137 137 <field id="SedpH" long_name="PH" unit="1" /> 138 138 <field id="SedCO3por" long_name="Bicarbonates" unit="mol/m3" /> 139 <field id="Sedligand" long_name="Ligands Concentration" unit="mol/m3" /> 139 <field id="Sedligand" long_name="Ligands" unit="mol/m3" /> 140 <field id="SaturCO3" long_name="CO3 Saturation" unit="-" /> 140 141 </field_group> 141 142 … … 153 154 <field id="FlxFe2" long_name="Fe2+ sediment flux" unit="mol/cm2/s" /> 154 155 <field id="dzdep" long_name="Sedimentation rate" unit="cm/s" /> 155 <field id="sflxclay" long_name="Clay sedimentation rate" unit="g/m2/s" /> 156 <field id="sflxcal" long_name="Calcite sedimentation rate" unit="mol/m2/s" /> 157 <field id="sflxbsi" long_name="BSi Sedimentation rate" unit="mol/m2/s" /> 158 <field id="sflxpoc" long_name="POC Sedimentation rate" unit="mol/m2/s" /> 159 <field id="sflxfeo" long_name="Fe(OH)3 Sedimentation rate" unit="mol/m2/s" /> 160 <field id="FlxBSi" long_name="BSi burial flux" unit="g/cm2/s" /> 161 <field id="FlxClay" long_name="Clay burial flux" unit="g/cm2/s" /> 162 <field id="FlxPOC" long_name="POC burial flux" unit="g/cm2/s" /> 163 <field id="FlxCaCO3" long_name="Calcite burial flux" unit="g/cm2/s" /> 164 <field id="FlxPOS" long_name="POS burial flux" unit="g/cm2/s" /> 165 <field id="FlxPOR" long_name="POR burial flux" unit="g/cm2/s" /> 166 <field id="FlxFeO" long_name="FeO burial flux" unit="g/cm2/s" /> 167 <field id="FlxFeS" long_name="FeS burial flux" unit="g/cm2/s" /> 156 <field id="sflxclay" long_name="Clay sedimentation rate" unit="g/cm2/s" /> 157 <field id="sflxbsi" long_name="BSi sedimentation rate" unit="g/cm2/s" /> 158 <field id="sflxpoc" long_name="POC sedimentation rate" unit="g/cm2/s" /> 159 <field id="sflxcal" long_name="Calcite sedimentation rate" unit="g/cm2/s" /> 160 <field id="FlxClay" long_name="Clay burial rate" unit="g/cm2/s" /> 161 <field id="FlxCaCO3" long_name="Calcite burial rate" unit="g/cm2/s" /> 162 <field id="FlxBSi" long_name="BSi burial rate" unit="g/cm2/s" /> 163 <field id="FlxPOC" long_name="POC burial rate" unit="g/cm2/s" /> 164 <field id="FlxFeO" long_name="Fe(OH)3 burial rate" unit="g/cm2/s" /> 165 <field id="FlxFeS" long_name="FeS burial rate" unit="g/cm2/s" /> 166 <field id="FlxPOR" long_name="POR burial rate" unit="g/cm2/s" /> 167 <field id="FlxPOS" long_name="POS burial rate" unit="g/cm2/s" /> 168 <field id="Flxtot" long_name="total burial flux" unit="g/cm2/s" /> 169 168 170 </field_group> 169 171 … … 215 217 <field id="SIZEP" long_name="Mean relative size of picophyto." unit="-" grid_ref="grid_T_3D" /> 216 218 <field id="SIZED" long_name="Mean relative size of diatoms" unit="-" grid_ref="grid_T_3D" /> 217 <field id="RASSD" long_name="Mean relative size of diatoms" unit="-" grid_ref="grid_T_3D" />218 <field id="RASSN" long_name="Mean relative size of diatoms" unit="-" grid_ref="grid_T_3D" />219 <field id="RASSP" long_name="Mean relative size of diatoms" unit="-" grid_ref="grid_T_3D" />220 219 <field id="Fe3" long_name="Iron III concentration" unit="nmol/m3" grid_ref="grid_T_3D" /> 221 220 <field id="FeL1" long_name="Complexed Iron concentration with L1" unit="nmol/m3" grid_ref="grid_T_3D" /> -
NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/cfgs/SHARED/namelist_sediment_ref
r14306 r14323 12 12 &nam_run ! Characteristics of the simulation 13 13 !----------------------------------------------------------------------- 14 nrseddt = 1! Nb of iterations for fast species14 nrseddt = 2 ! Nb of iterations for fast species 15 15 ln_sed_2way = .false. ! 2 way coupling with pisces 16 16 / … … 58 58 seddiag3d(1) = 'SedpH ' , 'pH ', '- ' 59 59 seddiag3d(2) = 'SedCO3por ' , 'Dissolved CO3 concentration ', 'mol/L' 60 seddiag3d(3) = 'Sedligand' , 'ligand concentration ', 'mol/L' 60 seddiag3d(3) = 'Sedligand ' , 'ligand concentration ', 'mol/L' 61 seddiag3d(4) = 'SaturCO3 ' , 'CO3 saturation ', '-' 61 62 seddiag2d(1) = 'FlxSi ' , 'Silicate flux ', 'mol/cm2/s' 62 63 seddiag2d(2) = 'FlxO2 ' , 'Dissolved Oxygen Flux ', 'mol/L' … … 99 100 rcorgl = 10. ! Reactivity for labile POC [an-1] 100 101 rcorgs = 0.1 ! Reactivity for semi-refractory POC [an-1] 101 rcorgr = 1.E-4 ! Reactivity for refractory POC [an-1]102 rcorgr = 5.E-4 ! Reactivity for refractory POC [an-1] 102 103 rcnh4 = 1E7 ! Reactivity for O2/NH4 [l.mol-1.an-1] 103 104 rch2s = 2E8 ! Reactivity for O2/H2S [l.mol-1.an-1] … … 106 107 rcfes = 1E6 ! Reactivity for FE2+/H2S [l.mol-1.an-1] 107 108 rcfeso = 2E7 ! Reactivity for FES/O2 [l.mol-1.an-1] 108 xksedo2 = 1.E-6 ! Half-saturation constant for oxic remin [mol/l]109 xksedo2 = 4.E-6 ! Half-saturation constant for oxic remin [mol/l] 109 110 xksedno3 = 10.E-6 ! Half-saturation constant for denitrification [mol/l] 110 111 xksedfeo = 0.007 ! Half-saturation constant for iron remin [%] … … 130 131 cn_sedrst_outdir = "." ! directory to which to write output sediment restarts 131 132 / 132 !-----------------------------------------------------------------------133 &nam_output ! parameters for outputing the sediment module134 !-----------------------------------------------------------------------135 ldefsedpis_avg = .true. ! write averaged output variables136 cn_sedwri_out = "output_sed.nc" ! name of the input restart file name of the sediment module137 nrpfsedpis_avg = 0 ! Frequency of the averaged outputs138 nwrtsedpis_avg = 24 ! Frequency of the averaged outputs139 ntssedpis_avg = 1 ! ???140 / -
NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/SED/par_sed.F90
r14276 r14323 61 61 INTEGER, PARAMETER :: & 62 62 jptrased = jpsol + jpwat , & 63 jpdia3dsed = 3, &63 jpdia3dsed = 4 , & 64 64 jpdia2dsed = 20 65 65 -
NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/SED/sed.F90
r10425 r14323 73 73 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: fromsed !: 74 74 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: tosed !: 75 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: rloss !:76 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: tokbot77 !78 75 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: temp !: temperature 79 76 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: salt !: salinity 80 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: press !: pressure81 77 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: raintg !: total massic flux rained in each cell (sum of sol. comp.) 82 78 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: fecratio !: Fe/C ratio in falling particles to the sediments 83 79 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: dzdep !: total thickness of solid material rained [cm] in each cell 84 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: zkbot !: total thickness of solid material rained [cm] in each cell80 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: zkbot !: Total water column depth 85 81 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: wacc !: total thickness of solid material rained [cm] in each cell 86 82 ! … … 91 87 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: vols3d !: ??? 92 88 93 94 89 !! Chemistry 95 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: densSW 96 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: borats 90 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: densSW 91 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: borats 97 92 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: calcon2 98 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: akbs 99 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: ak1s 100 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: ak2s 101 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: akws 102 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: ak12s 103 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: ak1ps 104 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: ak2ps 105 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: ak3ps 106 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: ak12ps 93 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: akbs 94 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: ak1s 95 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: ak2s 96 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: akws 97 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: ak12s 98 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: ak1ps 99 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: ak2ps 100 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: ak3ps 101 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: ak12ps 107 102 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: ak123ps 108 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: aksis 109 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: aksps 103 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: aksis 104 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: aksps 110 105 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: sieqs 111 106 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: aks3s … … 130 125 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: db !: bioturbation ceofficient 131 126 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: irrig !: bioturbation ceofficient 132 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: rdtsed !: sediment model time-step133 127 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: sedligand 128 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: saturco3 134 129 REAL(wp) :: dens !: density of solid material 135 130 !! Inputs / Outputs … … 157 152 !!------------------------------------------------------------------- 158 153 ! 159 ALLOCATE( trc_data(jpi,jpj,jpdta) , & 160 & epkbot(jpi,jpj), gdepbot(jpi,jpj) , & 161 & dz(jpksed) , por(jpksed) , por1(jpksed) , & 162 & volw(jpksed), vols(jpksed), rdtsed(jpksed) , & 163 & trcsedi (jpi,jpj,jpksed,jptrased) , & 164 & flxsedi3d(jpi,jpj,jpksed,jpdia3dsed) , & 165 & flxsedi2d(jpi,jpj,jpdia2dsed) , & 166 & mol_wgt(jpsol), STAT=sed_alloc ) 154 ALLOCATE( trc_data(jpi,jpj,jpdta), trcsedi (jpi,jpj,jpksed,jptrased) , & 155 & epkbot(jpi,jpj), gdepbot(jpi,jpj), dz(jpksed), por(jpksed) , & 156 & por1(jpksed), volw(jpksed), vols(jpksed), mol_wgt(jpsol) , & 157 & flxsedi2d(jpi,jpj,jpdia2dsed), flxsedi3d(jpi,jpj,jpksed,jpdia3dsed), STAT=sed_alloc ) 167 158 168 159 IF( sed_alloc /= 0 ) CALL ctl_stop( 'STOP', 'sed_alloc: failed to allocate arrays' ) -
NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/SED/sedadv.F90
r10425 r14323 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/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/SED/seddiff.F90
r10225 r14323 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/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/SED/seddsr.F90
r10362 r14323 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/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/SED/seddta.F90
r14306 r14323 49 49 50 50 REAL(wp), DIMENSION(jpoce) :: zdtap, zdtag 51 REAL(wp), DIMENSION(jpi,jpj) :: zwsbio4, zwsbio3 51 REAL(wp), DIMENSION(jpi,jpj) :: zwsbio4, zwsbio3, zddust 52 52 REAL(wp) :: zf0, zf1, zf2, zkapp, zratio, zdep 53 53 … … 76 76 ! conv2 = 1.0e+3 / ( 1.0e+4 * rsecday * 30. ) 77 77 conv2 = 1.0e+3 / 1.0e+4 78 rdtsed(2:jpksed) = dtsed / ( denssol * por1(2:jpksed) )79 78 ENDIF 80 79 … … 82 81 zdtap(:) = 0. 83 82 zdtag(:) = 0. 83 zddust(:,:) = 0.0 84 84 85 85 ! reading variables … … 157 157 ! zf1 = ( 0.02 * (reac_poc - reac_poc * 0.) + zkapp - reac_poc ) / ( reac_poc / 100. - reac_poc ) 158 158 ! zf1 = MIN(0.98, MAX(0., zf1 ) ) 159 zf1 = 0. 7159 zf1 = 0.48 160 160 zf2 = 0.02 161 161 zf0 = 1.0 - zf2 - zf1 … … 180 180 rainrm_dta(1:jpoce,jsfeo) = rainrm_dta(1:jpoce,jsclay) * mol_wgt(jsclay) / mol_wgt(jsfeo) * 0.035 * 0.5 * 0.333 181 181 rainrm_dta(1:jpoce,jsclay) = rainrm_dta(1:jpoce,jsclay) * (1.0 - 0.035 * 0.5 * 0.333 ) 182 CALL unpack_arr ( jpoce, zddust(1:jpi,1:jpj), iarroce(1:jpoce), wacc(1:jpoce) ) 183 zddust(:,:) = dust(:,:) + zddust(:,:) / ( rsecday * 365.0 ) * por1(2) * denssol / conv2 182 184 183 185 ! rainrm_dta(1:jpoce,jsclay) = 1.0E-4 * conv2 / mol_wgt(jsclay) … … 213 215 214 216 ! computation of dzdep = total thickness of solid material rained [cm] in each cell 215 dzdep(1:jpoce) = raintg(1:jpoce) * rdtsed(2)217 dzdep(1:jpoce) = raintg(1:jpoce) * dtsed / ( denssol * por1(2) ) 216 218 217 219 IF( lk_iomput ) THEN 218 IF( iom_use("sflxclay" ) ) CALL iom_put( "sflxclay", dust(:,:) * 1E3 / 1.E4 )220 IF( iom_use("sflxclay" ) ) CALL iom_put( "sflxclay", zddust(:,:) * 1E3 / 1.E4 ) 219 221 IF( iom_use("sflxcal" ) ) CALL iom_put( "sflxcal", trc_data(:,:,14) / 1.E4 ) 220 222 IF( iom_use("sflxbsi" ) ) CALL iom_put( "sflxbsi", trc_data(:,:,11) / 1.E4 ) -
NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/SED/sedini.F90
r14306 r14323 161 161 ALLOCATE( diff(jpoce,jpksed,jpwat ) ) ; ALLOCATE( irrig(jpoce, jpksed) ) 162 162 ALLOCATE( wacc(jpoce) ) ; ALLOCATE( fecratio(jpoce) ) 163 ALLOCATE( press(jpoce) ) ; ALLOCATE( densSW(jpoce) )164 163 ALLOCATE( hipor(jpoce,jpksed) ) ; ALLOCATE( co3por(jpoce,jpksed) ) 165 164 ALLOCATE( dz3d(jpoce,jpksed) ) ; ALLOCATE( volw3d(jpoce,jpksed) ) ; ALLOCATE( vols3d(jpoce,jpksed) ) 166 ALLOCATE( sedligand(jpoce, jpksed) ) 165 ALLOCATE( sedligand(jpoce, jpksed) ) ; ALLOCATE( saturco3(jpoce,jpksed) ) ; ALLOCATE( densSW(jpoce) ) 167 166 168 167 ! Initialization of global variables 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 press (: ) = 0. ; densSW (: ) = 0. ; db (:,:)= 0.175 hipor (:,: ) = 0. ; co3por (:,: ) = 0. ; irrig (:,:) = 0.176 dz3d (:,: ) = 0. ; volw3d (:,: ) = 0. ; vols3d (:,:) = 0.168 pwcp (:,:,:) = 0. ; pwcp0 (:,:,:) = 0. ; pwcp_dta (:,:) = 0. 169 solcp (:,:,:) = 0. ; solcp0 (:,:,:) = 0. ; rainrm_dta(:,:) = 0. 170 rainrm(:,: ) = 0. ; rainrg (:,: ) = 0. ; raintg (: ) = 0. 171 dzdep (: ) = 0. ; iarroce(: ) = 0 ; dzkbot (: ) = 0. 172 temp (: ) = 0. ; salt (: ) = 0. ; zkbot (: ) = 0. 173 db (:,: ) = 0. ; densSW (: ) = 0. 174 hipor (:,: ) = 0. ; co3por (:,: ) = 0. ; irrig (:,:) = 0. 175 dz3d (:,: ) = 0. ; volw3d (:,: ) = 0. ; vols3d (:,:) = 0. 177 176 fecratio(:) = 1E-5 178 177 sedligand(:,:) = 0.6E-9 179 178 180 179 ! Chemical variables 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) ) 180 ALLOCATE( akbs (jpoce) ) ; ALLOCATE( ak1s (jpoce) ) ; ALLOCATE( ak2s (jpoce) ) ; ALLOCATE( akws (jpoce) ) 181 ALLOCATE( ak1ps (jpoce) ) ; ALLOCATE( ak2ps (jpoce) ) ; ALLOCATE( ak3ps (jpoce) ) ; ALLOCATE( aksis (jpoce) ) 182 ALLOCATE( aksps (jpoce) ) ; ALLOCATE( ak12s (jpoce) ) ; ALLOCATE( ak12ps(jpoce) ) ; ALLOCATE( ak123ps(jpoce) ) 183 ALLOCATE( borats(jpoce) ) ; ALLOCATE( calcon2(jpoce) ) ; ALLOCATE( sieqs (jpoce) ) 185 184 ALLOCATE( aks3s(jpoce) ) ; ALLOCATE( akf3s(jpoce) ) ; ALLOCATE( sulfats(jpoce) ) 186 185 ALLOCATE( fluorids(jpoce) ) … … 193 192 194 193 ! Mass balance calculation 195 ALLOCATE( fromsed(jpoce, jpsol) ) ; ALLOCATE( tosed(jpoce, jpsol) ) ; ALLOCATE( rloss(jpoce, jpsol) ) 196 ALLOCATE( tokbot (jpoce, jpwat) ) 197 198 fromsed(:,:) = 0. ; tosed(:,:) = 0. ; rloss(:,:) = 0. ; tokbot(:,:) = 0. 194 ALLOCATE( fromsed(jpoce, jpsol) ) ; ALLOCATE( tosed(jpoce, jpsol) ) 195 196 fromsed(:,:) = 0. ; tosed(:,:) = 0. 199 197 200 198 ! Initialization of sediment geometry -
NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/SED/sedinorg.F90
r12360 r14323 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) 178 165 ENDDO 179 166 ENDDO … … 181 168 DO jk = 2, jpksed 182 169 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) ) 170 IF ( zundsat(ji,jk) > 0. ) THEN 171 zsolid1 = zvolc(ji,jk,jscal) * solcp(ji,jk,jscal) 172 zbeta = 1.0 - reac_cal * dtsed2 * zundsat(ji,jk) + reac_cal * dtsed2 * zsolid1 173 zgamma = -reac_cal * dtsed2 * zundsat(ji,jk) 174 zundsat(ji,jk) = ( -zbeta + SQRT( zbeta**2 - 4.0 * zgamma ) ) / ( 2.0 * reac_sil * znusil * dtsed2 ) 175 saturco3(ji,jk) = zundsat(ji,jk) 176 zreasat = reac_cal * dtsed2 * zundsat(ji,jk) * zsolid1 / ( 1. + reac_cal * dtsed2 * zundsat(ji,jk) ) 177 solcp(ji,jk,jscal) = solcp(ji,jk,jscal) - zreasat / zvolc(ji,jk,jscal) 178 ! For DIC 179 pwcp(ji,jk,jwdic) = pwcp(ji,jk,jwdic) + zreasat 180 ! For alkalinity 181 pwcp(ji,jk,jwalk) = pwcp(ji,jk,jwalk) + 2.0 * zreasat 182 ENDIF 186 183 END DO 187 184 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 185 297 186 IF( ln_timing ) CALL timing_stop('sed_inorg') -
NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/SED/sedstp.F90
r10222 r14323 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 … … 45 43 !!---------------------------------------------------------------------- 46 44 INTEGER, INTENT(in) :: kt ! number of iteration 47 INTEGER :: ji,jk,js,jn,jw48 45 !!---------------------------------------------------------------------- 49 IF( ln_timing ) CALL timing_start('sed_stp')46 IF( ln_timing ) CALL timing_start('sed_stp') 50 47 ! 51 48 CALL sed_rst_opn ( kt ) ! Open tracer restart file … … 65 62 ENDIF 66 63 67 CALL sed_btb( kt ) ! 1st pass of bioturbation at t+1/2 68 CALL sed_org( kt ) ! Organic related reactions and diffusion 69 CALL sed_inorg( kt ) ! Dissolution reaction 70 CALL sed_btb( kt ) ! 2nd pass of bioturbation at t+1 71 tokbot(:,:) = 0.0 72 DO jw = 1, jpwat 73 DO ji = 1, jpoce 74 tokbot(ji,jw) = pwcp(ji,1,jw) * 1.e-3 * dzkbot(ji) 75 END DO 76 ENDDO 64 CALL sed_btb( kt ) ! 1st pass of bioturbation at t+1/2 65 CALL sed_sol( kt ) ! Solute diffusion and reactions 66 CALL sed_btb( kt ) ! 2nd pass of bioturbation at t+1 77 67 CALL sed_adv( kt ) ! advection 78 68 CALL sed_co3( kt ) ! pH actualization for saving 79 ! This routine is commented out since it does not work at all80 CALL sed_mbc( kt ) ! cumulation for mass balance calculation81 82 69 IF (ln_sed_2way) CALL sed_sfc( kt ) ! Give back new bottom wat chem to tracer model 83 70 ENDIF -
NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/SED/sedwri.F90
r14306 r14323 90 90 CALL unpack_arr( jpoce, flxsedi3d(1:jpi,1:jpj,1:jpksed,3) , iarroce(1:jpoce), & 91 91 & sedligand(1:jpoce,1:jpksed) ) 92 93 CALL unpack_arr( jpoce, flxsedi3d(1:jpi,1:jpj,1:jpksed,4) , iarroce(1:jpoce), & 94 & saturco3(1:jpoce,1:jpksed) ) 95 92 96 93 97 ! flxsedi3d = 0. … … 97 101 DO ji = 1, jpoce 98 102 zflx(ji,jw) = ( pwcp(ji,1,jw) - pwcp_dta(ji,jw) ) & 99 & * 1.e3 / ( 1.e2 * dzkbot(ji) ) / 1.E4 / r2dttrc103 & * 1.e3 * ( 1.e-2 * dzkbot(ji) ) / 1.E4 / r2dttrc 100 104 ENDDO 101 105 ENDDO
Note: See TracChangeset
for help on using the changeset viewer.