Changeset 5864


Ignore:
Timestamp:
2015-11-05T17:38:47+01:00 (5 years ago)
Author:
gm
Message:

#1613: vvl by default: use r1_e1 instead of 1 / e1 in dynhpg.F90

File:
1 edited

Legend:

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

    r5845 r5864  
    234234            zcoef1 = zcoef0 * e3w_n(ji,jj,1) 
    235235            ! hydrostatic pressure gradient 
    236             zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) / e1u(ji,jj) 
    237             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) 
     236            zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj) 
     237            zhpj(ji,jj,1) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj) 
    240238            ! add to the general momentum trend 
    241239            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) 
     
    253251               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1)   & 
    254252                  &           + zcoef1 * (  ( rhd(ji+1,jj,jk)+rhd(ji+1,jj,jk-1) )   & 
    255                   &                       - ( 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) 
     253                  &                       - ( rhd(ji  ,jj,jk)+rhd(ji  ,jj,jk-1) )  ) * r1_e1u(ji,jj) 
    257254 
    258255               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1)   & 
    259256                  &           + zcoef1 * (  ( rhd(ji,jj+1,jk)+rhd(ji,jj+1,jk-1) )   & 
    260                   &                       - ( 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) 
     257                  &                       - ( rhd(ji,jj,  jk)+rhd(ji,jj  ,jk-1) )  ) * r1_e2v(ji,jj) 
    262258               ! add to the general momentum trend 
    263259               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) 
     
    305301            zcoef1 = zcoef0 * e3w_n(ji,jj,1) 
    306302            ! hydrostatic pressure gradient 
    307             zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj  ,1) - rhd(ji,jj,1) ) / e1u(ji,jj) 
    308             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) 
     303            zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj  ,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj) 
     304            zhpj(ji,jj,1) = zcoef1 * ( rhd(ji  ,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj) 
    311305            ! add to the general momentum trend 
    312306            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) 
     
    324318               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1)   & 
    325319                  &           + zcoef1 * (  ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) )   & 
    326                   &                       - ( 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) 
     320                  &                       - ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) )  ) * r1_e1u(ji,jj) 
    328321 
    329322               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1)   & 
    330323                  &           + zcoef1 * (  ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) )   & 
    331                   &                       - ( 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) 
     324                  &                       - ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) )  ) * r1_e2v(ji,jj) 
    333325               ! add to the general momentum trend 
    334326               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) 
     
    349341               ua  (ji,jj,iku) = ua(ji,jj,iku) - zhpi(ji,jj,iku)         ! subtract old value 
    350342               zhpi(ji,jj,iku) = zhpi(ji,jj,iku-1)                   &   ! compute the new one 
    351                   &            + 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) 
     343                  &            + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + gru(ji,jj) ) * r1_e1u(ji,jj) 
    353344               ua  (ji,jj,iku) = ua(ji,jj,iku) + zhpi(ji,jj,iku)         ! add the new one to the general momentum trend 
    354345            ENDIF 
     
    356347               va  (ji,jj,ikv) = va(ji,jj,ikv) - zhpj(ji,jj,ikv)         ! subtract old value 
    357348               zhpj(ji,jj,ikv) = zhpj(ji,jj,ikv-1)                   &   ! compute the new one 
    358                   &            + 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) 
     349                  &            + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + grv(ji,jj) ) * r1_e2v(ji,jj) 
    360350               va  (ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv)         ! add the new one to the general momentum trend 
    361351            ENDIF 
     
    411401         DO ji = fs_2, fs_jpim1   ! vector opt. 
    412402            ! hydrostatic pressure gradient along s-surfaces 
    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) ) ) 
     403            zhpi(ji,jj,1) = zcoef0 * r1_e1u(ji,jj) * (  e3w_n(ji+1,jj  ,1) * ( znad + rhd(ji+1,jj  ,1) )  & 
     404               &                                      - e3w_n(ji  ,jj  ,1) * ( znad + rhd(ji  ,jj  ,1) )  ) 
     405            zhpj(ji,jj,1) = zcoef0 * r1_e2v(ji,jj) * (  e3w_n(ji  ,jj+1,1) * ( znad + rhd(ji  ,jj+1,1) )  & 
     406               &                                      - e3w_n(ji  ,jj  ,1) * ( znad + rhd(ji  ,jj  ,1) )  ) 
    419407            ! s-coordinate pressure gradient correction 
    420408            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) 
     409               &           * ( gde3w_n(ji+1,jj,1) - gde3w_n(ji,jj,1) ) * r1_e1u(ji,jj) 
    423410            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) 
     411               &           * ( gde3w_n(ji,jj+1,1) - gde3w_n(ji,jj,1) ) * r1_e2v(ji,jj) 
    426412            ! add to the general momentum trend 
    427413            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap 
     
    435421            DO ji = fs_2, fs_jpim1   ! vector opt. 
    436422               ! hydrostatic pressure gradient along s-surfaces 
    437 !!gm               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 * r1_e1u(ji,jj)   & 
    438                zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj)   & 
     423               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 * r1_e1u(ji,jj)   & 
    439424                  &           * (  e3w_n(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad )   & 
    440425                  &              - 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)   & 
    442426               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 * r1_e2v(ji,jj)   & 
    443427                  &           * (  e3w_n(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad )   & 
     
    445429               ! s-coordinate pressure gradient correction 
    446430               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) 
     431                  &           * ( gde3w_n(ji+1,jj  ,jk) - gde3w_n(ji,jj,jk) ) * r1_e1u(ji,jj) 
    449432               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) 
     433                  &           * ( gde3w_n(ji  ,jj+1,jk) - gde3w_n(ji,jj,jk) ) * r1_e2v(ji,jj) 
    452434               ! add to the general momentum trend 
    453435               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap 
     
    561543            ! hydrostatic pressure gradient along s-surfaces and ice shelf pressure 
    562544            ! we assume ISF is in isostatic equilibrium 
    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)                                    & 
    565                &                                  * ( 2._wp * znad + rhd(ji+1,jj  ,iktp1i) + zrhdtop_oce(ji+1,jj  ) )   & 
    566                &                                  - 0.5_wp * e3w_n(ji  ,jj  ,ikt   )                                    & 
    567                &                                  * ( 2._wp * znad + rhd(ji  ,jj  ,ikt   ) + zrhdtop_oce(ji  ,jj  ) )   & 
    568                &                                  + ( ziceload(ji+1,jj) - ziceload(ji,jj))                              )  
    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)                                    & 
    571                &                                  * ( 2._wp * znad + rhd(ji  ,jj+1,iktp1j) + zrhdtop_oce(ji  ,jj+1) )   & 
    572                &                                  - 0.5_wp * e3w_n(ji  ,jj  ,ikt   )                                    &  
    573                &                                  * ( 2._wp * znad + rhd(ji  ,jj  ,ikt   ) + zrhdtop_oce(ji  ,jj  ) )   & 
    574                &                                  + ( ziceload(ji,jj+1) - ziceload(ji,jj) )                             )  
     545            zhpi(ji,jj,1) = zcoef0 * r1_e1u(ji,jj) * ( 0.5_wp * e3w_n(ji+1,jj  ,iktp1i)                                  & 
     546               &                                   * ( 2._wp * znad + rhd(ji+1,jj  ,iktp1i) + zrhdtop_oce(ji+1,jj  ) )   & 
     547               &                                   - 0.5_wp * e3w_n(ji  ,jj  ,ikt   )                                    & 
     548               &                                   * ( 2._wp * znad + rhd(ji  ,jj  ,ikt   ) + zrhdtop_oce(ji  ,jj  ) )   & 
     549               &                                   + ( ziceload(ji+1,jj) - ziceload(ji,jj) )                             )  
     550            zhpj(ji,jj,1) = zcoef0 * r1_e2v(ji,jj) * ( 0.5_wp * e3w_n(ji  ,jj+1,iktp1j)                                  & 
     551               &                                   * ( 2._wp * znad + rhd(ji  ,jj+1,iktp1j) + zrhdtop_oce(ji  ,jj+1) )   & 
     552               &                                   - 0.5_wp * e3w_n(ji  ,jj  ,ikt   )                                    &  
     553               &                                   * ( 2._wp * znad + rhd(ji  ,jj  ,ikt   ) + zrhdtop_oce(ji  ,jj  ) )   & 
     554               &                                   + ( ziceload(ji,jj+1) - ziceload(ji,jj) )                             )  
    575555            ! s-coordinate pressure gradient correction (=0 if z coordinate) 
    576556            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) 
     557               &           * ( gde3w_n(ji+1,jj,1) - gde3w_n(ji,jj,1) ) * r1_e1u(ji,jj) 
    579558            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) 
     559               &           * ( gde3w_n(ji,jj+1,1) - gde3w_n(ji,jj,1) ) * r1_e2v(ji,jj) 
    582560            ! add to the general momentum trend 
    583561            ua(ji,jj,1) = ua(ji,jj,1) + (zhpi(ji,jj,1) + zuap) * umask(ji,jj,1) 
     
    596574            IF ( iku .GT. 1 ) THEN 
    597575               ! case iku 
    598 !!gm               zhpi(ji,jj,iku)   =  zcoef0 * r1_e1u(ji,jj) * ze3wu                                            & 
    599                zhpi(ji,jj,iku)   =  zcoef0 / e1u(ji,jj) * ze3wu                                            & 
     576               zhpi(ji,jj,iku)   =  zcoef0 * r1_e1u(ji,jj) * ze3wu                                         & 
    600577                      &                                 * ( rhd    (ji+1,jj,iku) + rhd   (ji,jj,iku)       & 
    601578                      &                                   + SIGN(1._wp,ze3wu) * grui(ji,jj) + 2._wp * znad ) 
    602579               ! corrective term ( = 0 if z coordinate ) 
    603 !!gm               zuap              = -zcoef0 * ( arui(ji,jj) + 2._wp * znad ) * gzui(ji,jj) * r1_e1u(ji,jj) 
    604                zuap              = -zcoef0 * ( arui(ji,jj) + 2._wp * znad ) * gzui(ji,jj) / e1u(ji,jj) 
     580               zuap              = -zcoef0 * ( arui(ji,jj) + 2._wp * znad ) * gzui(ji,jj) * r1_e1u(ji,jj) 
    605581               ! zhpi will be added in interior loop 
    606582               ua(ji,jj,iku)     = ua(ji,jj,iku) + zuap 
     
    609585 
    610586               ! 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)                                                               &     
    612                zhpiint        =  zcoef0 / e1u(ji,jj)                                                               &     
     587               zhpiint        =  zcoef0 * r1_e1u(ji,jj)                                                               &     
    613588                  &           * (  e3w_n(ji+1,jj  ,iku+1) * ( (rhd(ji+1,jj,iku+1) + znad)                          & 
    614589                  &                                         + (rhd(ji+1,jj,iku  ) + znad) ) * tmask(ji+1,jj,iku)   & 
    615590                  &              - e3w_n(ji  ,jj  ,iku+1) * ( (rhd(ji  ,jj,iku+1) + znad)                          & 
    616591                  &                                         + (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  
    618                zhpi(ji,jj,iku+1) =  zcoef0 / e1u(ji,jj) * ge3rui(ji,jj) - zhpiint  
     592               zhpi(ji,jj,iku+1) =  zcoef0 * r1_e1u(ji,jj) * ge3rui(ji,jj) - zhpiint  
    619593            END IF 
    620594                
     
    624598            IF ( ikv .GT. 1 ) THEN 
    625599               ! case ikv 
    626 !!gm               zhpj(ji,jj,ikv)   =  zcoef0 * r1_e2v(ji,jj) * ze3wv                                            & 
    627                zhpj(ji,jj,ikv)   =  zcoef0 / e2v(ji,jj) * ze3wv                                            & 
    628                      &                                  * ( rhd(ji,jj+1,ikv) + rhd   (ji,jj,ikv)           & 
    629                      &                                    + SIGN(1._wp,ze3wv) * grvi(ji,jj) + 2._wp * znad ) 
     600               zhpj(ji,jj,ikv)   =  zcoef0 * r1_e2v(ji,jj) * ze3wv                                            & 
     601                     &                                     * ( rhd(ji,jj+1,ikv) + rhd   (ji,jj,ikv)           & 
     602                     &                                       + SIGN(1._wp,ze3wv) * grvi(ji,jj) + 2._wp * znad ) 
    630603               ! corrective term ( = 0 if z coordinate ) 
    631 !!gm               zvap              = -zcoef0 * ( arvi(ji,jj) + 2._wp * znad ) * gzvi(ji,jj) * r1_e2v(ji,jj) 
    632                zvap              = -zcoef0 * ( arvi(ji,jj) + 2._wp * znad ) * gzvi(ji,jj) / e2v(ji,jj) 
     604               zvap              = -zcoef0 * ( arvi(ji,jj) + 2._wp * znad ) * gzvi(ji,jj) * r1_e2v(ji,jj) 
    633605               ! zhpi will be added in interior loop 
    634606               va(ji,jj,ikv)      = va(ji,jj,ikv) + zvap 
     
    637609                
    638610               ! 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)                                                              & 
    640                zhpjint        =  zcoef0 / e2v(ji,jj)                                                              & 
     611               zhpjint        =  zcoef0 * r1_e2v(ji,jj)                                                           & 
    641612                  &           * (  e3w_n(ji  ,jj+1,ikv+1) * ( (rhd(ji,jj+1,ikv+1) + znad)                         & 
    642613                  &                                       + (rhd(ji,jj+1,ikv  ) + znad) ) * tmask(ji,jj+1,ikv)    & 
    643614                  &              - e3w_n(ji  ,jj  ,ikv+1) * ( (rhd(ji,jj  ,ikv+1) + znad)                         & 
    644615                  &                                       + (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 
    646                zhpj(ji,jj,ikv+1) =  zcoef0 / e2v(ji,jj) * ge3rvi(ji,jj) - zhpjint 
     616               zhpj(ji,jj,ikv+1) =  zcoef0 * r1_e2v(ji,jj) * ge3rvi(ji,jj) - zhpjint 
    647617            END IF 
    648618         END DO 
     
    655625      DO jj = 2, jpjm1 
    656626         DO ji = fs_2, fs_jpim1   ! vector opt. 
    657 !!gm useles ! 
    658             iku=miku(ji,jj); ikv=mikv(ji,jj) 
    659 !!gm 
    660627            DO jk = 2, jpkm1 
    661628               ! hydrostatic pressure gradient along s-surfaces 
    662629               ! zhpi is masked for the first wet cell  (contribution already done in the upper bloc) 
    663630               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)                                                                & 
     631                  &           + zcoef0 * r1_e1u(ji,jj)                                                                & 
    666632                  &                    * ( e3w_n(ji+1,jj,jk) * ( (rhd(ji+1,jj,jk  ) + znad)                           & 
    667633                  &                                            + (rhd(ji+1,jj,jk-1) + znad) ) * tmask(ji+1,jj,jk-1)   & 
     
    671637               ! corrective term, we mask this term for the first wet level beneath the ice shelf (contribution done in the upper bloc)  
    672638               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) 
     639                  &            * ( gde3w_n(ji+1,jj  ,jk) - gde3w_n(ji,jj,jk) ) * r1_e1u(ji,jj) * umask(ji,jj,jk-1) 
    675640               ua(ji,jj,jk) = ua(ji,jj,jk) + ( zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk) 
    676641 
     
    678643               ! zhpi is masked for the first wet cell  (contribution already done in the upper bloc) 
    679644               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)                                                                & 
     645                  &           + zcoef0 * r1_e2v(ji,jj)                                                                & 
    682646                  &                    * ( e3w_n(ji  ,jj+1,jk) * ( (rhd(ji,jj+1,jk  ) + znad)                           & 
    683647                  &                                              + (rhd(ji,jj+1,jk-1) + znad) ) * tmask(ji,jj+1,jk-1)   & 
     
    687651               ! corrective term, we mask this term for the first wet level beneath the ice shelf (contribution done in the upper bloc) 
    688652               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) 
     653                  &            * ( gde3w_n(ji  ,jj+1,jk) - gde3w_n(ji,jj,jk) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk-1) 
    691654               ! add to the general momentum trend 
    692655               va(ji,jj,jk) = va(ji,jj,jk) + ( zhpj(ji,jj,jk) + zvap ) * vmask(ji,jj,jk) 
     
    707670               ! remove old value (interior case) 
    708671               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) 
     672                     &                   * ( gde3w_n(ji+1,jj  ,iku) - gde3w_n(ji,jj,iku) ) * r1_e1u(ji,jj) 
    711673               ua(ji,jj,iku)   = ua(ji,jj,iku) - zhpi(ji,jj,iku) - zuap 
    712674               ! put new value 
    713675               ! -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) 
    715                zuap            = -zcoef0 * ( aru(ji,jj) + 2._wp * znad ) * gzu(ji,jj)  / e1u(ji,jj) 
    716                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)  
     676               zuap            = -zcoef0 * ( aru(ji,jj) + 2._wp * znad ) * gzu(ji,jj)  * r1_e1u(ji,jj) 
     677               zhpi(ji,jj,iku) =  zhpi(ji,jj,iku-1) + zcoef0 * r1_e1u(ji,jj) * ge3ru(ji,jj) - zpshpi(ji,jj)  
    718678               ua(ji,jj,iku)   =  ua(ji,jj,iku) + zhpi(ji,jj,iku) + zuap 
    719679            END IF 
     
    722682               ! remove old value (interior case) 
    723683               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) 
     684                     &                   * ( gde3w_n(ji  ,jj+1,ikv) - gde3w_n(ji,jj,ikv) ) * r1_e2v(ji,jj) 
    726685               va(ji,jj,ikv)   = va(ji,jj,ikv) - zhpj(ji,jj,ikv) - zvap 
    727686               ! put new value 
    728687               ! -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) 
    730                zvap            = -zcoef0 * ( arv(ji,jj) + 2._wp * znad ) * gzv(ji,jj)     / e2v(ji,jj) 
    731                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)    
     688               zvap            = -zcoef0 * ( arv(ji,jj) + 2._wp * znad ) * gzv(ji,jj) * r1_e2v(ji,jj) 
     689               zhpj(ji,jj,ikv) =  zhpj(ji,jj,ikv-1) + zcoef0 * r1_e2v(ji,jj) * ge3rv(ji,jj) - zpshpj(ji,jj)    
    733690               va(ji,jj,ikv)   =  va(ji,jj,ikv) + zhpj(ji,jj,ikv) + zvap 
    734691            END IF 
     
    938895      DO jj = 2, jpjm1 
    939896         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) 
    941             zhpi(ji,jj,1) = ( rho_k(ji+1,jj  ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) / e1u(ji,jj) 
    942             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) 
     897            zhpi(ji,jj,1) = ( rho_k(ji+1,jj  ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) * r1_e1u(ji,jj) 
     898            zhpj(ji,jj,1) = ( rho_k(ji  ,jj+1,1) - rho_k(ji,jj,1) - rho_j(ji,jj,1) ) * r1_e2v(ji,jj) 
    944899            ! add to the general momentum trend 
    945900            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) 
     
    957912               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1)                                & 
    958913                  &           + (  ( rho_k(ji+1,jj,jk) - rho_k(ji,jj,jk  ) )    & 
    959                   &              - ( 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) 
     914                  &              - ( rho_i(ji  ,jj,jk) - rho_i(ji,jj,jk-1) )  ) * r1_e1u(ji,jj) 
    961915               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1)                                & 
    962916                  &           + (  ( rho_k(ji,jj+1,jk) - rho_k(ji,jj,jk  ) )    & 
    963                   &               -( 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) 
     917                  &               -( rho_j(ji,jj  ,jk) - rho_j(ji,jj,jk-1) )  ) * r1_e2v(ji,jj) 
    965918               ! add to the general momentum trend 
    966919               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) 
     
    992945      !! 
    993946      INTEGER  ::   ji, jj, jk, jkk                 ! dummy loop indices 
    994       REAL(wp) ::   zcoef0, znad                    ! temporary scalars 
    995       !! 
     947      REAL(wp) ::   zcoef0, znad                    ! local scalars 
     948      ! 
    996949      !! The local variables for the correction term 
    997950      INTEGER  :: jk1, jis, jid, jjs, jjd 
     
    1004957      !!---------------------------------------------------------------------- 
    1005958      ! 
    1006       CALL wrk_alloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp ) 
    1007       CALL wrk_alloc( jpi,jpj,jpk, zdept, zrhh ) 
    1008       CALL wrk_alloc( jpi,jpj, zsshu_n, zsshv_n ) 
     959      CALL wrk_alloc( jpi,jpj,jpk,   zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp ) 
     960      CALL wrk_alloc( jpi,jpj,jpk,   zdept, zrhh ) 
     961      CALL wrk_alloc( jpi,jpj,       zsshu_n, zsshv_n ) 
    1009962      ! 
    1010963      IF( kt == nit000 ) THEN 
     
    1014967      ENDIF 
    1015968 
    1016       !!---------------------------------------------------------------------- 
    1017969      ! Local constant initialization 
    1018970      zcoef0 = - grav 
     
    10601012      ! constrained cubic spline interpolation 
    10611013      ! rho(z) = asp + bsp*z + csp*z^2 + dsp*z^3 
    1062       CALL cspline(fsp,xsp,asp,bsp,csp,dsp,polynomial_type) 
     1014      CALL cspline( fsp, xsp, asp, bsp, csp, dsp, polynomial_type ) 
    10631015 
    10641016      ! Integrate the hydrostatic pressure "zhpi(:,:,:)" at "T(ji,jj,1)" 
    10651017      DO jj = 2, jpj 
    10661018        DO ji = 2, jpi 
    1067           zrhdt1 = zrhh(ji,jj,1) - interp3(zdept(ji,jj,1),asp(ji,jj,1), & 
    1068                                          bsp(ji,jj,1),   csp(ji,jj,1), & 
    1069                                          dsp(ji,jj,1) ) * 0.25_wp * e3w_n(ji,jj,1) 
     1019          zrhdt1 = zrhh(ji,jj,1) - interp3( zdept(ji,jj,1), asp(ji,jj,1), bsp(ji,jj,1),  & 
     1020             &                                              csp(ji,jj,1), dsp(ji,jj,1) ) * 0.25_wp * e3w_n(ji,jj,1) 
    10701021 
    10711022          ! assuming linear profile across the top half surface layer 
     
    10781029        DO jj = 2, jpj 
    10791030          DO ji = 2, jpi 
    1080             zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) +                          & 
    1081                              integ_spline(zdept(ji,jj,jk-1), zdept(ji,jj,jk),& 
    1082                                     asp(ji,jj,jk-1),    bsp(ji,jj,jk-1), & 
    1083                                     csp(ji,jj,jk-1),    dsp(ji,jj,jk-1)) 
     1031            zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) +                                  & 
     1032               &             integ_spline( zdept(ji,jj,jk-1), zdept(ji,jj,jk),   & 
     1033               &                           asp  (ji,jj,jk-1), bsp  (ji,jj,jk-1), & 
     1034               &                           csp  (ji,jj,jk-1), dsp  (ji,jj,jk-1)  ) 
    10841035          END DO 
    10851036        END DO 
     
    10911042      DO jj = 2, jpjm1 
    10921043        DO ji = 2, jpim1 
     1044!!gm BUG ?    if it is ssh at u- & v-point then it should be: 
     1045!          zsshu_n(ji,jj) = (e1e2t(ji,jj) * sshn(ji,jj) + e1e2t(ji+1,jj) * sshn(ji+1,jj)) * & 
     1046!                         & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp  
     1047!          zsshv_n(ji,jj) = (e1e2t(ji,jj) * sshn(ji,jj) + e1e2t(ji,jj+1) * sshn(ji,jj+1)) * & 
     1048!                         & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp  
     1049!!gm not this: 
    10931050          zsshu_n(ji,jj) = (e1e2u(ji,jj) * sshn(ji,jj) + e1e2u(ji+1, jj) * sshn(ji+1,jj)) * & 
    10941051                         & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp  
     
    11081065        DO jj = 2, jpjm1 
    11091066          DO ji = 2, jpim1 
    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) 
     1067            zu(ji,jj,jk) = zu(ji,jj,jk-1) - e3u_n(ji,jj,jk) 
     1068            zv(ji,jj,jk) = zv(ji,jj,jk-1) - e3v_n(ji,jj,jk) 
    11121069          END DO 
    11131070        END DO 
     
    11261083        DO jj = 2, jpjm1 
    11271084          DO ji = 2, jpim1 
    1128             zu(ji,jj,jk) = min(zu(ji,jj,jk), max(-zdept(ji,jj,jk), -zdept(ji+1,jj,jk))) 
    1129             zu(ji,jj,jk) = max(zu(ji,jj,jk), min(-zdept(ji,jj,jk), -zdept(ji+1,jj,jk))) 
    1130             zv(ji,jj,jk) = min(zv(ji,jj,jk), max(-zdept(ji,jj,jk), -zdept(ji,jj+1,jk))) 
    1131             zv(ji,jj,jk) = max(zv(ji,jj,jk), min(-zdept(ji,jj,jk), -zdept(ji,jj+1,jk))) 
     1085            zu(ji,jj,jk) = MIN(  zu(ji,jj,jk) , MAX( -zdept(ji,jj,jk) , -zdept(ji+1,jj,jk) )  ) 
     1086            zu(ji,jj,jk) = MAX(  zu(ji,jj,jk) , MIN( -zdept(ji,jj,jk) , -zdept(ji+1,jj,jk) )  ) 
     1087            zv(ji,jj,jk) = MIN(  zv(ji,jj,jk) , MAX( -zdept(ji,jj,jk) , -zdept(ji,jj+1,jk) )  ) 
     1088            zv(ji,jj,jk) = MAX(  zv(ji,jj,jk) , MIN( -zdept(ji,jj,jk) , -zdept(ji,jj+1,jk) )  ) 
    11321089          END DO 
    11331090        END DO 
     
    11871144               ! update the momentum trends in u direction 
    11881145 
    1189 !!gm               zdpdx1 = zcoef0 * r1_e1u(ji,jj) * (zhpi(ji+1,jj,jk) - zhpi(ji,jj,jk)) 
    1190                zdpdx1 = zcoef0 / e1u(ji,jj) * (zhpi(ji+1,jj,jk) - zhpi(ji,jj,jk)) 
     1146               zdpdx1 = zcoef0 * r1_e1u(ji,jj) * ( zhpi(ji+1,jj,jk) - zhpi(ji,jj,jk) ) 
    11911147               IF( lk_vvl ) THEN 
    1192 !!gm                 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * & 
    1193                  zdpdx2 = zcoef0 / e1u(ji,jj) * & 
    1194                          ( REAL(jis-jid, wp) * (zpwes + zpwed) + (sshn(ji+1,jj)-sshn(ji,jj)) ) 
     1148                 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * & 
     1149                    &    ( REAL(jis-jid, wp) * (zpwes + zpwed) + (sshn(ji+1,jj)-sshn(ji,jj)) ) 
    11951150                ELSE 
    1196 !!gm                 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed) 
    1197                  zdpdx2 = zcoef0 / e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed) 
     1151                 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed) 
    11981152               ENDIF 
    1199  
    1200                ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * & 
    1201                &          umask(ji,jj,jk) * tmask(ji,jj,jk) * tmask(ji+1,jj,jk) 
     1153!!gm  Since umask(ji,:,:) = tmask(ji,:,:) * tmask(ji+1,:,:)  by definition 
     1154!!gm      in the line below only * umask(ji,jj,jk)  is needed !! 
     1155               ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk) * tmask(ji,jj,jk) * tmask(ji+1,jj,jk) 
    12021156            ENDIF 
    12031157 
     
    12471201               ! update the momentum trends in v direction 
    12481202 
    1249 !!gm               zdpdy1 = zcoef0 * r1_e2v(ji,jj) * (zhpi(ji,jj+1,jk) - zhpi(ji,jj,jk)) 
    1250                zdpdy1 = zcoef0 / e2v(ji,jj) * (zhpi(ji,jj+1,jk) - zhpi(ji,jj,jk)) 
     1203               zdpdy1 = zcoef0 * r1_e2v(ji,jj) * ( zhpi(ji,jj+1,jk) - zhpi(ji,jj,jk) ) 
    12511204               IF( lk_vvl ) THEN 
    1252 !!gm                   zdpdy2 = zcoef0 * r1_e2v(ji,jj) * & 
    1253                    zdpdy2 = zcoef0 / e2v(ji,jj) * & 
     1205                   zdpdy2 = zcoef0 * r1_e2v(ji,jj) * & 
    12541206                           ( REAL(jjs-jjd, wp) * (zpnss + zpnsd) + (sshn(ji,jj+1)-sshn(ji,jj)) ) 
    12551207               ELSE 
    1256 !!gm                   zdpdy2 = zcoef0 * r1_e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd ) 
    1257                    zdpdy2 = zcoef0 / e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd ) 
     1208                   zdpdy2 = zcoef0 * r1_e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd ) 
    12581209               ENDIF 
    1259  
    1260                va(ji,jj,jk) = va(ji,jj,jk) + (zdpdy1 + zdpdy2)*& 
    1261                &              vmask(ji,jj,jk)*tmask(ji,jj,jk)*tmask(ji,jj+1,jk) 
     1210!!gm  Since vmask(:,jj,:) = tmask(:,jj,:) * tmask(:,jj+1,:)  by definition 
     1211!!gm      in the line below only * vmask(ji,jj,jk)  is needed !! 
     1212               va(ji,jj,jk) = va(ji,jj,jk) + ( zdpdy1 + zdpdy2 ) * vmask(ji,jj,jk)*tmask(ji,jj,jk)*tmask(ji,jj+1,jk) 
    12621213            ENDIF 
    1263  
    1264  
    1265            END DO 
    1266         END DO 
    1267       END DO 
    1268       ! 
    1269       CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp ) 
    1270       CALL wrk_dealloc( jpi,jpj,jpk, zdept, zrhh ) 
    1271       CALL wrk_dealloc( jpi,jpj, zsshu_n, zsshv_n ) 
     1214               ! 
     1215            END DO 
     1216         END DO 
     1217      END DO 
     1218      ! 
     1219      CALL wrk_dealloc( jpi,jpj,jpk,   zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp ) 
     1220      CALL wrk_dealloc( jpi,jpj,jpk,   zdept, zrhh ) 
     1221      CALL wrk_dealloc( jpi,jpj,       zsshu_n, zsshv_n ) 
    12721222      ! 
    12731223   END SUBROUTINE hpg_prj 
    12741224 
    12751225 
    1276    SUBROUTINE cspline(fsp, xsp, asp, bsp, csp, dsp, polynomial_type) 
     1226   SUBROUTINE cspline( fsp, xsp, asp, bsp, csp, dsp, polynomial_type ) 
    12771227      !!---------------------------------------------------------------------- 
    12781228      !!                 ***  ROUTINE cspline  *** 
     
    12971247      REAL(wp) ::   zdf(size(fsp,3)) 
    12981248      !!---------------------------------------------------------------------- 
    1299  
     1249      ! 
     1250!!gm  WHAT !!!!!   THIS IS VERY DANGEROUS !!!!!   
    13001251      jpi   = size(fsp,1) 
    13011252      jpj   = size(fsp,2) 
    13021253      jpkm1 = size(fsp,3) - 1 
    1303  
    1304  
     1254      ! 
    13051255      IF (polynomial_type == 1) THEN     ! Constrained Cubic Spline 
    13061256         DO ji = 1, jpi 
     
    13351285 
    13361286               zdf(1)     = 1.5_wp * ( fsp(ji,jj,2) - fsp(ji,jj,1) ) / & 
    1337                           &          ( xsp(ji,jj,2) - xsp(ji,jj,1) ) -  0.5_wp * zdf(2) 
     1287                          &          ( xsp(ji,jj,2) - xsp(ji,jj,1) )           -  0.5_wp * zdf(2) 
    13381288               zdf(jpkm1) = 1.5_wp * ( fsp(ji,jj,jpkm1) - fsp(ji,jj,jpkm1-1) ) / & 
    1339                           &          ( xsp(ji,jj,jpkm1) - xsp(ji,jj,jpkm1-1) ) - & 
    1340                           & 0.5_wp * zdf(jpkm1 - 1) 
     1289                          &          ( xsp(ji,jj,jpkm1) - xsp(ji,jj,jpkm1-1) ) - 0.5_wp * zdf(jpkm1 - 1) 
    13411290 
    13421291               DO jk = 1, jpkm1 - 1 
     
    13611310         END DO 
    13621311 
    1363       ELSE IF (polynomial_type == 2) THEN     ! Linear 
     1312      ELSEIF ( polynomial_type == 2 ) THEN     ! Linear 
    13641313         DO ji = 1, jpi 
    13651314            DO jj = 1, jpj 
     
    13921341      !!                extrapolation is also permitted (no value limit) 
    13931342      !!---------------------------------------------------------------------- 
    1394       IMPLICIT NONE 
    13951343      REAL(wp), INTENT(in) ::  x, xl, xr, fl, fr 
    13961344      REAL(wp)             ::  f ! result of the interpolation (extrapolation) 
    13971345      REAL(wp)             ::  zdeltx 
    13981346      !!---------------------------------------------------------------------- 
    1399  
     1347      ! 
    14001348      zdeltx = xr - xl 
    1401       IF(abs(zdeltx) <= 10._wp * EPSILON(x)) THEN 
    1402         f = 0.5_wp * (fl + fr) 
     1349      IF( abs(zdeltx) <= 10._wp * EPSILON(x) ) THEN 
     1350         f = 0.5_wp * (fl + fr) 
    14031351      ELSE 
    1404         f = ( (x - xl ) * fr - ( x - xr ) * fl ) / zdeltx 
     1352         f = ( (x - xl ) * fr - ( x - xr ) * fl ) / zdeltx 
    14051353      ENDIF 
    1406  
     1354      ! 
    14071355   END FUNCTION interp1 
    14081356 
    14091357 
    1410    FUNCTION interp2(x, a, b, c, d)  RESULT(f) 
     1358   FUNCTION interp2( x, a, b, c, d )  RESULT(f) 
    14111359      !!---------------------------------------------------------------------- 
    14121360      !!                 ***  ROUTINE interp1  *** 
     
    14171365      !! 
    14181366      !!---------------------------------------------------------------------- 
    1419       IMPLICIT NONE 
    14201367      REAL(wp), INTENT(in) ::  x, a, b, c, d 
    14211368      REAL(wp)             ::  f ! value from the interpolation 
    14221369      !!---------------------------------------------------------------------- 
    1423  
     1370      ! 
    14241371      f = a + x* ( b + x * ( c + d * x ) ) 
    1425  
     1372      ! 
    14261373   END FUNCTION interp2 
    14271374 
    14281375 
    1429    FUNCTION interp3(x, a, b, c, d)  RESULT(f) 
     1376   FUNCTION interp3( x, a, b, c, d )  RESULT(f) 
    14301377      !!---------------------------------------------------------------------- 
    14311378      !!                 ***  ROUTINE interp1  *** 
     
    14371384      !! 
    14381385      !!---------------------------------------------------------------------- 
    1439       IMPLICIT NONE 
    14401386      REAL(wp), INTENT(in) ::  x, a, b, c, d 
    14411387      REAL(wp)             ::  f ! value from the interpolation 
    14421388      !!---------------------------------------------------------------------- 
    1443  
     1389      ! 
    14441390      f = b + x * ( 2._wp * c + 3._wp * d * x) 
    1445  
     1391      ! 
    14461392   END FUNCTION interp3 
    14471393 
    14481394 
    1449    FUNCTION integ_spline(xl, xr, a, b, c, d)  RESULT(f) 
     1395   FUNCTION integ_spline( xl, xr, a, b, c, d )  RESULT(f) 
    14501396      !!---------------------------------------------------------------------- 
    14511397      !!                 ***  ROUTINE interp1  *** 
     
    14561402      !! 
    14571403      !!---------------------------------------------------------------------- 
    1458       IMPLICIT NONE 
    14591404      REAL(wp), INTENT(in) ::  xl, xr, a, b, c, d 
    14601405      REAL(wp)             ::  za1, za2, za3 
    14611406      REAL(wp)             ::  f                   ! integration result 
    14621407      !!---------------------------------------------------------------------- 
    1463  
     1408      ! 
    14641409      za1 = 0.5_wp * b 
    14651410      za2 = c / 3.0_wp 
    14661411      za3 = 0.25_wp * d 
    1467  
     1412      ! 
    14681413      f  = xr * ( a + xr * ( za1 + xr * ( za2 + za3 * xr ) ) ) - & 
    14691414         & xl * ( a + xl * ( za1 + xl * ( za2 + za3 * xl ) ) ) 
    1470  
     1415      ! 
    14711416   END FUNCTION integ_spline 
    14721417 
Note: See TracChangeset for help on using the changeset viewer.