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 3294 for trunk/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Ignore:
Timestamp:
2012-01-28T17:44:18+01:00 (12 years ago)
Author:
rblod
Message:

Merge of 3.4beta into the trunk

Location:
trunk/NEMOGCM/NEMO/OPA_SRC/DYN
Files:
24 edited
1 copied

Legend:

Unmodified
Added
Removed

  • trunk/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90

    r2715 r3294  
    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 
    3231   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    3332   USE lib_mpp         ! MPP library 
     33   USE wrk_nemo        ! Memory Allocation 
     34   USE timing          ! Timing 
    3435 
    3536   IMPLICIT NONE 
    … …  
    8485      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
    8586      ! 
    86       REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zwu   ! specific 2D workspace 
    87       REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zwv   ! specific 2D workspace 
    88       ! 
    8987      INTEGER ::   ji, jj, jk, jl           ! dummy loop indices 
    9088      INTEGER ::   ii, ij, ijt, iju, ierr   ! local integer 
    9189      REAL(wp) ::  zraur, zdep              ! local scalar 
    92       !!---------------------------------------------------------------------- 
    93  
     90      REAL(wp), POINTER,  DIMENSION(:,:) ::   zwu   ! specific 2D workspace 
     91      REAL(wp), POINTER,  DIMENSION(:,:) ::   zwv   ! specific 2D workspace 
     92      !!---------------------------------------------------------------------- 
     93      ! 
     94      IF( nn_timing == 1 )  CALL timing_start('div_cur') 
     95      ! 
     96      CALL wrk_alloc( jpi  , jpj+2, zwu               ) 
     97      CALL wrk_alloc( jpi+4, jpj  , zwv, kjstart = -1 ) 
     98      ! 
    9499      IF( kt == nit000 ) THEN 
    95100         IF(lwp) WRITE(numout,*) 
    96101         IF(lwp) WRITE(numout,*) 'div_cur : horizontal velocity divergence and relative vorticity' 
    97102         IF(lwp) WRITE(numout,*) '~~~~~~~   NOT optimal for auto-tasking case' 
    98          ! 
    99          ALLOCATE( zwu( jpi, 1:jpj+2) , zwv(-1:jpi+2, jpj) , STAT=ierr ) 
    100          IF( lk_mpp    )   CALL mpp_sum( ierr ) 
    101          IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'div_cur : unable to allocate arrays' ) 
    102103      ENDIF 
    103104 
    … …  
    121122         END DO 
    122123 
    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          
    133124         IF( .NOT. AGRIF_Root() ) THEN 
    134125            IF ((nbondi ==  1).OR.(nbondi == 2)) hdivn(nlci-1 , :     ,jk) = 0.e0      ! east 
    … …  
    241232      CALL lbc_lnk( hdivn, 'T', 1. )   ;   CALL lbc_lnk( rotn , 'F', 1. )    ! lateral boundary cond. (no sign change) 
    242233      ! 
     234      CALL wrk_dealloc( jpi  , jpj+2, zwu               ) 
     235      CALL wrk_dealloc( jpi+4, jpj  , zwv, kjstart = -1 ) 
     236      ! 
     237      IF( nn_timing == 1 )  CALL timing_stop('div_cur') 
     238      ! 
    243239   END SUBROUTINE div_cur 
    244240    
    … …  
    278274      REAL(wp) ::   zraur, zdep   ! local scalars 
    279275      !!---------------------------------------------------------------------- 
    280  
     276      ! 
     277      IF( nn_timing == 1 )  CALL timing_start('div_cur') 
     278      ! 
    281279      IF( kt == nit000 ) THEN 
    282280         IF(lwp) WRITE(numout,*) 
    … …  
    304302         END DO   
    305303 
    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          
    316304         IF( .NOT. AGRIF_Root() ) THEN 
    317305            IF ((nbondi ==  1).OR.(nbondi == 2)) hdivn(nlci-1 , :     ,jk) = 0.e0      ! east 
    … …  
    340328      CALL lbc_lnk( hdivn, 'T', 1. )   ;   CALL lbc_lnk( rotn , 'F', 1. )     ! lateral boundary cond. (no sign change) 
    341329      ! 
     330      IF( nn_timing == 1 )  CALL timing_stop('div_cur') 
     331      ! 
    342332   END SUBROUTINE div_cur 
    343333    

  • trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv.F90

    r2715 r3294  
    1919   USE in_out_manager  ! I/O manager 
    2020   USE lib_mpp         ! MPP library 
     21   USE timing          ! Timing 
    2122 
    2223   IMPLICIT NONE 
    … …  
    5758      !!---------------------------------------------------------------------- 
    5859      ! 
     60      IF( nn_timing == 1 )  CALL timing_start('dyn_adv') 
     61      ! 
    5962      SELECT CASE ( nadv )                  ! compute advection trend and add it to general trend 
    6063      CASE ( 0 )      
    … …  
    7275                      CALL dyn_adv_ubs ( kt ) 
    7376      END SELECT 
     77      ! 
     78      IF( nn_timing == 1 )  CALL timing_stop('dyn_adv') 
    7479      ! 
    7580   END SUBROUTINE dyn_adv 

  • trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv_cen2.F90

    r2715 r3294  
    2020   USE lib_mpp        ! MPP library 
    2121   USE prtctl         ! Print control 
     22   USE wrk_nemo        ! Memory Allocation 
     23   USE timing          ! Timing 
    2224 
    2325   IMPLICIT NONE 
    … …  
    4749      !! ** Action  :   (ua,va) updated with the now vorticity term trend 
    4850      !!---------------------------------------------------------------------- 
    49       USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released 
    50       USE oce     , ONLY:   zfu   => ta       , zfv   => sa       ! (ta,sa) used as 3D workspace 
    51       USE wrk_nemo, ONLY:   zfu_t => wrk_3d_1 , zfv_t => wrk_3d_4 , zfu_uw =>wrk_3d_6   ! 3D workspaces 
    52       USE wrk_nemo, ONLY:   zfu_f => wrk_3d_2 , zfv_f => wrk_3d_5 , zfv_vw =>wrk_3d_7 
    53       USE wrk_nemo, ONLY:   zfw   => wrk_3d_3  
    54       ! 
    5551      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index 
    5652      ! 
    5753      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    5854      REAL(wp) ::   zbu, zbv     ! local scalars 
     55      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zfu_t, zfv_t, zfu_f, zfv_f, zfu_uw, zfv_vw, zfw 
     56      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zfu, zfv 
    5957      !!---------------------------------------------------------------------- 
    60  
     58      ! 
     59      IF( nn_timing == 1 )  CALL timing_start('dyn_adv_cen2') 
     60      ! 
     61      CALL wrk_alloc( jpi, jpj, jpk, zfu_t, zfv_t, zfu_f, zfv_f, zfu_uw, zfv_vw, zfu, zfv, zfw ) 
     62      ! 
    6163      IF( kt == nit000 .AND. lwp ) THEN 
    6264         WRITE(numout,*) 
    … …  
    6466         WRITE(numout,*) '~~~~~~~~~~~~' 
    6567      ENDIF 
    66  
    67       ! Check that global workspace arrays aren't already in use 
    68       IF( wrk_in_use(3, 1,2,3,4,5,6,7) ) THEN 
    69          CALL ctl_stop('dyn_adv_cen2 : requested workspace array unavailable')   ;   RETURN 
    70       ENDIF 
    71  
     68      ! 
    7269      IF( l_trddyn ) THEN           ! Save ua and va trends 
    7370         zfu_uw(:,:,:) = ua(:,:,:) 
    … …  
    162159         &                       tab3d_2=va, clinfo2=           ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    163160      ! 
    164       IF( wrk_not_released(3, 1,2,3,4,5,6,7) )   CALL ctl_stop('dyn_adv_cen2: failed to release workspace array') 
     161      CALL wrk_dealloc( jpi, jpj, jpk, zfu_t, zfv_t, zfu_f, zfv_f, zfu_uw, zfv_vw, zfu, zfv, zfw ) 
     162      ! 
     163      IF( nn_timing == 1 )  CALL timing_stop('dyn_adv_cen2') 
    165164      ! 
    166165   END SUBROUTINE dyn_adv_cen2 

  • trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv_ubs.F90

    r2715 r3294  
    2222   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
    2323   USE lib_mpp        ! MPP library 
     24   USE wrk_nemo        ! Memory Allocation 
     25   USE timing          ! Timing 
    2426 
    2527   IMPLICIT NONE 
    2628   PRIVATE 
    2729 
    28    REAL(wp), PARAMETER :: gamma1 = 1._wp/4._wp  ! =1/4 quick      ; =1/3  3rd order UBS 
     30   REAL(wp), PARAMETER :: gamma1 = 1._wp/3._wp  ! =1/4 quick      ; =1/3  3rd order UBS 
    2931   REAL(wp), PARAMETER :: gamma2 = 1._wp/8._wp  ! =0   2nd order  ; =1/8  4th order centred 
    3032 
    … …  
    6264      !!      before velocity (forward in time).  
    6365      !!      Default value (hard coded in the begining of the module) are  
    64       !!      gamma1=1/4 and gamma2=1/8. 
     66      !!      gamma1=1/3 and gamma2=1/8. 
    6567      !! 
    6668      !! ** Action : - (ua,va) updated with the 3D advective momentum trends 
    … …  
    6870      !! Reference : Shchepetkin & McWilliams, 2005, Ocean Modelling.  
    6971      !!---------------------------------------------------------------------- 
    70       USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released 
    71       USE oce     , ONLY:   zfu    => ta       , zfv    => sa      ! (ta,sa) used as 3D workspace 
    72       USE wrk_nemo, ONLY:   zfu_t  => wrk_3d_1 , zfv_t  =>wrk_3d_4 , zfu_uw =>wrk_3d_6   ! 3D workspace 
    73       USE wrk_nemo, ONLY:   zfu_f  => wrk_3d_2 , zfv_f  =>wrk_3d_5 , zfv_vw =>wrk_3d_7 
    74       USE wrk_nemo, ONLY:   zfw    => wrk_3d_3 
    75       USE wrk_nemo, ONLY:   zlu_uu => wrk_4d_1 , zlv_vv=>wrk_4d_3   ! 4D workspace 
    76       USE wrk_nemo, ONLY:   zlu_uv => wrk_4d_2 , zlv_vu=>wrk_4d_4 
    77       ! 
    7872      INTEGER, INTENT(in) ::   kt     ! ocean time-step index 
    7973      ! 
    … …  
    8175      REAL(wp) ::   zbu, zbv    ! temporary scalars 
    8276      REAL(wp) ::   zui, zvj, zfuj, zfvi, zl_u, zl_v   ! temporary scalars 
     77      REAL(wp), POINTER, DIMENSION(:,:,:  ) ::  zfu, zfv 
     78      REAL(wp), POINTER, DIMENSION(:,:,:  ) ::  zfu_t, zfv_t, zfu_f, zfv_f, zfu_uw, zfv_vw, zfw 
     79      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::  zlu_uu, zlv_vv, zlu_uv, zlv_vu 
    8380      !!---------------------------------------------------------------------- 
    84  
     81      ! 
     82      IF( nn_timing == 1 )  CALL timing_start('dyn_adv_ubs') 
     83      ! 
     84      CALL wrk_alloc( jpi, jpj, jpk,       zfu_t , zfv_t , zfu_f , zfv_f, zfu_uw, zfv_vw, zfu, zfv, zfw ) 
     85      CALL wrk_alloc( jpi, jpj, jpk, jpts, zlu_uu, zlv_vv, zlu_uv, zlv_vu                               ) 
     86      ! 
    8587      IF( kt == nit000 ) THEN 
    8688         IF(lwp) WRITE(numout,*) 
    … …  
    8890         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 
    8991      ENDIF 
    90  
    91       ! Check that required workspace arrays are not already in use 
    92       IF( wrk_in_use(3, 1,2,3,4,5,6,7) .OR. wrk_in_use(4, 1,2,3,4) ) THEN 
    93          CALL ctl_stop('dyn_adv_ubs: requested workspace array unavailable')   ;   RETURN 
    94       ENDIF 
    95  
     92      ! 
    9693      zfu_t(:,:,:) = 0._wp 
    9794      zfv_t(:,:,:) = 0._wp 
    … …  
    254251         &                       tab3d_2=va, clinfo2=           ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    255252      ! 
    256       IF( wrk_not_released(3, 1,2,3,4,5,6,7) .OR.   & 
    257           wrk_not_released(4, 1,2,3,4)       )   CALL ctl_stop('dyn_adv_ubs: failed to release workspace array') 
     253      CALL wrk_dealloc( jpi, jpj, jpk,       zfu_t , zfv_t , zfu_f , zfv_f, zfu_uw, zfv_vw, zfu, zfv, zfw ) 
     254      CALL wrk_dealloc( jpi, jpj, jpk, jpts, zlu_uu, zlv_vv, zlu_uv, zlv_vu                               ) 
     255      ! 
     256      IF( nn_timing == 1 )  CALL timing_stop('dyn_adv_ubs') 
    258257      ! 
    259258   END SUBROUTINE dyn_adv_ubs 

  • trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynbfr.F90

    r2715 r3294  
    55   !!============================================================================== 
    66   !! History :  3.2  ! 2008-11  (A. C. Coward)  Original code 
     7   !!            3.4  ! 2011-09  (H. Liu) Make it consistent with semi-implicit 
     8   !!                            Bottom friction (ln_bfrimp = .true.)  
    79   !!---------------------------------------------------------------------- 
    810 
    … …  
    1315   USE dom_oce         ! ocean space and time domain variables  
    1416   USE zdf_oce         ! ocean vertical physics variables 
     17   USE zdfbfr          ! ocean bottom friction variables 
    1518   USE trdmod          ! ocean active dynamics and tracers trends  
    1619   USE trdmod_oce      ! ocean variables trends 
    1720   USE in_out_manager  ! I/O manager 
    1821   USE prtctl          ! Print control 
     22   USE timing          ! Timing 
     23   USE wrk_nemo        ! Memory Allocation 
    1924 
    2025   IMPLICIT NONE 
    … …  
    4247      !! ** Action  :   (ua,va)   momentum trend increased by bottom friction trend 
    4348      !!--------------------------------------------------------------------- 
    44       USE oce, ONLY:   ztrduv => tsa   ! tsa used as 4D workspace 
    45       !! 
    4649      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
    4750      !!  
    … …  
    4952      INTEGER  ::   ikbu, ikbv   ! local integers 
    5053      REAL(wp) ::   zm1_2dt      ! local scalar 
     54      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdu, ztrdv 
    5155      !!--------------------------------------------------------------------- 
    5256      ! 
    53       zm1_2dt = - 1._wp / ( 2._wp * rdt ) 
     57      IF( nn_timing == 1 )  CALL timing_start('dyn_bfr') 
     58      ! 
     59      IF( .NOT.ln_bfrimp) THEN     ! only for explicit bottom friction form 
     60                                    ! implicit bfr is implemented in dynzdf_imp 
    5461 
    55       IF( l_trddyn )   THEN                      ! temporary save of ua and va trends 
    56          ztrduv(:,:,:,1) = ua(:,:,:) 
    57          ztrduv(:,:,:,2) = va(:,:,:) 
    58       ENDIF 
     62        zm1_2dt = - 1._wp / ( 2._wp * rdt ) 
     63 
     64        IF( l_trddyn )   THEN                      ! temporary save of ua and va trends 
     65           CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv ) 
     66           ztrdu(:,:,:) = ua(:,:,:) 
     67           ztrdv(:,:,:) = va(:,:,:) 
     68        ENDIF 
     69 
    5970 
    6071# if defined key_vectopt_loop 
    61       DO jj = 1, 1 
    62          DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling) 
     72        DO jj = 1, 1 
     73           DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling) 
    6374# else 
    64       DO jj = 2, jpjm1 
    65          DO ji = 2, jpim1 
     75        DO jj = 2, jpjm1 
     76           DO ji = 2, jpim1 
    6677# 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 
     78              ikbu = mbku(ji,jj)          ! deepest ocean u- & v-levels 
     79              ikbv = mbkv(ji,jj) 
     80              ! 
     81              ! Apply stability criteria on absolute value  : abs(bfr/e3) < 1/(2dt) => bfr/e3 > -1/(2dt) 
     82              ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + MAX(  bfrua(ji,jj) / fse3u(ji,jj,ikbu) , zm1_2dt  ) * ub(ji,jj,ikbu) 
     83              va(ji,jj,ikbv) = va(ji,jj,ikbv) + MAX(  bfrva(ji,jj) / fse3v(ji,jj,ikbv) , zm1_2dt  ) * vb(ji,jj,ikbv) 
     84           END DO 
     85        END DO 
    7586 
     87        ! 
     88        IF( l_trddyn )   THEN                      ! save the vertical diffusive trends for further diagnostics 
     89           ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
     90           ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
     91           CALL trd_mod( ztrdu(:,:,:), ztrdv(:,:,:), jpdyn_trd_bfr, 'DYN', kt ) 
     92           CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv ) 
     93        ENDIF 
     94        !                                          ! print mean trends (used for debugging) 
     95        IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' bfr  - Ua: ', mask1=umask,               & 
     96           &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
     97        ! 
     98      ENDIF     ! end explicit bottom friction 
    7699      ! 
    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' ) 
     100      IF( nn_timing == 1 )  CALL timing_stop('dyn_bfr') 
    85101      ! 
    86102   END SUBROUTINE dyn_bfr 

  • trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90

    r2715 r3294  
    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  (H. Liu) hpg_prj: Original code for s-coordinates 
     17   !!                 !           (A. Coward) 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 
    … …  
    3737   USE lbclnk          ! lateral boundary condition  
    3838   USE lib_mpp         ! MPP library 
     39   USE wrk_nemo        ! Memory Allocation 
     40   USE timing          ! Timing 
    3941 
    4042   IMPLICIT NONE 
    … …  
    4850   LOGICAL , PUBLIC ::   ln_hpg_zps    = .FALSE.   !: z-coordinate - partial steps (interpolation) 
    4951   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) 
    5252   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 
     53   LOGICAL , PUBLIC ::   ln_hpg_prj    = .FALSE.   !: s-coordinate (Pressure Jacobian scheme) 
    5554   LOGICAL , PUBLIC ::   ln_dynhpg_imp = .FALSE.   !: semi-implicite hpg flag 
    5655 
    57    INTEGER  ::   nhpg  =  0   ! = 0 to 6, type of pressure gradient scheme used ! (deduced from ln_hpg_... flags) 
     56   INTEGER  ::   nhpg  =  0   ! = 0 to 7, type of pressure gradient scheme used ! (deduced from ln_hpg_... flags) 
    5857 
    5958   !! * Substitutions 
    … …  
    7776      !!             - Save the trend (l_trddyn=T) 
    7877      !!---------------------------------------------------------------------- 
    79       USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released 
    80       USE wrk_nemo, ONLY:   ztrdu => wrk_3d_1 , ztrdv => wrk_3d_2   ! 3D workspace 
    81       !! 
    8278      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
    83       !!---------------------------------------------------------------------- 
    84       ! 
    85       IF( wrk_in_use(3, 1,2) ) THEN 
    86          CALL ctl_stop('dyn_hpg: requested workspace arrays are unavailable')   ;   RETURN 
    87       ENDIF 
     79      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdu, ztrdv 
     80      !!---------------------------------------------------------------------- 
     81      ! 
     82      IF( nn_timing == 1 )  CALL timing_start('dyn_hpg') 
    8883      ! 
    8984      IF( l_trddyn ) THEN                    ! Temporary saving of ua and va trends (l_trddyn) 
     85         CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv ) 
    9086         ztrdu(:,:,:) = ua(:,:,:)   
    9187         ztrdv(:,:,:) = va(:,:,:)  
    9288      ENDIF       
    9389      ! 
    94       SELECT CASE ( nhpg )      ! Hydrastatic pressure gradient computation 
     90      SELECT CASE ( nhpg )      ! Hydrostatic pressure gradient computation 
    9591      CASE (  0 )   ;   CALL hpg_zco    ( kt )      ! z-coordinate 
    9692      CASE (  1 )   ;   CALL hpg_zps    ( kt )      ! z-coordinate plus partial steps (interpolation) 
    9793      CASE (  2 )   ;   CALL hpg_sco    ( kt )      ! s-coordinate (standard jacobian formulation) 
    98       CASE (  3 )   ;   CALL hpg_hel    ( kt )      ! s-coordinate (helsinki modification) 
    99       CASE (  4 )   ;   CALL hpg_wdj    ( kt )      ! s-coordinate (weighted density jacobian) 
    100       CASE (  5 )   ;   CALL hpg_djc    ( kt )      ! s-coordinate (Density Jacobian with Cubic polynomial) 
    101       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) 
    10296      END SELECT 
    10397      ! 
    … …  
    106100         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    107101         CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_hpg, 'DYN', kt ) 
     102         CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv ) 
    108103      ENDIF           
    109104      ! 
    … …  
    111106         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    112107      ! 
    113       IF( wrk_not_released(3, 1,2) )   CALL ctl_stop('dyn_hpg: failed to release workspace arrays') 
     108      IF( nn_timing == 1 )  CALL timing_stop('dyn_hpg') 
    114109      ! 
    115110   END SUBROUTINE dyn_hpg 
    … …  
    128123      INTEGER ::   ioptio = 0      ! temporary integer 
    129124      !! 
    130       NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco, ln_hpg_hel,    & 
    131          &                 ln_hpg_wdj, ln_hpg_djc, ln_hpg_rot, rn_gamma  , ln_dynhpg_imp 
     125      NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco,     & 
     126         &                 ln_hpg_djc, ln_hpg_prj, ln_dynhpg_imp 
    132127      !!---------------------------------------------------------------------- 
    133128      ! 
    … …  
    143138         WRITE(numout,*) '      z-coord. - partial steps (interpolation)          ln_hpg_zps    = ', ln_hpg_zps 
    144139         WRITE(numout,*) '      s-coord. (standard jacobian formulation)          ln_hpg_sco    = ', ln_hpg_sco 
    145          WRITE(numout,*) '      s-coord. (helsinki modification)                  ln_hpg_hel    = ', ln_hpg_hel 
    146          WRITE(numout,*) '      s-coord. (weighted density jacobian)              ln_hpg_wdj    = ', ln_hpg_wdj 
    147140         WRITE(numout,*) '      s-coord. (Density Jacobian: Cubic polynomial)     ln_hpg_djc    = ', ln_hpg_djc 
    148          WRITE(numout,*) '      s-coord. (ROTated axes scheme)                    ln_hpg_rot    = ', ln_hpg_rot 
    149          WRITE(numout,*) '      weighting coeff. (wdj scheme)                     rn_gamma      = ', rn_gamma 
     141         WRITE(numout,*) '      s-coord. (Pressure Jacobian: Cubic polynomial)    ln_hpg_prj    = ', ln_hpg_prj 
    150142         WRITE(numout,*) '      time stepping: centered (F) or semi-implicit (T)  ln_dynhpg_imp = ', ln_dynhpg_imp 
    151143      ENDIF 
    152144      ! 
    153       IF( lk_vvl .AND. .NOT. ln_hpg_sco )   & 
    154          &   CALL ctl_stop('dyn_hpg_init : variable volume key_vvl require the standard jacobian formulation hpg_sco') 
     145      IF( ln_hpg_djc )   & 
     146         &   CALL ctl_stop('dyn_hpg_init : Density Jacobian: Cubic polynominal method & 
     147                           & currently disabled (bugs under investigation). Please select & 
     148                           & either  ln_hpg_sco or  ln_hpg_prj instead') 
     149      ! 
     150      IF( lk_vvl .AND. .NOT. (ln_hpg_sco.OR.ln_hpg_prj) )   & 
     151         &   CALL ctl_stop('dyn_hpg_init : variable volume key_vvl requires:& 
     152                           & the standard jacobian formulation hpg_sco or & 
     153                           & the pressure jacobian formulation hpg_prj') 
    155154      ! 
    156155      !                               ! Set nhpg from ln_hpg_... flags 
    … …  
    158157      IF( ln_hpg_zps )   nhpg = 1 
    159158      IF( ln_hpg_sco )   nhpg = 2 
    160       IF( ln_hpg_hel )   nhpg = 3 
    161       IF( ln_hpg_wdj )   nhpg = 4 
    162       IF( ln_hpg_djc )   nhpg = 5 
    163       IF( ln_hpg_rot )   nhpg = 6 
    164       ! 
    165       !                               ! Consitency check 
     159      IF( ln_hpg_djc )   nhpg = 3 
     160      IF( ln_hpg_prj )   nhpg = 4 
     161      ! 
     162      !                               ! Consistency check 
    166163      ioptio = 0  
    167164      IF( ln_hpg_zco )   ioptio = ioptio + 1 
    168165      IF( ln_hpg_zps )   ioptio = ioptio + 1 
    169166      IF( ln_hpg_sco )   ioptio = ioptio + 1 
    170       IF( ln_hpg_hel )   ioptio = ioptio + 1 
    171       IF( ln_hpg_wdj )   ioptio = ioptio + 1 
    172167      IF( ln_hpg_djc )   ioptio = ioptio + 1 
    173       IF( ln_hpg_rot )   ioptio = ioptio + 1 
     168      IF( ln_hpg_prj )   ioptio = ioptio + 1 
    174169      IF( ioptio /= 1 )   CALL ctl_stop( 'NO or several hydrostatic pressure gradient options used' ) 
    175170      ! 
    … …  
    193188      !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 
    194189      !!---------------------------------------------------------------------- 
    195       USE oce, ONLY:   zhpi => ta , zhpj => sa   ! (ta,sa) used as 3D workspace 
    196       !! 
    197190      INTEGER, INTENT(in) ::   kt    ! ocean time-step index 
    198191      !! 
    199192      INTEGER  ::   ji, jj, jk       ! dummy loop indices 
    200193      REAL(wp) ::   zcoef0, zcoef1   ! temporary scalars 
    201       !!---------------------------------------------------------------------- 
    202        
     194      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zhpi, zhpj  
     195      !!---------------------------------------------------------------------- 
     196      !   
     197      CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj ) 
     198      ! 
    203199      IF( kt == nit000 ) THEN 
    204200         IF(lwp) WRITE(numout,*) 
    … …  
    221217         END DO 
    222218      END DO 
     219 
    223220      ! 
    224221      ! interior value (2=<jk=<jpkm1) 
    … …  
    242239      END DO 
    243240      ! 
     241      CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj ) 
     242      ! 
    244243   END SUBROUTINE hpg_zco 
    245244 
    … …  
    253252      !! ** Action  : - Update (ua,va) with the now hydrastatic pressure trend 
    254253      !!----------------------------------------------------------------------  
    255       USE oce, ONLY:   zhpi => ta , zhpj => sa   ! (ta,sa) used as 3D workspace 
    256       !! 
    257254      INTEGER, INTENT(in) ::   kt    ! ocean time-step index 
    258255      !! 
    … …  
    260257      INTEGER  ::   iku, ikv                         ! temporary integers 
    261258      REAL(wp) ::   zcoef0, zcoef1, zcoef2, zcoef3   ! temporary scalars 
    262       !!---------------------------------------------------------------------- 
    263  
     259      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zhpi, zhpj  
     260      !!---------------------------------------------------------------------- 
     261      ! 
     262      CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj ) 
     263      ! 
    264264      IF( kt == nit000 ) THEN 
    265265         IF(lwp) WRITE(numout,*) 
    … …  
    267267         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   z-coordinate with partial steps - vector optimization' 
    268268      ENDIF 
     269 
    269270 
    270271      ! Local constant initialization 
    … …  
    284285      END DO 
    285286 
     287 
    286288      ! interior value (2=<jk=<jpkm1) 
    287289      DO jk = 2, jpkm1 
    … …  
    303305         END DO 
    304306      END DO 
     307 
    305308 
    306309      ! partial steps correction at the last level  (use gru & grv computed in zpshde.F90) 
    … …  
    333336      END DO 
    334337      ! 
     338      CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj ) 
     339      ! 
    335340   END SUBROUTINE hpg_zps 
    336341 
    … …  
    354359      !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 
    355360      !!---------------------------------------------------------------------- 
    356       USE oce, ONLY:   zhpi => ta , zhpj => sa   ! (ta,sa) used as 3D workspace 
    357       !! 
    358361      INTEGER, INTENT(in) ::   kt    ! ocean time-step index 
    359362      !! 
    360363      INTEGER  ::   ji, jj, jk                 ! dummy loop indices 
    361364      REAL(wp) ::   zcoef0, zuap, zvap, znad   ! temporary scalars 
    362       !!---------------------------------------------------------------------- 
    363  
     365      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zhpi, zhpj  
     366      !!---------------------------------------------------------------------- 
     367      ! 
     368      CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj ) 
     369      ! 
    364370      IF( kt == nit000 ) THEN 
    365371         IF(lwp) WRITE(numout,*) 
    … …  
    417423      END DO 
    418424      ! 
     425      CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj ) 
     426      ! 
    419427   END SUBROUTINE hpg_sco 
    420  
    421  
    422    SUBROUTINE hpg_hel( kt ) 
    423       !!--------------------------------------------------------------------- 
    424       !!                  ***  ROUTINE hpg_hel  *** 
    425       !! 
    426       !! ** Method  :   s-coordinate case. 
    427       !!      The now hydrostatic pressure gradient at a given level 
    428       !!      jk is computed by taking the vertical integral of the in-situ  
    429       !!      density gradient along the model level from the suface to that  
    430       !!      level. s-coordinates (ln_sco): a corrective term is added 
    431       !!      to the horizontal pressure gradient : 
    432       !!         zhpi = grav .....  + 1/e1u mi(rhd) di[ grav dep3w ] 
    433       !!         zhpj = grav .....  + 1/e2v mj(rhd) dj[ grav dep3w ] 
    434       !!      add it to the general momentum trend (ua,va). 
    435       !!         ua = ua - 1/e1u * zhpi 
    436       !!         va = va - 1/e2v * zhpj 
    437       !! 
    438       !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 
    439       !!             - Save the trend (l_trddyn=T) 
    440       !!---------------------------------------------------------------------- 
    441       USE oce, ONLY:   zhpi => ta , zhpj => sa   ! (ta,sa) used as 3D workspace 
    442       !! 
    443       INTEGER, INTENT(in) ::   kt    ! ocean time-step index 
    444       !! 
    445       INTEGER  ::   ji, jj, jk           ! dummy loop indices 
    446       REAL(wp) ::   zcoef0, zuap, zvap   ! temporary scalars 
    447       !!---------------------------------------------------------------------- 
    448  
    449       IF( kt == nit000 ) THEN 
    450          IF(lwp) WRITE(numout,*) 
    451          IF(lwp) WRITE(numout,*) 'dyn:hpg_hel : hydrostatic pressure gradient trend' 
    452          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate case, helsinki modified scheme' 
    453       ENDIF 
    454  
    455       ! Local constant initialization 
    456       zcoef0 = - grav * 0.5_wp 
    457   
    458       ! Surface value 
    459       DO jj = 2, jpjm1 
    460          DO ji = fs_2, fs_jpim1   ! vector opt. 
    461             ! hydrostatic pressure gradient along s-surfaces 
    462             zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( fse3t(ji+1,jj  ,1) * rhd(ji+1,jj  ,1)  & 
    463                &                                  - fse3t(ji  ,jj  ,1) * rhd(ji  ,jj  ,1) ) 
    464             zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( fse3t(ji  ,jj+1,1) * rhd(ji  ,jj+1,1)  & 
    465                &                                  - fse3t(ji  ,jj  ,1) * rhd(ji  ,jj  ,1) ) 
    466             ! s-coordinate pressure gradient correction 
    467             zuap = -zcoef0 * ( rhd   (ji+1,jj,1) + rhd   (ji,jj,1) )   & 
    468                &           * ( fsdept(ji+1,jj,1) - fsdept(ji,jj,1) ) / e1u(ji,jj) 
    469             zvap = -zcoef0 * ( rhd   (ji,jj+1,1) + rhd   (ji,jj,1) )   & 
    470                &           * ( fsdept(ji,jj+1,1) - fsdept(ji,jj,1) ) / e2v(ji,jj) 
    471             ! add to the general momentum trend 
    472             ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap 
    473             va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap 
    474          END DO 
    475       END DO 
    476       ! 
    477       ! interior value (2=<jk=<jpkm1) 
    478       DO jk = 2, jpkm1 
    479          DO jj = 2, jpjm1 
    480             DO ji = fs_2, fs_jpim1   ! vector opt. 
    481                ! hydrostatic pressure gradient along s-surfaces 
    482                zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) & 
    483                   &           +  zcoef0 / e1u(ji,jj) * ( fse3t(ji+1,jj,jk  ) * rhd(ji+1,jj,jk)     & 
    484                   &                                     -fse3t(ji  ,jj,jk  ) * rhd(ji  ,jj,jk)   ) & 
    485                   &           +  zcoef0 / e1u(ji,jj) * ( fse3t(ji+1,jj,jk-1) * rhd(ji+1,jj,jk-1)   & 
    486                   &                                     -fse3t(ji  ,jj,jk-1) * rhd(ji  ,jj,jk-1) ) 
    487                zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) & 
    488                   &           +  zcoef0 / e2v(ji,jj) * ( fse3t(ji,jj+1,jk  ) * rhd(ji,jj+1,jk)     & 
    489                   &                                     -fse3t(ji,jj  ,jk  ) * rhd(ji,jj,  jk)   ) & 
    490                   &           +  zcoef0 / e2v(ji,jj) * ( fse3t(ji,jj+1,jk-1) * rhd(ji,jj+1,jk-1)   & 
    491                   &                                     -fse3t(ji,jj  ,jk-1) * rhd(ji,jj,  jk-1) ) 
    492                ! s-coordinate pressure gradient correction 
    493                zuap = - zcoef0 * ( rhd   (ji+1,jj,jk) + rhd   (ji,jj,jk) )   & 
    494                   &            * ( fsdept(ji+1,jj,jk) - fsdept(ji,jj,jk) ) / e1u(ji,jj) 
    495                zvap = - zcoef0 * ( rhd   (ji,jj+1,jk) + rhd   (ji,jj,jk) )   & 
    496                   &            * ( fsdept(ji,jj+1,jk) - fsdept(ji,jj,jk) ) / e2v(ji,jj) 
    497                ! add to the general momentum trend 
    498                ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap 
    499                va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) + zvap 
    500             END DO 
    501          END DO 
    502       END DO 
    503       ! 
    504    END SUBROUTINE hpg_hel 
    505  
    506  
    507    SUBROUTINE hpg_wdj( kt ) 
    508       !!--------------------------------------------------------------------- 
    509       !!                  ***  ROUTINE hpg_wdj  *** 
    510       !! 
    511       !! ** Method  :   Weighted Density Jacobian (wdj) scheme (song 1998) 
    512       !!      The weighting coefficients from the namelist parameter rn_gamma 
    513       !!      (alpha=0.5-rn_gamma ; beta=1-alpha=0.5+rn_gamma 
    514       !! 
    515       !! Reference : Song, Mon. Wea. Rev., 126, 3213-3230, 1998. 
    516       !!---------------------------------------------------------------------- 
    517       USE oce, ONLY:   zhpi => ta , zhpj => sa   ! (ta,sa) used as 3D workspace 
    518       !! 
    519       INTEGER, INTENT(in) ::   kt    ! ocean time-step index 
    520       !! 
    521       INTEGER  ::   ji, jj, jk           ! dummy loop indices 
    522       REAL(wp) ::   zcoef0, zuap, zvap   ! temporary scalars 
    523       REAL(wp) ::   zalph , zbeta        !    "         " 
    524       !!---------------------------------------------------------------------- 
    525  
    526       IF( kt == nit000 ) THEN 
    527          IF(lwp) WRITE(numout,*) 
    528          IF(lwp) WRITE(numout,*) 'dyn:hpg_wdj : hydrostatic pressure gradient trend' 
    529          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   Weighted Density Jacobian' 
    530       ENDIF 
    531  
    532       ! Local constant initialization 
    533       zcoef0 = - grav * 0.5_wp 
    534       zalph  = 0.5_wp - rn_gamma    ! weighting coefficients (alpha=0.5-rn_gamma 
    535       zbeta  = 0.5_wp + rn_gamma    !                        (beta =1-alpha=0.5+rn_gamma 
    536  
    537       ! Surface value (no ponderation) 
    538       DO jj = 2, jpjm1 
    539          DO ji = fs_2, fs_jpim1   ! vector opt. 
    540             ! hydrostatic pressure gradient along s-surfaces 
    541             zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * (  fse3w(ji+1,jj  ,1) * rhd(ji+1,jj  ,1)   & 
    542                &                                   - fse3w(ji  ,jj  ,1) * rhd(ji  ,jj  ,1)  ) 
    543             zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * (  fse3w(ji  ,jj+1,1) * rhd(ji  ,jj+1,1)   & 
    544                &                                   - fse3w(ji  ,jj  ,1) * rhd(ji,  jj  ,1)  ) 
    545             ! s-coordinate pressure gradient correction 
    546             zuap = -zcoef0 * ( rhd   (ji+1,jj,1) + rhd   (ji,jj,1) )   & 
    547                &           * ( fsde3w(ji+1,jj,1) - fsde3w(ji,jj,1) ) / e1u(ji,jj) 
    548             zvap = -zcoef0 * ( rhd   (ji,jj+1,1) + rhd   (ji,jj,1) )   & 
    549                &           * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj) 
    550             ! add to the general momentum trend 
    551             ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap 
    552             va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap 
    553          END DO 
    554       END DO 
    555  
    556       ! Interior value (2=<jk=<jpkm1) (weighted with zalph & zbeta) 
    557       DO jk = 2, jpkm1 
    558          DO jj = 2, jpjm1 
    559             DO ji = fs_2, fs_jpim1   ! vector opt. 
    560                zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj)                            & 
    561                   &           * (   (            fsde3w(ji+1,jj,jk  ) + fsde3w(ji,jj,jk  )        & 
    562                   &                            - fsde3w(ji+1,jj,jk-1) - fsde3w(ji,jj,jk-1)    )   & 
    563                   &               * (  zalph * ( rhd   (ji+1,jj,jk-1) - rhd   (ji,jj,jk-1) )      & 
    564                   &                  + zbeta * ( rhd   (ji+1,jj,jk  ) - rhd   (ji,jj,jk  ) )  )   & 
    565                   &             -   (            rhd   (ji+1,jj,jk  ) + rhd   (ji,jj,jk  )        & 
    566                   &                           - rhd   (ji+1,jj,jk-1) - rhd   (ji,jj,jk-1)     )   & 
    567                   &               * (  zalph * ( fsde3w(ji+1,jj,jk-1) - fsde3w(ji,jj,jk-1) )      & 
    568                   &                  + zbeta * ( fsde3w(ji+1,jj,jk  ) - fsde3w(ji,jj,jk  ) )  )  ) 
    569                zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj)                            & 
    570                   &           * (   (           fsde3w(ji,jj+1,jk  ) + fsde3w(ji,jj,jk  )         & 
    571                   &                           - fsde3w(ji,jj+1,jk-1) - fsde3w(ji,jj,jk-1)     )   & 
    572                   &               * (  zalph * ( rhd   (ji,jj+1,jk-1) - rhd   (ji,jj,jk-1) )      & 
    573                   &                  + zbeta * ( rhd   (ji,jj+1,jk  ) - rhd   (ji,jj,jk  ) )  )   & 
    574                   &             -   (            rhd   (ji,jj+1,jk  ) + rhd   (ji,jj,jk  )        & 
    575                   &                            - rhd   (ji,jj+1,jk-1) - rhd   (ji,jj,jk-1)    )   & 
    576                   &               * (  zalph * ( fsde3w(ji,jj+1,jk-1) - fsde3w(ji,jj,jk-1) )      & 
    577                   &                  + zbeta * ( fsde3w(ji,jj+1,jk  ) - fsde3w(ji,jj,jk  ) )  )  ) 
    578                ! add to the general momentum trend 
    579                ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) 
    580                va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) 
    581             END DO 
    582          END DO 
    583       END DO 
    584       ! 
    585    END SUBROUTINE hpg_wdj 
    586  
    587428 
    588429   SUBROUTINE hpg_djc( kt ) 
    … …  
    594435      !! Reference: Shchepetkin and McWilliams, J. Geophys. Res., 108(C3), 3090, 2003 
    595436      !!---------------------------------------------------------------------- 
    596       USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released 
    597       USE oce     , ONLY:   zhpi  => ta        , zhpj => sa       ! (ta,sa) used as 3D workspace 
    598       USE wrk_nemo, ONLY:   drhox => wrk_3d_1  , dzx  => wrk_3d_2 
    599       USE wrk_nemo, ONLY:   drhou => wrk_3d_3  , dzu  => wrk_3d_4 , rho_i => wrk_3d_5 
    600       USE wrk_nemo, ONLY:   drhoy => wrk_3d_6  , dzy  => wrk_3d_7 
    601       USE wrk_nemo, ONLY:   drhov => wrk_3d_8  , dzv  => wrk_3d_9 , rho_j => wrk_3d_10 
    602       USE wrk_nemo, ONLY:   drhoz => wrk_3d_11 , dzz  => wrk_3d_12  
    603       USE wrk_nemo, ONLY:   drhow => wrk_3d_13 , dzw  => wrk_3d_14 
    604       USE wrk_nemo, ONLY:   rho_k => wrk_3d_15 
    605       !! 
    606437      INTEGER, INTENT(in) ::   kt    ! ocean time-step index 
    607438      !! 
    … …  
    610441      REAL(wp) ::   z1_10, cffu, cffx   !    "         " 
    611442      REAL(wp) ::   z1_12, cffv, cffy   !    "         " 
    612       !!---------------------------------------------------------------------- 
    613  
    614       IF( wrk_in_use(3, 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15) ) THEN 
    615          CALL ctl_stop('dyn:hpg_djc: requested workspace arrays unavailable')   ;   RETURN 
    616       ENDIF 
     443      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zhpi, zhpj  
     444      REAL(wp), POINTER, DIMENSION(:,:,:) ::  dzx, dzy, dzz, dzu, dzv, dzw 
     445      REAL(wp), POINTER, DIMENSION(:,:,:) ::  drhox, drhoy, drhoz, drhou, drhov, drhow 
     446      REAL(wp), POINTER, DIMENSION(:,:,:) ::  rho_i, rho_j, rho_k 
     447      !!---------------------------------------------------------------------- 
     448      ! 
     449      CALL wrk_alloc( jpi, jpj, jpk, dzx  , dzy  , dzz  , dzu  , dzv  , dzw   )  
     450      CALL wrk_alloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow )  
     451      CALL wrk_alloc( jpi, jpj, jpk, rho_i, rho_j, rho_k,  zhpi,  zhpj        )  
     452      ! 
    617453 
    618454      IF( kt == nit000 ) THEN 
    … …  
    811647      END DO 
    812648      ! 
    813       IF( wrk_not_released(3, 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15) )   & 
    814          CALL ctl_stop('dyn:hpg_djc: failed to release workspace arrays') 
     649      CALL wrk_dealloc( jpi, jpj, jpk, dzx  , dzy  , dzz  , dzu  , dzv  , dzw   )  
     650      CALL wrk_dealloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow )  
     651      CALL wrk_dealloc( jpi, jpj, jpk, rho_i, rho_j, rho_k,  zhpi,  zhpj        )  
    815652      ! 
    816653   END SUBROUTINE hpg_djc 
    817654 
    818655 
    819    SUBROUTINE hpg_rot( kt ) 
     656   SUBROUTINE hpg_prj( kt ) 
    820657      !!--------------------------------------------------------------------- 
    821       !!                  ***  ROUTINE hpg_rot  *** 
    822       !! 
    823       !! ** Method  :   rotated axes scheme (Thiem and Berntsen 2005) 
    824       !! 
    825       !! Reference: Thiem & Berntsen, Ocean Modelling, In press, 2005. 
    826       !!---------------------------------------------------------------------- 
    827       USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released 
    828       USE oce     , ONLY:   zhpi    => ta       , zhpj    => sa       ! (ta,sa) used as 3D workspace 
    829       USE wrk_nemo, ONLY:   zdistr  => wrk_2d_1 , zsina   => wrk_2d_2 , zcosa  => wrk_2d_3 
    830       USE wrk_nemo, ONLY:   zhpiorg => wrk_3d_1 , zhpirot => wrk_3d_2 
    831       USE wrk_nemo, ONLY:   zhpitra => wrk_3d_3 , zhpine  => wrk_3d_4 
    832       USE wrk_nemo, ONLY:   zhpjorg => wrk_3d_5 , zhpjrot => wrk_3d_6 
    833       USE wrk_nemo, ONLY:   zhpjtra => wrk_3d_7 , zhpjne  => wrk_3d_8 
    834       !! 
    835       INTEGER, INTENT(in) ::   kt    ! ocean time-step index 
    836       !! 
    837       INTEGER  ::   ji, jj, jk          ! dummy loop indices 
    838       REAL(wp) ::   zforg, zcoef0, zuap, zmskd1, zmskd1m   ! temporary scalar 
    839       REAL(wp) ::   zfrot        , zvap, zmskd2, zmskd2m   !    "         " 
    840       !!---------------------------------------------------------------------- 
    841  
    842       IF( wrk_in_use(2, 1,2,3)             .OR.   & 
    843           wrk_in_use(3, 1,2,3,4,5,6,7,8) ) THEN 
    844          CALL ctl_stop('dyn:hpg_rot: requested workspace arrays unavailable')   ;   RETURN 
    845       ENDIF 
    846  
     658      !!                  ***  ROUTINE hpg_prj  *** 
     659      !! 
     660      !! ** Method  :   s-coordinate case. 
     661      !!      A Pressure-Jacobian horizontal pressure gradient method 
     662      !!      based on the constrained cubic-spline interpolation for 
     663      !!      all vertical coordinate systems 
     664      !! 
     665      !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 
     666      !!             - Save the trend (l_trddyn=T) 
     667      !! 
     668      !!---------------------------------------------------------------------- 
     669      INTEGER, PARAMETER  :: polynomial_type = 1    ! 1: cubic spline, 2: linear 
     670      INTEGER, INTENT(in) ::   kt                   ! ocean time-step index 
     671      !! 
     672      INTEGER  ::   ji, jj, jk, jkk                 ! dummy loop indices 
     673      REAL(wp) ::   zcoef0, znad                    ! temporary scalars 
     674      !! 
     675      !! The local variables for the correction term 
     676      INTEGER  :: jk1, jis, jid, jjs, jjd 
     677      REAL(wp) :: zuijk, zvijk, zpwes, zpwed, zpnss, zpnsd, zdeps 
     678      REAL(wp) :: zrhdt1  
     679      REAL(wp) :: zdpdx1, zdpdx2, zdpdy1, zdpdy2 
     680      INTEGER  :: zbhitwe, zbhitns 
     681      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zdeptht, zrhh  
     682      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp 
     683      !!---------------------------------------------------------------------- 
     684      ! 
     685      CALL wrk_alloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp )  
     686      CALL wrk_alloc( jpi,jpj,jpk, zdeptht, zrhh )  
     687      ! 
    847688      IF( kt == nit000 ) THEN 
    848689         IF(lwp) WRITE(numout,*) 
    849          IF(lwp) WRITE(numout,*) 'dyn:hpg_rot : hydrostatic pressure gradient trend' 
    850          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate case, rotated axes scheme used' 
     690         IF(lwp) WRITE(numout,*) 'dyn:hpg_prj : hydrostatic pressure gradient trend' 
     691         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate case, cubic spline pressure Jacobian' 
    851692      ENDIF 
    852693 
    853       ! ------------------------------- 
    854       !  Local constant initialization 
    855       ! ------------------------------- 
    856       zcoef0 = - grav * 0.5_wp 
    857       zforg  = 0.95_wp 
    858       zfrot  = 1._wp - zforg 
    859  
    860       ! inverse of the distance between 2 diagonal T-points (defined at F-point) (here zcoef0/distance) 
    861       zdistr(:,:) = zcoef0 / SQRT( e1f(:,:)*e1f(:,:) + e2f(:,:)*e1f(:,:) ) 
    862  
    863       ! sinus and cosinus of diagonal angle at F-point 
    864       zsina(:,:) = ATAN2( e2f(:,:), e1f(:,:) ) 
    865       zcosa(:,:) = COS( zsina(:,:) ) 
    866       zsina(:,:) = SIN( zsina(:,:) ) 
    867  
    868       ! --------------- 
    869       !  Surface value 
    870       ! --------------- 
    871       ! compute and add to the general trend the pressure gradients along the axes 
    872       DO jj = 2, jpjm1 
    873          DO ji = fs_2, fs_jpim1   ! vector opt. 
    874             ! hydrostatic pressure gradient along s-surfaces 
    875             zhpiorg(ji,jj,1) = zcoef0 / e1u(ji,jj) * (  fse3t(ji+1,jj,1) * rhd(ji+1,jj,1)   & 
    876                &                                      - fse3t(ji  ,jj,1) * rhd(ji  ,jj,1)  ) 
    877             zhpjorg(ji,jj,1) = zcoef0 / e2v(ji,jj) * (  fse3t(ji,jj+1,1) * rhd(ji,jj+1,1)   & 
    878                &                                      - fse3t(ji,jj  ,1) * rhd(ji,jj  ,1)  ) 
    879             ! s-coordinate pressure gradient correction 
    880             zuap = -zcoef0 * ( rhd   (ji+1,jj  ,1) + rhd   (ji,jj,1) )   & 
    881                &           * ( fsdept(ji+1,jj  ,1) - fsdept(ji,jj,1) ) / e1u(ji,jj) 
    882             zvap = -zcoef0 * ( rhd   (ji  ,jj+1,1) + rhd   (ji,jj,1) )   & 
    883                &           * ( fsdept(ji  ,jj+1,1) - fsdept(ji,jj,1) ) / e2v(ji,jj) 
    884             ! add to the general momentum trend 
    885             ua(ji,jj,1) = ua(ji,jj,1) + zforg * ( zhpiorg(ji,jj,1) + zuap ) 
    886             va(ji,jj,1) = va(ji,jj,1) + zforg * ( zhpjorg(ji,jj,1) + zvap ) 
    887          END DO 
    888       END DO 
    889  
    890       ! compute the pressure gradients in the diagonal directions 
    891       DO jj = 1, jpjm1 
    892          DO ji = 1, fs_jpim1   ! vector opt. 
    893             zmskd1 = tmask(ji+1,jj+1,1) * tmask(ji  ,jj,1)      ! mask in the 1st diagnonal 
    894             zmskd2 = tmask(ji  ,jj+1,1) * tmask(ji+1,jj,1)      ! mask in the 2nd diagnonal 
    895             ! hydrostatic pressure gradient along s-surfaces 
    896             zhpitra(ji,jj,1) = zdistr(ji,jj) * zmskd1 * (  fse3t(ji+1,jj+1,1) * rhd(ji+1,jj+1,1)   & 
    897                &                                         - fse3t(ji  ,jj  ,1) * rhd(ji  ,jj  ,1)  ) 
    898             zhpjtra(ji,jj,1) = zdistr(ji,jj) * zmskd2 * (  fse3t(ji  ,jj+1,1) * rhd(ji  ,jj+1,1)   & 
    899                &                                         - fse3t(ji+1,jj  ,1) * rhd(ji+1,jj  ,1)  ) 
    900             ! s-coordinate pressure gradient correction 
    901             zuap = -zdistr(ji,jj) * zmskd1 * ( rhd   (ji+1,jj+1,1) + rhd   (ji  ,jj,1) )   & 
    902                &                           * ( fsdept(ji+1,jj+1,1) - fsdept(ji  ,jj,1) ) 
    903             zvap = -zdistr(ji,jj) * zmskd2 * ( rhd   (ji  ,jj+1,1) + rhd   (ji+1,jj,1) )   & 
    904                &                           * ( fsdept(ji  ,jj+1,1) - fsdept(ji+1,jj,1) ) 
    905             ! back rotation 
    906             zhpine(ji,jj,1) = zcosa(ji,jj) * ( zhpitra(ji,jj,1) + zuap )   & 
    907                &            - zsina(ji,jj) * ( zhpjtra(ji,jj,1) + zvap ) 
    908             zhpjne(ji,jj,1) = zsina(ji,jj) * ( zhpitra(ji,jj,1) + zuap )   & 
    909                &            + zcosa(ji,jj) * ( zhpjtra(ji,jj,1) + zvap ) 
    910          END DO 
    911       END DO 
    912  
    913       ! interpolate and add to the general trend the diagonal gradient 
    914       DO jj = 2, jpjm1 
    915          DO ji = fs_2, fs_jpim1   ! vector opt. 
    916             ! averaging 
    917             zhpirot(ji,jj,1) = 0.5 * ( zhpine(ji,jj,1) + zhpine(ji  ,jj-1,1) ) 
    918             zhpjrot(ji,jj,1) = 0.5 * ( zhpjne(ji,jj,1) + zhpjne(ji-1,jj  ,1) ) 
    919             ! add to the general momentum trend 
    920             ua(ji,jj,1) = ua(ji,jj,1) + zfrot * zhpirot(ji,jj,1)  
    921             va(ji,jj,1) = va(ji,jj,1) + zfrot * zhpjrot(ji,jj,1)  
    922          END DO 
    923       END DO 
    924  
    925       ! ----------------- 
    926       ! 2. interior value (2=<jk=<jpkm1) 
    927       ! ----------------- 
    928       ! compute and add to the general trend the pressure gradients along the axes 
    929       DO jk = 2, jpkm1 
    930          DO jj = 2, jpjm1 
    931             DO ji = fs_2, fs_jpim1   ! vector opt. 
    932                ! hydrostatic pressure gradient along s-surfaces 
    933                zhpiorg(ji,jj,jk) = zhpiorg(ji,jj,jk-1)                                                 & 
    934                   &              +  zcoef0 / e1u(ji,jj) * (  fse3t(ji+1,jj,jk  ) * rhd(ji+1,jj,jk  )   & 
    935                   &                                        - fse3t(ji  ,jj,jk  ) * rhd(ji  ,jj,jk  )   & 
    936                   &                                        + fse3t(ji+1,jj,jk-1) * rhd(ji+1,jj,jk-1)   & 
    937                   &                                        - fse3t(ji  ,jj,jk-1) * rhd(ji  ,jj,jk-1)  ) 
    938                zhpjorg(ji,jj,jk) = zhpjorg(ji,jj,jk-1)                                                 & 
    939                   &              +  zcoef0 / e2v(ji,jj) * (  fse3t(ji,jj+1,jk  ) * rhd(ji,jj+1,jk  )   & 
    940                   &                                        - fse3t(ji,jj  ,jk  ) * rhd(ji,jj,  jk  )   & 
    941                   &                                        + fse3t(ji,jj+1,jk-1) * rhd(ji,jj+1,jk-1)   & 
    942                   &                                        - fse3t(ji,jj  ,jk-1) * rhd(ji,jj,  jk-1)  ) 
    943                ! s-coordinate pressure gradient correction 
    944                zuap = - zcoef0 * ( rhd   (ji+1,jj  ,jk) + rhd   (ji,jj,jk) )   & 
    945                   &            * ( fsdept(ji+1,jj  ,jk) - fsdept(ji,jj,jk) ) / e1u(ji,jj) 
    946                zvap = - zcoef0 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji,jj,jk) )   & 
    947                   &            * ( fsdept(ji  ,jj+1,jk) - fsdept(ji,jj,jk) ) / e2v(ji,jj) 
    948                ! add to the general momentum trend 
    949                ua(ji,jj,jk) = ua(ji,jj,jk) + zforg*( zhpiorg(ji,jj,jk) + zuap ) 
    950                va(ji,jj,jk) = va(ji,jj,jk) + zforg*( zhpjorg(ji,jj,jk) + zvap ) 
     694      !!---------------------------------------------------------------------- 
     695      ! Local constant initialization 
     696      zcoef0 = - grav  
     697      znad = 0.0_wp 
     698      IF( lk_vvl ) znad = 1._wp 
     699 
     700      ! Clean 3-D work arrays 
     701      zhpi(:,:,:) = 0._wp 
     702      zrhh(:,:,:) = rhd(:,:,:) 
     703       
     704      ! Preparing vertical density profile "zrhh(:,:,:)" for hybrid-sco coordinate 
     705      DO jj = 1, jpj 
     706        DO ji = 1, jpi    
     707          jk = mbathy(ji,jj) 
     708          IF( jk <= 0 ) THEN; zrhh(ji,jj,:) = 0._wp 
     709          ELSE IF(jk == 1) THEN; zrhh(ji,jj, jk+1:jpk) = rhd(ji,jj,jk) 
     710          ELSE IF(jk < jpkm1) THEN 
     711             DO jkk = jk+1, jpk 
     712                zrhh(ji,jj,jkk) = interp1(fsde3w(ji,jj,jkk),   fsde3w(ji,jj,jkk-1), & 
     713                                         fsde3w(ji,jj,jkk-2), rhd(ji,jj,jkk-1), rhd(ji,jj,jkk-2)) 
     714             END DO  
     715          ENDIF 
     716        END DO 
     717      END DO 
     718 
     719      ! Transfer the depth of "T(:,:,:)" to vertical coordinate "zdeptht(:,:,:)" 
     720      DO jj = 1, jpj 
     721        DO ji = 1, jpi 
     722          zdeptht(ji,jj,1) = 0.5_wp * fse3w(ji,jj,1) 
     723          zdeptht(ji,jj,1) = zdeptht(ji,jj,1) - sshn(ji,jj) * znad 
     724          DO jk = 2, jpk 
     725             zdeptht(ji,jj,jk) = zdeptht(ji,jj,jk-1) + fse3w(ji,jj,jk) 
     726          END DO 
     727        END DO 
     728      END DO 
     729 
     730      DO jk = 1, jpkm1 
     731        DO jj = 1, jpj 
     732          DO ji = 1, jpi 
     733            fsp(ji,jj,jk) = zrhh(ji,jj,jk) 
     734            xsp(ji,jj,jk) = zdeptht(ji,jj,jk) 
     735          END DO 
     736        END DO 
     737      END DO 
     738 
     739      ! Construct the vertical density profile with the  
     740      ! constrained cubic spline interpolation 
     741      ! rho(z) = asp + bsp*z + csp*z^2 + dsp*z^3 
     742      CALL cspline(fsp,xsp,asp,bsp,csp,dsp,polynomial_type)       
     743 
     744      ! Integrate the hydrostatic pressure "zhpi(:,:,:)" at "T(ji,jj,1)" 
     745      DO jj = 2, jpj 
     746        DO ji = 2, jpi  
     747          zrhdt1 = zrhh(ji,jj,1) - interp3(zdeptht(ji,jj,1),asp(ji,jj,1), & 
     748                                         bsp(ji,jj,1),   csp(ji,jj,1), & 
     749                                         dsp(ji,jj,1) ) * 0.5_wp * zdeptht(ji,jj,1) 
     750          zrhdt1 = MAX(zrhdt1, 1000._wp - rau0)        ! no lighter than fresh water 
     751 
     752          ! assuming linear profile across the top half surface layer 
     753          zhpi(ji,jj,1) =  0.5_wp * fse3w(ji,jj,1) * zrhdt1   
     754        END DO 
     755      END DO 
     756 
     757      ! Calculate the pressure "zhpi(:,:,:)" at "T(ji,jj,2:jpkm1)" 
     758      DO jk = 2, jpkm1                                   
     759        DO jj = 2, jpj      
     760          DO ji = 2, jpi 
     761            zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) +                          & 
     762                             integ2(zdeptht(ji,jj,jk-1), zdeptht(ji,jj,jk),& 
     763                                    asp(ji,jj,jk-1),    bsp(ji,jj,jk-1), & 
     764                                    csp(ji,jj,jk-1),    dsp(ji,jj,jk-1)) 
     765          END DO 
     766        END DO 
     767      END DO 
     768 
     769      ! Z coordinate of U(ji,jj,1:jpkm1) and V(ji,jj,1:jpkm1) 
     770      DO jj = 2, jpjm1      
     771        DO ji = 2, jpim1   
     772          zu(ji,jj,1) = - ( fse3u(ji,jj,1) - sshu_n(ji,jj) * znad) 
     773          zv(ji,jj,1) = - ( fse3v(ji,jj,1) - sshv_n(ji,jj) * znad) 
     774        END DO 
     775      END DO 
     776 
     777      DO jk = 2, jpkm1                                   
     778        DO jj = 2, jpjm1      
     779          DO ji = 2, jpim1   
     780            zu(ji,jj,jk) = zu(ji,jj,jk-1)- fse3u(ji,jj,jk) 
     781            zv(ji,jj,jk) = zv(ji,jj,jk-1)- fse3v(ji,jj,jk) 
     782          END DO 
     783        END DO 
     784      END DO 
     785                
     786      DO jk = 1, jpkm1                                   
     787        DO jj = 2, jpjm1      
     788          DO ji = 2, jpim1   
     789            zu(ji,jj,jk) = zu(ji,jj,jk) + 0.5_wp * fse3u(ji,jj,jk) 
     790            zv(ji,jj,jk) = zv(ji,jj,jk) + 0.5_wp * fse3v(ji,jj,jk) 
     791          END DO 
     792        END DO 
     793      END DO 
     794 
     795      DO jk = 1, jpkm1                                   
     796        DO jj = 2, jpjm1      
     797          DO ji = 2, jpim1   
     798            zpwes = 0._wp; zpwed = 0._wp 
     799            zpnss = 0._wp; zpnsd = 0._wp 
     800            zuijk = zu(ji,jj,jk) 
     801            zvijk = zv(ji,jj,jk) 
     802 
     803            !!!!!     for u equation 
     804            IF( jk <= mbku(ji,jj) ) THEN 
     805               IF( -zdeptht(ji+1,jj,mbku(ji,jj)) >= -zdeptht(ji,jj,mbku(ji,jj)) ) THEN 
     806                 jis = ji + 1; jid = ji 
     807               ELSE 
     808                 jis = ji;     jid = ji +1 
     809               ENDIF 
     810 
     811               ! integrate the pressure on the shallow side 
     812               jk1 = jk  
     813               zbhitwe = 0 
     814               DO WHILE ( -zdeptht(jis,jj,jk1) > zuijk ) 
     815                 IF( jk1 == mbku(ji,jj) ) THEN 
     816                   zbhitwe = 1 
     817                   EXIT 
     818                 ENDIF 
     819                 zdeps = MIN(zdeptht(jis,jj,jk1+1), -zuijk) 
     820                 zpwes = zpwes +                                    &  
     821                      integ2(zdeptht(jis,jj,jk1), zdeps,            & 
     822                             asp(jis,jj,jk1),    bsp(jis,jj,jk1), & 
     823                             csp(jis,jj,jk1),    dsp(jis,jj,jk1)) 
     824                 jk1 = jk1 + 1 
     825               END DO 
     826             
     827               IF(zbhitwe == 1) THEN 
     828                 zuijk = -zdeptht(jis,jj,jk1) 
     829               ENDIF 
     830 
     831               ! integrate the pressure on the deep side 
     832               jk1 = jk  
     833               zbhitwe = 0 
     834               DO WHILE ( -zdeptht(jid,jj,jk1) < zuijk ) 
     835                 IF( jk1 == 1 ) THEN 
     836                   zbhitwe = 1 
     837                   EXIT 
     838                 ENDIF 
     839                 zdeps = MAX(zdeptht(jid,jj,jk1-1), -zuijk) 
     840                 zpwed = zpwed +                                        &  
     841                        integ2(zdeps,              zdeptht(jid,jj,jk1), & 
     842                               asp(jid,jj,jk1-1), bsp(jid,jj,jk1-1),  & 
     843                               csp(jid,jj,jk1-1), dsp(jid,jj,jk1-1) ) 
     844                 jk1 = jk1 - 1 
     845               END DO 
     846             
     847               IF( zbhitwe == 1 ) THEN 
     848                 zdeps = zdeptht(jid,jj,1) + MIN(zuijk, sshn(jid,jj)*znad) 
     849                 zrhdt1 = zrhh(jid,jj,1) - interp3(zdeptht(jid,jj,1), asp(jid,jj,1), & 
     850                                                 bsp(jid,jj,1),    csp(jid,jj,1), & 
     851                                                 dsp(jid,jj,1)) * zdeps 
     852                 zrhdt1 = MAX(zrhdt1, 1000._wp - rau0)        ! no lighter than fresh water 
     853                 zpwed  = zpwed + 0.5_wp * (zrhh(jid,jj,1) + zrhdt1) * zdeps 
     854               ENDIF 
     855 
     856               ! update the momentum trends in u direction 
     857 
     858               zdpdx1 = zcoef0 / e1u(ji,jj) * (zhpi(ji+1,jj,jk) - zhpi(ji,jj,jk)) 
     859               IF( lk_vvl ) THEN 
     860                 zdpdx2 = zcoef0 / e1u(ji,jj) * &  
     861                         ( REAL(jis-jid, wp) * (zpwes + zpwed) + (sshn(ji+1,jj)-sshn(ji,jj)) )  
     862                ELSE 
     863                 zdpdx2 = zcoef0 / e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed)  
     864               ENDIF 
     865 
     866               ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * & 
     867               &           umask(ji,jj,jk) * tmask(ji,jj,jk) * tmask(ji+1,jj,jk) 
     868            ENDIF 
     869   
     870            !!!!!     for v equation 
     871            IF( jk <= mbkv(ji,jj) ) THEN 
     872               IF( -zdeptht(ji,jj+1,mbkv(ji,jj)) >= -zdeptht(ji,jj,mbkv(ji,jj)) ) THEN 
     873                 jjs = jj + 1; jjd = jj 
     874               ELSE 
     875                 jjs = jj    ; jjd = jj + 1 
     876               ENDIF 
     877 
     878               ! integrate the pressure on the shallow side 
     879               jk1 = jk  
     880               zbhitns = 0 
     881               DO WHILE ( -zdeptht(ji,jjs,jk1) > zvijk ) 
     882                 IF( jk1 == mbkv(ji,jj) ) THEN 
     883                   zbhitns = 1 
     884                   EXIT 
     885                 ENDIF 
     886                 zdeps = MIN(zdeptht(ji,jjs,jk1+1), -zvijk) 
     887                 zpnss = zpnss +                                      &  
     888                        integ2(zdeptht(ji,jjs,jk1), zdeps,            & 
     889                               asp(ji,jjs,jk1),    bsp(ji,jjs,jk1), & 
     890                               csp(ji,jjs,jk1),    dsp(ji,jjs,jk1) ) 
     891                 jk1 = jk1 + 1 
     892               END DO 
     893             
     894               IF(zbhitns == 1) THEN 
     895                 zvijk = -zdeptht(ji,jjs,jk1) 
     896               ENDIF 
     897 
     898               ! integrate the pressure on the deep side 
     899               jk1 = jk  
     900               zbhitns = 0 
     901               DO WHILE ( -zdeptht(ji,jjd,jk1) < zvijk ) 
     902                 IF( jk1 == 1 ) THEN 
     903                   zbhitns = 1 
     904                   EXIT 
     905                 ENDIF 
     906                 zdeps = MAX(zdeptht(ji,jjd,jk1-1), -zvijk) 
     907                 zpnsd = zpnsd +                                        &  
     908                        integ2(zdeps,              zdeptht(ji,jjd,jk1), & 
     909                               asp(ji,jjd,jk1-1), bsp(ji,jjd,jk1-1), & 
     910                               csp(ji,jjd,jk1-1), dsp(ji,jjd,jk1-1) ) 
     911                 jk1 = jk1 - 1 
     912               END DO 
     913             
     914               IF( zbhitns == 1 ) THEN 
     915                 zdeps = zdeptht(ji,jjd,1) + MIN(zvijk, sshn(ji,jjd)*znad) 
     916                 zrhdt1 = zrhh(ji,jjd,1) - interp3(zdeptht(ji,jjd,1), asp(ji,jjd,1), & 
     917                                                 bsp(ji,jjd,1),    csp(ji,jjd,1), & 
     918                                                 dsp(ji,jjd,1) ) * zdeps 
     919                 zrhdt1 = MAX(zrhdt1, 1000._wp - rau0)        ! no lighter than fresh water 
     920                 zpnsd  = zpnsd + 0.5_wp * (zrhh(ji,jjd,1) + zrhdt1) * zdeps 
     921               ENDIF 
     922 
     923               ! update the momentum trends in v direction 
     924 
     925               zdpdy1 = zcoef0 / e2v(ji,jj) * (zhpi(ji,jj+1,jk) - zhpi(ji,jj,jk)) 
     926               IF( lk_vvl ) THEN 
     927                   zdpdy2 = zcoef0 / e2v(ji,jj) * & 
     928                           ( REAL(jjs-jjd, wp) * (zpnss + zpnsd) + (sshn(ji,jj+1)-sshn(ji,jj)) )  
     929               ELSE 
     930                   zdpdy2 = zcoef0 / e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd )  
     931               ENDIF 
     932 
     933               va(ji,jj,jk) = va(ji,jj,jk) + (zdpdy1 + zdpdy2)*& 
     934               &              vmask(ji,jj,jk)*tmask(ji,jj,jk)*tmask(ji,jj+1,jk) 
     935            ENDIF 
     936 
     937                     
     938           END DO 
     939        END DO 
     940      END DO 
     941      ! 
     942      CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp )  
     943      CALL wrk_dealloc( jpi,jpj,jpk, zdeptht, zrhh )  
     944      ! 
     945   END SUBROUTINE hpg_prj 
     946 
     947   SUBROUTINE cspline(fsp, xsp, asp, bsp, csp, dsp, polynomial_type) 
     948      !!---------------------------------------------------------------------- 
     949      !!                 ***  ROUTINE cspline  *** 
     950      !!        
     951      !! ** Purpose :   constrained cubic spline interpolation 
     952      !!           
     953      !! ** Method  :   f(x) = asp + bsp*x + csp*x^2 + dsp*x^3  
     954      !! Reference: CJC Kruger, Constrained Cubic Spline Interpoltation 
     955      !! 
     956      !!---------------------------------------------------------------------- 
     957      IMPLICIT NONE 
     958      REAL(wp), DIMENSION(:,:,:), INTENT(in)  :: fsp, xsp           ! value and coordinate 
     959      REAL(wp), DIMENSION(:,:,:), INTENT(out) :: asp, bsp, csp, dsp ! coefficients of  
     960                                                                    ! the interpoated function 
     961      INTEGER, INTENT(in) :: polynomial_type                        ! 1: cubic spline  
     962                                                                    ! 2: Linear 
     963 
     964      ! Local Variables       
     965      INTEGER  ::   ji, jj, jk                 ! dummy loop indices 
     966      INTEGER  ::   jpi, jpj, jpkm1 
     967      REAL(wp) ::   zdf1, zdf2, zddf1, zddf2, ztmp1, ztmp2, zdxtmp 
     968      REAL(wp) ::   zdxtmp1, zdxtmp2, zalpha 
     969      REAL(wp) ::   zdf(size(fsp,3)) 
     970      !!---------------------------------------------------------------------- 
     971 
     972      jpi   = size(fsp,1) 
     973      jpj   = size(fsp,2) 
     974      jpkm1 = size(fsp,3) - 1 
     975 
     976       
     977      IF (polynomial_type == 1) THEN     ! Constrained Cubic Spline 
     978         DO ji = 1, jpi 
     979            DO jj = 1, jpj 
     980           !!Fritsch&Butland's method, 1984 (preferred, but more computation)               
     981           !    DO jk = 2, jpkm1-1 
     982           !       zdxtmp1 = xsp(ji,jj,jk)   - xsp(ji,jj,jk-1)   
     983           !       zdxtmp2 = xsp(ji,jj,jk+1) - xsp(ji,jj,jk)   
     984           !       zdf1    = ( fsp(ji,jj,jk)   - fsp(ji,jj,jk-1) ) / zdxtmp1 
     985           !       zdf2    = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk)   ) / zdxtmp2 
     986           ! 
     987           !       zalpha = ( zdxtmp1 + 2._wp * zdxtmp2 ) / ( zdxtmp1 + zdxtmp2 ) / 3._wp 
     988           !      
     989           !       IF(zdf1 * zdf2 <= 0._wp) THEN 
     990           !           zdf(jk) = 0._wp 
     991           !       ELSE 
     992           !         zdf(jk) = zdf1 * zdf2 / ( ( 1._wp - zalpha ) * zdf1 + zalpha * zdf2 ) 
     993           !       ENDIF 
     994           !    END DO 
     995            
     996           !!Simply geometric average 
     997               DO jk = 2, jpkm1-1 
     998                  zdf1 = (fsp(ji,jj,jk) - fsp(ji,jj,jk-1)) / (xsp(ji,jj,jk) - xsp(ji,jj,jk-1)) 
     999                  zdf2 = (fsp(ji,jj,jk+1) - fsp(ji,jj,jk)) / (xsp(ji,jj,jk+1) - xsp(ji,jj,jk)) 
     1000             
     1001                  IF(zdf1 * zdf2 <= 0._wp) THEN 
     1002                     zdf(jk) = 0._wp 
     1003                  ELSE 
     1004                     zdf(jk) = 2._wp * zdf1 * zdf2 / (zdf1 + zdf2) 
     1005                  ENDIF 
     1006               END DO 
     1007            
     1008               zdf(1)     = 1.5_wp * ( fsp(ji,jj,2) - fsp(ji,jj,1) ) / & 
     1009                          &          ( xsp(ji,jj,2) - xsp(ji,jj,1) ) -  0.5_wp * zdf(2) 
     1010               zdf(jpkm1) = 1.5_wp * ( fsp(ji,jj,jpkm1) - fsp(ji,jj,jpkm1-1) ) / & 
     1011                          &          ( xsp(ji,jj,jpkm1) - xsp(ji,jj,jpkm1-1) ) - & 
     1012                          & 0.5_wp * zdf(jpkm1 - 1) 
     1013    
     1014               DO jk = 1, jpkm1 - 1 
     1015                 zdxtmp = xsp(ji,jj,jk+1) - xsp(ji,jj,jk)  
     1016                 ztmp1  = (zdf(jk+1) + 2._wp * zdf(jk)) / zdxtmp 
     1017                 ztmp2  =  6._wp * (fsp(ji,jj,jk+1) - fsp(ji,jj,jk)) / zdxtmp / zdxtmp 
     1018                 zddf1  = -2._wp * ztmp1 + ztmp2  
     1019                 ztmp1  = (2._wp * zdf(jk+1) + zdf(jk)) / zdxtmp 
     1020                 zddf2  =  2._wp * ztmp1 - ztmp2  
     1021       
     1022                 dsp(ji,jj,jk) = (zddf2 - zddf1) / 6._wp / zdxtmp 
     1023                 csp(ji,jj,jk) = ( xsp(ji,jj,jk+1) * zddf1 - xsp(ji,jj,jk)*zddf2 ) / 2._wp / zdxtmp 
     1024                 bsp(ji,jj,jk) = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk) ) / zdxtmp - &  
     1025                               & csp(ji,jj,jk) * ( xsp(ji,jj,jk+1) + xsp(ji,jj,jk) ) - & 
     1026                               & dsp(ji,jj,jk) * ((xsp(ji,jj,jk+1) + xsp(ji,jj,jk))**2 - & 
     1027                               &                   xsp(ji,jj,jk+1) * xsp(ji,jj,jk)) 
     1028                 asp(ji,jj,jk) = fsp(ji,jj,jk) - xsp(ji,jj,jk) * (bsp(ji,jj,jk) + & 
     1029                               &                (xsp(ji,jj,jk) * (csp(ji,jj,jk) + & 
     1030                               &                 dsp(ji,jj,jk) * xsp(ji,jj,jk)))) 
     1031               END DO 
    9511032            END DO 
    9521033         END DO 
    953       END DO 
    954  
    955       ! compute the pressure gradients in the diagonal directions 
    956       DO jk = 2, jpkm1 
    957          DO jj = 1, jpjm1 
    958             DO ji = 1, fs_jpim1   ! vector opt. 
    959                zmskd1  = tmask(ji+1,jj+1,jk  ) * tmask(ji  ,jj,jk  )      ! level jk   mask in the 1st diagnonal 
    960                zmskd1m = tmask(ji+1,jj+1,jk-1) * tmask(ji  ,jj,jk-1)      ! level jk-1    "               "      
    961                zmskd2  = tmask(ji  ,jj+1,jk  ) * tmask(ji+1,jj,jk  )      ! level jk   mask in the 2nd diagnonal 
    962                zmskd2m = tmask(ji  ,jj+1,jk-1) * tmask(ji+1,jj,jk-1)      ! level jk-1    "               "      
    963                ! hydrostatic pressure gradient along s-surfaces 
    964                zhpitra(ji,jj,jk) = zhpitra(ji,jj,jk-1)                                                       & 
    965                   &              + zdistr(ji,jj) * zmskd1  * ( fse3t(ji+1,jj+1,jk  ) * rhd(ji+1,jj+1,jk)     & 
    966                   &                                           -fse3t(ji  ,jj  ,jk  ) * rhd(ji  ,jj  ,jk) )   & 
    967                   &              + zdistr(ji,jj) * zmskd1m * ( fse3t(ji+1,jj+1,jk-1) * rhd(ji+1,jj+1,jk-1)   & 
    968                   &                                           -fse3t(ji  ,jj  ,jk-1) * rhd(ji  ,jj  ,jk-1) ) 
    969                zhpjtra(ji,jj,jk) = zhpjtra(ji,jj,jk-1)                                                       & 
    970                   &              + zdistr(ji,jj) * zmskd2  * ( fse3t(ji  ,jj+1,jk  ) * rhd(ji  ,jj+1,jk)     & 
    971                   &                                           -fse3t(ji+1,jj  ,jk  ) * rhd(ji+1,jj,  jk) )   & 
    972                   &              + zdistr(ji,jj) * zmskd2m * ( fse3t(ji  ,jj+1,jk-1) * rhd(ji  ,jj+1,jk-1)   & 
    973                   &                                           -fse3t(ji+1,jj  ,jk-1) * rhd(ji+1,jj,  jk-1) ) 
    974                ! s-coordinate pressure gradient correction 
    975                zuap = - zdistr(ji,jj) * zmskd1 * ( rhd   (ji+1,jj+1,jk) + rhd   (ji  ,jj,jk) )   & 
    976                   &                            * ( fsdept(ji+1,jj+1,jk) - fsdept(ji  ,jj,jk) ) 
    977                zvap = - zdistr(ji,jj) * zmskd2 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji+1,jj,jk) )   & 
    978                   &                            * ( fsdept(ji  ,jj+1,jk) - fsdept(ji+1,jj,jk) ) 
    979                ! back rotation 
    980                zhpine(ji,jj,jk) = zcosa(ji,jj) * ( zhpitra(ji,jj,jk) + zuap )   & 
    981                   &             - zsina(ji,jj) * ( zhpjtra(ji,jj,jk) + zvap ) 
    982                zhpjne(ji,jj,jk) = zsina(ji,jj) * ( zhpitra(ji,jj,jk) + zuap )   & 
    983                   &             + zcosa(ji,jj) * ( zhpjtra(ji,jj,jk) + zvap ) 
     1034  
     1035      ELSE IF (polynomial_type == 2) THEN     ! Linear 
     1036         DO ji = 1, jpi 
     1037            DO jj = 1, jpj 
     1038               DO jk = 1, jpkm1-1 
     1039                  zdxtmp =xsp(ji,jj,jk+1) - xsp(ji,jj,jk)  
     1040                  ztmp1 = fsp(ji,jj,jk+1) - fsp(ji,jj,jk) 
     1041    
     1042                  dsp(ji,jj,jk) = 0._wp 
     1043                  csp(ji,jj,jk) = 0._wp 
     1044                  bsp(ji,jj,jk) = ztmp1 / zdxtmp 
     1045                  asp(ji,jj,jk) = fsp(ji,jj,jk) - bsp(ji,jj,jk) * xsp(ji,jj,jk) 
     1046               END DO 
    9841047            END DO 
    9851048         END DO 
    986       END DO 
    987  
    988       ! interpolate and add to the general trend 
    989       DO jk = 2, jpkm1 
    990          DO jj = 2, jpjm1 
    991             DO ji = fs_2, fs_jpim1   ! vector opt. 
    992                ! averaging 
    993                zhpirot(ji,jj,jk) = 0.5 * ( zhpine(ji,jj,jk) + zhpine(ji  ,jj-1,jk) ) 
    994                zhpjrot(ji,jj,jk) = 0.5 * ( zhpjne(ji,jj,jk) + zhpjne(ji-1,jj  ,jk) ) 
    995                ! add to the general momentum trend 
    996                ua(ji,jj,jk) = ua(ji,jj,jk) + zfrot * zhpirot(ji,jj,jk)  
    997                va(ji,jj,jk) = va(ji,jj,jk) + zfrot * zhpjrot(ji,jj,jk)  
    998             END DO 
    999          END DO 
    1000       END DO 
    1001       ! 
    1002       IF( wrk_not_released(2, 1,2,3)           .OR.   & 
    1003           wrk_not_released(3, 1,2,3,4,5,6,7,8) )   CALL ctl_stop('dyn:hpg_rot: failed to release workspace arrays') 
    1004       ! 
    1005    END SUBROUTINE hpg_rot 
     1049 
     1050      ELSE 
     1051           CALL ctl_stop( 'invalid polynomial type in cspline' ) 
     1052      ENDIF 
     1053 
     1054       
     1055   END SUBROUTINE cspline 
     1056 
     1057 
     1058   FUNCTION interp1(x, xl, xr, fl, fr)  RESULT(f)  
     1059      !!---------------------------------------------------------------------- 
     1060      !!                 ***  ROUTINE interp1  *** 
     1061      !!        
     1062      !! ** Purpose :   1-d linear interpolation 
     1063      !!           
     1064      !! ** Method  :   
     1065      !!                interpolation is straight forward 
     1066      !!                extrapolation is also permitted (no value limit)  
     1067      !! 
     1068      !!---------------------------------------------------------------------- 
     1069      IMPLICIT NONE 
     1070      REAL(wp), INTENT(in) ::  x, xl, xr, fl, fr    
     1071      REAL(wp)             ::  f ! result of the interpolation (extrapolation) 
     1072      REAL(wp)             ::  zdeltx 
     1073      !!---------------------------------------------------------------------- 
     1074 
     1075      zdeltx = xr - xl 
     1076      IF(abs(zdeltx) <= 10._wp * EPSILON(x)) THEN 
     1077        f = 0.5_wp * (fl + fr) 
     1078      ELSE 
     1079        f = ( (x - xl ) * fr - ( x - xr ) * fl ) / zdeltx 
     1080      ENDIF 
     1081       
     1082   END FUNCTION interp1 
     1083 
     1084   FUNCTION interp2(x, a, b, c, d)  RESULT(f)  
     1085      !!---------------------------------------------------------------------- 
     1086      !!                 ***  ROUTINE interp1  *** 
     1087      !!        
     1088      !! ** Purpose :   1-d constrained cubic spline interpolation 
     1089      !!           
     1090      !! ** Method  :  cubic spline interpolation 
     1091      !! 
     1092      !!---------------------------------------------------------------------- 
     1093      IMPLICIT NONE 
     1094      REAL(wp), INTENT(in) ::  x, a, b, c, d    
     1095      REAL(wp)             ::  f ! value from the interpolation 
     1096      !!---------------------------------------------------------------------- 
     1097 
     1098      f = a + x* ( b + x * ( c + d * x ) )  
     1099 
     1100   END FUNCTION interp2 
     1101 
     1102 
     1103   FUNCTION interp3(x, a, b, c, d)  RESULT(f)  
     1104      !!---------------------------------------------------------------------- 
     1105      !!                 ***  ROUTINE interp1  *** 
     1106      !!        
     1107      !! ** Purpose :   Calculate the first order of deriavtive of 
     1108      !!                a cubic spline function y=a+b*x+c*x^2+d*x^3 
     1109      !!           
     1110      !! ** Method  :   f=dy/dx=b+2*c*x+3*d*x^2 
     1111      !! 
     1112      !!---------------------------------------------------------------------- 
     1113      IMPLICIT NONE 
     1114      REAL(wp), INTENT(in) ::  x, a, b, c, d    
     1115      REAL(wp)             ::  f ! value from the interpolation 
     1116      !!---------------------------------------------------------------------- 
     1117 
     1118      f = b + x * ( 2._wp * c + 3._wp * d * x) 
     1119 
     1120   END FUNCTION interp3 
     1121 
     1122    
     1123   FUNCTION integ2(xl, xr, a, b, c, d)  RESULT(f)  
     1124      !!---------------------------------------------------------------------- 
     1125      !!                 ***  ROUTINE interp1  *** 
     1126      !!        
     1127      !! ** Purpose :   1-d constrained cubic spline integration 
     1128      !!           
     1129      !! ** Method  :  integrate polynomial a+bx+cx^2+dx^3 from xl to xr  
     1130      !! 
     1131      !!---------------------------------------------------------------------- 
     1132      IMPLICIT NONE 
     1133      REAL(wp), INTENT(in) ::  xl, xr, a, b, c, d    
     1134      REAL(wp)             ::  za1, za2, za3       
     1135      REAL(wp)             ::  f                   ! integration result 
     1136      !!---------------------------------------------------------------------- 
     1137 
     1138      za1 = 0.5_wp * b  
     1139      za2 = c / 3.0_wp  
     1140      za3 = 0.25_wp * d  
     1141 
     1142      f  = xr * ( a + xr * ( za1 + xr * ( za2 + za3 * xr ) ) ) - & 
     1143         & xl * ( a + xl * ( za1 + xl * ( za2 + za3 * xl ) ) ) 
     1144 
     1145   END FUNCTION integ2 
     1146 
    10061147 
    10071148   !!====================================================================== 

  • trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90

    r2777 r3294  
    1919   USE lib_mpp         ! MPP library 
    2020   USE prtctl          ! Print control 
     21   USE wrk_nemo        ! Memory Allocation 
     22   USE timing          ! Timing 
    2123 
    2224   IMPLICIT NONE 
    … …  
    5254      !!             - save this trends (l_trddyn=T) for post-processing 
    5355      !!---------------------------------------------------------------------- 
    54       USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released 
    55       USE oce     , ONLY:   ztrdu => ta       , ztrdv => sa   ! (ta,sa) used as 3D workspace    
    56       USE wrk_nemo, ONLY:   zhke  => wrk_3d_1                 ! 3D workspace 
    57       !! 
    5856      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index 
    5957      !! 
    6058      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    6159      REAL(wp) ::   zu, zv       ! temporary scalars 
     60      REAL(wp), POINTER, DIMENSION(:,:,:) :: zhke 
     61      REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdu, ztrdv  
    6262      !!---------------------------------------------------------------------- 
    63  
    64       IF( wrk_in_use(3,1) ) THEN 
    65          CALL ctl_stop('dyn_keg: requested workspace array is unavailable')   ;   RETURN 
    66       ENDIF 
    67  
     63      ! 
     64      IF( nn_timing == 1 )  CALL timing_start('dyn_keg') 
     65      ! 
     66      CALL wrk_alloc( jpi, jpj, jpk, zhke ) 
     67      ! 
    6868      IF( kt == nit000 ) THEN 
    6969         IF(lwp) WRITE(numout,*) 
    … …  
    7373 
    7474      IF( l_trddyn ) THEN           ! Save ua and va trends 
     75         CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv ) 
    7576         ztrdu(:,:,:) = ua(:,:,:)  
    7677         ztrdv(:,:,:) = va(:,:,:)  
    … …  
    131132         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    132133         CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_keg, 'DYN', kt ) 
     134         CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv ) 
    133135      ENDIF 
    134136      ! 
    … …  
    136138         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    137139      ! 
    138       IF( wrk_not_released(3, 1) )   CALL ctl_stop('dyn_keg: failed to release workspace array') 
     140      CALL wrk_dealloc( jpi, jpj, jpk, zhke ) 
     141      ! 
     142      IF( nn_timing == 1 )  CALL timing_stop('dyn_keg') 
    139143      ! 
    140144   END SUBROUTINE dyn_keg 

  • trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf.F90

    r2715 r3294  
    2626   USE lib_mpp        ! distribued memory computing library 
    2727   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
     28   USE wrk_nemo        ! Memory Allocation 
     29   USE timing          ! Timing 
     30 
    2831 
    2932   IMPLICIT NONE 
    … …  
    5154      !! ** Purpose :   compute the lateral ocean dynamics physics. 
    5255      !!---------------------------------------------------------------------- 
    53       USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released 
    54       USE wrk_nemo, ONLY:   ztrdu => wrk_3d_1 , ztrdv => wrk_3d_2 
    5556      ! 
    5657      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
    57       !!---------------------------------------------------------------------- 
    58  
    59       IF( wrk_in_use(3, 1,2) ) THEN 
    60          CALL ctl_stop('dyn_ldf: requested workspace arrays unavailable')   ;   RETURN 
    61       ENDIF 
     58      ! 
     59      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdu, ztrdv 
     60      !!---------------------------------------------------------------------- 
     61      ! 
     62      IF( nn_timing == 1 )  CALL timing_start('dyn_ldf') 
    6263      ! 
    6364      IF( l_trddyn )   THEN                      ! temporary save of ta and sa trends 
     65         CALL wrk_alloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
    6466         ztrdu(:,:,:) = ua(:,:,:)  
    6567         ztrdv(:,:,:) = va(:,:,:)  
    … …  
    105107         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    106108         CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_ldf, 'DYN', kt ) 
     109         CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
    107110      ENDIF 
    108111      !                                          ! print sum trends (used for debugging) 
    … …  
    110113         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    111114      ! 
    112       IF( wrk_not_released(3, 1,2) )   CALL ctl_stop('dyn_ldf: failed to release workspace arrays') 
     115      IF( nn_timing == 1 )  CALL timing_stop('dyn_ldf') 
    113116      ! 
    114117   END SUBROUTINE dyn_ldf 
    … …  
    147150      IF( ln_dynldf_hor   )   ioptio = ioptio + 1 
    148151      IF( ln_dynldf_iso   )   ioptio = ioptio + 1 
    149       IF( ioptio /= 1 ) CALL ctl_stop( '          use only ONE direction (level/hor/iso)' ) 
     152      IF( ioptio > 1 ) CALL ctl_stop( '          use only ONE direction (level/hor/iso)' ) 
    150153 
    151154      !                                   ! Set nldf, the type of lateral diffusion, from ln_dynldf_... logicals 

  • trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilap.F90

    r2715 r3294  
    2323   USE trdmod_oce      ! ocean variables trends 
    2424   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
     25   USE wrk_nemo        ! Memory Allocation 
     26   USE timing          ! Timing 
    2527 
    2628   IMPLICIT NONE 
    … …  
    7476      !!               mixing trend. 
    7577      !!---------------------------------------------------------------------- 
    76       USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released 
    77       USE wrk_nemo, ONLY:   zcu => wrk_2d_1 , zcv => wrk_2d_2   ! 3D workspace 
    78       USE wrk_nemo, ONLY:   zuf => wrk_3d_3 , zut => wrk_3d_4   ! 3D workspace 
    79       USE wrk_nemo, ONLY:   zlu => wrk_3d_5 , zlv => wrk_3d_6 
    8078      ! 
    8179      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
    … …  
    8381      INTEGER  ::   ji, jj, jk                  ! dummy loop indices 
    8482      REAL(wp) ::   zua, zva, zbt, ze2u, ze2v   ! temporary scalar 
     83      REAL(wp), POINTER, DIMENSION(:,:  ) :: zcu, zcv 
     84      REAL(wp), POINTER, DIMENSION(:,:,:) :: zuf, zut, zlu, zlv 
    8585      !!---------------------------------------------------------------------- 
    86  
    87       IF( wrk_in_use(2, 1,2) .OR. wrk_in_use(3, 3,4,5,6) ) THEN 
    88          CALL ctl_stop('dyn_ldf_bilap : requested workspace arrays unavailable')   ;   RETURN 
    89       ENDIF 
    90  
     86      ! 
     87      IF( nn_timing == 1 )  CALL timing_start('dyn_ldf_bilap') 
     88      ! 
     89      CALL wrk_alloc( jpi, jpj,      zcu, zcv           ) 
     90      CALL wrk_alloc( jpi, jpj, jpk, zuf, zut, zlu, zlv )  
     91      ! 
    9192      IF( kt == nit000 .AND. lwp ) THEN 
    9293         WRITE(numout,*) 
    … …  
    207208      END DO                                           !   End of slab 
    208209      !                                                ! =============== 
    209       IF( wrk_not_released(2, 1,2)     .OR.   & 
    210           wrk_not_released(3, 3,4,5,6) )   CALL ctl_stop('dyn_ldf_bilap: failed to release workspace arrays') 
     210      CALL wrk_dealloc( jpi, jpj,      zcu, zcv           ) 
     211      CALL wrk_dealloc( jpi, jpj, jpk, zuf, zut, zlu, zlv )  
     212      ! 
     213      IF( nn_timing == 1 )  CALL timing_stop('dyn_ldf_bilap') 
    211214      ! 
    212215   END SUBROUTINE dyn_ldf_bilap 

  • trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilapg.F90

    r2715 r3294  
    2727   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    2828   USE prtctl          ! Print control 
     29   USE wrk_nemo        ! Memory Allocation 
     30   USE timing          ! Timing 
     31 
    2932 
    3033   IMPLICIT NONE 
    … …  
    8487      !!              - save the trend in (zwk3,zwk4) ('key_trddyn') 
    8588      !!---------------------------------------------------------------------- 
    86       USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released 
    87       USE wrk_nemo, ONLY:   zwk1 => wrk_3d_3 , zwk2 => wrk_3d_4   ! 3D workspace 
    88       USE oce     , ONLY:   zwk3 => ta       , zwk4 => sa         ! ta, sa used as 3D workspace    
    89       ! 
    9089      INTEGER, INTENT( in ) ::   kt           ! ocean time-step index 
    9190      ! 
    9291      INTEGER ::   ji, jj, jk                 ! dummy loop indices 
    93       !!---------------------------------------------------------------------- 
    94  
    95       IF( wrk_in_use(3, 3,4) ) THEN 
    96          CALL ctl_stop('dyn_ldf_bilapg: requested workspace arrays unavailable')   ;   RETURN 
    97       ENDIF 
    98  
     92      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwk1, zwk2, zwk3, zwk4 
     93      !!---------------------------------------------------------------------- 
     94      ! 
     95      IF( nn_timing == 1 )  CALL timing_start('dyn_ldf_bilapg') 
     96      ! 
     97      CALL wrk_alloc( jpi, jpj, jpk, zwk1, zwk2, zwk3, zwk4 )  
     98      ! 
    9999      IF( kt == nit000 ) THEN 
    100100         IF(lwp) WRITE(numout,*) 
    101101         IF(lwp) WRITE(numout,*) 'dyn_ldf_bilapg : horizontal biharmonic operator in s-coordinate' 
    102102         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~' 
    103          zwk1(:,:,:) = 0.e0   ;   zwk3(:,:,:) = 0.e0 
    104          zwk2(:,:,:) = 0.e0   ;   zwk4(:,:,:) = 0.e0 
    105103         !                                      ! allocate dyn_ldf_bilapg arrays 
    106104         IF( dyn_ldf_bilapg_alloc() /= 0 )   CALL ctl_stop('STOP', 'dyn_ldf_bilapg: failed to allocate arrays') 
    107105      ENDIF 
     106      ! 
     107      zwk1(:,:,:) = 0.e0   ;   zwk3(:,:,:) = 0.e0 
     108      zwk2(:,:,:) = 0.e0   ;   zwk4(:,:,:) = 0.e0 
    108109 
    109110      ! Laplacian of (ub,vb) multiplied by ahm 
    … …  
    129130      END DO 
    130131      ! 
    131       IF( wrk_not_released(3, 3,4) )   CALL ctl_stop('dyn_ldf_bilapg: failed to release workspace arrays') 
     132      CALL wrk_dealloc( jpi, jpj, jpk, zwk1, zwk2, zwk3, zwk4 )  
     133      ! 
     134      IF( nn_timing == 1 )  CALL timing_stop('dyn_ldf_bilapg') 
    132135      ! 
    133136   END SUBROUTINE dyn_ldf_bilapg 
    … …  
    175178      !!      'key_trddyn' defined: the trend is saved for diagnostics. 
    176179      !!---------------------------------------------------------------------- 
    177       USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released 
    178       USE wrk_nemo, ONLY:   ziut => wrk_2d_1 , zjuf  => wrk_2d_2 , zjvt  => wrk_2d_3 
    179       USE wrk_nemo, ONLY:   zivf => wrk_2d_4 , zdku  => wrk_2d_5 , zdk1u => wrk_2d_6 
    180       USE wrk_nemo, ONLY:   zdkv => wrk_2d_7 , zdk1v => wrk_2d_8 
    181180      !! 
    182181      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pu , pv    ! 1st call: before horizontal velocity  
    … …  
    192191      REAL(wp) ::   zbur, zbvr, zmkt, zmkf, zuav, zvav   !   -      - 
    193192      REAL(wp) ::   zuwslpi, zuwslpj, zvwslpi, zvwslpj   !   -      - 
    194       !!---------------------------------------------------------------------- 
    195  
    196       IF( wrk_in_use(2, 1,2,3,4,5,6,7,8) ) THEN 
    197          CALL ctl_stop('dyn:ldfguv: requested workspace arrays unavailable')   ;   RETURN 
    198       END IF 
     193      ! 
     194      REAL(wp), POINTER, DIMENSION(:,:) :: ziut, zjuf, zjvt, zivf, zdku, zdk1u, zdkv, zdk1v 
     195      !!---------------------------------------------------------------------- 
     196      ! 
     197      IF( nn_timing == 1 )  CALL timing_start('ldfguv') 
     198      ! 
     199      CALL wrk_alloc( jpi, jpj, ziut, zjuf, zjvt, zivf, zdku, zdk1u, zdkv, zdk1v )  
     200      ! 
    199201      !                               ! ********** !   ! =============== 
    200202      DO jk = 1, jpkm1                ! First step !   ! Horizontal slab 
    … …  
    452454      !                                                ! =============== 
    453455 
    454       IF( wrk_not_released(2, 1,2,3,4,5,6,7,8) )   CALL ctl_stop('dyn:ldfguv: failed to release workspace arrays') 
     456      CALL wrk_dealloc( jpi, jpj, ziut, zjuf, zjvt, zivf, zdku, zdk1u, zdkv, zdk1v )  
     457      ! 
     458      IF( nn_timing == 1 )  CALL timing_stop('ldfguv') 
    455459      ! 
    456460   END SUBROUTINE ldfguv 

  • trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90

    r2715 r3294  
    2929   USE lib_mpp         ! MPP library 
    3030   USE prtctl          ! Print control 
     31   USE wrk_nemo        ! Memory Allocation 
     32   USE timing          ! Timing 
    3133 
    3234   IMPLICIT NONE 
    … …  
    105107      !!      of the rotated operator in dynzdf module 
    106108      !!---------------------------------------------------------------------- 
    107       USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released 
    108       USE wrk_nemo, ONLY:   ziut  => wrk_2d_1 , zjuf  => wrk_2d_2 , zjvt => wrk_2d_3    ! 2D workspace 
    109       USE wrk_nemo, ONLY:   zivf  => wrk_2d_4 , zdku  => wrk_2d_5 , zdkv => wrk_2d_6    ! 2D workspace 
    110       USE wrk_nemo, ONLY:   zdk1u => wrk_2d_7 , zdk1v => wrk_2d_8 
    111109      ! 
    112110      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index 
    … …  
    117115      REAL(wp) ::   zcoef0, zcoef3, zcoef4, zmkt, zmkf               !   -      - 
    118116      REAL(wp) ::   zuav, zvav, zuwslpi, zuwslpj, zvwslpi, zvwslpj   !   -      - 
     117      ! 
     118      REAL(wp), POINTER, DIMENSION(:,:) :: ziut, zjuf, zjvt, zivf, zdku, zdk1u, zdkv, zdk1v 
    119119      !!---------------------------------------------------------------------- 
    120  
    121       IF( wrk_in_use(2, 1,2,3,4,5,6,7,8) ) THEN 
    122          CALL ctl_stop('dyn_ldf_iso: requested workspace arrays unavailable')   ;   RETURN 
    123       END IF 
    124  
     120      ! 
     121      IF( nn_timing == 1 )  CALL timing_start('dyn_ldf_iso') 
     122      ! 
     123      CALL wrk_alloc( jpi, jpj, ziut, zjuf, zjvt, zivf, zdku, zdk1u, zdkv, zdk1v )  
     124      ! 
    125125      IF( kt == nit000 ) THEN 
    126126         IF(lwp) WRITE(numout,*) 
    … …  
    427427      END DO                                           !   End of slab 
    428428      !                                                ! =============== 
    429  
    430       IF( wrk_not_released(2, 1,2,3,4,5,6,7,8) )   CALL ctl_stop('dyn_ldf_iso: failed to release workspace arrays') 
     429      CALL wrk_dealloc( jpi, jpj, ziut, zjuf, zjvt, zivf, zdku, zdk1u, zdkv, zdk1v )  
     430      ! 
     431      IF( nn_timing == 1 )  CALL timing_stop('dyn_ldf_iso') 
    431432      ! 
    432433   END SUBROUTINE dyn_ldf_iso 

  • trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_lap.F90

    r2715 r3294  
    2323   USE trdmod_oce      ! ocean variables trends 
    2424   USE ldfslp          ! iso-neutral slopes  
     25   USE timing          ! Timing 
    2526 
    2627   IMPLICIT NONE 
    … …  
    6869      !!---------------------------------------------------------------------- 
    6970      ! 
     71      IF( nn_timing == 1 )  CALL timing_start('dyn_ldf_lap') 
     72      ! 
    7073      IF( kt == nit000 ) THEN 
    7174         IF(lwp) WRITE(numout,*) 
    … …  
    9598      END DO                                           !   End of slab 
    9699      !                                                ! =============== 
     100      IF( nn_timing == 1 )  CALL timing_stop('dyn_ldf_lap') 
     101      ! 
    97102   END SUBROUTINE dyn_ldf_lap 
    98103 

  • trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90

    r2779 r3294  
    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) 
    4041   USE lib_mpp         ! MPP library 
     42   USE wrk_nemo        ! Memory Allocation 
    4143   USE prtctl          ! Print control 
    4244#if defined key_agrif 
    4345   USE agrif_opa_interp 
    4446#endif 
     47   USE timing          ! Timing 
    4548 
    4649   IMPLICIT NONE 
    … …  
    7780      !!              * Apply lateral boundary conditions on after velocity  
    7881      !!             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 
     82      !!             at the one-way open boundaries (lk_obc=T), 
    8183      !!             at the AGRIF zoom     boundaries (lk_agrif=T) 
    8284      !! 
    … …  
    9294      !!               un,vn   now horizontal velocity of next time-step 
    9395      !!---------------------------------------------------------------------- 
    94       USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released 
    95       USE oce     , ONLY:   ze3u_f => ta       , ze3v_f => sa       ! (ta,sa) used as 3D workspace 
    96       USE wrk_nemo, ONLY:   zs_t   => wrk_2d_1 , zs_u_1 => wrk_2d_2 , zs_v_1 => wrk_2d_3 
    97       ! 
    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        !   -      - 
     105      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ze3u_f, ze3v_f  
    107106      !!---------------------------------------------------------------------- 
    108  
    109       IF( wrk_in_use(2, 1,2,3) ) THEN 
    110          CALL ctl_stop('dyn_nxt: requested workspace arrays unavailable')   ;   RETURN 
    111       ENDIF 
    112  
     107      ! 
     108      IF( nn_timing == 1 )  CALL timing_start('dyn_nxt') 
     109      ! 
     110      CALL wrk_alloc( jpi,jpj,jpk, ze3u_f, ze3v_f ) 
     111      ! 
    113112      IF( kt == nit000 ) THEN 
    114113         IF(lwp) WRITE(numout,*) 
    … …  
    174173      ENDIF 
    175174      ! 
    176 # elif defined key_bdy  
     175# elif defined key_bdy 
    177176      !                                !* BDY open boundaries 
    178       IF( .NOT. lk_dynspg_flt ) THEN 
    179          CALL bdy_dyn_frs( kt ) 
    180 #  if ! defined key_vvl 
    181          ua_e(:,:) = 0.e0 
    182          va_e(:,:) = 0.e0 
    183          ! Set these variables for use in bdy_dyn_fla 
    184          hur_e(:,:) = hur(:,:) 
    185          hvr_e(:,:) = hvr(:,:) 
    186          DO jk = 1, jpkm1   !! Vertically integrated momentum trends 
    187             ua_e(:,:) = ua_e(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk) 
    188             va_e(:,:) = va_e(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk) 
    189          END DO 
    190          ua_e(:,:) = ua_e(:,:) * hur(:,:) 
    191          va_e(:,:) = va_e(:,:) * hvr(:,:) 
    192          DO jk = 1 , jpkm1 
    193             ua(:,:,jk) = ua(:,:,jk) - ua_e(:,:) 
    194             va(:,:,jk) = va(:,:,jk) - va_e(:,:) 
    195          END DO 
    196          CALL bdy_dta_fla( kt+1, 0,2*nn_baro) 
    197          CALL bdy_dyn_fla( sshn_b ) 
    198          CALL lbc_lnk( ua_e, 'U', -1. )   ! Boundary points should be updated 
    199          CALL lbc_lnk( va_e, 'V', -1. )   ! 
    200          DO jk = 1 , jpkm1 
    201             ua(:,:,jk) = ( ua(:,:,jk) + ua_e(:,:) ) * umask(:,:,jk) 
    202             va(:,:,jk) = ( va(:,:,jk) + va_e(:,:) ) * vmask(:,:,jk) 
    203          END DO 
    204 #  endif 
    205       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      ! 
    206182# endif 
    207183      ! 
    … …  
    238214         ELSE                             ! Variable volume ! 
    239215            !                             ! ================! 
    240             ! Before scale factor at t-points 
    241             ! ------------------------------- 
    242             DO jk = 1, jpkm1 
     216            ! 
     217            DO jk = 1, jpkm1                 ! Before scale factor at t-points 
    243218               fse3t_b(:,:,jk) = fse3t_n(:,:,jk)                                   & 
    244219                  &              + atfp * (  fse3t_b(:,:,jk) + fse3t_a(:,:,jk)     & 
    245                   &                         - 2.e0 * fse3t_n(:,:,jk)            ) 
    246             ENDDO 
    247             ! Add volume filter correction only at the first level of t-point scale factors 
    248             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 
    249223            fse3t_b(:,:,1) = fse3t_b(:,:,1) - zec * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1) 
    250             ! surface at t-points and inverse surface at (u/v)-points used in surface averaging computations 
    251             zs_t  (:,:) =       e1t(:,:) * e2t(:,:) 
    252             zs_u_1(:,:) = 0.5 / ( e1u(:,:) * e2u(:,:) ) 
    253             zs_v_1(:,:) = 0.5 / ( e1v(:,:) * e2v(:,:) ) 
    254224            ! 
    255             IF( ln_dynadv_vec ) THEN 
    256                ! Before scale factor at (u/v)-points 
    257                ! ----------------------------------- 
    258                ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points 
    259                DO jk = 1, jpkm1 
    260                   DO jj = 1, jpjm1 
    261                      DO ji = 1, jpim1 
    262                         zv_t_ij           = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk) 
    263                         zv_t_ip1j         = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk) 
    264                         zv_t_ijp1         = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk) 
    265                         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) ) 
    266                         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) ) 
    267                      END DO 
    268                   END DO 
    269                END DO 
    270                ! lateral boundary conditions 
    271                CALL lbc_lnk( fse3u_b(:,:,:), 'U', 1. ) 
    272                CALL lbc_lnk( fse3v_b(:,:,:), 'V', 1. ) 
    273                ! Add initial scale factor to scale factor anomaly 
    274                fse3u_b(:,:,:) = fse3u_b(:,:,:) + fse3u_0(:,:,:) 
    275                fse3v_b(:,:,:) = fse3v_b(:,:,:) + fse3v_0(:,:,:) 
    276                ! Leap-Frog - Asselin filter and swap: applied on velocity 
    277                ! ----------------------------------- 
    278                DO jk = 1, jpkm1 
    279                   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                      !                                                 -------- 
    280232                     DO ji = 1, jpi 
    281233                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) ) 
    … …  
    290242               END DO 
    291243               ! 
    292             ELSE 
    293                ! Temporary filered scale factor at (u/v)-points (will become before scale factor) 
    294                !----------------------------------------------- 
    295                ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points 
    296                DO jk = 1, jpkm1 
    297                   DO jj = 1, jpjm1 
    298                      DO ji = 1, jpim1 
    299                         zv_t_ij          = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk) 
    300                         zv_t_ip1j        = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk) 
    301                         zv_t_ijp1        = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk) 
    302                         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) ) 
    303                         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) ) 
    304                      END DO 
    305                   END DO 
    306                END DO 
    307                ! lateral boundary conditions 
    308                CALL lbc_lnk( ze3u_f, 'U', 1. ) 
    309                CALL lbc_lnk( ze3v_f, 'V', 1. ) 
    310                ! Add initial scale factor to scale factor anomaly 
    311                ze3u_f(:,:,:) = ze3u_f(:,:,:) + fse3u_0(:,:,:) 
    312                ze3v_f(:,:,:) = ze3v_f(:,:,:) + fse3v_0(:,:,:) 
    313                ! Leap-Frog - Asselin filter and swap: applied on thickness weighted velocity 
    314                ! -----------------------------------             =========================== 
    315                DO jk = 1, jpkm1 
    316                   DO jj = 1, jpj 
    317                      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                 !                              --------------------------- 
    318251                        zue3a = ua(ji,jj,jk) * fse3u_a(ji,jj,jk) 
    319252                        zve3a = va(ji,jj,jk) * fse3v_a(ji,jj,jk) 
    … …  
    323256                        zve3b = vb(ji,jj,jk) * fse3v_b(ji,jj,jk) 
    324257                        ! 
    325                         zuf  = ( zue3n + atfp * ( zue3b - 2.e0 * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk) 
    326                         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) 
    327260                        ! 
    328                         ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity 
     261                        ub(ji,jj,jk) = zuf                     ! ub <-- filtered velocity 
    329262                        vb(ji,jj,jk) = zvf 
    330                         un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua 
     263                        un(ji,jj,jk) = ua(ji,jj,jk)            ! un <-- ua 
    331264                        vn(ji,jj,jk) = va(ji,jj,jk) 
    332265                     END DO 
    333266                  END DO 
    334267               END DO 
    335                fse3u_b(:,:,:) = ze3u_f(:,:,:)                   ! e3u_b <-- filtered scale factor 
    336                fse3v_b(:,:,:) = ze3v_f(:,:,:) 
    337                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 
    338271               CALL lbc_lnk( vb, 'V', -1. ) 
    339272            ENDIF 
    … …  
    346279         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask ) 
    347280      !  
    348       IF( wrk_not_released(2, 1,2,3) )   CALL ctl_stop('dyn_nxt: failed to release workspace arrays') 
     281      CALL wrk_dealloc( jpi,jpj,jpk, ze3u_f, ze3v_f ) 
     282      ! 
     283      IF( nn_timing == 1 )  CALL timing_stop('dyn_nxt') 
    349284      ! 
    350285   END SUBROUTINE dyn_nxt 

  • trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90

    r2715 r3294  
    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 
    … …  
    2928   USE lib_mpp        ! MPP library 
    3029   USE solver          ! solver initialization 
     30   USE wrk_nemo        ! Memory Allocation 
     31   USE timing          ! Timing 
     32 
    3133 
    3234   IMPLICIT NONE 
    … …  
    7476      !!        of the physical meaning of the results.  
    7577      !!---------------------------------------------------------------------- 
    76       USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released 
    77       USE wrk_nemo, ONLY:   ztrdu => wrk_3d_4 , ztrdv => wrk_3d_5    ! 3D workspace 
    7878      ! 
    7979      INTEGER, INTENT(in   ) ::   kt       ! ocean time-step index 
    … …  
    8282      INTEGER  ::   ji, jj, jk                             ! dummy loop indices 
    8383      REAL(wp) ::   z2dt, zg_2                             ! temporary scalar 
    84       !!---------------------------------------------------------------------- 
    85  
    86       IF( wrk_in_use(3, 4,5) ) THEN 
    87          CALL ctl_stop('dyn_spg: requested workspace arrays unavailable')   ;   RETURN 
    88       ENDIF 
     84      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdu, ztrdv 
     85      !!---------------------------------------------------------------------- 
     86      ! 
     87      IF( nn_timing == 1 )  CALL timing_start('dyn_spg') 
     88      ! 
    8989 
    9090!!gm NOTA BENE : the dynspg_exp and dynspg_ts should be modified so that  
    … …  
    9494 
    9595      IF( l_trddyn )   THEN                      ! temporary save of ta and sa trends 
     96         CALL wrk_alloc( jpi, jpj, jpk, ztrdu, ztrdv )  
    9697         ztrdu(:,:,:) = ua(:,:,:) 
    9798         ztrdv(:,:,:) = va(:,:,:) 
    … …  
    149150         END SELECT 
    150151         CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_spg, 'DYN', kt ) 
     152         ! 
     153         CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv )  
    151154      ENDIF 
    152155      !                                          ! print mean trends (used for debugging) 
    … …  
    154157         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    155158      ! 
    156       IF( wrk_not_released(3, 4,5) )   CALL ctl_stop('dyn_spg: failed to release workspace arrays') 
     159      IF( nn_timing == 1 )  CALL timing_stop('dyn_spg') 
    157160      ! 
    158161   END SUBROUTINE dyn_spg 
    … …  
    168171      INTEGER ::   ioptio 
    169172      !!---------------------------------------------------------------------- 
    170  
     173      ! 
     174      IF( nn_timing == 1 )  CALL timing_start('dyn_spg_init') 
     175      ! 
    171176      IF(lwp) THEN             ! Control print 
    172177         WRITE(numout,*) 
    … …  
    221226         IF( .NOT.ln_dynadv_vec )   CALL ctl_stop( 'Flux form not implemented for this free surface formulation' ) 
    222227      ENDIF 
    223  
    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 
     228      ! 
     229      IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_init') 
    233230      ! 
    234231   END SUBROUTINE dyn_spg_init 

  • trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_exp.F90

    r2715 r3294  
    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 
    … …  
    2828   USE iom             ! I/O library 
    2929   USE restart         ! only for lrst_oce 
     30   USE timing          ! Timing 
     31 
    3032 
    3133   IMPLICIT NONE 
    … …  
    6567      INTEGER ::   ji, jj, jk   ! dummy loop indices 
    6668      !!---------------------------------------------------------------------- 
    67  
     69      ! 
     70      IF( nn_timing == 1 )  CALL timing_start('dyn_spg_exp') 
     71      ! 
    6872      IF( kt == nit000 ) THEN 
    6973         IF(lwp) WRITE(numout,*) 
    … …  
    100104      ENDIF 
    101105      ! 
     106      IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_exp') 
     107      ! 
    102108   END SUBROUTINE dyn_spg_exp 
    103109 

  • trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_flt.F90

    r2715 r3294  
    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 
    4242   USE lib_mpp         ! distributed memory computing library 
     43   USE wrk_nemo        ! Memory Allocation 
    4344   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    4445   USE prtctl          ! Print control 
    … …  
    4950   USE agrif_opa_interp 
    5051#endif 
     52   USE timing          ! Timing 
    5153 
    5254   IMPLICIT NONE 
    … …  
    103105      !! References : Roullet and Madec 1999, JGR. 
    104106      !!--------------------------------------------------------------------- 
    105       USE oce, ONLY:   zub   => ta , zvb   => sa   ! (ta,sa) used as workspace 
    106       !! 
    107107      INTEGER, INTENT(in   ) ::   kt       ! ocean time-step index 
    108108      INTEGER, INTENT(  out) ::   kindic   ! solver convergence flag (<0 if not converge) 
    … …  
    110110      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    111111      REAL(wp) ::   z2dt, z2dtg, zgcb, zbtd, ztdgu, ztdgv   ! local scalars 
     112      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zub, zvb 
    112113      !!---------------------------------------------------------------------- 
     114      ! 
     115      IF( nn_timing == 1 )  CALL timing_start('dyn_spg_flt') 
     116      ! 
     117      CALL wrk_alloc( jpi,jpj,jpk, zub, zvb ) 
    113118      ! 
    114119      IF( kt == nit000 ) THEN 
    … …  
    187192#endif 
    188193#if defined key_bdy 
    189       CALL bdy_dyn_frs( kt )       ! Update velocities on unstructured boundary using the Flow Relaxation Scheme 
    190       CALL bdy_vol( kt )           ! Correction of the barotropic component velocity to control the volume of the system 
     194      CALL bdy_dyn( kt )      ! Update velocities on each open boundary 
     195      CALL bdy_vol( kt )      ! Correction of the barotropic component velocity to control the volume of the system 
    191196#endif 
    192197#if defined key_agrif 
    … …  
    304309#if defined key_obc 
    305310            ! caution : grad D = 0 along open boundaries 
     311            ! Remark: The filtering force could be reduced here in the FRS zone 
     312            !         by multiplying spgu/spgv by (1-alpha) ??   
    306313            spgu(ji,jj) = z2dt * ztdgu * obcumask(ji,jj) 
    307314            spgv(ji,jj) = z2dt * ztdgv * obcvmask(ji,jj) 
    308315#elif defined key_bdy 
    309316            ! 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) ??   
    312317            spgu(ji,jj) = z2dt * ztdgu * bdyumask(ji,jj) 
    313             spgv(ji,jj) = z2dt * ztdgv * bdyvmask(ji,jj)            
     318            spgv(ji,jj) = z2dt * ztdgv * bdyvmask(ji,jj) 
    314319#else 
    315320            spgu(ji,jj) = z2dt * ztdgu 
    … …  
    345350      ! -------------------------------------------------- 
    346351      IF( lrst_oce ) CALL flt_rst( kt, 'WRITE' ) 
     352      ! 
     353      CALL wrk_dealloc( jpi,jpj,jpk, zub, zvb ) 
     354      ! 
     355      IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_flt') 
    347356      ! 
    348357   END SUBROUTINE dyn_spg_flt 

  • trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_oce.F90

    r2715 r3294  
    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   !!---------------------------------------------------------------------- 

  • trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r2724 r3294  
    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 
     36   USE sbctide 
     37   USE updtide 
    3738   USE lib_mpp         ! distributed memory computing library 
    3839   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    … …  
    4142   USE iom             ! IOM library 
    4243   USE restart         ! only for lrst_oce 
    43    USE zdf_oce 
     44   USE zdf_oce         ! Vertical diffusion 
     45   USE wrk_nemo        ! Memory Allocation 
     46   USE timing          ! Timing 
     47 
    4448 
    4549   IMPLICIT NONE 
    … …  
    107111      !! References : Griffies et al., (2003): A technical guide to MOM4. NOAA/GFDL 
    108112      !!--------------------------------------------------------------------- 
    109       USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released 
    110       USE wrk_nemo, ONLY:   zsshun_e => wrk_2d_1 , zsshb_e  => wrk_2d_2  , zhdiv => wrk_2d_3 
    111       USE wrk_nemo, ONLY:   zsshvn_e => wrk_2d_4 , zssh_sum => wrk_2d_5 
    112       USE wrk_nemo, ONLY:   zcu => wrk_2d_6  , zwx   => wrk_2d_7  , zua   => wrk_2d_8  , zbfru  => wrk_2d_9 
    113       USE wrk_nemo, ONLY:   zcv => wrk_2d_10 , zwy   => wrk_2d_11 , zva   => wrk_2d_12 , zbfrv  => wrk_2d_13 
    114       USE wrk_nemo, ONLY:   zun => wrk_2d_14 , zun_e => wrk_2d_15 , zub_e => wrk_2d_16 , zu_sum => wrk_2d_17 
    115       USE wrk_nemo, ONLY:   zvn => wrk_2d_18 , zvn_e => wrk_2d_19 , zvb_e => wrk_2d_20 , zv_sum => wrk_2d_21 
    116113      ! 
    117114      INTEGER, INTENT(in)  ::   kt   ! ocean time-step index 
    … …  
    119116      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices 
    120117      INTEGER  ::   icycle           ! local scalar 
    121       REAL(wp) ::   zraur, zcoef, z2dt_e, z2dt_b     ! local scalars 
    122       REAL(wp) ::   z1_8, zx1, zy1                   !   -      - 
    123       REAL(wp) ::   z1_4, zx2, zy2                   !   -      - 
    124       REAL(wp) ::   zu_spg, zu_cor, zu_sld, zu_asp   !   -      - 
    125       REAL(wp) ::   zv_spg, zv_cor, zv_sld, zv_asp   !   -      - 
     118      INTEGER  ::   ikbu, ikbv       ! local scalar 
     119      REAL(wp) ::   zraur, zcoef, z2dt_e, z1_2dt_b, z2dt_bf   ! local scalars 
     120      REAL(wp) ::   z1_8, zx1, zy1                            !   -      - 
     121      REAL(wp) ::   z1_4, zx2, zy2                            !   -      - 
     122      REAL(wp) ::   zu_spg, zu_cor, zu_sld, zu_asp            !   -      - 
     123      REAL(wp) ::   zv_spg, zv_cor, zv_sld, zv_asp            !   -      - 
     124      REAL(wp) ::   ua_btm, va_btm                            !   -      - 
     125      ! 
     126      REAL(wp), POINTER, DIMENSION(:,:) :: zsshun_e, zsshvn_e, zsshb_e, zssh_sum, zhdiv  
     127      REAL(wp), POINTER, DIMENSION(:,:) :: zua, zva, zun, zvn, zun_e, zvn_e, zub_e, zvb_e  
     128      REAL(wp), POINTER, DIMENSION(:,:) :: zcu, zcv, zwx, zwy, zbfru, zbfrv, zu_sum, zv_sum 
    126129      !!---------------------------------------------------------------------- 
    127  
    128       IF( wrk_in_use(2,  1, 2, 3, 4, 5, 6, 7, 8, 9,10,     & 
    129          &              11,12,13,14,15,16,17,18,19,20,21 ) ) THEN 
    130          CALL ctl_stop( 'dyn_spg_ts: requested workspace arrays unavailable' )   ;   RETURN 
    131       ENDIF 
    132  
     130      ! 
     131      IF( nn_timing == 1 )  CALL timing_start('dyn_spg_ts') 
     132      ! 
     133      CALL wrk_alloc( jpi, jpj, zsshun_e, zsshvn_e, zsshb_e, zssh_sum, zhdiv     ) 
     134      CALL wrk_alloc( jpi, jpj, zua, zva, zun, zvn, zun_e, zvn_e, zub_e, zvb_e   ) 
     135      CALL wrk_alloc( jpi, jpj, zcu, zcv, zwx, zwy, zbfru, zbfrv, zu_sum, zv_sum ) 
     136      ! 
    133137      IF( kt == nit000 ) THEN             !* initialisation 
    134138         ! 
    … …  
    147151         hvr_e (:,:) = hvr  (:,:) 
    148152         IF( ln_dynvor_een ) THEN 
    149             ftne(1,:) = 0.e0   ;   ftnw(1,:) = 0.e0   ;   ftse(1,:) = 0.e0   ;   ftsw(1,:) = 0.e0 
     153            ftne(1,:) = 0._wp   ;   ftnw(1,:) = 0._wp   ;   ftse(1,:) = 0._wp   ;   ftsw(1,:) = 0._wp 
    150154            DO jj = 2, jpj 
    151155               DO ji = fs_2, jpi   ! vector opt. 
    152                   ftne(ji,jj) = ( ff(ji-1,jj  ) + ff(ji  ,jj  ) + ff(ji  ,jj-1) ) / 3. 
    153                   ftnw(ji,jj) = ( ff(ji-1,jj-1) + ff(ji-1,jj  ) + ff(ji  ,jj  ) ) / 3. 
    154                   ftse(ji,jj) = ( ff(ji  ,jj  ) + ff(ji  ,jj-1) + ff(ji-1,jj-1) ) / 3. 
    155                   ftsw(ji,jj) = ( ff(ji  ,jj-1) + ff(ji-1,jj-1) + ff(ji-1,jj  ) ) / 3. 
     156                  ftne(ji,jj) = ( ff(ji-1,jj  ) + ff(ji  ,jj  ) + ff(ji  ,jj-1) ) / 3._wp 
     157                  ftnw(ji,jj) = ( ff(ji-1,jj-1) + ff(ji-1,jj  ) + ff(ji  ,jj  ) ) / 3._wp 
     158                  ftse(ji,jj) = ( ff(ji  ,jj  ) + ff(ji  ,jj-1) + ff(ji-1,jj-1) ) / 3._wp 
     159                  ftsw(ji,jj) = ( ff(ji  ,jj-1) + ff(ji-1,jj-1) + ff(ji-1,jj  ) ) / 3._wp 
    156160               END DO 
    157161            END DO 
    … …  
    160164      ENDIF 
    161165 
    162       !                                   !* Local constant initialization 
    163       z2dt_b = 2.0 * rdt                                    ! baroclinic time step 
    164       z1_8 = 0.5 * 0.25                                     ! coefficient for vorticity estimates 
    165       z1_4 = 0.5 * 0.5 
    166       zraur  = 1. / rau0                                    ! 1 / volumic mass 
    167       ! 
    168       zhdiv(:,:) = 0.e0                                     ! barotropic divergence 
    169       zu_sld = 0.e0   ;   zu_asp = 0.e0                     ! tides trends (lk_tide=F) 
    170       zv_sld = 0.e0   ;   zv_asp = 0.e0 
     166      !                                                     !* Local constant initialization 
     167      z1_2dt_b = 1._wp / ( 2.0_wp * rdt )                   ! reciprocal of baroclinic time step 
     168      IF( neuler == 0 .AND. kt == nit000 )   z1_2dt_b = 1.0_wp / rdt    ! reciprocal of baroclinic  
     169                                                                        ! time step (euler timestep) 
     170      z1_8     = 0.125_wp                                   ! coefficient for vorticity estimates 
     171      z1_4     = 0.25_wp         
     172      zraur    = 1._wp / rau0                               ! 1 / volumic mass 
     173      ! 
     174      zhdiv(:,:) = 0._wp                                    ! barotropic divergence 
     175      zu_sld = 0._wp   ;   zu_asp = 0._wp                   ! tides trends (lk_tide=F) 
     176      zv_sld = 0._wp   ;   zv_asp = 0._wp 
     177 
     178      IF( kt == nit000 .AND. neuler == 0) THEN              ! for implicit bottom friction 
     179        z2dt_bf = rdt 
     180      ELSE 
     181        z2dt_bf = 2.0_wp * rdt 
     182      ENDIF 
    171183 
    172184      ! ----------------------------------------------------------------------------- 
    … …  
    176188      !                                   !* e3*d/dt(Ua), e3*Ub, e3*Vn (Vertically integrated) 
    177189      !                                   ! -------------------------- 
    178       zua(:,:) = 0.e0   ;   zun(:,:) = 0.e0   ;   ub_b(:,:) = 0.e0 
    179       zva(:,:) = 0.e0   ;   zvn(:,:) = 0.e0   ;   vb_b(:,:) = 0.e0 
     190      zua(:,:) = 0._wp   ;   zun(:,:) = 0._wp   ;   ub_b(:,:) = 0._wp 
     191      zva(:,:) = 0._wp   ;   zvn(:,:) = 0._wp   ;   vb_b(:,:) = 0._wp 
    180192      ! 
    181193      DO jk = 1, jpkm1 
    … …  
    195207               ! 
    196208#if defined key_vvl 
    197                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)  
    198                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)    
     209               ub_b(ji,jj) = ub_b(ji,jj) + fse3u_b(ji,jj,jk)* ub(ji,jj,jk)   *umask(ji,jj,jk)  
     210               vb_b(ji,jj) = vb_b(ji,jj) + fse3v_b(ji,jj,jk)* vb(ji,jj,jk)   *vmask(ji,jj,jk) 
    199211#else 
    200212               ub_b(ji,jj) = ub_b(ji,jj) + fse3u_0(ji,jj,jk) * ub(ji,jj,jk)  * umask(ji,jj,jk) 
    … …  
    270282      DO jj = 2, jpjm1                             ! Remove coriolis term (and possibly spg) from barotropic trend 
    271283         DO ji = fs_2, fs_jpim1 
    272             zua(ji,jj) = zua(ji,jj) - zcu(ji,jj) 
    273             zva(ji,jj) = zva(ji,jj) - zcv(ji,jj) 
    274          END DO 
     284             zua(ji,jj) = zua(ji,jj) - zcu(ji,jj) 
     285             zva(ji,jj) = zva(ji,jj) - zcv(ji,jj) 
     286          END DO 
    275287      END DO 
    276288 
    … …  
    278290      !                                             ! Remove barotropic contribution of bottom friction  
    279291      !                                             ! from the barotropic transport trend 
    280       zcoef = -1. / z2dt_b 
     292      zcoef = -1._wp * z1_2dt_b 
     293 
     294      IF(ln_bfrimp) THEN 
     295      !                                   ! Remove the bottom stress trend from 3-D sea surface level gradient 
     296      !                                   ! and Coriolis forcing in case of 3D semi-implicit bottom friction  
     297        DO jj = 2, jpjm1          
     298           DO ji = fs_2, fs_jpim1 
     299              ikbu = mbku(ji,jj) 
     300              ikbv = mbkv(ji,jj) 
     301              ua_btm = zcu(ji,jj) * z2dt_bf * hur(ji,jj) * umask (ji,jj,ikbu) 
     302              va_btm = zcv(ji,jj) * z2dt_bf * hvr(ji,jj) * vmask (ji,jj,ikbv) 
     303 
     304              zua(ji,jj) = zua(ji,jj) - bfrua(ji,jj) * ua_btm 
     305              zva(ji,jj) = zva(ji,jj) - bfrva(ji,jj) * va_btm 
     306           END DO 
     307        END DO 
     308 
     309      ELSE 
     310 
    281311# if defined key_vectopt_loop 
    282       DO jj = 1, 1 
    283          DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling) 
     312        DO jj = 1, 1 
     313           DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling) 
    284314# else 
    285       DO jj = 2, jpjm1 
    286          DO ji = 2, jpim1 
     315        DO jj = 2, jpjm1 
     316           DO ji = 2, jpim1 
    287317# endif 
    288318            ! Apply stability criteria for bottom friction 
    289319            !RBbug for vvl and external mode we may need to use varying fse3 
    290320            !!gm  Rq: the bottom e3 present the smallest variation, the use of e3u_0 is not a big approx. 
    291             zbfru(ji,jj) = MAX(  bfrua(ji,jj) , fse3u(ji,jj,mbku(ji,jj)) * zcoef  ) 
    292             zbfrv(ji,jj) = MAX(  bfrva(ji,jj) , fse3v(ji,jj,mbkv(ji,jj)) * zcoef  ) 
    293          END DO 
    294       END DO 
    295  
    296       IF( lk_vvl ) THEN 
    297          DO jj = 2, jpjm1 
    298             DO ji = fs_2, fs_jpim1   ! vector opt. 
    299                zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj)   & 
    300                   &       / ( hu_0(ji,jj) + sshu_b(ji,jj) + 1.e0 - umask(ji,jj,1) ) 
    301                zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj)   & 
    302                   &       / ( hv_0(ji,jj) + sshv_b(ji,jj) + 1.e0 - vmask(ji,jj,1) ) 
    303             END DO 
    304          END DO 
    305       ELSE 
    306          DO jj = 2, jpjm1 
    307             DO ji = fs_2, fs_jpim1   ! vector opt. 
    308                zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj) * hur(ji,jj) 
    309                zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj) * hvr(ji,jj) 
    310             END DO 
    311          END DO 
    312       ENDIF 
    313  
     321              zbfru(ji,jj) = MAX(  bfrua(ji,jj) , fse3u(ji,jj,mbku(ji,jj)) * zcoef  ) 
     322              zbfrv(ji,jj) = MAX(  bfrva(ji,jj) , fse3v(ji,jj,mbkv(ji,jj)) * zcoef  ) 
     323           END DO 
     324        END DO 
     325 
     326        IF( lk_vvl ) THEN 
     327           DO jj = 2, jpjm1 
     328              DO ji = fs_2, fs_jpim1   ! vector opt. 
     329                 zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj)   & 
     330                    &       / ( hu_0(ji,jj) + sshu_b(ji,jj) + 1._wp - umask(ji,jj,1) ) 
     331                 zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj)   & 
     332                    &       / ( hv_0(ji,jj) + sshv_b(ji,jj) + 1._wp - vmask(ji,jj,1) ) 
     333              END DO 
     334           END DO 
     335        ELSE 
     336           DO jj = 2, jpjm1 
     337              DO ji = fs_2, fs_jpim1   ! vector opt. 
     338                 zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj) * hur(ji,jj) 
     339                 zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj) * hvr(ji,jj) 
     340              END DO 
     341           END DO 
     342        ENDIF 
     343      END IF    ! end (ln_bfrimp) 
     344 
     345                     
    314346      !                                   !* d/dt(Ua), Ub, Vn (Vertical mean velocity) 
    315347      !                                   ! --------------------------  
    … …  
    318350      ! 
    319351      IF( lk_vvl ) THEN 
    320          ub_b(:,:) = ub_b(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1.e0 - umask(:,:,1) ) 
    321          vb_b(:,:) = vb_b(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1.e0 - vmask(:,:,1) ) 
     352         ub_b(:,:) = ub_b(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1._wp - umask(:,:,1) ) 
     353         vb_b(:,:) = vb_b(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1._wp - vmask(:,:,1) ) 
    322354      ELSE 
    323355         ub_b(:,:) = ub_b(:,:) * hur(:,:) 
    … …  
    355387      ! set ssh corrections to 0 
    356388      ! ssh corrections are applied to normal velocities (Flather's algorithm) and averaged over the barotropic loop 
    357       IF( lp_obc_east  )   sshfoe_b(:,:) = 0.e0 
    358       IF( lp_obc_west  )   sshfow_b(:,:) = 0.e0 
    359       IF( lp_obc_south )   sshfos_b(:,:) = 0.e0 
    360       IF( lp_obc_north )   sshfon_b(:,:) = 0.e0 
     389      IF( lp_obc_east  )   sshfoe_b(:,:) = 0._wp 
     390      IF( lp_obc_west  )   sshfow_b(:,:) = 0._wp 
     391      IF( lp_obc_south )   sshfos_b(:,:) = 0._wp 
     392      IF( lp_obc_north )   sshfon_b(:,:) = 0._wp 
    361393#endif 
    362394 
    … …  
    367399         IF( jn == 1 )   z2dt_e = rdt / nn_baro 
    368400 
    369          !                                                !* Update the forcing (OBC, BDY and tides) 
     401         !                                                !* Update the forcing (BDY and tides) 
    370402         !                                                !  ------------------ 
    371403         IF( lk_obc )   CALL obc_dta_bt ( kt, jn   ) 
    372          IF( lk_bdy )   CALL bdy_dta_fla( kt, jn+1, icycle ) 
     404         IF( lk_bdy )   CALL bdy_dta ( kt, jit=jn, time_offset=+1 ) 
     405         IF ( ln_tide_pot ) CALL upd_tide( kt, jn ) 
    373406 
    374407         !                                                !* after ssh_e 
    375408         !                                                !  ----------- 
    376          DO jj = 2, jpjm1                                      ! Horizontal divergence of barotropic transports 
     409         DO jj = 2, jpjm1                                 ! Horizontal divergence of barotropic transports 
    377410            DO ji = fs_2, fs_jpim1   ! vector opt. 
    378411               zhdiv(ji,jj) = (   e2u(ji  ,jj) * zun_e(ji  ,jj) * hu_e(ji  ,jj)     & 
    … …  
    386419         !                                                     ! OBC : zhdiv must be zero behind the open boundary 
    387420!!  mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column 
    388          IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1  ) = 0.e0      ! east 
    389          IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1  ) = 0.e0      ! west 
    390          IF( lp_obc_north )   zhdiv(nin0  :nin1  ,njn0p1:njn1p1) = 0.e0      ! north 
    391          IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1  ) = 0.e0      ! south 
     421         IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1  ) = 0._wp      ! east 
     422         IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1  ) = 0._wp      ! west 
     423         IF( lp_obc_north )   zhdiv(nin0  :nin1  ,njn0p1:njn1p1) = 0._wp      ! north 
     424         IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1  ) = 0._wp      ! south 
    392425#endif 
    393426#if defined key_bdy 
    … …  
    403436         !                                                !* after barotropic velocities (vorticity scheme dependent) 
    404437         !                                                !  ---------------------------   
    405          zwx(:,:) = e2u(:,:) * zun_e(:,:) * hu_e(:,:)           ! now_e transport 
     438         zwx(:,:) = e2u(:,:) * zun_e(:,:) * hu_e(:,:)     ! now_e transport 
    406439         zwy(:,:) = e1v(:,:) * zvn_e(:,:) * hv_e(:,:) 
    407440         ! 
    … …  
    419452                     zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj) 
    420453                  ENDIF 
     454                  ! add tidal astronomical forcing 
     455                  IF ( ln_tide_pot ) THEN  
     456                  zu_spg = zu_spg + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj) 
     457                  zv_spg = zv_spg + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj) 
     458                  ENDIF 
    421459                  ! energy conserving formulation for planetary vorticity term 
    422460                  zy1 = ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj) 
    … …  
    427465                  zv_cor =-z1_4 * ( ff(ji-1,jj  ) * zx1 + ff(ji,jj) * zx2 ) * hvr_e(ji,jj) 
    428466                  ! after velocities with implicit bottom friction 
    429                   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)   & 
    430                      &         / ( 1.e0         - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 
    431                   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)   & 
    432                      &         / ( 1.e0         - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 
     467 
     468                  IF( ln_bfrimp ) THEN      ! implicit bottom friction 
     469                     !   A new method to implement the implicit bottom friction.  
     470                     !   H. Liu 
     471                     !   Sept 2011 
     472                     ua_e(ji,jj) = umask(ji,jj,1) * ( zub_e(ji,jj) +                                            & 
     473                      &                               z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp )            & 
     474                      &                               / ( 1._wp      - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) ) 
     475                     ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e *   zua(ji,jj)  ) * umask(ji,jj,1)    
     476                     ! 
     477                     va_e(ji,jj) = vmask(ji,jj,1) * ( zvb_e(ji,jj) +                                            & 
     478                      &                               z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp )            & 
     479                      &                               / ( 1._wp      - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) ) 
     480                     va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e *   zva(ji,jj)  ) * vmask(ji,jj,1)    
     481                     ! 
     482                  ELSE 
     483                     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)   & 
     484                      &           / ( 1._wp         - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 
     485                     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)   & 
     486                      &           / ( 1._wp         - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 
     487                  ENDIF 
    433488               END DO 
    434489            END DO 
    … …  
    447502                     zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj) 
    448503                  ENDIF 
     504                  ! add tidal astronomical forcing 
     505                  IF ( ln_tide_pot ) THEN 
     506                  zu_spg = zu_spg + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj) 
     507                  zv_spg = zv_spg + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj) 
     508                  ENDIF 
    449509                  ! enstrophy conserving formulation for planetary vorticity term 
    450510                  zy1 =   z1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) + zwy(ji,jj) + zwy(ji+1,jj  ) ) / e1u(ji,jj) 
    … …  
    453513                  zv_cor  = zx1 * ( ff(ji-1,jj  ) + ff(ji,jj) ) * hvr_e(ji,jj) 
    454514                  ! after velocities with implicit bottom friction 
    455                   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)   & 
    456                      &         / ( 1.e0         - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 
    457                   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)   & 
    458                      &         / ( 1.e0         - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 
     515                  IF( ln_bfrimp ) THEN 
     516                     !   A new method to implement the implicit bottom friction.  
     517                     !   H. Liu 
     518                     !   Sept 2011 
     519                     ua_e(ji,jj) = umask(ji,jj,1) * ( zub_e(ji,jj) +                                            & 
     520                      &                               z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp )            & 
     521                      &                               / ( 1._wp      - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) ) 
     522                     ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e *   zua(ji,jj)  ) * umask(ji,jj,1)    
     523                     ! 
     524                     va_e(ji,jj) = vmask(ji,jj,1) * ( zvb_e(ji,jj) +                                            & 
     525                      &                               z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp )            & 
     526                      &                               / ( 1._wp      - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) ) 
     527                     va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e *   zva(ji,jj)  ) * vmask(ji,jj,1)    
     528                     ! 
     529                  ELSE 
     530                     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)   & 
     531                     &            / ( 1._wp        - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 
     532                     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)   & 
     533                     &            / ( 1._wp        - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 
     534                  ENDIF 
    459535               END DO 
    460536            END DO 
    … …  
    473549                     zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj) 
    474550                  ENDIF 
     551                  ! add tidal astronomical forcing 
     552                  IF ( ln_tide_pot ) THEN 
     553                  zu_spg = zu_spg + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj) 
     554                  zv_spg = zv_spg + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj) 
     555                  ENDIF 
    475556                  ! energy/enstrophy conserving formulation for planetary vorticity term 
    476557                  zu_cor = + z1_4 / e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) + ftnw(ji+1,jj) * zwy(ji+1,jj  )   & 
    … …  
    479560                     &                           + ftnw(ji,jj  ) * zwx(ji-1,jj  ) + ftne(ji,jj  ) * zwx(ji  ,jj  ) ) * hvr_e(ji,jj) 
    480561                  ! after velocities with implicit bottom friction 
    481                   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)   & 
    482                      &         / ( 1.e0         - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 
    483                   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)   & 
    484                      &         / ( 1.e0         - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 
     562                  IF( ln_bfrimp ) THEN 
     563                     !   A new method to implement the implicit bottom friction.  
     564                     !   H. Liu 
     565                     !   Sept 2011 
     566                     ua_e(ji,jj) = umask(ji,jj,1) * ( zub_e(ji,jj) +                                            & 
     567                      &                               z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp )            & 
     568                      &                               / ( 1._wp      - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) ) 
     569                     ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e *   zua(ji,jj)  ) * umask(ji,jj,1)    
     570                     ! 
     571                     va_e(ji,jj) = vmask(ji,jj,1) * ( zvb_e(ji,jj) +                                            & 
     572                      &                               z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp )            & 
     573                      &                               / ( 1._wp      - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) ) 
     574                     va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e *   zva(ji,jj)  ) * vmask(ji,jj,1)    
     575                     ! 
     576                  ELSE 
     577                     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)   & 
     578                     &            / ( 1._wp        - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 
     579                     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)   & 
     580                     &            / ( 1._wp        - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 
     581                  ENDIF 
    485582               END DO 
    486583            END DO 
    487584            !  
    488585         ENDIF 
    489          !                                                !* domain lateral boundary 
    490          !                                                !  ----------------------- 
    491          !                                                      ! Flather's boundary condition for the barotropic loop : 
    492          !                                                      !         - Update sea surface height on each open boundary 
    493          !                                                      !         - Correct the velocity 
    494  
     586         !                                                     !* domain lateral boundary 
     587         !                                                     !  ----------------------- 
     588 
     589                                                               ! OBC open boundaries 
    495590         IF( lk_obc               )   CALL obc_fla_ts ( ua_e, va_e, sshn_e, ssha_e ) 
    496          IF( lk_bdy .OR. ln_tides )   CALL bdy_dyn_fla( sshn_e )  
     591 
     592                                                               ! BDY open boundaries 
     593#if defined key_bdy 
     594         pssh => sshn_e 
     595         phur => hur_e 
     596         phvr => hvr_e 
     597         pu2d => ua_e 
     598         pv2d => va_e 
     599 
     600         IF( lk_bdy )   CALL bdy_dyn2d( kt )  
     601#endif 
     602 
    497603         ! 
    498604         CALL lbc_lnk( ua_e  , 'U', -1. )                      ! local domain boundaries  
    … …  
    526632            DO jj = 1, jpjm1                                    ! Sea Surface Height at u- & v-points 
    527633               DO ji = 1, fs_jpim1   ! Vector opt. 
    528                   zsshun_e(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )       & 
    529                      &                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn_e(ji  ,jj)    & 
    530                      &                    + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn_e(ji+1,jj) ) 
    531                   zsshvn_e(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )       & 
    532                      &                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn_e(ji,jj  )    & 
    533                      &                    + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn_e(ji,jj+1) ) 
     634                  zsshun_e(ji,jj) = 0.5_wp * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )       & 
     635                     &                     * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn_e(ji  ,jj)    & 
     636                     &                     +  e1t(ji+1,jj) * e2t(ji+1,jj) * sshn_e(ji+1,jj) ) 
     637                  zsshvn_e(ji,jj) = 0.5_wp * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )       & 
     638                     &                     * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn_e(ji,jj  )    & 
     639                     &                     +  e1t(ji,jj+1) * e2t(ji,jj+1) * sshn_e(ji,jj+1) ) 
    534640               END DO 
    535641            END DO 
    … …  
    539645            hu_e (:,:) = hu_0(:,:) + zsshun_e(:,:)              ! Ocean depth at U- and V-points 
    540646            hv_e (:,:) = hv_0(:,:) + zsshvn_e(:,:) 
    541             hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1.e0 - umask(:,:,1) ) 
    542             hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1.e0 - vmask(:,:,1) ) 
     647            hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1._wp - umask(:,:,1) ) 
     648            hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1._wp - vmask(:,:,1) ) 
    543649            ! 
    544650         ENDIF 
    … …  
    559665      ! 
    560666      !                                   !* Time average ==> after barotropic u, v, ssh 
    561       zcoef =  1.e0 / ( 2 * nn_baro  + 1 )  
     667      zcoef =  1._wp / ( 2 * nn_baro  + 1 )  
    562668      zu_sum(:,:) = zcoef * zu_sum  (:,:)  
    563669      zv_sum(:,:) = zcoef * zv_sum  (:,:)  
    … …  
    565671      !                                   !* update the general momentum trend 
    566672      DO jk=1,jpkm1 
    567          ua(:,:,jk) = ua(:,:,jk) + ( zu_sum(:,:) - ub_b(:,:) ) / z2dt_b 
    568          va(:,:,jk) = va(:,:,jk) + ( zv_sum(:,:) - vb_b(:,:) ) / z2dt_b 
     673         ua(:,:,jk) = ua(:,:,jk) + ( zu_sum(:,:) - ub_b(:,:) ) * z1_2dt_b 
     674         va(:,:,jk) = va(:,:,jk) + ( zv_sum(:,:) - vb_b(:,:) ) * z1_2dt_b 
    569675      END DO 
    570676      un_b  (:,:) =  zu_sum(:,:)  
    … …  
    575681      IF( lrst_oce )   CALL ts_rst( kt, 'WRITE' ) 
    576682      ! 
    577       ! 
    578       IF( wrk_not_released(2,  1, 2, 3, 4, 5, 6, 7, 8, 9,10,         & 
    579          &                    11,12,13,14,15,16,17,18,19,20,21) )    & 
    580          CALL ctl_stop('dyn_spg_ts: failed to release workspace arrays') 
     683      CALL wrk_dealloc( jpi, jpj, zsshun_e, zsshvn_e, zsshb_e, zssh_sum, zhdiv     ) 
     684      CALL wrk_dealloc( jpi, jpj, zua, zva, zun, zvn, zun_e, zvn_e, zub_e, zvb_e   ) 
     685      CALL wrk_dealloc( jpi, jpj, zcu, zcv, zwx, zwy, zbfru, zbfrv, zu_sum, zv_sum ) 
     686      ! 
     687      IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_ts') 
    581688      ! 
    582689   END SUBROUTINE dyn_spg_ts 
    … …  
    600707            CALL iom_get( numror, jpdom_autoglo, 'vn_b'  , vn_b  (:,:) )   ! from barotropic loop 
    601708         ELSE 
    602             un_b (:,:) = 0.e0 
    603             vn_b (:,:) = 0.e0 
     709            un_b (:,:) = 0._wp 
     710            vn_b (:,:) = 0._wp 
    604711            ! vertical sum 
    605712            IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll 
    … …  
    622729         ! Vertically integrated velocity (before) 
    623730         IF (neuler/=0) THEN 
    624             ub_b (:,:) = 0.e0 
    625             vb_b (:,:) = 0.e0 
     731            ub_b (:,:) = 0._wp 
     732            vb_b (:,:) = 0._wp 
    626733 
    627734            ! vertical sum 
    … …  
    641748 
    642749            IF( lk_vvl ) THEN 
    643                ub_b (:,:) = ub_b(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1.e0 - umask(:,:,1) ) 
    644                vb_b (:,:) = vb_b(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1.e0 - vmask(:,:,1) ) 
     750               ub_b (:,:) = ub_b(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1._wp - umask(:,:,1) ) 
     751               vb_b (:,:) = vb_b(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1._wp - vmask(:,:,1) ) 
    645752            ELSE 
    646753               ub_b(:,:) = ub_b(:,:) * hur(:,:) 

  • trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90

    r2715 r3294  
    2727   USE oce            ! ocean dynamics and tracers 
    2828   USE dom_oce        ! ocean space and time domain 
     29   USE dommsk         ! ocean mask 
    2930   USE dynadv         ! momentum advection (use ln_dynadv_vec value) 
    3031   USE trdmod         ! ocean dynamics trends  
    … …  
    3334   USE prtctl         ! Print control 
    3435   USE in_out_manager ! I/O manager 
    35    USE lib_mpp 
     36   USE lib_mpp        ! MPP library 
     37   USE wrk_nemo       ! Memory Allocation 
     38   USE timing         ! Timing 
     39 
    3640 
    3741   IMPLICIT NONE 
    … …  
    7175      !!               and planetary vorticity trends) ('key_trddyn') 
    7276      !!---------------------------------------------------------------------- 
    73       USE oce, ONLY:   ztrdu => ta , ztrdv => sa   ! (ta,sa) used as 3D workspace 
    74       ! 
    7577      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index 
    76       !!---------------------------------------------------------------------- 
     78      ! 
     79      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdu, ztrdv 
     80      !!---------------------------------------------------------------------- 
     81      ! 
     82      IF( nn_timing == 1 )  CALL timing_start('dyn_vor') 
     83      ! 
     84      IF( l_trddyn )   CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv ) 
    7785      ! 
    7886      !                                          ! vorticity term  
    … …  
    175183         &                     tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    176184      ! 
     185      IF( l_trddyn )   CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv ) 
     186      ! 
     187      IF( nn_timing == 1 )  CALL timing_stop('dyn_vor') 
     188      ! 
    177189   END SUBROUTINE dyn_vor 
    178190 
    … …  
    204216      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. 
    205217      !!---------------------------------------------------------------------- 
    206       USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released 
    207       USE wrk_nemo, ONLY:   zwx => wrk_2d_1 , zwy => wrk_2d_2 , zwz => wrk_2d_3     ! 2D workspace 
    208218      ! 
    209219      INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index 
    … …  
    215225      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    216226      REAL(wp) ::   zx1, zy1, zfact2, zx2, zy2   ! local scalars 
    217       !!---------------------------------------------------------------------- 
    218  
    219       IF( wrk_in_use(2, 1,2,3) ) THEN 
    220          CALL ctl_stop('dyn:vor_ene: requested workspace arrays unavailable')   ;   RETURN 
    221       ENDIF 
    222  
     227      REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz 
     228      !!---------------------------------------------------------------------- 
     229      ! 
     230      IF( nn_timing == 1 )  CALL timing_start('vor_ene') 
     231      ! 
     232      CALL wrk_alloc( jpi, jpj, zwx, zwy, zwz )  
     233      ! 
    223234      IF( kt == nit000 ) THEN 
    224235         IF(lwp) WRITE(numout,*) 
    … …  
    284295      END DO                                           !   End of slab 
    285296      !                                                ! =============== 
    286       IF( wrk_not_released(2, 1,2,3) )   CALL ctl_stop('dyn:vor_ene: failed to release workspace arrays') 
     297      CALL wrk_dealloc( jpi, jpj, zwx, zwy, zwz )  
     298      ! 
     299      IF( nn_timing == 1 )  CALL timing_stop('vor_ene') 
    287300      ! 
    288301   END SUBROUTINE vor_ene 
    … …  
    320333      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. 
    321334      !!---------------------------------------------------------------------- 
    322       USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released 
    323       USE wrk_nemo, ONLY:   zwx => wrk_2d_4 , zwy => wrk_2d_5 , zwz => wrk_2d_6 , zww => wrk_2d_7   ! 2D workspace 
    324335      ! 
    325336      INTEGER, INTENT(in) ::   kt   ! ocean timestep index 
    … …  
    328339      REAL(wp) ::   zfact1, zua, zcua, zx1, zy1   ! local scalars 
    329340      REAL(wp) ::   zfact2, zva, zcva, zx2, zy2   !   -      - 
    330       !!---------------------------------------------------------------------- 
    331  
    332       IF( wrk_in_use(2, 4,5,6,7) ) THEN 
    333          CALL ctl_stop('dyn:vor_mix: requested workspace arrays unavailable')   ;   RETURN 
    334       ENDIF 
    335  
     341      REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz, zww 
     342      !!---------------------------------------------------------------------- 
     343      ! 
     344      IF( nn_timing == 1 )  CALL timing_start('vor_mix') 
     345      ! 
     346      CALL wrk_alloc( jpi, jpj, zwx, zwy, zwz, zww )  
     347      ! 
    336348      IF( kt == nit000 ) THEN 
    337349         IF(lwp) WRITE(numout,*) 
    … …  
    404416      END DO                                           !   End of slab 
    405417      !                                                ! =============== 
    406       IF( wrk_not_released(2, 4,5,6,7) )   CALL ctl_stop('dyn:vor_mix: failed to release workspace arrays') 
     418      CALL wrk_dealloc( jpi, jpj, zwx, zwy, zwz, zww )  
     419      ! 
     420      IF( nn_timing == 1 )  CALL timing_stop('vor_mix') 
    407421      ! 
    408422   END SUBROUTINE vor_mix 
    … …  
    435449      !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. 
    436450      !!---------------------------------------------------------------------- 
    437       USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released 
    438       USE wrk_nemo, ONLY:   zwx => wrk_2d_4, zwy => wrk_2d_5, zwz => wrk_2d_6    ! 2D workspace 
    439451      ! 
    440452      INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index 
    … …  
    446458      INTEGER  ::   ji, jj, jk           ! dummy loop indices 
    447459      REAL(wp) ::   zfact1, zuav, zvau   ! temporary scalars 
    448       !!---------------------------------------------------------------------- 
    449        
    450       IF( wrk_in_use(2, 4,5,6) ) THEN 
    451          CALL ctl_stop('dyn:vor_ens: requested workspace arrays unavailable')   ;   RETURN 
    452       END IF 
    453  
     460      REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz, zww 
     461      !!---------------------------------------------------------------------- 
     462      ! 
     463      IF( nn_timing == 1 )  CALL timing_start('vor_ens') 
     464      ! 
     465      CALL wrk_alloc( jpi, jpj, zwx, zwy, zwz )  
     466      ! 
    454467      IF( kt == nit000 ) THEN 
    455468         IF(lwp) WRITE(numout,*) 
    … …  
    523536      END DO                                           !   End of slab 
    524537      !                                                ! =============== 
    525       IF( wrk_not_released(2, 4,5,6) )   CALL ctl_stop('dyn:vor_ens: failed to release workspace arrays') 
     538      CALL wrk_dealloc( jpi, jpj, zwx, zwy, zwz )  
     539      ! 
     540      IF( nn_timing == 1 )  CALL timing_stop('vor_ens') 
    526541      ! 
    527542   END SUBROUTINE vor_ens 
    … …  
    547562      !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36 
    548563      !!---------------------------------------------------------------------- 
    549       USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released 
    550       USE wrk_nemo, ONLY:   zwx  => wrk_2d_1 , zwy  => wrk_2d_2 ,  zwz => wrk_2d_3     ! 2D workspace 
    551       USE wrk_nemo, ONLY:   ztnw => wrk_2d_4 , ztne => wrk_2d_5  
    552       USE wrk_nemo, ONLY:   ztsw => wrk_2d_6 , ztse => wrk_2d_7 
    553 #if defined key_vvl 
    554       USE wrk_nemo, ONLY:   ze3f => wrk_3d_1                                           ! 3D workspace (lk_vvl=T) 
    555 #endif 
    556564      ! 
    557565      INTEGER , INTENT(in   )                         ::   kt     ! ocean time-step index 
    … …  
    561569      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pva    ! total v-trend 
    562570      !! 
    563       INTEGER  ::   ji, jj, jk         ! dummy loop indices 
    564       INTEGER  ::   ierr               ! local integer 
    565       REAL(wp) ::   zfac12, zua, zva   ! local scalars 
     571      INTEGER  ::   ji, jj, jk                                    ! dummy loop indices 
     572      INTEGER  ::   ierr                                          ! local integer 
     573      REAL(wp) ::   zfac12, zua, zva                              ! local scalars 
     574      !                                                           !  3D workspace  
     575      REAL(wp), POINTER    , DIMENSION(:,:  )         :: zwx, zwy, zwz 
     576      REAL(wp), POINTER    , DIMENSION(:,:  )         :: ztnw, ztne, ztsw, ztse 
     577#if defined key_vvl 
     578      REAL(wp), POINTER    , DIMENSION(:,:,:)         :: ze3f     !  3D workspace (lk_vvl=T) 
     579#endif 
    566580#if ! defined key_vvl 
    567       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), SAVE ::  ze3f     ! lk_vvl=F, ze3f=1/e3f saved one for all 
     581      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), SAVE   :: ze3f     ! lk_vvl=F, ze3f=1/e3f saved one for all 
    568582#endif 
    569583      !!---------------------------------------------------------------------- 
    570  
    571       IF( wrk_in_use(2, 1,2,3,4,5,6,7) .OR. wrk_in_use(3, 1) ) THEN 
    572          CALL ctl_stop('dyn:vor_een: requested workspace arrays unavailable')   ;   RETURN 
    573       ENDIF 
    574  
     584      ! 
     585      IF( nn_timing == 1 )  CALL timing_start('vor_een') 
     586      ! 
     587      CALL wrk_alloc( jpi, jpj,      zwx , zwy , zwz        )  
     588      CALL wrk_alloc( jpi, jpj,      ztnw, ztne, ztsw, ztse )  
     589#if defined key_vvl 
     590      CALL wrk_alloc( jpi, jpj, jpk, ze3f                   ) 
     591#endif 
     592      ! 
    575593      IF( kt == nit000 ) THEN 
    576594         IF(lwp) WRITE(numout,*) 
    … …  
    670688      END DO                                           !   End of slab 
    671689      !                                                ! =============== 
    672       IF( wrk_not_released(2, 1,2,3,4,5,6,7) .OR.   & 
    673           wrk_not_released(3, 1)             )   CALL ctl_stop('dyn:vor_een: failed to release workspace arrays') 
     690      CALL wrk_dealloc( jpi, jpj,      zwx , zwy , zwz        )  
     691      CALL wrk_dealloc( jpi, jpj,      ztnw, ztne, ztsw, ztse )  
     692#if defined key_vvl 
     693      CALL wrk_dealloc( jpi, jpj, jpk, ze3f                   ) 
     694#endif 
     695      ! 
     696      IF( nn_timing == 1 )  CALL timing_stop('vor_een') 
    674697      ! 
    675698   END SUBROUTINE vor_een 
    … …  
    684707      !!---------------------------------------------------------------------- 
    685708      INTEGER ::   ioptio          ! local integer 
     709      INTEGER ::   ji, jj, jk      ! dummy loop indices 
    686710      !! 
    687711      NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix, ln_dynvor_een 
    … …  
    700724         WRITE(numout,*) '           mixed enstrophy/energy conserving scheme   ln_dynvor_mix = ', ln_dynvor_mix 
    701725         WRITE(numout,*) '           enstrophy and energy conserving scheme     ln_dynvor_een = ', ln_dynvor_een 
     726      ENDIF 
     727 
     728      ! If energy, enstrophy or mixed advection of momentum in vector form change the value for masks 
     729      ! at angles with three ocean points and one land point 
     730      IF( ln_vorlat .AND. ( ln_dynvor_ene .OR. ln_dynvor_ens .OR. ln_dynvor_mix ) ) THEN 
     731         DO jk = 1, jpk 
     732            DO jj = 2, jpjm1 
     733               DO ji = 2, jpim1 
     734                  IF( tmask(ji,jj,jk)+tmask(ji+1,jj,jk)+tmask(ji,jj+1,jk)+tmask(ji+1,jj+1,jk) == 3._wp ) & 
     735                      fmask(ji,jj,jk) = 1._wp 
     736               END DO 
     737            END DO 
     738         END DO 
     739          ! 
     740          CALL lbc_lnk( fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask 
     741          ! 
    702742      ENDIF 
    703743 

  • trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90

    r2715 r3294  
    2121   USE lib_mpp         ! MPP library 
    2222   USE prtctl         ! Print control 
     23   USE wrk_nemo        ! Memory Allocation 
     24   USE timing          ! Timing 
    2325 
    2426   IMPLICIT NONE 
    … …  
    5254      !! ** Action  : - Update (ua,va) with the vert. momentum adv. trends 
    5355      !!              - Save the trends in (ztrdu,ztrdv) ('key_trddyn') 
    54      !!---------------------------------------------------------------------- 
    55       USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released 
    56       USE wrk_nemo, ONLY:   zww   => wrk_2d_1                        ! 2D workspace 
    57       USE oce     , ONLY:   zwuw  => ta       , zwvw  => sa          ! (ta,sa) used as 3D workspace 
    58       USE wrk_nemo, ONLY:   ztrdu => wrk_3d_1 , ztrdv => wrk_3d_2    ! 3D workspace 
    59       ! 
     56      !!---------------------------------------------------------------------- 
    6057      INTEGER, INTENT(in) ::   kt   ! ocean time-step inedx 
    6158      ! 
    6259      INTEGER  ::   ji, jj, jk      ! dummy loop indices 
    6360      REAL(wp) ::   zua, zva        ! temporary scalars 
     61      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwuw , zwvw 
     62      REAL(wp), POINTER, DIMENSION(:,:  ) ::  zww 
     63      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdu, ztrdv 
    6464      !!---------------------------------------------------------------------- 
    65        
    66       IF( wrk_in_use(2, 1) .OR. wrk_in_use(3, 1,2) ) THEN 
    67          CALL ctl_stop('dyn_zad: requested workspace arrays unavailable')   ;   RETURN 
    68       ENDIF 
    69  
     65      ! 
     66      IF( nn_timing == 1 )  CALL timing_start('dyn_zad') 
     67      ! 
     68      CALL wrk_alloc( jpi,jpj, zww )  
     69      CALL wrk_alloc( jpi,jpj,jpk, zwuw , zwvw )  
     70      ! 
    7071      IF( kt == nit000 ) THEN 
    7172         IF(lwp)WRITE(numout,*) 
    … …  
    7475 
    7576      IF( l_trddyn )   THEN         ! Save ua and va trends 
     77         CALL wrk_alloc( jpi, jpj, jpk, ztrdu, ztrdv )  
    7678         ztrdu(:,:,:) = ua(:,:,:)  
    7779         ztrdv(:,:,:) = va(:,:,:)  
    … …  
    117119         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    118120         CALL trd_mod(ztrdu, ztrdv, jpdyn_trd_zad, 'DYN', kt) 
     121         CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv )  
    119122      ENDIF 
    120123      !                             ! Control print 
    … …  
    122125         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    123126      ! 
    124       IF( wrk_not_released(2, 1)   .OR.   & 
    125           wrk_not_released(3, 1,2) )   CALL ctl_stop('dyn_zad: failed to release workspace arrays') 
     127      CALL wrk_dealloc( jpi,jpj, zww )  
     128      CALL wrk_dealloc( jpi,jpj,jpk, zwuw , zwvw )  
     129      ! 
     130      IF( nn_timing == 1 )  CALL timing_stop('dyn_zad') 
    126131      ! 
    127132   END SUBROUTINE dyn_zad 

  • trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf.F90

    r2715 r3294  
    2525   USE lib_mpp         ! MPP library 
    2626   USE prtctl          ! Print control 
     27   USE wrk_nemo        ! Memory Allocation 
     28   USE timing          ! Timing 
    2729 
    2830   IMPLICIT NONE 
    … …  
    5355      !! ** Purpose :   compute the vertical ocean dynamics physics. 
    5456      !!--------------------------------------------------------------------- 
    55       USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released 
    56       USE wrk_nemo, ONLY:   ztrdu => wrk_3d_1 , ztrdv => wrk_3d_2    ! 3D workspace 
    5757      !! 
    5858      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index 
     59      ! 
     60      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdu, ztrdv 
    5961      !!--------------------------------------------------------------------- 
    60  
    61       IF( wrk_in_use(3, 1,2) ) THEN 
    62          CALL ctl_stop('dyn_zdf: requested workspace arrays unavailable')   ;   RETURN 
    63       END IF 
     62      ! 
     63      IF( nn_timing == 1 )  CALL timing_start('dyn_zdf') 
     64      ! 
    6465      !                                          ! set time step 
    6566      IF( neuler == 0 .AND. kt == nit000     ) THEN   ;   r2dt =      rdt   ! = rdtra (restart with Euler time stepping) 
    … …  
    6869 
    6970      IF( l_trddyn )   THEN                      ! temporary save of ta and sa trends 
     71         CALL wrk_alloc( jpi, jpj, jpk, ztrdu, ztrdv )  
    7072         ztrdu(:,:,:) = ua(:,:,:) 
    7173         ztrdv(:,:,:) = va(:,:,:) 
    … …  
    9092         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
    9193         CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_zdf, 'DYN', kt ) 
     94         ! 
     95         CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv )  
    9296      ENDIF 
    9397      !                                          ! print mean trends (used for debugging) 
    … …  
    9599            &                    tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    96100      ! 
    97       IF( wrk_not_released(3, 1,2) )   CALL ctl_stop('dyn_zdf: failed to release workspace arrays') 
     101      IF( nn_timing == 1 )  CALL timing_stop('dyn_zdf') 
    98102      ! 
    99103   END SUBROUTINE dyn_zdf 

  • trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_exp.F90

    r2715 r3294  
    2222   USE in_out_manager  ! I/O manager 
    2323   USE lib_mpp         ! MPP library 
     24   USE wrk_nemo        ! Memory Allocation 
     25   USE timing          ! Timing 
     26 
    2427 
    2528   IMPLICIT NONE 
    … …  
    5457      !! ** Action : - Update (ua,va) with the vertical diffusive trend 
    5558      !!--------------------------------------------------------------------- 
    56       USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released 
    57       USE oce     , ONLY:   zwx => ta       , zwy => sa         ! (ta,sa) used as 3D workspace 
    58       USE wrk_nemo, ONLY:   zwz => wrk_3d_1 , zww => wrk_3d_2   ! 3D workspace 
    59       ! 
    6059      INTEGER , INTENT(in) ::   kt     ! ocean time-step index 
    6160      REAL(wp), INTENT(in) ::   p2dt   ! time-step  
    … …  
    6362      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices 
    6463      REAL(wp) ::   zrau0r, zlavmr, zua, zva   ! local scalars 
     64      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwx, zwy, zwz, zww 
    6565      !!---------------------------------------------------------------------- 
    66  
    67       IF( wrk_in_use(3, 1,2) ) THEN 
    68          CALL ctl_stop('dyn_zdf_exp: requested workspace arrays unavailable')   ;   RETURN 
    69       ENDIF 
    70  
     66      ! 
     67      IF( nn_timing == 1 )  CALL timing_start('dyn_zdf_exp') 
     68      ! 
     69      CALL wrk_alloc( jpi,jpj,jpk, zwx, zwy, zwz, zww )  
     70      ! 
    7171      IF( kt == nit000 .AND. lwp ) THEN 
    7272         WRITE(numout,*) 
    … …  
    120120      END DO                           ! End of time splitting 
    121121      ! 
    122       IF( wrk_not_released(3, 1,2) )   CALL ctl_stop('dyn_zdf_exp: failed to release workspace arrays') 
     122      CALL wrk_dealloc( jpi,jpj,jpk, zwx, zwy, zwz, zww )  
     123      ! 
     124      IF( nn_timing == 1 )  CALL timing_stop('dyn_zdf_exp') 
    123125      ! 
    124126   END SUBROUTINE dyn_zdf_exp 

  • trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90

    r2715 r3294  
    88   !!   NEMO     0.5  !  2002-08  (G. Madec)  F90: Free form and module 
    99   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  Forcing averaged over 2 time steps 
     10   !!            3.4  !  2012-01  (H. Liu) Semi-implicit bottom friction 
    1011   !!---------------------------------------------------------------------- 
    1112 
    … …  
    2021   USE in_out_manager  ! I/O manager 
    2122   USE lib_mpp         ! MPP library 
     23   USE zdfbfr          ! Bottom friction setup 
     24   USE wrk_nemo        ! Memory Allocation 
     25   USE timing          ! Timing 
    2226 
    2327   IMPLICIT NONE 
    … …  
    5458      !! ** Action : - Update (ua,va) arrays with the after vertical diffusive mixing trend. 
    5559      !!--------------------------------------------------------------------- 
    56       USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released 
    57       USE oce     , ONLY:  zwd  => ta       , zws   => sa   ! (ta,sa) used as 3D workspace 
    58       USE wrk_nemo, ONLY:   zwi => wrk_3d_3                 ! 3D workspace 
    59       !! 
    60       INTEGER , INTENT(in) ::   kt    ! ocean time-step index 
     60      INTEGER , INTENT(in) ::  kt     ! ocean time-step index 
    6161      REAL(wp), INTENT(in) ::  p2dt   ! vertical profile of tracer time-step 
    6262      !! 
    6363      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     64      INTEGER  ::   ikbu, ikbv   ! local integers 
    6465      REAL(wp) ::   z1_p2dt, zcoef, zzwi, zzws, zrhs   ! local scalars 
    6566      !!---------------------------------------------------------------------- 
    6667 
    67       IF( wrk_in_use(3, 3) ) THEN 
    68          CALL ctl_stop('dyn_zdf_imp: requested workspace array unavailable')   ;   RETURN 
    69       END IF 
    70  
     68      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwi, zwd, zws 
     69      REAL(wp), POINTER, DIMENSION(:,:)   ::  zavmu, zavmv 
     70      !!---------------------------------------------------------------------- 
     71      ! 
     72      IF( nn_timing == 1 )  CALL timing_start('dyn_zdf_imp') 
     73      ! 
     74      CALL wrk_alloc( jpi,jpj,jpk, zwi, zwd, zws )  
     75      CALL wrk_alloc( jpi,jpj, zavmu, zavmv )  
     76      ! 
    7177      IF( kt == nit000 ) THEN 
    7278         IF(lwp) WRITE(numout,*) 
    … …  
    7985      z1_p2dt = 1._wp / p2dt      ! inverse of the timestep 
    8086 
    81       ! 1. Vertical diffusion on u 
     87      ! 1. Apply semi-implicit bottom friction 
     88      ! -------------------------------------- 
     89      ! Only needed for semi-implicit bottom friction setup. The explicit 
     90      ! bottom friction has been included in "u(v)a" which act as the R.H.S 
     91      ! column vector of the tri-diagonal matrix equation 
     92      ! 
     93 
     94      IF( ln_bfrimp ) THEN 
     95# if defined key_vectopt_loop 
     96      DO jj = 1, 1 
     97         DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling) 
     98# else 
     99      DO jj = 2, jpjm1 
     100         DO ji = 2, jpim1 
     101# endif 
     102            ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points  
     103            ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points) 
     104            zavmu(ji,jj) = avmu(ji,jj,ikbu+1) 
     105            zavmv(ji,jj) = avmv(ji,jj,ikbv+1) 
     106            avmu(ji,jj,ikbu+1) = -bfrua(ji,jj) * fse3uw(ji,jj,ikbu+1)  
     107            avmv(ji,jj,ikbv+1) = -bfrva(ji,jj) * fse3vw(ji,jj,ikbv+1) 
     108         END DO 
     109      END DO 
     110      ENDIF 
     111 
     112      ! 2. Vertical diffusion on u 
    82113      ! --------------------------- 
    83114      ! Matrix and second member construction 
    84115      ! bottom boundary condition: both zwi and zws must be masked as avmu can take 
    85       ! non zero value at the ocean bottom depending on the bottom friction 
    86       ! used but the bottom velocities have already been updated with the bottom 
    87       ! friction velocity in dyn_bfr using values from the previous timestep. There 
    88       ! is no need to include these in the implicit calculation. 
     116      ! non zero value at the ocean bottom depending on the bottom friction used. 
    89117      ! 
    90118      DO jk = 1, jpkm1        ! Matrix 
    … …  
    168196 
    169197 
    170       ! 2. Vertical diffusion on v 
     198      ! 3. Vertical diffusion on v 
    171199      ! --------------------------- 
    172200      ! Matrix and second member construction 
    173201      ! bottom boundary condition: both zwi and zws must be masked as avmv can take 
    174       ! non zero value at the ocean bottom depending on the bottom friction 
    175       ! used but the bottom velocities have already been updated with the bottom 
    176       ! friction velocity in dyn_bfr using values from the previous timestep. There 
    177       ! is no need to include these in the implicit calculation. 
     202      ! non zero value at the ocean bottom depending on the bottom friction used 
    178203      ! 
    179204      DO jk = 1, jpkm1        ! Matrix 
    … …  
    255280         END DO 
    256281      END DO 
    257       ! 
    258       IF( wrk_not_released(3, 3) )   CALL ctl_stop('dyn_zdf_imp: failed to release workspace array') 
     282 
     283      !! restore bottom layer avmu(v)  
     284      IF( ln_bfrimp ) THEN 
     285# if defined key_vectopt_loop 
     286      DO jj = 1, 1 
     287         DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling) 
     288# else 
     289      DO jj = 2, jpjm1 
     290         DO ji = 2, jpim1 
     291# endif 
     292            ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points  
     293            ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points) 
     294            avmu(ji,jj,ikbu+1) = zavmu(ji,jj) 
     295            avmv(ji,jj,ikbv+1) = zavmv(ji,jj) 
     296         END DO 
     297      END DO 
     298      ENDIF 
     299      ! 
     300      CALL wrk_dealloc( jpi,jpj,jpk, zwi, zwd, zws)  
     301      CALL wrk_dealloc( jpi,jpj, zavmu, zavmv)  
     302      ! 
     303      IF( nn_timing == 1 )  CALL timing_stop('dyn_zdf_imp') 
    259304      ! 
    260305   END SUBROUTINE dyn_zdf_imp 

  • trunk/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90

    r2715 r3294  
    3939   USE asminc          ! Assimilation increment 
    4040#endif 
     41   USE wrk_nemo        ! Memory Allocation 
     42   USE timing          ! Timing 
    4143 
    4244   IMPLICIT NONE 
    … …  
    7577      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling. 
    7678      !!---------------------------------------------------------------------- 
    77       USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released 
    78       USE oce     , ONLY:   z3d   => ta                           ! ta used as 3D workspace 
    79       USE wrk_nemo, ONLY:   zhdiv => wrk_2d_1 , z2d => wrk_2d_2   ! 2D workspace 
    80       ! 
    8179      INTEGER, INTENT(in) ::   kt   ! time step 
    8280      ! 
    8381      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    8482      REAL(wp) ::   zcoefu, zcoefv, zcoeff, z2dt, z1_2dt, z1_rau0   ! local scalars 
    85       !!---------------------------------------------------------------------- 
    86  
    87       IF( wrk_in_use(2, 1,2) ) THEN 
    88          CALL ctl_stop('ssh_wzv: requested workspace arrays unavailable')   ;   RETURN 
    89       ENDIF 
    90  
     83      REAL(wp), POINTER, DIMENSION(:,:  ) ::  z2d, zhdiv 
     84      REAL(wp), POINTER, DIMENSION(:,:,:) ::  z3d 
     85      !!---------------------------------------------------------------------- 
     86      ! 
     87      IF( nn_timing == 1 )  CALL timing_start('ssh_wzv') 
     88      ! 
     89      CALL wrk_alloc( jpi, jpj, z2d, zhdiv )  
     90      ! 
    9191      IF( kt == nit000 ) THEN 
    9292         ! 
    … …  
    182182#if defined key_bdy 
    183183      ssha(:,:) = ssha(:,:) * bdytmask(:,:) 
    184       CALL lbc_lnk( ssha, 'T', 1. )  
     184      CALL lbc_lnk( ssha, 'T', 1. )                 ! absolutly compulsory !! (jmm) 
    185185#endif 
    186186 
    … …  
    230230      IF( lk_diaar5 ) THEN                            ! vertical mass transport & its square value 
    231231         ! Caution: in the VVL case, it only correponds to the baroclinic mass transport. 
     232         CALL wrk_alloc( jpi,jpj,jpk, z3d ) 
    232233         z2d(:,:) = rau0 * e1t(:,:) * e2t(:,:) 
    233234         DO jk = 1, jpk 
    … …  
    236237         CALL iom_put( "w_masstr" , z3d                     )   
    237238         CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) ) 
     239         CALL wrk_dealloc( jpi,jpj,jpk, z3d ) 
    238240      ENDIF 
    239241      ! 
    240242      IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha  - : ', mask1=tmask, ovlap=1 ) 
    241243      ! 
    242       IF( wrk_not_released(2, 1,2) )   CALL ctl_stop('ssh_wzv: failed to release workspace arrays') 
     244      CALL wrk_dealloc( jpi, jpj, z2d, zhdiv )  
     245      ! 
     246      IF( nn_timing == 1 )  CALL timing_stop('ssh_wzv') 
    243247      ! 
    244248   END SUBROUTINE ssh_wzv 
    … …  
    269273      REAL(wp) ::   zec      ! temporary scalar 
    270274      !!---------------------------------------------------------------------- 
    271  
     275      ! 
     276      IF( nn_timing == 1 )  CALL timing_start('ssh_nxt') 
     277      ! 
    272278      IF( kt == nit000 ) THEN 
    273279         IF(lwp) WRITE(numout,*) 
    … …  
    354360      IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb  - : ', mask1=tmask, ovlap=1 ) 
    355361      ! 
     362      IF( nn_timing == 1 )  CALL timing_stop('ssh_nxt') 
     363      ! 
    356364   END SUBROUTINE ssh_nxt 
    357365 
Note: See TracChangeset for help on using the changeset viewer.