Changeset 15297
- Timestamp:
- 2021-09-28T11:20:56+02:00 (20 months ago)
- Location:
- NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP
- Files:
-
- 1 added
- 22 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/par_sed.F90
r15127 r15297 19 19 20 20 INTEGER, PARAMETER :: jpdta = 18 21 22 ! Vertical sediment geometry 23 INTEGER, PUBLIC :: & 24 jpksed = 11 , & 25 jpksedm1 = 10 21 26 22 27 ! sediment tracer species … … 56 61 INTEGER, PARAMETER :: & 57 62 jptrased = jpsol + jpwat , & 58 jpvode = jptrased - 11 , &63 jpvode = jptrased - 11 , & 59 64 jpdia3dsed = 3 , & 60 jpdia2dsed = 2 165 jpdia2dsed = 22 61 66 62 67 END MODULE par_sed -
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/sed.F90
r15127 r15297 9 9 USE par_sed 10 10 USE oce_sed 11 USE trc, ONLY : profsed12 11 USE in_out_manager 13 12 … … 43 42 INTEGER , PUBLIC :: nitsed000 44 43 INTEGER , PUBLIC :: nitsedend 45 REAL(wp), PUBLIC :: denssol !: density of solid material46 44 LOGICAL , PUBLIC :: lrst_sed !: logical to control the trc restart write 47 45 LOGICAL , PUBLIC :: ln_rst_sed = .TRUE. !: initialisation from a restart file or not … … 56 54 CHARACTER(len = 80) , PUBLIC :: cn_sedrst_out !: suffix of pass. tracer restart name (output) 57 55 CHARACTER(len = 256), PUBLIC :: cn_sedrst_outdir !: restart output directory 56 INTEGER, PUBLIC :: nrosorder !: order of the rosenbrock method 57 REAL(wp), PUBLIC :: rosatol !: Tolerance for absolute error 58 REAL(wp), PUBLIC :: rosrtol !: Tolerance for relative error 58 59 59 60 ! 60 61 REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: pwcp, pwcpa 61 62 REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: solcp, solcpa 62 REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: Jacobian63 63 REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: diff 64 64 … … 66 66 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: pwcp_dta !: pore water data at given time-step 67 67 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: rainrm_dta !: rain data at at initial time 68 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: rainrm !: rain data at given time-step69 68 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: rainrg !: rain of each solid component in [g/(cm**2.s)] 70 69 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: fromsed !: 71 70 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: tosed !: 72 71 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: rearatpom !: 72 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: apluss, aminuss !: 73 73 ! 74 74 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: temp !: temperature 75 75 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: salt !: salinity 76 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: raintg !: total massic flux rained in each cell (sum of sol. comp.)77 76 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: fecratio !: Fe/C ratio in falling particles to the sediments 78 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: dzdep 77 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: dzdep, slatit, slongit !: total thickness of solid material rained [cm] in each cell 79 78 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: zkbot !: total thickness of solid material rained [cm] in each cell 80 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: wacc , wacc1!: total thickness of solid material rained [cm] in each cell79 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: wacc !: total thickness of solid material rained [cm] in each cell 81 80 ! 82 81 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: hipor !: [h+] in mol/kg*densSW … … 86 85 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: vols3d !: ??? 87 86 REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: volc 87 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: dens_sol !: Density of each solid fraction 88 88 89 89 … … 121 121 INTEGER , PUBLIC, DIMENSION(: ), ALLOCATABLE :: jsvode, isvode !: Computation of 1D array of sediments points 122 122 123 ! Vertical sediment geometry124 INTEGER, PUBLIC, SAVE :: jpksed125 INTEGER, PUBLIC, SAVE :: jpksedm1126 127 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: profsedw !: depth of middle of each layer128 129 123 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: epkbot !: ocean bottom layer thickness 130 124 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: gdepbot !: Depth of the sediment … … 137 131 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: db !: bioturbation ceofficient 138 132 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: irrig !: bioturbation ceofficient 139 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: radsnh4, radsfe2133 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: radssol, rads1sol 140 134 REAL(wp), PUBLIC, DIMENSION(:,: ), ALLOCATABLE :: saturco3 141 135 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: rdtsed !: sediment model time-step 142 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: stepros136 REAL(wp), PUBLIC, DIMENSION(: ), ALLOCATABLE :: rstepros !: Number of iteration of rosenbrock method 143 137 REAL(wp) :: dens !: density of solid material 144 138 !! Inputs / Outputs … … 163 157 !!------------------------------------------------------------------- 164 158 ! 165 ALLOCATE( profsed(jpksed) , profsedw(jpksed),trc_data(jpi,jpj,jpdta), &159 ALLOCATE( trc_data(jpi,jpj,jpdta) , & 166 160 & epkbot(jpi,jpj), gdepbot(jpi,jpj) , & 167 161 & dz(jpksed) , por(jpksed) , por1(jpksed) , & -
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/sedadv.F90
r15075 r15297 15 15 16 16 PUBLIC sed_adv 17 PUBLIC sed_adv_alloc18 19 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: dvolsp, dvolsm20 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: c2por, ckpor21 17 22 18 REAL(wp) :: cpor … … 49 45 ! * local variables 50 46 INTEGER :: ji, jk, js 51 INTEGER :: jn , ntimes, nztime, ikwneg47 INTEGER :: jn 52 48 53 49 REAL(wp), DIMENSION(jpoce,jpksed,jpsol+jpads) :: zsolcp 54 REAL(wp), DIMENSION(jpksed,jpsol+jpads) :: zsolcpno 55 REAL(wp), DIMENSION(jpksed) :: zfilled, zfull, zfromup, zempty 56 REAL(wp), DIMENSION(jpoce,jpksed) :: zgap, zwb 57 REAL(wp), DIMENSION(jpoce,jpsol+jpads) :: zrainrf 58 REAL(wp), DIMENSION(: ), ALLOCATABLE :: zraipush 59 60 REAL(wp) :: zkwnup, zkwnlo, zfrac, zfromce, zrest 50 REAL(wp) :: solfu, zfilled, zwb, fulsed, uebers, seddef 61 51 62 52 !------------------------------------------------------------------------ 63 64 53 65 54 IF( ln_timing ) CALL timing_start('sed_adv') … … 71 60 WRITE(numsed,*) ' ' 72 61 ENDIF 73 por1clay = denssol * por1(jpksed) * dz(jpksed)74 cpor = por1(jpksed) / por1(2)75 DO jk = 2, jpksed76 c2por(jk) = por1(2) / por1(jk)77 ckpor(jk) = por1(jpksed) / por1(jk)78 ENDDO79 DO jk = jpksedm1, 2, -180 dvolsp(jk) = vols(jk+1) / vols(jk)81 ENDDO82 DO jk = 3, jpksed83 dvolsm(jk) = vols(jk-1) / vols(jk)84 ENDDO85 62 ENDIF 86 63 … … 100 77 fromsed(:,:) = 0. 101 78 tosed (:,:) = 0. 102 ikwneg = 1103 nztime = jpksed104 79 105 ALLOCATE( zraipush(nztime) ) 106 80 solfu = 0.0 81 DO jk = 2, jpksed 82 solfu = solfu + dz(jk) * por1(jk) 83 END DO 107 84 ! Initiate gap 108 85 !-------------- 109 zgap(:,:) = 0.110 DO js = 1, jpsol111 zgap(:,:) = zgap(:,:) + zsolcp(:,:,js)112 ENDDO113 zgap(:,:) = 1. - zgap(:,:)114 86 115 87 ! Initiate burial rates 116 88 !----------------------- 117 zwb(:,:) = 0. 118 DO jk = 2, jpksed 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 131 DO ji = 1, jpoce 132 IF( raintg(ji) /= 0. ) zrainrf(ji,js) = rainrg(ji,js) / raintg(ji) 133 END DO 134 ENDDO 135 136 ! Computation of full and empty solid fraction in each layer 137 ! for all 'burial' case 138 !---------------------------------------------------------- 89 DO ji = 1, jpoce 90 DO jk = 2, jpksed-1 91 zfilled = 0._wp 92 DO js = 1, jpsol 93 zfilled = zfilled + zsolcp(ji,jk,js) / dens_sol(js) 94 END DO 95 zwb = MAX(0._wp, (zfilled - 1._wp) / (zfilled + rtrn) ) 139 96 140 97 98 DO js = 1, jpsol + jpads 99 uebers = zwb * zsolcp(ji,jk,js) 100 zsolcp(ji,jk,js) = zsolcp(ji,jk,js) - uebers 101 zsolcp(ji,jk+1,js) = zsolcp(ji,jk+1,js) + uebers * dz(jk) * por1(jk) / ( dz(jk+1) * por1(jk+1) ) 102 END DO 103 END DO 104 105 zfilled = 0._wp 106 DO js = 1, jpsol 107 zfilled = zfilled + zsolcp(ji,jpksed,js) / dens_sol(js) 108 END DO 109 zwb = MAX(0._wp, (zfilled - 1._wp) / (zfilled + rtrn) ) 110 DO js = 1, jpsol + jpads 111 uebers = zwb * zsolcp(ji,jpksed,js) 112 zsolcp(ji,jpksed,js) = zsolcp(ji,jpksed,js) - uebers 113 tosed(ji,js) = uebers * dz(jpksed) * por1(jpksed) 114 END DO 115 END DO 116 141 117 DO ji = 1, jpoce 118 fulsed = 0._wp 119 DO jk = 2, jpksed 120 zfilled = 0._wp 121 DO js = 1, jpsol 122 zfilled = zfilled + zsolcp(ji,jk,js) / dens_sol(js) 123 END DO 124 fulsed = fulsed + zfilled * dz(jk) * por1(jk) 125 END DO 142 126 143 ! computation of total weight fraction in sediment 144 !------------------------------------------------- 145 zfilled(:) = 0. 146 DO js = 1, jpsol 147 DO jk = 2, jpksed 148 zfilled(jk) = zfilled(jk) + zsolcp(ji,jk,js) 127 seddef = solfu - fulsed 128 129 zsolcp(ji,jpksed,jsclay) = zsolcp(ji,jpksed,jsclay) + seddef / ( por1(jpksed) * dz(jpksed) ) & 130 & * dens_sol(jsclay) 131 fromsed(ji,jsclay) = seddef * dens_sol(jsclay) 132 133 DO jk = jpksed, 3, -1 134 zfilled = 0._wp 135 DO js = 1, jpsol 136 zfilled = zfilled + zsolcp(ji,jk,js) / dens_sol(js) 149 137 END DO 150 ENDDO151 152 DO js = 1, jpsol + jpads153 DO jk = 2, jpksed154 zsolcp no(jk,js) = zsolcp(ji,jk,js) / zfilled(jk)138 zwb = MAX(0._wp, (zfilled - 1._wp) / (zfilled + rtrn) ) 139 DO js = 1, jpsol + jpads 140 uebers = zwb * zsolcp(ji,jk,js) 141 zsolcp(ji,jk,js) = zsolcp(ji,jk,js) - uebers 142 zsolcp(ji,jk-1,js) = zsolcp(ji,jk-1,js) + uebers * dz(jk) * por1(jk) / ( dz(jk-1) * por1(jk-1) ) 155 143 END DO 156 END DO144 END DO 157 145 158 ! burial 3 cases: 159 ! zwb > 0 ==> rain > total rection loss 160 ! zwb = 0 ==> rain = 0 161 ! zwb < 0 ==> rain > 0 and rain < total reaction loss 162 !---------------------------------------------------------------- 163 164 IF( zwb(ji,jpksed) > 0. ) THEN 165 166 zfull (jpksed) = zfilled(jpksed) 167 zempty(jpksed) = 1. - zfull(jpksed) 168 DO jk = jpksedm1, 2, -1 169 zfull (jk) = zfilled(jk) - zempty(jk+1) * dvolsp(jk) 170 zempty(jk) = 1. - zfull(jk) 171 ENDDO 172 173 ! Computation of solid sediment species 174 !-------------------------------------- 175 ! push entire sediment column downward to account rest of rain 176 DO js = 1, jpsol + jpads 177 DO jk = jpksed, 3, -1 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) 182 183 DO jk = 2, jpksed 184 zsolcpno(jk,js) = zsolcp(ji,jk,js) 185 END DO 186 ENDDO 187 188 zrest = zwb(ji,jpksed) * cpor 189 ! what is remaining is less than dz(2) 190 IF( zrest <= dz(2) ) THEN 191 192 zfromup(2) = zrest / dz(2) 193 DO jk = 3, jpksed 194 zfromup(jk) = zwb(ji,jpksed) * ckpor(jk) / dz(jk) 195 ENDDO 196 DO js = 1, jpsol + jpads 197 zfromce = 1. - zfromup(2) 198 zsolcp(ji,2,js) = zfromce * zsolcpno(2,js) + zfromup(2) * zrainrf(ji,js) 199 DO jk = 3, jpksed 200 zfromce = 1. - zfromup(jk) 201 zsolcp(ji,jk,js) = zfromce * zsolcpno(jk,js) + zfromup(jk) * zsolcpno(jk-1,js) 202 ENDDO 203 fromsed(ji,js) = 0. 204 ! quantities to push in deeper sediment 205 tosed (ji,js) = zsolcpno(jpksed,js) & 206 & * zwb(ji,jpksed) * denssol * por1(jpksed) 207 ENDDO 208 209 ELSE ! what is remaining is great than dz(2) 210 211 ntimes = INT( zrest / dz(2) ) + 1 212 IF( ntimes > nztime ) CALL ctl_stop( 'STOP', 'sed_adv : rest too large ' ) 213 zraipush(1) = dz(2) 214 zrest = zrest - zraipush(1) 215 DO jn = 2, ntimes 216 IF( zrest >= dz(2) ) THEN 217 zraipush(jn) = dz(2) 218 zrest = zrest - zraipush(jn) 219 ELSE 220 zraipush(jn) = zrest 221 zrest = 0. 222 ENDIF 223 ENDDO 224 225 DO jn = 1, ntimes 226 DO js = 1, jpsol + jpads 227 DO jk = 2, jpksed 228 zsolcpno(jk,js) = zsolcp(ji,jk,js) 229 END DO 230 ENDDO 231 232 zfromup(2) = zraipush(jn) / dz(2) 233 DO jk = 3, jpksed 234 zfromup(jk) = ( zraipush(jn) / dz(jk) ) * c2por(jk) 235 ENDDO 236 237 DO js = 1, jpsol + jpads 238 zfromce = 1. - zfromup(2) 239 zsolcp(ji,2,js) = zfromce * zsolcpno(2,js) + zfromup(2) * zrainrf(ji,js) 240 DO jk = 3, jpksed 241 zfromce = 1. - zfromup(jk) 242 zsolcp(ji,jk,js) = zfromce * zsolcpno(jk,js) + zfromup(jk) * zsolcpno(jk-1,js) 243 ENDDO 244 fromsed(ji,js) = 0. 245 tosed (ji,js) = tosed(ji,js) + zsolcpno(jpksed,js) * zraipush(jn) & 246 & * denssol * por1(2) 247 ENDDO 248 ENDDO 249 250 ENDIF 251 252 ELSE IF( raintg(ji) < eps ) THEN ! rain = 0 253 254 zfull (2) = zfilled(2) 255 zempty(2) = 1. - zfull(2) 256 DO jk = 3, jpksed 257 zfull (jk) = zfilled(jk) - zempty(jk-1) * dvolsm(jk) 258 zempty(jk) = 1. - zfull(jk) 259 ENDDO 260 261 ! fill boxes with weight fraction from underlying box 262 DO js = 1, jpsol + jpads 263 DO jk = 2, jpksedm1 264 zsolcp(ji,jk,js) = zfull(jk) * zsolcpno(jk,js) + zempty(jk) * zsolcpno(jk+1,js) 265 END DO 266 zsolcp(ji,jpksed,js) = zsolcpno(jpksed,js) * zfull(jpksed) 267 tosed (ji,js) = 0. 268 fromsed(ji,js) = 0. 269 ENDDO 270 ! for the last layer, one make go up clay 271 zsolcp(ji,jpksed,jsclay) = zsolcp(ji,jpksed,jsclay) + zempty(jpksed) * 1. 272 fromsed(ji,jsclay) = zempty(jpksed) * 1. * por1clay 273 ELSE ! rain > 0 and rain < total reaction loss 274 275 DO jk = 2, jpksed 276 zfull (jk) = zfilled(jk) 277 zempty(jk) = 1. - zfull(jk) 278 ENDDO 279 280 ! Determination of indice of layer - ikwneg - where advection is reversed 281 !------------------------------------------------------------------------ 282 iflag: DO jk = 2, jpksed 283 IF( zwb(ji,jk) < 0. ) THEN 284 ikwneg = jk 285 EXIT iflag 286 ENDIF 287 ENDDO iflag 288 289 ! computation of zfull and zempty 290 ! 3 cases : a/ ikwneg=2, b/ikwneg=3...jpksedm1, c/ikwneg=jpksed 291 !------------------------------------------------------------- 292 IF( ikwneg == 2 ) THEN ! advection is reversed in the first sediment layer 293 zkwnup = dtsed * raintg(ji) / ( denssol * por1(ikwneg) * dz(ikwneg) ) 294 zkwnlo = ABS( zwb(ji,ikwneg) ) / dz(ikwneg) 295 zfull (ikwneg+1) = zfilled(ikwneg+1) - zkwnlo * dvolsm(ikwneg+1) 296 zempty(ikwneg+1) = 1. - zfull(ikwneg+1) 297 DO jk = ikwneg+2, jpksed 298 zfull (jk) = zfilled(jk) - zempty(jk-1) * dvolsm(jk) 299 zempty(jk) = 1. - zfull(jk) 300 ENDDO 301 DO js = 1, jpsol + jpads 302 zsolcp(ji,2,js) = zfull(2) * zsolcpno(2,js)+ zkwnlo * zsolcpno(3,js) & 303 & + zkwnup * zrainrf(ji,js) 304 DO jk = 3, jpksedm1 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) 308 tosed(ji,js) = 0. 309 fromsed(ji,js) = 0. 310 ENDDO 311 zsolcp(ji,jpksed,jsclay) = zsolcp(ji,jpksed,jsclay) + zempty(jpksed) * 1. 312 fromsed(ji,jsclay) = zempty(jpksed) * 1. * por1clay 313 314 ELSE IF( ikwneg == jpksed ) THEN 315 316 zkwnup = ABS( zwb(ji,ikwneg-1) ) * dvolsm(ikwneg) / dz(ikwneg) 317 zkwnlo = ABS( zwb(ji,ikwneg) ) / dz(ikwneg) 318 zfull (ikwneg-1) = zfilled(ikwneg-1) - zkwnup * dvolsp(ikwneg-1) 319 zempty(ikwneg-1) = 1. - zfull(ikwneg-1) 320 DO jk = ikwneg-2, 2, -1 321 zfull (jk) = zfilled(jk) - zempty(jk+1) * dvolsp(jk) 322 zempty(jk) = 1. - zfull(jk) 323 ENDDO 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 328 DO jk = jpksedm1, 3, -1 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) & 332 & + zkwnup * zsolcpno(jpksedm1,js) 333 tosed(ji,js) = 0. 334 fromsed(ji,js) = 0. 335 ENDDO 336 zsolcp(ji,jpksed,jsclay) = zsolcp(ji,jpksed,jsclay) + zkwnlo * 1. 337 fromsed(ji,jsclay) = zkwnlo * 1.* por1clay 338 ELSE IF( ikwneg > 2 .AND. ikwneg < jpksed-1 ) THEN 339 340 zkwnup = ABS( zwb(ji,ikwneg-1) ) * por1(ikwneg-1) / ( dz(ikwneg) * por1(ikwneg) ) 341 zkwnlo = ABS( zwb(ji,ikwneg) ) / dz(ikwneg) 342 343 IF( ikwneg > 3 ) THEN 344 345 zfull (ikwneg-1) = zfilled(ikwneg-1) - zkwnup * dvolsp(ikwneg-1) 346 zempty(ikwneg-1) = 1. - zfull(ikwneg-1) 347 DO jk = ikwneg-2, 2, -1 348 zfull (jk) = zfilled(jk) - zempty(jk+1) * dvolsp(jk) 349 zempty(jk) = 1. - zfull(jk) 350 ENDDO 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 355 DO jk = ikwneg-1, 3, -1 356 zsolcp(ji,jk,js) = zfull(jk) * zsolcpno(jk,js) + zempty(jk) * zsolcpno(jk-1,js) 357 ENDDO 358 ENDDO 359 ELSE ! ikw = 3 360 361 zfull (2) = zfilled(2) - zkwnup * dvolsm(3) 362 zempty(2) = 1. - zfull(2) 363 DO js = 1, jpsol + jpads 364 zsolcp(ji,2,js) = zfull(2) * zsolcpno(2,js) + zempty(2) * zrainrf(ji,js) 365 ENDDO 366 ENDIF 367 368 IF( ikwneg < jpksedm1) THEN 369 370 zfull (ikwneg+1) = zfilled(ikwneg+1) - zkwnlo * dvolsm(ikwneg+1) 371 zempty(ikwneg+1) = 1. - zfull(ikwneg+1) 372 DO jk = ikwneg+2, jpksed 373 zfull (jk) = zfilled(jk) - zempty(jk-1) * dvolsm(jk) 374 zempty(jk) = 1. - zfull(jk) 375 ENDDO 376 DO js = 1, jpsol + jpads 377 DO jk = ikwneg+1, jpksedm1 378 zsolcp(ji,jk,js) = zfull(jk) * zsolcpno(jk,js) + zempty(jk) * zsolcpno(jk+1,js) 379 ENDDO 380 zsolcp(ji,jpksed,js) = zfull(jpksed) * zsolcpno(jpksed,js) 381 ENDDO 382 zsolcp(ji,jpksed,jsclay) = zsolcp(ji,jpksed,jsclay) + zempty(jpksed) * 1. 383 ELSE 384 385 zfull (jpksed) = zfilled(jpksed) - zkwnlo * dvolsm(jpksed) 386 zempty(jpksed) = 1. - zfull(jpksed) 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. 391 ENDIF ! jpksedm1 392 393 ! ikwneg = jpksedm1 ; ikwneg+1 = jpksed ; ikwneg-1 = jpksed - 2 394 DO js = 1, jpsol + jpads 395 zsolcp(ji,ikwneg,js) = zfull(ikwneg) * zsolcpno(ikwneg ,js) & 396 & + zkwnup * zsolcpno(ikwneg-1,js) & 397 & + zkwnlo * zsolcpno(ikwneg+1,js) 398 tosed (ji,js) = 0. 399 fromsed(ji,js) = 0. 400 ENDDO 401 fromsed(ji,jsclay) = zempty(jpksed) * 1. * por1clay 402 403 ENDIF ! ikwneg(ji) = 2 404 ENDIF ! zwb > 0 405 ENDDO ! ji = 1, jpoce 146 END DO 406 147 407 148 DO js = 1, jpsol … … 409 150 END DO 410 151 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) 152 pwcp(:,jk,jwnh4) = (pwcp(:,jk,jwnh4) + zsolcp(:,jk,jpsol+1) * por1(jk) / por(jk) ) * radssol(jk,jwnh4) 153 IF (jpoce == 146 .and. slatit(145) > 0.) write(0,*) 'plante advection ',pwcp(145,jk,jwfe2)*1E6,zsolcp(145,jk,jpsol+2)*1E6, & 154 & (pwcp(145,jk,jwfe2) + zsolcp(145,jk,jpsol+2) * por1(jk) / por(jk) ) * radssol(jk,jwfe2) * 1E6 155 pwcp(:,jk,jwfe2) = (pwcp(:,jk,jwfe2) + zsolcp(:,jk,jpsol+2) * por1(jk) / por(jk) ) * radssol(jk,jwfe2) 413 156 END DO 414 157 415 rainrm(:,:) = 0.416 158 rainrg(:,:) = 0. 417 raintg(:) = 0.418 419 DEALLOCATE( zraipush )420 159 421 160 IF( ln_timing ) CALL timing_stop('sed_adv') … … 424 163 425 164 426 INTEGER FUNCTION sed_adv_alloc()427 !!----------------------------------------------------------------------428 !! *** ROUTINE p4z_prod_alloc ***429 !!----------------------------------------------------------------------430 ALLOCATE( dvolsp(jpksed), dvolsm(jpksed), c2por(jpksed), &431 & ckpor(jpksed) , STAT = sed_adv_alloc )432 !433 IF( sed_adv_alloc /= 0 ) CALL ctl_stop( 'STOP', 'sed_adv_alloc : failed to allocate arrays.' )434 !435 END FUNCTION sed_adv_alloc436 437 165 END MODULE sedadv -
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/sedarr.F90
r15075 r15297 51 51 jjd = ( tab_ind(jn) - 1 ) / jpi + 1 52 52 tab1d(jn) = tab2d(jid, jjd) 53 ! IF (mig(jid) == 150 .and. mjg(jjd) == 136) write(0,*) 'plante indices ',jn,ndim1d,slatit(jn),slongit(jn) 53 54 END DO 54 55 -
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/sedbtb.F90
r15127 r15297 6 6 !! * Modules used 7 7 USE sed ! sediment global variable 8 USE sed_oce 8 9 USE sedmat ! linear system of equations 9 10 USE sedinorg … … 38 39 ! * local variables 39 40 INTEGER :: ji, jk, js 40 REAL(wp), DIMENSION(jpoce,jpksed m1,jpsol+2) :: zsol,zrearat ! solution41 REAL(wp), DIMENSION(jpoce,jpksed,jpsol) :: zrearat ! solution 41 42 REAL(wp) :: zsolid1, zsolid2, zsolid3, zsumtot, zlimo2 42 43 !------------------------------------------------------------------------ … … 52 53 ! Initializations 53 54 !---------------- 54 zsol(:,:,:) = 0.55 55 zrearat = 0. 56 56 … … 62 62 ! right hand side of coefficient matrix 63 63 !-------------------------------------- 64 DO js = 1, jpsol 65 DO jk = 1, jpksedm1 66 zsol(:,jk,js) = solcp(:,jk+1,js) 67 END DO 68 END DO 69 70 DO jk = 1, jpksedm1 71 zsol(:,jk,jpsol+1) = pwcp(:,jk+1,jwnh4) * adsnh4 72 zsol(:,jk,jpsol+2) = pwcp(:,jk+1,jwfe2) * adsfe2 73 END DO 74 75 CALL sed_mat_btb( jpksedm1, jpsol+2, zsol, zrearat(:,:,:), dtsed ) 76 77 ! store solution of the tridiagonal system 78 !------------------------ 79 DO js = 1, jpsol 80 DO jk = 1, jpksedm1 81 solcp(:,jk+1,js) = zsol(:,jk,js) 82 END DO 83 ENDDO 84 85 ! NH4 86 DO jk = 1, jpksedm1 87 pwcp(:,jk+1,jwnh4) = (pwcp(:,jk+1,jwnh4) + zsol(:,jk,jpsol+1) * por1(jk+1) / por(jk+1) ) * radsnh4(jk+1) 88 END DO 89 90 ! Fe2 91 DO jk = 1, jpksedm1 92 pwcp(:,jk+1,jwfe2) = (pwcp(:,jk+1,jwfe2) + zsol(:,jk,jpsol+2) * por1(jk+1) / por(jk+1) ) * radsfe2(jk+1) 93 END DO 64 CALL sed_mat_btbi( jpsol, solcp, zrearat(:,:,:), dtsed ) 94 65 95 66 DO ji = 1, jpoce … … 100 71 zsolid2 = volc(ji,jk,jspos) * solcp(ji,jk,jspos) 101 72 zsolid3 = volc(ji,jk,jspor) * solcp(ji,jk,jspor) 102 rearatpom(ji,jk) = ( reac_pocl * zsolid1 + reac_pocs * zsolid2 + reac_pocr * zsolid3 ) * dtsed103 zsumtot = zsumtot + rearatpom(ji,jk) / dtsed* volw3d(ji,jk) * 1.e-3 * 86400. * 365. * 1E373 rearatpom(ji,jk) = ( reac_pocl * zsolid1 + reac_pocs * zsolid2 + reac_pocr * zsolid3 ) 74 zsumtot = zsumtot + rearatpom(ji,jk) * volw3d(ji,jk) * 1.e-3 * 86400. * 365. * 1E3 104 75 END DO 105 76 … … 107 78 ! This parameterization is taken from Archer et al. (2002) 108 79 ! -------------------------------------------------------------- 109 zlimo2 = pwcp(ji,1,jwoxy) / (pwcp(ji,1,jwoxy) + 20.E-6)80 zlimo2 = max(0.01, pwcp(ji,1,jwoxy) / (pwcp(ji,1,jwoxy) + 20.E-6) ) 110 81 db(ji,:) = dbiot * zsumtot**0.85 * zlimo2 / (365.0 * 86400.0) 111 82 … … 115 86 IF (ln_btbz) THEN 116 87 DO jk = 1, jpksed 117 IF (profsedw(jk) < 2.0 * dbtbzsc) THEN118 88 db(ji,jk) = db(ji,jk) * exp( -(profsedw(jk) / dbtbzsc)**2 ) 119 ELSE120 db(ji,jk) = 0.121 ENDIF122 89 END DO 123 90 ELSE … … 142 109 ! --------------------------------------------------------- 143 110 CALL sed_inorg( kt ) 144 CALL sed_org( kt )145 111 146 112 IF( ln_timing ) CALL timing_stop('sed_btb') -
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/sedchem.F90
r15075 r15297 149 149 CALL sed_chem_cst 150 150 ELSE 151 DO_2D( 1, 1, 1, 1)151 DO_2D( 0, 0, 0, 0 ) 152 152 ikt = mbkt(ji,jj) 153 153 IF ( tmask(ji,jj,ikt) == 1 ) THEN -
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/seddsr.F90
r15127 r15297 7 7 !! * Modules used 8 8 USE sed ! sediment global variable 9 USE sed_oce 9 10 USE sedini 10 11 USE lib_mpp ! distribued memory computing library … … 23 24 CONTAINS 24 25 25 SUBROUTINE sed_dsr( ji )26 SUBROUTINE sed_dsr( accmask ) 26 27 !!---------------------------------------------------------------------- 27 28 !! *** ROUTINE sed_dsr *** … … 45 46 !!---------------------------------------------------------------------- 46 47 !! Arguments 47 INTEGER, INTENT(in) :: ji ! number of iteration48 48 ! --- local variables 49 INTEGER :: jk, js, jw, jn ! dummy looop indices 50 51 REAL(wp), DIMENSION(jpksed) ::zlimo2, zlimno3, zlimso4, zlimfeo ! undersaturation ; indice jpwatp1 is for calcite 49 INTEGER, DIMENSION(jpoce), INTENT(in) :: accmask 50 INTEGER :: ji, jk, js, jw, jn ! dummy looop indices 51 52 REAL(wp), DIMENSION(jpoce,jpksed) ::zlimo2, zlimno3, zlimso4, zlimfeo ! undersaturation ; indice jpwatp1 is for calcite 52 53 REAL(wp) :: zreasat 53 54 !! … … 59 60 !---------------------- 60 61 61 zlimo2 (: ) = 0. ; zlimno3(:) = 0.62 zlimso4(: ) = 0.62 zlimo2 (:,:) = 0. ; zlimno3(:,:) = 0. 63 zlimso4(:,:) = 0. 63 64 64 65 !---------------------------------------------------------- … … 73 74 ! -------------------------------------------------------------- 74 75 DO jk = 2, jpksed 75 zreasat = rearatpom(ji,jk) 76 ! For alkalinity 77 pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) + zreasat * ( srno3 - 2.* spo4r ) 78 ! For Phosphate (in mol/l) 79 pwcpa(ji,jk,jwpo4) = pwcpa(ji,jk,jwpo4) + zreasat * spo4r 80 81 ! For iron (in mol/l) 82 pwcpa(ji,jk,jwfe2) = pwcpa(ji,jk,jwfe2) + fecratio(ji) * zreasat * radsfe2(jk) 76 DO ji = 1, jpoce 77 IF (accmask(ji) == 0 ) THEN 78 zreasat = rearatpom(ji,jk) 79 ! For alkalinity 80 pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) + zreasat * ( srno3 - 2.* spo4r ) 81 ! For Phosphate (in mol/l) 82 pwcpa(ji,jk,jwpo4) = pwcpa(ji,jk,jwpo4) + zreasat * spo4r 83 ! For Ammonium 84 pwcpa(ji,jk,jwnh4) = pwcpa(ji,jk,jwnh4) + zreasat * srno3 * radssol(jk,jwnh4) 85 ! For iron (in mol/l) 86 pwcpa(ji,jk,jwfe2) = pwcpa(ji,jk,jwfe2) + fecratio(ji) * zreasat * radssol(jk,jwfe2) 87 ENDIF 88 END DO 83 89 ENDDO 84 90 85 91 ! Computes SMS of oxygen 86 92 DO jk = 2, jpksed 87 zlimo2(jk) = pwcp(ji,jk,jwoxy) / ( pwcp(ji,jk,jwoxy) + xksedo2 ) 88 zreasat = rearatpom(ji,jk) * zlimo2(jk) 89 ! Acid Silicic 90 pwcpa(ji,jk,jwoxy) = pwcpa(ji,jk,jwoxy) - so2ut * zreasat 93 DO ji = 1, jpoce 94 IF (accmask(ji) == 0 ) THEN 95 zlimo2(ji,jk) = pwcp(ji,jk,jwoxy) / ( pwcp(ji,jk,jwoxy) + xksedo2 ) 96 zlimo2(ji,jk) = MAX( 0., zlimo2(ji,jk) ) 97 zreasat = rearatpom(ji,jk) * zlimo2(ji,jk) 98 ! Acid Silicic 99 pwcpa(ji,jk,jwoxy) = pwcpa(ji,jk,jwoxy) - so2ut * zreasat 100 ENDIF 101 END DO 91 102 ENDDO 92 103 … … 95 106 !-------------------------------------------------------------------- 96 107 DO jk = 2, jpksed 97 zlimno3(jk) = ( 1.0 - zlimo2(jk) ) * pwcp(ji,jk,jwno3) / ( pwcp(ji,jk,jwno3) + xksedno3 ) 98 zreasat = rearatpom(ji,jk) * zlimno3(jk) 99 ! For nitrates 100 pwcpa(ji,jk,jwno3) = pwcpa(ji,jk,jwno3) - zreasat * srDnit 101 ! For alkalinity 102 pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) + zreasat * srDnit 108 DO ji = 1, jpoce 109 IF (accmask(ji) == 0 ) THEN 110 zlimno3(ji,jk) = ( 1.0 - zlimo2(ji,jk) ) * pwcp(ji,jk,jwno3) / ( pwcp(ji,jk,jwno3) + xksedno3 ) 111 zlimno3(ji,jk) = MAX(0., zlimno3(ji,jk) ) 112 zreasat = rearatpom(ji,jk) * zlimno3(ji,jk) * srDnit 113 ! For nitrates 114 pwcpa(ji,jk,jwno3) = pwcpa(ji,jk,jwno3) - zreasat 115 ! For alkalinity 116 pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) + zreasat 117 ENDIF 118 END DO 103 119 ENDDO 104 120 … … 109 125 110 126 DO jk = 2, jpksed 111 zlimfeo(jk) = ( 1.0 - zlimno3(jk) - zlimo2(jk) ) * solcp(ji,jk,jsfeo) / ( solcp(ji,jk,jsfeo) + xksedfeo ) 112 zreasat = rearatpom(ji,jk) * zlimfeo(jk) 113 ! For FEOH 114 solcpa(ji,jk,jsfeo) = solcpa(ji,jk,jsfeo) - 4.0 * zreasat / volc(ji,jk,jsfeo) 115 ! For PO4 116 pwcpa(ji,jk,jwpo4) = pwcpa(ji,jk,jwpo4) + zreasat * 4.0 * redfep 117 ! For alkalinity 118 pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) + 8.0 * zreasat 119 ! Iron 120 pwcpa(ji,jk,jwfe2) = pwcpa(ji,jk,jwfe2) + zreasat * 4.0 * radsfe2(jk) 127 DO ji = 1, jpoce 128 IF (accmask(ji) == 0 ) THEN 129 zlimfeo(ji,jk) = ( 1.0 - zlimno3(ji,jk) - zlimo2(ji,jk) ) * solcp(ji,jk,jsfeo) / ( solcp(ji,jk,jsfeo) + xksedfeo ) 130 zlimfeo(ji,jk) = MAX(0., zlimfeo(ji,jk) ) 131 zreasat = 4.0 * rearatpom(ji,jk) * zlimfeo(ji,jk) 132 ! For FEOH 133 solcpa(ji,jk,jsfeo) = solcpa(ji,jk,jsfeo) - zreasat / volc(ji,jk,jsfeo) 134 ! For PO4 135 pwcpa(ji,jk,jwpo4) = pwcpa(ji,jk,jwpo4) + zreasat * redfep 136 ! For alkalinity 137 pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) + 2.0 * zreasat 138 ! Iron 139 pwcpa(ji,jk,jwfe2) = pwcpa(ji,jk,jwfe2) + zreasat * radssol(jk,jwfe2) 140 ENDIF 141 END DO 121 142 ENDDO 122 143 … … 126 147 127 148 DO jk = 2, jpksed 128 zlimso4(jk) = ( 1.0 - zlimno3(jk) - zlimo2(jk) - zlimfeo(jk) ) * pwcp(ji,jk,jwso4) / ( pwcp(ji,jk,jwso4) + xksedso4 ) 129 zreasat = rearatpom(ji,jk) * zlimso4(jk) 130 ! For sulfur 131 pwcpa(ji,jk,jwso4) = pwcpa(ji,jk,jwso4) - 0.5 * zreasat 132 pwcpa(ji,jk,jwh2s) = pwcpa(ji,jk,jwh2s) + 0.5 * zreasat 133 ! For alkalinity 134 pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) + zreasat 149 DO ji = 1, jpoce 150 IF (accmask(ji) == 0 ) THEN 151 zlimso4(ji,jk) = ( 1.0 - zlimno3(ji,jk) - zlimo2(ji,jk) - zlimfeo(ji,jk) ) * pwcp(ji,jk,jwso4) / ( pwcp(ji,jk,jwso4) + xksedso4 ) 152 zlimso4(ji,jk) = MAX(0., zlimso4(ji,jk) ) 153 zreasat = 0.5 * rearatpom(ji,jk) * zlimso4(ji,jk) 154 ! For sulfur 155 pwcpa(ji,jk,jwso4) = pwcpa(ji,jk,jwso4) - zreasat 156 pwcpa(ji,jk,jwh2s) = pwcpa(ji,jk,jwh2s) + zreasat 157 ! For alkalinity 158 pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) + 2.0 * zreasat 159 ENDIF 160 END DO 135 161 ENDDO 136 162 … … 138 164 ! ------------------------- 139 165 140 call sed_dsr_redoxb( ji)166 call sed_dsr_redoxb( accmask ) 141 167 142 168 IF( ln_timing ) CALL timing_stop('sed_dsr') … … 144 170 END SUBROUTINE sed_dsr 145 171 146 SUBROUTINE sed_dsr_redoxb( ji)172 SUBROUTINE sed_dsr_redoxb( accmask ) 147 173 !!---------------------------------------------------------------------- 148 174 !! *** ROUTINE sed_dsr_redox *** … … 155 181 !! Arguments 156 182 ! --- local variables 157 INTEGER, INTENT(IN) :: ji158 INTEGER :: j ni, jnj, jk ! dummy looop indices183 INTEGER, DIMENSION(jpoce), INTENT(in) :: accmask 184 INTEGER :: ji, jni, jnj, jk ! dummy looop indices 159 185 160 186 REAL(wp) :: zalpha, zexcess, zh2seq, zsedfer, zreasat … … 165 191 166 192 DO jk = 2, jpksed 167 zalpha = ( pwcp(ji,jk,jwfe2) - pwcp(ji,jk,jwlgw) ) * 1E9 168 zsedfer = ( zalpha + SQRT( zalpha**2 + 1.E-10 ) ) /2.0 * 1E-9 169 170 ! First pass of the scheme. At the end, it is 1st order 171 ! ----------------------------------------------------- 172 ! Fe (both adsorbed and solute) + O2 173 zreasat = reac_fe2 * dtsed * zsedfer / radsfe2(jk) * pwcp(ji,jk,jwoxy) 174 pwcpa(ji,jk,jwfe2) = pwcpa(ji,jk,jwfe2) - zreasat * radsfe2(jk) 175 pwcpa(ji,jk,jwoxy) = pwcpa(ji,jk,jwoxy) - 0.25 * zreasat 176 pwcpa(ji,jk,jwpo4) = pwcpa(ji,jk,jwpo4) - redfep * zreasat 177 pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) - 2.0 * zreasat 178 solcpa(ji,jk,jsfeo) = solcpa(ji,jk,jsfeo) + zreasat / volc(ji,jk,jsfeo) 179 180 ! H2S + O2 181 zreasat = reac_h2s * dtsed * pwcp(ji,jk,jwoxy) * pwcp(ji,jk,jwh2s) 182 pwcpa(ji,jk,jwh2s) = pwcpa(ji,jk,jwh2s) - zreasat 183 pwcpa(ji,jk,jwoxy) = pwcpa(ji,jk,jwoxy) - 2.0 * zreasat 184 pwcpa(ji,jk,jwso4) = pwcpa(ji,jk,jwso4) + zreasat 185 pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) - 2.0 * zreasat 186 187 ! NH4 + O2 188 zreasat = reac_nh4 * dtsed * pwcp(ji,jk,jwnh4) / radsnh4(jk) * pwcp(ji,jk,jwoxy) 189 pwcpa(ji,jk,jwnh4) = pwcpa(ji,jk,jwnh4) - zreasat * radsnh4(jk) 190 pwcpa(ji,jk,jwoxy) = pwcpa(ji,jk,jwoxy) - 2.0 * zreasat 191 pwcpa(ji,jk,jwno3) = pwcpa(ji,jk,jwno3) + zreasat 192 pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) - 2.0 * zreasat 193 194 ! FeS - O2 195 zreasat = reac_feso * dtsed * solcp(ji,jk,jsfes) * volc(ji,jk,jsfes) * pwcp(ji,jk,jwoxy) 196 solcpa(ji,jk,jsfes) = solcpa(ji,jk,jsfes) - zreasat / volc(ji,jk,jsfes) 197 pwcpa(ji,jk,jwoxy) = pwcpa(ji,jk,jwoxy) - 2.0 * zreasat 198 pwcpa(ji,jk,jwfe2) = pwcpa(ji,jk,jwfe2) + zreasat * radsfe2(jk) 199 pwcpa(ji,jk,jwso4) = pwcpa(ji,jk,jwso4) + zreasat 200 201 ! FEOH + H2S 202 zreasat = reac_feh2s * dtsed * solcp(ji,jk,jsfeo) * volc(ji,jk,jsfeo) * pwcp(ji,jk,jwh2s) 203 solcpa(ji,jk,jsfeo) = solcpa(ji,jk,jsfeo) - 8.0 * zreasat / volc(ji,jk,jsfeo) 204 pwcpa(ji,jk,jwh2s) = pwcpa(ji,jk,jwh2s) - zreasat 205 pwcpa(ji,jk,jwfe2) = pwcpa(ji,jk,jwfe2) + 8.0 * zreasat * radsfe2(jk) 206 pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) + 14.0 * zreasat 207 pwcpa(ji,jk,jwso4) = pwcpa(ji,jk,jwso4) + zreasat 208 pwcpa(ji,jk,jwpo4) = pwcpa(ji,jk,jwpo4) + 8.0 * redfep * zreasat 209 210 ! Fe + H2S 211 zh2seq = MAX(rtrn, 2.1E-3 * hipor(ji,jk)) 212 zexcess = pwcp(ji,jk,jwh2s) * zsedfer / zh2seq - 1.0 213 IF ( zexcess >= 0.0 ) THEN 214 zreasat = reac_fesp * zexcess * dtsed 215 pwcpa (ji,jk,jwfe2) = pwcpa(ji,jk,jwfe2) - zreasat * radsfe2(jk) 216 solcpa(ji,jk,jsfes) = solcpa(ji,jk,jsfes) + zreasat / volc(ji,jk,jsfes) 217 pwcpa(ji,jk,jwh2s) = pwcpa(ji,jk,jwh2s) - zreasat 218 ELSE 219 zreasat = reac_fesd * zexcess * dtsed * solcp(ji,jk,jsfes) * volc(ji,jk,jsfes) 220 pwcpa (ji,jk,jwfe2) = pwcpa(ji,jk,jwfe2) - zreasat * radsfe2(jk) 221 solcpa(ji,jk,jsfes) = solcpa(ji,jk,jsfes) + zreasat / volc(ji,jk,jsfes) 222 pwcpa(ji,jk,jwh2s) = pwcpa(ji,jk,jwh2s) - zreasat 223 ENDIF 224 ! For alkalinity 225 pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) - zreasat * 2.0 193 DO ji = 1, jpoce 194 IF (accmask(ji) == 0 ) THEN 195 zalpha = ( pwcp(ji,jk,jwfe2) - pwcp(ji,jk,jwlgw) ) * 1E9 196 zsedfer = ( zalpha + SQRT( zalpha**2 + 1.E-10 ) ) /2.0 * 1E-9 197 198 ! First pass of the scheme. At the end, it is 1st order 199 ! ----------------------------------------------------- 200 ! Fe (both adsorbed and solute) + O2 201 zreasat = reac_fe2 * zsedfer / radssol(jk,jwfe2) * pwcp(ji,jk,jwoxy) 202 pwcpa(ji,jk,jwfe2) = pwcpa(ji,jk,jwfe2) - zreasat * radssol(jk,jwfe2) 203 pwcpa(ji,jk,jwoxy) = pwcpa(ji,jk,jwoxy) - 0.25 * zreasat 204 pwcpa(ji,jk,jwpo4) = pwcpa(ji,jk,jwpo4) - redfep * zreasat 205 pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) - 2.0 * zreasat 206 solcpa(ji,jk,jsfeo) = solcpa(ji,jk,jsfeo) + zreasat / volc(ji,jk,jsfeo) 207 208 ! H2S + O2 209 zreasat = reac_h2s * pwcp(ji,jk,jwoxy) * pwcp(ji,jk,jwh2s) 210 pwcpa(ji,jk,jwh2s) = pwcpa(ji,jk,jwh2s) - zreasat 211 pwcpa(ji,jk,jwoxy) = pwcpa(ji,jk,jwoxy) - 2.0 * zreasat 212 pwcpa(ji,jk,jwso4) = pwcpa(ji,jk,jwso4) + zreasat 213 pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) - 2.0 * zreasat 214 215 ! NH4 + O2 216 zreasat = reac_nh4 * pwcp(ji,jk,jwnh4) / radssol(jk,jwnh4) * pwcp(ji,jk,jwoxy) 217 pwcpa(ji,jk,jwnh4) = pwcpa(ji,jk,jwnh4) - zreasat * radssol(jk,jwnh4) 218 pwcpa(ji,jk,jwoxy) = pwcpa(ji,jk,jwoxy) - 2.0 * zreasat 219 pwcpa(ji,jk,jwno3) = pwcpa(ji,jk,jwno3) + zreasat 220 pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) - 2.0 * zreasat 221 222 ! FeS - O2 223 zreasat = reac_feso * solcp(ji,jk,jsfes) * volc(ji,jk,jsfes) * pwcp(ji,jk,jwoxy) 224 solcpa(ji,jk,jsfes) = solcpa(ji,jk,jsfes) - zreasat / volc(ji,jk,jsfes) 225 pwcpa(ji,jk,jwoxy) = pwcpa(ji,jk,jwoxy) - 2.0 * zreasat 226 pwcpa(ji,jk,jwfe2) = pwcpa(ji,jk,jwfe2) + zreasat * radssol(jk,jwfe2) 227 pwcpa(ji,jk,jwso4) = pwcpa(ji,jk,jwso4) + zreasat 228 229 ! FEOH + H2S 230 zreasat = reac_feh2s * solcp(ji,jk,jsfeo) * volc(ji,jk,jsfeo) * pwcp(ji,jk,jwh2s) 231 solcpa(ji,jk,jsfeo) = solcpa(ji,jk,jsfeo) - 8.0 * zreasat / volc(ji,jk,jsfeo) 232 pwcpa(ji,jk,jwh2s) = pwcpa(ji,jk,jwh2s) - zreasat 233 pwcpa(ji,jk,jwfe2) = pwcpa(ji,jk,jwfe2) + 8.0 * zreasat * radssol(jk,jwfe2) 234 pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) + 14.0 * zreasat 235 pwcpa(ji,jk,jwso4) = pwcpa(ji,jk,jwso4) + zreasat 236 pwcpa(ji,jk,jwpo4) = pwcpa(ji,jk,jwpo4) + 8.0 * redfep * zreasat 237 238 ! Fe + H2S 239 zh2seq = MAX(rtrn, 2.1E-3 * hipor(ji,jk)) 240 zexcess = pwcp(ji,jk,jwh2s) * zsedfer / zh2seq - 1.0 241 IF ( zexcess >= 0.0 ) THEN 242 zreasat = reac_fesp * zexcess 243 pwcpa (ji,jk,jwfe2) = pwcpa(ji,jk,jwfe2) - zreasat * radssol(jk,jwfe2) 244 solcpa(ji,jk,jsfes) = solcpa(ji,jk,jsfes) + zreasat / volc(ji,jk,jsfes) 245 pwcpa(ji,jk,jwh2s) = pwcpa(ji,jk,jwh2s) - zreasat 246 ELSE 247 zreasat = reac_fesd * zexcess * solcp(ji,jk,jsfes) * volc(ji,jk,jsfes) 248 pwcpa (ji,jk,jwfe2) = pwcpa(ji,jk,jwfe2) - zreasat * radssol(jk,jwfe2) 249 solcpa(ji,jk,jsfes) = solcpa(ji,jk,jsfes) + zreasat / volc(ji,jk,jsfes) 250 pwcpa(ji,jk,jwh2s) = pwcpa(ji,jk,jwh2s) - zreasat 251 ENDIF 252 ! For alkalinity 253 pwcpa(ji,jk,jwalk) = pwcpa(ji,jk,jwalk) - zreasat * 2.0 254 ENDIF 255 END DO 226 256 END DO 227 257 -
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/seddsrjac.F90
r15127 r15297 7 7 !! * Modules used 8 8 USE sed ! sediment global variable 9 USE sed_oce 9 10 USE sedini 10 11 USE seddsr … … 22 23 CONTAINS 23 24 24 SUBROUTINE sed_dsrjac( ji )25 SUBROUTINE sed_dsrjac( NEQ, NROWPD, jacvode, accmask ) 25 26 !!---------------------------------------------------------------------- 26 27 !! *** ROUTINE sed_dsr *** … … 44 45 !!---------------------------------------------------------------------- 45 46 !! Arguments 46 INTEGER, INTENT(in) :: ji ! number of iteration47 47 ! --- local variables 48 INTEGER :: jni, jnj, jk, js, jw, jn ! dummy looop indices 49 50 REAL(wp), DIMENSION(jpksed) ::zlimo2, zlimno3, zlimso4, zlimfeo ! undersaturation ; indice jpwatp1 is for calcite 51 REAL(wp) :: zvolw, zreasat, zlimtmpo2, zlimtmpno3 48 INTEGER, INTENT(in) :: NEQ, NROWPD 49 INTEGER, DIMENSION(jpoce), INTENT(in) :: accmask 50 REAL, DIMENSION(jpoce,NROWPD,NEQ), INTENT(inout) :: jacvode 51 INTEGER :: ji, jni, jnj, jnij, jk, js, jw, jn ! dummy looop indices 52 53 REAL(wp), DIMENSION(jpoce, jpksed) ::zlimo2, zlimno3, zlimso4, zlimfeo ! undersaturation ; indice jpwatp1 is for calcite 54 REAL(wp), DIMENSION(jpoce, jpksed) ::zlimdo2, zlimdno3, zlimdso4, zlimdfeo ! undersaturation ; indice jpwatp1 is for calcite 55 REAL(wp) :: zvolw, zreasat, zlimtmpo2, zlimtmpno3, zlimtmpfeo, zlimtmpso4 52 56 !! 53 57 !!---------------------------------------------------------------------- … … 58 62 !---------------------- 59 63 60 zlimo2 (: ) = 0. ; zlimno3(:) = 0.61 zlimso4(: ) = 0.64 zlimo2 (:,:) = 0. ; zlimno3(:,:) = 0. 65 zlimso4(:,:) = 0. 62 66 63 67 !---------------------------------------------------------- … … 67 71 ! Computes SMS of oxygen 68 72 DO jk = 2, jpksed 69 zlimo2(jk) = pwcp(ji,jk,jwoxy) / ( pwcp(ji,jk,jwoxy) + xksedo2 ) 70 ! Acid Silicic 71 jni = (jk-1)*jpvode+isvode(jwoxy) 72 Jacobian(ji,jni,jni) = Jacobian(ji,jni,jni) - so2ut * rearatpom(ji,jk) * xksedo2 & 73 & / ( pwcp(ji,jk,jwoxy) + xksedo2 )**2 73 DO ji = 1, jpoce 74 IF (accmask(ji) == 0) THEN 75 zlimo2(ji,jk) = pwcp(ji,jk,jwoxy) / ( pwcp(ji,jk,jwoxy) + xksedo2 ) 76 zlimdo2(ji,jk) = xksedo2 / ( pwcp(ji,jk,jwoxy) + xksedo2 )**2 77 ! Acid Silicic 78 jni = (jk-1)*jpvode+isvode(jwoxy) 79 jnij = jpvode + 1 80 jacvode(ji,jnij,jni) = jacvode(ji,jnij,jni) - so2ut * rearatpom(ji,jk) * zlimdo2(ji,jk) 81 ENDIF 82 END DO 74 83 ENDDO 75 84 … … 78 87 !-------------------------------------------------------------------- 79 88 DO jk = 2, jpksed 80 zlimno3(jk) = ( 1.0 - zlimo2(jk) ) * pwcp(ji,jk,jwno3) / ( pwcp(ji,jk,jwno3) + xksedno3 ) 81 ! For nitrates 82 jni = (jk-1)*jpvode+isvode(jwno3) 83 jnj = (jk-1)*jpvode+isvode(jwoxy) 84 Jacobian(ji,jni,jni) = Jacobian(ji,jni,jni) - srDnit * rearatpom(ji,jk) * (1.0 - zlimo2(jk) ) & 85 & * xksedno3 / ( pwcp(ji,jk,jwno3) + xksedno3 )**2 86 Jacobian(ji,jni,jnj) = Jacobian(ji,jni,jnj) + srDnit * rearatpom(ji,jk) * pwcp(ji,jk,jwno3) & 87 & / ( pwcp(ji,jk,jwno3) + xksedno3 ) * xksedo2 / ( pwcp(ji,jk,jwoxy) + xksedo2 )**2 89 DO ji = 1, jpoce 90 IF (accmask(ji) == 0) THEN 91 92 zlimno3(ji,jk) = pwcp(ji,jk,jwno3) / ( pwcp(ji,jk,jwno3) + xksedno3 ) 93 zlimdno3(ji,jk) = xksedno3 / ( pwcp(ji,jk,jwno3) + xksedno3 )**2 94 ! For nitrates 95 jni = (jk-1)*jpvode+isvode(jwno3) 96 jnj = (jk-1)*jpvode+isvode(jwoxy) 97 jnij = jpvode + 1 98 jacvode(ji,jnij,jni) = jacvode(ji,jnij,jni) - srDnit * rearatpom(ji,jk) * (1.0 - zlimo2(ji,jk) ) & 99 & * zlimdno3(ji,jk) 100 jnij = jni - jnj + jpvode + 1 101 jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) + srDnit * rearatpom(ji,jk) * zlimno3(ji,jk) * zlimdo2(ji,jk) 102 ENDIF 103 END DO 88 104 ENDDO 89 105 90 91 106 !-------------------------------------------------------------------- 92 107 ! Begining POC iron reduction … … 94 109 95 110 DO jk = 2, jpksed 96 zlimfeo(jk) = ( 1.0 - zlimno3(jk) - zlimo2(jk) ) * solcp(ji,jk,jsfeo) / ( solcp(ji,jk,jsfeo) + xksedfeo ) 97 ! For FEOH 98 jni = (jk-1)*jpvode+isvode(jpwat+jsfeo) 99 Jacobian(ji,jni,jni) = Jacobian(ji,jni,jni) - 4.0 * rearatpom(ji,jk) / volc(ji,jk,jsfeo) & 100 & * ( 1.0 - zlimno3(jk) - zlimo2(jk) ) * xksedfeo / ( solcp(ji,jk,jsfeo) + xksedfeo )**2 101 jnj = (jk-1)*jpvode+isvode(jwoxy) 102 zlimtmpo2 = solcp(ji,jk,jsfeo) / ( solcp(ji,jk,jsfeo) + xksedfeo ) * xksedo2 & 103 & / ( pwcp(ji,jk,jwoxy) + xksedo2 )**2 * ( 1.0 - pwcp(ji,jk,jwno3) / ( pwcp(ji,jk,jwno3) + xksedno3 ) ) 104 Jacobian(ji,jni,jnj) = Jacobian(ji,jni,jnj) + 4.0 * rearatpom(ji,jk) / volc(ji,jk,jsfeo) * zlimtmpo2 105 jnj = (jk-1)*jpvode+isvode(jwno3) 106 zlimtmpno3 = solcp(ji,jk,jsfeo) / ( solcp(ji,jk,jsfeo) + xksedfeo ) * xksedno3 & 107 & / ( pwcp(ji,jk,jwno3) + xksedno3 )**2 * ( 1.0 - zlimo2(jk) ) 108 Jacobian(ji,jni,jnj) = Jacobian(ji,jni,jnj) + 4.0 * rearatpom(ji,jk) / volc(ji,jk,jsfeo) * zlimtmpno3 109 ! Iron 110 zreasat = rearatpom(ji,jk) * 4.0 / ( 1.0 + adsfe2 * por1(jk) ) 111 jni = (jk-1)*jpvode+isvode(jwfe2) 112 jnj = (jk-1)*jpvode+isvode(jpwat+jsfeo) 113 Jacobian(ji,jni,jnj) = Jacobian(ji,jni,jnj) + zreasat & 114 & * ( 1.0 - zlimno3(jk) - zlimo2(jk) ) * xksedfeo / ( solcp(ji,jk,jsfeo) + xksedfeo )**2 115 jnj = (jk-1)*jpvode+isvode(jwoxy) 116 Jacobian(ji,jni,jnj) = Jacobian(ji,jni,jnj) - zreasat * zlimtmpo2 117 jnj = (jk-1)*jpvode+isvode(jwno3) 118 Jacobian(ji,jni,jnj) = Jacobian(ji,jni,jnj) - zreasat * zlimtmpno3 111 DO ji = 1, jpoce 112 IF (accmask(ji) == 0) THEN 113 zlimfeo(ji,jk) = solcp(ji,jk,jsfeo) / ( solcp(ji,jk,jsfeo) + xksedfeo ) 114 zlimdfeo(ji,jk) = xksedfeo / ( solcp(ji,jk,jsfeo) + xksedfeo )**2 115 ! For FEOH 116 jni = (jk-1)*jpvode+isvode(jpwat+jsfeo) 117 jnij = jpvode + 1 118 zlimtmpfeo = ( 1.0 - zlimno3(ji,jk) ) * ( 1.0 - zlimo2(ji,jk) ) * zlimdfeo(ji,jk) 119 jacvode(ji,jnij,jni) = jacvode(ji,jnij,jni) - 4.0 * rearatpom(ji,jk) / volc(ji,jk,jsfeo) * zlimtmpfeo 120 jnj = (jk-1)*jpvode+isvode(jwoxy) 121 jnij = jni - jnj + jpvode + 1 122 zlimtmpo2 = zlimfeo(ji,jk) * zlimdo2(ji,jk) * ( 1.0 - zlimno3(ji,jk) ) 123 jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) + 4.0 * rearatpom(ji,jk) / volc(ji,jk,jsfeo) * zlimtmpo2 124 jnj = (jk-1)*jpvode+isvode(jwno3) 125 jnij = jni - jnj + jpvode + 1 126 zlimtmpno3 = zlimfeo(ji,jk) * zlimdno3(ji,jk) * ( 1.0 - zlimo2(ji,jk) ) 127 jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) + 4.0 * rearatpom(ji,jk) / volc(ji,jk,jsfeo) * zlimtmpno3 128 ! Iron 129 zreasat = rearatpom(ji,jk) * 4.0 * radssol(jk,jwfe2) 130 jni = (jk-1)*jpvode+isvode(jwfe2) 131 jnj = (jk-1)*jpvode+isvode(jpwat+jsfeo) 132 jnij = jni - jnj + jpvode + 1 133 jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) + zreasat * zlimtmpfeo 134 jnj = (jk-1)*jpvode+isvode(jwoxy) 135 jnij = jni - jnj + jpvode + 1 136 jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) - zreasat * zlimtmpo2 137 jnj = (jk-1)*jpvode+isvode(jwno3) 138 jnij = jni - jnj + jpvode + 1 139 jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) - zreasat * zlimtmpno3 140 ENDIF 141 END DO 119 142 ENDDO 120 143 121 ! write(*,*) 'plante sedfeo1 ',jacobian(ji,15,9) / dtsed122 123 144 !-------------------------------------------------------------------- 124 145 ! Begining POC denitrification and NO3- diffusion … … 126 147 127 148 DO jk = 2, jpksed 128 zlimso4(jk) = ( 1.0 - zlimno3(jk) - zlimo2(jk) - zlimfeo(jk) ) * pwcp(ji,jk,jwso4) / ( pwcp(ji,jk,jwso4) + xksedso4 ) 129 ! For sulfur 130 jni = (jk-1) * jpvode + isvode(jwso4) 131 Jacobian(ji,jni,jni) = Jacobian(ji,jni,jni) - 0.5 * rearatpom(ji,jk) * ( 1.0 - zlimno3(jk) & 132 & - zlimo2(jk) - zlimfeo(jk) ) * xksedso4 / ( pwcp(ji,jk,jwso4) + xksedso4 )**2 133 jnj = (jk-1) * jpvode + isvode(jwoxy) 134 zlimtmpo2 = pwcp(ji,jk,jwso4) / ( pwcp(ji,jk,jwso4) + xksedso4 ) * xksedo2 & 135 & / (pwcp(ji,jk,jwoxy) + xksedo2 )**2 * ( 1.0 - pwcp(ji,jk,jwno3) / ( xksedno3 + pwcp(ji,jk,jwno3) ) ) & 136 & * ( 1.0 - solcp(ji,jk,jsfeo) / ( solcp(ji,jk,jsfeo) + xksedfeo ) ) 137 Jacobian(ji,jni,jnj) = Jacobian(ji,jni,jnj) + 0.5 * rearatpom(ji,jk) * zlimtmpo2 138 jnj = (jk-1) * jpvode + isvode(jwno3) 139 zlimtmpno3 = pwcp(ji,jk,jwso4) / ( pwcp(ji,jk,jwso4) + xksedso4 ) * xksedno3 & 140 & / ( pwcp(ji,jk,jwno3) + xksedno3 )**2 * ( 1.0 - zlimo2(jk) ) * ( 1.0 - solcp(ji,jk,jsfeo) & 141 & / ( solcp(ji,jk,jsfeo) + xksedfeo ) ) 142 Jacobian(ji,jni,jnj) = Jacobian(ji,jni,jnj) + 0.5 * rearatpom(ji,jk) * zlimtmpno3 143 jnj = (jk-1) * jpvode + isvode(jpwat+jsfeo) 144 Jacobian(ji,jni,jnj) = Jacobian(ji,jni,jnj) + 0.5 * rearatpom(ji,jk) * pwcp(ji,jk,jwso4) & 145 & / ( pwcp(ji,jk,jwso4) + xksedso4 ) * xksedfeo / ( solcp(ji,jk,jsfeo) + xksedfeo )**2 & 146 & * ( 1.0 - zlimo2(jk) - zlimno3(jk) ) 147 jni = (jk-1) * jpvode + isvode(jwh2s) 148 jnj = (jk-1) * jpvode + isvode(jwso4) 149 Jacobian(ji,jni,jnj) = Jacobian(ji,jni,jnj) + 0.5 * rearatpom(ji,jk) * ( 1.0 - zlimno3(jk) & 150 & - zlimo2(jk) - zlimfeo(jk) ) * xksedso4 / ( pwcp(ji,jk,jwso4) + xksedso4 )**2 151 jnj = (jk-1) * jpvode + isvode(jwoxy) 152 Jacobian(ji,jni,jnj) = Jacobian(ji,jni,jnj) - 0.5 * rearatpom(ji,jk) * zlimtmpo2 153 jnj = (jk-1) * jpvode + isvode(jwno3) 154 Jacobian(ji,jni,jnj) = Jacobian(ji,jni,jnj) - 0.5 * rearatpom(ji,jk) * zlimtmpno3 155 jnj = (jk-1) * jpvode + isvode(jpwat+jsfeo) 156 Jacobian(ji,jni,jnj) = Jacobian(ji,jni,jnj) - 0.5 * rearatpom(ji,jk) * pwcp(ji,jk,jwso4) & 157 & / ( pwcp(ji,jk,jwso4) + xksedso4 ) * xksedfeo / ( solcp(ji,jk,jsfeo) + xksedfeo )**2 & 158 & * ( 1.0 - zlimo2(jk) - zlimno3(jk) ) 149 DO ji = 1, jpoce 150 IF (accmask(ji) == 0) THEN 151 zlimso4(ji,jk) = pwcp(ji,jk,jwso4) / ( pwcp(ji,jk,jwso4) + xksedso4 ) 152 zlimdso4(ji,jk) = xksedso4 / ( pwcp(ji,jk,jwso4) + xksedso4 )**2 153 ! For sulfur 154 jni = (jk-1) * jpvode + isvode(jwso4) 155 jnij = jpvode + 1 156 zlimtmpso4 = ( 1.0 - zlimno3(ji,jk) ) * ( 1.0 - zlimo2(ji,jk) ) * ( 1.0 - zlimfeo(ji,jk) ) * zlimdso4(ji,jk) 157 jacvode(ji,jnij,jni) = jacvode(ji,jnij,jni) - 0.5 * rearatpom(ji,jk) * zlimtmpso4 158 jnj = (jk-1) * jpvode + isvode(jwoxy) 159 jnij = jni - jnj + jpvode + 1 160 zlimtmpo2 = zlimso4(ji,jk) * zlimdo2(ji,jk) * ( 1.0 - zlimno3(ji,jk) ) & 161 & * ( 1.0 - zlimfeo(ji,jk) ) 162 jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) + 0.5 * rearatpom(ji,jk) * zlimtmpo2 163 jnj = (jk-1) * jpvode + isvode(jwno3) 164 jnij = jni - jnj + jpvode + 1 165 zlimtmpno3 = zlimso4(ji,jk) * ( 1.0 - zlimo2(ji,jk) ) * zlimdno3(ji,jk) & 166 & * ( 1.0 - zlimfeo(ji,jk) ) 167 jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) + 0.5 * rearatpom(ji,jk) * zlimtmpno3 168 jnj = (jk-1) * jpvode + isvode(jpwat+jsfeo) 169 jnij = jni - jnj + jpvode + 1 170 zlimtmpfeo = zlimso4(ji,jk) * ( 1.0 - zlimo2(ji,jk) ) * ( 1.0 - zlimno3(ji,jk) ) * zlimdfeo(ji,jk) 171 jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) + 0.5 * rearatpom(ji,jk) * zlimtmpfeo 172 jni = (jk-1) * jpvode + isvode(jwh2s) 173 jnj = (jk-1) * jpvode + isvode(jwso4) 174 jnij = jni - jnj + jpvode + 1 175 jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) + 0.5 * rearatpom(ji,jk) * zlimtmpso4 176 jnj = (jk-1) * jpvode + isvode(jwoxy) 177 jnij = jni - jnj + jpvode + 1 178 jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) - 0.5 * rearatpom(ji,jk) * zlimtmpo2 179 jnj = (jk-1) * jpvode + isvode(jwno3) 180 jnij = jni - jnj + jpvode + 1 181 jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) - 0.5 * rearatpom(ji,jk) * zlimtmpno3 182 jnj = (jk-1) * jpvode + isvode(jpwat+jsfeo) 183 jnij = jni - jnj + jpvode + 1 184 jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) - 0.5 * rearatpom(ji,jk) * zlimtmpfeo 185 ENDIF 186 END DO 159 187 ENDDO 160 188 … … 162 190 ! ------------------------- 163 191 164 call sed_dsr_redoxbjac( ji)192 call sed_dsr_redoxbjac( NEQ, NROWPD, jacvode, accmask ) 165 193 166 194 IF( ln_timing ) CALL timing_stop('sed_dsrjac') … … 168 196 END SUBROUTINE sed_dsrjac 169 197 170 SUBROUTINE sed_dsr_redoxbjac( ji)198 SUBROUTINE sed_dsr_redoxbjac( NEQ, NROWPD, jacvode, accmask ) 171 199 !!---------------------------------------------------------------------- 172 200 !! *** ROUTINE sed_dsr_redox *** … … 178 206 !!---------------------------------------------------------------------- 179 207 !! Arguments 208 INTEGER, INTENT(in) :: NEQ, NROWPD 209 REAL, DIMENSION(jpoce,NROWPD,NEQ), INTENT(inout) :: jacvode 210 INTEGER, DIMENSION(jpoce), INTENT(in) :: accmask 180 211 ! --- local variables 181 INTEGER, INTENT(IN) :: ji 182 INTEGER :: jni, jnj, jk ! dummy looop indices 212 INTEGER :: ji, jni, jnj, jnij, jk ! dummy looop indices 183 213 184 214 REAL(wp) :: zalpha, zexcess, zh2seq, zsedfer, zdsedfer … … 189 219 190 220 DO jk = 2, jpksed 191 zalpha = ( pwcp(ji,jk,jwfe2) - pwcp(ji,jk,jwlgw) ) * 1E9 192 zsedfer = ( zalpha + SQRT( zalpha**2 + 1.E-10 ) ) /2.0 * 1E-9 193 zdsedfer = zsedfer*1E9 / SQRT( zalpha**2 +1E-10 ) 194 195 196 ! First pass of the scheme. At the end, it is 1st order 197 ! ----------------------------------------------------- 198 ! Fe (both adsorbed and solute) + O2 199 jni = (jk-1) * jpvode + isvode(jwfe2) 200 Jacobian(ji,jni,jni) = Jacobian(ji,jni,jni) - reac_fe2 * dtsed * pwcp(ji,jk,jwoxy) * zdsedfer 201 jnj = (jk-1) * jpvode + isvode(jwoxy) 202 Jacobian(ji,jni,jnj) = Jacobian(ji,jni,jnj) - reac_fe2 * dtsed * zsedfer 203 jni = (jk-1) * jpvode + isvode(jwoxy) 204 Jacobian(ji,jni,jni) = Jacobian(ji,jni,jni) - 0.25 * reac_fe2 * dtsed * zsedfer / radsfe2(jk) 205 jnj = (jk-1) * jpvode + isvode(jwfe2) 206 Jacobian(ji,jni,jnj) = Jacobian(ji,jni,jnj) - 0.25 * reac_fe2 * dtsed / radsfe2(jk) * pwcp(ji,jk,jwoxy) * zdsedfer 207 jni = (jk-1) * jpvode + isvode(jpwat+jsfeo) 208 jnj = (jk-1) * jpvode + isvode(jwfe2) 209 Jacobian(ji,jni,jnj) = Jacobian(ji,jni,jnj) + reac_fe2 * dtsed / radsfe2(jk) * pwcp(ji,jk,jwoxy) & 210 & * zdsedfer / volc(ji,jk,jsfeo) 211 jnj = (jk-1) * jpvode + isvode(jwoxy) 212 Jacobian(ji,jni,jnj) = Jacobian(ji,jni,jnj) + reac_fe2 * dtsed / radsfe2(jk) * zsedfer & 213 & / volc(ji,jk,jsfeo) 214 215 216 ! H2S + O2 217 jni = (jk-1) * jpvode + isvode(jwh2s) 218 Jacobian(ji,jni,jni) = Jacobian(ji,jni,jni) - reac_h2s * dtsed * pwcp(ji,jk,jwoxy) 219 jnj = (jk-1) * jpvode + isvode(jwoxy) 220 Jacobian(ji,jni,jnj) = Jacobian(ji,jni,jnj) - reac_h2s * dtsed * pwcp(ji,jk,jwh2s) 221 jni = (jk-1) * jpvode + isvode(jwoxy) 222 Jacobian(ji,jni,jni) = Jacobian(ji,jni,jni) - 2.0 * reac_h2s * dtsed * pwcp(ji,jk,jwh2s) 223 jnj = (jk-1) * jpvode + isvode(jwh2s) 224 Jacobian(ji,jni,jnj) = Jacobian(ji,jni,jnj) - 2.0 * reac_h2s * dtsed * pwcp(ji,jk,jwoxy) 225 jni = (jk-1) * jpvode + isvode(jwso4) 226 jnj = (jk-1) * jpvode + isvode(jwh2s) 227 Jacobian(ji,jni,jnj) = Jacobian(ji,jni,jnj) + reac_h2s * dtsed * pwcp(ji,jk,jwoxy) 228 jnj = (jk-1) * jpvode + isvode(jwoxy) 229 Jacobian(ji,jni,jnj) = Jacobian(ji,jni,jnj) + reac_h2s * dtsed * pwcp(ji,jk,jwh2s) 230 231 ! NH4 + O2 232 jni = (jk-1) * jpvode + isvode(jwnh4) 233 Jacobian(ji,jni,jni) = Jacobian(ji,jni,jni) - reac_nh4 * dtsed * pwcp(ji,jk,jwoxy) 234 jnj = (jk-1) * jpvode + isvode(jwoxy) 235 Jacobian(ji,jni,jnj) = Jacobian(ji,jni,jnj) - reac_nh4 * dtsed * pwcp(ji,jk,jwnh4) 236 jni = (jk-1) * jpvode + isvode(jwoxy) 237 Jacobian(ji,jni,jni) = Jacobian(ji,jni,jni) - 2.0 * reac_nh4 * dtsed * pwcp(ji,jk,jwnh4) / radsnh4(jk) 238 jnj = (jk-1) * jpvode + isvode(jwnh4) 239 Jacobian(ji,jni,jnj) = Jacobian(ji,jni,jnj) - 2.0 * reac_nh4 * dtsed * pwcp(ji,jk,jwoxy) / radsnh4(jk) 240 jni = (jk-1) * jpvode + isvode(jwno3) 241 jnj = (jk-1) * jpvode + isvode(jwoxy) 242 Jacobian(ji,jni,jnj) = Jacobian(ji,jni,jnj) + reac_nh4 * dtsed * pwcp(ji,jk,jwnh4) / radsnh4(jk) 243 jnj = (jk-1) * jpvode + isvode(jwnh4) 244 Jacobian(ji,jni,jnj) = Jacobian(ji,jni,jnj) + reac_nh4 * dtsed * pwcp(ji,jk,jwoxy) / radsnh4(jk) 245 246 ! FeS - O2 247 jni = (jk-1) * jpvode + isvode(jpwat+jsfes) 248 Jacobian(ji,jni,jni) = Jacobian(ji,jni,jni) - reac_feso * dtsed * pwcp(ji,jk,jwoxy) 249 jnj = (jk-1) * jpvode + isvode(jwoxy) 250 Jacobian(ji,jni,jnj) = Jacobian(ji,jni,jnj) - reac_feso * dtsed * solcp(ji,jk,jsfes) 251 jni = (jk-1) * jpvode + isvode(jwoxy) 252 Jacobian(ji,jni,jni) = Jacobian(ji,jni,jni) - 2.0 * reac_feso * dtsed * solcp(ji,jk,jsfes) & 253 & * volc(ji,jk,jsfes) 254 jnj = (jk-1) * jpvode + isvode(jpwat+jsfes) 255 Jacobian(ji,jni,jnj) = Jacobian(ji,jni,jnj) - 2.0 * reac_feso * dtsed * pwcp(ji,jk,jwoxy) & 256 & * volc(ji,jk,jsfes) 257 jni = (jk-1) * jpvode + isvode(jwfe2) 258 jnj = (jk-1) * jpvode + isvode(jwoxy) 259 Jacobian(ji,jni,jnj) = Jacobian(ji,jni,jnj) + reac_feso * dtsed * solcp(ji,jk,jsfes) & 260 & * volc(ji,jk,jsfes) * radsfe2(jk) 261 jnj = (jk-1) * jpvode + isvode(jpwat+jsfes) 262 Jacobian(ji,jni,jnj) = Jacobian(ji,jni,jnj) + reac_feso * dtsed * pwcp(ji,jk,jwoxy) & 263 & * volc(ji,jk,jsfes) * radsfe2(jk) 264 jni = (jk-1) * jpvode + isvode(jwso4) 265 jnj = (jk-1) * jpvode + isvode(jwoxy) 266 Jacobian(ji,jni,jnj) = Jacobian(ji,jni,jnj) + reac_feso * dtsed * solcp(ji,jk,jsfes) & 267 & * volc(ji,jk,jsfes) 268 jnj = (jk-1) * jpvode + isvode(jpwat+jsfes) 269 Jacobian(ji,jni,jnj) = Jacobian(ji,jni,jnj) + reac_feso * dtsed * pwcp(ji,jk,jwoxy) & 270 & * volc(ji,jk,jsfes) 271 272 ! FEOH + H2S 273 jni = (jk-1) * jpvode + isvode(jpwat+jsfeo) 274 Jacobian(ji,jni,jni) = Jacobian(ji,jni,jni) - 8.0 * reac_feh2s * dtsed * pwcp(ji,jk,jwh2s) 275 jnj = (jk-1) * jpvode + isvode(jwh2s) 276 Jacobian(ji,jni,jnj) = Jacobian(ji,jni,jnj) - 8.0 * reac_feh2s * dtsed * solcp(ji,jk,jsfeo) 277 jni = (jk-1) * jpvode + isvode(jwh2s) 278 Jacobian(ji,jni,jni) = Jacobian(ji,jni,jni) - reac_feh2s * dtsed * solcp(ji,jk,jsfeo) & 279 & * volc(ji,jk,jsfeo) 280 jnj = (jk-1) * jpvode + isvode(jpwat+jsfeo) 281 Jacobian(ji,jni,jnj) = Jacobian(ji,jni,jnj) - reac_feh2s * dtsed * pwcp(ji,jk,jwh2s) & 282 & * volc(ji,jk,jsfeo) 283 jni = (jk-1) * jpvode + isvode(jwfe2) 284 jnj = (jk-1) * jpvode + isvode(jwh2s) 285 Jacobian(ji,jni,jnj) = Jacobian(ji,jni,jnj) + 8.0 * reac_feh2s * dtsed * solcp(ji,jk,jsfeo) & 286 & * volc(ji,jk,jsfeo) * radsfe2(jk) 287 jnj = (jk-1) * jpvode + isvode(jpwat+jsfeo) 288 Jacobian(ji,jni,jnj) = Jacobian(ji,jni,jnj) + 8.0 * reac_feh2s * dtsed * pwcp(ji,jk,jwh2s) & 289 & * volc(ji,jk,jsfeo) * radsfe2(jk) 290 jni = (jk-1) * jpvode + isvode(jwso4) 291 jnj = (jk-1) * jpvode + isvode(jwh2s) 292 Jacobian(ji,jni,jnj) = Jacobian(ji,jni,jnj) + reac_feh2s * dtsed * solcp(ji,jk,jsfeo) & 293 & * volc(ji,jk,jsfeo) 294 jnj = (jk-1) * jpvode + isvode(jpwat+jsfeo) 295 Jacobian(ji,jni,jnj) = Jacobian(ji,jni,jnj) + reac_feh2s * dtsed * pwcp(ji,jk,jwh2s) & 296 & * volc(ji,jk,jsfeo) 297 298 ! Fe + H2S 299 zh2seq = MAX(rtrn, 2.1E-3 * hipor(ji,jk)) 300 zexcess = pwcp(ji,jk,jwh2s) * zsedfer / zh2seq - 1.0 301 IF ( zexcess >= 0.0 ) THEN 302 jni = (jk-1) * jpvode + isvode(jwfe2) 303 Jacobian(ji,jni,jni) = Jacobian(ji,jni,jni) - reac_fesp * dtsed * pwcp(ji,jk,jwh2s) * zdsedfer / zh2seq * radsfe2(jk) 304 jnj = (jk-1) * jpvode + isvode(jwh2s) 305 Jacobian(ji,jni,jnj) = Jacobian(ji,jni,jnj) - reac_fesp * dtsed * zsedfer / zh2seq * radsfe2(jk) 306 jni = (jk-1) * jpvode + isvode(jpwat+jsfes) 307 jnj = (jk-1) * jpvode + isvode(jwfe2) 308 Jacobian(ji,jni,jnj) = Jacobian(ji,jni,jnj) + reac_fesp * dtsed * pwcp(ji,jk,jwh2s) / zh2seq & 309 & * zdsedfer / volc(ji,jk,jsfes) 310 jnj = (jk-1) * jpvode + isvode(jwh2s) 311 Jacobian(ji,jni,jnj) = Jacobian(ji,jni,jnj) + reac_fesp * dtsed * zsedfer / zh2seq & 312 & / volc(ji,jk,jsfes) 313 jni = (jk-1) * jpvode + isvode(jwh2s) 314 Jacobian(ji,jni,jni) = Jacobian(ji,jni,jni) - reac_fesp * dtsed * zsedfer / zh2seq 315 jnj = (jk-1) * jpvode + isvode(jwfe2) 316 Jacobian(ji,jni,jnj) = Jacobian(ji,jni,jnj) - reac_fesp * dtsed * pwcp(ji,jk,jwh2s) * zdsedfer / zh2seq 317 ELSE 318 jni = (jk-1) * jpvode + isvode(jwfe2) 319 Jacobian(ji,jni,jni) = Jacobian(ji,jni,jni) - reac_fesd * dtsed * pwcp(ji,jk,jwh2s) / zh2seq & 320 & * zdsedfer * radsfe2(jk) * solcp(ji,jk,jsfes) * volc(ji,jk,jsfes) 321 jnj = (jk-1) * jpvode + isvode(jwh2s) 322 Jacobian(ji,jni,jnj) = Jacobian(ji,jni,jnj) - reac_fesd * dtsed * zsedfer / zh2seq * radsfe2(jk) & 323 & * solcp(ji,jk,jsfes) * volc(ji,jk,jsfes) 324 jnj = (jk-1) * jpvode + isvode(jpwat+jsfes) 325 Jacobian(ji,jni,jnj) = Jacobian(ji,jni,jnj) - reac_fesd * dtsed * zexcess * radsfe2(jk) * volc(ji,jk,jsfes) 326 jni = (jk-1) * jpvode + isvode(jpwat+jsfes) 327 Jacobian(ji,jni,jni) = Jacobian(ji,jni,jni) + reac_fesd * dtsed * zexcess 328 jnj = (jk-1) * jpvode + isvode(jwfe2) 329 Jacobian(ji,jni,jnj) = Jacobian(ji,jni,jnj) + reac_fesd * dtsed * solcp(ji,jk,jsfes) & 330 & * zdsedfer * pwcp(ji,jk,jwh2s) / zh2seq 331 jnj = (jk-1) * jpvode + isvode(jwh2s) 332 Jacobian(ji,jni,jnj) = Jacobian(ji,jni,jnj) + reac_fesd * dtsed * solcp(ji,jk,jsfes) & 333 & * zsedfer / zh2seq 334 jni = (jk-1) * jpvode + isvode(jwh2s) 335 Jacobian(ji,jni,jni) = Jacobian(ji,jni,jni) - reac_fesd * dtsed * zsedfer / zh2seq & 336 & * solcp(ji,jk,jsfes) * volc(ji,jk,jsfes) 337 jnj = (jk-1) * jpvode + isvode(jwfe2) 338 Jacobian(ji,jni,jnj) = Jacobian(ji,jni,jnj) - reac_fesd * dtsed * pwcp(ji,jk,jwh2s) / zh2seq & 339 & * zdsedfer * solcp(ji,jk,jsfes) * volc(ji,jk,jsfes) 340 jnj = (jk-1) * jpvode + isvode(jpwat+jsfes) 341 Jacobian(ji,jni,jnj) = Jacobian(ji,jni,jnj) - reac_fesd * dtsed * zexcess * radsfe2(jk) * volc(ji,jk,jsfes) 342 ENDIF 221 DO ji = 1, jpoce 222 IF (accmask(ji) == 0) THEN 223 zalpha = ( pwcp(ji,jk,jwfe2) - pwcp(ji,jk,jwlgw) ) * 1E9 224 zsedfer = ( zalpha + SQRT( zalpha**2 + 1.E-10 ) ) /2.0 * 1E-9 225 zdsedfer = zsedfer*1E9 / SQRT( zalpha**2 +1E-10 ) 226 227 ! First pass of the scheme. At the end, it is 1st order 228 ! ----------------------------------------------------- 229 ! Fe (both adsorbed and solute) + O2 230 jni = (jk-1) * jpvode + isvode(jwfe2) 231 jnij = jpvode + 1 232 jacvode(ji,jnij,jni) = jacvode(ji,jnij,jni) - reac_fe2 * pwcp(ji,jk,jwoxy) * zdsedfer 233 jnj = (jk-1) * jpvode + isvode(jwoxy) 234 jnij = jni - jnj + jpvode + 1 235 jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) - reac_fe2 * zsedfer 236 jni = (jk-1) * jpvode + isvode(jwoxy) 237 jnij = jpvode + 1 238 jacvode(ji,jnij,jni) = jacvode(ji,jnij,jni) - 0.25 * reac_fe2 * zsedfer / radssol(jk,jwfe2) 239 jnj = (jk-1) * jpvode + isvode(jwfe2) 240 jnij = jni - jnj + jpvode + 1 241 jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) - 0.25 * reac_fe2 / radssol(jk,jwfe2) * pwcp(ji,jk,jwoxy) * zdsedfer 242 jni = (jk-1) * jpvode + isvode(jpwat+jsfeo) 243 jnj = (jk-1) * jpvode + isvode(jwfe2) 244 jnij = jni - jnj + jpvode + 1 245 jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) + reac_fe2 / radssol(jk,jwfe2) * pwcp(ji,jk,jwoxy) & 246 & * zdsedfer / volc(ji,jk,jsfeo) 247 jnj = (jk-1) * jpvode + isvode(jwoxy) 248 jnij = jni - jnj + jpvode + 1 249 jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) + reac_fe2 / radssol(jk,jwfe2) * zsedfer & 250 & / volc(ji,jk,jsfeo) 251 252 ! H2S + O2 253 jni = (jk-1) * jpvode + isvode(jwh2s) 254 jnij = jpvode + 1 255 jacvode(ji,jnij,jni) = jacvode(ji,jnij,jni) - reac_h2s * pwcp(ji,jk,jwoxy) 256 jnj = (jk-1) * jpvode + isvode(jwoxy) 257 jnij = jni - jnj + jpvode + 1 258 jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) - reac_h2s * pwcp(ji,jk,jwh2s) 259 jni = (jk-1) * jpvode + isvode(jwoxy) 260 jnij = jpvode + 1 261 jacvode(ji,jnij,jni) = jacvode(ji,jnij,jni) - 2.0 * reac_h2s * pwcp(ji,jk,jwh2s) 262 jnj = (jk-1) * jpvode + isvode(jwh2s) 263 jnij = jni - jnj + jpvode + 1 264 jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) - 2.0 * reac_h2s * pwcp(ji,jk,jwoxy) 265 jni = (jk-1) * jpvode + isvode(jwso4) 266 jnj = (jk-1) * jpvode + isvode(jwh2s) 267 jnij = jni - jnj + jpvode + 1 268 jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) + reac_h2s * pwcp(ji,jk,jwoxy) 269 jnj = (jk-1) * jpvode + isvode(jwoxy) 270 jnij = jni - jnj + jpvode + 1 271 jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) + reac_h2s * pwcp(ji,jk,jwh2s) 272 273 ! NH4 + O2 274 jni = (jk-1) * jpvode + isvode(jwnh4) 275 jnij = jpvode + 1 276 jacvode(ji,jnij,jni) = jacvode(ji,jnij,jni) - reac_nh4 * pwcp(ji,jk,jwoxy) 277 jnj = (jk-1) * jpvode + isvode(jwoxy) 278 jnij = jni - jnj + jpvode + 1 279 jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) - reac_nh4 * pwcp(ji,jk,jwnh4) 280 jni = (jk-1) * jpvode + isvode(jwoxy) 281 jnij = jpvode + 1 282 jacvode(ji,jnij,jni) = jacvode(ji,jnij,jni) - 2.0 * reac_nh4 * pwcp(ji,jk,jwnh4) / radssol(jk,jwnh4) 283 jnj = (jk-1) * jpvode + isvode(jwnh4) 284 jnij = jni - jnj + jpvode + 1 285 jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) - 2.0 * reac_nh4 * pwcp(ji,jk,jwoxy) / radssol(jk,jwnh4) 286 jni = (jk-1) * jpvode + isvode(jwno3) 287 jnj = (jk-1) * jpvode + isvode(jwoxy) 288 jnij = jni - jnj + jpvode + 1 289 jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) + reac_nh4 * pwcp(ji,jk,jwnh4) / radssol(jk,jwnh4) 290 jnj = (jk-1) * jpvode + isvode(jwnh4) 291 jnij = jni - jnj + jpvode + 1 292 jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) + reac_nh4 * pwcp(ji,jk,jwoxy) / radssol(jk,jwnh4) 293 294 ! FeS - O2 295 jni = (jk-1) * jpvode + isvode(jpwat+jsfes) 296 jnij = jpvode + 1 297 jacvode(ji,jnij,jni) = jacvode(ji,jnij,jni) - reac_feso * pwcp(ji,jk,jwoxy) 298 jnj = (jk-1) * jpvode + isvode(jwoxy) 299 jnij = jni - jnj + jpvode + 1 300 jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) - reac_feso * solcp(ji,jk,jsfes) 301 jni = (jk-1) * jpvode + isvode(jwoxy) 302 jnij = jpvode + 1 303 jacvode(ji,jnij,jni) = jacvode(ji,jnij,jni) - 2.0 * reac_feso * solcp(ji,jk,jsfes) & 304 & * volc(ji,jk,jsfes) 305 jnj = (jk-1) * jpvode + isvode(jpwat+jsfes) 306 jnij = jni - jnj + jpvode + 1 307 jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) - 2.0 * reac_feso * pwcp(ji,jk,jwoxy) & 308 & * volc(ji,jk,jsfes) 309 jni = (jk-1) * jpvode + isvode(jwfe2) 310 jnj = (jk-1) * jpvode + isvode(jwoxy) 311 jnij = jni - jnj + jpvode + 1 312 jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) + reac_feso * solcp(ji,jk,jsfes) & 313 & * volc(ji,jk,jsfes) * radssol(jk,jwfe2) 314 jnj = (jk-1) * jpvode + isvode(jpwat+jsfes) 315 jnij = jni - jnj + jpvode + 1 316 jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) + reac_feso * pwcp(ji,jk,jwoxy) & 317 & * volc(ji,jk,jsfes) * radssol(jk,jwfe2) 318 jni = (jk-1) * jpvode + isvode(jwso4) 319 jnj = (jk-1) * jpvode + isvode(jwoxy) 320 jnij = jni - jnj + jpvode + 1 321 jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) + reac_feso * solcp(ji,jk,jsfes) & 322 & * volc(ji,jk,jsfes) 323 jnj = (jk-1) * jpvode + isvode(jpwat+jsfes) 324 jnij = jni - jnj + jpvode + 1 325 jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) + reac_feso * pwcp(ji,jk,jwoxy) & 326 & * volc(ji,jk,jsfes) 327 328 ! FEOH + H2S 329 jni = (jk-1) * jpvode + isvode(jpwat+jsfeo) 330 jnij = jpvode + 1 331 jacvode(ji,jnij,jni) = jacvode(ji,jnij,jni) - 8.0 * reac_feh2s * pwcp(ji,jk,jwh2s) 332 jnj = (jk-1) * jpvode + isvode(jwh2s) 333 jnij = jni - jnj + jpvode + 1 334 jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) - 8.0 * reac_feh2s * solcp(ji,jk,jsfeo) 335 jni = (jk-1) * jpvode + isvode(jwh2s) 336 jnij = jpvode + 1 337 jacvode(ji,jnij,jni) = jacvode(ji,jnij,jni) - reac_feh2s * solcp(ji,jk,jsfeo) & 338 & * volc(ji,jk,jsfeo) 339 jnj = (jk-1) * jpvode + isvode(jpwat+jsfeo) 340 jnij = jni - jnj + jpvode + 1 341 jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) - reac_feh2s * pwcp(ji,jk,jwh2s) & 342 & * volc(ji,jk,jsfeo) 343 jni = (jk-1) * jpvode + isvode(jwfe2) 344 jnj = (jk-1) * jpvode + isvode(jwh2s) 345 jnij = jni - jnj + jpvode + 1 346 jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) + 8.0 * reac_feh2s * solcp(ji,jk,jsfeo) & 347 & * volc(ji,jk,jsfeo) * radssol(jk,jwfe2) 348 jnj = (jk-1) * jpvode + isvode(jpwat+jsfeo) 349 jnij = jni - jnj + jpvode + 1 350 jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) + 8.0 * reac_feh2s * pwcp(ji,jk,jwh2s) & 351 & * volc(ji,jk,jsfeo) * radssol(jk,jwfe2) 352 jni = (jk-1) * jpvode + isvode(jwso4) 353 jnj = (jk-1) * jpvode + isvode(jwh2s) 354 jnij = jni - jnj + jpvode + 1 355 jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) + reac_feh2s * solcp(ji,jk,jsfeo) & 356 & * volc(ji,jk,jsfeo) 357 jnj = (jk-1) * jpvode + isvode(jpwat+jsfeo) 358 jnij = jni - jnj + jpvode + 1 359 jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) + reac_feh2s * pwcp(ji,jk,jwh2s) & 360 & * volc(ji,jk,jsfeo) 361 362 ! Fe + H2S 363 zh2seq = MAX(rtrn, 2.1E-3 * hipor(ji,jk)) 364 zexcess = pwcp(ji,jk,jwh2s) * zsedfer / zh2seq - 1.0 365 IF ( zexcess >= 0.0 ) THEN 366 jni = (jk-1) * jpvode + isvode(jwfe2) 367 jnij = jpvode + 1 368 jacvode(ji,jnij,jni) = jacvode(ji,jnij,jni) - reac_fesp * pwcp(ji,jk,jwh2s) * zdsedfer / zh2seq * radssol(jk,jwfe2) 369 jnj = (jk-1) * jpvode + isvode(jwh2s) 370 jnij = jni - jnj + jpvode + 1 371 jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) - reac_fesp * zsedfer / zh2seq * radssol(jk,jwfe2) 372 jni = (jk-1) * jpvode + isvode(jpwat+jsfes) 373 jnj = (jk-1) * jpvode + isvode(jwfe2) 374 jnij = jni - jnj + jpvode + 1 375 jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) + reac_fesp * pwcp(ji,jk,jwh2s) / zh2seq & 376 & * zdsedfer / volc(ji,jk,jsfes) 377 jnj = (jk-1) * jpvode + isvode(jwh2s) 378 jnij = jni - jnj + jpvode + 1 379 jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) + reac_fesp * zsedfer / zh2seq & 380 & / volc(ji,jk,jsfes) 381 jni = (jk-1) * jpvode + isvode(jwh2s) 382 jnij = jpvode + 1 383 jacvode(ji,jnij,jni) = jacvode(ji,jnij,jni) - reac_fesp * zsedfer / zh2seq 384 jnj = (jk-1) * jpvode + isvode(jwfe2) 385 jnij = jni - jnj + jpvode + 1 386 jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) - reac_fesp * pwcp(ji,jk,jwh2s) * zdsedfer / zh2seq 387 ELSE 388 jni = (jk-1) * jpvode + isvode(jwfe2) 389 jnij = jpvode + 1 390 jacvode(ji,jnij,jni) = jacvode(ji,jnij,jni) - reac_fesd * pwcp(ji,jk,jwh2s) / zh2seq & 391 & * zdsedfer * radssol(jk,jwfe2) * solcp(ji,jk,jsfes) * volc(ji,jk,jsfes) 392 jnj = (jk-1) * jpvode + isvode(jwh2s) 393 jnij = jni - jnj + jpvode + 1 394 jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) - reac_fesd * zsedfer / zh2seq * radssol(jk,jwfe2) & 395 & * solcp(ji,jk,jsfes) * volc(ji,jk,jsfes) 396 jnj = (jk-1) * jpvode + isvode(jpwat+jsfes) 397 jnij = jni - jnj + jpvode + 1 398 jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) - reac_fesd * zexcess * radssol(jk,jwfe2) * volc(ji,jk,jsfes) 399 jni = (jk-1) * jpvode + isvode(jpwat+jsfes) 400 jnij = jpvode + 1 401 jacvode(ji,jnij,jni) = jacvode(ji,jnij,jni) + reac_fesd * zexcess 402 jnj = (jk-1) * jpvode + isvode(jwfe2) 403 jnij = jni - jnj + jpvode + 1 404 jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) + reac_fesd * solcp(ji,jk,jsfes) & 405 & * zdsedfer * pwcp(ji,jk,jwh2s) / zh2seq 406 jnj = (jk-1) * jpvode + isvode(jwh2s) 407 jnij = jni - jnj + jpvode + 1 408 jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) + reac_fesd * solcp(ji,jk,jsfes) & 409 & * zsedfer / zh2seq 410 jni = (jk-1) * jpvode + isvode(jwh2s) 411 jnij = jpvode + 1 412 jacvode(ji,jnij,jni) = jacvode(ji,jnij,jni) - reac_fesd * zsedfer / zh2seq & 413 & * solcp(ji,jk,jsfes) * volc(ji,jk,jsfes) 414 jnj = (jk-1) * jpvode + isvode(jwfe2) 415 jnij = jni - jnj + jpvode + 1 416 jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) - reac_fesd * pwcp(ji,jk,jwh2s) / zh2seq & 417 & * zdsedfer * solcp(ji,jk,jsfes) * volc(ji,jk,jsfes) 418 jnj = (jk-1) * jpvode + isvode(jpwat+jsfes) 419 jnij = jni - jnj + jpvode + 1 420 jacvode(ji,jnij,jnj) = jacvode(ji,jnij,jnj) - reac_fesd * zexcess * volc(ji,jk,jsfes) 421 ENDIF 422 ENDIF 423 END DO 343 424 END DO 344 425 345 ! write(*,*) 'plante sedfeo2 ',jacobian(ji,15,9) / dtsed346 347 426 IF( ln_timing ) CALL timing_stop('sed_dsr_redoxbjac') 348 427 -
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/seddta.F90
r15075 r15297 8 8 USE sed 9 9 USE sedarr 10 USE sedini 10 11 USE phycst, ONLY : rday 11 12 USE iom … … 19 20 20 21 !! * Module variables 21 REAL(wp) :: rsecday ! number of second per a day22 22 REAL(wp) :: conv2 ! [kg/m2/month]-->[g/cm2/s] ( 1 month has 30 days ) 23 23 … … 78 78 IF (lwp) WRITE(numsed,*) ' sed_dta : Sediment fields' 79 79 dtsed = rDt_trc 80 rsecday = 60.* 60. * 24.81 80 conv2 = 1.0e+3 / 1.0e+4 82 81 ENDIF … … 96 95 ! ----------------------------------------------------------- 97 96 IF (ln_sediment_offline) THEN 98 DO_2D( 1, 1, 1, 1)97 DO_2D( 0, 0, 0, 0 ) 99 98 ikt = mbkt(ji,jj) 100 99 zwsbio4(ji,jj) = wsbio2 / rday … … 102 101 END_2D 103 102 ELSE 104 DO_2D( 1, 1, 1, 1)103 DO_2D( 0, 0, 0, 0 ) 105 104 ikt = mbkt(ji,jj) 106 105 zdep = e3t(ji,jj,ikt,Kmm) / rDt_trc … … 111 110 112 111 trc_data(:,:,:) = 0. 113 DO_2D( 1, 1, 1, 1)112 DO_2D( 0, 0, 0, 0 ) 114 113 ikt = mbkt(ji,jj) 115 IF ( tmask(ji,jj,ikt) == 1 ) THEN116 trc_data(ji,jj,jwsil) 117 trc_data(ji,jj,jwoxy) 118 trc_data(ji,jj,jwdic) 119 trc_data(ji,jj,jwno3) = tr(ji,jj,ikt,jpno3,Kbb) / 7.625120 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.625123 trc_data(ji,jj,jwh2s) 124 trc_data(ji,jj,jwso4) 125 trc_data(ji,jj,jwfe2) 126 trc_data(ji,jj,jwlgw) 127 trc_data(ji,jj,12 ) = MIN(tr(ji,jj,ikt,jpgsi,Kbb), 1E-4) * zwsbio4(ji,jj) * 1E3128 trc_data(ji,jj,13 ) = MIN(tr(ji,jj,ikt,jppoc,Kbb), 1E-4) * zwsbio3(ji,jj) * 1E3129 trc_data(ji,jj,14 ) = MIN(tr(ji,jj,ikt,jpgoc,Kbb), 1E-4) * zwsbio4(ji,jj) * 1E3130 trc_data(ji,jj,15) = MIN(tr(ji,jj,ikt,jpcal,Kbb), 1E-4) * zwsbio4(ji,jj) * 1E3131 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 ) )114 IF ( tmask(ji,jj,ikt) == 1.0 ) THEN 115 trc_data(ji,jj,jwsil) = tr(ji,jj,ikt,jpsil,Kbb) 116 trc_data(ji,jj,jwoxy) = tr(ji,jj,ikt,jpoxy,Kbb) 117 trc_data(ji,jj,jwdic) = tr(ji,jj,ikt,jpdic,Kbb) 118 trc_data(ji,jj,jwno3) = tr(ji,jj,ikt,jpno3,Kbb) * redNo3 / redC 119 trc_data(ji,jj,jwpo4) = tr(ji,jj,ikt,jppo4,Kbb) / redC 120 trc_data(ji,jj,jwalk) = tr(ji,jj,ikt,jptal,Kbb) 121 trc_data(ji,jj,jwnh4) = tr(ji,jj,ikt,jpnh4,Kbb) * redNo3 / redC 122 trc_data(ji,jj,jwh2s) = 0.0 123 trc_data(ji,jj,jwso4) = 0.14 * ts(ji,jj,ikt,jp_sal,Kmm) / 1.80655 / 96.062 124 trc_data(ji,jj,jwfe2) = tr(ji,jj,ikt,jpfer,Kbb) 125 trc_data(ji,jj,jwlgw) = 1E-9 126 trc_data(ji,jj,12 ) = MIN(tr(ji,jj,ikt,jpgsi,Kbb), 1E-4) * zwsbio4(ji,jj) * 1E3 127 trc_data(ji,jj,13 ) = MIN(tr(ji,jj,ikt,jppoc,Kbb), 1E-4) * zwsbio3(ji,jj) * 1E3 128 trc_data(ji,jj,14 ) = MIN(tr(ji,jj,ikt,jpgoc,Kbb), 1E-4) * zwsbio4(ji,jj) * 1E3 129 trc_data(ji,jj,15) = MIN(tr(ji,jj,ikt,jpcal,Kbb), 1E-4) * zwsbio4(ji,jj) * 1E3 130 trc_data(ji,jj,16) = ts(ji,jj,ikt,jp_tem,Kmm) 131 trc_data(ji,jj,17) = ts(ji,jj,ikt,jp_sal,Kmm) 132 trc_data(ji,jj,18 ) = ( tr(ji,jj,ikt,jpsfe,Kbb) * zwsbio3(ji,jj) + tr(ji,jj,ikt,jpbfe,Kbb) & 133 & * zwsbio4(ji,jj) ) * 1E3 / ( trc_data(ji,jj,13 ) + trc_data(ji,jj,14 ) + rtrn ) 134 trc_data(ji,jj,18 ) = MIN(1E-3, trc_data(ji,jj,18 ) ) 136 135 ENDIF 137 136 END_2D … … 148 147 CALL pack_arr ( jpoce, rainrm_dta(1:jpoce,jsopal), trc_data(1:jpi,1:jpj,12), iarroce(1:jpoce) ) 149 148 rainrm_dta(1:jpoce,jsopal) = rainrm_dta(1:jpoce,jsopal) * 1e-4 149 150 150 ! Sinking fluxes for POC in mol.m-2.s-1 ; conversion in mol.cm-2.s-1 151 151 CALL pack_arr ( jpoce, zdtap(1:jpoce), trc_data(1:jpi,1:jpj,13) , iarroce(1:jpoce) ) … … 153 153 DO ji = 1, jpoce 154 154 zzf2 = 2E-2 155 zzf1 = 0. 5155 zzf1 = 0.3 156 156 zzf0 = 1.0 - zzf1 - zzf2 157 157 zf0s = zzf0 … … 169 169 CALL pack_arr ( jpoce, rainrm_dta(1:jpoce,jscal), trc_data(1:jpi,1:jpj,15), iarroce(1:jpoce) ) 170 170 rainrm_dta(1:jpoce,jscal) = rainrm_dta(1:jpoce,jscal) * 1e-4 171 171 172 ! vector temperature [°C] and salinity 172 173 CALL pack_arr ( jpoce, temp(1:jpoce), trc_data(1:jpi,1:jpj,16), iarroce(1:jpoce) ) … … 178 179 CALL pack_arr ( jpoce, rainrm_dta(1:jpoce,jsclay), dust(1:jpi,1:jpj), iarroce(1:jpoce) ) 179 180 rainrm_dta(1:jpoce,jsclay) = rainrm_dta(1:jpoce,jsclay) * conv2 / mol_wgt(jsclay) & 180 & + wacc(1:jpoce) * por1(2) * dens sol / mol_wgt(jsclay) / ( rsecday * 365.0 )181 rainrm_dta(1:jpoce,jsfeo) = rainrm_dta(1:jpoce,jsclay) * mol_wgt(jsclay) / mol_wgt(jsfeo) * 0.035 * 0.5 * 0.333182 rainrm_dta(1:jpoce,jsclay) = rainrm_dta(1:jpoce,jsclay) * ( 1.0 - 0.035 * 0.5 * 0.333 )181 & + wacc(1:jpoce) * por1(2) * dens_sol(jsclay) / mol_wgt(jsclay) / ( rday * 365.0 ) 182 rainrm_dta(1:jpoce,jsfeo) = rainrm_dta(1:jpoce,jsclay) * mol_wgt(jsclay) / mol_wgt(jsfeo) * 0.035 * 0.5 183 rainrm_dta(1:jpoce,jsclay) = rainrm_dta(1:jpoce,jsclay) * ( 1.0 - 0.035 * 0.5 ) 183 184 CALL unpack_arr ( jpoce, zddust(1:jpi,1:jpj), iarroce(1:jpoce), wacc(1:jpoce) ) 184 zddust(:,:) = dust(:,:) + zddust(:,:) / ( r secday * 365.0 ) * por1(2) * denssol/ conv2185 zddust(:,:) = dust(:,:) + zddust(:,:) / ( rday * 365.0 ) * por1(2) * dens_sol(jsclay) / conv2 185 186 186 187 ! rainrm_dta(1:jpoce,jsclay) = 1.0E-4 * conv2 / mol_wgt(jsclay) … … 193 194 194 195 ! sediment pore water at 1st layer (k=1) 195 DO jw = 1, jpwat 196 pwcp(1:jpoce,1,jw) = pwcp_dta(1:jpoce,jw) 197 ENDDO 198 199 ! rain 200 DO js = 1, jpsol 201 rainrm(1:jpoce,js) = rainrm_dta(1:jpoce,js) 202 ENDDO 196 pwcp(1:jpoce,1,1:jpwat) = pwcp_dta(1:jpoce,1:jpwat) 203 197 204 198 ! Calculation of raintg of each sol. comp.: rainrm in [g/(cm**2.s)] 205 199 DO js = 1, jpsol 206 rainrg(1:jpoce,js) = rainrm (1:jpoce,js) *mol_wgt(js)200 rainrg(1:jpoce,js) = rainrm_dta(1:jpoce,js) * mol_wgt(js) 207 201 ENDDO 208 202 209 ! Calculation of raintg = total massic flux rained in each cell (sum of sol. comp.)210 raintg(:) = 0.203 ! computation of dzdep = total thickness of solid material rained [cm] in each cell 204 dzdep(:) = 0. 211 205 DO js = 1, jpsol 212 raintg(1:jpoce) = raintg(1:jpoce) + rainrg(1:jpoce,js) 213 ENDDO 214 215 ! computation of dzdep = total thickness of solid material rained [cm] in each cell 216 dzdep(1:jpoce) = raintg(1:jpoce) * dtsed / ( denssol * por1(2) ) 206 dzdep(1:jpoce) = dzdep(1:jpoce) + rainrg(1:jpoce,js) * dtsed / ( dens_sol(js) * por1(2) ) 207 END DO 217 208 218 209 IF( lk_iomput ) THEN -
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/sedfunc.F90
r15127 r15297 8 8 !! * Modules used 9 9 USE sed ! sediment global variable 10 USE sed_oce 10 11 USE sedini 11 12 USE seddsr … … 24 25 CONTAINS 25 26 26 SUBROUTINE sed_func( NEQ, JINDEX, T, X, fval0)27 SUBROUTINE sed_func( NEQ, X, fval0, accmask ) 27 28 !!---------------------------------------------------------------------- 28 29 !! *** ROUTINE sed_sol *** … … 49 50 !!---------------------------------------------------------------------- 50 51 !! Arguments 51 INTEGER, INTENT(in) :: NEQ , JINDEX52 REAL(wp), INTENT(in) :: T53 REAL, DIMENSION( NEQ), INTENT(in) :: X54 REAL, DIMENSION( NEQ), INTENT(out) :: fval052 INTEGER, INTENT(in) :: NEQ 53 INTEGER, DIMENSION(jpoce), INTENT(in) :: accmask 54 REAL, DIMENSION(jpoce,NEQ), INTENT(in) :: X 55 REAL, DIMENSION(jpoce,NEQ), INTENT(out) :: fval0 55 56 ! --- local variables 56 57 INTEGER :: ji, jk, js, jn ! dummy looop indices … … 60 61 IF( ln_timing ) CALL timing_start('sed_func') 61 62 ! 62 ji = JINDEX 63 pwcpa(ji,:,:) = 0. 64 solcpa(ji,:,:) = 0. 63 pwcpa(:,:,:) = 0. 64 solcpa(:,:,:) = 0. 65 65 66 66 do jn = 1, NEQ … … 68 68 js = jarr(jn,2) 69 69 IF (js <= jpwat) THEN 70 pwcp( ji,jk,js) = X(jn) * 1E-670 pwcp(:,jk,js) = X(:,jn) * 1E-6 71 71 ELSE 72 solcp( ji,jk,js-jpwat) = X(jn) * 1E-672 solcp(:,jk,js-jpwat) = X(:,jn) * 1E-6 73 73 ENDIF 74 74 END DO 75 75 76 CALL sed_dsr ( ji) ! Redox reactions76 CALL sed_dsr( accmask ) ! Redox reactions 77 77 ! Computes diffusive fluxes 78 78 DO jn = 1, jpvode 79 79 js = jsvode(jn) 80 IF (js <= jpwat) CALL sed_mat_dsr( j i, js, dtsed)80 IF (js <= jpwat) CALL sed_mat_dsr( js, accmask ) 81 81 END DO 82 call sed_mat_btb( jwnh4, accmask ) 83 call sed_mat_btb( jwfe2, accmask ) 82 84 83 85 do jn = 1, NEQ … … 85 87 js = jarr(jn,2) 86 88 IF (js <= jpwat) THEN 87 fval0( jn) = pwcpa(ji,jk,js) * 1E6 / dtsed89 fval0(:,jn) = pwcpa(:,jk,js) * 1E6 88 90 ELSE 89 fval0( jn) = solcpa(ji,jk,js-jpwat) * 1E6 / dtsed91 fval0(:,jn) = solcpa(:,jk,js-jpwat) * 1E6 90 92 ENDIF 91 93 END DO -
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/sedini.F90
r15127 r15297 10 10 !! * Modules used 11 11 USE sed ! sediment global variable 12 USE sed_oce 12 13 USE sedarr 13 14 USE sedadv … … 36 37 37 38 REAL(wp) :: & 38 rcopal = 40. , & !: reactivity for si [l.mol-1.an-1] 39 dcoef = 8.e-6 !: diffusion coefficient (*por) [cm**2/s] 39 rcopal = 40. !: reactivity for si [l.mol-1.an-1] 40 40 41 41 REAL(wp), PUBLIC :: & 42 redO2 = 13 3. , & !: Redfield coef for Oxygen42 redO2 = 138. , & !: Redfield coef for Oxygen 43 43 redNo3 = 16. , & !: Redfield coef for Nitrate 44 44 redPo4 = 1. , & !: Redfield coef for Phosphate 45 redC = 1 22. , & !: Redfield coef for Carbon45 redC = 117. , & !: Redfield coef for Carbon 46 46 redfep = 0.175 , & !: Ratio for iron bound phosphorus 47 47 rcorgl = 50. , & !: reactivity for POC/O2 [l.mol-1.an-1] … … 103 103 !! ! 06-07 (C. Ethe) Re-organization 104 104 !!---------------------------------------------------------------------- 105 INTEGER :: ji, jj, js, jn, jk, ikt, ierr 105 INTEGER :: ji, jj, js, jn, jk, ikt, ierr 106 REAL(wp) :: ztmp1, ztmp2 106 107 !!---------------------------------------------------------------------- 107 108 … … 129 130 ! Allocate SEDIMENT arrays 130 131 ierr = sed_alloc() 131 ierr = ierr + sed_ adv_alloc()132 ierr = ierr + sed_oce_alloc() 132 133 IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'sed_ini: unable to allocate sediment model arrays' ) 133 134 … … 135 136 epkbot(:,:) = 0. 136 137 gdepbot(:,:) = 0. 137 DO_2D( 1, 1, 1, 1)138 DO_2D( 0, 0, 0, 0 ) 138 139 ikt = mbkt(ji,jj) 139 140 IF( tmask(ji,jj,ikt) == 1 ) epkbot(ji,jj) = e3t_0(ji,jj,ikt) … … 155 156 ALLOCATE( pwcpa(jpoce,jpksed,jpwat) ) ; ALLOCATE( solcpa(jpoce,jpksed,jpsol) ) 156 157 ALLOCATE( solcp(jpoce,jpksed,jpsol) ) ; ALLOCATE( rainrm_dta(jpoce,jpsol) ) 157 ALLOCATE( rainr m(jpoce,jpsol) ) ; ALLOCATE( rainrg(jpoce,jpsol) ) ; ALLOCATE( raintg(jpoce) )158 ALLOCATE( rainrg(jpoce,jpsol) ) 158 159 ALLOCATE( dzdep(jpoce) ) ; ALLOCATE( iarroce(jpoce) ) ; ALLOCATE( dzkbot(jpoce) ) 160 ALLOCATE( slatit(jpoce) ) ; ALLOCATE( slongit(jpoce) ) 159 161 ALLOCATE( zkbot(jpoce) ) ; ALLOCATE( db(jpoce,jpksed) ) 160 162 ALLOCATE( temp(jpoce) ) ; ALLOCATE( salt(jpoce) ) … … 165 167 ALLOCATE( dz3d(jpoce,jpksed) ) ; ALLOCATE( volw3d(jpoce,jpksed) ) ; ALLOCATE( vols3d(jpoce,jpksed) ) 166 168 ALLOCATE( rearatpom(jpoce, jpksed) ) ; ALLOCATE( volc(jpoce,jpksed,jpsol) ) 167 ALLOCATE( Jacobian(jpoce, jpvode*jpksed, jpvode*jpksed) ) 168 ALLOCATE( radsfe2(jpksed) ) ; ALLOCATE( radsnh4(jpksed) ) 169 ALLOCATE( wacc1(jpoce) ) 169 ALLOCATE( radssol(jpksed, jpwat) ) ; ALLOCATE( rads1sol(jpksed, jpwat) ) 170 ALLOCATE( apluss(jpoce, jpksed) ) ; ALLOCATE( aminuss(jpoce,jpksed) ) 170 171 171 172 ! Initialization of global variables … … 173 174 pwcpa (:,:,:) = 0. ; solcpa(:,:,:) = 0. 174 175 solcp (:,:,:) = 0. ; rainrm_dta(:,:) = 0. 175 rainr m(:,: ) = 0. ; rainrg(:,: ) = 0. ; raintg (: ) = 0.176 rainrg(:,: ) = 0. 176 177 dzdep (: ) = 0. ; iarroce(: ) = 0 ; dzkbot (: ) = 0. 177 178 temp (: ) = 0. ; salt (: ) = 0. ; zkbot (: ) = 0. … … 180 181 dz3d (:,: ) = 0. ; volw3d (:,: ) = 0. ; vols3d (:,:) = 0. 181 182 fecratio(:) = 1E-5 ; rearatpom(:,:) = 0. 182 radsfe2(:) = 1.0 ; radsnh4(:) = 1.0 183 radssol(:,:) = 1.0 ; rads1sol(:,:) = 0. 184 apluss(:,:) = 0.0 ; aminuss(:,:) = 0.0 183 185 184 186 ! Chemical variables … … 210 212 ln_sediment_offline = .FALSE. 211 213 214 !--------------------------------------------- 215 ! Molecular weight [g/mol] for solid species 216 !--------------------------------------------- 217 218 ! opal=sio2*0.4(h20)=28+2*16+0.4*(2+16) 219 !--------------------------------------- 220 mol_wgt(jsopal) = 28. + 2. * 16. + 0.4 * ( 2. + 16. ) 221 222 ! clay 223 ! some kind of Illit (according to Pape) 224 ! K0.58(Al 1.38 Fe(III)0.37Fe(II)0.04Mg0.34)[(OH)2|(Si3.41Al0.59)O10] 225 !-------------------------------------------------------------------- 226 mol_wgt(jsclay) = 0.58 * 39. + 1.38 * 27. + ( 0.37 + 0.04 ) * 56.+ & 227 & 0.34 * 24. + 2. * ( 16. + 1. ) + 3.41 * 38. + & 228 & 0.59 * 27. + 10. * 16. 229 230 mol_wgt(jsfeo) = 55.0 + 3.0 * ( 16.0 + 1.0) 231 232 233 mol_wgt(jsfes) = 55.0 + 32.0 234 235 ! for chemistry Poc : C(122)H(244)O(86)N(16)P(1) 236 ! But den sity of Poc is an Hydrated material (= POC + 30H2O) 237 ! So C(122)H(355)O(120)N(16)P(1) 238 !------------------------------------------------------------ 239 mol_wgt(jspoc) = ( redC * 12. + 355. + 120. * 16. + redNo3 * 14. + 31. ) / redC 240 mol_wgt(jspos) = mol_wgt(jspoc) 241 mol_wgt(jspor) = mol_wgt(jspoc) 242 243 ! CaCO3 244 !--------- 245 mol_wgt(jscal) = 40. + 12. + 3. * 16. 246 247 ! Density of solid material in sediment [g/cm**3] 248 !------------------------------------------------ 249 ALLOCATE( dens_sol(jpsol) ) 250 dens_sol(jsclay) = 2.6 251 dens_sol(jscal) = 2.7 252 dens_sol(jsopal) = 2.1 253 dens_sol(jsfeo) = 3.4 254 dens_sol(jsfes) = 4.8 255 dens_sol(jspoc) = 1.0 256 dens_sol(jspos) = 1.0 257 dens_sol(jspor) = 1.0 258 259 ! Accumulation rate from Burwicz et al. (2011). This is used to 260 ! compute the flux of clays and minerals 261 ! -------------------------------------------------------------- 262 DO ji = 1, jpoce 263 ztmp1 = 0.0117 * 3.0 / ( 1.0 + ( zkbot(ji) / 200.)**3 ) 264 ztmp2 = 0.006 / ( 1.0 + ( zkbot(ji) / 5000.)**10 ) 265 wacc(ji) = ztmp2+ztmp1 266 END DO 267 212 268 ! Vertical profile of of the adsorption factor for adsorbed species 213 269 ! ----------------------------------------------------------------- 214 radsfe2(:) = 1.0 / ( 1.0 + adsfe2 * por1(:) / por(:) ) 215 radsnh4(:) = 1.0 / ( 1.0 + adsnh4 * por1(:) / por(:) ) 270 radssol(:,jwfe2) = 1.0 / ( 1.0 + adsfe2 * por1(:) / por(:) ) 271 radssol(:,jwnh4) = 1.0 / ( 1.0 + adsnh4 * por1(:) / por(:) ) 272 rads1sol(:,jwnh4) = adsnh4 * por1(:) / ( por(:) + adsnh4 * por1(:) ) 273 rads1sol(:,jwfe2) = adsfe2 * por1(:) / ( por(:) + adsfe2 * por1(:) ) 216 274 217 275 ! Initialization of the array for non linear solving 218 276 ! -------------------------------------------------- 219 277 220 ALLOCATE( jarr(jp vode*jpksed,2) )278 ALLOCATE( jarr(jpksed*jpvode,2) ) 221 279 ALLOCATE( jsvode(jpvode), isvode(jptrased) ) 222 280 jsvode(1) = jwoxy ; jsvode(2) = jwno3 ; jsvode(3) = jwnh4 … … 234 292 END DO 235 293 236 ALLOCATE( stepros(jpoce) )294 ALLOCATE( rstepros(jpoce) ) 237 295 238 296 #if defined key_sed_off … … 277 335 ! Computation of 1D array of sediments points 278 336 indoce = 0 279 DO_2D( 1, 1, 1, 1)337 DO_2D( 0, 0, 0, 0 ) 280 338 IF ( epkbot(ji,jj) > 0. ) THEN 281 339 indoce = indoce + 1 … … 300 358 !------------------------------------------------ 301 359 CALL pack_arr ( jpoce, dzkbot(1:jpoce), epkbot(1:jpi,1:jpj), iarroce(1:jpoce) ) 360 #if defined key_sed_off 361 dzkbot(1:jpoce) = 1.E8 362 #else 302 363 dzkbot(1:jpoce) = dzkbot(1:jpoce) * 1.e+2 364 #endif 303 365 CALL pack_arr ( jpoce, zkbot(1:jpoce), gdepbot(1:jpi,1:jpj), iarroce(1:jpoce) ) 366 CALL pack_arr ( jpoce, slatit(1:jpoce), gphit(1:jpi,1:jpj), iarroce(1:jpoce) ) 367 CALL pack_arr ( jpoce, slongit(1:jpoce), glamt(1:jpi,1:jpj), iarroce(1:jpoce) ) 304 368 305 369 ! Geometry and constants … … 332 396 por(1) = 1.0 333 397 DO jk = 2, jpksed 334 por(jk) = porinf + ( porsurf-porinf) * exp(-rhox * (profsed(jk)) )398 por(jk) = porinf + ( porsurf-porinf) * exp(-rhox * profsed(jk) ) 335 399 END DO 336 400 … … 356 420 dz3d(1:jpoce,1) = dz(1) 357 421 358 !---------------------------------------------359 ! Molecular weight [g/mol] for solid species360 !---------------------------------------------361 362 ! opal=sio2*0.4(h20)=28+2*16+0.4*(2+16)363 !---------------------------------------364 mol_wgt(jsopal) = 28. + 2. * 16. + 0.4 * ( 2. + 16. )365 366 ! clay367 ! some kind of Illit (according to Pape)368 ! K0.58(Al 1.38 Fe(III)0.37Fe(II)0.04Mg0.34)[(OH)2|(Si3.41Al0.59)O10]369 !--------------------------------------------------------------------370 mol_wgt(jsclay) = 0.58 * 39. + 1.38 * 27. + ( 0.37 + 0.04 ) * 56.+ &371 & 0.34 * 24. + 2. * ( 16. + 1. ) + 3.41 * 38. + &372 & 0.59 * 27. + 10. * 16.373 374 mol_wgt(jsfeo) = 55.0 + 3.0 * ( 16.0 + 1.0)375 376 mol_wgt(jsfes) = 55.0 + 32.0377 378 ! for chemistry Poc : C(122)H(244)O(86)N(16)P(1)379 ! But den sity of Poc is an Hydrated material (= POC + 30H2O)380 ! So C(122)H(355)O(120)N(16)P(1)381 !------------------------------------------------------------382 mol_wgt(jspoc) = ( 122. * 12. + 355. + 120. * 16.+ &383 & 16. * 14. + 31. ) / 122.384 mol_wgt(jspos) = mol_wgt(jspoc)385 mol_wgt(jspor) = mol_wgt(jspoc)386 387 ! CaCO3388 !---------389 mol_wgt(jscal) = 40. + 12. + 3. * 16.390 391 ! Density of solid material in sediment [g/cm**3]392 !------------------------------------------------393 denssol = 2.6394 395 ! Accumulation rate from Burwicz et al. (2011). This is used to396 ! compute the flux of clays and minerals397 ! --------------------------------------------------------------398 DO ji = 1, jpoce399 ! ztmp1 = 0.117 / ( 1.0 + ( zkbot(ji) / 200.)**3 )400 ztmp1 = 0.401 ztmp2 = 0.006 / ( 1.0 + ( zkbot(ji) / 4000.)**10 )402 wacc(ji) = ztmp2+ztmp1403 END DO404 405 406 ! Initialization of time step as function of porosity [cm**2/s]407 !------------------------------------------------------------------408 422 END SUBROUTINE sed_ini_geom 409 423 … … 435 449 TYPE(PSED) , DIMENSION(jpdia2dsed) :: seddiag2d 436 450 437 NAMELIST/nam_run/ln_sed_2way 451 NAMELIST/nam_run/ln_sed_2way,nrosorder,rosatol,rosrtol 438 452 NAMELIST/nam_geom/jpksed, sedzmin, sedhmax, sedkth, sedacr, porsurf, porinf, rhox 439 453 NAMELIST/nam_trased/sedsol, sedwat 440 454 NAMELIST/nam_diased/seddiag3d, seddiag2d 441 NAMELIST/nam_inorg/rcopal, dcoef,rccal, ratligc, rcligc455 NAMELIST/nam_inorg/rcopal, rccal, ratligc, rcligc 442 456 NAMELIST/nam_poc/redO2, redNo3, redPo4, redC, redfep, rcorgl, rcorgs, & 443 457 & rcorgr, rcnh4, rch2s, rcfe2, rcfeh2s, rcfeso, rcfesp, & … … 482 496 WRITE(numsed,*) ' namelist nam_run' 483 497 WRITE(numsed,*) ' 2-way coupling between PISCES and Sed ln_sed_2way = ', ln_sed_2way 498 WRITE(numsed,*) ' Order of the Rosenbrock method (2,3,4) = ', nrosorder 499 WRITE(numsed,*) ' Tolerance for absolute error = ', rosatol 500 WRITE(numsed,*) ' Tolerance for relative order = ', rosrtol 484 501 ENDIF 485 502 … … 595 612 WRITE(numsed,*) ' namelist nam_inorg' 596 613 WRITE(numsed,*) ' reactivity for Si rcopal = ', rcopal 597 WRITE(numsed,*) ' diff. coef for por. dcoef = ', dcoef598 614 WRITE(numsed,*) ' reactivity for calcite rccal = ', rccal 599 615 WRITE(numsed,*) ' L/C ratio in POC ratligc = ', ratligc … … 659 675 reac_fesd = rcfesd / ryear 660 676 661 662 677 ! reactivity rc in [l.mol-1.s-1] 663 678 reac_cal = rccal / ryear -
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/sedinitrc.F90
r15127 r15297 16 16 USE sedchem 17 17 USE sedarr 18 USE sed_oce 18 19 USE lib_mpp ! distribued memory computing library 19 20 -
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/sedinorg.F90
r15127 r15297 7 7 !! * Modules used 8 8 USE sed ! sediment global variable 9 USE sed_oce 9 10 USE sedini 11 USE sedmat 10 12 USE lib_mpp ! distribued memory computing library 11 13 USE lib_fortran … … 45 47 ! --- local variables 46 48 INTEGER :: ji,jk ! dummy looop indices 47 REAL(wp) :: zsieq48 REAL(wp) :: zsolid1, zreasat 49 REAL(wp), DIMENSION(jpoce) :: zsieq, reac_silf 50 REAL(wp) :: zsolid1, zreasat, zco3sat 49 51 REAL(wp) :: zsatur, zsatur2, znusil, zsolcpcl, zsolcpsi, zexcess 52 REAL(wp), DIMENSION(jpoce,jpksed) :: zundsat, zrearat, psms 50 53 !! 51 54 !!---------------------------------------------------------------------- … … 71 74 END DO 72 75 zsolcpsi = MAX( zsolcpsi, rtrn ) 73 zsieq = sieqs(ji) * MAX(0.25, 1.0 - (0.045 * zsolcpcl / zsolcpsi )**0.58 ) 76 zsieq(ji) = sieqs(ji) * MAX(0.2, 1.0 - (0.045 * zsolcpcl / zsolcpsi )**0.58 ) 77 reac_silf(ji) = reac_sil * ( 0.05 + 0.055 * ( 1.64 * ( zsolcpcl / zsolcpsi + 0.01 ) )**(-0.75) ) / 1.25 78 END DO 74 79 75 !----------------------------------------------------------76 ! 5. Beginning of Pore Water diffusion and solid reaction77 !---------------------------------------------------------78 80 79 !----------------------------------------------------------------------------- 80 ! For jk=2,jpksed, and for couple 81 ! 1 : jwsil/jsopal ( SI/Opal ) 82 ! 2 : jsclay/jsclay ( clay/clay ) 83 ! 3 : jwoxy/jspoc ( O2/POC ) 84 ! reaction rate is a function of solid=concentration in solid reactif in [mol/l] 85 ! and undersaturation in [mol/l]. 86 ! Solid weight fractions should be in ie [mol/l]) 87 ! second member and solution are in zundsat variable 88 !------------------------------------------------------------------------- 81 DO ji = 1, jpoce 89 82 DO jk = 2, jpksed 90 83 zsolid1 = volc(ji,jk,jsopal) * solcp(ji,jk,jsopal) 91 zsatur = MAX(0., ( zsieq - pwcp(ji,jk,jwsil) ) / zsieq)84 zsatur = MAX(0., ( zsieq(ji) - pwcp(ji,jk,jwsil) ) / zsieq(ji) ) 92 85 zsatur2 = (1.0 + temp(ji) / 400.0 )**37 93 86 znusil = ( 0.225 * ( 1.0 + temp(ji) / 15.) * zsatur + 0.775 * zsatur2 * zsatur**9.25 ) 94 solcp(ji,jk,jsopal) = solcp(ji,jk,jsopal) - reac_sil * znusil * dtsed * solcp(ji,jk,jsopal)95 pwcp(ji,jk,jwsil) = pwcp(ji,jk,jwsil) + reac_sil * znusil * dtsed * zsolid187 solcp(ji,jk,jsopal) = solcp(ji,jk,jsopal) - reac_silf(ji) * znusil * dtsed * solcp(ji,jk,jsopal) 88 pwcp(ji,jk,jwsil) = pwcp(ji,jk,jwsil) + reac_silf(ji) * znusil * dtsed * zsolid1 96 89 END DO 97 90 END DO … … 104 97 ! *densSW(l)**2 converts aksps [mol2/kg sol2] into [mol2/l2] to get [undsat] in [mol/l] 105 98 DO ji = 1, jpoce 106 saturco3(ji,:) = 1.0 - co3por(ji,:) * calcon2(ji) / ( aksps(ji) * densSW(ji) * densSW(ji) + rtrn ) 107 DO jk = 2, jpksed 99 zco3sat = aksps(ji) * densSW(ji) * densSW(ji) / ( calcon2(ji) + rtrn ) 100 saturco3(ji,:) = 1.0 - co3por(ji,:) / ( rtrn + zco3sat ) 101 DO jk = 1, jpksed 108 102 zsolid1 = volc(ji,jk,jscal) * solcp(ji,jk,jscal) 109 zexcess = MAX( 0., saturco3(ji,jk) ) 110 zreasat = reac_cal * dtsed * zexcess * zsolid1 111 solcp(ji,jk,jscal) = solcp(ji,jk,jscal) - zreasat / volc(ji,jk,jscal) 103 zundsat(ji,jk) = MAX( 0., zco3sat - co3por(ji,jk) ) 104 zrearat(ji,jk) = ( reac_cal * zsolid1 / ( zco3sat + rtrn ) ) / & 105 & ( 1. + reac_cal * dtsed * zundsat(ji,jk) / ( zco3sat + rtrn ) ) 106 END DO 107 END DO 108 109 psms(:,:) = 0.0 110 ! solves tridiagonal system 111 CALL sed_mat_dsri( jwdic, -zrearat(:,:), psms(:,:), dtsed, zundsat ) 112 113 ! New solid concentration values (jk=2 to jksed) for cacO3 114 DO jk = 2, jpksed 115 DO ji = 1, jpoce 116 zreasat = zrearat(ji,jk) * dtsed * zundsat(ji,jk) / ( volc(ji,jk,jscal) + rtrn ) 117 solcp(ji,jk,jscal) = solcp(ji,jk,jscal) - zreasat 118 zreasat = zrearat(ji,jk) * dtsed * zundsat(ji,jk) 112 119 ! For DIC 113 120 pwcp(ji,jk,jwdic) = pwcp(ji,jk,jwdic) + zreasat 114 121 ! For alkalinity 115 pwcp(ji,jk,jwalk) = pwcp(ji,jk,jwalk) + 2.0 * zreasat 116 END DO 117 END DO 122 pwcp(ji,jk,jwalk) = pwcp(ji,jk,jwalk) + 2.* zreasat 123 ENDDO 124 ENDDO 125 118 126 119 127 IF( ln_timing ) CALL timing_stop('sed_inorg') -
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/sedjac.F90
r15127 r15297 8 8 !! * Modules used 9 9 USE sed ! sediment global variable 10 USE sed_oce 10 11 USE sedini 11 12 USE seddsrjac … … 24 25 CONTAINS 25 26 26 SUBROUTINE sed_jac( NEQ, JINDEX, T, X, jacvode, NROWPD)27 SUBROUTINE sed_jac( NEQ, X, jacvode, NROWPD, accmask ) 27 28 !!---------------------------------------------------------------------- 28 29 !! *** ROUTINE sed_sol *** … … 49 50 !!---------------------------------------------------------------------- 50 51 !! Arguments 51 INTEGER, INTENT(in) :: NEQ, NROWPD , JINDEX52 REAL(wp), INTENT(in) :: T53 REAL, DIMENSION( NEQ), INTENT(in) :: X54 REAL, DIMENSION( NROWPD,NEQ), INTENT(out) :: jacvode52 INTEGER, INTENT(in) :: NEQ, NROWPD 53 INTEGER, DIMENSION(jpoce), INTENT(in) :: accmask 54 REAL, DIMENSION(jpoce,NEQ), INTENT(in) :: X 55 REAL, DIMENSION(jpoce,NROWPD,NEQ), INTENT(out) :: jacvode 55 56 ! --- local variables 56 INTEGER :: j i, jk, js, jn, jni, jnj! dummy looop indices57 INTEGER :: jk, js, jn ! dummy looop indices 57 58 !! 58 59 !!---------------------------------------------------------------------- … … 60 61 IF( ln_timing ) CALL timing_start('sed_jac') 61 62 ! 62 ji = jindex63 Jacobian(ji,:,:) = 0.64 63 jacvode = 0. 65 64 … … 68 67 js = jarr(jn,2) 69 68 IF (js <= jpwat) THEN 70 pwcp( ji,jk,js) = X(jn) * 1E-669 pwcp(:,jk,js) = X(:,jn) * 1E-6 71 70 ELSE 72 solcp( ji,jk,js-jpwat) = X(jn) * 1E-671 solcp(:,jk,js-jpwat) = X(:,jn) * 1E-6 73 72 ENDIF 74 73 END DO 75 74 76 CALL sed_dsrjac ( ji ) ! Redox reactions 75 CALL sed_dsrjac( NEQ, NROWPD, jacvode, accmask ) ! Redox reactions 76 77 77 ! Computes diffusive fluxes 78 78 DO jn = 1, jpvode 79 79 js = jsvode(jn) 80 IF (js <= jpwat) CALL sed_mat_dsrjac( j i, js, dtsed)80 IF (js <= jpwat) CALL sed_mat_dsrjac( js, NEQ, NROWPD, jacvode, accmask ) 81 81 END DO 82 83 do jni = 1, NEQ 84 DO jnj = 1, NEQ 85 IF (Jacobian(ji,jni,jnj) .NE. 0.) THEN 86 jacvode(jni-jnj + jpvode +1 , jnj) = Jacobian(ji,jni,jnj) / dtsed 87 ENDIF 88 END DO 89 END DO 82 CALL sed_mat_btbjac( jwnh4, NEQ, NROWPD, jacvode, accmask ) 83 CALL sed_mat_btbjac( jwfe2, NEQ, NROWPD, jacvode, accmask ) 90 84 91 85 IF( ln_timing ) CALL timing_stop('sed_jac') -
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/sedmat.F90
r15115 r15297 18 18 PUBLIC sed_mat_dsri 19 19 PUBLIC sed_mat_btb 20 21 INTEGER, PARAMETER :: nmax = 6022 20 PUBLIC sed_mat_btbjac 21 PUBLIC sed_mat_btbi 22 PUBLIC sed_mat_coef 23 23 24 24 !! $Id$ 25 25 CONTAINS 26 26 27 SUBROUTINE sed_mat_ dsr( ji, nvar, dtsed_in )27 SUBROUTINE sed_mat_coef 28 28 !!--------------------------------------------------------------------- 29 !! *** ROUTINE sed_mat_ dsr***29 !! *** ROUTINE sed_mat_coef *** 30 30 !! 31 31 !! ** Purpose : solves tridiagonal system of linear equations … … 48 48 !!---------------------------------------------------------------------- 49 49 !! * Arguments 50 INTEGER , INTENT(in) :: ji, nvar ! number of variable 51 52 REAL(wp), INTENT(in) :: dtsed_in 53 50 INTEGER :: ji 51 52 !---Local declarations 53 INTEGER :: jk 54 REAL(wp) :: aplus, aminus, dxplus, dxminus 55 !---------------------------------------------------------------------- 56 57 IF( ln_timing ) CALL timing_start('sed_mat_coef') 58 59 ! Computation left hand side of linear system of 60 ! equations for dissolution reaction 61 !--------------------------------------------- 62 ! first sediment level 63 DO ji = 1, jpoce 64 aplus = ( por(1) + por(2) ) * 0.5 65 dxplus = ( dz3d(ji,1) + dz3d(ji,2) ) / 2. 66 apluss(ji,1) = ( 1.0 / ( volw3d(ji,1) ) ) * aplus / dxplus 67 68 DO jk = 2, jpksed - 1 69 aminus = ( por(jk-1) + por(jk) ) * 0.5 70 dxminus = ( dz3d(ji,jk-1) + dz3d(ji,jk) ) / 2. 71 72 aplus = ( por(jk+1) + por(jk) ) * 0.5 73 dxplus = ( dz3d(ji,jk) + dz3d(ji,jk+1) ) / 2 74 ! 75 aminuss(ji,jk) = ( 1.0 / volw3d(ji,jk) ) * aminus / dxminus 76 apluss (ji,jk) = ( 1.0 / volw3d(ji,jk) ) * aplus / dxplus 77 END DO 78 79 aminus = ( por(jpksed-1) + por(jpksed) ) * 0.5 80 dxminus = ( dz3d(ji,jpksed-1) + dz3d(ji,jpksed) ) / 2. 81 aminuss(ji,jpksed) = ( 1.0 / volw3d(ji,jpksed) ) * aminus / dxminus 82 83 END DO 84 85 IF( ln_timing ) CALL timing_stop('sed_mat_coef') 86 87 END SUBROUTINE sed_mat_coef 88 89 SUBROUTINE sed_mat_dsr( nvar, accmask ) 90 !!--------------------------------------------------------------------- 91 !! *** ROUTINE sed_mat_dsr *** 92 !! 93 !! ** Purpose : solves tridiagonal system of linear equations 94 !! 95 !! ** Method : 96 !! 1 - computes left hand side of linear system of equations 97 !! for dissolution reaction 98 !! For mass balance in kbot+sediment : 99 !! dz3d (:,1) = dz(1) = 0.5 cm 100 !! volw3d(:,1) = dzkbot ( see sedini.F90 ) 101 !! dz(2) = 0.3 cm 102 !! dz3d(:,2) = 0.3 + dzdep ( see seddsr.F90 ) 103 !! volw3d(:,2) and vols3d(l,2) are thickened ( see seddsr.F90 ) 104 !! 105 !! 2 - forward/backward substitution. 106 !! 107 !! History : 108 !! ! 04-10 (N. Emprin, M. Gehlen ) original 109 !! ! 06-04 (C. Ethe) Module Re-organization 110 !!---------------------------------------------------------------------- 111 !! * Arguments 112 INTEGER , INTENT(in) :: nvar ! number of variable 113 INTEGER, DIMENSION(jpoce) :: accmask 114 INTEGER :: ji 115 54 116 !---Local declarations 55 117 INTEGER :: jk, jn 56 118 REAL(wp), DIMENSION(jpksed) :: za, zb, zc 57 119 58 REAL(wp) :: aplus,aminus59 120 REAL(wp) :: rplus,rminus 60 REAL(wp) :: dxplus,dxminus61 62 121 !---------------------------------------------------------------------- 63 122 … … 67 126 ! equations for dissolution reaction 68 127 !--------------------------------------------- 69 70 71 128 jn = nvar 72 129 ! first sediment level 73 aplus = ( por(1) + por(2) ) * 0.574 dxplus = ( dz3d(ji,1) + dz3d(ji,2) ) / 2.75 rplus = ( dtsed_in / ( volw3d(ji,1) ) ) * diff(ji,1,jn) * aplus / dxplus76 77 za(1) = 0.78 zb(1) = rplus79 zc(1) = -rplus130 DO ji = 1, jpoce 131 IF (accmask(ji) == 0) THEN 132 rplus = apluss(ji,1) * diff(ji,1,jn) * radssol(1,jn) 133 134 za(1) = 0. 135 zb(1) = rplus 136 zc(1) = -rplus 80 137 81 DO jk = 2, jpksed - 1 82 aminus = ( por(jk-1) + por(jk) ) * 0.5 83 dxminus = ( dz3d(ji,jk-1) + dz3d(ji,jk) ) / 2. 84 85 aplus = ( por(jk+1) + por(jk) ) * 0.5 86 dxplus = ( dz3d(ji,jk) + dz3d(ji,jk+1) ) / 2 87 ! 88 rminus = ( dtsed_in / volw3d(ji,jk) ) * diff(ji,jk-1,jn) * aminus / dxminus 89 rplus = ( dtsed_in / volw3d(ji,jk) ) * diff(ji,jk,jn) * aplus / dxplus 90 ! 91 za(jk) = -rminus 92 zb(jk) = rminus + rplus 93 zc(jk) = -rplus 138 DO jk = 2, jpksed - 1 139 rminus = aminuss(ji,jk) * diff(ji,jk-1,jn) * radssol(jk,jn) 140 rplus = apluss (ji,jk) * diff(ji,jk,jn) * radssol(jk,jn) 141 ! 142 za(jk) = -rminus 143 zb(jk) = rminus + rplus 144 zc(jk) = -rplus 145 END DO 146 147 rminus = aminuss(ji,jpksed) * diff(ji,jpksed-1,jn) * radssol(jpksed,jn) 148 ! 149 za(jpksed) = -rminus 150 zb(jpksed) = rminus 151 zc(jpksed) = 0. 152 153 ! solves tridiagonal system of linear equations 154 ! ----------------------------------------------- 155 156 pwcpa(ji,1,jn) = pwcpa(ji,1,jn) - ( zc(1) * pwcp(ji,2,jn) + zb(1) * pwcp(ji,1,jn) ) 157 DO jk = 2, jpksed - 1 158 pwcpa(ji,jk,jn) = pwcpa(ji,jk,jn) - ( zc(jk) * pwcp(ji,jk+1,jn) + za(jk) * pwcp(ji,jk-1,jn) & 159 & + zb(jk) * pwcp(ji,jk,jn) ) 160 ENDDO 161 pwcpa(ji,jpksed,jn) = pwcpa(ji,jpksed,jn) - ( za(jpksed) * pwcp(ji,jpksed-1,jn) & 162 & + zb(jpksed) * pwcp(ji,jpksed,jn) ) 163 164 ENDIF 94 165 END DO 95 166 96 aminus = ( por(jpksed-1) + por(jpksed) ) * 0.597 dxminus = ( dz3d(ji,jpksed-1) + dz3d(ji,jpksed) ) / 2.98 rminus = ( dtsed_in / volw3d(ji,jpksed) ) * diff(ji,jpksed-1,jn) * aminus / dxminus99 !100 za(jpksed) = -rminus101 zb(jpksed) = rminus102 zc(jpksed) = 0.103 104 ! solves tridiagonal system of linear equations105 ! -----------------------------------------------106 107 pwcpa(ji,1,jn) = pwcpa(ji,1,jn) - ( zc(1) * pwcp(ji,2,jn) + zb(1) * pwcp(ji,1,jn) )108 DO jk = 2, jpksed - 1109 pwcpa(ji,jk,jn) = pwcpa(ji,jk,jn) - ( zc(jk) * pwcp(ji,jk+1,jn) + za(jk) * pwcp(ji,jk-1,jn) &110 & + zb(jk) * pwcp(ji,jk,jn) )111 ENDDO112 pwcpa(ji,jpksed,jn) = pwcpa(ji,jpksed,jn) - ( za(jk) * pwcp(ji,jk-1,jn) + zb(jk) * pwcp(ji,jk,jn) )113 114 167 IF( ln_timing ) CALL timing_stop('sed_mat_dsr') 115 168 116 169 END SUBROUTINE sed_mat_dsr 117 170 118 SUBROUTINE sed_mat_dsrjac( ji, nvar, dtsed_in)171 SUBROUTINE sed_mat_dsrjac( nvar, NEQ, NROWPD, jacvode, accmask ) 119 172 !!--------------------------------------------------------------------- 120 173 !! *** ROUTINE sed_mat_dsrjac *** … … 140 193 !!---------------------------------------------------------------------- 141 194 !! * Arguments 142 INTEGER , INTENT(in) :: ji, nvar! number of variable143 144 REAL(wp), INTENT(in) :: dtsed_in195 INTEGER , INTENT(in) :: nvar, NEQ, NROWPD ! number of variable 196 REAL, DIMENSION(jpoce,NROWPD,NEQ), INTENT(inout) :: jacvode 197 INTEGER, DIMENSION(jpoce), INTENT(in) :: accmask 145 198 146 199 !---Local declarations 147 INTEGER :: j k, jn, jnn200 INTEGER :: ji,jk, jn, jnn, jni, jnj ,jnij 148 201 REAL(wp), DIMENSION(jpksed) :: za, zb, zc 149 202 150 REAL(wp) :: aplus,aminus151 203 REAL(wp) :: rplus,rminus 152 REAL(wp) :: dxplus,dxminus153 154 204 !---------------------------------------------------------------------- 155 205 … … 159 209 ! equations for dissolution reaction 160 210 !--------------------------------------------- 161 162 163 211 jn = nvar 164 212 ! first sediment level 165 aplus = ( por(1) + por(2) ) *0.5 166 dxplus = ( dz3d(ji,1) + dz3d(ji,2) ) / 2. 167 rplus = ( dtsed_in / ( volw3d(ji,1) ) ) * diff(ji,1,jn) * aplus / dxplus 168 169 za(1) = 0. 170 zb(1) = rplus 171 zc(1) = -rplus 172 173 DO jk = 2, jpksed - 1 174 aminus = ( por(jk-1) + por(jk) ) * 0.5 175 dxminus = ( dz3d(ji,jk-1) + dz3d(ji,jk) ) / 2. 176 177 aplus = ( por(jk+1) + por(jk) ) * 0.5 178 dxplus = ( dz3d(ji,jk) + dz3d(ji,jk+1) ) / 2 179 ! 180 rminus = ( dtsed_in / volw3d(ji,jk) ) * diff(ji,jk-1,jn) * aminus / dxminus 181 rplus = ( dtsed_in / volw3d(ji,jk) ) * diff(ji,jk,jn) * aplus / dxplus 182 ! 183 za(jk) = -rminus 184 zb(jk) = rminus + rplus 185 zc(jk) = -rplus 213 DO ji = 1, jpoce 214 IF (accmask(ji) == 0 ) THEN 215 rplus = apluss(ji,1) * diff(ji,1,jn) * radssol(1,jn) 216 217 za(1) = 0. 218 zb(1) = rplus 219 zc(1) = -rplus 220 221 DO jk = 2, jpksed - 1 222 rminus = aminuss(ji,jk) * diff(ji,jk-1,jn) * radssol(jk,jn) 223 rplus = apluss (ji,jk) * diff(ji,jk,jn) * radssol(jk,jn) 224 ! 225 za(jk) = -rminus 226 zb(jk) = rminus + rplus 227 zc(jk) = -rplus 228 END DO 229 230 rminus = aminuss(ji,jpksed) * diff(ji,jpksed-1,jn) * radssol(jpksed,jn) 231 ! 232 za(jpksed) = -rminus 233 zb(jpksed) = rminus 234 zc(jpksed) = 0. 235 236 ! solves tridiagonal system of linear equations 237 238 jnn = isvode(jn) 239 jnij = jpvode + 1 240 jacvode(ji, jnij, jnn) = jacvode(ji,jnij,jnn) - zb(1) 241 jnj = jpvode + jnn 242 jnij = jnn - jnj + jpvode + 1 243 jacvode(ji, jnij, jnj) = jacvode(ji, jnij, jnj) -zc(1) 244 DO jk = 2, jpksed - 1 245 jni = (jk-1) * jpvode + jnn 246 jnj = (jk-2) * jpvode + jnn 247 jnij = jni - jnj + jpvode + 1 248 jacvode(ji, jnij, jnj) = jacvode(ji, jnij, jnj) - za(jk) 249 jnj = (jk-1) * jpvode + jnn 250 jnij = jni - jnj + jpvode + 1 251 jacvode(ji, jnij, jnj) = jacvode(ji, jnij, jnj) - zb(jk) 252 jnj = (jk) * jpvode + jnn 253 jnij = jni - jnj + jpvode + 1 254 jacvode(ji, jnij, jnj) = jacvode(ji, jnij, jnj) - zc(jk) 255 END DO 256 jni = (jpksed-1) * jpvode + jnn 257 jnj = (jpksed-2) * jpvode + jnn 258 jnij = jni - jnj + jpvode + 1 259 jacvode(ji, jnij, jnj) = jacvode(ji, jnij, jnj) - za(jpksed) 260 jnij = jpvode + 1 261 jacvode(ji, jnij, jni) = jacvode(ji, jnij, jni) - zb(jpksed) 262 ENDIF 186 263 END DO 187 264 188 aminus = ( por(jpksed-1) + por(jpksed) ) * 0.5189 dxminus = ( dz3d(ji,jpksed-1) + dz3d(ji,jpksed) ) / 2.190 rminus = ( dtsed_in / volw3d(ji,jpksed) ) * diff(ji,jpksed-1,jn) * aminus / dxminus191 !192 za(jpksed) = -rminus193 zb(jpksed) = rminus194 zc(jpksed) = 0.195 196 ! solves tridiagonal system of linear equations197 198 IF (isvode(jn) > 0) THEN199 jnn = isvode(jn)200 Jacobian(ji, jnn, jnn) = Jacobian(ji,jnn,jnn) - zb(1)201 Jacobian(ji, jnn, jpvode + jnn) = Jacobian(ji, jnn, jpvode + jnn) -zc(1)202 DO jk = 2, jpksed - 1203 Jacobian(ji, (jk-1) * jpvode + jnn, (jk-2) * jpvode + jnn) = Jacobian(ji, (jk-1) * jpvode + jnn, (jk-2) * jpvode + jnn) - za(jk)204 Jacobian(ji, (jk-1) * jpvode + jnn, (jk-1) * jpvode + jnn) = Jacobian(ji, (jk-1) * jpvode + jnn, (jk-1) * jpvode + jnn) - zb(jk)205 Jacobian(ji, (jk-1) * jpvode + jnn, (jk) * jpvode + jnn) = Jacobian(ji, (jk-1) * jpvode + jnn, (jk) * jpvode + jnn) - zc(jk)206 END DO207 Jacobian(ji, (jpksed-1) * jpvode + jnn, (jpksed-2) * jpvode + jnn) = Jacobian(ji, (jpksed-1) * jpvode + jnn, (jpksed-2) * jpvode + jnn) - za(jpksed)208 Jacobian(ji, (jpksed-1) * jpvode + jnn, (jpksed-1) * jpvode + jnn) = Jacobian(ji, (jpksed-1) * jpvode + jnn, (jpksed-1) * jpvode + jnn) - zb(jpksed)209 ENDIF210 211 265 IF( ln_timing ) CALL timing_stop('sed_mat_dsrjac') 212 266 213 267 END SUBROUTINE sed_mat_dsrjac 214 268 215 216 217 218 SUBROUTINE sed_mat_btb( nlev, nvar, psol, preac, dtsed_in ) 269 SUBROUTINE sed_mat_btbi( nvar, psol, preac, dtsed_in ) 270 !!--------------------------------------------------------------------- 271 !! *** ROUTINE sed_mat_btb *** 272 !! 273 !! ** Purpose : solves tridiagonal system of linear equations 274 !! 275 !! ** Method : 276 !! 1 - computes left hand side of linear system of equations 277 !! for dissolution reaction 278 !! 279 !! 280 !! 2 - forward/backward substitution. 281 !! 282 !! History : 283 !! ! 04-10 (N. Emprin, M. Gehlen ) original 284 !! ! 06-04 (C. Ethe) Module Re-organization 285 !!---------------------------------------------------------------------- 286 !! * Arguments 287 INTEGER , INTENT(in) :: nvar ! number of sediment levels 288 289 REAL(wp), DIMENSION(jpoce,jpksed,nvar), INTENT(inout) :: & 290 psol, preac 291 292 REAL(wp), INTENT(in) :: dtsed_in 293 294 !---Local declarations 295 INTEGER :: & 296 ji, jk, jn 297 298 REAL(wp) :: & 299 aplus,aminus , & 300 rplus,rminus , & 301 dxplus,dxminus 302 303 REAL(wp), DIMENSION(jpksed) :: za, zb, zc 304 REAL(wp), DIMENSION(jpksed) :: zr, zgamm 305 REAL(wp) :: zbet 306 307 !---------------------------------------------------------------------- 308 309 ! Computation left hand side of linear system of 310 ! equations for dissolution reaction 311 !--------------------------------------------- 312 IF( ln_timing ) CALL timing_start('sed_mat_btbi') 313 314 ! first sediment level 315 DO ji = 1, jpoce 316 aplus = ( por1(2) + por1(3) ) / 2.0 317 dxplus = ( dz(2) + dz(3) ) / 2. 318 rplus = ( dtsed_in / vols(2) ) * db(ji,2) * aplus / dxplus 319 za(2) = 0. 320 zb(2) = 1. + rplus 321 zc(2) = -rplus 322 323 DO jk = 3, jpksed - 1 324 aminus = ( por1(jk-1) + por1(jk) ) * 0.5 325 aminus = ( ( vols(jk-1) / dz(jk-1) ) + ( vols(jk) / dz(jk) ) ) / 2. 326 dxminus = ( dz(jk-1) + dz(jk) ) / 2. 327 rminus = ( dtsed_in / vols(jk) ) * db(ji,jk-1) * aminus / dxminus 328 ! 329 aplus = ( por1(jk) + por1(jk+1) ) * 0.5 330 dxplus = ( dz(jk) + dz(jk+1) ) / 2. 331 rplus = ( dtsed_in / vols(jk) ) * db(ji,jk) * aplus / dxplus 332 ! 333 za(jk) = -rminus 334 zb(jk) = 1. + rminus + rplus 335 zc(jk) = -rplus 336 337 ENDDO 338 339 aminus = ( por1(jpksed-1) + por1(jpksed) ) * 0.5 340 dxminus = ( dz(jpksed-1) + dz(jpksed) ) / 2. 341 rminus = ( dtsed_in / vols(jpksed) ) * db(ji,jpksed-1) * aminus / dxminus 342 ! 343 za(jpksed) = -rminus 344 zb(jpksed) = 1. + rminus 345 zc(jpksed) = 0. 346 347 ! solves tridiagonal system of linear equations 348 ! ----------------------------------------------- 349 DO jn = 1, nvar 350 zr(:) = psol(ji,:,jn) 351 zbet = zb(2) - preac(ji,2,jn) * dtsed_in 352 psol(ji,2,jn) = zr(2) / zbet 353 ! 354 DO jk = 3, jpksed 355 zgamm(jk) = zc(jk-1) / zbet 356 zbet = zb(jk) - preac(ji,jk,jn) * dtsed_in - za(jk) * zgamm(jk) 357 psol(ji,jk,jn) = ( zr(jk) - za(jk) * psol(ji,jk-1,jn) ) / zbet 358 ENDDO 359 ! 360 DO jk = jpksed - 1, 2, -1 361 psol(ji,jk,jn) = psol(ji,jk,jn) - zgamm(jk+1) * psol(ji,jk+1,jn) 362 ENDDO 363 END DO 364 END DO 365 ! 366 IF( ln_timing ) CALL timing_stop('sed_mat_btbi') 367 368 END SUBROUTINE sed_mat_btbi 369 370 371 SUBROUTINE sed_mat_btb( nvar, accmask ) 219 372 !!--------------------------------------------------------------------- 220 373 !! *** ROUTINE sed_mat_btb *** … … 234 387 !! * Arguments 235 388 INTEGER , INTENT(in) :: & 236 nlev, nvar ! number of sediment levels 237 238 REAL(wp), DIMENSION(jpoce,nlev,nvar), INTENT(inout) :: & 239 psol, preac 240 241 REAL(wp), INTENT(in) :: dtsed_in 389 nvar ! number of sediment levels 390 INTEGER, DIMENSION(jpoce) :: accmask 242 391 243 392 !---Local declarations 244 INTEGER :: & 245 ji, jk, jn 393 INTEGER :: ji, jk, jn 246 394 247 395 REAL(wp) :: & … … 250 398 dxplus,dxminus 251 399 252 REAL(wp), DIMENSION(nlev) :: za, zb, zc 253 REAL(wp), DIMENSION(jpoce,nlev) :: zr 254 REAL(wp), DIMENSION(nmax) :: zgamm 255 REAL(wp) :: zbet 400 REAL(wp), DIMENSION(jpksed) :: za, zb, zc 256 401 257 402 !---------------------------------------------------------------------- 258 403 259 ! Computation left hand side of linear system of 260 ! equations for dissolution reaction 261 !--------------------------------------------- 262 263 404 ! Computation left hand side of linear system of 405 ! equations for dissolution reaction 406 !--------------------------------------------- 264 407 IF( ln_timing ) CALL timing_start('sed_mat_btb') 265 408 266 ! first sediment level 409 ! first sediment level 410 jn = nvar 267 411 DO ji = 1, jpoce 412 IF (accmask(ji) == 0) THEN 268 413 aplus = ( por1(2) + por1(3) ) / 2.0 269 414 dxplus = ( dz(2) + dz(3) ) / 2. 270 rplus = ( dtsed_in / vols(2) ) * db(ji,2) * aplus / dxplus271 272 za( 1) = 0.273 zb( 1) = 1. +rplus274 zc( 1) = -rplus275 276 DO jk = 2, nlev- 1277 aminus = ( por1(jk ) + por1(jk+1) ) * 0.5278 aminus = ( ( vols(jk ) / dz(jk) ) + ( vols(jk+1) / dz(jk+1) ) ) / 2.279 dxminus = ( dz(jk ) + dz(jk+1) ) / 2.280 rminus = ( dtsed_in / vols(jk+1) ) * db(ji,jk) * aminus / dxminus415 rplus = ( 1.0 / vols(2) ) * db(ji,2) * aplus / dxplus * rads1sol(2,jn) 416 417 za(2) = 0. 418 zb(2) = rplus 419 zc(2) = -rplus 420 421 DO jk = 3, jpksed - 1 422 aminus = ( por1(jk-1) + por1(jk) ) * 0.5 423 aminus = ( ( vols(jk-1) / dz(jk-1) ) + ( vols(jk) / dz(jk) ) ) / 2. 424 dxminus = ( dz(jk-1) + dz(jk) ) / 2. 425 rminus = ( 1.0 / vols(jk) ) * db(ji,jk-1) * aminus / dxminus * rads1sol(jk,jn) 281 426 ! 282 aplus = ( por1(jk +1) + por1(jk+2) ) * 0.5283 dxplus = ( dz(jk +1) + dz(jk+2) ) / 2.284 rplus = ( dtsed_in / vols(jk+1) ) * db(ji,jk+1) * aplus / dxplus427 aplus = ( por1(jk) + por1(jk+1) ) * 0.5 428 dxplus = ( dz(jk) + dz(jk+1) ) / 2. 429 rplus = ( 1.0 / vols(jk) ) * db(ji,jk) * aplus / dxplus * rads1sol(jk,jn) 285 430 ! 286 431 za(jk) = -rminus 287 zb(jk) = 1. +rminus + rplus432 zb(jk) = rminus + rplus 288 433 zc(jk) = -rplus 289 434 290 435 ENDDO 291 436 292 aminus = ( por1( nlev) + por1(nlev+1) ) * 0.5293 dxminus = ( dz( nlev) + dz(nlev+1) ) / 2.294 rminus = ( dtsed_in / vols(nlev+1) ) * db(ji,nlev) * aminus / dxminus437 aminus = ( por1(jpksed-1) + por1(jpksed) ) * 0.5 438 dxminus = ( dz(jpksed-1) + dz(jpksed) ) / 2. 439 rminus = ( 1.0 / vols(jpksed) ) * db(ji,jpksed-1) * aminus / dxminus * rads1sol(jpksed,jn) 295 440 ! 296 za( nlev) = -rminus297 zb( nlev) = 1. +rminus298 zc( nlev) = 0.441 za(jpksed) = -rminus 442 zb(jpksed) = rminus 443 zc(jpksed) = 0. 299 444 300 445 ! solves tridiagonal system of linear equations 301 446 ! ----------------------------------------------- 302 DO jn = 1, nvar 303 DO jk = 1, nlev 304 zr(ji,jk) = psol(ji,jk,jn) 305 END DO 306 zbet = zb(1) - preac(ji,1,jn) * dtsed_in 307 psol(ji,1,jn) = zr(ji,1) / zbet 308 ! 309 DO jk = 2, nlev 310 zgamm(jk) = zc(jk-1) / zbet 311 zbet = zb(jk) - preac(ji,jk,jn) * dtsed_in - za(jk) * zgamm(jk) 312 psol(ji,jk,jn) = ( zr(ji,jk) - za(jk) * psol(ji,jk-1,jn) ) / zbet 313 ENDDO 314 ! 315 DO jk = nlev - 1, 1, -1 316 psol(ji,jk,jn) = psol(ji,jk,jn) - zgamm(jk+1) * psol(ji,jk+1,jn) 317 ENDDO 318 END DO 447 pwcpa(ji,2,jn) = pwcpa(ji,2,jn) - ( zc(2) * pwcp(ji,3,jn) + zb(2) * pwcp(ji,2,jn) ) 448 DO jk = 3, jpksed-1 449 pwcpa(ji,jk,jn) = pwcpa(ji,jk,jn) - ( zc(jk) * pwcp(ji,jk+1,jn) + za(jk) * pwcp(ji,jk-1,jn) & 450 & + zb(jk) * pwcp(ji,jk,jn) ) 451 ENDDO 452 pwcpa(ji,jpksed,jn) = pwcpa(ji,jpksed,jn) - ( za(jpksed) * pwcp(ji,jpksed-1,jn) & 453 & + zb(jpksed) * pwcp(ji,jpksed,jn) ) 454 ! 455 ENDIF 319 456 END DO 320 457 ! 321 458 IF( ln_timing ) CALL timing_stop('sed_mat_btb') 322 323 459 324 460 END SUBROUTINE sed_mat_btb 325 461 326 SUBROUTINE sed_mat_dsri( nvar, preac, psms, dtsed_in ) 462 SUBROUTINE sed_mat_btbjac( nvar, NEQ, NROWPD, jacvode, accmask ) 463 !!--------------------------------------------------------------------- 464 !! *** ROUTINE sed_mat_btb *** 465 !! 466 !! ** Purpose : solves tridiagonal system of linear equations 467 !! 468 !! ** Method : 469 !! 1 - computes left hand side of linear system of equations 470 !! for dissolution reaction 471 !! 472 !! 2 - forward/backward substitution. 473 !! 474 !! History : 475 !! ! 04-10 (N. Emprin, M. Gehlen ) original 476 !! ! 06-04 (C. Ethe) Module Re-organization 477 !!---------------------------------------------------------------------- 478 !! * Arguments 479 INTEGER , INTENT(in) :: nvar, NEQ, NROWPD ! number of variable 480 REAL, DIMENSION(jpoce,NROWPD,NEQ), INTENT(inout) :: jacvode 481 INTEGER, DIMENSION(jpoce), INTENT(in) :: accmask 482 483 !---Local declarations 484 INTEGER :: ji, jk, jn, jnn, jni, jnj ,jnij 485 486 REAL(wp) :: & 487 aplus,aminus , & 488 rplus,rminus , & 489 dxplus,dxminus 490 491 REAL(wp), DIMENSION(jpksed) :: za, zb, zc 492 493 !---------------------------------------------------------------------- 494 495 ! Computation left hand side of linear system of 496 ! equations for dissolution reaction 497 !--------------------------------------------- 498 IF( ln_timing ) CALL timing_start('sed_mat_btbjac') 499 500 ! first sediment level 501 jn = nvar 502 DO ji = 1, jpoce 503 IF (accmask(ji) == 0) THEN 504 aplus = ( por1(2) + por1(3) ) / 2.0 505 dxplus = ( dz(2) + dz(3) ) / 2. 506 rplus = ( 1.0 / vols(2) ) * db(ji,2) * aplus / dxplus * rads1sol(2,jn) 507 508 za(2) = 0. 509 zb(2) = rplus 510 zc(2) = -rplus 511 512 DO jk = 3, jpksed - 1 513 aminus = ( por1(jk-1) + por1(jk) ) * 0.5 514 aminus = ( ( vols(jk-1) / dz(jk-1) ) + ( vols(jk) / dz(jk) ) ) / 2. 515 dxminus = ( dz(jk-1) + dz(jk) ) / 2. 516 rminus = ( 1.0 / vols(jk) ) * db(ji,jk-1) * aminus / dxminus * rads1sol(jk,jn) 517 ! 518 aplus = ( por1(jk) + por1(jk+1) ) * 0.5 519 dxplus = ( dz(jk) + dz(jk+1) ) / 2. 520 rplus = ( 1.0 / vols(jk) ) * db(ji,jk) * aplus / dxplus * rads1sol(jk,jn) 521 ! 522 za(jk) = -rminus 523 zb(jk) = rminus + rplus 524 zc(jk) = -rplus 525 526 ENDDO 527 528 aminus = ( por1(jpksed-1) + por1(jpksed) ) * 0.5 529 dxminus = ( dz(jpksed-1) + dz(jpksed) ) / 2. 530 rminus = ( 1.0 / vols(jpksed) ) * db(ji,jpksed-1) * aminus / dxminus * rads1sol(jpksed,jn) 531 ! 532 za(jpksed) = -rminus 533 zb(jpksed) = rminus 534 zc(jpksed) = 0. 535 536 ! solves tridiagonal system of linear equations 537 ! ----------------------------------------------- 538 jnn = isvode(jn) 539 jni = jpvode + jnn 540 jnij = jpvode + 1 541 jacvode(ji, jnij, jni) = jacvode(ji,jnij,jni) - zb(2) 542 jnj = 2 * jpvode + jnn 543 jnij = jni - jnj + jpvode + 1 544 jacvode(ji, jnij, jnj) = jacvode(ji, jnij, jnj) -zc(2) 545 DO jk = 3, jpksed-1 546 jni = (jk-1) * jpvode + jnn 547 jnj = (jk-2) * jpvode + jnn 548 jnij = jni - jnj + jpvode + 1 549 jacvode(ji, jnij, jnj) = jacvode(ji, jnij, jnj) - za(jk) 550 jnj = (jk-1) * jpvode + jnn 551 jnij = jni - jnj + jpvode + 1 552 jacvode(ji, jnij, jnj) = jacvode(ji, jnij, jnj) - zb(jk) 553 jnj = (jk) * jpvode + jnn 554 jnij = jni - jnj + jpvode + 1 555 jacvode(ji, jnij, jnj) = jacvode(ji, jnij, jnj) - zc(jk) 556 ENDDO 557 jni = (jpksed-1) * jpvode + jnn 558 jnj = (jpksed-2) * jpvode + jnn 559 jnij = jni - jnj + jpvode + 1 560 jacvode(ji, jnij, jnj) = jacvode(ji, jnij, jnj) - za(jpksed) 561 jnij = jpvode + 1 562 jacvode(ji, jnij, jni) = jacvode(ji, jnij, jni) - zb(jpksed) 563 ! 564 ENDIF 565 END DO 566 ! 567 IF( ln_timing ) CALL timing_stop('sed_mat_btbjac') 568 569 END SUBROUTINE sed_mat_btbjac 570 571 572 SUBROUTINE sed_mat_dsri( nvar, preac, psms, dtsed_in, psol ) 327 573 !!--------------------------------------------------------------------- 328 574 !! *** ROUTINE sed_mat_dsr *** … … 352 598 REAL(wp), DIMENSION(jpoce,jpksed), INTENT(in ) :: preac ! reaction rates 353 599 REAL(wp), DIMENSION(jpoce,jpksed), INTENT(in ) :: psms ! reaction rates 600 REAL(wp), DIMENSION(jpoce,jpksed), INTENT(inout) :: psol ! reaction rates 354 601 REAL(wp), INTENT(in) :: dtsed_in 355 602 … … 357 604 INTEGER :: ji, jk, jn 358 605 REAL(wp), DIMENSION(jpoce,jpksed) :: za, zb, zc, zr 359 REAL(wp), DIMENSION(jpoce) :: zbet 360 REAL(wp), DIMENSION(jpoce,nmax) :: zgamm 361 362 REAL(wp) :: aplus,aminus 606 REAL(wp), DIMENSION(jpoce) :: zbet 607 REAL(wp), DIMENSION(jpoce,jpksed) :: zgamm 608 363 609 REAL(wp) :: rplus,rminus 364 REAL(wp) :: dxplus,dxminus365 366 610 !---------------------------------------------------------------------- 367 611 … … 371 615 ! equations for dissolution reaction 372 616 !--------------------------------------------- 373 374 375 617 jn = nvar 376 618 ! first sediment level 377 619 DO ji = 1, jpoce 378 aplus = ( por(1) + por(2) ) * 0.5 379 dxplus = ( dz3d(ji,1) + dz3d(ji,2) ) / 2. 380 rplus = ( dtsed_in / ( volw3d(ji,1) ) ) * diff(ji,1,jn) * aplus / dxplus 620 rplus = dtsed_in * apluss(ji,1) * diff(ji,1,jn) * radssol(1,jn) 381 621 382 622 za(ji,1) = 0. … … 387 627 DO jk = 2, jpksed - 1 388 628 DO ji = 1, jpoce 389 aminus = ( por(jk-1) + por(jk) ) * 0.5 390 dxminus = ( dz3d(ji,jk-1) + dz3d(ji,jk) ) / 2. 391 392 aplus = ( por(jk+1) + por(jk) ) * 0.5 393 dxplus = ( dz3d(ji,jk) + dz3d(ji,jk+1) ) / 2 394 ! 395 rminus = ( dtsed_in / volw3d(ji,jk) ) * diff(ji,jk-1,jn) * aminus / dxminus 396 rplus = ( dtsed_in / volw3d(ji,jk) ) * diff(ji,jk,jn) * aplus /dxplus 629 rminus = dtsed_in * aminuss(ji,jk) * diff(ji,jk-1,jn) * radssol(jk,jn) 630 rplus = dtsed_in * apluss (ji,jk) * diff(ji,jk,jn) * radssol(jk,jn) 397 631 ! 398 632 za(ji,jk) = -rminus … … 403 637 404 638 DO ji = 1, jpoce 405 aminus = ( por(jpksed-1) + por(jpksed) ) *0.5 406 dxminus = ( dz3d(ji,jpksed-1) + dz3d(ji,jpksed) ) / 2. 407 rminus = ( dtsed_in / volw3d(ji,jpksed) ) * diff(ji,jpksed-1,jn) * aminus/ dxminus 639 rminus = dtsed_in * aminuss(ji,jpksed) * diff(ji,jpksed-1,jn) * radssol(jpksed,jn) 408 640 ! 409 641 za(ji,jpksed) = -rminus … … 416 648 ! ----------------------------------------------- 417 649 418 zr (:,:) = p wcp(:,:,jn) + psms(:,:)419 zb (:,:) = zb(:,:) - preac(:,:) 650 zr (:,:) = psol(:,:) + psms(:,:) * dtsed_in 651 zb (:,:) = zb(:,:) - preac(:,:) * dtsed_in 420 652 zbet(: ) = zb(:,1) 421 p wcp(:,1,jn) = zr(:,1) / zbet(:)653 psol(:,1) = zr(:,1) / zbet(:) 422 654 423 655 ! … … 426 658 zgamm(ji,jk) = zc(ji,jk-1) / zbet(ji) 427 659 zbet(ji) = zb(ji,jk) - za(ji,jk) * zgamm(ji,jk) 428 p wcp(ji,jk,jn) = ( zr(ji,jk) - za(ji,jk) * pwcp(ji,jk-1,jn) ) / zbet(ji)660 psol(ji,jk) = ( zr(ji,jk) - za(ji,jk) * psol(ji,jk-1) ) / zbet(ji) 429 661 END DO 430 662 ENDDO … … 432 664 DO jk = jpksed - 1, 1, -1 433 665 DO ji = 1, jpoce 434 p wcp(ji,jk,jn) = pwcp(ji,jk,jn) - zgamm(ji,jk+1) * pwcp(ji,jk+1,jn)666 psol(ji,jk) = psol(ji,jk) - zgamm(ji,jk+1) * psol(ji,jk+1) 435 667 END DO 436 668 ENDDO -
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/sedorg.F90
r15127 r15297 7 7 !! * Modules used 8 8 USE sed ! sediment global variable 9 USE sed_oce 9 10 USE sedini 10 11 USE sedmat … … 59 60 preac(:,:) = 0.0_wp 60 61 psms (:,:) = rearatpom(:,:) 61 CALL sed_mat_dsri( jwdic, preac, psms, dtsed )62 CALL sed_mat_dsri( jwdic, preac, psms, dtsed, pwcp(:,:,jwdic) ) 62 63 63 64 ! Silicate in pore water 64 65 psms (:,:) = 0.0_wp 65 CALL sed_mat_dsri( jwsil, preac, psms, dtsed )66 CALL sed_mat_dsri( jwsil, preac, psms, dtsed, pwcp(:,:,jwsil) ) 66 67 67 68 ! Iron ligands in pore water 68 69 psms (:,:) = ratligc * rearatpom(:,:) 69 preac(:,:) = -reac_ligc * dtsed70 CALL sed_mat_dsri( jwlgw, preac, psms, dtsed )70 preac(:,:) = -reac_ligc 71 CALL sed_mat_dsri( jwlgw, preac, psms, dtsed, pwcp(:,:,jwlgw) ) 71 72 72 73 IF( ln_timing ) CALL timing_stop('sed_org') -
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/sedrst.F90
r15075 r15297 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=3) :: cdcomp 44 45 !!---------------------------------------------------------------------- 45 46 ! … … 79 80 IF(lwp) WRITE(numsed,*) & 80 81 ' open sed restart.output NetCDF file: ',TRIM(clpath)//clname 81 CALL iom_open( TRIM(clpath)//TRIM(clname), numrsw, ldwrt = .TRUE., kdlev = jpksed ) 82 cdcomp ='SED' 83 CALL iom_open( TRIM(clpath)//TRIM(clname), numrsw, ldwrt = .TRUE., kdlev = jpksed, cdcomp = cdcomp ) 82 84 lrst_sed = .TRUE. 83 85 ENDIF -
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/sedsol.F90
r15127 r15297 8 8 !! * Modules used 9 9 USE sed ! sediment global variable 10 USE sed_oce 10 11 USE sedini 11 12 USE sedfunc … … 15 16 USE sedco3 16 17 USE sedmat 17 # if ! defined key_agrif 18 USE trosk 19 #endif 18 USE sedorg 19 USE tros 20 20 USE lib_mpp ! distribued memory computing library 21 21 USE lib_fortran … … 57 57 ! --- local variables 58 58 INTEGER :: ji, jk, js, jw, jn, neq ! dummy looop indices 59 REAL(wp), DIMENSION( jpoce, jpvode * jpksed ) :: ZXIN, FVAL 59 60 REAL(wp), DIMENSION(jpoce,jpksed) :: preac 60 # if ! defined key_agrif61 REAL(wp), DIMENSION( jpvode * jpksed ) :: zxin62 61 INTEGER :: JINDEX, ITOL, IJAC, MLJAC, IMAX 63 INTEGER :: MUJAC,LE1, LJAC, L IWORK, LWORK62 INTEGER :: MUJAC,LE1, LJAC, LWORK 64 63 INTEGER :: IDID, NMAXSTP, ROSM 65 REAL(wp) :: X, XEND, H, STEPMIN 66 REAL(wp), DIMENSION(1) :: RTOL, ATOL 64 REAL(wp) :: X, XEND 65 REAL(wp),DIMENSION(jpoce) :: H 66 INTEGER, DIMENSION(jpoce) :: accmask 67 REAL(wp), DIMENSION(jpvode * jpksed) :: RTOL, ATOL 67 68 REAL(wp), POINTER :: WORK(:) 68 INTEGER , POINTER :: IWORK(:) 69 INTEGER, DIMENSION(3) :: ISTAT 70 REAL(wp), DIMENSION(2) :: RSTAT 71 #endif 69 INTEGER, DIMENSION(jpoce,3) :: ISTAT 70 REAL(wp), DIMENSION(jpoce,2) :: RSTAT 72 71 !! 73 72 !!---------------------------------------------------------------------- … … 81 80 ENDIF 82 81 ! ! 83 dens_mol_wgt(1:jpsol) = denssol / mol_wgt(1:jpsol)84 stepros(:) = dtsed85 !86 82 ENDIF 87 83 … … 108 104 DO js = 1, jpsol 109 105 DO ji = 1, jpoce 110 IF (raintg(ji) .ne. 0) THEN 111 solcp(ji,2,js) = solcp(ji,2,js) + & 112 & ( rainrg(ji,js) / raintg(ji) ) * ( dzdep(ji) / dz3d(ji,2) ) 106 IF (dzdep(ji) .ne. 0) THEN 107 solcp(ji,2,js) = solcp(ji,2,js) + rainrg(ji,js) * dtsed / ( por1(2) * dz3d(ji,2) ) 113 108 ! rainrm are temporary cancel 114 rainrm(ji,js) = 0.115 109 ENDIF 116 110 END DO … … 148 142 diff(:,:,jwso4) = diff(:,:,jwso4) * ( 1.0 + irrig(:,:) ) 149 143 diff(:,:,jwlgw) = diff(:,:,jwlgw) * ( 1.0 + irrig(:,:) ) 150 DO jk = 1, jpksed 151 diff(:,jk,jwfe2) = diff(:,jk,jwfe2) * ( 1.0 + 0.1 * irrig(:,jk) ) * radsfe2(jk) 152 diff(:,jk,jwnh4) = diff(:,jk,jwnh4) * ( 1.0 + irrig(:,jk) ) * radsnh4(jk) 153 END DO 144 diff(:,:,jwnh4) = diff(:,:,jwnh4) * ( 1.0 + irrig(:,:) ) 145 diff(:,:,jwfe2) = diff(:,:,jwfe2) * ( 1.0 + 0.1 * irrig(:,:) ) 154 146 155 147 ! Conversion of volume units … … 157 149 DO js = 1, jpsol 158 150 DO jk = 1, jpksed 159 volc(:,jk,js) = ( vols3d(:,jk) * dens_mol_wgt(js) ) / &151 volc(:,jk,js) = ( vols3d(:,jk) / mol_wgt(js) ) / & 160 152 & ( volw3d(:,jk) * 1.e-3 ) 161 153 ENDDO 162 154 ENDDO 163 155 156 ! Compute coefficients commonly used in diffusion 157 CALL sed_mat_coef 164 158 ! Apply bioturbation and compute the impact of the slow SMS on species 165 159 CALL sed_btb( kt ) 166 160 ! Recompute pH after bioturbation and slow chemistry 167 161 CALL sed_co3( kt ) 168 # if ! defined key_agrif 162 169 163 ! The following part deals with the stiff ODEs 170 164 ! This is the expensive part of the code and should be carefully … … 176 170 ! Brent, Powell's hybrid method, ... 177 171 ! --------------------------------------------------------------------- 178 ROSM = 2179 172 NEQ = jpvode * jpksed 180 173 XEND = dtsed 181 RTOL = 0.005182 ATOL = 0.005183 ITOL = 0174 RTOL = rosrtol 175 ATOL = rosatol 176 ITOL = 1 184 177 IJAC = 1 178 DO jn = 1, NEQ 179 js = jarr(jn,2) 180 IF (js == jwfe2) ATOL(jn) = rosatol / 100.0 181 END DO 185 182 MLJAC = jpvode 186 183 MUJAC = jpvode 187 184 LE1 = 2*MLJAC+MUJAC+1 188 185 LJAC = MLJAC+MUJAC+1 189 LIWORK = NEQ + 2190 186 LWORK = NEQ*(LJAC+LE1+8) + 5 191 ALLOCATE(WORK(LWORK), IWORK(LIWORK) ) 192 NMAXSTP = 0 193 STEPMIN = dtsed 194 195 DO ji = 1, jpoce 196 X = 0.0 197 H = stepros(ji) 198 WORK = 0. 199 IWORK = 0 200 IWORK(2) = 6 201 202 ! Put all the species in one local array (nb of tracers * vertical 203 ! dimension 204 jindex = ji 205 DO jn = 1, NEQ 206 jk = jarr(jn,1) 207 js = jarr(jn,2) 208 IF (js <= jpwat) THEN 209 zxin(jn) = pwcp(ji,jk,js) * 1E6 210 ELSE 211 zxin(jn) = solcp(ji,jk,js-jpwat) * 1E6 212 ENDIF 213 END DO 214 215 ! Set options for VODE : banded matrix. SParse option is much more 216 ! expensive except if one computes the sparse Jacobian explicitly 217 ! To speed up the computation, one way is to reduce ATOL and RTOL 218 ! which may be a good option at the beginning of the simulations 219 ! during the spin up 220 ! ---------------------------------------------------------------- 221 CALL ROSK(ROSM, NEQ,JINDEX,X,zxin,XEND,H, & 222 RTOL,ATOL,ITOL, sed_jac,IJAC,MLJAC,MUJAC, & 223 WORK,LWORK,IWORK,LIWORK,IDID,ISTAT,RSTAT) 224 225 DO jn = 1, NEQ 226 jk = jarr(jn,1) 227 js = jarr(jn,2) 228 IF (js <= jpwat) THEN 229 pwcp(ji,jk,js) = zxin(jn) * 1E-6 230 ELSE 231 solcp(ji,jk,js-jpwat) = zxin(jn) * 1E-6 232 ENDIF 233 END DO 234 NMAXSTP = MAX(NMAXSTP,ISTAT(3)) 235 STEPMIN = MIN(STEPMIN, RSTAT(1) ) 236 stepros(ji) = RSTAT(1) 237 END DO 238 239 DEALLOCATE(WORK, IWORK) 240 #endif 187 ALLOCATE(WORK(LWORK) ) 188 189 X = 0.0 190 H(:) = dtsed 191 WORK = 0. 192 193 ! Put all the species in one local array (nb of tracers * vertical 194 ! dimension 195 DO jn = 1, NEQ 196 jk = jarr(jn,1) 197 js = jarr(jn,2) 198 IF (js <= jpwat) THEN 199 zxin(:,jn) = pwcp(:,jk,js) * 1E6 200 ELSE 201 zxin(:,jn) = solcp(:,jk,js-jpwat) * 1E6 202 ENDIF 203 END DO 204 205 ! Set options for VODE : banded matrix. SParse option is much more 206 ! expensive except if one computes the sparse Jacobian explicitly 207 ! To speed up the computation, one way is to reduce ATOL and RTOL 208 ! which may be a good option at the beginning of the simulations 209 ! during the spin up 210 ! ---------------------------------------------------------------- 211 CALL ROSK(NROSORDER, NEQ,X,zxin,XEND,H,RTOL,ATOL,ITOL, sed_jac, & 212 & MLJAC,MUJAC,WORK,LWORK,IDID,ISTAT,RSTAT) 213 214 accmask(:) = 0.0 215 CALL sed_func( NEQ, ZXIN, FVAL, ACCMASK ) 216 217 DO jn = 1, NEQ 218 jk = jarr(jn,1) 219 js = jarr(jn,2) 220 IF (js <= jpwat) THEN 221 pwcp(:,jk,js) = zxin(:,jn) * 1E-6 222 ELSE 223 solcp(:,jk,js-jpwat) = zxin(:,jn) * 1E-6 224 ENDIF 225 END DO 226 rstepros(:) = ISTAT(:,3) 227 228 DEALLOCATE( WORK ) 241 229 242 230 preac(:,:) = 0. 243 CALL sed_mat_dsri( jwalk, preac, pwcpa(:,:,jwalk), dtsed ) 244 CALL sed_mat_dsri( jwpo4, preac, pwcpa(:,:,jwpo4), dtsed ) 231 CALL sed_mat_dsri( jwpo4, preac, pwcpa(:,:,jwpo4), dtsed, pwcp(:,:,jwpo4) ) 232 CALL sed_mat_dsri( jwalk, preac, pwcpa(:,:,jwalk), dtsed, pwcp(:,:,jwalk) ) 233 234 CALL sed_org( kt ) 245 235 246 236 !---------------------------------- … … 256 246 ENDDO 257 247 258 !----------------------------------------------------------------------- 259 ! 2/ Det of new rainrg taking account of the new weight fraction 260 ! obtained 261 ! in dz3d(2) after diffusion/reaction (react/diffu are also in 262 ! dzdep!) 263 ! This new rain (rgntg rm) will be used in advection/burial routine 264 !------------------------------------------------------------------------ 265 DO js = 1, jpsol 266 rainrg(:,js) = raintg(:) * solcp(:,2,js) 267 rainrm(:,js) = rainrg(:,js) / mol_wgt(js) 268 ENDDO 269 270 ! New raintg 271 raintg(:) = 0. 272 DO js = 1, jpsol 273 raintg(:) = raintg(:) + rainrg(:,js) 274 ENDDO 248 ! 2. Change of previous solid fractions (due to volum changes) for k=2 249 !--------------------------------------------------------------------- 250 DO js = 1, jpsol 251 solcp(:,2,js) = solcp(:,2,js) * dz3d(:,2) / dz(2) 252 END DO 275 253 276 254 !-------------------------------- … … 285 263 END SUBROUTINE sed_sol 286 264 265 SUBROUTINE JACFAC 266 267 END SUBROUTINE JACFAC 268 287 269 END MODULE sedsol -
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/sedwri.F90
r15075 r15297 121 121 ENDDO 122 122 123 DO jn = 1, jpdia2dsed - 1123 DO jn = 1, jpdia2dsed - 2 124 124 CALL unpack_arr( jpoce, flxsedi2d(1:jpi,1:jpj,jn), iarroce(1:jpoce), zflx(1:jpoce,jn) ) 125 125 END DO 126 126 127 127 zflx(:,1) = dzdep(:) / dtsed 128 CALL unpack_arr( jpoce, flxsedi2d(1:jpi,1:jpj,jpdia2dsed ), iarroce(1:jpoce), zflx(1:jpoce,1) )128 CALL unpack_arr( jpoce, flxsedi2d(1:jpi,1:jpj,jpdia2dsed-1), iarroce(1:jpoce), zflx(1:jpoce,1) ) 129 129 130 CALL unpack_arr( jpoce, flxsedi2d(1:jpi,1:jpj,jpdia2dsed), iarroce(1:jpoce), rstepros(1:jpoce) ) 130 131 ! 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 )132 ! CALL lbc_lnk( 'sedwri', trcsedi(:,:,:,:), 'T', 1._wp ) 133 ! CALL lbc_lnk( 'sedwri', flxsedi3d(:,:,:,:), 'T', 1._wp ) 134 ! CALL lbc_lnk( 'sedwri', flxsedi2d(:,:,:), 'T', 1._wp ) 134 135 135 136 ! Start writing data -
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/trcdmp_sed.F90
r15127 r15297 93 93 CALL trc_dta( kt, jl, ztrcdta ) ! read tracer data at nit000 94 94 ! 95 DO_2D( 1, 1, 1, 1)95 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 96 96 ikt = mbkt(ji,jj) 97 97 tr(ji,jj,ikt,jn,Kbb) = ztrcdta(ji,jj,ikt) + ( tr(ji,jj,ikt,jn,Kbb) - ztrcdta(ji,jj,ikt) ) & -
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/trosk.F90
r15127 r15297 1 MODULE trosk 2 # if ! defined key_agrif 1 MODULE TROS 3 2 !**************************************************************** 4 3 !* NUMERICAL SOLUTION OF A STIFF SYSTEM OF FIRST 0RDER ORDINARY * 5 4 !* DIFFERENTIAL EQUATIONS Y'=F(X,Y) BY ROSENBROCK METHOD. * 6 5 !* ------------------------------------------------------------ * 7 !* SAMPLE RUN: *8 !* Example #1: *9 !* (Solve set of differential equations (N=2): *10 !* F(1) = Y(1) * Y(2) + COS(X) - HALF * SIN(TWO * X) *11 !* F(2) = Y(1) * Y(1) + Y(2) * Y(2) - (ONE + SIN(X)) *12 !* Find values of F(1), F(2) at X=1.5). *13 !* *14 !* SOLUTION AT X= 1.50000000000000 *15 !* Y(1) = 0.12359935E+01 *16 !* Y(2) = -0.10494372E+00 *17 !* *18 !* LAST STEP SIZE = 4.150113101356574E-002 *19 !* ERROR CODE = 1 *20 !* *21 !* Example #2: *22 !* (Solve set of differential equations (N=5): *23 !* F(1) = Y(2) *24 !* F(2) = Y(3) *25 !* F(3) = Y(4) *26 !* F(4) = Y(5) *27 !* F(5) = (45.d0 * Y(3) * Y(4) * Y(5) - *28 !* 40.d0 * Y(4) * Y(4) * Y(4)) / (NINE * Y(3) * Y(3)) *29 !* Find values of F(1), F(2), ..., F(5) at X=1.5). *30 !* *31 !* SOLUTION AT X= 1.50000000000000 *32 !* Y(1) = 0.43639610E+01 *33 !* Y(2) = 0.40000000E+01 *34 !* Y(3) = 0.28284271E+01 *35 !* Y(4) = 0.14790900E-10 *36 !* Y(5) = -0.37712362E+01 *37 !* *38 !* LAST STEP SIZE = 3.825256949194526E-003 *39 !* ERROR CODE = 1 *40 6 !* ------------------------------------------------------------ * 41 7 !* Ref.: From Numath Library By Tuan Dang Trong in Fortran 77 * … … 45 11 !* (www.jpmoreau.fr) * 46 12 !**************************************************************** 13 USE timing 14 USE in_out_manager, ONLY : ln_timing ! I/O manager 47 15 USE sed 48 16 USE sedfunc 49 17 50 INTEGER, PRIVATE :: JINDIC 51 52 PUBLIC rosk 18 IMPLICIT NONE 19 PRIVATE 20 21 PUBLIC ROSK 22 23 INTEGER, ALLOCATABLE, SAVE, DIMENSION(:) :: NFCN, NJAC, NSTEP, NACCPT, NREJCT, NDEC, NSOL 24 53 25 54 26 !define example #1 55 27 INTERFACE 56 SUBROUTINE JAC(NEQ, JINDEX,X,Y,DFY,LDFY)28 SUBROUTINE JAC(NEQ,Y,DFY,LDFY,ACCMASK) 57 29 INTEGER, PARAMETER :: WP = KIND(1.0D0) 58 INTEGER NEQ, JINDEX, LDFY 59 REAL(WP), DIMENSION(LDFY,NEQ) :: DFY 60 INTENT(IN) :: NEQ,JINDEX,X,Y,LDFY 61 INTENT(OUT) :: DFY 30 INTEGER, INTENT(IN) :: NEQ, LDFY 31 REAL(WP), DIMENSION(:,:,:), INTENT(OUT) :: DFY 32 INTEGER , DIMENSION(:), INTENT(IN) :: ACCMASK 62 33 END SUBROUTINE JAC 63 34 END INTERFACE … … 67 38 68 39 !********************************************************************** 69 SUBROUTINE rosk(ROSM,N,JINDEX,X,Y,XEND,H, & 70 RTOL,ATOL,ITOL, & 71 JAC ,IJAC,MLJAC,MUJAC, & 72 WORK,LWORK,IWORK,LIWORK,IDID,ISTAT,RSTAT) 40 SUBROUTINE ROSK(ROSM,N,X,Y,XEND,H, RTOL,ATOL,ITOL, & 41 & JAC, MLJAC, MUJAC, WORK,LWORK,IDID,ISTAT,RSTAT) 73 42 ! --------------------------------------------------------------------- 74 43 ! NUMERICAL SOLUTION OF A STIFF (OR DIFFERENTIAL ALGEBRAIC) … … 128 97 ! JAC NAME (EXTERNAL) OF THE SUBROUTINE WHICH COMPUTES 129 98 ! THE PARTIAL DERIVATIVES OF F(X,Y) WITH RESPECT TO Y 130 ! (THIS ROUTINE IS ONLY CALLED IF IJAC=1; SUPPLY131 ! A DUMMY SUBROUTINE IN THE CASE IJAC=0).132 99 ! FOR IJAC=1, THIS SUBROUTINE MUST HAVE THE FORM: 133 100 ! SUBROUTINE JAC(N,X,Y,DFY,LDFY) … … 136 103 ! LDFY, THE COLUMN-LENGTH OF THE ARRAY, IS 137 104 ! FURNISHED BY THE CALLING PROGRAM. 138 ! IF (MLJAC.EQ.N) THE JACOBIAN IS SUPPOSED TO 139 ! BE FULL AND THE PARTIAL DERIVATIVES ARE 140 ! STORED IN DFY AS 141 ! DFY(I,J) = PARTIAL F(I) / PARTIAL Y(J) 142 ! ELSE, THE JACOBIAN IS TAKEN AS BANDED AND 105 ! THE JACOBIAN IS TAKEN AS BANDED AND 143 106 ! THE PARTIAL DERIVATIVES ARE STORED 144 107 ! DIAGONAL-WISE AS 145 108 ! DFY(I-J+MUJAC+1,J) = PARTIAL F(I) / PARTIAL Y(J). 146 109 ! 147 ! IJAC SWITCH FOR THE COMPUTATION OF THE JACOBIAN:148 ! IJAC=0: JACOBIAN IS COMPUTED INTERNALLY BY FINITE149 ! DIFFERENCES, SUBROUTINE "JAC" IS NEVER CALLED.150 ! IJAC=1: JACOBIAN IS SUPPLIED BY SUBROUTINE JAC.151 !152 110 ! MLJAC SWITCH FOR THE BANDED STRUCTURE OF THE JACOBIAN: 153 ! MLJAC=N: JACOBIAN IS A FULL MATRIX. THE LINEAR154 ! ALGEBRA IS DONE BY FULL-MATRIX GAUSS-ELIMINATION.155 111 ! 0<=MLJAC<N: MLJAC IS THE LOWER BANDWITH OF JACOBIAN 156 112 ! MATRIX (>= NUMBER OF NON-ZERO DIAGONALS BELOW … … 179 135 ! LWORK DECLARED LENGHT OF ARRAY "WORK". 180 136 ! 181 ! IWORK INTEGER WORKING SPACE OF LENGTH "LIWORK".182 ! "LIWORK" MUST BE AT LEAST N+2.183 !184 ! LIWORK DECLARED LENGHT OF ARRAY "IWORK".185 !186 137 ! ---------------------------------------------------------------------- 187 138 ! … … 193 144 ! FOR ZERO INPUT, THE CODE CHOOSES DEFAULT VALUES: 194 145 ! 195 ! IWORK(1) THIS IS THE MAXIMAL NUMBER OF ALLOWED STEPS.196 ! THE DEFAULT VALUE (FOR IWORK(1)=0) IS 100000.197 !198 ! IWORK(2) SWITCH FOR THE CHOICE OF THE COEFFICIENTS199 ! IF IWORK(2).EQ.1 METHOD OF SHAMPINE200 ! IF IWORK(2).EQ.2 METHOD GRK4T OF KAPS-RENTROP201 ! IF IWORK(2).EQ.3 METHOD GRK4A OF KAPS-RENTROP202 ! IF IWORK(2).EQ.4 METHOD OF VAN VELDHUIZEN (GAMMA=1/2)203 ! IF IWORK(2).EQ.5 METHOD OF VAN VELDHUIZEN ("D-STABLE")204 ! IF IWORK(2).EQ.6 AN L-STABLE METHOD205 ! THE DEFAULT VALUE (FOR IWORK(2)=0) IS IWORK(2)=2.206 !207 146 ! WORK(1) UROUND, THE ROUNDING UNIT, DEFAULT 1.D-16. 208 147 ! … … 237 176 ! DECLARATIONS 238 177 ! *** *** *** *** *** *** *** *** *** *** *** *** *** 239 IMPLICIT REAL(wp) (A-H,O-Z) 240 INTEGER :: JINDEX, ROSM 241 DIMENSION Y(N),ATOL(1),RTOL(1),WORK(LWORK),IWORK(LIWORK) 242 INTEGER , DIMENSION(3) :: ISTAT 243 REAL(wp), DIMENSION(2) :: RSTAT 244 LOGICAL JBAND,ARRET 178 INTEGER, INTENT(in) :: ROSM, N, ITOL, MLJAC, MUJAC, LWORK 179 REAL(wp), DIMENSION(1), INTENT(in) :: ATOL, RTOL 180 INTEGER, INTENT(inout) :: IDID 181 INTEGER , DIMENSION(jpoce,3), INTENT(out) :: ISTAT 182 REAL(wp), DIMENSION(jpoce,2), INTENT(out) :: RSTAT 183 184 INTEGER :: NMAX, LDJAC, LDE, IEYNEW, IEDY1, IEDY, IEAK1 185 INTEGER :: IEAK2, IEAK3, IEAK4, IEFX, IEJAC, IEE 186 INTEGER :: ISTORE 187 REAL(wp) :: UROUND, HMAX, FAC1, FAC2, FACREJ, XEND, X 188 REAL(wp), DIMENSION(jpoce) :: H 189 REAL(wp), DIMENSION(jpoce, N) :: Y 190 REAL(wp), DIMENSION(LWORK) :: WORK 191 LOGICAL ARRET 245 192 EXTERNAL JAC 246 COMMON/STAT/NFCN,NJAC,NSTEP,NACCPT,NREJCT,NDEC,NSOL247 193 ! -------------------------------------------------------------------- 248 194 ! --- COMMON STAT CAN BE USED FOR STATISTICS … … 261 207 ! SETTING THE PARAMETERS 262 208 ! *** *** *** *** *** *** *** 209 210 IF ( ln_timing ) CALL timing_start('rosk') 211 212 ALLOCATE (NFCN(jpoce), NJAC(jpoce), NSTEP(jpoce), NACCPT(jpoce), NREJCT(jpoce), NDEC(jpoce), NSOL(jpoce)) 213 263 214 NFCN=0 264 215 NJAC=0 … … 269 220 NSOL=0 270 221 ARRET=.FALSE. 271 JINDIC = JINDEX272 222 ! -------- NMAX , THE MAXIMAL NUMBER OF STEPS ----- 273 IF(IWORK(1).EQ.0)THEN 274 NMAX=100000 223 NMAX = 100000 224 ! -------- UROUND SMALLEST NUMBER SATISFYING 1.D0+UROUND>1.D0 225 IF(WORK(1) == 0.0)THEN 226 UROUND = 1.E-16 275 227 ELSE 276 NMAX=IWORK(1) 277 IF(NMAX.LE.0)THEN 278 WRITE(NUMSED,*)' WRONG INPUT IWORK(1)=',IWORK(1) 279 ARRET=.TRUE. 280 END IF 281 END IF 282 ! -------- METH COEFFICIENTS OF THE METHOD 283 IF(IWORK(2).EQ.0)THEN 284 METH=2 285 ELSE 286 METH=IWORK(2) 287 IF(METH.LE.0.OR.METH.GE.7)THEN 288 WRITE(NUMSED,*)' CURIOUS INPUT IWORK(2)=',IWORK(2) 289 ARRET=.TRUE. 290 END IF 291 END IF 292 ! -------- UROUND SMALLEST NUMBER SATISFYING 1.D0+UROUND>1.D0 293 IF(WORK(1).EQ.0.D0)THEN 294 UROUND=1.D-16 295 ELSE 296 UROUND=WORK(1) 297 IF(UROUND.LE.1.D-14.OR.UROUND.GE.1.D0)THEN 228 UROUND = WORK(1) 229 IF(UROUND <= 1.E-14 .OR. UROUND >= 1.0)THEN 298 230 WRITE(NUMSED,*)' COEFFICIENTS HAVE 16 DIGITS, UROUND=',WORK(1) 299 231 ARRET=.TRUE. … … 301 233 END IF 302 234 ! -------- MAXIMAL STEP SIZE 303 IF(WORK(2) .EQ.0.D0)THEN304 HMAX =XEND-X235 IF(WORK(2) == 0.0)THEN 236 HMAX = XEND-X 305 237 ELSE 306 HMAX =WORK(2)238 HMAX = WORK(2) 307 239 END IF 308 240 ! ------- FAC1,FAC2 PARAMETERS FOR STEP SIZE SELECTION 309 IF(WORK(3) .EQ.0.D0)THEN310 FAC1 =5.D0241 IF(WORK(3) == 0.0)THEN 242 FAC1 = 5.0_wp 311 243 ELSE 312 FAC1 =1.D0/WORK(3)244 FAC1 = 1.0/WORK(3) 313 245 END IF 314 IF(WORK(4) .EQ.0.D0)THEN315 FAC2 =1.D0/6.0D0246 IF(WORK(4) == 0.0)THEN 247 FAC2 = 1.0_wp / 6.0_wp 316 248 ELSE 317 FAC2 =1.D0/WORK(4)249 FAC2 = 1.0_wp / WORK(4) 318 250 END IF 319 251 ! ------- FACREJ FOR THE HUMP 320 IF(WORK(5) .EQ.0.D0)THEN321 FACREJ =0.1D0252 IF(WORK(5) == 0.0)THEN 253 FACREJ = 0.1_wp 322 254 ELSE 323 FACREJ=WORK(5) 324 END IF 325 ! --------- CHECK IF TOLERANCES ARE O.K. 326 IF (ITOL.EQ.0) THEN 327 IF (ATOL(1).LE.0.D0.OR.RTOL(1).LE.10.D0*UROUND) THEN 328 WRITE (NUMSED,*) ' TOLERANCES ARE TOO SMALL' 329 ARRET=.TRUE. 330 END IF 331 ELSE 332 DO 15 I=1,N 333 IF (ATOL(I).LE.0.D0.OR.RTOL(I).LE.10.D0*UROUND) THEN 334 WRITE (NUMSED,*) ' TOLERANCES(',I,') ARE TOO SMALL' 335 ARRET=.TRUE. 336 END IF 337 15 CONTINUE 255 FACREJ = WORK(5) 338 256 END IF 339 257 ! *** *** *** *** *** *** *** *** *** *** *** *** *** 340 258 ! COMPUTATION OF ARRAY ENTRIES 341 259 ! *** *** *** *** *** *** *** *** *** *** *** *** *** 342 ! ---- AUTONOMOUS, IMPLICIT, BANDED OR NOT ?343 JBAND=MLJAC.NE.N344 ARRET=.FALSE.345 260 ! -------- COMPUTATION OF THE ROW-DIMENSIONS OF THE 2-ARRAYS --- 346 261 ! -- JACOBIAN 347 IF(JBAND)THEN 348 LDJAC=MLJAC+MUJAC+1 349 ELSE 350 LDJAC=N 351 END IF 352 ! -- MATRIX E FOR LINEAR ALGEBRA 353 IF(JBAND)THEN 354 LDE=2*MLJAC+MUJAC+1 355 ELSE 356 LDE=N 357 END IF 262 LDJAC=MLJAC+MUJAC+1 263 LDE=2*MLJAC+MUJAC+1 358 264 ! ------- PREPARE THE ENTRY-POINTS FOR THE ARRAYS IN WORK ----- 359 265 IEYNEW=6 … … 366 272 IEFX =IEAK4+N 367 273 IEJAC=IEFX +N 368 IEMAS=IEJAC+N*LDJAC 369 IEE =IEMAS 274 IEE =IEJAC+N*LDJAC 370 275 ! ------ TOTAL STORAGE REQUIREMENT ----------- 371 276 ISTORE=IEE+N*LDE-1 372 IF(ISTORE .GT.LWORK)THEN277 IF(ISTORE > LWORK)THEN 373 278 WRITE(NUMSED,*)' INSUFFICIENT STORAGE FOR WORK, MIN. LWORK=',ISTORE 374 ARRET=.TRUE.375 END IF376 ! ------- ENTRY POINTS FOR INTEGER WORKSPACE -----377 IEIP=3378 ! --------- TOTAL REQUIREMENT ---------------379 ISTORE=IEIP+N-1380 IF(ISTORE.GT.LIWORK)THEN381 WRITE(NUMSED,*)' INSUFF. STORAGE FOR IWORK, MIN. LIWORK=',ISTORE382 279 ARRET=.TRUE. 383 280 END IF … … 387 284 RETURN 388 285 END IF 286 389 287 ! -------- CALL TO CORE INTEGRATOR ------------ 390 IF (ROSM == 1) THEN391 CALL RO4COR(N,X,Y,XEND,HMAX,H,RTOL,ATOL,ITOL,JAC, IJAC,&288 IF (ROSM == 4) THEN 289 CALL RO4COR(N,X,Y,XEND,HMAX,H,RTOL,ATOL,ITOL,JAC, & 392 290 MLJAC,MUJAC,IDID, & 393 NMAX,UROUND,METH,FAC1,FAC2,FACREJ,JBAND, & 394 LDJAC,LDE,WORK(IEYNEW),WORK(IEDY1),WORK(IEDY), & 395 WORK(IEAK1),WORK(IEAK2),WORK(IEAK3),WORK(IEAK4), & 396 WORK(IEFX),WORK(IEJAC),WORK(IEE),IWORK(IEIP),RSTAT) 397 ELSE IF (ROSM == 2) THEN 398 ! -------- CALL TO CORE INTEGRATOR ------------ 399 CALL RO3COR(N,X,Y,XEND,HMAX,H,RTOL,ATOL,ITOL,JAC,IJAC, & 291 NMAX,UROUND,FAC1,FAC2,FACREJ, & 292 LDJAC,LDE,RSTAT ) 293 ELSE IF (ROSM == 3) THEN 294 CALL RO3COR(N,X,Y,XEND,HMAX,H,RTOL,ATOL,ITOL,JAC, & 400 295 MLJAC,MUJAC,IDID, & 401 NMAX,UROUND,METH,FAC1,FAC2,FACREJ,JBAND, & 402 LDJAC,LDE,WORK(IEYNEW),WORK(IEDY1),WORK(IEDY), & 403 WORK(IEAK1),WORK(IEAK2),WORK(IEAK3),WORK(IEAK4), & 404 WORK(IEFX),WORK(IEJAC),WORK(IEE),IWORK(IEIP),RSTAT) 405 ELSE IF (ROSM == 3) THEN 406 CALL RO2COR(N,X,Y,XEND,HMAX,H,RTOL,ATOL,ITOL,JAC,IJAC, & 407 MLJAC,MUJAC,IDID, & 408 NMAX,UROUND,METH,FAC1,FAC2,FACREJ,JBAND, & 409 LDJAC,LDE,WORK(IEYNEW),WORK(IEDY1),WORK(IEDY), & 410 WORK(IEAK1),WORK(IEAK2),WORK(IEAK3),WORK(IEAK4), & 411 WORK(IEFX),WORK(IEJAC),WORK(IEE),IWORK(IEIP),RSTAT) 296 NMAX,UROUND,FAC1,FAC2,FACREJ, & 297 LDJAC,LDE,RSTAT ) 412 298 ENDIF 413 299 ! ----------- RETURN ----------- 414 300 415 ISTAT(1) = NFCN 416 ISTAT(2) = NJAC 417 ISTAT(3) = NSTEP 301 ISTAT(:,1) = NFCN(:) 302 ISTAT(:,2) = NJAC(:) 303 ISTAT(:,3) = NSTEP(:) 304 305 DEALLOCATE (NFCN, NJAC, NSTEP, NACCPT, NREJCT, NDEC, NSOL ) 306 307 IF ( ln_timing ) CALL timing_stop('rosk') 418 308 419 309 RETURN 420 310 421 END SUBROUTINE rosk 422 423 SUBROUTINE RO4COR(N,X,Y,XEND,HMAX,H,RTOL,ATOL,ITOL,JAC, & 424 IJAC,MLJAC,MUJAC,IDID, & 425 NMAX,UROUND,METH,FAC1,FAC2,FACREJ,BANDED, & 426 LFJAC,LE,YNEW,DY1,DY,AK1,AK2,AK3,AK4,FX,FJAC,E,IP,RSTAT) 311 END SUBROUTINE ROSK 312 313 SUBROUTINE RO3COR(N,X,Y,XEND,HMAX,H,RTOL,ATOL,ITOL,JAC, & 314 MLJAC,MUJAC,IDID, NMAX,UROUND,FAC1,FAC2,FACREJ, & 315 LFJAC,LE,RSTAT) 427 316 ! ---------------------------------------------------------- 428 317 ! CORE INTEGRATOR FOR ROS4 … … 432 321 ! DECLARATIONS 433 322 ! ---------------------------------------------------------- 434 IMPLICIT REAL*8 (A-H,O-Z) 435 REAL*8 Y(N),YNEW(N),DY1(N),DY(N),AK1(N), & 436 AK2(N),AK3(N),AK4(N),FX(N), & 437 FJAC(LFJAC,N),E(LE,N),ATOL(1),RTOL(1) 438 INTEGER IP(N) 439 LOGICAL REJECT,RJECT2,BANDED 440 REAL(wp), DIMENSION(2) :: RSTAT 441 COMMON/STAT/NFCN,NJAC,NSTEP,NACCPT,NREJCT,NDEC,NSOL 323 IMPLICIT REAL(wp) (A-H,O-Z) 324 IMPLICIT INTEGER (I-N) 325 326 REAL(wp) :: ATOL(1),RTOL(1) 327 REAL(wp), DIMENSION(jpoce,N) :: Y, YNEW, DY1, DY, AK1, AK2, AK3, AK4, FX 328 REAL(wp), DIMENSION(jpoce,LFJAC,N) :: FJAC 329 REAL(wp), DIMENSION(jpoce, LE, N) :: E 330 REAL(wp), DIMENSION(jpoce) :: H, HNEW, HMAXN, XI, FAC 331 REAL(wp), DIMENSION(jpoce) :: HC21, HC31, HC32, HC41, HC42, HC43 332 REAL(wp), DIMENSION(jpoce) :: XXOLD, HOPT, ERR 333 INTEGER, DIMENSION(jpoce,N) :: IP 334 LOGICAL, DIMENSION(jpoce) :: REJECT,RJECT2 335 INTEGER, DIMENSION(jpoce) :: ACCMASK, ENDMASK, ERRMASK 336 REAL(wp), DIMENSION(jpoce,2) :: RSTAT 337 338 IF ( ln_timing ) CALL timing_start('ro3cor') 442 339 443 340 ! ---- PREPARE BANDWIDTHS ----- 444 IF (BANDED) THEN 445 MLE=MLJAC 446 MUE=MUJAC 447 MBJAC=MLJAC+MUJAC+1 448 MDIAG=MLE+MUE+1 449 END IF 341 MLE=MLJAC 342 MUE=MUJAC 343 MBJAC=MLJAC+MUJAC+1 344 MDIAG=MLE+MUE+1 450 345 ! *** *** *** *** *** *** *** 451 346 ! INITIALISATIONS 452 347 ! *** *** *** *** *** *** *** 453 POSNEG= DSIGN(1.D0,XEND-X)454 IF (METH.EQ.1) CALL SHAMP (A21,A31,A32,C21,C31,C32,C41,C42,C43, &348 POSNEG=SIGN(1.0,XEND-X) 349 CALL RODAS3 (A21,A31,A32,A41,A42,A43,C21,C31,C32,C41,C42,C43, & 455 350 B1,B2,B3,B4,E1,E2,E3,E4,GAMMA) 456 IF (METH.EQ.2) CALL GRK4T (A21,A31,A32,C21,C31,C32,C41,C42,C43, &457 B1,B2,B3,B4,E1,E2,E3,E4,GAMMA)458 IF (METH.EQ.3) CALL GRK4A (A21,A31,A32,C21,C31,C32,C41,C42,C43, &459 B1,B2,B3,B4,E1,E2,E3,E4,GAMMA)460 IF (METH.EQ.4) CALL VELDS (A21,A31,A32,C21,C31,C32,C41,C42,C43, &461 B1,B2,B3,B4,E1,E2,E3,E4,GAMMA)462 IF (METH.EQ.5) CALL VELDD (A21,A31,A32,C21,C31,C32,C41,C42,C43, &463 B1,B2,B3,B4,E1,E2,E3,E4,GAMMA)464 IF (METH.EQ.6) CALL LSTAB (A21,A31,A32,C21,C31,C32,C41,C42,C43, &465 B1,B2,B3,B4,E1,E2,E3,E4,GAMMA)466 351 467 352 ! --- INITIAL PREPARATIONS 468 HMAXN=DMIN1(DABS(HMAX),DABS(XEND-X)) 469 H=DMIN1(DMAX1(1.D-10,DABS(H)),HMAXN) 470 H=DSIGN(H,POSNEG) 471 REJECT=.FALSE. 472 NSING=0 473 IRTRN=1 474 XXOLD=X 475 476 IF (IRTRN.LT.0) GOTO 79 353 DO ji = 1, jpoce 354 HMAXN(ji)=MIN(ABS(HMAX),ABS(XEND-X)) 355 H(ji)=MIN(MAX(1.E-10,ABS(H(ji))),HMAXN(ji)) 356 H(ji)=SIGN(H(ji),POSNEG) 357 REJECT(ji)=.FALSE. 358 XXOLD(ji)=X 359 XI(ji) = X 360 END DO 361 IRTRN = 1 362 ERRMASK(:) = 0 363 ENDMASK(:) = 0 364 ACCMASK(:) = ENDMASK(:) 365 366 IF (IRTRN < 0) GOTO 79 477 367 ! --- BASIC INTEGRATION STEP 478 1 IF (NSTEP.GT.NMAX.OR.X+.1D0*H.EQ.X.OR.ABS(H).LE.UROUND) GOTO 79 479 IF ((X-XEND)*POSNEG+UROUND.GT.0.D0) THEN 480 H=HOPT 481 IDID=1 482 RETURN 483 END IF 484 HOPT=H 485 IF ((X+H-XEND)*POSNEG.GT.0.D0) H=XEND-X 486 487 CALL sed_func(N,JINDIC,X,Y,DY1) 488 489 NFCN=NFCN+1 368 1 CONTINUE 369 DO JI = 1, jpoce 370 IF (NSTEP(ji) > NMAX .OR. XI(ji)+0.1*H(ji) == XI(ji) .OR. ABS(H(ji)) <= UROUND) ERRMASK(ji) = 1 371 IF ((XI(ji)-XEND)*POSNEG+UROUND > 0.0) THEN 372 H(ji)=HOPT(ji) 373 ENDMASK(JI) = 1 374 END IF 375 IF ( ENDMASK(ji) == 0 ) THEN 376 HOPT(JI)=H(JI) 377 IF ((XI(ji)+H(ji)-XEND)*POSNEG > 0.0) H(ji)=XEND-XI(ji) 378 ENDIF 379 END DO 380 381 ACCMASK(:) = ENDMASK(:) 382 383 IF ( COUNT( ENDMASK(:) == 1 ) == jpoce ) THEN 384 IF ( ln_timing ) CALL timing_stop('ro3cor') 385 RETURN 386 ENDIF 387 IF ( COUNT( ERRMASK(:) == 1 ) > 0 ) GOTO 79 388 389 CALL sed_func(N,Y,DY1,ACCMASK) 390 391 NFCN(:)=NFCN(:)+1 392 490 393 ! *** *** *** *** *** *** *** 491 394 ! COMPUTATION OF THE JACOBIAN 492 395 ! *** *** *** *** *** *** *** 493 NJAC=NJAC+1 494 IF (IJAC.EQ.0) THEN 495 ! --- COMPUTE JACOBIAN MATRIX NUMERICALLY 496 IF (BANDED) THEN 497 ! --- JACOBIAN IS BANDED 498 MUJACP=MUJAC+1 499 MD=MIN(MBJAC,N) 500 DO 16 K=1,MD 501 J=K 502 12 AK2(J)=Y(J) 503 AK3(J)=DSQRT(UROUND*MAX(1.D-5,ABS(Y(J)))) 504 Y(J)=Y(J)+AK3(J) 505 J=J+MD 506 IF (J.LE.N) GOTO 12 507 CALL sed_func(N,JINDIC,X,Y,AK1) 508 J=K 509 LBEG=MAX(1,J-MUJAC) 510 14 LEND=MIN(N,J+MLJAC) 511 Y(J)=AK2(J) 512 MUJACJ=MUJACP-J 513 DO 15 L=LBEG,LEND 514 15 FJAC(L+MUJACJ,J)=(AK1(L)-DY1(L))/AK3(J) 515 J=J+MD 516 LBEG=LEND+1 517 IF (J.LE.N) GOTO 14 518 16 CONTINUE 519 ELSE 520 ! --- JACOBIAN IS FULL 521 DO 18 I=1,N 522 YSAFE=Y(I) 523 DELT=DSQRT(UROUND*MAX(1.D-5,ABS(YSAFE))) 524 Y(I)=YSAFE+DELT 525 CALL sed_func(N,JINDIC,X,Y,AK1) 526 DO 17 J=1,N 527 17 FJAC(J,I)=(AK1(J)-DY1(J))/DELT 528 18 Y(I)=YSAFE 529 MLJAC=N 530 END IF 531 ELSE 532 ! --- COMPUTE JACOBIAN MATRIX ANALYTICALLY 533 CALL JAC(N,JINDIC,X,Y,FJAC,LFJAC) 534 END IF 396 NJAC(:)=NJAC(:)+1 397 CALL JAC(N,Y,FJAC,LFJAC,ACCMASK) 535 398 2 CONTINUE 399 536 400 ! *** *** *** *** *** *** *** 537 401 ! COMPUTE THE STAGES 538 402 ! *** *** *** *** *** *** *** 539 NDEC=NDEC+1 540 HC21=C21/H 541 HC31=C31/H 542 HC32=C32/H 543 HC41=C41/H 544 HC42=C42/H 545 HC43=C43/H 546 FAC=1.D0/(H*GAMMA) 547 IF (BANDED) THEN 403 DO ji = 1, jpoce 404 IF (ACCMASK(ji) == 0) THEN 405 NDEC(ji)=NDEC(ji)+1 406 HC21(ji)=C21/H(ji) 407 HC31(ji)=C31/H(ji) 408 HC32(ji)=C32/H(ji) 409 HC41(ji)=C41/H(ji) 410 HC42(ji)=C42/H(ji) 411 HC43(ji)=C43/H(ji) 412 FAC(ji)=1.0/(H(ji)*GAMMA) 548 413 ! --- THE MATRIX E (B=IDENTITY, JACOBIAN A BANDED MATRIX) 549 DO 601 J=1,N 550 I1=MAX0(1,MUJAC+2-J) 551 I2=MIN0(MBJAC,N+MUJAC+1-J) 552 DO 600 I=I1,I2 553 600 E(I+MLE,J)=-FJAC(I,J) 554 601 E(MDIAG,J)=E(MDIAG,J)+FAC 555 CALL DECB(N,LE,E,MLE,MUE,IP,INFO) 556 IF (INFO.NE.0) GOTO 80 414 DO 601 J=1,N 415 I1=MAX(1,MUJAC+2-J) 416 I2=MIN(MBJAC,N+MUJAC+1-J) 417 DO 600 I=I1,I2 418 600 E(ji,I+MLE,J)=-FJAC(ji,I,J) 419 601 E(ji,MDIAG,J)=E(ji,MDIAG,J)+FAC(ji) 420 ENDIF 421 END DO 422 CALL DECB(N,LE,E,MLE,MUE,IP,INFO,ACCMASK) 423 557 424 ! --- THIS PART COMPUTES THE STAGES IN THE CASE WHERE 558 425 ! --- 1) THE DIFFERENTIAL EQUATION IS IN EXPLICIT FORM 559 426 ! --- 2) THE JACOBIAN OF THE PROBLEM IS A BANDED MATRIX 560 427 ! --- 3) THE DIFFERENTIAL EQUATION IS AUTONOMOUS 561 DO 603 I=1,N 562 603 AK1(I)=DY1(I) 563 CALL SOLB(N,LE,E,MLE,MUE,AK1,IP) 564 DO 610 I=1,N 565 610 YNEW(I)=Y(I)+A21*AK1(I) 566 CALL sed_func(N,JINDIC,X,YNEW,DY) 567 DO 611 I=1,N 568 611 AK2(I)=DY(I)+HC21*AK1(I) 569 CALL SOLB(N,LE,E,MLE,MUE,AK2,IP) 570 DO 620 I=1,N 571 620 YNEW(I)=Y(I)+A31*AK1(I)+A32*AK2(I) 572 CALL sed_func(N,JINDIC,X,YNEW,DY) 573 DO 621 I=1,N 574 621 AK3(I)=DY(I)+HC31*AK1(I)+HC32*AK2(I) 575 CALL SOLB(N,LE,E,MLE,MUE,AK3,IP) 576 DO 631 I=1,N 577 631 AK4(I)=DY(I)+HC41*AK1(I)+HC42*AK2(I)+HC43*AK3(I) 578 CALL SOLB(N,LE,E,MLE,MUE,AK4,IP) 579 ELSE 580 ! --- THE MATRIX E (B=IDENTITY, JACOBIAN A FULL MATRIX) 581 DO 801 J=1,N 582 DO 800 I=1,N 583 800 E(I,J)=-FJAC(I,J) 584 801 E(J,J)=E(J,J)+FAC 585 CALL DEC(N,LE,E,IP,INFO) 586 IF (INFO.NE.0) GOTO 80 587 ! --- THIS PART COMPUTES THE STAGES IN THE CASE WHERE 588 ! --- 1) THE DIFFERENTIAL EQUATION IS IN EXPLICIT FORM 589 ! --- 2) THE JACOBIAN OF THE PROBLEM IS A FULL MATRIX 590 ! --- 3) THE DIFFERENTIAL EQUATION IS AUTONOMOUS 591 DO 803 I=1,N 592 803 AK1(I)=DY1(I) 593 CALL SOL(N,LE,E,AK1,IP) 594 DO 810 I=1,N 595 810 YNEW(I)=Y(I)+A21*AK1(I) 596 CALL sed_func(N,JINDIC,X,YNEW,DY) 597 DO 811 I=1,N 598 811 AK2(I)=DY(I)+HC21*AK1(I) 599 CALL SOL(N,LE,E,AK2,IP) 600 DO 820 I=1,N 601 820 YNEW(I)=Y(I)+A31*AK1(I)+A32*AK2(I) 602 CALL sed_func(N,JINDIC,X,YNEW,DY) 603 DO 821 I=1,N 604 821 AK3(I)=DY(I)+HC31*AK1(I)+HC32*AK2(I) 605 CALL SOL(N,LE,E,AK3,IP) 606 DO 831 I=1,N 607 831 AK4(I)=DY(I)+HC41*AK1(I)+HC42*AK2(I)+HC43*AK3(I) 608 CALL SOL(N,LE,E,AK4,IP) 609 END IF 610 NSOL=NSOL+4 611 NFCN=NFCN+2 428 DO ji = 1, jpoce 429 IF (ACCMASK(ji) == 0) THEN 430 AK1(ji,1:N)=DY1(ji,1:N) 431 ENDIF 432 END DO 433 CALL SOLB(N,LE,E,MLE,MUE,AK1,IP,ACCMASK) 434 DO ji = 1, jpoce 435 IF (ACCMASK(ji) == 0) THEN 436 AK2(ji,1:N)=DY1(ji,1:N)+HC21(ji)*AK1(ji,1:N) 437 ENDIF 438 END DO 439 CALL SOLB(N,LE,E,MLE,MUE,AK2,IP,ACCMASK) 440 DO ji = 1, jpoce 441 IF (ACCMASK(ji) == 0) THEN 442 YNEW(ji,1:N)=Y(ji,1:N)+A31*AK1(ji,1:N)+A32*AK2(ji,1:N) 443 ENDIF 444 END DO 445 CALL sed_func(N,YNEW,DY,ACCMASK) 446 DO ji = 1, jpoce 447 IF (ACCMASK(ji) == 0) THEN 448 AK3(ji,1:N)=DY(ji,1:N)+HC31(ji)*AK1(ji,1:N)+HC32(ji)*AK2(ji,1:N) 449 ENDIF 450 END DO 451 CALL SOLB(N,LE,E,MLE,MUE,AK3,IP,ACCMASK) 452 DO ji = 1, jpoce 453 IF (ACCMASK(ji) == 0) THEN 454 YNEW(ji,1:N)=Y(ji,1:N)+A41*AK1(ji,1:N)+A42*AK2(ji,1:N)+A43*AK3(ji,1:N) 455 ENDIF 456 END DO 457 CALL sed_func(N,YNEW,DY,ACCMASK) 458 DO ji = 1, jpoce 459 IF (ACCMASK(ji) == 0) THEN 460 DO I = 1, N 461 AK4(ji,I)=DY(ji,I)+HC41(ji)*AK1(ji,I)+HC42(ji)*AK2(ji,I)+HC43(ji)*AK3(ji,I) 462 END DO 463 ENDIF 464 END DO 465 CALL SOLB(N,LE,E,MLE,MUE,AK4,IP,ACCMASK) 466 467 DO ji = 1, jpoce 468 IF (ACCMASK(ji) == 0) THEN 469 NSOL(ji) = NSOL(ji)+4 470 NFCN(ji) = NFCN(ji)+2 612 471 ! *** *** *** *** *** *** *** 613 472 ! ERROR ESTIMATION 614 473 ! *** *** *** *** *** *** *** 615 NSTEP=NSTEP+1474 NSTEP(ji)=NSTEP(ji)+1 616 475 ! ------------ NEW SOLUTION --------------- 617 DO 240 I=1,N 618 240 YNEW(I)=Y(I)+B1*AK1(I)+B2*AK2(I)+B3*AK3(I)+B4*AK4(I) 476 DO I = 1, N 477 YNEW(ji,I)=Y(ji,I)+B1*AK1(ji,I)+B2*AK2(ji,I)+B3*AK3(ji,I)+B4*AK4(ji,I) 478 END DO 619 479 ! ------------ COMPUTE ERROR ESTIMATION ---------------- 620 ERR=0.D0 621 DO 300 I=1,N 622 S=E1*AK1(I)+E2*AK2(I)+E3*AK3(I)+E4*AK4(I) 623 IF (ITOL.EQ.0) THEN 624 SK=ATOL(1)+RTOL(1)*DMAX1(DABS(Y(I)),DABS(YNEW(I))) 625 ELSE 626 SK=ATOL(I)+RTOL(I)*DMAX1(DABS(Y(I)),DABS(YNEW(I))) 627 END IF 628 300 ERR=ERR+(S/SK)**2 629 ERR=DSQRT(ERR/N) 480 ERR(JI) = 0.0_wp 481 DO I = 1, N 482 S = E1*AK1(ji,I)+E2*AK2(ji,I)+E3*AK3(ji,I)+E4*AK4(ji,I) 483 IF (ITOL == 0) THEN 484 SK = ATOL(1)+RTOL(1)*MAX(ABS(Y(ji,I)),ABS(YNEW(ji,I))) 485 ELSE 486 SK = ATOL(I)+RTOL(I)*MAX(ABS(Y(ji,I)),ABS(YNEW(ji,I))) 487 END IF 488 ERR(ji) = ERR(ji)+(S/SK)**2 489 END DO 490 ERR(ji) = SQRT(ERR(ji)/N) 630 491 ! --- COMPUTATION OF HNEW 631 492 ! --- WE REQUIRE .2<=HNEW/H<=6. 632 FAC=DMAX1(FAC2,DMIN1(FAC1,(ERR)**.25D0/.9D0))633 HNEW=H/FAC493 FAC(ji) = MAX(FAC2,MIN(FAC1,(ERR(ji))**.333D0/.9D0)) 494 HNEW(ji) = H(ji)/FAC(ji) 634 495 ! *** *** *** *** *** *** *** 635 496 ! IS THE ERROR SMALL ENOUGH ? 636 497 ! *** *** *** *** *** *** *** 637 RJECT2= .TRUE.638 IF (ERR.LE.1.D0) THEN498 RJECT2(ji) = .TRUE. 499 IF (ERR(ji) <= 1.0) THEN 639 500 ! --- STEP IS ACCEPTED 640 NACCPT=NACCPT+1 641 DO 44 I=1,N 642 44 Y(I)=YNEW(I) 643 XXOLD=X 644 X=X+H 645 RSTAT(2) = H 646 IF (IRTRN.LT.0) GOTO 79 647 IF (DABS(HNEW).GT.HMAXN) HNEW=POSNEG*HMAXN 648 IF (REJECT) HNEW=POSNEG*DMIN1(DABS(HNEW),DABS(H)) 649 REJECT=.FALSE. 650 RJECT2=.FALSE. 651 IF (NACCPT == 1) RSTAT(1) = (H+HNEW)/2.0 652 H=HNEW 653 GOTO 1 654 ELSE 501 NACCPT(ji) = NACCPT(ji)+1 502 Y(ji,1:N) = YNEW(ji,1:N) 503 XXOLD(ji) = XI(ji) 504 XI(ji) = XI(ji)+H(ji) 505 RSTAT(ji,2) = H(ji) 506 IF (IRTRN < 0) GOTO 79 507 IF (ABS(HNEW(ji)) > HMAXN(ji)) HNEW(ji)=POSNEG*HMAXN(ji) 508 IF (REJECT(ji)) HNEW(ji)=POSNEG*MIN(ABS(HNEW(ji)),ABS(H(ji))) 509 REJECT(ji) = .FALSE. 510 RJECT2(ji) = .FALSE. 511 IF (NACCPT(ji) == 1) RSTAT(ji,1) = (H(ji)+HNEW(ji))/2.0 512 H(ji) = HNEW(ji) 513 ACCMASK(ji) = 1 514 ELSE 655 515 ! --- STEP IS REJECTED 656 IF (RJECT2) HNEW=H*FACREJ 657 IF (REJECT) RJECT2=.TRUE. 658 REJECT=.TRUE. 659 H=HNEW 660 IF (NACCPT.GE.1) NREJCT=NREJCT+1 661 GOTO 2 662 END IF 516 IF (RJECT2(ji)) HNEW(ji) = H(ji)*FACREJ 517 IF (REJECT(ji)) RJECT2(ji) = .TRUE. 518 REJECT(ji) = .TRUE. 519 H(ji)=HNEW(ji) 520 IF (NACCPT(ji) >= 1) NREJCT(ji) = NREJCT(ji)+1 521 ACCMASK(ji) = 0 522 END IF 523 ENDIF 524 END DO 525 IF (COUNT( ACCMASK(:) == 0 ) > 0 ) GOTO 2 526 GOTO 1 663 527 ! --- EXIT 664 80 WRITE (NUMSED,*) ' MATRIX E IS SINGULAR, INFO = ',INFO 665 NSING=NSING+1 666 IF (NSING.GE.5) GOTO 79 667 H=H*0.5D0 668 GOTO 2 669 79 WRITE(NUMSED,979)X,H 670 979 FORMAT(' EXIT OF ROS4 AT X=',D16.7,' H=',D16.7) 528 79 CONTINUE 529 979 FORMAT(' EXIT OF ROS3 AT X=',D16.7,' H=',D16.7) 671 530 IDID=-1 531 532 IF ( ln_timing ) CALL timing_stop('ro3cor') 533 672 534 RETURN 673 535 674 END SUBROUTINE RO4COR 675 676 677 SUBROUTINE RO3COR(N,X,Y,XEND,HMAX,H,RTOL,ATOL,ITOL,JAC, & 678 IJAC,MLJAC,MUJAC,IDID, & 679 NMAX,UROUND,METH,FAC1,FAC2,FACREJ,BANDED, & 680 LFJAC,LE,YNEW,DY1,DY,AK1,AK2,AK3,AK4,FX,FJAC,E,IP,RSTAT) 536 END SUBROUTINE RO3COR 537 538 SUBROUTINE RO4COR(N,X,Y,XEND,HMAX,H,RTOL,ATOL,ITOL,JAC, & 539 MLJAC,MUJAC,IDID, NMAX,UROUND,FAC1,FAC2,FACREJ, & 540 LFJAC,LE,RSTAT) 681 541 ! ---------------------------------------------------------- 682 542 ! CORE INTEGRATOR FOR ROS4 683 543 ! PARAMETERS SAME AS IN ROS4 WITH WORKSPACE ADDED 684 544 ! ---------------------------------------------------------- 545 ! ---------------------------------------------------------- 685 546 ! DECLARATIONS 686 547 ! ---------------------------------------------------------- 687 IMPLICIT REAL*8 (A-H,O-Z) 688 REAL*8 Y(N),YNEW(N),DY1(N),DY(N),AK1(N), & 689 AK2(N),AK3(N),AK4(N),FX(N), & 690 FJAC(LFJAC,N),E(LE,N),ATOL(1),RTOL(1) 691 INTEGER IP(N) 692 LOGICAL REJECT,RJECT2,BANDED 693 REAL(wp), DIMENSION(2) :: RSTAT 694 COMMON/STAT/NFCN,NJAC,NSTEP,NACCPT,NREJCT,NDEC,NSOL 695 548 IMPLICIT REAL(wp) (A-H,O-Z) 549 IMPLICIT INTEGER (I-N) 550 551 REAL(wp) :: ATOL(1),RTOL(1) 552 REAL(wp), DIMENSION(jpoce,N) :: Y, YNEW, DY1, DY, AK1, AK2, AK3, AK4, FX 553 REAL(wp), DIMENSION(jpoce,N) :: AK5, AK6 554 REAL(wp), DIMENSION(jpoce,LFJAC,N) :: FJAC 555 REAL(wp), DIMENSION(jpoce, LE, N) :: E 556 REAL(wp), DIMENSION(jpoce) :: H, HNEW, HMAXN, XI, FAC 557 REAL(wp), DIMENSION(jpoce) :: HC21, HC31, HC32, HC41, HC42, HC43 558 REAL(wp), DIMENSION(jpoce) :: HC51, HC52, HC53, HC54, HC61, HC62 559 REAL(wp), DIMENSION(jpoce) :: HC63, HC64, HC65 560 REAL(wp), DIMENSION(jpoce) :: XXOLD, HOPT, ERR 561 INTEGER, DIMENSION(jpoce,N) :: IP 562 LOGICAL, DIMENSION(jpoce) :: REJECT,RJECT2 563 INTEGER, DIMENSION(jpoce) :: ACCMASK, ENDMASK, ERRMASK 564 REAL(wp), DIMENSION(jpoce,2) :: RSTAT 696 565 ! ---- PREPARE BANDWIDTHS ----- 697 IF (BANDED) THEN 698 MLE=MLJAC 699 MUE=MUJAC 700 MBJAC=MLJAC+MUJAC+1 701 MDIAG=MLE+MUE+1 702 END IF 566 MLE = MLJAC 567 MUE = MUJAC 568 MBJAC = MLJAC+MUJAC+1 569 MDIAG = MLE+MUE+1 703 570 ! *** *** *** *** *** *** *** 704 571 ! INITIALISATIONS 705 572 ! *** *** *** *** *** *** *** 706 POSNEG=DSIGN(1.D0,XEND-X) 707 CALL ROS3PR(A21,A31,A32,C21,C31,C32,B1,B2,B3,E1,E2,E3,GAMMA) 573 POSNEG = SIGN(1.0,XEND-X) 574 CALL RODAS4(A21,A31,A32,A41,A42,A43,A51,A52,A53,A54,A61,A62,A63, & 575 A64,A65,C21,C31,C32,C41,C42,C43,C51,C52,C53,C54,C61,C62,C63, & 576 C64,C65,B1,B2,B3,B4,B5,B6,E1,E2,E3,E4,E5,E6,GAMMA) 708 577 709 578 ! --- INITIAL PREPARATIONS 710 HMAXN=DMIN1(DABS(HMAX),DABS(XEND-X)) 711 H=DMIN1(DMAX1(1.D-10,DABS(H)),HMAXN) 712 H=DSIGN(H,POSNEG) 713 REJECT=.FALSE. 714 NSING=0 715 IRTRN=1 716 XXOLD=X 717 718 IF (IRTRN.LT.0) GOTO 79 579 DO ji = 1, jpoce 580 HMAXN(ji) = MIN(ABS(HMAX),ABS(XEND-X)) 581 H(ji) = MIN(MAX(1.E-10,ABS(H(ji))),HMAXN(ji)) 582 H(ji) = SIGN(H(ji),POSNEG) 583 REJECT(ji) = .FALSE. 584 XXOLD(ji) = X 585 XI(ji) = X 586 END DO 587 IRTRN = 1 588 ERRMASK(:) = 0 589 ENDMASK(:) = 0 590 ACCMASK(:) = ENDMASK(:) 591 592 IF (IRTRN < 0) GOTO 79 719 593 ! --- BASIC INTEGRATION STEP 720 1 IF (NSTEP.GT.NMAX.OR.X+.1D0*H.EQ.X.OR.ABS(H).LE.UROUND) GOTO 79 721 IF ((X-XEND)*POSNEG+UROUND.GT.0.D0) THEN 722 H=HOPT 723 IDID=1 724 RETURN 725 END IF 726 HOPT=H 727 IF ((X+H-XEND)*POSNEG.GT.0.D0) H=XEND-X 728 729 CALL sed_func(N,JINDIC,X,Y,DY1) 730 731 NFCN=NFCN+1 594 1 CONTINUE 595 DO JI = 1, jpoce 596 IF (NSTEP(ji) > NMAX .OR. XI(ji)+0.1*H(ji) == XI(ji) .OR. ABS(H(ji)) <= UROUND) THEN 597 ERRMASK(ji) = 1 598 ENDIF 599 IF ((XI(ji)-XEND)*POSNEG+UROUND > 0.0) THEN 600 H(ji) = HOPT(ji) 601 ENDMASK(JI) = 1 602 END IF 603 IF ( ENDMASK(ji) == 0 ) THEN 604 HOPT(JI) = H(JI) 605 IF ((XI(ji)+H(ji)-XEND)*POSNEG > 0.0) H(ji)=XEND-XI(ji) 606 ENDIF 607 END DO 608 609 ACCMASK(:) = ENDMASK(:) 610 611 IF ( COUNT( ENDMASK(:) == 1 ) == jpoce ) RETURN 612 IF ( COUNT( ERRMASK(:) == 1 ) > 0 ) GOTO 79 613 614 CALL sed_func(N,Y,DY1,ACCMASK) 615 616 NFCN(:) = NFCN(:) + 1 732 617 ! *** *** *** *** *** *** *** 733 618 ! COMPUTATION OF THE JACOBIAN 734 619 ! *** *** *** *** *** *** *** 735 NJAC=NJAC+1 736 IF (IJAC.EQ.0) THEN 737 ! --- COMPUTE JACOBIAN MATRIX NUMERICALLY 738 IF (BANDED) THEN 739 ! --- JACOBIAN IS BANDED 740 MUJACP=MUJAC+1 741 MD=MIN(MBJAC,N) 742 DO 16 K=1,MD 743 J=K 744 12 AK2(J)=Y(J) 745 AK3(J)=DSQRT(UROUND*MAX(1.D-5,ABS(Y(J)))) 746 Y(J)=Y(J)+AK3(J) 747 J=J+MD 748 IF (J.LE.N) GOTO 12 749 CALL sed_func(N,JINDIC,X,Y,AK1) 750 J=K 751 LBEG=MAX(1,J-MUJAC) 752 14 LEND=MIN(N,J+MLJAC) 753 Y(J)=AK2(J) 754 MUJACJ=MUJACP-J 755 DO 15 L=LBEG,LEND 756 15 FJAC(L+MUJACJ,J)=(AK1(L)-DY1(L))/AK3(J) 757 J=J+MD 758 LBEG=LEND+1 759 IF (J.LE.N) GOTO 14 760 16 CONTINUE 761 ELSE 762 ! --- JACOBIAN IS FULL 763 DO 18 I=1,N 764 YSAFE=Y(I) 765 DELT=DSQRT(UROUND*MAX(1.D-5,ABS(YSAFE))) 766 Y(I)=YSAFE+DELT 767 CALL sed_func(N,JINDIC,X,Y,AK1) 768 DO 17 J=1,N 769 17 FJAC(J,I)=(AK1(J)-DY1(J))/DELT 770 18 Y(I)=YSAFE 771 MLJAC=N 772 END IF 773 ELSE 774 ! --- COMPUTE JACOBIAN MATRIX ANALYTICALLY 775 CALL JAC(N,JINDIC,X,Y,FJAC,LFJAC) 776 END IF 620 NJAC(:) = NJAC(:)+1 621 CALL JAC(N,Y,FJAC,LFJAC,jpoce,ACCMASK) 777 622 2 CONTINUE 778 623 ! *** *** *** *** *** *** *** 779 624 ! COMPUTE THE STAGES 780 625 ! *** *** *** *** *** *** *** 781 NDEC=NDEC+1 782 HC21=C21/H 783 HC31=C31/H 784 HC32=C32/H 785 FAC=1.D0/(H*GAMMA) 786 IF (BANDED) THEN 626 DO ji = 1, jpoce 627 IF (ACCMASK(ji) == 0) THEN 628 NDEC(ji) = NDEC(ji)+1 629 HC21(ji) = C21/H(ji) 630 HC31(ji) = C31/H(ji) 631 HC32(ji) = C32/H(ji) 632 HC41(ji) = C41/H(ji) 633 HC42(ji) = C42/H(ji) 634 HC43(ji) = C43/H(ji) 635 HC51(ji) = C51/H(ji) 636 HC52(ji) = C52/H(ji) 637 HC53(ji) = C53/H(ji) 638 HC54(ji) = C54/H(ji) 639 HC61(ji) = C61/H(ji) 640 HC62(ji) = C62/H(ji) 641 HC63(ji) = C63/H(ji) 642 HC64(ji) = C64/H(ji) 643 HC65(ji) = C65/H(ji) 644 645 FAC(ji) = 1.0/(H(ji)*GAMMA) 787 646 ! --- THE MATRIX E (B=IDENTITY, JACOBIAN A BANDED MATRIX) 788 DO 601 J=1,N 789 I1=MAX0(1,MUJAC+2-J) 790 I2=MIN0(MBJAC,N+MUJAC+1-J) 791 DO 600 I=I1,I2 792 600 E(I+MLE,J)=-FJAC(I,J) 793 601 E(MDIAG,J)=E(MDIAG,J)+FAC 794 CALL DECB(N,LE,E,MLE,MUE,IP,INFO) 795 IF (INFO.NE.0) GOTO 80 647 DO 601 J=1,N 648 I1=MAX0(1,MUJAC+2-J) 649 I2=MIN0(MBJAC,N+MUJAC+1-J) 650 DO 600 I=I1,I2 651 600 E(ji,I+MLE,J)=-FJAC(ji,I,J) 652 601 E(ji,MDIAG,J)=E(ji,MDIAG,J)+FAC(ji) 653 ENDIF 654 END DO 655 CALL DECB(N,LE,E,MLE,MUE,IP,INFO,ACCMASK) 656 796 657 ! --- THIS PART COMPUTES THE STAGES IN THE CASE WHERE 797 658 ! --- 1) THE DIFFERENTIAL EQUATION IS IN EXPLICIT FORM 798 659 ! --- 2) THE JACOBIAN OF THE PROBLEM IS A BANDED MATRIX 799 660 ! --- 3) THE DIFFERENTIAL EQUATION IS AUTONOMOUS 800 DO 603 I=1,N 801 603 AK1(I)=DY1(I) 802 CALL SOLB(N,LE,E,MLE,MUE,AK1,IP) 803 DO 610 I=1,N 804 610 YNEW(I)=Y(I)+A21*AK1(I) 805 CALL sed_func(N,JINDIC,X,YNEW,DY) 806 DO 611 I=1,N 807 611 AK2(I)=DY(I)+HC21*AK1(I) 808 CALL SOLB(N,LE,E,MLE,MUE,AK2,IP) 809 DO 620 I=1,N 810 620 YNEW(I)=Y(I)+A31*AK1(I)+A32*AK2(I) 811 CALL sed_func(N,JINDIC,X,YNEW,DY) 812 DO 621 I=1,N 813 621 AK3(I)=DY(I)+HC31*AK1(I)+HC32*AK2(I) 814 CALL SOLB(N,LE,E,MLE,MUE,AK3,IP) 815 ELSE 816 ! --- THE MATRIX E (B=IDENTITY, JACOBIAN A FULL MATRIX) 817 DO 801 J=1,N 818 DO 800 I=1,N 819 800 E(I,J)=-FJAC(I,J) 820 801 E(J,J)=E(J,J)+FAC 821 CALL DEC(N,LE,E,IP,INFO) 822 IF (INFO.NE.0) GOTO 80 823 ! --- THIS PART COMPUTES THE STAGES IN THE CASE WHERE 824 ! --- 1) THE DIFFERENTIAL EQUATION IS IN EXPLICIT FORM 825 ! --- 2) THE JACOBIAN OF THE PROBLEM IS A FULL MATRIX 826 ! --- 3) THE DIFFERENTIAL EQUATION IS AUTONOMOUS 827 DO 803 I=1,N 828 803 AK1(I)=DY1(I) 829 CALL SOL(N,LE,E,AK1,IP) 830 DO 810 I=1,N 831 810 YNEW(I)=Y(I)+A21*AK1(I) 832 CALL sed_func(N,JINDIC,X,YNEW,DY) 833 DO 811 I=1,N 834 811 AK2(I)=DY(I)+HC21*AK1(I) 835 CALL SOL(N,LE,E,AK2,IP) 836 DO 820 I=1,N 837 820 YNEW(I)=Y(I)+A31*AK1(I)+A32*AK2(I) 838 CALL sed_func(N,JINDIC,X,YNEW,DY) 839 DO 821 I=1,N 840 821 AK3(I)=DY(I)+HC31*AK1(I)+HC32*AK2(I) 841 CALL SOL(N,LE,E,AK3,IP) 842 END IF 843 NSOL=NSOL+3 844 NFCN=NFCN+2 661 DO ji = 1, jpoce 662 IF (ACCMASK(ji) == 0) THEN 663 AK1(ji,1:N) = DY1(ji,1:N) 664 ENDIF 665 END DO 666 CALL SOLB(N,LE,E,MLE,MUE,AK1,IP,ACCMASK) 667 DO ji = 1, jpoce 668 IF (ACCMASK(ji) == 0) THEN 669 YNEW(ji,1:N)=Y(ji,1:N)+A21*AK1(ji,1:N) 670 ENDIF 671 END DO 672 CALL sed_func(N,YNEW,DY,ACCMASK) 673 DO ji = 1, jpoce 674 IF (ACCMASK(ji) == 0) THEN 675 AK2(ji,1:N)=DY(ji,1:N)+HC21(ji)*AK1(ji,1:N) 676 ENDIF 677 END DO 678 CALL SOLB(N,LE,E,MLE,MUE,AK2,IP,ACCMASK) 679 DO ji = 1, jpoce 680 IF (ACCMASK(ji) == 0) THEN 681 YNEW(ji,1:N)=Y(ji,1:N)+A31*AK1(ji,1:N)+A32*AK2(ji,1:N) 682 ENDIF 683 END DO 684 CALL sed_func(N,YNEW,DY,ACCMASK) 685 DO ji = 1, jpoce 686 IF (ACCMASK(ji) == 0) THEN 687 AK3(ji,1:N)=DY(ji,1:N)+HC31(ji)*AK1(ji,1:N)+HC32(ji)*AK2(ji,1:N) 688 ENDIF 689 END DO 690 CALL SOLB(N,LE,E,MLE,MUE,AK3,IP,ACCMASK) 691 DO ji = 1, jpoce 692 IF (ACCMASK(ji) == 0) THEN 693 DO I = 1, N 694 YNEW(ji,I)=Y(ji,I)+A41*AK1(ji,I)+A42*AK2(ji,I)+A43*AK3(ji,I) 695 END DO 696 ENDIF 697 END DO 698 CALL sed_func(N,YNEW,DY,ACCMASK) 699 DO ji = 1, jpoce 700 IF (ACCMASK(ji) == 0) THEN 701 DO I = 1, N 702 AK4(ji,I)=DY(ji,I)+HC41(ji)*AK1(ji,I)+HC42(ji)*AK2(ji,I)+HC43(ji)*AK3(ji,I) 703 END DO 704 ENDIF 705 END DO 706 CALL SOLB(N,LE,E,MLE,MUE,AK4,IP,ACCMASK) 707 DO ji = 1, jpoce 708 IF (ACCMASK(ji) == 0) THEN 709 DO I = 1, N 710 YNEW(ji,I)=Y(ji,I)+A51*AK1(ji,I)+A52*AK2(ji,I)+A53*AK3(ji,I)+A54*AK4(ji,I) 711 END DO 712 ENDIF 713 END DO 714 CALL sed_func(N,YNEW,DY,ACCMASK) 715 DO ji = 1, jpoce 716 IF (ACCMASK(ji) == 0) THEN 717 DO I = 1, N 718 AK5(ji,I)=DY(ji,I)+HC51(ji)*AK1(ji,I)+HC52(ji)*AK2(ji,I)+HC53(ji)*AK3(ji,I)+HC54(ji)*AK4(ji,I) 719 END DO 720 ENDIF 721 END DO 722 CALL SOLB(N,LE,E,MLE,MUE,AK5,IP,ACCMASK) 723 DO ji = 1, jpoce 724 IF (ACCMASK(ji) == 0) THEN 725 DO I = 1, N 726 YNEW(ji,I)=Y(ji,I)+A61*AK1(ji,I)+A62*AK2(ji,I)+A63*AK3(ji,I)+A64*AK4(ji,I)+A65*AK5(ji,I) 727 END DO 728 ENDIF 729 END DO 730 CALL sed_func(N,YNEW,DY,ACCMASK) 731 DO ji = 1, jpoce 732 IF (ACCMASK(ji) == 0) THEN 733 DO I = 1, N 734 AK6(ji,I)=DY(ji,I)+HC61(ji)*AK1(ji,I)+HC62(ji)*AK2(ji,I)+HC63(ji)*AK3(ji,I)+HC64(ji)*AK4(ji,I) & 735 & + HC65(ji)*AK5(ji,I) 736 END DO 737 ENDIF 738 END DO 739 CALL SOLB(N,LE,E,MLE,MUE,AK6,IP,ACCMASK) 740 741 DO ji = 1, jpoce 742 IF (ACCMASK(ji) == 0) THEN 743 NSOL(ji) = NSOL(ji) + 6 744 NFCN(ji) = NFCN(ji) + 5 845 745 ! *** *** *** *** *** *** *** 846 746 ! ERROR ESTIMATION 847 747 ! *** *** *** *** *** *** *** 848 NSTEP=NSTEP+1748 NSTEP(ji) = NSTEP(ji)+1 849 749 ! ------------ NEW SOLUTION --------------- 850 DO 240 I=1,N851 240 YNEW(I)=Y(I)+B1*AK1(I)+B2*AK2(I)+B3*AK3(I)750 DO 240 I = 1, N 751 240 YNEW(ji,I)=Y(ji,I)+B1*AK1(ji,I)+B2*AK2(ji,I)+B3*AK3(ji,I)+B4*AK4(ji,I)+B5*AK5(ji,I)+B6*AK6(ji,I) 852 752 ! ------------ COMPUTE ERROR ESTIMATION ---------------- 853 ERR=0.D0854 DO 300 I=1,N855 S=E1*AK1(I)+E2*AK2(I)+E3*AK3(I)856 IF (ITOL.EQ.0) THEN857 SK=ATOL(1)+RTOL(1)*DMAX1(DABS(Y(I)),DABS(YNEW(I)))858 ELSE859 SK=ATOL(I)+RTOL(I)*DMAX1(DABS(Y(I)),DABS(YNEW(I)))860 END IF861 300 ERR=ERR+(S/SK)**2862 ERR=DSQRT(ERR/N)753 ERR(JI) = 0.0_wp 754 DO 300 I = 1, N 755 S = E1*AK1(ji,I)+E2*AK2(ji,I)+E3*AK3(ji,I)+E4*AK4(ji,I)+E5*AK5(ji,I)+E6*AK6(ji,I) 756 IF (ITOL == 0) THEN 757 SK = ATOL(1)+RTOL(1)*MAX(ABS(Y(ji,I)),ABS(YNEW(ji,I))) 758 ELSE 759 SK = ATOL(I)+RTOL(I)*MAX(ABS(Y(ji,I)),ABS(YNEW(ji,I))) 760 END IF 761 300 ERR(ji) = ERR(ji)+(S/SK)**2 762 ERR(ji) = SQRT(ERR(ji)/N) 863 763 ! --- COMPUTATION OF HNEW 864 764 ! --- WE REQUIRE .2<=HNEW/H<=6. 865 FAC=DMAX1(FAC2,DMIN1(FAC1,(ERR)**.33D0/0.9D0))866 HNEW=H/FAC765 FAC(ji) = MAX(FAC2,MIN(FAC1,(ERR(ji))**0.25/0.9)) 766 HNEW(ji) = H(ji)/FAC(ji) 867 767 ! *** *** *** *** *** *** *** 868 768 ! IS THE ERROR SMALL ENOUGH ? 869 769 ! *** *** *** *** *** *** *** 870 RJECT2= .TRUE.871 IF (ERR.LE.1.D0) THEN770 RJECT2(ji) = .TRUE. 771 IF (ERR(ji) <= 1.0) THEN 872 772 ! --- STEP IS ACCEPTED 873 NACCPT=NACCPT+1 874 DO 44 I=1,N 875 44 Y(I)=YNEW(I) 876 XXOLD=X 877 X=X+H 878 RSTAT(2) = H 879 IF (IRTRN.LT.0) GOTO 79 880 IF (DABS(HNEW).GT.HMAXN) HNEW=POSNEG*HMAXN 881 IF (REJECT) HNEW=POSNEG*DMIN1(DABS(HNEW),DABS(H)) 882 REJECT=.FALSE. 883 RJECT2=.FALSE. 884 IF (NACCPT == 1) RSTAT(1) = (H+HNEW)/2.0 885 H=HNEW 886 GOTO 1 887 ELSE 773 NACCPT(ji) = NACCPT(ji) + 1 774 Y(ji,1:N) = YNEW(ji,1:N) 775 XXOLD(ji) = XI(ji) 776 XI(ji) = XI(ji)+H(ji) 777 RSTAT(ji,2) = H(ji) 778 IF (IRTRN < 0) GOTO 79 779 IF (ABS(HNEW(ji)) > HMAXN(ji)) HNEW(ji)=POSNEG*HMAXN(ji) 780 IF (REJECT(ji)) HNEW(ji)=POSNEG*MIN(ABS(HNEW(ji)),ABS(H(ji))) 781 REJECT(ji) = .FALSE. 782 RJECT2(ji) = .FALSE. 783 IF (NACCPT(ji) == 1) RSTAT(ji,1) = (H(ji)+HNEW(ji))/2.0 784 H(ji) = HNEW(ji) 785 ACCMASK(ji) = 1 786 ELSE 888 787 ! --- STEP IS REJECTED 889 IF (RJECT2) HNEW=H*FACREJ 890 IF (REJECT) RJECT2=.TRUE. 891 REJECT=.TRUE. 892 H=HNEW 893 IF (NACCPT.GE.1) NREJCT=NREJCT+1 894 GOTO 2 895 END IF 788 IF (RJECT2(ji)) HNEW(ji)=H(ji)*FACREJ 789 IF (REJECT(ji)) RJECT2(ji)=.TRUE. 790 REJECT(ji) = .TRUE. 791 H(ji) = HNEW(ji) 792 IF (NACCPT(ji) >= 1) NREJCT(ji)=NREJCT(ji)+1 793 ACCMASK(ji) = 0 794 END IF 795 ENDIF 796 END DO 797 IF (COUNT( ACCMASK(:) == 0 ) > 0 ) GOTO 2 798 GOTO 1 896 799 ! --- EXIT 897 80 WRITE (NUMSED,*) ' MATRIX E IS SINGULAR, INFO = ',INFO 898 NSING=NSING+1 899 IF (NSING.GE.5) GOTO 79 900 H=H*0.5D0 901 GOTO 2 902 79 WRITE(NUMSED,979)X,H 800 79 CONTINUE 903 801 979 FORMAT(' EXIT OF ROS4 AT X=',D16.7,' H=',D16.7) 904 802 IDID=-1 905 803 RETURN 906 804 907 END SUBROUTINE RO3COR 908 909 ! 910 SUBROUTINE RO2COR(N,X,Y,XEND,HMAX,H,RTOL,ATOL,ITOL,JAC, & 911 IJAC,MLJAC,MUJAC,IDID, & 912 NMAX,UROUND,METH,FAC1,FAC2,FACREJ,BANDED, & 913 LFJAC,LE,YNEW,DY1,DY,AK1,AK2,AK3,AK4,FX,FJAC,E,IP,RSTAT) 914 ! ---------------------------------------------------------- 915 ! CORE INTEGRATOR FOR ROS4 916 ! PARAMETERS SAME AS IN ROS4 WITH WORKSPACE ADDED 917 ! ---------------------------------------------------------- 918 ! DECLARATIONS 919 ! ---------------------------------------------------------- 920 IMPLICIT REAL*8 (A-H,O-Z) 921 REAL*8 Y(N),YNEW(N),DY1(N),DY(N),AK1(N), & 922 AK2(N),AK3(N),AK4(N),FX(N), & 923 FJAC(LFJAC,N),E(LE,N),ATOL(1),RTOL(1) 924 INTEGER IP(N) 925 LOGICAL REJECT,RJECT2,BANDED 926 REAL(wp), DIMENSION(2) :: RSTAT 927 COMMON/STAT/NFCN,NJAC,NSTEP,NACCPT,NREJCT,NDEC,NSOL 928 929 ! ---- PREPARE BANDWIDTHS ----- 930 IF (BANDED) THEN 931 MLE=MLJAC 932 MUE=MUJAC 933 MBJAC=MLJAC+MUJAC+1 934 MDIAG=MLE+MUE+1 935 END IF 936 ! *** *** *** *** *** *** *** 937 ! INITIALISATIONS 938 ! *** *** *** *** *** *** *** 939 POSNEG=DSIGN(1.D0,XEND-X) 940 CALL RO2COEF (A21,C21,B1,B2,E1,E2,GAMMA) 941 ! --- INITIAL PREPARATIONS 942 HMAXN=DMIN1(DABS(HMAX),DABS(XEND-X)) 943 H=DMIN1(DMAX1(1.D-10,DABS(H)),HMAXN) 944 H=DSIGN(H,POSNEG) 945 REJECT=.FALSE. 946 NSING=0 947 IRTRN=1 948 XXOLD=X 949 950 IF (IRTRN.LT.0) GOTO 79 951 ! --- BASIC INTEGRATION STEP 952 1 IF (NSTEP.GT.NMAX.OR.X+.1D0*H.EQ.X.OR.ABS(H).LE.UROUND) GOTO 79 953 IF ((X-XEND)*POSNEG+UROUND.GT.0.D0) THEN 954 H=HOPT 955 IDID=1 956 RETURN 957 END IF 958 HOPT=H 959 IF ((X+H-XEND)*POSNEG.GT.0.D0) H=XEND-X 960 961 CALL sed_func(N,JINDIC,X,Y,DY1) 962 963 NFCN=NFCN+1 964 ! *** *** *** *** *** *** *** 965 ! COMPUTATION OF THE JACOBIAN 966 ! *** *** *** *** *** *** *** 967 NJAC=NJAC+1 968 IF (IJAC.EQ.0) THEN 969 ! --- COMPUTE JACOBIAN MATRIX NUMERICALLY 970 IF (BANDED) THEN 971 ! --- JACOBIAN IS BANDED 972 MUJACP=MUJAC+1 973 MD=MIN(MBJAC,N) 974 DO 16 K=1,MD 975 J=K 976 12 AK2(J)=Y(J) 977 AK3(J)=DSQRT(UROUND*MAX(1.D-5,ABS(Y(J)))) 978 Y(J)=Y(J)+AK3(J) 979 J=J+MD 980 IF (J.LE.N) GOTO 12 981 CALL sed_func(N,JINDIC,X,Y,AK1) 982 J=K 983 LBEG=MAX(1,J-MUJAC) 984 14 LEND=MIN(N,J+MLJAC) 985 Y(J)=AK2(J) 986 MUJACJ=MUJACP-J 987 DO 15 L=LBEG,LEND 988 15 FJAC(L+MUJACJ,J)=(AK1(L)-DY1(L))/AK3(J) 989 J=J+MD 990 LBEG=LEND+1 991 IF (J.LE.N) GOTO 14 992 16 CONTINUE 993 ELSE 994 ! --- JACOBIAN IS FULL 995 DO 18 I=1,N 996 YSAFE=Y(I) 997 DELT=DSQRT(UROUND*MAX(1.D-5,ABS(YSAFE))) 998 Y(I)=YSAFE+DELT 999 CALL sed_func(N,JINDIC,X,Y,AK1) 1000 DO 17 J=1,N 1001 17 FJAC(J,I)=(AK1(J)-DY1(J))/DELT 1002 18 Y(I)=YSAFE 1003 MLJAC=N 1004 END IF 1005 ELSE 1006 ! --- COMPUTE JACOBIAN MATRIX ANALYTICALLY 1007 CALL JAC(N,JINDIC,X,Y,FJAC,LFJAC) 1008 END IF 1009 2 CONTINUE 1010 ! *** *** *** *** *** *** *** 1011 ! COMPUTE THE STAGES 1012 ! *** *** *** *** *** *** *** 1013 NDEC=NDEC+1 1014 HC21=C21/H 1015 FAC=1.D0/(H*GAMMA) 1016 IF (BANDED) THEN 1017 ! --- THE MATRIX E (B=IDENTITY, JACOBIAN A BANDED MATRIX) 1018 DO 601 J=1,N 1019 I1=MAX0(1,MUJAC+2-J) 1020 I2=MIN0(MBJAC,N+MUJAC+1-J) 1021 DO 600 I=I1,I2 1022 600 E(I+MLE,J)=-FJAC(I,J) 1023 601 E(MDIAG,J)=E(MDIAG,J)+FAC 1024 CALL DECB(N,LE,E,MLE,MUE,IP,INFO) 1025 IF (INFO.NE.0) GOTO 80 1026 ! --- THIS PART COMPUTES THE STAGES IN THE CASE WHERE 1027 ! --- 1) THE DIFFERENTIAL EQUATION IS IN EXPLICIT FORM 1028 ! --- 2) THE JACOBIAN OF THE PROBLEM IS A BANDED MATRIX 1029 ! --- 3) THE DIFFERENTIAL EQUATION IS AUTONOMOUS 1030 DO 603 I=1,N 1031 603 AK1(I)=DY1(I) 1032 CALL SOLB(N,LE,E,MLE,MUE,AK1,IP) 1033 DO 610 I=1,N 1034 610 YNEW(I)=Y(I)+A21*AK1(I) 1035 CALL sed_func(N,JINDIC,X,YNEW,DY) 1036 DO 611 I=1,N 1037 611 AK2(I)=DY(I)+HC21*AK1(I) 1038 CALL SOLB(N,LE,E,MLE,MUE,AK2,IP) 1039 ELSE 1040 ! --- THE MATRIX E (B=IDENTITY, JACOBIAN A FULL MATRIX) 1041 DO 801 J=1,N 1042 DO 800 I=1,N 1043 800 E(I,J)=-FJAC(I,J) 1044 801 E(J,J)=E(J,J)+FAC 1045 CALL DEC(N,LE,E,IP,INFO) 1046 IF (INFO.NE.0) GOTO 80 1047 ! --- THIS PART COMPUTES THE STAGES IN THE CASE WHERE 1048 ! --- 1) THE DIFFERENTIAL EQUATION IS IN EXPLICIT FORM 1049 ! --- 2) THE JACOBIAN OF THE PROBLEM IS A FULL MATRIX 1050 ! --- 3) THE DIFFERENTIAL EQUATION IS AUTONOMOUS 1051 DO 803 I=1,N 1052 803 AK1(I)=DY1(I) 1053 CALL SOL(N,LE,E,AK1,IP) 1054 DO 810 I=1,N 1055 810 YNEW(I)=Y(I)+A21*AK1(I) 1056 CALL sed_func(N,JINDIC,X,YNEW,DY) 1057 DO 811 I=1,N 1058 811 AK2(I)=DY(I)+HC21*AK1(I) 1059 CALL SOL(N,LE,E,AK2,IP) 1060 END IF 1061 NSOL=NSOL+2 1062 NFCN=NFCN+1 1063 ! *** *** *** *** *** *** *** 1064 ! ERROR ESTIMATION 1065 ! *** *** *** *** *** *** *** 1066 NSTEP=NSTEP+1 1067 ! ------------ NEW SOLUTION --------------- 1068 DO 240 I=1,N 1069 240 YNEW(I)=Y(I)+B1*AK1(I)+B2*AK2(I) 1070 ! ------------ COMPUTE ERROR ESTIMATION ---------------- 1071 ERR=0.D0 1072 DO 300 I=1,N 1073 S=E1*AK1(I)+E2*AK2(I) 1074 IF (ITOL.EQ.0) THEN 1075 SK=ATOL(1)+RTOL(1)*DMAX1(DABS(Y(I)),DABS(YNEW(I))) 1076 ELSE 1077 SK=ATOL(I)+RTOL(I)*DMAX1(DABS(Y(I)),DABS(YNEW(I))) 1078 END IF 1079 300 ERR=ERR+(S/SK)**2 1080 ERR=DSQRT(ERR/N) 1081 ! --- COMPUTATION OF HNEW 1082 ! --- WE REQUIRE .2<=HNEW/H<=6. 1083 FAC=DMAX1(FAC2,DMIN1(FAC1,(ERR)**.25D0/.9D0)) 1084 HNEW=H/FAC 1085 ! *** *** *** *** *** *** *** 1086 ! IS THE ERROR SMALL ENOUGH ? 1087 ! *** *** *** *** *** *** *** 1088 RJECT2 = .TRUE. 1089 IF (ERR.LE.1.D0) THEN 1090 ! --- STEP IS ACCEPTED 1091 NACCPT=NACCPT+1 1092 DO 44 I=1,N 1093 44 Y(I)=YNEW(I) 1094 XXOLD=X 1095 X=X+H 1096 RSTAT(2) = H 1097 IF (IRTRN.LT.0) GOTO 79 1098 IF (DABS(HNEW).GT.HMAXN) HNEW=POSNEG*HMAXN 1099 IF (REJECT) HNEW=POSNEG*DMIN1(DABS(HNEW),DABS(H)) 1100 REJECT=.FALSE. 1101 RJECT2=.FALSE. 1102 IF (NACCPT == 1) RSTAT(1) = (H+HNEW)/2.0 1103 H=HNEW 1104 GOTO 1 1105 ELSE 1106 ! --- STEP IS REJECTED 1107 IF (RJECT2) HNEW=H*FACREJ 1108 IF (REJECT) RJECT2=.TRUE. 1109 REJECT=.TRUE. 1110 H=HNEW 1111 IF (NACCPT.GE.1) NREJCT=NREJCT+1 1112 GOTO 2 1113 END IF 1114 ! --- EXIT 1115 80 WRITE (NUMSED,*) ' MATRIX E IS SINGULAR, INFO = ',INFO 1116 NSING=NSING+1 1117 IF (NSING.GE.5) GOTO 79 1118 H=H*0.5D0 1119 GOTO 2 1120 79 WRITE(NUMSED,979)X,H 1121 979 FORMAT(' EXIT OF ROS2 AT X=',D16.7,' H=',D16.7) 1122 IDID=-1 805 END SUBROUTINE RO4COR 806 807 SUBROUTINE RODAS3 (A21,A31,A32,A41,A42,A43,C21,C31,C32,C41,C42,C43, & 808 B1,B2,B3,B4,E1,E2,E3,E4,GAMMA) 809 810 REAL(wp), INTENT(out) :: A21, A31, A32, A41, A42, A43, C21, C31 811 REAL(wp), INTENT(out) :: C32, C41, C42, C43, B1, B2, B3, B4, E1 812 REAL(wp), INTENT(out) :: E2, E3, E4, GAMMA 813 814 A21= 0.0 815 A31= 2.0 816 A32= 0.0 817 A41= 2.0 818 A42= 0.0 819 A43= 1.0 820 C21= 4.0 821 C31= 1.0 822 C32=-1.0 823 C41= 1.0 824 C42=-1.0 825 C43=-8.0/3.0 826 B1= 2.0 827 B2= 0.0 828 B3= 1.0 829 B4= 1.0 830 E1= 0.0 831 E2= 0.0 832 E3= 0.0 833 E4= 1.0 834 GAMMA= 0.5 1123 835 RETURN 1124 1125 END SUBROUTINE RO2COR 1126 1127 SUBROUTINE SHAMP (A21,A31,A32,C21,C31,C32,C41,C42,C43, & 1128 B1,B2,B3,B4,E1,E2,E3,E4,GAMMA) 1129 IMPLICIT REAL*8 (A-H,O-Z) 1130 A21=2.D0 1131 A31=48.D0/25.D0 1132 A32=6.D0/25.D0 1133 C21=-8.D0 1134 C31=372.D0/25.D0 1135 C32=12.D0/5.D0 1136 C41=-112.D0/125.D0 1137 C42=-54.D0/125.D0 1138 C43=-2.D0/5.D0 1139 B1=19.D0/9.D0 1140 B2=1.D0/2.D0 1141 B3=25.D0/108.D0 1142 B4=125.D0/108.D0 1143 E1=17.D0/54.D0 1144 E2=7.D0/36.D0 1145 E3=0.D0 1146 E4=125.D0/108.D0 1147 GAMMA=.5D0 836 END SUBROUTINE RODAS3 837 838 SUBROUTINE RODAS4(A21,A31,A32,A41,A42,A43,A51,A52,A53,A54,A61,A62,A63, & 839 A64,A65,C21,C31,C32,C41,C42,C43,C51,C52,C53,C54,C61,C62,C63, & 840 C64,C65,B1,B2,B3,B4,B5,B6,E1,E2,E3,E4,E5,E6,GAMMA) 841 842 REAL(wp), INTENT(out) :: A21,A31,A32,A41,A42,A43,A51,A52,A53,A54,A61 843 REAL(wp), INTENT(out) :: A62,A63,A64,A65,C21,C31,C32,C41,C42,C43,C51 844 REAL(wp), INTENT(out) :: C52,C53,C54,C61,C62,C63,C64,C65,B1,B2,B3,B4,B5 845 REAL(wp), INTENT(out) :: B6,E1,E2,E3,E4,E5,E6,GAMMA 846 847 A21 = 0.1544000000000000E+01 848 A31 = 0.9466785280815826 849 A32 = 0.2557011698983284 850 A41 = 0.3314825187068521E+01 851 A42 = 0.2896124015972201E+01 852 A43 = 0.9986419139977817 853 A51 = 0.1221224509226641E+01 854 A52 = 0.6019134481288629E+01 855 A53 = 0.1253708332932087E+02 856 A54 =-0.6878860361058950 857 A61 = A51 858 A62 = A52 859 A63 = A53 860 A64 = A54 861 A65 = 1.0 862 C21 =-0.5668800000000000E+01 863 C31 =-0.2430093356833875E+01 864 C32 =-0.2063599157091915 865 C41 =-0.1073529058151375 866 C42 =-0.9594562251023355E+01 867 C43 =-0.2047028614809616E+02 868 C51 = 0.7496443313967647E+01 869 C52 =-0.1024680431464352E+02 870 C53 =-0.3399990352819905E+02 871 C54 = 0.1170890893206160E+02 872 C61 = 0.8083246795921522E+01 873 C62 =-0.7981132988064893E+01 874 C63 =-0.3152159432874371E+02 875 C64 = 0.1631930543123136E+02 876 C65 =-0.6058818238834054E+01 877 B1 = A51 878 B2 = A52 879 B3 = A53 880 B4 = A54 881 B5 = 1.0 882 B6 = 1.0 883 E1 = 0.0 884 E2 = 0.0 885 E3 = 0.0 886 E4 = 0.0 887 E5 = 0.0 888 E6 = 1.0 889 GAMMA= 0.25 1148 890 RETURN 1149 END SUBROUTINE SHAMP 1150 ! 1151 SUBROUTINE GRK4A (A21,A31,A32,C21,C31,C32,C41,C42,C43, & 1152 B1,B2,B3,B4,E1,E2,E3,E4,GAMMA) 1153 IMPLICIT REAL*8 (A-H,O-Z) 1154 A21= 0.1108860759493671D+01 1155 A31= 0.2377085261983360D+01 1156 A32= 0.1850114988899692D+00 1157 C21=-0.4920188402397641D+01 1158 C31= 0.1055588686048583D+01 1159 C32= 0.3351817267668938D+01 1160 C41= 0.3846869007049313D+01 1161 C42= 0.3427109241268180D+01 1162 C43=-0.2162408848753263D+01 1163 B1= 0.1845683240405840D+01 1164 B2= 0.1369796894360503D+00 1165 B3= 0.7129097783291559D+00 1166 B4= 0.6329113924050632D+00 1167 E1= 0.4831870177201765D-01 1168 E2=-0.6471108651049505D+00 1169 E3= 0.2186876660500240D+00 1170 E4=-0.6329113924050632D+00 1171 GAMMA= 0.3950000000000000D+00 1172 RETURN 1173 1174 END SUBROUTINE GRK4A 1175 ! 1176 SUBROUTINE GRK4T (A21,A31,A32,C21,C31,C32,C41,C42,C43, & 1177 B1,B2,B3,B4,E1,E2,E3,E4,GAMMA) 1178 IMPLICIT REAL*8 (A-H,O-Z) 1179 A21= 0.2000000000000000D+01 1180 A31= 0.4524708207373116D+01 1181 A32= 0.4163528788597648D+01 1182 C21=-0.5071675338776316D+01 1183 C31= 0.6020152728650786D+01 1184 C32= 0.1597506846727117D+00 1185 C41=-0.1856343618686113D+01 1186 C42=-0.8505380858179826D+01 1187 C43=-0.2084075136023187D+01 1188 B1= 0.3957503746640777D+01 1189 B2= 0.4624892388363313D+01 1190 B3= 0.6174772638750108D+00 1191 B4= 0.1282612945269037D+01 1192 E1= 0.2302155402932996D+01 1193 E2= 0.3073634485392623D+01 1194 E3=-0.8732808018045032D+00 1195 E4=-0.1282612945269037D+01 1196 GAMMA= 0.2310000000000000D+00 1197 RETURN 1198 1199 END SUBROUTINE GRK4T 1200 ! 1201 SUBROUTINE VELDS (A21,A31,A32,C21,C31,C32,C41,C42,C43, & 1202 B1,B2,B3,B4,E1,E2,E3,E4,GAMMA) 1203 ! --- METHOD GIVEN BY VAN VELDHUIZEN 1204 IMPLICIT REAL*8 (A-H,O-Z) 1205 A21= 0.2000000000000000D+01 1206 A31= 0.1750000000000000D+01 1207 A32= 0.2500000000000000D+00 1208 C21=-0.8000000000000000D+01 1209 C31=-0.8000000000000000D+01 1210 C32=-0.1000000000000000D+01 1211 C41= 0.5000000000000000D+00 1212 C42=-0.5000000000000000D+00 1213 C43= 0.2000000000000000D+01 1214 B1= 0.1333333333333333D+01 1215 B2= 0.6666666666666667D+00 1216 B3=-0.1333333333333333D+01 1217 B4= 0.1333333333333333D+01 1218 E1=-0.3333333333333333D+00 1219 E2=-0.3333333333333333D+00 1220 E3=-0.0000000000000000D+00 1221 E4=-0.1333333333333333D+01 1222 GAMMA= 0.5000000000000000D+00 1223 RETURN 1224 END SUBROUTINE VELDS 1225 ! 1226 SUBROUTINE VELDD (A21,A31,A32,C21,C31,C32,C41,C42,C43, & 1227 B1,B2,B3,B4,E1,E2,E3,E4,GAMMA) 1228 ! --- METHOD GIVEN BY VAN VELDHUIZEN 1229 IMPLICIT REAL*8 (A-H,O-Z) 1230 A21= 0.2000000000000000D+01 1231 A31= 0.4812234362695436D+01 1232 A32= 0.4578146956747842D+01 1233 C21=-0.5333333333333331D+01 1234 C31= 0.6100529678848254D+01 1235 C32= 0.1804736797378427D+01 1236 C41=-0.2540515456634749D+01 1237 C42=-0.9443746328915205D+01 1238 C43=-0.1988471753215993D+01 1239 B1= 0.4289339254654537D+01 1240 B2= 0.5036098482851414D+01 1241 B3= 0.6085736420673917D+00 1242 B4= 0.1355958941201148D+01 1243 E1= 0.2175672787531755D+01 1244 E2= 0.2950911222575741D+01 1245 E3=-0.7859744544887430D+00 1246 E4=-0.1355958941201148D+01 1247 GAMMA= 0.2257081148225682D+00 1248 RETURN 1249 END SUBROUTINE VELDD 1250 ! 1251 SUBROUTINE LSTAB (A21,A31,A32,C21,C31,C32,C41,C42,C43, & 1252 B1,B2,B3,B4,E1,E2,E3,E4,GAMMA) 1253 ! --- AN L-STABLE METHOD 1254 IMPLICIT REAL*8 (A-H,O-Z) 1255 A21= 0.2000000000000000D+01 1256 A31= 0.1867943637803922D+01 1257 A32= 0.2344449711399156D+00 1258 C21=-0.7137615036412310D+01 1259 C31= 0.2580708087951457D+01 1260 C32= 0.6515950076447975D+00 1261 C41=-0.2137148994382534D+01 1262 C42=-0.3214669691237626D+00 1263 C43=-0.6949742501781779D+00 1264 B1= 0.2255570073418735D+01 1265 B2= 0.2870493262186792D+00 1266 B3= 0.4353179431840180D+00 1267 B4= 0.1093502252409163D+01 1268 E1=-0.2815431932141155D+00 1269 E2=-0.7276199124938920D-01 1270 E3=-0.1082196201495311D+00 1271 E4=-0.1093502252409163D+01 1272 GAMMA= 0.5728200000000000D+00 1273 RETURN 1274 1275 END SUBROUTINE LSTAB 1276 1277 SUBROUTINE ROS3PR(A21,A31,A32,C21,C31,C32,B1,B2,B3,E1,E2,E3,GAMMA) 1278 IMPLICIT REAL*8 (A-H,O-Z) 1279 A21=1.0 1280 A31=1.0 1281 A32=0.0 1282 C21=-0.10156171083877702091975600115545E+01 1283 C31= 0.40759956452537699824805835358067E+01 1284 C32= 0.92076794298330791242156818474003E+01 1285 B1= 0.1E+01 1286 B2= 0.61697947043828245592553615689730E+01 1287 B3= -0.42772256543218573326238373806514 1288 E1= 0.5 1289 E2=-0.29079558716805469821718236208017E+01 1290 E3= 0.22354069897811569627360909276199 1291 GAMMA=0.43586652150845899941601945119356 1292 RETURN 1293 END SUBROUTINE ROS3PR 1294 1295 SUBROUTINE RO2COEF (A21,C21,B1,B2,E1,E2,GAMMA) 1296 IMPLICIT REAL*8 (A-H,O-Z) 1297 A21=1.D0 1298 C21=-2.D0 1299 B1=3./2. 1300 B2=1./2. 1301 E1=1./2. 1302 E2=1./2. 1303 GAMMA=1.0 + 1.0/SQRT(2.) 1304 RETURN 1305 END SUBROUTINE RO2COEF 1306 ! 1307 SUBROUTINE DEC (N, NDIM, A, IP, IER) 1308 ! VERSION REAL DOUBLE PRECISION 1309 INTEGER N,NDIM,IP,IER,NM1,K,KP1,M,I,J 1310 DOUBLE PRECISION A,T 1311 DIMENSION A(NDIM,N), IP(N) 1312 !----------------------------------------------------------------------- 1313 ! MATRIX TRIANGULARIZATION BY GAUSSIAN ELIMINATION. 1314 ! INPUT.. 1315 ! N = ORDER OF MATRIX. 1316 ! NDIM = DECLARED DIMENSION OF ARRAY A . 1317 ! A = MATRIX TO BE TRIANGULARIZED. 1318 ! OUTPUT.. 1319 ! A(I,J), I.LE.J = UPPER TRIANGULAR FACTOR, U . 1320 ! A(I,J), I.GT.J = MULTIPLIERS = LOWER TRIANGULAR FACTOR, I - L. 1321 ! IP(K), K.LT.N = INDEX OF K-TH PIVOT ROW. 1322 ! IP(N) = (-1)**(NUMBER OF INTERCHANGES) OR O . 1323 ! IER = 0 IF MATRIX A IS NONSINGULAR, OR K IF FOUND TO BE 1324 ! SINGULAR AT STAGE K. 1325 ! USE SOL TO OBTAIN SOLUTION OF LINEAR SYSTEM. 1326 ! DETERM(A) = IP(N)*A(1,1)*A(2,2)*...*A(N,N). 1327 ! IF IP(N)=O, A IS SINGULAR, SOL WILL DIVIDE BY ZERO. 1328 ! 1329 ! REFERENCE.. 1330 ! C. B. MOLER, ALGORITHM 423, LINEAR EQUATION SOLVER, 1331 ! C.A.C.M. 15 (1972), P. 274. 1332 !----------------------------------------------------------------------- 1333 IER = 0 1334 IP(N) = 1 1335 IF (N .EQ. 1) GO TO 70 1336 NM1 = N - 1 1337 DO 60 K = 1,NM1 1338 KP1 = K + 1 1339 M = K 1340 DO 10 I = KP1,N 1341 IF (DABS(A(I,K)) .GT. DABS(A(M,K))) M = I 1342 10 CONTINUE 1343 IP(K) = M 1344 T = A(M,K) 1345 IF (M .EQ. K) GO TO 20 1346 IP(N) = -IP(N) 1347 A(M,K) = A(K,K) 1348 A(K,K) = T 1349 20 CONTINUE 1350 IF (T .EQ. 0.D0) GO TO 80 1351