Ignore:
Timestamp:
2020-07-30T17:05:58+02:00 (8 weeks ago)
Author:
techene
Message:

hydrostatic pressure gradient is computed with density anomaly when possible

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG/src/OCE/DYN/dynhpg.F90

    r13295 r13364  
    150150      !! 
    151151      INTEGER  ::   ji, jj, jk, ikt    ! dummy loop indices      ISF 
    152       REAL(wp) ::   znad 
    153152      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  zts_top, zrhd   ! hypothesys on isf density 
    154153      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::  zrhdtop_isf    ! density at bottom of ISF 
     
    245244      INTEGER  ::   ji, jj, jk       ! dummy loop indices 
    246245      REAL(wp) ::   zcoef0, zcoef1   ! temporary scalars 
    247       REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zhpi, zhpj 
     246      REAL(wp), DIMENSION(jpi,jpj) ::  zhpi, zhpj 
    248247      !!---------------------------------------------------------------------- 
    249248      ! 
     
    253252         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   z-coordinate case ' 
    254253      ENDIF 
    255  
    256       zcoef0 = - grav * 0.5_wp      ! Local constant initialization 
    257  
    258       ! Surface value 
    259       DO_2D( 0, 0, 0, 0 ) 
     254      ! 
     255      zcoef0 = - grav * 0.5_wp            ! Local constant initialization 
     256      ! 
     257      DO_2D( 0, 0, 0, 0 )                 ! Surface value 
    260258         zcoef1 = zcoef0 * e3w(ji,jj,1,Kmm) 
    261          ! hydrostatic pressure gradient 
    262          zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj) 
    263          zhpj(ji,jj,1) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj) 
    264          ! add to the general momentum trend 
    265          puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + zhpi(ji,jj,1) 
    266          pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + zhpj(ji,jj,1) 
     259         !                                   ! hydrostatic pressure gradient 
     260         zhpi(ji,jj) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj) 
     261         zhpj(ji,jj) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj) 
     262         !                                   ! add to the general momentum trend 
     263         puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + zhpi(ji,jj) 
     264         pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + zhpj(ji,jj) 
    267265      END_2D 
    268  
    269       ! 
    270       ! interior value (2=<jk=<jpkm1) 
    271       DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     266      ! 
     267      DO_3D( 0, 0, 0, 0, 2, jpkm1 )        ! interior value (2=<jk=<jpkm1) 
    272268         zcoef1 = zcoef0 * e3w(ji,jj,jk,Kmm) 
    273          ! hydrostatic pressure gradient 
    274          zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1)   & 
    275             &           + zcoef1 * (  ( rhd(ji+1,jj,jk)+rhd(ji+1,jj,jk-1) )    & 
    276             &                       - ( rhd(ji  ,jj,jk)+rhd(ji  ,jj,jk-1) )  ) * r1_e1u(ji,jj) 
    277  
    278          zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1)   & 
    279             &           + zcoef1 * (  ( rhd(ji,jj+1,jk)+rhd(ji,jj+1,jk-1) )    & 
    280             &                       - ( rhd(ji,jj,  jk)+rhd(ji,jj  ,jk-1) )  ) * r1_e2v(ji,jj) 
    281          ! add to the general momentum trend 
    282          puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zhpi(ji,jj,jk) 
    283          pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zhpj(ji,jj,jk) 
     269         !                                   ! hydrostatic pressure gradient 
     270         zhpi(ji,jj) = zhpi(ji,jj) + zcoef1 * (  ( rhd(ji+1,jj,jk)+rhd(ji+1,jj,jk-1) )  & 
     271            &                                  - ( rhd(ji  ,jj,jk)+rhd(ji  ,jj,jk-1) )  ) * r1_e1u(ji,jj) 
     272 
     273         zhpj(ji,jj) = zhpj(ji,jj) + zcoef1 * (  ( rhd(ji,jj+1,jk)+rhd(ji,jj+1,jk-1) )  & 
     274            &                                  - ( rhd(ji,jj,  jk)+rhd(ji,jj  ,jk-1) )  ) * r1_e2v(ji,jj) 
     275         !                                   ! add to the general momentum trend 
     276         puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zhpi(ji,jj) 
     277         pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zhpj(ji,jj) 
    284278      END_3D 
    285279      ! 
     
    390384      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv    ! ocean velocities and RHS of momentum equation 
    391385      !! 
    392       INTEGER  ::   ji, jj, jk, jii, jjj                 ! dummy loop indices 
    393       REAL(wp) ::   zcoef0, zuap, zvap, znad, ztmp       ! temporary scalars 
    394       LOGICAL  ::   ll_tmp1, ll_tmp2                     ! local logical variables 
     386      INTEGER  ::   ji, jj, jk, jii, jjj           ! dummy loop indices 
     387      REAL(wp) ::   zcoef0, zuap, zvap, ztmp       ! local scalars 
     388      LOGICAL  ::   ll_tmp1, ll_tmp2               ! local logical variables 
    395389      REAL(wp), DIMENSION(jpi,jpj,jpk)      ::   zhpi, zhpj 
    396390      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zcpx, zcpy   !W/D pressure filter 
     
    406400      ! 
    407401      zcoef0 = - grav * 0.5_wp 
    408       IF ( ln_linssh ) THEN   ;   znad = 0._wp         ! Fixed    volume: density anomaly 
    409       ELSE                    ;   znad = 1._wp         ! Variable volume: density 
    410       ENDIF 
    411402      ! 
    412403      IF( ln_wd_il ) THEN 
     
    450441        CALL lbc_lnk_multi( 'dynhpg', zcpx, 'U', 1.0_wp, zcpy, 'V', 1.0_wp ) 
    451442      END IF 
    452  
    453       ! Surface value 
    454       DO_2D( 0, 0, 0, 0 ) 
    455          ! hydrostatic pressure gradient along s-surfaces 
    456          zhpi(ji,jj,1) =   & 
    457             &  zcoef0 * (  e3w(ji+1,jj  ,1,Kmm) * ( znad + rhd(ji+1,jj  ,1) )    & 
    458             &            - e3w(ji  ,jj  ,1,Kmm) * ( znad + rhd(ji  ,jj  ,1) )  ) & 
    459             &           * r1_e1u(ji,jj) 
    460          zhpj(ji,jj,1) =   & 
    461             &  zcoef0 * (  e3w(ji  ,jj+1,1,Kmm) * ( znad + rhd(ji  ,jj+1,1) )    & 
    462             &            - e3w(ji  ,jj  ,1,Kmm) * ( znad + rhd(ji  ,jj  ,1) )  ) & 
    463             &           * r1_e2v(ji,jj) 
    464          ! s-coordinate pressure gradient correction 
    465          zuap = -zcoef0 * ( rhd    (ji+1,jj,1) + rhd    (ji,jj,1) + 2._wp * znad )   & 
     443      ! 
     444      DO_2D( 0, 0, 0, 0 )              ! Surface value 
     445         !                                   ! hydrostatic pressure gradient along s-surfaces 
     446         zhpi(ji,jj,1) = zcoef0 * r1_e1u(ji,jj)                      & 
     447            &          * (  e3w(ji+1,jj  ,1,Kmm) * rhd(ji+1,jj  ,1)  & 
     448            &             - e3w(ji  ,jj  ,1,Kmm) * rhd(ji  ,jj  ,1)  ) 
     449         zhpj(ji,jj,1) = zcoef0 * r1_e2v(ji,jj)                      & 
     450            &          * (  e3w(ji  ,jj+1,1,Kmm) * rhd(ji  ,jj+1,1)  & 
     451            &             - e3w(ji  ,jj  ,1,Kmm) * rhd(ji  ,jj  ,1)  ) 
     452         !                                   ! s-coordinate pressure gradient correction 
     453         zuap = -zcoef0 * ( rhd    (ji+1,jj,1) + rhd    (ji,jj,1) )   & 
    466454            &           * ( gde3w(ji+1,jj,1) - gde3w(ji,jj,1) ) * r1_e1u(ji,jj) 
    467          zvap = -zcoef0 * ( rhd    (ji,jj+1,1) + rhd    (ji,jj,1) + 2._wp * znad )   & 
     455         zvap = -zcoef0 * ( rhd    (ji,jj+1,1) + rhd    (ji,jj,1) )   & 
    468456            &           * ( gde3w(ji,jj+1,1) - gde3w(ji,jj,1) ) * r1_e2v(ji,jj) 
    469457         ! 
     
    474462            zvap = zvap * zcpy(ji,jj) 
    475463         ENDIF 
    476          ! 
    477          ! add to the general momentum trend 
     464         !                                   ! add to the general momentum trend 
    478465         puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + zhpi(ji,jj,1) + zuap 
    479466         pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + zhpj(ji,jj,1) + zvap 
    480467      END_2D 
    481  
    482       ! interior value (2=<jk=<jpkm1) 
    483       DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
    484          ! hydrostatic pressure gradient along s-surfaces 
    485          zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 * r1_e1u(ji,jj)   & 
    486             &           * (  e3w(ji+1,jj,jk,Kmm) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad )   & 
    487             &              - e3w(ji  ,jj,jk,Kmm) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) + 2*znad )  ) 
    488          zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 * r1_e2v(ji,jj)   & 
    489             &           * (  e3w(ji,jj+1,jk,Kmm) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad )   & 
    490             &              - e3w(ji,jj  ,jk,Kmm) * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) + 2*znad )  ) 
    491          ! s-coordinate pressure gradient correction 
    492          zuap = -zcoef0 * ( rhd    (ji+1,jj  ,jk) + rhd    (ji,jj,jk) + 2._wp * znad )   & 
     468      ! 
     469      DO_3D( 0, 0, 0, 0, 2, jpkm1 )    ! interior value (2=<jk=<jpkm1) 
     470         !                                   ! hydrostatic pressure gradient along s-surfaces 
     471         zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 * r1_e1u(ji,jj)                         & 
     472            &           * (  e3w(ji+1,jj,jk,Kmm) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) )  & 
     473            &              - e3w(ji  ,jj,jk,Kmm) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) )  ) 
     474         zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 * r1_e2v(ji,jj)                         & 
     475            &           * (  e3w(ji,jj+1,jk,Kmm) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) )  & 
     476            &              - e3w(ji,jj  ,jk,Kmm) * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) )  ) 
     477         !                                   ! s-coordinate pressure gradient correction 
     478         zuap = -zcoef0 * ( rhd  (ji+1,jj  ,jk) + rhd  (ji,jj,jk) ) & 
    493479            &           * ( gde3w(ji+1,jj  ,jk) - gde3w(ji,jj,jk) ) * r1_e1u(ji,jj) 
    494          zvap = -zcoef0 * ( rhd    (ji  ,jj+1,jk) + rhd    (ji,jj,jk) + 2._wp * znad )  & 
     480         zvap = -zcoef0 * ( rhd  (ji  ,jj+1,jk) + rhd  (ji,jj,jk) ) & 
    495481            &           * ( gde3w(ji  ,jj+1,jk) - gde3w(ji,jj,jk) ) * r1_e2v(ji,jj) 
    496482         ! 
     
    535521      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv    ! ocean velocities and RHS of momentum equation 
    536522      !! 
    537       INTEGER  ::   ji, jj, jk, ikt, iktp1i, iktp1j   ! dummy loop indices 
    538       REAL(wp) ::   zcoef0, zuap, zvap, znad          ! temporary scalars 
     523      INTEGER  ::   ji, jj, jk             ! dummy loop indices 
     524      INTEGER  ::   ikt ,  ikti1,  iktj1   ! local integer 
     525      REAL(wp) ::   ze3w, ze3wi1, ze3wj1   ! local scalars 
     526      REAL(wp) ::   zcoef0, zuap, zvap     !   -      - 
    539527      REAL(wp), DIMENSION(jpi,jpj,jpk ) ::  zhpi, zhpj 
    540528      REAL(wp), DIMENSION(jpi,jpj,jpts) ::  zts_top 
     
    543531      ! 
    544532      zcoef0 = - grav * 0.5_wp   ! Local constant initialization 
    545       ! 
    546       znad=1._wp                 ! To use density and not density anomaly 
    547533      ! 
    548534      !                          ! iniitialised to 0. zhpi zhpi  
     
    560546      CALL eos( zts_top, risfdep, zrhdtop_oce ) 
    561547 
    562 !==================================================================================      
    563 !===== Compute surface value =====================================================  
    564 !================================================================================== 
     548      !                     !===========================! 
     549      !                     !=====  surface value  =====! 
     550      !                     !===========================! 
    565551      DO_2D( 0, 0, 0, 0 ) 
    566          ikt    = mikt(ji,jj) 
    567          iktp1i = mikt(ji+1,jj) 
    568          iktp1j = mikt(ji,jj+1) 
    569          ! hydrostatic pressure gradient along s-surfaces and ice shelf pressure 
    570          ! we assume ISF is in isostatic equilibrium 
    571          zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( 0.5_wp * e3w(ji+1,jj,iktp1i,Kmm)                                    & 
    572             &                                    * ( 2._wp * znad + rhd(ji+1,jj,iktp1i) + zrhdtop_oce(ji+1,jj) )   & 
    573             &                                  - 0.5_wp * e3w(ji,jj,ikt,Kmm)                                         & 
    574             &                                    * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) )          & 
    575             &                                  + ( risfload(ji+1,jj) - risfload(ji,jj))                            )  
    576          zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * e3w(ji,jj+1,iktp1j,Kmm)                                    & 
    577             &                                    * ( 2._wp * znad + rhd(ji,jj+1,iktp1j) + zrhdtop_oce(ji,jj+1) )   & 
    578             &                                  - 0.5_wp * e3w(ji,jj,ikt,Kmm)                                         &  
    579             &                                    * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) )          & 
    580             &                                  + ( risfload(ji,jj+1) - risfload(ji,jj))                            )  
    581          ! s-coordinate pressure gradient correction (=0 if z coordinate) 
    582          zuap = -zcoef0 * ( rhd    (ji+1,jj,1) + rhd    (ji,jj,1) + 2._wp * znad )   & 
     552         ikt   = mikt(ji  ,jj  )   ;   ze3w   = e3w(ji  ,jj  ,ikt  ,Kmm) 
     553         ikti1 = mikt(ji+1,jj  )   ;   ze3wi1 = e3w(ji+1,jj  ,ikti1,Kmm) 
     554         iktj1 = mikt(ji  ,jj+1)   ;   ze3wj1 = e3w(ji  ,jj+1,iktj1,Kmm) 
     555         !                          ! hydrostatic pressure gradient along s-surfaces and ice shelf pressure 
     556         !                          ! we assume ISF is in isostatic equilibrium 
     557         zhpi(ji,jj,1) = zcoef0 * r1_e1u(ji,jj) * (   risfload(ji+1,jj) - risfload(ji,jj)  & 
     558            &                                       + 0.5_wp * ( ze3wi1 * ( rhd(ji+1,jj,ikti1) + zrhdtop_oce(ji+1,jj) )     & 
     559            &                                                  - ze3w   * ( rhd(ji  ,jj,ikt  ) + zrhdtop_oce(ji  ,jj) ) )   ) 
     560         zhpj(ji,jj,1) = zcoef0 * r1_e2v(ji,jj) * (   risfload(ji,jj+1) - risfload(ji,jj)  & 
     561            &                                       + 0.5_wp * ( ze3wj1 * ( rhd(ji,jj+1,iktj1) + zrhdtop_oce(ji,jj+1) )      & 
     562            &                                                  - ze3w   * ( rhd(ji,jj  ,ikt  ) + zrhdtop_oce(ji,jj  ) ) )   )  
     563         !                          ! s-coordinate pressure gradient correction (=0 if z coordinate) 
     564         zuap = -zcoef0 * ( rhd    (ji+1,jj,1) + rhd    (ji,jj,1) )   & 
    583565            &           * ( gde3w(ji+1,jj,1) - gde3w(ji,jj,1) ) * r1_e1u(ji,jj) 
    584          zvap = -zcoef0 * ( rhd    (ji,jj+1,1) + rhd    (ji,jj,1) + 2._wp * znad )   & 
     566         zvap = -zcoef0 * ( rhd    (ji,jj+1,1) + rhd    (ji,jj,1) )   & 
    585567            &           * ( gde3w(ji,jj+1,1) - gde3w(ji,jj,1) ) * r1_e2v(ji,jj) 
    586          ! add to the general momentum trend 
     568         !                          ! add to the general momentum trend 
    587569         puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + (zhpi(ji,jj,1) + zuap) * umask(ji,jj,1) 
    588570         pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + (zhpj(ji,jj,1) + zvap) * vmask(ji,jj,1) 
    589571      END_2D 
    590 !==================================================================================      
    591 !===== Compute interior value =====================================================  
    592 !================================================================================== 
    593       ! interior value (2=<jk=<jpkm1) 
     572      !    
     573      !                     !=============================! 
     574      !                     !=====  interior values  =====! 
     575      !                     !=============================! 
    594576      DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
    595          ! hydrostatic pressure gradient along s-surfaces 
     577         ze3w   = e3w(ji  ,jj  ,jk,Kmm) 
     578         ze3wi1 = e3w(ji+1,jj  ,jk,Kmm) 
     579         ze3wj1 = e3w(ji  ,jj+1,jk,Kmm) 
     580         !                          ! hydrostatic pressure gradient along s-surfaces 
    596581         zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj)   & 
    597             &           * (  e3w(ji+1,jj,jk,Kmm)                   & 
    598             &                  * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) * wmask(ji+1,jj,jk)   & 
    599             &              - e3w(ji  ,jj,jk,Kmm)                   & 
    600             &                  * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) + 2*znad ) * wmask(ji  ,jj,jk)   ) 
     582            &           * (  ze3wi1 * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) ) * wmask(ji+1,jj,jk)   & 
     583            &              - ze3w   * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) ) * wmask(ji  ,jj,jk)   ) 
    601584         zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj)   & 
    602             &           * (  e3w(ji,jj+1,jk,Kmm)                   & 
    603             &                  * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) * wmask(ji,jj+1,jk)   & 
    604             &              - e3w(ji,jj  ,jk,Kmm)                   & 
    605             &                  * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) + 2*znad ) * wmask(ji,jj  ,jk)   ) 
    606          ! s-coordinate pressure gradient correction 
    607          zuap = -zcoef0 * ( rhd   (ji+1,jj  ,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   & 
     585            &           * (  ze3wj1 * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) ) * wmask(ji,jj+1,jk)   & 
     586            &              - ze3w   * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) ) * wmask(ji,jj  ,jk)   ) 
     587         !                          ! s-coordinate pressure gradient correction 
     588         zuap = -zcoef0 * ( rhd   (ji+1,jj  ,jk) + rhd   (ji,jj,jk) )   & 
    608589            &           * ( gde3w(ji+1,jj  ,jk) - gde3w(ji,jj,jk) ) / e1u(ji,jj) 
    609          zvap = -zcoef0 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   & 
     590         zvap = -zcoef0 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji,jj,jk) )   & 
    610591            &           * ( gde3w(ji  ,jj+1,jk) - gde3w(ji,jj,jk) ) / e2v(ji,jj) 
    611          ! add to the general momentum trend 
     592         !                          ! add to the general momentum trend 
    612593         puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk) 
    613594         pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + (zhpj(ji,jj,jk) + zvap) * vmask(ji,jj,jk) 
     
    634615      REAL(wp) ::   z1_12, cffv, cffy   !    "         " 
    635616      LOGICAL  ::   ll_tmp1, ll_tmp2    ! local logical variables 
     617      REAL(wp), DIMENSION(jpi,jpj)     ::   zpgu, zpgv   ! 2D workspace 
    636618      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zhpi, zhpj 
    637619      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   dzx, dzy, dzz, dzu, dzv, dzw 
     
    826808      END_3D 
    827809      CALL lbc_lnk_multi( 'dynhpg', rho_k, 'W', 1.0_wp, rho_i, 'U', 1.0_wp, rho_j, 'V', 1.0_wp ) 
    828  
     810      ! 
     811      ! --------------- 
     812      !  Surface pressure gradient to be removed 
     813      ! --------------- 
     814      DO_2D( 0, 0, 0, 0 ) 
     815         zpgu(ji,jj) = - grav * ( ssh(ji+1,jj,Kmm) - ssh(ji,jj,Kmm) ) * r1_e1u(ji,jj) 
     816         zpgv(ji,jj) = - grav * ( ssh(ji,jj+1,Kmm) - ssh(ji,jj,Kmm) ) * r1_e2v(ji,jj) 
     817      END_2D 
     818      ! 
    829819      ! --------------- 
    830820      !  Surface value 
     
    838828         ENDIF 
    839829         ! add to the general momentum trend 
    840          puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + zhpi(ji,jj,1) 
    841          pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + zhpj(ji,jj,1) 
     830         puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + zhpi(ji,jj,1) - zpgu(ji,jj) 
     831         pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + zhpj(ji,jj,1) - zpgv(ji,jj) 
    842832      END_2D 
    843833 
     
    858848         ENDIF 
    859849         ! add to the general momentum trend 
    860          puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zhpi(ji,jj,jk) 
    861          pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zhpj(ji,jj,jk) 
     850         puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zhpi(ji,jj,jk) - zpgu(ji,jj) 
     851         pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zhpj(ji,jj,jk) - zpgv(ji,jj) 
    862852      END_3D 
    863853      ! 
     
    892882      REAL(wp) :: zrhdt1 
    893883      REAL(wp) :: zdpdx1, zdpdx2, zdpdy1, zdpdy2 
     884      REAL(wp), DIMENSION(jpi,jpj)     ::   zpgu, zpgv   ! 2D workspace 
     885      REAL(wp), DIMENSION(jpi,jpj)     ::   zsshu_n, zsshv_n 
    894886      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdept, zrhh 
    895887      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp 
    896       REAL(wp), DIMENSION(jpi,jpj)   ::   zsshu_n, zsshv_n 
    897888      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zcpx, zcpy   !W/D pressure filter 
    898889      !!---------------------------------------------------------------------- 
     
    907898      zcoef0 = - grav 
    908899      znad = 1._wp 
    909       IF( ln_linssh )   znad = 0._wp 
    910  
     900      IF( ln_linssh )   znad = 1._wp 
     901      ! 
     902      ! --------------- 
     903      !  Surface pressure gradient to be removed 
     904      ! --------------- 
     905      DO_2D( 0, 0, 0, 0 ) 
     906         zpgu(ji,jj) = - grav * ( ssh(ji+1,jj,Kmm) - ssh(ji,jj,Kmm) ) * r1_e1u(ji,jj) 
     907         zpgv(ji,jj) = - grav * ( ssh(ji,jj+1,Kmm) - ssh(ji,jj,Kmm) ) * r1_e2v(ji,jj) 
     908      END_2D 
     909      ! 
    911910      IF( ln_wd_il ) THEN 
    912911         ALLOCATE( zcpx(jpi,jpj) , zcpy(jpi,jpj) ) 
     
    974973      ! Transfer the depth of "T(:,:,:)" to vertical coordinate "zdept(:,:,:)" 
    975974      DO_2D( 1, 1, 1, 1 ) 
    976          zdept(ji,jj,1) = 0.5_wp * e3w(ji,jj,1,Kmm) - ssh(ji,jj,Kmm) * znad 
     975         zdept(ji,jj,1) = 0.5_wp * e3w(ji,jj,1,Kmm) - ssh(ji,jj,Kmm) 
    977976      END_2D 
    978977 
     
    10251024 
    10261025      DO_2D( 0, 0, 0, 0 ) 
    1027        zu(ji,jj,1) = - ( e3u(ji,jj,1,Kmm) - zsshu_n(ji,jj) * znad)  
    1028        zv(ji,jj,1) = - ( e3v(ji,jj,1,Kmm) - zsshv_n(ji,jj) * znad) 
     1026       zu(ji,jj,1) = - ( e3u(ji,jj,1,Kmm) - zsshu_n(ji,jj) )  
     1027       zv(ji,jj,1) = - ( e3v(ji,jj,1,Kmm) - zsshv_n(ji,jj) ) 
    10291028      END_2D 
    10301029 
     
    11081107            zdpdx2 = zdpdx2 * zcpx(ji,jj) * wdrampu(ji,jj) 
    11091108         ENDIF 
    1110          puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk)  
     1109         puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (zdpdx1 + zdpdx2 - zpgu(ji,jj)) * umask(ji,jj,jk)  
    11111110      ENDIF 
    11121111 
     
    11681167         ENDIF 
    11691168 
    1170          pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + (zdpdy1 + zdpdy2) * vmask(ji,jj,jk) 
     1169         pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + (zdpdy1 + zdpdy2 - zpgv(ji,jj)) * vmask(ji,jj,jk) 
    11711170      ENDIF 
    11721171         ! 
Note: See TracChangeset for help on using the changeset viewer.