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 5845 for branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90 – NEMO

Ignore:
Timestamp:
2015-10-31T08:40:45+01:00 (8 years ago)
Author:
gm
Message:

#1613: vvl by default: suppression of domzgr_substitute.h90

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90

    r5836 r5845  
    6363 
    6464   !! * Substitutions 
    65 #  include "domzgr_substitute.h90" 
    6665#  include "vectopt_loop_substitute.h90" 
    6766   !!---------------------------------------------------------------------- 
     
    214213      !!---------------------------------------------------------------------- 
    215214      INTEGER, INTENT(in) ::   kt    ! ocean time-step index 
    216       !! 
     215      ! 
    217216      INTEGER  ::   ji, jj, jk       ! dummy loop indices 
    218217      REAL(wp) ::   zcoef0, zcoef1   ! temporary scalars 
     
    233232      DO jj = 2, jpjm1 
    234233         DO ji = fs_2, fs_jpim1   ! vector opt. 
    235             zcoef1 = zcoef0 * fse3w(ji,jj,1) 
     234            zcoef1 = zcoef0 * e3w_n(ji,jj,1) 
    236235            ! hydrostatic pressure gradient 
    237236            zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) / e1u(ji,jj) 
    238237            zhpj(ji,jj,1) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) / e2v(ji,jj) 
     238!!gm            zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj) 
     239!!gm            zhpj(ji,jj,1) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj) 
    239240            ! add to the general momentum trend 
    240241            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) 
     
    248249         DO jj = 2, jpjm1 
    249250            DO ji = fs_2, fs_jpim1   ! vector opt. 
    250                zcoef1 = zcoef0 * fse3w(ji,jj,jk) 
     251               zcoef1 = zcoef0 * e3w_n(ji,jj,jk) 
    251252               ! hydrostatic pressure gradient 
    252253               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1)   & 
    253254                  &           + zcoef1 * (  ( rhd(ji+1,jj,jk)+rhd(ji+1,jj,jk-1) )   & 
    254255                  &                       - ( rhd(ji  ,jj,jk)+rhd(ji  ,jj,jk-1) )  ) / e1u(ji,jj) 
     256!!gm                  &                       - ( rhd(ji  ,jj,jk)+rhd(ji  ,jj,jk-1) )  ) * r1_e1u(ji,jj) 
    255257 
    256258               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1)   & 
    257259                  &           + zcoef1 * (  ( rhd(ji,jj+1,jk)+rhd(ji,jj+1,jk-1) )   & 
    258260                  &                       - ( rhd(ji,jj,  jk)+rhd(ji,jj  ,jk-1) )  ) / e2v(ji,jj) 
     261!!gm                  &                       - ( rhd(ji,jj,  jk)+rhd(ji,jj  ,jk-1) )  ) * r1_e2v(ji,jj) 
    259262               ! add to the general momentum trend 
    260263               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) 
     
    300303      DO jj = 2, jpjm1 
    301304         DO ji = fs_2, fs_jpim1   ! vector opt. 
    302             zcoef1 = zcoef0 * fse3w(ji,jj,1) 
     305            zcoef1 = zcoef0 * e3w_n(ji,jj,1) 
    303306            ! hydrostatic pressure gradient 
    304307            zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj  ,1) - rhd(ji,jj,1) ) / e1u(ji,jj) 
    305308            zhpj(ji,jj,1) = zcoef1 * ( rhd(ji  ,jj+1,1) - rhd(ji,jj,1) ) / e2v(ji,jj) 
     309!!gm            zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj  ,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj) 
     310!!gm            zhpj(ji,jj,1) = zcoef1 * ( rhd(ji  ,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj) 
    306311            ! add to the general momentum trend 
    307312            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) 
     
    315320         DO jj = 2, jpjm1 
    316321            DO ji = fs_2, fs_jpim1   ! vector opt. 
    317                zcoef1 = zcoef0 * fse3w(ji,jj,jk) 
     322               zcoef1 = zcoef0 * e3w_n(ji,jj,jk) 
    318323               ! hydrostatic pressure gradient 
    319324               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1)   & 
    320325                  &           + zcoef1 * (  ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) )   & 
    321326                  &                       - ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) )  ) / e1u(ji,jj) 
     327!!gm                  &                       - ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) )  ) * r1_e1u(ji,jj) 
    322328 
    323329               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1)   & 
    324330                  &           + zcoef1 * (  ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) )   & 
    325331                  &                       - ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) )  ) / e2v(ji,jj) 
     332!!gm                  &                       - ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) )  ) * r1_e2v(ji,jj) 
    326333               ! add to the general momentum trend 
    327334               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) 
     
    337344            iku = mbku(ji,jj) 
    338345            ikv = mbkv(ji,jj) 
    339             zcoef2 = zcoef0 * MIN( fse3w(ji,jj,iku), fse3w(ji+1,jj  ,iku) ) 
    340             zcoef3 = zcoef0 * MIN( fse3w(ji,jj,ikv), fse3w(ji  ,jj+1,ikv) ) 
     346            zcoef2 = zcoef0 * MIN( e3w_n(ji,jj,iku), e3w_n(ji+1,jj  ,iku) ) 
     347            zcoef3 = zcoef0 * MIN( e3w_n(ji,jj,ikv), e3w_n(ji  ,jj+1,ikv) ) 
    341348            IF( iku > 1 ) THEN            ! on i-direction (level 2 or more) 
    342349               ua  (ji,jj,iku) = ua(ji,jj,iku) - zhpi(ji,jj,iku)         ! subtract old value 
    343350               zhpi(ji,jj,iku) = zhpi(ji,jj,iku-1)                   &   ! compute the new one 
    344351                  &            + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + gru(ji,jj) ) / e1u(ji,jj) 
     352!!gm                  &            + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + gru(ji,jj) ) * r1_e1u(ji,jj) 
    345353               ua  (ji,jj,iku) = ua(ji,jj,iku) + zhpi(ji,jj,iku)         ! add the new one to the general momentum trend 
    346354            ENDIF 
     
    349357               zhpj(ji,jj,ikv) = zhpj(ji,jj,ikv-1)                   &   ! compute the new one 
    350358                  &            + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + grv(ji,jj) ) / e2v(ji,jj) 
     359!!gm                  &            + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + grv(ji,jj) ) * r1_e2v(ji,jj) 
    351360               va  (ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv)         ! add the new one to the general momentum trend 
    352361            ENDIF 
     
    402411         DO ji = fs_2, fs_jpim1   ! vector opt. 
    403412            ! hydrostatic pressure gradient along s-surfaces 
    404             zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( fse3w(ji+1,jj  ,1) * ( znad + rhd(ji+1,jj  ,1) )   & 
    405                &                                  - fse3w(ji  ,jj  ,1) * ( znad + rhd(ji  ,jj  ,1) ) ) 
    406             zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( fse3w(ji  ,jj+1,1) * ( znad + rhd(ji  ,jj+1,1) )   & 
    407                &                                  - fse3w(ji  ,jj  ,1) * ( znad + rhd(ji  ,jj  ,1) ) ) 
     413!!gm            zhpi(ji,jj,1) = zcoef0 * r1_e1u(ji,jj) * ( e3w_n(ji+1,jj  ,1) * ( znad + rhd(ji+1,jj  ,1) )   & 
     414            zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( e3w_n(ji+1,jj  ,1) * ( znad + rhd(ji+1,jj  ,1) )   & 
     415               &                                     - e3w_n(ji  ,jj  ,1) * ( znad + rhd(ji  ,jj  ,1) ) ) 
     416!!gm            zhpj(ji,jj,1) = zcoef0 * r1_e2v(ji,jj) * ( e3w_n(ji  ,jj+1,1) * ( znad + rhd(ji  ,jj+1,1) )   & 
     417            zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( e3w_n(ji  ,jj+1,1) * ( znad + rhd(ji  ,jj+1,1) )   & 
     418               &                                     - e3w_n(ji  ,jj  ,1) * ( znad + rhd(ji  ,jj  ,1) ) ) 
    408419            ! s-coordinate pressure gradient correction 
    409             zuap = -zcoef0 * ( rhd   (ji+1,jj,1) + rhd   (ji,jj,1) + 2._wp * znad )   & 
    410                &           * ( fsde3w(ji+1,jj,1) - fsde3w(ji,jj,1) ) / e1u(ji,jj) 
    411             zvap = -zcoef0 * ( rhd   (ji,jj+1,1) + rhd   (ji,jj,1) + 2._wp * znad )   & 
    412                &           * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj) 
     420            zuap = -zcoef0 * ( rhd    (ji+1,jj,1) + rhd    (ji,jj,1) + 2._wp * znad )   & 
     421               &           * ( gde3w_n(ji+1,jj,1) - gde3w_n(ji,jj,1) ) / e1u(ji,jj) 
     422!!gm               &           * ( gde3w_n(ji+1,jj,1) - gde3w_n(ji,jj,1) ) * r1_e1u(ji,jj) 
     423            zvap = -zcoef0 * ( rhd    (ji,jj+1,1) + rhd    (ji,jj,1) + 2._wp * znad )   & 
     424               &           * ( gde3w_n(ji,jj+1,1) - gde3w_n(ji,jj,1) ) / e2v(ji,jj) 
     425!!gm               &           * ( gde3w_n(ji,jj+1,1) - gde3w_n(ji,jj,1) ) * r1_e2v(ji,jj) 
    413426            ! add to the general momentum trend 
    414427            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap 
     
    422435            DO ji = fs_2, fs_jpim1   ! vector opt. 
    423436               ! hydrostatic pressure gradient along s-surfaces 
     437!!gm               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 * r1_e1u(ji,jj)   & 
    424438               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj)   & 
    425                   &           * (  fse3w(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad )   & 
    426                   &              - fse3w(ji  ,jj,jk) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) + 2*znad )  ) 
    427                zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj)   & 
    428                   &           * (  fse3w(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad )   & 
    429                   &              - fse3w(ji,jj  ,jk) * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) + 2*znad )  ) 
     439                  &           * (  e3w_n(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad )   & 
     440                  &              - e3w_n(ji  ,jj,jk) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) + 2*znad )  ) 
     441!!gm               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj)   & 
     442               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 * r1_e2v(ji,jj)   & 
     443                  &           * (  e3w_n(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad )   & 
     444                  &              - e3w_n(ji,jj  ,jk) * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) + 2*znad )  ) 
    430445               ! s-coordinate pressure gradient correction 
    431                zuap = -zcoef0 * ( rhd   (ji+1,jj  ,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   & 
    432                   &           * ( fsde3w(ji+1,jj  ,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj) 
    433                zvap = -zcoef0 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   & 
    434                   &           * ( fsde3w(ji  ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) 
     446               zuap = -zcoef0 * ( rhd    (ji+1,jj  ,jk) + rhd    (ji,jj,jk) + 2._wp * znad )   & 
     447                  &           * ( gde3w_n(ji+1,jj  ,jk) - gde3w_n(ji,jj,jk) ) / e1u(ji,jj) 
     448!!gm                  &           * ( gde3w_n(ji+1,jj  ,jk) - gde3w_n(ji,jj,jk) ) * r1_e1u(ji,jj) 
     449               zvap = -zcoef0 * ( rhd    (ji  ,jj+1,jk) + rhd    (ji,jj,jk) + 2._wp * znad )   & 
     450                  &           * ( gde3w_n(ji  ,jj+1,jk) - gde3w_n(ji,jj,jk) ) / e2v(ji,jj) 
     451!!gm                  &           * ( gde3w_n(ji  ,jj+1,jk) - gde3w_n(ji,jj,jk) ) * r1_e2v(ji,jj) 
    435452               ! add to the general momentum trend 
    436453               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap 
     
    528545         DO ji = 1, jpi   ! vector opt. 
    529546            ikt=mikt(ji,jj) 
    530             ziceload(ji,jj) = ziceload(ji,jj) + (znad + rhd(ji,jj,1) ) * fse3w(ji,jj,1) * (1._wp - tmask(ji,jj,1)) 
     547            ziceload(ji,jj) = ziceload(ji,jj) + (znad + rhd(ji,jj,1) ) * e3w_n(ji,jj,1) * (1._wp - tmask(ji,jj,1)) 
    531548            DO jk=2,ikt-1 
    532                ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + rhd(ji,jj,jk-1) + rhd(ji,jj,jk)) * fse3w(ji,jj,jk) & 
     549               ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + rhd(ji,jj,jk-1) + rhd(ji,jj,jk)) * e3w_n(ji,jj,jk) & 
    533550                  &                              * (1._wp - tmask(ji,jj,jk)) 
    534551            END DO 
     
    544561            ! hydrostatic pressure gradient along s-surfaces and ice shelf pressure 
    545562            ! we assume ISF is in isostatic equilibrium 
    546             zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( 0.5_wp * fse3w(ji+1,jj  ,iktp1i)                                    & 
     563!!gm            zhpi(ji,jj,1) = zcoef0 * r1_e1u(ji,jj) * ( 0.5_wp * e3w_n(ji+1,jj  ,iktp1i)                                    & 
     564            zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( 0.5_wp * e3w_n(ji+1,jj  ,iktp1i)                                    & 
    547565               &                                  * ( 2._wp * znad + rhd(ji+1,jj  ,iktp1i) + zrhdtop_oce(ji+1,jj  ) )   & 
    548                &                                  - 0.5_wp * fse3w(ji  ,jj  ,ikt   )                                    & 
     566               &                                  - 0.5_wp * e3w_n(ji  ,jj  ,ikt   )                                    & 
    549567               &                                  * ( 2._wp * znad + rhd(ji  ,jj  ,ikt   ) + zrhdtop_oce(ji  ,jj  ) )   & 
    550568               &                                  + ( ziceload(ji+1,jj) - ziceload(ji,jj))                              )  
    551             zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * fse3w(ji  ,jj+1,iktp1j)                                    & 
     569!!gm            zhpj(ji,jj,1) = zcoef0 * r1_e2v(ji,jj) * ( 0.5_wp * e3w_n(ji  ,jj+1,iktp1j)                                    & 
     570            zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * e3w_n(ji  ,jj+1,iktp1j)                                    & 
    552571               &                                  * ( 2._wp * znad + rhd(ji  ,jj+1,iktp1j) + zrhdtop_oce(ji  ,jj+1) )   & 
    553                &                                  - 0.5_wp * fse3w(ji  ,jj  ,ikt   )                                    &  
     572               &                                  - 0.5_wp * e3w_n(ji  ,jj  ,ikt   )                                    &  
    554573               &                                  * ( 2._wp * znad + rhd(ji  ,jj  ,ikt   ) + zrhdtop_oce(ji  ,jj  ) )   & 
    555574               &                                  + ( ziceload(ji,jj+1) - ziceload(ji,jj) )                             )  
    556575            ! s-coordinate pressure gradient correction (=0 if z coordinate) 
    557             zuap = -zcoef0 * ( rhd   (ji+1,jj,1) + rhd   (ji,jj,1) + 2._wp * znad )   & 
    558                &           * ( fsde3w(ji+1,jj,1) - fsde3w(ji,jj,1) ) / e1u(ji,jj) 
    559             zvap = -zcoef0 * ( rhd   (ji,jj+1,1) + rhd   (ji,jj,1) + 2._wp * znad )   & 
    560                &           * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj) 
     576            zuap = -zcoef0 * ( rhd    (ji+1,jj,1) + rhd    (ji,jj,1) + 2._wp * znad )   & 
     577               &           * ( gde3w_n(ji+1,jj,1) - gde3w_n(ji,jj,1) ) / e1u(ji,jj) 
     578!!gm               &           * ( gde3w_n(ji+1,jj,1) - gde3w_n(ji,jj,1) ) * r1_e1u(ji,jj) 
     579            zvap = -zcoef0 * ( rhd    (ji,jj+1,1) + rhd    (ji,jj,1) + 2._wp * znad )   & 
     580               &           * ( gde3w_n(ji,jj+1,1) - gde3w_n(ji,jj,1) ) / e2v(ji,jj) 
     581!!gm               &           * ( gde3w_n(ji,jj+1,1) - gde3w_n(ji,jj,1) ) * r1_e2v(ji,jj) 
    561582            ! add to the general momentum trend 
    562583            ua(ji,jj,1) = ua(ji,jj,1) + (zhpi(ji,jj,1) + zuap) * umask(ji,jj,1) 
     
    569590      DO jj = 2, jpjm1 
    570591         DO ji = fs_2, fs_jpim1   ! vector opt. 
    571             iku = miku(ji,jj) ;  
    572             zpshpi(ji,jj)=0.0_wp ; zpshpj(ji,jj)=0.0_wp 
     592            iku = miku(ji,jj) 
     593            zpshpi(ji,jj) = 0._wp   ;   zpshpj(ji,jj) = 0._wp 
    573594            ze3wu  = (gdepw_0(ji+1,jj,iku+1) - gdept_0(ji+1,jj,iku)) - (gdepw_0(ji,jj,iku+1) - gdept_0(ji,jj,iku)) 
    574595            ! u direction 
    575596            IF ( iku .GT. 1 ) THEN 
    576597               ! case iku 
     598!!gm               zhpi(ji,jj,iku)   =  zcoef0 * r1_e1u(ji,jj) * ze3wu                                            & 
    577599               zhpi(ji,jj,iku)   =  zcoef0 / e1u(ji,jj) * ze3wu                                            & 
    578600                      &                                 * ( rhd    (ji+1,jj,iku) + rhd   (ji,jj,iku)       & 
    579601                      &                                   + SIGN(1._wp,ze3wu) * grui(ji,jj) + 2._wp * znad ) 
    580602               ! corrective term ( = 0 if z coordinate ) 
     603!!gm               zuap              = -zcoef0 * ( arui(ji,jj) + 2._wp * znad ) * gzui(ji,jj) * r1_e1u(ji,jj) 
    581604               zuap              = -zcoef0 * ( arui(ji,jj) + 2._wp * znad ) * gzui(ji,jj) / e1u(ji,jj) 
    582605               ! zhpi will be added in interior loop 
     
    586609 
    587610               ! case iku + 1 (remove the zphi term added in the interior loop and compute the one corrected for zps) 
     611!!gm               zhpiint        =  zcoef0 * r1_e1u(ji,jj)                                                               &     
    588612               zhpiint        =  zcoef0 / e1u(ji,jj)                                                               &     
    589                   &           * (  fse3w(ji+1,jj  ,iku+1) * ( (rhd(ji+1,jj,iku+1) + znad)                          & 
     613                  &           * (  e3w_n(ji+1,jj  ,iku+1) * ( (rhd(ji+1,jj,iku+1) + znad)                          & 
    590614                  &                                         + (rhd(ji+1,jj,iku  ) + znad) ) * tmask(ji+1,jj,iku)   & 
    591                   &              - fse3w(ji  ,jj  ,iku+1) * ( (rhd(ji  ,jj,iku+1) + znad)                          & 
     615                  &              - e3w_n(ji  ,jj  ,iku+1) * ( (rhd(ji  ,jj,iku+1) + znad)                          & 
    592616                  &                                         + (rhd(ji  ,jj,iku  ) + znad) ) * tmask(ji  ,jj,iku)   ) 
     617!!gm               zhpi(ji,jj,iku+1) =  zcoef0 * r1_e1u(ji,jj) * ge3rui(ji,jj) - zhpiint  
    593618               zhpi(ji,jj,iku+1) =  zcoef0 / e1u(ji,jj) * ge3rui(ji,jj) - zhpiint  
    594619            END IF 
     
    599624            IF ( ikv .GT. 1 ) THEN 
    600625               ! case ikv 
     626!!gm               zhpj(ji,jj,ikv)   =  zcoef0 * r1_e2v(ji,jj) * ze3wv                                            & 
    601627               zhpj(ji,jj,ikv)   =  zcoef0 / e2v(ji,jj) * ze3wv                                            & 
    602628                     &                                  * ( rhd(ji,jj+1,ikv) + rhd   (ji,jj,ikv)           & 
    603629                     &                                    + SIGN(1._wp,ze3wv) * grvi(ji,jj) + 2._wp * znad ) 
    604630               ! corrective term ( = 0 if z coordinate ) 
     631!!gm               zvap              = -zcoef0 * ( arvi(ji,jj) + 2._wp * znad ) * gzvi(ji,jj) * r1_e2v(ji,jj) 
    605632               zvap              = -zcoef0 * ( arvi(ji,jj) + 2._wp * znad ) * gzvi(ji,jj) / e2v(ji,jj) 
    606633               ! zhpi will be added in interior loop 
     
    610637                
    611638               ! case ikv + 1 (remove the zphj term added in the interior loop and compute the one corrected for zps) 
     639!!gm               zhpjint        =  zcoef0 * r1_e2v(ji,jj)                                                              & 
    612640               zhpjint        =  zcoef0 / e2v(ji,jj)                                                              & 
    613                   &           * (  fse3w(ji  ,jj+1,ikv+1) * ( (rhd(ji,jj+1,ikv+1) + znad)                         & 
     641                  &           * (  e3w_n(ji  ,jj+1,ikv+1) * ( (rhd(ji,jj+1,ikv+1) + znad)                         & 
    614642                  &                                       + (rhd(ji,jj+1,ikv  ) + znad) ) * tmask(ji,jj+1,ikv)    & 
    615                   &              - fse3w(ji  ,jj  ,ikv+1) * ( (rhd(ji,jj  ,ikv+1) + znad)                         & 
     643                  &              - e3w_n(ji  ,jj  ,ikv+1) * ( (rhd(ji,jj  ,ikv+1) + znad)                         & 
    616644                  &                                       + (rhd(ji,jj  ,ikv  ) + znad) ) * tmask(ji,jj  ,ikv)  ) 
     645!!gm               zhpj(ji,jj,ikv+1) =  zcoef0 * r1_e2v(ji,jj) * ge3rvi(ji,jj) - zhpjint 
    617646               zhpj(ji,jj,ikv+1) =  zcoef0 / e2v(ji,jj) * ge3rvi(ji,jj) - zhpjint 
    618647            END IF 
     
    626655      DO jj = 2, jpjm1 
    627656         DO ji = fs_2, fs_jpim1   ! vector opt. 
     657!!gm useles ! 
    628658            iku=miku(ji,jj); ikv=mikv(ji,jj) 
     659!!gm 
    629660            DO jk = 2, jpkm1 
    630661               ! hydrostatic pressure gradient along s-surfaces 
    631662               ! zhpi is masked for the first wet cell  (contribution already done in the upper bloc) 
    632                zhpi(ji,jj,jk) = zhpi(ji,jj,jk) + zhpi(ji,jj,jk-1)                                                              & 
    633                   &                            + zcoef0 / e1u(ji,jj)                                                           & 
    634                   &                                     * ( fse3w(ji+1,jj  ,jk) * ( (rhd(ji+1,jj,jk  ) + znad)                 & 
    635                   &                                                     + (rhd(ji+1,jj,jk-1) + znad) ) * tmask(ji+1,jj,jk-1)   & 
    636                   &                                       - fse3w(ji  ,jj  ,jk) * ( (rhd(ji  ,jj,jk  ) + znad)                 & 
    637                   &                                                     + (rhd(ji  ,jj,jk-1) + znad) ) * tmask(ji  ,jj,jk-1)   )  
     663               zhpi(ji,jj,jk) = zhpi(ji,jj,jk) + zhpi(ji,jj,jk-1)                                                     & 
     664!!gm                  &           + zcoef0 * r1_e1u(ji,jj)                                                                & 
     665                  &           + zcoef0 / e1u(ji,jj)                                                                & 
     666                  &                    * ( e3w_n(ji+1,jj,jk) * ( (rhd(ji+1,jj,jk  ) + znad)                           & 
     667                  &                                            + (rhd(ji+1,jj,jk-1) + znad) ) * tmask(ji+1,jj,jk-1)   & 
     668                  &                      - e3w_n(ji  ,jj,jk) * ( (rhd(ji  ,jj,jk  ) + znad)                           & 
     669                  &                                            + (rhd(ji  ,jj,jk-1) + znad) ) * tmask(ji  ,jj,jk-1)   )  
    638670               ! s-coordinate pressure gradient correction 
    639671               ! corrective term, we mask this term for the first wet level beneath the ice shelf (contribution done in the upper bloc)  
    640                zuap = - zcoef0 * ( rhd   (ji+1,jj  ,jk) + rhd   (ji,jj,jk) + 2._wp * znad )                    & 
    641                   &            * ( fsde3w(ji+1,jj  ,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj) * umask(ji,jj,jk-1) 
     672               zuap = - zcoef0 * ( rhd    (ji+1,jj  ,jk) + rhd    (ji,jj,jk) + 2._wp * znad )                    & 
     673                  &            * ( gde3w_n(ji+1,jj  ,jk) - gde3w_n(ji,jj,jk) ) / e1u(ji,jj) * umask(ji,jj,jk-1) 
     674!!gm                  &            * ( gde3w_n(ji+1,jj  ,jk) - gde3w_n(ji,jj,jk) ) * r1_e1u(ji,jj) * umask(ji,jj,jk-1) 
    642675               ua(ji,jj,jk) = ua(ji,jj,jk) + ( zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk) 
    643676 
    644677               ! hydrostatic pressure gradient along s-surfaces 
    645678               ! zhpi is masked for the first wet cell  (contribution already done in the upper bloc) 
    646                zhpj(ji,jj,jk) = zhpj(ji,jj,jk) + zhpj(ji,jj,jk-1)                                                              & 
    647                   &                            + zcoef0 / e2v(ji,jj)                                                           & 
    648                   &                                     * ( fse3w(ji  ,jj+1,jk) * ( (rhd(ji,jj+1,jk  ) + znad)                 & 
    649                   &                                                     + (rhd(ji,jj+1,jk-1) + znad) ) * tmask(ji,jj+1,jk-1)   & 
    650                   &                                       - fse3w(ji  ,jj  ,jk) * ( (rhd(ji,jj  ,jk  ) + znad)                 & 
    651                   &                                                     + (rhd(ji,jj  ,jk-1) + znad) ) * tmask(ji,jj  ,jk-1)   ) 
     679               zhpj(ji,jj,jk) = zhpj(ji,jj,jk) + zhpj(ji,jj,jk-1)                                                     & 
     680!!gm                  &           + zcoef0 * r1_e2v(ji,jj)                                                                & 
     681                  &           + zcoef0 / e2v(ji,jj)                                                                & 
     682                  &                    * ( e3w_n(ji  ,jj+1,jk) * ( (rhd(ji,jj+1,jk  ) + znad)                           & 
     683                  &                                              + (rhd(ji,jj+1,jk-1) + znad) ) * tmask(ji,jj+1,jk-1)   & 
     684                  &                      - e3w_n(ji  ,jj  ,jk) * ( (rhd(ji,jj  ,jk  ) + znad)                           & 
     685                  &                                              + (rhd(ji,jj  ,jk-1) + znad) ) * tmask(ji,jj  ,jk-1)   ) 
    652686               ! s-coordinate pressure gradient correction 
    653687               ! corrective term, we mask this term for the first wet level beneath the ice shelf (contribution done in the upper bloc) 
    654                zvap = - zcoef0 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji,jj,jk) + 2._wp * znad )                     & 
    655                   &            * ( fsde3w(ji  ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) * vmask(ji,jj,jk-1) 
     688               zvap = - zcoef0 * ( rhd    (ji  ,jj+1,jk) + rhd    (ji,jj,jk) + 2._wp * znad )                    & 
     689                  &            * ( gde3w_n(ji  ,jj+1,jk) - gde3w_n(ji,jj,jk) ) / e2v(ji,jj) * vmask(ji,jj,jk-1) 
     690!!gm                  &            * ( gde3w_n(ji  ,jj+1,jk) - gde3w_n(ji,jj,jk) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk-1) 
    656691               ! add to the general momentum trend 
    657692               va(ji,jj,jk) = va(ji,jj,jk) + ( zhpj(ji,jj,jk) + zvap ) * vmask(ji,jj,jk) 
     
    671706            IF (iku .GT. 1) THEN 
    672707               ! remove old value (interior case) 
    673                zuap            = -zcoef0 * ( rhd   (ji+1,jj  ,iku) + rhd   (ji,jj,iku) + 2._wp * znad )   & 
    674                      &                   * ( fsde3w(ji+1,jj  ,iku) - fsde3w(ji,jj,iku) ) / e1u(ji,jj) 
     708               zuap            = -zcoef0 * ( rhd    (ji+1,jj  ,iku) + rhd    (ji,jj,iku) + 2._wp * znad )   & 
     709                     &                   * ( gde3w_n(ji+1,jj  ,iku) - gde3w_n(ji,jj,iku) ) / e1u(ji,jj) 
     710!!gm                     &                   * ( gde3w_n(ji+1,jj  ,iku) - gde3w_n(ji,jj,iku) ) * r1_e1u(ji,jj) 
    675711               ua(ji,jj,iku)   = ua(ji,jj,iku) - zhpi(ji,jj,iku) - zuap 
    676712               ! put new value 
    677713               ! -zpshpi to avoid double contribution of the partial step in the top layer  
     714!!gm               zuap            = -zcoef0 * ( aru(ji,jj) + 2._wp * znad ) * gzu(ji,jj)  * r1_e1u(ji,jj) 
    678715               zuap            = -zcoef0 * ( aru(ji,jj) + 2._wp * znad ) * gzu(ji,jj)  / e1u(ji,jj) 
    679716               zhpi(ji,jj,iku) =  zhpi(ji,jj,iku-1) + zcoef0 / e1u(ji,jj) * ge3ru(ji,jj) - zpshpi(ji,jj)  
     717!!gm               zhpi(ji,jj,iku) =  zhpi(ji,jj,iku-1) + zcoef0 * r1_e1u(ji,jj) * ge3ru(ji,jj) - zpshpi(ji,jj)  
    680718               ua(ji,jj,iku)   =  ua(ji,jj,iku) + zhpi(ji,jj,iku) + zuap 
    681719            END IF 
     
    683721            IF (ikv .GT. 1) THEN 
    684722               ! remove old value (interior case) 
    685                zvap            = -zcoef0 * ( rhd   (ji  ,jj+1,ikv) + rhd   (ji,jj,ikv) + 2._wp * znad )   & 
    686                      &                   * ( fsde3w(ji  ,jj+1,ikv) - fsde3w(ji,jj,ikv) )   / e2v(ji,jj) 
     723               zvap            = -zcoef0 * ( rhd    (ji  ,jj+1,ikv) + rhd    (ji,jj,ikv) + 2._wp * znad )   & 
     724                     &                   * ( gde3w_n(ji  ,jj+1,ikv) - gde3w_n(ji,jj,ikv) )   / e2v(ji,jj) 
     725!!gm                     &                   * ( gde3w_n(ji  ,jj+1,ikv) - gde3w_n(ji,jj,ikv) )   * r1_e2v(ji,jj) 
    687726               va(ji,jj,ikv)   = va(ji,jj,ikv) - zhpj(ji,jj,ikv) - zvap 
    688727               ! put new value 
    689728               ! -zpshpj to avoid double contribution of the partial step in the top layer  
     729!!gm               zvap            = -zcoef0 * ( arv(ji,jj) + 2._wp * znad ) * gzv(ji,jj)     * r1_e2v(ji,jj) 
    690730               zvap            = -zcoef0 * ( arv(ji,jj) + 2._wp * znad ) * gzv(ji,jj)     / e2v(ji,jj) 
    691731               zhpj(ji,jj,ikv) =  zhpj(ji,jj,ikv-1) + zcoef0 / e2v(ji,jj) * ge3rv(ji,jj) - zpshpj(ji,jj)    
     732!!gm               zhpj(ji,jj,ikv) =  zhpj(ji,jj,ikv-1) + zcoef0 * r1_e2v(ji,jj) * ge3rv(ji,jj) - zpshpj(ji,jj)    
    692733               va(ji,jj,ikv)   =  va(ji,jj,ikv) + zhpj(ji,jj,ikv) + zvap 
    693734            END IF 
     
    750791         DO jj = 2, jpjm1 
    751792            DO ji = fs_2, fs_jpim1   ! vector opt. 
    752                drhoz(ji,jj,jk) = rhd   (ji  ,jj  ,jk) - rhd   (ji,jj,jk-1) 
    753                dzz  (ji,jj,jk) = fsde3w(ji  ,jj  ,jk) - fsde3w(ji,jj,jk-1) 
    754                drhox(ji,jj,jk) = rhd   (ji+1,jj  ,jk) - rhd   (ji,jj,jk  ) 
    755                dzx  (ji,jj,jk) = fsde3w(ji+1,jj  ,jk) - fsde3w(ji,jj,jk  ) 
    756                drhoy(ji,jj,jk) = rhd   (ji  ,jj+1,jk) - rhd   (ji,jj,jk  ) 
    757                dzy  (ji,jj,jk) = fsde3w(ji  ,jj+1,jk) - fsde3w(ji,jj,jk  ) 
     793               drhoz(ji,jj,jk) = rhd    (ji  ,jj  ,jk) - rhd    (ji,jj,jk-1) 
     794               dzz  (ji,jj,jk) = gde3w_n(ji  ,jj  ,jk) - gde3w_n(ji,jj,jk-1) 
     795               drhox(ji,jj,jk) = rhd    (ji+1,jj  ,jk) - rhd    (ji,jj,jk  ) 
     796               dzx  (ji,jj,jk) = gde3w_n(ji+1,jj  ,jk) - gde3w_n(ji,jj,jk  ) 
     797               drhoy(ji,jj,jk) = rhd    (ji  ,jj+1,jk) - rhd    (ji,jj,jk  ) 
     798               dzy  (ji,jj,jk) = gde3w_n(ji  ,jj+1,jk) - gde3w_n(ji,jj,jk  ) 
    758799            END DO 
    759800         END DO 
     
    837878      !------------------------------------------------------------- 
    838879 
    839 !!bug gm   :  e3w-de3w = 0.5*e3w  ....  and de3w(2)-de3w(1)=e3w(2) ....   to be verified 
    840 !          true if de3w is really defined as the sum of the e3w scale factors as, it seems to me, it should be 
     880!!bug gm   :  e3w-gde3w = 0.5*e3w  ....  and gde3w(2)-gde3w(1)=e3w(2) ....   to be verified 
     881!          true if gde3w is really defined as the sum of the e3w scale factors as, it seems to me, it should be 
    841882 
    842883      DO jj = 2, jpjm1 
    843884         DO ji = fs_2, fs_jpim1   ! vector opt. 
    844             rho_k(ji,jj,1) = -grav * ( fse3w(ji,jj,1) - fsde3w(ji,jj,1) )               & 
    845                &                   * (  rhd(ji,jj,1)                                    & 
    846                &                     + 0.5_wp * ( rhd(ji,jj,2) - rhd(ji,jj,1) )         & 
    847                &                              * ( fse3w (ji,jj,1) - fsde3w(ji,jj,1) )   & 
    848                &                              / ( fsde3w(ji,jj,2) - fsde3w(ji,jj,1) )  ) 
     885            rho_k(ji,jj,1) = -grav * ( e3w_n(ji,jj,1) - gde3w_n(ji,jj,1) )               & 
     886               &                   * (  rhd(ji,jj,1)                                     & 
     887               &                     + 0.5_wp * ( rhd    (ji,jj,2) - rhd    (ji,jj,1) )  & 
     888               &                              * ( e3w_n  (ji,jj,1) - gde3w_n(ji,jj,1) )  & 
     889               &                              / ( gde3w_n(ji,jj,2) - gde3w_n(ji,jj,1) )  ) 
    849890         END DO 
    850891      END DO 
     
    857898            DO ji = fs_2, fs_jpim1   ! vector opt. 
    858899 
    859                rho_k(ji,jj,jk) = zcoef0 * ( rhd   (ji,jj,jk) + rhd   (ji,jj,jk-1) )                                   & 
    860                   &                     * ( fsde3w(ji,jj,jk) - fsde3w(ji,jj,jk-1) )                                   & 
    861                   &            - grav * z1_10 * (                                                                     & 
    862                   &     ( drhow (ji,jj,jk) - drhow (ji,jj,jk-1) )                                                     & 
    863                   &   * ( fsde3w(ji,jj,jk) - fsde3w(ji,jj,jk-1) - z1_12 * ( dzw  (ji,jj,jk) + dzw  (ji,jj,jk-1) ) )   & 
    864                   &   - ( dzw   (ji,jj,jk) - dzw   (ji,jj,jk-1) )                                                     & 
    865                   &   * ( rhd   (ji,jj,jk) - rhd   (ji,jj,jk-1) - z1_12 * ( drhow(ji,jj,jk) + drhow(ji,jj,jk-1) ) )   & 
     900               rho_k(ji,jj,jk) = zcoef0 * ( rhd    (ji,jj,jk) + rhd    (ji,jj,jk-1) )                                   & 
     901                  &                     * ( gde3w_n(ji,jj,jk) - gde3w_n(ji,jj,jk-1) )                                   & 
     902                  &            - grav * z1_10 * (                                                                   & 
     903                  &     ( drhow  (ji,jj,jk) - drhow (ji,jj,jk-1) )                                                     & 
     904                  &   * ( gde3w_n(ji,jj,jk) - gde3w_n(ji,jj,jk-1) - z1_12 * ( dzw  (ji,jj,jk) + dzw  (ji,jj,jk-1) ) )   & 
     905                  &   - ( dzw    (ji,jj,jk) - dzw    (ji,jj,jk-1) )                                                     & 
     906                  &   * ( rhd    (ji,jj,jk) - rhd    (ji,jj,jk-1) - z1_12 * ( drhow(ji,jj,jk) + drhow(ji,jj,jk-1) ) )   & 
    866907                  &                             ) 
    867908 
    868                rho_i(ji,jj,jk) = zcoef0 * ( rhd   (ji+1,jj,jk) + rhd   (ji,jj,jk) )                                   & 
    869                   &                     * ( fsde3w(ji+1,jj,jk) - fsde3w(ji,jj,jk) )                                   & 
    870                   &            - grav* z1_10 * (                                                                      & 
    871                   &     ( drhou (ji+1,jj,jk) - drhou (ji,jj,jk) )                                                     & 
    872                   &   * ( fsde3w(ji+1,jj,jk) - fsde3w(ji,jj,jk) - z1_12 * ( dzu  (ji+1,jj,jk) + dzu  (ji,jj,jk) ) )   & 
    873                   &   - ( dzu   (ji+1,jj,jk) - dzu   (ji,jj,jk) )                                                     & 
    874                   &   * ( rhd   (ji+1,jj,jk) - rhd   (ji,jj,jk) - z1_12 * ( drhou(ji+1,jj,jk) + drhou(ji,jj,jk) ) )   & 
     909               rho_i(ji,jj,jk) = zcoef0 * ( rhd    (ji+1,jj,jk) + rhd    (ji,jj,jk) )                                   & 
     910                  &                     * ( gde3w_n(ji+1,jj,jk) - gde3w_n(ji,jj,jk) )                                   & 
     911                  &            - grav* z1_10 * (                                                                    & 
     912                  &     ( drhou  (ji+1,jj,jk) - drhou (ji,jj,jk) )                                                     & 
     913                  &   * ( gde3w_n(ji+1,jj,jk) - gde3w_n(ji,jj,jk) - z1_12 * ( dzu  (ji+1,jj,jk) + dzu  (ji,jj,jk) ) )   & 
     914                  &   - ( dzu    (ji+1,jj,jk) - dzu    (ji,jj,jk) )                                                     & 
     915                  &   * ( rhd    (ji+1,jj,jk) - rhd    (ji,jj,jk) - z1_12 * ( drhou(ji+1,jj,jk) + drhou(ji,jj,jk) ) )   & 
    875916                  &                            ) 
    876917 
    877                rho_j(ji,jj,jk) = zcoef0 * ( rhd   (ji,jj+1,jk) + rhd   (ji,jj,jk) )                                   & 
    878                   &                     * ( fsde3w(ji,jj+1,jk) - fsde3w(ji,jj,jk) )                                   & 
    879                   &            - grav* z1_10 * (                                                                      & 
    880                   &     ( drhov (ji,jj+1,jk) - drhov (ji,jj,jk) )                                                     & 
    881                   &   * ( fsde3w(ji,jj+1,jk) - fsde3w(ji,jj,jk) - z1_12 * ( dzv  (ji,jj+1,jk) + dzv  (ji,jj,jk) ) )   & 
    882                   &   - ( dzv   (ji,jj+1,jk) - dzv   (ji,jj,jk) )                                                     & 
    883                   &   * ( rhd   (ji,jj+1,jk) - rhd   (ji,jj,jk) - z1_12 * ( drhov(ji,jj+1,jk) + drhov(ji,jj,jk) ) )   & 
     918               rho_j(ji,jj,jk) = zcoef0 * ( rhd    (ji,jj+1,jk) + rhd    (ji,jj,jk) )                                 & 
     919                  &                     * ( gde3w_n(ji,jj+1,jk) - gde3w_n(ji,jj,jk) )                                   & 
     920                  &            - grav* z1_10 * (                                                                    & 
     921                  &     ( drhov  (ji,jj+1,jk) - drhov (ji,jj,jk) )                                                     & 
     922                  &   * ( gde3w_n(ji,jj+1,jk) - gde3w_n(ji,jj,jk) - z1_12 * ( dzv  (ji,jj+1,jk) + dzv  (ji,jj,jk) ) )   & 
     923                  &   - ( dzv    (ji,jj+1,jk) - dzv    (ji,jj,jk) )                                                     & 
     924                  &   * ( rhd    (ji,jj+1,jk) - rhd    (ji,jj,jk) - z1_12 * ( drhov(ji,jj+1,jk) + drhov(ji,jj,jk) ) )   & 
    884925                  &                            ) 
    885926 
     
    897938      DO jj = 2, jpjm1 
    898939         DO ji = fs_2, fs_jpim1   ! vector opt. 
     940!!gm            zhpi(ji,jj,1) = ( rho_k(ji+1,jj  ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) * r1_e1u(ji,jj) 
    899941            zhpi(ji,jj,1) = ( rho_k(ji+1,jj  ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) / e1u(ji,jj) 
    900942            zhpj(ji,jj,1) = ( rho_k(ji  ,jj+1,1) - rho_k(ji,jj,1) - rho_j(ji,jj,1) ) / e2v(ji,jj) 
     943!!gm            zhpj(ji,jj,1) = ( rho_k(ji  ,jj+1,1) - rho_k(ji,jj,1) - rho_j(ji,jj,1) ) * r1_e2v(ji,jj) 
    901944            ! add to the general momentum trend 
    902945            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) 
     
    915958                  &           + (  ( rho_k(ji+1,jj,jk) - rho_k(ji,jj,jk  ) )    & 
    916959                  &              - ( rho_i(ji  ,jj,jk) - rho_i(ji,jj,jk-1) )  ) / e1u(ji,jj) 
     960!!gm                  &              - ( rho_i(ji  ,jj,jk) - rho_i(ji,jj,jk-1) )  ) * r1_e1u(ji,jj) 
    917961               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1)                                & 
    918962                  &           + (  ( rho_k(ji,jj+1,jk) - rho_k(ji,jj,jk  ) )    & 
    919963                  &               -( rho_j(ji,jj  ,jk) - rho_j(ji,jj,jk-1) )  ) / e2v(ji,jj) 
     964!!gm                  &               -( rho_j(ji,jj  ,jk) - rho_j(ji,jj,jk-1) )  ) * r1_e2v(ji,jj) 
    920965               ! add to the general momentum trend 
    921966               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) 
     
    9831028        DO ji = 1, jpi 
    9841029          jk = mbathy(ji,jj) 
    985           IF( jk <= 0 ) THEN; zrhh(ji,jj,:) = 0._wp 
    986           ELSE IF(jk == 1) THEN; zrhh(ji,jj, jk+1:jpk) = rhd(ji,jj,jk) 
    987           ELSE IF(jk < jpkm1) THEN 
     1030          IF(     jk <=  0   ) THEN   ;   zrhh(ji,jj,    :   ) = 0._wp 
     1031          ELSEIF( jk ==  1   ) THEN   ;   zrhh(ji,jj,jk+1:jpk) = rhd(ji,jj,jk) 
     1032          ELSEIF( jk < jpkm1 ) THEN 
    9881033             DO jkk = jk+1, jpk 
    989                 zrhh(ji,jj,jkk) = interp1(fsde3w(ji,jj,jkk),   fsde3w(ji,jj,jkk-1), & 
    990                                          fsde3w(ji,jj,jkk-2), rhd(ji,jj,jkk-1), rhd(ji,jj,jkk-2)) 
     1034                zrhh(ji,jj,jkk) = interp1(gde3w_n(ji,jj,jkk  ), gde3w_n(ji,jj,jkk-1),  & 
     1035                   &                      gde3w_n(ji,jj,jkk-2), rhd    (ji,jj,jkk-1), rhd(ji,jj,jkk-2)) 
    9911036             END DO 
    9921037          ENDIF 
     
    9971042      DO jj = 1, jpj 
    9981043         DO ji = 1, jpi 
    999             zdept(ji,jj,1) = 0.5_wp * fse3w(ji,jj,1) - sshn(ji,jj) * znad 
     1044            zdept(ji,jj,1) = 0.5_wp * e3w_n(ji,jj,1) - sshn(ji,jj) * znad 
    10001045         END DO 
    10011046      END DO 
     
    10041049         DO jj = 1, jpj 
    10051050            DO ji = 1, jpi 
    1006                zdept(ji,jj,jk) = zdept(ji,jj,jk-1) + fse3w(ji,jj,jk) 
     1051               zdept(ji,jj,jk) = zdept(ji,jj,jk-1) + e3w_n(ji,jj,jk) 
    10071052            END DO 
    10081053         END DO 
     
    10221067          zrhdt1 = zrhh(ji,jj,1) - interp3(zdept(ji,jj,1),asp(ji,jj,1), & 
    10231068                                         bsp(ji,jj,1),   csp(ji,jj,1), & 
    1024                                          dsp(ji,jj,1) ) * 0.25_wp * fse3w(ji,jj,1) 
     1069                                         dsp(ji,jj,1) ) * 0.25_wp * e3w_n(ji,jj,1) 
    10251070 
    10261071          ! assuming linear profile across the top half surface layer 
    1027           zhpi(ji,jj,1) =  0.5_wp * fse3w(ji,jj,1) * zrhdt1 
     1072          zhpi(ji,jj,1) =  0.5_wp * e3w_n(ji,jj,1) * zrhdt1 
    10281073        END DO 
    10291074      END DO 
     
    10551100      DO jj = 2, jpjm1 
    10561101        DO ji = 2, jpim1 
    1057           zu(ji,jj,1) = - ( fse3u(ji,jj,1) - zsshu_n(ji,jj) * znad)  
    1058           zv(ji,jj,1) = - ( fse3v(ji,jj,1) - zsshv_n(ji,jj) * znad) 
     1102          zu(ji,jj,1) = - ( e3u_n(ji,jj,1) - zsshu_n(ji,jj) * znad)  
     1103          zv(ji,jj,1) = - ( e3v_n(ji,jj,1) - zsshv_n(ji,jj) * znad) 
    10591104        END DO 
    10601105      END DO 
     
    10631108        DO jj = 2, jpjm1 
    10641109          DO ji = 2, jpim1 
    1065             zu(ji,jj,jk) = zu(ji,jj,jk-1)- fse3u(ji,jj,jk) 
    1066             zv(ji,jj,jk) = zv(ji,jj,jk-1)- fse3v(ji,jj,jk) 
     1110            zu(ji,jj,jk) = zu(ji,jj,jk-1)- e3u_n(ji,jj,jk) 
     1111            zv(ji,jj,jk) = zv(ji,jj,jk-1)- e3v_n(ji,jj,jk) 
    10671112          END DO 
    10681113        END DO 
     
    10721117        DO jj = 2, jpjm1 
    10731118          DO ji = 2, jpim1 
    1074             zu(ji,jj,jk) = zu(ji,jj,jk) + 0.5_wp * fse3u(ji,jj,jk) 
    1075             zv(ji,jj,jk) = zv(ji,jj,jk) + 0.5_wp * fse3v(ji,jj,jk) 
     1119            zu(ji,jj,jk) = zu(ji,jj,jk) + 0.5_wp * e3u_n(ji,jj,jk) 
     1120            zv(ji,jj,jk) = zv(ji,jj,jk) + 0.5_wp * e3v_n(ji,jj,jk) 
    10761121          END DO 
    10771122        END DO 
     
    11421187               ! update the momentum trends in u direction 
    11431188 
     1189!!gm               zdpdx1 = zcoef0 * r1_e1u(ji,jj) * (zhpi(ji+1,jj,jk) - zhpi(ji,jj,jk)) 
    11441190               zdpdx1 = zcoef0 / e1u(ji,jj) * (zhpi(ji+1,jj,jk) - zhpi(ji,jj,jk)) 
    11451191               IF( lk_vvl ) THEN 
     1192!!gm                 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * & 
    11461193                 zdpdx2 = zcoef0 / e1u(ji,jj) * & 
    11471194                         ( REAL(jis-jid, wp) * (zpwes + zpwed) + (sshn(ji+1,jj)-sshn(ji,jj)) ) 
    11481195                ELSE 
     1196!!gm                 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed) 
    11491197                 zdpdx2 = zcoef0 / e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed) 
    11501198               ENDIF 
     
    11991247               ! update the momentum trends in v direction 
    12001248 
     1249!!gm               zdpdy1 = zcoef0 * r1_e2v(ji,jj) * (zhpi(ji,jj+1,jk) - zhpi(ji,jj,jk)) 
    12011250               zdpdy1 = zcoef0 / e2v(ji,jj) * (zhpi(ji,jj+1,jk) - zhpi(ji,jj,jk)) 
    12021251               IF( lk_vvl ) THEN 
     1252!!gm                   zdpdy2 = zcoef0 * r1_e2v(ji,jj) * & 
    12031253                   zdpdy2 = zcoef0 / e2v(ji,jj) * & 
    12041254                           ( REAL(jjs-jjd, wp) * (zpnss + zpnsd) + (sshn(ji,jj+1)-sshn(ji,jj)) ) 
    12051255               ELSE 
     1256!!gm                   zdpdy2 = zcoef0 * r1_e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd ) 
    12061257                   zdpdy2 = zcoef0 / e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd ) 
    12071258               ENDIF 
Note: See TracChangeset for help on using the changeset viewer.