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 6808 for branches/NERC/dev_r5549_BDY_ZEROGRAD/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90 – NEMO

Ignore:
Timestamp:
2016-07-19T10:38:35+02:00 (8 years ago)
Author:
jamesharle
Message:

merge with trunk@6232 for consistency with SSB code

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/NERC/dev_r5549_BDY_ZEROGRAD/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90

    r5120 r6808  
    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 
     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   !   -      - 
    7274      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwi, zwd, zws 
    7375      !!---------------------------------------------------------------------- 
     
    8284         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ' 
    8385         ! 
    84          IF( lk_vvl ) THEN   ;    r_vvl = 1._wp       ! Variable volume indicator 
    85          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 
    8688         ENDIF 
    8789      ENDIF 
    88  
    89       ! 0. Local constant initialization 
    90       ! -------------------------------- 
    91       z1_p2dt = 1._wp / p2dt      ! inverse of the timestep 
    92  
    93       ! 1. Apply semi-implicit bottom friction 
    94       ! -------------------------------------- 
     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      ! 
    95109      ! Only needed for semi-implicit bottom friction setup. The explicit 
    96110      ! bottom friction has been included in "u(v)a" which act as the R.H.S 
    97111      ! column vector of the tri-diagonal matrix equation 
    98112      ! 
    99  
    100113      IF( ln_bfrimp ) THEN 
    101114         DO jj = 2, jpjm1 
     
    103116               ikbu = mbku(ji,jj)       ! ocean bottom level at u- and v-points  
    104117               ikbv = mbkv(ji,jj)       ! (deepest ocean u- and v-points) 
    105                avmu(ji,jj,ikbu+1) = -bfrua(ji,jj) * fse3uw(ji,jj,ikbu+1) 
    106                avmv(ji,jj,ikbv+1) = -bfrva(ji,jj) * fse3vw(ji,jj,ikbv+1) 
     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) 
    107120            END DO 
    108121         END DO 
     
    112125                  ikbu = miku(ji,jj)       ! ocean top level at u- and v-points  
    113126                  ikbv = mikv(ji,jj)       ! (first wet ocean u- and v-points) 
    114                   IF (ikbu .GE. 2) avmu(ji,jj,ikbu) = -tfrua(ji,jj) * fse3uw(ji,jj,ikbu) 
    115                   IF (ikbv .GE. 2) avmv(ji,jj,ikbv) = -tfrva(ji,jj) * fse3vw(ji,jj,ikbv) 
     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) 
    116129               END DO 
    117130            END DO 
    118131         END IF 
    119132      ENDIF 
    120  
    121 #if defined key_dynspg_ts 
    122       IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN      ! applied on velocity 
    123          DO jk = 1, jpkm1 
    124             ua(:,:,jk) = ( ub(:,:,jk) + p2dt * ua(:,:,jk) ) * umask(:,:,jk) 
    125             va(:,:,jk) = ( vb(:,:,jk) + p2dt * va(:,:,jk) ) * vmask(:,:,jk) 
    126          END DO 
    127       ELSE                                            ! applied on thickness weighted velocity 
    128          DO jk = 1, jpkm1 
    129             ua(:,:,jk) = (          ub(:,:,jk) * fse3u_b(:,:,jk)      & 
    130                &           + p2dt * ua(:,:,jk) * fse3u_n(:,:,jk)  )   & 
    131                &                               / fse3u_a(:,:,jk) * umask(:,:,jk) 
    132             va(:,:,jk) = (          vb(:,:,jk) * fse3v_b(:,:,jk)      & 
    133                &           + p2dt * va(:,:,jk) * fse3v_n(:,:,jk)  )   & 
    134                &                               / fse3v_a(:,:,jk) * vmask(:,:,jk) 
    135          END DO 
    136       ENDIF 
    137  
    138       IF ( ln_bfrimp .AND.lk_dynspg_ts ) THEN 
    139          ! remove barotropic velocities: 
    140          DO jk = 1, jpkm1 
    141             ua(:,:,jk) = (ua(:,:,jk) - ua_b(:,:)) * umask(:,:,jk) 
    142             va(:,:,jk) = (va(:,:,jk) - va_b(:,:)) * vmask(:,:,jk) 
    143          END DO 
    144          ! Add bottom/top stress due to barotropic component only: 
    145          DO jj = 2, jpjm1         
     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 
    146145            DO ji = fs_2, fs_jpim1   ! vector opt. 
    147146               ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points  
    148147               ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points) 
    149                ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,ikbu) + r_vvl   * fse3u_a(ji,jj,ikbu) 
    150                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) 
    151150               ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * bfrua(ji,jj) * ua_b(ji,jj) / ze3ua 
    152151               va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * bfrva(ji,jj) * va_b(ji,jj) / ze3va 
    153152            END DO 
    154153         END DO 
    155          IF ( ln_isfcav ) THEN 
     154         IF( ln_isfcav ) THEN    ! Ocean cavities (ISF) 
    156155            DO jj = 2, jpjm1         
    157156               DO ji = fs_2, fs_jpim1   ! vector opt. 
    158157                  ikbu = miku(ji,jj)         ! top ocean level at u- and v-points  
    159158                  ikbv = mikv(ji,jj)         ! (first wet ocean u- and v-points) 
    160                   ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,ikbu) + r_vvl   * fse3u_a(ji,jj,ikbu) 
    161                   ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,ikbv) + r_vvl   * fse3v_a(ji,jj,ikbv) 
     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) 
    162161                  ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * tfrua(ji,jj) * ua_b(ji,jj) / ze3ua 
    163162                  va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * tfrva(ji,jj) * va_b(ji,jj) / ze3va 
     
    166165         END IF 
    167166      ENDIF 
    168 #endif 
    169  
    170       ! 2. Vertical diffusion on u 
    171       ! --------------------------- 
     167      ! 
     168      !              !==  Vertical diffusion on u  ==! 
     169      ! 
    172170      ! Matrix and second member construction 
    173171      ! bottom boundary condition: both zwi and zws must be masked as avmu can take 
     
    177175         DO jj = 2, jpjm1  
    178176            DO ji = fs_2, fs_jpim1   ! vector opt. 
    179                ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,jk) + r_vvl   * fse3u_a(ji,jj,jk)   ! after scale factor at T-point 
    180                zcoef = - p2dt / ze3ua       
    181                zzwi          = zcoef * avmu  (ji,jj,jk  ) / fse3uw(ji,jj,jk  ) 
    182                zwi(ji,jj,jk) = zzwi  * wumask(ji,jj,jk  ) 
    183                zzws          = zcoef * avmu  (ji,jj,jk+1) / fse3uw(ji,jj,jk+1)  
    184                zws(ji,jj,jk) = zzws  * wumask(ji,jj,jk+1) 
     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) 
    185182               zwd(ji,jj,jk) = 1._wp - zzwi - zzws 
    186183            END DO 
     
    209206      !----------------------------------------------------------------------- 
    210207      ! 
    211       !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  == 
    212       DO jk = 2, jpkm1 
     208      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  == 
    213209         DO jj = 2, jpjm1    
    214210            DO ji = fs_2, fs_jpim1   ! vector opt. 
     
    218214      END DO 
    219215      ! 
    220       DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  == 
    221          DO ji = fs_2, fs_jpim1   ! vector opt. 
    222 #if defined key_dynspg_ts 
    223             ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,1) + r_vvl   * fse3u_a(ji,jj,1)  
     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)  
    224219            ua(ji,jj,1) = ua(ji,jj,1) + p2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   & 
    225220               &                                      / ( ze3ua * rau0 ) * umask(ji,jj,1)  
    226 #else 
    227             ua(ji,jj,1) = ub(ji,jj,1) & 
    228                &                   + p2dt *(ua(ji,jj,1) +  0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   & 
    229                &                                      / ( fse3u(ji,jj,1) * rau0     ) * umask(ji,jj,1) )  
    230 #endif 
    231221         END DO 
    232222      END DO 
     
    234224         DO jj = 2, jpjm1 
    235225            DO ji = fs_2, fs_jpim1 
    236 #if defined key_dynspg_ts 
    237                zrhs = ua(ji,jj,jk)   ! zrhs=right hand side 
    238 #else 
    239                zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk) 
    240 #endif 
    241                ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1) 
    242             END DO 
    243          END DO 
    244       END DO 
    245       ! 
    246       DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk  == 
     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  ==! 
    247232         DO ji = fs_2, fs_jpim1   ! vector opt. 
    248233            ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 
     
    256241         END DO 
    257242      END DO 
    258  
    259 #if ! defined key_dynspg_ts 
    260       ! Normalization to obtain the general momentum trend ua 
    261       DO jk = 1, jpkm1 
    262          DO jj = 2, jpjm1    
    263             DO ji = fs_2, fs_jpim1   ! vector opt. 
    264                ua(ji,jj,jk) = ( ua(ji,jj,jk) - ub(ji,jj,jk) ) * z1_p2dt 
    265             END DO 
    266          END DO 
    267       END DO 
    268 #endif 
    269  
    270       ! 3. Vertical diffusion on v 
    271       ! --------------------------- 
     243      ! 
     244      !              !==  Vertical diffusion on v  ==! 
     245      ! 
    272246      ! Matrix and second member construction 
    273247      ! bottom boundary condition: both zwi and zws must be masked as avmv can take 
     
    277251         DO jj = 2, jpjm1    
    278252            DO ji = fs_2, fs_jpim1   ! vector opt. 
    279                ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,jk)  + r_vvl * fse3v_a(ji,jj,jk)   ! after scale factor at T-point 
    280                zcoef = - p2dt / ze3va 
    281                zzwi          = zcoef * avmv (ji,jj,jk  ) / fse3vw(ji,jj,jk  ) 
    282                zwi(ji,jj,jk) =  zzwi * wvmask(ji,jj,jk) 
    283                zzws          = zcoef * avmv (ji,jj,jk+1) / fse3vw(ji,jj,jk+1) 
    284                zws(ji,jj,jk) =  zzws * wvmask(ji,jj,jk+1) 
     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) 
    285258               zwd(ji,jj,jk) = 1._wp - zzwi - zzws 
    286259            END DO 
     
    309282      !----------------------------------------------------------------------- 
    310283      ! 
    311       !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  == 
    312       DO jk = 2, jpkm1         
     284      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  == 
    313285         DO jj = 2, jpjm1    
    314286            DO ji = fs_2, fs_jpim1   ! vector opt. 
     
    318290      END DO 
    319291      ! 
    320       DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  == 
    321          DO ji = fs_2, fs_jpim1   ! vector opt. 
    322 #if defined key_dynspg_ts             
    323             ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,1) + r_vvl   * fse3v_a(ji,jj,1)  
     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)  
    324295            va(ji,jj,1) = va(ji,jj,1) + p2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   & 
    325296               &                                      / ( ze3va * rau0 )  
    326 #else 
    327             va(ji,jj,1) = vb(ji,jj,1) & 
    328                &                   + p2dt *(va(ji,jj,1) +  0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   & 
    329                &                                                       / ( fse3v(ji,jj,1) * rau0     )  ) 
    330 #endif 
    331297         END DO 
    332298      END DO 
     
    334300         DO jj = 2, jpjm1 
    335301            DO ji = fs_2, fs_jpim1   ! vector opt. 
    336 #if defined key_dynspg_ts 
    337                zrhs = va(ji,jj,jk)   ! zrhs=right hand side 
    338 #else 
    339                zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk) 
    340 #endif 
    341                va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1) 
    342             END DO 
    343          END DO 
    344       END DO 
    345       ! 
    346       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  ==! 
    347308         DO ji = fs_2, fs_jpim1   ! vector opt. 
    348309            va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 
     
    356317         END DO 
    357318      END DO 
    358  
    359       ! Normalization to obtain the general momentum trend va 
    360 #if ! defined key_dynspg_ts 
    361       DO jk = 1, jpkm1 
    362          DO jj = 2, jpjm1    
    363             DO ji = fs_2, fs_jpim1   ! vector opt. 
    364                va(ji,jj,jk) = ( va(ji,jj,jk) - vb(ji,jj,jk) ) * z1_p2dt 
    365             END DO 
    366          END DO 
    367       END DO 
    368 #endif 
    369  
     319       
    370320      ! J. Chanut: Lines below are useless ? 
    371321      !! restore bottom layer avmu(v)  
     322      !!gm  I almost sure it is !!!! 
    372323      IF( ln_bfrimp ) THEN 
    373324        DO jj = 2, jpjm1 
     
    375326              ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points  
    376327              ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points) 
    377               avmu(ji,jj,ikbu+1) = 0.e0 
    378               avmv(ji,jj,ikbv+1) = 0.e0 
     328              avmu(ji,jj,ikbu+1) = 0._wp 
     329              avmv(ji,jj,ikbv+1) = 0._wp 
    379330           END DO 
    380331        END DO 
     
    384335                 ikbu = miku(ji,jj)         ! ocean top level at u- and v-points  
    385336                 ikbv = mikv(ji,jj)         ! (first wet ocean u- and v-points) 
    386                  IF (ikbu > 1) avmu(ji,jj,ikbu) = 0.e0 
    387                  IF (ikbv > 1) avmv(ji,jj,ikbv) = 0.e0 
     337                 IF( ikbu > 1 )   avmu(ji,jj,ikbu) = 0._wp 
     338                 IF( ikbv > 1 )   avmv(ji,jj,ikbv) = 0._wp 
    388339              END DO 
    389340           END DO 
    390         END IF 
    391       ENDIF 
    392       ! 
    393       CALL wrk_dealloc( jpi,jpj,jpk, zwi, zwd, zws)  
    394       ! 
    395       IF( nn_timing == 1 )  CALL timing_stop('dyn_zdf_imp') 
     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') 
    396347      ! 
    397348   END SUBROUTINE dyn_zdf_imp 
Note: See TracChangeset for help on using the changeset viewer.