Changeset 15450
- Timestamp:
- 2021-10-27T16:32:08+02:00 (2 years ago)
- Location:
- NEMO/trunk/src/TOP/PISCES
- Files:
-
- 5 added
- 1 deleted
- 20 edited
-
SED/par_sed.F90 (modified) (2 diffs)
-
SED/sed.F90 (modified) (10 diffs)
-
SED/sedadv.F90 (modified) (5 diffs)
-
SED/sedarr.F90 (modified) (2 diffs)
-
SED/sedbtb.F90 (modified) (3 diffs)
-
SED/sedchem.F90 (modified) (1 diff)
-
SED/seddiff.F90 (modified) (2 diffs)
-
SED/seddsr.F90 (modified) (7 diffs)
-
SED/seddsrjac.F90 (added)
-
SED/seddta.F90 (modified) (13 diffs)
-
SED/sedfunc.F90 (added)
-
SED/sedini.F90 (modified) (32 diffs)
-
SED/sedinitrc.F90 (modified) (5 diffs)
-
SED/sedinorg.F90 (modified) (4 diffs)
-
SED/sedjac.F90 (added)
-
SED/sedmat.F90 (modified) (5 diffs)
-
SED/sedmbc.F90 (deleted)
-
SED/sedorg.F90 (modified) (4 diffs)
-
SED/sedrst.F90 (modified) (13 diffs)
-
SED/sedsfc.F90 (modified) (3 diffs)
-
SED/sedsol.F90 (added)
-
SED/sedstp.F90 (modified) (5 diffs)
-
SED/sedwri.F90 (modified) (6 diffs)
-
SED/trcdmp_sed.F90 (modified) (1 diff)
-
SED/trosk.F90 (added)
-
trcini_pisces.F90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/TOP/PISCES/SED/par_sed.F90
r15033 r15450 18 18 jp_sal => jp_sal !: indice of salintity 19 19 20 INTEGER, P UBLIC, PARAMETER :: jpdta = 1720 INTEGER, PARAMETER :: jpdta = 18 21 21 22 22 ! Vertical sediment geometry 23 INTEGER, PUBLIC :: jpksed = 11 24 INTEGER, PUBLIC :: jpksedm1 23 INTEGER, PUBLIC :: & 24 jpksed = 11 , & 25 jpksedm1 = 10 25 26 26 27 ! sediment tracer species 27 INTEGER, P UBLIC, PARAMETER :: &28 INTEGER, PARAMETER :: & 28 29 jpsol = 8, & !: number of solid component 29 jpwat = 10, & !: number of pore water component 30 jpwat = 11, & !: number of pore water component 31 jpads = 2 , & !! number adsorbed species 30 32 jpwatp1 = jpwat +1, & 31 33 jpsol1 = jpsol - 1 … … 33 35 34 36 ! pore water components 35 INTEGER, PUBLIC, PARAMETER :: & 36 jwsil = 1, & !: silic acid 37 jwoxy = 2, & !: oxygen 38 jwdic = 3, & !: dissolved inorganic carbon 39 jwno3 = 4, & !: nitrate 40 jwpo4 = 5, & !: phosphate 41 jwalk = 6, & !: alkalinity 42 jwnh4 = 7, & !: Ammonium 43 jwh2s = 8, & !: Sulfate 44 jwso4 = 9, & !: H2S 45 jwfe2 = 10 !: Fe2+ 37 INTEGER, PARAMETER :: & 38 jwoxy = 1, & !: oxygen 39 jwno3 = 2, & !: nitrate 40 jwpo4 = 3, & !: phosphate 41 jwnh4 = 4, & !: Ammonium 42 jwh2s = 5, & !: Sulfate 43 jwso4 = 6, & !: H2S 44 jwfe2 = 7, & !: Fe2+ 45 jwalk = 8, & !: Alkalinity 46 jwlgw = 9, & !: Alkalinity 47 jwdic = 10, & !: DIC 48 jwsil = 11 !: Silicate 46 49 47 50 ! solid components 48 INTEGER, P UBLIC, PARAMETER :: &49 js opal = 1, & !: opal sediment50 js clay = 2, & !: clay51 INTEGER, PARAMETER :: & 52 jsfeo = 1, & !: iron hydroxides 53 jsfes = 2, & !: FeS 51 54 jspoc = 3, & !: organic carbon 52 js cal = 4, & !: calcite53 jspo s = 5, & !: semi-refPOC54 js por = 6, & !: refractory POC55 js feo = 7, & !: iron hydroxides56 js fes = 8 !: FeS55 jspos = 4, & !: semi-ref POC 56 jspor = 5, & !: refractory POC 57 jscal = 6, & !: Calcite 58 jsopal = 7, & !: Opal 59 jsclay = 8 !: clay 57 60 58 INTEGER, P UBLIC, PARAMETER :: &61 INTEGER, PARAMETER :: & 59 62 jptrased = jpsol + jpwat , & 60 jpdia3dsed = 2 , & 61 jpdia2dsed = 12 63 jpvode = jptrased - 11 , & 64 jpdia3dsed = 3 , & 65 jpdia2dsed = 22 62 66 63 67 END MODULE par_sed -
NEMO/trunk/src/TOP/PISCES/SED/sed.F90
r14086 r15450 7 7 !! ! 06-12 (C. Ethe) Orignal 8 8 !!---------------------------------------------------------------------- 9 USE par_sed 9 10 USE oce_sed 10 11 USE in_out_manager … … 27 28 REAL(wp), PUBLIC :: reac_fe2 !: reactivity of Fe2+ in [l.mol-1.s-1] 28 29 REAL(wp), PUBLIC :: reac_feh2s !: reactivity of Fe2+ in [l.mol-1.s-1] 29 REAL(wp), PUBLIC :: reac_fes !: reactivity of Fe with H2S in [l.mol-1.s-1]30 30 REAL(wp), PUBLIC :: reac_feso !: reactivity of FeS with O2 in [l.mol-1.s-1] 31 REAL(wp), PUBLIC :: reac_fesp !: precipitation of FeS [mol.l-1.s-1] 32 REAL(wp), PUBLIC :: reac_fesd !: Dissolution of FeS [s-1] 31 33 REAL(wp), PUBLIC :: reac_cal !: reactivity of cal in [l.mol-1.s-1] 32 34 REAL(wp), PUBLIC :: adsnh4 !: adsorption coefficient of NH4 35 REAL(wp), PUBLIC :: adsfe2 !: adsorption coefficient of Fe2 33 36 REAL(wp), PUBLIC :: ratligc !: C/L ratio in POC 34 37 REAL(wp), PUBLIC :: so2ut … … 37 40 REAL(wp), PUBLIC :: srDnit 38 41 REAL(wp), PUBLIC :: dtsed !: sedimentation time step 39 REAL(wp), PUBLIC :: dtsed2 !: sedimentation time step40 42 INTEGER , PUBLIC :: nitsed000 41 43 INTEGER , PUBLIC :: nitsedend 42 INTEGER, PUBLIC :: nrseddt43 REAL , PUBLIC :: sedmask44 REAL(wp), PUBLIC :: denssol !: density of solid material45 44 LOGICAL , PUBLIC :: lrst_sed !: logical to control the trc restart write 46 45 LOGICAL , PUBLIC :: ln_rst_sed = .TRUE. !: initialisation from a restart file or not … … 55 54 CHARACTER(len = 80) , PUBLIC :: cn_sedrst_out !: suffix of pass. tracer restart name (output) 56 55 CHARACTER(len = 256), PUBLIC :: cn_sedrst_outdir !: restart output directory 56 INTEGER, PUBLIC :: nrosorder !: order of the rosenbrock method 57 REAL(wp), PUBLIC :: rosatol !: Tolerance for absolute error 58 REAL(wp), PUBLIC :: rosrtol !: Tolerance for relative error 57 59 58 60 ! 59 REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: pwcp !: pore water sediment data at given time-step 60 REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: pwcp0 !: pore water sediment data at initial time 61 REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: solcp !: solid sediment data at given time-step 62 REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: solcp0 !: solid sediment data at initial time 63 REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: diff 61 REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: pwcp, pwcpa 62 REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: solcp, solcpa 63 REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: diff 64 64 65 65 !! * Shared module variables 66 66 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: pwcp_dta !: pore water data at given time-step 67 67 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: rainrm_dta !: rain data at at initial time 68 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: rainrm !: rain data at given time-step69 68 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: rainrg !: rain of each solid component in [g/(cm**2.s)] 70 69 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: fromsed !: 71 70 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: tosed !: 72 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: r loss!:73 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: tokbot71 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: rearatpom !: 72 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: apluss, aminuss !: 74 73 ! 75 74 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: temp !: temperature 76 75 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: salt !: salinity 77 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: press !: pressure78 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: raintg !: total massic flux rained in each cell (sum of sol. comp.)79 76 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: fecratio !: Fe/C ratio in falling particles to the sediments 80 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: dzdep 77 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: dzdep, slatit, slongit !: total thickness of solid material rained [cm] in each cell 81 78 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: zkbot !: total thickness of solid material rained [cm] in each cell 82 79 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: wacc !: total thickness of solid material rained [cm] in each cell … … 87 84 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: volw3d !: ??? 88 85 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: vols3d !: ??? 86 REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: volc 87 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: dens_sol !: Density of each solid fraction 89 88 90 89 … … 104 103 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: ak123ps 105 104 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: aksis 105 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: aknh3 106 106 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: aksps 107 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: akh2s 107 108 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: sieqs 108 109 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: aks3s … … 117 118 INTEGER , PUBLIC, SAVE :: jpoce, indoce !: Ocean points ( number/indices ) 118 119 INTEGER , PUBLIC, DIMENSION(: ), ALLOCATABLE :: iarroce !: Computation of 1D array of sediments points 120 INTEGER , PUBLIC, DIMENSION(:, : ), ALLOCATABLE :: jarr !: Computation of 1D array of sediments points 121 INTEGER , PUBLIC, DIMENSION(: ), ALLOCATABLE :: jsvode, isvode !: Computation of 1D array of sediments points 122 119 123 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: epkbot !: ocean bottom layer thickness 120 124 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: gdepbot !: Depth of the sediment … … 127 131 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: db !: bioturbation ceofficient 128 132 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: irrig !: bioturbation ceofficient 133 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: radssol, rads1sol 134 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: saturco3 129 135 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: rdtsed !: sediment model time-step 130 REAL(wp), PUBLIC, DIMENSION(: ,: ), ALLOCATABLE :: sedligand136 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: rstepros !: Number of iteration of rosenbrock method 131 137 REAL(wp) :: dens !: density of solid material 132 138 !! Inputs / Outputs … … 138 144 CHARACTER( len = 20 ), DIMENSION(jpdia2dsed) :: seddia2d, seddia2u 139 145 ! 140 REAL(wp), PUBLIC, DIMENSION(:,:,:,:), ALLOCATABLE :: trcsedi141 REAL(wp), PUBLIC, DIMENSION(:,:,:,:), ALLOCATABLE :: flxsedi3d142 REAL(wp), PUBLIC, DIMENSION(:,:,: ), ALLOCATABLE :: flxsedi2d143 146 144 147 INTEGER, PUBLIC :: numsed = 27 ! units … … 158 161 & dz(jpksed) , por(jpksed) , por1(jpksed) , & 159 162 & volw(jpksed), vols(jpksed), rdtsed(jpksed) , & 160 & trcsedi (jpi,jpj,jpksed,jptrased) , &161 & flxsedi3d(jpi,jpj,jpksed,jpdia3dsed) , &162 & flxsedi2d(jpi,jpj,jpdia2dsed) , &163 163 & mol_wgt(jpsol), STAT=sed_alloc ) 164 164 -
NEMO/trunk/src/TOP/PISCES/SED/sedadv.F90
r10425 r15450 15 15 16 16 PUBLIC sed_adv 17 PUBLIC sed_adv_alloc18 19 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: dvolsp, dvolsm20 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: c2por, ckpor21 17 22 18 REAL(wp) :: cpor … … 49 45 ! * local variables 50 46 INTEGER :: ji, jk, js 51 INTEGER :: jn , ntimes, nztime, ikwneg47 INTEGER :: jn 52 48 53 REAL(wp), DIMENSION(jpksed,jpsol) :: zsolcpno 54 REAL(wp), DIMENSION(jpksed) :: zfilled, zfull, zfromup, zempty 55 REAL(wp), DIMENSION(jpoce,jpksed) :: zgap, zwb 56 REAL(wp), DIMENSION(jpoce,jpsol) :: zrainrf 57 REAL(wp), DIMENSION(: ), ALLOCATABLE :: zraipush 58 59 REAL(wp) :: zkwnup, zkwnlo, zfrac, zfromce, zrest, sumtot, zsumtot1 49 REAL(wp), DIMENSION(jpoce,jpksed,jpsol+jpads) :: zsolcp 50 REAL(wp) :: solfu, zfilled, zwb, fulsed, uebers, seddef 60 51 61 52 !------------------------------------------------------------------------ 62 63 53 64 54 IF( ln_timing ) CALL timing_start('sed_adv') … … 70 60 WRITE(numsed,*) ' ' 71 61 ENDIF 72 por1clay = denssol * por1(jpksed) * dz(jpksed)73 cpor = por1(jpksed) / por1(2)74 DO jk = 2, jpksed75 c2por(jk) = por1(2) / por1(jk)76 ckpor(jk) = por1(jpksed) / por1(jk)77 ENDDO78 DO jk = jpksedm1, 2, -179 dvolsp(jk) = vols(jk+1) / vols(jk)80 ENDDO81 DO jk = 3, jpksed82 dvolsm(jk) = vols(jk-1) / vols(jk)83 ENDDO84 62 ENDIF 63 64 ! Allocation of the temporary arrays 65 ! ---------------------------------- 66 zsolcp(:,:,:) = 0._wp 67 DO js = 1, jpsol 68 zsolcp(:,:,js) = solcp(:,:,js) 69 END DO 70 DO jk = 2, jpksed 71 zsolcp(:,jk,jpsol+1) = pwcp(:,jk,jwnh4) * adsnh4 72 zsolcp(:,jk,jpsol+2) = pwcp(:,jk,jwfe2) * adsfe2 73 END DO 85 74 86 75 ! Initialization of data for mass balance calculation … … 88 77 fromsed(:,:) = 0. 89 78 tosed (:,:) = 0. 90 rloss (:,:) = 0.91 ikwneg = 192 nztime = jpksed93 79 94 ALLOCATE( zraipush(nztime) ) 95 80 solfu = 0.0 81 DO jk = 2, jpksed 82 solfu = solfu + dz(jk) * por1(jk) 83 END DO 96 84 ! Initiate gap 97 85 !-------------- 98 zgap(:,:) = 0.99 DO js = 1, jpsol100 DO jk = 1, jpksed101 DO ji = 1, jpoce102 zgap(ji,jk) = zgap(ji,jk) + solcp(ji,jk,js)103 END DO104 ENDDO105 ENDDO106 107 zgap(1:jpoce,1:jpksed) = 1. - zgap(1:jpoce,1:jpksed)108 86 109 87 ! Initiate burial rates 110 88 !----------------------- 111 zwb(:,:) = 0.112 DO jk = 2, jpksed113 zfrac = dtsed / ( denssol * por1(jk) )114 DO ji = 1, jpoce115 zwb(ji,jk) = zfrac * raintg(ji)116 END DO117 ENDDO89 DO ji = 1, jpoce 90 DO jk = 2, jpksed-1 91 zfilled = 0._wp 92 DO js = 1, jpsol 93 zfilled = zfilled + zsolcp(ji,jk,js) / dens_sol(js) 94 END DO 95 zwb = MAX(0._wp, (zfilled - 1._wp) / (zfilled + rtrn) ) 118 96 119 97 120 DO ji = 1, jpoce 121 zwb(ji,2) = zwb(ji,2) - zgap(ji,2) * dz(2) 122 ENDDO 98 DO js = 1, jpsol + jpads 99 uebers = zwb * zsolcp(ji,jk,js) 100 zsolcp(ji,jk,js) = zsolcp(ji,jk,js) - uebers 101 zsolcp(ji,jk+1,js) = zsolcp(ji,jk+1,js) + uebers * dz(jk) * por1(jk) / ( dz(jk+1) * por1(jk+1) ) 102 END DO 103 END DO 123 104 124 DO jk = 3, jpksed 125 zfrac = por1(jk-1) / por1(jk) 126 DO ji = 1, jpoce 127 zwb(ji,jk) = zwb(ji,jk-1) * zfrac - zgap(ji,jk) * dz(jk) 105 zfilled = 0._wp 106 DO js = 1, jpsol 107 zfilled = zfilled + zsolcp(ji,jpksed,js) / dens_sol(js) 128 108 END DO 129 ENDDO 130 131 zrainrf(:,:) = 0. 132 DO ji = 1, jpoce 133 IF( raintg(ji) /= 0. ) & 134 & zrainrf(ji,:) = rainrg(ji,:) / raintg(ji) 135 ENDDO 136 137 138 ! Computation of full and empty solid fraction in each layer 139 ! for all 'burial' case 140 !---------------------------------------------------------- 141 109 zwb = MAX(0._wp, (zfilled - 1._wp) / (zfilled + rtrn) ) 110 DO js = 1, jpsol + jpads 111 uebers = zwb * zsolcp(ji,jpksed,js) 112 zsolcp(ji,jpksed,js) = zsolcp(ji,jpksed,js) - uebers 113 tosed(ji,js) = uebers * dz(jpksed) * por1(jpksed) 114 END DO 115 END DO 142 116 143 117 DO ji = 1, jpoce 118 fulsed = 0._wp 119 DO jk = 2, jpksed 120 zfilled = 0._wp 121 DO js = 1, jpsol 122 zfilled = zfilled + zsolcp(ji,jk,js) / dens_sol(js) 123 END DO 124 fulsed = fulsed + zfilled * dz(jk) * por1(jk) 125 END DO 144 126 145 ! computation of total weight fraction in sediment 146 !------------------------------------------------- 147 zfilled(:) = 0. 148 DO js = 1, jpsol 149 DO jk = 2, jpksed 150 zfilled(jk) = zfilled(jk) + solcp(ji,jk,js) 151 ENDDO 152 ENDDO 153 154 DO js = 1, jpsol 155 DO jk = 2, jpksed 156 zsolcpno(jk,js) = solcp(ji,jk,js) / zfilled(jk) 157 ENDDO 158 ENDDO 127 seddef = solfu - fulsed 159 128 160 ! burial 3 cases: 161 ! zwb > 0 ==> rain > total rection loss 162 ! zwb = 0 ==> rain = 0 163 ! zwb < 0 ==> rain > 0 and rain < total reaction loss 164 !---------------------------------------------------------------- 129 zsolcp(ji,jpksed,jsclay) = zsolcp(ji,jpksed,jsclay) + seddef / ( por1(jpksed) * dz(jpksed) ) & 130 & * dens_sol(jsclay) 131 fromsed(ji,jsclay) = seddef * dens_sol(jsclay) 165 132 166 IF( zwb(ji,jpksed) > 0. ) THEN 133 DO jk = jpksed, 3, -1 134 zfilled = 0._wp 135 DO js = 1, jpsol 136 zfilled = zfilled + zsolcp(ji,jk,js) / dens_sol(js) 137 END DO 138 zwb = MAX(0._wp, (zfilled - 1._wp) / (zfilled + rtrn) ) 139 DO js = 1, jpsol + jpads 140 uebers = zwb * zsolcp(ji,jk,js) 141 zsolcp(ji,jk,js) = zsolcp(ji,jk,js) - uebers 142 zsolcp(ji,jk-1,js) = zsolcp(ji,jk-1,js) + uebers * dz(jk) * por1(jk) / ( dz(jk-1) * por1(jk-1) ) 143 END DO 144 END DO 167 145 168 zfull (jpksed) = zfilled(jpksed) 169 zempty(jpksed) = 1. - zfull(jpksed) 170 DO jk = jpksedm1, 2, -1 171 zfull (jk) = zfilled(jk) 172 zfull (jk) = zfull(jk) - zempty(jk+1) * dvolsp(jk) 173 zempty(jk) = 1. - zfull(jk) 174 ENDDO 146 END DO 175 147 176 ! Computation of solid sediment species 177 !-------------------------------------- 178 ! push entire sediment column downward to account rest of rain 179 DO js = 1, jpsol 180 DO jk = jpksed, 3, -1 181 solcp(ji,jk,js) = zfull(jk) * zsolcpno(jk,js) + zempty(jk) * zsolcpno(jk-1,js) 182 ENDDO 148 DO js = 1, jpsol 149 solcp(:,:,js) = zsolcp(:,:,js) 150 END DO 151 DO jk = 2, jpksed 152 pwcp(:,jk,jwnh4) = (pwcp(:,jk,jwnh4) + zsolcp(:,jk,jpsol+1) * por1(jk) / por(jk) ) * radssol(jk,jwnh4) 153 IF (jpoce == 146 .and. slatit(145) > 0.) write(0,*) 'plante advection ',pwcp(145,jk,jwfe2)*1E6,zsolcp(145,jk,jpsol+2)*1E6, & 154 & (pwcp(145,jk,jwfe2) + zsolcp(145,jk,jpsol+2) * por1(jk) / por(jk) ) * radssol(jk,jwfe2) * 1E6 155 pwcp(:,jk,jwfe2) = (pwcp(:,jk,jwfe2) + zsolcp(:,jk,jpsol+2) * por1(jk) / por(jk) ) * radssol(jk,jwfe2) 156 END DO 183 157 184 solcp(ji,2,js) = zfull(2) * zsolcpno(2,js) + zempty(2) * zrainrf(ji,js)185 186 DO jk = 2, jpksed187 zsolcpno(jk,js) = solcp(ji,jk,js)188 END DO189 ENDDO190 191 zrest = zwb(ji,jpksed) * cpor192 ! what is remaining is less than dz(2)193 IF( zrest <= dz(2) ) THEN194 195 zfromup(2) = zrest / dz(2)196 DO jk = 3, jpksed197 zfromup(jk) = zwb(ji,jpksed) * ckpor(jk) / dz(jk)198 ENDDO199 DO js = 1, jpsol200 zfromce = 1. - zfromup(2)201 solcp(ji,2,js) = zfromce * zsolcpno(2,js) + zfromup(2) * zrainrf(ji,js)202 DO jk = 3, jpksed203 zfromce = 1. - zfromup(jk)204 solcp(ji,jk,js) = zfromce * zsolcpno(jk,js) + zfromup(jk) * zsolcpno(jk-1,js)205 ENDDO206 fromsed(ji,js) = 0.207 ! quantities to push in deeper sediment208 tosed (ji,js) = zsolcpno(jpksed,js) &209 & * zwb(ji,jpksed) * denssol * por1(jpksed)210 ENDDO211 212 ELSE ! what is remaining is great than dz(2)213 214 ntimes = INT( zrest / dz(2) ) + 1215 IF( ntimes > nztime ) CALL ctl_stop( 'STOP', 'sed_adv : rest too large ' )216 zraipush(1) = dz(2)217 zrest = zrest - zraipush(1)218 DO jn = 2, ntimes219 IF( zrest >= dz(2) ) THEN220 zraipush(jn) = dz(2)221 zrest = zrest - zraipush(jn)222 ELSE223 zraipush(jn) = zrest224 zrest = 0.225 ENDIF226 ENDDO227 228 DO jn = 1, ntimes229 DO js = 1, jpsol230 DO jk = 2, jpksed231 zsolcpno(jk,js) = solcp(ji,jk,js)232 END DO233 ENDDO234 235 zfromup(2) = zraipush(jn) / dz(2)236 DO jk = 3, jpksed237 zfromup(jk) = ( zraipush(jn) / dz(jk) ) * c2por(jk)238 ENDDO239 240 DO js = 1, jpsol241 zfromce = 1. - zfromup(2)242 solcp(ji,2,js) = zfromce * zsolcpno(2,js) + zfromup(2) * zrainrf(ji,js)243 DO jk = 3, jpksed244 zfromce = 1. - zfromup(jk)245 solcp(ji,jk,js) = zfromce * zsolcpno(jk,js) + zfromup(jk) * zsolcpno(jk-1,js)246 ENDDO247 fromsed(ji,js) = 0.248 tosed (ji,js) = tosed(ji,js) + zsolcpno(jpksed,js) * zraipush(jn) &249 & * denssol * por1(2)250 ENDDO251 ENDDO252 253 ENDIF254 255 ELSE IF( raintg(ji) < eps ) THEN ! rain = 0256 !! Nadia rloss(:,:) = rainrm(:,:) bug ??????257 258 rloss(ji,1:jpsol) = rainrm(ji,1:jpsol)259 260 zfull (2) = zfilled(2)261 zempty(2) = 1. - zfull(2)262 DO jk = 3, jpksed263 zfull (jk) = zfilled(jk)264 zfull (jk) = zfull (jk) - zempty(jk-1) * dvolsm(jk)265 zempty(jk) = 1. - zfull(jk)266 ENDDO267 268 ! fill boxes with weight fraction from underlying box269 DO js = 1, jpsol270 DO jk = 2, jpksedm1271 solcp(ji,jk,js) = zfull(jk) * zsolcpno(jk,js) + zempty(jk) * zsolcpno(jk+1,js)272 END DO273 solcp(ji,jpksed,js) = zsolcpno(jpksed,js) * zfull(jpksed)274 tosed (ji,js) = 0.275 fromsed(ji,js) = 0.276 ENDDO277 ! for the last layer, one make go up clay278 solcp(ji,jpksed,jsclay) = solcp(ji,jpksed,jsclay) + zempty(jpksed) * 1.279 fromsed(ji,jsclay) = zempty(jpksed) * 1. * por1clay280 ELSE ! rain > 0 and rain < total reaction loss281 282 283 DO jk = 2, jpksed284 zfull (jk) = zfilled(jk)285 zempty(jk) = 1. - zfull(jk)286 ENDDO287 288 ! Determination of indice of layer - ikwneg - where advection is reversed289 !------------------------------------------------------------------------290 iflag: DO jk = 2, jpksed291 IF( zwb(ji,jk) < 0. ) THEN292 ikwneg = jk293 EXIT iflag294 ENDIF295 ENDDO iflag296 297 ! computation of zfull and zempty298 ! 3 cases : a/ ikwneg=2, b/ikwneg=3...jpksedm1, c/ikwneg=jpksed299 !-------------------------------------------------------------300 IF( ikwneg == 2 ) THEN ! advection is reversed in the first sediment layer301 302 zkwnup = rdtsed(ikwneg) * raintg(ji) / dz(ikwneg)303 zkwnlo = ABS( zwb(ji,ikwneg) ) / dz(ikwneg)304 zfull (ikwneg+1) = zfilled(ikwneg+1) - zkwnlo * dvolsm(ikwneg+1)305 zempty(ikwneg+1) = 1. - zfull(ikwneg+1)306 DO jk = ikwneg+2, jpksed307 zfull (jk) = zfilled(jk) - zempty(jk-1) * dvolsm(jk)308 zempty(jk) = 1. - zfull(jk)309 ENDDO310 DO js = 1, jpsol311 solcp(ji,2,js) = zfull(2) * zsolcpno(2,js)+ zkwnlo * zsolcpno(3,js) &312 & + zkwnup * zrainrf(ji,js)313 DO jk = 3, jpksedm1314 solcp(ji,jk,js) = zfull(jk) * zsolcpno(jk,js) + zempty(jk) * zsolcpno(jk+1,js)315 ENDDO316 solcp(ji,jpksed,js) = zfull(jpksed) * zsolcpno(jpksed,js)317 tosed(ji,js) = 0.318 fromsed(ji,js) = 0.319 ENDDO320 solcp(ji,jpksed,jsclay) = solcp(ji,jpksed,jsclay) + zempty(jpksed) * 1.321 !! C. Heinze fromsed(ji,jsclay) = zempty(jpksed) * 1. * denssol * por1(jpksed) / mol_wgt(jsclay)322 fromsed(ji,jsclay) = zempty(jpksed) * 1. * por1clay323 324 ELSE IF( ikwneg == jpksed ) THEN325 326 zkwnup = ABS( zwb(ji,ikwneg-1) ) * dvolsm(ikwneg) / dz(ikwneg)327 zkwnlo = ABS( zwb(ji,ikwneg) ) / dz(ikwneg)328 zfull (ikwneg-1) = zfilled(ikwneg-1) - zkwnup * dvolsp(ikwneg-1)329 zempty(ikwneg-1) = 1. - zfull(ikwneg-1)330 DO jk = ikwneg-2, 2, -1331 zfull (jk) = zfilled(jk) - zempty(jk+1) * dvolsp(jk)332 zempty(jk) = 1. - zfull(jk)333 ENDDO334 DO js = 1, jpsol335 solcp(ji,2,js) = zfull(2) * zsolcpno(2,js) + zempty(2) * zrainrf(ji,js)336 ENDDO337 DO js = 1, jpsol338 DO jk = jpksedm1, 3, -1339 solcp(ji,jk,js) = zfull(jk) * zsolcpno(jk,js) + zempty(jk) * zsolcpno(jk-1,js)340 ENDDO341 solcp(ji,jpksed,js) = zfull(jpksed) * zsolcpno(jpksed,js) &342 & + zkwnup * zsolcpno(jpksedm1,js)343 tosed(ji,js) = 0.344 fromsed(ji,js) = 0.345 ENDDO346 solcp(ji,jpksed,jsclay) = solcp(ji,jpksed,jsclay) + zkwnlo * 1.347 ! Heinze fromsed(ji,jsclay) = zkwnlo * 1. * denssol * por1(jpksed) / mol_wgt(jsclay)348 fromsed(ji,jsclay) = zkwnlo * 1.* por1clay349 ELSE ! 2 < ikwneg(ji) <= jpksedm1350 351 zkwnup = ABS( zwb(ji,ikwneg-1) ) * por1(ikwneg-1) / ( dz(ikwneg) * por1(ikwneg) )352 zkwnlo = ABS( zwb(ji,ikwneg) ) / dz(ikwneg)353 354 IF( ikwneg > 3 ) THEN355 356 zfull (ikwneg-1) = zfilled(ikwneg-1) - zkwnup * dvolsp(ikwneg-1)357 zempty(ikwneg-1) = 1. - zfull(ikwneg-1)358 DO jk = ikwneg-2, 2, -1359 zfull (jk) = zfilled(jk) - zempty(jk+1) * dvolsp(jk)360 zempty(jk) = 1. - zfull(jk)361 ENDDO362 DO js = 1, jpsol363 solcp(ji,2,js) = zfull(2) * zsolcpno(2,js) + zempty(2) * zrainrf(ji,js)364 ENDDO365 DO js = 1, jpsol366 DO jk = ikwneg-1, 3, -1367 solcp(ji,jk,js) = zfull(jk) * zsolcpno(jk,js) + zempty(jk) * zsolcpno(jk-1,js)368 ENDDO369 ENDDO370 ELSE ! ikw = 3371 372 373 zfull (2) = zfilled(2) - zkwnup * dvolsm(3)374 zempty(2) = 1. - zfull(2)375 DO js = 1, jpsol376 solcp(ji,2,js) = zfull(2) * zsolcpno(2,js) + zempty(2) * zrainrf(ji,js)377 ENDDO378 ENDIF379 380 IF( ikwneg < jpksedm1) THEN381 382 zfull (ikwneg+1) = zfilled(ikwneg+1) - zkwnlo * dvolsm(ikwneg+1)383 zempty(ikwneg+1) = 1. - zfull(ikwneg+1)384 DO jk = ikwneg+2, jpksed385 zfull (jk) = zfilled(jk) - zempty(jk-1) * dvolsm(jk)386 zempty(jk) = 1. - zfull(jk)387 ENDDO388 DO js = 1, jpsol389 DO jk = ikwneg+1, jpksedm1390 solcp(ji,jk,js) = zfull(jk) * zsolcpno(jk,js) + zempty(jk) * zsolcpno(jk+1,js)391 ENDDO392 solcp(ji,jpksed,js) = zfull(jpksed) * zsolcpno(jpksed,js)393 ENDDO394 solcp(ji,jpksed,jsclay) = solcp(ji,jpksed,jsclay) + zempty(jpksed) * 1.395 ELSE396 397 zfull (jpksed) = zfilled(jpksed) - zkwnlo * dvolsm(jpksed)398 zempty(jpksed) = 1. - zfull(jpksed)399 DO js = 1, jpsol400 solcp(ji,jpksed,js) = zfull(jpksed) * zsolcpno(jpksed,js)401 ENDDO402 solcp(ji,jpksed,jsclay) = solcp(ji,jpksed,jsclay) + zempty(jpksed) * 1.403 ENDIF ! jpksedm1404 405 ! ikwneg = jpksedm1 ; ikwneg+1 = jpksed ; ikwneg-1 = jpksed - 2406 DO js = 1, jpsol407 solcp(ji,ikwneg,js) = zfull(ikwneg) * zsolcpno(ikwneg ,js) &408 & + zkwnup * zsolcpno(ikwneg-1,js) &409 & + zkwnlo * zsolcpno(ikwneg+1,js)410 tosed (ji,js) = 0.411 fromsed(ji,js) = 0.412 ENDDO413 ! Heinze fromsed(ji,jsclay) = zempty * 1. * denssol * por1(jpksed) / mol_wgt(jsclay)414 fromsed(ji,jsclay) = zempty(jpksed) * 1. * por1clay415 416 ENDIF ! ikwneg(ji) = 2417 ENDIF ! zwb > 0418 ENDDO ! ji = 1, jpoce419 420 rainrm(:,:) = 0.421 158 rainrg(:,:) = 0. 422 raintg(:) = 0.423 424 DEALLOCATE( zraipush )425 159 426 160 IF( ln_timing ) CALL timing_stop('sed_adv') … … 429 163 430 164 431 INTEGER FUNCTION sed_adv_alloc()432 !!----------------------------------------------------------------------433 !! *** ROUTINE p4z_prod_alloc ***434 !!----------------------------------------------------------------------435 ALLOCATE( dvolsp(jpksed), dvolsm(jpksed), c2por(jpksed), &436 & ckpor(jpksed) , STAT = sed_adv_alloc )437 !438 IF( sed_adv_alloc /= 0 ) CALL ctl_stop( 'STOP', 'sed_adv_alloc : failed to allocate arrays.' )439 !440 END FUNCTION sed_adv_alloc441 442 165 END MODULE sedadv -
NEMO/trunk/src/TOP/PISCES/SED/sedarr.F90
r14086 r15450 10 10 !!---------------------------------------------------------------------- 11 11 !! * Modules used 12 USE par_oce13 12 USE par_sed 14 USE in_out_manager, ONLY : ln_timing15 USE timing13 USE dom_oce 14 USE sed 16 15 17 16 IMPLICIT NONE … … 52 51 jjd = ( tab_ind(jn) - 1 ) / jpi + 1 53 52 tab1d(jn) = tab2d(jid, jjd) 53 ! IF (mig(jid) == 150 .and. mjg(jjd) == 136) write(0,*) 'plante indices ',jn,ndim1d,slatit(jn),slongit(jn) 54 54 END DO 55 55 -
NEMO/trunk/src/TOP/PISCES/SED/sedbtb.F90
r10222 r15450 6 6 !! * Modules used 7 7 USE sed ! sediment global variable 8 USE sed_oce 8 9 USE sedmat ! linear system of equations 10 USE sedinorg 11 USE sedorg 12 USE sedini 9 13 USE lib_mpp ! distribued memory computing library 10 14 … … 31 35 !! ! 06-04 (C. Ethe) Re-organization 32 36 !!---------------------------------------------------------------------- 33 !!* Arguments 34 INTEGER, INTENT(in) :: kt ! time step 35 37 !! * Arguments 38 INTEGER, INTENT(in) :: kt ! time step 36 39 ! * local variables 37 40 INTEGER :: ji, jk, js 38 REAL(wp), DIMENSION(jpoce,jpksedm1,jpsol) :: zsol ! solution 41 REAL(wp), DIMENSION(jpoce,jpksed,jpsol) :: zrearat ! solution 42 REAL(wp) :: zsolid1, zsolid2, zsolid3, zsumtot, zlimo2 39 43 !------------------------------------------------------------------------ 40 44 … … 42 46 43 47 IF( kt == nitsed000 ) THEN 44 IF (lwp) WRITE(numsed,*) ' sed_btb : Bioturbation'48 IF (lwp) WRITE(numsed,*) ' sed_btb : bioturbation of solid and adsorbed species ' 45 49 IF (lwp) WRITE(numsed,*) ' ' 46 50 ENDIF 47 51 52 48 53 ! Initializations 49 54 !---------------- 50 zsol(:,:,:) = 0. 55 zrearat = 0. 56 57 ! Remineralization rates of the different POC pools 58 zrearat(:,:,jspoc) = -reac_pocl 59 zrearat(:,:,jspos) = -reac_pocs 60 zrearat(:,:,jspor) = -reac_pocr 51 61 52 62 ! right hand side of coefficient matrix 53 63 !-------------------------------------- 54 DO js = 1, jpsol 55 DO jk = 1, jpksedm1 56 DO ji = 1, jpoce 57 zsol(ji,jk,js) = solcp(ji,jk+1,js) 58 ENDDO 59 ENDDO 60 ENDDO 64 CALL sed_mat_btbi( jpksed, jpsol, solcp, zrearat(:,:,:), dtsed ) 61 65 62 CALL sed_mat( jpsol, jpoce, jpksedm1, zsol, dtsed / 2.0 )66 DO ji = 1, jpoce 63 67 68 zsumtot = 0. 69 DO jk = 2, jpksed 70 zsolid1 = volc(ji,jk,jspoc) * solcp(ji,jk,jspoc) 71 zsolid2 = volc(ji,jk,jspos) * solcp(ji,jk,jspos) 72 zsolid3 = volc(ji,jk,jspor) * solcp(ji,jk,jspor) 73 rearatpom(ji,jk) = ( reac_pocl * zsolid1 + reac_pocs * zsolid2 + reac_pocr * zsolid3 ) 74 zsumtot = zsumtot + rearatpom(ji,jk) * volw3d(ji,jk) * 1.e-3 * 86400. * 365. * 1E3 75 END DO 64 76 65 ! store solution of the tridiagonal system 66 !------------------------ 67 DO js = 1, jpsol 68 DO jk = 1, jpksedm1 69 DO ji = 1, jpoce 70 solcp(ji,jk+1,js) = zsol(ji,jk,js) 71 ENDDO 72 ENDDO 73 ENDDO 77 ! 4/ Computation of the bioturbation coefficient 78 ! This parameterization is taken from Archer et al. (2002) 79 ! -------------------------------------------------------------- 80 zlimo2 = max(0.01, pwcp(ji,1,jwoxy) / (pwcp(ji,1,jwoxy) + 20.E-6) ) 81 db(ji,:) = dbiot * zsumtot**0.85 * zlimo2 / (365.0 * 86400.0) 82 83 ! ------------------------------------------------------ 84 ! Vertical variations of the bioturbation coefficient 85 ! ------------------------------------------------------ 86 IF (ln_btbz) THEN 87 DO jk = 1, jpksed 88 db(ji,jk) = db(ji,jk) * exp( -(profsedw(jk) / dbtbzsc)**2 ) 89 END DO 90 ELSE 91 DO jk = 1, jpksed 92 IF (profsedw(jk) > dbtbzsc) THEN 93 db(ji,jk) = 0.0 94 ENDIF 95 END DO 96 ENDIF 97 98 ! Computation of the bioirrigation factor (from Archer, MUDS model) 99 irrig(ji,:) = 0.0 100 IF (ln_irrig) THEN 101 DO jk = 1, jpksed 102 irrig(ji,jk) = ( 7.63752 - 7.4465 * exp( -0.89603 * zsumtot ) ) * zlimo2 103 irrig(ji,jk) = irrig(ji,jk) * exp( -(profsedw(jk) / xirrzsc) ) 104 END DO 105 ENDIF 106 END DO 107 108 ! CALL inorganic and organic slow redow/chemistry processes 109 ! --------------------------------------------------------- 110 CALL sed_inorg( kt ) 74 111 75 112 IF( ln_timing ) CALL timing_stop('sed_btb') -
NEMO/trunk/src/TOP/PISCES/SED/sedchem.F90
r15237 r15450 139 139 CALL sed_chem_cst 140 140 ELSE 141 DO_2D( 1, 1, 1, 1)141 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 142 142 ikt = mbkt(ji,jj) 143 143 IF ( tmask(ji,jj,ikt) == 1 ) THEN -
NEMO/trunk/src/TOP/PISCES/SED/seddiff.F90
r10225 r15450 7 7 !! * Modules used 8 8 USE sed ! sediment global variable 9 USE sed_oce10 9 USE sedmat ! linear system of equations 11 10 USE sedini … … 42 41 !! Arguments 43 42 INTEGER, INTENT(in) :: kt, knt ! number of iteration 44 ! --- local variables 45 INTEGER :: ji, jk, js ! dummy looop indices 43 ! 46 44 47 REAL(wp), DIMENSION(jpoce,jpksed) :: zrearat1, zrearat2 ! reaction rate in pore water48 !!49 !!----------------------------------------------------------------------50 51 IF( ln_timing ) CALL timing_start('sed_diff')52 !53 IF( kt == nitsed000 .AND. knt == 1 ) THEN54 IF (lwp) THEN55 WRITE(numsed,*) ' sed_diff : pore-water diffusion '56 WRITE(numsed,*) ' '57 ENDIF58 ENDIF59 60 ! Initializations61 !----------------------62 zrearat1(:,:) = 0.63 zrearat2(:,:) = 0.64 65 !---------------------------66 ! Solves PO4 diffusion67 !----------------------------68 69 ! solves tridiagonal system70 CALL sed_mat( jwpo4, jpoce, jpksed, zrearat1, zrearat2, pwcp(:,:,jwpo4), dtsed2 / 2.0 )71 72 !---------------------------73 ! Solves NH4 diffusion74 !----------------------------75 76 ! solves tridiagonal system77 CALL sed_mat( jwnh4, jpoce, jpksed, zrearat1, zrearat2, pwcp(:,:,jwnh4), dtsed2 / 2.0 )78 79 !---------------------------80 ! Solves Fe2+ diffusion81 !----------------------------82 83 ! solves tridiagonal system84 CALL sed_mat( jwfe2, jpoce, jpksed, zrearat1, zrearat2, pwcp(:,:,jwfe2), dtsed2 / 2.0 )85 86 !---------------------------87 ! Solves H2S diffusion88 !----------------------------89 90 ! solves tridiagonal system91 CALL sed_mat( jwh2s, jpoce, jpksed, zrearat1, zrearat2, pwcp(:,:,jwh2s), dtsed2 / 2.0 )92 93 !---------------------------94 ! Solves SO4 diffusion95 !----------------------------96 97 ! solves tridiagonal system98 CALL sed_mat( jwso4, jpoce, jpksed, zrearat1, zrearat2, pwcp(:,:,jwso4), dtsed2 / 2.0 )99 100 !---------------------------101 ! Solves O2 diffusion102 !----------------------------103 104 ! solves tridiagonal system105 CALL sed_mat( jwoxy, jpoce, jpksed, zrearat1, zrearat2, pwcp(:,:,jwoxy), dtsed2 / 2.0 )106 107 !---------------------------108 ! Solves NO3 diffusion109 !----------------------------110 111 ! solves tridiagonal system112 CALL sed_mat( jwno3, jpoce, jpksed, zrearat1, zrearat2, pwcp(:,:,jwno3), dtsed2 / 2.0 )113 114 CALL sed_mat( jwdic, jpoce, jpksed, zrearat1, zrearat2, sedligand(:,:), dtsed2 / 2.0 )115 116 IF( ln_timing ) CALL timing_stop('sed_diff')117 !118 45 END SUBROUTINE sed_diff 119 46 -
NEMO/trunk/src/TOP/PISCES/SED/seddsr.F90
r10362 r15450 19 19 !! * Module variables 20 20 21 REAL(wp) :: zadsnh422 21 REAL(wp), DIMENSION(jpsol), PUBLIC :: dens_mol_wgt ! molecular density 23 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zvolc ! temp. variables24 25 22 26 23 !! $Id$ 27 24 CONTAINS 28 25 29 SUBROUTINE sed_dsr( kt, knt )26 SUBROUTINE sed_dsr( accmask ) 30 27 !!---------------------------------------------------------------------- 31 28 !! *** ROUTINE sed_dsr *** … … 49 46 !!---------------------------------------------------------------------- 50 47 !! Arguments 51 INTEGER, INTENT(in) :: kt, knt ! number of iteration52 48 ! --- local variables 49 INTEGER, DIMENSION(jpoce), INTENT(in) :: accmask 53 50 INTEGER :: ji, jk, js, jw, jn ! dummy looop indices 54 51 55 REAL(wp), DIMENSION(jpoce,jpksed) :: zrearat1, zrearat2, zrearat3 ! reaction rate in pore water 56 REAL(wp), DIMENSION(jpoce,jpksed) :: zundsat ! undersaturation ; indice jpwatp1 is for calcite 57 REAL(wp), DIMENSION(jpoce,jpksed) :: zkpoc, zkpos, zkpor, zlimo2, zlimno3, zlimso4, zlimfeo ! undersaturation ; indice jpwatp1 is for calcite 58 REAL(wp), DIMENSION(jpoce) :: zsumtot 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 52 REAL(wp), DIMENSION(jpoce,jpksed) ::zlimo2, zlimno3, zlimso4, zlimfeo ! undersaturation ; indice jpwatp1 is for calcite 53 REAL(wp) :: zreasat 62 54 !! 63 55 !!---------------------------------------------------------------------- … … 65 57 IF( ln_timing ) CALL timing_start('sed_dsr') 66 58 ! 67 IF( kt == nitsed000 .AND. knt == 1 ) THEN68 IF (lwp) THEN69 WRITE(numsed,*) ' sed_dsr : Dissolution reaction '70 WRITE(numsed,*) ' '71 ENDIF72 ENDIF73 74 59 ! Initializations 75 60 !---------------------- 76 61 77 zrearat1(:,:) = 0. ; zundsat(:,:) = 0. ; zkpoc(:,:) = 0. 78 zlimo2 (:,:) = 0. ; zlimno3(:,:) = 0. ; zrearat2(:,:) = 0. 79 zlimso4(:,:) = 0. ; zkpor(:,:) = 0. ; zrearat3(:,:) = 0. 80 zkpos (:,:) = 0. 81 zsumtot(:) = rtrn 62 zlimo2 (:,:) = 0. ; zlimno3(:,:) = 0. 63 zlimso4(:,:) = 0. 82 64 83 ALLOCATE( zvolc(jpoce, jpksed, jpsol) )84 zvolc(:,:,:) = 0.85 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 97 ! Conversion of volume units98 !----------------------------99 DO js = 1, jpsol100 DO jk = 1, jpksed101 DO ji = 1, jpoce102 zvolc(ji,jk,js) = ( vols3d(ji,jk) * dens_mol_wgt(js) ) / &103 & ( volw3d(ji,jk) * 1.e-3 )104 ENDDO105 ENDDO106 ENDDO107 108 65 !---------------------------------------------------------- 109 66 ! 5. Beginning of solid reaction 110 67 !--------------------------------------------------------- 111 68 112 ! Definition of reaction rates [rearat]=sans dim 113 ! For jk=1 no reaction (pure water without solid) for each solid compo 114 zrearat1(:,:) = 0. 115 zrearat2(:,:) = 0. 116 zrearat3(:,:) = 0. 117 118 zundsat(:,:) = pwcp(:,:,jwoxy) 119 120 DO jk = 2, jpksed 121 DO ji = 1, jpoce 122 zlimo2(ji,jk) = 1.0 / ( zundsat(ji,jk) + xksedo2 ) 123 zsolid1 = zvolc(ji,jk,jspoc) * solcp(ji,jk,jspoc) 124 zsolid2 = zvolc(ji,jk,jspos) * solcp(ji,jk,jspos) 125 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 186 ! For DIC 187 pwcp(ji,jk,jwdic) = pwcp(ji,jk,jwdic) + zreasat 188 zsumtot(ji) = zsumtot(ji) + zreasat / dtsed2 * volw3d(ji,jk) * 1.e-3 * 86400. * 365. * 1E3 189 ! For Phosphate (in mol/l) 190 pwcp(ji,jk,jwpo4) = pwcp(ji,jk,jwpo4) + zreasat * spo4r 191 ! For iron (in mol/l) 192 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 ) 195 ! Ammonium 196 pwcp(ji,jk,jwnh4) = pwcp(ji,jk,jwnh4) + zreasat * srno3 * zadsnh4 197 ! Ligands 198 sedligand(ji,jk) = sedligand(ji,jk) + ratligc * zreasat - reac_ligc * sedligand(ji,jk) 199 ENDDO 69 ! Computes the SMS of the species which do not affect the remin 70 ! processes (Alkalinity, PO4, NH4, and the release of Fe from 71 ! biogenic Fe 72 ! In the following, all loops start from jpk = 2 as reactions 73 ! only occur in the sediment 74 ! -------------------------------------------------------------- 75 DO jk = 2, jpksed 76 DO ji = 1, jpoce 77 IF (accmask(ji) == 0 ) THEN 78 zreasat = rearatpom(ji,jk) 79 ! For alkalinity 80 pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) + zreasat * ( srno3 - 2.* spo4r ) 81 ! For Phosphate (in mol/l) 82 pwcpa(ji,jk,jwpo4) = pwcpa(ji,jk,jwpo4) + zreasat * spo4r 83 ! For Ammonium 84 pwcpa(ji,jk,jwnh4) = pwcpa(ji,jk,jwnh4) + zreasat * srno3 * radssol(jk,jwnh4) 85 ! For iron (in mol/l) 86 pwcpa(ji,jk,jwfe2) = pwcpa(ji,jk,jwfe2) + fecratio(ji) * zreasat * radssol(jk,jwfe2) 87 ENDIF 88 END DO 89 ENDDO 90 91 ! Computes SMS of oxygen 92 DO jk = 2, jpksed 93 DO ji = 1, jpoce 94 IF (accmask(ji) == 0 ) THEN 95 zlimo2(ji,jk) = pwcp(ji,jk,jwoxy) / ( pwcp(ji,jk,jwoxy) + xksedo2 ) 96 zlimo2(ji,jk) = MAX( 0., zlimo2(ji,jk) ) 97 zreasat = rearatpom(ji,jk) * zlimo2(ji,jk) 98 ! Acid Silicic 99 pwcpa(ji,jk,jwoxy) = pwcpa(ji,jk,jwoxy) - so2ut * zreasat 100 ENDIF 101 END DO 102 ENDDO 103 104 !-------------------------------------------------------------------- 105 ! Denitrification 106 !-------------------------------------------------------------------- 107 DO jk = 2, jpksed 108 DO ji = 1, jpoce 109 IF (accmask(ji) == 0 ) THEN 110 zlimno3(ji,jk) = ( 1.0 - zlimo2(ji,jk) ) * pwcp(ji,jk,jwno3) / ( pwcp(ji,jk,jwno3) + xksedno3 ) 111 zlimno3(ji,jk) = MAX(0., zlimno3(ji,jk) ) 112 zreasat = rearatpom(ji,jk) * zlimno3(ji,jk) * srDnit 113 ! For nitrates 114 pwcpa(ji,jk,jwno3) = pwcpa(ji,jk,jwno3) - zreasat 115 ! For alkalinity 116 pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) + zreasat 117 ENDIF 118 END DO 119 ENDDO 120 121 122 !-------------------------------------------------------------------- 123 ! Begining POC iron reduction 124 !-------------------------------------------------------------------- 125 126 DO jk = 2, jpksed 127 DO ji = 1, jpoce 128 IF (accmask(ji) == 0 ) THEN 129 zlimfeo(ji,jk) = ( 1.0 - zlimno3(ji,jk) - zlimo2(ji,jk) ) * solcp(ji,jk,jsfeo) / ( solcp(ji,jk,jsfeo) + xksedfeo ) 130 zlimfeo(ji,jk) = MAX(0., zlimfeo(ji,jk) ) 131 zreasat = 4.0 * rearatpom(ji,jk) * zlimfeo(ji,jk) 132 ! For FEOH 133 solcpa(ji,jk,jsfeo) = solcpa(ji,jk,jsfeo) - zreasat / volc(ji,jk,jsfeo) 134 ! For PO4 135 pwcpa(ji,jk,jwpo4) = pwcpa(ji,jk,jwpo4) + zreasat * redfep 136 ! For alkalinity 137 pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) + 2.0 * zreasat 138 ! Iron 139 pwcpa(ji,jk,jwfe2) = pwcpa(ji,jk,jwfe2) + zreasat * radssol(jk,jwfe2) 140 ENDIF 141 END DO 200 142 ENDDO 201 143 202 144 !-------------------------------------------------------------------- 203 145 ! Begining POC denitrification and NO3- diffusion 204 ! (indice n°5 for couple POC/NO3- ie solcp(:,:,jspoc)/pwcp(:,:,jwno3)) 205 !-------------------------------------------------------------------- 206 207 zrearat1(:,:) = 0. 208 zrearat2(:,:) = 0. 209 zrearat3(:,:) = 0. 210 211 zundsat(:,:) = pwcp(:,:,jwno3) 212 213 DO jk = 2, jpksed 214 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 277 ! 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 289 ! 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 293 ENDDO 294 ENDDO 295 296 !-------------------------------------------------------------------- 297 ! Begining POC iron reduction 298 ! (indice n�5 for couple POFe(OH)3 ie solcp(:,:,jspoc)/pwcp(:,:,jsfeo)) 299 !-------------------------------------------------------------------- 300 301 zrearat1(:,:) = 0. 302 zrearat2(:,:) = 0. 303 zrearat3(:,:) = 0. 304 305 zundsat(:,:) = solcp(:,:,jsfeo) 306 307 DO jk = 2, jpksed 308 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) 375 ! 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 380 ! 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 386 ! For alkalinity 387 pwcp(ji,jk,jwalk) = pwcp(ji,jk,jwalk) + zreasat * ( srno3 * zadsnh4 - 2.* spo4r ) + 8.0 * zreasat 388 ! Ammonium 389 pwcp(ji,jk,jwnh4) = pwcp(ji,jk,jwnh4) + zreasat * srno3 * zadsnh4 390 pwcp(ji,jk,jwfe2) = pwcp(ji,jk,jwfe2) + zreasat * 4.0 391 ENDDO 392 ENDDO 393 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) 404 405 DO jk = 2, jpksed 406 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 474 ! 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 487 ! 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 491 ENDDO 492 ENDDO 493 494 ! Oxydation of the reduced products. Here only ammonium and ODU is accounted for 495 ! There are two options here: A simple time splitting scheme and a modified 496 ! Patankar scheme 497 ! ------------------------------------------------------------------------------ 498 499 call sed_dsr_redoxb 500 501 ! -------------------------------------------------------------- 502 ! 4/ Computation of the bioturbation coefficient 503 ! This parameterization is taken from Archer et al. (2002) 504 ! -------------------------------------------------------------- 505 506 DO ji = 1, jpoce 507 db(ji,:) = dbiot * zsumtot(ji) * pwcp(ji,1,jwoxy) / (pwcp(ji,1,jwoxy) + 20.E-6) 508 END DO 509 510 ! ------------------------------------------------------ 511 ! Vertical variations of the bioturbation coefficient 512 ! ------------------------------------------------------ 513 IF (ln_btbz) THEN 514 DO ji = 1, jpoce 515 db(ji,:) = db(ji,:) * exp( -(profsedw(:) / dbtbzsc)**2 ) / (365.0 * 86400.0) 516 END DO 517 ELSE 518 DO jk = 1, jpksed 519 IF (profsedw(jk) > dbtbzsc) THEN 520 db(:,jk) = 0.0 521 ENDIF 522 END DO 523 ENDIF 524 525 IF (ln_irrig) THEN 526 DO jk = 1, jpksed 527 DO ji = 1, jpoce 528 irrig(ji,jk) = ( 7.63752 - 7.4465 * exp( -0.89603 * zsumtot(ji) ) ) * pwcp(ji,1,jwoxy) & 529 & / (pwcp(ji,1,jwoxy) + 20.E-6) 530 irrig(ji,jk) = irrig(ji,jk) * exp( -(profsedw(jk) / xirrzsc) ) 531 END DO 532 END DO 533 ELSE 534 irrig(:,:) = 0.0 535 ENDIF 536 537 DEALLOCATE( zvolc ) 146 !-------------------------------------------------------------------- 147 148 DO jk = 2, jpksed 149 DO ji = 1, jpoce 150 IF (accmask(ji) == 0 ) THEN 151 zlimso4(ji,jk) = ( 1.0 - zlimno3(ji,jk) - zlimo2(ji,jk) - zlimfeo(ji,jk) ) * pwcp(ji,jk,jwso4) / ( pwcp(ji,jk,jwso4) + xksedso4 ) 152 zlimso4(ji,jk) = MAX(0., zlimso4(ji,jk) ) 153 zreasat = 0.5 * rearatpom(ji,jk) * zlimso4(ji,jk) 154 ! For sulfur 155 pwcpa(ji,jk,jwso4) = pwcpa(ji,jk,jwso4) - zreasat 156 pwcpa(ji,jk,jwh2s) = pwcpa(ji,jk,jwh2s) + zreasat 157 ! For alkalinity 158 pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) + 2.0 * zreasat 159 ENDIF 160 END DO 161 ENDDO 162 163 ! Secondary redox reactions 164 ! ------------------------- 165 166 call sed_dsr_redoxb( accmask ) 538 167 539 168 IF( ln_timing ) CALL timing_stop('sed_dsr') … … 541 170 END SUBROUTINE sed_dsr 542 171 543 SUBROUTINE sed_dsr_redoxb 172 SUBROUTINE sed_dsr_redoxb( accmask ) 544 173 !!---------------------------------------------------------------------- 545 174 !! *** ROUTINE sed_dsr_redox *** … … 547 176 !! ** Purpose : computes secondary redox reactions 548 177 !! 549 !! ** Methode : It uses Strand splitter algorithm proposed by550 !! Nguyen et al. (2009) and modified by Wang et al. (2018)551 !! Basically, each equation is solved analytically when552 !! feasible, otherwise numerically at t+1/2. Then553 !! the last equation is solved at t+1. The other equations554 !! are then solved at t+1 starting in the reverse order.555 !! Ideally, it's better to start from the fastest reaction556 !! to the slowest and then reverse the order to finish up557 !! with the fastest one. But random order works well also.558 !! The scheme is second order, positive and mass559 !! conserving. It works well for stiff systems.560 !!561 178 !! History : 562 179 !! ! 18-08 (O. Aumont) Original code … … 564 181 !! Arguments 565 182 ! --- local variables 566 INTEGER :: ji, jk, jn ! dummy looop indices567 568 REAL, DIMENSION(6) :: zsedtrn, zsedtra 569 REAL(wp) :: zalpha, z beta, zgamma, zdelta, zepsi, zsedfer183 INTEGER, DIMENSION(jpoce), INTENT(in) :: accmask 184 INTEGER :: ji, jni, jnj, jk ! dummy looop indices 185 186 REAL(wp) :: zalpha, zexcess, zh2seq, zsedfer, zreasat 570 187 !! 571 188 !!---------------------------------------------------------------------- … … 573 190 IF( ln_timing ) CALL timing_start('sed_dsr_redoxb') 574 191 575 DO ji = 1, jpoce 576 DO jk = 2, jpksed 577 zsedtrn(1) = pwcp(ji,jk,jwoxy) 578 zsedtrn(2) = MAX(0., pwcp(ji,jk,jwh2s) ) 579 zsedtrn(3) = pwcp(ji,jk,jwnh4) 580 zsedtrn(4) = MAX(0., pwcp(ji,jk,jwfe2) - sedligand(ji,jk) ) 581 zsedfer = MIN(0., pwcp(ji,jk,jwfe2) - sedligand(ji,jk) ) 582 zsedtrn(5) = solcp(ji,jk,jsfeo) * zvolc(ji,jk,jsfeo) 583 zsedtrn(6) = solcp(ji,jk,jsfes) * zvolc(ji,jk,jsfes) 584 zsedtra(:) = zsedtrn(:) 585 586 ! First pass of the scheme. At the end, it is 1st order 587 ! ----------------------------------------------------- 588 ! Fe + O2 589 zalpha = zsedtra(1) - 0.25 * zsedtra(4) 590 zbeta = zsedtra(4) + zsedtra(5) 591 zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4) 592 zdelta = pwcp(ji,jk,jwpo4) - redfep * zsedtra(4) 593 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 ) & 597 & + zalpha * exp( reac_fe2 * zalpha * dtsed2 / 2. ) ) 598 ENDIF 599 zsedtra(1) = zalpha + 0.25 * zsedtra(4) 600 zsedtra(5) = zbeta - zsedtra(4) 601 pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(4) 602 pwcp(ji,jk,jwpo4) = zdelta + redfep * zsedtra(4) 603 ! H2S + O2 604 zalpha = zsedtra(1) - 2.0 * zsedtra(2) 605 zbeta = pwcp(ji,jk,jwso4) + zsedtra(2) 606 zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(2) 607 IF ( zalpha == 0. ) THEN 608 zsedtra(2) = zsedtra(2) / ( 1.0 + zsedtra(2) * reac_h2s * dtsed2 / 2.0 ) 609 ELSE 610 zsedtra(2) = ( zsedtra(2) * zalpha ) / ( 2.0 * zsedtra(2) * ( exp( reac_h2s * zalpha * dtsed2 / 2. ) - 1.0 ) & 611 & + zalpha * exp( reac_h2s * zalpha * dtsed2 / 2. ) ) 612 ENDIF 613 zsedtra(1) = zalpha + 2.0 * zsedtra(2) 614 pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(2) 615 pwcp(ji,jk,jwso4) = zbeta - zsedtra(2) 616 ! 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) 627 ! FeS - O2 628 zalpha = zsedtra(1) - 2.0 * zsedtra(6) 629 zbeta = zsedtra(4) + zsedtra(6) 630 zgamma = pwcp(ji,jk,jwso4) + zsedtra(6) 631 IF ( zalpha == 0. ) THEN 632 zsedtra(6) = zsedtra(6) / ( 1.0 + zsedtra(6) * reac_feso * dtsed2 / 2.0 ) 633 ELSE 634 zsedtra(6) = ( zsedtra(6) * zalpha ) / ( 2.0 * zsedtra(6) * ( exp( reac_feso * zalpha * dtsed2 / 2. ) - 1.0 ) & 635 & + zalpha * exp( reac_feso * zalpha * dtsed2 /2. ) ) 636 ENDIF 637 zsedtra(1) = zalpha + 2.0 * zsedtra(6) 638 zsedtra(4) = zbeta - zsedtra(6) 639 pwcp(ji,jk,jwso4) = zgamma - zsedtra(6) 640 ! ! Fe - H2S 641 zalpha = zsedtra(2) - zsedtra(4) 642 zbeta = zsedtra(4) + zsedtra(6) 643 zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4) 644 IF ( zalpha == 0. ) THEN 645 zsedtra(4) = zsedtra(4) / ( 1.0 + zsedtra(4) * reac_fes * dtsed2 / 2.0 ) 646 ELSE 647 zsedtra(4) = ( zsedtra(4) * zalpha ) / ( zsedtra(4) * ( exp( reac_fes * zalpha * dtsed2 / 2. ) - 1.0 ) & 648 & + zalpha * exp( reac_fes * zalpha * dtsed2 /2. ) ) 649 ENDIF 650 zsedtra(2) = zalpha + zsedtra(4) 651 zsedtra(6) = zbeta - zsedtra(4) 652 pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(4) 653 ! FEOH + H2S 654 zalpha = zsedtra(5) - 2.0 * zsedtra(2) 655 zbeta = zsedtra(5) + zsedtra(4) 656 zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4) 657 zdelta = pwcp(ji,jk,jwso4) + zsedtra(2) 658 zepsi = pwcp(ji,jk,jwpo4) + redfep * zsedtra(5) 659 IF ( zalpha == 0. ) THEN 660 zsedtra(2) = zsedtra(2) / ( 1.0 + zsedtra(2) * reac_feh2s * dtsed2 ) 661 ELSE 662 zsedtra(2) = ( zsedtra(2) * zalpha ) / ( 2.0 * zsedtra(2) * ( exp( reac_feh2s * zalpha * dtsed2 ) - 1.0 ) & 663 & + zalpha * exp( reac_feh2s * zalpha * dtsed2 ) ) 664 ENDIF 665 zsedtra(5) = zalpha + 2.0 * zsedtra(2) 666 zsedtra(4) = zbeta - zsedtra(5) 667 pwcp(ji,jk,jwso4) = zdelta - zsedtra(2) 668 pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(4) 669 pwcp(ji,jk,jwpo4) = zepsi - redfep * zsedtra(5) 670 ! Fe - H2S 671 zalpha = zsedtra(2) - zsedtra(4) 672 zbeta = zsedtra(4) + zsedtra(6) 673 zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4) 674 IF ( zalpha == 0. ) THEN 675 zsedtra(4) = zsedtra(4) / ( 1.0 + zsedtra(4) * reac_fes * dtsed2 / 2.0 ) 676 ELSE 677 zsedtra(4) = ( zsedtra(4) * zalpha ) / ( zsedtra(4) * ( exp( reac_fes * zalpha * dtsed2 / 2. ) - 1.0 ) & 678 & + zalpha * exp( reac_fes * zalpha * dtsed2 /2. ) ) 679 ENDIF 680 zsedtra(2) = zalpha + zsedtra(4) 681 zsedtra(6) = zbeta - zsedtra(4) 682 pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(4) 683 ! FeS - O2 684 zalpha = zsedtra(1) - 2.0 * zsedtra(6) 685 zbeta = zsedtra(4) + zsedtra(6) 686 zgamma = pwcp(ji,jk,jwso4) + zsedtra(6) 687 IF (zalpha == 0.) THEN 688 zsedtra(6) = zsedtra(6) / ( 1.0 + zsedtra(6) * reac_feso * dtsed2 / 2. ) 689 ELSE 690 zsedtra(6) = ( zsedtra(6) * zalpha ) / ( 2.0 * zsedtra(6) * ( exp( reac_feso * zalpha * dtsed2 / 2. ) - 1.0 ) & 691 & + zalpha * exp( reac_feso * zalpha * dtsed2 /2. ) ) 692 ENDIF 693 zsedtra(1) = zalpha + 2.0 * zsedtra(6) 694 zsedtra(4) = zbeta - zsedtra(6) 695 pwcp(ji,jk,jwso4) = zgamma - zsedtra(6) 696 ! NH4 + O2 697 zalpha = zsedtra(1) - 2.0 * zsedtra(3) 698 zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(3) 699 IF (zalpha == 0.) THEN 700 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) 707 ! H2S + O2 708 zalpha = zsedtra(1) - 2.0 * zsedtra(2) 709 zbeta = pwcp(ji,jk,jwso4) + zsedtra(2) 710 zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(2) 711 IF ( zalpha == 0. ) THEN 712 zsedtra(2) = zsedtra(2) / ( 1.0 + zsedtra(2) * reac_h2s * dtsed2 / 2.0 ) 713 ELSE 714 zsedtra(2) = ( zsedtra(2) * zalpha ) / ( 2.0 * zsedtra(2) * ( exp( reac_h2s * zalpha * dtsed2 / 2. ) - 1.0 ) & 715 & + zalpha * exp( reac_h2s * zalpha * dtsed2 / 2. ) ) 716 ENDIF 717 zsedtra(1) = zalpha + 2.0 * zsedtra(2) 718 pwcp(ji,jk,jwso4) = zbeta - zsedtra(2) 719 pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(2) 720 ! Fe + O2 721 zalpha = zsedtra(1) - 0.25 * zsedtra(4) 722 zbeta = zsedtra(4) + zsedtra(5) 723 zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4) 724 zdelta = pwcp(ji,jk,jwpo4) - redfep * zsedtra(4) 725 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 ) & 729 & + zalpha * exp( reac_fe2 * zalpha * dtsed2 / 2. ) ) 730 ENDIF 731 zsedtra(1) = zalpha + 0.25 * zsedtra(4) 732 zsedtra(5) = zbeta - zsedtra(4) 733 pwcp(ji,jk,jwpo4) = zdelta + redfep * zsedtra(4) 734 pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(4) 735 pwcp(ji,jk,jwoxy) = zsedtra(1) 736 pwcp(ji,jk,jwh2s) = zsedtra(2) 737 pwcp(ji,jk,jwnh4) = zsedtra(3) 738 pwcp(ji,jk,jwfe2) = zsedtra(4) + sedligand(ji,jk) + zsedfer 739 pwcp(ji,jk,jwno3) = pwcp(ji,jk,jwno3) + ( zsedtrn(3) - pwcp(ji,jk,jwnh4) ) 740 solcp(ji,jk,jsfeo) = zsedtra(5) / zvolc(ji,jk,jsfeo) 741 solcp(ji,jk,jsfes) = zsedtra(6) / zvolc(ji,jk,jsfes) 742 END DO 743 END DO 192 DO jk = 2, jpksed 193 DO ji = 1, jpoce 194 IF (accmask(ji) == 0 ) THEN 195 zalpha = ( pwcp(ji,jk,jwfe2) - pwcp(ji,jk,jwlgw) ) * 1E9 196 zsedfer = ( zalpha + SQRT( zalpha**2 + 1.E-10 ) ) /2.0 * 1E-9 197 198 ! First pass of the scheme. At the end, it is 1st order 199 ! ----------------------------------------------------- 200 ! Fe (both adsorbed and solute) + O2 201 zreasat = reac_fe2 * zsedfer / radssol(jk,jwfe2) * pwcp(ji,jk,jwoxy) 202 pwcpa(ji,jk,jwfe2) = pwcpa(ji,jk,jwfe2) - zreasat * radssol(jk,jwfe2) 203 pwcpa(ji,jk,jwoxy) = pwcpa(ji,jk,jwoxy) - 0.25 * zreasat 204 pwcpa(ji,jk,jwpo4) = pwcpa(ji,jk,jwpo4) - redfep * zreasat 205 pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) - 2.0 * zreasat 206 solcpa(ji,jk,jsfeo) = solcpa(ji,jk,jsfeo) + zreasat / volc(ji,jk,jsfeo) 207 208 ! H2S + O2 209 zreasat = reac_h2s * pwcp(ji,jk,jwoxy) * pwcp(ji,jk,jwh2s) 210 pwcpa(ji,jk,jwh2s) = pwcpa(ji,jk,jwh2s) - zreasat 211 pwcpa(ji,jk,jwoxy) = pwcpa(ji,jk,jwoxy) - 2.0 * zreasat 212 pwcpa(ji,jk,jwso4) = pwcpa(ji,jk,jwso4) + zreasat 213 pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) - 2.0 * zreasat 214 215 ! NH4 + O2 216 zreasat = reac_nh4 * pwcp(ji,jk,jwnh4) / radssol(jk,jwnh4) * pwcp(ji,jk,jwoxy) 217 pwcpa(ji,jk,jwnh4) = pwcpa(ji,jk,jwnh4) - zreasat * radssol(jk,jwnh4) 218 pwcpa(ji,jk,jwoxy) = pwcpa(ji,jk,jwoxy) - 2.0 * zreasat 219 pwcpa(ji,jk,jwno3) = pwcpa(ji,jk,jwno3) + zreasat 220 pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) - 2.0 * zreasat 221 222 ! FeS - O2 223 zreasat = reac_feso * solcp(ji,jk,jsfes) * volc(ji,jk,jsfes) * pwcp(ji,jk,jwoxy) 224 solcpa(ji,jk,jsfes) = solcpa(ji,jk,jsfes) - zreasat / volc(ji,jk,jsfes) 225 pwcpa(ji,jk,jwoxy) = pwcpa(ji,jk,jwoxy) - 2.0 * zreasat 226 pwcpa(ji,jk,jwfe2) = pwcpa(ji,jk,jwfe2) + zreasat * radssol(jk,jwfe2) 227 pwcpa(ji,jk,jwso4) = pwcpa(ji,jk,jwso4) + zreasat 228 229 ! FEOH + H2S 230 zreasat = reac_feh2s * solcp(ji,jk,jsfeo) * volc(ji,jk,jsfeo) * pwcp(ji,jk,jwh2s) 231 solcpa(ji,jk,jsfeo) = solcpa(ji,jk,jsfeo) - 8.0 * zreasat / volc(ji,jk,jsfeo) 232 pwcpa(ji,jk,jwh2s) = pwcpa(ji,jk,jwh2s) - zreasat 233 pwcpa(ji,jk,jwfe2) = pwcpa(ji,jk,jwfe2) + 8.0 * zreasat * radssol(jk,jwfe2) 234 pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) + 14.0 * zreasat 235 pwcpa(ji,jk,jwso4) = pwcpa(ji,jk,jwso4) + zreasat 236 pwcpa(ji,jk,jwpo4) = pwcpa(ji,jk,jwpo4) + 8.0 * redfep * zreasat 237 238 ! Fe + H2S 239 zh2seq = MAX(rtrn, 2.1E-3 * hipor(ji,jk)) 240 zexcess = pwcp(ji,jk,jwh2s) * zsedfer / zh2seq - 1.0 241 IF ( zexcess >= 0.0 ) THEN 242 zreasat = reac_fesp * zexcess 243 pwcpa (ji,jk,jwfe2) = pwcpa(ji,jk,jwfe2) - zreasat * radssol(jk,jwfe2) 244 solcpa(ji,jk,jsfes) = solcpa(ji,jk,jsfes) + zreasat / volc(ji,jk,jsfes) 245 pwcpa(ji,jk,jwh2s) = pwcpa(ji,jk,jwh2s) - zreasat 246 ELSE 247 zreasat = reac_fesd * zexcess * solcp(ji,jk,jsfes) * volc(ji,jk,jsfes) 248 pwcpa (ji,jk,jwfe2) = pwcpa(ji,jk,jwfe2) - zreasat * radssol(jk,jwfe2) 249 solcpa(ji,jk,jsfes) = solcpa(ji,jk,jsfes) + zreasat / volc(ji,jk,jsfes) 250 pwcpa(ji,jk,jwh2s) = pwcpa(ji,jk,jwh2s) - zreasat 251 ENDIF 252 ! For alkalinity 253 pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) - zreasat * 2.0 254 ENDIF 255 END DO 256 END DO 744 257 745 258 IF( ln_timing ) CALL timing_stop('sed_dsr_redoxb') -
NEMO/trunk/src/TOP/PISCES/SED/seddta.F90
r15237 r15450 8 8 USE sed 9 9 USE sedarr 10 USE par_pisces10 USE sedini 11 11 USE phycst, ONLY : rday 12 12 USE iom … … 20 20 21 21 !! * Module variables 22 REAL(wp) :: rsecday ! number of second per a day23 22 REAL(wp) :: conv2 ! [kg/m2/month]-->[g/cm2/s] ( 1 month has 30 days ) 24 23 … … 26 25 # include "do_loop_substitute.h90" 27 26 # include "domzgr_substitute.h90" 28 !! $Id$ 27 29 28 CONTAINS 30 29 … … 47 46 48 47 !! Arguments 49 INTEGER, INTENT( in) :: kt! time-step50 INTEGER, INTENT( in) :: Kbb, Kmm! time level indices48 INTEGER, INTENT( in ) :: kt ! time-step 49 INTEGER, INTENT( in ) :: Kbb, Kmm ! time level indices 51 50 52 51 !! * Local declarations … … 54 53 55 54 REAL(wp), DIMENSION(jpoce) :: zdtap, zdtag 56 REAL(wp), DIMENSION(jpi,jpj) :: zwsbio4, zwsbio3 55 REAL(wp), DIMENSION(jpi,jpj) :: zwsbio4, zwsbio3, zddust 57 56 REAL(wp) :: zf0, zf1, zf2, zkapp, zratio, zdep 57 REAL(wp) :: zzf0, zf0s, zf0b, zzf1, zf1s, zf1b, zzf2, zf2s, zf2b 58 58 59 59 !---------------------------------------------------------------------- … … 78 78 IF (lwp) WRITE(numsed,*) ' sed_dta : Sediment fields' 79 79 dtsed = rDt_trc 80 rsecday = 60.* 60. * 24.81 ! conv2 = 1.0e+3 / ( 1.0e+4 * rsecday * 30. )82 80 conv2 = 1.0e+3 / 1.0e+4 83 rdtsed(2:jpksed) = dtsed / ( denssol * por1(2:jpksed) )84 81 ENDIF 85 82 … … 87 84 zdtap(:) = 0. 88 85 zdtag(:) = 0. 86 zddust(:,:) = 0.0 89 87 90 88 ! reading variables … … 97 95 ! ----------------------------------------------------------- 98 96 IF (ln_sediment_offline) THEN 99 DO_2D( 1, 1, 1, 1)97 DO_2D( 0, 0, 0, 0 ) 100 98 ikt = mbkt(ji,jj) 101 99 zwsbio4(ji,jj) = wsbio2 / rday … … 103 101 END_2D 104 102 ELSE 105 DO_2D( 1, 1, 1, 1)103 DO_2D( 0, 0, 0, 0 ) 106 104 ikt = mbkt(ji,jj) 107 105 zdep = e3t(ji,jj,ikt,Kmm) / rDt_trc … … 112 110 113 111 trc_data(:,:,:) = 0. 114 DO_2D( 1, 1, 1, 1)112 DO_2D( 0, 0, 0, 0 ) 115 113 ikt = mbkt(ji,jj) 116 IF ( tmask(ji,jj,ikt) == 1 ) THEN 117 trc_data(ji,jj,1) = tr(ji,jj,ikt,jpsil,Kbb) 118 trc_data(ji,jj,2) = tr(ji,jj,ikt,jpoxy,Kbb) 119 trc_data(ji,jj,3) = tr(ji,jj,ikt,jpdic,Kbb) 120 trc_data(ji,jj,4) = tr(ji,jj,ikt,jpno3,Kbb) / 7.625 121 trc_data(ji,jj,5) = tr(ji,jj,ikt,jppo4,Kbb) / 122. 122 trc_data(ji,jj,6) = tr(ji,jj,ikt,jptal,Kbb) 123 trc_data(ji,jj,7) = tr(ji,jj,ikt,jpnh4,Kbb) / 7.625 124 trc_data(ji,jj,8) = 0.0 125 trc_data(ji,jj,9) = 28.0E-3 126 trc_data(ji,jj,10) = tr(ji,jj,ikt,jpfer,Kbb) 127 trc_data(ji,jj,11 ) = MIN(tr(ji,jj,ikt,jpgsi,Kbb), 1E-4) * zwsbio4(ji,jj) * 1E3 128 trc_data(ji,jj,12 ) = MIN(tr(ji,jj,ikt,jppoc,Kbb), 1E-4) * zwsbio3(ji,jj) * 1E3 129 trc_data(ji,jj,13 ) = MIN(tr(ji,jj,ikt,jpgoc,Kbb), 1E-4) * zwsbio4(ji,jj) * 1E3 130 trc_data(ji,jj,14) = MIN(tr(ji,jj,ikt,jpcal,Kbb), 1E-4) * zwsbio4(ji,jj) * 1E3 131 trc_data(ji,jj,15) = ts(ji,jj,ikt,jp_tem,Kmm) 132 trc_data(ji,jj,16) = ts(ji,jj,ikt,jp_sal,Kmm) 133 trc_data(ji,jj,17 ) = ( tr(ji,jj,ikt,jpsfe,Kbb) * zwsbio3(ji,jj) + tr(ji,jj,ikt,jpbfe,Kbb) & 134 & * zwsbio4(ji,jj) ) * 1E3 / ( trc_data(ji,jj,12 ) + trc_data(ji,jj,13 ) + rtrn ) 135 trc_data(ji,jj,17 ) = MIN(1E-3, trc_data(ji,jj,17 ) ) 114 IF ( tmask(ji,jj,ikt) == 1.0 ) THEN 115 trc_data(ji,jj,jwsil) = tr(ji,jj,ikt,jpsil,Kbb) 116 trc_data(ji,jj,jwoxy) = tr(ji,jj,ikt,jpoxy,Kbb) 117 trc_data(ji,jj,jwdic) = tr(ji,jj,ikt,jpdic,Kbb) 118 trc_data(ji,jj,jwno3) = tr(ji,jj,ikt,jpno3,Kbb) * redNo3 / redC 119 trc_data(ji,jj,jwpo4) = tr(ji,jj,ikt,jppo4,Kbb) / redC 120 trc_data(ji,jj,jwalk) = tr(ji,jj,ikt,jptal,Kbb) 121 trc_data(ji,jj,jwnh4) = tr(ji,jj,ikt,jpnh4,Kbb) * redNo3 / redC 122 trc_data(ji,jj,jwh2s) = 0.0 123 trc_data(ji,jj,jwso4) = 0.14 * ts(ji,jj,ikt,jp_sal,Kmm) / 1.80655 / 96.062 124 trc_data(ji,jj,jwfe2) = tr(ji,jj,ikt,jpfer,Kbb) 125 trc_data(ji,jj,jwlgw) = 1E-9 126 trc_data(ji,jj,12 ) = MIN(tr(ji,jj,ikt,jpgsi,Kbb), 1E-4) * zwsbio4(ji,jj) * 1E3 127 trc_data(ji,jj,13 ) = MIN(tr(ji,jj,ikt,jppoc,Kbb), 1E-4) * zwsbio3(ji,jj) * 1E3 128 trc_data(ji,jj,14 ) = MIN(tr(ji,jj,ikt,jpgoc,Kbb), 1E-4) * zwsbio4(ji,jj) * 1E3 129 trc_data(ji,jj,15) = MIN(tr(ji,jj,ikt,jpcal,Kbb), 1E-4) * zwsbio4(ji,jj) * 1E3 130 trc_data(ji,jj,16) = ts(ji,jj,ikt,jp_tem,Kmm) 131 trc_data(ji,jj,17) = ts(ji,jj,ikt,jp_sal,Kmm) 132 trc_data(ji,jj,18 ) = ( tr(ji,jj,ikt,jpsfe,Kbb) * zwsbio3(ji,jj) + tr(ji,jj,ikt,jpbfe,Kbb) & 133 & * zwsbio4(ji,jj) ) * 1E3 / ( trc_data(ji,jj,13 ) + trc_data(ji,jj,14 ) + rtrn ) 134 trc_data(ji,jj,18 ) = MIN(1E-3, trc_data(ji,jj,18 ) ) 136 135 ENDIF 137 136 END_2D … … 142 141 CALL pack_arr ( jpoce, pwcp_dta(1:jpoce,jw), trc_data(1:jpi,1:jpj,jw), iarroce(1:jpoce) ) 143 142 END DO 143 144 144 ! Solid components : 145 145 !----------------------- 146 146 ! Sinking fluxes for OPAL in mol.m-2.s-1 ; conversion in mol.cm-2.s-1 147 CALL pack_arr ( jpoce, rainrm_dta(1:jpoce,jsopal), trc_data(1:jpi,1:jpj,1 1), iarroce(1:jpoce) )147 CALL pack_arr ( jpoce, rainrm_dta(1:jpoce,jsopal), trc_data(1:jpi,1:jpj,12), iarroce(1:jpoce) ) 148 148 rainrm_dta(1:jpoce,jsopal) = rainrm_dta(1:jpoce,jsopal) * 1e-4 149 149 150 ! Sinking fluxes for POC in mol.m-2.s-1 ; conversion in mol.cm-2.s-1 150 CALL pack_arr ( jpoce, zdtap(1:jpoce), trc_data(1:jpi,1:jpj,1 2) , iarroce(1:jpoce) )151 CALL pack_arr ( jpoce, zdtag(1:jpoce), trc_data(1:jpi,1:jpj,1 3) , iarroce(1:jpoce) )151 CALL pack_arr ( jpoce, zdtap(1:jpoce), trc_data(1:jpi,1:jpj,13) , iarroce(1:jpoce) ) 152 CALL pack_arr ( jpoce, zdtag(1:jpoce), trc_data(1:jpi,1:jpj,14) , iarroce(1:jpoce) ) 152 153 DO ji = 1, jpoce 153 ! zkapp = MIN( (1.0 - 0.02 ) * reac_poc, 3731.0 * max(100.0, zkbot(ji) )**(-1.011) / ( 365.0 * 24.0 * 3600.0 ) ) 154 ! zkapp = MIN( 0.98 * reac_poc, 100.0 * max(100.0, zkbot(ji) )**(-0.6) / ( 365.0 * 24.0 * 3600.0 ) ) 155 ! zratio = ( ( 1.0 - 0.02 ) * reac_poc + 0.02 * reac_poc * 0. - zkapp) / ( ( 0.02 - 1.0 ) * reac_poc / 100. - 0.02 * reac_poc * 0. + zkapp ) 156 ! zf1 = ( 0.02 * (reac_poc - reac_poc * 0.) + zkapp - reac_poc ) / ( reac_poc / 100. - reac_poc ) 157 ! zf1 = MIN(0.98, MAX(0., zf1 ) ) 158 zf1 = 0.48 159 zf0 = 1.0 - 0.02 - zf1 160 zf2 = 0.02 161 rainrm_dta(ji,jspoc) = ( zdtap(ji) + zdtag(ji) ) * 1e-4 * zf0 162 rainrm_dta(ji,jspos) = ( zdtap(ji) + zdtag(ji) ) * 1e-4 * zf1 163 rainrm_dta(ji,jspor) = ( zdtap(ji) + zdtag(ji) ) * 1e-4 * zf2 154 zzf2 = 2E-2 155 zzf1 = 0.3 156 zzf0 = 1.0 - zzf1 - zzf2 157 zf0s = zzf0 158 zf1s = zzf1 159 zf2s = 1.0 - zf1s - zf0s 160 zf0b = zzf0 161 zf1b = zzf1 162 zf2b = 1.0 - zf1b - zf0b 163 rainrm_dta(ji,jspoc) = ( zdtap(ji) * zf0s + zdtag(ji) * zf0b ) * 1e-4 164 rainrm_dta(ji,jspos) = ( zdtap(ji) * zf1s + zdtag(ji) * zf1b ) * 1e-4 165 rainrm_dta(ji,jspor) = ( zdtap(ji) * zf2s + zdtag(ji) * zf2b ) * 1e-4 164 166 END DO 167 165 168 ! Sinking fluxes for Calcite in mol.m-2.s-1 ; conversion in mol.cm-2.s-1 166 CALL pack_arr ( jpoce, rainrm_dta(1:jpoce,jscal), trc_data(1:jpi,1:jpj,1 4), iarroce(1:jpoce) )169 CALL pack_arr ( jpoce, rainrm_dta(1:jpoce,jscal), trc_data(1:jpi,1:jpj,15), iarroce(1:jpoce) ) 167 170 rainrm_dta(1:jpoce,jscal) = rainrm_dta(1:jpoce,jscal) * 1e-4 168 ! vector temperature [�C] and salinity 169 CALL pack_arr ( jpoce, temp(1:jpoce), trc_data(1:jpi,1:jpj,15), iarroce(1:jpoce) ) 170 CALL pack_arr ( jpoce, salt(1:jpoce), trc_data(1:jpi,1:jpj,16), iarroce(1:jpoce) ) 171 172 ! vector temperature [°C] and salinity 173 CALL pack_arr ( jpoce, temp(1:jpoce), trc_data(1:jpi,1:jpj,16), iarroce(1:jpoce) ) 174 CALL pack_arr ( jpoce, salt(1:jpoce), trc_data(1:jpi,1:jpj,17), iarroce(1:jpoce) ) 171 175 172 176 ! Clay rain rate in [mol/(cm**2.s)] … … 175 179 CALL pack_arr ( jpoce, rainrm_dta(1:jpoce,jsclay), dust(1:jpi,1:jpj), iarroce(1:jpoce) ) 176 180 rainrm_dta(1:jpoce,jsclay) = rainrm_dta(1:jpoce,jsclay) * conv2 / mol_wgt(jsclay) & 177 & + wacc(1:jpoce) * por1(2) * denssol / mol_wgt(jsclay) / ( rsecday * 365.0 ) 178 rainrm_dta(1:jpoce,jsclay) = rainrm_dta(1:jpoce,jsclay) * 0.965 179 rainrm_dta(1:jpoce,jsfeo) = rainrm_dta(1:jpoce,jsclay) * mol_wgt(jsclay) / mol_wgt(jsfeo) * 0.035 / 0.965 181 & + wacc(1:jpoce) * por1(2) * dens_sol(jsclay) / mol_wgt(jsclay) / ( rday * 365.0 ) 182 rainrm_dta(1:jpoce,jsfeo) = rainrm_dta(1:jpoce,jsclay) * mol_wgt(jsclay) / mol_wgt(jsfeo) * 0.035 * 0.5 183 rainrm_dta(1:jpoce,jsclay) = rainrm_dta(1:jpoce,jsclay) * ( 1.0 - 0.035 * 0.5 ) 184 CALL unpack_arr ( jpoce, zddust(1:jpi,1:jpj), iarroce(1:jpoce), wacc(1:jpoce) ) 185 zddust(:,:) = dust(:,:) + zddust(:,:) / ( rday * 365.0 ) * por1(2) * dens_sol(jsclay) / conv2 186 180 187 ! rainrm_dta(1:jpoce,jsclay) = 1.0E-4 * conv2 / mol_wgt(jsclay) 181 188 … … 184 191 185 192 ! Fe/C ratio in sinking particles that fall to the sediments 186 CALL pack_arr ( jpoce, fecratio(1:jpoce), trc_data(1:jpi,1:jpj,17), iarroce(1:jpoce) ) 187 188 sedligand(:,1) = 1.E-9 193 CALL pack_arr ( jpoce, fecratio(1:jpoce), trc_data(1:jpi,1:jpj,18), iarroce(1:jpoce) ) 189 194 190 195 ! sediment pore water at 1st layer (k=1) 191 DO jw = 1, jpwat 192 pwcp(1:jpoce,1,jw) = pwcp_dta(1:jpoce,jw) 193 ENDDO 194 195 ! rain 196 DO js = 1, jpsol 197 rainrm(1:jpoce,js) = rainrm_dta(1:jpoce,js) 198 ENDDO 196 pwcp(1:jpoce,1,1:jpwat) = pwcp_dta(1:jpoce,1:jpwat) 199 197 200 198 ! Calculation of raintg of each sol. comp.: rainrm in [g/(cm**2.s)] 201 199 DO js = 1, jpsol 202 rainrg(1:jpoce,js) = rainrm (1:jpoce,js) *mol_wgt(js)200 rainrg(1:jpoce,js) = rainrm_dta(1:jpoce,js) * mol_wgt(js) 203 201 ENDDO 204 202 205 ! Calculation of raintg = total massic flux rained in each cell (sum of sol. comp.)206 raintg(:) = 0.203 ! computation of dzdep = total thickness of solid material rained [cm] in each cell 204 dzdep(:) = 0. 207 205 DO js = 1, jpsol 208 raintg(1:jpoce) = raintg(1:jpoce) + rainrg(1:jpoce,js) 209 ENDDO 210 211 ! computation of dzdep = total thickness of solid material rained [cm] in each cell 212 dzdep(1:jpoce) = raintg(1:jpoce) * rdtsed(2) 206 dzdep(1:jpoce) = dzdep(1:jpoce) + rainrg(1:jpoce,js) * dtsed / ( dens_sol(js) * por1(2) ) 207 END DO 213 208 214 209 IF( lk_iomput ) THEN 215 IF( iom_use("sflxclay" ) ) CALL iom_put( "sflxclay", dust(:,:) * conv2 * 1E4 )216 IF( iom_use("sflxcal" ) ) CALL iom_put( "sflxcal", trc_data(:,:,1 3))217 IF( iom_use("sflxbsi" ) ) CALL iom_put( "sflxbsi", trc_data(:,:,1 0))218 IF( iom_use("sflxpoc" ) ) CALL iom_put( "sflxpoc", trc_data(:,:,11) + trc_data(:,:,12))210 IF( iom_use("sflxclay" ) ) CALL iom_put( "sflxclay", zddust(:,:) * 1E3 / 1.E4 ) 211 IF( iom_use("sflxcal" ) ) CALL iom_put( "sflxcal", trc_data(:,:,15) / 1.E4 ) 212 IF( iom_use("sflxbsi" ) ) CALL iom_put( "sflxbsi", trc_data(:,:,12) / 1.E4 ) 213 IF( iom_use("sflxpoc" ) ) CALL iom_put( "sflxpoc", ( trc_data(:,:,13) + trc_data(:,:,14) ) / 1.E4 ) 219 214 ENDIF 220 215 -
NEMO/trunk/src/TOP/PISCES/SED/sedini.F90
r14086 r15450 6 6 7 7 !!---------------------------------------------------------------------- 8 !! sed_ini t: initialization, namelist read, and parameters control8 !! sed_ini : initialization, namelist read, and parameters control 9 9 !!---------------------------------------------------------------------- 10 10 !! * Modules used 11 USE par_trc ! need jptra, number of passive tracers12 USE par_sed13 11 USE sed ! sediment global variable 14 12 USE sed_oce … … 24 22 PRIVATE 25 23 26 !! * Substitutions27 # include "do_loop_substitute.h90"28 24 !! Module variables 25 REAL(wp), PUBLIC :: sedmask 26 29 27 REAL(wp) :: & 30 28 sedzmin = 0.3 , & !: Minimum vertical spacing … … 39 37 40 38 REAL(wp) :: & 41 rcopal = 40. , & !: reactivity for si [l.mol-1.an-1] 42 dcoef = 8.e-6 !: diffusion coefficient (*por) [cm**2/s] 39 rcopal = 40. !: reactivity for si [l.mol-1.an-1] 43 40 44 41 REAL(wp), PUBLIC :: & 45 redO2 = 1 72. , & !: Redfield coef for Oxygen42 redO2 = 138. , & !: Redfield coef for Oxygen 46 43 redNo3 = 16. , & !: Redfield coef for Nitrate 47 44 redPo4 = 1. , & !: Redfield coef for Phosphate 48 redC = 1 22. , & !: Redfield coef for Carbon45 redC = 117. , & !: Redfield coef for Carbon 49 46 redfep = 0.175 , & !: Ratio for iron bound phosphorus 50 47 rcorgl = 50. , & !: reactivity for POC/O2 [l.mol-1.an-1] … … 55 52 rcfe2 = 5.E8 , & !: reactivity for O2/Fe2+ [l.mol-1.an-1] 56 53 rcfeh2s = 1.E4 , & !: Reactivity for FEOH/H2S [l.mol-1.an-1] 57 rcfes = 1.E5 , & !: Reactivity for FE2+/H2S [l.mol-1.an-1]58 54 rcfeso = 3.E5 , & !: Reactivity for FES/O2 [l.mol-1.an-1] 55 rcfesp = 5E-6 , & !: Precipitation of FeS [mol/l-1.an-1] 56 rcfesd = 1E-3 , & !: Dissolution of FeS [an-1] 59 57 xksedo2 = 5E-6 , & !: half-sturation constant for oxic remin. 60 58 xksedno3 = 5E-6 , & !: half-saturation constant for denitrification … … 73 71 74 72 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 / 73 DATA diff1/ 1.104E-5, 9.78E-6, 3.58E-6, 9.8E-6, 9.73E-6, 5.0E-6, 3.31E-6, 4.81E-6, 4.81E-6, 4.81E-6, 4.59E-6 / 74 76 75 77 76 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.5E-7, 3.89E-7, 3.06E-7, 2.5E-7, 1.5E-7 / 79 80 77 DATA diff2/ 4.47E-7, 3.89E-7, 1.77E-7, 3.89E-7, 3.06E-7, 2.5E-7, 1.5E-7, 2.51E-7, 2.51E-7, 2.51E-7, 1.74E-7 / 81 78 82 79 !! * Routine accessibility 83 PUBLIC sed_init ! routine called by opa.F90 84 80 PUBLIC sed_ini ! routine called by opa.F90 81 82 !! * Substitutions 83 # include "do_loop_substitute.h90" 85 84 !! $Id$ 86 85 CONTAINS 87 86 88 87 89 SUBROUTINE sed_ini t88 SUBROUTINE sed_ini 90 89 !!---------------------------------------------------------------------- 91 !! *** ROUTINE sed_ini t***90 !! *** ROUTINE sed_ini *** 92 91 !! 93 92 !! ** Purpose : Initialization of sediment module … … 104 103 !! ! 06-07 (C. Ethe) Re-organization 105 104 !!---------------------------------------------------------------------- 106 INTEGER :: ji, jj, ikt, ierr 105 INTEGER :: ji, jj, js, jn, jk, ikt, ierr 106 REAL(wp) :: ztmp1, ztmp2 107 107 !!---------------------------------------------------------------------- 108 109 108 110 109 ! Reading namelist.sed variables … … 122 121 ENDIF 123 122 124 IF(lwp) WRITE(numsed,*) ' sed_ini t: Initialization of sediment module '123 IF(lwp) WRITE(numsed,*) ' sed_ini : Initialization of sediment module ' 125 124 IF(lwp) WRITE(numsed,*) ' ' 126 125 127 126 ! Read sediment Namelist 128 127 !------------------------- 129 CALL sed_ini t_nam128 CALL sed_ini_nam 130 129 131 130 ! Allocate SEDIMENT arrays 132 131 ierr = sed_alloc() 133 132 ierr = ierr + sed_oce_alloc() 134 ierr = ierr + sed_adv_alloc()135 133 IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'sed_ini: unable to allocate sediment model arrays' ) 136 134 137 135 ! Determination of sediments number of points and allocate global variables 138 136 epkbot(:,:) = 0. 139 DO_2D( 1, 1, 1, 1 ) 137 gdepbot(:,:) = 0. 138 DO_2D( 0, 0, 0, 0 ) 140 139 ikt = mbkt(ji,jj) 141 IF( tmask(ji,jj,ikt) == 1 ) epkbot(ji,jj) = e3t_ 1d(ikt)142 gdepbot(ji,jj) = gdepw_0(ji,jj,ikt )140 IF( tmask(ji,jj,ikt) == 1 ) epkbot(ji,jj) = e3t_0(ji,jj,ikt) 141 gdepbot(ji,jj) = gdepw_0(ji,jj,ikt+1) 143 142 END_2D 144 143 … … 154 153 155 154 ! Allocate memory size of global variables 156 ALLOCATE( pwcp (jpoce,jpksed,jpwat) ) ; ALLOCATE( pwcp0 (jpoce,jpksed,jpwat) ) ; ALLOCATE( pwcp_dta (jpoce,jpwat) ) 157 ALLOCATE( solcp(jpoce,jpksed,jpsol) ) ; ALLOCATE( solcp0(jpoce,jpksed,jpsol) ) ; ALLOCATE( rainrm_dta(jpoce,jpsol) ) 158 ALLOCATE( rainrm(jpoce,jpsol) ) ; ALLOCATE( rainrg(jpoce,jpsol) ) ; ALLOCATE( raintg(jpoce) ) 155 ALLOCATE( pwcp (jpoce,jpksed,jpwat) ) ; ALLOCATE( pwcp_dta (jpoce,jpwat) ) 156 ALLOCATE( pwcpa(jpoce,jpksed,jpwat) ) ; ALLOCATE( solcpa(jpoce,jpksed,jpsol) ) 157 ALLOCATE( solcp(jpoce,jpksed,jpsol) ) ; ALLOCATE( rainrm_dta(jpoce,jpsol) ) 158 ALLOCATE( rainrg(jpoce,jpsol) ) 159 159 ALLOCATE( dzdep(jpoce) ) ; ALLOCATE( iarroce(jpoce) ) ; ALLOCATE( dzkbot(jpoce) ) 160 ALLOCATE( slatit(jpoce) ) ; ALLOCATE( slongit(jpoce) ) 160 161 ALLOCATE( zkbot(jpoce) ) ; ALLOCATE( db(jpoce,jpksed) ) 161 162 ALLOCATE( temp(jpoce) ) ; ALLOCATE( salt(jpoce) ) 162 163 ALLOCATE( diff(jpoce,jpksed,jpwat ) ) ; ALLOCATE( irrig(jpoce, jpksed) ) 163 164 ALLOCATE( wacc(jpoce) ) ; ALLOCATE( fecratio(jpoce) ) 164 ALLOCATE( press(jpoce) ) ; ALLOCATE( densSW(jpoce) )165 ALLOCATE( densSW(jpoce) ) ; ALLOCATE( saturco3(jpoce,jpksed) ) 165 166 ALLOCATE( hipor(jpoce,jpksed) ) ; ALLOCATE( co3por(jpoce,jpksed) ) 166 167 ALLOCATE( dz3d(jpoce,jpksed) ) ; ALLOCATE( volw3d(jpoce,jpksed) ) ; ALLOCATE( vols3d(jpoce,jpksed) ) 167 ALLOCATE( sedligand(jpoce, jpksed) ) 168 ALLOCATE( rearatpom(jpoce, jpksed) ) ; ALLOCATE( volc(jpoce,jpksed,jpsol) ) 169 ALLOCATE( radssol(jpksed, jpwat) ) ; ALLOCATE( rads1sol(jpksed, jpwat) ) 170 ALLOCATE( apluss(jpoce, jpksed) ) ; ALLOCATE( aminuss(jpoce,jpksed) ) 168 171 169 172 ! 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 pwcp (:,:,:) = 0. ; pwcp_dta (:,:) = 0. 174 pwcpa (:,:,:) = 0. ; solcpa(:,:,:) = 0. 175 solcp (:,:,:) = 0. ; rainrm_dta(:,:) = 0. 176 rainrg(:,: ) = 0. 173 177 dzdep (: ) = 0. ; iarroce(: ) = 0 ; dzkbot (: ) = 0. 174 178 temp (: ) = 0. ; salt (: ) = 0. ; zkbot (: ) = 0. 175 press (: ) = 0. ; densSW (: ) = 0. ; db(:,:) = 0.179 densSW (: ) = 0. ; db (:,:) = 0. 176 180 hipor (:,: ) = 0. ; co3por (:,: ) = 0. ; irrig (:,:) = 0. 177 181 dz3d (:,: ) = 0. ; volw3d (:,: ) = 0. ; vols3d (:,:) = 0. 178 fecratio(:) = 1E-5 179 sedligand(:,:) = 0.6E-9 182 fecratio(:) = 1E-5 ; rearatpom(:,:) = 0. 183 radssol(:,:) = 1.0 ; rads1sol(:,:) = 0. 184 apluss(:,:) = 0.0 ; aminuss(:,:) = 0.0 180 185 181 186 ! Chemical variables … … 185 190 ALLOCATE( borats(jpoce) ) ; ALLOCATE( calcon2(jpoce) ) ; ALLOCATE( sieqs (jpoce) ) 186 191 ALLOCATE( aks3s(jpoce) ) ; ALLOCATE( akf3s(jpoce) ) ; ALLOCATE( sulfats(jpoce) ) 187 ALLOCATE( fluorids(jpoce) ) 192 ALLOCATE( fluorids(jpoce) ) ; ALLOCATE( akh2s(jpoce) ) ; ALLOCATE( aknh3(jpoce) ) 188 193 189 194 akbs (:) = 0. ; ak1s (:) = 0. ; ak2s (:) = 0. ; akws (:) = 0. 190 195 ak1ps (:) = 0. ; ak2ps (:) = 0. ; ak3ps (:) = 0. ; aksis (:) = 0. 191 196 aksps (:) = 0. ; ak12s (:) = 0. ; ak12ps(:) = 0. ; ak123ps(:) = 0. 192 borats(:) = 0. ; calcon2(:) = 0. ; sieqs (:) = 0. 197 borats(:) = 0. ; calcon2(:) = 0. ; sieqs (:) = 0. ; akh2s (:) = 0. 193 198 aks3s(:) = 0. ; akf3s(:) = 0. ; sulfats(:) = 0. ; fluorids(:) = 0. 199 aknh3(:) = 0. 194 200 195 201 ! 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. 202 ALLOCATE( fromsed(jpoce, jpsol+jpads) ) ; ALLOCATE( tosed(jpoce, jpsol+jpads) ) 203 204 fromsed(:,:) = 0. ; tosed(:,:) = 0. 200 205 201 206 ! Initialization of sediment geometry 202 207 !------------------------------------ 203 CALL sed_ini t_geom208 CALL sed_ini_geom 204 209 205 210 ! Offline specific mode 206 211 ! --------------------- 207 212 ln_sediment_offline = .FALSE. 213 214 !--------------------------------------------- 215 ! Molecular weight [g/mol] for solid species 216 !--------------------------------------------- 217 218 ! opal=sio2*0.4(h20)=28+2*16+0.4*(2+16) 219 !--------------------------------------- 220 mol_wgt(jsopal) = 28. + 2. * 16. + 0.4 * ( 2. + 16. ) 221 222 ! clay 223 ! some kind of Illit (according to Pape) 224 ! K0.58(Al 1.38 Fe(III)0.37Fe(II)0.04Mg0.34)[(OH)2|(Si3.41Al0.59)O10] 225 !-------------------------------------------------------------------- 226 mol_wgt(jsclay) = 0.58 * 39. + 1.38 * 27. + ( 0.37 + 0.04 ) * 56.+ & 227 & 0.34 * 24. + 2. * ( 16. + 1. ) + 3.41 * 38. + & 228 & 0.59 * 27. + 10. * 16. 229 230 mol_wgt(jsfeo) = 55.0 + 3.0 * ( 16.0 + 1.0) 231 232 233 mol_wgt(jsfes) = 55.0 + 32.0 234 235 ! for chemistry Poc : C(122)H(244)O(86)N(16)P(1) 236 ! But den sity of Poc is an Hydrated material (= POC + 30H2O) 237 ! So C(122)H(355)O(120)N(16)P(1) 238 !------------------------------------------------------------ 239 mol_wgt(jspoc) = ( redC * 12. + 355. + 120. * 16. + redNo3 * 14. + 31. ) / redC 240 mol_wgt(jspos) = mol_wgt(jspoc) 241 mol_wgt(jspor) = mol_wgt(jspoc) 242 243 ! CaCO3 244 !--------- 245 mol_wgt(jscal) = 40. + 12. + 3. * 16. 246 247 ! Density of solid material in sediment [g/cm**3] 248 !------------------------------------------------ 249 ALLOCATE( dens_sol(jpsol) ) 250 dens_sol(jsclay) = 2.6 251 dens_sol(jscal) = 2.7 252 dens_sol(jsopal) = 2.1 253 dens_sol(jsfeo) = 3.4 254 dens_sol(jsfes) = 4.8 255 dens_sol(jspoc) = 1.0 256 dens_sol(jspos) = 1.0 257 dens_sol(jspor) = 1.0 258 259 ! Accumulation rate from Burwicz et al. (2011). This is used to 260 ! compute the flux of clays and minerals 261 ! -------------------------------------------------------------- 262 DO ji = 1, jpoce 263 ztmp1 = 0.0117 * 3.0 / ( 1.0 + ( zkbot(ji) / 200.)**3 ) 264 ztmp2 = 0.006 / ( 1.0 + ( zkbot(ji) / 5000.)**10 ) 265 wacc(ji) = ztmp2+ztmp1 266 END DO 267 268 ! Vertical profile of of the adsorption factor for adsorbed species 269 ! ----------------------------------------------------------------- 270 radssol(:,jwfe2) = 1.0 / ( 1.0 + adsfe2 * por1(:) / por(:) ) 271 radssol(:,jwnh4) = 1.0 / ( 1.0 + adsnh4 * por1(:) / por(:) ) 272 rads1sol(:,jwnh4) = adsnh4 * por1(:) / ( por(:) + adsnh4 * por1(:) ) 273 rads1sol(:,jwfe2) = adsfe2 * por1(:) / ( por(:) + adsfe2 * por1(:) ) 274 275 ! Initialization of the array for non linear solving 276 ! -------------------------------------------------- 277 278 ALLOCATE( jarr(jpksed*jpvode,2) ) 279 ALLOCATE( jsvode(jpvode), isvode(jptrased) ) 280 jsvode(1) = jwoxy ; jsvode(2) = jwno3 ; jsvode(3) = jwnh4 281 jsvode(4) = jwh2s ; jsvode(5) = jwso4 ; jsvode(6) = jwfe2 282 jsvode(7) = jpwat+jsfeo ; jsvode(8) = jpwat+jsfes 283 isvode(jwoxy) = 1 ; isvode(jwno3) = 2 ; isvode(jwnh4) = 3 284 isvode(jwh2s) = 4 ; isvode(jwso4) = 5 ; isvode(jwfe2) = 6 285 isvode(jpwat+jsfeo) = 7 ; isvode(jpwat+jsfes) = 8 286 DO js = 1, jpvode 287 DO jk = 1, jpksed 288 jn = (jk-1) * jpvode + js 289 jarr(jn,1) = jk 290 jarr(jn,2) = jsvode(js) 291 END DO 292 END DO 293 294 ALLOCATE( rstepros(jpoce) ) 208 295 209 296 #if defined key_sed_off … … 218 305 ! -------------------------------- 219 306 IF (ln_sediment_offline) THEN 220 CALL trc_dta_ini(jptra)221 307 CALL trc_dmp_sed_ini 222 308 ENDIF 223 309 224 END SUBROUTINE sed_ini t225 226 SUBROUTINE sed_ini t_geom310 END SUBROUTINE sed_ini 311 312 SUBROUTINE sed_ini_geom 227 313 !!---------------------------------------------------------------------- 228 !! *** ROUTINE sed_ini t_geom ***314 !! *** ROUTINE sed_ini_geom *** 229 315 !! 230 316 !! ** Purpose : Initialization of sediment geometry … … 244 330 !---------------------------------------------------------- 245 331 246 IF(lwp) WRITE(numsed,*) ' sed_ini t_geom : Initialization of sediment geometry '332 IF(lwp) WRITE(numsed,*) ' sed_ini_geom : Initialization of sediment geometry ' 247 333 IF(lwp) WRITE(numsed,*) ' ' 248 334 249 335 ! Computation of 1D array of sediments points 250 336 indoce = 0 251 DO_2D( 1, 1, 1, 1)337 DO_2D( 0, 0, 0, 0 ) 252 338 IF ( epkbot(ji,jj) > 0. ) THEN 253 339 indoce = indoce + 1 … … 272 358 !------------------------------------------------ 273 359 CALL pack_arr ( jpoce, dzkbot(1:jpoce), epkbot(1:jpi,1:jpj), iarroce(1:jpoce) ) 360 #if defined key_sed_off 361 dzkbot(1:jpoce) = 1.E8 362 #else 274 363 dzkbot(1:jpoce) = dzkbot(1:jpoce) * 1.e+2 364 #endif 275 365 CALL pack_arr ( jpoce, zkbot(1:jpoce), gdepbot(1:jpi,1:jpj), iarroce(1:jpoce) ) 366 CALL pack_arr ( jpoce, slatit(1:jpoce), gphit(1:jpi,1:jpj), iarroce(1:jpoce) ) 367 CALL pack_arr ( jpoce, slongit(1:jpoce), glamt(1:jpi,1:jpj), iarroce(1:jpoce) ) 276 368 277 369 ! Geometry and constants … … 285 377 zsur = - za0 - za1 * sedacr * LOG( COSH( (1-sedkth) / sedacr ) ) 286 378 379 dz(1) = 0.1 287 380 profsedw(1) = 0.0 288 profsed(1) = -dz(1) / 2.381 profsed(1) = -dz(1) / 2. 289 382 DO jk = 2, jpksed 290 383 zw = REAL( jk , wp ) … … 292 385 profsed(jk) = ( zsur + za0 * zt + za1 * sedacr * LOG ( COSH( (zt-sedkth) / sedacr ) ) ) 293 386 profsedw(jk) = ( zsur + za0 * zw + za1 * sedacr * LOG ( COSH( (zw-sedkth) / sedacr ) ) ) 294 END DO295 296 dz(1) = 0.1297 DO jk = 2, jpksed298 387 dz(jk) = profsedw(jk) - profsedw(jk-1) 299 388 END DO 300 389 301 DO jk = 1, jpksed 302 DO ji = 1, jpoce 303 dz3d(ji,jk) = dz(jk) 304 END DO 305 ENDDO 390 DO ji = 1, jpoce 391 dz3d(ji,:) = dz(:) 392 END DO 306 393 307 394 ! Porosity profile [0] … … 309 396 por(1) = 1.0 310 397 DO jk = 2, jpksed 311 por(jk) = porinf + ( porsurf-porinf) * exp(-rhox * (profsed(jk)) )398 por(jk) = porinf + ( porsurf-porinf) * exp(-rhox * profsed(jk) ) 312 399 END DO 313 400 … … 333 420 dz3d(1:jpoce,1) = dz(1) 334 421 335 !--------------------------------------------- 336 ! Molecular weight [g/mol] for solid species 337 !--------------------------------------------- 338 339 ! opal=sio2*0.4(h20)=28+2*16+0.4*(2+16) 340 !--------------------------------------- 341 mol_wgt(jsopal) = 28. + 2. * 16. + 0.4 * ( 2. + 16. ) 342 343 ! clay 344 ! some kind of Illit (according to Pape) 345 ! K0.58(Al 1.38 Fe(III)0.37Fe(II)0.04Mg0.34)[(OH)2|(Si3.41Al0.59)O10] 346 !-------------------------------------------------------------------- 347 mol_wgt(jsclay) = 0.58 * 39. + 1.38 * 27. + ( 0.37 + 0.04 ) * 56.+ & 348 & 0.34 * 24. + 2. * ( 16. + 1. ) + 3.41 * 38. + & 349 & 0.59 * 27. + 10. * 16. 350 351 mol_wgt(jsfeo) = 55.0 + 3.0 * ( 16.0 + 1.0) 352 353 mol_wgt(jsfes) = 55.0 + 32.0 354 355 ! for chemistry Poc : C(122)H(244)O(86)N(16)P(1) 356 ! But den sity of Poc is an Hydrated material (= POC + 30H2O) 357 ! So C(122)H(355)O(120)N(16)P(1) 358 !------------------------------------------------------------ 359 mol_wgt(jspoc) = ( 122. * 12. + 355. + 120. * 16.+ & 360 & 16. * 14. + 31. ) / 122. 361 mol_wgt(jspos) = mol_wgt(jspoc) 362 mol_wgt(jspor) = mol_wgt(jspoc) 363 364 ! CaCO3 365 !--------- 366 mol_wgt(jscal) = 40. + 12. + 3. * 16. 367 368 ! Density of solid material in sediment [g/cm**3] 369 !------------------------------------------------ 370 denssol = 2.6 371 372 ! Initialization of diffusion coefficient as function of porosity [cm**2/s] 373 !-------------------------------------------------------------------- 374 ! DO jn = 1, jpsol 375 ! DO jk = 1, jpksed 376 ! DO ji = 1, jpoce 377 ! diff(ji,jk,jn) = dcoef / ( 1.0 - 2.0 * log(por(jk)) ) 378 ! END DO 379 ! END DO 380 ! END DO 381 382 ! Accumulation rate from Burwicz et al. (2011). This is used to 383 ! compute the flux of clays and minerals 384 ! -------------------------------------------------------------- 385 DO ji = 1, jpoce 386 ztmp1 = 0.117 / ( 1.0 + ( zkbot(ji) / 200.)**3 ) 387 ztmp2 = 0.006 / ( 1.0 + ( zkbot(ji) / 4000.)**10 ) 388 wacc(ji) = ztmp1 + ztmp2 389 END DO 390 391 392 ! Initialization of time step as function of porosity [cm**2/s] 393 !------------------------------------------------------------------ 394 END SUBROUTINE sed_init_geom 395 396 SUBROUTINE sed_init_nam 422 END SUBROUTINE sed_ini_geom 423 424 SUBROUTINE sed_ini_nam 397 425 !!---------------------------------------------------------------------- 398 !! *** ROUTINE sed_ini t_nam ***426 !! *** ROUTINE sed_ini_nam *** 399 427 !! 400 428 !! ** Purpose : Initialization of sediment geometry … … 405 433 !!---------------------------------------------------------------------- 406 434 407 CHARACTER(:), ALLOCATABLE :: numnamsed_ref !! Character buffer for referencenamelist sediment408 CHARACTER(:), ALLOCATABLE :: numnamsed_cfg !! Character buffer for configurationnamelist sediment435 INTEGER :: numnamsed_ref = -1 !! Logical units for namelist sediment 436 INTEGER :: numnamsed_cfg = -1 !! Logical units for namelist sediment 409 437 INTEGER :: ios ! Local integer output status for namelist read 410 438 CHARACTER(LEN=20) :: clname … … 421 449 TYPE(PSED) , DIMENSION(jpdia2dsed) :: seddiag2d 422 450 423 NAMELIST/nam_run/ nrseddt,ln_sed_2way451 NAMELIST/nam_run/ln_sed_2way,nrosorder,rosatol,rosrtol 424 452 NAMELIST/nam_geom/jpksed, sedzmin, sedhmax, sedkth, sedacr, porsurf, porinf, rhox 425 453 NAMELIST/nam_trased/sedsol, sedwat 426 454 NAMELIST/nam_diased/seddiag3d, seddiag2d 427 NAMELIST/nam_inorg/rcopal, dcoef,rccal, ratligc, rcligc455 NAMELIST/nam_inorg/rcopal, rccal, ratligc, rcligc 428 456 NAMELIST/nam_poc/redO2, redNo3, redPo4, redC, redfep, rcorgl, rcorgs, & 429 & rcorgr, rcnh4, rch2s, rcfe2, rcfeh2s, rcfes , rcfeso,&430 & xksedo2, xksedno3, xksedfeo, xksedso4431 NAMELIST/nam_btb/dbiot, ln_btbz, dbtbzsc, adsnh4, ln_irrig, xirrzsc457 & rcorgr, rcnh4, rch2s, rcfe2, rcfeh2s, rcfeso, rcfesp, & 458 & rcfesd, xksedo2, xksedno3, xksedfeo, xksedso4 459 NAMELIST/nam_btb/dbiot, ln_btbz, dbtbzsc, adsnh4, adsfe2, ln_irrig, xirrzsc 432 460 NAMELIST/nam_rst/ln_rst_sed, cn_sedrst_indir, cn_sedrst_outdir, cn_sedrst_in, cn_sedrst_out 433 461 … … 435 463 !------------------------------------------------------- 436 464 437 IF(lwp) WRITE(numsed,*) ' sed_ini t_nam : Read namelists '465 IF(lwp) WRITE(numsed,*) ' sed_ini_nam : Read namelists ' 438 466 IF(lwp) WRITE(numsed,*) ' ' 439 467 … … 449 477 !--------------------------------- 450 478 clname = 'namelist_sediment' 451 IF(lwp) WRITE(numsed,*) ' sed_ini t_nam : read SEDIMENT namelist'479 IF(lwp) WRITE(numsed,*) ' sed_ini_nam : read SEDIMENT namelist' 452 480 IF(lwp) WRITE(numsed,*) ' ~~~~~~~~~~~~~~' 453 CALL load_nml( numnamsed_ref, TRIM( clname )//'_ref', numout, lwm)454 CALL load_nml( numnamsed_cfg, TRIM( clname )//'_cfg', numout, lwm)481 CALL ctl_opn( numnamsed_ref, TRIM( clname )//'_ref', 'OLD' , 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. ) 482 CALL ctl_opn( numnamsed_cfg, TRIM( clname )//'_cfg', 'OLD' , 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. ) 455 483 456 484 nitsed000 = nittrc000 457 485 nitsedend = nitend 458 486 ! Namelist nam_run 487 REWIND( numnamsed_ref ) ! Namelist nam_run in reference namelist : Pisces variables 459 488 READ ( numnamsed_ref, nam_run, IOSTAT = ios, ERR = 901) 460 489 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_run in reference namelist' ) 461 490 491 REWIND( numnamsed_cfg ) ! Namelist nam_run in reference namelist : Pisces variables 462 492 READ ( numnamsed_cfg, nam_run, IOSTAT = ios, ERR = 902) 463 493 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_run in configuration namelist' ) … … 465 495 IF (lwp) THEN 466 496 WRITE(numsed,*) ' namelist nam_run' 467 WRITE(numsed,*) ' Nb of iterations for fast species nrseddt = ', nrseddt468 497 WRITE(numsed,*) ' 2-way coupling between PISCES and Sed ln_sed_2way = ', ln_sed_2way 498 WRITE(numsed,*) ' Order of the Rosenbrock method (2,3,4) = ', nrosorder 499 WRITE(numsed,*) ' Tolerance for absolute error = ', rosatol 500 WRITE(numsed,*) ' Tolerance for relative order = ', rosrtol 469 501 ENDIF 470 502 471 503 IF ( ln_p5z .AND. ln_sed_2way ) CALL ctl_stop( '2 ways coupling with sediment cannot be activated with PISCES-QUOTA' ) 472 504 505 REWIND( numnamsed_ref ) ! Namelist nam_geom in reference namelist : Pisces variables 473 506 READ ( numnamsed_ref, nam_geom, IOSTAT = ios, ERR = 903) 474 507 903 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_geom in reference namelist' ) 475 508 509 REWIND( numnamsed_cfg ) ! Namelist nam_geom in reference namelist : Pisces variables 476 510 READ ( numnamsed_cfg, nam_geom, IOSTAT = ios, ERR = 904) 477 511 904 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_geom in configuration namelist' ) … … 492 526 dtsed = rDt_trc 493 527 528 REWIND( numnamsed_ref ) ! Namelist nam_trased in reference namelist : Pisces variables 494 529 READ ( numnamsed_ref, nam_trased, IOSTAT = ios, ERR = 905) 495 530 905 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_trased in reference namelist' ) 496 531 532 REWIND( numnamsed_cfg ) ! Namelist nam_trased in reference namelist : Pisces variables 497 533 READ ( numnamsed_cfg, nam_trased, IOSTAT = ios, ERR = 906) 498 534 906 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_trased in configuration namelist' ) … … 523 559 ENDIF 524 560 561 REWIND( numnamsed_ref ) ! Namelist nam_diased in reference namelist : Pisces variables 525 562 READ ( numnamsed_ref, nam_diased, IOSTAT = ios, ERR = 907) 526 563 907 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diased in reference namelist' ) 527 564 565 REWIND( numnamsed_cfg ) ! Namelist nam_diased in reference namelist : Pisces variables 528 566 READ ( numnamsed_cfg, nam_diased, IOSTAT = ios, ERR = 908) 529 567 908 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diased in configuration namelist' ) … … 563 601 ! Inorganic chemistry parameters 564 602 !---------------------------------- 603 REWIND( numnamsed_ref ) ! Namelist nam_inorg in reference namelist : Pisces variables 565 604 READ ( numnamsed_ref, nam_inorg, IOSTAT = ios, ERR = 909) 566 605 909 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_inorg in reference namelist' ) 567 606 607 REWIND( numnamsed_cfg ) ! Namelist nam_inorg in reference namelist : Pisces variables 568 608 READ ( numnamsed_cfg, nam_inorg, IOSTAT = ios, ERR = 910) 569 609 910 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_inorg in configuration namelist' ) … … 572 612 WRITE(numsed,*) ' namelist nam_inorg' 573 613 WRITE(numsed,*) ' reactivity for Si rcopal = ', rcopal 574 WRITE(numsed,*) ' diff. coef for por. dcoef = ', dcoef575 614 WRITE(numsed,*) ' reactivity for calcite rccal = ', rccal 576 615 WRITE(numsed,*) ' L/C ratio in POC ratligc = ', ratligc … … 587 626 ! Additional parameter linked to POC/O2/No3/Po4 588 627 !---------------------------------------------- 628 REWIND( numnamsed_ref ) ! Namelist nam_poc in reference namelist : Pisces variables 589 629 READ ( numnamsed_ref, nam_poc, IOSTAT = ios, ERR = 911) 590 630 911 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_poc in reference namelist' ) 591 631 632 REWIND( numnamsed_cfg ) ! Namelist nam_poc in reference namelist : Pisces variables 592 633 READ ( numnamsed_cfg, nam_poc, IOSTAT = ios, ERR = 912) 593 634 912 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_poc in configuration namelist' ) … … 607 648 WRITE(numsed,*) ' reactivity for Fe2+ rcfe2 = ', rcfe2 608 649 WRITE(numsed,*) ' reactivity for FeOH/H2S rcfeh2s = ', rcfeh2s 609 WRITE(numsed,*) ' reactivity for Fe2+/H2S rcfes = ', rcfes610 650 WRITE(numsed,*) ' reactivity for FeS/O2 rcfeso = ', rcfeso 651 WRITE(numsed,*) ' Precipitation of FeS rcfesp = ', rcfesp 652 WRITE(numsed,*) ' Dissolution of FeS rcfesd = ', rcfesd 611 653 WRITE(numsed,*) ' Half-sat. cste for oxic remin xksedo2 = ', xksedo2 612 654 WRITE(numsed,*) ' Half-sat. cste for denit. xksedno3 = ', xksedno3 … … 629 671 reac_fe2 = rcfe2 / ryear 630 672 reac_feh2s = rcfeh2s/ ryear 631 reac_fes = rcfes / ryear632 673 reac_feso = rcfeso / ryear 674 reac_fesp = rcfesp / ryear 675 reac_fesd = rcfesd / ryear 633 676 634 677 ! reactivity rc in [l.mol-1.s-1] … … 637 680 ! Bioturbation parameter 638 681 !------------------------ 682 REWIND( numnamsed_ref ) ! Namelist nam_btb in reference namelist : Pisces variables 639 683 READ ( numnamsed_ref, nam_btb, IOSTAT = ios, ERR = 913) 640 684 913 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_btb in reference namelist' ) 641 685 686 REWIND( numnamsed_cfg ) ! Namelist nam_btb in reference namelist : Pisces variables 642 687 READ ( numnamsed_cfg, nam_btb, IOSTAT = ios, ERR = 914) 643 688 914 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_btb in configuration namelist' ) … … 649 694 WRITE(numsed,*) ' coefficient for btb attenuation dbtbzsc = ', dbtbzsc 650 695 WRITE(numsed,*) ' Adsorption coefficient of NH4 adsnh4 = ', adsnh4 696 WRITE(numsed,*) ' Adsorption coefficient of Fe2 adsfe2 = ', adsfe2 651 697 WRITE(numsed,*) ' Bioirrigation in sediment ln_irrig = ', ln_irrig 652 698 WRITE(numsed,*) ' coefficient for irrig attenuation xirrzsc = ', xirrzsc … … 656 702 ! Initial value (t=0) for sediment pore water and solid components 657 703 !---------------------------------------------------------------- 704 REWIND( numnamsed_ref ) ! Namelist nam_rst in reference namelist : Pisces variables 658 705 READ ( numnamsed_ref, nam_rst, IOSTAT = ios, ERR = 915) 659 706 915 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_rst in reference namelist' ) 660 707 708 REWIND( numnamsed_cfg ) ! Namelist nam_rst in reference namelist : Pisces variables 661 709 READ ( numnamsed_cfg, nam_rst, IOSTAT = ios, ERR = 916) 662 710 916 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_rst in configuration namelist' ) … … 667 715 WRITE(numsed,*) ' ' 668 716 ENDIF 669 nn_dtsed = 1 670 671 672 END SUBROUTINE sed_init_nam 717 718 CLOSE( numnamsed_cfg ) 719 CLOSE( numnamsed_ref ) 720 721 END SUBROUTINE sed_ini_nam 673 722 674 723 END MODULE sedini -
NEMO/trunk/src/TOP/PISCES/SED/sedinitrc.F90
r12377 r15450 10 10 !! * Modules used 11 11 USE sed ! sediment global variable 12 USE sed_oce13 12 USE sedini 14 13 USE seddta … … 17 16 USE sedchem 18 17 USE sedarr 18 USE sed_oce 19 19 USE lib_mpp ! distribued memory computing library 20 20 … … 50 50 !! ! 06-07 (C. Ethe) Re-organization 51 51 !!---------------------------------------------------------------------- 52 INTEGER, INTENT(in) :: Kbb, Kmm ! time level indices 52 INTEGER, INTENT(in) :: Kbb, Kmm ! time level indices 53 53 54 INTEGER :: ji, jj, ikt 54 55 !!---------------------------------------------------------------------- … … 86 87 !! ! 06-07 (C. Ethe) original 87 88 !!---------------------------------------------------------------------- 88 INTEGER, INTENT(in) :: Kbb, Kmm ! time level indices89 89 90 INTEGER, INTENT(in) :: Kbb, Kmm ! time level indices 91 90 92 ! local variables 91 INTEGER :: & 92 ji, jk, zhipor 93 INTEGER :: ji, jk, zhipor 93 94 94 95 !-------------------------------------------------------------------- … … 134 135 ! Initialization of chemical constants 135 136 CALL sed_chem ( nitsed000 ) 136 137 ! Stores initial sediment data for mass balance calculation138 pwcp0 (1:jpoce,1:jpksed,1:jpwat ) = pwcp (1:jpoce,1:jpksed,1:jpwat )139 solcp0(1:jpoce,1:jpksed,1:jpsol ) = solcp(1:jpoce,1:jpksed,1:jpsol)140 137 141 138 ! Conversion of [h+] in mol/Kg to get it in mol/l ( multiplication by density) -
NEMO/trunk/src/TOP/PISCES/SED/sedinorg.F90
r12839 r15450 8 8 USE sed ! sediment global variable 9 9 USE sed_oce 10 USE sedmat ! linear system of equations11 USE sedco3 ! carbonate ion and proton concentration12 10 USE sedini 13 USE sed dsr11 USE sedmat 14 12 USE lib_mpp ! distribued memory computing library 15 13 USE lib_fortran … … 23 21 CONTAINS 24 22 25 SUBROUTINE sed_inorg( kt ) 23 SUBROUTINE sed_inorg( kt ) 26 24 !!---------------------------------------------------------------------- 27 25 !! *** ROUTINE sed_inorg *** … … 46 44 !!---------------------------------------------------------------------- 47 45 !! Arguments 48 INTEGER, INTENT(in) :: kt ! number of iteration46 INTEGER, INTENT(in) :: kt ! time step 49 47 ! --- local variables 50 INTEGER :: ji, jk, js, jw ! dummy looop indices 51 REAL(wp), DIMENSION(jpoce,jpksed) :: zrearat1, zrearat2 ! reaction rate in pore water 52 REAL(wp), DIMENSION(jpoce,jpksed) :: zundsat ! undersaturation ; indice jpwatp1 is for calcite 53 REAL(wp), DIMENSION(jpoce) :: zco3eq 54 REAL(wp), DIMENSION(jpoce,jpksed,jpsol) :: zvolc ! temp. variables 55 REAL(wp), DIMENSION(jpoce) :: zsieq 56 REAL(wp) :: zsolid1, zvolw, zreasat 57 REAL(wp) :: zsatur, zsatur2, znusil, zsolcpcl, zsolcpsi 48 INTEGER :: ji,jk ! dummy looop indices 49 REAL(wp), DIMENSION(jpoce) :: zsieq, reac_silf 50 REAL(wp) :: zsolid1, zreasat, zco3sat 51 REAL(wp) :: zsatur, zsatur2, znusil, zsolcpcl, zsolcpsi, zexcess 52 REAL(wp), DIMENSION(jpoce,jpksed) :: zundsat, zrearat, psms 58 53 !! 59 54 !!---------------------------------------------------------------------- 60 55 61 56 IF( ln_timing ) CALL timing_start('sed_inorg') 57 58 IF( kt == nitsed000 ) THEN 59 IF (lwp) WRITE(numsed,*) ' sed_inorg : Dissolution of CaCO3 and BSi ' 60 IF (lwp) WRITE(numsed,*) ' ' 61 ENDIF 62 62 ! 63 IF( kt == nitsed000 ) THEN 64 IF (lwp) THEN 65 WRITE(numsed,*) ' sed_inorg : Dissolution reaction ' 66 WRITE(numsed,*) ' ' 67 ENDIF 68 ! ! 69 ENDIF 63 DO ji = 1, jpoce 64 ! ----------------------------------------------- 65 ! Computation of Si solubility 66 ! Param of Ridgwell et al. 2002 67 ! ----------------------------------------------- 70 68 71 ! Initializations72 !----------------------73 74 zrearat1(:,:) = 0. ; zundsat(:,:) = 0.75 zrearat2(:,:) = 0. ; zrearat2(:,:) = 0.76 zco3eq(:) = rtrn77 zvolc(:,:,:) = 0.78 79 ! -----------------------------------------------80 ! Computation of Si solubility81 ! Param of Ridgwell et al. 200282 ! -----------------------------------------------83 84 DO ji = 1, jpoce85 69 zsolcpcl = 0.0 86 70 zsolcpsi = 0.0 87 71 DO jk = 1, jpksed 88 zsolcpsi = zsolcpsi + solcp(ji,jk,jsopal) * dz(jk)89 zsolcpcl = zsolcpcl + solcp(ji,jk,jsclay) * dz(jk)72 zsolcpsi = zsolcpsi + solcp(ji,jk,jsopal) * vols3d(ji,jk) 73 zsolcpcl = zsolcpcl + solcp(ji,jk,jsclay) * vols3d(ji,jk) 90 74 END DO 91 75 zsolcpsi = MAX( zsolcpsi, rtrn ) 92 zsieq(ji) = sieqs(ji) * MAX(0.2 5, 1.0 - (0.045 * zsolcpcl / zsolcpsi )**0.58 )93 zsieq(ji) = MAX( rtrn, sieqs(ji) )76 zsieq(ji) = sieqs(ji) * MAX(0.2, 1.0 - (0.045 * zsolcpcl / zsolcpsi )**0.58 ) 77 reac_silf(ji) = reac_sil * ( 0.05 + 0.055 * ( 1.64 * ( zsolcpcl / zsolcpsi + 0.01 ) )**(-0.75) ) / 1.25 94 78 END DO 95 79 96 DO js = 1, jpsol97 DO jk = 1, jpksed98 DO ji = 1, jpoce99 zvolc(ji,jk,js) = ( vols3d(ji,jk) * dens_mol_wgt(js) ) / &100 & ( volw3d(ji,jk) * 1.e-3 )101 ENDDO102 ENDDO103 ENDDO104 105 !----------------------------------------------------------106 ! 5. Beginning of Pore Water diffusion and solid reaction107 !---------------------------------------------------------108 80 109 !-----------------------------------------------------------------------------110 ! For jk=2,jpksed, and for couple111 ! 1 : jwsil/jsopal ( SI/Opal )112 ! 2 : jsclay/jsclay ( clay/clay )113 ! 3 : jwoxy/jspoc ( O2/POC )114 ! reaction rate is a function of solid=concentration in solid reactif in [mol/l]115 ! and undersaturation in [mol/l].116 ! Solid weight fractions should be in ie [mol/l])117 ! second member and solution are in zundsat variable118 !-------------------------------------------------------------------------119 120 DO jk = 1, jpksed121 DO ji = 1, jpoce122 ! For Silicic Acid and clay123 zundsat(ji,jk) = zsieq(ji) - pwcp(ji,jk,jwsil)124 ENDDO125 ENDDO126 127 ! Definition of reaction rates [rearat]=sans dim128 ! For jk=1 no reaction (pure water without solid) for each solid compo129 81 DO ji = 1, jpoce 130 zrearat1(ji,:) = 0. 131 zrearat2(ji,:) = 0. 132 ENDDO 133 134 ! left hand side of coefficient matrix 135 DO jk = 2, jpksed 136 DO ji = 1, jpoce 137 zsolid1 = zvolc(ji,jk,jsopal) * solcp(ji,jk,jsopal) 138 zsatur = MAX(0., zundsat(ji,jk) / zsieq(ji) ) 82 DO jk = 2, jpksed 83 zsolid1 = volc(ji,jk,jsopal) * solcp(ji,jk,jsopal) 84 zsatur = MAX(0., ( zsieq(ji) - pwcp(ji,jk,jwsil) ) / zsieq(ji) ) 139 85 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) 161 ENDDO 162 ENDDO 86 znusil = ( 0.225 * ( 1.0 + temp(ji) / 15.) * zsatur + 0.775 * zsatur2 * zsatur**9.25 ) 87 solcp(ji,jk,jsopal) = solcp(ji,jk,jsopal) - reac_silf(ji) * znusil * dtsed * solcp(ji,jk,jsopal) 88 pwcp(ji,jk,jwsil) = pwcp(ji,jk,jwsil) + reac_silf(ji) * znusil * dtsed * zsolid1 89 END DO 90 END DO 163 91 164 92 !--------------------------------------------------------------- … … 167 95 168 96 ! computes co3por from the updated pwcp concentrations (note [co3por] = mol/l) 169 170 CALL sed_co3( kt )171 172 97 ! *densSW(l)**2 converts aksps [mol2/kg sol2] into [mol2/l2] to get [undsat] in [mol/l] 173 DO jk = 1, jpksed 174 DO ji = 1, jpoce 175 zco3eq(ji) = aksps(ji) * densSW(ji) * densSW(ji) / ( calcon2(ji) + rtrn ) 176 zco3eq(ji) = MAX( rtrn, zco3eq(ji) ) 177 zundsat(ji,jk) = MAX(0., zco3eq(ji) - co3por(ji,jk) ) 178 ENDDO 179 ENDDO 180 181 DO jk = 2, jpksed 182 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) ) 98 DO ji = 1, jpoce 99 zco3sat = aksps(ji) * densSW(ji) * densSW(ji) / ( calcon2(ji) + rtrn ) 100 saturco3(ji,:) = 1.0 - co3por(ji,:) / ( rtrn + zco3sat ) 101 DO jk = 1, jpksed 102 zsolid1 = volc(ji,jk,jscal) * solcp(ji,jk,jscal) 103 zundsat(ji,jk) = MAX( 0., zco3sat - co3por(ji,jk) ) 104 zrearat(ji,jk) = ( reac_cal * zsolid1 / ( zco3sat + rtrn ) ) / & 105 & ( 1. + reac_cal * dtsed * zundsat(ji,jk) / ( zco3sat + rtrn ) ) 186 106 END DO 187 107 END DO 188 108 109 psms(:,:) = 0.0 189 110 ! solves tridiagonal system 190 CALL sed_mat ( jwdic, jpoce, jpksed, zrearat1, zrearat2, zundsat, dtsed)111 CALL sed_mat_dsri( jpksed, jwdic, -zrearat(:,:), psms(:,:), dtsed, zundsat ) 191 112 192 113 ! New solid concentration values (jk=2 to jksed) for cacO3 193 114 DO jk = 2, jpksed 194 115 DO ji = 1, jpoce 195 zreasat = zrearat 1(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jscal)116 zreasat = zrearat(ji,jk) * dtsed * zundsat(ji,jk) / ( volc(ji,jk,jscal) + rtrn ) 196 117 solcp(ji,jk,jscal) = solcp(ji,jk,jscal) - zreasat 118 zreasat = zrearat(ji,jk) * dtsed * zundsat(ji,jk) 119 ! For DIC 120 pwcp(ji,jk,jwdic) = pwcp(ji,jk,jwdic) + zreasat 121 ! For alkalinity 122 pwcp(ji,jk,jwalk) = pwcp(ji,jk,jwalk) + 2.* zreasat 197 123 ENDDO 198 124 ENDDO 199 125 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 126 297 127 IF( ln_timing ) CALL timing_stop('sed_inorg') -
NEMO/trunk/src/TOP/PISCES/SED/sedmat.F90
r10222 r15450 14 14 PRIVATE 15 15 16 PUBLIC sed_mat 17 18 INTERFACE sed_mat 19 MODULE PROCEDURE sed_mat_dsr, sed_mat_btb 20 END INTERFACE 21 22 INTEGER, PARAMETER :: nmax = 30 23 16 PUBLIC sed_mat_dsr 17 PUBLIC sed_mat_dsrjac 18 PUBLIC sed_mat_dsri 19 PUBLIC sed_mat_btb 20 PUBLIC sed_mat_btbjac 21 PUBLIC sed_mat_btbi 22 PUBLIC sed_mat_coef 24 23 25 24 !! $Id$ 26 25 CONTAINS 27 26 28 SUBROUTINE sed_mat_ dsr( nvar, ndim, nlev, preac, psms, psol, dtsed_in)27 SUBROUTINE sed_mat_coef( nksed ) 29 28 !!--------------------------------------------------------------------- 30 !! *** ROUTINE sed_mat_ dsr***29 !! *** ROUTINE sed_mat_coef *** 31 30 !! 32 31 !! ** Purpose : solves tridiagonal system of linear equations … … 49 48 !!---------------------------------------------------------------------- 50 49 !! * Arguments 51 INTEGER , INTENT(in) :: nvar ! number of variable 52 INTEGER , INTENT(in) :: ndim ! number of points 53 INTEGER , INTENT(in) :: nlev ! number of sediment levels 54 55 REAL(wp), DIMENSION(ndim,nlev), INTENT(in ) :: preac ! reaction rates 56 REAL(wp), DIMENSION(ndim,nlev), INTENT(in ) :: psms ! reaction rates 57 REAL(wp), DIMENSION(ndim,nlev), INTENT(inout) :: psol ! solution ( undersaturation values ) 58 REAL(wp), INTENT(in) :: dtsed_in 59 50 INTEGER, INTENT(in) :: nksed 51 60 52 !---Local declarations 61 INTEGER :: ji, jk, jn 62 REAL(wp), DIMENSION(ndim,nlev) :: za, zb, zc, zr 63 REAL(wp), DIMENSION(ndim) :: zbet 64 REAL(wp), DIMENSION(ndim,nmax) :: zgamm 65 66 REAL(wp) :: aplus,aminus 67 REAL(wp) :: rplus,rminus 68 REAL(wp) :: dxplus,dxminus 69 53 INTEGER :: ji, jk 54 REAL(wp) :: aplus, aminus, dxplus, dxminus 70 55 !---------------------------------------------------------------------- 71 56 72 IF( ln_timing ) CALL timing_start('sed_mat_ dsr')57 IF( ln_timing ) CALL timing_start('sed_mat_coef') 73 58 74 59 ! Computation left hand side of linear system of 75 60 ! equations for dissolution reaction 76 61 !--------------------------------------------- 77 78 62 ! first sediment level 63 DO ji = 1, jpoce 64 aplus = ( por(1) + por(2) ) * 0.5 65 dxplus = ( dz3d(ji,1) + dz3d(ji,2) ) / 2. 66 apluss(ji,1) = ( 1.0 / ( volw3d(ji,1) ) ) * aplus / dxplus 67 68 DO jk = 2, nksed - 1 69 aminus = ( por(jk-1) + por(jk) ) * 0.5 70 dxminus = ( dz3d(ji,jk-1) + dz3d(ji,jk) ) / 2. 71 72 aplus = ( por(jk+1) + por(jk) ) * 0.5 73 dxplus = ( dz3d(ji,jk) + dz3d(ji,jk+1) ) / 2 74 ! 75 aminuss(ji,jk) = ( 1.0 / volw3d(ji,jk) ) * aminus / dxminus 76 apluss (ji,jk) = ( 1.0 / volw3d(ji,jk) ) * aplus / dxplus 77 END DO 78 79 aminus = ( por(nksed-1) + por(nksed) ) * 0.5 80 dxminus = ( dz3d(ji,nksed-1) + dz3d(ji,nksed) ) / 2. 81 aminuss(ji,nksed) = ( 1.0 / volw3d(ji,nksed) ) * aminus / dxminus 82 83 END DO 84 85 IF( ln_timing ) CALL timing_stop('sed_mat_coef') 86 87 END SUBROUTINE sed_mat_coef 88 89 SUBROUTINE sed_mat_dsr( nksed, nvar, accmask ) 90 !!--------------------------------------------------------------------- 91 !! *** ROUTINE sed_mat_dsr *** 92 !! 93 !! ** Purpose : solves tridiagonal system of linear equations 94 !! 95 !! ** Method : 96 !! 1 - computes left hand side of linear system of equations 97 !! for dissolution reaction 98 !! For mass balance in kbot+sediment : 99 !! dz3d (:,1) = dz(1) = 0.5 cm 100 !! volw3d(:,1) = dzkbot ( see sedini.F90 ) 101 !! dz(2) = 0.3 cm 102 !! dz3d(:,2) = 0.3 + dzdep ( see seddsr.F90 ) 103 !! volw3d(:,2) and vols3d(l,2) are thickened ( see seddsr.F90 ) 104 !! 105 !! 2 - forward/backward substitution. 106 !! 107 !! History : 108 !! ! 04-10 (N. Emprin, M. Gehlen ) original 109 !! ! 06-04 (C. Ethe) Module Re-organization 110 !!---------------------------------------------------------------------- 111 !! * Arguments 112 INTEGER , INTENT(in) :: nvar, nksed ! number of variable 113 INTEGER, DIMENSION(jpoce) :: accmask 114 INTEGER :: ji 115 116 !---Local declarations 117 INTEGER :: jk, jn 118 REAL(wp), DIMENSION(nksed) :: za, zb, zc 119 120 REAL(wp) :: rplus,rminus 121 !---------------------------------------------------------------------- 122 123 IF( ln_timing ) CALL timing_start('sed_mat_dsr') 124 125 ! Computation left hand side of linear system of 126 ! equations for dissolution reaction 127 !--------------------------------------------- 79 128 jn = nvar 80 129 ! first sediment level 81 DO ji = 1, ndim 82 aplus = ( ( volw3d(ji,1) / ( dz3d(ji,1) ) ) + & 83 ( volw3d(ji,2) / ( dz3d(ji,2) ) ) ) / 2. 84 dxplus = ( dz3d(ji,1) + dz3d(ji,2) ) / 2. 85 rplus = ( dtsed_in / ( volw3d(ji,1) ) ) * diff(ji,1,jn) * aplus / dxplus 130 DO ji = 1, jpoce 131 IF (accmask(ji) == 0) THEN 132 rplus = apluss(ji,1) * diff(ji,1,jn) * radssol(1,jn) 133 134 za(1) = 0. 135 zb(1) = rplus 136 zc(1) = -rplus 137 138 DO jk = 2, nksed - 1 139 rminus = aminuss(ji,jk) * diff(ji,jk-1,jn) * radssol(jk,jn) 140 rplus = apluss (ji,jk) * diff(ji,jk,jn) * radssol(jk,jn) 141 ! 142 za(jk) = -rminus 143 zb(jk) = rminus + rplus 144 zc(jk) = -rplus 145 END DO 146 147 rminus = aminuss(ji,nksed) * diff(ji,nksed-1,jn) * radssol(nksed,jn) 148 ! 149 za(nksed) = -rminus 150 zb(nksed) = rminus 151 zc(nksed) = 0. 152 153 ! solves tridiagonal system of linear equations 154 ! ----------------------------------------------- 155 156 pwcpa(ji,1,jn) = pwcpa(ji,1,jn) - ( zc(1) * pwcp(ji,2,jn) + zb(1) * pwcp(ji,1,jn) ) 157 DO jk = 2, nksed - 1 158 pwcpa(ji,jk,jn) = pwcpa(ji,jk,jn) - ( zc(jk) * pwcp(ji,jk+1,jn) + za(jk) * pwcp(ji,jk-1,jn) & 159 & + zb(jk) * pwcp(ji,jk,jn) ) 160 ENDDO 161 pwcpa(ji,nksed,jn) = pwcpa(ji,nksed,jn) - ( za(nksed) * pwcp(ji,nksed-1,jn) & 162 & + zb(nksed) * pwcp(ji,nksed,jn) ) 163 164 ENDIF 165 END DO 166 167 IF( ln_timing ) CALL timing_stop('sed_mat_dsr') 168 169 END SUBROUTINE sed_mat_dsr 170 171 SUBROUTINE sed_mat_dsrjac( nksed, nvar, NEQ, NROWPD, jacvode, accmask ) 172 !!--------------------------------------------------------------------- 173 !! *** ROUTINE sed_mat_dsrjac *** 174 !! 175 !! ** Purpose : solves tridiagonal system of linear equations 176 !! 177 !! ** Method : 178 !! 1 - computes left hand side of linear system of equations 179 !! for dissolution reaction 180 !! For mass balance in kbot+sediment : 181 !! dz3d (:,1) = dz(1) = 0.5 cm 182 !! volw3d(:,1) = dzkbot ( see sedini.F90 ) 183 !! dz(2) = 0.3 cm 184 !! dz3d(:,2) = 0.3 + dzdep ( see seddsr.F90 ) 185 !! volw3d(:,2) and vols3d(l,2) are thickened ( see 186 !seddsr.F90 ) 187 !! 188 !! 2 - forward/backward substitution. 189 !! 190 !! History : 191 !! ! 04-10 (N. Emprin, M. Gehlen ) original 192 !! ! 06-04 (C. Ethe) Module Re-organization 193 !!---------------------------------------------------------------------- 194 !! * Arguments 195 INTEGER , INTENT(in) :: nvar, nksed, NEQ, NROWPD ! number of variable 196 REAL, DIMENSION(jpoce,NROWPD,NEQ), INTENT(inout) :: jacvode 197 INTEGER, DIMENSION(jpoce), INTENT(in) :: accmask 198 199 !---Local declarations 200 INTEGER :: ji,jk, jn, jnn, jni, jnj ,jnij 201 REAL(wp), DIMENSION(nksed) :: za, zb, zc 202 203 REAL(wp) :: rplus,rminus 204 !---------------------------------------------------------------------- 205 206 IF( ln_timing ) CALL timing_start('sed_mat_dsrjac') 207 208 ! Computation left hand side of linear system of 209 ! equations for dissolution reaction 210 !--------------------------------------------- 211 jn = nvar 212 ! first sediment level 213 DO ji = 1, jpoce 214 IF (accmask(ji) == 0 ) THEN 215 rplus = apluss(ji,1) * diff(ji,1,jn) * radssol(1,jn) 216 217 za(1) = 0. 218 zb(1) = rplus 219 zc(1) = -rplus 220 221 DO jk = 2, nksed - 1 222 rminus = aminuss(ji,jk) * diff(ji,jk-1,jn) * radssol(jk,jn) 223 rplus = apluss (ji,jk) * diff(ji,jk,jn) * radssol(jk,jn) 224 ! 225 za(jk) = -rminus 226 zb(jk) = rminus + rplus 227 zc(jk) = -rplus 228 END DO 229 230 rminus = aminuss(ji,nksed) * diff(ji,nksed-1,jn) * radssol(nksed,jn) 231 ! 232 za(nksed) = -rminus 233 zb(nksed) = rminus 234 zc(nksed) = 0. 235 236 ! solves tridiagonal system of linear equations 237 238 jnn = isvode(jn) 239 jnij = jpvode + 1 240 jacvode(ji, jnij, jnn) = jacvode(ji,jnij,jnn) - zb(1) 241 jnj = jpvode + jnn 242 jnij = jnn - jnj + jpvode + 1 243 jacvode(ji, jnij, jnj) = jacvode(ji, jnij, jnj) -zc(1) 244 DO jk = 2, nksed - 1 245 jni = (jk-1) * jpvode + jnn 246 jnj = (jk-2) * jpvode + jnn 247 jnij = jni - jnj + jpvode + 1 248 jacvode(ji, jnij, jnj) = jacvode(ji, jnij, jnj) - za(jk) 249 jnj = (jk-1) * jpvode + jnn 250 jnij = jni - jnj + jpvode + 1 251 jacvode(ji, jnij, jnj) = jacvode(ji, jnij, jnj) - zb(jk) 252 jnj = (jk) * jpvode + jnn 253 jnij = jni - jnj + jpvode + 1 254 jacvode(ji, jnij, jnj) = jacvode(ji, jnij, jnj) - zc(jk) 255 END DO 256 jni = (nksed-1) * jpvode + jnn 257 jnj = (nksed-2) * jpvode + jnn 258 jnij = jni - jnj + jpvode + 1 259 jacvode(ji, jnij, jnj) = jacvode(ji, jnij, jnj) - za(nksed) 260 jnij = jpvode + 1 261 jacvode(ji, jnij, jni) = jacvode(ji, jnij, jni) - zb(nksed) 262 ENDIF 263 END DO 264 265 IF( ln_timing ) CALL timing_stop('sed_mat_dsrjac') 266 267 END SUBROUTINE sed_mat_dsrjac 268 269 SUBROUTINE sed_mat_btbi( nksed, nvar, psol, preac, dtsed_in ) 270 !!--------------------------------------------------------------------- 271 !! *** ROUTINE sed_mat_btb *** 272 !! 273 !! ** Purpose : solves tridiagonal system of linear equations 274 !! 275 !! ** Method : 276 !! 1 - computes left hand side of linear system of equations 277 !! for dissolution reaction 278 !! 279 !! 280 !! 2 - forward/backward substitution. 281 !! 282 !! History : 283 !! ! 04-10 (N. Emprin, M. Gehlen ) original 284 !! ! 06-04 (C. Ethe) Module Re-organization 285 !!---------------------------------------------------------------------- 286 !! * Arguments 287 INTEGER , INTENT(in) :: nksed, nvar ! number of sediment levels 288 289 REAL(wp), DIMENSION(jpoce,nksed,nvar), INTENT(inout) :: & 290 psol, preac 291 292 REAL(wp), INTENT(in) :: dtsed_in 293 294 !---Local declarations 295 INTEGER :: & 296 ji, jk, jn 297 298 REAL(wp) :: & 299 aplus,aminus , & 300 rplus,rminus , & 301 dxplus,dxminus 302 303 REAL(wp), DIMENSION(nksed) :: za, zb, zc 304 REAL(wp), DIMENSION(nksed) :: zr, zgamm 305 REAL(wp) :: zbet 306 307 !---------------------------------------------------------------------- 308 309 ! Computation left hand side of linear system of 310 ! equations for dissolution reaction 311 !--------------------------------------------- 312 IF( ln_timing ) CALL timing_start('sed_mat_btbi') 313 314 ! first sediment level 315 DO ji = 1, jpoce 316 aplus = ( por1(2) + por1(3) ) / 2.0 317 dxplus = ( dz(2) + dz(3) ) / 2. 318 rplus = ( dtsed_in / vols(2) ) * db(ji,2) * aplus / dxplus 319 za(2) = 0. 320 zb(2) = 1. + rplus 321 zc(2) = -rplus 322 323 DO jk = 3, nksed - 1 324 aminus = ( por1(jk-1) + por1(jk) ) * 0.5 325 aminus = ( ( vols(jk-1) / dz(jk-1) ) + ( vols(jk) / dz(jk) ) ) / 2. 326 dxminus = ( dz(jk-1) + dz(jk) ) / 2. 327 rminus = ( dtsed_in / vols(jk) ) * db(ji,jk-1) * aminus / dxminus 328 ! 329 aplus = ( por1(jk) + por1(jk+1) ) * 0.5 330 dxplus = ( dz(jk) + dz(jk+1) ) / 2. 331 rplus = ( dtsed_in / vols(jk) ) * db(ji,jk) * aplus / dxplus 332 ! 333 za(jk) = -rminus 334 zb(jk) = 1. + rminus + rplus 335 zc(jk) = -rplus 336 337 ENDDO 338 339 aminus = ( por1(nksed-1) + por1(nksed) ) * 0.5 340 dxminus = ( dz(nksed-1) + dz(nksed) ) / 2. 341 rminus = ( dtsed_in / vols(nksed) ) * db(ji,nksed-1) * aminus / dxminus 342 ! 343 za(nksed) = -rminus 344 zb(nksed) = 1. + rminus 345 zc(nksed) = 0. 346 347 ! solves tridiagonal system of linear equations 348 ! ----------------------------------------------- 349 DO jn = 1, nvar 350 zr(:) = psol(ji,:,jn) 351 zbet = zb(2) - preac(ji,2,jn) * dtsed_in 352 psol(ji,2,jn) = zr(2) / zbet 353 ! 354 DO jk = 3, nksed 355 zgamm(jk) = zc(jk-1) / zbet 356 zbet = zb(jk) - preac(ji,jk,jn) * dtsed_in - za(jk) * zgamm(jk) 357 psol(ji,jk,jn) = ( zr(jk) - za(jk) * psol(ji,jk-1,jn) ) / zbet 358 ENDDO 359 ! 360 DO jk = nksed - 1, 2, -1 361 psol(ji,jk,jn) = psol(ji,jk,jn) - zgamm(jk+1) * psol(ji,jk+1,jn) 362 ENDDO 363 END DO 364 END DO 365 ! 366 IF( ln_timing ) CALL timing_stop('sed_mat_btbi') 367 368 END SUBROUTINE sed_mat_btbi 369 370 371 SUBROUTINE sed_mat_btb( nksed, nvar, accmask ) 372 !!--------------------------------------------------------------------- 373 !! *** ROUTINE sed_mat_btb *** 374 !! 375 !! ** Purpose : solves tridiagonal system of linear equations 376 !! 377 !! ** Method : 378 !! 1 - computes left hand side of linear system of equations 379 !! for dissolution reaction 380 !! 381 !! 2 - forward/backward substitution. 382 !! 383 !! History : 384 !! ! 04-10 (N. Emprin, M. Gehlen ) original 385 !! ! 06-04 (C. Ethe) Module Re-organization 386 !!---------------------------------------------------------------------- 387 !! * Arguments 388 INTEGER , INTENT(in) :: & 389 nvar, nksed ! number of sediment levels 390 INTEGER, DIMENSION(jpoce) :: accmask 391 392 !---Local declarations 393 INTEGER :: ji, jk, jn 394 395 REAL(wp) :: & 396 aplus,aminus , & 397 rplus,rminus , & 398 dxplus,dxminus 399 400 REAL(wp), DIMENSION(nksed) :: za, zb, zc 401 402 !---------------------------------------------------------------------- 403 404 ! Computation left hand side of linear system of 405 ! equations for dissolution reaction 406 !--------------------------------------------- 407 IF( ln_timing ) CALL timing_start('sed_mat_btb') 408 409 ! first sediment level 410 jn = nvar 411 DO ji = 1, jpoce 412 IF (accmask(ji) == 0) THEN 413 aplus = ( por1(2) + por1(3) ) / 2.0 414 dxplus = ( dz(2) + dz(3) ) / 2. 415 rplus = ( 1.0 / vols(2) ) * db(ji,2) * aplus / dxplus * rads1sol(2,jn) 416 417 za(2) = 0. 418 zb(2) = rplus 419 zc(2) = -rplus 420 421 DO jk = 3, nksed - 1 422 aminus = ( por1(jk-1) + por1(jk) ) * 0.5 423 aminus = ( ( vols(jk-1) / dz(jk-1) ) + ( vols(jk) / dz(jk) ) ) / 2. 424 dxminus = ( dz(jk-1) + dz(jk) ) / 2. 425 rminus = ( 1.0 / vols(jk) ) * db(ji,jk-1) * aminus / dxminus * rads1sol(jk,jn) 426 ! 427 aplus = ( por1(jk) + por1(jk+1) ) * 0.5 428 dxplus = ( dz(jk) + dz(jk+1) ) / 2. 429 rplus = ( 1.0 / vols(jk) ) * db(ji,jk) * aplus / dxplus * rads1sol(jk,jn) 430 ! 431 za(jk) = -rminus 432 zb(jk) = rminus + rplus 433 zc(jk) = -rplus 434 435 ENDDO 436 437 aminus = ( por1(nksed-1) + por1(nksed) ) * 0.5 438 dxminus = ( dz(nksed-1) + dz(nksed) ) / 2. 439 rminus = ( 1.0 / vols(nksed) ) * db(ji,nksed-1) * aminus / dxminus * rads1sol(nksed,jn) 440 ! 441 za(nksed) = -rminus 442 zb(nksed) = rminus 443 zc(nksed) = 0. 444 445 ! solves tridiagonal system of linear equations 446 ! ----------------------------------------------- 447 pwcpa(ji,2,jn) = pwcpa(ji,2,jn) - ( zc(2) * pwcp(ji,3,jn) + zb(2) * pwcp(ji,2,jn) ) 448 DO jk = 3, nksed-1 449 pwcpa(ji,jk,jn) = pwcpa(ji,jk,jn) - ( zc(jk) * pwcp(ji,jk+1,jn) + za(jk) * pwcp(ji,jk-1,jn) & 450 & + zb(jk) * pwcp(ji,jk,jn) ) 451 ENDDO 452 pwcpa(ji,nksed,jn) = pwcpa(ji,nksed,jn) - ( za(nksed) * pwcp(ji,nksed-1,jn) & 453 & + zb(nksed) * pwcp(ji,nksed,jn) ) 454 ! 455 ENDIF 456 END DO 457 ! 458 IF( ln_timing ) CALL timing_stop('sed_mat_btb') 459 460 END SUBROUTINE sed_mat_btb 461 462 SUBROUTINE sed_mat_btbjac( nksed, nvar, NEQ, NROWPD, jacvode, accmask ) 463 !!--------------------------------------------------------------------- 464 !! *** ROUTINE sed_mat_btb *** 465 !! 466 !! ** Purpose : solves tridiagonal system of linear equations 467 !! 468 !! ** Method : 469 !! 1 - computes left hand side of linear system of equations 470 !! for dissolution reaction 471 !! 472 !! 2 - forward/backward substitution. 473 !! 474 !! History : 475 !! ! 04-10 (N. Emprin, M. Gehlen ) original 476 !! ! 06-04 (C. Ethe) Module Re-organization 477 !!---------------------------------------------------------------------- 478 !! * Arguments 479 INTEGER , INTENT(in) :: nvar, nksed, NEQ, NROWPD ! number of variable 480 REAL, DIMENSION(jpoce,NROWPD,NEQ), INTENT(inout) :: jacvode 481 INTEGER, DIMENSION(jpoce), INTENT(in) :: accmask 482 483 !---Local declarations 484 INTEGER :: ji, jk, jn, jnn, jni, jnj ,jnij 485 486 REAL(wp) :: & 487 aplus,aminus , & 488 rplus,rminus , & 489 dxplus,dxminus 490 491 REAL(wp), DIMENSION(nksed) :: za, zb, zc 492 493 !---------------------------------------------------------------------- 494 495 ! Computation left hand side of linear system of 496 ! equations for dissolution reaction 497 !--------------------------------------------- 498 IF( ln_timing ) CALL timing_start('sed_mat_btbjac') 499 500 ! first sediment level 501 jn = nvar 502 DO ji = 1, jpoce 503 IF (accmask(ji) == 0) THEN 504 aplus = ( por1(2) + por1(3) ) / 2.0 505 dxplus = ( dz(2) + dz(3) ) / 2. 506 rplus = ( 1.0 / vols(2) ) * db(ji,2) * aplus / dxplus * rads1sol(2,jn) 507 508 za(2) = 0. 509 zb(2) = rplus 510 zc(2) = -rplus 511 512 DO jk = 3, nksed - 1 513 aminus = ( por1(jk-1) + por1(jk) ) * 0.5 514 aminus = ( ( vols(jk-1) / dz(jk-1) ) + ( vols(jk) / dz(jk) ) ) / 2. 515 dxminus = ( dz(jk-1) + dz(jk) ) / 2. 516 rminus = ( 1.0 / vols(jk) ) * db(ji,jk-1) * aminus / dxminus * rads1sol(jk,jn) 517 ! 518 aplus = ( por1(jk) + por1(jk+1) ) * 0.5 519 dxplus = ( dz(jk) + dz(jk+1) ) / 2. 520 rplus = ( 1.0 / vols(jk) ) * db(ji,jk) * aplus / dxplus * rads1sol(jk,jn) 521 ! 522 za(jk) = -rminus 523 zb(jk) = rminus + rplus 524 zc(jk) = -rplus 525 526 ENDDO 527 528 aminus = ( por1(nksed-1) + por1(nksed) ) * 0.5 529 dxminus = ( dz(nksed-1) + dz(nksed) ) / 2. 530 rminus = ( 1.0 / vols(nksed) ) * db(ji,nksed-1) * aminus / dxminus * rads1sol(nksed,jn) 531 ! 532 za(nksed) = -rminus 533 zb(nksed) = rminus 534 zc(nksed) = 0. 535 536 ! solves tridiagonal system of linear equations 537 ! ----------------------------------------------- 538 jnn = isvode(jn) 539 jni = jpvode + jnn 540 jnij = jpvode + 1 541 jacvode(ji, jnij, jni) = jacvode(ji,jnij,jni) - zb(2) 542 jnj = 2 * jpvode + jnn 543 jnij = jni - jnj + jpvode + 1 544 jacvode(ji, jnij, jnj) = jacvode(ji, jnij, jnj) -zc(2) 545 DO jk = 3, nksed-1 546 jni = (jk-1) * jpvode + jnn 547 jnj = (jk-2) * jpvode + jnn 548 jnij = jni - jnj + jpvode + 1 549 jacvode(ji, jnij, jnj) = jacvode(ji, jnij, jnj) - za(jk) 550 jnj = (jk-1) * jpvode + jnn 551 jnij = jni - jnj + jpvode + 1 552 jacvode(ji, jnij, jnj) = jacvode(ji, jnij, jnj) - zb(jk) 553 jnj = (jk) * jpvode + jnn 554 jnij = jni - jnj + jpvode + 1 555 jacvode(ji, jnij, jnj) = jacvode(ji, jnij, jnj) - zc(jk) 556 ENDDO 557 jni = (nksed-1) * jpvode + jnn 558 jnj = (nksed-2) * jpvode + jnn 559 jnij = jni - jnj + jpvode + 1 560 jacvode(ji, jnij, jnj) = jacvode(ji, jnij, jnj) - za(nksed) 561 jnij = jpvode + 1 562 jacvode(ji, jnij, jni) = jacvode(ji, jnij, jni) - zb(nksed) 563 ! 564 ENDIF 565 END DO 566 ! 567 IF( ln_timing ) CALL timing_stop('sed_mat_btbjac') 568 569 END SUBROUTINE sed_mat_btbjac 570 571 572 SUBROUTINE sed_mat_dsri( nksed, nvar, preac, psms, dtsed_in, psol ) 573 !!--------------------------------------------------------------------- 574 !! *** ROUTINE sed_mat_dsr *** 575 !! 576 !! ** Purpose : solves tridiagonal system of linear equations 577 !! 578 !! ** Method : 579 !! 1 - computes left hand side of linear system of equations 580 !! for dissolution reaction 581 !! For mass balance in kbot+sediment : 582 !! dz3d (:,1) = dz(1) = 0.5 cm 583 !! volw3d(:,1) = dzkbot ( see sedini.F90 ) 584 !! dz(2) = 0.3 cm 585 !! dz3d(:,2) = 0.3 + dzdep ( see seddsr.F90 ) 586 !! volw3d(:,2) and vols3d(l,2) are thickened ( see 587 !seddsr.F90 ) 588 !! 589 !! 2 - forward/backward substitution. 590 !! 591 !! History : 592 !! ! 04-10 (N. Emprin, M. Gehlen ) original 593 !! ! 06-04 (C. Ethe) Module Re-organization 594 !!---------------------------------------------------------------------- 595 !! * Arguments 596 INTEGER , INTENT(in) :: nksed, nvar ! number of variable 597 598 REAL(wp), DIMENSION(jpoce,nksed), INTENT(in ) :: preac ! reaction rates 599 REAL(wp), DIMENSION(jpoce,nksed), INTENT(in ) :: psms ! reaction rates 600 REAL(wp), DIMENSION(jpoce,nksed), INTENT(inout) :: psol ! reaction rates 601 REAL(wp), INTENT(in) :: dtsed_in 602 603 !---Local declarations 604 INTEGER :: ji, jk, jn 605 REAL(wp), DIMENSION(jpoce,nksed) :: za, zb, zc, zr 606 REAL(wp), DIMENSION(jpoce) :: zbet 607 REAL(wp), DIMENSION(jpoce,nksed) :: zgamm 608 609 REAL(wp) :: rplus,rminus 610 !---------------------------------------------------------------------- 611 612 IF( ln_timing ) CALL timing_start('sed_mat_dsri') 613 614 ! Computation left hand side of linear system of 615 ! equations for dissolution reaction 616 !--------------------------------------------- 617 jn = nvar 618 ! first sediment level 619 DO ji = 1, jpoce 620 rplus = dtsed_in * apluss(ji,1) * diff(ji,1,jn) * radssol(1,jn) 86 621 87 622 za(ji,1) = 0. … … 89 624 zc(ji,1) = -rplus 90 625 ENDDO 91 92 DO jk = 2, nlev - 1 93 DO ji = 1, ndim 94 aminus = ( ( volw3d(ji,jk-1) / ( dz3d(ji,jk-1) ) ) + & 95 & ( volw3d(ji,jk ) / ( dz3d(ji,jk ) ) ) ) / 2. 96 dxminus = ( dz3d(ji,jk-1) + dz3d(ji,jk) ) / 2. 97 98 aplus = ( ( volw3d(ji,jk ) / ( dz3d(ji,jk ) ) ) + & 99 & ( volw3d(ji,jk+1) / ( dz3d(ji,jk+1) ) ) ) / 2. 100 dxplus = ( dz3d(ji,jk) + dz3d(ji,jk+1) ) / 2 101 ! 102 rminus = ( dtsed_in / volw3d(ji,jk) ) * diff(ji,jk-1,jn) * aminus / dxminus 103 rplus = ( dtsed_in / volw3d(ji,jk) ) * diff(ji,jk,jn) * aplus / dxplus 626 627 DO jk = 2, nksed - 1 628 DO ji = 1, jpoce 629 rminus = dtsed_in * aminuss(ji,jk) * diff(ji,jk-1,jn) * radssol(jk,jn) 630 rplus = dtsed_in * apluss (ji,jk) * diff(ji,jk,jn) * radssol(jk,jn) 104 631 ! 105 632 za(ji,jk) = -rminus 106 zb(ji,jk) = 1. + rminus + rplus 633 zb(ji,jk) = 1. + rminus + rplus 107 634 zc(ji,jk) = -rplus 108 635 END DO 109 636 END DO 110 637 111 DO ji = 1, ndim 112 aminus = ( ( volw3d(ji,nlev-1) / dz3d(ji,nlev-1) ) + & 113 & ( volw3d(ji,nlev) / dz3d(ji,nlev) ) ) / 2. 114 dxminus = ( dz3d(ji,nlev-1) + dz3d(ji,nlev) ) / 2. 115 rminus = ( dtsed_in / volw3d(ji,nlev) ) * diff(ji,nlev-1,jn) * aminus / dxminus 638 DO ji = 1, jpoce 639 rminus = dtsed_in * aminuss(ji,nksed) * diff(ji,nksed-1,jn) * radssol(nksed,jn) 116 640 ! 117 za(ji,n lev) = -rminus118 zb(ji,n lev) = 1. + rminus119 zc(ji,n lev) = 0.641 za(ji,nksed) = -rminus 642 zb(ji,nksed) = 1. + rminus 643 zc(ji,nksed) = 0. 120 644 END DO 121 645 … … 124 648 ! ----------------------------------------------- 125 649 126 zr (:,:) = psol(:,:) + psms(:,:) 127 zb (:,:) = zb(:,:) + preac(:,:)650 zr (:,:) = psol(:,:) + psms(:,:) * dtsed_in 651 zb (:,:) = zb(:,:) - preac(:,:) * dtsed_in 128 652 zbet(: ) = zb(:,1) 129 653 psol(:,1) = zr(:,1) / zbet(:) 130 654 131 655 ! 132 DO jk = 2, n lev133 DO ji = 1, ndim656 DO jk = 2, nksed 657 DO ji = 1, jpoce 134 658 zgamm(ji,jk) = zc(ji,jk-1) / zbet(ji) 135 659 zbet(ji) = zb(ji,jk) - za(ji,jk) * zgamm(ji,jk) … … 138 662 ENDDO 139 663 ! 140 DO jk = n lev- 1, 1, -1141 DO ji = 1, ndim664 DO jk = nksed - 1, 1, -1 665 DO ji = 1, jpoce 142 666 psol(ji,jk) = psol(ji,jk) - zgamm(ji,jk+1) * psol(ji,jk+1) 143 667 END DO 144 668 ENDDO 145 669 146 IF( ln_timing ) CALL timing_stop('sed_mat_dsr') 147 148 149 END SUBROUTINE sed_mat_dsr 150 151 SUBROUTINE sed_mat_btb( nvar, ndim, nlev, psol, dtsed_in ) 152 !!--------------------------------------------------------------------- 153 !! *** ROUTINE sed_mat_btb *** 154 !! 155 !! ** Purpose : solves tridiagonal system of linear equations 156 !! 157 !! ** Method : 158 !! 1 - computes left hand side of linear system of equations 159 !! for dissolution reaction 160 !! 161 !! 2 - forward/backward substitution. 162 !! 163 !! History : 164 !! ! 04-10 (N. Emprin, M. Gehlen ) original 165 !! ! 06-04 (C. Ethe) Module Re-organization 166 !!---------------------------------------------------------------------- 167 !! * Arguments 168 INTEGER , INTENT(in) :: & 169 nvar , & ! number of variables 170 ndim , & ! number of points 171 nlev ! number of sediment levels 172 173 REAL(wp), DIMENSION(ndim,nlev,nvar), INTENT(inout) :: & 174 psol ! solution 175 176 REAL(wp), INTENT(in) :: dtsed_in 177 178 !---Local declarations 179 INTEGER :: & 180 ji, jk, jn 181 182 REAL(wp) :: & 183 aplus,aminus , & 184 rplus,rminus , & 185 dxplus,dxminus 186 187 REAL(wp), DIMENSION(nlev) :: za, zb, zc 188 REAL(wp), DIMENSION(ndim,nlev) :: zr 189 REAL(wp), DIMENSION(nmax) :: zgamm 190 REAL(wp) :: zbet 191 192 193 !---------------------------------------------------------------------- 194 195 ! Computation left hand side of linear system of 196 ! equations for dissolution reaction 197 !--------------------------------------------- 198 199 200 IF( ln_timing ) CALL timing_start('sed_mat_btb') 201 202 ! first sediment level 203 DO ji = 1, ndim 204 aplus = ( ( vols(2) / dz(2) ) + ( vols(3) / dz(3) ) ) / 2. 205 dxplus = ( dz(2) + dz(3) ) / 2. 206 rplus = ( dtsed_in / vols(2) ) * db(ji,2) * aplus / dxplus 207 208 za(1) = 0. 209 zb(1) = 1. + rplus 210 zc(1) = -rplus 211 212 213 DO jk = 2, nlev - 1 214 aminus = ( ( vols(jk) / dz(jk) ) + ( vols(jk+1) / dz(jk+1) ) ) / 2. 215 dxminus = ( dz(jk) + dz(jk+1) ) / 2. 216 rminus = ( dtsed_in / vols(jk+1) ) * db(ji,jk) * aminus / dxminus 217 ! 218 aplus = ( ( vols(jk+1) / dz(jk+1 ) ) + ( vols(jk+2) / dz(jk+2) ) ) / 2. 219 dxplus = ( dz(jk+1) + dz(jk+2) ) / 2. 220 rplus = ( dtsed_in / vols(jk+1) ) * db(ji,jk+1) * aplus / dxplus 221 ! 222 za(jk) = -rminus 223 zb(jk) = 1. + rminus + rplus 224 zc(jk) = -rplus 225 ENDDO 226 227 aminus = ( ( vols(nlev) / dz(nlev) ) + ( vols(nlev+1) / dz(nlev+1) ) ) / 2. 228 dxminus = ( dz(nlev) + dz(nlev+1) ) / 2. 229 rminus = ( dtsed_in / vols(nlev+1) ) * db(ji,nlev) * aminus / dxminus 230 ! 231 za(nlev) = -rminus 232 zb(nlev) = 1. + rminus 233 zc(nlev) = 0. 234 235 236 ! solves tridiagonal system of linear equations 237 ! ----------------------------------------------- 238 DO jn = 1, nvar 239 240 DO jk = 1, nlev 241 zr (ji,jk) = psol(ji,jk,jn) 242 END DO 243 zbet = zb(1) 244 psol(ji,1,jn) = zr(ji,1) / zbet 245 ! 246 DO jk = 2, nlev 247 zgamm(jk) = zc(jk-1) / zbet 248 zbet = zb(jk) - za(jk) * zgamm(jk) 249 psol(ji,jk,jn) = ( zr(ji,jk) - za(jk) * psol(ji,jk-1,jn) ) / zbet 250 ENDDO 251 ! 252 DO jk = nlev - 1, 1, -1 253 psol(ji,jk,jn) = psol(ji,jk,jn) - zgamm(jk+1) * psol(ji,jk+1,jn) 254 ENDDO 255 256 ENDDO 257 258 END DO 259 ! 260 IF( ln_timing ) CALL timing_stop('sed_mat_btb') 261 262 263 END SUBROUTINE sed_mat_btb 670 IF( ln_timing ) CALL timing_stop('sed_mat_dsri') 671 672 673 END SUBROUTINE sed_mat_dsri 674 264 675 265 676 END MODULE sedmat -
NEMO/trunk/src/TOP/PISCES/SED/sedorg.F90
r10225 r15450 1 1 MODULE sedorg 2 2 !!====================================================================== 3 !! *** MODULE sed dsr***4 !! Sediment : dissolution and reaction in pore water related5 !! related to organic matter3 !! *** MODULE sedinorg *** 4 !! Sediment : dissolution and reaction in pore water of 5 !! inorganic species 6 6 !!===================================================================== 7 7 !! * Modules used … … 9 9 USE sed_oce 10 10 USE sedini 11 USE seddiff 12 USE seddsr 11 USE sedmat 13 12 USE lib_mpp ! distribued memory computing library 14 13 USE lib_fortran … … 19 18 PUBLIC sed_org 20 19 21 !! * Module variables22 23 REAL(wp) :: zadsnh424 25 20 !! $Id: seddsr.F90 5215 2015-04-15 16:11:56Z nicolasmartin $ 26 21 CONTAINS 27 22 28 SUBROUTINE sed_org( kt ) 23 SUBROUTINE sed_org( kt ) 29 24 !!---------------------------------------------------------------------- 30 25 !! *** ROUTINE sed_org *** 31 26 !! 32 !! ** Purpose : computes pore water di ffusion and reaction27 !! ** Purpose : computes pore water dissolution and reaction 33 28 !! 34 !! ** Methode : Computation of the redox reactions in sediment. 35 !! The main redox reactions are solved in sed_dsr whereas 36 !! the secondary reactions are solved in sed_dsr_redoxb. 37 !! A strand spliting approach is being used here (see 38 !! sed_dsr_redoxb for more information). 39 !! Diffusive fluxes are computed in sed_diff 29 !! ** Methode : implicit simultaneous computation of undersaturation 30 !! resulting from diffusive pore water transport and chemical 31 !! pore water reactions. Solid material is consumed according 32 !! to redissolution and remineralisation 33 !! 34 !! ** Remarks : 35 !! - undersaturation : deviation from saturation concentration 36 !! - reaction rate : sink of undersaturation from dissolution 37 !! of solid material 40 38 !! 41 39 !! History : … … 43 41 !! ! 04-10 (N. Emprin, M. Gehlen ) f90 44 42 !! ! 06-04 (C. Ethe) Re-organization 45 !! ! 19-08 (O. Aumont) Debugging and improvement of the model. 46 !! The original method is replaced by a 47 !! Strand splitting method which deals 48 !! well with stiff reactions. 43 !! ! 19-08 (O. Aumont) Debugging and improvement of the model 49 44 !!---------------------------------------------------------------------- 50 45 !! Arguments 51 INTEGER, INTENT(in) :: kt46 INTEGER, INTENT(in) :: kt ! time step 52 47 ! --- local variables 53 INTEGER :: ji, jk, js, jw, jnt ! dummy looop indices 54 REAL(wp) :: zadsnh4 48 REAL(wp), DIMENSION(jpoce, jpksed) :: psms, preac 55 49 !! 56 50 !!---------------------------------------------------------------------- 57 51 58 52 IF( ln_timing ) CALL timing_start('sed_org') 53 54 IF( kt == nitsed000 ) THEN 55 IF (lwp) WRITE(numsed,*) ' sed_org : solute species which do not experience redox reactions ' 56 IF (lwp) WRITE(numsed,*) ' ' 57 ENDIF 59 58 ! 60 IF( kt == nitsed000 ) THEN 61 IF (lwp) THEN 62 WRITE(numsed,*) ' sed_org : Organic degradation related reactions and diffusion' 63 WRITE(numsed,*) ' ' 64 ENDIF 65 ! ! 66 dens_mol_wgt(1:jpsol) = denssol / mol_wgt(1:jpsol) 67 ! 68 ENDIF 59 ! DIC in pore water 60 preac(:,:) = 0.0_wp 61 psms (:,:) = rearatpom(:,:) 62 CALL sed_mat_dsri( jpksed, jwdic, preac, psms, dtsed, pwcp(:,:,jwdic) ) 63 64 ! Silicate in pore water 65 psms (:,:) = 0.0_wp 66 CALL sed_mat_dsri( jpksed, jwsil, preac, psms, dtsed, pwcp(:,:,jwsil) ) 69 67 70 dtsed2 = dtsed / REAL( nrseddt, wp ) 71 72 ! 1. Change of geometry 73 ! Increase of dz3d(2) thickness : dz3d(2) = dz3d(2)+dzdep 74 ! Warning : no change for dz(2) 75 !--------------------------------------------------------- 76 dz3d(1:jpoce,2) = dz3d(1:jpoce,2) + dzdep(1:jpoce) 77 78 ! New values for volw3d(:,2) and vols3d(:,2) 79 ! Warning : no change neither for volw(2) nor vols(2) 80 !------------------------------------------------------ 81 volw3d(1:jpoce,2) = dz3d(1:jpoce,2) * por(2) 82 vols3d(1:jpoce,2) = dz3d(1:jpoce,2) * por1(2) 83 84 ! 2. Change of previous solid fractions (due to volum changes) for k=2 85 !--------------------------------------------------------------------- 86 87 DO js = 1, jpsol 88 DO ji = 1, jpoce 89 solcp(ji,2,js) = solcp(ji,2,js) * dz(2) / dz3d(ji,2) 90 ENDDO 91 END DO 92 93 ! 3. New solid fractions (including solid rain fractions) for k=2 94 !------------------------------------------------------------------ 95 DO js = 1, jpsol 96 DO ji = 1, jpoce 97 IF (raintg(ji) .ne. 0) THEN 98 solcp(ji,2,js) = solcp(ji,2,js) + & 99 & ( rainrg(ji,js) / raintg(ji) ) * ( dzdep(ji) / dz3d(ji,2) ) 100 ! rainrm are temporary cancel 101 rainrm(ji,js) = 0. 102 ENDIF 103 END DO 104 ENDDO 105 106 ! 4. Adjustment of bottom water concen.(pwcp(1)): 107 ! We impose that pwcp(2) is constant. Including dzdep in dz3d(:,2) we assume 108 ! that dzdep has got a porosity of por(2). So pore water volum of jk=2 increase. 109 ! To keep pwcp(2) cste we must compensate this "increase" by a slight adjusment 110 ! of bottom water concentration. 111 ! This adjustment is compensate at the end of routine 112 !------------------------------------------------------------- 113 DO jw = 1, jpwat 114 DO ji = 1, jpoce 115 pwcp(ji,1,jw) = pwcp(ji,1,jw) - & 116 & pwcp(ji,2,jw) * dzdep(ji) * por(2) / ( dzkbot(ji) + rtrn ) 117 END DO 118 ENDDO 119 120 zadsnh4 = 1.0 / ( 1.0 + adsnh4 ) 121 122 ! -------------------------------------------------- 123 ! Computation of the diffusivities 124 ! -------------------------------------------------- 125 126 DO js = 1, jpwat 127 DO jk = 1, jpksed 128 DO ji = 1, jpoce 129 diff(ji,jk,js) = ( diff1(js) + diff2(js) * temp(ji) ) / ( 1.0 - 2.0 * log( por(jk) ) ) 130 END DO 131 END DO 132 END DO 133 134 ! Impact of bioirrigation and adsorption on diffusion 135 ! --------------------------------------------------- 136 137 diff(:,:,jwnh4) = diff(:,:,jwnh4) * ( 1.0 + irrig(:,:) ) * zadsnh4 138 diff(:,:,jwsil) = diff(:,:,jwsil) * ( 1.0 + irrig(:,:) ) 139 diff(:,:,jwoxy) = diff(:,:,jwoxy) * ( 1.0 + irrig(:,:) ) 140 diff(:,:,jwdic) = diff(:,:,jwdic) * ( 1.0 + irrig(:,:) ) 141 diff(:,:,jwno3) = diff(:,:,jwno3) * ( 1.0 + irrig(:,:) ) 142 diff(:,:,jwpo4) = diff(:,:,jwpo4) * ( 1.0 + irrig(:,:) ) 143 diff(:,:,jwalk) = diff(:,:,jwalk) * ( 1.0 + irrig(:,:) ) 144 diff(:,:,jwh2s) = diff(:,:,jwh2s) * ( 1.0 + irrig(:,:) ) 145 diff(:,:,jwso4) = diff(:,:,jwso4) * ( 1.0 + irrig(:,:) ) 146 diff(:,:,jwfe2) = diff(:,:,jwfe2) * ( 1.0 + 0.2 * irrig(:,:) ) 147 148 DO jnt = 1, nrseddt 149 CALL sed_diff( kt, jnt ) ! 1st pass in diffusion to get values at t+1/2 150 CALL sed_dsr ( kt, jnt ) ! Dissolution reaction 151 CALL sed_diff( kt, jnt ) ! 2nd pass in diffusion to get values at t+1 152 END DO 153 68 ! Iron ligands in pore water 69 psms (:,:) = ratligc * rearatpom(:,:) 70 preac(:,:) = -reac_ligc 71 CALL sed_mat_dsri( jpksed, jwlgw, preac, psms, dtsed, pwcp(:,:,jwlgw) ) 154 72 155 73 IF( ln_timing ) CALL timing_stop('sed_org') -
NEMO/trunk/src/TOP/PISCES/SED/sedrst.F90
r14239 r15450 42 42 CHARACTER(LEN=50) :: clname ! trc output restart file name 43 43 CHARACTER(LEN=256) :: clpath ! full path to ocean output restart file 44 CHARACTER(LEN= 52) :: clpname ! trc output restart file name including AGRIF44 CHARACTER(LEN=3) :: cdcomp 45 45 !!---------------------------------------------------------------------- 46 46 ! … … 65 65 66 66 IF( .NOT. ln_rst_list .AND. nn_stock == -1 ) RETURN ! we will never do any restart 67 68 67 ! to get better performances with NetCDF format: 69 ! we open and define the tracer restart file one tracer time step before writing the data (-> at nitrst - 1)70 ! except if we write tracer restart files every tracer time step or if a tracer restart file was writen at nitend - 171 IF( kt == nitrst - 2*nn_dtsed .OR. nn_stock == nn_dtsed .OR. ( kt == nitend - nn_dtsed.AND. .NOT. lrst_sed ) ) THEN68 ! we open and define the tracer restart file one tracer time step before writing the data (-> at nitrst - 2*nn_dttrc + 1) 69 ! except if we write tracer restart files every tracer time step or if a tracer restart file was writen at nitend - 2*nn_dttrc + 1 70 IF( kt == nitrst - 1 .OR. nn_stock == 1 .OR. ( kt == nitend - 1 .AND. .NOT. lrst_sed ) ) THEN 72 71 ! beware of the format used to write kt (default is i8.8, that should be large enough) 73 72 IF( nitrst > 1.0e9 ) THEN ; WRITE(clkt,* ) nitrst … … 81 80 IF(lwp) WRITE(numsed,*) & 82 81 ' open sed restart.output NetCDF file: ',TRIM(clpath)//clname 83 IF(.NOT.lwxios) THEN 84 CALL iom_open( TRIM(clpath)//TRIM(clname), numrsw, ldwrt = .TRUE., kdlev = jpksed, cdcomp = 'SED' ) 85 ELSE 86 #if defined key_xios 87 cw_sedrst_cxt = "rstws_"//TRIM(ADJUSTL(clkt)) 88 IF( TRIM(Agrif_CFixed()) == '0' ) THEN 89 clpname = clname 90 ELSE 91 clpname = TRIM(Agrif_CFixed())//"_"//clname 92 ENDIF 93 numrsw = iom_xios_setid(TRIM(clpath)//TRIM(clpname)) 94 CALL iom_init( cw_sedrst_cxt, kdid = numrsw, ld_closedef = .FALSE. ) 95 #else 96 CALL ctl_stop( 'Can not use XIOS in trc_rst_opn' ) 97 #endif 98 ENDIF 99 82 cdcomp ='SED' 83 CALL iom_open( TRIM(clpath)//TRIM(clname), numrsw, ldwrt = .TRUE., kdlev = jpksed, cdcomp = cdcomp ) 100 84 lrst_sed = .TRUE. 101 85 ENDIF … … 204 188 & zdta2(1:jpi,1:jpj,1:jpksed), iarroce(1:jpoce) ) 205 189 206 cltra = "sedligand"207 IF( iom_varid( numrsr, TRIM(cltra) , ldstop = .FALSE. ) > 0 ) THEN208 CALL iom_get( numrsr, jpdom_auto, TRIM(cltra), zdta2(:,:,:) )209 ELSE210 zdta2(:,:,:) = 0.0211 ENDIF212 213 CALL pack_arr( jpoce, sedligand(1:jpoce,1:jpksed), &214 & zdta2(1:jpi,1:jpj,1:jpksed), iarroce(1:jpoce) )215 190 IF( ln_timing ) CALL timing_stop('sed_rst_read') 216 191 … … 247 222 IF(lwp) WRITE(numsed,*) '~~~~~~~~~' 248 223 249 250 trcsedi(:,:,:,:) = 0.0251 flxsedi3d(:,:,:,:) = 0.0252 224 zdta(:,:) = 1.0 253 225 zdta2(:,:,:) = 0.0 254 255 226 256 227 !! 1. WRITE in nutwrs 257 228 !! ------------------ 258 !zinfo(1) = REAL( kt)259 CALL iom_rstput( kt, nitrst, numrsw, 'kt', REAL( kt , wp))229 zinfo(1) = REAL( kt) 230 CALL iom_rstput( kt, nitrst, numrsw, 'kt', zinfo ) 260 231 261 232 ! Back to 2D geometry 262 233 DO jn = 1, jpsol 263 CALL unpack_arr( jpoce, trcsedi(1:jpi,1:jpj,1:jpksed,jn) , iarroce(1:jpoce), &234 CALL unpack_arr( jpoce, zdta2(1:jpi,1:jpj,1:jpksed) , iarroce(1:jpoce), & 264 235 & solcp(1:jpoce,1:jpksed,jn ) ) 236 cltra = TRIM(sedtrcd(jn)) 237 CALL iom_rstput( kt, nitrst, numrsw, TRIM(cltra), zdta2(:,:,:) ) 265 238 END DO 266 239 267 240 DO jn = 1, jpwat 268 CALL unpack_arr( jpoce, trcsedi(1:jpi,1:jpj,1:jpksed,jpsol+jn) , iarroce(1:jpoce), &241 CALL unpack_arr( jpoce, zdta2(1:jpi,1:jpj,1:jpksed) , iarroce(1:jpoce), & 269 242 & pwcp(1:jpoce,1:jpksed,jn ) ) 243 cltra = TRIM(sedtrcd(jpsol+jn)) 244 CALL iom_rstput( kt, nitrst, numrsw, TRIM(cltra), zdta2(:,:,:) ) 270 245 END DO 271 246 ! pH … … 276 251 ENDDO 277 252 278 CALL unpack_arr( jpoce, flxsedi3d(1:jpi,1:jpj,1:jpksed,1) , iarroce(1:jpoce), &253 CALL unpack_arr( jpoce, zdta2(1:jpi,1:jpj,1:jpksed) , iarroce(1:jpoce), & 279 254 & zdta(1:jpoce,1:jpksed) ) 255 cltra = TRIM(seddia3d(1)) 256 CALL iom_rstput( kt, nitrst, numrsw, TRIM(cltra), zdta2(:,:,:) ) 280 257 281 CALL unpack_arr( jpoce, flxsedi3d(1:jpi,1:jpj,1:jpksed,2) , iarroce(1:jpoce), &258 CALL unpack_arr( jpoce, zdta2(1:jpi,1:jpj,1:jpksed) , iarroce(1:jpoce), & 282 259 & co3por(1:jpoce,1:jpksed) ) 260 cltra = TRIM(seddia3d(2)) 261 CALL iom_rstput( kt, nitrst, numrsw, TRIM(cltra), zdta2(:,:,:) ) 283 262 284 263 ! prognostic variables 285 264 ! -------------------- 286 265 287 DO jn = 1, jptrased288 cltra = TRIM(sedtrcd(jn))289 CALL iom_rstput( kt, nitrst, numrsw, TRIM(cltra), trcsedi(:,:,:,jn) )290 ENDDO291 292 DO jn = 1, 2293 cltra = TRIM(seddia3d(jn))294 CALL iom_rstput( kt, nitrst, numrsw, TRIM(cltra), flxsedi3d(:,:,:,jn) )295 ENDDO296 297 266 CALL unpack_arr( jpoce, zdta2(1:jpi,1:jpj,1:jpksed) , iarroce(1:jpoce), & 298 267 & db(1:jpoce,1:jpksed) ) … … 307 276 CALL iom_rstput( kt, nitrst, numrsw, TRIM(cltra), zdta2(:,:,:) ) 308 277 309 CALL unpack_arr( jpoce, zdta2(1:jpi,1:jpj,1:jpksed) , iarroce(1:jpoce), &310 & sedligand(1:jpoce,1:jpksed) )311 312 cltra = "sedligand"313 CALL iom_rstput( kt, nitrst, numrsw, TRIM(cltra), zdta2(:,:,:) )314 315 278 IF( kt == nitrst ) THEN 316 IF(.NOT.lwxios) THEN 317 CALL iom_close( numrsw ) ! close the restart file (only at last time step) 318 ELSE 319 CALL iom_context_finalize( cw_sedrst_cxt ) 320 iom_file(numrsw)%nfid = 0 321 numrsw = 0 322 ENDIF 279 CALL iom_close( numrsw ) ! close the restart file (only at last time step) 323 280 IF( l_offline .AND. ln_rst_list ) THEN 324 281 nrst_lst = nrst_lst + 1 … … 351 308 !! In both those options, the exact duration of the experiment 352 309 !! since the beginning (cumulated duration of all previous restart runs) 353 !! is not stored in the restart and is assumed to be (nittrc000-1)*r n_Dt.310 !! is not stored in the restart and is assumed to be (nittrc000-1)*rdt. 354 311 !! This is valid is the time step has remained constant. 355 312 !! … … 363 320 REAL(wp) :: zkt, zrdttrc1 364 321 REAL(wp) :: zndastp 365 CHARACTER(len = 82) :: clpname366 322 367 323 ! Time domain : restart … … 375 331 376 332 IF( ln_rst_sed ) THEN 377 lxios_sini = .FALSE.378 333 CALL iom_open( TRIM(cn_sedrst_indir)//'/'//cn_sedrst_in, numrsr ) 379 380 IF( lrxios) THEN381 cr_sedrst_cxt = 'sed_rst'382 IF(lwp) WRITE(numout,*) 'Enable restart reading by XIOS for SED'383 ! IF( TRIM(Agrif_CFixed()) == '0' ) THEN384 ! clpname = cn_sedrst_in385 ! ELSE386 ! clpname = TRIM(Agrif_CFixed())//"_"//cn_sedrst_in387 ! ENDIF388 CALL iom_init( cr_sedrst_cxt, kdid = numrsr, ld_closedef = .TRUE. )389 ENDIF390 334 CALL iom_get ( numrsr, 'kt', zkt ) ! last time-step of previous run 335 391 336 IF(lwp) THEN 392 337 WRITE(numsed,*) ' *** Info read in restart : ' … … 401 346 ENDIF 402 347 ! Control of date 403 IF( nittrc000 - NINT( zkt ) /= nn_dtsed.AND. nn_rstsed /= 0 ) &348 IF( nittrc000 - NINT( zkt ) /= 1 .AND. nn_rstsed /= 0 ) & 404 349 & CALL ctl_stop( ' ===>>>> : problem with nittrc000 for the restart', & 405 350 & ' verify the restart file or rerun with nn_rsttr = 0 (namelist)' ) … … 414 359 ELSE 415 360 ndastp = ndate0 - 1 ! ndate0 read in the namelist in dom_nam 416 adatrj = ( REAL( nittrc000-1, wp ) * r n_Dt ) / rday361 adatrj = ( REAL( nittrc000-1, wp ) * rdt ) / rday 417 362 ! note this is wrong if time step has changed during run 418 363 ENDIF … … 435 380 IF(lwp) WRITE(numsed,*) 'trc_wri : write the TOP restart file (NetCDF) at it= ', kt, ' date= ', ndastp 436 381 IF(lwp) WRITE(numsed,*) '~~~~~~~' 437 IF( lwxios ) CALL iom_init_closedef(cw_sedrst_cxt)438 382 ENDIF 439 383 CALL iom_rstput( kt, nitrst, numrsw, 'kt' , REAL( kt , wp) ) ! time-step 440 384 CALL iom_rstput( kt, nitrst, numrsw, 'ndastp' , REAL( ndastp, wp) ) ! date 441 CALL iom_rstput( kt, nitrst, numrsw, 'adatrj' , adatrj ) ! number of elapsed days since442 ! 385 CALL iom_rstput( kt, nitrst, numrsw, 'adatrj' , adatrj ) ! number of elapsed days since 386 ! ! the begining of the run [s] 443 387 ENDIF 444 388 -
NEMO/trunk/src/TOP/PISCES/SED/sedsfc.F90
r13295 r15450 13 13 !! * Substitutions 14 14 # include "do_loop_substitute.h90" 15 !! $Id$ 15 # include "domzgr_substitute.h90" 16 16 17 CONTAINS 17 18 … … 28 29 !!* Arguments 29 30 INTEGER, INTENT(in) :: kt ! time step 30 INTEGER, INTENT(in) :: Kbb ! time index31 INTEGER, INTENT(in) :: Kbb ! time level indices 31 32 32 33 ! * local variables … … 46 47 CALL unpack_arr ( jpoce, trc_data(1:jpi,1:jpj,7), iarroce(1:jpoce), pwcp(1:jpoce,1,jwnh4) ) 47 48 CALL unpack_arr ( jpoce, trc_data(1:jpi,1:jpj,8), iarroce(1:jpoce), pwcp(1:jpoce,1,jwfe2) ) 49 CALL unpack_arr ( jpoce, trc_data(1:jpi,1:jpj,9), iarroce(1:jpoce), pwcp(1:jpoce,1,jwlgw) ) 48 50 49 50 DO_2D( 1, 1, 1, 1 ) 51 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 51 52 ikt = mbkt(ji,jj) 52 53 IF ( tmask(ji,jj,ikt) == 1 ) THEN -
NEMO/trunk/src/TOP/PISCES/SED/sedstp.F90
r13970 r15450 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 12 USE sedbtb ! bioturbation 10 USE sedsol ! Organic reactions and diffusion 13 11 USE sedadv ! vertical advection 14 USE sedmbc ! mass balance calculation15 12 USE sedsfc ! sediment surface data 16 13 USE sedrst ! restart 17 14 USE sedwri ! outputs 15 USE sedini 18 16 USE trcdmp_sed 19 17 USE lib_mpp ! distribued memory computing library … … 25 23 !! * Routine accessibility 26 24 PUBLIC sed_stp ! called by step.F90 25 26 !! * Substitutions 27 # include "do_loop_substitute.h90" 28 # include "domzgr_substitute.h90" 27 29 28 30 !! $Id$ … … 44 46 !! ! 06-04 (C. Ethe) Re-organization 45 47 !!---------------------------------------------------------------------- 46 INTEGER, INTENT(in) :: kt ! number of iteration 47 INTEGER, INTENT(in) :: Kbb, Kmm, Krhs ! time level indices 48 INTEGER :: ji,jk,js,jn,jw 48 INTEGER, INTENT(in) :: kt ! number of iteration 49 INTEGER, INTENT(in) :: Kbb, Kmm, Krhs ! time level indices 50 51 INTEGER :: ji,jk,js,jn,jw,jkmax,jsmax 49 52 !!---------------------------------------------------------------------- 50 IF( ln_timing ) CALL timing_start('sed_stp')53 IF( ln_timing ) CALL timing_start('sed_stp') 51 54 ! 52 55 CALL sed_rst_opn ( kt ) ! Open tracer restart file … … 56 59 57 60 dtsed = rDt_trc 58 ! dtsed2 = dtsed59 61 IF (kt /= nitsed000) THEN 60 CALL sed_dta( kt, Kbb, Kmm ) 62 CALL sed_dta( kt, Kbb, Kmm ) ! Load Data for bot. wat. Chem and fluxes 61 63 ENDIF 62 64 … … 66 68 ENDIF 67 69 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 70 CALL sed_sol( kt ) ! Solute diffusion and reactions 78 71 CALL sed_adv( kt ) ! advection 79 72 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 73 83 IF (ln_sed_2way) CALL sed_sfc( kt, Kbb ) 74 IF (ln_sed_2way) CALL sed_sfc( kt, Kbb ) ! Give back new bottom wat chem to tracer model 84 75 ENDIF 85 76 CALL sed_wri( kt ) ! outputs 86 77 IF( kt == nitsed000 ) THEN 87 78 CALL iom_close( numrsr ) ! close input tracer restart file 88 IF(lrxios) CALL iom_context_finalize( cr_sedrst_cxt ) 89 ! IF(lwm) CALL FLUSH( numont ) ! flush namelist output 79 ! IF(lwm) CALL FLUSH( numont ) ! flush namelist output 90 80 ENDIF 91 81 IF( lrst_sed ) CALL sed_rst_wri( kt ) ! restart file output 92 82 93 IF( kt == nitsedend ) CLOSE( numsed )83 IF( kt == nitsedend ) CLOSE( numsed ) 94 84 95 IF( ln_timing ) CALL timing_stop('sed_stp')85 IF( ln_timing ) CALL timing_stop('sed_stp') 96 86 97 87 END SUBROUTINE sed_stp -
NEMO/trunk/src/TOP/PISCES/SED/sedwri.F90
r14086 r15450 4 4 !! Sediment diagnostics : write sediment output files 5 5 !!====================================================================== 6 USE par_sed7 6 USE sed 8 7 USE sedarr 9 8 USE lib_mpp ! distribued memory computing library 10 9 USE iom 10 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 11 11 12 12 IMPLICIT NONE … … 38 38 CHARACTER(len = 20) :: cltra 39 39 REAL(wp) :: zrate 40 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdta, zflx 40 REAL(wp), DIMENSION(jpoce, jpksed) :: zdta 41 REAL(wp), DIMENSION(jpoce, jptrased+1) :: zflx 42 REAL(wp), DIMENSION(jpi, jpj, jpksed, jptrased) :: trcsedi 43 REAL(wp), DIMENSION(jpi, jpj, jpksed, jpdia3dsed) :: flxsedi3d 44 REAL(wp), DIMENSION(jpi, jpj, jpdia2dsed) :: flxsedi2d 41 45 42 46 !!------------------------------------------------------------------- … … 54 58 IF (lwp) WRITE(numsed,*) ' ' 55 59 56 ALLOCATE( zdta(jpoce,jpksed) ) ; ALLOCATE( zflx(jpoce,jpwatp1) )57 58 60 ! Initialize variables 59 61 ! -------------------- … … 88 90 CALL unpack_arr( jpoce, flxsedi3d(1:jpi,1:jpj,1:jpksed,2) , iarroce(1:jpoce), & 89 91 & co3por(1:jpoce,1:jpksed) ) 92 93 CALL unpack_arr( jpoce, flxsedi3d(1:jpi,1:jpj,1:jpksed,3) , iarroce(1:jpoce), & 94 & saturco3(1:jpoce,1:jpksed) ) 95 90 96 91 97 ! flxsedi3d = 0. … … 95 101 DO ji = 1, jpoce 96 102 zflx(ji,jw) = ( pwcp(ji,1,jw) - pwcp_dta(ji,jw) ) & 97 & * 1.e3 / 1.e2 * dzkbot(ji) / rDt_trc 103 & * 1.e3 * ( 1.e-2 * dzkbot(ji) ) / 1.E4 / rDt_trc 104 ENDDO 105 ENDDO 106 107 ! Calculation of fluxes g/cm2/s 108 DO js = 1, jpsol 109 zrate = 1.0 / rDt_trc 110 DO ji = 1, jpoce 111 zflx(ji,jpwat+js) = zflx(ji,jpwat+js) + ( tosed(ji,js) - fromsed(ji,js) ) * zrate 98 112 ENDDO 99 113 ENDDO … … 101 115 ! Calculation of accumulation rate per dt 102 116 DO js = 1, jpsol 103 zrate = 1.0 / ( denssol * por1(jpksed) ) /rDt_trc117 zrate = 1.0 / rDt_trc 104 118 DO ji = 1, jpoce 105 zflx(ji,jp watp1) = zflx(ji,jpwatp1) + ( tosed(ji,js) - fromsed(ji,js) ) * zrate119 zflx(ji,jptrased+1) = zflx(ji,jptrased+1) + ( tosed(ji,js) - fromsed(ji,js) ) * zrate 106 120 ENDDO 107 121 ENDDO 108 122 109 DO jn = 1, jpdia2dsed - 1123 DO jn = 1, jpdia2dsed - 2 110 124 CALL unpack_arr( jpoce, flxsedi2d(1:jpi,1:jpj,jn), iarroce(1:jpoce), zflx(1:jpoce,jn) ) 111 125 END DO 126 112 127 zflx(:,1) = dzdep(:) / dtsed 113 CALL unpack_arr( jpoce, flxsedi2d(1:jpi,1:jpj,jpdia2dsed ), iarroce(1:jpoce), zflx(1:jpoce,1) )128 CALL unpack_arr( jpoce, flxsedi2d(1:jpi,1:jpj,jpdia2dsed-1), iarroce(1:jpoce), zflx(1:jpoce,1) ) 114 129 115 ! Start writing data 116 ! --------------------- 117 DO jn = 1, jptrased 118 cltra = sedtrcd(jn) ! short title for 3D diagnostic 119 CALL iom_put( cltra, trcsedi(:,:,:,jn) ) 120 END DO 130 CALL unpack_arr( jpoce, flxsedi2d(1:jpi,1:jpj,jpdia2dsed), iarroce(1:jpoce), rstepros(1:jpoce) ) 131 ! 132 ! CALL lbc_lnk( 'sedwri', trcsedi(:,:,:,:), 'T', 1._wp ) 133 ! CALL lbc_lnk( 'sedwri', flxsedi3d(:,:,:,:), 'T', 1._wp ) 134 ! CALL lbc_lnk( 'sedwri', flxsedi2d(:,:,:), 'T', 1._wp ) 121 135 122 DO jn = 1, jpdia3dsed 123 cltra = seddia3d(jn) ! short title for 3D diagnostic 124 CALL iom_put( cltra, flxsedi3d(:,:,:,jn) ) 125 END DO 136 ! Start writing data 137 ! --------------------- 138 DO jn = 1, jptrased 139 cltra = sedtrcd(jn) ! short title for 3D diagnostic 140 CALL iom_put( cltra, trcsedi(:,:,:,jn) ) 141 END DO 126 142 127 DO jn = 1, jpdia2dsed128 cltra = seddia2d(jn) ! short title for 2D diagnostic129 CALL iom_put( cltra, flxsedi2d(:,:,jn) )130 143 DO jn = 1, jpdia3dsed 144 cltra = seddia3d(jn) ! short title for 3D diagnostic 145 CALL iom_put( cltra, flxsedi3d(:,:,:,jn) ) 146 END DO 131 147 132 133 DEALLOCATE( zdta ) ; DEALLOCATE( zflx ) 148 DO jn = 1, jpdia2dsed 149 cltra = seddia2d(jn) ! short title for 2D diagnostic 150 CALL iom_put( cltra, flxsedi2d(:,:,jn) ) 151 END DO 134 152 135 153 IF( ln_timing ) CALL timing_stop('sed_wri') -
NEMO/trunk/src/TOP/PISCES/SED/trcdmp_sed.F90
r15023 r15450 93 93 CALL trc_dta( kt, jl, ztrcdta ) ! read tracer data at nit000 94 94 ! 95 DO_2D( 1, 1, 1, 1)95 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 96 96 ikt = mbkt(ji,jj) 97 97 tr(ji,jj,ikt,jn,Kbb) = ztrcdta(ji,jj,ikt) + ( tr(ji,jj,ikt,jn,Kbb) - ztrcdta(ji,jj,ikt) ) & -
NEMO/trunk/src/TOP/PISCES/trcini_pisces.F90
r14086 r15450 280 280 ! Initialization of the sediment model 281 281 IF( ln_sediment) & 282 & CALL sed_ini t! Initialization of the sediment model282 & CALL sed_ini ! Initialization of the sediment model 283 283 284 284 CALL p4z_sed_init ! loss of organic matter in the sediments
Note: See TracChangeset
for help on using the changeset viewer.
