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 2024 for branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/trazdf_imp.F90 – NEMO

Ignore:
Timestamp:
2010-07-29T12:57:35+02:00 (14 years ago)
Author:
cetlod
Message:

Merge of active and passive tracer advection/diffusion modules, see ticket:693

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/trazdf_imp.F90

    r1517 r2024  
    22   !!====================================================================== 
    33   !!                 ***  MODULE  trazdf_imp  *** 
    4    !! Ocean active tracers:  vertical component of the tracer mixing trend 
     4   !! Ocean tracers:  vertical component of the tracer mixing trend 
    55   !!====================================================================== 
    66   !! History :  OPA  !  1990-10  (B. Blanke)  Original code 
     
    1414   !!            2.0  !  2006-11  (G. Madec) New step reorganisation 
    1515   !!            3.2  !  2009-03  (G. Madec)  heat and salt content trends 
     16   !!            3.3  !  2010-06  (C. Ethe, G. Madec) Merge TRA-TRC 
    1617   !!---------------------------------------------------------------------- 
    1718   
     
    2526   USE ldftra_oce      ! ocean active tracers: lateral physics 
    2627   USE ldfslp          ! lateral physics: slope of diffusion 
    27    USE trdmod          ! ocean active tracers trends  
    28    USE trdmod_oce      ! ocean variables trends 
    2928   USE zdfddm          ! ocean vertical physics: double diffusion 
    3029   USE in_out_manager  ! I/O manager 
    3130   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    32    USE prtctl          ! Print control 
    3331   USE domvvl          ! variable volume 
    3432   USE ldftra          ! lateral mixing type 
     
    5048   !!---------------------------------------------------------------------- 
    5149CONTAINS 
    52     
    53    SUBROUTINE tra_zdf_imp( kt, p2dt ) 
     50  
     51   SUBROUTINE tra_zdf_imp( kt    , cdtype, p2dt,    & 
     52      &                    ptrab , ptraa , kjpt     ) 
    5453      !!---------------------------------------------------------------------- 
    5554      !!                  ***  ROUTINE tra_zdf_imp  *** 
     
    7170      !!                  associated with the lateral mixing, through the 
    7271      !!                  update of avt) 
    73       !!      The vertical diffusion of tracers (t & s) is given by: 
     72      !!      The vertical diffusion of the tracer t is given by: 
    7473      !!             difft = dz( avt dz(t) ) = 1/e3t dk+1( avt/e3w dk(t) ) 
    7574      !!      It is computed using a backward time scheme (t=ta). 
     
    7877      !!      Add this trend to the general trend ta,sa : 
    7978      !!         ta = ta + dz( avt dz(t) ) 
    80       !!         (sa = sa + dz( avs dz(t) ) if lk_zdfddm=T ) 
     79      !!         if lk_zdfddm=T, use avs for salinity or for passive tracers 
     80      !!         (sa = sa + dz( avs dz(t) )  
    8181      !! 
    8282      !!      Third part: recover avt resulting from the vertical physics 
    8383      !!      ==========  alone, for further diagnostics (for example to 
    8484      !!                  compute the turbocline depth in zdfmxl.F90). 
    85       !!         avt = zavt 
     85      !!         if lk_zdfddm=T, use avt = zavt 
    8686      !!         (avs = zavs if lk_zdfddm=T ) 
    8787      !! 
    88       !! ** Action  : - Update (ta,sa) with before vertical diffusion trend 
     88      !! ** Action  : - Update (ta) with before vertical diffusion trend 
    8989      !! 
    9090      !!--------------------------------------------------------------------- 
     91      !! * Modules used 
    9192      USE oce    , ONLY :   zwd   => ua   ! ua used as workspace 
    9293      USE oce    , ONLY :   zws   => va   ! va  -          - 
    93       !! 
    94       INTEGER                 , INTENT(in) ::   kt     ! ocean time-step index 
    95       REAL(wp), DIMENSION(jpk), INTENT(in) ::   p2dt   ! vertical profile of tracer time-step 
    96       !! 
    97       INTEGER  ::   ji, jj, jk            ! dummy loop indices 
    98       REAL(wp) ::   zavi, zrhs, znvvl     ! temporary scalars 
    99       REAL(wp) ::   ze3tb, ze3tn, ze3ta   ! variable vertical scale factors 
    100       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwi, zwt, zavsi   ! workspace arrays 
     94      !! * Arguments 
     95      INTEGER         , INTENT(in   )                                ::   kt             ! ocean time-step index 
     96      CHARACTER(len=3), INTENT(in   )                                ::   cdtype         ! =TRA or TRC (tracer indicator) 
     97      INTEGER         , INTENT(in   )                                ::   kjpt            ! number of tracers 
     98      REAL(wp)        , INTENT(in   ), DIMENSION(jpk)                ::   p2dt        ! vertical profile of tracer time-step 
     99      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk,kjpt)   ::   ptrab          ! before and now tracer fields 
     100      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)   ::   ptraa          ! tracer trend  
     101      !! 
     102      INTEGER  ::  ji, jj, jk, jn        ! dummy loop indices 
     103      REAL(wp) ::  zavi, zrhs, znvvl     ! temporary scalars 
     104      REAL(wp) ::  ze3tb, ze3tn, ze3ta   ! variable vertical scale factors 
     105      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwi, zwt   ! workspace arrays 
    101106      !!--------------------------------------------------------------------- 
    102107 
     
    107112         zavi = 0.e0      ! avoid warning at compilation phase when lk_ldfslp=F 
    108113      ENDIF 
    109  
     114      ! 
    110115      ! I. Local initialization 
    111116      ! ----------------------- 
    112       zwd  (1,:, : ) = 0.e0     ;     zwd  (jpi,:,:) = 0.e0 
    113       zws  (1,:, : ) = 0.e0     ;     zws  (jpi,:,:) = 0.e0 
    114       zwi  (1,:, : ) = 0.e0     ;     zwi  (jpi,:,:) = 0.e0 
    115       zwt  (1,:, : ) = 0.e0     ;     zwt  (jpi,:,:) = 0.e0 
    116       zavsi(1,:, : ) = 0.e0     ;     zavsi(jpi,:,:) = 0.e0 
    117       zwt  (:,:,jpk) = 0.e0     ;     zwt  ( : ,:,1) = 0.e0 
    118       zavsi(:,:,jpk) = 0.e0     ;     zavsi( : ,:,1) = 0.e0 
     117      zwd(1,:, : ) = 0.e0     ;     zwd(jpi,:,:) = 0.e0 
     118      zws(1,:, : ) = 0.e0     ;     zws(jpi,:,:) = 0.e0 
     119      zwi(1,:, : ) = 0.e0     ;     zwi(jpi,:,:) = 0.e0 
     120      zwt(1,:, : ) = 0.e0     ;     zwt(jpi,:,:) = 0.e0 
     121      zwt(:,:,jpk) = 0.e0     ;     zwt( : ,:,1) = 0.e0 
    119122 
    120123      ! I.1 Variable volume : to take into account vertical variable vertical scale factors 
     
    130133      !     dk[ avt dk[ (t,s) ] ] diffusive trends 
    131134 
    132  
     135      ! 
    133136      ! II.0 Matrix construction 
    134137      ! ------------------------ 
    135  
     138      DO jn = 1, kjpt 
     139         ! 
     140         !  Matrix construction 
     141         ! ------------------------ 
     142         IF( cdtype == 'TRA' .AND. jn == jp_tem )  THEN  
    136143#if defined key_ldfslp 
    137       ! update and save of avt (and avs if double diffusive mixing) 
    138       IF( l_traldf_rot ) THEN 
    139          DO jk = 2, jpkm1 
     144            ! update and save of avt (and avs if double diffusive mixing) 
     145            IF( l_traldf_rot ) THEN 
     146               DO jk = 2, jpkm1 
     147                  DO jj = 2, jpjm1 
     148                     DO ji = fs_2, fs_jpim1   ! vector opt. 
     149                        zavi = fsahtw(ji,jj,jk)                       &   ! vertical mixing coef. due to lateral mixing 
     150                           & * (  wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   & 
     151                           &    + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)  ) 
     152                        zwt(ji,jj,jk) = avt(ji,jj,jk) + zavi              ! zwt=avt+zavi (total vertical mixing coef. on temperature) 
     153                     END DO 
     154                  END DO 
     155               END DO 
     156            ELSE                         ! no rotation but key_ldfslp defined 
     157               zwt  (:,:,:) = avt(:,:,:) 
     158            ENDIF 
     159#else 
     160            ! No isopycnal diffusion 
     161            zwt(:,:,:) = avt(:,:,:)            
     162#endif 
     163            ! Diagonal, inferior, superior  (including the bottom boundary condition via avt masked) 
     164            DO jk = 1, jpkm1 
     165               DO jj = 2, jpjm1 
     166                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     167                     ze3ta =  ( 1. - znvvl ) +        znvvl   * fse3t_a(ji,jj,jk)   ! after scale factor at T-point 
     168                     ze3tn =         znvvl   + ( 1. - znvvl ) * fse3t_n(ji,jj,jk)   ! now   scale factor at T-point 
     169                     zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk  ) / ( ze3tn * fse3w(ji,jj,jk  ) ) 
     170                     zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) ) 
     171                     zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk) 
     172                 END DO 
     173               END DO 
     174            END DO 
     175            ! Surface boudary conditions 
    140176            DO jj = 2, jpjm1 
    141177               DO ji = fs_2, fs_jpim1   ! vector opt. 
    142                   zavi = fsahtw(ji,jj,jk)                       &   ! vertical mixing coef. due to lateral mixing 
    143                      & * (  wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   & 
    144                      &    + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)  ) 
    145                   zwt(ji,jj,jk) = avt(ji,jj,jk) + zavi              ! zwt=avt+zavi (total vertical mixing coef. on temperature) 
    146 # if defined key_zdfddm 
    147                   zavsi(ji,jj,jk) = fsavs(ji,jj,jk) + zavi          ! dd mixing: zavsi = total vertical mixing coef. on salinity 
    148 # endif 
    149                END DO 
    150             END DO 
    151          END DO 
    152       ELSE                         ! no rotation but key_ldfslp defined 
    153          zwt  (:,:,:) = avt(:,:,:) 
    154 # if defined key_zdfddm 
    155          zavsi(:,:,:) = avs(:,:,:)       ! avs /= avt (double diffusion mixing) 
    156 # endif 
    157       ENDIF    
     178                 ze3ta = ( 1. - znvvl ) +  znvvl * fse3t_a(ji,jj,1)    ! after scale factor at T-point 
     179                 zwi(ji,jj,1) = 0.e0 
     180                 zwd(ji,jj,1) = ze3ta - zws(ji,jj,1) 
     181               END DO 
     182            END DO 
     183            ! 
     184         ELSE IF( ( cdtype == 'TRA' .AND. jn == jp_sal ) .OR. ( cdtype == 'TRC' .AND. jn == 1 ) ) THEN 
     185#if defined key_ldfslp 
     186            ! update and save of avt (and avs if double diffusive mixing) 
     187            IF( l_traldf_rot ) THEN 
     188               DO jk = 2, jpkm1 
     189                  DO jj = 2, jpjm1 
     190                     DO ji = fs_2, fs_jpim1   ! vector opt. 
     191                        zavi = fsahtw(ji,jj,jk)                       &   ! vertical mixing coef. due to lateral mixing 
     192                           & * (  wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   & 
     193                           &    + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)  ) 
     194                        zwt(ji,jj,jk) = fsavs(ji,jj,jk) + zavi              ! zwt=avt+zavi (total vertical mixing coef. on temperature) 
     195                     END DO 
     196                  END DO 
     197               END DO 
     198            ELSE                         ! no rotation but key_ldfslp defined 
     199               zwt  (:,:,:) = fsavs(:,:,:) 
     200            ENDIF 
    158201#else 
    159       ! No isopycnal diffusion 
    160       zwt(:,:,:) = avt(:,:,:) 
    161 # if defined key_zdfddm 
    162       zavsi(:,:,:) = avs(:,:,:) 
    163 # endif 
    164  
     202            ! No isopycnal diffusion 
     203            zwt(:,:,:) = fsavs(:,:,:)            
    165204#endif 
    166  
    167       ! Diagonal, inferior, superior  (including the bottom boundary condition via avt masked) 
    168       DO jk = 1, jpkm1 
    169          DO jj = 2, jpjm1 
    170             DO ji = fs_2, fs_jpim1   ! vector opt. 
    171                ze3ta =  ( 1. - znvvl )   &                            ! after scale factor at T-point 
    172                   &   +        znvvl    * fse3t_a(ji,jj,jk)  
    173                ze3tn =    znvvl          &                            ! now   scale factor at T-point 
    174                   &   + ( 1. - znvvl ) * fse3t_n(ji,jj,jk) 
    175                zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk  ) / ( ze3tn * fse3w(ji,jj,jk  ) ) 
    176                zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) ) 
    177                zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk) 
    178             END DO 
    179          END DO 
    180       END DO 
    181  
    182       ! Surface boudary conditions 
    183       DO jj = 2, jpjm1 
    184          DO ji = fs_2, fs_jpim1   ! vector opt. 
    185             ze3ta = ( 1. - znvvl )        &                           ! after scale factor at T-point 
    186                &   +       znvvl      * fse3t_a(ji,jj,1)  
    187             zwi(ji,jj,1) = 0.e0 
    188             zwd(ji,jj,1) = ze3ta - zws(ji,jj,1) 
    189          END DO 
    190       END DO 
    191  
    192  
    193       ! II.1. Vertical diffusion on t 
    194       ! --------------------------- 
    195  
    196       !! Matrix inversion from the first level 
    197       !!---------------------------------------------------------------------- 
    198       !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk ) 
    199       ! 
    200       !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 ) 
    201       !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 ) 
    202       !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 ) 
    203       !        (        ...               )( ...  ) ( ...  ) 
    204       !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk ) 
    205       ! 
    206       !   m is decomposed in the product of an upper and lower triangular matrix 
    207       !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi 
    208       !   The second member is in 2d array zwy 
    209       !   The solution is in 2d array zwx 
    210       !   The 3d arry zwt is a work space array 
    211       !   zwy is used and then used as a work space array : its value is modified! 
    212  
    213       ! first recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k) 
    214       DO jj = 2, jpjm1 
    215          DO ji = fs_2, fs_jpim1 
    216             zwt(ji,jj,1) = zwd(ji,jj,1) 
    217          END DO 
    218       END DO 
    219       DO jk = 2, jpkm1 
     205            ! Diagonal, inferior, superior  (including the bottom boundary condition via avt masked) 
     206            DO jk = 1, jpkm1 
     207               DO jj = 2, jpjm1 
     208                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     209                     ze3ta =  ( 1. - znvvl ) +        znvvl   * fse3t_a(ji,jj,jk)   ! after scale factor at T-point 
     210                     ze3tn =         znvvl   + ( 1. - znvvl ) * fse3t_n(ji,jj,jk)   ! now   scale factor at T-point 
     211                     zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk  ) / ( ze3tn * fse3w(ji,jj,jk  ) ) 
     212                     zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) ) 
     213                     zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk) 
     214                 END DO 
     215               END DO 
     216            END DO 
     217            ! Surface boudary conditions 
     218            DO jj = 2, jpjm1 
     219               DO ji = fs_2, fs_jpim1   ! vector opt. 
     220                 ze3ta = ( 1. - znvvl ) +  znvvl * fse3t_a(ji,jj,1)    ! after scale factor at T-point 
     221                 zwi(ji,jj,1) = 0.e0 
     222                 zwd(ji,jj,1) = ze3ta - zws(ji,jj,1) 
     223               END DO 
     224            END DO 
     225            ! 
     226         END IF 
     227 
     228         ! II.1. Vertical diffusion on tracer 
     229         ! --------------------------- 
     230 
     231         !! Matrix inversion from the first level 
     232         !!---------------------------------------------------------------------- 
     233         !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk ) 
     234         ! 
     235         !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 ) 
     236         !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 ) 
     237         !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 ) 
     238         !        (        ...               )( ...  ) ( ...  ) 
     239         !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk ) 
     240         ! 
     241         !   m is decomposed in the product of an upper and lower triangular matrix 
     242         !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi 
     243         !   The second member is in 2d array zwy 
     244         !   The solution is in 2d array zwx 
     245         !   The 3d arry zwt is a work space array 
     246         !   zwy is used and then used as a work space array : its value is modified! 
     247          
     248         ! first recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k) 
    220249         DO jj = 2, jpjm1 
    221250            DO ji = fs_2, fs_jpim1 
    222                zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1)  /zwt(ji,jj,jk-1) 
    223             END DO 
    224          END DO 
    225       END DO 
    226  
    227       ! second recurrence:    Zk = Yk - Ik / Tk-1  Zk-1 
    228       DO jj = 2, jpjm1 
    229          DO ji = fs_2, fs_jpim1 
    230             ze3tb = ( 1. - znvvl ) + znvvl * fse3t_b(ji,jj,1) 
    231             ze3tn = ( 1. - znvvl ) + znvvl * fse3t(ji,jj,1) 
    232             ta(ji,jj,1) = ze3tb * tb(ji,jj,1) + p2dt(1) * ze3tn * ta(ji,jj,1) 
    233          END DO 
    234       END DO 
    235       DO jk = 2, jpkm1 
     251               zwt(ji,jj,1) = zwd(ji,jj,1) 
     252            END DO 
     253         END DO 
     254         DO jk = 2, jpkm1 
     255            DO jj = 2, jpjm1 
     256               DO ji = fs_2, fs_jpim1 
     257                  zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1)  /zwt(ji,jj,jk-1) 
     258               END DO 
     259            END DO 
     260         END DO 
     261          
     262         ! second recurrence:    Zk = Yk - Ik / Tk-1  Zk-1 
    236263         DO jj = 2, jpjm1 
    237264            DO ji = fs_2, fs_jpim1 
    238                ze3tb = ( 1. - znvvl ) + znvvl * fse3t_b(ji,jj,jk) 
    239                ze3tn = ( 1. - znvvl ) + znvvl * fse3t  (ji,jj,jk) 
    240                zrhs = ze3tb * tb(ji,jj,jk) + p2dt(jk) * ze3tn * ta(ji,jj,jk)   ! zrhs=right hand side  
    241                ta(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * ta(ji,jj,jk-1) 
    242             END DO 
    243          END DO 
    244       END DO 
    245  
    246       ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk 
    247       ! Save the masked temperature after in ta 
    248       ! (c a u t i o n:  temperature not its trend, Leap-frog scheme done it will not be done in tranxt) 
    249       DO jj = 2, jpjm1 
    250          DO ji = fs_2, fs_jpim1 
    251             ta(ji,jj,jpkm1) = ta(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1) 
    252          END DO 
    253       END DO 
    254       DO jk = jpk-2, 1, -1 
     265               ze3tb = ( 1. - znvvl ) + znvvl * fse3t_b(ji,jj,1) 
     266               ze3tn = ( 1. - znvvl ) + znvvl * fse3t(ji,jj,1) 
     267               ptraa(ji,jj,1,jn) = ze3tb * ptrab(ji,jj,1,jn) + p2dt(1) * ze3tn * ptraa(ji,jj,1,jn) 
     268            END DO 
     269         END DO 
     270         DO jk = 2, jpkm1 
     271            DO jj = 2, jpjm1 
     272               DO ji = fs_2, fs_jpim1 
     273                  ze3tb = ( 1. - znvvl ) + znvvl * fse3t_b(ji,jj,jk) 
     274                  ze3tn = ( 1. - znvvl ) + znvvl * fse3t  (ji,jj,jk) 
     275                  zrhs = ze3tb * ptrab(ji,jj,jk,jn) + p2dt(jk) * ze3tn * ptraa(ji,jj,jk,jn)   ! zrhs=right hand side  
     276                  ptraa(ji,jj,jk,jn) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * ptraa(ji,jj,jk-1,jn) 
     277               END DO 
     278            END DO 
     279         END DO 
     280          
     281         ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk 
     282         ! Save the masked temperature after in ta 
     283         ! (c a u t i o n:  temperature not its trend, Leap-frog scheme done it will not be done in tranxt) 
    255284         DO jj = 2, jpjm1 
    256285            DO ji = fs_2, fs_jpim1 
    257                ta(ji,jj,jk) = ( ta(ji,jj,jk) - zws(ji,jj,jk) * ta(ji,jj,jk+1) ) / zwt(ji,jj,jk) * tmask(ji,jj,jk) 
    258             END DO 
    259          END DO 
    260       END DO 
    261  
    262       ! II.2 Vertical diffusion on salinity 
    263       ! ----------------------------------- 
    264  
    265 #if defined key_zdfddm 
    266       ! Rebuild the Matrix as avt /= avs 
    267  
    268       ! Diagonal, inferior, superior  (including the bottom boundary condition via avs masked) 
    269       DO jk = 1, jpkm1 
    270          DO jj = 2, jpjm1 
    271             DO ji = fs_2, fs_jpim1   ! vector opt. 
    272                ze3ta =  ( 1. - znvvl )                        &         ! after scale factor at T-point 
    273                   &   +        znvvl   * fse3t_a(ji,jj,jk)            
    274                ze3tn =    znvvl                               &         ! now   scale factor at T-point 
    275                   &   + ( 1. - znvvl ) * fse3t_n(ji,jj,jk) 
    276                zwi(ji,jj,jk) = - p2dt(jk) * zavsi(ji,jj,jk  ) / ( ze3tn * fse3w(ji,jj,jk  ) ) 
    277                zws(ji,jj,jk) = - p2dt(jk) * zavsi(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) ) 
    278                zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk) 
    279             END DO 
    280          END DO 
    281       END DO 
    282  
    283       ! Surface boudary conditions 
    284       DO jj = 2, jpjm1 
    285          DO ji = fs_2, fs_jpim1   ! vector opt. 
    286             ze3ta = ( 1. - znvvl ) + znvvl * fse3t_a(ji,jj,1) 
    287             zwi(ji,jj,1) = 0.e0 
    288             zwd(ji,jj,1) = ze3ta - zws(ji,jj,1) 
    289          END DO 
    290       END DO 
    291 #endif 
    292  
    293  
    294       !! Matrix inversion from the first level 
    295       !!---------------------------------------------------------------------- 
    296       !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk ) 
    297       ! 
    298       !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 ) 
    299       !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 ) 
    300       !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 ) 
    301       !        (        ...               )( ...  ) ( ...  ) 
    302       !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk ) 
    303       ! 
    304       !   m is decomposed in the product of an upper and lower triangular 
    305       !   matrix 
    306       !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi 
    307       !   The second member is in 2d array zwy 
    308       !   The solution is in 2d array zwx 
    309       !   The 3d arry zwt is a work space array 
    310       !   zwy is used and then used as a work space array : its value is modified! 
    311  
    312       ! first recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k) 
    313       DO jj = 2, jpjm1 
    314          DO ji = fs_2, fs_jpim1 
    315             zwt(ji,jj,1) = zwd(ji,jj,1) 
    316          END DO 
    317       END DO 
    318       DO jk = 2, jpkm1 
    319          DO jj = 2, jpjm1 
    320             DO ji = fs_2, fs_jpim1 
    321                zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1)  /zwt(ji,jj,jk-1) 
    322             END DO 
    323          END DO 
    324       END DO 
    325  
    326       ! second recurrence:    Zk = Yk - Ik / Tk-1  Zk-1 
    327       DO jj = 2, jpjm1 
    328          DO ji = fs_2, fs_jpim1 
    329             ze3tb = ( 1. - znvvl )   &                                           ! before scale factor at T-point 
    330                &   +  znvvl       * fse3t_b(ji,jj,1) 
    331             ze3tn = ( 1. - znvvl ) + znvvl * fse3t  (ji,jj,1)                    ! now    scale factor at T-point 
    332             sa(ji,jj,1) = ze3tb * sb(ji,jj,1) + p2dt(1) * ze3tn * sa(ji,jj,1) 
    333          END DO 
    334       END DO 
    335       DO jk = 2, jpkm1 
    336          DO jj = 2, jpjm1 
    337             DO ji = fs_2, fs_jpim1 
    338                ze3tb = ( 1. - znvvl )   &                                        ! before scale factor at T-point 
    339                   &   +  znvvl       * fse3t_b(ji,jj,jk) 
    340                ze3tn = ( 1. - znvvl ) + znvvl * fse3t  (ji,jj,jk)                ! now    scale factor at T-point 
    341                zrhs = ze3tb * sb(ji,jj,jk) + p2dt(jk) * ze3tn * sa(ji,jj,jk)     ! zrhs=right hand side 
    342                sa(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *sa(ji,jj,jk-1) 
    343             END DO 
    344          END DO 
    345       END DO 
    346  
    347       ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk 
    348       ! Save the masked temperature after in ta 
    349       ! (c a u t i o n:  temperature not its trend, Leap-frog scheme done it will not be done in tranxt) 
    350       DO jj = 2, jpjm1 
    351          DO ji = fs_2, fs_jpim1 
    352             sa(ji,jj,jpkm1) = sa(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1) 
    353          END DO 
    354       END DO 
    355       DO jk = jpk-2, 1, -1 
    356          DO jj = 2, jpjm1 
    357             DO ji = fs_2, fs_jpim1 
    358                sa(ji,jj,jk) = ( sa(ji,jj,jk) - zws(ji,jj,jk) * sa(ji,jj,jk+1) ) / zwt(ji,jj,jk) * tmask(ji,jj,jk) 
    359             END DO 
    360          END DO 
     286               ptraa(ji,jj,jpkm1,jn) = ptraa(ji,jj,jpkm1,jn) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1) 
     287            END DO 
     288         END DO 
     289         DO jk = jpk-2, 1, -1 
     290            DO jj = 2, jpjm1 
     291               DO ji = fs_2, fs_jpim1 
     292                  ptraa(ji,jj,jk,jn) = ( ptraa(ji,jj,jk,jn) - zws(ji,jj,jk) * ptraa(ji,jj,jk+1,jn) ) & 
     293                  &                    / zwt(ji,jj,jk) * tmask(ji,jj,jk) 
     294               END DO 
     295            END DO 
     296         END DO 
     297         ! 
    361298      END DO 
    362299      ! 
Note: See TracChangeset for help on using the changeset viewer.