Changeset 15075
- Timestamp:
- 2021-07-02T17:15:57+02:00 (3 years ago)
- Location:
- NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src
- Files:
-
- 5 added
- 23 edited
-
OFF/nemogcm.F90 (modified) (2 diffs)
-
TOP/PISCES/P4Z/p4zbc.F90 (modified) (4 diffs)
-
TOP/PISCES/P4Z/p4zrem.F90 (modified) (2 diffs)
-
TOP/PISCES/SED/par_sed.F90 (modified) (2 diffs)
-
TOP/PISCES/SED/sed.F90 (modified) (11 diffs)
-
TOP/PISCES/SED/sedadv.F90 (modified) (18 diffs)
-
TOP/PISCES/SED/sedarr.F90 (modified) (1 diff)
-
TOP/PISCES/SED/sedbtb.F90 (modified) (5 diffs)
-
TOP/PISCES/SED/sedchem.F90 (modified) (19 diffs)
-
TOP/PISCES/SED/seddiff.F90 (modified) (1 diff)
-
TOP/PISCES/SED/seddsr.F90 (modified) (7 diffs)
-
TOP/PISCES/SED/seddsrjac.F90 (added)
-
TOP/PISCES/SED/seddta.F90 (modified) (8 diffs)
-
TOP/PISCES/SED/sedfunc.F90 (added)
-
TOP/PISCES/SED/sedini.F90 (modified) (33 diffs)
-
TOP/PISCES/SED/sedinitrc.F90 (modified) (3 diffs)
-
TOP/PISCES/SED/sedinorg.F90 (modified) (5 diffs)
-
TOP/PISCES/SED/sedjac.F90 (added)
-
TOP/PISCES/SED/sedmat.F90 (modified) (9 diffs)
-
TOP/PISCES/SED/sedorg.F90 (added)
-
TOP/PISCES/SED/sedrst.F90 (modified) (13 diffs)
-
TOP/PISCES/SED/sedsfc.F90 (modified) (3 diffs)
-
TOP/PISCES/SED/sedsol.F90 (modified) (13 diffs)
-
TOP/PISCES/SED/sedstp.F90 (modified) (5 diffs)
-
TOP/PISCES/SED/sedwri.F90 (modified) (6 diffs)
-
TOP/PISCES/SED/trosk.F90 (added)
-
TOP/PISCES/trcini_pisces.F90 (modified) (1 diff)
-
TOP/trcstp.F90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OFF/nemogcm.F90
r14448 r15075 123 123 ! 124 124 DO WHILE ( istp <= nitend .AND. nstop == 0 ) !== OFF time-stepping ==! 125 126 125 IF( ln_timing ) THEN 127 126 zstptiming = MPI_Wtime() … … 153 152 Naa = Nrhs 154 153 ! 155 # if ! defined key_qco154 # if ! defined key_qco 156 155 IF( .NOT.ln_linssh ) CALL dta_dyn_sf_interp( istp, Nnn ) ! calculate now grid parameters 157 # endif156 # endif 158 157 159 158 #else 160 159 CALL dta_dyn_sed( istp, Nnn ) ! Interpolation of the dynamical fields 161 160 CALL trc_stp ( istp, Nbb, Nnn, Nrhs, Naa ) ! time-stepping 161 ! Swap time levels 162 Nrhs = Nbb 163 Nbb = Nnn 164 Nnn = Naa 165 Naa = Nrhs 162 166 #endif 163 167 CALL stp_ctl ( istp ) ! Time loop: control and print -
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/P4Z/p4zbc.F90
r14786 r15075 24 24 LOGICAL , PUBLIC :: ln_ironsed !: boolean for Fe input from sediments 25 25 LOGICAL , PUBLIC :: ln_hydrofe !: boolean for Fe input from hydrothermal vents 26 LOGICAL , PUBLIC :: ln_dust_inp !: boolean for Fe input from hydrothermal vents 26 27 REAL(wp), PUBLIC :: sedfeinput !: Coastal release of Iron 27 28 REAL(wp), PUBLIC :: icefeinput !: Iron concentration in sea ice … … 97 98 98 99 zwdust = 0.03 / ( wdust / rday ) / ( 250. * rday ) 99 DO_3D( 1, 1, 1, 1, 2, jpkm1 ) 100 zdustdep = dust(ji,jj) * zwdust * rfact * EXP( -gdept(ji,jj,jk,Kmm) /( 250. * wdust ) ) 101 ! 102 tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) + zdustdep * mfrac / mMass_Fe 103 tr(ji,jj,jk,jppo4,Krhs) = tr(ji,jj,jk,jppo4,Krhs) + zdustdep * 1.e-3 / mMass_P 104 tr(ji,jj,jk,jpsil,Krhs) = tr(ji,jj,jk,jpsil,Krhs) + zdustdep * 0.269 / mMass_Si 105 END_3D 100 101 ! Atmospheric input of Iron dissolves in the water column 102 IF ( ln_trc_sbc(jpfer) ) THEN 103 DO_3D( 1, 1, 1, 1, 2, jpkm1 ) 104 zdustdep = dust(ji,jj) * zwdust * rfact * EXP( -gdept(ji,jj,jk,Kmm) /( 250. * wdust ) ) 105 ! 106 tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) + zdustdep * mfrac / mMass_Fe 107 END_3D 108 109 IF( lk_iomput ) THEN 110 ! surface downward dust depo of iron 111 jl = n_trc_indsbc(jpfer) 112 CALL iom_put( "Irondep", ( rf_trsfac(jl) * sf_trcsbc(jl)%fnow(:,:,1) / rn_sbc_time ) * 1.e+3 * tmask(:,:,1) ) 113 114 ENDIF 115 116 ENDIF 117 118 ! Atmospheric input of PO4 dissolves in the water column 119 IF ( ln_trc_sbc(jppo4) ) THEN 120 DO_3D( 1, 1, 1, 1, 2, jpkm1 ) 121 zdustdep = dust(ji,jj) * zwdust * rfact * EXP( -gdept(ji,jj,jk,Kmm) /( 250. * wdust ) ) 122 ! 123 tr(ji,jj,jk,jppo4,Krhs) = tr(ji,jj,jk,jppo4,Krhs) + zdustdep * 1.e-3 / mMass_P 124 END_3D 125 ENDIF 126 127 ! Atmospheric input of Si dissolves in the water column 128 IF ( ln_trc_sbc(jpsil) ) THEN 129 DO_3D( 1, 1, 1, 1, 2, jpkm1 ) 130 zdustdep = dust(ji,jj) * zwdust * rfact * EXP( -gdept(ji,jj,jk,Kmm) /( 250. * wdust ) ) 131 ! 132 tr(ji,jj,jk,jpsil,Krhs) = tr(ji,jj,jk,jpsil,Krhs) + zdustdep * 0.269 / mMass_Si 133 END_3D 134 ENDIF 135 106 136 ! 107 137 IF( lk_iomput ) THEN 108 ! surface downward dust depo of iron109 jl = n_trc_indsbc(jpfer)110 CALL iom_put( "Irondep", ( rf_trsfac(jl) * sf_trcsbc(jl)%fnow(:,:,1) / rn_sbc_time ) * 1.e+3 * tmask(:,:,1) )111 112 138 ! dust concentration at surface 113 139 CALL iom_put( "pdust" , dust(:,:) / ( wdust / rday ) * tmask(:,:,1) ) ! dust concentration at surface … … 225 251 !! 226 252 NAMELIST/nampisbc/cn_dir, sn_dust, sn_ironsed, sn_hydrofe, & 227 & ln_ironsed, ln_ironice, ln_hydrofe, &228 & sedfeinput, distcoast, icefeinput, wdust, mfrac, &253 & ln_ironsed, ln_ironice, ln_hydrofe, ln_dust_inp, & 254 & sedfeinput, distcoast, icefeinput, wdust, mfrac, & 229 255 & hratio, lgw_rath 230 256 !!---------------------------------------------------------------------- … … 267 293 END IF 268 294 269 ll_bc = ( ln_trcbc .AND. lltrcbc ) .OR. ln_hydrofe .OR. ln_ironsed .OR. ln_ironice 270 ll_dust = ln_ trc_sbc(jpfer)295 ll_bc = ( ln_trcbc .AND. lltrcbc ) .OR. ln_hydrofe .OR. ln_ironsed .OR. ln_ironice .OR. ln_dust_inp 296 ll_dust = ln_dust_inp 271 297 ll_ndepo = ln_trc_sbc(jpno3) .OR. ln_trc_sbc(jpnh4) 272 298 ll_river = ln_trc_cbc(jpno3) -
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/P4Z/p4zrem.F90
r14963 r15075 104 104 zdepmin = MIN( 1., zdep / gdept(ji,jj,jk,Kmm) ) 105 105 zdepbac (ji,jj,jk) = zdepmin**0.683 * ztempbac(ji,jj) 106 zdepeff(ji,jj,jk) = zdepeff(ji,jj,jk) * zdepmin**0.3 106 ! zdepeff(ji,jj,jk) = zdepeff(ji,jj,jk) * zdepmin**0.3 107 zdepeff(ji,jj,jk) = zdepeff(ji,jj,jk) * zdepmin**0.6 107 108 ENDIF 108 109 END_3D … … 187 188 ! is treated here. The GGE of bacteria supposed to be equal to 188 189 ! 0.33. This is hard-coded. 189 tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) - zbactfer*0.1 2190 tr(ji,jj,jk,jpsfe,Krhs) = tr(ji,jj,jk,jpsfe,Krhs) + zbactfer*0. 10190 tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) - zbactfer*0.1 191 tr(ji,jj,jk,jpsfe,Krhs) = tr(ji,jj,jk,jpsfe,Krhs) + zbactfer*0.08 191 192 tr(ji,jj,jk,jpbfe,Krhs) = tr(ji,jj,jk,jpbfe,Krhs) + zbactfer*0.02 192 zfebact(ji,jj,jk) = zbactfer * 0.1 2193 zfebact(ji,jj,jk) = zbactfer * 0.1 193 194 blim(ji,jj,jk) = xlimbacl(ji,jj,jk) * zdepbac(ji,jj,jk) / 1.e-6 194 195 END_3D -
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/par_sed.F90
r14385 r15075 20 20 jp_sal => jp_sal !: indice of salintity 21 21 22 INTEGER, P UBLIC, PARAMETER :: jpdta = 1722 INTEGER, PARAMETER :: jpdta = 18 23 23 24 24 ! Vertical sediment geometry 25 INTEGER, PUBLIC :: jpksed = 11 26 INTEGER, PUBLIC :: jpksedm1 25 INTEGER, PUBLIC :: & 26 jpksed = 11 , & 27 jpksedm1 = 10 27 28 28 29 ! sediment tracer species 29 INTEGER, P UBLIC, PARAMETER :: &30 INTEGER, PARAMETER :: & 30 31 jpsol = 8, & !: number of solid component 31 jpwat = 10, & !: number of pore water component 32 jpwat = 11, & !: number of pore water component 33 jpads = 2 , & !! number adsorbed species 32 34 jpwatp1 = jpwat +1, & 33 35 jpsol1 = jpsol - 1 … … 35 37 36 38 ! pore water components 37 INTEGER, PUBLIC, PARAMETER :: & 38 jwsil = 1, & !: silic acid 39 jwoxy = 2, & !: oxygen 40 jwdic = 3, & !: dissolved inorganic carbon 41 jwno3 = 4, & !: nitrate 42 jwpo4 = 5, & !: phosphate 43 jwalk = 6, & !: alkalinity 44 jwnh4 = 7, & !: Ammonium 45 jwh2s = 8, & !: Sulfate 46 jwso4 = 9, & !: H2S 47 jwfe2 = 10 !: Fe2+ 39 INTEGER, PARAMETER :: & 40 jwoxy = 1, & !: oxygen 41 jwno3 = 2, & !: nitrate 42 jwpo4 = 3, & !: phosphate 43 jwnh4 = 4, & !: Ammonium 44 jwh2s = 5, & !: Sulfate 45 jwso4 = 6, & !: H2S 46 jwfe2 = 7, & !: Fe2+ 47 jwalk = 8, & !: Alkalinity 48 jwlgw = 9, & !: Alkalinity 49 jwdic = 10, & !: DIC 50 jwsil = 11 !: Silicate 48 51 49 52 ! solid components 50 INTEGER, P UBLIC, PARAMETER :: &51 js opal = 1, & !: opal sediment52 js clay = 2, & !: clay53 INTEGER, PARAMETER :: & 54 jsfeo = 1, & !: iron hydroxides 55 jsfes = 2, & !: FeS 53 56 jspoc = 3, & !: organic carbon 54 js cal = 4, & !: calcite55 jspo s = 5, & !: semi-refPOC56 js por = 6, & !: refractory POC57 js feo = 7, & !: iron hydroxides58 js fes = 8 !: FeS57 jspos = 4, & !: semi-ref POC 58 jspor = 5, & !: refractory POC 59 jscal = 6, & !: Calcite 60 jsopal = 7, & !: Opal 61 jsclay = 8 !: clay 59 62 60 INTEGER, P UBLIC, PARAMETER :: &63 INTEGER, PARAMETER :: & 61 64 jptrased = jpsol + jpwat , & 62 jpdia3dsed = 4 , & 63 jpdia2dsed = 20 65 jpvode = jptrased - 11 , & 66 jpdia3dsed = 3 , & 67 jpdia2dsed = 21 64 68 65 69 END MODULE par_sed -
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/sed.F90
r14385 r15075 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 44 REAL(wp), PUBLIC :: denssol !: density of solid material 45 45 LOGICAL , PUBLIC :: lrst_sed !: logical to control the trc restart write … … 57 57 58 58 ! 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 59 REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: pwcp, pwcpa 60 REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: solcp, solcpa 61 REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: Jacobian 62 REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: diff 64 63 65 64 !! * Shared module variables … … 70 69 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: fromsed !: 71 70 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: tosed !: 71 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: rearatpom !: 72 ! 72 73 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: temp !: temperature 73 74 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: salt !: salinity … … 75 76 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: fecratio !: Fe/C ratio in falling particles to the sediments 76 77 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: dzdep !: total thickness of solid material rained [cm] in each cell 77 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: zkbot !: Total water column depth78 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: wacc 78 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: zkbot !: total thickness of solid material rained [cm] in each cell 79 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: wacc, wacc1 !: total thickness of solid material rained [cm] in each cell 79 80 ! 80 81 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: hipor !: [h+] in mol/kg*densSW … … 83 84 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: volw3d !: ??? 84 85 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: vols3d !: ??? 86 REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: volc 87 85 88 86 89 !! Chemistry 87 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: densSW 88 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: borats 90 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: densSW 91 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: borats 89 92 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: calcon2 90 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: akbs 91 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: ak1s 92 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: ak2s 93 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: akws 94 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: ak12s 95 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: ak1ps 96 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: ak2ps 97 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: ak3ps 98 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: ak12ps 93 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: akbs 94 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: ak1s 95 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: ak2s 96 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: akws 97 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: ak12s 98 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: ak1ps 99 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: ak2ps 100 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: ak3ps 101 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: ak12ps 99 102 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: ak123ps 100 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: aksis 101 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: aksps 103 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: aksis 104 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: aknh3 105 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: aksps 106 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: akh2s 102 107 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: sieqs 103 108 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: aks3s … … 112 117 INTEGER , PUBLIC, SAVE :: jpoce, indoce !: Ocean points ( number/indices ) 113 118 INTEGER , PUBLIC, DIMENSION(: ), ALLOCATABLE :: iarroce !: Computation of 1D array of sediments points 119 INTEGER , PUBLIC, DIMENSION(:, : ), ALLOCATABLE :: jarr !: Computation of 1D array of sediments points 120 INTEGER , PUBLIC, DIMENSION(: ), ALLOCATABLE :: jsvode, isvode !: Computation of 1D array of sediments points 121 114 122 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: epkbot !: ocean bottom layer thickness 115 123 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: gdepbot !: Depth of the sediment … … 122 130 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: db !: bioturbation ceofficient 123 131 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: irrig !: bioturbation ceofficient 124 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: sedligand 125 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: saturco3 132 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: radsnh4, radsfe2 133 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: saturco3 134 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: rdtsed !: sediment model time-step 135 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: stepros 126 136 REAL(wp) :: dens !: density of solid material 127 137 !! Inputs / Outputs … … 133 143 CHARACTER( len = 20 ), DIMENSION(jpdia2dsed) :: seddia2d, seddia2u 134 144 ! 135 REAL(wp), PUBLIC, DIMENSION(:,:,:,:), ALLOCATABLE :: trcsedi136 REAL(wp), PUBLIC, DIMENSION(:,:,:,:), ALLOCATABLE :: flxsedi3d137 REAL(wp), PUBLIC, DIMENSION(:,:,: ), ALLOCATABLE :: flxsedi2d138 145 139 146 INTEGER, PUBLIC :: numsed = 27 ! units … … 149 156 !!------------------------------------------------------------------- 150 157 ! 151 ALLOCATE( trc_data(jpi,jpj,jpdta), trcsedi (jpi,jpj,jpksed,jptrased) , & 152 & epkbot(jpi,jpj), gdepbot(jpi,jpj), dz(jpksed), por(jpksed) , & 153 & por1(jpksed), volw(jpksed), vols(jpksed), mol_wgt(jpsol) , & 154 & flxsedi2d(jpi,jpj,jpdia2dsed), flxsedi3d(jpi,jpj,jpksed,jpdia3dsed), STAT=sed_alloc ) 158 ALLOCATE( trc_data(jpi,jpj,jpdta) , & 159 & epkbot(jpi,jpj), gdepbot(jpi,jpj) , & 160 & dz(jpksed) , por(jpksed) , por1(jpksed) , & 161 & volw(jpksed), vols(jpksed), rdtsed(jpksed) , & 162 & mol_wgt(jpsol), STAT=sed_alloc ) 155 163 156 164 IF( sed_alloc /= 0 ) CALL ctl_stop( 'STOP', 'sed_alloc: failed to allocate arrays' ) -
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/sedadv.F90
r14385 r15075 51 51 INTEGER :: jn, ntimes, nztime, ikwneg 52 52 53 REAL(wp), DIMENSION(jpksed,jpsol) :: zsolcpno 53 REAL(wp), DIMENSION(jpoce,jpksed,jpsol+jpads) :: zsolcp 54 REAL(wp), DIMENSION(jpksed,jpsol+jpads) :: zsolcpno 54 55 REAL(wp), DIMENSION(jpksed) :: zfilled, zfull, zfromup, zempty 55 56 REAL(wp), DIMENSION(jpoce,jpksed) :: zgap, zwb 56 REAL(wp), DIMENSION(jpoce,jpsol ) :: zrainrf57 REAL(wp), DIMENSION(jpoce,jpsol+jpads) :: zrainrf 57 58 REAL(wp), DIMENSION(: ), ALLOCATABLE :: zraipush 58 59 59 REAL(wp) :: zkwnup, zkwnlo, zfrac, zfromce, zrest , sumtot, zsumtot160 REAL(wp) :: zkwnup, zkwnlo, zfrac, zfromce, zrest 60 61 61 62 !------------------------------------------------------------------------ … … 84 85 ENDIF 85 86 87 ! Allocation of the temporary arrays 88 ! ---------------------------------- 89 zsolcp(:,:,:) = 0._wp 90 DO js = 1, jpsol 91 zsolcp(:,:,js) = solcp(:,:,js) 92 END DO 93 DO jk = 2, jpksed 94 zsolcp(:,jk,jpsol+1) = pwcp(:,jk,jwnh4) * adsnh4 95 zsolcp(:,jk,jpsol+2) = pwcp(:,jk,jwfe2) * adsfe2 96 END DO 97 86 98 ! Initialization of data for mass balance calculation 87 99 !--------------------------------------------------- … … 97 109 zgap(:,:) = 0. 98 110 DO js = 1, jpsol 99 DO jk = 1, jpksed 100 DO ji = 1, jpoce 101 zgap(ji,jk) = zgap(ji,jk) + solcp(ji,jk,js) 102 END DO 103 ENDDO 111 zgap(:,:) = zgap(:,:) + zsolcp(:,:,js) 104 112 ENDDO 105 106 zgap(1:jpoce,1:jpksed) = 1. - zgap(1:jpoce,1:jpksed) 113 zgap(:,:) = 1. - zgap(:,:) 107 114 108 115 ! Initiate burial rates … … 111 118 DO jk = 2, jpksed 112 119 zfrac = dtsed / ( denssol * por1(jk) ) 120 zwb(:,jk) = zfrac * raintg(:) 121 ENDDO 122 zwb(:,2) = zwb(:,2) - zgap(:,2) * dz(2) 123 124 DO jk = 3, jpksed 125 zfrac = por1(jk-1) / por1(jk) 126 zwb(:,jk) = zwb(:,jk-1) * zfrac - zgap(:,jk) * dz(jk) 127 ENDDO 128 129 zrainrf(:,:) = 0. 130 DO js = 1, jpsol 113 131 DO ji = 1, jpoce 114 zwb(ji,jk) = zfrac *raintg(ji)132 IF( raintg(ji) /= 0. ) zrainrf(ji,js) = rainrg(ji,js) / raintg(ji) 115 133 END DO 116 134 ENDDO 117 118 119 DO ji = 1, jpoce120 zwb(ji,2) = zwb(ji,2) - zgap(ji,2) * dz(2)121 ENDDO122 123 DO jk = 3, jpksed124 zfrac = por1(jk-1) / por1(jk)125 DO ji = 1, jpoce126 zwb(ji,jk) = zwb(ji,jk-1) * zfrac - zgap(ji,jk) * dz(jk)127 END DO128 ENDDO129 130 zrainrf(:,:) = 0.131 DO ji = 1, jpoce132 IF( raintg(ji) /= 0. ) &133 & zrainrf(ji,:) = rainrg(ji,:) / raintg(ji)134 ENDDO135 136 135 137 136 ! Computation of full and empty solid fraction in each layer … … 147 146 DO js = 1, jpsol 148 147 DO jk = 2, jpksed 149 zfilled(jk) = zfilled(jk) + solcp(ji,jk,js)150 END DO148 zfilled(jk) = zfilled(jk) + zsolcp(ji,jk,js) 149 END DO 151 150 ENDDO 152 151 153 DO js = 1, jpsol 152 DO js = 1, jpsol + jpads 154 153 DO jk = 2, jpksed 155 zsolcpno(jk,js) = solcp(ji,jk,js) / zfilled(jk)156 END DO154 zsolcpno(jk,js) = zsolcp(ji,jk,js) / zfilled(jk) 155 END DO 157 156 ENDDO 158 157 … … 168 167 zempty(jpksed) = 1. - zfull(jpksed) 169 168 DO jk = jpksedm1, 2, -1 170 zfull (jk) = zfilled(jk) 171 zfull (jk) = zfull(jk) - zempty(jk+1) * dvolsp(jk) 169 zfull (jk) = zfilled(jk) - zempty(jk+1) * dvolsp(jk) 172 170 zempty(jk) = 1. - zfull(jk) 173 171 ENDDO … … 176 174 !-------------------------------------- 177 175 ! push entire sediment column downward to account rest of rain 178 DO js = 1, jpsol 176 DO js = 1, jpsol + jpads 179 177 DO jk = jpksed, 3, -1 180 solcp(ji,jk,js) = zfull(jk) * zsolcpno(jk,js) + zempty(jk) * zsolcpno(jk-1,js)181 ENDDO 182 183 solcp(ji,2,js) = zfull(2) * zsolcpno(2,js) + zempty(2) * zrainrf(ji,js)178 zsolcp(ji,jk,js) = zfull(jk) * zsolcpno(jk,js) + zempty(jk) * zsolcpno(jk-1,js) 179 ENDDO 180 181 zsolcp(ji,2,js) = zfull(2) * zsolcpno(2,js) + zempty(2) * zrainrf(ji,js) 184 182 185 183 DO jk = 2, jpksed 186 zsolcpno(jk,js) = solcp(ji,jk,js)184 zsolcpno(jk,js) = zsolcp(ji,jk,js) 187 185 END DO 188 186 ENDDO … … 196 194 zfromup(jk) = zwb(ji,jpksed) * ckpor(jk) / dz(jk) 197 195 ENDDO 198 DO js = 1, jpsol 196 DO js = 1, jpsol + jpads 199 197 zfromce = 1. - zfromup(2) 200 solcp(ji,2,js) = zfromce * zsolcpno(2,js) + zfromup(2) * zrainrf(ji,js)198 zsolcp(ji,2,js) = zfromce * zsolcpno(2,js) + zfromup(2) * zrainrf(ji,js) 201 199 DO jk = 3, jpksed 202 200 zfromce = 1. - zfromup(jk) 203 solcp(ji,jk,js) = zfromce * zsolcpno(jk,js) + zfromup(jk) * zsolcpno(jk-1,js)201 zsolcp(ji,jk,js) = zfromce * zsolcpno(jk,js) + zfromup(jk) * zsolcpno(jk-1,js) 204 202 ENDDO 205 203 fromsed(ji,js) = 0. … … 226 224 227 225 DO jn = 1, ntimes 228 DO js = 1, jpsol 226 DO js = 1, jpsol + jpads 229 227 DO jk = 2, jpksed 230 zsolcpno(jk,js) = solcp(ji,jk,js)228 zsolcpno(jk,js) = zsolcp(ji,jk,js) 231 229 END DO 232 230 ENDDO … … 237 235 ENDDO 238 236 239 DO js = 1, jpsol 237 DO js = 1, jpsol + jpads 240 238 zfromce = 1. - zfromup(2) 241 solcp(ji,2,js) = zfromce * zsolcpno(2,js) + zfromup(2) * zrainrf(ji,js)239 zsolcp(ji,2,js) = zfromce * zsolcpno(2,js) + zfromup(2) * zrainrf(ji,js) 242 240 DO jk = 3, jpksed 243 241 zfromce = 1. - zfromup(jk) 244 solcp(ji,jk,js) = zfromce * zsolcpno(jk,js) + zfromup(jk) * zsolcpno(jk-1,js)242 zsolcp(ji,jk,js) = zfromce * zsolcpno(jk,js) + zfromup(jk) * zsolcpno(jk-1,js) 245 243 ENDDO 246 244 fromsed(ji,js) = 0. … … 257 255 zempty(2) = 1. - zfull(2) 258 256 DO jk = 3, jpksed 259 zfull (jk) = zfilled(jk) 260 zfull (jk) = zfull (jk) - zempty(jk-1) * dvolsm(jk) 257 zfull (jk) = zfilled(jk) - zempty(jk-1) * dvolsm(jk) 261 258 zempty(jk) = 1. - zfull(jk) 262 259 ENDDO 263 260 264 261 ! fill boxes with weight fraction from underlying box 265 DO js = 1, jpsol 262 DO js = 1, jpsol + jpads 266 263 DO jk = 2, jpksedm1 267 solcp(ji,jk,js) = zfull(jk) * zsolcpno(jk,js) + zempty(jk) * zsolcpno(jk+1,js)264 zsolcp(ji,jk,js) = zfull(jk) * zsolcpno(jk,js) + zempty(jk) * zsolcpno(jk+1,js) 268 265 END DO 269 solcp(ji,jpksed,js) = zsolcpno(jpksed,js) * zfull(jpksed)266 zsolcp(ji,jpksed,js) = zsolcpno(jpksed,js) * zfull(jpksed) 270 267 tosed (ji,js) = 0. 271 268 fromsed(ji,js) = 0. 272 269 ENDDO 273 270 ! for the last layer, one make go up clay 274 solcp(ji,jpksed,jsclay) =solcp(ji,jpksed,jsclay) + zempty(jpksed) * 1.271 zsolcp(ji,jpksed,jsclay) = zsolcp(ji,jpksed,jsclay) + zempty(jpksed) * 1. 275 272 fromsed(ji,jsclay) = zempty(jpksed) * 1. * por1clay 276 273 ELSE ! rain > 0 and rain < total reaction loss 277 278 274 279 275 DO jk = 2, jpksed … … 295 291 !------------------------------------------------------------- 296 292 IF( ikwneg == 2 ) THEN ! advection is reversed in the first sediment layer 297 298 293 zkwnup = dtsed * raintg(ji) / ( denssol * por1(ikwneg) * dz(ikwneg) ) 299 294 zkwnlo = ABS( zwb(ji,ikwneg) ) / dz(ikwneg) … … 304 299 zempty(jk) = 1. - zfull(jk) 305 300 ENDDO 306 DO js = 1, jpsol 307 solcp(ji,2,js) = zfull(2) * zsolcpno(2,js)+ zkwnlo * zsolcpno(3,js) &301 DO js = 1, jpsol + jpads 302 zsolcp(ji,2,js) = zfull(2) * zsolcpno(2,js)+ zkwnlo * zsolcpno(3,js) & 308 303 & + zkwnup * zrainrf(ji,js) 309 304 DO jk = 3, jpksedm1 310 solcp(ji,jk,js) = zfull(jk) * zsolcpno(jk,js) + zempty(jk) * zsolcpno(jk+1,js)311 ENDDO 312 solcp(ji,jpksed,js) = zfull(jpksed) * zsolcpno(jpksed,js)305 zsolcp(ji,jk,js) = zfull(jk) * zsolcpno(jk,js) + zempty(jk) * zsolcpno(jk+1,js) 306 ENDDO 307 zsolcp(ji,jpksed,js) = zfull(jpksed) * zsolcpno(jpksed,js) 313 308 tosed(ji,js) = 0. 314 309 fromsed(ji,js) = 0. 315 310 ENDDO 316 solcp(ji,jpksed,jsclay) = solcp(ji,jpksed,jsclay) + zempty(jpksed) * 1. 317 !! C. Heinze fromsed(ji,jsclay) = zempty(jpksed) * 1. * denssol * por1(jpksed) / mol_wgt(jsclay) 311 zsolcp(ji,jpksed,jsclay) = zsolcp(ji,jpksed,jsclay) + zempty(jpksed) * 1. 318 312 fromsed(ji,jsclay) = zempty(jpksed) * 1. * por1clay 319 313 … … 328 322 zempty(jk) = 1. - zfull(jk) 329 323 ENDDO 330 DO js = 1, jpsol 331 solcp(ji,2,js) = zfull(2) * zsolcpno(2,js) + zempty(2) * zrainrf(ji,js)332 ENDDO 333 DO js = 1, jpsol 324 DO js = 1, jpsol + jpads 325 zsolcp(ji,2,js) = zfull(2) * zsolcpno(2,js) + zempty(2) * zrainrf(ji,js) 326 ENDDO 327 DO js = 1, jpsol + jpads 334 328 DO jk = jpksedm1, 3, -1 335 solcp(ji,jk,js) = zfull(jk) * zsolcpno(jk,js) + zempty(jk) * zsolcpno(jk-1,js)336 ENDDO 337 solcp(ji,jpksed,js) = zfull(jpksed) * zsolcpno(jpksed,js) &329 zsolcp(ji,jk,js) = zfull(jk) * zsolcpno(jk,js) + zempty(jk) * zsolcpno(jk-1,js) 330 ENDDO 331 zsolcp(ji,jpksed,js) = zfull(jpksed) * zsolcpno(jpksed,js) & 338 332 & + zkwnup * zsolcpno(jpksedm1,js) 339 333 tosed(ji,js) = 0. 340 334 fromsed(ji,js) = 0. 341 335 ENDDO 342 solcp(ji,jpksed,jsclay) = solcp(ji,jpksed,jsclay) + zkwnlo * 1. 343 ! Heinze fromsed(ji,jsclay) = zkwnlo * 1. * denssol * por1(jpksed) / mol_wgt(jsclay) 336 zsolcp(ji,jpksed,jsclay) = zsolcp(ji,jpksed,jsclay) + zkwnlo * 1. 344 337 fromsed(ji,jsclay) = zkwnlo * 1.* por1clay 345 338 ELSE IF( ikwneg > 2 .AND. ikwneg < jpksed-1 ) THEN … … 356 349 zempty(jk) = 1. - zfull(jk) 357 350 ENDDO 358 DO js = 1, jpsol 359 solcp(ji,2,js) = zfull(2) * zsolcpno(2,js) + zempty(2) * zrainrf(ji,js)360 ENDDO 361 DO js = 1, jpsol 351 DO js = 1, jpsol + jpads 352 zsolcp(ji,2,js) = zfull(2) * zsolcpno(2,js) + zempty(2) * zrainrf(ji,js) 353 ENDDO 354 DO js = 1, jpsol + jpads 362 355 DO jk = ikwneg-1, 3, -1 363 solcp(ji,jk,js) = zfull(jk) * zsolcpno(jk,js) + zempty(jk) * zsolcpno(jk-1,js)356 zsolcp(ji,jk,js) = zfull(jk) * zsolcpno(jk,js) + zempty(jk) * zsolcpno(jk-1,js) 364 357 ENDDO 365 358 ENDDO 366 359 ELSE ! ikw = 3 367 368 360 369 361 zfull (2) = zfilled(2) - zkwnup * dvolsm(3) 370 362 zempty(2) = 1. - zfull(2) 371 DO js = 1, jpsol 372 solcp(ji,2,js) = zfull(2) * zsolcpno(2,js) + zempty(2) * zrainrf(ji,js)363 DO js = 1, jpsol + jpads 364 zsolcp(ji,2,js) = zfull(2) * zsolcpno(2,js) + zempty(2) * zrainrf(ji,js) 373 365 ENDDO 374 366 ENDIF … … 382 374 zempty(jk) = 1. - zfull(jk) 383 375 ENDDO 384 DO js = 1, jpsol 376 DO js = 1, jpsol + jpads 385 377 DO jk = ikwneg+1, jpksedm1 386 solcp(ji,jk,js) = zfull(jk) * zsolcpno(jk,js) + zempty(jk) * zsolcpno(jk+1,js)378 zsolcp(ji,jk,js) = zfull(jk) * zsolcpno(jk,js) + zempty(jk) * zsolcpno(jk+1,js) 387 379 ENDDO 388 solcp(ji,jpksed,js) = zfull(jpksed) * zsolcpno(jpksed,js)389 ENDDO 390 solcp(ji,jpksed,jsclay) =solcp(ji,jpksed,jsclay) + zempty(jpksed) * 1.380 zsolcp(ji,jpksed,js) = zfull(jpksed) * zsolcpno(jpksed,js) 381 ENDDO 382 zsolcp(ji,jpksed,jsclay) = zsolcp(ji,jpksed,jsclay) + zempty(jpksed) * 1. 391 383 ELSE 392 384 393 385 zfull (jpksed) = zfilled(jpksed) - zkwnlo * dvolsm(jpksed) 394 386 zempty(jpksed) = 1. - zfull(jpksed) 395 DO js = 1, jpsol 396 solcp(ji,jpksed,js) = zfull(jpksed) * zsolcpno(jpksed,js)397 ENDDO 398 solcp(ji,jpksed,jsclay) =solcp(ji,jpksed,jsclay) + zempty(jpksed) * 1.387 DO js = 1, jpsol + jpads 388 zsolcp(ji,jpksed,js) = zfull(jpksed) * zsolcpno(jpksed,js) 389 ENDDO 390 zsolcp(ji,jpksed,jsclay) = zsolcp(ji,jpksed,jsclay) + zempty(jpksed) * 1. 399 391 ENDIF ! jpksedm1 400 392 401 393 ! ikwneg = jpksedm1 ; ikwneg+1 = jpksed ; ikwneg-1 = jpksed - 2 402 DO js = 1, jpsol 403 solcp(ji,ikwneg,js) =zfull(ikwneg) * zsolcpno(ikwneg ,js) &394 DO js = 1, jpsol + jpads 395 zsolcp(ji,ikwneg,js) = zfull(ikwneg) * zsolcpno(ikwneg ,js) & 404 396 & + zkwnup * zsolcpno(ikwneg-1,js) & 405 397 & + zkwnlo * zsolcpno(ikwneg+1,js) … … 407 399 fromsed(ji,js) = 0. 408 400 ENDDO 409 ! Heinze fromsed(ji,jsclay) = zempty * 1. * denssol * por1(jpksed) / mol_wgt(jsclay)410 401 fromsed(ji,jsclay) = zempty(jpksed) * 1. * por1clay 411 402 … … 413 404 ENDIF ! zwb > 0 414 405 ENDDO ! ji = 1, jpoce 406 407 DO js = 1, jpsol 408 solcp(:,:,js) = zsolcp(:,:,js) 409 END DO 410 DO jk = 2, jpksed 411 pwcp(:,jk,jwnh4) = (pwcp(:,jk,jwnh4) + zsolcp(:,jk,jpsol+1) * por1(jk) / por(jk) ) * radsnh4(jk) 412 pwcp(:,jk,jwfe2) = (pwcp(:,jk,jwfe2) + zsolcp(:,jk,jpsol+2) * por1(jk) / por(jk) ) * radsfe2(jk) 413 END DO 415 414 416 415 rainrm(:,:) = 0. -
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/sedarr.F90
r14786 r15075 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 -
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/sedbtb.F90
r10222 r15075 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,jpksedm1,jpsol+2) :: zsol, 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 51 47 52 48 53 ! Initializations 49 54 !---------------- 50 55 zsol(:,:,:) = 0. 56 zrearat = 0. 57 58 ! Remineralization rates of the different POC pools 59 zrearat(:,:,jspoc) = -reac_pocl 60 zrearat(:,:,jspos) = -reac_pocs 61 zrearat(:,:,jspor) = -reac_pocr 51 62 52 63 ! right hand side of coefficient matrix … … 54 65 DO js = 1, jpsol 55 66 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 67 zsol(:,jk,js) = solcp(:,jk+1,js) 68 END DO 69 END DO 61 70 62 CALL sed_mat( jpsol, jpoce, jpksedm1, zsol, dtsed / 2.0 ) 71 DO jk = 1, jpksedm1 72 zsol(:,jk,jpsol+1) = pwcp(:,jk+1,jwnh4) * adsnh4 73 zsol(:,jk,jpsol+2) = pwcp(:,jk+1,jwfe2) * adsfe2 74 END DO 63 75 76 CALL sed_mat_btb( jpksedm1, jpsol+2, zsol, zrearat(:,:,:), dtsed ) 64 77 65 78 ! store solution of the tridiagonal system … … 67 80 DO js = 1, jpsol 68 81 DO jk = 1, jpksedm1 69 DO ji = 1, jpoce 70 solcp(ji,jk+1,js) = zsol(ji,jk,js) 71 ENDDO 72 ENDDO 82 solcp(:,jk+1,js) = zsol(:,jk,js) 83 END DO 73 84 ENDDO 85 86 ! NH4 87 DO jk = 1, jpksedm1 88 pwcp(:,jk+1,jwnh4) = (pwcp(:,jk+1,jwnh4) + zsol(:,jk,jpsol+1) * por1(jk+1) / por(jk+1) ) * radsnh4(jk+1) 89 END DO 90 91 ! Fe2 92 DO jk = 1, jpksedm1 93 pwcp(:,jk+1,jwfe2) = (pwcp(:,jk+1,jwfe2) + zsol(:,jk,jpsol+2) * por1(jk+1) / por(jk+1) ) * radsfe2(jk+1) 94 END DO 95 96 DO ji = 1, jpoce 97 98 zsumtot = 0. 99 DO jk = 2, jpksed 100 zsolid1 = volc(ji,jk,jspoc) * solcp(ji,jk,jspoc) 101 zsolid2 = volc(ji,jk,jspos) * solcp(ji,jk,jspos) 102 zsolid3 = volc(ji,jk,jspor) * solcp(ji,jk,jspor) 103 rearatpom(ji,jk) = ( reac_pocl * zsolid1 + reac_pocs * zsolid2 + reac_pocr * zsolid3 ) * dtsed 104 zsumtot = zsumtot + rearatpom(ji,jk) / dtsed * volw3d(ji,jk) * 1.e-3 * 86400. * 365. * 1E3 105 END DO 106 107 ! 4/ Computation of the bioturbation coefficient 108 ! This parameterization is taken from Archer et al. (2002) 109 ! -------------------------------------------------------------- 110 zlimo2 = pwcp(ji,1,jwoxy) / (pwcp(ji,1,jwoxy) + 20.E-6) 111 db(ji,:) = dbiot * zsumtot**0.85 * zlimo2 / (365.0 * 86400.0) 112 113 ! ------------------------------------------------------ 114 ! Vertical variations of the bioturbation coefficient 115 ! ------------------------------------------------------ 116 IF (ln_btbz) THEN 117 DO jk = 1, jpksed 118 IF (profsedw(jk) < 2.0 * dbtbzsc) THEN 119 db(ji,jk) = db(ji,jk) * exp( -(profsedw(jk) / dbtbzsc)**2 ) 120 ELSE 121 db(ji,jk) = 0. 122 ENDIF 123 END DO 124 ELSE 125 DO jk = 1, jpksed 126 IF (profsedw(jk) > dbtbzsc) THEN 127 db(ji,jk) = 0.0 128 ENDIF 129 END DO 130 ENDIF 131 132 ! Computation of the bioirrigation factor (from Archer, MUDS model) 133 irrig(ji,:) = 0.0 134 IF (ln_irrig) THEN 135 DO jk = 1, jpksed 136 irrig(ji,jk) = ( 7.63752 - 7.4465 * exp( -0.89603 * zsumtot ) ) * zlimo2 137 irrig(ji,jk) = irrig(ji,jk) * exp( -(profsedw(jk) / xirrzsc) ) 138 END DO 139 ENDIF 140 END DO 141 142 ! CALL inorganic and organic slow redow/chemistry processes 143 ! --------------------------------------------------------- 144 CALL sed_inorg( kt ) 145 CALL sed_org( kt ) 74 146 75 147 IF( ln_timing ) CALL timing_stop('sed_btb') -
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/sedchem.F90
r14086 r15075 6 6 !!====================================================================== 7 7 !! modules used 8 USE par_sed9 8 USE sed ! sediment global variable 10 9 USE sedarr … … 24 23 REAL(wp), PARAMETER :: pp_rdel_ah_target = 1.E-4_wp 25 24 26 !! * Substitutions27 # include "do_loop_substitute.h90"28 25 !! * Module variables 29 26 REAL(wp) :: & … … 55 52 REAL(wp) :: devk19 = -26.57 56 53 REAL(wp) :: devk110 = -29.48 54 REAL(wp) :: devk120 = -14.8 55 REAL(wp) :: devk130 = -26.43 57 56 ! 58 57 REAL(wp) :: devk20 = 0.1271 … … 67 66 REAL(wp) :: devk29 = 0.2020 68 67 REAL(wp) :: devk210 = 0.1622 68 REAL(wp) :: devk220 = 0.0020 69 REAL(wp) :: devk230 = 0.0889 69 70 ! 70 71 REAL(wp) :: devk30 = 0. … … 79 80 REAL(wp) :: devk39 = -3.042e-3 80 81 REAL(wp) :: devk310 = -2.6080e-3 82 REAL(wp) :: devk320 = -0.4E-3 83 REAL(wp) :: devk330 = -0.905E-3 81 84 ! 82 85 REAL(wp) :: devk40 = -3.08E-3 … … 91 94 REAL(wp) :: devk49 = -4.08e-3 92 95 REAL(wp) :: devk410 = -2.84e-3 96 REAL(wp) :: devk420 = 2.89E-3 97 REAL(wp) :: devk430 = -5.03E-3 93 98 ! 94 99 REAL(wp) :: devk50 = 0.0877E-3 … … 103 108 REAL(wp) :: devk59 = 0.0714e-3 104 109 REAL(wp) :: devk510 = 0.0 110 REAL(wp) :: devk520 = 0.054E-3 111 REAL(wp) :: devk530 = 0.0814E-3 112 113 !! * Substitutions 114 # include "do_loop_substitute.h90" 105 115 106 116 !! $Id$ … … 180 190 ztc = temp(ji) 181 191 ztc2 = ztc * ztc 182 ! zqtt = ztkel * 0.01183 192 zsal = salt(ji) 184 193 zsal15 = SQRT( zsal ) * zsal … … 232 241 p_alkcb = pwcp(ji,jk,jwalk) / densSW(ji) 233 242 p_dictot = pwcp(ji,jk,jwdic) / densSW(ji) 234 p_bortot = borats(ji) / densSW(ji)243 p_bortot = borats(ji) 235 244 IF (p_alkcb <= 0.) THEN 236 245 p_hini(ji,jk) = 1.e-3 … … 283 292 DO jk = 1, jpksed 284 293 p_alknw_inf(:,jk) = -pwcp(:,jk,jwpo4) / densSW(:) 285 p_alknw_sup(:,jk) = (2. * pwcp(:,jk,jwdic) + 2. * pwcp(:,jk,jwpo4) + pwcp(:,jk,jwsil) &286 & + borats(:) ) / densSW(:)294 p_alknw_sup(:,jk) = (2. * pwcp(:,jk,jwdic) + 2. * pwcp(:,jk,jwpo4) + pwcp(:,jk,jwsil) ) & 295 & / densSW(:) + borats(:) 287 296 END DO 288 297 … … 313 322 REAL(wp) :: znumer_so4, zdnumer_so4, zdenom_so4, zalk_so4, zdalk_so4 314 323 REAL(wp) :: znumer_flu, zdnumer_flu, zdenom_flu, zalk_flu, zdalk_flu 324 REAL(wp) :: znumer_h2s, zdnumer_h2s, zdenom_h2s, zalk_h2s, zdalk_h2s 325 REAL(wp) :: znumer_nh3, zdnumer_nh3, zdenom_nh3, zalk_nh3, zdalk_nh3 315 326 REAL(wp) :: zalk_wat, zdalk_wat 316 REAL(wp) :: zfact, p_alktot, zdic, zbot, zpt, zst, zft, zsit 327 REAL(wp) :: zfact, p_alktot, zdic, zbot, zpt, zst, zft, zsit, zh2s, znh3 317 328 LOGICAL :: l_exitnow 318 329 REAL(wp), PARAMETER :: pz_exp_threshold = 1.0 … … 363 374 364 375 p_alktot = pwcp(ji,jk,jwalk) / densSW(ji) 365 zdic 366 zbot = borats(ji) / densSW(ji)367 zpt =pwcp(ji,jk,jwpo4) / densSW(ji)376 zdic = pwcp(ji,jk,jwdic) / densSW(ji) 377 zbot = borats(ji) 378 zpt = pwcp(ji,jk,jwpo4) / densSW(ji) 368 379 zsit = pwcp(ji,jk,jwsil) / densSW(ji) 369 zst = sulfats(ji) 370 zft = fluorids(ji) 380 zst = pwcp(ji,jk,jwso4) / densSW(ji) 381 zft = fluorids(ji) 382 zh2s = pwcp(ji,jk,jwh2s) / densSW(ji) 383 znh3 = pwcp(ji,jk,jwnh4) / densSW(ji) 371 384 aphscale = 1. + sulfats(ji)/aks3s(ji) 372 385 zh = zhi(ji,jk) … … 383 396 384 397 ! B(OH)3 - B(OH)4 : n=1, m=0 385 znumer_bor = akbs(ji)386 zdenom_bor = akbs(ji) + zh387 zalk_bor = zbot * (znumer_bor/zdenom_bor)398 znumer_bor = akbs(ji) 399 zdenom_bor = akbs(ji) + zh 400 zalk_bor = zbot * (znumer_bor/zdenom_bor) 388 401 zdnumer_bor = akbs(ji) 389 402 zdalk_bor = -zbot*(zdnumer_bor/zdenom_bor**2) … … 410 423 zdalk_sil = -zsit * (zdnumer_sil/zdenom_sil**2) 411 424 425 ! H2S - HS : n=1i, m=1 426 znumer_h2s = akh2s(ji) 427 zdenom_h2s = akh2s(ji) + zh 428 zalk_h2s = zh2s * (znumer_h2s/zdenom_h2s) 429 zdnumer_h2s = akh2s(ji) 430 zdalk_h2s = -zh2s * (zdnumer_h2s/zdenom_h2s**2) 431 432 ! NH4 - NH3 : n=1i, m=1 433 znumer_nh3 = aknh3(ji) 434 zdenom_nh3 = aknh3(ji) + zh 435 zalk_nh3 = znh3 * (znumer_nh3/zdenom_nh3) 436 zdnumer_nh3 = aknh3(ji) 437 zdalk_nh3 = -znh3 * (zdnumer_nh3/zdenom_nh3**2) 438 412 439 ! HSO4 - SO4 : n=1, m=1 413 440 aphscale = 1.0 + zst/aks3s(ji) … … 432 459 zeqn = zalk_dic + zalk_bor + zalk_po4 + zalk_sil & 433 460 & + zalk_so4 + zalk_flu & 434 & + zalk_wat - p_alktot461 & + zalk_wat + zalk_h2s + zalk_nh3 - p_alktot 435 462 436 463 zalka = p_alktot - (zalk_bor + zalk_po4 + zalk_sil & 437 & + zalk_so4 + zalk_flu + zalk_wat )464 & + zalk_so4 + zalk_flu + zalk_wat + zalk_h2s + zalk_nh3) 438 465 439 466 zdeqndh = zdalk_dic + zdalk_bor + zdalk_po4 + zdalk_sil & 440 & + zdalk_so4 + zdalk_flu + zdalk_wat 467 & + zdalk_so4 + zdalk_flu + zdalk_wat + zdalk_h2s + zdalk_nh3 441 468 442 469 ! Adapt bracketing interval … … 565 592 REAL(wp) :: zckb , zck1 , zck2 , zckw , zak1 , zak2 , zakb , zaksp0, zakw 566 593 REAL(wp) :: zck1p, zck2p, zck3p, zcksi, zak1p, zak2p, zak3p, zaksi 567 REAL(wp) :: zst , zft , zcks , zck f , zaksp1594 REAL(wp) :: zst , zft , zcks , zckhs, zckf , zaksp1, zcknh3 568 595 REAL(wp) :: total2free, free2SWS, total2SWS, SWS2total 569 596 !!--------------------------------------------------------------------- … … 630 657 & + LOG(1.0 - 0.001005 * zsal)) 631 658 659 ! DISSOCIATION CONSTANT FOR H2S on free H scale 660 zckhs = EXP(-13275.3 * ztr + 225.838 - 34.6435 * zlogt & 661 & + 0.3449 * zsqrt - 0.0274 * zsal ) 662 663 ! DISSOCIATION CONSTANT FOR NH3 on free H scale 664 zcknh3 = EXP(-6285.33 * ztr - 0.25444 + 0.0001635 * ztkel & 665 & + (0.46532 - 123.7184 * ztr) * zsqrt & 666 & + (-0.01992 + 3.17556 * ztr) * zsal ) 667 632 668 ! DISSOCIATION CONSTANT FOR FLUORIDES on free H scale (Dickson and Riley 79) 633 669 zckf = EXP( 1590.2*ztr - 12.641 + 1.525*zisqrt & … … 751 787 aksis(ji) = zaksi * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 752 788 789 zbuf1 = - ( devk120 + devk220 * ztc + devk320 * ztc * ztc ) 790 zbuf2 = 0.5 * ( devk420 + devk520 * ztc ) 791 akh2s(ji) = zckhs * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 792 793 zbuf1 = - ( devk130 + devk230 * ztc + devk330 * ztc * ztc ) 794 zbuf2 = 0.5 * ( devk430 + devk530 * ztc ) 795 aknh3(ji) = zcknh3 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 ) 796 753 797 ! Convert to total scale 754 798 ak1s(ji) = ak1s(ji) * SWS2total … … 761 805 aksis(ji) = aksis(ji) * SWS2total 762 806 akf3s(ji) = akf3s(ji) / total2free 807 akh2s(ji) = akh2s(ji) * SWS2total 808 aknh3(ji) = aknh3(ji) * SWS2total 763 809 764 810 ! APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE -
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/seddiff.F90
r14385 r15075 42 42 !! Arguments 43 43 INTEGER, INTENT(in) :: kt, knt ! number of iteration 44 ! --- local variables 45 INTEGER :: ji, jk, jw ! dummy looop indices 44 ! 46 45 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 solute diffusion in pore waters67 ! --------------------------------------68 69 ! Solves tridiagonal system70 DO jw = 1, jpwat71 CALL sed_mat( jw, jpoce, jpksed, zrearat1, zrearat2, pwcp(:,:,jw), dtsed2 / 2.0 )72 END DO73 74 CALL sed_mat( jwdic, jpoce, jpksed, zrearat1, zrearat2, sedligand(:,:), dtsed2 / 2.0 )75 76 IF( ln_timing ) CALL timing_stop('sed_diff')77 !78 46 END SUBROUTINE sed_diff 79 47 -
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/seddsr.F90
r14385 r15075 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( ji ) 30 27 !!---------------------------------------------------------------------- 31 28 !! *** ROUTINE sed_dsr *** … … 49 46 !!---------------------------------------------------------------------- 50 47 !! Arguments 51 INTEGER, INTENT(in) :: kt, knt! number of iteration48 INTEGER, INTENT(in) :: ji ! number of iteration 52 49 ! --- local variables 53 INTEGER :: ji, jk, js, jw, jn ! dummy looop indices 54 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) :: 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) :: zgamma, zbeta, zlimtmp 50 INTEGER :: jk, js, jw, jn ! dummy looop indices 51 52 REAL(wp), DIMENSION(jpksed) ::zlimo2, zlimno3, zlimso4, zlimfeo ! undersaturation ; indice jpwatp1 is for calcite 53 REAL(wp) :: zreasat 61 54 !! 62 55 !!---------------------------------------------------------------------- … … 64 57 IF( ln_timing ) CALL timing_start('sed_dsr') 65 58 ! 66 IF( kt == nitsed000 .AND. knt == 1 ) THEN67 IF (lwp) THEN68 WRITE(numsed,*) ' sed_dsr : Dissolution reaction '69 WRITE(numsed,*) ' '70 ENDIF71 ENDIF72 73 59 ! Initializations 74 60 !---------------------- 75 61 76 zrearat1(:,:) = 0. ; zundsat(:,:) = 0. 77 zlimo2 (:,:) = 0. ; zlimno3(:,:) = 0. ; zrearat2(:,:) = 0. 78 zlimso4(:,:) = 0. ; zrearat3(:,:) = 0. 79 zsumtot(:) = rtrn 62 zlimo2 (:) = 0. ; zlimno3(:) = 0. 63 zlimso4(:) = 0. 80 64 81 ALLOCATE( zvolc(jpoce, jpksed, jpsol) )82 zvolc(:,:,:) = 0.83 zadsnh4 = 1.0 / ( 1.0 + adsnh4 )84 85 ! Conversion of volume units86 !----------------------------87 DO js = 1, jpsol88 DO jk = 1, jpksed89 DO ji = 1, jpoce90 zvolc(ji,jk,js) = ( vols3d(ji,jk) * dens_mol_wgt(js) ) / &91 & ( volw3d(ji,jk) * 1.e-3 )92 ENDDO93 ENDDO94 ENDDO95 96 65 !---------------------------------------------------------- 97 66 ! 5. Beginning of solid reaction 98 67 !--------------------------------------------------------- 99 68 100 ! Definition of reaction rates [rearat]=sans dim 101 ! For jk=1 no reaction (pure water without solid) for each solid compo 102 zrearat1(:,:) = 0. 103 zrearat2(:,:) = 0. 104 zrearat3(:,:) = 0. 105 106 zundsat(:,:) = MAX( pwcp(:,:,jwoxy) - rtrn, 0. ) 107 108 DO jk = 2, jpksed 109 DO ji = 1, jpoce 110 zsolid1 = zvolc(ji,jk,jspoc) * solcp(ji,jk,jspoc) 111 zsolid2 = zvolc(ji,jk,jspos) * solcp(ji,jk,jspos) 112 zsolid3 = zvolc(ji,jk,jspor) * solcp(ji,jk,jspor) 113 zrearat1(ji,jk) = ( reac_pocl * dtsed2 * zsolid1 ) / ( 1. + reac_pocl * dtsed2 ) 114 zrearat2(ji,jk) = ( reac_pocs * dtsed2 * zsolid2 ) / ( 1. + reac_pocs * dtsed2 ) 115 zrearat3(ji,jk) = ( reac_pocr * dtsed2 * zsolid3 ) / ( 1. + reac_pocr * dtsed2 ) 116 solcp(ji,jk,jspoc) = solcp(ji,jk,jspoc) - zrearat1(ji,jk) / zvolc(ji,jk,jspoc) 117 solcp(ji,jk,jspos) = solcp(ji,jk,jspos) - zrearat2(ji,jk) / zvolc(ji,jk,jspos) 118 solcp(ji,jk,jspor) = solcp(ji,jk,jspor) - zrearat3(ji,jk) / zvolc(ji,jk,jspor) 119 zreasat = zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) 120 121 ! For DIC 122 pwcp(ji,jk,jwdic) = pwcp(ji,jk,jwdic) + zreasat 123 zsumtot(ji) = zsumtot(ji) + zreasat / dtsed2 * volw3d(ji,jk) * 1.e-3 * 86400. * 365. * 1E3 124 125 ! For alkalinity 126 pwcp(ji,jk,jwalk) = pwcp(ji,jk,jwalk) + zreasat * ( srno3 - 2.* spo4r ) 127 128 ! For Phosphate (in mol/l) 129 pwcp(ji,jk,jwpo4) = pwcp(ji,jk,jwpo4) + zreasat * spo4r 130 131 ! For iron (in mol/l) 132 pwcp(ji,jk,jwfe2) = pwcp(ji,jk,jwfe2) + fecratio(ji) * zreasat 133 134 ! Ammonium 135 pwcp(ji,jk,jwnh4) = pwcp(ji,jk,jwnh4) + zreasat * srno3 * zadsnh4 136 137 ! Ligands 138 sedligand(ji,jk) = sedligand(ji,jk) + ratligc * zreasat 139 140 ENDDO 141 ENDDO 142 143 ! left hand side of coefficient matrix 144 DO jk = 2, jpksed 145 DO ji = 1, jpoce 146 zreasat = zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) 147 zbeta = xksedo2 - pwcp(ji,jk,jwoxy) + so2ut * zreasat 148 zgamma = - xksedo2 * pwcp(ji,jk,jwoxy) 149 zundsat(ji,jk) = ( - zbeta + SQRT( zbeta**2 - 4.0 * zgamma ) ) / 2.0 150 zlimo2(ji,jk) = zundsat(ji,jk) / ( zundsat(ji,jk) + xksedo2 ) 151 152 ! Oxygen 153 pwcp(ji,jk,jwoxy) = zundsat(ji,jk) + MIN( pwcp(ji,jk,jwoxy), rtrn ) 154 zreasat = zreasat * zlimo2(ji,jk) 155 156 ! Ligands 157 sedligand(ji,jk) = sedligand(ji,jk) - reac_ligc * sedligand(ji,jk) 158 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 zreasat = rearatpom(ji,jk) 77 ! For alkalinity 78 pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) + zreasat * ( srno3 - 2.* spo4r ) 79 ! For Phosphate (in mol/l) 80 pwcpa(ji,jk,jwpo4) = pwcpa(ji,jk,jwpo4) + zreasat * spo4r 81 82 ! For iron (in mol/l) 83 pwcpa(ji,jk,jwfe2) = pwcpa(ji,jk,jwfe2) + fecratio(ji) * zreasat * radsfe2(jk) 84 ENDDO 85 86 ! Computes SMS of oxygen 87 DO jk = 2, jpksed 88 zlimo2(jk) = pwcp(ji,jk,jwoxy) / ( pwcp(ji,jk,jwoxy) + xksedo2 ) 89 zreasat = rearatpom(ji,jk) * zlimo2(jk) 90 ! Acid Silicic 91 pwcpa(ji,jk,jwoxy) = pwcpa(ji,jk,jwoxy) - so2ut * zreasat 92 ENDDO 93 94 !-------------------------------------------------------------------- 95 ! Denitrification 96 !-------------------------------------------------------------------- 97 DO jk = 2, jpksed 98 zlimno3(jk) = ( 1.0 - zlimo2(jk) ) * pwcp(ji,jk,jwno3) / ( pwcp(ji,jk,jwno3) + xksedno3 ) 99 zreasat = rearatpom(ji,jk) * zlimno3(jk) 100 ! For nitrates 101 pwcpa(ji,jk,jwno3) = pwcpa(ji,jk,jwno3) - zreasat * srDnit 102 ! For alkalinity 103 pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) + zreasat * srDnit 104 ENDDO 105 106 107 !-------------------------------------------------------------------- 108 ! Begining POC iron reduction 109 !-------------------------------------------------------------------- 110 111 DO jk = 2, jpksed 112 zlimfeo(jk) = ( 1.0 - zlimno3(jk) - zlimo2(jk) ) * solcp(ji,jk,jsfeo) / ( solcp(ji,jk,jsfeo) + xksedfeo ) 113 zreasat = rearatpom(ji,jk) * zlimfeo(jk) 114 ! For FEOH 115 solcpa(ji,jk,jsfeo) = solcpa(ji,jk,jsfeo) - 4.0 * zreasat / volc(ji,jk,jsfeo) 116 ! For PO4 117 pwcpa(ji,jk,jwpo4) = pwcpa(ji,jk,jwpo4) + zreasat * 4.0 * redfep 118 ! For alkalinity 119 pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) + 8.0 * zreasat 120 ! Iron 121 pwcpa(ji,jk,jwfe2) = pwcpa(ji,jk,jwfe2) + zreasat * 4.0 * radsfe2(jk) 159 122 ENDDO 160 123 161 124 !-------------------------------------------------------------------- 162 125 ! Begining POC denitrification and NO3- diffusion 163 ! (indice n°5 for couple POC/NO3- ie solcp(:,:,jspoc)/pwcp(:,:,jwno3)) 164 !-------------------------------------------------------------------- 165 166 zundsat(:,:) = MAX(0., pwcp(:,:,jwno3) - rtrn) 167 168 DO jk = 2, jpksed 169 DO ji = 1, jpoce 170 zlimtmp = ( 1.0 - zlimo2(ji,jk) ) 171 zreasat = zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) 172 zbeta = xksedno3 - pwcp(ji,jk,jwno3) + srDnit * zreasat * zlimtmp 173 zgamma = - xksedno3 * pwcp(ji,jk,jwno3) 174 zundsat(ji,jk) = ( - zbeta + SQRT( zbeta**2 - 4.0 * zgamma ) ) / 2.0 175 zlimno3(ji,jk) = zlimtmp * zundsat(ji,jk) / ( zundsat(ji,jk) + xksedno3 ) 176 177 ! For nitrates 178 pwcp(ji,jk,jwno3) = zundsat(ji,jk) + MIN(pwcp(ji,jk,jwno3), rtrn) 179 zreasat = zreasat * zlimno3(ji,jk) 180 181 ! For alkalinity 182 pwcp(ji,jk,jwalk) = pwcp(ji,jk,jwalk) + zreasat * srDnit 183 ENDDO 184 ENDDO 185 186 !-------------------------------------------------------------------- 187 ! Begining POC iron reduction 188 ! (indice n�5 for couple POFe(OH)3 ie solcp(:,:,jspoc)/pwcp(:,:,jsfeo)) 189 !-------------------------------------------------------------------- 190 191 zundsat(:,:) = MAX(0., solcp(:,:,jsfeo) - rtrn) 192 193 DO jk = 2, jpksed 194 DO ji = 1, jpoce 195 zlimtmp = 1.0 - zlimo2(ji,jk) - zlimno3(ji,jk) 196 zreasat = ( zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) ) / zvolc(ji,jk,jsfeo) 197 zbeta = xksedfeo - solcp(ji,jk,jsfeo) + 4.0 * zreasat * zlimtmp 198 zgamma = -xksedfeo * solcp(ji,jk,jsfeo) 199 zundsat(ji,jk) = ( - zbeta + SQRT( zbeta**2 - 4.0 * zgamma ) ) / 2.0 200 zlimfeo(ji,jk) = zlimtmp * zundsat(ji,jk) / ( zundsat(ji,jk) + xksedfeo ) 201 zreasat = ( zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) ) * zlimfeo(ji,jk) 202 203 ! For FEOH 204 solcp(ji,jk,jsfeo) = zundsat(ji,jk) + MIN(solcp(ji,jk,jsfeo), rtrn) 205 206 ! For Phosphate (in mol/l) 207 pwcp(ji,jk,jwpo4) = pwcp(ji,jk,jwpo4) + zreasat * 4.0 * redfep 208 209 ! For alkalinity 210 pwcp(ji,jk,jwalk) = pwcp(ji,jk,jwalk) + 8.0 * zreasat 211 212 ! Iron 213 pwcp(ji,jk,jwfe2) = pwcp(ji,jk,jwfe2) + zreasat * 4.0 214 ENDDO 215 ENDDO 216 217 ! !-------------------------------------------------------------------- 218 ! ! Begining POC denitrification and NO3- diffusion 219 ! ! (indice n�5 for couple POC/NO3- ie solcp(:,:,jspoc)/pwcp(:,:,jwno3)) 220 ! !-------------------------------------------------------------------- 221 ! 222 zundsat(:,:) = MAX(0., pwcp(:,:,jwso4) - rtrn) 223 224 DO jk = 2, jpksed 225 DO ji = 1, jpoce 226 zlimtmp = 1.0 - zlimo2(ji,jk) - zlimno3(ji,jk) - zlimfeo(ji,jk) 227 zreasat = zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) 228 zbeta = xksedso4 - pwcp(ji,jk,jwso4) + 0.5 * zreasat * zlimtmp 229 zgamma = - xksedso4 * pwcp(ji,jk,jwso4) 230 zundsat(ji,jk) = ( - zbeta + SQRT( zbeta**2 - 4.0 * zgamma ) ) / 2.0 231 zlimso4(ji,jk) = zlimtmp * zundsat(ji,jk) / ( zundsat(ji,jk) + xksedso4 ) 232 233 ! For sulfur 234 pwcp(ji,jk,jwso4) = zundsat(ji,jk) + MIN(pwcp(ji,jk,jwso4), rtrn) 235 zreasat = zreasat * zlimso4(ji,jk) 236 pwcp(ji,jk,jwh2s) = pwcp(ji,jk,jwh2s) + 0.5 * zreasat 237 238 ! For alkalinity 239 pwcp(ji,jk,jwalk) = pwcp(ji,jk,jwalk) + zreasat 240 ENDDO 241 ENDDO 242 243 ! Oxydation of the reduced products. Here only ammonium and ODU is accounted for 244 ! There are two options here: A simple time splitting scheme and a modified 245 ! Patankar scheme 246 ! ------------------------------------------------------------------------------ 247 call sed_dsr_redoxb 248 249 ! -------------------------------------------------------------- 250 ! 4/ Computation of the bioturbation coefficient 251 ! This parameterization is taken from Archer et al. (2002) 252 ! -------------------------------------------------------------- 253 254 DO ji = 1, jpoce 255 db(ji,:) = dbiot * zsumtot(ji) * pwcp(ji,1,jwoxy) / (pwcp(ji,1,jwoxy) + 20.E-6) 256 END DO 257 258 ! ------------------------------------------------------ 259 ! Vertical variations of the bioturbation coefficient 260 ! ------------------------------------------------------ 261 IF (ln_btbz) THEN 262 DO ji = 1, jpoce 263 db(ji,:) = db(ji,:) * exp( -(profsedw(:) / dbtbzsc)**2 ) / (365.0 * 86400.0) 264 END DO 265 ELSE 266 DO jk = 1, jpksed 267 IF (profsedw(jk) > dbtbzsc) THEN 268 db(:,jk) = 0.0 269 ENDIF 270 END DO 271 ENDIF 272 273 IF (ln_irrig) THEN 274 DO jk = 1, jpksed 275 DO ji = 1, jpoce 276 irrig(ji,jk) = ( 7.63752 - 7.4465 * exp( -0.89603 * zsumtot(ji) ) ) * pwcp(ji,1,jwoxy) & 277 & / (pwcp(ji,1,jwoxy) + 20.E-6) 278 irrig(ji,jk) = irrig(ji,jk) * exp( -(profsedw(jk) / xirrzsc) ) 279 END DO 280 END DO 281 ELSE 282 irrig(:,:) = 0.0 283 ENDIF 284 285 DEALLOCATE( zvolc ) 126 !-------------------------------------------------------------------- 127 128 DO jk = 2, jpksed 129 zlimso4(jk) = ( 1.0 - zlimno3(jk) - zlimo2(jk) - zlimfeo(jk) ) * pwcp(ji,jk,jwso4) / ( pwcp(ji,jk,jwso4) + xksedso4 ) 130 zreasat = rearatpom(ji,jk) * zlimso4(jk) 131 ! For sulfur 132 pwcpa(ji,jk,jwso4) = pwcpa(ji,jk,jwso4) - 0.5 * zreasat 133 pwcpa(ji,jk,jwh2s) = pwcpa(ji,jk,jwh2s) + 0.5 * zreasat 134 ! For alkalinity 135 pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) + zreasat 136 ENDDO 137 138 ! Secondary redox reactions 139 ! ------------------------- 140 141 call sed_dsr_redoxb( ji ) 286 142 287 143 IF( ln_timing ) CALL timing_stop('sed_dsr') … … 289 145 END SUBROUTINE sed_dsr 290 146 291 SUBROUTINE sed_dsr_redoxb 147 SUBROUTINE sed_dsr_redoxb(ji) 292 148 !!---------------------------------------------------------------------- 293 149 !! *** ROUTINE sed_dsr_redox *** … … 295 151 !! ** Purpose : computes secondary redox reactions 296 152 !! 297 !! ** Methode : It uses Strand splitter algorithm proposed by298 !! Nguyen et al. (2009) and modified by Wang et al. (2018)299 !! Basically, each equation is solved analytically when300 !! feasible, otherwise numerically at t+1/2. Then301 !! the last equation is solved at t+1. The other equations302 !! are then solved at t+1 starting in the reverse order.303 !! Ideally, it's better to start from the fastest reaction304 !! to the slowest and then reverse the order to finish up305 !! with the fastest one. But random order works well also.306 !! The scheme is second order, positive and mass307 !! conserving. It works well for stiff systems.308 !!309 153 !! History : 310 154 !! ! 18-08 (O. Aumont) Original code … … 312 156 !! Arguments 313 157 ! --- local variables 314 INTEGER :: ji, jk, jn, jw ! dummy looop indices315 316 REAL, DIMENSION(6) :: zsedtrn, zsedtra 317 REAL(wp) :: zalpha, z beta, zgamma, zdelta, zepsi, zsedfer158 INTEGER, INTENT(IN) :: ji 159 INTEGER :: jni, jnj, jk ! dummy looop indices 160 161 REAL(wp) :: zalpha, zexcess, zh2seq, zsedfer, zreasat 318 162 !! 319 163 !!---------------------------------------------------------------------- … … 321 165 IF( ln_timing ) CALL timing_start('sed_dsr_redoxb') 322 166 323 DO ji = 1, jpoce 324 DO jk = 2, jpksed 325 zsedtrn(1) = pwcp(ji,jk,jwoxy) 326 zsedtrn(2) = MAX(0., pwcp(ji,jk,jwh2s) ) 327 zsedtrn(3) = pwcp(ji,jk,jwnh4) 328 zsedtrn(4) = MAX(0., pwcp(ji,jk,jwfe2) - sedligand(ji,jk) ) 329 zsedfer = MIN(0., pwcp(ji,jk,jwfe2) - sedligand(ji,jk) ) 330 zsedtrn(5) = solcp(ji,jk,jsfeo) * zvolc(ji,jk,jsfeo) 331 zsedtrn(6) = solcp(ji,jk,jsfes) * zvolc(ji,jk,jsfes) 332 zsedtra(:) = MAX(0., zsedtrn(:) - rtrn ) 333 334 ! First pass of the scheme. At the end, it is 1st order 335 ! ----------------------------------------------------- 336 ! Fe + O2 337 zalpha = zsedtra(1) - 0.25 * zsedtra(4) 338 zbeta = zsedtra(4) + zsedtra(5) 339 zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4) 340 zdelta = pwcp(ji,jk,jwpo4) - redfep * zsedtra(4) 341 IF ( zalpha == 0. ) THEN 342 zsedtra(4) = zsedtra(4) / ( 1.0 + 0.25 * zsedtra(4) * reac_fe2 * dtsed2 / 2.0 ) 343 ELSE 344 zsedtra(4) = ( zsedtra(4) * zalpha ) / ( 0.25 * zsedtra(4) & 345 & * ( exp( reac_fe2 * zalpha * dtsed2 / 2. ) - 1.0 ) & 346 & + zalpha * exp( reac_fe2 * zalpha * dtsed2 / 2. ) ) 347 ENDIF 348 zsedtra(1) = zalpha + 0.25 * zsedtra(4) 349 zsedtra(5) = zbeta - zsedtra(4) 350 pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(4) 351 pwcp(ji,jk,jwpo4) = zdelta + redfep * zsedtra(4) 352 ! H2S + O2 353 zalpha = zsedtra(1) - 2.0 * zsedtra(2) 354 zbeta = pwcp(ji,jk,jwso4) + zsedtra(2) 355 zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(2) 356 IF ( zalpha == 0. ) THEN 357 zsedtra(2) = zsedtra(2) / ( 1.0 + 2.0 * zsedtra(2) * reac_h2s * dtsed2 / 2.0 ) 358 ELSE 359 zsedtra(2) = ( zsedtra(2) * zalpha ) / ( 2.0 * zsedtra(2) * ( exp( reac_h2s * zalpha * dtsed2 / 2. ) - 1.0 ) & 360 & + zalpha * exp( reac_h2s * zalpha * dtsed2 / 2. ) ) 361 ENDIF 362 zsedtra(1) = zalpha + 2.0 * zsedtra(2) 363 pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(2) 364 pwcp(ji,jk,jwso4) = zbeta - zsedtra(2) 365 ! NH4 + O2 366 zalpha = zsedtra(1) - 2.0 * zsedtra(3) / zadsnh4 367 zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(3) /zadsnh4 368 IF ( zalpha == 0. ) THEN 369 zsedtra(3) = zsedtra(3) / ( 1.0 + 2.0 * zsedtra(3) * reac_nh4 / zadsnh4 * dtsed2 / 2.0 ) 370 ELSE 371 zsedtra(3) = ( zsedtra(3) * zalpha * zadsnh4 ) / ( 2.0 * zsedtra(3) * ( exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 / 2. ) - 1.0 ) & 372 & + zalpha * zadsnh4 * exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 /2. ) ) 373 ENDIF 374 zsedtra(1) = zalpha + 2.0 * zsedtra(3) / zadsnh4 375 pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(3) / zadsnh4 376 ! FeS - O2 377 zalpha = zsedtra(1) - 2.0 * zsedtra(6) 378 zbeta = zsedtra(4) + zsedtra(6) 379 zgamma = pwcp(ji,jk,jwso4) + zsedtra(6) 380 IF ( zalpha == 0. ) THEN 381 zsedtra(6) = zsedtra(6) / ( 1.0 + 2.0 * zsedtra(6) * reac_feso * dtsed2 / 2.0 ) 382 ELSE 383 zsedtra(6) = ( zsedtra(6) * zalpha ) / ( 2.0 * zsedtra(6) * ( exp( reac_feso * zalpha * dtsed2 / 2. ) - 1.0 ) & 384 & + zalpha * exp( reac_feso * zalpha * dtsed2 /2. ) ) 385 ENDIF 386 zsedtra(1) = zalpha + 2.0 * zsedtra(6) 387 zsedtra(4) = zbeta - zsedtra(6) 388 pwcp(ji,jk,jwso4) = zgamma - zsedtra(6) 389 ! ! Fe - H2S 390 zalpha = zsedtra(2) - zsedtra(4) 391 zbeta = zsedtra(4) + zsedtra(6) 392 zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4) 393 IF ( zalpha == 0. ) THEN 394 zsedtra(4) = zsedtra(4) / ( 1.0 + zsedtra(4) * reac_fes * dtsed2 / 2.0 ) 395 ELSE 396 zsedtra(4) = ( zsedtra(4) * zalpha ) / ( zsedtra(4) * ( exp( reac_fes * zalpha * dtsed2 / 2. ) - 1.0 ) & 397 & + zalpha * exp( reac_fes * zalpha * dtsed2 /2. ) ) 398 ENDIF 399 zsedtra(2) = zalpha + zsedtra(4) 400 zsedtra(6) = zbeta - zsedtra(4) 401 pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(4) 402 ! FEOH + H2S 403 zalpha = zsedtra(5) - 2.0 * zsedtra(2) 404 zbeta = zsedtra(5) + zsedtra(4) 405 zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4) 406 zdelta = pwcp(ji,jk,jwso4) + zsedtra(2) 407 zepsi = pwcp(ji,jk,jwpo4) + redfep * zsedtra(5) 408 IF ( zalpha == 0. ) THEN 409 zsedtra(2) = zsedtra(2) / ( 1.0 + 2.0 * zsedtra(2) * reac_feh2s * dtsed2 ) 410 ELSE 411 zsedtra(2) = ( zsedtra(2) * zalpha ) / ( 2.0 * zsedtra(2) * ( exp( reac_feh2s * zalpha * dtsed2 ) - 1.0 ) & 412 & + zalpha * exp( reac_feh2s * zalpha * dtsed2 ) ) 413 ENDIF 414 zsedtra(5) = zalpha + 2.0 * zsedtra(2) 415 zsedtra(4) = zbeta - zsedtra(5) 416 pwcp(ji,jk,jwso4) = zdelta - zsedtra(2) 417 pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(4) 418 pwcp(ji,jk,jwpo4) = zepsi - redfep * zsedtra(5) 419 ! Fe - H2S 420 zalpha = zsedtra(2) - zsedtra(4) 421 zbeta = zsedtra(4) + zsedtra(6) 422 zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4) 423 IF ( zalpha == 0. ) THEN 424 zsedtra(4) = zsedtra(4) / ( 1.0 + zsedtra(4) * reac_fes * dtsed2 / 2.0 ) 425 ELSE 426 zsedtra(4) = ( zsedtra(4) * zalpha ) / ( zsedtra(4) * ( exp( reac_fes * zalpha * dtsed2 / 2. ) - 1.0 ) & 427 & + zalpha * exp( reac_fes * zalpha * dtsed2 /2. ) ) 428 ENDIF 429 zsedtra(2) = zalpha + zsedtra(4) 430 zsedtra(6) = zbeta - zsedtra(4) 431 pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(4) 432 ! FeS - O2 433 zalpha = zsedtra(1) - 2.0 * zsedtra(6) 434 zbeta = zsedtra(4) + zsedtra(6) 435 zgamma = pwcp(ji,jk,jwso4) + zsedtra(6) 436 IF (zalpha == 0.) THEN 437 zsedtra(6) = zsedtra(6) / ( 1.0 + 2.0 * zsedtra(6) * reac_feso * dtsed2 / 2. ) 438 ELSE 439 zsedtra(6) = ( zsedtra(6) * zalpha ) / ( 2.0 * zsedtra(6) * ( exp( reac_feso * zalpha * dtsed2 / 2. ) - 1.0 ) & 440 & + zalpha * exp( reac_feso * zalpha * dtsed2 /2. ) ) 441 ENDIF 442 zsedtra(1) = zalpha + 2.0 * zsedtra(6) 443 zsedtra(4) = zbeta - zsedtra(6) 444 pwcp(ji,jk,jwso4) = zgamma - zsedtra(6) 445 446 ! NH4 + O2 447 zalpha = zsedtra(1) - 2.0 * zsedtra(3) / zadsnh4 448 zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(3) / zadsnh4 449 IF ( zalpha == 0. ) THEN 450 zsedtra(3) = zsedtra(3) / ( 1.0 + 2.0 * zsedtra(3) * reac_nh4 / zadsnh4 * dtsed2 / 2.0 ) 451 ELSE 452 zsedtra(3) = ( zsedtra(3) * zalpha * zadsnh4 ) / ( 2.0 * zsedtra(3) * ( exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 / 2. ) - 1.0 ) & 453 & + zalpha * zadsnh4 * exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 /2. ) ) 454 ENDIF 455 zsedtra(1) = zalpha + 2.0 * zsedtra(3) / zadsnh4 456 pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(3) / zadsnh4 457 ! H2S + O2 458 zalpha = zsedtra(1) - 2.0 * zsedtra(2) 459 zbeta = pwcp(ji,jk,jwso4) + zsedtra(2) 460 zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(2) 461 IF ( zalpha == 0. ) THEN 462 zsedtra(2) = zsedtra(2) / ( 1.0 + 2.0 * zsedtra(2) * reac_h2s * dtsed2 / 2.0 ) 463 ELSE 464 zsedtra(2) = ( zsedtra(2) * zalpha ) / ( 2.0 * zsedtra(2) * ( exp( reac_h2s * zalpha * dtsed2 / 2. ) - 1.0 ) & 465 & + zalpha * exp( reac_h2s * zalpha * dtsed2 / 2. ) ) 466 ENDIF 467 zsedtra(1) = zalpha + 2.0 * zsedtra(2) 468 pwcp(ji,jk,jwso4) = zbeta - zsedtra(2) 469 pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(2) 470 471 ! Fe + O2 472 zalpha = zsedtra(1) - 0.25 * zsedtra(4) 473 zbeta = zsedtra(4) + zsedtra(5) 474 zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4) 475 zdelta = pwcp(ji,jk,jwpo4) - redfep * zsedtra(4) 476 IF ( zalpha == 0. ) THEN 477 zsedtra(4) = zsedtra(4) / ( 1.0 + 0.25 * zsedtra(4) * reac_fe2 * dtsed2 / 2.0 ) 478 ELSE 479 zsedtra(4) = ( zsedtra(4) * zalpha ) / ( 0.25 * zsedtra(4) & 480 & * ( exp( reac_fe2 * zalpha * dtsed2 / 2. ) - 1.0 ) & 481 & + zalpha * exp( reac_fe2 * zalpha * dtsed2 / 2. ) ) 482 ENDIF 483 zsedtra(1) = zalpha + 0.25 * zsedtra(4) 484 zsedtra(5) = zbeta - zsedtra(4) 485 pwcp(ji,jk,jwpo4) = zdelta + redfep * zsedtra(4) 486 pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(4) 487 488 ! Update the concentrations after the secondary reactions 489 zsedtra(:) = zsedtra(:) + MIN( zsedtrn(:), rtrn ) 490 pwcp(ji,jk,jwoxy) = zsedtra(1) 491 pwcp(ji,jk,jwh2s) = zsedtra(2) 492 pwcp(ji,jk,jwnh4) = zsedtra(3) 493 pwcp(ji,jk,jwfe2) = zsedtra(4) + sedligand(ji,jk) + zsedfer 494 pwcp(ji,jk,jwno3) = pwcp(ji,jk,jwno3) + ( zsedtrn(3) - pwcp(ji,jk,jwnh4) ) 495 solcp(ji,jk,jsfeo) = zsedtra(5) / zvolc(ji,jk,jsfeo) 496 solcp(ji,jk,jsfes) = zsedtra(6) / zvolc(ji,jk,jsfes) 497 END DO 498 END DO 167 DO jk = 2, jpksed 168 zalpha = ( pwcp(ji,jk,jwfe2) - pwcp(ji,jk,jwlgw) ) * 1E9 169 zsedfer = ( zalpha + SQRT( zalpha**2 + 1.E-10 ) ) /2.0 * 1E-9 170 171 ! First pass of the scheme. At the end, it is 1st order 172 ! ----------------------------------------------------- 173 ! Fe (both adsorbed and solute) + O2 174 zreasat = reac_fe2 * dtsed * zsedfer / radsfe2(jk) * pwcp(ji,jk,jwoxy) 175 pwcpa(ji,jk,jwfe2) = pwcpa(ji,jk,jwfe2) - zreasat * radsfe2(jk) 176 pwcpa(ji,jk,jwoxy) = pwcpa(ji,jk,jwoxy) - 0.25 * zreasat 177 pwcpa(ji,jk,jwpo4) = pwcpa(ji,jk,jwpo4) - redfep * zreasat 178 pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) - 2.0 * zreasat 179 solcpa(ji,jk,jsfeo) = solcpa(ji,jk,jsfeo) + zreasat / volc(ji,jk,jsfeo) 180 181 ! H2S + O2 182 zreasat = reac_h2s * dtsed * pwcp(ji,jk,jwoxy) * pwcp(ji,jk,jwh2s) 183 pwcpa(ji,jk,jwh2s) = pwcpa(ji,jk,jwh2s) - zreasat 184 pwcpa(ji,jk,jwoxy) = pwcpa(ji,jk,jwoxy) - 2.0 * zreasat 185 pwcpa(ji,jk,jwso4) = pwcpa(ji,jk,jwso4) + zreasat 186 pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) - 2.0 * zreasat 187 188 ! NH4 + O2 189 zreasat = reac_nh4 * dtsed * pwcp(ji,jk,jwnh4) / radsnh4(jk) * pwcp(ji,jk,jwoxy) 190 pwcpa(ji,jk,jwnh4) = pwcpa(ji,jk,jwnh4) - zreasat * radsnh4(jk) 191 pwcpa(ji,jk,jwoxy) = pwcpa(ji,jk,jwoxy) - 2.0 * zreasat 192 pwcpa(ji,jk,jwno3) = pwcpa(ji,jk,jwno3) + zreasat 193 pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) - 2.0 * zreasat 194 195 ! FeS - O2 196 zreasat = reac_feso * dtsed * solcp(ji,jk,jsfes) * volc(ji,jk,jsfes) * pwcp(ji,jk,jwoxy) 197 solcpa(ji,jk,jsfes) = solcpa(ji,jk,jsfes) - zreasat / volc(ji,jk,jsfes) 198 pwcpa(ji,jk,jwoxy) = pwcpa(ji,jk,jwoxy) - 2.0 * zreasat 199 pwcpa(ji,jk,jwfe2) = pwcpa(ji,jk,jwfe2) + zreasat * radsfe2(jk) 200 pwcpa(ji,jk,jwso4) = pwcpa(ji,jk,jwso4) + zreasat 201 202 ! FEOH + H2S 203 zreasat = reac_feh2s * dtsed * solcp(ji,jk,jsfeo) * volc(ji,jk,jsfeo) * pwcp(ji,jk,jwh2s) 204 solcpa(ji,jk,jsfeo) = solcpa(ji,jk,jsfeo) - 8.0 * zreasat / volc(ji,jk,jsfeo) 205 pwcpa(ji,jk,jwh2s) = pwcpa(ji,jk,jwh2s) - zreasat 206 pwcpa(ji,jk,jwfe2) = pwcpa(ji,jk,jwfe2) + 8.0 * zreasat * radsfe2(jk) 207 pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) + 14.0 * zreasat 208 pwcpa(ji,jk,jwso4) = pwcpa(ji,jk,jwso4) + zreasat 209 pwcpa(ji,jk,jwpo4) = pwcpa(ji,jk,jwpo4) + 8.0 * redfep * zreasat 210 211 ! Fe + H2S 212 zh2seq = MAX(rtrn, 2.1E-3 * hipor(ji,jk)) 213 zexcess = pwcp(ji,jk,jwh2s) * zsedfer / zh2seq - 1.0 214 IF ( zexcess >= 0.0 ) THEN 215 zreasat = reac_fesp * zexcess * dtsed 216 pwcpa (ji,jk,jwfe2) = pwcpa(ji,jk,jwfe2) - zreasat * radsfe2(jk) 217 solcpa(ji,jk,jsfes) = solcpa(ji,jk,jsfes) + zreasat / volc(ji,jk,jsfes) 218 pwcpa(ji,jk,jwh2s) = pwcpa(ji,jk,jwh2s) - zreasat 219 ELSE 220 zreasat = reac_fesd * zexcess * dtsed * solcp(ji,jk,jsfes) * volc(ji,jk,jsfes) 221 pwcpa (ji,jk,jwfe2) = pwcpa(ji,jk,jwfe2) - zreasat * radsfe2(jk) 222 solcpa(ji,jk,jsfes) = solcpa(ji,jk,jsfes) + zreasat / volc(ji,jk,jsfes) 223 pwcpa(ji,jk,jwh2s) = pwcpa(ji,jk,jwh2s) - zreasat 224 ENDIF 225 ! For alkalinity 226 pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) - zreasat * 2.0 227 END DO 499 228 500 229 IF( ln_timing ) CALL timing_stop('sed_dsr_redoxb') -
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/seddta.F90
r14385 r15075 25 25 # include "do_loop_substitute.h90" 26 26 # include "domzgr_substitute.h90" 27 !! $Id$ 27 28 28 CONTAINS 29 29 … … 46 46 47 47 !! Arguments 48 INTEGER, INTENT( in) :: kt! time-step49 INTEGER, INTENT( in) :: Kbb, Kmm! time level indices48 INTEGER, INTENT( in ) :: kt ! time-step 49 INTEGER, INTENT( in ) :: Kbb, Kmm ! time level indices 50 50 51 51 !! * Local declarations … … 55 55 REAL(wp), DIMENSION(jpi,jpj) :: zwsbio4, zwsbio3, zddust 56 56 REAL(wp) :: zf0, zf1, zf2, zkapp, zratio, zdep 57 REAL(wp) :: zzf0, zf0s, zf0b, zzf1, zf1s, zf1b, zzf2, zf2s, zf2b 57 58 58 59 !---------------------------------------------------------------------- … … 78 79 dtsed = rDt_trc 79 80 rsecday = 60.* 60. * 24. 80 ! conv2 = 1.0e+3 / ( 1.0e+4 * rsecday * 30. )81 81 conv2 = 1.0e+3 / 1.0e+4 82 82 ENDIF … … 114 114 ikt = mbkt(ji,jj) 115 115 IF ( tmask(ji,jj,ikt) == 1 ) THEN 116 trc_data(ji,jj,1) = tr(ji,jj,ikt,jpsil,Kbb) 117 trc_data(ji,jj,2) = tr(ji,jj,ikt,jpoxy,Kbb) 118 trc_data(ji,jj,3) = tr(ji,jj,ikt,jpdic,Kbb) 119 trc_data(ji,jj,4) = tr(ji,jj,ikt,jpno3,Kbb) / 7.625 120 trc_data(ji,jj,5) = tr(ji,jj,ikt,jppo4,Kbb) / 122. 121 trc_data(ji,jj,6) = tr(ji,jj,ikt,jptal,Kbb) 122 trc_data(ji,jj,7) = tr(ji,jj,ikt,jpnh4,Kbb) / 7.625 123 trc_data(ji,jj,8) = 0.0 124 trc_data(ji,jj,9) = 28.0E-3 125 trc_data(ji,jj,10) = tr(ji,jj,ikt,jpfer,Kbb) 126 trc_data(ji,jj,11 ) = MIN(tr(ji,jj,ikt,jpgsi,Kbb), 1E-4) * zwsbio4(ji,jj) * 1E3 127 trc_data(ji,jj,12 ) = MIN(tr(ji,jj,ikt,jppoc,Kbb), 1E-4) * zwsbio3(ji,jj) * 1E3 128 trc_data(ji,jj,13 ) = MIN(tr(ji,jj,ikt,jpgoc,Kbb), 1E-4) * zwsbio4(ji,jj) * 1E3 129 trc_data(ji,jj,14) = MIN(tr(ji,jj,ikt,jpcal,Kbb), 1E-4) * zwsbio4(ji,jj) * 1E3 130 trc_data(ji,jj,15) = ts(ji,jj,ikt,jp_tem,Kmm) 131 trc_data(ji,jj,16) = ts(ji,jj,ikt,jp_sal,Kmm) 132 trc_data(ji,jj,17 ) = ( tr(ji,jj,ikt,jpsfe,Kbb) * zwsbio3(ji,jj) + tr(ji,jj,ikt,jpbfe,Kbb) & 133 & * zwsbio4(ji,jj) ) * 1E3 / ( trc_data(ji,jj,12 ) + trc_data(ji,jj,13 ) + rtrn ) 134 trc_data(ji,jj,17 ) = MIN(1E-3, trc_data(ji,jj,17 ) ) 116 trc_data(ji,jj,jwsil) = tr(ji,jj,ikt,jpsil,Kbb) 117 trc_data(ji,jj,jwoxy) = tr(ji,jj,ikt,jpoxy,Kbb) 118 trc_data(ji,jj,jwdic) = tr(ji,jj,ikt,jpdic,Kbb) 119 trc_data(ji,jj,jwno3) = tr(ji,jj,ikt,jpno3,Kbb) / 7.625 120 trc_data(ji,jj,jwpo4) = tr(ji,jj,ikt,jppo4,Kbb) / 122. 121 trc_data(ji,jj,jwalk) = tr(ji,jj,ikt,jptal,Kbb) 122 trc_data(ji,jj,jwnh4) = tr(ji,jj,ikt,jpnh4,Kbb) / 7.625 123 trc_data(ji,jj,jwh2s) = 0.0 124 trc_data(ji,jj,jwso4) = 0.14 * ts(ji,jj,ikt,jp_sal,Kmm) / 1.80655 / 96.062 125 trc_data(ji,jj,jwfe2) = tr(ji,jj,ikt,jpfer,Kbb) 126 trc_data(ji,jj,jwlgw) = 1E-9 127 trc_data(ji,jj,12 ) = MIN(tr(ji,jj,ikt,jpgsi,Kbb), 1E-4) * zwsbio4(ji,jj) * 1E3 128 trc_data(ji,jj,13 ) = MIN(tr(ji,jj,ikt,jppoc,Kbb), 1E-4) * zwsbio3(ji,jj) * 1E3 129 trc_data(ji,jj,14 ) = MIN(tr(ji,jj,ikt,jpgoc,Kbb), 1E-4) * zwsbio4(ji,jj) * 1E3 130 trc_data(ji,jj,15) = MIN(tr(ji,jj,ikt,jpcal,Kbb), 1E-4) * zwsbio4(ji,jj) * 1E3 131 trc_data(ji,jj,16) = ts(ji,jj,ikt,jp_tem,Kmm) 132 trc_data(ji,jj,17) = ts(ji,jj,ikt,jp_sal,Kmm) 133 trc_data(ji,jj,18 ) = ( tr(ji,jj,ikt,jpsfe,Kbb) * zwsbio3(ji,jj) + tr(ji,jj,ikt,jpbfe,Kbb) & 134 & * zwsbio4(ji,jj) ) * 1E3 / ( trc_data(ji,jj,13 ) + trc_data(ji,jj,14 ) + rtrn ) 135 trc_data(ji,jj,18 ) = MIN(1E-3, trc_data(ji,jj,18 ) ) 135 136 ENDIF 136 137 END_2D … … 141 142 CALL pack_arr ( jpoce, pwcp_dta(1:jpoce,jw), trc_data(1:jpi,1:jpj,jw), iarroce(1:jpoce) ) 142 143 END DO 144 143 145 ! Solid components : 144 146 !----------------------- 145 147 ! Sinking fluxes for OPAL in mol.m-2.s-1 ; conversion in mol.cm-2.s-1 146 CALL pack_arr ( jpoce, rainrm_dta(1:jpoce,jsopal), trc_data(1:jpi,1:jpj,1 1), iarroce(1:jpoce) )148 CALL pack_arr ( jpoce, rainrm_dta(1:jpoce,jsopal), trc_data(1:jpi,1:jpj,12), iarroce(1:jpoce) ) 147 149 rainrm_dta(1:jpoce,jsopal) = rainrm_dta(1:jpoce,jsopal) * 1e-4 148 150 ! Sinking fluxes for POC in mol.m-2.s-1 ; conversion in mol.cm-2.s-1 149 CALL pack_arr ( jpoce, zdtap(1:jpoce), trc_data(1:jpi,1:jpj,1 2) , iarroce(1:jpoce) )150 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) ) 151 153 DO ji = 1, jpoce 152 ! zkapp = MIN( (1.0 - 0.02 ) * reac_poc, 3731.0 * max(100.0, zkbot(ji) )**(-1.011) / ( 365.0 * 24.0 * 3600.0 ) ) 153 ! zkapp = MIN( 0.98 * reac_poc, 100.0 * max(100.0, zkbot(ji) )**(-0.6) / ( 365.0 * 24.0 * 3600.0 ) ) 154 ! 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 ) 155 ! zf1 = ( 0.02 * (reac_poc - reac_poc * 0.) + zkapp - reac_poc ) / ( reac_poc / 100. - reac_poc ) 156 ! zf1 = MIN(0.98, MAX(0., zf1 ) ) 157 zf1 = 0.48 158 zf2 = 0.02 159 zf0 = 1.0 - zf2 - zf1 160 rainrm_dta(ji,jspoc) = ( zdtap(ji) + zdtag(ji) ) * 1e-4 * zf0 161 rainrm_dta(ji,jspos) = ( zdtap(ji) + zdtag(ji) ) * 1e-4 * zf1 162 rainrm_dta(ji,jspor) = ( zdtap(ji) + zdtag(ji) ) * 1e-4 * zf2 154 zzf2 = 2E-2 155 zzf1 = 0.5 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 163 166 END DO 164 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 salinity169 CALL pack_arr ( jpoce, temp(1:jpoce), trc_data(1:jpi,1:jpj,1 5), iarroce(1:jpoce) )170 CALL pack_arr ( jpoce, salt(1:jpoce), trc_data(1:jpi,1:jpj,1 6), iarroce(1:jpoce) )171 ! vector temperature [°C] and salinity 172 CALL pack_arr ( jpoce, temp(1:jpoce), trc_data(1:jpi,1:jpj,16), iarroce(1:jpoce) ) 173 CALL pack_arr ( jpoce, salt(1:jpoce), trc_data(1:jpi,1:jpj,17), iarroce(1:jpoce) ) 171 174 172 175 ! Clay rain rate in [mol/(cm**2.s)] … … 187 190 188 191 ! Fe/C ratio in sinking particles that fall to the sediments 189 CALL pack_arr ( jpoce, fecratio(1:jpoce), trc_data(1:jpi,1:jpj,17), iarroce(1:jpoce) ) 190 191 sedligand(:,1) = 1.E-9 192 CALL pack_arr ( jpoce, fecratio(1:jpoce), trc_data(1:jpi,1:jpj,18), iarroce(1:jpoce) ) 192 193 193 194 ! sediment pore water at 1st layer (k=1) … … 217 218 IF( lk_iomput ) THEN 218 219 IF( iom_use("sflxclay" ) ) CALL iom_put( "sflxclay", zddust(:,:) * 1E3 / 1.E4 ) 219 IF( iom_use("sflxcal" ) ) CALL iom_put( "sflxcal", trc_data(:,:,1 4) / 1.E4 )220 IF( iom_use("sflxbsi" ) ) CALL iom_put( "sflxbsi", trc_data(:,:,1 1) / 1.E4 )221 IF( iom_use("sflxpoc" ) ) CALL iom_put( "sflxpoc", ( trc_data(:,:,1 2) + trc_data(:,:,13) ) / 1.E4 )220 IF( iom_use("sflxcal" ) ) CALL iom_put( "sflxcal", trc_data(:,:,15) / 1.E4 ) 221 IF( iom_use("sflxbsi" ) ) CALL iom_put( "sflxbsi", trc_data(:,:,12) / 1.E4 ) 222 IF( iom_use("sflxpoc" ) ) CALL iom_put( "sflxpoc", ( trc_data(:,:,13) + trc_data(:,:,14) ) / 1.E4 ) 222 223 ENDIF 223 224 -
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/sedini.F90
r14385 r15075 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 … … 55 53 rcfe2 = 5.E8 , & !: reactivity for O2/Fe2+ [l.mol-1.an-1] 56 54 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 55 rcfeso = 3.E5 , & !: Reactivity for FES/O2 [l.mol-1.an-1] 56 rcfesp = 5E-6 , & !: Precipitation of FeS [mol/l-1.an-1] 57 rcfesd = 1E-3 , & !: Dissolution of FeS [an-1] 59 58 xksedo2 = 5E-6 , & !: half-sturation constant for oxic remin. 60 59 xksedno3 = 5E-6 , & !: half-saturation constant for denitrification … … 73 72 74 73 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.81E-6, 9.8E-6, 9.73E-6, 5.0E-6, 3.31E-6 / 74 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 / 75 76 76 77 77 REAL(wp), DIMENSION(jpwat), PUBLIC :: diff2 78 DATA diff2/1.74E-7, 4.47E-7, 2.51E-7, 3.89E-7, 1.77E-7, 2.51E-7, 3.89E-7, 3.06E-7, 2.5E-7, 1.5E-7 / 79 80 78 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 79 82 80 !! * Routine accessibility 83 PUBLIC sed_init ! routine called by opa.F90 84 81 PUBLIC sed_ini ! routine called by opa.F90 82 83 !! * Substitutions 84 # include "do_loop_substitute.h90" 85 85 !! $Id$ 86 86 CONTAINS 87 87 88 88 89 SUBROUTINE sed_ini t89 SUBROUTINE sed_ini 90 90 !!---------------------------------------------------------------------- 91 !! *** ROUTINE sed_ini t***91 !! *** ROUTINE sed_ini *** 92 92 !! 93 93 !! ** Purpose : Initialization of sediment module … … 104 104 !! ! 06-07 (C. Ethe) Re-organization 105 105 !!---------------------------------------------------------------------- 106 INTEGER :: ji, jj, ikt, ierr106 INTEGER :: ji, jj, js, jn, jk, ikt, ierr 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 … … 137 136 ! Determination of sediments number of points and allocate global variables 138 137 epkbot(:,:) = 0. 138 gdepbot(:,:) = 0. 139 139 DO_2D( 1, 1, 1, 1 ) 140 140 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 )141 IF( tmask(ji,jj,ikt) == 1 ) epkbot(ji,jj) = e3t_0(ji,jj,ikt) 142 gdepbot(ji,jj) = gdepw_0(ji,jj,ikt+1) 143 143 END_2D 144 144 … … 154 154 155 155 ! 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) ) 156 ALLOCATE( pwcp (jpoce,jpksed,jpwat) ) ; ALLOCATE( pwcp_dta (jpoce,jpwat) ) 157 ALLOCATE( pwcpa(jpoce,jpksed,jpwat) ) ; ALLOCATE( solcpa(jpoce,jpksed,jpsol) ) 158 ALLOCATE( solcp(jpoce,jpksed,jpsol) ) ; ALLOCATE( rainrm_dta(jpoce,jpsol) ) 158 159 ALLOCATE( rainrm(jpoce,jpsol) ) ; ALLOCATE( rainrg(jpoce,jpsol) ) ; ALLOCATE( raintg(jpoce) ) 159 160 ALLOCATE( dzdep(jpoce) ) ; ALLOCATE( iarroce(jpoce) ) ; ALLOCATE( dzkbot(jpoce) ) … … 162 163 ALLOCATE( diff(jpoce,jpksed,jpwat ) ) ; ALLOCATE( irrig(jpoce, jpksed) ) 163 164 ALLOCATE( wacc(jpoce) ) ; ALLOCATE( fecratio(jpoce) ) 165 ALLOCATE( densSW(jpoce) ) ; ALLOCATE( saturco3(jpoce,jpksed) ) 164 166 ALLOCATE( hipor(jpoce,jpksed) ) ; ALLOCATE( co3por(jpoce,jpksed) ) 165 167 ALLOCATE( dz3d(jpoce,jpksed) ) ; ALLOCATE( volw3d(jpoce,jpksed) ) ; ALLOCATE( vols3d(jpoce,jpksed) ) 166 ALLOCATE( sedligand(jpoce, jpksed) ) ; ALLOCATE( saturco3(jpoce,jpksed) ) ; ALLOCATE( densSW(jpoce) ) 168 ALLOCATE( rearatpom(jpoce, jpksed) ) ; ALLOCATE( volc(jpoce,jpksed,jpsol) ) 169 ALLOCATE( Jacobian(jpoce, jpvode*jpksed, jpvode*jpksed) ) 170 ALLOCATE( radsfe2(jpksed) ) ; ALLOCATE( radsnh4(jpksed) ) 171 ALLOCATE( wacc1(jpoce) ) 167 172 168 173 ! Initialization of global variables 169 pwcp (:,:,:) = 0. ; pwcp0 (:,:,:) = 0. ; pwcp_dta (:,:) = 0. 170 solcp (:,:,:) = 0. ; solcp0 (:,:,:) = 0. ; rainrm_dta(:,:) = 0. 171 rainrm(:,: ) = 0. ; rainrg (:,: ) = 0. ; raintg (: ) = 0. 172 dzdep (: ) = 0. ; iarroce(: ) = 0 ; dzkbot (: ) = 0. 173 temp (: ) = 0. ; salt (: ) = 0. ; zkbot (: ) = 0. 174 db (:,: ) = 0. ; densSW (: ) = 0. 175 hipor (:,: ) = 0. ; co3por (:,: ) = 0. ; irrig (:,:) = 0. 176 dz3d (:,: ) = 0. ; volw3d (:,: ) = 0. ; vols3d (:,:) = 0. 177 fecratio(:) = 1E-5 178 sedligand(:,:) = 0.6E-9 174 pwcp (:,:,:) = 0. ; pwcp_dta (:,:) = 0. 175 pwcpa (:,:,:) = 0. ; solcpa(:,:,:) = 0. 176 solcp (:,:,:) = 0. ; rainrm_dta(:,:) = 0. 177 rainrm(:,: ) = 0. ; rainrg(:,: ) = 0. ; raintg (: ) = 0. 178 dzdep (: ) = 0. ; iarroce(: ) = 0 ; dzkbot (: ) = 0. 179 temp (: ) = 0. ; salt (: ) = 0. ; zkbot (: ) = 0. 180 densSW (: ) = 0. ; db (:,:) = 0. 181 hipor (:,: ) = 0. ; co3por (:,: ) = 0. ; irrig (:,:) = 0. 182 dz3d (:,: ) = 0. ; volw3d (:,: ) = 0. ; vols3d (:,:) = 0. 183 fecratio(:) = 1E-5 ; rearatpom(:,:) = 0. 184 radsfe2(:) = 1.0 ; radsnh4(:) = 1.0 179 185 180 186 ! Chemical variables 181 ALLOCATE( akbs (jpoce) ) ; ALLOCATE( ak1s (jpoce) ) ; ALLOCATE( ak2s (jpoce) ) ; ALLOCATE( akws (jpoce) ) 182 ALLOCATE( ak1ps (jpoce) ) ; ALLOCATE( ak2ps (jpoce) ) ; ALLOCATE( ak3ps (jpoce) ) ; ALLOCATE( aksis (jpoce) ) 183 ALLOCATE( aksps (jpoce) ) ; ALLOCATE( ak12s (jpoce) ) ; ALLOCATE( ak12ps(jpoce) ) ; ALLOCATE( ak123ps(jpoce) ) 184 ALLOCATE( borats(jpoce) ) ; ALLOCATE( calcon2(jpoce) ) ; ALLOCATE( sieqs (jpoce) ) 187 ALLOCATE( akbs (jpoce) ) ; ALLOCATE( ak1s (jpoce) ) ; ALLOCATE( ak2s (jpoce) ) ; ALLOCATE( akws (jpoce) ) 188 ALLOCATE( ak1ps (jpoce) ) ; ALLOCATE( ak2ps (jpoce) ) ; ALLOCATE( ak3ps (jpoce) ) ; ALLOCATE( aksis (jpoce) ) 189 ALLOCATE( aksps (jpoce) ) ; ALLOCATE( ak12s (jpoce) ) ; ALLOCATE( ak12ps(jpoce) ) ; ALLOCATE( ak123ps(jpoce) ) 190 ALLOCATE( borats(jpoce) ) ; ALLOCATE( calcon2(jpoce) ) ; ALLOCATE( sieqs (jpoce) ) 185 191 ALLOCATE( aks3s(jpoce) ) ; ALLOCATE( akf3s(jpoce) ) ; ALLOCATE( sulfats(jpoce) ) 186 ALLOCATE( fluorids(jpoce) ) 192 ALLOCATE( fluorids(jpoce) ) ; ALLOCATE( akh2s(jpoce) ) ; ALLOCATE( aknh3(jpoce) ) 187 193 188 194 akbs (:) = 0. ; ak1s (:) = 0. ; ak2s (:) = 0. ; akws (:) = 0. 189 195 ak1ps (:) = 0. ; ak2ps (:) = 0. ; ak3ps (:) = 0. ; aksis (:) = 0. 190 196 aksps (:) = 0. ; ak12s (:) = 0. ; ak12ps(:) = 0. ; ak123ps(:) = 0. 191 borats(:) = 0. ; calcon2(:) = 0. ; sieqs (:) = 0. 197 borats(:) = 0. ; calcon2(:) = 0. ; sieqs (:) = 0. ; akh2s (:) = 0. 192 198 aks3s(:) = 0. ; akf3s(:) = 0. ; sulfats(:) = 0. ; fluorids(:) = 0. 199 aknh3(:) = 0. 193 200 194 201 ! Mass balance calculation 195 ALLOCATE( fromsed(jpoce, jpsol ) ) ; ALLOCATE( tosed(jpoce, jpsol) )202 ALLOCATE( fromsed(jpoce, jpsol+jpads) ) ; ALLOCATE( tosed(jpoce, jpsol+jpads) ) 196 203 197 204 fromsed(:,:) = 0. ; tosed(:,:) = 0. … … 199 206 ! Initialization of sediment geometry 200 207 !------------------------------------ 201 CALL sed_ini t_geom208 CALL sed_ini_geom 202 209 203 210 ! Offline specific mode 204 211 ! --------------------- 205 212 ln_sediment_offline = .FALSE. 213 214 ! Vertical profile of of the adsorption factor for adsorbed species 215 ! ----------------------------------------------------------------- 216 radsfe2(:) = 1.0 / ( 1.0 + adsfe2 * por1(:) / por(:) ) 217 radsnh4(:) = 1.0 / ( 1.0 + adsnh4 * por1(:) / por(:) ) 218 219 ! Initialization of the array for non linear solving 220 ! -------------------------------------------------- 221 222 ALLOCATE( jarr(jpvode*jpksed,2) ) 223 ALLOCATE( jsvode(jpvode), isvode(jptrased) ) 224 jsvode(1) = jwoxy ; jsvode(2) = jwno3 ; jsvode(3) = jwnh4 225 jsvode(4) = jwh2s ; jsvode(5) = jwso4 ; jsvode(6) = jwfe2 226 jsvode(7) = jpwat+jsfeo ; jsvode(8) = jpwat+jsfes 227 isvode(jwoxy) = 1 ; isvode(jwno3) = 2 ; isvode(jwnh4) = 3 228 isvode(jwh2s) = 4 ; isvode(jwso4) = 5 ; isvode(jwfe2) = 6 229 isvode(jpwat+jsfeo) = 7 ; isvode(jpwat+jsfes) = 8 230 DO js = 1, jpvode 231 DO jk = 1, jpksed 232 jn = (jk-1) * jpvode + js 233 jarr(jn,1) = jk 234 jarr(jn,2) = jsvode(js) 235 END DO 236 END DO 237 238 ALLOCATE( stepros(jpoce) ) 206 239 207 240 #if defined key_sed_off … … 219 252 ENDIF 220 253 221 END SUBROUTINE sed_ini t222 223 SUBROUTINE sed_ini t_geom254 END SUBROUTINE sed_ini 255 256 SUBROUTINE sed_ini_geom 224 257 !!---------------------------------------------------------------------- 225 !! *** ROUTINE sed_ini t_geom ***258 !! *** ROUTINE sed_ini_geom *** 226 259 !! 227 260 !! ** Purpose : Initialization of sediment geometry … … 241 274 !---------------------------------------------------------- 242 275 243 IF(lwp) WRITE(numsed,*) ' sed_ini t_geom : Initialization of sediment geometry '276 IF(lwp) WRITE(numsed,*) ' sed_ini_geom : Initialization of sediment geometry ' 244 277 IF(lwp) WRITE(numsed,*) ' ' 245 278 … … 282 315 zsur = - za0 - za1 * sedacr * LOG( COSH( (1-sedkth) / sedacr ) ) 283 316 317 dz(1) = 0.1 284 318 profsedw(1) = 0.0 285 profsed(1) = -dz(1) / 2.319 profsed(1) = -dz(1) / 2. 286 320 DO jk = 2, jpksed 287 321 zw = REAL( jk , wp ) … … 289 323 profsed(jk) = ( zsur + za0 * zt + za1 * sedacr * LOG ( COSH( (zt-sedkth) / sedacr ) ) ) 290 324 profsedw(jk) = ( zsur + za0 * zw + za1 * sedacr * LOG ( COSH( (zw-sedkth) / sedacr ) ) ) 291 END DO292 293 dz(1) = 0.1294 DO jk = 2, jpksed295 325 dz(jk) = profsedw(jk) - profsedw(jk-1) 296 326 END DO 297 327 298 DO jk = 1, jpksed 299 DO ji = 1, jpoce 300 dz3d(ji,jk) = dz(jk) 301 END DO 302 ENDDO 328 DO ji = 1, jpoce 329 dz3d(ji,:) = dz(:) 330 END DO 303 331 304 332 ! Porosity profile [0] … … 371 399 ! -------------------------------------------------------------- 372 400 DO ji = 1, jpoce 373 ztmp1 = 0.117 / ( 1.0 + ( zkbot(ji) / 200.)**3 ) 401 ! ztmp1 = 0.117 / ( 1.0 + ( zkbot(ji) / 200.)**3 ) 402 ztmp1 = 0. 374 403 ztmp2 = 0.006 / ( 1.0 + ( zkbot(ji) / 4000.)**10 ) 375 wacc(ji) = (ztmp1 + ztmp2)/10.0 376 wacc(ji) = ztmp2 / 10.0 404 wacc(ji) = ztmp2+ztmp1 377 405 END DO 378 406 … … 380 408 ! Initialization of time step as function of porosity [cm**2/s] 381 409 !------------------------------------------------------------------ 382 END SUBROUTINE sed_ini t_geom383 384 SUBROUTINE sed_ini t_nam410 END SUBROUTINE sed_ini_geom 411 412 SUBROUTINE sed_ini_nam 385 413 !!---------------------------------------------------------------------- 386 !! *** ROUTINE sed_ini t_nam ***414 !! *** ROUTINE sed_ini_nam *** 387 415 !! 388 416 !! ** Purpose : Initialization of sediment geometry … … 393 421 !!---------------------------------------------------------------------- 394 422 395 CHARACTER(:), ALLOCATABLE :: numnamsed_ref !! Character buffer for referencenamelist sediment396 CHARACTER(:), ALLOCATABLE :: numnamsed_cfg !! Character buffer for configurationnamelist sediment423 INTEGER :: numnamsed_ref = -1 !! Logical units for namelist sediment 424 INTEGER :: numnamsed_cfg = -1 !! Logical units for namelist sediment 397 425 INTEGER :: ios ! Local integer output status for namelist read 398 426 CHARACTER(LEN=20) :: clname … … 409 437 TYPE(PSED) , DIMENSION(jpdia2dsed) :: seddiag2d 410 438 411 NAMELIST/nam_run/ nrseddt,ln_sed_2way439 NAMELIST/nam_run/ln_sed_2way 412 440 NAMELIST/nam_geom/jpksed, sedzmin, sedhmax, sedkth, sedacr, porsurf, porinf, rhox 413 441 NAMELIST/nam_trased/sedsol, sedwat … … 415 443 NAMELIST/nam_inorg/rcopal, dcoef, rccal, ratligc, rcligc 416 444 NAMELIST/nam_poc/redO2, redNo3, redPo4, redC, redfep, rcorgl, rcorgs, & 417 & rcorgr, rcnh4, rch2s, rcfe2, rcfeh2s, rcfes , rcfeso,&418 & xksedo2, xksedno3, xksedfeo, xksedso4419 NAMELIST/nam_btb/dbiot, ln_btbz, dbtbzsc, adsnh4, ln_irrig, xirrzsc445 & rcorgr, rcnh4, rch2s, rcfe2, rcfeh2s, rcfeso, rcfesp, & 446 & rcfesd, xksedo2, xksedno3, xksedfeo, xksedso4 447 NAMELIST/nam_btb/dbiot, ln_btbz, dbtbzsc, adsnh4, adsfe2, ln_irrig, xirrzsc 420 448 NAMELIST/nam_rst/ln_rst_sed, cn_sedrst_indir, cn_sedrst_outdir, cn_sedrst_in, cn_sedrst_out 421 449 … … 423 451 !------------------------------------------------------- 424 452 425 IF(lwp) WRITE(numsed,*) ' sed_ini t_nam : Read namelists '453 IF(lwp) WRITE(numsed,*) ' sed_ini_nam : Read namelists ' 426 454 IF(lwp) WRITE(numsed,*) ' ' 427 455 … … 437 465 !--------------------------------- 438 466 clname = 'namelist_sediment' 439 IF(lwp) WRITE(numsed,*) ' sed_ini t_nam : read SEDIMENT namelist'467 IF(lwp) WRITE(numsed,*) ' sed_ini_nam : read SEDIMENT namelist' 440 468 IF(lwp) WRITE(numsed,*) ' ~~~~~~~~~~~~~~' 441 CALL load_nml( numnamsed_ref, TRIM( clname )//'_ref', numout, lwm)442 CALL load_nml( numnamsed_cfg, TRIM( clname )//'_cfg', numout, lwm)469 CALL ctl_opn( numnamsed_ref, TRIM( clname )//'_ref', 'OLD' , 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. ) 470 CALL ctl_opn( numnamsed_cfg, TRIM( clname )//'_cfg', 'OLD' , 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. ) 443 471 444 472 nitsed000 = nittrc000 445 473 nitsedend = nitend 446 474 ! Namelist nam_run 475 REWIND( numnamsed_ref ) ! Namelist nam_run in reference namelist : Pisces variables 447 476 READ ( numnamsed_ref, nam_run, IOSTAT = ios, ERR = 901) 448 477 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_run in reference namelist' ) 449 478 479 REWIND( numnamsed_cfg ) ! Namelist nam_run in reference namelist : Pisces variables 450 480 READ ( numnamsed_cfg, nam_run, IOSTAT = ios, ERR = 902) 451 481 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_run in configuration namelist' ) … … 453 483 IF (lwp) THEN 454 484 WRITE(numsed,*) ' namelist nam_run' 455 WRITE(numsed,*) ' Nb of iterations for fast species nrseddt = ', nrseddt456 485 WRITE(numsed,*) ' 2-way coupling between PISCES and Sed ln_sed_2way = ', ln_sed_2way 457 486 ENDIF … … 459 488 IF ( ln_p5z .AND. ln_sed_2way ) CALL ctl_stop( '2 ways coupling with sediment cannot be activated with PISCES-QUOTA' ) 460 489 490 REWIND( numnamsed_ref ) ! Namelist nam_geom in reference namelist : Pisces variables 461 491 READ ( numnamsed_ref, nam_geom, IOSTAT = ios, ERR = 903) 462 492 903 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_geom in reference namelist' ) 463 493 494 REWIND( numnamsed_cfg ) ! Namelist nam_geom in reference namelist : Pisces variables 464 495 READ ( numnamsed_cfg, nam_geom, IOSTAT = ios, ERR = 904) 465 496 904 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_geom in configuration namelist' ) … … 480 511 dtsed = rDt_trc 481 512 513 REWIND( numnamsed_ref ) ! Namelist nam_trased in reference namelist : Pisces variables 482 514 READ ( numnamsed_ref, nam_trased, IOSTAT = ios, ERR = 905) 483 515 905 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_trased in reference namelist' ) 484 516 517 REWIND( numnamsed_cfg ) ! Namelist nam_trased in reference namelist : Pisces variables 485 518 READ ( numnamsed_cfg, nam_trased, IOSTAT = ios, ERR = 906) 486 519 906 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_trased in configuration namelist' ) … … 511 544 ENDIF 512 545 546 REWIND( numnamsed_ref ) ! Namelist nam_diased in reference namelist : Pisces variables 513 547 READ ( numnamsed_ref, nam_diased, IOSTAT = ios, ERR = 907) 514 548 907 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diased in reference namelist' ) 515 549 550 REWIND( numnamsed_cfg ) ! Namelist nam_diased in reference namelist : Pisces variables 516 551 READ ( numnamsed_cfg, nam_diased, IOSTAT = ios, ERR = 908) 517 552 908 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diased in configuration namelist' ) … … 551 586 ! Inorganic chemistry parameters 552 587 !---------------------------------- 588 REWIND( numnamsed_ref ) ! Namelist nam_inorg in reference namelist : Pisces variables 553 589 READ ( numnamsed_ref, nam_inorg, IOSTAT = ios, ERR = 909) 554 590 909 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_inorg in reference namelist' ) 555 591 592 REWIND( numnamsed_cfg ) ! Namelist nam_inorg in reference namelist : Pisces variables 556 593 READ ( numnamsed_cfg, nam_inorg, IOSTAT = ios, ERR = 910) 557 594 910 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_inorg in configuration namelist' ) … … 575 612 ! Additional parameter linked to POC/O2/No3/Po4 576 613 !---------------------------------------------- 614 REWIND( numnamsed_ref ) ! Namelist nam_poc in reference namelist : Pisces variables 577 615 READ ( numnamsed_ref, nam_poc, IOSTAT = ios, ERR = 911) 578 616 911 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_poc in reference namelist' ) 579 617 618 REWIND( numnamsed_cfg ) ! Namelist nam_poc in reference namelist : Pisces variables 580 619 READ ( numnamsed_cfg, nam_poc, IOSTAT = ios, ERR = 912) 581 620 912 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_poc in configuration namelist' ) … … 595 634 WRITE(numsed,*) ' reactivity for Fe2+ rcfe2 = ', rcfe2 596 635 WRITE(numsed,*) ' reactivity for FeOH/H2S rcfeh2s = ', rcfeh2s 597 WRITE(numsed,*) ' reactivity for Fe2+/H2S rcfes = ', rcfes598 636 WRITE(numsed,*) ' reactivity for FeS/O2 rcfeso = ', rcfeso 637 WRITE(numsed,*) ' Precipitation of FeS rcfesp = ', rcfesp 638 WRITE(numsed,*) ' Dissolution of FeS rcfesd = ', rcfesd 599 639 WRITE(numsed,*) ' Half-sat. cste for oxic remin xksedo2 = ', xksedo2 600 640 WRITE(numsed,*) ' Half-sat. cste for denit. xksedno3 = ', xksedno3 … … 617 657 reac_fe2 = rcfe2 / ryear 618 658 reac_feh2s = rcfeh2s/ ryear 619 reac_fes = rcfes / ryear620 659 reac_feso = rcfeso / ryear 660 reac_fesp = rcfesp / ryear 661 reac_fesd = rcfesd / ryear 662 621 663 622 664 ! reactivity rc in [l.mol-1.s-1] … … 625 667 ! Bioturbation parameter 626 668 !------------------------ 669 REWIND( numnamsed_ref ) ! Namelist nam_btb in reference namelist : Pisces variables 627 670 READ ( numnamsed_ref, nam_btb, IOSTAT = ios, ERR = 913) 628 671 913 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_btb in reference namelist' ) 629 672 673 REWIND( numnamsed_cfg ) ! Namelist nam_btb in reference namelist : Pisces variables 630 674 READ ( numnamsed_cfg, nam_btb, IOSTAT = ios, ERR = 914) 631 675 914 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_btb in configuration namelist' ) … … 637 681 WRITE(numsed,*) ' coefficient for btb attenuation dbtbzsc = ', dbtbzsc 638 682 WRITE(numsed,*) ' Adsorption coefficient of NH4 adsnh4 = ', adsnh4 683 WRITE(numsed,*) ' Adsorption coefficient of Fe2 adsfe2 = ', adsfe2 639 684 WRITE(numsed,*) ' Bioirrigation in sediment ln_irrig = ', ln_irrig 640 685 WRITE(numsed,*) ' coefficient for irrig attenuation xirrzsc = ', xirrzsc … … 644 689 ! Initial value (t=0) for sediment pore water and solid components 645 690 !---------------------------------------------------------------- 691 REWIND( numnamsed_ref ) ! Namelist nam_rst in reference namelist : Pisces variables 646 692 READ ( numnamsed_ref, nam_rst, IOSTAT = ios, ERR = 915) 647 693 915 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_rst in reference namelist' ) 648 694 695 REWIND( numnamsed_cfg ) ! Namelist nam_rst in reference namelist : Pisces variables 649 696 READ ( numnamsed_cfg, nam_rst, IOSTAT = ios, ERR = 916) 650 697 916 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_rst in configuration namelist' ) … … 655 702 WRITE(numsed,*) ' ' 656 703 ENDIF 657 nn_dtsed = 1 658 659 660 END SUBROUTINE sed_init_nam 704 705 CLOSE( numnamsed_cfg ) 706 CLOSE( numnamsed_ref ) 707 708 END SUBROUTINE sed_ini_nam 661 709 662 710 END MODULE sedini -
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/sedinitrc.F90
r12377 r15075 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/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/sedinorg.F90
r14385 r15075 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 seddsr14 11 USE lib_mpp ! distribued memory computing library 15 12 USE lib_fortran … … 23 20 CONTAINS 24 21 25 SUBROUTINE sed_inorg( kt , knt )22 SUBROUTINE sed_inorg( kt ) 26 23 !!---------------------------------------------------------------------- 27 24 !! *** ROUTINE sed_inorg *** … … 46 43 !!---------------------------------------------------------------------- 47 44 !! Arguments 48 INTEGER, INTENT(in) :: kt, knt ! number of iteration45 INTEGER, INTENT(in) :: kt ! time step 49 46 ! --- 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, zbeta, zgamma 57 REAL(wp) :: zsatur, zsatur2, znusil, zsolcpcl, zsolcpsi 47 INTEGER :: ji,jk ! dummy looop indices 48 REAL(wp) :: zsieq 49 REAL(wp) :: zsolid1, zreasat 50 REAL(wp) :: zsatur, zsatur2, znusil, zsolcpcl, zsolcpsi, zexcess 58 51 !! 59 52 !!---------------------------------------------------------------------- 60 53 61 54 IF( ln_timing ) CALL timing_start('sed_inorg') 55 56 IF( kt == nitsed000 ) THEN 57 IF (lwp) WRITE(numsed,*) ' sed_inorg : Dissolution of CaCO3 and BSi ' 58 IF (lwp) WRITE(numsed,*) ' ' 59 ENDIF 62 60 ! 63 IF( kt == nitsed000 .AND. knt == 1 ) THEN 64 IF (lwp) THEN 65 WRITE(numsed,*) ' sed_inorg : Dissolution reaction ' 66 WRITE(numsed,*) ' ' 67 ENDIF 68 ! ! 69 ENDIF 61 DO ji = 1, jpoce 62 ! ----------------------------------------------- 63 ! Computation of Si solubility 64 ! Param of Ridgwell et al. 2002 65 ! ----------------------------------------------- 70 66 71 ! Initializations72 !----------------------73 74 zrearat1(:,:) = 0. ; zundsat(:,:) = 0.75 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 67 zsolcpcl = 0.0 86 68 zsolcpsi = 0.0 … … 90 72 END DO 91 73 zsolcpsi = MAX( zsolcpsi, rtrn ) 92 zsieq(ji) = sieqs(ji) * MAX(0.25, 1.0 - (0.045 * zsolcpcl / zsolcpsi )**0.58 ) 93 zsieq(ji) = MAX( rtrn, sieqs(ji) ) 74 zsieq = sieqs(ji) * MAX(0.25, 1.0 - (0.045 * zsolcpcl / zsolcpsi )**0.58 ) 75 76 !---------------------------------------------------------- 77 ! 5. Beginning of Pore Water diffusion and solid reaction 78 !--------------------------------------------------------- 79 80 !----------------------------------------------------------------------------- 81 ! For jk=2,jpksed, and for couple 82 ! 1 : jwsil/jsopal ( SI/Opal ) 83 ! 2 : jsclay/jsclay ( clay/clay ) 84 ! 3 : jwoxy/jspoc ( O2/POC ) 85 ! reaction rate is a function of solid=concentration in solid reactif in [mol/l] 86 ! and undersaturation in [mol/l]. 87 ! Solid weight fractions should be in ie [mol/l]) 88 ! second member and solution are in zundsat variable 89 !------------------------------------------------------------------------- 90 DO jk = 2, jpksed 91 zsolid1 = volc(ji,jk,jsopal) * solcp(ji,jk,jsopal) 92 zsatur = MAX(0., ( zsieq - pwcp(ji,jk,jwsil) ) / zsieq ) 93 zsatur2 = (1.0 + temp(ji) / 400.0 )**37 94 znusil = ( 0.225 * ( 1.0 + temp(ji) / 15.) * zsatur + 0.775 * zsatur2 * zsatur**9.25 ) 95 solcp(ji,jk,jsopal) = solcp(ji,jk,jsopal) - reac_sil * znusil * dtsed * solcp(ji,jk,jsopal) 96 pwcp(ji,jk,jwsil) = pwcp(ji,jk,jwsil) + reac_sil * znusil * dtsed * zsolid1 97 END DO 94 98 END DO 95 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 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 DO ji = 1, jpoce130 zrearat1(ji,:) = 0.131 zrearat2(ji,:) = 0.132 ENDDO133 134 ! left hand side of coefficient matrix135 DO jk = 2, jpksed136 DO ji = 1, jpoce137 IF ( zundsat(ji,jk) > 0. ) THEN138 zsolid1 = zvolc(ji,jk,jsopal) * solcp(ji,jk,jsopal)139 zsatur = zundsat(ji,jk) / zsieq(ji)140 zsatur2 = (1.0 + temp(ji) / 400.0 )**37141 znusil = ( 0.225 * ( 1.0 + temp(ji) / 15.) + 0.775 * zsatur2 * zsatur**8.25 ) / zsieq(ji)142 zbeta = 1.0 - reac_sil * znusil * dtsed2 * zundsat(ji,jk) + reac_sil * znusil * dtsed2 * zsolid1143 zgamma = - reac_sil * znusil * dtsed2 * zundsat(ji,jk)144 zundsat(ji,jk) = ( - zbeta + SQRT( zbeta**2 - 4.0 * zgamma ) ) / ( 2.0 * reac_sil * znusil * dtsed2 )145 solcp(ji,jk,jsopal ) = zsolid1 / ( 1. + reac_sil * znusil * dtsed2 * zundsat(ji,jk) ) / zvolc(ji,jk,jsopal)146 pwcp(ji,jk,jwsil) = zsieq(ji) - zundsat(ji,jk)147 ENDIF148 ENDDO149 ENDDO150 99 151 100 !--------------------------------------------------------------- … … 154 103 155 104 ! computes co3por from the updated pwcp concentrations (note [co3por] = mol/l) 156 157 CALL sed_co3( kt )158 159 105 ! *densSW(l)**2 converts aksps [mol2/kg sol2] into [mol2/l2] to get [undsat] in [mol/l] 160 DO jk = 1, jpksed 161 DO ji = 1, jpoce 162 zco3eq(ji) = aksps(ji) * densSW(ji) * densSW(ji) / ( calcon2(ji) + rtrn ) 163 zco3eq(ji) = MAX( rtrn, zco3eq(ji) ) 164 zundsat(ji,jk) = ( zco3eq(ji) - co3por(ji,jk) ) / zco3eq(ji) 165 saturco3(ji,jk) = zundsat(ji,jk) 166 ENDDO 167 ENDDO 168 169 DO jk = 2, jpksed 170 DO ji = 1, jpoce 171 IF ( zundsat(ji,jk) > 0. ) THEN 172 zsolid1 = zvolc(ji,jk,jscal) * solcp(ji,jk,jscal) 173 zbeta = 1.0 - reac_cal * dtsed2 * zundsat(ji,jk) + reac_cal * dtsed2 * zsolid1 174 zgamma = -reac_cal * dtsed2 * zundsat(ji,jk) 175 zundsat(ji,jk) = ( -zbeta + SQRT( zbeta**2 - 4.0 * zgamma ) ) / ( 2.0 * reac_cal * dtsed2 ) 176 saturco3(ji,jk) = zundsat(ji,jk) 177 zreasat = reac_cal * dtsed2 * zundsat(ji,jk) * zsolid1 / ( 1. + reac_cal * dtsed2 * zundsat(ji,jk) ) 178 solcp(ji,jk,jscal) = solcp(ji,jk,jscal) - zreasat / zvolc(ji,jk,jscal) 179 ! For DIC 180 pwcp(ji,jk,jwdic) = pwcp(ji,jk,jwdic) + zreasat 181 ! For alkalinity 182 pwcp(ji,jk,jwalk) = pwcp(ji,jk,jwalk) + 2.0 * zreasat 183 ENDIF 106 DO ji = 1, jpoce 107 saturco3(ji,:) = 1.0 - co3por(ji,:) * calcon2(ji) / ( aksps(ji) * densSW(ji) * densSW(ji) + rtrn ) 108 DO jk = 2, jpksed 109 zsolid1 = volc(ji,jk,jscal) * solcp(ji,jk,jscal) 110 zexcess = MAX( 0., saturco3(ji,jk) ) 111 zreasat = reac_cal * dtsed * zexcess * zsolid1 112 solcp(ji,jk,jscal) = solcp(ji,jk,jscal) - zreasat / volc(ji,jk,jscal) 113 ! For DIC 114 pwcp(ji,jk,jwdic) = pwcp(ji,jk,jwdic) + zreasat 115 ! For alkalinity 116 pwcp(ji,jk,jwalk) = pwcp(ji,jk,jwalk) + 2.0 * zreasat 184 117 END DO 185 118 END DO -
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/sedmat.F90
r10222 r15075 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 16 PUBLIC sed_mat_dsr 17 PUBLIC sed_mat_dsrjac 18 PUBLIC sed_mat_dsri 19 PUBLIC sed_mat_btb 20 21 INTEGER, PARAMETER :: nmax = 60 23 22 24 23 … … 26 25 CONTAINS 27 26 28 SUBROUTINE sed_mat_dsr( nvar, ndim, nlev, preac, psms, psol, dtsed_in )27 SUBROUTINE sed_mat_dsr( ji, nvar, dtsed_in ) 29 28 !!--------------------------------------------------------------------- 30 29 !! *** ROUTINE sed_mat_dsr *** … … 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 ) 50 INTEGER , INTENT(in) :: ji, nvar ! number of variable 51 58 52 REAL(wp), INTENT(in) :: dtsed_in 59 53 60 54 !---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 55 INTEGER :: jk, jn 56 REAL(wp), DIMENSION(jpksed) :: za, zb, zc 65 57 66 58 REAL(wp) :: aplus,aminus … … 79 71 jn = nvar 80 72 ! first sediment level 81 DO ji = 1, ndim 82 aplus = ( ( volw3d(ji,1) / ( dz3d(ji,1) ) ) + & 73 aplus = ( ( volw3d(ji,1) / ( dz3d(ji,1) ) ) + & 83 74 ( 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 86 87 za(ji,1) = 0. 88 zb(ji,1) = 1. + rplus 89 zc(ji,1) = -rplus 90 ENDDO 75 dxplus = ( dz3d(ji,1) + dz3d(ji,2) ) / 2. 76 rplus = ( dtsed_in / ( volw3d(ji,1) ) ) * diff(ji,1,jn) * aplus / dxplus 77 78 za(1) = 0. 79 zb(1) = rplus 80 zc(1) = -rplus 91 81 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 104 ! 105 za(ji,jk) = -rminus 106 zb(ji,jk) = 1. + rminus + rplus 107 zc(ji,jk) = -rplus 108 END DO 82 DO jk = 2, jpksed - 1 83 aminus = ( ( volw3d(ji,jk-1) / ( dz3d(ji,jk-1) ) ) + & 84 & ( volw3d(ji,jk ) / ( dz3d(ji,jk ) ) ) ) / 2. 85 dxminus = ( dz3d(ji,jk-1) + dz3d(ji,jk) ) / 2. 86 87 aplus = ( ( volw3d(ji,jk ) / ( dz3d(ji,jk ) ) ) + & 88 & ( volw3d(ji,jk+1) / ( dz3d(ji,jk+1) ) ) ) / 2. 89 dxplus = ( dz3d(ji,jk) + dz3d(ji,jk+1) ) / 2 90 ! 91 rminus = ( dtsed_in / volw3d(ji,jk) ) * diff(ji,jk-1,jn) * aminus / dxminus 92 rplus = ( dtsed_in / volw3d(ji,jk) ) * diff(ji,jk,jn) * aplus / dxplus 93 ! 94 za(jk) = -rminus 95 zb(jk) = rminus + rplus 96 zc(jk) = -rplus 109 97 END DO 110 98 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 116 ! 117 za(ji,nlev) = -rminus 118 zb(ji,nlev) = 1. + rminus 119 zc(ji,nlev) = 0. 120 END DO 121 99 aminus = ( ( volw3d(ji,jpksed-1) / dz3d(ji,jpksed-1) ) + & 100 & ( volw3d(ji,jpksed) / dz3d(ji,jpksed) ) ) / 2. 101 dxminus = ( dz3d(ji,jpksed-1) + dz3d(ji,jpksed) ) / 2. 102 rminus = ( dtsed_in / volw3d(ji,jpksed) ) * diff(ji,jpksed-1,jn) * aminus / dxminus 103 ! 104 za(jpksed) = -rminus 105 zb(jpksed) = rminus 106 zc(jpksed) = 0. 122 107 123 108 ! solves tridiagonal system of linear equations 124 109 ! ----------------------------------------------- 125 110 126 zr (:,:) = psol(:,:) + psms(:,:) 127 zb (:,:) = zb(:,:) + preac(:,:) 128 zbet(: ) = zb(:,1) 129 psol(:,1) = zr(:,1) / zbet(:) 130 131 ! 132 DO jk = 2, nlev 133 DO ji = 1, ndim 134 zgamm(ji,jk) = zc(ji,jk-1) / zbet(ji) 135 zbet(ji) = zb(ji,jk) - za(ji,jk) * zgamm(ji,jk) 136 psol(ji,jk) = ( zr(ji,jk) - za(ji,jk) * psol(ji,jk-1) ) / zbet(ji) 111 pwcpa(ji,1,jn) = pwcpa(ji,1,jn) - ( zc(1) * pwcp(ji,2,jn) + zb(1) * pwcp(ji,1,jn) ) 112 DO jk = 2, jpksed - 1 113 pwcpa(ji,jk,jn) = pwcpa(ji,jk,jn) - ( zc(jk) * pwcp(ji,jk+1,jn) + za(jk) * pwcp(ji,jk-1,jn) & 114 & + zb(jk) * pwcp(ji,jk,jn) ) 115 ENDDO 116 pwcpa(ji,jpksed,jn) = pwcpa(ji,jpksed,jn) - ( za(jk) * pwcp(ji,jk-1,jn) + zb(jk) * pwcp(ji,jk,jn) ) 117 118 IF( ln_timing ) CALL timing_stop('sed_mat_dsr') 119 120 END SUBROUTINE sed_mat_dsr 121 122 SUBROUTINE sed_mat_dsrjac( ji, nvar, dtsed_in ) 123 !!--------------------------------------------------------------------- 124 !! *** ROUTINE sed_mat_dsrjac *** 125 !! 126 !! ** Purpose : solves tridiagonal system of linear equations 127 !! 128 !! ** Method : 129 !! 1 - computes left hand side of linear system of equations 130 !! for dissolution reaction 131 !! For mass balance in kbot+sediment : 132 !! dz3d (:,1) = dz(1) = 0.5 cm 133 !! volw3d(:,1) = dzkbot ( see sedini.F90 ) 134 !! dz(2) = 0.3 cm 135 !! dz3d(:,2) = 0.3 + dzdep ( see seddsr.F90 ) 136 !! volw3d(:,2) and vols3d(l,2) are thickened ( see 137 !seddsr.F90 ) 138 !! 139 !! 2 - forward/backward substitution. 140 !! 141 !! History : 142 !! ! 04-10 (N. Emprin, M. Gehlen ) original 143 !! ! 06-04 (C. Ethe) Module Re-organization 144 !!---------------------------------------------------------------------- 145 !! * Arguments 146 INTEGER , INTENT(in) :: ji, nvar ! number of variable 147 148 REAL(wp), INTENT(in) :: dtsed_in 149 150 !---Local declarations 151 INTEGER :: jk, jn, jnn 152 REAL(wp), DIMENSION(jpksed) :: za, zb, zc 153 154 REAL(wp) :: aplus,aminus 155 REAL(wp) :: rplus,rminus 156 REAL(wp) :: dxplus,dxminus 157 158 !---------------------------------------------------------------------- 159 160 IF( ln_timing ) CALL timing_start('sed_mat_dsrjac') 161 162 ! Computation left hand side of linear system of 163 ! equations for dissolution reaction 164 !--------------------------------------------- 165 166 167 jn = nvar 168 ! first sediment level 169 aplus = ( ( volw3d(ji,1) / ( dz3d(ji,1) ) ) + & 170 ( volw3d(ji,2) / ( dz3d(ji,2) ) ) ) / 2. 171 dxplus = ( dz3d(ji,1) + dz3d(ji,2) ) / 2. 172 rplus = ( dtsed_in / ( volw3d(ji,1) ) ) * diff(ji,1,jn) * aplus / dxplus 173 174 za(1) = 0. 175 zb(1) = rplus 176 zc(1) = -rplus 177 178 DO jk = 2, jpksed - 1 179 aminus = ( ( volw3d(ji,jk-1) / ( dz3d(ji,jk-1) ) ) + & 180 & ( volw3d(ji,jk ) / ( dz3d(ji,jk ) ) ) ) / 2. 181 dxminus = ( dz3d(ji,jk-1) + dz3d(ji,jk) ) / 2. 182 183 aplus = ( ( volw3d(ji,jk ) / ( dz3d(ji,jk ) ) ) + & 184 & ( volw3d(ji,jk+1) / ( dz3d(ji,jk+1) ) ) ) / 2. 185 dxplus = ( dz3d(ji,jk) + dz3d(ji,jk+1) ) / 2 186 ! 187 rminus = ( dtsed_in / volw3d(ji,jk) ) * diff(ji,jk-1,jn) * aminus / dxminus 188 rplus = ( dtsed_in / volw3d(ji,jk) ) * diff(ji,jk,jn) * aplus / dxplus 189 ! 190 za(jk) = -rminus 191 zb(jk) = rminus + rplus 192 zc(jk) = -rplus 193 END DO 194 195 aminus = ( ( volw3d(ji,jpksed-1) / dz3d(ji,jpksed-1) ) + & 196 & ( volw3d(ji,jpksed) / dz3d(ji,jpksed) ) ) / 2. 197 dxminus = ( dz3d(ji,jpksed-1) + dz3d(ji,jpksed) ) / 2. 198 rminus = ( dtsed_in / volw3d(ji,jpksed) ) * diff(ji,jpksed-1,jn) * aminus / dxminus 199 ! 200 za(jpksed) = -rminus 201 zb(jpksed) = rminus 202 zc(jpksed) = 0. 203 204 ! solves tridiagonal system of linear equations 205 206 IF (isvode(jn) > 0) THEN 207 jnn = isvode(jn) 208 Jacobian(ji, jnn, jnn) = Jacobian(ji,jnn,jnn) - zb(1) 209 Jacobian(ji, jnn, jpvode + jnn) = Jacobian(ji, jnn, jpvode + jnn) -zc(1) 210 DO jk = 2, jpksed - 1 211 Jacobian(ji, (jk-1) * jpvode + jnn, (jk-2) * jpvode + jnn) = Jacobian(ji, (jk-1) * jpvode + jnn, (jk-2) * jpvode + jnn) - za(jk) 212 Jacobian(ji, (jk-1) * jpvode + jnn, (jk-1) * jpvode + jnn) = Jacobian(ji, (jk-1) * jpvode + jnn, (jk-1) * jpvode + jnn) - zb(jk) 213 Jacobian(ji, (jk-1) * jpvode + jnn, (jk) * jpvode + jnn) = Jacobian(ji, (jk-1) * jpvode + jnn, (jk) * jpvode + jnn) - zc(jk) 137 214 END DO 138 ENDDO 139 ! 140 DO jk = nlev - 1, 1, -1 141 DO ji = 1,ndim 142 psol(ji,jk) = psol(ji,jk) - zgamm(ji,jk+1) * psol(ji,jk+1) 143 END DO 144 ENDDO 145 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 ) 215 Jacobian(ji, (jpksed-1) * jpvode + jnn, (jpksed-2) * jpvode + jnn) = Jacobian(ji, (jpksed-1) * jpvode + jnn, (jpksed-2) * jpvode + jnn) - za(jpksed) 216 Jacobian(ji, (jpksed-1) * jpvode + jnn, (jpksed-1) * jpvode + jnn) = Jacobian(ji, (jpksed-1) * jpvode + jnn, (jpksed-1) * jpvode + jnn) - zb(jpksed) 217 ENDIF 218 219 IF( ln_timing ) CALL timing_stop('sed_mat_dsrjac') 220 221 END SUBROUTINE sed_mat_dsrjac 222 223 224 225 226 SUBROUTINE sed_mat_btb( nlev, nvar, psol, preac, dtsed_in ) 152 227 !!--------------------------------------------------------------------- 153 228 !! *** ROUTINE sed_mat_btb *** … … 167 242 !! * Arguments 168 243 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 244 nlev, nvar ! number of sediment levels 245 246 REAL(wp), DIMENSION(jpoce,nlev,nvar), INTENT(inout) :: & 247 psol, preac 175 248 176 249 REAL(wp), INTENT(in) :: dtsed_in … … 181 254 182 255 REAL(wp) :: & 183 aplus,aminus , & 256 aplus,aminus , & 184 257 rplus,rminus , & 185 dxplus,dxminus 258 dxplus,dxminus 186 259 187 260 REAL(wp), DIMENSION(nlev) :: za, zb, zc 188 REAL(wp), DIMENSION( ndim,nlev) :: zr189 REAL(wp), DIMENSION(nmax) :: zgamm 261 REAL(wp), DIMENSION(jpoce,nlev) :: zr 262 REAL(wp), DIMENSION(nmax) :: zgamm 190 263 REAL(wp) :: zbet 191 264 192 193 265 !---------------------------------------------------------------------- 194 266 … … 201 273 202 274 ! first sediment level 203 DO ji = 1, ndim275 DO ji = 1, jpoce 204 276 aplus = ( ( vols(2) / dz(2) ) + ( vols(3) / dz(3) ) ) / 2. 205 277 dxplus = ( dz(2) + dz(3) ) / 2. 206 rplus = ( dtsed_in / vols(2) ) * db(ji,2) * aplus / dxplus 278 rplus = ( dtsed_in / vols(2) ) * db(ji,2) * aplus / dxplus 207 279 208 280 za(1) = 0. 209 zb(1) = 1. + rplus 281 zb(1) = 1. + rplus 210 282 zc(1) = -rplus 211 283 212 213 284 DO jk = 2, nlev - 1 214 285 aminus = ( ( vols(jk) / dz(jk) ) + ( vols(jk+1) / dz(jk+1) ) ) / 2. … … 221 292 ! 222 293 za(jk) = -rminus 223 zb(jk) = 1. + rminus + rplus 294 zb(jk) = 1. + rminus + rplus 224 295 zc(jk) = -rplus 296 225 297 ENDDO 226 298 227 299 aminus = ( ( vols(nlev) / dz(nlev) ) + ( vols(nlev+1) / dz(nlev+1) ) ) / 2. 228 300 dxminus = ( dz(nlev) + dz(nlev+1) ) / 2. … … 233 305 zc(nlev) = 0. 234 306 235 236 307 ! solves tridiagonal system of linear equations 237 308 ! ----------------------------------------------- 238 309 DO jn = 1, nvar 239 240 310 DO jk = 1, nlev 241 zr (ji,jk)= psol(ji,jk,jn)311 zr(ji,jk) = psol(ji,jk,jn) 242 312 END DO 243 zbet = zb(1) 313 zbet = zb(1) - preac(ji,1,jn) * dtsed_in 244 314 psol(ji,1,jn) = zr(ji,1) / zbet 245 315 ! 246 316 DO jk = 2, nlev 247 317 zgamm(jk) = zc(jk-1) / zbet 248 zbet = zb(jk) - za(jk) * zgamm(jk)318 zbet = zb(jk) - preac(ji,jk,jn) * dtsed_in - za(jk) * zgamm(jk) 249 319 psol(ji,jk,jn) = ( zr(ji,jk) - za(jk) * psol(ji,jk-1,jn) ) / zbet 250 320 ENDDO 251 321 ! 252 322 DO jk = nlev - 1, 1, -1 253 323 psol(ji,jk,jn) = psol(ji,jk,jn) - zgamm(jk+1) * psol(ji,jk+1,jn) 254 324 ENDDO 255 256 ENDDO 257 258 END DO 259 ! 260 IF( ln_timing ) CALL timing_stop('sed_mat_btb') 325 END DO 326 END DO 327 ! 328 IF( ln_timing ) CALL timing_stop('sed_mat_btb') 261 329 262 330 263 331 END SUBROUTINE sed_mat_btb 264 332 333 SUBROUTINE sed_mat_dsri( nvar, preac, psms, dtsed_in ) 334 !!--------------------------------------------------------------------- 335 !! *** ROUTINE sed_mat_dsr *** 336 !! 337 !! ** Purpose : solves tridiagonal system of linear equations 338 !! 339 !! ** Method : 340 !! 1 - computes left hand side of linear system of equations 341 !! for dissolution reaction 342 !! For mass balance in kbot+sediment : 343 !! dz3d (:,1) = dz(1) = 0.5 cm 344 !! volw3d(:,1) = dzkbot ( see sedini.F90 ) 345 !! dz(2) = 0.3 cm 346 !! dz3d(:,2) = 0.3 + dzdep ( see seddsr.F90 ) 347 !! volw3d(:,2) and vols3d(l,2) are thickened ( see 348 !seddsr.F90 ) 349 !! 350 !! 2 - forward/backward substitution. 351 !! 352 !! History : 353 !! ! 04-10 (N. Emprin, M. Gehlen ) original 354 !! ! 06-04 (C. Ethe) Module Re-organization 355 !!---------------------------------------------------------------------- 356 !! * Arguments 357 INTEGER , INTENT(in) :: nvar ! number of variable 358 359 REAL(wp), DIMENSION(jpoce,jpksed), INTENT(in ) :: preac ! reaction rates 360 REAL(wp), DIMENSION(jpoce,jpksed), INTENT(in ) :: psms ! reaction rates 361 REAL(wp), INTENT(in) :: dtsed_in 362 363 !---Local declarations 364 INTEGER :: ji, jk, jn 365 REAL(wp), DIMENSION(jpoce,jpksed) :: za, zb, zc, zr 366 REAL(wp), DIMENSION(jpoce) :: zbet 367 REAL(wp), DIMENSION(jpoce,nmax) :: zgamm 368 369 REAL(wp) :: aplus,aminus 370 REAL(wp) :: rplus,rminus 371 REAL(wp) :: dxplus,dxminus 372 373 !---------------------------------------------------------------------- 374 375 IF( ln_timing ) CALL timing_start('sed_mat_dsri') 376 377 ! Computation left hand side of linear system of 378 ! equations for dissolution reaction 379 !--------------------------------------------- 380 381 382 jn = nvar 383 ! first sediment level 384 DO ji = 1, jpoce 385 aplus = ( ( volw3d(ji,1) / ( dz3d(ji,1) ) ) + & 386 ( volw3d(ji,2) / ( dz3d(ji,2) ) ) ) / 2. 387 dxplus = ( dz3d(ji,1) + dz3d(ji,2) ) / 2. 388 rplus = ( dtsed_in / ( volw3d(ji,1) ) ) * diff(ji,1,jn) * aplus / dxplus 389 390 za(ji,1) = 0. 391 zb(ji,1) = 1. + rplus 392 zc(ji,1) = -rplus 393 ENDDO 394 395 DO jk = 2, jpksed - 1 396 DO ji = 1, jpoce 397 aminus = ( ( volw3d(ji,jk-1) / ( dz3d(ji,jk-1) ) ) + & 398 & ( volw3d(ji,jk ) / ( dz3d(ji,jk ) ) ) ) / 2. 399 dxminus = ( dz3d(ji,jk-1) + dz3d(ji,jk) ) / 2. 400 401 aplus = ( ( volw3d(ji,jk ) / ( dz3d(ji,jk ) ) ) + & 402 & ( volw3d(ji,jk+1) / ( dz3d(ji,jk+1) ) ) ) / 2. 403 dxplus = ( dz3d(ji,jk) + dz3d(ji,jk+1) ) / 2 404 ! 405 rminus = ( dtsed_in / volw3d(ji,jk) ) * diff(ji,jk-1,jn) * aminus / dxminus 406 rplus = ( dtsed_in / volw3d(ji,jk) ) * diff(ji,jk,jn) * aplus /dxplus 407 ! 408 za(ji,jk) = -rminus 409 zb(ji,jk) = 1. + rminus + rplus 410 zc(ji,jk) = -rplus 411 END DO 412 END DO 413 414 DO ji = 1, jpoce 415 aminus = ( ( volw3d(ji,jpksed-1) / dz3d(ji,jpksed-1) ) + & 416 & ( volw3d(ji,jpksed) / dz3d(ji,jpksed) ) ) / 2. 417 dxminus = ( dz3d(ji,jpksed-1) + dz3d(ji,jpksed) ) / 2. 418 rminus = ( dtsed_in / volw3d(ji,jpksed) ) * diff(ji,jpksed-1,jn) * aminus/ dxminus 419 ! 420 za(ji,jpksed) = -rminus 421 zb(ji,jpksed) = 1. + rminus 422 zc(ji,jpksed) = 0. 423 END DO 424 425 426 ! solves tridiagonal system of linear equations 427 ! ----------------------------------------------- 428 429 zr (:,:) = pwcp(:,:,jn) + psms(:,:) 430 zb (:,:) = zb(:,:) - preac(:,:) 431 zbet(: ) = zb(:,1) 432 pwcp(:,1,jn) = zr(:,1) / zbet(:) 433 434 ! 435 DO jk = 2, jpksed 436 DO ji = 1, jpoce 437 zgamm(ji,jk) = zc(ji,jk-1) / zbet(ji) 438 zbet(ji) = zb(ji,jk) - za(ji,jk) * zgamm(ji,jk) 439 pwcp(ji,jk,jn) = ( zr(ji,jk) - za(ji,jk) * pwcp(ji,jk-1,jn) ) / zbet(ji) 440 END DO 441 ENDDO 442 ! 443 DO jk = jpksed - 1, 1, -1 444 DO ji = 1, jpoce 445 pwcp(ji,jk,jn) = pwcp(ji,jk,jn) - zgamm(ji,jk+1) * pwcp(ji,jk+1,jn) 446 END DO 447 ENDDO 448 449 IF( ln_timing ) CALL timing_stop('sed_mat_dsri') 450 451 452 END SUBROUTINE sed_mat_dsri 453 454 265 455 END MODULE sedmat -
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/sedrst.F90
r14239 r15075 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 AGRIF45 44 !!---------------------------------------------------------------------- 46 45 ! … … 65 64 66 65 IF( .NOT. ln_rst_list .AND. nn_stock == -1 ) RETURN ! we will never do any restart 67 68 66 ! 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 ) ) THEN67 ! we open and define the tracer restart file one tracer time step before writing the data (-> at nitrst - 2*nn_dttrc + 1) 68 ! 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 69 IF( kt == nitrst - 1 .OR. nn_stock == 1 .OR. ( kt == nitend - 1 .AND. .NOT. lrst_sed ) ) THEN 72 70 ! beware of the format used to write kt (default is i8.8, that should be large enough) 73 71 IF( nitrst > 1.0e9 ) THEN ; WRITE(clkt,* ) nitrst … … 81 79 IF(lwp) WRITE(numsed,*) & 82 80 ' 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 81 CALL iom_open( TRIM(clpath)//TRIM(clname), numrsw, ldwrt = .TRUE., kdlev = jpksed ) 100 82 lrst_sed = .TRUE. 101 83 ENDIF … … 204 186 & zdta2(1:jpi,1:jpj,1:jpksed), iarroce(1:jpoce) ) 205 187 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 188 IF( ln_timing ) CALL timing_stop('sed_rst_read') 216 189 … … 247 220 IF(lwp) WRITE(numsed,*) '~~~~~~~~~' 248 221 249 250 trcsedi(:,:,:,:) = 0.0251 flxsedi3d(:,:,:,:) = 0.0252 222 zdta(:,:) = 1.0 253 223 zdta2(:,:,:) = 0.0 254 255 224 256 225 !! 1. WRITE in nutwrs 257 226 !! ------------------ 258 !zinfo(1) = REAL( kt)259 CALL iom_rstput( kt, nitrst, numrsw, 'kt', REAL( kt , wp))227 zinfo(1) = REAL( kt) 228 CALL iom_rstput( kt, nitrst, numrsw, 'kt', zinfo ) 260 229 261 230 ! Back to 2D geometry 262 231 DO jn = 1, jpsol 263 CALL unpack_arr( jpoce, trcsedi(1:jpi,1:jpj,1:jpksed,jn) , iarroce(1:jpoce), &232 CALL unpack_arr( jpoce, zdta2(1:jpi,1:jpj,1:jpksed) , iarroce(1:jpoce), & 264 233 & solcp(1:jpoce,1:jpksed,jn ) ) 234 cltra = TRIM(sedtrcd(jn)) 235 CALL iom_rstput( kt, nitrst, numrsw, TRIM(cltra), zdta2(:,:,:) ) 265 236 END DO 266 237 267 238 DO jn = 1, jpwat 268 CALL unpack_arr( jpoce, trcsedi(1:jpi,1:jpj,1:jpksed,jpsol+jn) , iarroce(1:jpoce), &239 CALL unpack_arr( jpoce, zdta2(1:jpi,1:jpj,1:jpksed) , iarroce(1:jpoce), & 269 240 & pwcp(1:jpoce,1:jpksed,jn ) ) 241 cltra = TRIM(sedtrcd(jpsol+jn)) 242 CALL iom_rstput( kt, nitrst, numrsw, TRIM(cltra), zdta2(:,:,:) ) 270 243 END DO 271 244 ! pH … … 276 249 ENDDO 277 250 278 CALL unpack_arr( jpoce, flxsedi3d(1:jpi,1:jpj,1:jpksed,1) , iarroce(1:jpoce), &251 CALL unpack_arr( jpoce, zdta2(1:jpi,1:jpj,1:jpksed) , iarroce(1:jpoce), & 279 252 & zdta(1:jpoce,1:jpksed) ) 253 cltra = TRIM(seddia3d(1)) 254 CALL iom_rstput( kt, nitrst, numrsw, TRIM(cltra), zdta2(:,:,:) ) 280 255 281 CALL unpack_arr( jpoce, flxsedi3d(1:jpi,1:jpj,1:jpksed,2) , iarroce(1:jpoce), &256 CALL unpack_arr( jpoce, zdta2(1:jpi,1:jpj,1:jpksed) , iarroce(1:jpoce), & 282 257 & co3por(1:jpoce,1:jpksed) ) 258 cltra = TRIM(seddia3d(2)) 259 CALL iom_rstput( kt, nitrst, numrsw, TRIM(cltra), zdta2(:,:,:) ) 283 260 284 261 ! prognostic variables 285 262 ! -------------------- 286 263 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 264 CALL unpack_arr( jpoce, zdta2(1:jpi,1:jpj,1:jpksed) , iarroce(1:jpoce), & 298 265 & db(1:jpoce,1:jpksed) ) … … 307 274 CALL iom_rstput( kt, nitrst, numrsw, TRIM(cltra), zdta2(:,:,:) ) 308 275 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 276 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 277 CALL iom_close( numrsw ) ! close the restart file (only at last time step) 323 278 IF( l_offline .AND. ln_rst_list ) THEN 324 279 nrst_lst = nrst_lst + 1 … … 351 306 !! In both those options, the exact duration of the experiment 352 307 !! 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.308 !! is not stored in the restart and is assumed to be (nittrc000-1)*rdt. 354 309 !! This is valid is the time step has remained constant. 355 310 !! … … 363 318 REAL(wp) :: zkt, zrdttrc1 364 319 REAL(wp) :: zndastp 365 CHARACTER(len = 82) :: clpname366 320 367 321 ! Time domain : restart … … 375 329 376 330 IF( ln_rst_sed ) THEN 377 lxios_sini = .FALSE.378 331 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 332 CALL iom_get ( numrsr, 'kt', zkt ) ! last time-step of previous run 333 391 334 IF(lwp) THEN 392 335 WRITE(numsed,*) ' *** Info read in restart : ' … … 401 344 ENDIF 402 345 ! Control of date 403 IF( nittrc000 - NINT( zkt ) /= nn_dtsed.AND. nn_rstsed /= 0 ) &346 IF( nittrc000 - NINT( zkt ) /= 1 .AND. nn_rstsed /= 0 ) & 404 347 & CALL ctl_stop( ' ===>>>> : problem with nittrc000 for the restart', & 405 348 & ' verify the restart file or rerun with nn_rsttr = 0 (namelist)' ) … … 414 357 ELSE 415 358 ndastp = ndate0 - 1 ! ndate0 read in the namelist in dom_nam 416 adatrj = ( REAL( nittrc000-1, wp ) * r n_Dt ) / rday359 adatrj = ( REAL( nittrc000-1, wp ) * rdt ) / rday 417 360 ! note this is wrong if time step has changed during run 418 361 ENDIF … … 435 378 IF(lwp) WRITE(numsed,*) 'trc_wri : write the TOP restart file (NetCDF) at it= ', kt, ' date= ', ndastp 436 379 IF(lwp) WRITE(numsed,*) '~~~~~~~' 437 IF( lwxios ) CALL iom_init_closedef(cw_sedrst_cxt)438 380 ENDIF 439 381 CALL iom_rstput( kt, nitrst, numrsw, 'kt' , REAL( kt , wp) ) ! time-step 440 382 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 ! 383 CALL iom_rstput( kt, nitrst, numrsw, 'adatrj' , adatrj ) ! number of elapsed days since 384 ! ! the begining of the run [s] 443 385 ENDIF 444 386 -
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/sedsfc.F90
r13295 r15075 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 51 -
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/sedsol.F90
r14385 r15075 10 10 USE sed_oce 11 11 USE sedini 12 USE sed diff12 USE sedfunc 13 13 USE seddsr 14 USE sedinorg 14 USE sedjac 15 USE sedbtb 16 USE sedco3 17 USE sedmat 18 USE TROS_F90 15 19 USE lib_mpp ! distribued memory computing library 16 20 USE lib_fortran … … 20 24 21 25 PUBLIC sed_sol 22 23 !! * Module variables24 26 25 27 !! $Id: sedsol.F90 5215 2015-04-15 16:11:56Z nicolasmartin $ … … 53 55 INTEGER, INTENT(in) :: kt 54 56 ! --- local variables 55 INTEGER :: ji, jk, js, jw, jnt ! dummy looop indices 56 REAL(wp) :: zadsnh4 57 INTEGER :: ji, jk, js, jw, jn, neq ! dummy looop indices 58 REAL(wp), DIMENSION( jpvode * jpksed ) :: ZXIN 59 REAL(wp), DIMENSION(jpoce,jpksed) :: preac 60 INTEGER :: JINDEX, ITOL, IJAC, MLJAC, IMAX 61 INTEGER :: MUJAC,LE1, LJAC, LIWORK, LWORK 62 INTEGER :: IDID, NMAXSTP, ROSM 63 REAL(wp) :: X, XEND, H, STEPMIN 64 REAL(wp), DIMENSION(1) :: RTOL, ATOL 65 REAL(wp), POINTER :: WORK(:) 66 INTEGER , POINTER :: IWORK(:) 67 INTEGER, DIMENSION(3) :: ISTAT 68 REAL(wp), DIMENSION(2) :: RSTAT 57 69 !! 58 70 !!---------------------------------------------------------------------- … … 67 79 ! ! 68 80 dens_mol_wgt(1:jpsol) = denssol / mol_wgt(1:jpsol) 81 stepros(:) = dtsed 69 82 ! 70 83 ENDIF 71 72 dtsed2 = dtsed / REAL( nrseddt, wp )73 84 74 85 ! 1. Change of geometry … … 86 97 ! 2. Change of previous solid fractions (due to volum changes) for k=2 87 98 !--------------------------------------------------------------------- 88 89 DO js = 1, jpsol 90 DO ji = 1, jpoce 91 solcp(ji,2,js) = solcp(ji,2,js) * dz(2) / dz3d(ji,2) 92 ENDDO 99 DO js = 1, jpsol 100 solcp(:,2,js) = solcp(:,2,js) * dz(2) / dz3d(:,2) 93 101 END DO 94 102 … … 114 122 !------------------------------------------------------------- 115 123 DO jw = 1, jpwat 116 DO ji = 1, jpoce 117 pwcp(ji,1,jw) = pwcp(ji,1,jw) - & 118 & pwcp(ji,2,jw) * dzdep(ji) * por(2) / ( dzkbot(ji) + rtrn ) 119 END DO 120 ENDDO 121 122 zadsnh4 = 1.0 / ( 1.0 + adsnh4 ) 124 pwcp(:,1,jw) = pwcp(:,1,jw) - pwcp(:,2,jw) * dzdep(:) * por(2) / dzkbot(:) 125 ENDDO 123 126 124 127 ! -------------------------------------------------- 125 128 ! Computation of the diffusivities 126 129 ! -------------------------------------------------- 127 128 130 DO js = 1, jpwat 129 131 DO jk = 1, jpksed 130 DO ji = 1, jpoce 131 diff(ji,jk,js) = ( diff1(js) + diff2(js) * temp(ji) ) / ( 1.0 - 2.0 * log( por(jk) ) ) 132 END DO 132 diff(:,jk,js) = ( diff1(js) + diff2(js) * temp(:) ) / ( 1.0 - 2.0 * log( por(jk) ) ) 133 133 END DO 134 134 END DO … … 136 136 ! Impact of bioirrigation and adsorption on diffusion 137 137 ! --------------------------------------------------- 138 139 diff(:,:,jwnh4) = diff(:,:,jwnh4) * ( 1.0 + irrig(:,:) ) * zadsnh4140 138 diff(:,:,jwsil) = diff(:,:,jwsil) * ( 1.0 + irrig(:,:) ) 141 139 diff(:,:,jwoxy) = diff(:,:,jwoxy) * ( 1.0 + irrig(:,:) ) … … 146 144 diff(:,:,jwh2s) = diff(:,:,jwh2s) * ( 1.0 + irrig(:,:) ) 147 145 diff(:,:,jwso4) = diff(:,:,jwso4) * ( 1.0 + irrig(:,:) ) 148 diff(:,:,jwfe2) = diff(:,:,jwfe2) * ( 1.0 + 0.2 * irrig(:,:) ) 149 150 DO jnt = 1, nrseddt 151 CALL sed_diff( kt, jnt ) ! 1st pass in diffusion to get values at t+1/2 152 CALL sed_dsr ( kt, jnt ) ! Redox reactions 153 CALL sed_inorg( kt, jnt ) ! Inorganic reactions (dissolution) 154 CALL sed_diff ( kt, jnt ) ! 2nd pass in diffusion to get values at t+1 155 END DO 146 diff(:,:,jwlgw) = diff(:,:,jwlgw) * ( 1.0 + irrig(:,:) ) 147 DO jk = 1, jpksed 148 diff(:,jk,jwfe2) = diff(:,jk,jwfe2) * ( 1.0 + 0.1 * irrig(:,jk) ) * radsfe2(jk) 149 diff(:,jk,jwnh4) = diff(:,jk,jwnh4) * ( 1.0 + irrig(:,jk) ) * radsnh4(jk) 150 END DO 151 152 ! Conversion of volume units 153 !---------------------------- 154 DO js = 1, jpsol 155 DO jk = 1, jpksed 156 volc(:,jk,js) = ( vols3d(:,jk) * dens_mol_wgt(js) ) / & 157 & ( volw3d(:,jk) * 1.e-3 ) 158 ENDDO 159 ENDDO 160 161 ! Apply bioturbation and compute the impact of the slow SMS on species 162 CALL sed_btb( kt ) 163 ! Recompute pH after bioturbation and slow chemistry 164 CALL sed_co3( kt ) 165 166 ! The following part deals with the stiff ODEs 167 ! This is the expensive part of the code and should be carefully 168 ! chosen. We use the DVODE solver after many trials to find a cheap 169 ! way to solve the ODEs. This is not necessarily the most efficient 170 ! but this is the one that was not too much of a pain to code and that 171 ! was the most precise and quick. 172 ! The ones I tried : operator splitting (Strang), hybrid spectral methods 173 ! Brent, Powell's hybrid method, ... 174 ! --------------------------------------------------------------------- 175 ROSM = 2 176 NEQ = jpvode * jpksed 177 XEND = dtsed 178 RTOL = 0.005 179 ATOL = 0.005 180 ITOL = 0 181 IJAC = 1 182 MLJAC = jpvode 183 MUJAC = jpvode 184 LE1 = 2*MLJAC+MUJAC+1 185 LJAC = MLJAC+MUJAC+1 186 LIWORK = NEQ + 2 187 LWORK = NEQ*(LJAC+LE1+8) + 5 188 ALLOCATE(WORK(LWORK), IWORK(LIWORK) ) 189 NMAXSTP = 0 190 STEPMIN = dtsed 191 192 DO ji = 1, jpoce 193 X = 0.0 194 H = stepros(ji) 195 WORK = 0. 196 IWORK = 0 197 IWORK(2) = 6 198 199 ! Put all the species in one local array (nb of tracers * vertical 200 ! dimension 201 jindex = ji 202 DO jn = 1, NEQ 203 jk = jarr(jn,1) 204 js = jarr(jn,2) 205 IF (js <= jpwat) THEN 206 zxin(jn) = pwcp(ji,jk,js) * 1E6 207 ELSE 208 zxin(jn) = solcp(ji,jk,js-jpwat) * 1E6 209 ENDIF 210 END DO 211 212 ! Set options for VODE : banded matrix. SParse option is much more 213 ! expensive except if one computes the sparse Jacobian explicitly 214 ! To speed up the computation, one way is to reduce ATOL and RTOL 215 ! which may be a good option at the beginning of the simulations 216 ! during the spin up 217 ! ---------------------------------------------------------------- 218 CALL ROSK(ROSM, NEQ,JINDEX,X,zxin,XEND,H, & 219 RTOL,ATOL,ITOL, sed_jac,IJAC,MLJAC,MUJAC, & 220 WORK,LWORK,IWORK,LIWORK,IDID,ISTAT,RSTAT) 221 222 DO jn = 1, NEQ 223 jk = jarr(jn,1) 224 js = jarr(jn,2) 225 IF (js <= jpwat) THEN 226 pwcp(ji,jk,js) = zxin(jn) * 1E-6 227 ELSE 228 solcp(ji,jk,js-jpwat) = zxin(jn) * 1E-6 229 ENDIF 230 END DO 231 NMAXSTP = MAX(NMAXSTP,ISTAT(3)) 232 STEPMIN = MIN(STEPMIN, RSTAT(1) ) 233 stepros(ji) = RSTAT(1) 234 END DO 235 236 DEALLOCATE(WORK, IWORK) 237 238 preac(:,:) = 0. 239 CALL sed_mat_dsri( jwalk, preac, pwcpa(:,:,jwalk), dtsed ) 240 CALL sed_mat_dsri( jwpo4, preac, pwcpa(:,:,jwpo4), dtsed ) 156 241 157 242 !---------------------------------- … … 164 249 !-------------------------------------------------------------------- 165 250 DO jw = 1, jpwat 166 DO ji = 1, jpoce 167 pwcp(ji,1,jw) = pwcp(ji,1,jw) + & 168 & pwcp(ji,2,jw) * dzdep(ji) * por(2) / dzkbot(ji) 169 END DO 251 pwcp(:,1,jw) = pwcp(:,1,jw) + pwcp(:,2,jw) * dzdep(:) * por(2) / dzkbot(:) 170 252 ENDDO 171 253 … … 178 260 !------------------------------------------------------------------------ 179 261 DO js = 1, jpsol 180 DO ji = 1, jpoce 181 rainrg(ji,js) = raintg(ji) * solcp(ji,2,js) 182 rainrm(ji,js) = rainrg(ji,js) / mol_wgt(js) 183 END DO 262 rainrg(:,js) = raintg(:) * solcp(:,2,js) 263 rainrm(:,js) = rainrg(:,js) / mol_wgt(js) 184 264 ENDDO 185 265 … … 187 267 raintg(:) = 0. 188 268 DO js = 1, jpsol 189 DO ji = 1, jpoce 190 raintg(ji) = raintg(ji) + rainrg(ji,js) 191 END DO 269 raintg(:) = raintg(:) + rainrg(:,js) 192 270 ENDDO 193 271 … … 195 273 ! 3/ back to initial geometry 196 274 !-------------------------------- 197 DO ji = 1, jpoce 198 dz3d (ji,2) = dz(2) 199 volw3d(ji,2) = dz3d(ji,2) * por(2) 200 vols3d(ji,2) = dz3d(ji,2) * por1(2) 201 ENDDO 275 dz3d (:,2) = dz(2) 276 volw3d(:,2) = dz3d(:,2) * por(2) 277 vols3d(:,2) = dz3d(:,2) * por1(2) 202 278 203 279 IF( ln_timing ) CALL timing_stop('sed_sol') … … 205 281 END SUBROUTINE sed_sol 206 282 283 SUBROUTINE JACFAC 284 285 END SUBROUTINE JACFAC 286 207 287 END MODULE sedsol -
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/sedstp.F90
r14385 r15075 9 9 USE sedco3 ! carbonate in sediment pore water 10 10 USE sedsol ! Organic reactions and diffusion 11 USE sedbtb ! bioturbation12 11 USE sedadv ! vertical advection 13 12 USE sedsfc ! sediment surface data 14 13 USE sedrst ! restart 15 14 USE sedwri ! outputs 15 USE sedini 16 16 USE trcdmp_sed 17 17 USE lib_mpp ! distribued memory computing library … … 23 23 !! * Routine accessibility 24 24 PUBLIC sed_stp ! called by step.F90 25 26 !! * Substitutions 27 # include "do_loop_substitute.h90" 28 # include "domzgr_substitute.h90" 25 29 26 30 !! $Id$ … … 42 46 !! ! 06-04 (C. Ethe) Re-organization 43 47 !!---------------------------------------------------------------------- 44 INTEGER, INTENT(in) :: kt ! number of iteration 45 INTEGER, INTENT(in) :: Kbb, Kmm, Krhs ! time level indices 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 46 52 !!---------------------------------------------------------------------- 47 53 IF( ln_timing ) CALL timing_start('sed_stp') … … 53 59 54 60 dtsed = rDt_trc 55 ! dtsed2 = dtsed56 61 IF (kt /= nitsed000) THEN 57 CALL sed_dta( kt, Kbb, Kmm ) 62 CALL sed_dta( kt, Kbb, Kmm ) ! Load Data for bot. wat. Chem and fluxes 58 63 ENDIF 59 64 … … 63 68 ENDIF 64 69 65 CALL sed_btb( kt ) ! 1st pass of bioturbation at t+1/266 70 CALL sed_sol( kt ) ! Solute diffusion and reactions 67 CALL sed_btb( kt ) ! 2nd pass of bioturbation at t+168 71 CALL sed_adv( kt ) ! advection 69 72 CALL sed_co3( kt ) ! pH actualization for saving 70 IF (ln_sed_2way) CALL sed_sfc( kt, Kbb ) ! Give back new bottom wat chem to tracer model 73 74 IF (ln_sed_2way) CALL sed_sfc( kt, Kbb ) ! Give back new bottom wat chem to tracer model 71 75 ENDIF 72 76 CALL sed_wri( kt ) ! outputs 73 77 IF( kt == nitsed000 ) THEN 74 78 CALL iom_close( numrsr ) ! close input tracer restart file 75 IF(lrxios) CALL iom_context_finalize( cr_sedrst_cxt ) 76 ! IF(lwm) CALL FLUSH( numont ) ! flush namelist output 79 ! IF(lwm) CALL FLUSH( numont ) ! flush namelist output 77 80 ENDIF 78 81 IF( lrst_sed ) CALL sed_rst_wri( kt ) ! restart file output 79 82 80 IF( kt == nitsedend ) CLOSE( numsed )83 IF( kt == nitsedend ) CLOSE( numsed ) 81 84 82 IF( ln_timing ) CALL timing_stop('sed_stp')85 IF( ln_timing ) CALL timing_stop('sed_stp') 83 86 84 87 END SUBROUTINE sed_stp -
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/sedwri.F90
r14385 r15075 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,jptrased+1) )57 58 60 ! Initialize variables 59 61 ! -------------------- … … 90 92 91 93 CALL unpack_arr( jpoce, flxsedi3d(1:jpi,1:jpj,1:jpksed,3) , iarroce(1:jpoce), & 92 & sedligand(1:jpoce,1:jpksed) )93 94 CALL unpack_arr( jpoce, flxsedi3d(1:jpi,1:jpj,1:jpksed,4) , iarroce(1:jpoce), &95 94 & saturco3(1:jpoce,1:jpksed) ) 96 95 … … 129 128 CALL unpack_arr( jpoce, flxsedi2d(1:jpi,1:jpj,jpdia2dsed), iarroce(1:jpoce), zflx(1:jpoce,1) ) 130 129 130 ! 131 CALL lbc_lnk( 'sedwri', trcsedi(:,:,:,:), 'T', 1._wp ) 132 CALL lbc_lnk( 'sedwri', flxsedi3d(:,:,:,:), 'T', 1._wp ) 133 CALL lbc_lnk( 'sedwri', flxsedi2d(:,:,:), 'T', 1._wp ) 134 131 135 ! Start writing data 132 136 ! --------------------- … … 146 150 END DO 147 151 148 149 DEALLOCATE( zdta ) ; DEALLOCATE( zflx )150 151 152 IF( ln_timing ) CALL timing_stop('sed_wri') 152 153 -
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/trcini_pisces.F90
r14786 r15075 293 293 ! Initialization of the sediment model 294 294 IF( ln_sediment) & 295 & CALL sed_ini t! Initialization of the sediment model295 & CALL sed_ini ! Initialization of the sediment model 296 296 297 297 CALL p4z_sed_init ! loss of organic matter in the sediments -
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/trcstp.F90
r14086 r15075 102 102 CALL trc_wri ( kt, Kmm ) ! output of passive tracers with iom I/O manager 103 103 CALL trc_sms ( kt, Kbb, Kmm, Krhs ) ! tracers: sinks and sources 104 #if ! defined key_sed_off 104 105 CALL trc_trp ( kt, Kbb, Kmm, Krhs, Kaa ) ! transport of passive tracers 106 #endif 105 107 ! 106 108 ! Note passive tracers have been time-filtered in trc_trp but the time level
Note: See TracChangeset
for help on using the changeset viewer.
