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 6225 for branches/2014/dev_r4704_NOC5_MPP_BDY_UPDATE/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90 – NEMO

Ignore:
Timestamp:
2016-01-08T10:35:19+01:00 (8 years ago)
Author:
jamesharle
Message:

Update MPP_BDY_UPDATE branch to be consistent with head of trunk

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/2014/dev_r4704_NOC5_MPP_BDY_UPDATE/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90

    r4370 r6225  
    22   !!====================================================================== 
    33   !!                    ***  MODULE  dynzdf_imp  *** 
    4    !! Ocean dynamics:  vertical component(s) of the momentum mixing trend 
     4   !! Ocean dynamics:  vertical component(s) of the momentum mixing trend, implicit scheme 
    55   !!====================================================================== 
    66   !! History :  OPA  !  1990-10  (B. Blanke)  Original code 
    … …  
    1212 
    1313   !!---------------------------------------------------------------------- 
    14    !!   dyn_zdf_imp  : update the momentum trend with the vertical diffusion using a implicit time-stepping 
    15    !!---------------------------------------------------------------------- 
    16    USE oce             ! ocean dynamics and tracers 
    17    USE dom_oce         ! ocean space and time domain 
    18    USE domvvl          ! variable volume 
    19    USE sbc_oce         ! surface boundary condition: ocean 
    20    USE zdf_oce         ! ocean vertical physics 
    21    USE phycst          ! physical constants 
    22    USE in_out_manager  ! I/O manager 
    23    USE lib_mpp         ! MPP library 
    24    USE zdfbfr          ! Bottom friction setup 
    25    USE wrk_nemo        ! Memory Allocation 
    26    USE timing          ! Timing 
    27    USE dynadv          ! dynamics: vector invariant versus flux form 
    28    USE dynspg_oce, ONLY: lk_dynspg_ts 
     14   !!   dyn_zdf_imp   : compute the vertical diffusion using a implicit scheme 
     15   !!                   together with the Leap-Frog time integration. 
     16   !!---------------------------------------------------------------------- 
     17   USE oce            ! ocean dynamics and tracers 
     18   USE phycst         ! physical constants 
     19   USE dom_oce        ! ocean space and time domain 
     20   USE domvvl         ! variable volume 
     21   USE sbc_oce        ! surface boundary condition: ocean 
     22   USE dynadv   , ONLY: ln_dynadv_vec ! Momentum advection form 
     23   USE zdf_oce        ! ocean vertical physics 
     24   USE zdfbfr         ! Bottom friction setup 
     25   ! 
     26   USE in_out_manager ! I/O manager 
     27   USE lib_mpp        ! MPP library 
     28   USE wrk_nemo       ! Memory Allocation 
     29   USE timing         ! Timing 
    2930 
    3031   IMPLICIT NONE 
    … …  
    3334   PUBLIC   dyn_zdf_imp   ! called by step.F90 
    3435 
    35    REAL(wp) ::  r_vvl     ! variable volume indicator, =1 if lk_vvl=T, =0 otherwise  
     36   REAL(wp) ::  r_vvl     ! non-linear free surface indicator: =0 if ln_linssh=T, =1 otherwise  
    3637 
    3738   !! * Substitutions 
    38 #  include "domzgr_substitute.h90" 
    3939#  include "vectopt_loop_substitute.h90" 
    4040   !!---------------------------------------------------------------------- 
    … …  
    5050      !!                    
    5151      !! ** Purpose :   Compute the trend due to the vert. momentum diffusion 
    52       !!      and the surface forcing, and add it to the general trend of  
    53       !!      the momentum equations. 
     52      !!              together with the Leap-Frog time stepping using an  
     53      !!              implicit scheme. 
    5454      !! 
    55       !! ** Method  :   The vertical momentum mixing trend is given by : 
    56       !!             dz( avmu dz(u) ) = 1/e3u dk+1( avmu/e3uw dk(ua) ) 
    57       !!      backward time stepping 
    58       !!      Surface boundary conditions: wind stress input (averaged over kt-1/2 & kt+1/2) 
    59       !!      Bottom boundary conditions : bottom stress (cf zdfbfr.F) 
    60       !!      Add this trend to the general trend ua : 
    61       !!         ua = ua + dz( avmu dz(u) ) 
     55      !! ** Method  :  - Leap-Frog time stepping on all trends but the vertical mixing 
     56      !!         ua =         ub + 2*dt *       ua             vector form or linear free surf. 
     57      !!         ua = ( e3u_b*ub + 2*dt * e3u_n*ua ) / e3u_a   otherwise 
     58      !!               - update the after velocity with the implicit vertical mixing. 
     59      !!      This requires to solver the following system:  
     60      !!         ua = ua + 1/e3u_a dk+1[ avmu / e3uw_a dk[ua] ] 
     61      !!      with the following surface/top/bottom boundary condition: 
     62      !!      surface: wind stress input (averaged over kt-1/2 & kt+1/2) 
     63      !!      top & bottom : top stress (iceshelf-ocean) & bottom stress (cf zdfbfr.F) 
    6264      !! 
    63       !! ** Action : - Update (ua,va) arrays with the after vertical diffusive mixing trend. 
     65      !! ** Action :   (ua,va) after velocity  
    6466      !!--------------------------------------------------------------------- 
    6567      INTEGER , INTENT(in) ::  kt     ! ocean time-step index 
    6668      REAL(wp), INTENT(in) ::  p2dt   ! vertical profile of tracer time-step 
    67       !! 
    68       INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    69       INTEGER  ::   ikbu, ikbv   ! local integers 
    70       REAL(wp) ::   z1_p2dt, zcoef, zzwi, zzws, zrhs   ! local scalars 
    71       REAL(wp) ::   ze3ua, ze3va 
    72       !!---------------------------------------------------------------------- 
    73  
     69      ! 
     70      INTEGER  ::   ji, jj, jk    ! dummy loop indices 
     71      INTEGER  ::   ikbu, ikbv    ! local integers 
     72      REAL(wp) ::   zzwi, ze3ua   ! local scalars 
     73      REAL(wp) ::   zzws, ze3va   !   -      - 
    7474      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwi, zwd, zws 
    7575      !!---------------------------------------------------------------------- 
    … …  
    8484         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ' 
    8585         ! 
    86          IF( lk_vvl ) THEN   ;    r_vvl = 1._wp       ! Variable volume indicator 
    87          ELSE                ;    r_vvl = 0._wp        
     86         If( ln_linssh ) THEN   ;    r_vvl = 0._wp    ! non-linear free surface indicator 
     87         ELSE                   ;    r_vvl = 1._wp 
    8888         ENDIF 
    8989      ENDIF 
    90  
    91       ! 0. Local constant initialization 
    92       ! -------------------------------- 
    93       z1_p2dt = 1._wp / p2dt      ! inverse of the timestep 
    94  
    95       ! 1. Apply semi-implicit bottom friction 
    96       ! -------------------------------------- 
     90      ! 
     91      !              !==  Time step dynamics  ==! 
     92      ! 
     93      IF( ln_dynadv_vec .OR. ln_linssh ) THEN      ! applied on velocity 
     94         DO jk = 1, jpkm1 
     95            ua(:,:,jk) = ( ub(:,:,jk) + p2dt * ua(:,:,jk) ) * umask(:,:,jk) 
     96            va(:,:,jk) = ( vb(:,:,jk) + p2dt * va(:,:,jk) ) * vmask(:,:,jk) 
     97         END DO 
     98      ELSE                                         ! applied on thickness weighted velocity 
     99         DO jk = 1, jpkm1 
     100            ua(:,:,jk) = (         e3u_b(:,:,jk) * ub(:,:,jk)  & 
     101               &          + p2dt * e3u_n(:,:,jk) * ua(:,:,jk)  ) / e3u_a(:,:,jk) * umask(:,:,jk) 
     102            va(:,:,jk) = (         e3v_b(:,:,jk) * vb(:,:,jk)  & 
     103               &          + p2dt * e3v_n(:,:,jk) * va(:,:,jk)  ) / e3v_a(:,:,jk) * vmask(:,:,jk) 
     104         END DO 
     105      ENDIF 
     106      ! 
     107      !              !==  Apply semi-implicit bottom friction  ==! 
     108      ! 
    97109      ! Only needed for semi-implicit bottom friction setup. The explicit 
    98110      ! bottom friction has been included in "u(v)a" which act as the R.H.S 
    99111      ! column vector of the tri-diagonal matrix equation 
    100112      ! 
    101  
    102113      IF( ln_bfrimp ) THEN 
    103 # if defined key_vectopt_loop 
    104          DO jj = 1, 1 
    105             DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling) 
    106 # else 
    107114         DO jj = 2, jpjm1 
    108115            DO ji = 2, jpim1 
    109 # endif 
    110116               ikbu = mbku(ji,jj)       ! ocean bottom level at u- and v-points  
    111117               ikbv = mbkv(ji,jj)       ! (deepest ocean u- and v-points) 
    112                avmu(ji,jj,ikbu+1) = -bfrua(ji,jj) * fse3uw(ji,jj,ikbu+1) 
    113                avmv(ji,jj,ikbv+1) = -bfrva(ji,jj) * fse3vw(ji,jj,ikbv+1) 
    114             END DO 
    115          END DO 
    116       ENDIF 
    117  
    118 #if defined key_dynspg_ts 
    119       IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN      ! applied on velocity 
    120          DO jk = 1, jpkm1 
    121             ua(:,:,jk) = ( ub(:,:,jk) + p2dt * ua(:,:,jk) ) * umask(:,:,jk) 
    122             va(:,:,jk) = ( vb(:,:,jk) + p2dt * va(:,:,jk) ) * vmask(:,:,jk) 
    123          END DO 
    124       ELSE                                            ! applied on thickness weighted velocity 
    125          DO jk = 1, jpkm1 
    126             ua(:,:,jk) = (          ub(:,:,jk) * fse3u_b(:,:,jk)      & 
    127                &           + p2dt * ua(:,:,jk) * fse3u_n(:,:,jk)  )   & 
    128                &                               / fse3u_a(:,:,jk) * umask(:,:,jk) 
    129             va(:,:,jk) = (          vb(:,:,jk) * fse3v_b(:,:,jk)      & 
    130                &           + p2dt * va(:,:,jk) * fse3v_n(:,:,jk)  )   & 
    131                &                               / fse3v_a(:,:,jk) * vmask(:,:,jk) 
    132          END DO 
    133       ENDIF 
    134  
    135       IF ( ln_bfrimp .AND.lk_dynspg_ts ) THEN 
    136          ! remove barotropic velocities: 
    137          DO jk = 1, jpkm1 
    138             ua(:,:,jk) = (ua(:,:,jk) - ua_b(:,:)) * umask(:,:,jk) 
    139             va(:,:,jk) = (va(:,:,jk) - va_b(:,:)) * vmask(:,:,jk) 
    140          ENDDO 
    141          ! Add bottom stress due to barotropic component only: 
    142          DO jj = 2, jpjm1         
     118               avmu(ji,jj,ikbu+1) = -bfrua(ji,jj) * e3uw_n(ji,jj,ikbu+1) 
     119               avmv(ji,jj,ikbv+1) = -bfrva(ji,jj) * e3vw_n(ji,jj,ikbv+1) 
     120            END DO 
     121         END DO 
     122         IF ( ln_isfcav ) THEN 
     123            DO jj = 2, jpjm1 
     124               DO ji = 2, jpim1 
     125                  ikbu = miku(ji,jj)       ! ocean top level at u- and v-points  
     126                  ikbv = mikv(ji,jj)       ! (first wet ocean u- and v-points) 
     127                  IF( ikbu >= 2 )   avmu(ji,jj,ikbu) = -tfrua(ji,jj) * e3uw_n(ji,jj,ikbu) 
     128                  IF( ikbv >= 2 )   avmv(ji,jj,ikbv) = -tfrva(ji,jj) * e3vw_n(ji,jj,ikbv) 
     129               END DO 
     130            END DO 
     131         END IF 
     132      ENDIF 
     133      ! 
     134      ! With split-explicit free surface, barotropic stress is treated explicitly 
     135      ! Update velocities at the bottom. 
     136      ! J. Chanut: The bottom stress is computed considering after barotropic velocities, which does  
     137      !            not lead to the effective stress seen over the whole barotropic loop.  
     138      ! G. Madec : in linear free surface, e3u_a = e3u_n = e3u_0, so systematic use of e3u_a 
     139      IF( ln_bfrimp .AND. ln_dynspg_ts ) THEN 
     140         DO jk = 1, jpkm1        ! remove barotropic velocities 
     141            ua(:,:,jk) = ( ua(:,:,jk) - ua_b(:,:) ) * umask(:,:,jk) 
     142            va(:,:,jk) = ( va(:,:,jk) - va_b(:,:) ) * vmask(:,:,jk) 
     143         END DO 
     144         DO jj = 2, jpjm1        ! Add bottom/top stress due to barotropic component only 
    143145            DO ji = fs_2, fs_jpim1   ! vector opt. 
    144146               ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points  
    145147               ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points) 
    146                ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,ikbu) + r_vvl   * fse3u_a(ji,jj,ikbu) 
    147                ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,ikbv) + r_vvl   * fse3v_a(ji,jj,ikbv) 
     148               ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,ikbu) + r_vvl * e3u_a(ji,jj,ikbu) 
     149               ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikbv) + r_vvl * e3v_a(ji,jj,ikbv) 
    148150               ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * bfrua(ji,jj) * ua_b(ji,jj) / ze3ua 
    149151               va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * bfrva(ji,jj) * va_b(ji,jj) / ze3va 
    150152            END DO 
    151153         END DO 
    152       ENDIF 
    153 #endif 
    154  
    155       ! 2. Vertical diffusion on u 
    156       ! --------------------------- 
     154         IF( ln_isfcav ) THEN    ! Ocean cavities (ISF) 
     155            DO jj = 2, jpjm1         
     156               DO ji = fs_2, fs_jpim1   ! vector opt. 
     157                  ikbu = miku(ji,jj)         ! top ocean level at u- and v-points  
     158                  ikbv = mikv(ji,jj)         ! (first wet ocean u- and v-points) 
     159                  ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,ikbu) + r_vvl * e3u_a(ji,jj,ikbu) 
     160                  ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikbv) + r_vvl * e3v_a(ji,jj,ikbv) 
     161                  ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * tfrua(ji,jj) * ua_b(ji,jj) / ze3ua 
     162                  va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * tfrva(ji,jj) * va_b(ji,jj) / ze3va 
     163               END DO 
     164            END DO 
     165         END IF 
     166      ENDIF 
     167      ! 
     168      !              !==  Vertical diffusion on u  ==! 
     169      ! 
    157170      ! Matrix and second member construction 
    158171      ! bottom boundary condition: both zwi and zws must be masked as avmu can take 
    … …  
    162175         DO jj = 2, jpjm1  
    163176            DO ji = fs_2, fs_jpim1   ! vector opt. 
    164                ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,jk) + r_vvl   * fse3u_a(ji,jj,jk)   ! after scale factor at T-point 
    165                zcoef = - p2dt / ze3ua       
    166                zzwi          = zcoef * avmu (ji,jj,jk  ) / fse3uw(ji,jj,jk  ) 
    167                zwi(ji,jj,jk) = zzwi  * umask(ji,jj,jk) 
    168                zzws          = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1) 
    169                zws(ji,jj,jk) = zzws  * umask(ji,jj,jk+1) 
    170                zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws 
     177               ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at T-point 
     178               zzwi = - p2dt * avmu(ji,jj,jk  ) / ( ze3ua * e3uw_n(ji,jj,jk  ) ) 
     179               zzws = - p2dt * avmu(ji,jj,jk+1) / ( ze3ua * e3uw_n(ji,jj,jk+1) ) 
     180               zwi(ji,jj,jk) = zzwi * wumask(ji,jj,jk  ) 
     181               zws(ji,jj,jk) = zzws * wumask(ji,jj,jk+1) 
     182               zwd(ji,jj,jk) = 1._wp - zzwi - zzws 
    171183            END DO 
    172184         END DO 
    … …  
    202214      END DO 
    203215      ! 
    204       DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  == 
    205          DO ji = fs_2, fs_jpim1   ! vector opt. 
    206             ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,1) + r_vvl   * fse3u_a(ji,jj,1)  
    207 #if defined key_dynspg_ts 
     216      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==! 
     217         DO ji = fs_2, fs_jpim1   ! vector opt. 
     218            ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1)  
    208219            ua(ji,jj,1) = ua(ji,jj,1) + p2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   & 
    209                &                                      / ( ze3ua * rau0 )  
    210 #else 
    211             ua(ji,jj,1) = ub(ji,jj,1) + p2dt *(ua(ji,jj,1) +  0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   & 
    212                &                                                     / ( fse3u(ji,jj,1) * rau0     ) )  
    213 #endif 
     220               &                                      / ( ze3ua * rau0 ) * umask(ji,jj,1)  
    214221         END DO 
    215222      END DO 
    216223      DO jk = 2, jpkm1 
    217          DO jj = 2, jpjm1    
    218             DO ji = fs_2, fs_jpim1   ! vector opt. 
    219 #if defined key_dynspg_ts 
    220                zrhs = ua(ji,jj,jk)   ! zrhs=right hand side 
    221 #else 
    222                zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk) 
    223 #endif 
    224                ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1) 
    225             END DO 
    226          END DO 
    227       END DO 
    228       ! 
    229       DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk  == 
     224         DO jj = 2, jpjm1 
     225            DO ji = fs_2, fs_jpim1 
     226               ua(ji,jj,jk) = ua(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1) 
     227            END DO 
     228         END DO 
     229      END DO 
     230      ! 
     231      DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk  ==! 
    230232         DO ji = fs_2, fs_jpim1   ! vector opt. 
    231233            ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 
    … …  
    233235      END DO 
    234236      DO jk = jpk-2, 1, -1 
    235          DO jj = 2, jpjm1    
    236             DO ji = fs_2, fs_jpim1   ! vector opt. 
     237         DO jj = 2, jpjm1 
     238            DO ji = fs_2, fs_jpim1 
    237239               ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk) 
    238240            END DO 
    239241         END DO 
    240242      END DO 
    241  
    242 #if ! defined key_dynspg_ts 
    243       ! Normalization to obtain the general momentum trend ua 
    244       DO jk = 1, jpkm1 
    245          DO jj = 2, jpjm1    
    246             DO ji = fs_2, fs_jpim1   ! vector opt. 
    247                ua(ji,jj,jk) = ( ua(ji,jj,jk) - ub(ji,jj,jk) ) * z1_p2dt 
    248             END DO 
    249          END DO 
    250       END DO 
    251 #endif 
    252  
    253       ! 3. Vertical diffusion on v 
    254       ! --------------------------- 
     243      ! 
     244      !              !==  Vertical diffusion on v  ==! 
     245      ! 
    255246      ! Matrix and second member construction 
    256247      ! bottom boundary condition: both zwi and zws must be masked as avmv can take 
    … …  
    260251         DO jj = 2, jpjm1    
    261252            DO ji = fs_2, fs_jpim1   ! vector opt. 
    262                ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,jk)  + r_vvl * fse3v_a(ji,jj,jk)   ! after scale factor at T-point 
    263                zcoef = - p2dt / ze3va 
    264                zzwi          = zcoef * avmv (ji,jj,jk  ) / fse3vw(ji,jj,jk  ) 
    265                zwi(ji,jj,jk) =  zzwi * vmask(ji,jj,jk) 
    266                zzws          = zcoef * avmv (ji,jj,jk+1) / fse3vw(ji,jj,jk+1) 
    267                zws(ji,jj,jk) =  zzws * vmask(ji,jj,jk+1) 
    268                zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws 
     253               ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at T-point 
     254               zzwi = - p2dt * avmv (ji,jj,jk  ) / ( ze3va * e3vw_n(ji,jj,jk  ) ) 
     255               zzws = - p2dt * avmv (ji,jj,jk+1) / ( ze3va * e3vw_n(ji,jj,jk+1) ) 
     256               zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk  ) 
     257               zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1) 
     258               zwd(ji,jj,jk) = 1._wp - zzwi - zzws 
    269259            END DO 
    270260         END DO 
    … …  
    300290      END DO 
    301291      ! 
    302       DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  == 
    303          DO ji = fs_2, fs_jpim1   ! vector opt. 
    304             ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,1) + r_vvl   * fse3v_a(ji,jj,1)  
    305 #if defined key_dynspg_ts             
     292      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==! 
     293         DO ji = fs_2, fs_jpim1   ! vector opt.           
     294            ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1)  
    306295            va(ji,jj,1) = va(ji,jj,1) + p2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   & 
    307296               &                                      / ( ze3va * rau0 )  
    308 #else 
    309             va(ji,jj,1) = vb(ji,jj,1) + p2dt *(va(ji,jj,1) +  0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   & 
    310                &                                                       / ( fse3v(ji,jj,1) * rau0     )  ) 
    311 #endif 
    312297         END DO 
    313298      END DO 
    … …  
    315300         DO jj = 2, jpjm1 
    316301            DO ji = fs_2, fs_jpim1   ! vector opt. 
    317 #if defined key_dynspg_ts 
    318                zrhs = va(ji,jj,jk)   ! zrhs=right hand side 
    319 #else 
    320                zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk) 
    321 #endif 
    322                va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1) 
    323             END DO 
    324          END DO 
    325       END DO 
    326       ! 
    327       DO jj = 2, jpjm1        !==  third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  == 
     302               va(ji,jj,jk) = va(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1) 
     303            END DO 
     304         END DO 
     305      END DO 
     306      ! 
     307      DO jj = 2, jpjm1        !==  third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  ==! 
    328308         DO ji = fs_2, fs_jpim1   ! vector opt. 
    329309            va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 
    … …  
    331311      END DO 
    332312      DO jk = jpk-2, 1, -1 
    333          DO jj = 2, jpjm1    
    334             DO ji = fs_2, fs_jpim1   ! vector opt. 
     313         DO jj = 2, jpjm1 
     314            DO ji = fs_2, fs_jpim1 
    335315               va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk) 
    336316            END DO 
    337317         END DO 
    338318      END DO 
    339  
    340       ! Normalization to obtain the general momentum trend va 
    341 #if ! defined key_dynspg_ts 
    342       DO jk = 1, jpkm1 
    343          DO jj = 2, jpjm1    
    344             DO ji = fs_2, fs_jpim1   ! vector opt. 
    345                va(ji,jj,jk) = ( va(ji,jj,jk) - vb(ji,jj,jk) ) * z1_p2dt 
    346             END DO 
    347          END DO 
    348       END DO 
    349 #endif 
    350  
     319       
    351320      ! J. Chanut: Lines below are useless ? 
    352321      !! restore bottom layer avmu(v)  
     322      !!gm  I almost sure it is !!!! 
    353323      IF( ln_bfrimp ) THEN 
    354 # if defined key_vectopt_loop 
    355       DO jj = 1, 1 
    356          DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling) 
    357 # else 
    358       DO jj = 2, jpjm1 
    359          DO ji = 2, jpim1 
    360 # endif 
    361             ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points  
    362             ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points) 
    363             avmu(ji,jj,ikbu+1) = 0.e0 
    364             avmv(ji,jj,ikbv+1) = 0.e0 
    365          END DO 
    366       END DO 
    367       ENDIF 
    368       ! 
    369       CALL wrk_dealloc( jpi,jpj,jpk, zwi, zwd, zws)  
    370       ! 
    371       IF( nn_timing == 1 )  CALL timing_stop('dyn_zdf_imp') 
     324        DO jj = 2, jpjm1 
     325           DO ji = 2, jpim1 
     326              ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points  
     327              ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points) 
     328              avmu(ji,jj,ikbu+1) = 0._wp 
     329              avmv(ji,jj,ikbv+1) = 0._wp 
     330           END DO 
     331        END DO 
     332        IF (ln_isfcav) THEN 
     333           DO jj = 2, jpjm1 
     334              DO ji = 2, jpim1 
     335                 ikbu = miku(ji,jj)         ! ocean top level at u- and v-points  
     336                 ikbv = mikv(ji,jj)         ! (first wet ocean u- and v-points) 
     337                 IF( ikbu > 1 )   avmu(ji,jj,ikbu) = 0._wp 
     338                 IF( ikbv > 1 )   avmv(ji,jj,ikbv) = 0._wp 
     339              END DO 
     340           END DO 
     341        ENDIF 
     342      ENDIF 
     343      ! 
     344      CALL wrk_dealloc( jpi,jpj,jpk,   zwi, zwd, zws)  
     345      ! 
     346      IF( nn_timing == 1 )   CALL timing_stop('dyn_zdf_imp') 
    372347      ! 
    373348   END SUBROUTINE dyn_zdf_imp 
Note: See TracChangeset for help on using the changeset viewer.