source: NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/sedsol.F90 @ 15075

MODULE sedsol
   !!======================================================================
   !!              ***  MODULE  sedsol  ***
   !!    Sediment : dissolution and reaction in pore water related 
   !!    related to organic matter
   !!    Diffusion of solutes in pore water
   !!=====================================================================
   !! * Modules used
   USE sed     ! sediment global variable
   USE sed_oce
   USE sedini
   USE sedfunc
   USE seddsr
   USE sedjac
   USE sedbtb
   USE sedco3
   USE sedmat
   USE TROS_F90
   USE lib_mpp         ! distribued memory computing library
   USE lib_fortran

   IMPLICIT NONE
   PRIVATE

   PUBLIC sed_sol

   !! $Id: sedsol.F90 5215 2015-04-15 16:11:56Z nicolasmartin $
CONTAINS
   
   SUBROUTINE sed_sol( kt ) 
      !!----------------------------------------------------------------------
      !!                   ***  ROUTINE sed_sol  ***
      !! 
      !!  ** Purpose :  computes pore water diffusion and reactions
      !!
      !!  ** Methode :  Computation of the redox and dissolution reactions 
      !!                in the sediment.
      !!                The main redox reactions are solved in sed_dsr whereas
      !!                the secondary reactions are solved in sed_dsr_redoxb.
      !!                Inorganic dissolution is solved in sed_inorg
      !!                A strand spliting approach is being used here (see 
      !!                sed_dsr_redoxb for more information). 
      !!                Diffusive fluxes are computed in sed_diff
      !!
      !!   History :
      !!        !  98-08 (E. Maier-Reimer, Christoph Heinze )  Original code
      !!        !  04-10 (N. Emprin, M. Gehlen ) f90
      !!        !  06-04 (C. Ethe)  Re-organization
      !!        !  19-08 (O. Aumont) Debugging and improvement of the model.
      !!                             The original method is replaced by a 
      !!                             Strand splitting method which deals 
      !!                             well with stiff reactions.
      !!----------------------------------------------------------------------
      !! Arguments
      INTEGER, INTENT(in) ::   kt
      ! --- local variables
      INTEGER  :: ji, jk, js, jw, jn, neq   ! dummy looop indices
      REAL(wp), DIMENSION( jpvode * jpksed ) :: ZXIN
      REAL(wp), DIMENSION(jpoce,jpksed) :: preac
      INTEGER :: JINDEX, ITOL, IJAC, MLJAC, IMAX
      INTEGER :: MUJAC,LE1, LJAC, LIWORK, LWORK
      INTEGER :: IDID, NMAXSTP, ROSM
      REAL(wp) :: X, XEND, H, STEPMIN
      REAL(wp), DIMENSION(1) :: RTOL, ATOL
      REAL(wp), POINTER :: WORK(:)
      INTEGER , POINTER :: IWORK(:)
      INTEGER, DIMENSION(3)   :: ISTAT
      REAL(wp), DIMENSION(2)   :: RSTAT
      !!
      !!----------------------------------------------------------------------

      IF( ln_timing )  CALL timing_start('sed_sol')
!
      IF( kt == nitsed000 ) THEN
         IF (lwp) THEN
            WRITE(numsed,*) ' sed_sol : Organic/inorganic degradation related reactions and diffusion'
            WRITE(numsed,*) ' '
         ENDIF
