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 3116 for branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Ignore:
Timestamp:
2011-11-15T21:55:40+01:00 (13 years ago)
Author:
cetlod
Message:

dev_NEMO_MERGE_2011: add in changes dev_NOC_UKMO_MERGE developments

Location:
branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/DYN
Files:
11 edited
1 copied

Legend:

Unmodified
Added
Removed

  • branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90

    r2715 r3116  
    2727   USE sbc_oce, ONLY : ln_rnf   ! surface boundary condition: ocean 
    2828   USE sbcrnf          ! river runoff  
    29    USE obc_oce         ! ocean lateral open boundary condition 
    3029   USE cla             ! cross land advection             (cla_div routine) 
    3130   USE in_out_manager  ! I/O manager 
    … …  
    121120         END DO 
    122121 
    123 #if defined key_obc 
    124          IF( Agrif_Root() ) THEN 
    125             ! open boundaries (div must be zero behind the open boundary) 
    126             !  mpp remark: The zeroing of hdivn can probably be extended to 1->jpi/jpj for the correct row/column 
    127             IF( lp_obc_east  )   hdivn(nie0p1:nie1p1,nje0  :nje1  ,jk) = 0.e0      ! east 
    128             IF( lp_obc_west  )   hdivn(niw0  :niw1  ,njw0  :njw1  ,jk) = 0.e0      ! west 
    129             IF( lp_obc_north )   hdivn(nin0  :nin1  ,njn0p1:njn1p1,jk) = 0.e0      ! north 
    130             IF( lp_obc_south )   hdivn(nis0  :nis1  ,njs0  :njs1  ,jk) = 0.e0      ! south 
    131          ENDIF 
    132 #endif          
    133122         IF( .NOT. AGRIF_Root() ) THEN 
    134123            IF ((nbondi ==  1).OR.(nbondi == 2)) hdivn(nlci-1 , :     ,jk) = 0.e0      ! east 
    … …  
    304293         END DO   
    305294 
    306 #if defined key_obc 
    307          IF( Agrif_Root() ) THEN 
    308             ! open boundaries (div must be zero behind the open boundary) 
    309             !  mpp remark: The zeroing of hdivn can probably be extended to 1->jpi/jpj for the correct row/column 
    310             IF( lp_obc_east  )   hdivn(nie0p1:nie1p1,nje0  :nje1  ,jk) = 0.e0      ! east 
    311             IF( lp_obc_west  )   hdivn(niw0  :niw1  ,njw0  :njw1  ,jk) = 0.e0      ! west 
    312             IF( lp_obc_north )   hdivn(nin0  :nin1  ,njn0p1:njn1p1,jk) = 0.e0      ! north 
    313             IF( lp_obc_south )   hdivn(nis0  :nis1  ,njs0  :njs1  ,jk) = 0.e0      ! south 
    314          ENDIF 
    315 #endif          
    316295         IF( .NOT. AGRIF_Root() ) THEN 
    317296            IF ((nbondi ==  1).OR.(nbondi == 2)) hdivn(nlci-1 , :     ,jk) = 0.e0      ! east 

  • branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/DYN/dynbfr.F90

    r2715 r3116  
    1313   USE dom_oce         ! ocean space and time domain variables  
    1414   USE zdf_oce         ! ocean vertical physics variables 
     15   USE zdfbfr          ! ocean bottom friction variables 
    1516   USE trdmod          ! ocean active dynamics and tracers trends  
    1617   USE trdmod_oce      ! ocean variables trends 
    … …  
    5152      !!--------------------------------------------------------------------- 
    5253      ! 
    53       zm1_2dt = - 1._wp / ( 2._wp * rdt ) 
     54      IF( .not. ln_bfrimp) THEN     ! only for explicit bottom friction form 
     55                                    ! implicit bfr is implemented in dynzdf_imp 
     56                                    ! H. Liu,  Sept. 2011 
    5457 
    55       IF( l_trddyn )   THEN                      ! temporary save of ua and va trends 
    56          ztrduv(:,:,:,1) = ua(:,:,:) 
    57          ztrduv(:,:,:,2) = va(:,:,:) 
    58       ENDIF 
     58        zm1_2dt = - 1._wp / ( 2._wp * rdt ) 
     59 
     60        IF( l_trddyn )   THEN                      ! temporary save of ua and va trends 
     61           ztrduv(:,:,:,1) = ua(:,:,:) 
     62           ztrduv(:,:,:,2) = va(:,:,:) 
     63        ENDIF 
     64 
    5965 
    6066# if defined key_vectopt_loop 
    61       DO jj = 1, 1 
    62          DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling) 
     67        DO jj = 1, 1 
     68           DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling) 
    6369# else 
    64       DO jj = 2, jpjm1 
    65          DO ji = 2, jpim1 
     70        DO jj = 2, jpjm1 
     71           DO ji = 2, jpim1 
    6672# endif 
    67             ikbu = mbku(ji,jj)          ! deepest ocean u- & v-levels 
    68             ikbv = mbkv(ji,jj) 
    69             ! 
    70             ! Apply stability criteria on absolute value  : abs(bfr/e3) < 1/(2dt) => bfr/e3 > -1/(2dt) 
    71             ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + MAX(  bfrua(ji,jj) / fse3u(ji,jj,ikbu) , zm1_2dt  ) * ub(ji,jj,ikbu) 
    72             va(ji,jj,ikbv) = va(ji,jj,ikbv) + MAX(  bfrva(ji,jj) / fse3v(ji,jj,ikbv) , zm1_2dt  ) * vb(ji,jj,ikbv) 
    73          END DO 
    74       END DO 
     73              ikbu = mbku(ji,jj)          ! deepest ocean u- & v-levels 
     74              ikbv = mbkv(ji,jj) 
     75              ! 
     76              ! Apply stability criteria on absolute value  : abs(bfr/e3) < 1/(2dt) => bfr/e3 > -1/(2dt) 
     77              ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + MAX(  bfrua(ji,jj) / fse3u(ji,jj,ikbu) , zm1_2dt  ) * ub(ji,jj,ikbu) 
     78              va(ji,jj,ikbv) = va(ji,jj,ikbv) + MAX(  bfrva(ji,jj) / fse3v(ji,jj,ikbv) , zm1_2dt  ) * vb(ji,jj,ikbv) 
     79           END DO 
     80        END DO 
    7581 
    76       ! 
    77       IF( l_trddyn )   THEN                      ! save the vertical diffusive trends for further diagnostics 
    78          ztrduv(:,:,:,1) = ua(:,:,:) - ztrduv(:,:,:,1) 
    79          ztrduv(:,:,:,2) = va(:,:,:) - ztrduv(:,:,:,2) 
    80          CALL trd_mod( ztrduv(:,:,:,1), ztrduv(:,:,:,2), jpdyn_trd_bfr, 'DYN', kt ) 
    81       ENDIF 
    82       !                                          ! print mean trends (used for debugging) 
    83       IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' bfr  - Ua: ', mask1=umask,               & 
    84          &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    85       ! 
     82        ! 
     83        IF( l_trddyn )   THEN                      ! save the vertical diffusive trends for further diagnostics 
     84           ztrduv(:,:,:,1) = ua(:,:,:) - ztrduv(:,:,:,1) 
     85           ztrduv(:,:,:,2) = va(:,:,:) - ztrduv(:,:,:,2) 
     86           CALL trd_mod( ztrduv(:,:,:,1), ztrduv(:,:,:,2), jpdyn_trd_bfr, 'DYN', kt ) 
     87        ENDIF 
     88        !                                          ! print mean trends (used for debugging) 
     89        IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' bfr  - Ua: ', mask1=umask,               & 
     90           &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
     91        ! 
     92      ENDIF     ! end explicit bottom friction 
    8693   END SUBROUTINE dyn_bfr 
    8794 

  • branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90

    r2977 r3116  
    1414   !!             -   !  2005-11  (G. Madec) style & small optimisation 
    1515   !!            3.3  !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase 
     16   !!            3.4  !  2011-11  (A. Coward, H. Liu) introduction of prj scheme;  
     17   !!                 !           suppression of hel, wdj and rot options 
    1618   !!---------------------------------------------------------------------- 
    1719 
    … …  
    2325   !!       hpg_zps  : z-coordinate plus partial steps (interpolation) 
    2426   !!       hpg_sco  : s-coordinate (standard jacobian formulation) 
    25    !!       hpg_hel  : s-coordinate (helsinki modification) 
    26    !!       hpg_wdj  : s-coordinate (weighted density jacobian) 
    2727   !!       hpg_djc  : s-coordinate (Density Jacobian with Cubic polynomial) 
    28    !!       hpg_rot  : s-coordinate (ROTated axes scheme) 
     28   !!       hpg_prj  : s-coordinate (Pressure Jacobian with Cubic polynomial) 
    2929   !!---------------------------------------------------------------------- 
    3030   USE oce             ! ocean dynamics and tracers 
    … …  
    4848   LOGICAL , PUBLIC ::   ln_hpg_zps    = .FALSE.   !: z-coordinate - partial steps (interpolation) 
    4949   LOGICAL , PUBLIC ::   ln_hpg_sco    = .FALSE.   !: s-coordinate (standard jacobian formulation) 
    50    LOGICAL , PUBLIC ::   ln_hpg_hel    = .FALSE.   !: s-coordinate (helsinki modification) 
    51    LOGICAL , PUBLIC ::   ln_hpg_wdj    = .FALSE.   !: s-coordinate (weighted density jacobian) 
    5250   LOGICAL , PUBLIC ::   ln_hpg_djc    = .FALSE.   !: s-coordinate (Density Jacobian with Cubic polynomial) 
    53    LOGICAL , PUBLIC ::   ln_hpg_rot    = .FALSE.   !: s-coordinate (ROTated axes scheme) 
    54    REAL(wp), PUBLIC ::   rn_gamma      = 0._wp     !: weighting coefficient 
     51   LOGICAL , PUBLIC ::   ln_hpg_prj    = .FALSE.   !: s-coordinate (Pressure Jacobian scheme) 
    5552   LOGICAL , PUBLIC ::   ln_dynhpg_imp = .FALSE.   !: semi-implicite hpg flag 
    5653 
    57    INTEGER  ::   nhpg  =  0   ! = 0 to 6, type of pressure gradient scheme used ! (deduced from ln_hpg_... flags) 
     54   INTEGER  ::   nhpg  =  0   ! = 0 to 7, type of pressure gradient scheme used ! (deduced from ln_hpg_... flags) 
    5855 
    5956   !! * Substitutions 
    … …  
    9188      ENDIF       
    9289      ! 
    93       SELECT CASE ( nhpg )      ! Hydrastatic pressure gradient computation 
     90      SELECT CASE ( nhpg )      ! Hydrostatic pressure gradient computation 
    9491      CASE (  0 )   ;   CALL hpg_zco    ( kt )      ! z-coordinate 
    9592      CASE (  1 )   ;   CALL hpg_zps    ( kt )      ! z-coordinate plus partial steps (interpolation) 
    9693      CASE (  2 )   ;   CALL hpg_sco    ( kt )      ! s-coordinate (standard jacobian formulation) 
    97       CASE (  3 )   ;   CALL hpg_hel    ( kt )      ! s-coordinate (helsinki modification) 
    98       CASE (  4 )   ;   CALL hpg_wdj    ( kt )      ! s-coordinate (weighted density jacobian) 
    99       CASE (  5 )   ;   CALL hpg_djc    ( kt )      ! s-coordinate (Density Jacobian with Cubic polynomial) 
    100       CASE (  6 )   ;   CALL hpg_rot    ( kt )      ! s-coordinate (ROTated axes scheme) 
     94      CASE (  3 )   ;   CALL hpg_djc    ( kt )      ! s-coordinate (Density Jacobian with Cubic polynomial) 
     95      CASE (  4 )   ;   CALL hpg_prj    ( kt )      ! s-coordinate (Pressure Jacobian scheme) 
    10196      END SELECT 
    10297      ! 
    … …  
    125120      INTEGER ::   ioptio = 0      ! temporary integer 
    126121      !! 
    127       NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco, ln_hpg_hel,    & 
    128          &                 ln_hpg_wdj, ln_hpg_djc, ln_hpg_rot, rn_gamma  , ln_dynhpg_imp 
     122      NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco,     & 
     123         &                 ln_hpg_djc, ln_hpg_prj, ln_dynhpg_imp 
    129124      !!---------------------------------------------------------------------- 
    130125      ! 
    … …  
    140135         WRITE(numout,*) '      z-coord. - partial steps (interpolation)          ln_hpg_zps    = ', ln_hpg_zps 
    141136         WRITE(numout,*) '      s-coord. (standard jacobian formulation)          ln_hpg_sco    = ', ln_hpg_sco 
    142          WRITE(numout,*) '      s-coord. (helsinki modification)                  ln_hpg_hel    = ', ln_hpg_hel 
    143          WRITE(numout,*) '      s-coord. (weighted density jacobian)              ln_hpg_wdj    = ', ln_hpg_wdj 
    144137         WRITE(numout,*) '      s-coord. (Density Jacobian: Cubic polynomial)     ln_hpg_djc    = ', ln_hpg_djc 
    145          WRITE(numout,*) '      s-coord. (ROTated axes scheme)                    ln_hpg_rot    = ', ln_hpg_rot 
    146          WRITE(numout,*) '      weighting coeff. (wdj scheme)                     rn_gamma      = ', rn_gamma 
     138         WRITE(numout,*) '      s-coord. (Pressure Jacobian: Cubic polynomial)    ln_hpg_prj    = ', ln_hpg_prj 
    147139         WRITE(numout,*) '      time stepping: centered (F) or semi-implicit (T)  ln_dynhpg_imp = ', ln_dynhpg_imp 
    148140      ENDIF 
    149141      ! 
    150       IF( lk_vvl .AND. .NOT. ln_hpg_sco )   & 
    151          &   CALL ctl_stop('dyn_hpg_init : variable volume key_vvl require the standard jacobian formulation hpg_sco') 
     142      IF( lk_vvl .AND. .NOT. (ln_hpg_sco.OR.ln_hpg_prj) )   & 
     143         &   CALL ctl_stop('dyn_hpg_init : variable volume key_vvl requires:& 
     144                           & the standard jacobian formulation hpg_sco or & 
     145                           & the pressure jacobian formulation hpg_prj') 
    152146      ! 
    153147      !                               ! Set nhpg from ln_hpg_... flags 
    … …  
    155149      IF( ln_hpg_zps )   nhpg = 1 
    156150      IF( ln_hpg_sco )   nhpg = 2 
    157       IF( ln_hpg_hel )   nhpg = 3 
    158       IF( ln_hpg_wdj )   nhpg = 4 
    159       IF( ln_hpg_djc )   nhpg = 5 
    160       IF( ln_hpg_rot )   nhpg = 6 
    161       ! 
    162       !                               ! Consitency check 
     151      IF( ln_hpg_djc )   nhpg = 3 
     152      IF( ln_hpg_prj )   nhpg = 4 
     153      ! 
     154      !                               ! Consistency check 
    163155      ioptio = 0  
    164156      IF( ln_hpg_zco )   ioptio = ioptio + 1 
    165157      IF( ln_hpg_zps )   ioptio = ioptio + 1 
    166158      IF( ln_hpg_sco )   ioptio = ioptio + 1 
    167       IF( ln_hpg_hel )   ioptio = ioptio + 1 
    168       IF( ln_hpg_wdj )   ioptio = ioptio + 1 
    169159      IF( ln_hpg_djc )   ioptio = ioptio + 1 
    170       IF( ln_hpg_rot )   ioptio = ioptio + 1 
     160      IF( ln_hpg_prj )   ioptio = ioptio + 1 
    171161      IF( ioptio /= 1 )   CALL ctl_stop( 'NO or several hydrostatic pressure gradient options used' ) 
    172162      ! 
    … …  
    433423   END SUBROUTINE hpg_sco 
    434424 
    435  
    436    SUBROUTINE hpg_hel( kt ) 
    437       !!--------------------------------------------------------------------- 
    438       !!                  ***  ROUTINE hpg_hel  *** 
    439       !! 
    440       !! ** Method  :   s-coordinate case. 
    441       !!      The now hydrostatic pressure gradient at a given level 
    442       !!      jk is computed by taking the vertical integral of the in-situ  
    443       !!      density gradient along the model level from the suface to that  
    444       !!      level. s-coordinates (ln_sco): a corrective term is added 
    445       !!      to the horizontal pressure gradient : 
    446       !!         zhpi = grav .....  + 1/e1u mi(rhd) di[ grav dep3w ] 
    447       !!         zhpj = grav .....  + 1/e2v mj(rhd) dj[ grav dep3w ] 
    448       !!      add it to the general momentum trend (ua,va). 
    449       !!         ua = ua - 1/e1u * zhpi 
    450       !!         va = va - 1/e2v * zhpj 
    451       !! 
    452       !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 
    453       !!             - Save the trend (l_trddyn=T) 
    454       !!---------------------------------------------------------------------- 
    455       USE oce, ONLY:   tsa                         ! (tsa) used as 2 3D workspace 
    456       !! 
    457       INTEGER, INTENT(in) ::   kt    ! ocean time-step index 
    458       !! 
    459       INTEGER  ::   ji, jj, jk           ! dummy loop indices 
    460       REAL(wp) ::   zcoef0, zuap, zvap   ! temporary scalars 
    461       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zhpi, zhpj  
    462       !!---------------------------------------------------------------------- 
    463  
    464       zhpi => tsa(:,:,:,1)  
    465       zhpj => tsa(:,:,:,2)  
    466       ! 
    467       IF( kt == nit000 ) THEN 
    468          IF(lwp) WRITE(numout,*) 
    469          IF(lwp) WRITE(numout,*) 'dyn:hpg_hel : hydrostatic pressure gradient trend' 
    470          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate case, helsinki modified scheme' 
    471       ENDIF 
    472  
    473       ! Local constant initialization 
    474       zcoef0 = - grav * 0.5_wp 
    475   
    476       ! Surface value 
    477       DO jj = 2, jpjm1 
    478          DO ji = fs_2, fs_jpim1   ! vector opt. 
    479             ! hydrostatic pressure gradient along s-surfaces 
    480             zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( fse3t(ji+1,jj  ,1) * rhd(ji+1,jj  ,1)  & 
    481                &                                  - fse3t(ji  ,jj  ,1) * rhd(ji  ,jj  ,1) ) 
    482             zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( fse3t(ji  ,jj+1,1) * rhd(ji  ,jj+1,1)  & 
    483                &                                  - fse3t(ji  ,jj  ,1) * rhd(ji  ,jj  ,1) ) 
    484             ! s-coordinate pressure gradient correction 
    485             zuap = -zcoef0 * ( rhd   (ji+1,jj,1) + rhd   (ji,jj,1) )   & 
    486                &           * ( fsdept(ji+1,jj,1) - fsdept(ji,jj,1) ) / e1u(ji,jj) 
    487             zvap = -zcoef0 * ( rhd   (ji,jj+1,1) + rhd   (ji,jj,1) )   & 
    488                &           * ( fsdept(ji,jj+1,1) - fsdept(ji,jj,1) ) / e2v(ji,jj) 
    489             ! add to the general momentum trend 
    490             ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap 
    491             va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap 
    492          END DO 
    493       END DO 
    494       ! 
    495       ! interior value (2=<jk=<jpkm1) 
    496       DO jk = 2, jpkm1 
    497          DO jj = 2, jpjm1 
    498             DO ji = fs_2, fs_jpim1   ! vector opt. 
    499                ! hydrostatic pressure gradient along s-surfaces 
    500                zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) & 
    501                   &           +  zcoef0 / e1u(ji,jj) * ( fse3t(ji+1,jj,jk  ) * rhd(ji+1,jj,jk)     & 
    502                   &                                     -fse3t(ji  ,jj,jk  ) * rhd(ji  ,jj,jk)   ) & 
    503                   &           +  zcoef0 / e1u(ji,jj) * ( fse3t(ji+1,jj,jk-1) * rhd(ji+1,jj,jk-1)   & 
    504                   &                                     -fse3t(ji  ,jj,jk-1) * rhd(ji  ,jj,jk-1) ) 
    505                zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) & 
    506                   &           +  zcoef0 / e2v(ji,jj) * ( fse3t(ji,jj+1,jk  ) * rhd(ji,jj+1,jk)     & 
    507                   &                                     -fse3t(ji,jj  ,jk  ) * rhd(ji,jj,  jk)   ) & 
    508                   &           +  zcoef0 / e2v(ji,jj) * ( fse3t(ji,jj+1,jk-1) * rhd(ji,jj+1,jk-1)   & 
    509                   &                                     -fse3t(ji,jj  ,jk-1) * rhd(ji,jj,  jk-1) ) 
    510                ! s-coordinate pressure gradient correction 
    511                zuap = - zcoef0 * ( rhd   (ji+1,jj,jk) + rhd   (ji,jj,jk) )   & 
    512                   &            * ( fsdept(ji+1,jj,jk) - fsdept(ji,jj,jk) ) / e1u(ji,jj) 
    513                zvap = - zcoef0 * ( rhd   (ji,jj+1,jk) + rhd   (ji,jj,jk) )   & 
    514                   &            * ( fsdept(ji,jj+1,jk) - fsdept(ji,jj,jk) ) / e2v(ji,jj) 
    515                ! add to the general momentum trend 
    516                ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap 
    517                va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) + zvap 
    518             END DO 
    519          END DO 
    520       END DO 
    521       ! 
    522    END SUBROUTINE hpg_hel 
    523  
    524  
    525    SUBROUTINE hpg_wdj( kt ) 
    526       !!--------------------------------------------------------------------- 
    527       !!                  ***  ROUTINE hpg_wdj  *** 
    528       !! 
    529       !! ** Method  :   Weighted Density Jacobian (wdj) scheme (song 1998) 
    530       !!      The weighting coefficients from the namelist parameter rn_gamma 
    531       !!      (alpha=0.5-rn_gamma ; beta=1-alpha=0.5+rn_gamma 
    532       !! 
    533       !! Reference : Song, Mon. Wea. Rev., 126, 3213-3230, 1998. 
    534       !!---------------------------------------------------------------------- 
    535       USE oce, ONLY:   tsa                         ! (tsa) used as 2 3D workspace 
    536       !! 
    537       INTEGER, INTENT(in) ::   kt    ! ocean time-step index 
    538       !! 
    539       INTEGER  ::   ji, jj, jk           ! dummy loop indices 
    540       REAL(wp) ::   zcoef0, zuap, zvap   ! temporary scalars 
    541       REAL(wp) ::   zalph , zbeta        !    "         " 
    542       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zhpi, zhpj  
    543       !!---------------------------------------------------------------------- 
    544       ! 
    545       zhpi => tsa(:,:,:,1)  
    546       zhpj => tsa(:,:,:,2)  
    547       ! 
    548       IF( kt == nit000 ) THEN 
    549          IF(lwp) WRITE(numout,*) 
    550          IF(lwp) WRITE(numout,*) 'dyn:hpg_wdj : hydrostatic pressure gradient trend' 
    551          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   Weighted Density Jacobian' 
    552       ENDIF 
    553  
    554       ! Local constant initialization 
    555       zcoef0 = - grav * 0.5_wp 
    556       zalph  = 0.5_wp - rn_gamma    ! weighting coefficients (alpha=0.5-rn_gamma 
    557       zbeta  = 0.5_wp + rn_gamma    !                        (beta =1-alpha=0.5+rn_gamma 
    558  
    559       ! Surface value (no ponderation) 
    560       DO jj = 2, jpjm1 
    561          DO ji = fs_2, fs_jpim1   ! vector opt. 
    562             ! hydrostatic pressure gradient along s-surfaces 
    563             zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * (  fse3w(ji+1,jj  ,1) * rhd(ji+1,jj  ,1)   & 
    564                &                                   - fse3w(ji  ,jj  ,1) * rhd(ji  ,jj  ,1)  ) 
    565             zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * (  fse3w(ji  ,jj+1,1) * rhd(ji  ,jj+1,1)   & 
    566                &                                   - fse3w(ji  ,jj  ,1) * rhd(ji,  jj  ,1)  ) 
    567             ! s-coordinate pressure gradient correction 
    568             zuap = -zcoef0 * ( rhd   (ji+1,jj,1) + rhd   (ji,jj,1) )   & 
    569                &           * ( fsde3w(ji+1,jj,1) - fsde3w(ji,jj,1) ) / e1u(ji,jj) 
    570             zvap = -zcoef0 * ( rhd   (ji,jj+1,1) + rhd   (ji,jj,1) )   & 
    571                &           * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj) 
    572             ! add to the general momentum trend 
    573             ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap 
    574             va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap 
    575          END DO 
    576       END DO 
    577  
    578       ! Interior value (2=<jk=<jpkm1) (weighted with zalph & zbeta) 
    579       DO jk = 2, jpkm1 
    580          DO jj = 2, jpjm1 
    581             DO ji = fs_2, fs_jpim1   ! vector opt. 
    582                zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj)                            & 
    583                   &           * (   (            fsde3w(ji+1,jj,jk  ) + fsde3w(ji,jj,jk  )        & 
    584                   &                            - fsde3w(ji+1,jj,jk-1) - fsde3w(ji,jj,jk-1)    )   & 
    585                   &               * (  zalph * ( rhd   (ji+1,jj,jk-1) - rhd   (ji,jj,jk-1) )      & 
    586                   &                  + zbeta * ( rhd   (ji+1,jj,jk  ) - rhd   (ji,jj,jk  ) )  )   & 
    587                   &             -   (            rhd   (ji+1,jj,jk  ) + rhd   (ji,jj,jk  )        & 
    588                   &                           - rhd   (ji+1,jj,jk-1) - rhd   (ji,jj,jk-1)     )   & 
    589                   &               * (  zalph * ( fsde3w(ji+1,jj,jk-1) - fsde3w(ji,jj,jk-1) )      & 
    590                   &                  + zbeta * ( fsde3w(ji+1,jj,jk  ) - fsde3w(ji,jj,jk  ) )  )  ) 
    591                zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj)                            & 
    592                   &           * (   (           fsde3w(ji,jj+1,jk  ) + fsde3w(ji,jj,jk  )         & 
    593                   &                           - fsde3w(ji,jj+1,jk-1) - fsde3w(ji,jj,jk-1)     )   & 
    594                   &               * (  zalph * ( rhd   (ji,jj+1,jk-1) - rhd   (ji,jj,jk-1) )      & 
    595                   &                  + zbeta * ( rhd   (ji,jj+1,jk  ) - rhd   (ji,jj,jk  ) )  )   & 
    596                   &             -   (            rhd   (ji,jj+1,jk  ) + rhd   (ji,jj,jk  )        & 
    597                   &                            - rhd   (ji,jj+1,jk-1) - rhd   (ji,jj,jk-1)    )   & 
    598                   &               * (  zalph * ( fsde3w(ji,jj+1,jk-1) - fsde3w(ji,jj,jk-1) )      & 
    599                   &                  + zbeta * ( fsde3w(ji,jj+1,jk  ) - fsde3w(ji,jj,jk  ) )  )  ) 
    600                ! add to the general momentum trend 
    601                ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) 
    602                va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) 
    603             END DO 
    604          END DO 
    605       END DO 
    606       ! 
    607    END SUBROUTINE hpg_wdj 
    608  
    609  
    610425   SUBROUTINE hpg_djc( kt ) 
    611426      !!--------------------------------------------------------------------- 
    … …  
    843658 
    844659 
    845    SUBROUTINE hpg_rot( kt ) 
     660   SUBROUTINE hpg_prj( kt ) 
    846661      !!--------------------------------------------------------------------- 
    847       !!                  ***  ROUTINE hpg_rot  *** 
    848       !! 
    849       !! ** Method  :   rotated axes scheme (Thiem and Berntsen 2005) 
    850       !! 
    851       !! Reference: Thiem & Berntsen, Ocean Modelling, In press, 2005. 
    852       !!---------------------------------------------------------------------- 
     662      !!                  ***  ROUTINE hpg_prj  *** 
     663      !! 
     664      !! ** Method  :   s-coordinate case. 
     665      !!      A Pressure-Jacobian horizontal pressure gradient method 
     666      !!      based on the constrained cubic-spline interpolation for 
     667      !!      all vertical coordinate systems 
     668      !! 
     669      !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 
     670      !!             - Save the trend (l_trddyn=T) 
     671      !! 
     672      !!---------------------------------------------------------------------- 
     673 
    853674      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released 
    854675      USE oce     , ONLY:   tsa                          ! (tsa) used as 2 3D workspace 
    855       USE wrk_nemo, ONLY:   zdistr  => wrk_2d_1 , zsina   => wrk_2d_2 , zcosa  => wrk_2d_3 
    856       USE wrk_nemo, ONLY:   zhpiorg => wrk_3d_1 , zhpirot => wrk_3d_2 
    857       USE wrk_nemo, ONLY:   zhpitra => wrk_3d_3 , zhpine  => wrk_3d_4 
    858       USE wrk_nemo, ONLY:   zhpjorg => wrk_3d_5 , zhpjrot => wrk_3d_6 
    859       USE wrk_nemo, ONLY:   zhpjtra => wrk_3d_7 , zhpjne  => wrk_3d_8 
    860       !! 
    861       INTEGER, INTENT(in) ::   kt    ! ocean time-step index 
    862       !! 
    863       INTEGER  ::   ji, jj, jk          ! dummy loop indices 
    864       REAL(wp) ::   zforg, zcoef0, zuap, zmskd1, zmskd1m   ! temporary scalar 
    865       REAL(wp) ::   zfrot        , zvap, zmskd2, zmskd2m   !    "         " 
    866       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zhpi, zhpj  
    867       !!---------------------------------------------------------------------- 
    868  
    869       IF( wrk_in_use(2, 1,2,3)             .OR.   & 
    870           wrk_in_use(3, 1,2,3,4,5,6,7,8) ) THEN 
    871          CALL ctl_stop('dyn:hpg_rot: requested workspace arrays unavailable')   ;   RETURN 
     676      USE wrk_nemo, ONLY:   zhpi => wrk_3d_3  
     677      USE wrk_nemo, ONLY:   zu => wrk_3d_4  
     678      USE wrk_nemo, ONLY:   zv => wrk_3d_5 
     679      USE wrk_nemo, ONLY:   sp => wrk_3d_6  
     680      USE wrk_nemo, ONLY:   sp => wrk_3d_7 
     681      USE wrk_nemo, ONLY:   sp => wrk_3d_8 
     682      USE wrk_nemo, ONLY:   sp => wrk_3d_9 
     683      USE wrk_nemo, ONLY:   sp => wrk_3d_10 
     684      USE wrk_nemo, ONLY:   sp => wrk_3d_11 
     685      !! 
     686      !!---------------------------------------------------------------------- 
     687      !! 
     688      INTEGER, PARAMETER  :: polynomial_type = 1    ! 1: cubic spline, 2: linear 
     689      INTEGER, INTENT(in) ::   kt                   ! ocean time-step index 
     690      !! 
     691      INTEGER  ::   ji, jj, jk, jkk                 ! dummy loop indices 
     692      REAL(wp) ::   zcoef0, znad                    ! temporary scalars 
     693      !! 
     694      !! The local variables for the correction term 
     695      INTEGER  :: jk1, jis, jid, jjs, jjd 
     696      REAL(wp) :: zuijk, zvijk, zpwes, zpwed, zpnss, zpnsd, zdeps 
     697      REAL(wp) :: zrhdt1  
     698      REAL(wp) :: zdpdx1, zdpdx2, zdpdy1, zdpdy2 
     699      INTEGER  :: zbhitwe, zbhitns 
     700      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zdeptht, zrhh  
     701      !!---------------------------------------------------------------------- 
     702 
     703      IF( wrk_in_use(3, 3,4,5,6,7,8,9,10,11) ) THEN 
     704         CALL ctl_stop('dyn:hpg_prj: requested workspace arrays unavailable')   ;   RETURN 
    872705      ENDIF 
    873706      ! 
    874       zhpi => tsa(:,:,:,1)  
    875       zhpj => tsa(:,:,:,2)  
     707      zdeptht => tsa(:,:,:,1)  
     708      zrhh    => tsa(:,:,:,2)  
    876709 
    877710      IF( kt == nit000 ) THEN 
    878711         IF(lwp) WRITE(numout,*) 
    879          IF(lwp) WRITE(numout,*) 'dyn:hpg_rot : hydrostatic pressure gradient trend' 
    880          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate case, rotated axes scheme used' 
     712         IF(lwp) WRITE(numout,*) 'dyn:hpg_prj : hydrostatic pressure gradient trend' 
     713         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate case, cubic spline pressure Jacobian' 
    881714      ENDIF 
    882715 
    883       ! ------------------------------- 
    884       !  Local constant initialization 
    885       ! ------------------------------- 
    886       zcoef0 = - grav * 0.5_wp 
    887       zforg  = 0.95_wp 
    888       zfrot  = 1._wp - zforg 
    889  
    890       ! inverse of the distance between 2 diagonal T-points (defined at F-point) (here zcoef0/distance) 
    891       zdistr(:,:) = zcoef0 / SQRT( e1f(:,:)*e1f(:,:) + e2f(:,:)*e1f(:,:) ) 
    892  
    893       ! sinus and cosinus of diagonal angle at F-point 
    894       zsina(:,:) = ATAN2( e2f(:,:), e1f(:,:) ) 
    895       zcosa(:,:) = COS( zsina(:,:) ) 
    896       zsina(:,:) = SIN( zsina(:,:) ) 
    897  
    898       ! --------------- 
    899       !  Surface value 
    900       ! --------------- 
    901       ! compute and add to the general trend the pressure gradients along the axes 
    902       DO jj = 2, jpjm1 
    903          DO ji = fs_2, fs_jpim1   ! vector opt. 
    904             ! hydrostatic pressure gradient along s-surfaces 
    905             zhpiorg(ji,jj,1) = zcoef0 / e1u(ji,jj) * (  fse3t(ji+1,jj,1) * rhd(ji+1,jj,1)   & 
    906                &                                      - fse3t(ji  ,jj,1) * rhd(ji  ,jj,1)  ) 
    907             zhpjorg(ji,jj,1) = zcoef0 / e2v(ji,jj) * (  fse3t(ji,jj+1,1) * rhd(ji,jj+1,1)   & 
    908                &                                      - fse3t(ji,jj  ,1) * rhd(ji,jj  ,1)  ) 
    909             ! s-coordinate pressure gradient correction 
    910             zuap = -zcoef0 * ( rhd   (ji+1,jj  ,1) + rhd   (ji,jj,1) )   & 
    911                &           * ( fsdept(ji+1,jj  ,1) - fsdept(ji,jj,1) ) / e1u(ji,jj) 
    912             zvap = -zcoef0 * ( rhd   (ji  ,jj+1,1) + rhd   (ji,jj,1) )   & 
    913                &           * ( fsdept(ji  ,jj+1,1) - fsdept(ji,jj,1) ) / e2v(ji,jj) 
    914             ! add to the general momentum trend 
    915             ua(ji,jj,1) = ua(ji,jj,1) + zforg * ( zhpiorg(ji,jj,1) + zuap ) 
    916             va(ji,jj,1) = va(ji,jj,1) + zforg * ( zhpjorg(ji,jj,1) + zvap ) 
    917          END DO 
    918       END DO 
    919  
    920       ! compute the pressure gradients in the diagonal directions 
    921       DO jj = 1, jpjm1 
    922          DO ji = 1, fs_jpim1   ! vector opt. 
    923             zmskd1 = tmask(ji+1,jj+1,1) * tmask(ji  ,jj,1)      ! mask in the 1st diagnonal 
    924             zmskd2 = tmask(ji  ,jj+1,1) * tmask(ji+1,jj,1)      ! mask in the 2nd diagnonal 
    925             ! hydrostatic pressure gradient along s-surfaces 
    926             zhpitra(ji,jj,1) = zdistr(ji,jj) * zmskd1 * (  fse3t(ji+1,jj+1,1) * rhd(ji+1,jj+1,1)   & 
    927                &                                         - fse3t(ji  ,jj  ,1) * rhd(ji  ,jj  ,1)  ) 
    928             zhpjtra(ji,jj,1) = zdistr(ji,jj) * zmskd2 * (  fse3t(ji  ,jj+1,1) * rhd(ji  ,jj+1,1)   & 
    929                &                                         - fse3t(ji+1,jj  ,1) * rhd(ji+1,jj  ,1)  ) 
    930             ! s-coordinate pressure gradient correction 
    931             zuap = -zdistr(ji,jj) * zmskd1 * ( rhd   (ji+1,jj+1,1) + rhd   (ji  ,jj,1) )   & 
    932                &                           * ( fsdept(ji+1,jj+1,1) - fsdept(ji  ,jj,1) ) 
    933             zvap = -zdistr(ji,jj) * zmskd2 * ( rhd   (ji  ,jj+1,1) + rhd   (ji+1,jj,1) )   & 
    934                &                           * ( fsdept(ji  ,jj+1,1) - fsdept(ji+1,jj,1) ) 
    935             ! back rotation 
    936             zhpine(ji,jj,1) = zcosa(ji,jj) * ( zhpitra(ji,jj,1) + zuap )   & 
    937                &            - zsina(ji,jj) * ( zhpjtra(ji,jj,1) + zvap ) 
    938             zhpjne(ji,jj,1) = zsina(ji,jj) * ( zhpitra(ji,jj,1) + zuap )   & 
    939                &            + zcosa(ji,jj) * ( zhpjtra(ji,jj,1) + zvap ) 
    940          END DO 
    941       END DO 
    942  
    943       ! interpolate and add to the general trend the diagonal gradient 
    944       DO jj = 2, jpjm1 
    945          DO ji = fs_2, fs_jpim1   ! vector opt. 
    946             ! averaging 
    947             zhpirot(ji,jj,1) = 0.5 * ( zhpine(ji,jj,1) + zhpine(ji  ,jj-1,1) ) 
    948             zhpjrot(ji,jj,1) = 0.5 * ( zhpjne(ji,jj,1) + zhpjne(ji-1,jj  ,1) ) 
    949             ! add to the general momentum trend 
    950             ua(ji,jj,1) = ua(ji,jj,1) + zfrot * zhpirot(ji,jj,1)  
    951             va(ji,jj,1) = va(ji,jj,1) + zfrot * zhpjrot(ji,jj,1)  
    952          END DO 
    953       END DO 
    954  
    955       ! ----------------- 
    956       ! 2. interior value (2=<jk=<jpkm1) 
    957       ! ----------------- 
    958       ! compute and add to the general trend the pressure gradients along the axes 
    959       DO jk = 2, jpkm1 
    960          DO jj = 2, jpjm1 
    961             DO ji = fs_2, fs_jpim1   ! vector opt. 
    962                ! hydrostatic pressure gradient along s-surfaces 
    963                zhpiorg(ji,jj,jk) = zhpiorg(ji,jj,jk-1)                                                 & 
    964                   &              +  zcoef0 / e1u(ji,jj) * (  fse3t(ji+1,jj,jk  ) * rhd(ji+1,jj,jk  )   & 
    965                   &                                        - fse3t(ji  ,jj,jk  ) * rhd(ji  ,jj,jk  )   & 
    966                   &                                        + fse3t(ji+1,jj,jk-1) * rhd(ji+1,jj,jk-1)   & 
    967                   &                                        - fse3t(ji  ,jj,jk-1) * rhd(ji  ,jj,jk-1)  ) 
    968                zhpjorg(ji,jj,jk) = zhpjorg(ji,jj,jk-1)                                                 & 
    969                   &              +  zcoef0 / e2v(ji,jj) * (  fse3t(ji,jj+1,jk  ) * rhd(ji,jj+1,jk  )   & 
    970                   &                                        - fse3t(ji,jj  ,jk  ) * rhd(ji,jj,  jk  )   & 
    971                   &                                        + fse3t(ji,jj+1,jk-1) * rhd(ji,jj+1,jk-1)   & 
    972                   &                                        - fse3t(ji,jj  ,jk-1) * rhd(ji,jj,  jk-1)  ) 
    973                ! s-coordinate pressure gradient correction 
    974                zuap = - zcoef0 * ( rhd   (ji+1,jj  ,jk) + rhd   (ji,jj,jk) )   & 
    975                   &            * ( fsdept(ji+1,jj  ,jk) - fsdept(ji,jj,jk) ) / e1u(ji,jj) 
    976                zvap = - zcoef0 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji,jj,jk) )   & 
    977                   &            * ( fsdept(ji  ,jj+1,jk) - fsdept(ji,jj,jk) ) / e2v(ji,jj) 
    978                ! add to the general momentum trend 
    979                ua(ji,jj,jk) = ua(ji,jj,jk) + zforg*( zhpiorg(ji,jj,jk) + zuap ) 
    980                va(ji,jj,jk) = va(ji,jj,jk) + zforg*( zhpjorg(ji,jj,jk) + zvap ) 
    981             END DO 
    982          END DO 
    983       END DO 
    984  
    985       ! compute the pressure gradients in the diagonal directions 
    986       DO jk = 2, jpkm1 
    987          DO jj = 1, jpjm1 
    988             DO ji = 1, fs_jpim1   ! vector opt. 
    989                zmskd1  = tmask(ji+1,jj+1,jk  ) * tmask(ji  ,jj,jk  )      ! level jk   mask in the 1st diagnonal 
    990                zmskd1m = tmask(ji+1,jj+1,jk-1) * tmask(ji  ,jj,jk-1)      ! level jk-1    "               "      
    991                zmskd2  = tmask(ji  ,jj+1,jk  ) * tmask(ji+1,jj,jk  )      ! level jk   mask in the 2nd diagnonal 
    992                zmskd2m = tmask(ji  ,jj+1,jk-1) * tmask(ji+1,jj,jk-1)      ! level jk-1    "               "      
    993                ! hydrostatic pressure gradient along s-surfaces 
    994                zhpitra(ji,jj,jk) = zhpitra(ji,jj,jk-1)                                                       & 
    995                   &              + zdistr(ji,jj) * zmskd1  * ( fse3t(ji+1,jj+1,jk  ) * rhd(ji+1,jj+1,jk)     & 
    996                   &                                           -fse3t(ji  ,jj  ,jk  ) * rhd(ji  ,jj  ,jk) )   & 
    997                   &              + zdistr(ji,jj) * zmskd1m * ( fse3t(ji+1,jj+1,jk-1) * rhd(ji+1,jj+1,jk-1)   & 
    998                   &                                           -fse3t(ji  ,jj  ,jk-1) * rhd(ji  ,jj  ,jk-1) ) 
    999                zhpjtra(ji,jj,jk) = zhpjtra(ji,jj,jk-1)                                                       & 
    1000                   &              + zdistr(ji,jj) * zmskd2  * ( fse3t(ji  ,jj+1,jk  ) * rhd(ji  ,jj+1,jk)     & 
    1001                   &                                           -fse3t(ji+1,jj  ,jk  ) * rhd(ji+1,jj,  jk) )   & 
    1002                   &              + zdistr(ji,jj) * zmskd2m * ( fse3t(ji  ,jj+1,jk-1) * rhd(ji  ,jj+1,jk-1)   & 
    1003                   &                                           -fse3t(ji+1,jj  ,jk-1) * rhd(ji+1,jj,  jk-1) ) 
    1004                ! s-coordinate pressure gradient correction 
    1005                zuap = - zdistr(ji,jj) * zmskd1 * ( rhd   (ji+1,jj+1,jk) + rhd   (ji  ,jj,jk) )   & 
    1006                   &                            * ( fsdept(ji+1,jj+1,jk) - fsdept(ji  ,jj,jk) ) 
    1007                zvap = - zdistr(ji,jj) * zmskd2 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji+1,jj,jk) )   & 
    1008                   &                            * ( fsdept(ji  ,jj+1,jk) - fsdept(ji+1,jj,jk) ) 
    1009                ! back rotation 
    1010                zhpine(ji,jj,jk) = zcosa(ji,jj) * ( zhpitra(ji,jj,jk) + zuap )   & 
    1011                   &             - zsina(ji,jj) * ( zhpjtra(ji,jj,jk) + zvap ) 
    1012                zhpjne(ji,jj,jk) = zsina(ji,jj) * ( zhpitra(ji,jj,jk) + zuap )   & 
    1013                   &             + zcosa(ji,jj) * ( zhpjtra(ji,jj,jk) + zvap ) 
    1014             END DO 
    1015          END DO 
    1016       END DO 
    1017  
    1018       ! interpolate and add to the general trend 
    1019       DO jk = 2, jpkm1 
    1020          DO jj = 2, jpjm1 
    1021             DO ji = fs_2, fs_jpim1   ! vector opt. 
    1022                ! averaging 
    1023                zhpirot(ji,jj,jk) = 0.5 * ( zhpine(ji,jj,jk) + zhpine(ji  ,jj-1,jk) ) 
    1024                zhpjrot(ji,jj,jk) = 0.5 * ( zhpjne(ji,jj,jk) + zhpjne(ji-1,jj  ,jk) ) 
    1025                ! add to the general momentum trend 
    1026                ua(ji,jj,jk) = ua(ji,jj,jk) + zfrot * zhpirot(ji,jj,jk)  
    1027                va(ji,jj,jk) = va(ji,jj,jk) + zfrot * zhpjrot(ji,jj,jk)  
    1028             END DO 
    1029          END DO 
    1030       END DO 
    1031       ! 
    1032       IF( wrk_not_released(2, 1,2,3)           .OR.   & 
    1033           wrk_not_released(3, 1,2,3,4,5,6,7,8) )   CALL ctl_stop('dyn:hpg_rot: failed to release workspace arrays') 
    1034       ! 
    1035    END SUBROUTINE hpg_rot 
     716      !!---------------------------------------------------------------------- 
     717      ! Local constant initialization 
     718      zcoef0 = - grav  
     719      znad = 0.0_wp 
     720      IF( lk_vvl ) znad = 1._wp 
     721 
     722      ! Clean 3-D work arrays 
     723      zhpi(:,:,:) = 0._wp 
     724      zrhh(:,:,:) = rhd(:,:,:) 
     725       
     726      ! Preparing vertical density profile for hybrid-sco coordinate 
     727      DO jj = 1, jpj 
     728        DO ji = 1, jpi    
     729          jk = mbathy(ji,jj) 
     730          IF( jk <= 0 ) THEN; zrhh(ji,jj,:) = 0._wp 
     731          ELSE IF(jk == 1) THEN; zrhh(ji,jj, jk+1:jpk) = rhd(ji,jj,jk) 
     732          ELSE IF(jk < jpkm1) THEN 
     733             DO jkk = jk+1, jpk 
     734                zrhh(ji,jj,jkk) = interp1(fsde3w(ji,jj,jkk),   fsde3w(ji,jj,jkk-1), & 
     735                                         fsde3w(ji,jj,jkk-2), rhd(ji,jj,jkk-1), rhd(ji,jj,jkk-2)) 
     736             END DO  
     737          ENDIF 
     738        END DO 
     739      END DO 
     740 
     741      DO jj = 1, jpj 
     742        DO ji = 1, jpi 
     743          zdeptht(ji,jj,1) = 0.5_wp * fse3w(ji,jj,1) 
     744          zdeptht(ji,jj,1) = zdeptht(ji,jj,1) - sshn(ji,jj) * znad 
     745          DO jk = 2, jpk 
     746             zdeptht(ji,jj,jk) = zdeptht(ji,jj,jk-1) + fse3w(ji,jj,jk) 
     747          END DO 
     748        END DO 
     749      END DO 
     750 
     751      DO jk = 1, jpkm1 
     752        DO jj = 1, jpj 
     753          DO ji = 1, jpi 
     754            fsp(ji,jj,jk) = zrhh(ji,jj,jk) 
     755            xsp(ji,jj,jk) = zdeptht(ji,jj,jk) 
     756          END DO 
     757        END DO 
     758      END DO 
     759 
     760      ! Construct the vertical density profile with the  
     761      ! constrained cubic spline interpolation 
     762      CALL cspline(fsp,xsp,asp,bsp,csp,dsp,polynomial_type)       
     763 
     764      ! Calculate the hydrostatic pressure at T(ji,jj,1) 
     765      DO jj = 2, jpj 
     766        DO ji = 2, jpi  
     767          zrhdt1 = zrhh(ji,jj,1) - interp3(zdeptht(ji,jj,1),asp(ji,jj,1), & 
     768                                         bsp(ji,jj,1),   csp(ji,jj,1), & 
     769                                         dsp(ji,jj,1) ) * 0.5_wp * zdeptht(ji,jj,1) 
     770          zrhdt1 = MAX(zrhdt1, 1000._wp - rau0)        ! no lighter than fresh water 
     771 
     772          ! assuming linear profile across the top half surface layer 
     773          zhpi(ji,jj,1) =  0.5_wp * fse3w(ji,jj,1) * zrhdt1   
     774        END DO 
     775      END DO 
     776 
     777      ! Calculate the pressure at T(ji,jj,2:jpkm1) 
     778      DO jk = 2, jpkm1                                   
     779        DO jj = 2, jpj      
     780          DO ji = 2, jpi 
     781            zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) +                          & 
     782                             integ2(zdeptht(ji,jj,jk-1), zdeptht(ji,jj,jk),& 
     783                                    asp(ji,jj,jk-1),    bsp(ji,jj,jk-1), & 
     784                                    csp(ji,jj,jk-1),    dsp(ji,jj,jk-1)) 
     785          END DO 
     786        END DO 
     787      END DO 
     788 
     789      ! Z coordinate of U(ji,jj,1:jpkm1) and V(ji,jj,1:jpkm1) 
     790      DO jj = 2, jpjm1      
     791        DO ji = 2, jpim1   
     792          zu(ji,jj,1) = - ( fse3u(ji,jj,1) - sshu_n(ji,jj) * znad) 
     793          zv(ji,jj,1) = - ( fse3v(ji,jj,1) - sshv_n(ji,jj) * znad) 
     794        END DO 
     795      END DO 
     796 
     797      DO jk = 2, jpkm1                                   
     798        DO jj = 2, jpjm1      
     799          DO ji = 2, jpim1   
     800            zu(ji,jj,jk) = zu(ji,jj,jk-1)- fse3u(ji,jj,jk) 
     801            zv(ji,jj,jk) = zv(ji,jj,jk-1)- fse3v(ji,jj,jk) 
     802          END DO 
     803        END DO 
     804      END DO 
     805                
     806      DO jk = 1, jpkm1                                   
     807        DO jj = 2, jpjm1      
     808          DO ji = 2, jpim1   
     809            zu(ji,jj,jk) = zu(ji,jj,jk) + 0.5_wp * fse3u(ji,jj,jk) 
     810            zv(ji,jj,jk) = zv(ji,jj,jk) + 0.5_wp * fse3v(ji,jj,jk) 
     811          END DO 
     812        END DO 
     813      END DO 
     814 
     815      DO jk = 1, jpkm1                                   
     816        DO jj = 2, jpjm1      
     817          DO ji = 2, jpim1   
     818            zpwes = 0._wp; zpwed = 0._wp 
     819            zpnss = 0._wp; zpnsd = 0._wp 
     820            zuijk = zu(ji,jj,jk) 
     821            zvijk = zv(ji,jj,jk) 
     822 
     823            !!!!!     for u equation 
     824            IF( jk <= mbku(ji,jj) ) THEN 
     825               IF( -zdeptht(ji+1,jj,mbku(ji,jj)) >= -zdeptht(ji,jj,mbku(ji,jj)) ) THEN 
     826                 jis = ji + 1; jid = ji 
     827               ELSE 
     828                 jis = ji;     jid = ji +1 
     829               ENDIF 
     830 
     831               ! integrate the pressure on the shallow side 
     832               jk1 = jk  
     833               zbhitwe = 0 
     834               DO WHILE ( -zdeptht(jis,jj,jk1) > zuijk ) 
     835                 IF( jk1 == mbku(ji,jj) ) THEN 
     836                   zbhitwe = 1 
     837                   EXIT 
     838                 ENDIF 
     839                 zdeps = MIN(zdeptht(jis,jj,jk1+1), -zuijk) 
     840                 zpwes = zpwes +                                    &  
     841                      integ2(zdeptht(jis,jj,jk1), zdeps,            & 
     842                             asp(jis,jj,jk1),    bsp(jis,jj,jk1), & 
     843                             csp(jis,jj,jk1),    dsp(jis,jj,jk1)) 
     844                 jk1 = jk1 + 1 
     845               END DO 
     846             
     847               IF(zbhitwe == 1) THEN 
     848                 zuijk = -zdeptht(jis,jj,jk1) 
     849               ENDIF 
     850 
     851               ! integrate the pressure on the deep side 
     852               jk1 = jk  
     853               zbhitwe = 0 
     854               DO WHILE ( -zdeptht(jid,jj,jk1) < zuijk ) 
     855                 IF( jk1 == 1 ) THEN 
     856                   zbhitwe = 1 
     857                   EXIT 
     858                 ENDIF 
     859                 zdeps = MAX(zdeptht(jid,jj,jk1-1), -zuijk) 
     860                 zpwed = zpwed +                                        &  
     861                        integ2(zdeps,              zdeptht(jid,jj,jk1), & 
     862                               asp(jid,jj,jk1-1), bsp(jid,jj,jk1-1),  & 
     863                               csp(jid,jj,jk1-1), dsp(jid,jj,jk1-1) ) 
     864                 jk1 = jk1 - 1 
     865               END DO 
     866             
     867               IF( zbhitwe == 1 ) THEN 
     868                 zdeps = zdeptht(jid,jj,1) + MIN(zuijk, sshn(jid,jj)*znad) 
     869                 zrhdt1 = zrhh(jid,jj,1) - interp3(zdeptht(jid,jj,1), asp(jid,jj,1), & 
     870                                                 bsp(jid,jj,1),    csp(jid,jj,1), & 
     871                                                 dsp(jid,jj,1)) * zdeps 
     872                 zrhdt1 = MAX(zrhdt1, 1000._wp - rau0)        ! no lighter than fresh water 
     873                 zpwed  = zpwed + 0.5_wp * (zrhh(jid,jj,1) + zrhdt1) * zdeps 
     874               ENDIF 
     875 
     876               ! update the momentum trends in u direction 
     877 
     878               zdpdx1 = zcoef0 / e1u(ji,jj) * (zhpi(ji+1,jj,jk) - zhpi(ji,jj,jk)) 
     879               IF( lk_vvl ) THEN 
     880                 zdpdx2 = zcoef0 / e1u(ji,jj) * &  
     881                         ( REAL(jis-jid, wp) * (zpwes + zpwed) + (sshn(ji+1,jj)-sshn(ji,jj)) )  
     882                ELSE 
     883                 zdpdx2 = zcoef0 / e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed)  
     884               ENDIF 
     885 
     886               ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * & 
     887               &           umask(ji,jj,jk) * tmask(ji,jj,jk) * tmask(ji+1,jj,jk) 
     888            ENDIF 
     889   
     890            !!!!!     for v equation 
     891            IF( jk <= mbkv(ji,jj) ) THEN 
     892               IF( -zdeptht(ji,jj+1,mbkv(ji,jj)) >= -zdeptht(ji,jj,mbkv(ji,jj)) ) THEN 
     893                 jjs = jj + 1; jjd = jj 
     894               ELSE 
     895                 jjs = jj    ; jjd = jj + 1 
     896               ENDIF 
     897 
     898               ! integrate the pressure on the shallow side 
     899               jk1 = jk  
     900               zbhitns = 0 
     901               DO WHILE ( -zdeptht(ji,jjs,jk1) > zvijk ) 
     902                 IF( jk1 == mbkv(ji,jj) ) THEN 
     903                   zbhitns = 1 
     904                   EXIT 
     905                 ENDIF 
     906                 zdeps = MIN(zdeptht(ji,jjs,jk1+1), -zvijk) 
     907                 zpnss = zpnss +                                      &  
     908                        integ2(zdeptht(ji,jjs,jk1), zdeps,            & 
     909                               asp(ji,jjs,jk1),    bsp(ji,jjs,jk1), & 
     910                               csp(ji,jjs,jk1),    dsp(ji,jjs,jk1) ) 
     911                 jk1 = jk1 + 1 
     912               END DO 
     913             
     914               IF(zbhitns == 1) THEN 
     915                 zvijk = -zdeptht(ji,jjs,jk1) 
     916               ENDIF 
     917 
     918               ! integrate the pressure on the deep side 
     919               jk1 = jk  
     920               zbhitns = 0 
     921               DO WHILE ( -zdeptht(ji,jjd,jk1) < zvijk ) 
     922                 IF( jk1 == 1 ) THEN 
     923                   zbhitns = 1 
     924                   EXIT 
     925                 ENDIF 
     926                 zdeps = MAX(zdeptht(ji,jjd,jk1-1), -zvijk) 
     927                 zpnsd = zpnsd +                                        &  
     928                        integ2(zdeps,              zdeptht(ji,jjd,jk1), & 
     929                               asp(ji,jjd,jk1-1), bsp(ji,jjd,jk1-1), & 
     930                               csp(ji,jjd,jk1-1), dsp(ji,jjd,jk1-1) ) 
     931                 jk1 = jk1 - 1 
     932               END DO 
     933             
     934               IF( zbhitns == 1 ) THEN 
     935                 zdeps = zdeptht(ji,jjd,1) + MIN(zvijk, sshn(ji,jjd)*znad) 
     936                 zrhdt1 = zrhh(ji,jjd,1) - interp3(zdeptht(ji,jjd,1), asp(ji,jjd,1), & 
     937                                                 bsp(ji,jjd,1),    csp(ji,jjd,1), & 
     938                                                 dsp(ji,jjd,1) ) * zdeps 
     939                 zrhdt1 = MAX(zrhdt1, 1000._wp - rau0)        ! no lighter than fresh water 
     940                 zpnsd  = zpnsd + 0.5_wp * (zrhh(ji,jjd,1) + zrhdt1) * zdeps 
     941               ENDIF 
     942 
     943               ! update the momentum trends in v direction 
     944 
     945               zdpdy1 = zcoef0 / e2v(ji,jj) * (zhpi(ji,jj+1,jk) - zhpi(ji,jj,jk)) 
     946               IF( lk_vvl ) THEN 
     947                   zdpdy2 = zcoef0 / e2v(ji,jj) * & 
     948                           ( REAL(jjs-jjd, wp) * (zpnss + zpnsd) + (sshn(ji,jj+1)-sshn(ji,jj)) )  
     949               ELSE 
     950                   zdpdy2 = zcoef0 / e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd )  
     951               ENDIF 
     952 
     953               va(ji,jj,jk) = va(ji,jj,jk) + (zdpdy1 + zdpdy2)*& 
     954               &              vmask(ji,jj,jk)*tmask(ji,jj,jk)*tmask(ji,jj+1,jk) 
     955            ENDIF 
     956 
     957                     
     958           END DO 
     959        END DO 
     960      END DO 
     961 
     962      ! 
     963      IF( wrk_not_released(3, 3,4,5,6,7,8,9,10,11) )   & 
     964         CALL ctl_stop('dyn:hpg_prj: failed to release workspace arrays') 
     965      ! 
     966   END SUBROUTINE hpg_prj 
     967 
     968   SUBROUTINE cspline(fsp, xsp, asp, bsp, csp, dsp, polynomial_type) 
     969      !!---------------------------------------------------------------------- 
     970      !!                 ***  ROUTINE cspline  *** 
     971      !!        
     972      !! ** Purpose :   constrained cubic spline interpolation 
     973      !!           
     974      !! ** Method  :   f(x) = asp + bsp*x + csp*x^2 + dsp*x^3  
     975      !! Reference: K.W. Brodlie, A review of mehtods for curve and function 
     976      !!                          drawing, 1980 
     977      !! 
     978      !!---------------------------------------------------------------------- 
     979      IMPLICIT NONE 
     980      REAL(wp), DIMENSION(:,:,:), INTENT(in)  :: fsp, xsp           ! value and coordinate 
     981      REAL(wp), DIMENSION(:,:,:), INTENT(out) :: asp, bsp, csp, dsp ! coefficients of  
     982                                                                    ! the interpoated function 
     983      INTEGER, INTENT(in) :: polynomial_type                        ! 1: cubic spline  
     984                                                                    ! 2: Linear 
     985 
     986      ! Local Variables       
     987      INTEGER  ::   ji, jj, jk                 ! dummy loop indices 
     988      INTEGER  ::   jpi, jpj, jpkm1 
     989      REAL(wp) ::   zdf1, zdf2, zddf1, zddf2, ztmp1, ztmp2, zdxtmp 
     990      REAL(wp) ::   zdxtmp1, zdxtmp2, zalpha 
     991      REAL(wp) ::   zdf(size(fsp,3)) 
     992      !!---------------------------------------------------------------------- 
     993 
     994      jpi   = size(fsp,1) 
     995      jpj   = size(fsp,2) 
     996      jpkm1 = size(fsp,3) - 1 
     997 
     998      ! Clean output arrays 
     999      asp = 0.0_wp 
     1000      bsp = 0.0_wp 
     1001      csp = 0.0_wp 
     1002      dsp = 0.0_wp 
     1003       
     1004      DO ji = 1, jpi 
     1005        DO jj = 1, jpj 
     1006          IF (polynomial_type == 1) THEN     ! Constrained Cubic Spline 
     1007             DO jk = 2, jpkm1-1 
     1008                zdxtmp1 = xsp(ji,jj,jk)   - xsp(ji,jj,jk-1)   
     1009                zdxtmp2 = xsp(ji,jj,jk+1) - xsp(ji,jj,jk)   
     1010                zdf1    = ( fsp(ji,jj,jk)   - fsp(ji,jj,jk-1) ) / zdxtmp1 
     1011                zdf2    = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk)   ) / zdxtmp2 
     1012    
     1013                zalpha = ( zdxtmp1 + 2._wp * zdxtmp2 ) / ( zdxtmp1 + zdxtmp2 ) / 3._wp 
     1014               
     1015                IF(zdf1 * zdf2 <= 0._wp) THEN 
     1016                    zdf(jk) = 0._wp 
     1017                ELSE 
     1018                  zdf(jk) = zdf1 * zdf2 / ( ( 1._wp - zalpha ) * zdf1 + zalpha * zdf2 ) 
     1019                ENDIF 
     1020             END DO 
     1021 
     1022             zdf(1)     = 1.5_wp * ( fsp(ji,jj,2) - fsp(ji,jj,1) ) / & 
     1023                        &          ( xsp(ji,jj,2) - xsp(ji,jj,1) ) -  0.5_wp * zdf(2) 
     1024             zdf(jpkm1) = 1.5_wp * ( fsp(ji,jj,jpkm1) - fsp(ji,jj,jpkm1-1) ) / & 
     1025                        &          ( xsp(ji,jj,jpkm1) - xsp(ji,jj,jpkm1-1) ) - & 
     1026                        & 0.5_wp * zdf(jpkm1 - 1) 
     1027    
     1028             DO jk = 1, jpkm1 - 1 
     1029                zdxtmp = xsp(ji,jj,jk+1) - xsp(ji,jj,jk)  
     1030                ztmp1  = (zdf(jk+1) + 2._wp * zdf(jk)) / zdxtmp 
     1031                ztmp2  =  6._wp * (fsp(ji,jj,jk+1) - fsp(ji,jj,jk)) / zdxtmp / zdxtmp 
     1032                zddf1  = -2._wp * ztmp1 + ztmp2  
     1033                ztmp1  = (2._wp * zdf(jk+1) + zdf(jk)) / zdxtmp 
     1034                zddf2  =  2._wp * ztmp1 - ztmp2  
     1035    
     1036                dsp(ji,jj,jk) = (zddf2 - zddf1) / 6._wp / zdxtmp 
     1037                csp(ji,jj,jk) = ( xsp(ji,jj,jk+1) * zddf1 - xsp(ji,jj,jk)*zddf2 ) / 2._wp / zdxtmp 
     1038                bsp(ji,jj,jk) = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk) ) / zdxtmp - &  
     1039                              & csp(ji,jj,jk) * ( xsp(ji,jj,jk+1) + xsp(ji,jj,jk) ) - & 
     1040                              & dsp(ji,jj,jk) * ( xsp(ji,jj,jk+1)**2 + & 
     1041                              &                   xsp(ji,jj,jk+1) * xsp(ji,jj,jk) + & 
     1042                              &                   xsp(ji,jj,jk)**2 ) 
     1043                asp(ji,jj,jk) = fsp(ji,jj,jk) - bsp(ji,jj,jk) * xsp(ji,jj,jk) - & 
     1044                              &                 csp(ji,jj,jk) * xsp(ji,jj,jk)**2 - & 
     1045                              &                 dsp(ji,jj,jk) * xsp(ji,jj,jk)**3 
     1046             END DO 
     1047  
     1048          ELSE IF (polynomial_type == 2) THEN     ! Linear 
     1049  
     1050             DO jk = 1, jpkm1-1 
     1051                zdxtmp =xsp(ji,jj,jk+1) - xsp(ji,jj,jk)  
     1052                ztmp1 = fsp(ji,jj,jk+1) - fsp(ji,jj,jk) 
     1053    
     1054                dsp(ji,jj,jk) = 0._wp 
     1055                csp(ji,jj,jk) = 0._wp 
     1056                bsp(ji,jj,jk) = ztmp1 / zdxtmp 
     1057                asp(ji,jj,jk) = fsp(ji,jj,jk) - bsp(ji,jj,jk) * xsp(ji,jj,jk) 
     1058             END DO 
     1059 
     1060          ELSE 
     1061             CALL ctl_stop( 'invalid polynomial type in cspline' ) 
     1062          ENDIF 
     1063 
     1064        END DO 
     1065      END DO 
     1066       
     1067   END SUBROUTINE cspline 
     1068 
     1069 
     1070   FUNCTION interp1(x, xl, xr, fl, fr)  RESULT(f)  
     1071      !!---------------------------------------------------------------------- 
     1072      !!                 ***  ROUTINE interp1  *** 
     1073      !!        
     1074      !! ** Purpose :   1-d linear interpolation 
     1075      !!           
     1076      !! ** Method  :   
     1077      !!                interpolation is straight forward 
     1078      !!                extrapolation is also permitted (no value limit)  
     1079      !! 
     1080      !! H.Liu, Jan 2009,  POL  
     1081      !!---------------------------------------------------------------------- 
     1082      IMPLICIT NONE 
     1083      REAL(wp), INTENT(in) ::  x, xl, xr, fl, fr    
     1084      REAL(wp)             ::  f ! result of the interpolation (extrapolation) 
     1085      REAL(wp)             ::  zdeltx 
     1086      !!---------------------------------------------------------------------- 
     1087 
     1088      zdeltx = xr - xl 
     1089      IF(abs(zdeltx) <= 10._wp * EPSILON(x)) THEN 
     1090        f = 0.5_wp * (fl + fr) 
     1091      ELSE 
     1092        f = ( (x - xl ) * fr - ( x - xr ) * fl ) / zdeltx 
     1093      ENDIF 
     1094       
     1095   END FUNCTION interp1 
     1096 
     1097   FUNCTION interp2(x, a, b, c, d)  RESULT(f)  
     1098      !!---------------------------------------------------------------------- 
     1099      !!                 ***  ROUTINE interp1  *** 
     1100      !!        
     1101      !! ** Purpose :   1-d constrained cubic spline interpolation 
     1102      !!           
     1103      !! ** Method  :  cubic spline interpolation 
     1104      !! 
     1105      !!---------------------------------------------------------------------- 
     1106      IMPLICIT NONE 
     1107      REAL(wp), INTENT(in) ::  x, a, b, c, d    
     1108      REAL(wp)             ::  f ! value from the interpolation 
     1109      !!---------------------------------------------------------------------- 
     1110 
     1111      f = a + x* ( b + x * ( c + d * x ) )  
     1112 
     1113   END FUNCTION interp2 
     1114 
     1115 
     1116   FUNCTION interp3(x, a, b, c, d)  RESULT(f)  
     1117      !!---------------------------------------------------------------------- 
     1118      !!                 ***  ROUTINE interp1  *** 
     1119      !!        
     1120      !! ** Purpose :   Calculate the first order of deriavtive of 
     1121      !!                a cubic spline function y=a+b*x+c*x^2+d*x^3 
     1122      !!           
     1123      !! ** Method  :   f=dy/dx=b+2*c*x+3*d*x^2 
     1124      !! 
     1125      !!---------------------------------------------------------------------- 
     1126      IMPLICIT NONE 
     1127      REAL(wp), INTENT(in) ::  x, a, b, c, d    
     1128      REAL(wp)             ::  f ! value from the interpolation 
     1129      !!---------------------------------------------------------------------- 
     1130 
     1131      f = b + x * ( 2._wp * c + 3._wp * d * x) 
     1132 
     1133   END FUNCTION interp3 
     1134 
     1135    
     1136   FUNCTION integ2(xl, xr, a, b, c, d)  RESULT(f)  
     1137      !!---------------------------------------------------------------------- 
     1138      !!                 ***  ROUTINE interp1  *** 
     1139      !!        
     1140      !! ** Purpose :   1-d constrained cubic spline integration 
     1141      !!           
     1142      !! ** Method  :  integrate polynomial a+bx+cx^2+dx^3 from xl to xr  
     1143      !! 
     1144      !!---------------------------------------------------------------------- 
     1145      IMPLICIT NONE 
     1146      REAL(wp), INTENT(in) ::  xl, xr, a, b, c, d    
     1147      REAL(wp)             ::  za1, za2, za3       
     1148      REAL(wp)             ::  f                   ! integration result 
     1149      !!---------------------------------------------------------------------- 
     1150 
     1151      za1 = 0.5_wp * b  
     1152      za2 = c / 3.0_wp  
     1153      za3 = 0.25_wp * d  
     1154 
     1155      f  = xr * ( a + xr * ( za1 + xr * ( za2 + za3 * xr ) ) ) - & 
     1156         & xl * ( a + xl * ( za1 + xl * ( za2 + za3 * xl ) ) ) 
     1157 
     1158   END FUNCTION integ2 
     1159 
    10361160 
    10371161   !!====================================================================== 

  • branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90

    r2977 r3116  
    3333   USE obcdyn_bt       ! 2D open boundary condition for momentum (obc_dyn_bt routine) 
    3434   USE obcvol          ! ocean open boundary condition (obc_vol routines) 
    35    USE bdy_oce         ! unstructured open boundary conditions 
    36    USE bdydta          ! unstructured open boundary conditions 
    37    USE bdydyn          ! unstructured open boundary conditions 
     35   USE bdy_oce         ! ocean open boundary conditions 
     36   USE bdydta          ! ocean open boundary conditions 
     37   USE bdydyn          ! ocean open boundary conditions 
     38   USE bdyvol          ! ocean open boundary condition (bdy_vol routines) 
    3839   USE in_out_manager  ! I/O manager 
    3940   USE lbclnk          ! lateral boundary condition (or mpp link) 
    … …  
    7778      !!              * Apply lateral boundary conditions on after velocity  
    7879      !!             at the local domain boundaries through lbc_lnk call, 
    79       !!             at the radiative open boundaries (lk_obc=T), 
    80       !!             at the relaxed   open boundaries (lk_bdy=T), and 
     80      !!             at the one-way open boundaries (lk_obc=T), 
    8181      !!             at the AGRIF zoom     boundaries (lk_agrif=T) 
    8282      !! 
    … …  
    9292      !!               un,vn   now horizontal velocity of next time-step 
    9393      !!---------------------------------------------------------------------- 
    94       USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released 
    9594      USE oce     , ONLY:   tsa             ! tsa used as 2 3D workspace 
    96       USE wrk_nemo, ONLY:   zs_t   => wrk_2d_1 , zs_u_1 => wrk_2d_2 , zs_v_1 => wrk_2d_3 
    9795      ! 
    9896      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index 
    9997      ! 
    10098      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     99      INTEGER  ::   iku, ikv     ! local integers 
    101100#if ! defined key_dynspg_flt 
    102101      REAL(wp) ::   z2dt         ! temporary scalar 
    103102#endif 
    104       REAL(wp) ::   zue3a, zue3n, zue3b, zuf    ! local scalars 
    105       REAL(wp) ::   zve3a, zve3n, zve3b, zvf    !   -      - 
    106       REAL(wp) ::   zec, zv_t_ij, zv_t_ip1j, zv_t_ijp1 
     103      REAL(wp) ::   zue3a, zue3n, zue3b, zuf, zec   ! local scalars 
     104      REAL(wp) ::   zve3a, zve3n, zve3b, zvf        !   -      - 
    107105      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ze3u_f, ze3v_f  
    108106      !!---------------------------------------------------------------------- 
    109107 
    110       IF( wrk_in_use(2, 1,2,3) ) THEN 
    111          CALL ctl_stop('dyn_nxt: requested workspace arrays unavailable')   ;   RETURN 
    112       ENDIF 
    113108      ! 
    114109      ze3u_f => tsa(:,:,:,1)  
    … …  
    178173      ENDIF 
    179174      ! 
    180 # elif defined key_bdy  
     175# elif defined key_bdy 
    181176      !                                !* BDY open boundaries 
    182       IF( .NOT. lk_dynspg_flt ) THEN 
    183          CALL bdy_dyn_frs( kt ) 
    184 #  if ! defined key_vvl 
    185          ua_e(:,:) = 0.e0 
    186          va_e(:,:) = 0.e0 
    187          ! Set these variables for use in bdy_dyn_fla 
    188          hur_e(:,:) = hur(:,:) 
    189          hvr_e(:,:) = hvr(:,:) 
    190          DO jk = 1, jpkm1   !! Vertically integrated momentum trends 
    191             ua_e(:,:) = ua_e(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk) 
    192             va_e(:,:) = va_e(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk) 
    193          END DO 
    194          ua_e(:,:) = ua_e(:,:) * hur(:,:) 
    195          va_e(:,:) = va_e(:,:) * hvr(:,:) 
    196          DO jk = 1 , jpkm1 
    197             ua(:,:,jk) = ua(:,:,jk) - ua_e(:,:) 
    198             va(:,:,jk) = va(:,:,jk) - va_e(:,:) 
    199          END DO 
    200          CALL bdy_dta_fla( kt+1, 0,2*nn_baro) 
    201          CALL bdy_dyn_fla( sshn_b ) 
    202          CALL lbc_lnk( ua_e, 'U', -1. )   ! Boundary points should be updated 
    203          CALL lbc_lnk( va_e, 'V', -1. )   ! 
    204          DO jk = 1 , jpkm1 
    205             ua(:,:,jk) = ( ua(:,:,jk) + ua_e(:,:) ) * umask(:,:,jk) 
    206             va(:,:,jk) = ( va(:,:,jk) + va_e(:,:) ) * vmask(:,:,jk) 
    207          END DO 
    208 #  endif 
    209       ENDIF 
     177      IF( lk_dynspg_exp ) CALL bdy_dyn( kt ) 
     178      IF( lk_dynspg_ts )  CALL bdy_dyn( kt, dyn3d_only=.true. ) 
     179 
     180!!$   Do we need a call to bdy_vol here?? 
     181      ! 
    210182# endif 
    211183      ! 
    … …  
    242214         ELSE                             ! Variable volume ! 
    243215            !                             ! ================! 
    244             ! Before scale factor at t-points 
    245             ! ------------------------------- 
    246             DO jk = 1, jpkm1 
     216            ! 
     217            DO jk = 1, jpkm1                 ! Before scale factor at t-points 
    247218               fse3t_b(:,:,jk) = fse3t_n(:,:,jk)                                   & 
    248219                  &              + atfp * (  fse3t_b(:,:,jk) + fse3t_a(:,:,jk)     & 
    249                   &                         - 2.e0 * fse3t_n(:,:,jk)            ) 
    250             ENDDO 
    251             ! Add volume filter correction only at the first level of t-point scale factors 
    252             zec = atfp * rdt / rau0 
     220                  &                         - 2._wp * fse3t_n(:,:,jk)            ) 
     221            END DO 
     222            zec = atfp * rdt / rau0          ! Add filter correction only at the 1st level of t-point scale factors 
    253223            fse3t_b(:,:,1) = fse3t_b(:,:,1) - zec * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1) 
    254             ! surface at t-points and inverse surface at (u/v)-points used in surface averaging computations 
    255             zs_t  (:,:) =       e1t(:,:) * e2t(:,:) 
    256             zs_u_1(:,:) = 0.5 / ( e1u(:,:) * e2u(:,:) ) 
    257             zs_v_1(:,:) = 0.5 / ( e1v(:,:) * e2v(:,:) ) 
    258224            ! 
    259             IF( ln_dynadv_vec ) THEN 
    260                ! Before scale factor at (u/v)-points 
    261                ! ----------------------------------- 
    262                ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points 
    263                DO jk = 1, jpkm1 
    264                   DO jj = 1, jpjm1 
    265                      DO ji = 1, jpim1 
    266                         zv_t_ij           = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk) 
    267                         zv_t_ip1j         = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk) 
    268                         zv_t_ijp1         = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk) 
    269                         fse3u_b(ji,jj,jk) = umask(ji,jj,jk) * ( zs_u_1(ji,jj) * ( zv_t_ij + zv_t_ip1j ) - fse3u_0(ji,jj,jk) ) 
    270                         fse3v_b(ji,jj,jk) = vmask(ji,jj,jk) * ( zs_v_1(ji,jj) * ( zv_t_ij + zv_t_ijp1 ) - fse3v_0(ji,jj,jk) ) 
    271                      END DO 
    272                   END DO 
    273                END DO 
    274                ! lateral boundary conditions 
    275                CALL lbc_lnk( fse3u_b(:,:,:), 'U', 1. ) 
    276                CALL lbc_lnk( fse3v_b(:,:,:), 'V', 1. ) 
    277                ! Add initial scale factor to scale factor anomaly 
    278                fse3u_b(:,:,:) = fse3u_b(:,:,:) + fse3u_0(:,:,:) 
    279                fse3v_b(:,:,:) = fse3v_b(:,:,:) + fse3v_0(:,:,:) 
    280                ! Leap-Frog - Asselin filter and swap: applied on velocity 
    281                ! ----------------------------------- 
    282                DO jk = 1, jpkm1 
    283                   DO jj = 1, jpj 
     225            IF( ln_dynadv_vec ) THEN         ! vector invariant form (no thickness weighted calulation) 
     226               ! 
     227               !                                      ! before scale factors at u- & v-pts (computed from fse3t_b) 
     228               CALL dom_vvl_2( kt, fse3u_b(:,:,:), fse3v_b(:,:,:) ) 
     229               ! 
     230               DO jk = 1, jpkm1                       ! Leap-Frog - Asselin filter and swap: applied on velocity 
     231                  DO jj = 1, jpj                      !                                                 -------- 
    284232                     DO ji = 1, jpi 
    285233                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) ) 
    … …  
    294242               END DO 
    295243               ! 
    296             ELSE 
    297                ! Temporary filered scale factor at (u/v)-points (will become before scale factor) 
    298                !----------------------------------------------- 
    299                ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points 
    300                DO jk = 1, jpkm1 
    301                   DO jj = 1, jpjm1 
    302                      DO ji = 1, jpim1 
    303                         zv_t_ij          = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk) 
    304                         zv_t_ip1j        = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk) 
    305                         zv_t_ijp1        = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk) 
    306                         ze3u_f(ji,jj,jk) = umask(ji,jj,jk) * ( zs_u_1(ji,jj) * ( zv_t_ij + zv_t_ip1j ) - fse3u_0(ji,jj,jk) ) 
    307                         ze3v_f(ji,jj,jk) = vmask(ji,jj,jk) * ( zs_v_1(ji,jj) * ( zv_t_ij + zv_t_ijp1 ) - fse3v_0(ji,jj,jk) ) 
    308                      END DO 
    309                   END DO 
    310                END DO 
    311                ! lateral boundary conditions 
    312                CALL lbc_lnk( ze3u_f, 'U', 1. ) 
    313                CALL lbc_lnk( ze3v_f, 'V', 1. ) 
    314                ! Add initial scale factor to scale factor anomaly 
    315                ze3u_f(:,:,:) = ze3u_f(:,:,:) + fse3u_0(:,:,:) 
    316                ze3v_f(:,:,:) = ze3v_f(:,:,:) + fse3v_0(:,:,:) 
    317                ! Leap-Frog - Asselin filter and swap: applied on thickness weighted velocity 
    318                ! -----------------------------------             =========================== 
    319                DO jk = 1, jpkm1 
    320                   DO jj = 1, jpj 
    321                      DO ji = 1, jpim1 
     244            ELSE                             ! flux form (thickness weighted calulation) 
     245               ! 
     246               CALL dom_vvl_2( kt, ze3u_f, ze3v_f )   ! before scale factors at u- & v-pts (computed from fse3t_b) 
     247               ! 
     248               DO jk = 1, jpkm1                       ! Leap-Frog - Asselin filter and swap:  
     249                  DO jj = 1, jpj                      !                   applied on thickness weighted velocity 
     250                     DO ji = 1, jpim1                 !                              --------------------------- 
    322251                        zue3a = ua(ji,jj,jk) * fse3u_a(ji,jj,jk) 
    323252                        zve3a = va(ji,jj,jk) * fse3v_a(ji,jj,jk) 
    … …  
    327256                        zve3b = vb(ji,jj,jk) * fse3v_b(ji,jj,jk) 
    328257                        ! 
    329                         zuf  = ( zue3n + atfp * ( zue3b - 2.e0 * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk) 
    330                         zvf  = ( zve3n + atfp * ( zve3b - 2.e0 * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk) 
     258                        zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk) 
     259                        zvf = ( zve3n + atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk) 
    331260                        ! 
    332                         ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity 
     261                        ub(ji,jj,jk) = zuf                     ! ub <-- filtered velocity 
    333262                        vb(ji,jj,jk) = zvf 
    334                         un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua 
     263                        un(ji,jj,jk) = ua(ji,jj,jk)            ! un <-- ua 
    335264                        vn(ji,jj,jk) = va(ji,jj,jk) 
    336265                     END DO 
    337266                  END DO 
    338267               END DO 
    339                fse3u_b(:,:,:) = ze3u_f(:,:,:)                   ! e3u_b <-- filtered scale factor 
    340                fse3v_b(:,:,:) = ze3v_f(:,:,:) 
    341                CALL lbc_lnk( ub, 'U', -1. )                     ! lateral boundary conditions 
     268               fse3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1)      ! e3u_b <-- filtered scale factor 
     269               fse3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1) 
     270               CALL lbc_lnk( ub, 'U', -1. )                    ! lateral boundary conditions 
    342271               CALL lbc_lnk( vb, 'V', -1. ) 
    343272            ENDIF 
    … …  
    350279         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask ) 
    351280      !  
    352       IF( wrk_not_released(2, 1,2,3) )   CALL ctl_stop('dyn_nxt: failed to release workspace arrays') 
    353       ! 
    354281   END SUBROUTINE dyn_nxt 
    355282 

  • branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90

    r2715 r3116  
    1515   USE dom_oce        ! ocean space and time domain variables 
    1616   USE phycst         ! physical constants 
    17    USE obc_oce        ! ocean open boundary conditions 
    1817   USE sbc_oce        ! surface boundary condition: ocean 
    1918   USE sbcapr         ! surface boundary condition: atmospheric pressure 
    … …  
    222221      ENDIF 
    223222 
    224 #if defined key_obc 
    225       !                        ! Conservation of ocean volume (key_dynspg_flt) 
    226       IF( lk_dynspg_flt )   ln_vol_cst = .true. 
    227  
    228       !                        ! Application of Flather's algorithm at open boundaries 
    229       IF( lk_dynspg_flt )   ln_obc_fla = .false. 
    230       IF( lk_dynspg_exp )   ln_obc_fla = .true. 
    231       IF( lk_dynspg_ts  )   ln_obc_fla = .true. 
    232 #endif 
    233223      ! 
    234224   END SUBROUTINE dyn_spg_init 

  • branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_exp.F90

    r2715 r3116  
    2121   USE phycst          ! physical constants 
    2222   USE obc_par         ! open boundary condition parameters 
    23    USE obcdta          ! open boundary condition data     (obc_dta_bt routine) 
     23   USE obcdta          ! open boundary condition data     (bdy_dta_bt routine) 
    2424   USE in_out_manager  ! I/O manager 
    2525   USE lib_mpp         ! distributed memory computing library 

  • branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_flt.F90

    r2977 r3116  
    2626   USE sbc_oce         ! surface boundary condition: ocean 
    2727   USE obc_oce         ! Lateral open boundary condition 
     28   USE bdy_oce         ! Lateral open boundary condition 
    2829   USE sol_oce         ! ocean elliptic solver 
    2930   USE phycst          ! physical constants 
    … …  
    3334   USE solpcg          ! preconditionned conjugate gradient solver 
    3435   USE solsor          ! Successive Over-relaxation solver 
    35    USE obcdyn          ! ocean open boundary condition (obc_dyn routines) 
    36    USE obcvol          ! ocean open boundary condition (obc_vol routines) 
    37    USE bdy_oce         ! Unstructured open boundaries condition 
    38    USE bdydyn          ! Unstructured open boundaries condition (bdy_dyn routine)  
    39    USE bdyvol          ! Unstructured open boundaries condition (bdy_vol routine) 
     36   USE obcdyn          ! ocean open boundary condition on dynamics 
     37   USE obcvol          ! ocean open boundary condition (obc_vol routine) 
     38   USE bdydyn          ! ocean open boundary condition on dynamics 
     39   USE bdyvol          ! ocean open boundary condition (bdy_vol routine) 
    4040   USE cla             ! cross land advection 
    4141   USE in_out_manager  ! I/O manager 
    … …  
    191191#endif 
    192192#if defined key_bdy 
    193       CALL bdy_dyn_frs( kt )       ! Update velocities on unstructured boundary using the Flow Relaxation Scheme 
    194       CALL bdy_vol( kt )           ! Correction of the barotropic component velocity to control the volume of the system 
     193      CALL bdy_dyn( kt )      ! Update velocities on each open boundary 
     194      CALL bdy_vol( kt )      ! Correction of the barotropic component velocity to control the volume of the system 
    195195#endif 
    196196#if defined key_agrif 
    … …  
    308308#if defined key_obc 
    309309            ! caution : grad D = 0 along open boundaries 
     310            ! Remark: The filtering force could be reduced here in the FRS zone 
     311            !         by multiplying spgu/spgv by (1-alpha) ??   
    310312            spgu(ji,jj) = z2dt * ztdgu * obcumask(ji,jj) 
    311313            spgv(ji,jj) = z2dt * ztdgv * obcvmask(ji,jj) 
    312314#elif defined key_bdy 
    313315            ! caution : grad D = 0 along open boundaries 
    314             ! Remark: The filtering force could be reduced here in the FRS zone 
    315             !         by multiplying spgu/spgv by (1-alpha) ??   
    316316            spgu(ji,jj) = z2dt * ztdgu * bdyumask(ji,jj) 
    317             spgv(ji,jj) = z2dt * ztdgv * bdyvmask(ji,jj)            
     317            spgv(ji,jj) = z2dt * ztdgv * bdyvmask(ji,jj) 
    318318#else 
    319319            spgu(ji,jj) = z2dt * ztdgu 

  • branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_oce.F90

    r2715 r3116  
    3434 
    3535  !                                                                         !!! Time splitting scheme (key_dynspg_ts)  
    36    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sshn_e, ssha_e   ! sea surface heigth (now, after, average) 
    37    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ua_e  , va_e     ! barotropic velocities (after) 
    38    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hu_e  , hv_e     ! now ocean depth ( = Ho+sshn_e ) 
    39    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hur_e , hvr_e    ! inverse of hu_e and hv_e 
    40    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sshn_b       ! before field without time-filter 
     36   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:), TARGET ::   sshn_e, ssha_e   ! sea surface heigth (now, after, average) 
     37   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:), TARGET ::   ua_e  , va_e     ! barotropic velocities (after) 
     38   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   hu_e  , hv_e     ! now ocean depth ( = Ho+sshn_e ) 
     39   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:), TARGET ::   hur_e , hvr_e    ! inverse of hu_e and hv_e 
     40   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:), TARGET ::   sshn_b           ! before field without time-filter 
    4141 
    4242   !!---------------------------------------------------------------------- 

  • branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r3104 r3116  
    2525   USE domvvl          ! variable volume 
    2626   USE zdfbfr          ! bottom friction 
    27    USE obcdta          ! open boundary condition data      
    28    USE obcfla          ! Flather open boundary condition   
    2927   USE dynvor          ! vorticity term 
    3028   USE obc_oce         ! Lateral open boundary condition 
    3129   USE obc_par         ! open boundary condition parameters 
    32    USE bdy_oce         ! unstructured open boundaries 
    33    USE bdy_par         ! unstructured open boundaries 
    34    USE bdydta          ! unstructured open boundaries 
    35    USE bdydyn          ! unstructured open boundaries 
    36    USE bdytides        ! tidal forcing at unstructured open boundaries. 
     30   USE obcdta          ! open boundary condition data      
     31   USE obcfla          ! Flather open boundary condition   
     32   USE bdy_par         ! for lk_bdy 
     33   USE bdy_oce         ! Lateral open boundary condition 
     34   USE bdydta          ! open boundary condition data      
     35   USE bdydyn2d        ! open boundary conditions on barotropic variables 
    3736   USE sbctide 
    3837   USE updtide 
    … …  
    121120      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices 
    122121      INTEGER  ::   icycle           ! local scalar 
    123       REAL(wp) ::   zraur, zcoef, z2dt_e, z2dt_b     ! local scalars 
    124       REAL(wp) ::   z1_8, zx1, zy1                   !   -      - 
    125       REAL(wp) ::   z1_4, zx2, zy2                   !   -      - 
    126       REAL(wp) ::   zu_spg, zu_cor, zu_sld, zu_asp   !   -      - 
    127       REAL(wp) ::   zv_spg, zv_cor, zv_sld, zv_asp   !   -      - 
     122      INTEGER  ::   ikbu, ikbv       ! local scalar 
     123      REAL(wp) ::   zraur, zcoef, z2dt_e, z1_2dt_b, z2dt_bf   ! local scalars 
     124      REAL(wp) ::   z1_8, zx1, zy1                            !   -      - 
     125      REAL(wp) ::   z1_4, zx2, zy2                            !   -      - 
     126      REAL(wp) ::   zu_spg, zu_cor, zu_sld, zu_asp            !   -      - 
     127      REAL(wp) ::   zv_spg, zv_cor, zv_sld, zv_asp            !   -      - 
     128      REAL(wp) ::   ua_btm, va_btm                            !   -      - 
    128129      !!---------------------------------------------------------------------- 
    129130 
    … …  
    149150         hvr_e (:,:) = hvr  (:,:) 
    150151         IF( ln_dynvor_een ) THEN 
    151             ftne(1,:) = 0.e0   ;   ftnw(1,:) = 0.e0   ;   ftse(1,:) = 0.e0   ;   ftsw(1,:) = 0.e0 
     152            ftne(1,:) = 0._wp   ;   ftnw(1,:) = 0._wp   ;   ftse(1,:) = 0._wp   ;   ftsw(1,:) = 0._wp 
    152153            DO jj = 2, jpj 
    153154               DO ji = fs_2, jpi   ! vector opt. 
    154                   ftne(ji,jj) = ( ff(ji-1,jj  ) + ff(ji  ,jj  ) + ff(ji  ,jj-1) ) / 3. 
    155                   ftnw(ji,jj) = ( ff(ji-1,jj-1) + ff(ji-1,jj  ) + ff(ji  ,jj  ) ) / 3. 
    156                   ftse(ji,jj) = ( ff(ji  ,jj  ) + ff(ji  ,jj-1) + ff(ji-1,jj-1) ) / 3. 
    157                   ftsw(ji,jj) = ( ff(ji  ,jj-1) + ff(ji-1,jj-1) + ff(ji-1,jj  ) ) / 3. 
     155                  ftne(ji,jj) = ( ff(ji-1,jj  ) + ff(ji  ,jj  ) + ff(ji  ,jj-1) ) / 3._wp 
     156                  ftnw(ji,jj) = ( ff(ji-1,jj-1) + ff(ji-1,jj  ) + ff(ji  ,jj  ) ) / 3._wp 
     157                  ftse(ji,jj) = ( ff(ji  ,jj  ) + ff(ji  ,jj-1) + ff(ji-1,jj-1) ) / 3._wp 
     158                  ftsw(ji,jj) = ( ff(ji  ,jj-1) + ff(ji-1,jj-1) + ff(ji-1,jj  ) ) / 3._wp 
    158159               END DO 
    159160            END DO 
    … …  
    162163      ENDIF 
    163164 
    164       !                                   !* Local constant initialization 
    165       z2dt_b = 2.0 * rdt                                    ! baroclinic time step 
    166       z1_8 = 0.5 * 0.25                                     ! coefficient for vorticity estimates 
    167       z1_4 = 0.5 * 0.5 
    168       zraur  = 1. / rau0                                    ! 1 / volumic mass 
    169       ! 
    170       zhdiv(:,:) = 0.e0                                     ! barotropic divergence 
    171       zu_sld = 0.e0   ;   zu_asp = 0.e0                     ! tides trends (lk_tide=F) 
    172       zv_sld = 0.e0   ;   zv_asp = 0.e0 
     165      !                                                     !* Local constant initialization 
     166      z1_2dt_b = 1._wp / ( 2.0_wp * rdt )                   ! reciprocal of baroclinic time step 
     167      IF( neuler == 0 .AND. kt == nit000 )   z1_2dt_b = 1.0_wp / rdt    ! reciprocal of baroclinic  
     168                                                                        ! time step (euler timestep) 
     169      z1_8     = 0.125_wp                                   ! coefficient for vorticity estimates 
     170      z1_4     = 0.25_wp         
     171      zraur    = 1._wp / rau0                               ! 1 / volumic mass 
     172      ! 
     173      zhdiv(:,:) = 0._wp                                    ! barotropic divergence 
     174      zu_sld = 0._wp   ;   zu_asp = 0._wp                   ! tides trends (lk_tide=F) 
     175      zv_sld = 0._wp   ;   zv_asp = 0._wp 
     176 
     177      IF( kt == nit000 .AND. neuler == 0) THEN              ! for implicit bottom friction 
     178        z2dt_bf = rdt 
     179      ELSE 
     180        z2dt_bf = 2.0_wp * rdt 
     181      ENDIF 
    173182 
    174183      ! ----------------------------------------------------------------------------- 
    … …  
    178187      !                                   !* e3*d/dt(Ua), e3*Ub, e3*Vn (Vertically integrated) 
    179188      !                                   ! -------------------------- 
    180       zua(:,:) = 0.e0   ;   zun(:,:) = 0.e0   ;   ub_b(:,:) = 0.e0 
    181       zva(:,:) = 0.e0   ;   zvn(:,:) = 0.e0   ;   vb_b(:,:) = 0.e0 
     189      zua(:,:) = 0._wp   ;   zun(:,:) = 0._wp   ;   ub_b(:,:) = 0._wp 
     190      zva(:,:) = 0._wp   ;   zvn(:,:) = 0._wp   ;   vb_b(:,:) = 0._wp 
    182191      ! 
    183192      DO jk = 1, jpkm1 
    … …  
    197206               ! 
    198207#if defined key_vvl 
    199                ub_b(ji,jj) = ub_b(ji,jj) + (fse3u_0(ji,jj,jk)*(1.+sshu_b(ji,jj)*muu(ji,jj,jk)))* ub(ji,jj,jk)  
    200                vb_b(ji,jj) = vb_b(ji,jj) + (fse3v_0(ji,jj,jk)*(1.+sshv_b(ji,jj)*muv(ji,jj,jk)))* vb(ji,jj,jk)    
     208               ub_b(ji,jj) = ub_b(ji,jj) + fse3u_b(ji,jj,jk)* ub(ji,jj,jk)   *umask(ji,jj,jk)  
     209               vb_b(ji,jj) = vb_b(ji,jj) + fse3v_b(ji,jj,jk)* vb(ji,jj,jk)   *vmask(ji,jj,jk) 
    201210#else 
    202211               ub_b(ji,jj) = ub_b(ji,jj) + fse3u_0(ji,jj,jk) * ub(ji,jj,jk)  * umask(ji,jj,jk) 
    … …  
    272281      DO jj = 2, jpjm1                             ! Remove coriolis term (and possibly spg) from barotropic trend 
    273282         DO ji = fs_2, fs_jpim1 
    274             zua(ji,jj) = zua(ji,jj) - zcu(ji,jj) 
    275             zva(ji,jj) = zva(ji,jj) - zcv(ji,jj) 
    276          END DO 
     283             zua(ji,jj) = zua(ji,jj) - zcu(ji,jj) 
     284             zva(ji,jj) = zva(ji,jj) - zcv(ji,jj) 
     285          END DO 
    277286      END DO 
    278287 
    … …  
    280289      !                                             ! Remove barotropic contribution of bottom friction  
    281290      !                                             ! from the barotropic transport trend 
    282       zcoef = -1. / z2dt_b 
     291      zcoef = -1._wp * z1_2dt_b 
     292 
     293      IF(ln_bfrimp) THEN 
     294      !                                   ! Remove the bottom stress trend from 3-D sea surface level gradient 
     295      !                                   ! and Coriolis forcing in case of 3D semi-implicit bottom friction  
     296        DO jj = 2, jpjm1          
     297           DO ji = fs_2, fs_jpim1 
     298              ikbu = mbku(ji,jj) 
     299              ikbv = mbkv(ji,jj) 
     300              ua_btm = zcu(ji,jj) * z2dt_bf * hur(ji,jj) * umask (ji,jj,ikbu) 
     301              va_btm = zcv(ji,jj) * z2dt_bf * hvr(ji,jj) * vmask (ji,jj,ikbv) 
     302 
     303              zua(ji,jj) = zua(ji,jj) - bfrua(ji,jj) * ua_btm 
     304              zva(ji,jj) = zva(ji,jj) - bfrva(ji,jj) * va_btm 
     305           END DO 
     306        END DO 
     307 
     308      ELSE 
     309 
    283310# if defined key_vectopt_loop 
    284       DO jj = 1, 1 
    285          DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling) 
     311        DO jj = 1, 1 
     312           DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling) 
    286313# else 
    287       DO jj = 2, jpjm1 
    288          DO ji = 2, jpim1 
     314        DO jj = 2, jpjm1 
     315           DO ji = 2, jpim1 
    289316# endif 
    290317            ! Apply stability criteria for bottom friction 
    291318            !RBbug for vvl and external mode we may need to use varying fse3 
    292319            !!gm  Rq: the bottom e3 present the smallest variation, the use of e3u_0 is not a big approx. 
    293             zbfru(ji,jj) = MAX(  bfrua(ji,jj) , fse3u(ji,jj,mbku(ji,jj)) * zcoef  ) 
    294             zbfrv(ji,jj) = MAX(  bfrva(ji,jj) , fse3v(ji,jj,mbkv(ji,jj)) * zcoef  ) 
    295          END DO 
    296       END DO 
    297  
    298       IF( lk_vvl ) THEN 
    299          DO jj = 2, jpjm1 
    300             DO ji = fs_2, fs_jpim1   ! vector opt. 
    301                zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj)   & 
    302                   &       / ( hu_0(ji,jj) + sshu_b(ji,jj) + 1.e0 - umask(ji,jj,1) ) 
    303                zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj)   & 
    304                   &       / ( hv_0(ji,jj) + sshv_b(ji,jj) + 1.e0 - vmask(ji,jj,1) ) 
    305             END DO 
    306          END DO 
    307       ELSE 
    308          DO jj = 2, jpjm1 
    309             DO ji = fs_2, fs_jpim1   ! vector opt. 
    310                zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj) * hur(ji,jj) 
    311                zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj) * hvr(ji,jj) 
    312             END DO 
    313          END DO 
    314       ENDIF 
    315  
     320              zbfru(ji,jj) = MAX(  bfrua(ji,jj) , fse3u(ji,jj,mbku(ji,jj)) * zcoef  ) 
     321              zbfrv(ji,jj) = MAX(  bfrva(ji,jj) , fse3v(ji,jj,mbkv(ji,jj)) * zcoef  ) 
     322           END DO 
     323        END DO 
     324 
     325        IF( lk_vvl ) THEN 
     326           DO jj = 2, jpjm1 
     327              DO ji = fs_2, fs_jpim1   ! vector opt. 
     328                 zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj)   & 
     329                    &       / ( hu_0(ji,jj) + sshu_b(ji,jj) + 1._wp - umask(ji,jj,1) ) 
     330                 zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj)   & 
     331                    &       / ( hv_0(ji,jj) + sshv_b(ji,jj) + 1._wp - vmask(ji,jj,1) ) 
     332              END DO 
     333           END DO 
     334        ELSE 
     335           DO jj = 2, jpjm1 
     336              DO ji = fs_2, fs_jpim1   ! vector opt. 
     337                 zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj) * hur(ji,jj) 
     338                 zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj) * hvr(ji,jj) 
     339              END DO 
     340           END DO 
     341        ENDIF 
     342      END IF    ! end (ln_bfrimp) 
     343 
     344                     
    316345      !                                   !* d/dt(Ua), Ub, Vn (Vertical mean velocity) 
    317346      !                                   ! --------------------------  
    … …  
    320349      ! 
    321350      IF( lk_vvl ) THEN 
    322          ub_b(:,:) = ub_b(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1.e0 - umask(:,:,1) ) 
    323          vb_b(:,:) = vb_b(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1.e0 - vmask(:,:,1) ) 
     351         ub_b(:,:) = ub_b(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1._wp - umask(:,:,1) ) 
     352         vb_b(:,:) = vb_b(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1._wp - vmask(:,:,1) ) 
    324353      ELSE 
    325354         ub_b(:,:) = ub_b(:,:) * hur(:,:) 
    … …  
    357386      ! set ssh corrections to 0 
    358387      ! ssh corrections are applied to normal velocities (Flather's algorithm) and averaged over the barotropic loop 
    359       IF( lp_obc_east  )   sshfoe_b(:,:) = 0.e0 
    360       IF( lp_obc_west  )   sshfow_b(:,:) = 0.e0 
    361       IF( lp_obc_south )   sshfos_b(:,:) = 0.e0 
    362       IF( lp_obc_north )   sshfon_b(:,:) = 0.e0 
     388      IF( lp_obc_east  )   sshfoe_b(:,:) = 0._wp 
     389      IF( lp_obc_west  )   sshfow_b(:,:) = 0._wp 
     390      IF( lp_obc_south )   sshfos_b(:,:) = 0._wp 
     391      IF( lp_obc_north )   sshfon_b(:,:) = 0._wp 
    363392#endif 
    364393 
    … …  
    369398         IF( jn == 1 )   z2dt_e = rdt / nn_baro 
    370399 
    371          !                                                !* Update the forcing (OBC, BDY and tides) 
     400         !                                                !* Update the forcing (BDY and tides) 
    372401         !                                                !  ------------------ 
    373402         IF( lk_obc )   CALL obc_dta_bt ( kt, jn   ) 
    374          IF( lk_bdy )   CALL bdy_dta_fla( kt, jn+1, icycle ) 
     403         IF( lk_bdy )   CALL bdy_dta ( kt, jit=jn, time_offset=+1 ) 
    375404         IF ( ln_tide_pot ) CALL upd_tide( kt, jn ) 
    376405 
    377406         !                                                !* after ssh_e 
    378407         !                                                !  ----------- 
    379          DO jj = 2, jpjm1                                      ! Horizontal divergence of barotropic transports 
     408         DO jj = 2, jpjm1                                 ! Horizontal divergence of barotropic transports 
    380409            DO ji = fs_2, fs_jpim1   ! vector opt. 
    381410               zhdiv(ji,jj) = (   e2u(ji  ,jj) * zun_e(ji  ,jj) * hu_e(ji  ,jj)     & 
    … …  
    389418         !                                                     ! OBC : zhdiv must be zero behind the open boundary 
    390419!!  mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column 
    391          IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1  ) = 0.e0      ! east 
    392          IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1  ) = 0.e0      ! west 
    393          IF( lp_obc_north )   zhdiv(nin0  :nin1  ,njn0p1:njn1p1) = 0.e0      ! north 
    394          IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1  ) = 0.e0      ! south 
     420         IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1  ) = 0._wp      ! east 
     421         IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1  ) = 0._wp      ! west 
     422         IF( lp_obc_north )   zhdiv(nin0  :nin1  ,njn0p1:njn1p1) = 0._wp      ! north 
     423         IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1  ) = 0._wp      ! south 
    395424#endif 
    396425#if defined key_bdy 
    … …  
    406435         !                                                !* after barotropic velocities (vorticity scheme dependent) 
    407436         !                                                !  ---------------------------   
    408          zwx(:,:) = e2u(:,:) * zun_e(:,:) * hu_e(:,:)           ! now_e transport 
     437         zwx(:,:) = e2u(:,:) * zun_e(:,:) * hu_e(:,:)     ! now_e transport 
    409438         zwy(:,:) = e1v(:,:) * zvn_e(:,:) * hv_e(:,:) 
    410439         ! 
    … …  
    435464                  zv_cor =-z1_4 * ( ff(ji-1,jj  ) * zx1 + ff(ji,jj) * zx2 ) * hvr_e(ji,jj) 
    436465                  ! after velocities with implicit bottom friction 
    437                   ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj) ) ) * umask(ji,jj,1)   & 
    438                      &         / ( 1.e0         - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 
    439                   va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj) ) ) * vmask(ji,jj,1)   & 
    440                      &         / ( 1.e0         - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 
     466 
     467                  IF( ln_bfrimp ) THEN      ! implicit bottom friction 
     468                     !   A new method to implement the implicit bottom friction.  
     469                     !   H. Liu 
     470                     !   Sept 2011 
     471                     ua_e(ji,jj) = umask(ji,jj,1) * ( zub_e(ji,jj) +                                            & 
     472                      &                               z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp )            & 
     473                      &                               / ( 1._wp      - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) ) 
     474                     ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e *   zua(ji,jj)  ) * umask(ji,jj,1)    
     475                     ! 
     476                     va_e(ji,jj) = vmask(ji,jj,1) * ( zvb_e(ji,jj) +                                            & 
     477                      &                               z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp )            & 
     478                      &                               / ( 1._wp      - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) ) 
     479                     va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e *   zva(ji,jj)  ) * vmask(ji,jj,1)    
     480                     ! 
     481                  ELSE 
     482                     ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj) ) ) * umask(ji,jj,1)   & 
     483                      &           / ( 1._wp         - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 
     484                     va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj) ) ) * vmask(ji,jj,1)   & 
     485                      &           / ( 1._wp         - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 
     486                  ENDIF 
    441487               END DO 
    442488            END DO 
    … …  
    466512                  zv_cor  = zx1 * ( ff(ji-1,jj  ) + ff(ji,jj) ) * hvr_e(ji,jj) 
    467513                  ! after velocities with implicit bottom friction 
    468                   ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj) ) ) * umask(ji,jj,1)   & 
    469                      &         / ( 1.e0         - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 
    470                   va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj) ) ) * vmask(ji,jj,1)   & 
    471                      &         / ( 1.e0         - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 
     514                  IF( ln_bfrimp ) THEN 
     515                     !   A new method to implement the implicit bottom friction.  
     516                     !   H. Liu 
     517                     !   Sept 2011 
     518                     ua_e(ji,jj) = umask(ji,jj,1) * ( zub_e(ji,jj) +                                            & 
     519                      &                               z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp )            & 
     520                      &                               / ( 1._wp      - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) ) 
     521                     ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e *   zua(ji,jj)  ) * umask(ji,jj,1)    
     522                     ! 
     523                     va_e(ji,jj) = vmask(ji,jj,1) * ( zvb_e(ji,jj) +                                            & 
     524                      &                               z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp )            & 
     525                      &                               / ( 1._wp      - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) ) 
     526                     va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e *   zva(ji,jj)  ) * vmask(ji,jj,1)    
     527                     ! 
     528                  ELSE 
     529                     ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj) ) ) * umask(ji,jj,1)   & 
     530                     &            / ( 1._wp        - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 
     531                     va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj) ) ) * vmask(ji,jj,1)   & 
     532                     &            / ( 1._wp        - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 
     533                  ENDIF 
    472534               END DO 
    473535            END DO 
    … …  
    497559                     &                           + ftnw(ji,jj  ) * zwx(ji-1,jj  ) + ftne(ji,jj  ) * zwx(ji  ,jj  ) ) * hvr_e(ji,jj) 
    498560                  ! after velocities with implicit bottom friction 
    499                   ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj) ) ) * umask(ji,jj,1)   & 
    500                      &         / ( 1.e0         - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 
    501                   va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj) ) ) * vmask(ji,jj,1)   & 
    502                      &         / ( 1.e0         - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 
     561                  IF( ln_bfrimp ) THEN 
     562                     !   A new method to implement the implicit bottom friction.  
     563                     !   H. Liu 
     564                     !   Sept 2011 
     565                     ua_e(ji,jj) = umask(ji,jj,1) * ( zub_e(ji,jj) +                                            & 
     566                      &                               z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp )            & 
     567                      &                               / ( 1._wp      - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) ) 
     568                     ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e *   zua(ji,jj)  ) * umask(ji,jj,1)    
     569                     ! 
     570                     va_e(ji,jj) = vmask(ji,jj,1) * ( zvb_e(ji,jj) +                                            & 
     571                      &                               z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp )            & 
     572                      &                               / ( 1._wp      - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) ) 
     573                     va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e *   zva(ji,jj)  ) * vmask(ji,jj,1)    
     574                     ! 
     575                  ELSE 
     576                     ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj) ) ) * umask(ji,jj,1)   & 
     577                     &            / ( 1._wp        - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 
     578                     va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj) ) ) * vmask(ji,jj,1)   & 
     579                     &            / ( 1._wp        - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 
     580                  ENDIF 
    503581               END DO 
    504582            END DO 
    505583            !  
    506584         ENDIF 
    507          !                                                !* domain lateral boundary 
    508          !                                                !  ----------------------- 
    509          !                                                      ! Flather's boundary condition for the barotropic loop : 
    510          !                                                      !         - Update sea surface height on each open boundary 
    511          !                                                      !         - Correct the velocity 
    512  
     585         !                                                     !* domain lateral boundary 
     586         !                                                     !  ----------------------- 
     587 
     588                                                               ! OBC open boundaries 
    513589         IF( lk_obc               )   CALL obc_fla_ts ( ua_e, va_e, sshn_e, ssha_e ) 
    514          IF( lk_bdy .OR. ln_tides )   CALL bdy_dyn_fla( sshn_e )  
     590 
     591                                                               ! BDY open boundaries 
     592#if defined key_bdy 
     593         pssh => sshn_e 
     594         phur => hur_e 
     595         phvr => hvr_e 
     596         pu2d => ua_e 
     597         pv2d => va_e 
     598 
     599         IF( lk_bdy )   CALL bdy_dyn2d( kt )  
     600#endif 
     601 
    515602         ! 
    516603         CALL lbc_lnk( ua_e  , 'U', -1. )                      ! local domain boundaries  
    … …  
    544631            DO jj = 1, jpjm1                                    ! Sea Surface Height at u- & v-points 
    545632               DO ji = 1, fs_jpim1   ! Vector opt. 
    546                   zsshun_e(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )       & 
    547                      &                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn_e(ji  ,jj)    & 
    548                      &                    + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn_e(ji+1,jj) ) 
    549                   zsshvn_e(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )       & 
    550                      &                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn_e(ji,jj  )    & 
    551                      &                    + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn_e(ji,jj+1) ) 
     633                  zsshun_e(ji,jj) = 0.5_wp * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )       & 
     634                     &                     * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn_e(ji  ,jj)    & 
     635                     &                     +  e1t(ji+1,jj) * e2t(ji+1,jj) * sshn_e(ji+1,jj) ) 
     636                  zsshvn_e(ji,jj) = 0.5_wp * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )       & 
     637                     &                     * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn_e(ji,jj  )    & 
     638                     &                     +  e1t(ji,jj+1) * e2t(ji,jj+1) * sshn_e(ji,jj+1) ) 
    552639               END DO 
    553640            END DO 
    … …  
    557644            hu_e (:,:) = hu_0(:,:) + zsshun_e(:,:)              ! Ocean depth at U- and V-points 
    558645            hv_e (:,:) = hv_0(:,:) + zsshvn_e(:,:) 
    559             hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1.e0 - umask(:,:,1) ) 
    560             hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1.e0 - vmask(:,:,1) ) 
     646            hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1._wp - umask(:,:,1) ) 
     647            hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1._wp - vmask(:,:,1) ) 
    561648            ! 
    562649         ENDIF 
    … …  
    577664      ! 
    578665      !                                   !* Time average ==> after barotropic u, v, ssh 
    579       zcoef =  1.e0 / ( 2 * nn_baro  + 1 )  
     666      zcoef =  1._wp / ( 2 * nn_baro  + 1 )  
    580667      zu_sum(:,:) = zcoef * zu_sum  (:,:)  
    581668      zv_sum(:,:) = zcoef * zv_sum  (:,:)  
    … …  
    583670      !                                   !* update the general momentum trend 
    584671      DO jk=1,jpkm1 
    585          ua(:,:,jk) = ua(:,:,jk) + ( zu_sum(:,:) - ub_b(:,:) ) / z2dt_b 
    586          va(:,:,jk) = va(:,:,jk) + ( zv_sum(:,:) - vb_b(:,:) ) / z2dt_b 
     672         ua(:,:,jk) = ua(:,:,jk) + ( zu_sum(:,:) - ub_b(:,:) ) * z1_2dt_b 
     673         va(:,:,jk) = va(:,:,jk) + ( zv_sum(:,:) - vb_b(:,:) ) * z1_2dt_b 
    587674      END DO 
    588675      un_b  (:,:) =  zu_sum(:,:)  
    … …  
    618705            CALL iom_get( numror, jpdom_autoglo, 'vn_b'  , vn_b  (:,:) )   ! from barotropic loop 
    619706         ELSE 
    620             un_b (:,:) = 0.e0 
    621             vn_b (:,:) = 0.e0 
     707            un_b (:,:) = 0._wp 
     708            vn_b (:,:) = 0._wp 
    622709            ! vertical sum 
    623710            IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll 
    … …  
    640727         ! Vertically integrated velocity (before) 
    641728         IF (neuler/=0) THEN 
    642             ub_b (:,:) = 0.e0 
    643             vb_b (:,:) = 0.e0 
     729            ub_b (:,:) = 0._wp 
     730            vb_b (:,:) = 0._wp 
    644731 
    645732            ! vertical sum 
    … …  
    659746 
    660747            IF( lk_vvl ) THEN 
    661                ub_b (:,:) = ub_b(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1.e0 - umask(:,:,1) ) 
    662                vb_b (:,:) = vb_b(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1.e0 - vmask(:,:,1) ) 
     748               ub_b (:,:) = ub_b(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1._wp - umask(:,:,1) ) 
     749               vb_b (:,:) = vb_b(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1._wp - vmask(:,:,1) ) 
    663750            ELSE 
    664751               ub_b(:,:) = ub_b(:,:) * hur(:,:) 

  • branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90

    r2977 r3116  
    2020   USE in_out_manager  ! I/O manager 
    2121   USE lib_mpp         ! MPP library 
     22   USE zdfbfr          ! bottom friction setup 
    2223 
    2324   IMPLICIT NONE 
    … …  
    6162      REAL(wp), INTENT(in) ::  p2dt   ! vertical profile of tracer time-step 
    6263      !! 
    63       INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    64       REAL(wp) ::   z1_p2dt, zcoef, zzwi, zzws, zrhs   ! local scalars 
     64      INTEGER  ::   ji, jj, jk     ! dummy loop indices 
     65      INTEGER  ::   ikbum1, ikbvm1 ! local variable 
     66      REAL(wp) ::   z1_p2dt, z2dtf, zcoef, zzwi, zzws, zrhs ! local scalars 
     67 
     68      !! * Local variables for implicit bottom friction.    H. Liu 
     69      REAL(wp) ::   zbfru, zbfrv  
     70      REAL(wp) ::   zbfr_imp = 0._wp                        ! toggle (SAVE'd by assignment)  
    6571      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwd, zws 
     72      !!---------------------------------------------------------------------- 
    6673      !!---------------------------------------------------------------------- 
    6774 
    … …  
    7784         IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator' 
    7885         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ' 
     86         IF(ln_bfrimp) zbfr_imp = 1._wp 
    7987      ENDIF 
    8088 
    … …  
    8492 
    8593      ! 1. Vertical diffusion on u 
     94 
     95      ! Vertical diffusion on u&v 
    8696      ! --------------------------- 
    8797      ! Matrix and second member construction 
    88       ! bottom boundary condition: both zwi and zws must be masked as avmu can take 
    89       ! non zero value at the ocean bottom depending on the bottom friction 
    90       ! used but the bottom velocities have already been updated with the bottom 
    91       ! friction velocity in dyn_bfr using values from the previous timestep. There 
    92       ! is no need to include these in the implicit calculation. 
    93       ! 
    94       DO jk = 1, jpkm1        ! Matrix 
    95          DO jj = 2, jpjm1  
    96             DO ji = fs_2, fs_jpim1   ! vector opt. 
     98      !! bottom boundary condition: both zwi and zws must be masked as avmu can take 
     99      !! non zero value at the ocean bottom depending on the bottom friction 
     100      !! used but the bottom velocities have already been updated with the bottom 
     101      !! friction velocity in dyn_bfr using values from the previous timestep. There 
     102      !! is no need to include these in the implicit calculation. 
     103 
     104      ! The code has been modified here to implicitly implement bottom 
     105      ! friction: u(v)mask is not necessary here anymore.  
     106      ! H. Liu, April 2010. 
     107 
     108      ! 1. Vertical diffusion on u 
     109      DO jj = 2, jpjm1  
     110         DO ji = fs_2, fs_jpim1   ! vector opt. 
     111            ikbum1 = mbku(ji,jj) 
     112               zbfru = bfrua(ji,jj) 
     113 
     114            DO jk = 1, ikbum1 
    97115               zcoef = - p2dt / fse3u(ji,jj,jk) 
    98                zzwi          = zcoef * avmu (ji,jj,jk  ) / fse3uw(ji,jj,jk  ) 
    99                zwi(ji,jj,jk) = zzwi  * umask(ji,jj,jk) 
    100                zzws          = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1) 
    101                zws(ji,jj,jk) = zzws  * umask(ji,jj,jk+1) 
    102                zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws 
    103             END DO 
    104          END DO 
    105       END DO 
    106       DO jj = 2, jpjm1        ! Surface boudary conditions 
    107          DO ji = fs_2, fs_jpim1   ! vector opt. 
     116               zwi(ji,jj,jk) = zcoef * avmu(ji,jj,jk  ) / fse3uw(ji,jj,jk  ) 
     117               zws(ji,jj,jk) = zcoef * avmu(ji,jj,jk+1) / fse3uw(ji,jj,jk+1) 
     118               zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zws(ji,jj,jk) 
     119            END DO 
     120 
     121      ! Surface boundary conditions 
    108122            zwi(ji,jj,1) = 0._wp 
    109123            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) 
    110          END DO 
    111       END DO 
     124 
     125      ! Bottom boundary conditions  ! H. Liu, May, 2010 
     126!           !commented out to be consistent with v3.2, h.liu 
     127!           z2dtf = p2dt * zbfru / fse3u(ji,jj,ikbum1) * 2.0_wp * zbfr_imp 
     128            z2dtf = p2dt * zbfru / fse3u(ji,jj,ikbum1) * 1.0_wp * zbfr_imp 
     129            zws(ji,jj,ikbum1) = 0._wp 
     130            zwd(ji,jj,ikbum1) = 1._wp - zwi(ji,jj,ikbum1) - z2dtf  
    112131 
    113132      ! Matrix inversion starting from the first level 
    … …  
    125144      !   The solution (the after velocity) is in ua 
    126145      !----------------------------------------------------------------------- 
    127       ! 
    128       DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  == 
    129          DO jj = 2, jpjm1    
    130             DO ji = fs_2, fs_jpim1   ! vector opt. 
     146 
     147      ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k) 
     148            DO jk = 2, ikbum1 
    131149               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 
    132150            END DO 
    133          END DO 
    134       END DO 
    135       ! 
    136       DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  == 
    137          DO ji = fs_2, fs_jpim1   ! vector opt. 
    138             ua(ji,jj,1) = ub(ji,jj,1) + p2dt * (  ua(ji,jj,1) + 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   & 
    139                &                                                       / ( fse3u(ji,jj,1) * rau0       )  ) 
    140          END DO 
    141       END DO 
    142       DO jk = 2, jpkm1 
    143          DO jj = 2, jpjm1    
    144             DO ji = fs_2, fs_jpim1   ! vector opt. 
     151 
     152      ! second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1 
     153            z2dtf = 0.5_wp * p2dt / ( fse3u(ji,jj,1) * rau0 ) 
     154            ua(ji,jj,1) = ub(ji,jj,1) + p2dt * ua(ji,jj,1) + z2dtf * (utau_b(ji,jj) + utau(ji,jj)) 
     155           DO jk = 2, ikbum1    
    145156               zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk)   ! zrhs=right hand side 
    146157               ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1) 
    147158            END DO 
    148          END DO 
    149       END DO 
    150       ! 
    151       DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk  == 
    152          DO ji = fs_2, fs_jpim1   ! vector opt. 
    153             ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 
    154          END DO 
    155       END DO 
    156       DO jk = jpk-2, 1, -1 
    157          DO jj = 2, jpjm1    
    158             DO ji = fs_2, fs_jpim1   ! vector opt. 
    159                ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk) 
     159 
     160 
     161      ! third recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk 
     162            ua(ji,jj,ikbum1) = ua(ji,jj,ikbum1) / zwd(ji,jj,ikbum1) 
     163            DO jk = ikbum1-1, 1, -1 
     164               ua(ji,jj,jk) =( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk) 
    160165            END DO 
    161166         END DO 
    … …  
    174179      ! 2. Vertical diffusion on v 
    175180      ! --------------------------- 
    176       ! Matrix and second member construction 
    177       ! bottom boundary condition: both zwi and zws must be masked as avmv can take 
    178       ! non zero value at the ocean bottom depending on the bottom friction 
    179       ! used but the bottom velocities have already been updated with the bottom 
    180       ! friction velocity in dyn_bfr using values from the previous timestep. There 
    181       ! is no need to include these in the implicit calculation. 
    182       ! 
    183       DO jk = 1, jpkm1        ! Matrix 
     181 
     182      DO ji = fs_2, fs_jpim1   ! vector opt. 
    184183         DO jj = 2, jpjm1    
    185             DO ji = fs_2, fs_jpim1   ! vector opt. 
     184            ikbvm1 = mbkv(ji,jj) 
     185               zbfrv = bfrva(ji,jj) 
     186 
     187            DO jk = 1, ikbvm1 
    186188               zcoef = -p2dt / fse3v(ji,jj,jk) 
    187                zzwi          = zcoef * avmv (ji,jj,jk  ) / fse3vw(ji,jj,jk  ) 
    188                zwi(ji,jj,jk) =  zzwi * vmask(ji,jj,jk) 
    189                zzws          = zcoef * avmv (ji,jj,jk+1) / fse3vw(ji,jj,jk+1) 
    190                zws(ji,jj,jk) =  zzws * vmask(ji,jj,jk+1) 
    191                zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws 
    192             END DO 
    193          END DO 
    194       END DO 
    195       DO jj = 2, jpjm1        ! Surface boudary conditions 
    196          DO ji = fs_2, fs_jpim1   ! vector opt. 
     189               zwi(ji,jj,jk) = zcoef * avmv(ji,jj,jk  ) / fse3vw(ji,jj,jk  ) 
     190               zws(ji,jj,jk) = zcoef * avmv(ji,jj,jk+1) / fse3vw(ji,jj,jk+1) 
     191               zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zws(ji,jj,jk) 
     192            END DO 
     193 
     194      ! Surface boundary conditions 
    197195            zwi(ji,jj,1) = 0._wp 
    198196            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) 
    199          END DO 
    200       END DO 
     197 
     198      ! Bottom boundary conditions  ! H. Liu, May, 2010 
     199!           !commented out to be consistent with v3.2, h.liu 
     200!           z2dtf = p2dt * zbfrv / fse3v(ji,jj,ikbvm1) * 2.0_wp * zbfr_imp 
     201            z2dtf = p2dt * zbfrv / fse3v(ji,jj,ikbvm1) * 1.0_wp * zbfr_imp 
     202            zws(ji,jj,ikbvm1) = 0._wp 
     203            zwd(ji,jj,ikbvm1) = 1._wp - zwi(ji,jj,ikbvm1) - z2dtf  
    201204 
    202205      ! Matrix inversion 
    … …  
    210213      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk ) 
    211214      ! 
    212       !   m is decomposed in the product of an upper and lower triangular matrix 
     215      !   m is decomposed in the product of an upper and lower triangular 
     216      !   matrix 
    213217      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi 
    214218      !   The solution (after velocity) is in 2d array va 
    215219      !----------------------------------------------------------------------- 
    216       ! 
    217       DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  == 
    218          DO jj = 2, jpjm1    
    219             DO ji = fs_2, fs_jpim1   ! vector opt. 
     220 
     221      ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k) 
     222            DO jk = 2, ikbvm1 
    220223               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 
    221224            END DO 
    222          END DO 
    223       END DO 
    224       ! 
    225       DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  == 
    226          DO ji = fs_2, fs_jpim1   ! vector opt. 
    227             va(ji,jj,1) = vb(ji,jj,1) + p2dt * (  va(ji,jj,1) + 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   & 
    228                &                                                       / ( fse3v(ji,jj,1) * rau0       )  ) 
    229          END DO 
    230       END DO 
    231       DO jk = 2, jpkm1 
    232          DO jj = 2, jpjm1 
    233             DO ji = fs_2, fs_jpim1   ! vector opt. 
     225 
     226      ! second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1 
     227            z2dtf = 0.5_wp * p2dt / ( fse3v(ji,jj,1)*rau0 ) 
     228            va(ji,jj,1) = vb(ji,jj,1) + p2dt * va(ji,jj,1) + z2dtf * (vtau_b(ji,jj) + vtau(ji,jj)) 
     229            DO jk = 2, ikbvm1 
    234230               zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk)   ! zrhs=right hand side 
    235231               va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1) 
    236232            END DO 
    237          END DO 
    238       END DO 
    239       ! 
    240       DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  == 
    241          DO ji = fs_2, fs_jpim1   ! vector opt. 
    242             va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 
    243          END DO 
    244       END DO 
    245       DO jk = jpk-2, 1, -1 
    246          DO jj = 2, jpjm1    
    247             DO ji = fs_2, fs_jpim1   ! vector opt. 
    248                va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk) 
    249             END DO 
     233 
     234      ! third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk 
     235            va(ji,jj,ikbvm1) = va(ji,jj,ikbvm1) / zwd(ji,jj,ikbvm1) 
     236 
     237            DO jk = ikbvm1-1, 1, -1 
     238               va(ji,jj,jk) =( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk) 
     239            END DO 
     240 
    250241         END DO 
    251242      END DO 
    … …  
    262253      IF( wrk_not_released(3, 3) )   CALL ctl_stop('dyn_zdf_imp: failed to release workspace array') 
    263254      ! 
     255 
    264256   END SUBROUTINE dyn_zdf_imp 
    265257 

  • branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90

    r2977 r3116  
    183183#if defined key_bdy 
    184184      ssha(:,:) = ssha(:,:) * bdytmask(:,:) 
    185       CALL lbc_lnk( ssha, 'T', 1. )  
     185      CALL lbc_lnk( ssha, 'T', 1. )                 ! absolutly compulsory !! (jmm) 
    186186#endif 
    187187 
Note: See TracChangeset for help on using the changeset viewer.