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 15548 for NEMO/branches/2021/ticket2632_r14588_theta_sbcblk/src/TOP/PISCES/SED/sedinorg.F90 – NEMO

Ignore:
Timestamp:
2021-11-28T18:59:49+01:00 (3 years ago)
Author:
gsamson
Message:

update branch to the head of the trunk (r15547); ticket #2632

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

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2021/ticket2632_r14588_theta_sbcblk

    • Property svn:externals

      •  

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

  • NEMO/branches/2021/ticket2632_r14588_theta_sbcblk/src/TOP/PISCES/SED/sedinorg.F90

    r12839 r15548  
    88   USE sed     ! sediment global variable 
    99   USE sed_oce 
    10    USE sedmat  ! linear system of equations 
    11    USE sedco3  ! carbonate ion and proton concentration  
    1210   USE sedini 
    13    USE seddsr 
     11   USE sedmat 
    1412   USE lib_mpp         ! distribued memory computing library 
    1513   USE lib_fortran 
    … …  
    2321CONTAINS 
    2422    
    25    SUBROUTINE sed_inorg( kt )  
     23   SUBROUTINE sed_inorg( kt ) 
    2624      !!---------------------------------------------------------------------- 
    2725      !!                   ***  ROUTINE sed_inorg  *** 
    … …  
    4644      !!---------------------------------------------------------------------- 
    4745      !! Arguments 
    48       INTEGER, INTENT(in) ::   kt       ! number of iteration 
     46      INTEGER, INTENT(in)  :: kt   ! time step 
    4947      ! --- local variables 
    50       INTEGER :: ji, jk, js, jw         ! dummy looop indices 
    51       REAL(wp), DIMENSION(jpoce,jpksed) :: zrearat1, zrearat2    ! reaction rate in pore water 
    52       REAL(wp), DIMENSION(jpoce,jpksed) :: zundsat    ! undersaturation ; indice jpwatp1 is for calcite    
    53       REAL(wp), DIMENSION(jpoce) :: zco3eq 
    54       REAL(wp), DIMENSION(jpoce,jpksed,jpsol) :: zvolc    ! temp. variables 
    55       REAL(wp), DIMENSION(jpoce) :: zsieq 
    56       REAL(wp)  ::  zsolid1, zvolw, zreasat 
    57       REAL(wp)  ::  zsatur, zsatur2, znusil, zsolcpcl, zsolcpsi 
     48      INTEGER   ::  ji,jk          ! dummy looop indices 
     49      REAL(wp), DIMENSION(jpoce) ::  zsieq, reac_silf 
     50      REAL(wp)  ::  zsolid1, zreasat, zco3sat 
     51      REAL(wp)  ::  zsatur, zsatur2, znusil, zsolcpcl, zsolcpsi, zexcess 
     52      REAL(wp), DIMENSION(jpoce,jpksed) :: zundsat, zrearat, psms 
    5853      !! 
    5954      !!---------------------------------------------------------------------- 
    6055 
    6156      IF( ln_timing )  CALL timing_start('sed_inorg') 
     57 
     58      IF( kt == nitsed000 ) THEN 
     59         IF (lwp) WRITE(numsed,*) ' sed_inorg : Dissolution of CaCO3 and BSi  ' 
     60         IF (lwp) WRITE(numsed,*) ' ' 
     61      ENDIF 
    6262! 
    63       IF( kt == nitsed000 ) THEN 
    64          IF (lwp) THEN 
    65             WRITE(numsed,*) ' sed_inorg : Dissolution reaction ' 
    66             WRITE(numsed,*) ' ' 
    67          ENDIF 
    68 !         !  
    69       ENDIF 
     63      DO ji = 1, jpoce 
     64         ! ----------------------------------------------- 
     65         ! Computation of Si solubility 
     66         ! Param of Ridgwell et al. 2002 
     67         ! ----------------------------------------------- 
    7068 
    71      ! Initializations 
    72      !---------------------- 
    73        
    74       zrearat1(:,:) = 0.    ;   zundsat(:,:)  = 0.  
    75       zrearat2(:,:) = 0.    ;   zrearat2(:,:) = 0. 
    76       zco3eq(:)     = rtrn 
    77       zvolc(:,:,:)  = 0. 
    78  
    79       ! ----------------------------------------------- 
    80       ! Computation of Si solubility 
    81       ! Param of Ridgwell et al. 2002 
    82       ! ----------------------------------------------- 
    83  
    84       DO ji = 1, jpoce 
    8569         zsolcpcl = 0.0 
    8670         zsolcpsi = 0.0 
    8771         DO jk = 1, jpksed 
    88             zsolcpsi = zsolcpsi + solcp(ji,jk,jsopal) * dz(jk) 
    89             zsolcpcl = zsolcpcl + solcp(ji,jk,jsclay) * dz(jk) 
     72            zsolcpsi = zsolcpsi + solcp(ji,jk,jsopal) * vols3d(ji,jk) 
     73            zsolcpcl = zsolcpcl + solcp(ji,jk,jsclay) * vols3d(ji,jk) 
    9074         END DO 
    9175         zsolcpsi = MAX( zsolcpsi, rtrn ) 
    92          zsieq(ji) = sieqs(ji) * MAX(0.25, 1.0 - (0.045 * zsolcpcl / zsolcpsi )**0.58 ) 
    93          zsieq(ji) = MAX( rtrn, sieqs(ji) ) 
     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  
    9478      END DO 
    9579 
    96       DO js = 1, jpsol 
    97          DO jk = 1, jpksed 
    98             DO ji = 1, jpoce     
    99                zvolc(ji,jk,js) = ( vols3d(ji,jk) * dens_mol_wgt(js) ) /  & 
    100                   &              ( volw3d(ji,jk) * 1.e-3  )      
    101             ENDDO 
    102          ENDDO 
    103       ENDDO 
    104  
    105       !---------------------------------------------------------- 
    106       ! 5.  Beginning of  Pore Water diffusion and solid reaction 
    107       !--------------------------------------------------------- 
    10880       
    109       !----------------------------------------------------------------------------- 
    110       ! For jk=2,jpksed, and for couple  
    111       !  1 : jwsil/jsopal  ( SI/Opal ) 
    112       !  2 : jsclay/jsclay ( clay/clay )  
    113       !  3 : jwoxy/jspoc   ( O2/POC ) 
    114       !  reaction rate is a function of solid=concentration in solid reactif in [mol/l]  
    115       !  and undersaturation in [mol/l]. 
    116       !  Solid weight fractions should be in ie [mol/l]) 
    117       !  second member and solution are in zundsat variable 
    118       !------------------------------------------------------------------------- 
    119  
    120       DO jk = 1, jpksed 
    121          DO ji = 1, jpoce 
    122             ! For Silicic Acid and clay 
    123             zundsat(ji,jk) = zsieq(ji) - pwcp(ji,jk,jwsil) 
    124          ENDDO 
    125       ENDDO 
    126        
    127       ! Definition of reaction rates [rearat]=sans dim  
    128       ! For jk=1 no reaction (pure water without solid) for each solid compo 
    12981      DO ji = 1, jpoce 
    130          zrearat1(ji,:) = 0. 
    131          zrearat2(ji,:) = 0. 
    132       ENDDO 
    133  
    134       ! left hand side of coefficient matrix 
    135       DO jk = 2, jpksed 
    136          DO ji = 1, jpoce 
    137             zsolid1 = zvolc(ji,jk,jsopal) * solcp(ji,jk,jsopal) 
    138             zsatur = MAX(0., zundsat(ji,jk) / zsieq(ji) ) 
     82         DO jk = 2, jpksed 
     83            zsolid1 = volc(ji,jk,jsopal) * solcp(ji,jk,jsopal) 
     84            zsatur = MAX(0., ( zsieq(ji) - pwcp(ji,jk,jwsil) ) / zsieq(ji) ) 
    13985            zsatur2 = (1.0 + temp(ji) / 400.0 )**37 
    140             znusil = ( 0.225 * ( 1.0 + temp(ji) / 15.) + 0.775 * zsatur2 * zsatur**2.25 ) / zsieq(ji) 
    141             zrearat1(ji,jk)  = ( reac_sil * znusil * dtsed * zsolid1 ) / & 
    142                &                ( 1. + reac_sil * znusil * dtsed * zundsat(ji,jk) ) 
    143          ENDDO 
    144       ENDDO 
    145  
    146       CALL sed_mat( jwsil, jpoce, jpksed, zrearat1, zrearat2, zundsat, dtsed ) 
    147  
    148       ! New solid concentration values (jk=2 to jksed) for each couple  
    149       DO jk = 2, jpksed 
    150          DO ji = 1, jpoce 
    151             zreasat = zrearat1(ji,jk) * zundsat(ji,jk) / ( zvolc(ji,jk,jsopal) ) 
    152             solcp(ji,jk,jsopal) = solcp(ji,jk,jsopal) - zreasat 
    153          ENDDO 
    154       ENDDO 
    155  
    156       ! New pore water concentrations     
    157       DO jk = 1, jpksed 
    158          DO ji = 1, jpoce 
    159             ! Acid Silicic  
    160             pwcp(ji,jk,jwsil)  = zsieq(ji) - zundsat(ji,jk) 
    161          ENDDO 
    162       ENDDO 
     86            znusil = ( 0.225 * ( 1.0 + temp(ji) / 15.) * zsatur + 0.775 * zsatur2 * zsatur**9.25 ) 
     87            solcp(ji,jk,jsopal) = solcp(ji,jk,jsopal) - reac_silf(ji) * znusil * dtsed * solcp(ji,jk,jsopal) 
     88            pwcp(ji,jk,jwsil) = pwcp(ji,jk,jwsil) + reac_silf(ji) * znusil * dtsed * zsolid1 
     89         END DO 
     90      END DO 
    16391 
    16492      !--------------------------------------------------------------- 
    … …  
    16795 
    16896      ! computes co3por from the updated pwcp concentrations (note [co3por] = mol/l) 
    169  
    170       CALL sed_co3( kt ) 
    171  
    17297      ! *densSW(l)**2 converts aksps [mol2/kg sol2] into [mol2/l2] to get [undsat] in [mol/l] 
    173       DO jk = 1, jpksed 
    174          DO ji = 1, jpoce 
    175             zco3eq(ji)     = aksps(ji) * densSW(ji) * densSW(ji) / ( calcon2(ji) + rtrn ) 
    176             zco3eq(ji)     = MAX( rtrn, zco3eq(ji) )  
    177             zundsat(ji,jk) = MAX(0., zco3eq(ji) - co3por(ji,jk) ) 
    178          ENDDO 
    179       ENDDO 
    180  
    181       DO jk = 2, jpksed 
    182          DO ji = 1, jpoce 
    183             zsolid1 = zvolc(ji,jk,jscal) * solcp(ji,jk,jscal) 
    184             zrearat1(ji,jk) = ( reac_cal * dtsed * zsolid1 / zco3eq(ji) ) / & 
    185                   &               ( 1. + reac_cal * dtsed * zundsat(ji,jk) / zco3eq(ji) ) 
     98      DO ji = 1, jpoce 
     99         zco3sat = aksps(ji) * densSW(ji) * densSW(ji) / ( calcon2(ji) + rtrn ) 
     100         saturco3(ji,:) = 1.0 - co3por(ji,:) / ( rtrn + zco3sat ) 
     101         DO jk = 1, jpksed 
     102            zsolid1 = volc(ji,jk,jscal) * solcp(ji,jk,jscal) 
     103            zundsat(ji,jk) = MAX( 0., zco3sat - co3por(ji,jk) ) 
     104            zrearat(ji,jk) = ( reac_cal * zsolid1 / ( zco3sat + rtrn ) ) / & 
     105            &                ( 1. + reac_cal * dtsed * zundsat(ji,jk) / ( zco3sat + rtrn ) ) 
    186106         END DO 
    187107      END DO 
    188108 
     109      psms(:,:) = 0.0 
    189110      ! solves tridiagonal system 
    190       CALL sed_mat( jwdic, jpoce, jpksed, zrearat1, zrearat2, zundsat, dtsed ) 
     111      CALL sed_mat_dsri( jpksed, jwdic, -zrearat(:,:), psms(:,:), dtsed, zundsat ) 
    191112 
    192113      ! New solid concentration values (jk=2 to jksed) for cacO3 
    193114      DO jk = 2, jpksed 
    194115         DO ji = 1, jpoce 
    195             zreasat = zrearat1(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jscal) 
     116            zreasat = zrearat(ji,jk) * dtsed * zundsat(ji,jk) / ( volc(ji,jk,jscal) + rtrn ) 
    196117            solcp(ji,jk,jscal) = solcp(ji,jk,jscal) - zreasat 
     118            zreasat = zrearat(ji,jk) * dtsed * zundsat(ji,jk) 
     119            ! For DIC 
     120            pwcp(ji,jk,jwdic)  = pwcp(ji,jk,jwdic) + zreasat 
     121            ! For alkalinity 
     122            pwcp(ji,jk,jwalk)  = pwcp(ji,jk,jwalk) + 2.* zreasat 
    197123         ENDDO 
    198124      ENDDO 
    199125 
    200       ! New dissolved concentrations 
    201       DO jk = 1, jpksed 
    202          DO ji = 1, jpoce 
    203             zreasat = zrearat1(ji,jk) * zundsat(ji,jk)     
    204             ! For DIC 
    205             pwcp(ji,jk,jwdic)  = pwcp(ji,jk,jwdic) + zreasat 
    206             ! For alkalinity 
    207             pwcp(ji,jk,jwalk)  = pwcp(ji,jk,jwalk) + 2.0 * zreasat  
    208          ENDDO 
    209       ENDDO 
    210  
    211       !------------------------------------------------- 
    212       ! Beginning DIC, Alkalinity  
    213       !------------------------------------------------- 
    214        
    215       DO jk = 1, jpksed 
    216          DO ji = 1, jpoce       
    217             zundsat(ji,jk)   = pwcp(ji,jk,jwdic) 
    218             zrearat1(ji,jk)  = 0. 
    219          ENDDO 
    220       ENDDO 
    221  
    222       ! solves tridiagonal system 
    223       CALL sed_mat( jwdic, jpoce, jpksed, zrearat1, zrearat2, zundsat, dtsed ) 
    224  
    225      ! New dissolved concentrations       
    226       DO jk = 1, jpksed 
    227          DO ji = 1, jpoce                       
    228             pwcp(ji,jk,jwdic) = zundsat(ji,jk) 
    229          ENDDO 
    230       ENDDO             
    231  
    232       !------------------------------------------------- 
    233       ! Beginning DIC, Alkalinity  
    234       !------------------------------------------------- 
    235  
    236       DO jk = 1, jpksed 
    237          DO ji = 1, jpoce 
    238             zundsat(ji,jk) = pwcp(ji,jk,jwalk) 
    239             zrearat1(ji,jk) = 0. 
    240          ENDDO 
    241       ENDDO 
    242 ! 
    243 !      ! solves tridiagonal system 
    244       CALL sed_mat( jwalk, jpoce, jpksed, zrearat1, zrearat2, zundsat, dtsed ) 
    245 ! 
    246 !      ! New dissolved concentrations       
    247       DO jk = 1, jpksed 
    248          DO ji = 1, jpoce 
    249             pwcp(ji,jk,jwalk) = zundsat(ji,jk) 
    250          ENDDO 
    251       ENDDO 
    252        
    253       !---------------------------------- 
    254       !   Back to initial geometry 
    255       !----------------------------- 
    256        
    257       !--------------------------------------------------------------------- 
    258       !   1/ Compensation for ajustement of the bottom water concentrations 
    259       !      (see note n° 1 about *por(2)) 
    260       !-------------------------------------------------------------------- 
    261       DO jw = 1, jpwat 
    262          DO ji = 1, jpoce 
    263             pwcp(ji,1,jw) = pwcp(ji,1,jw) + & 
    264                &            pwcp(ji,2,jw) * dzdep(ji) * por(2) / dzkbot(ji) 
    265          END DO 
    266       ENDDO 
    267        
    268       !----------------------------------------------------------------------- 
    269       !    2/ Det of new rainrg taking account of the new weight fraction obtained  
    270       !      in dz3d(2) after diffusion/reaction (react/diffu are also in dzdep!) 
    271       !      This new rain (rgntg rm) will be used in advection/burial routine 
    272       !------------------------------------------------------------------------ 
    273       DO js = 1, jpsol 
    274          DO ji = 1, jpoce 
    275             rainrg(ji,js) = raintg(ji) * solcp(ji,2,js) 
    276             rainrm(ji,js) = rainrg(ji,js) / mol_wgt(js) 
    277          END DO 
    278       ENDDO 
    279  
    280       !  New raintg 
    281       raintg(:) = 0. 
    282       DO js = 1, jpsol 
    283          DO ji = 1, jpoce 
    284             raintg(ji) = raintg(ji) + rainrg(ji,js) 
    285          END DO 
    286       ENDDO 
    287        
    288       !-------------------------------- 
    289       !    3/ back to initial geometry 
    290       !-------------------------------- 
    291       DO ji = 1, jpoce 
    292          dz3d  (ji,2) = dz(2) 
    293          volw3d(ji,2) = dz3d(ji,2) * por(2) 
    294          vols3d(ji,2) = dz3d(ji,2) * por1(2) 
    295       ENDDO 
    296126 
    297127      IF( ln_timing )  CALL timing_stop('sed_inorg') 
Note: See TracChangeset for help on using the changeset viewer.