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 8882 for branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf.F90 – NEMO

Ignore:
Timestamp:
2017-12-01T18:44:09+01:00 (6 years ago)
Author:
flavoni
Message:

dev_CNRS_2017 branch: merged dev_r7881_ENHANCE09_RK3 with trunk r8864

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf.F90

    r8698 r8882  
    44   !! Ocean active tracers:  vertical component of the tracer mixing trend 
    55   !!============================================================================== 
    6    !! History :  1.0  ! 2005-11  (G. Madec)  Original code 
    7    !!            3.0  ! 2008-01  (C. Ethe, G. Madec)  merge TRC-TRA 
     6   !! History :  1.0  !  2005-11  (G. Madec)  Original code 
     7   !!            3.0  !  2008-01  (C. Ethe, G. Madec)  merge TRC-TRA 
     8   !!            4.0  !  2017-06  (G. Madec)  remove explict time-stepping option 
    89   !!---------------------------------------------------------------------- 
    910 
    1011   !!---------------------------------------------------------------------- 
    1112   !!   tra_zdf       : Update the tracer trend with the vertical diffusion 
    12    !!   tra_zdf_init  : initialisation of the computation 
    1313   !!---------------------------------------------------------------------- 
    1414   USE oce            ! ocean dynamics and tracers variables 
     
    2020   USE ldftra         ! lateral diffusion: eddy diffusivity 
    2121   USE ldfslp         ! lateral diffusion: iso-neutral slope  
    22    USE trazdf_exp     ! vertical diffusion: explicit (tra_zdf_exp routine) 
    23    USE trazdf_imp     ! vertical diffusion: implicit (tra_zdf_imp routine) 
    2422   USE trd_oce        ! trends: ocean variables 
    2523   USE trdtra         ! trends: tracer trend manager 
     
    2927   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
    3028   USE lib_mpp        ! MPP library 
    31    USE wrk_nemo       ! Memory allocation 
    3229   USE timing         ! Timing 
    3330 
     
    3532   PRIVATE 
    3633 
    37    PUBLIC   tra_zdf        ! routine called by step.F90 
    38    PUBLIC   tra_zdf_init   ! routine called by nemogcm.F90 
    39  
    40    INTEGER ::   nzdf = 0   ! type vertical diffusion algorithm used (defined from ln_zdf...  namlist logicals) 
     34   PUBLIC   tra_zdf       ! called by step.F90 
     35   PUBLIC   tra_zdf_imp   ! called by trczdf.F90 
    4136 
    4237   !! * Substitutions 
    43 #  include "zdfddm_substitute.h90" 
    4438#  include "vectopt_loop_substitute.h90" 
    4539   !!---------------------------------------------------------------------- 
    46    !! NEMO/OPA 3.7 , NEMO Consortium (2015) 
     40   !! NEMO/OPA 4.0 , NEMO Consortium (2017) 
    4741   !! $Id$ 
    4842   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    5650      !! ** Purpose :   compute the vertical ocean tracer physics. 
    5751      !!--------------------------------------------------------------------- 
    58       INTEGER, INTENT( in ) ::   kt      ! ocean time-step index 
    59       ! 
    60       INTEGER  ::   jk                   ! Dummy loop indices 
    61       REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztrdt, ztrds   ! 3D workspace 
    62       !!--------------------------------------------------------------------- 
    63       ! 
    64       IF( nn_timing == 1 )  CALL timing_start('tra_zdf') 
    65       ! 
    66       IF( neuler == 0 .AND. kt == nit000 ) THEN     ! at nit000 
    67          r2dt =  rdt                          ! = rdt (restarting with Euler time stepping) 
    68       ELSEIF( kt <= nit000 + 1) THEN                ! at nit000 or nit000+1 
    69          r2dt = 2. * rdt                      ! = 2 rdt (leapfrog) 
    70       ENDIF 
    71       ! 
    72       IF( l_trdtra )   THEN                    !* Save ta and sa trends 
    73          CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds ) 
     52      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     53      ! 
     54      INTEGER  ::   jk   ! Dummy loop indices 
     55      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdt, ztrds   ! 3D workspace 
     56      !!--------------------------------------------------------------------- 
     57      ! 
     58      IF( ln_timing )   CALL timing_start('tra_zdf') 
     59      ! 
     60      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dt =      rdt   ! at nit000, =   rdt (restarting with Euler time stepping) 
     61      ELSEIF( kt <= nit000 + 1           ) THEN   ;   r2dt = 2. * rdt   ! otherwise, = 2 rdt (leapfrog) 
     62      ENDIF 
     63      ! 
     64      IF( l_trdtra )   THEN                  !* Save ta and sa trends 
     65         ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) ) 
    7466         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 
    7567         ztrds(:,:,:) = tsa(:,:,:,jp_sal) 
    7668      ENDIF 
    7769      ! 
    78       SELECT CASE ( nzdf )                       ! compute lateral mixing trend and add it to the general trend 
    79       CASE ( 0 )    ;    CALL tra_zdf_exp( kt, nit000, 'TRA', r2dt, nn_zdfexp, tsb, tsa, jpts )  !   explicit scheme  
    80       CASE ( 1 )    ;    CALL tra_zdf_imp( kt, nit000, 'TRA', r2dt,            tsb, tsa, jpts )  !   implicit scheme  
    81       END SELECT 
     70      !                                      !* compute lateral mixing trend and add it to the general trend 
     71      CALL tra_zdf_imp( kt, nit000, 'TRA', r2dt, tsb, tsa, jpts )  
     72 
    8273!!gm WHY here !   and I don't like that ! 
    8374      ! DRAKKAR SSS control { 
     
    9081         DO jk = 1, jpkm1 
    9182            ztrdt(:,:,jk) = ( ( tsa(:,:,jk,jp_tem)*e3t_a(:,:,jk) - tsb(:,:,jk,jp_tem)*e3t_b(:,:,jk) ) & 
    92                  & / (e3t_n(:,:,jk)*r2dt) ) - ztrdt(:,:,jk) 
     83               &          / (e3t_n(:,:,jk)*r2dt) ) - ztrdt(:,:,jk) 
    9384            ztrds(:,:,jk) = ( ( tsa(:,:,jk,jp_sal)*e3t_a(:,:,jk) - tsb(:,:,jk,jp_sal)*e3t_b(:,:,jk) ) & 
    94                  & / (e3t_n(:,:,jk)*r2dt) ) - ztrds(:,:,jk) 
     85              &          / (e3t_n(:,:,jk)*r2dt) ) - ztrds(:,:,jk) 
    9586         END DO 
    9687!!gm this should be moved in trdtra.F90 and done on all trends 
     
    10091         CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdf, ztrdt ) 
    10192         CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdf, ztrds ) 
    102          CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds ) 
     93         DEALLOCATE( ztrdt , ztrds ) 
    10394      ENDIF 
    10495      !                                          ! print mean trends (used for debugging) 
     
    10697         &                       tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 
    10798      ! 
    108       IF( nn_timing == 1 )  CALL timing_stop('tra_zdf') 
     99      IF( ln_timing )   CALL timing_stop('tra_zdf') 
    109100      ! 
    110101   END SUBROUTINE tra_zdf 
    111102 
    112  
    113    SUBROUTINE tra_zdf_init 
     103  
     104   SUBROUTINE tra_zdf_imp( kt, kit000, cdtype, p2dt, ptb, pta, kjpt )  
    114105      !!---------------------------------------------------------------------- 
    115       !!                 ***  ROUTINE tra_zdf_init  *** 
    116       !! 
    117       !! ** Purpose :   Choose the vertical mixing scheme 
    118       !! 
    119       !! ** Method  :   Set nzdf from ln_zdfexp 
    120       !!      nzdf = 0   explicit (time-splitting) scheme (ln_zdfexp=T) 
    121       !!           = 1   implicit (euler backward) scheme (ln_zdfexp=F) 
    122       !!      NB: rotation of lateral mixing operator or TKE & GLS schemes, 
    123       !!          an implicit scheme is required. 
    124       !!---------------------------------------------------------------------- 
    125       USE zdftke 
    126       USE zdfgls 
    127       !!---------------------------------------------------------------------- 
    128       ! 
    129       ! Choice from ln_zdfexp already read in namelist in zdfini module 
    130       IF( ln_zdfexp ) THEN   ;   nzdf = 0           ! use explicit scheme 
    131       ELSE                   ;   nzdf = 1           ! use implicit scheme 
    132       ENDIF 
    133       ! 
    134       ! Force implicit schemes 
    135       IF( lk_zdftke .OR. lk_zdfgls   )   nzdf = 1   ! TKE, or GLS physics 
    136       IF( ln_traldf_iso              )   nzdf = 1   ! iso-neutral lateral physics 
    137       IF( ln_traldf_hor .AND. ln_sco )   nzdf = 1   ! horizontal lateral physics in s-coordinate 
    138       IF( ln_zdfexp .AND. nzdf == 1 )   CALL ctl_stop( 'tra_zdf : If using the rotation of lateral mixing operator',   & 
    139             &                         ' GLS or TKE scheme, the implicit scheme is required, set ln_zdfexp = .false.' ) 
    140             ! 
    141       IF(lwp) THEN 
    142          WRITE(numout,*) 
    143          WRITE(numout,*) 'tra_zdf_init : vertical tracer physics scheme' 
    144          WRITE(numout,*) '~~~~~~~~~~~' 
    145          IF( nzdf ==  0 )   WRITE(numout,*) '      ===>>   Explicit time-splitting scheme' 
    146          IF( nzdf ==  1 )   WRITE(numout,*) '      ===>>   Implicit (euler backward) scheme' 
    147       ENDIF 
    148       ! 
    149    END SUBROUTINE tra_zdf_init 
     106      !!                  ***  ROUTINE tra_zdf_imp  *** 
     107      !! 
     108      !! ** Purpose :   Compute the after tracer through a implicit computation 
     109      !!     of the vertical tracer diffusion (including the vertical component  
     110      !!     of lateral mixing (only for 2nd order operator, for fourth order  
     111      !!     it is already computed and add to the general trend in traldf)  
     112      !! 
     113      !! ** Method  :  The vertical diffusion of a tracer ,t , is given by: 
     114      !!          difft = dz( avt dz(t) ) = 1/e3t dk+1( avt/e3w dk(t) ) 
     115      !!      It is computed using a backward time scheme (t=after field) 
     116      !!      which provide directly the after tracer field. 
     117      !!      If ln_zdfddm=T, use avs for salinity or for passive tracers 
     118      !!      Surface and bottom boundary conditions: no diffusive flux on 
     119      !!      both tracers (bottom, applied through the masked field avt). 
     120      !!      If iso-neutral mixing, add to avt the contribution due to lateral mixing. 
     121      !! 
     122      !! ** Action  : - pta  becomes the after tracer 
     123      !!--------------------------------------------------------------------- 
     124      INTEGER                              , INTENT(in   ) ::   kt       ! ocean time-step index 
     125      INTEGER                              , INTENT(in   ) ::   kit000   ! first time step index 
     126      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator) 
     127      INTEGER                              , INTENT(in   ) ::   kjpt     ! number of tracers 
     128      REAL(wp)                             , INTENT(in   ) ::   p2dt     ! tracer time-step 
     129      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb      ! before and now tracer fields 
     130      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta      ! in: tracer trend ; out: after tracer field 
     131      ! 
     132      INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices 
     133      REAL(wp) ::  zrhs             ! local scalars 
     134      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zwi, zwt, zwd, zws 
     135      !!--------------------------------------------------------------------- 
     136      ! 
     137      IF( ln_timing )   CALL timing_start('tra_zdf_imp') 
     138      ! 
     139      IF( kt == kit000 )  THEN 
     140         IF(lwp)WRITE(numout,*) 
     141         IF(lwp)WRITE(numout,*) 'tra_zdf_imp : implicit vertical mixing on ', cdtype 
     142         IF(lwp)WRITE(numout,*) '~~~~~~~~~~~ ' 
     143      ENDIF 
     144      !                                               ! ============= ! 
     145      DO jn = 1, kjpt                                 !  tracer loop  ! 
     146         !                                            ! ============= ! 
     147         !  Matrix construction 
     148         ! -------------------- 
     149         ! Build matrix if temperature or salinity (only in double diffusion case) or first passive tracer 
     150         ! 
     151         IF(  ( cdtype == 'TRA' .AND. ( jn == jp_tem .OR. ( jn == jp_sal .AND. ln_zdfddm ) ) ) .OR.   & 
     152            & ( cdtype == 'TRC' .AND. jn == 1 )  )  THEN 
     153            ! 
     154            ! vertical mixing coef.: avt for temperature, avs for salinity and passive tracers 
     155            IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN   ;   zwt(:,:,2:jpk) = avt(:,:,2:jpk) 
     156            ELSE                                            ;   zwt(:,:,2:jpk) = avs(:,:,2:jpk) 
     157            ENDIF 
     158            zwt(:,:,1) = 0._wp 
     159            ! 
     160            IF( l_ldfslp ) THEN            ! isoneutral diffusion: add the contribution  
     161               IF( ln_traldf_msc  ) THEN     ! MSC iso-neutral operator  
     162                  DO jk = 2, jpkm1 
     163                     DO jj = 2, jpjm1 
     164                        DO ji = fs_2, fs_jpim1   ! vector opt. 
     165                           zwt(ji,jj,jk) = zwt(ji,jj,jk) + akz(ji,jj,jk)   
     166                        END DO 
     167                     END DO 
     168                  END DO 
     169               ELSE                          ! standard or triad iso-neutral operator 
     170                  DO jk = 2, jpkm1 
     171                     DO jj = 2, jpjm1 
     172                        DO ji = fs_2, fs_jpim1   ! vector opt. 
     173                           zwt(ji,jj,jk) = zwt(ji,jj,jk) + ah_wslp2(ji,jj,jk) 
     174                        END DO 
     175                     END DO 
     176                  END DO 
     177               ENDIF 
     178            ENDIF 
     179            ! 
     180            ! Diagonal, lower (i), upper (s)  (including the bottom boundary condition since avt is masked) 
     181            DO jk = 1, jpkm1 
     182               DO jj = 2, jpjm1 
     183                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     184!!gm BUG  I think, use e3w_a instead of e3w_n, not sure of that 
     185                     zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk  ) / e3w_n(ji,jj,jk  ) 
     186                     zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / e3w_n(ji,jj,jk+1) 
     187                     zwd(ji,jj,jk) = e3t_a(ji,jj,jk) - zwi(ji,jj,jk) - zws(ji,jj,jk) 
     188                 END DO 
     189               END DO 
     190            END DO 
     191            ! 
     192            !! Matrix inversion from the first level 
     193            !!---------------------------------------------------------------------- 
     194            !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk ) 
     195            ! 
     196            !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 ) 
     197            !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 ) 
     198            !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 ) 
     199            !        (        ...               )( ...  ) ( ...  ) 
     200            !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk ) 
     201            ! 
     202            !   m is decomposed in the product of an upper and lower triangular matrix. 
     203            !   The 3 diagonal terms are in 3d arrays: zwd, zws, zwi. 
     204            !   Suffices i,s and d indicate "inferior" (below diagonal), diagonal 
     205            !   and "superior" (above diagonal) components of the tridiagonal system. 
     206            !   The solution will be in the 4d array pta. 
     207            !   The 3d array zwt is used as a work space array. 
     208            !   En route to the solution pta is used a to evaluate the rhs and then  
     209            !   used as a work space array: its value is modified. 
     210            ! 
     211            DO jj = 2, jpjm1        !* 1st recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k) 
     212               DO ji = fs_2, fs_jpim1            ! done one for all passive tracers (so included in the IF instruction) 
     213                  zwt(ji,jj,1) = zwd(ji,jj,1) 
     214               END DO 
     215            END DO 
     216            DO jk = 2, jpkm1 
     217               DO jj = 2, jpjm1 
     218                  DO ji = fs_2, fs_jpim1 
     219                     zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1) 
     220                  END DO 
     221               END DO 
     222            END DO 
     223            ! 
     224         ENDIF  
     225         !          
     226         DO jj = 2, jpjm1           !* 2nd recurrence:    Zk = Yk - Ik / Tk-1  Zk-1 
     227            DO ji = fs_2, fs_jpim1 
     228               pta(ji,jj,1,jn) = e3t_b(ji,jj,1) * ptb(ji,jj,1,jn) + p2dt * e3t_n(ji,jj,1) * pta(ji,jj,1,jn) 
     229            END DO 
     230         END DO 
     231         DO jk = 2, jpkm1 
     232            DO jj = 2, jpjm1 
     233               DO ji = fs_2, fs_jpim1 
     234                  zrhs = e3t_b(ji,jj,jk) * ptb(ji,jj,jk,jn) + p2dt * e3t_n(ji,jj,jk) * pta(ji,jj,jk,jn)   ! zrhs=right hand side 
     235                  pta(ji,jj,jk,jn) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pta(ji,jj,jk-1,jn) 
     236               END DO 
     237            END DO 
     238         END DO 
     239         ! 
     240         DO jj = 2, jpjm1           !* 3d recurrence:    Xk = (Zk - Sk Xk+1 ) / Tk   (result is the after tracer) 
     241            DO ji = fs_2, fs_jpim1 
     242               pta(ji,jj,jpkm1,jn) = pta(ji,jj,jpkm1,jn) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1) 
     243            END DO 
     244         END DO 
     245         DO jk = jpk-2, 1, -1 
     246            DO jj = 2, jpjm1 
     247               DO ji = fs_2, fs_jpim1 
     248                  pta(ji,jj,jk,jn) = ( pta(ji,jj,jk,jn) - zws(ji,jj,jk) * pta(ji,jj,jk+1,jn) )   & 
     249                     &             / zwt(ji,jj,jk) * tmask(ji,jj,jk) 
     250               END DO 
     251            END DO 
     252         END DO 
     253         !                                            ! ================= ! 
     254      END DO                                          !  end tracer loop  ! 
     255      !                                               ! ================= ! 
     256      ! 
     257      IF( ln_timing )   CALL timing_stop('tra_zdf_imp') 
     258      ! 
     259   END SUBROUTINE tra_zdf_imp 
    150260 
    151261   !!============================================================================== 
Note: See TracChangeset for help on using the changeset viewer.