!         ! 
         dens_mol_wgt(1:jpsol) = denssol / mol_wgt(1:jpsol)
         stepros(:) = dtsed
         ! 
      ENDIF

      ! 1. Change of geometry
      !    Increase of dz3d(2) thickness : dz3d(2) = dz3d(2)+dzdep
      !    Warning : no change for dz(2)
      !---------------------------------------------------------
      dz3d(1:jpoce,2) = dz3d(1:jpoce,2) + dzdep(1:jpoce)

      ! New values for volw3d(:,2) and vols3d(:,2)
      ! Warning : no change neither for volw(2) nor  vols(2)
      !------------------------------------------------------
      volw3d(1:jpoce,2) = dz3d(1:jpoce,2) * por(2)
      vols3d(1:jpoce,2) = dz3d(1:jpoce,2) * por1(2)

      ! 2. Change of previous solid fractions (due to volum changes) for k=2
      !---------------------------------------------------------------------
      DO js = 1, jpsol
         solcp(:,2,js) = solcp(:,2,js) * dz(2) / dz3d(:,2)
      END DO

      ! 3. New solid fractions (including solid rain fractions) for k=2
      !------------------------------------------------------------------
      DO js = 1, jpsol
         DO ji = 1, jpoce
            IF (raintg(ji) .ne. 0) THEN
               solcp(ji,2,js) = solcp(ji,2,js) + &
               &           ( rainrg(ji,js) / raintg(ji) ) * ( dzdep(ji) / dz3d(ji,2) )
               ! rainrm are temporary cancel
               rainrm(ji,js) = 0.
            ENDIF
         END DO
      ENDDO

      ! 4.  Adjustment of bottom water concen.(pwcp(1)):
      !     We impose that pwcp(2) is constant. Including dzdep in dz3d(:,2) we assume
      !     that dzdep has got a porosity of por(2). So pore water volum of jk=2 increase.
      !     To keep pwcp(2) cste we must compensate this "increase" by a slight adjusment
      !     of bottom water concentration.
      !     This adjustment is compensate at the end of routine
      !-------------------------------------------------------------
      DO jw = 1, jpwat
         pwcp(:,1,jw) = pwcp(:,1,jw) - pwcp(:,2,jw) * dzdep(:) * por(2) / dzkbot(:)
      ENDDO

      ! --------------------------------------------------
      ! Computation of the diffusivities
      ! --------------------------------------------------
      DO js = 1, jpwat
         DO jk = 1, jpksed
            diff(:,jk,js) = ( diff1(js) + diff2(js) * temp(:) ) / ( 1.0 - 2.0 * log( por(jk) ) )
         END DO
      END DO

      ! Impact of bioirrigation and adsorption on diffusion
      ! ---------------------------------------------------
      diff(:,:,jwsil) = diff(:,:,jwsil) * ( 1.0 + irrig(:,:) )
      diff(:,:,jwoxy) = diff(:,:,jwoxy) * ( 1.0 + irrig(:,:) )
      diff(:,:,jwdic) = diff(:,:,jwdic) * ( 1.0 + irrig(:,:) )
      diff(:,:,jwno3) = diff(:,:,jwno3) * ( 1.0 + irrig(:,:) )
      diff(:,:,jwpo4) = diff(:,:,jwpo4) * ( 1.0 + irrig(:,:) )
      diff(:,:,jwalk) = diff(:,:,jwalk) * ( 1.0 + irrig(:,:) )
      diff(:,:,jwh2s) = diff(:,:,jwh2s) * ( 1.0 + irrig(:,:) )
      diff(:,:,jwso4) = diff(:,:,jwso4) * ( 1.0 + irrig(:,:) )
      diff(:,:,jwlgw) = diff(:,:,jwlgw) * ( 1.0 + irrig(:,:) )
      DO jk = 1, jpksed
         diff(:,jk,jwfe2) = diff(:,jk,jwfe2) * ( 1.0 + 0.1 * irrig(:,jk) ) * radsfe2(jk)
         diff(:,jk,jwnh4) = diff(:,jk,jwnh4) * ( 1.0 + irrig(:,jk) ) * radsnh4(jk)
      END DO

      ! Conversion of volume units
      !----------------------------
      DO js = 1, jpsol
         DO jk = 1, jpksed
            volc(:,jk,js) = ( vols3d(:,jk) * dens_mol_wgt(js) ) /  &
            &                 ( volw3d(:,jk) * 1.e-3 )
         ENDDO
      ENDDO

      ! Apply bioturbation and compute the impact of the slow SMS on species
      CALL sed_btb( kt )
      ! Recompute pH after bioturbation and slow chemistry
      CALL sed_co3( kt )

      ! The following part deals with the stiff ODEs
      ! This is the expensive part of the code and should be carefully
      ! chosen. We use the DVODE solver after many trials to find a cheap 
      ! way to solve the ODEs. This is not necessarily the most efficient 
      ! but this is the one that was not too much of a pain to code and that
      ! was the most precise and quick.
      ! The ones I tried : operator splitting (Strang), hybrid spectral methods
      ! Brent, Powell's hybrid method, ...
      ! ---------------------------------------------------------------------
      ROSM = 2
      NEQ  = jpvode * jpksed
      XEND = dtsed 
      RTOL = 0.005
      ATOL = 0.005
      ITOL = 0
      IJAC = 1
      MLJAC = jpvode
      MUJAC = jpvode
      LE1  = 2*MLJAC+MUJAC+1
      LJAC = MLJAC+MUJAC+1
      LIWORK = NEQ + 2
      LWORK = NEQ*(LJAC+LE1+8) + 5
      ALLOCATE(WORK(LWORK), IWORK(LIWORK) )
      NMAXSTP = 0
      STEPMIN = dtsed

      DO ji = 1, jpoce
         X    = 0.0
         H    = stepros(ji)
         WORK  = 0.
         IWORK = 0
         IWORK(2) = 6

         ! Put all the species in one local array (nb of tracers * vertical
         ! dimension
         jindex = ji
         DO jn = 1, NEQ
            jk = jarr(jn,1)
            js = jarr(jn,2)
            IF (js <= jpwat) THEN
               zxin(jn) = pwcp(ji,jk,js) * 1E6
            ELSE
               zxin(jn) = solcp(ji,jk,js-jpwat) * 1E6
            ENDIF
         END DO

         ! Set options for VODE : banded matrix. SParse option is much more
         ! expensive except if one computes the sparse Jacobian explicitly
         ! To speed up the computation, one way is to reduce ATOL and RTOL
         ! which may be a good option at the beginning of the simulations 
         ! during the spin up
         ! ----------------------------------------------------------------
          CALL ROSK(ROSM, NEQ,JINDEX,X,zxin,XEND,H,  &
          RTOL,ATOL,ITOL, sed_jac,IJAC,MLJAC,MUJAC,  &
          WORK,LWORK,IWORK,LIWORK,IDID,ISTAT,RSTAT)

         DO jn = 1, NEQ
            jk = jarr(jn,1)
            js = jarr(jn,2)
            IF (js <= jpwat) THEN
               pwcp(ji,jk,js) = zxin(jn) * 1E-6
            ELSE
               solcp(ji,jk,js-jpwat) = zxin(jn) * 1E-6
            ENDIF
         END DO
         NMAXSTP = MAX(NMAXSTP,ISTAT(3))
         STEPMIN = MIN(STEPMIN, RSTAT(1) )
         stepros(ji) = RSTAT(1)
      END DO

      DEALLOCATE(WORK, IWORK)

      preac(:,:) = 0.
      CALL sed_mat_dsri( jwalk, preac, pwcpa(:,:,jwalk), dtsed )
      CALL sed_mat_dsri( jwpo4, preac, pwcpa(:,:,jwpo4), dtsed )

      !----------------------------------
      !   Back to initial geometry
      !-----------------------------

      !---------------------------------------------------------------------
      !   1/ Compensation for ajustement of the bottom water concentrations
      !      (see note n� 1 about *por(2))
      !--------------------------------------------------------------------
      DO jw = 1, jpwat
         pwcp(:,1,jw) = pwcp(:,1,jw) + pwcp(:,2,jw) * dzdep(:) * por(2) / dzkbot(:)
      ENDDO

      !-----------------------------------------------------------------------
      !    2/ Det of new rainrg taking account of the new weight fraction
      !    obtained 
      !      in dz3d(2) after diffusion/reaction (react/diffu are also in
      !      dzdep!)
      !      This new rain (rgntg rm) will be used in advection/burial routine
      !------------------------------------------------------------------------
      DO js = 1, jpsol
         rainrg(:,js) = raintg(:) * solcp(:,2,js)
         rainrm(:,js) = rainrg(:,js) / mol_wgt(js)
      ENDDO

      !  New raintg
      raintg(:) = 0.
      DO js = 1, jpsol
         raintg(:) = raintg(:) + rainrg(:,js)
      ENDDO

      !--------------------------------
      !    3/ back to initial geometry
      !--------------------------------
      dz3d  (:,2) = dz(2)
      volw3d(:,2) = dz3d(:,2) * por(2)
      vols3d(:,2) = dz3d(:,2) * por1(2)

      IF( ln_timing )  CALL timing_stop('sed_sol')
!      
   END SUBROUTINE sed_sol

   SUBROUTINE JACFAC
  
   END SUBROUTINE JACFAC

END MODULE sedsol