New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 15574 for NEMO/branches/2021/dev_r14318_RK3_stage1/src/TOP/PISCES/SED/sedmat.F90 – NEMO

Ignore:
Timestamp:
2021-12-03T20:32:50+01:00 (3 years ago)
Author:
techene
Message:

#2605 #2715 trunk merged into dev_r14318_RK3_stage1

Location:
NEMO/branches/2021/dev_r14318_RK3_stage1
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2021/dev_r14318_RK3_stage1

    • Property svn:externals

      •  

        old new  
        99 
        1010# SETTE 
        11 ^/utils/CI/sette@14244        sette 
         11^/utils/CI/sette@HEAD        sette 
         12 

  • NEMO/branches/2021/dev_r14318_RK3_stage1/src/TOP/PISCES/SED/sedmat.F90

    r10222 r15574  
    1414   PRIVATE 
    1515 
    16    PUBLIC sed_mat  
    17  
    18    INTERFACE sed_mat 
    19       MODULE PROCEDURE sed_mat_dsr, sed_mat_btb 
    20    END INTERFACE 
    21  
    22    INTEGER, PARAMETER :: nmax = 30  
    23  
     16   PUBLIC sed_mat_dsr  
     17   PUBLIC sed_mat_dsrjac 
     18   PUBLIC sed_mat_dsri 
     19   PUBLIC sed_mat_btb 
     20   PUBLIC sed_mat_btbjac 
     21   PUBLIC sed_mat_btbi 
     22   PUBLIC sed_mat_coef 
    2423 
    2524   !! $Id$ 
    2625 CONTAINS 
    2726 
    28     SUBROUTINE sed_mat_dsr( nvar, ndim, nlev, preac, psms, psol, dtsed_in ) 
     27    SUBROUTINE sed_mat_coef( nksed ) 
    2928       !!--------------------------------------------------------------------- 
    30        !!                  ***  ROUTINE sed_mat_dsr  *** 
     29       !!                  ***  ROUTINE sed_mat_coef  *** 
    3130       !! 
    3231       !! ** Purpose :  solves tridiagonal system of linear equations  
    … …  
    4948       !!---------------------------------------------------------------------- 
    5049       !! * Arguments 
    51        INTEGER , INTENT(in) ::  nvar  ! number of variable 
    52        INTEGER , INTENT(in) ::  ndim  ! number of points 
    53        INTEGER , INTENT(in) ::  nlev  ! number of sediment levels 
    54  
    55        REAL(wp), DIMENSION(ndim,nlev), INTENT(in   ) :: preac  ! reaction rates 
    56        REAL(wp), DIMENSION(ndim,nlev), INTENT(in   ) :: psms  ! reaction rates 
    57        REAL(wp), DIMENSION(ndim,nlev), INTENT(inout) :: psol   ! solution ( undersaturation values ) 
    58        REAL(wp), INTENT(in) ::  dtsed_in 
    59   
     50       INTEGER, INTENT(in) :: nksed 
     51 
    6052       !---Local declarations 
    61        INTEGER  ::  ji, jk, jn 
    62        REAL(wp), DIMENSION(ndim,nlev) :: za, zb, zc, zr 
    63        REAL(wp), DIMENSION(ndim)      :: zbet 
    64        REAL(wp), DIMENSION(ndim,nmax) :: zgamm 
    65  
    66        REAL(wp) ::  aplus,aminus   
    67        REAL(wp) ::  rplus,rminus    
    68        REAL(wp) ::  dxplus,dxminus 
    69  
     53       INTEGER  ::  ji, jk 
     54       REAL(wp) ::  aplus, aminus, dxplus, dxminus 
    7055       !---------------------------------------------------------------------- 
    7156 
    72        IF( ln_timing )  CALL timing_start('sed_mat_dsr') 
     57       IF( ln_timing )  CALL timing_start('sed_mat_coef') 
    7358 
    7459       ! Computation left hand side of linear system of  
    7560       ! equations for dissolution reaction 
    7661       !--------------------------------------------- 
    77  
    78  
     62       ! first sediment level           
     63       DO ji = 1, jpoce 
     64          aplus  = ( por(1) + por(2) ) * 0.5 
     65          dxplus = ( dz3d(ji,1) + dz3d(ji,2) ) / 2. 
     66          apluss(ji,1) = ( 1.0 / ( volw3d(ji,1) ) ) * aplus / dxplus 
     67 
     68          DO jk = 2, nksed - 1 
     69             aminus  = ( por(jk-1) + por(jk) ) * 0.5 
     70             dxminus = ( dz3d(ji,jk-1) + dz3d(ji,jk) ) / 2. 
     71 
     72             aplus   = ( por(jk+1) + por(jk) ) * 0.5 
     73             dxplus  = ( dz3d(ji,jk) + dz3d(ji,jk+1) ) / 2 
     74             ! 
     75             aminuss(ji,jk) = ( 1.0 / volw3d(ji,jk) ) * aminus / dxminus 
     76             apluss (ji,jk) = ( 1.0 / volw3d(ji,jk) ) * aplus / dxplus 
     77          END DO 
     78 
     79          aminus  = ( por(nksed-1) + por(nksed) ) * 0.5 
     80          dxminus = ( dz3d(ji,nksed-1) + dz3d(ji,nksed) ) / 2. 
     81          aminuss(ji,nksed)  = ( 1.0 / volw3d(ji,nksed) ) * aminus / dxminus 
     82 
     83       END DO 
     84 
     85       IF( ln_timing )  CALL timing_stop('sed_mat_coef') 
     86 
     87    END SUBROUTINE sed_mat_coef 
     88 
     89    SUBROUTINE sed_mat_dsr( nksed, nvar, accmask ) 
     90       !!--------------------------------------------------------------------- 
     91       !!                  ***  ROUTINE sed_mat_dsr  *** 
     92       !! 
     93       !! ** Purpose :  solves tridiagonal system of linear equations  
     94       !! 
     95       !! ** Method  :  
     96       !!        1 - computes left hand side of linear system of equations 
     97       !!            for dissolution reaction 
     98       !!         For mass balance in kbot+sediment : 
     99       !!              dz3d  (:,1) = dz(1) = 0.5 cm 
     100       !!              volw3d(:,1) = dzkbot ( see sedini.F90 )  
     101       !!              dz(2)       = 0.3 cm  
     102       !!              dz3d(:,2)   = 0.3 + dzdep   ( see seddsr.F90 )      
     103       !!              volw3d(:,2) and vols3d(l,2) are thickened ( see seddsr.F90 )  
     104       !! 
     105       !!         2 - forward/backward substitution.  
     106       !! 
     107       !!   History : 
     108       !!        !  04-10 (N. Emprin, M. Gehlen ) original 
     109       !!        !  06-04 (C. Ethe)  Module Re-organization 
     110       !!---------------------------------------------------------------------- 
     111       !! * Arguments 
     112       INTEGER , INTENT(in) ::  nvar, nksed  ! number of variable 
     113       INTEGER, DIMENSION(jpoce) :: accmask 
     114       INTEGER :: ji 
     115 
     116       !---Local declarations 
     117       INTEGER  ::  jk, jn 
     118       REAL(wp), DIMENSION(nksed) :: za, zb, zc 
     119 
     120       REAL(wp) ::  rplus,rminus    
     121       !---------------------------------------------------------------------- 
     122 
     123       IF( ln_timing )  CALL timing_start('sed_mat_dsr') 
     124 
     125       ! Computation left hand side of linear system of  
     126       ! equations for dissolution reaction 
     127       !--------------------------------------------- 
    79128       jn = nvar 
    80129       ! first sediment level           
    81        DO ji = 1, ndim 
    82           aplus  = ( ( volw3d(ji,1) / ( dz3d(ji,1) ) ) + & 
    83                         ( volw3d(ji,2) / ( dz3d(ji,2) ) ) ) / 2. 
    84           dxplus = ( dz3d(ji,1) + dz3d(ji,2) ) / 2. 
    85           rplus  = ( dtsed_in / ( volw3d(ji,1) ) ) * diff(ji,1,jn) * aplus / dxplus  
     130       DO ji = 1, jpoce 
     131          IF (accmask(ji) == 0) THEN 
     132             rplus  = apluss(ji,1) * diff(ji,1,jn) * radssol(1,jn) 
     133 
     134             za(1) = 0. 
     135             zb(1) = rplus 
     136             zc(1) = -rplus 
     137  
     138             DO jk = 2, nksed - 1 
     139                rminus  = aminuss(ji,jk) * diff(ji,jk-1,jn) * radssol(jk,jn) 
     140                rplus   = apluss (ji,jk) * diff(ji,jk,jn) * radssol(jk,jn) 
     141                !      
     142                za(jk) = -rminus 
     143                zb(jk) = rminus + rplus  
     144                zc(jk) = -rplus 
     145             END DO 
     146 
     147             rminus  = aminuss(ji,nksed) * diff(ji,nksed-1,jn) * radssol(nksed,jn) 
     148             ! 
     149             za(nksed) = -rminus 
     150             zb(nksed) = rminus 
     151             zc(nksed) = 0. 
     152 
     153             ! solves tridiagonal system of linear equations  
     154             ! ----------------------------------------------- 
     155 
     156             pwcpa(ji,1,jn) = pwcpa(ji,1,jn) - ( zc(1) * pwcp(ji,2,jn) + zb(1) * pwcp(ji,1,jn) ) 
     157             DO jk = 2, nksed - 1 
     158                pwcpa(ji,jk,jn) =  pwcpa(ji,jk,jn) - ( zc(jk) * pwcp(ji,jk+1,jn) + za(jk) * pwcp(ji,jk-1,jn)    & 
     159                &                  + zb(jk) * pwcp(ji,jk,jn) ) 
     160             ENDDO 
     161             pwcpa(ji,nksed,jn) = pwcpa(ji,nksed,jn) - ( za(nksed) * pwcp(ji,nksed-1,jn)    & 
     162             &                     + zb(nksed) * pwcp(ji,nksed,jn) ) 
     163 
     164          ENDIF 
     165       END DO 
     166 
     167       IF( ln_timing )  CALL timing_stop('sed_mat_dsr') 
     168 
     169    END SUBROUTINE sed_mat_dsr 
     170 
     171    SUBROUTINE sed_mat_dsrjac( nksed, nvar, NEQ, NROWPD, jacvode, accmask ) 
     172       !!--------------------------------------------------------------------- 
     173       !!                  ***  ROUTINE sed_mat_dsrjac  *** 
     174       !! 
     175       !! ** Purpose :  solves tridiagonal system of linear equations  
     176       !! 
     177       !! ** Method  :  
     178       !!        1 - computes left hand side of linear system of equations 
     179       !!            for dissolution reaction 
     180       !!         For mass balance in kbot+sediment : 
     181       !!              dz3d  (:,1) = dz(1) = 0.5 cm 
     182       !!              volw3d(:,1) = dzkbot ( see sedini.F90 )  
     183       !!              dz(2)       = 0.3 cm  
     184       !!              dz3d(:,2)   = 0.3 + dzdep   ( see seddsr.F90 )      
     185       !!              volw3d(:,2) and vols3d(l,2) are thickened ( see 
     186       !seddsr.F90 )  
     187       !! 
     188       !!         2 - forward/backward substitution.  
     189       !! 
     190       !!   History : 
     191       !!        !  04-10 (N. Emprin, M. Gehlen ) original 
     192       !!        !  06-04 (C. Ethe)  Module Re-organization 
     193       !!---------------------------------------------------------------------- 
     194       !! * Arguments 
     195       INTEGER , INTENT(in) ::  nvar, nksed, NEQ, NROWPD  ! number of variable 
     196       REAL, DIMENSION(jpoce,NROWPD,NEQ), INTENT(inout) :: jacvode 
     197       INTEGER, DIMENSION(jpoce), INTENT(in) :: accmask 
     198 
     199       !---Local declarations 
     200       INTEGER  ::  ji,jk, jn, jnn, jni, jnj ,jnij 
     201       REAL(wp), DIMENSION(nksed) :: za, zb, zc 
     202 
     203       REAL(wp) ::  rplus,rminus 
     204       !---------------------------------------------------------------------- 
     205 
     206       IF( ln_timing )  CALL timing_start('sed_mat_dsrjac') 
     207 
     208       ! Computation left hand side of linear system of  
     209       ! equations for dissolution reaction 
     210       !--------------------------------------------- 
     211       jn = nvar 
     212       ! first sediment level           
     213       DO ji = 1, jpoce 
     214          IF (accmask(ji) == 0 ) THEN 
     215             rplus  = apluss(ji,1) * diff(ji,1,jn) * radssol(1,jn) 
     216 
     217             za(1) = 0. 
     218             zb(1) = rplus 
     219             zc(1) = -rplus 
     220 
     221             DO jk = 2, nksed - 1 
     222                rminus  = aminuss(ji,jk) * diff(ji,jk-1,jn) * radssol(jk,jn) 
     223                rplus   = apluss (ji,jk) * diff(ji,jk,jn) * radssol(jk,jn) 
     224                !      
     225                za(jk) = -rminus 
     226                zb(jk) = rminus + rplus 
     227                zc(jk) = -rplus 
     228             END DO 
     229 
     230             rminus  = aminuss(ji,nksed) * diff(ji,nksed-1,jn) * radssol(nksed,jn) 
     231             ! 
     232             za(nksed) = -rminus 
     233             zb(nksed) = rminus 
     234             zc(nksed) = 0. 
     235 
     236             ! solves tridiagonal system of linear equations  
     237 
     238             jnn = isvode(jn) 
     239             jnij = jpvode + 1 
     240             jacvode(ji, jnij, jnn) = jacvode(ji,jnij,jnn) - zb(1) 
     241             jnj = jpvode + jnn 
     242             jnij = jnn - jnj + jpvode + 1 
     243             jacvode(ji, jnij, jnj) = jacvode(ji, jnij, jnj) -zc(1) 
     244             DO jk = 2, nksed - 1 
     245                jni = (jk-1) * jpvode + jnn 
     246                jnj = (jk-2) * jpvode + jnn 
     247                jnij = jni - jnj + jpvode + 1 
     248                jacvode(ji, jnij, jnj) = jacvode(ji, jnij, jnj) - za(jk) 
     249                jnj = (jk-1) * jpvode + jnn 
     250                jnij = jni - jnj + jpvode + 1 
     251                jacvode(ji, jnij, jnj) = jacvode(ji, jnij, jnj) - zb(jk) 
     252                jnj = (jk) * jpvode + jnn 
     253                jnij = jni - jnj + jpvode + 1 
     254                jacvode(ji, jnij, jnj)   = jacvode(ji, jnij, jnj) - zc(jk) 
     255             END DO 
     256             jni = (nksed-1) * jpvode + jnn 
     257             jnj = (nksed-2) * jpvode + jnn 
     258             jnij = jni - jnj + jpvode + 1 
     259             jacvode(ji, jnij, jnj) = jacvode(ji, jnij, jnj) - za(nksed) 
     260             jnij = jpvode + 1 
     261             jacvode(ji, jnij, jni) = jacvode(ji, jnij, jni) - zb(nksed) 
     262          ENDIF 
     263       END DO 
     264 
     265       IF( ln_timing )  CALL timing_stop('sed_mat_dsrjac') 
     266 
     267    END SUBROUTINE sed_mat_dsrjac 
     268 
     269    SUBROUTINE sed_mat_btbi( nksed, nvar, psol, preac, dtsed_in ) 
     270       !!--------------------------------------------------------------------- 
     271       !!                  ***  ROUTINE sed_mat_btb  *** 
     272       !! 
     273       !! ** Purpose :  solves tridiagonal system of linear equations  
     274       !! 
     275       !! ** Method  :  
     276       !!        1 - computes left hand side of linear system of equations 
     277       !!            for dissolution reaction 
     278       !! 
     279       !! 
     280       !!         2 - forward/backward substitution.  
     281       !! 
     282       !!   History : 
     283       !!        !  04-10 (N. Emprin, M. Gehlen ) original 
     284       !!        !  06-04 (C. Ethe)  Module Re-organization 
     285       !!---------------------------------------------------------------------- 
     286       !! * Arguments 
     287       INTEGER , INTENT(in) :: nksed, nvar      ! number of sediment levels 
     288 
     289      REAL(wp), DIMENSION(jpoce,nksed,nvar), INTENT(inout) :: & 
     290          psol, preac 
     291 
     292      REAL(wp), INTENT(in) :: dtsed_in 
     293 
     294       !---Local declarations 
     295       INTEGER  ::  & 
     296          ji, jk, jn 
     297 
     298       REAL(wp) ::  & 
     299          aplus,aminus   ,  & 
     300          rplus,rminus   ,  & 
     301          dxplus,dxminus 
     302 
     303       REAL(wp), DIMENSION(nksed)    :: za, zb, zc 
     304       REAL(wp), DIMENSION(nksed)    :: zr, zgamm 
     305       REAL(wp) ::  zbet 
     306 
     307       !---------------------------------------------------------------------- 
     308 
     309      ! Computation left hand side of linear system of  
     310      ! equations for dissolution reaction 
     311      !--------------------------------------------- 
     312      IF( ln_timing )  CALL timing_start('sed_mat_btbi') 
     313 
     314      ! first sediment level           
     315      DO ji = 1, jpoce 
     316         aplus  = ( por1(2) + por1(3) ) / 2.0 
     317         dxplus = ( dz(2) + dz(3) ) / 2. 
     318         rplus  = ( dtsed_in / vols(2) ) * db(ji,2) * aplus / dxplus 
     319         za(2) = 0. 
     320         zb(2) = 1. + rplus 
     321         zc(2) = -rplus 
     322 
     323         DO jk = 3, nksed - 1 
     324            aminus  = ( por1(jk-1) + por1(jk) ) * 0.5 
     325            aminus  = ( ( vols(jk-1) / dz(jk-1) ) + ( vols(jk) / dz(jk) ) ) / 2. 
     326            dxminus = ( dz(jk-1) + dz(jk) ) / 2. 
     327            rminus  = ( dtsed_in / vols(jk) ) * db(ji,jk-1) * aminus / dxminus 
     328            ! 
     329            aplus   = ( por1(jk) + por1(jk+1) ) * 0.5 
     330            dxplus  = ( dz(jk) + dz(jk+1) ) / 2. 
     331            rplus   = ( dtsed_in / vols(jk) ) * db(ji,jk) * aplus / dxplus 
     332            !      
     333            za(jk) = -rminus 
     334            zb(jk) = 1. + rminus + rplus 
     335            zc(jk) = -rplus 
     336 
     337         ENDDO 
     338 
     339         aminus = ( por1(nksed-1) + por1(nksed) ) * 0.5 
     340         dxminus = ( dz(nksed-1) + dz(nksed) ) / 2. 
     341         rminus  = ( dtsed_in / vols(nksed) ) * db(ji,nksed-1) * aminus / dxminus 
     342         ! 
     343         za(nksed) = -rminus 
     344         zb(nksed) = 1. + rminus 
     345         zc(nksed) = 0. 
     346 
     347         ! solves tridiagonal system of linear equations  
     348         ! -----------------------------------------------     
     349         DO jn = 1, nvar 
     350            zr(:) = psol(ji,:,jn) 
     351            zbet     = zb(2) - preac(ji,2,jn) * dtsed_in 
     352            psol(ji,2,jn) = zr(2) / zbet 
     353            !  
     354            DO jk = 3, nksed 
     355               zgamm(jk) =  zc(jk-1) / zbet 
     356               zbet      =  zb(jk) - preac(ji,jk,jn) * dtsed_in - za(jk) * zgamm(jk) 
     357               psol(ji,jk,jn) = ( zr(jk) - za(jk) * psol(ji,jk-1,jn) ) / zbet 
     358            ENDDO 
     359            !  
     360            DO jk = nksed - 1, 2, -1 
     361               psol(ji,jk,jn) = psol(ji,jk,jn) - zgamm(jk+1) * psol(ji,jk+1,jn) 
     362            ENDDO 
     363         END DO 
     364      END DO 
     365      ! 
     366      IF( ln_timing )  CALL timing_stop('sed_mat_btbi') 
     367 
     368    END SUBROUTINE sed_mat_btbi 
     369 
     370 
     371    SUBROUTINE sed_mat_btb( nksed, nvar, accmask ) 
     372       !!--------------------------------------------------------------------- 
     373       !!                  ***  ROUTINE sed_mat_btb  *** 
     374       !! 
     375       !! ** Purpose :  solves tridiagonal system of linear equations  
     376       !! 
     377       !! ** Method  :  
     378       !!        1 - computes left hand side of linear system of equations 
     379       !!            for dissolution reaction 
     380       !! 
     381       !!         2 - forward/backward substitution.  
     382       !! 
     383       !!   History : 
     384       !!        !  04-10 (N. Emprin, M. Gehlen ) original 
     385       !!        !  06-04 (C. Ethe)  Module Re-organization 
     386       !!---------------------------------------------------------------------- 
     387       !! * Arguments 
     388       INTEGER , INTENT(in) :: & 
     389          nvar, nksed     ! number of sediment levels 
     390       INTEGER, DIMENSION(jpoce) :: accmask 
     391 
     392       !---Local declarations 
     393       INTEGER  ::  ji, jk, jn 
     394 
     395       REAL(wp) ::  & 
     396          aplus,aminus   ,  & 
     397          rplus,rminus   ,  & 
     398          dxplus,dxminus 
     399 
     400       REAL(wp), DIMENSION(nksed)      :: za, zb, zc 
     401 
     402       !---------------------------------------------------------------------- 
     403 
     404      ! Computation left hand side of linear system of  
     405      ! equations for dissolution reaction 
     406      !--------------------------------------------- 
     407      IF( ln_timing )  CALL timing_start('sed_mat_btb') 
     408 
     409      ! first sediment level           
     410      jn = nvar 
     411      DO ji = 1, jpoce 
     412         IF (accmask(ji) == 0) THEN 
     413         aplus  = ( por1(2) + por1(3) ) / 2.0 
     414         dxplus = ( dz(2) + dz(3) ) / 2. 
     415         rplus  = ( 1.0 / vols(2) ) * db(ji,2) * aplus / dxplus * rads1sol(2,jn) 
     416 
     417         za(2) = 0. 
     418         zb(2) = rplus  
     419         zc(2) = -rplus 
     420 
     421         DO jk = 3, nksed - 1 
     422            aminus  = ( por1(jk-1) + por1(jk) ) * 0.5 
     423            aminus  = ( ( vols(jk-1) / dz(jk-1) ) + ( vols(jk) / dz(jk) ) ) / 2. 
     424            dxminus = ( dz(jk-1) + dz(jk) ) / 2. 
     425            rminus  = ( 1.0 / vols(jk) ) * db(ji,jk-1) * aminus / dxminus * rads1sol(jk,jn) 
     426            ! 
     427            aplus   = ( por1(jk) + por1(jk+1) ) * 0.5 
     428            dxplus  = ( dz(jk) + dz(jk+1) ) / 2. 
     429            rplus   = ( 1.0 / vols(jk) ) * db(ji,jk) * aplus / dxplus * rads1sol(jk,jn) 
     430            !      
     431            za(jk) = -rminus 
     432            zb(jk) = rminus + rplus 
     433            zc(jk) = -rplus 
     434 
     435         ENDDO 
     436 
     437         aminus = ( por1(nksed-1) + por1(nksed) ) * 0.5 
     438         dxminus = ( dz(nksed-1) + dz(nksed) ) / 2. 
     439         rminus  = ( 1.0 / vols(nksed) ) * db(ji,nksed-1) * aminus / dxminus * rads1sol(nksed,jn) 
     440         ! 
     441         za(nksed) = -rminus 
     442         zb(nksed) = rminus 
     443         zc(nksed) = 0. 
     444 
     445         ! solves tridiagonal system of linear equations  
     446         ! -----------------------------------------------     
     447         pwcpa(ji,2,jn) = pwcpa(ji,2,jn) - ( zc(2) * pwcp(ji,3,jn) + zb(2) * pwcp(ji,2,jn) ) 
     448         DO jk = 3, nksed-1 
     449            pwcpa(ji,jk,jn) =  pwcpa(ji,jk,jn) - ( zc(jk) * pwcp(ji,jk+1,jn) + za(jk) * pwcp(ji,jk-1,jn)    & 
     450            &                  + zb(jk) * pwcp(ji,jk,jn) ) 
     451         ENDDO 
     452         pwcpa(ji,nksed,jn) = pwcpa(ji,nksed,jn) - ( za(nksed) * pwcp(ji,nksed-1,jn)    & 
     453         &                     + zb(nksed) * pwcp(ji,nksed,jn) ) 
     454         !  
     455         ENDIF 
     456      END DO 
     457      ! 
     458      IF( ln_timing )  CALL timing_stop('sed_mat_btb') 
     459        
     460    END SUBROUTINE sed_mat_btb 
     461 
     462    SUBROUTINE sed_mat_btbjac( nksed, nvar, NEQ, NROWPD, jacvode, accmask ) 
     463       !!--------------------------------------------------------------------- 
     464       !!                  ***  ROUTINE sed_mat_btb  *** 
     465       !! 
     466       !! ** Purpose :  solves tridiagonal system of linear equations  
     467       !! 
     468       !! ** Method  :  
     469       !!        1 - computes left hand side of linear system of equations 
     470       !!            for dissolution reaction 
     471       !! 
     472       !!         2 - forward/backward substitution.  
     473       !! 
     474       !!   History : 
     475       !!        !  04-10 (N. Emprin, M. Gehlen ) original 
     476       !!        !  06-04 (C. Ethe)  Module Re-organization 
     477       !!---------------------------------------------------------------------- 
     478       !! * Arguments 
     479       INTEGER , INTENT(in) ::  nvar, nksed, NEQ, NROWPD  ! number of variable 
     480       REAL, DIMENSION(jpoce,NROWPD,NEQ), INTENT(inout) :: jacvode 
     481       INTEGER, DIMENSION(jpoce), INTENT(in) :: accmask 
     482 
     483       !---Local declarations 
     484       INTEGER  ::  ji, jk, jn, jnn, jni, jnj ,jnij 
     485 
     486       REAL(wp) ::  & 
     487          aplus,aminus   ,  & 
     488          rplus,rminus   ,  & 
     489          dxplus,dxminus 
     490 
     491       REAL(wp), DIMENSION(nksed)      :: za, zb, zc 
     492 
     493       !---------------------------------------------------------------------- 
     494 
     495      ! Computation left hand side of linear system of  
     496      ! equations for dissolution reaction 
     497      !--------------------------------------------- 
     498      IF( ln_timing )  CALL timing_start('sed_mat_btbjac') 
     499 
     500      ! first sediment level           
     501      jn = nvar 
     502      DO ji = 1, jpoce 
     503         IF (accmask(ji) == 0) THEN 
     504         aplus  = ( por1(2) + por1(3) ) / 2.0 
     505         dxplus = ( dz(2) + dz(3) ) / 2. 
     506         rplus  = ( 1.0 / vols(2) ) * db(ji,2) * aplus / dxplus * rads1sol(2,jn) 
     507 
     508         za(2) = 0. 
     509         zb(2) = rplus 
     510         zc(2) = -rplus 
     511 
     512         DO jk = 3, nksed - 1 
     513            aminus  = ( por1(jk-1) + por1(jk) ) * 0.5 
     514            aminus  = ( ( vols(jk-1) / dz(jk-1) ) + ( vols(jk) / dz(jk) ) ) / 2. 
     515            dxminus = ( dz(jk-1) + dz(jk) ) / 2. 
     516            rminus  = ( 1.0 / vols(jk) ) * db(ji,jk-1) * aminus / dxminus * rads1sol(jk,jn) 
     517            ! 
     518            aplus   = ( por1(jk) + por1(jk+1) ) * 0.5 
     519            dxplus  = ( dz(jk) + dz(jk+1) ) / 2. 
     520            rplus   = ( 1.0 / vols(jk) ) * db(ji,jk) * aplus / dxplus * rads1sol(jk,jn) 
     521            !      
     522            za(jk) = -rminus 
     523            zb(jk) = rminus + rplus 
     524            zc(jk) = -rplus 
     525 
     526         ENDDO 
     527 
     528         aminus = ( por1(nksed-1) + por1(nksed) ) * 0.5 
     529         dxminus = ( dz(nksed-1) + dz(nksed) ) / 2. 
     530         rminus  = ( 1.0 / vols(nksed) ) * db(ji,nksed-1) * aminus / dxminus * rads1sol(nksed,jn) 
     531         ! 
     532         za(nksed) = -rminus 
     533         zb(nksed) = rminus 
     534         zc(nksed) = 0. 
     535 
     536         ! solves tridiagonal system of linear equations  
     537         ! -----------------------------------------------     
     538         jnn = isvode(jn) 
     539         jni = jpvode + jnn 
     540         jnij = jpvode + 1 
     541         jacvode(ji, jnij, jni) = jacvode(ji,jnij,jni) - zb(2) 
     542         jnj = 2 * jpvode + jnn 
     543         jnij = jni - jnj + jpvode + 1 
     544         jacvode(ji, jnij, jnj) = jacvode(ji, jnij, jnj) -zc(2) 
     545         DO jk = 3, nksed-1 
     546            jni = (jk-1) * jpvode + jnn 
     547            jnj = (jk-2) * jpvode + jnn 
     548            jnij = jni - jnj + jpvode + 1 
     549            jacvode(ji, jnij, jnj) = jacvode(ji, jnij, jnj) - za(jk) 
     550            jnj = (jk-1) * jpvode + jnn 
     551            jnij = jni - jnj + jpvode + 1 
     552            jacvode(ji, jnij, jnj) = jacvode(ji, jnij, jnj) - zb(jk) 
     553            jnj = (jk) * jpvode + jnn 
     554            jnij = jni - jnj + jpvode + 1 
     555            jacvode(ji, jnij, jnj)   = jacvode(ji, jnij, jnj) - zc(jk) 
     556         ENDDO 
     557         jni = (nksed-1) * jpvode + jnn 
     558         jnj = (nksed-2) * jpvode + jnn 
     559         jnij = jni - jnj + jpvode + 1 
     560         jacvode(ji, jnij, jnj) = jacvode(ji, jnij, jnj) - za(nksed) 
     561         jnij = jpvode + 1 
     562         jacvode(ji, jnij, jni) = jacvode(ji, jnij, jni) - zb(nksed) 
     563         !  
     564         ENDIF 
     565      END DO 
     566      ! 
     567      IF( ln_timing )  CALL timing_stop('sed_mat_btbjac') 
     568 
     569    END SUBROUTINE sed_mat_btbjac 
     570 
     571 
     572    SUBROUTINE sed_mat_dsri( nksed, nvar, preac, psms, dtsed_in, psol ) 
     573       !!--------------------------------------------------------------------- 
     574       !!                  ***  ROUTINE sed_mat_dsr  *** 
     575       !! 
     576       !! ** Purpose :  solves tridiagonal system of linear equations  
     577       !! 
     578       !! ** Method  :  
     579       !!        1 - computes left hand side of linear system of equations 
     580       !!            for dissolution reaction 
     581       !!         For mass balance in kbot+sediment : 
     582       !!              dz3d  (:,1) = dz(1) = 0.5 cm 
     583       !!              volw3d(:,1) = dzkbot ( see sedini.F90 )  
     584       !!              dz(2)       = 0.3 cm  
     585       !!              dz3d(:,2)   = 0.3 + dzdep   ( see seddsr.F90 )      
     586       !!              volw3d(:,2) and vols3d(l,2) are thickened ( see 
     587       !seddsr.F90 )  
     588       !! 
     589       !!         2 - forward/backward substitution.  
     590       !! 
     591       !!   History : 
     592       !!        !  04-10 (N. Emprin, M. Gehlen ) original 
     593       !!        !  06-04 (C. Ethe)  Module Re-organization 
     594       !!---------------------------------------------------------------------- 
     595       !! * Arguments 
     596       INTEGER , INTENT(in) ::  nksed, nvar  ! number of variable 
     597 
     598       REAL(wp), DIMENSION(jpoce,nksed), INTENT(in   ) :: preac  ! reaction rates 
     599       REAL(wp), DIMENSION(jpoce,nksed), INTENT(in   ) :: psms  ! reaction rates 
     600       REAL(wp), DIMENSION(jpoce,nksed), INTENT(inout) :: psol  ! reaction rates 
     601       REAL(wp), INTENT(in) ::  dtsed_in 
     602 
     603       !---Local declarations 
     604       INTEGER  ::  ji, jk, jn 
     605       REAL(wp), DIMENSION(jpoce,nksed) :: za, zb, zc, zr 
     606       REAL(wp), DIMENSION(jpoce)        :: zbet 
     607       REAL(wp), DIMENSION(jpoce,nksed) :: zgamm 
     608 
     609       REAL(wp) ::  rplus,rminus 
     610       !---------------------------------------------------------------------- 
     611 
     612       IF( ln_timing )  CALL timing_start('sed_mat_dsri') 
     613 
     614       ! Computation left hand side of linear system of  
     615       ! equations for dissolution reaction 
     616       !--------------------------------------------- 
     617       jn = nvar 
     618       ! first sediment level           
     619       DO ji = 1, jpoce 
     620          rplus  = dtsed_in * apluss(ji,1) * diff(ji,1,jn) * radssol(1,jn) 
    86621 
    87622          za(ji,1) = 0. 
    … …  
    89624          zc(ji,1) = -rplus 
    90625       ENDDO 
    91   
    92        DO jk = 2, nlev - 1 
    93           DO ji = 1, ndim 
    94              aminus  = ( ( volw3d(ji,jk-1) / ( dz3d(ji,jk-1) ) ) + & 
    95              &        ( volw3d(ji,jk  ) / ( dz3d(ji,jk  ) ) ) ) / 2. 
    96              dxminus = ( dz3d(ji,jk-1) + dz3d(ji,jk) ) / 2. 
    97  
    98              aplus   = ( ( volw3d(ji,jk  ) / ( dz3d(ji,jk  ) ) ) + & 
    99              &        ( volw3d(ji,jk+1) / ( dz3d(ji,jk+1) ) ) ) / 2. 
    100              dxplus  = ( dz3d(ji,jk) + dz3d(ji,jk+1) ) / 2 
    101                 ! 
    102              rminus  = ( dtsed_in / volw3d(ji,jk) ) * diff(ji,jk-1,jn) * aminus / dxminus 
    103              rplus   = ( dtsed_in / volw3d(ji,jk) ) * diff(ji,jk,jn)   * aplus / dxplus 
     626 
     627       DO jk = 2, nksed - 1 
     628          DO ji = 1, jpoce 
     629             rminus  = dtsed_in * aminuss(ji,jk) * diff(ji,jk-1,jn) * radssol(jk,jn) 
     630             rplus   = dtsed_in * apluss (ji,jk) * diff(ji,jk,jn) * radssol(jk,jn) 
    104631                !      
    105632             za(ji,jk) = -rminus 
    106              zb(ji,jk) = 1. + rminus + rplus  
     633             zb(ji,jk) = 1. + rminus + rplus 
    107634             zc(ji,jk) = -rplus 
    108635          END DO 
    109636       END DO 
    110637 
    111        DO ji = 1, ndim 
    112           aminus  = ( ( volw3d(ji,nlev-1) / dz3d(ji,nlev-1)  ) + & 
    113           &        ( volw3d(ji,nlev)  / dz3d(ji,nlev) ) ) / 2. 
    114           dxminus = ( dz3d(ji,nlev-1) + dz3d(ji,nlev) ) / 2. 
    115           rminus  = ( dtsed_in / volw3d(ji,nlev) ) * diff(ji,nlev-1,jn) * aminus / dxminus 
     638       DO ji = 1, jpoce 
     639          rminus  = dtsed_in * aminuss(ji,nksed) * diff(ji,nksed-1,jn) * radssol(nksed,jn) 
    116640          ! 
    117           za(ji,nlev) = -rminus 
    118           zb(ji,nlev) = 1. + rminus 
    119           zc(ji,nlev) = 0. 
     641          za(ji,nksed) = -rminus 
     642          zb(ji,nksed) = 1. + rminus 
     643          zc(ji,nksed) = 0. 
    120644       END DO 
    121645 
    … …  
    124648       ! ----------------------------------------------- 
    125649 
    126        zr  (:,:) = psol(:,:) + psms(:,:) 
    127        zb  (:,:) = zb(:,:) + preac(:,:) 
     650       zr  (:,:) = psol(:,:) + psms(:,:) * dtsed_in 
     651       zb  (:,:) = zb(:,:) - preac(:,:) * dtsed_in 
    128652       zbet(:  ) = zb(:,1) 
    129653       psol(:,1) = zr(:,1) / zbet(:) 
    130654 
    131655          !  
    132        DO jk = 2, nlev 
    133           DO ji = 1, ndim 
     656       DO jk = 2, nksed 
     657          DO ji = 1, jpoce 
    134658             zgamm(ji,jk) =  zc(ji,jk-1) / zbet(ji) 
    135659             zbet(ji)     =  zb(ji,jk) - za(ji,jk) * zgamm(ji,jk) 
    … …  
    138662       ENDDO 
    139663          !  
    140        DO jk = nlev - 1, 1, -1 
    141           DO ji = 1,ndim 
     664       DO jk = nksed - 1, 1, -1 
     665          DO ji = 1, jpoce 
    142666             psol(ji,jk) = psol(ji,jk) - zgamm(ji,jk+1) * psol(ji,jk+1) 
    143667          END DO 
    144668       ENDDO 
    145669 
    146        IF( ln_timing )  CALL timing_stop('sed_mat_dsr') 
    147  
    148  
    149     END SUBROUTINE sed_mat_dsr 
    150  
    151     SUBROUTINE sed_mat_btb( nvar, ndim, nlev, psol, dtsed_in ) 
    152        !!--------------------------------------------------------------------- 
    153        !!                  ***  ROUTINE sed_mat_btb  *** 
    154        !! 
    155        !! ** Purpose :  solves tridiagonal system of linear equations  
    156        !! 
    157        !! ** Method  :  
    158        !!        1 - computes left hand side of linear system of equations 
    159        !!            for dissolution reaction 
    160        !! 
    161        !!         2 - forward/backward substitution.  
    162        !! 
    163        !!   History : 
    164        !!        !  04-10 (N. Emprin, M. Gehlen ) original 
    165        !!        !  06-04 (C. Ethe)  Module Re-organization 
    166        !!---------------------------------------------------------------------- 
    167        !! * Arguments 
    168        INTEGER , INTENT(in) :: & 
    169           nvar , &  ! number of variables 
    170           ndim , &  ! number of points 
    171           nlev      ! number of sediment levels 
    172  
    173       REAL(wp), DIMENSION(ndim,nlev,nvar), INTENT(inout) :: & 
    174           psol      ! solution 
    175  
    176       REAL(wp), INTENT(in) :: dtsed_in 
    177  
    178        !---Local declarations 
    179        INTEGER  ::  & 
    180           ji, jk, jn 
    181  
    182        REAL(wp) ::  & 
    183           aplus,aminus   ,  &  
    184           rplus,rminus   ,  & 
    185           dxplus,dxminus  
    186  
    187        REAL(wp), DIMENSION(nlev)      :: za, zb, zc 
    188        REAL(wp), DIMENSION(ndim,nlev) :: zr 
    189        REAL(wp), DIMENSION(nmax)      :: zgamm  
    190        REAL(wp) ::  zbet 
    191  
    192   
    193        !---------------------------------------------------------------------- 
    194  
    195        ! Computation left hand side of linear system of  
    196        ! equations for dissolution reaction 
    197        !--------------------------------------------- 
    198  
    199  
    200       IF( ln_timing )  CALL timing_start('sed_mat_btb') 
    201  
    202        ! first sediment level           
    203       DO ji = 1, ndim 
    204          aplus  = ( ( vols(2) / dz(2) ) + ( vols(3) / dz(3) ) ) / 2. 
    205          dxplus = ( dz(2) + dz(3) ) / 2. 
    206          rplus  = ( dtsed_in / vols(2) ) * db(ji,2) * aplus / dxplus  
    207  
    208          za(1) = 0. 
    209          zb(1) = 1. + rplus 
    210          zc(1) = -rplus 
    211  
    212               
    213          DO jk = 2, nlev - 1 
    214             aminus  = ( ( vols(jk) / dz(jk) ) + ( vols(jk+1) / dz(jk+1) ) ) / 2. 
    215             dxminus = ( dz(jk) + dz(jk+1) ) / 2. 
    216             rminus  = ( dtsed_in / vols(jk+1) ) * db(ji,jk) * aminus / dxminus 
    217             ! 
    218             aplus   = ( ( vols(jk+1) / dz(jk+1  ) ) + ( vols(jk+2) / dz(jk+2) ) ) / 2. 
    219             dxplus  = ( dz(jk+1) + dz(jk+2) ) / 2. 
    220             rplus   = ( dtsed_in / vols(jk+1) ) * db(ji,jk+1) * aplus / dxplus 
    221             !      
    222             za(jk) = -rminus 
    223             zb(jk) = 1. + rminus + rplus  
    224             zc(jk) = -rplus 
    225          ENDDO 
    226   
    227          aminus  = ( ( vols(nlev) / dz(nlev) ) + ( vols(nlev+1) / dz(nlev+1) ) ) / 2. 
    228          dxminus = ( dz(nlev) + dz(nlev+1) ) / 2. 
    229          rminus  = ( dtsed_in / vols(nlev+1) ) * db(ji,nlev) * aminus / dxminus 
    230          ! 
    231          za(nlev) = -rminus 
    232          zb(nlev) = 1. + rminus 
    233          zc(nlev) = 0. 
    234  
    235  
    236          ! solves tridiagonal system of linear equations  
    237          ! -----------------------------------------------     
    238          DO jn = 1, nvar 
    239            
    240             DO jk = 1, nlev 
    241                zr  (ji,jk)    = psol(ji,jk,jn) 
    242             END DO 
    243             zbet          = zb(1) 
    244             psol(ji,1,jn) = zr(ji,1) / zbet 
    245             !  
    246             DO jk = 2, nlev 
    247                zgamm(jk) =  zc(jk-1) / zbet 
    248                zbet      =  zb(jk) - za(jk) * zgamm(jk) 
    249                psol(ji,jk,jn) = ( zr(ji,jk) - za(jk) * psol(ji,jk-1,jn) ) / zbet 
    250             ENDDO 
    251             !  
    252             DO jk = nlev - 1, 1, -1 
    253                 psol(ji,jk,jn) = psol(ji,jk,jn) - zgamm(jk+1) * psol(ji,jk+1,jn) 
    254             ENDDO 
    255  
    256          ENDDO 
    257  
    258        END DO 
    259        ! 
    260        IF( ln_timing )  CALL timing_stop('sed_mat_btb') 
    261  
    262         
    263     END SUBROUTINE sed_mat_btb 
     670       IF( ln_timing )  CALL timing_stop('sed_mat_dsri') 
     671 
     672 
     673    END SUBROUTINE sed_mat_dsri 
     674 
    264675 
    265676 END MODULE sedmat 
Note: See TracChangeset for help on using the changeset viewer.