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/dynldf_iso.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/dynldf_iso.F90

    r4990 r6808  
    22   !!====================================================================== 
    33   !!                     ***  MODULE  dynldf_iso  *** 
    4    !! Ocean dynamics:  lateral viscosity trend 
     4   !! Ocean dynamics:   lateral viscosity trend (rotated laplacian operator) 
    55   !!====================================================================== 
    66   !! History :  OPA  !  97-07  (G. Madec)  Original code 
     
    88   !!             -   !  2004-08  (C. Talandier) New trends organization 
    99   !!            2.0  !  2005-11  (G. Madec)  s-coordinate: horizontal diffusion 
     10   !!            3.7  !  2014-01  (F. Lemarie, G. Madec)  restructuration/simplification of ahm specification, 
     11   !!                 !                                   add velocity dependent coefficient and optional read in file 
    1012   !!---------------------------------------------------------------------- 
    11 #if defined key_ldfslp   ||   defined key_esopa 
    12    !!---------------------------------------------------------------------- 
    13    !!   'key_ldfslp'                      slopes of the direction of mixing 
     13 
    1414   !!---------------------------------------------------------------------- 
    1515   !!   dyn_ldf_iso  : update the momentum trend with the horizontal part 
     
    1919   USE oce             ! ocean dynamics and tracers 
    2020   USE dom_oce         ! ocean space and time domain 
    21    USE ldfdyn_oce      ! ocean dynamics lateral physics 
    22    USE ldftra_oce      ! ocean tracer   lateral physics 
     21   USE ldfdyn          ! lateral diffusion: eddy viscosity coef. 
     22   USE ldftra          ! lateral physics: eddy diffusivity 
    2323   USE zdf_oce         ! ocean vertical physics 
    2424   USE ldfslp          ! iso-neutral slopes  
    2525   ! 
    26    USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    2726   USE in_out_manager  ! I/O manager 
    2827   USE lib_mpp         ! MPP library 
     28   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    2929   USE prtctl          ! Print control 
    3030   USE wrk_nemo        ! Memory Allocation 
     
    4141 
    4242   !! * Substitutions 
    43 #  include "domzgr_substitute.h90" 
    44 #  include "ldfdyn_substitute.h90" 
    4543#  include "vectopt_loop_substitute.h90" 
    4644   !!---------------------------------------------------------------------- 
     
    8381      !!      horizontal fluxes associated with the rotated lateral mixing: 
    8482      !!      u-component: 
    85       !!         ziut = ( ahmt + ahmb0 ) e2t * e3t / e1t  di[ ub ] 
    86       !!               -      ahmt       e2t * mi-1(uslp) dk[ mi(mk(ub)) ] 
    87       !!         zjuf = ( ahmf + ahmb0 ) e1f * e3f / e2f  dj[ ub ] 
    88       !!               -      ahmf       e1f * mi(vslp)   dk[ mj(mk(ub)) ] 
     83      !!         ziut = ( ahmt + rn_ahm_b ) e2t * e3t / e1t  di[ ub ] 
     84      !!               -  ahmt              e2t * mi-1(uslp) dk[ mi(mk(ub)) ] 
     85      !!         zjuf = ( ahmf + rn_ahm_b ) e1f * e3f / e2f  dj[ ub ] 
     86      !!               -  ahmf              e1f * mi(vslp)   dk[ mj(mk(ub)) ] 
    8987      !!      v-component: 
    90       !!         zivf = ( ahmf + ahmb0 ) e2t * e3t / e1t  di[ vb ] 
    91       !!               -      ahmf       e2t * mj(uslp)   dk[ mi(mk(vb)) ] 
    92       !!         zjvt = ( ahmt + ahmb0 ) e1f * e3f / e2f  dj[ ub ] 
    93       !!               -      ahmt       e1f * mj-1(vslp) dk[ mj(mk(vb)) ] 
     88      !!         zivf = ( ahmf + rn_ahm_b ) e2t * e3t / e1t  di[ vb ] 
     89      !!               -  ahmf              e2t * mj(uslp)   dk[ mi(mk(vb)) ] 
     90      !!         zjvt = ( ahmt + rn_ahm_b ) e1f * e3f / e2f  dj[ ub ] 
     91      !!               -  ahmt              e1f * mj-1(vslp) dk[ mj(mk(vb)) ] 
    9492      !!      take the horizontal divergence of the fluxes: 
    9593      !!         diffu = 1/(e1u*e2u*e3u) {  di  [ ziut ] + dj-1[ zjuf ]  } 
     
    106104      !!      of the rotated operator in dynzdf module 
    107105      !!---------------------------------------------------------------------- 
    108       ! 
    109106      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index 
    110107      ! 
    111108      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    112109      REAL(wp) ::   zabe1, zabe2, zcof1, zcof2                       ! local scalars 
    113       REAL(wp) ::   zmskt, zmskf, zbu, zbv, zuah, zvah               !   -      - 
     110      REAL(wp) ::   zmskt, zmskf                                     !   -      - 
    114111      REAL(wp) ::   zcoef0, zcoef3, zcoef4, zmkt, zmkf               !   -      - 
    115112      REAL(wp) ::   zuav, zvav, zuwslpi, zuwslpj, zvwslpi, zvwslpj   !   -      - 
     
    130127      ENDIF 
    131128 
    132       ! s-coordinate: Iso-level diffusion on tracer, but geopotential level diffusion on momentum 
     129!!gm bug is dyn_ldf_iso called before tra_ldf_iso ....   <<<<<===== TO BE CHECKED 
     130      ! s-coordinate: Iso-level diffusion on momentum but not on tracer 
    133131      IF( ln_dynldf_hor .AND. ln_traldf_iso ) THEN 
    134132         ! 
     
    136134            DO jj = 2, jpjm1 
    137135               DO ji = 2, jpim1 
    138                   uslp (ji,jj,jk) = -1./e1u(ji,jj) * ( fsdept_b(ji+1,jj,jk) - fsdept_b(ji ,jj ,jk) ) * umask(ji,jj,jk) 
    139                   vslp (ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept_b(ji,jj+1,jk) - fsdept_b(ji ,jj ,jk) ) * vmask(ji,jj,jk) 
    140                   wslpi(ji,jj,jk) = -1./e1t(ji,jj) * ( fsdepw_b(ji+1,jj,jk) - fsdepw_b(ji-1,jj,jk) ) * tmask(ji,jj,jk) * 0.5 
    141                   wslpj(ji,jj,jk) = -1./e2t(ji,jj) * ( fsdepw_b(ji,jj+1,jk) - fsdepw_b(ji,jj-1,jk) ) * tmask(ji,jj,jk) * 0.5 
     136                  uslp (ji,jj,jk) = - ( gdept_b(ji+1,jj,jk) - gdept_b(ji ,jj ,jk) ) * r1_e1u(ji,jj) * umask(ji,jj,jk) 
     137                  vslp (ji,jj,jk) = - ( gdept_b(ji,jj+1,jk) - gdept_b(ji ,jj ,jk) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk) 
     138                  wslpi(ji,jj,jk) = - ( gdepw_b(ji+1,jj,jk) - gdepw_b(ji-1,jj,jk) ) * r1_e1t(ji,jj) * tmask(ji,jj,jk) * 0.5 
     139                  wslpj(ji,jj,jk) = - ( gdepw_b(ji,jj+1,jk) - gdepw_b(ji,jj-1,jk) ) * r1_e2t(ji,jj) * tmask(ji,jj,jk) * 0.5 
    142140               END DO 
    143141            END DO 
     
    184182            DO jj = 2, jpjm1 
    185183               DO ji = fs_2, jpi   ! vector opt. 
    186                   zabe1 = (fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * MIN( fse3u(ji,jj,jk), fse3u(ji-1,jj,jk) ) / e1t(ji,jj) 
    187  
    188                   zmskt = 1./MAX(  umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)   & 
    189                      &           + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ), 1. ) 
    190  
    191                   zcof1 = - aht0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) 
     184                  zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj) * MIN( e3u_n(ji,jj,jk), e3u_n(ji-1,jj,jk) ) * r1_e1t(ji,jj) 
     185 
     186                  zmskt = 1._wp / MAX(   umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)     & 
     187                     &                 + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ) , 1._wp ) 
     188 
     189                  zcof1 = - rn_aht_0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) 
    192190    
     191                  ziut(ji,jj) = (  zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) )    & 
     192                     &           + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)      & 
     193                     &                      +zdk1u(ji,jj) + zdku (ji-1,jj) )  ) * tmask(ji,jj,jk) 
     194               END DO 
     195            END DO 
     196         ELSE                   ! other coordinate system (zco or sco) : e3t 
     197            DO jj = 2, jpjm1 
     198               DO ji = fs_2, jpi   ! vector opt. 
     199                  zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj) * e3t_n(ji,jj,jk) * r1_e1t(ji,jj) 
     200 
     201                  zmskt = 1._wp / MAX(   umask(ji-1,jj,jk  ) + umask(ji,jj,jk+1)     & 
     202                     &                 + umask(ji-1,jj,jk+1) + umask(ji,jj,jk  ) , 1._wp ) 
     203 
     204                  zcof1 = - rn_aht_0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) 
     205 
    193206                  ziut(ji,jj) = (  zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) )   & 
    194207                     &           + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)     & 
     
    196209               END DO 
    197210            END DO 
    198          ELSE                   ! other coordinate system (zco or sco) : e3t 
    199             DO jj = 2, jpjm1 
    200                DO ji = fs_2, jpi   ! vector opt. 
    201                   zabe1 = (fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj) 
    202  
    203                   zmskt = 1./MAX(  umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)   & 
    204                      &           + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ), 1. ) 
    205  
    206                   zcof1 = - aht0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) 
    207  
    208                   ziut(ji,jj) = (  zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) )   & 
    209                      &           + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)     & 
    210                      &                      +zdk1u(ji,jj) + zdku (ji-1,jj) )  ) * tmask(ji,jj,jk) 
    211                END DO 
    212             END DO 
    213211         ENDIF 
    214212 
     
    216214         DO jj = 1, jpjm1 
    217215            DO ji = 1, fs_jpim1   ! vector opt. 
    218                zabe2 = ( fsahmf(ji,jj,jk) + ahmb0 ) * e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj) 
    219  
    220                zmskf = 1./MAX(  umask(ji,jj+1,jk  )+umask(ji,jj,jk+1)   & 
    221                   &           + umask(ji,jj+1,jk+1)+umask(ji,jj,jk  ), 1. ) 
    222  
    223                zcof2 = - aht0 * e1f(ji,jj) * zmskf * 0.5  * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) ) 
     216               zabe2 = ( ahmf(ji,jj,jk) + rn_ahm_b ) * e1f(ji,jj) * e3f_n(ji,jj,jk) * r1_e2f(ji,jj) 
     217 
     218               zmskf = 1._wp / MAX(   umask(ji,jj+1,jk  )+umask(ji,jj,jk+1)     & 
     219                  &                 + umask(ji,jj+1,jk+1)+umask(ji,jj,jk  ) , 1._wp ) 
     220 
     221               zcof2 = - rn_aht_0 * e1f(ji,jj) * zmskf * 0.5  * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) ) 
    224222 
    225223               zjuf(ji,jj) = (  zabe2 * ( ub(ji,jj+1,jk) - ub(ji,jj,jk) )   & 
     
    237235         DO jj = 2, jpjm1 
    238236            DO ji = 1, fs_jpim1   ! vector opt. 
    239                zabe1 = ( fsahmf(ji,jj,jk) + ahmb0 ) * e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj) 
    240  
    241                zmskf = 1./MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)   & 
    242                   &           + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk  ), 1. ) 
    243  
    244                zcof1 = - aht0 * e2f(ji,jj) * zmskf * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) ) 
    245  
    246                zivf(ji,jj) = (  zabe1 * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) )   & 
    247                   &           + zcof1 * ( zdkv (ji,jj) + zdk1v(ji+1,jj)     & 
    248                   &                      +zdk1v(ji,jj) + zdkv (ji+1,jj) )  ) * fmask(ji,jj,jk) 
     237               zabe1 = ( ahmf(ji,jj,jk) + rn_ahm_b ) * e2f(ji,jj) * e3f_n(ji,jj,jk) * r1_e1f(ji,jj) 
     238 
     239               zmskf = 1._wp / MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)     & 
     240                  &                + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk  ) , 1._wp ) 
     241 
     242               zcof1 = - rn_aht_0 * e2f(ji,jj) * zmskf * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) ) 
     243 
     244               zivf(ji,jj) = (  zabe1 * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) )    & 
     245                  &           + zcof1 * (  zdkv (ji,jj) + zdk1v(ji+1,jj)      & 
     246                  &                      + zdk1v(ji,jj) + zdkv (ji+1,jj) )  ) * fmask(ji,jj,jk) 
    249247            END DO 
    250248         END DO 
     
    254252            DO jj = 2, jpj 
    255253               DO ji = 1, fs_jpim1   ! vector opt. 
    256                   zabe2 = (fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * MIN( fse3v(ji,jj,jk), fse3v(ji,jj-1,jk) ) / e2t(ji,jj) 
     254                  zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj) * MIN( e3v_n(ji,jj,jk), e3v_n(ji,jj-1,jk) ) * r1_e2t(ji,jj) 
     255 
     256                  zmskt = 1._wp / MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)     & 
     257                     &                + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ) , 1._wp ) 
     258 
     259                  zcof2 = - rn_aht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) ) 
     260 
     261                  zjvt(ji,jj) = (  zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) )    & 
     262                     &           + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj)      & 
     263                     &                      +zdk1v(ji,jj-1) + zdkv (ji,jj) )  ) * tmask(ji,jj,jk) 
     264               END DO 
     265            END DO 
     266         ELSE                   ! other coordinate system (zco or sco) : e3t 
     267            DO jj = 2, jpj 
     268               DO ji = 1, fs_jpim1   ! vector opt. 
     269                  zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj) * e3t_n(ji,jj,jk) * r1_e2t(ji,jj) 
    257270 
    258271                  zmskt = 1./MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   & 
    259272                     &           + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1. ) 
    260273 
    261                   zcof2 = - aht0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) ) 
     274                  zcof2 = - rn_aht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) ) 
    262275 
    263276                  zjvt(ji,jj) = (  zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) )   & 
     
    266279               END DO 
    267280            END DO 
    268          ELSE                   ! other coordinate system (zco or sco) : e3t 
    269             DO jj = 2, jpj 
    270                DO ji = 1, fs_jpim1   ! vector opt. 
    271                   zabe2 = (fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj) 
    272  
    273                   zmskt = 1./MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   & 
    274                      &           + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1. ) 
    275  
    276                   zcof2 = - aht0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) ) 
    277  
    278                   zjvt(ji,jj) = (  zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) )   & 
    279                      &           + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj)     & 
    280                      &                      +zdk1v(ji,jj-1) + zdkv (ji,jj) )  ) * tmask(ji,jj,jk) 
    281                END DO 
    282             END DO 
    283281         ENDIF 
    284282 
     
    286284         ! Second derivative (divergence) and add to the general trend 
    287285         ! ----------------------------------------------------------- 
    288  
    289286         DO jj = 2, jpjm1 
    290             DO ji = 2, jpim1          !! Question vectop possible??? !!bug 
    291                ! volume elements 
    292                zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) 
    293                zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) 
    294                ! horizontal component of isopycnal momentum diffusive trends 
    295                zuah =( ziut (ji+1,jj) - ziut (ji,jj  ) +   & 
    296                   &    zjuf (ji  ,jj) - zjuf (ji,jj-1)  ) / zbu 
    297                zvah =( zivf (ji,jj  ) - zivf (ji-1,jj) +   & 
    298                   &    zjvt (ji,jj+1) - zjvt (ji,jj  )  ) / zbv 
    299                ! add the trends to the general trends 
    300                ua (ji,jj,jk) = ua (ji,jj,jk) + zuah 
    301                va (ji,jj,jk) = va (ji,jj,jk) + zvah 
     287            DO ji = 2, jpim1          !!gm Question vectop possible??? !!bug 
     288               ua(ji,jj,jk) = ua(ji,jj,jk) + (  ziut(ji+1,jj) - ziut(ji,jj  )    & 
     289                  &                           + zjuf(ji  ,jj) - zjuf(ji,jj-1)  ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk) 
     290               va(ji,jj,jk) = va(ji,jj,jk) + (  zivf(ji,jj  ) - zivf(ji-1,jj)    & 
     291                  &                           + zjvt(ji,jj+1) - zjvt(ji,jj  )  ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk) 
    302292            END DO 
    303293         END DO 
     
    358348         DO jk = 2, jpkm1 
    359349            DO ji = 2, jpim1 
    360                zcoef0= 0.5 * aht0 * umask(ji,jj,jk) 
    361  
     350               zcoef0= 0.5 * rn_aht_0 * umask(ji,jj,jk) 
     351               ! 
    362352               zuwslpi = zcoef0 * ( wslpi(ji+1,jj,jk) + wslpi(ji,jj,jk) ) 
    363353               zuwslpj = zcoef0 * ( wslpj(ji+1,jj,jk) + wslpj(ji,jj,jk) ) 
    364  
     354               ! 
    365355               zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji+1,jj,jk-1)   & 
    366356                             + tmask(ji,jj,jk  )+tmask(ji+1,jj,jk  ), 1. ) 
    367                zmkf = 1./MAX(  fmask(ji,jj-1,jk-1)+fmask(ji,jj,jk-1)   & 
    368                              + fmask(ji,jj-1,jk  )+fmask(ji,jj,jk  ), 1. ) 
     357               zmkf = 1./MAX(  fmask(ji,jj-1,jk-1) + fmask(ji,jj,jk-1)   & 
     358                             + fmask(ji,jj-1,jk  ) + fmask(ji,jj,jk  ), 1. ) 
    369359 
    370360               zcoef3 = - e2u(ji,jj) * zmkt * zuwslpi 
     
    376366                                       +zdj1u(ji,jk  ) + zdju (ji  ,jk  ) ) 
    377367               ! update avmu (add isopycnal vertical coefficient to avmu) 
    378                ! Caution: zcoef0 include aht0, so divided by aht0 to obtain slp^2 * aht0 
    379                avmu(ji,jj,jk) = avmu(ji,jj,jk) + ( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) / aht0 
     368               ! Caution: zcoef0 include rn_aht_0, so divided by rn_aht_0 to obtain slp^2 * rn_aht_0 
     369               avmu(ji,jj,jk) = avmu(ji,jj,jk) + ( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) / rn_aht_0 
    380370            END DO 
    381371         END DO 
     
    384374         DO jk = 2, jpkm1 
    385375            DO ji = 2, jpim1 
    386                zcoef0= 0.5 * aht0 * vmask(ji,jj,jk) 
     376               zcoef0 = 0.5 * rn_aht_0 * vmask(ji,jj,jk) 
    387377 
    388378               zvwslpi = zcoef0 * ( wslpi(ji,jj+1,jk) + wslpi(ji,jj,jk) ) 
     
    398388               ! vertical flux on v field 
    399389               zfvw(ji,jk) = zcoef3 * ( zdiv (ji,jk-1) + zdiv (ji-1,jk-1)     & 
    400                                        +zdiv (ji,jk  ) + zdiv (ji-1,jk  ) )   & 
    401                            + zcoef4 * ( zdjv (ji,jk-1) + zdj1v(ji  ,jk-1)     & 
    402                                        +zdjv (ji,jk  ) + zdj1v(ji  ,jk  ) ) 
     390                  &                    +zdiv (ji,jk  ) + zdiv (ji-1,jk  ) )   & 
     391                  &        + zcoef4 * ( zdjv (ji,jk-1) + zdj1v(ji  ,jk-1)     & 
     392                  &                    +zdjv (ji,jk  ) + zdj1v(ji  ,jk  ) ) 
    403393               ! update avmv (add isopycnal vertical coefficient to avmv) 
    404                ! Caution: zcoef0 include aht0, so divided by aht0 to obtain slp^2 * aht0 
    405                avmv(ji,jj,jk) = avmv(ji,jj,jk) + ( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) / aht0 
     394               ! Caution: zcoef0 include rn_aht_0, so divided by rn_aht_0 to obtain slp^2 * rn_aht_0 
     395               avmv(ji,jj,jk) = avmv(ji,jj,jk) + ( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) / rn_aht_0 
    406396            END DO 
    407397         END DO 
     
    412402         DO jk = 1, jpkm1 
    413403            DO ji = 2, jpim1 
    414                ! volume elements 
    415                zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) 
    416                zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) 
    417                ! part of the k-component of isopycnal momentum diffusive trends 
    418                zuav = ( zfuw(ji,jk) - zfuw(ji,jk+1) ) / zbu 
    419                zvav = ( zfvw(ji,jk) - zfvw(ji,jk+1) ) / zbv 
    420                ! add the trends to the general trends 
    421                ua(ji,jj,jk) = ua(ji,jj,jk) + zuav 
    422                va(ji,jj,jk) = va(ji,jj,jk) + zvav 
     404               ua(ji,jj,jk) = ua(ji,jj,jk) + ( zfuw(ji,jk) - zfuw(ji,jk+1) ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk) 
     405               va(ji,jj,jk) = va(ji,jj,jk) + ( zfvw(ji,jk) - zfvw(ji,jk+1) ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk) 
    423406            END DO 
    424407         END DO 
     
    432415   END SUBROUTINE dyn_ldf_iso 
    433416 
    434 # else 
    435    !!---------------------------------------------------------------------- 
    436    !!   Dummy module                           NO rotation of mixing tensor 
    437    !!---------------------------------------------------------------------- 
    438 CONTAINS 
    439    SUBROUTINE dyn_ldf_iso( kt )               ! Empty routine 
    440       INTEGER, INTENT(in) :: kt 
    441       WRITE(*,*) 'dyn_ldf_iso: You should not have seen this print! error?', kt 
    442    END SUBROUTINE dyn_ldf_iso 
    443 #endif 
    444  
    445417   !!====================================================================== 
    446418END MODULE dynldf_iso 
Note: See TracChangeset for help on using the changeset viewer.