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 15115 – NEMO

Changeset 15115


Ignore:
Timestamp:
2021-07-12T17:00:27+02:00 (3 years ago)
Author:
aumont
Message:

bug fixes in the sediment module

Location:
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OFF/nemogcm.F90

    r15075 r15115  
    159159                                CALL dta_dyn_sed( istp,      Nnn      )       ! Interpolation of the dynamical fields 
    160160                                CALL trc_stp    ( istp, Nbb, Nnn, Nrhs, Naa ) ! time-stepping 
    161         ! Swap time levels 
    162          Nrhs = Nbb 
    163          Nbb  = Nnn 
    164          Nnn  = Naa 
    165          Naa  = Nrhs 
     161         ! Swap time levels 
     162         Nnn = Nbb 
     163         Naa = Nbb 
    166164#endif 
    167165         CALL stp_ctl    ( istp )             ! Time loop: control and print 

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

    r15075 r15115  
    7171       jn = nvar 
    7272       ! first sediment level           
    73        aplus  = ( ( volw3d(ji,1) / ( dz3d(ji,1) ) ) + & 
    74                         ( volw3d(ji,2) / ( dz3d(ji,2) ) ) ) / 2. 
     73       aplus  = ( por(1) + por(2) ) * 0.5 
    7574       dxplus = ( dz3d(ji,1) + dz3d(ji,2) ) / 2. 
    7675       rplus  = ( dtsed_in / ( volw3d(ji,1) ) ) * diff(ji,1,jn) * aplus / dxplus  
    … …  
    8180  
    8281       DO jk = 2, jpksed - 1 
    83           aminus  = ( ( volw3d(ji,jk-1) / ( dz3d(ji,jk-1) ) ) + & 
    84           &        ( volw3d(ji,jk  ) / ( dz3d(ji,jk  ) ) ) ) / 2. 
     82          aminus  = ( por(jk-1) + por(jk) ) * 0.5 
    8583          dxminus = ( dz3d(ji,jk-1) + dz3d(ji,jk) ) / 2. 
    8684 
    87           aplus   = ( ( volw3d(ji,jk  ) / ( dz3d(ji,jk  ) ) ) + & 
    88           &        ( volw3d(ji,jk+1) / ( dz3d(ji,jk+1) ) ) ) / 2. 
     85          aplus   = ( por(jk+1) + por(jk) ) * 0.5 
    8986          dxplus  = ( dz3d(ji,jk) + dz3d(ji,jk+1) ) / 2 
    9087          ! 
    … …  
    9794       END DO 
    9895 
    99        aminus  = ( ( volw3d(ji,jpksed-1) / dz3d(ji,jpksed-1)  ) + & 
    100        &        ( volw3d(ji,jpksed)  / dz3d(ji,jpksed) ) ) / 2. 
     96       aminus  = ( por(jpksed-1) + por(jpksed) ) * 0.5 
    10197       dxminus = ( dz3d(ji,jpksed-1) + dz3d(ji,jpksed) ) / 2. 
    10298       rminus  = ( dtsed_in / volw3d(ji,jpksed) ) * diff(ji,jpksed-1,jn) * aminus / dxminus 
    … …  
    167163       jn = nvar 
    168164       ! first sediment level           
    169        aplus  = ( ( volw3d(ji,1) / ( dz3d(ji,1) ) ) + & 
    170                         ( volw3d(ji,2) / ( dz3d(ji,2) ) ) ) / 2. 
     165       aplus  = ( por(1) + por(2) ) *0.5 
    171166       dxplus = ( dz3d(ji,1) + dz3d(ji,2) ) / 2. 
    172167       rplus  = ( dtsed_in / ( volw3d(ji,1) ) ) * diff(ji,1,jn) * aplus / dxplus 
    … …  
    177172 
    178173       DO jk = 2, jpksed - 1 
    179           aminus  = ( ( volw3d(ji,jk-1) / ( dz3d(ji,jk-1) ) ) + & 
    180           &        ( volw3d(ji,jk  ) / ( dz3d(ji,jk  ) ) ) ) / 2. 
     174          aminus  = ( por(jk-1) + por(jk) ) * 0.5 
    181175          dxminus = ( dz3d(ji,jk-1) + dz3d(ji,jk) ) / 2. 
    182176 
    183           aplus   = ( ( volw3d(ji,jk  ) / ( dz3d(ji,jk  ) ) ) + & 
    184           &        ( volw3d(ji,jk+1) / ( dz3d(ji,jk+1) ) ) ) / 2. 
     177          aplus   = ( por(jk+1) + por(jk) ) * 0.5 
    185178          dxplus  = ( dz3d(ji,jk) + dz3d(ji,jk+1) ) / 2 
    186179          ! 
    … …  
    193186       END DO 
    194187 
    195        aminus  = ( ( volw3d(ji,jpksed-1) / dz3d(ji,jpksed-1)  ) + & 
    196        &        ( volw3d(ji,jpksed)  / dz3d(ji,jpksed) ) ) / 2. 
     188       aminus  = ( por(jpksed-1) + por(jpksed) ) * 0.5 
    197189       dxminus = ( dz3d(ji,jpksed-1) + dz3d(ji,jpksed) ) / 2. 
    198190       rminus  = ( dtsed_in / volw3d(ji,jpksed) ) * diff(ji,jpksed-1,jn) * aminus / dxminus 
    … …  
    274266       ! first sediment level           
    275267      DO ji = 1, jpoce 
    276          aplus  = ( ( vols(2) / dz(2) ) + ( vols(3) / dz(3) ) ) / 2. 
     268         aplus  = ( por1(2) + por1(3) ) / 2.0 
    277269         dxplus = ( dz(2) + dz(3) ) / 2. 
    278270         rplus  = ( dtsed_in / vols(2) ) * db(ji,2) * aplus / dxplus 
    … …  
    283275 
    284276         DO jk = 2, nlev - 1 
     277            aminus  = ( por1(jk) + por1(jk+1) ) * 0.5 
    285278            aminus  = ( ( vols(jk) / dz(jk) ) + ( vols(jk+1) / dz(jk+1) ) ) / 2. 
    286279            dxminus = ( dz(jk) + dz(jk+1) ) / 2. 
    287280            rminus  = ( dtsed_in / vols(jk+1) ) * db(ji,jk) * aminus / dxminus 
    288281            ! 
    289             aplus   = ( ( vols(jk+1) / dz(jk+1  ) ) + ( vols(jk+2) / dz(jk+2) ) ) / 2. 
     282            aplus   = ( por1(jk+1) + por1(jk+2) ) * 0.5 
    290283            dxplus  = ( dz(jk+1) + dz(jk+2) ) / 2. 
    291284            rplus   = ( dtsed_in / vols(jk+1) ) * db(ji,jk+1) * aplus / dxplus 
    … …  
    297290         ENDDO 
    298291 
    299          aminus  = ( ( vols(nlev) / dz(nlev) ) + ( vols(nlev+1) / dz(nlev+1) ) ) / 2. 
     292         aminus = ( por1(nlev) + por1(nlev+1) ) * 0.5 
    300293         dxminus = ( dz(nlev) + dz(nlev+1) ) / 2. 
    301294         rminus  = ( dtsed_in / vols(nlev+1) ) * db(ji,nlev) * aminus / dxminus 
    … …  
    383376       ! first sediment level           
    384377       DO ji = 1, jpoce 
    385           aplus  = ( ( volw3d(ji,1) / ( dz3d(ji,1) ) ) + & 
    386                         ( volw3d(ji,2) / ( dz3d(ji,2) ) ) ) / 2. 
     378          aplus  = ( por(1) + por(2) ) * 0.5 
    387379          dxplus = ( dz3d(ji,1) + dz3d(ji,2) ) / 2. 
    388380          rplus  = ( dtsed_in / ( volw3d(ji,1) ) ) * diff(ji,1,jn) * aplus / dxplus 
    … …  
    395387       DO jk = 2, jpksed - 1 
    396388          DO ji = 1, jpoce 
    397              aminus  = ( ( volw3d(ji,jk-1) / ( dz3d(ji,jk-1) ) ) + & 
    398              &        ( volw3d(ji,jk  ) / ( dz3d(ji,jk  ) ) ) ) / 2. 
     389             aminus  = ( por(jk-1) + por(jk) ) * 0.5 
    399390             dxminus = ( dz3d(ji,jk-1) + dz3d(ji,jk) ) / 2. 
    400391 
    401              aplus   = ( ( volw3d(ji,jk  ) / ( dz3d(ji,jk  ) ) ) + & 
    402              &        ( volw3d(ji,jk+1) / ( dz3d(ji,jk+1) ) ) ) / 2. 
     392             aplus   = ( por(jk+1) + por(jk) ) * 0.5  
    403393             dxplus  = ( dz3d(ji,jk) + dz3d(ji,jk+1) ) / 2 
    404394                ! 
    … …  
    413403 
    414404       DO ji = 1, jpoce 
    415           aminus  = ( ( volw3d(ji,jpksed-1) / dz3d(ji,jpksed-1)  ) + & 
    416           &        ( volw3d(ji,jpksed)  / dz3d(ji,jpksed) ) ) / 2. 
     405          aminus  = ( por(jpksed-1) + por(jpksed) ) *0.5 
    417406          dxminus = ( dz3d(ji,jpksed-1) + dz3d(ji,jpksed) ) / 2. 
    418407          rminus  = ( dtsed_in / volw3d(ji,jpksed) ) * diff(ji,jpksed-1,jn) * aminus/ dxminus 

  • NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/sedsol.F90

    r15075 r15115  
    1616   USE sedco3 
    1717   USE sedmat 
    18    USE TROS_F90 
     18   USE tros 
    1919   USE lib_mpp         ! distribued memory computing library 
    2020   USE lib_fortran 

  • NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/trosk.F90

    r15075 r15115  
    1 MODULE TROS_F90 
     1MODULE TROS 
    22!**************************************************************** 
    33!* NUMERICAL SOLUTION OF A STIFF SYSTEM OF FIRST 0RDER ORDINARY * 
    … …  
    4444!*                                (www.jpmoreau.fr)             * 
    4545!**************************************************************** 
     46  USE sed 
    4647  USE sedfunc 
    4748 
    48   INTEGER, PARAMETER, PRIVATE :: WP = KIND(1.0D0) 
    4949  INTEGER, PRIVATE :: JINDIC 
    5050 
    … …  
    275275         NMAX=IWORK(1) 
    276276         IF(NMAX.LE.0)THEN 
    277             WRITE(6,*)' WRONG INPUT IWORK(1)=',IWORK(1) 
     277            WRITE(NUMSED,*)' WRONG INPUT IWORK(1)=',IWORK(1) 
    278278            ARRET=.TRUE. 
    279279         END IF 
    … …  
    285285         METH=IWORK(2) 
    286286         IF(METH.LE.0.OR.METH.GE.7)THEN 
    287             WRITE(6,*)' CURIOUS INPUT IWORK(2)=',IWORK(2) 
     287            WRITE(NUMSED,*)' CURIOUS INPUT IWORK(2)=',IWORK(2) 
    288288            ARRET=.TRUE. 
    289289         END IF 
    … …  
    295295         UROUND=WORK(1) 
    296296         IF(UROUND.LE.1.D-14.OR.UROUND.GE.1.D0)THEN 
    297             WRITE(6,*)' COEFFICIENTS HAVE 16 DIGITS, UROUND=',WORK(1) 
     297            WRITE(NUMSED,*)' COEFFICIENTS HAVE 16 DIGITS, UROUND=',WORK(1) 
    298298            ARRET=.TRUE. 
    299299         END IF 
    … …  
    325325      IF (ITOL.EQ.0) THEN 
    326326          IF (ATOL(1).LE.0.D0.OR.RTOL(1).LE.10.D0*UROUND) THEN 
    327               WRITE (6,*) ' TOLERANCES ARE TOO SMALL' 
     327              WRITE (NUMSED,*) ' TOLERANCES ARE TOO SMALL' 
    328328              ARRET=.TRUE. 
    329329          END IF 
    … …  
    331331          DO 15 I=1,N 
    332332          IF (ATOL(I).LE.0.D0.OR.RTOL(I).LE.10.D0*UROUND) THEN 
    333               WRITE (6,*) ' TOLERANCES(',I,') ARE TOO SMALL' 
     333              WRITE (NUMSED,*) ' TOLERANCES(',I,') ARE TOO SMALL' 
    334334              ARRET=.TRUE. 
    335335          END IF 
    … …  
    370370      ISTORE=IEE+N*LDE-1 
    371371      IF(ISTORE.GT.LWORK)THEN 
    372          WRITE(6,*)' INSUFFICIENT STORAGE FOR WORK, MIN. LWORK=',ISTORE 
     372         WRITE(NUMSED,*)' INSUFFICIENT STORAGE FOR WORK, MIN. LWORK=',ISTORE 
    373373         ARRET=.TRUE. 
    374374      END IF 
    … …  
    378378      ISTORE=IEIP+N-1 
    379379      IF(ISTORE.GT.LIWORK)THEN 
    380          WRITE(6,*)' INSUFF. STORAGE FOR IWORK, MIN. LIWORK=',ISTORE 
     380         WRITE(NUMSED,*)' INSUFF. STORAGE FOR IWORK, MIN. LIWORK=',ISTORE 
    381381         ARRET=.TRUE. 
    382382      END IF 
    … …  
    661661      END IF 
    662662! --- EXIT 
    663   80  WRITE (6,*) ' MATRIX E IS SINGULAR, INFO = ',INFO 
     663  80  WRITE (NUMSED,*) ' MATRIX E IS SINGULAR, INFO = ',INFO 
    664664      NSING=NSING+1 
    665665      IF (NSING.GE.5) GOTO 79 
    666666      H=H*0.5D0 
    667667      GOTO 2 
    668   79  WRITE(6,979)X,H 
     668  79  WRITE(NUMSED,979)X,H 
    669669 979  FORMAT(' EXIT OF ROS4 AT X=',D16.7,'   H=',D16.7) 
    670670      IDID=-1 
    … …  
    862862! --- COMPUTATION OF HNEW 
    863863! --- WE REQUIRE .2<=HNEW/H<=6. 
    864       FAC=DMAX1(FAC2,DMIN1(FAC1,(ERR)**.33D0/.9D0)) 
     864      FAC=DMAX1(FAC2,DMIN1(FAC1,(ERR)**.33D0/0.9D0)) 
    865865      HNEW=H/FAC 
    866866! *** *** *** *** *** *** *** 
    … …  
    894894      END IF 
    895895! --- EXIT 
    896   80  WRITE (6,*) ' MATRIX E IS SINGULAR, INFO = ',INFO 
     896  80  WRITE (NUMSED,*) ' MATRIX E IS SINGULAR, INFO = ',INFO 
    897897      NSING=NSING+1 
    898898      IF (NSING.GE.5) GOTO 79 
    899899      H=H*0.5D0 
    900900      GOTO 2 
    901   79  WRITE(6,979)X,H 
     901  79  WRITE(NUMSED,979)X,H 
    902902 979  FORMAT(' EXIT OF ROS4 AT X=',D16.7,'   H=',D16.7) 
    903903      IDID=-1 
    … …  
    11121112      END IF 
    11131113! --- EXIT 
    1114   80  WRITE (6,*) ' MATRIX E IS SINGULAR, INFO = ',INFO 
     1114  80  WRITE (NUMSED,*) ' MATRIX E IS SINGULAR, INFO = ',INFO 
    11151115      NSING=NSING+1 
    11161116      IF (NSING.GE.5) GOTO 79 
    11171117      H=H*0.5D0 
    11181118      GOTO 2 
    1119   79  WRITE(6,979)X,H 
     1119  79  WRITE(NUMSED,979)X,H 
    11201120 979  FORMAT(' EXIT OF ROS2 AT X=',D16.7,'   H=',D16.7) 
    11211121      IDID=-1 
    … …  
    15511551      END SUBROUTINE SOLB 
    15521552 
    1553 END MODULE TROS_F90 
     1553END MODULE TROS 
Note: See TracChangeset for help on using the changeset viewer.