Changeset 10928


Ignore:
Timestamp:
2019-05-03T17:44:56+02:00 (18 months ago)
Author:
davestorkey
Message:

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : Finish converting DYN module, including some updates to previously processed modules, but excluding dynnxt.F90 (which needs to be completely rewritten) and wet_dry.F90 - I need to talk to Enda about this.

Location:
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src
Files:
16 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM/istate.F90

    r10922 r10928  
    5151CONTAINS 
    5252 
    53    SUBROUTINE istate_init( Kbb, Kmm ) 
     53   SUBROUTINE istate_init( Kbb, Kmm, Kaa ) 
    5454      !!---------------------------------------------------------------------- 
    5555      !!                   ***  ROUTINE istate_init  *** 
     
    5757      !! ** Purpose :   Initialization of the dynamics and tracer fields. 
    5858      !!---------------------------------------------------------------------- 
    59       INTEGER, INTENT(in) ::   Kbb, Kmm   ! ocean time level indices 
     59      INTEGER, INTENT( in )  ::  Kbb, Kmm, Kaa   ! ocean time level indices 
     60      ! 
    6061      INTEGER ::   ji, jj, jk   ! dummy loop indices 
    6162!!gm see comment further down 
     
    7778      rhd  (:,:,:  ) = 0._wp   ;   rhop (:,:,:  ) = 0._wp      ! set one for all to 0 at level jpk 
    7879      rn2b (:,:,:  ) = 0._wp   ;   rn2  (:,:,:  ) = 0._wp      ! set one for all to 0 at levels 1 and jpk 
    79       tsa  (:,:,:,:) = 0._wp                                   ! set one for all to 0 at level jpk 
     80      ts  (:,:,:,:,Kaa) = 0._wp                                   ! set one for all to 0 at level jpk 
    8081      rab_b(:,:,:,:) = 0._wp   ;   rab_n(:,:,:,:) = 0._wp      ! set one for all to 0 at level jpk 
    8182#if defined key_agrif 
    82       ua   (:,:,:  ) = 0._wp   ! used in agrif_oce_sponge at initialization 
    83       va   (:,:,:  ) = 0._wp   ! used in agrif_oce_sponge at initialization     
     83      uu   (:,:,:  ,Kaa) = 0._wp   ! used in agrif_oce_sponge at initialization 
     84      vv   (:,:,:  ,Kaa) = 0._wp   ! used in agrif_oce_sponge at initialization     
    8485#endif 
    8586 
     
    9899         ! 
    99100         IF( ln_tsd_init ) THEN                
    100             CALL dta_tsd( nit000, tsb )       ! read 3D T and S data at nit000 
     101            CALL dta_tsd( nit000, ts(:,:,:,:,Kbb) )       ! read 3D T and S data at nit000 
    101102            ! 
    102             sshb(:,:)   = 0._wp               ! set the ocean at rest 
     103            ssh(:,:,Kbb)   = 0._wp               ! set the ocean at rest 
    103104            IF( ll_wd ) THEN 
    104                sshb(:,:) =  -ssh_ref  ! Added in 30 here for bathy that adds 30 as Iterative test CEOD  
     105               ssh(:,:,Kbb) =  -ssh_ref  ! Added in 30 here for bathy that adds 30 as Iterative test CEOD  
    105106               ! 
    106107               ! Apply minimum wetdepth criterion 
     
    108109               DO jj = 1,jpj 
    109110                  DO ji = 1,jpi 
    110                      IF( ht_0(ji,jj) + sshb(ji,jj)  < rn_wdmin1 ) THEN 
    111                         sshb(ji,jj) = tmask(ji,jj,1)*( rn_wdmin1 - (ht_0(ji,jj)) ) 
     111                     IF( ht_0(ji,jj) + ssh(ji,jj,Kbb)  < rn_wdmin1 ) THEN 
     112                        ssh(ji,jj,Kbb) = tmask(ji,jj,1)*( rn_wdmin1 - (ht_0(ji,jj)) ) 
    112113                     ENDIF 
    113114                  END DO 
    114115               END DO  
    115116            ENDIF  
    116             ub  (:,:,:) = 0._wp 
    117             vb  (:,:,:) = 0._wp   
     117            uu  (:,:,:,Kbb) = 0._wp 
     118            vv  (:,:,:,Kbb) = 0._wp   
    118119            ! 
    119120         ELSE                                 ! user defined initial T and S 
    120             CALL usr_def_istate( gdept_b, tmask, tsb, ub, vb, sshb  )          
     121            CALL usr_def_istate( gdept(:,:,:,Kbb), tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )          
    121122         ENDIF 
    122          tsn  (:,:,:,:) = tsb (:,:,:,:)       ! set now values from to before ones 
    123          sshn (:,:)     = sshb(:,:)    
    124          un   (:,:,:)   = ub  (:,:,:) 
    125          vn   (:,:,:)   = vb  (:,:,:) 
    126          hdivn(:,:,jpk) = 0._wp               ! bottom divergence set one for 0 to zero at jpk level 
    127          CALL div_hor( 0, Kmm )               ! compute interior hdivn value   
    128 !!gm                                    hdivn(:,:,:) = 0._wp 
     123         ts  (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)       ! set now values from to before ones 
     124         ssh (:,:,Kmm)     = ssh(:,:,Kbb)    
     125         uu   (:,:,:,Kmm)   = uu  (:,:,:,Kbb) 
     126         vv   (:,:,:,Kmm)   = vv  (:,:,:,Kbb) 
     127         hdiv(:,:,jpk) = 0._wp               ! bottom divergence set one for 0 to zero at jpk level 
     128         CALL div_hor( 0, Kmm )              ! compute interior hdiv value   
     129!!gm                                    hdiv(:,:,:) = 0._wp 
    129130 
    130131!!gm POTENTIAL BUG : 
    131 !!gm  ISSUE :  if sshb /= 0  then, in non linear free surface, the e3._n, e3._b should be recomputed 
     132!!gm  ISSUE :  if ssh(:,:,Kbb) /= 0  then, in non linear free surface, the e3._n, e3._b should be recomputed 
    132133!!             as well as gdept and gdepw....   !!!!!  
    133134!!      ===>>>>   probably a call to domvvl initialisation here.... 
     
    139140!            ALLOCATE( zuvd(jpi,jpj,jpk,2) ) 
    140141!            CALL dta_uvd( nit000, zuvd ) 
    141 !            ub(:,:,:) = zuvd(:,:,:,1) ;  un(:,:,:) = ub(:,:,:) 
    142 !            vb(:,:,:) = zuvd(:,:,:,2) ;  vn(:,:,:) = vb(:,:,:) 
     142!            uu(:,:,:,Kbb) = zuvd(:,:,:,1) ;  uu(:,:,:,Kmm) = uu(:,:,:,Kbb) 
     143!            vv(:,:,:,Kbb) = zuvd(:,:,:,2) ;  vv(:,:,:,Kmm) = vv(:,:,:,Kbb) 
    143144!            DEALLOCATE( zuvd ) 
    144145!         ENDIF 
    145146         ! 
    146147!!gm This is to be changed !!!! 
    147 !         ! - ML - sshn could be modified by istate_eel, so that initialization of e3t_b is done here 
     148!         ! - ML - ssh(:,:,Kmm) could be modified by istate_eel, so that initialization of e3t(:,:,:,Kbb) is done here 
    148149!         IF( .NOT.ln_linssh ) THEN 
    149150!            DO jk = 1, jpk 
    150 !               e3t_b(:,:,jk) = e3t_n(:,:,jk) 
     151!               e3t(:,:,jk,Kbb) = e3t(:,:,jk,Kmm) 
    151152!            END DO 
    152153!         ENDIF 
     
    158159      ! Do it whatever the free surface method, these arrays being eventually used 
    159160      ! 
    160       un_b(:,:) = 0._wp   ;   vn_b(:,:) = 0._wp 
    161       ub_b(:,:) = 0._wp   ;   vb_b(:,:) = 0._wp 
     161      uu_b(:,:,Kmm) = 0._wp   ;   vv_b(:,:,Kmm) = 0._wp 
     162      uu_b(:,:,Kbb) = 0._wp   ;   vv_b(:,:,Kbb) = 0._wp 
    162163      ! 
    163 !!gm  the use of umsak & vmask is not necessary below as un, vn, ub, vb are always masked 
     164!!gm  the use of umsak & vmask is not necessary below as uu(:,:,:,Kmm), vv(:,:,:,Kmm), uu(:,:,:,Kbb), vv(:,:,:,Kbb) are always masked 
    164165      DO jk = 1, jpkm1 
    165166         DO jj = 1, jpj 
    166167            DO ji = 1, jpi 
    167                un_b(ji,jj) = un_b(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk) 
    168                vn_b(ji,jj) = vn_b(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk) 
     168               uu_b(ji,jj,Kmm) = uu_b(ji,jj,Kmm) + e3u(ji,jj,jk,Kmm) * uu(ji,jj,jk,Kmm) * umask(ji,jj,jk) 
     169               vv_b(ji,jj,Kmm) = vv_b(ji,jj,Kmm) + e3v(ji,jj,jk,Kmm) * vv(ji,jj,jk,Kmm) * vmask(ji,jj,jk) 
    169170               ! 
    170                ub_b(ji,jj) = ub_b(ji,jj) + e3u_b(ji,jj,jk) * ub(ji,jj,jk) * umask(ji,jj,jk) 
    171                vb_b(ji,jj) = vb_b(ji,jj) + e3v_b(ji,jj,jk) * vb(ji,jj,jk) * vmask(ji,jj,jk) 
     171               uu_b(ji,jj,Kbb) = uu_b(ji,jj,Kbb) + e3u(ji,jj,jk,Kbb) * uu(ji,jj,jk,Kbb) * umask(ji,jj,jk) 
     172               vv_b(ji,jj,Kbb) = vv_b(ji,jj,Kbb) + e3v(ji,jj,jk,Kbb) * vv(ji,jj,jk,Kbb) * vmask(ji,jj,jk) 
    172173            END DO 
    173174         END DO 
    174175      END DO 
    175176      ! 
    176       un_b(:,:) = un_b(:,:) * r1_hu_n(:,:) 
    177       vn_b(:,:) = vn_b(:,:) * r1_hv_n(:,:) 
     177      uu_b(:,:,Kmm) = uu_b(:,:,Kmm) * r1_hu_n(:,:) 
     178      vv_b(:,:,Kmm) = vv_b(:,:,Kmm) * r1_hv_n(:,:) 
    178179      ! 
    179       ub_b(:,:) = ub_b(:,:) * r1_hu_b(:,:) 
    180       vb_b(:,:) = vb_b(:,:) * r1_hv_b(:,:) 
     180      uu_b(:,:,Kbb) = uu_b(:,:,Kbb) * r1_hu_b(:,:) 
     181      vv_b(:,:,Kbb) = vv_b(:,:,Kbb) * r1_hv_b(:,:) 
    181182      ! 
    182183   END SUBROUTINE istate_init 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/divhor.F90

    r10922 r10928  
    5555      !! 
    5656      !! ** Method  :   the now divergence is computed as : 
    57       !!         hdivn = 1/(e1e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] ) 
     57      !!         hdiv = 1/(e1e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] ) 
    5858      !!      and correct with runoff inflow (div_rnf) and cross land flow (div_cla)  
    5959      !! 
    60       !! ** Action  : - update hdivn, the now horizontal divergence 
     60      !! ** Action  : - update hdiv, the now horizontal divergence 
    6161      !!---------------------------------------------------------------------- 
    6262      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
    63       INTEGER, INTENT(in) ::   Kmm  ! ocean time-level index 
     63      INTEGER, INTENT(in) ::   Kmm  ! ocean time level index 
    6464      ! 
    6565      INTEGER  ::   ji, jj, jk    ! dummy loop indices 
     
    7878         DO jj = 2, jpjm1 
    7979            DO ji = fs_2, fs_jpim1   ! vector opt. 
    80                hdivn(ji,jj,jk) = (  e2u(ji  ,jj) * e3u_n(ji  ,jj,jk) * un(ji  ,jj,jk)      & 
    81                   &               - e2u(ji-1,jj) * e3u_n(ji-1,jj,jk) * un(ji-1,jj,jk)      & 
    82                   &               + e1v(ji,jj  ) * e3v_n(ji,jj  ,jk) * vn(ji,jj  ,jk)      & 
    83                   &               - e1v(ji,jj-1) * e3v_n(ji,jj-1,jk) * vn(ji,jj-1,jk)  )   & 
    84                   &            * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
     80               hdiv(ji,jj,jk) = (  e2u(ji  ,jj) * e3u(ji  ,jj,jk,Kmm) * uu(ji  ,jj,jk,Kmm)      & 
     81                  &               - e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) * uu(ji-1,jj,jk,Kmm)      & 
     82                  &               + e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm) * vv(ji,jj  ,jk,Kmm)      & 
     83                  &               - e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) * vv(ji,jj-1,jk,Kmm)  )   & 
     84                  &            * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
    8585            END DO   
    8686         END DO   
     
    8888#if defined key_agrif 
    8989      IF( .NOT. Agrif_Root() ) THEN 
    90          IF( nbondi == -1 .OR. nbondi == 2 )   hdivn(   2   ,  :   ,:) = 0._wp      ! west 
    91          IF( nbondi ==  1 .OR. nbondi == 2 )   hdivn( nlci-1,  :   ,:) = 0._wp      ! east 
    92          IF( nbondj == -1 .OR. nbondj == 2 )   hdivn(   :   ,  2   ,:) = 0._wp      ! south 
    93          IF( nbondj ==  1 .OR. nbondj == 2 )   hdivn(   :   ,nlcj-1,:) = 0._wp      ! north 
     90         IF( nbondi == -1 .OR. nbondi == 2 )   hdiv(   2   ,  :   ,:) = 0._wp      ! west 
     91         IF( nbondi ==  1 .OR. nbondi == 2 )   hdiv( nlci-1,  :   ,:) = 0._wp      ! east 
     92         IF( nbondj == -1 .OR. nbondj == 2 )   hdiv(   :   ,  2   ,:) = 0._wp      ! south 
     93         IF( nbondj ==  1 .OR. nbondj == 2 )   hdiv(   :   ,nlcj-1,:) = 0._wp      ! north 
    9494      ENDIF 
    9595#endif 
    9696      ! 
    97       IF( ln_rnf )   CALL sbc_rnf_div( hdivn, Kmm )                     !==  runoffs    ==!   (update hdivn field) 
     97      IF( ln_rnf )   CALL sbc_rnf_div( hdiv, Kmm )                     !==  runoffs    ==!   (update hdiv field) 
    9898      ! 
    9999#if defined key_asminc  
    100       IF( ln_sshinc .AND. ln_asmiau )   CALL ssh_asm_div( kt, hdivn )   !==  SSH assimilation  ==!   (update hdivn field) 
     100      IF( ln_sshinc .AND. ln_asmiau )   CALL ssh_asm_div( kt, hdiv )   !==  SSH assimilation  ==!   (update hdiv field) 
    101101      !  
    102102#endif 
    103       IF( ln_isf )   CALL sbc_isf_div( hdivn, Kmm )                     !==  ice shelf  ==!   (update hdivn field) 
     103      IF( ln_isf )   CALL sbc_isf_div( hdiv, Kmm )                     !==  ice shelf  ==!   (update hdiv field) 
    104104      ! 
    105       IF( ln_iscpl .AND. ln_hsb )   CALL iscpl_div( hdivn )             !==  ice sheet  ==!   (update hdivn field) 
     105      IF( ln_iscpl .AND. ln_hsb )   CALL iscpl_div( hdiv ) !==  ice sheet  ==!   (update hdiv field) 
    106106      ! 
    107       CALL lbc_lnk( 'divhor', hdivn, 'T', 1. )   !   (no sign change) 
     107      CALL lbc_lnk( 'divhor', hdiv, 'T', 1. )   !   (no sign change) 
    108108      ! 
    109109      IF( ln_timing )   CALL timing_stop('div_hor') 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynadv_cen2.F90

    r10877 r10928  
    141141      ENDIF 
    142142      !                                   ! Control print 
    143       IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' cen2 adv - Ua: ', mask1=umask,   & 
    144          &                       tab3d_2=va, clinfo2=           ' Va: ', mask2=vmask, clinfo3='dyn' ) 
     143      IF(ln_ctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' cen2 adv - Ua: ', mask1=umask,   & 
     144         &                       tab3d_2=pvv(:,:,:,Krhs), clinfo2=           ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    145145      ! 
    146146   END SUBROUTINE dyn_adv_cen2 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynadv_ubs.F90

    r10877 r10928  
    234234      ENDIF 
    235235      !                                         ! Control print 
    236       IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' ubs2 adv - Ua: ', mask1=umask,   & 
    237          &                       tab3d_2=va, clinfo2=           ' Va: ', mask2=vmask, clinfo3='dyn' ) 
     236      IF(ln_ctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' ubs2 adv - Ua: ', mask1=umask,   & 
     237         &                       tab3d_2=pvv(:,:,:,Krhs), clinfo2=           ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    238238      ! 
    239239   END SUBROUTINE dyn_adv_ubs 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynhpg.F90

    r10491 r10928  
    8181CONTAINS 
    8282 
    83    SUBROUTINE dyn_hpg( kt ) 
     83   SUBROUTINE dyn_hpg( kt, Kmm, puu, pvv, Krhs ) 
    8484      !!--------------------------------------------------------------------- 
    8585      !!                  ***  ROUTINE dyn_hpg  *** 
     
    8888      !!              using the scheme defined in the namelist 
    8989      !! 
    90       !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 
     90      !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend 
    9191      !!             - send trends to trd_dyn for futher diagnostics (l_trddyn=T) 
    9292      !!---------------------------------------------------------------------- 
    93       INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     93      INTEGER                             , INTENT( in )  ::  kt          ! ocean time-step index 
     94      INTEGER                             , INTENT( in )  ::  Kmm, Krhs   ! ocean time level indices 
     95      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv    ! ocean velocities and RHS of momentum equation 
     96      ! 
    9497      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdu, ztrdv 
    9598      !!---------------------------------------------------------------------- 
     
    97100      IF( ln_timing )   CALL timing_start('dyn_hpg') 
    98101      ! 
    99       IF( l_trddyn ) THEN                    ! Temporary saving of ua and va trends (l_trddyn) 
     102      IF( l_trddyn ) THEN                    ! Temporary saving of puu(:,:,:,Krhs) and pvv(:,:,:,Krhs) trends (l_trddyn) 
    100103         ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) ) 
    101          ztrdu(:,:,:) = ua(:,:,:) 
    102          ztrdv(:,:,:) = va(:,:,:) 
     104         ztrdu(:,:,:) = puu(:,:,:,Krhs) 
     105         ztrdv(:,:,:) = pvv(:,:,:,Krhs) 
    103106      ENDIF 
    104107      ! 
    105108      SELECT CASE ( nhpg )      ! Hydrostatic pressure gradient computation 
    106       CASE ( np_zco )   ;   CALL hpg_zco    ( kt )      ! z-coordinate 
    107       CASE ( np_zps )   ;   CALL hpg_zps    ( kt )      ! z-coordinate plus partial steps (interpolation) 
    108       CASE ( np_sco )   ;   CALL hpg_sco    ( kt )      ! s-coordinate (standard jacobian formulation) 
    109       CASE ( np_djc )   ;   CALL hpg_djc    ( kt )      ! s-coordinate (Density Jacobian with Cubic polynomial) 
    110       CASE ( np_prj )   ;   CALL hpg_prj    ( kt )      ! s-coordinate (Pressure Jacobian scheme) 
    111       CASE ( np_isf )   ;   CALL hpg_isf    ( kt )      ! s-coordinate similar to sco modify for ice shelf 
     109      CASE ( np_zco )   ;   CALL hpg_zco    ( kt, Kmm, puu, pvv, Krhs )  ! z-coordinate 
     110      CASE ( np_zps )   ;   CALL hpg_zps    ( kt, Kmm, puu, pvv, Krhs )  ! z-coordinate plus partial steps (interpolation) 
     111      CASE ( np_sco )   ;   CALL hpg_sco    ( kt, Kmm, puu, pvv, Krhs )  ! s-coordinate (standard jacobian formulation) 
     112      CASE ( np_djc )   ;   CALL hpg_djc    ( kt, Kmm, puu, pvv, Krhs )  ! s-coordinate (Density Jacobian with Cubic polynomial) 
     113      CASE ( np_prj )   ;   CALL hpg_prj    ( kt, Kmm, puu, pvv, Krhs )  ! s-coordinate (Pressure Jacobian scheme) 
     114      CASE ( np_isf )   ;   CALL hpg_isf    ( kt, Kmm, puu, pvv, Krhs )  ! s-coordinate similar to sco modify for ice shelf 
    112115      END SELECT 
    113116      ! 
    114117      IF( l_trddyn ) THEN      ! save the hydrostatic pressure gradient trends for momentum trend diagnostics 
    115          ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    116          ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
     118         ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:) 
     119         ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:) 
    117120         CALL trd_dyn( ztrdu, ztrdv, jpdyn_hpg, kt ) 
    118121         DEALLOCATE( ztrdu , ztrdv ) 
    119122      ENDIF 
    120123      ! 
    121       IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' hpg  - Ua: ', mask1=umask,   & 
    122          &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
     124      IF(ln_ctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' hpg  - Ua: ', mask1=umask,   & 
     125         &                       tab3d_2=pvv(:,:,:,Krhs), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    123126      ! 
    124127      IF( ln_timing )   CALL timing_stop('dyn_hpg') 
     
    127130 
    128131 
    129    SUBROUTINE dyn_hpg_init 
     132   SUBROUTINE dyn_hpg_init( Kmm ) 
    130133      !!---------------------------------------------------------------------- 
    131134      !!                 ***  ROUTINE dyn_hpg_init  *** 
     
    137140      !!      with the type of vertical coordinate used (zco, zps, sco) 
    138141      !!---------------------------------------------------------------------- 
     142      INTEGER, INTENT( in ) :: Kmm   ! ocean time level index 
     143      ! 
    139144      INTEGER ::   ioptio = 0      ! temporary integer 
    140145      INTEGER ::   ios             ! Local integer output status for namelist read 
     
    228233         ! 
    229234         DO jk = 1, jpk                   !- compute density of the water displaced by the ice shelf  
    230             CALL eos( zts_top(:,:,:), gdept_n(:,:,jk), zrhd(:,:,jk) ) 
     235            CALL eos( zts_top(:,:,:), gdept(:,:,jk,Kmm), zrhd(:,:,jk) ) 
    231236         END DO 
    232237         ! 
     
    239244            DO ji = 1, jpi                      ! divided by 2 later 
    240245               ikt = mikt(ji,jj) 
    241                ziceload(ji,jj) = ziceload(ji,jj) + (znad + zrhd(ji,jj,1) ) * e3w_n(ji,jj,1) * (1._wp - tmask(ji,jj,1)) 
     246               ziceload(ji,jj) = ziceload(ji,jj) + (znad + zrhd(ji,jj,1) ) * e3w(ji,jj,1,Kmm) * (1._wp - tmask(ji,jj,1)) 
    242247               DO jk = 2, ikt-1 
    243                   ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhd(ji,jj,jk-1) + zrhd(ji,jj,jk)) * e3w_n(ji,jj,jk) & 
     248                  ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhd(ji,jj,jk-1) + zrhd(ji,jj,jk)) * e3w(ji,jj,jk,Kmm) & 
    244249                     &                              * (1._wp - tmask(ji,jj,jk)) 
    245250               END DO 
    246251               IF (ikt  >=  2) ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + zrhd(ji,jj,ikt-1)) & 
    247                   &                                              * ( risfdep(ji,jj) - gdept_n(ji,jj,ikt-1) ) 
     252                  &                                              * ( risfdep(ji,jj) - gdept(ji,jj,ikt-1,Kmm) ) 
    248253            END DO 
    249254         END DO 
     
    256261 
    257262 
    258    SUBROUTINE hpg_zco( kt ) 
     263   SUBROUTINE hpg_zco( kt, Kmm, puu, pvv, Krhs ) 
    259264      !!--------------------------------------------------------------------- 
    260265      !!                  ***  ROUTINE hpg_zco  *** 
     
    266271      !!      level:    zhpi = grav ..... 
    267272      !!                zhpj = grav ..... 
    268       !!      add it to the general momentum trend (ua,va). 
    269       !!            ua = ua - 1/e1u * zhpi 
    270       !!            va = va - 1/e2v * zhpj 
    271       !! 
    272       !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 
    273       !!---------------------------------------------------------------------- 
    274       INTEGER, INTENT(in) ::   kt    ! ocean time-step index 
     273      !!      add it to the general momentum trend (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)). 
     274      !!            puu(:,:,:,Krhs) = puu(:,:,:,Krhs) - 1/e1u * zhpi 
     275      !!            pvv(:,:,:,Krhs) = pvv(:,:,:,Krhs) - 1/e2v * zhpj 
     276      !! 
     277      !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend 
     278      !!---------------------------------------------------------------------- 
     279      INTEGER                             , INTENT( in )  ::  kt          ! ocean time-step index 
     280      INTEGER                             , INTENT( in )  ::  Kmm, Krhs   ! ocean time level indices 
     281      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv    ! ocean velocities and RHS of momentum equation 
    275282      ! 
    276283      INTEGER  ::   ji, jj, jk       ! dummy loop indices 
     
    290297      DO jj = 2, jpjm1 
    291298         DO ji = fs_2, fs_jpim1   ! vector opt. 
    292             zcoef1 = zcoef0 * e3w_n(ji,jj,1) 
     299            zcoef1 = zcoef0 * e3w(ji,jj,1,Kmm) 
    293300            ! hydrostatic pressure gradient 
    294301            zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj) 
    295302            zhpj(ji,jj,1) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj) 
    296303            ! add to the general momentum trend 
    297             ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) 
    298             va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) 
     304            puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + zhpi(ji,jj,1) 
     305            pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + zhpj(ji,jj,1) 
    299306         END DO 
    300307      END DO 
     
    305312         DO jj = 2, jpjm1 
    306313            DO ji = fs_2, fs_jpim1   ! vector opt. 
    307                zcoef1 = zcoef0 * e3w_n(ji,jj,jk) 
     314               zcoef1 = zcoef0 * e3w(ji,jj,jk,Kmm) 
    308315               ! hydrostatic pressure gradient 
    309316               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1)   & 
     
    315322                  &                       - ( rhd(ji,jj,  jk)+rhd(ji,jj  ,jk-1) )  ) * r1_e2v(ji,jj) 
    316323               ! add to the general momentum trend 
    317                ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) 
    318                va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) 
     324               puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zhpi(ji,jj,jk) 
     325               pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zhpj(ji,jj,jk) 
    319326            END DO 
    320327         END DO 
     
    324331 
    325332 
    326    SUBROUTINE hpg_zps( kt ) 
     333   SUBROUTINE hpg_zps( kt, Kmm, puu, pvv, Krhs ) 
    327334      !!--------------------------------------------------------------------- 
    328335      !!                 ***  ROUTINE hpg_zps  *** 
     
    330337      !! ** Method  :   z-coordinate plus partial steps case.  blahblah... 
    331338      !! 
    332       !! ** Action  : - Update (ua,va) with the now hydrastatic pressure trend 
    333       !!---------------------------------------------------------------------- 
    334       INTEGER, INTENT(in) ::   kt    ! ocean time-step index 
     339      !! ** Action  : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend 
     340      !!---------------------------------------------------------------------- 
     341      INTEGER                             , INTENT( in )  ::  kt          ! ocean time-step index 
     342      INTEGER                             , INTENT( in )  ::  Kmm, Krhs   ! ocean time level indices 
     343      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv    ! ocean velocities and RHS of momentum equation 
    335344      !! 
    336345      INTEGER  ::   ji, jj, jk                       ! dummy loop indices 
     
    347356 
    348357      ! Partial steps: bottom before horizontal gradient of t, s, rd at the last ocean level 
    349 !jc      CALL zps_hde    ( kt, jpts, tsn, gtsu, gtsv, rhd, gru , grv ) 
     358!jc      CALL zps_hde    ( kt, jpts, ts(:,:,:,:,Kmm), gtsu, gtsv, rhd, gru , grv ) 
    350359 
    351360      ! Local constant initialization 
     
    355364      DO jj = 2, jpjm1 
    356365         DO ji = fs_2, fs_jpim1   ! vector opt. 
    357             zcoef1 = zcoef0 * e3w_n(ji,jj,1) 
     366            zcoef1 = zcoef0 * e3w(ji,jj,1,Kmm) 
    358367            ! hydrostatic pressure gradient 
    359368            zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj  ,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj) 
    360369            zhpj(ji,jj,1) = zcoef1 * ( rhd(ji  ,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj) 
    361370            ! add to the general momentum trend 
    362             ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) 
    363             va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) 
     371            puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + zhpi(ji,jj,1) 
     372            pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + zhpj(ji,jj,1) 
    364373         END DO 
    365374      END DO 
     
    369378         DO jj = 2, jpjm1 
    370379            DO ji = fs_2, fs_jpim1   ! vector opt. 
    371                zcoef1 = zcoef0 * e3w_n(ji,jj,jk) 
     380               zcoef1 = zcoef0 * e3w(ji,jj,jk,Kmm) 
    372381               ! hydrostatic pressure gradient 
    373382               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1)   & 
     
    379388                  &                       - ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) )  ) * r1_e2v(ji,jj) 
    380389               ! add to the general momentum trend 
    381                ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) 
    382                va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) 
     390               puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zhpi(ji,jj,jk) 
     391               pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zhpj(ji,jj,jk) 
    383392            END DO 
    384393         END DO 
     
    390399            iku = mbku(ji,jj) 
    391400            ikv = mbkv(ji,jj) 
    392             zcoef2 = zcoef0 * MIN( e3w_n(ji,jj,iku), e3w_n(ji+1,jj  ,iku) ) 
    393             zcoef3 = zcoef0 * MIN( e3w_n(ji,jj,ikv), e3w_n(ji  ,jj+1,ikv) ) 
     401            zcoef2 = zcoef0 * MIN( e3w(ji,jj,iku,Kmm), e3w(ji+1,jj  ,iku,Kmm) ) 
     402            zcoef3 = zcoef0 * MIN( e3w(ji,jj,ikv,Kmm), e3w(ji  ,jj+1,ikv,Kmm) ) 
    394403            IF( iku > 1 ) THEN            ! on i-direction (level 2 or more) 
    395                ua  (ji,jj,iku) = ua(ji,jj,iku) - zhpi(ji,jj,iku)         ! subtract old value 
     404               puu  (ji,jj,iku,Krhs) = puu(ji,jj,iku,Krhs) - zhpi(ji,jj,iku)         ! subtract old value 
    396405               zhpi(ji,jj,iku) = zhpi(ji,jj,iku-1)                   &   ! compute the new one 
    397406                  &            + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + gru(ji,jj) ) * r1_e1u(ji,jj) 
    398                ua  (ji,jj,iku) = ua(ji,jj,iku) + zhpi(ji,jj,iku)         ! add the new one to the general momentum trend 
     407               puu  (ji,jj,iku,Krhs) = puu(ji,jj,iku,Krhs) + zhpi(ji,jj,iku)         ! add the new one to the general momentum trend 
    399408            ENDIF 
    400409            IF( ikv > 1 ) THEN            ! on j-direction (level 2 or more) 
    401                va  (ji,jj,ikv) = va(ji,jj,ikv) - zhpj(ji,jj,ikv)         ! subtract old value 
     410               pvv  (ji,jj,ikv,Krhs) = pvv(ji,jj,ikv,Krhs) - zhpj(ji,jj,ikv)         ! subtract old value 
    402411               zhpj(ji,jj,ikv) = zhpj(ji,jj,ikv-1)                   &   ! compute the new one 
    403412                  &            + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + grv(ji,jj) ) * r1_e2v(ji,jj) 
    404                va  (ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv)         ! add the new one to the general momentum trend 
     413               pvv  (ji,jj,ikv,Krhs) = pvv(ji,jj,ikv,Krhs) + zhpj(ji,jj,ikv)         ! add the new one to the general momentum trend 
    405414            ENDIF 
    406415         END DO 
     
    410419 
    411420 
    412    SUBROUTINE hpg_sco( kt ) 
     421   SUBROUTINE hpg_sco( kt, Kmm, puu, pvv, Krhs ) 
    413422      !!--------------------------------------------------------------------- 
    414423      !!                  ***  ROUTINE hpg_sco  *** 
     
    422431      !!         zhpi = grav .....  + 1/e1u mi(rhd) di[ grav dep3w ] 
    423432      !!         zhpj = grav .....  + 1/e2v mj(rhd) dj[ grav dep3w ] 
    424       !!      add it to the general momentum trend (ua,va). 
    425       !!         ua = ua - 1/e1u * zhpi 
    426       !!         va = va - 1/e2v * zhpj 
    427       !! 
    428       !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 
    429       !!---------------------------------------------------------------------- 
    430       INTEGER, INTENT(in) ::   kt    ! ocean time-step index 
     433      !!      add it to the general momentum trend (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)). 
     434      !!         puu(:,:,:,Krhs) = puu(:,:,:,Krhs) - 1/e1u * zhpi 
     435      !!         pvv(:,:,:,Krhs) = pvv(:,:,:,Krhs) - 1/e2v * zhpj 
     436      !! 
     437      !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend 
     438      !!---------------------------------------------------------------------- 
     439      INTEGER                             , INTENT( in )  ::  kt          ! ocean time-step index 
     440      INTEGER                             , INTENT( in )  ::  Kmm, Krhs   ! ocean time level indices 
     441      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv    ! ocean velocities and RHS of momentum equation 
    431442      !! 
    432443      INTEGER  ::   ji, jj, jk, jii, jjj                 ! dummy loop indices 
     
    453464        DO jj = 2, jpjm1 
    454465           DO ji = 2, jpim1  
    455              ll_tmp1 = MIN(  sshn(ji,jj)               ,  sshn(ji+1,jj) ) >                & 
     466             ll_tmp1 = MIN(  ssh(ji,jj,Kmm)               ,  ssh(ji+1,jj,Kmm) ) >                & 
    456467                  &    MAX( -ht_0(ji,jj)               , -ht_0(ji+1,jj) ) .AND.            & 
    457                   &    MAX(  sshn(ji,jj) +  ht_0(ji,jj),  sshn(ji+1,jj) + ht_0(ji+1,jj) )  & 
     468                  &    MAX(  ssh(ji,jj,Kmm) +  ht_0(ji,jj),  ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) )  & 
    458469                  &                                                       > rn_wdmin1 + rn_wdmin2 
    459              ll_tmp2 = ( ABS( sshn(ji,jj)              -  sshn(ji+1,jj) ) > 1.E-12 ) .AND. (       & 
    460                   &    MAX(   sshn(ji,jj)              ,  sshn(ji+1,jj) ) >                & 
     470             ll_tmp2 = ( ABS( ssh(ji,jj,Kmm)              -  ssh(ji+1,jj,Kmm) ) > 1.E-12 ) .AND. (       & 
     471                  &    MAX(   ssh(ji,jj,Kmm)              ,  ssh(ji+1,jj,Kmm) ) >                & 
    461472                  &    MAX(  -ht_0(ji,jj)              , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 
    462473 
     
    464475               zcpx(ji,jj) = 1.0_wp 
    465476             ELSE IF(ll_tmp2) THEN 
    466                ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here 
    467                zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) & 
    468                            &    / (sshn(ji+1,jj) - sshn(ji  ,jj)) ) 
     477               ! no worries about  ssh(ji+1,jj,Kmm) -  ssh(ji  ,jj,Kmm) = 0, it won't happen ! here 
     478               zcpx(ji,jj) = ABS( (ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) & 
     479                           &    / (ssh(ji+1,jj,Kmm) - ssh(ji  ,jj,Kmm)) ) 
    469480             ELSE 
    470481               zcpx(ji,jj) = 0._wp 
    471482             END IF 
    472483       
    473              ll_tmp1 = MIN(  sshn(ji,jj)              ,  sshn(ji,jj+1) ) >                & 
     484             ll_tmp1 = MIN(  ssh(ji,jj,Kmm)              ,  ssh(ji,jj+1,Kmm) ) >                & 
    474485                  &    MAX( -ht_0(ji,jj)              , -ht_0(ji,jj+1) ) .AND.            & 
    475                   &    MAX(  sshn(ji,jj) + ht_0(ji,jj),  sshn(ji,jj+1) + ht_0(ji,jj+1) )  & 
     486                  &    MAX(  ssh(ji,jj,Kmm) + ht_0(ji,jj),  ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) )  & 
    476487                  &                                                      > rn_wdmin1 + rn_wdmin2 
    477              ll_tmp2 = ( ABS( sshn(ji,jj)             -  sshn(ji,jj+1) ) > 1.E-12 ) .AND. (        & 
    478                   &    MAX(   sshn(ji,jj)             ,  sshn(ji,jj+1) ) >                & 
     488             ll_tmp2 = ( ABS( ssh(ji,jj,Kmm)             -  ssh(ji,jj+1,Kmm) ) > 1.E-12 ) .AND. (        & 
     489                  &    MAX(   ssh(ji,jj,Kmm)             ,  ssh(ji,jj+1,Kmm) ) >                & 
    479490                  &    MAX(  -ht_0(ji,jj)             , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 
    480491 
     
    482493               zcpy(ji,jj) = 1.0_wp 
    483494             ELSE IF(ll_tmp2) THEN 
    484                ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here 
    485                zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) & 
    486                            &    / (sshn(ji,jj+1) - sshn(ji,jj  )) ) 
     495               ! no worries about  ssh(ji,jj+1,Kmm) -  ssh(ji,jj  ,Kmm) = 0, it won't happen ! here 
     496               zcpy(ji,jj) = ABS( (ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) & 
     497                           &    / (ssh(ji,jj+1,Kmm) - ssh(ji,jj  ,Kmm)) ) 
    487498             ELSE 
    488499               zcpy(ji,jj) = 0._wp 
     
    497508         DO ji = fs_2, fs_jpim1   ! vector opt. 
    498509            ! hydrostatic pressure gradient along s-surfaces 
    499             zhpi(ji,jj,1) = zcoef0 * (  e3w_n(ji+1,jj  ,1) * ( znad + rhd(ji+1,jj  ,1) )    & 
    500                &                      - e3w_n(ji  ,jj  ,1) * ( znad + rhd(ji  ,jj  ,1) )  ) * r1_e1u(ji,jj) 
    501             zhpj(ji,jj,1) = zcoef0 * (  e3w_n(ji  ,jj+1,1) * ( znad + rhd(ji  ,jj+1,1) )    & 
    502                &                      - e3w_n(ji  ,jj  ,1) * ( znad + rhd(ji  ,jj  ,1) )  ) * r1_e2v(ji,jj) 
     510            zhpi(ji,jj,1) = zcoef0 * (  e3w(ji+1,jj  ,1,Kmm) * ( znad + rhd(ji+1,jj  ,1) )    & 
     511               &                      - e3w(ji  ,jj  ,1,Kmm) * ( znad + rhd(ji  ,jj  ,1) )  ) * r1_e1u(ji,jj) 
     512            zhpj(ji,jj,1) = zcoef0 * (  e3w(ji  ,jj+1,1,Kmm) * ( znad + rhd(ji  ,jj+1,1) )    & 
     513               &                      - e3w(ji  ,jj  ,1,Kmm) * ( znad + rhd(ji  ,jj  ,1) )  ) * r1_e2v(ji,jj) 
    503514            ! s-coordinate pressure gradient correction 
    504515            zuap = -zcoef0 * ( rhd    (ji+1,jj,1) + rhd    (ji,jj,1) + 2._wp * znad )   & 
    505                &           * ( gde3w_n(ji+1,jj,1) - gde3w_n(ji,jj,1) ) * r1_e1u(ji,jj) 
     516               &           * ( gde3w(ji+1,jj,1) - gde3w(ji,jj,1) ) * r1_e1u(ji,jj) 
    506517            zvap = -zcoef0 * ( rhd    (ji,jj+1,1) + rhd    (ji,jj,1) + 2._wp * znad )   & 
    507                &           * ( gde3w_n(ji,jj+1,1) - gde3w_n(ji,jj,1) ) * r1_e2v(ji,jj) 
     518               &           * ( gde3w(ji,jj+1,1) - gde3w(ji,jj,1) ) * r1_e2v(ji,jj) 
    508519            ! 
    509520            IF( ln_wd_il ) THEN 
     
    515526            ! 
    516527            ! add to the general momentum trend 
    517             ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap 
    518             va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap 
     528            puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + zhpi(ji,jj,1) + zuap 
     529            pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + zhpj(ji,jj,1) + zvap 
    519530         END DO 
    520531      END DO 
     
    526537               ! hydrostatic pressure gradient along s-surfaces 
    527538               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 * r1_e1u(ji,jj)   & 
    528                   &           * (  e3w_n(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad )   & 
    529                   &              - e3w_n(ji  ,jj,jk) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) + 2*znad )  ) 
     539                  &           * (  e3w(ji+1,jj,jk,Kmm) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad )   & 
     540                  &              - e3w(ji  ,jj,jk,Kmm) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) + 2*znad )  ) 
    530541               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 * r1_e2v(ji,jj)   & 
    531                   &           * (  e3w_n(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad )   & 
    532                   &              - e3w_n(ji,jj  ,jk) * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) + 2*znad )  ) 
     542                  &           * (  e3w(ji,jj+1,jk,Kmm) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad )   & 
     543                  &              - e3w(ji,jj  ,jk,Kmm) * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) + 2*znad )  ) 
    533544               ! s-coordinate pressure gradient correction 
    534545               zuap = -zcoef0 * ( rhd    (ji+1,jj  ,jk) + rhd    (ji,jj,jk) + 2._wp * znad )   & 
    535                   &           * ( gde3w_n(ji+1,jj  ,jk) - gde3w_n(ji,jj,jk) ) * r1_e1u(ji,jj) 
     546                  &           * ( gde3w(ji+1,jj  ,jk) - gde3w(ji,jj,jk) ) * r1_e1u(ji,jj) 
    536547               zvap = -zcoef0 * ( rhd    (ji  ,jj+1,jk) + rhd    (ji,jj,jk) + 2._wp * znad )   & 
    537                   &           * ( gde3w_n(ji  ,jj+1,jk) - gde3w_n(ji,jj,jk) ) * r1_e2v(ji,jj) 
     548                  &           * ( gde3w(ji  ,jj+1,jk) - gde3w(ji,jj,jk) ) * r1_e2v(ji,jj) 
    538549               ! 
    539550               IF( ln_wd_il ) THEN 
     
    545556               ! 
    546557               ! add to the general momentum trend 
    547                ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap 
    548                va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) + zvap 
     558               puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zhpi(ji,jj,jk) + zuap 
     559               pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zhpj(ji,jj,jk) + zvap 
    549560            END DO 
    550561         END DO 
     
    556567 
    557568 
    558    SUBROUTINE hpg_isf( kt ) 
     569   SUBROUTINE hpg_isf( kt, Kmm, puu, pvv, Krhs ) 
    559570      !!--------------------------------------------------------------------- 
    560571      !!                  ***  ROUTINE hpg_isf  *** 
     
    568579      !!         zhpi = grav .....  + 1/e1u mi(rhd) di[ grav dep3w ] 
    569580      !!         zhpj = grav .....  + 1/e2v mj(rhd) dj[ grav dep3w ] 
    570       !!      add it to the general momentum trend (ua,va). 
    571       !!         ua = ua - 1/e1u * zhpi 
    572       !!         va = va - 1/e2v * zhpj 
     581      !!      add it to the general momentum trend (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)). 
     582      !!         puu(:,:,:,Krhs) = puu(:,:,:,Krhs) - 1/e1u * zhpi 
     583      !!         pvv(:,:,:,Krhs) = pvv(:,:,:,Krhs) - 1/e2v * zhpj 
    573584      !!      iceload is added and partial cell case are added to the top and bottom 
    574585      !!       
    575       !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 
    576       !!---------------------------------------------------------------------- 
    577       INTEGER, INTENT(in) ::   kt    ! ocean time-step index 
     586      !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend 
     587      !!---------------------------------------------------------------------- 
     588      INTEGER                             , INTENT( in )  ::  kt          ! ocean time-step index 
     589      INTEGER                             , INTENT( in )  ::  Kmm, Krhs   ! ocean time level indices 
     590      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv    ! ocean velocities and RHS of momentum equation 
    578591      !! 
    579592      INTEGER  ::   ji, jj, jk, ikt, iktp1i, iktp1j   ! dummy loop indices 
     
    596609        DO jj = 1, jpj 
    597610          ikt = mikt(ji,jj) 
    598           zts_top(ji,jj,1) = tsn(ji,jj,ikt,1) 
    599           zts_top(ji,jj,2) = tsn(ji,jj,ikt,2) 
     611          zts_top(ji,jj,1) = ts(ji,jj,ikt,1,Kmm) 
     612          zts_top(ji,jj,2) = ts(ji,jj,ikt,2,Kmm) 
    600613        END DO 
    601614      END DO 
     
    612625            ! hydrostatic pressure gradient along s-surfaces and ice shelf pressure 
    613626            ! we assume ISF is in isostatic equilibrium 
    614             zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( 0.5_wp * e3w_n(ji+1,jj,iktp1i)                                    & 
     627            zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( 0.5_wp * e3w(ji+1,jj,iktp1i,Kmm)                                    & 
    615628               &                                    * ( 2._wp * znad + rhd(ji+1,jj,iktp1i) + zrhdtop_oce(ji+1,jj) )   & 
    616                &                                  - 0.5_wp * e3w_n(ji,jj,ikt)                                         & 
     629               &                                  - 0.5_wp * e3w(ji,jj,ikt,Kmm)                                         & 
    617630               &                                    * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) )          & 
    618631               &                                  + ( riceload(ji+1,jj) - riceload(ji,jj))                            )  
    619             zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * e3w_n(ji,jj+1,iktp1j)                                    & 
     632            zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * e3w(ji,jj+1,iktp1j,Kmm)                                    & 
    620633               &                                    * ( 2._wp * znad + rhd(ji,jj+1,iktp1j) + zrhdtop_oce(ji,jj+1) )   & 
    621                &                                  - 0.5_wp * e3w_n(ji,jj,ikt)                                         &  
     634               &                                  - 0.5_wp * e3w(ji,jj,ikt,Kmm)                                         &  
    622635               &                                    * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) )          & 
    623636               &                                  + ( riceload(ji,jj+1) - riceload(ji,jj))                            )  
    624637            ! s-coordinate pressure gradient correction (=0 if z coordinate) 
    625638            zuap = -zcoef0 * ( rhd    (ji+1,jj,1) + rhd    (ji,jj,1) + 2._wp * znad )   & 
    626                &           * ( gde3w_n(ji+1,jj,1) - gde3w_n(ji,jj,1) ) * r1_e1u(ji,jj) 
     639               &           * ( gde3w(ji+1,jj,1) - gde3w(ji,jj,1) ) * r1_e1u(ji,jj) 
    627640            zvap = -zcoef0 * ( rhd    (ji,jj+1,1) + rhd    (ji,jj,1) + 2._wp * znad )   & 
    628                &           * ( gde3w_n(ji,jj+1,1) - gde3w_n(ji,jj,1) ) * r1_e2v(ji,jj) 
     641               &           * ( gde3w(ji,jj+1,1) - gde3w(ji,jj,1) ) * r1_e2v(ji,jj) 
    629642            ! add to the general momentum trend 
    630             ua(ji,jj,1) = ua(ji,jj,1) + (zhpi(ji,jj,1) + zuap) * umask(ji,jj,1) 
    631             va(ji,jj,1) = va(ji,jj,1) + (zhpj(ji,jj,1) + zvap) * vmask(ji,jj,1) 
     643            puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + (zhpi(ji,jj,1) + zuap) * umask(ji,jj,1) 
     644            pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + (zhpj(ji,jj,1) + zvap) * vmask(ji,jj,1) 
    632645         END DO 
    633646      END DO 
     
    641654               ! hydrostatic pressure gradient along s-surfaces 
    642655               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj)   & 
    643                   &           * (  e3w_n(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) * wmask(ji+1,jj,jk)   & 
    644                   &              - e3w_n(ji  ,jj,jk) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) + 2*znad ) * wmask(ji  ,jj,jk)   ) 
     656                  &           * (  e3w(ji+1,jj,jk,Kmm) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) * wmask(ji+1,jj,jk)   & 
     657                  &              - e3w(ji  ,jj,jk,Kmm) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) + 2*znad ) * wmask(ji  ,jj,jk)   ) 
    645658               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj)   & 
    646                   &           * (  e3w_n(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) * wmask(ji,jj+1,jk)   & 
    647                   &              - e3w_n(ji,jj  ,jk) * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) + 2*znad ) * wmask(ji,jj  ,jk)   ) 
     659                  &           * (  e3w(ji,jj+1,jk,Kmm) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) * wmask(ji,jj+1,jk)   & 
     660                  &              - e3w(ji,jj  ,jk,Kmm) * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) + 2*znad ) * wmask(ji,jj  ,jk)   ) 
    648661               ! s-coordinate pressure gradient correction 
    649662               zuap = -zcoef0 * ( rhd   (ji+1,jj  ,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   & 
    650                   &           * ( gde3w_n(ji+1,jj  ,jk) - gde3w_n(ji,jj,jk) ) / e1u(ji,jj) 
     663                  &           * ( gde3w(ji+1,jj  ,jk) - gde3w(ji,jj,jk) ) / e1u(ji,jj) 
    651664               zvap = -zcoef0 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   & 
    652                   &           * ( gde3w_n(ji  ,jj+1,jk) - gde3w_n(ji,jj,jk) ) / e2v(ji,jj) 
     665                  &           * ( gde3w(ji  ,jj+1,jk) - gde3w(ji,jj,jk) ) / e2v(ji,jj) 
    653666               ! add to the general momentum trend 
    654                ua(ji,jj,jk) = ua(ji,jj,jk) + (zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk) 
    655                va(ji,jj,jk) = va(ji,jj,jk) + (zhpj(ji,jj,jk) + zvap) * vmask(ji,jj,jk) 
     667               puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk) 
     668               pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + (zhpj(ji,jj,jk) + zvap) * vmask(ji,jj,jk) 
    656669            END DO 
    657670         END DO 
     
    661674 
    662675 
    663    SUBROUTINE hpg_djc( kt ) 
     676   SUBROUTINE hpg_djc( kt, Kmm, puu, pvv, Krhs ) 
    664677      !!--------------------------------------------------------------------- 
    665678      !!                  ***  ROUTINE hpg_djc  *** 
     
    669682      !! Reference: Shchepetkin and McWilliams, J. Geophys. Res., 108(C3), 3090, 2003 
    670683      !!---------------------------------------------------------------------- 
    671       INTEGER, INTENT(in) ::   kt    ! ocean time-step index 
     684      INTEGER                             , INTENT( in )  ::  kt          ! ocean time-step index 
     685      INTEGER                             , INTENT( in )  ::  Kmm, Krhs   ! ocean time level indices 
     686      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv    ! ocean velocities and RHS of momentum equation 
    672687      !! 
    673688      INTEGER  ::   ji, jj, jk          ! dummy loop indices 
     
    687702        DO jj = 2, jpjm1 
    688703           DO ji = 2, jpim1  
    689              ll_tmp1 = MIN(  sshn(ji,jj)              ,  sshn(ji+1,jj) ) >                & 
     704             ll_tmp1 = MIN(  ssh(ji,jj,Kmm)              ,  ssh(ji+1,jj,Kmm) ) >                & 
    690705                  &    MAX( -ht_0(ji,jj)              , -ht_0(ji+1,jj) ) .AND.            & 
    691                   &    MAX(  sshn(ji,jj) + ht_0(ji,jj),  sshn(ji+1,jj) + ht_0(ji+1,jj) )  & 
     706                  &    MAX(  ssh(ji,jj,Kmm) + ht_0(ji,jj),  ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) )  & 
    692707                  &                                                      > rn_wdmin1 + rn_wdmin2 
    693              ll_tmp2 = ( ABS( sshn(ji,jj)             -  sshn(ji+1,jj) ) > 1.E-12 ) .AND. (        & 
    694                   &    MAX(   sshn(ji,jj)             ,  sshn(ji+1,jj) ) >                & 
     708             ll_tmp2 = ( ABS( ssh(ji,jj,Kmm)             -  ssh(ji+1,jj,Kmm) ) > 1.E-12 ) .AND. (        & 
     709                  &    MAX(   ssh(ji,jj,Kmm)             ,  ssh(ji+1,jj,Kmm) ) >                & 
    695710                  &    MAX(  -ht_0(ji,jj)             , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 
    696711             IF(ll_tmp1) THEN 
    697712               zcpx(ji,jj) = 1.0_wp 
    698713             ELSE IF(ll_tmp2) THEN 
    699                ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here 
    700                zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) & 
    701                            &    / (sshn(ji+1,jj) - sshn(ji  ,jj)) ) 
     714               ! no worries about  ssh(ji+1,jj,Kmm) -  ssh(ji  ,jj,Kmm) = 0, it won't happen ! here 
     715               zcpx(ji,jj) = ABS( (ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) & 
     716                           &    / (ssh(ji+1,jj,Kmm) - ssh(ji  ,jj,Kmm)) ) 
    702717             ELSE 
    703718               zcpx(ji,jj) = 0._wp 
    704719             END IF 
    705720       
    706              ll_tmp1 = MIN(  sshn(ji,jj)              ,  sshn(ji,jj+1) ) >                & 
     721             ll_tmp1 = MIN(  ssh(ji,jj,Kmm)              ,  ssh(ji,jj+1,Kmm) ) >                & 
    707722                  &    MAX( -ht_0(ji,jj)              , -ht_0(ji,jj+1) ) .AND.            & 
    708                   &    MAX(  sshn(ji,jj) + ht_0(ji,jj),  sshn(ji,jj+1) + ht_0(ji,jj+1) )  & 
     723                  &    MAX(  ssh(ji,jj,Kmm) + ht_0(ji,jj),  ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) )  & 
    709724                  &                                                      > rn_wdmin1 + rn_wdmin2 
    710              ll_tmp2 = ( ABS( sshn(ji,jj)             -  sshn(ji,jj+1) ) > 1.E-12 ) .AND. (        & 
    711                   &    MAX(   sshn(ji,jj)             ,  sshn(ji,jj+1) ) >                & 
     725             ll_tmp2 = ( ABS( ssh(ji,jj,Kmm)             -  ssh(ji,jj+1,Kmm) ) > 1.E-12 ) .AND. (        & 
     726                  &    MAX(   ssh(ji,jj,Kmm)             ,  ssh(ji,jj+1,Kmm) ) >                & 
    712727                  &    MAX(  -ht_0(ji,jj)             , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 
    713728 
     
    715730               zcpy(ji,jj) = 1.0_wp 
    716731             ELSE IF(ll_tmp2) THEN 
    717                ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here 
    718                zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) & 
    719                            &    / (sshn(ji,jj+1) - sshn(ji,jj  )) ) 
     732               ! no worries about  ssh(ji,jj+1,Kmm) -  ssh(ji,jj  ,Kmm) = 0, it won't happen ! here 
     733               zcpy(ji,jj) = ABS( (ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) & 
     734                           &    / (ssh(ji,jj+1,Kmm) - ssh(ji,jj  ,Kmm)) ) 
    720735             ELSE 
    721736               zcpy(ji,jj) = 0._wp 
     
    747762            DO ji = fs_2, fs_jpim1   ! vector opt. 
    748763               drhoz(ji,jj,jk) = rhd    (ji  ,jj  ,jk) - rhd    (ji,jj,jk-1) 
    749                dzz  (ji,jj,jk) = gde3w_n(ji  ,jj  ,jk) - gde3w_n(ji,jj,jk-1) 
     764               dzz  (ji,jj,jk) = gde3w(ji  ,jj  ,jk) - gde3w(ji,jj,jk-1) 
    750765               drhox(ji,jj,jk) = rhd    (ji+1,jj  ,jk) - rhd    (ji,jj,jk  ) 
    751                dzx  (ji,jj,jk) = gde3w_n(ji+1,jj  ,jk) - gde3w_n(ji,jj,jk  ) 
     766               dzx  (ji,jj,jk) = gde3w(ji+1,jj  ,jk) - gde3w(ji,jj,jk  ) 
    752767               drhoy(ji,jj,jk) = rhd    (ji  ,jj+1,jk) - rhd    (ji,jj,jk  ) 
    753                dzy  (ji,jj,jk) = gde3w_n(ji  ,jj+1,jk) - gde3w_n(ji,jj,jk  ) 
     768               dzy  (ji,jj,jk) = gde3w(ji  ,jj+1,jk) - gde3w(ji,jj,jk  ) 
    754769            END DO 
    755770         END DO 
     
    838853      DO jj = 2, jpjm1 
    839854         DO ji = fs_2, fs_jpim1   ! vector opt. 
    840             rho_k(ji,jj,1) = -grav * ( e3w_n(ji,jj,1) - gde3w_n(ji,jj,1) )               & 
     855            rho_k(ji,jj,1) = -grav * ( e3w(ji,jj,1,Kmm) - gde3w(ji,jj,1) )               & 
    841856               &                   * (  rhd(ji,jj,1)                                     & 
    842857               &                     + 0.5_wp * ( rhd    (ji,jj,2) - rhd    (ji,jj,1) )  & 
    843                &                              * ( e3w_n  (ji,jj,1) - gde3w_n(ji,jj,1) )  & 
    844                &                              / ( gde3w_n(ji,jj,2) - gde3w_n(ji,jj,1) )  ) 
     858               &                              * ( e3w  (ji,jj,1,Kmm) - gde3w(ji,jj,1) )  & 
     859               &                              / ( gde3w(ji,jj,2) - gde3w(ji,jj,1) )  ) 
    845860         END DO 
    846861      END DO 
     
    854869 
    855870               rho_k(ji,jj,jk) = zcoef0 * ( rhd    (ji,jj,jk) + rhd    (ji,jj,jk-1) )                                   & 
    856                   &                     * ( gde3w_n(ji,jj,jk) - gde3w_n(ji,jj,jk-1) )                                   & 
     871                  &                     * ( gde3w(ji,jj,jk) - gde3w(ji,jj,jk-1) )                                   & 
    857872                  &            - grav * z1_10 * (                                                                   & 
    858873                  &     ( drhow  (ji,jj,jk) - drhow  (ji,jj,jk-1) )                                                     & 
    859                   &   * ( gde3w_n(ji,jj,jk) - gde3w_n(ji,jj,jk-1) - z1_12 * ( dzw  (ji,jj,jk) + dzw  (ji,jj,jk-1) ) )   & 
     874                  &   * ( gde3w(ji,jj,jk) - gde3w(ji,jj,jk-1) - z1_12 * ( dzw  (ji,jj,jk) + dzw  (ji,jj,jk-1) ) )   & 
    860875                  &   - ( dzw    (ji,jj,jk) - dzw    (ji,jj,jk-1) )                                                     & 
    861876                  &   * ( rhd    (ji,jj,jk) - rhd    (ji,jj,jk-1) - z1_12 * ( drhow(ji,jj,jk) + drhow(ji,jj,jk-1) ) )   & 
     
    863878 
    864879               rho_i(ji,jj,jk) = zcoef0 * ( rhd    (ji+1,jj,jk) + rhd    (ji,jj,jk) )                                   & 
    865                   &                     * ( gde3w_n(ji+1,jj,jk) - gde3w_n(ji,jj,jk) )                                   & 
     880                  &                     * ( gde3w(ji+1,jj,jk) - gde3w(ji,jj,jk) )                                   & 
    866881                  &            - grav* z1_10 * (                                                                    & 
    867882                  &     ( drhou  (ji+1,jj,jk) - drhou  (ji,jj,jk) )                                                     & 
    868                   &   * ( gde3w_n(ji+1,jj,jk) - gde3w_n(ji,jj,jk) - z1_12 * ( dzu  (ji+1,jj,jk) + dzu  (ji,jj,jk) ) )   & 
     883                  &   * ( gde3w(ji+1,jj,jk) - gde3w(ji,jj,jk) - z1_12 * ( dzu  (ji+1,jj,jk) + dzu  (ji,jj,jk) ) )   & 
    869884                  &   - ( dzu    (ji+1,jj,jk) - dzu    (ji,jj,jk) )                                                     & 
    870885                  &   * ( rhd    (ji+1,jj,jk) - rhd    (ji,jj,jk) - z1_12 * ( drhou(ji+1,jj,jk) + drhou(ji,jj,jk) ) )   & 
     
    872887 
    873888               rho_j(ji,jj,jk) = zcoef0 * ( rhd    (ji,jj+1,jk) + rhd    (ji,jj,jk) )                                 & 
    874                   &                     * ( gde3w_n(ji,jj+1,jk) - gde3w_n(ji,jj,jk) )                                   & 
     889                  &                     * ( gde3w(ji,jj+1,jk) - gde3w(ji,jj,jk) )                                   & 
    875890                  &            - grav* z1_10 * (                                                                    & 
    876891                  &     ( drhov  (ji,jj+1,jk) - drhov  (ji,jj,jk) )                                                     & 
    877                   &   * ( gde3w_n(ji,jj+1,jk) - gde3w_n(ji,jj,jk) - z1_12 * ( dzv  (ji,jj+1,jk) + dzv  (ji,jj,jk) ) )   & 
     892                  &   * ( gde3w(ji,jj+1,jk) - gde3w(ji,jj,jk) - z1_12 * ( dzv  (ji,jj+1,jk) + dzv  (ji,jj,jk) ) )   & 
    878893                  &   - ( dzv    (ji,jj+1,jk) - dzv    (ji,jj,jk) )                                                     & 
    879894                  &   * ( rhd    (ji,jj+1,jk) - rhd    (ji,jj,jk) - z1_12 * ( drhov(ji,jj+1,jk) + drhov(ji,jj,jk) ) )   & 
     
    897912            ENDIF 
    898913            ! add to the general momentum trend 
    899             ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) 
    900             va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) 
     914            puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + zhpi(ji,jj,1) 
     915            pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + zhpj(ji,jj,1) 
    901916         END DO 
    902917      END DO 
     
    920935               ENDIF 
    921936               ! add to the general momentum trend 
    922                ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) 
    923                va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) 
     937               puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zhpi(ji,jj,jk) 
     938               pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zhpj(ji,jj,jk) 
    924939            END DO 
    925940         END DO 
     
    931946 
    932947 
    933    SUBROUTINE hpg_prj( kt ) 
     948   SUBROUTINE hpg_prj( kt, Kmm, puu, pvv, Krhs ) 
    934949      !!--------------------------------------------------------------------- 
    935950      !!                  ***  ROUTINE hpg_prj  *** 
     
    940955      !!      all vertical coordinate systems 
    941956      !! 
    942       !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 
     957      !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend 
    943958      !!---------------------------------------------------------------------- 
    944959      INTEGER, PARAMETER  :: polynomial_type = 1    ! 1: cubic spline, 2: linear 
    945       INTEGER, INTENT(in) ::   kt                   ! ocean time-step index 
     960      INTEGER                             , INTENT( in )  ::  kt          ! ocean time-step index 
     961      INTEGER                             , INTENT( in )  ::  Kmm, Krhs   ! ocean time level indices 
     962      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv    ! ocean velocities and RHS of momentum equation 
    946963      !! 
    947964      INTEGER  ::   ji, jj, jk, jkk                 ! dummy loop indices 
     
    975992         DO jj = 2, jpjm1 
    976993           DO ji = 2, jpim1  
    977              ll_tmp1 = MIN(  sshn(ji,jj)              ,  sshn(ji+1,jj) ) >                & 
     994             ll_tmp1 = MIN(  ssh(ji,jj,Kmm)              ,  ssh(ji+1,jj,Kmm) ) >                & 
    978995                  &    MAX( -ht_0(ji,jj)              , -ht_0(ji+1,jj) ) .AND.            & 
    979                   &    MAX(  sshn(ji,jj) + ht_0(ji,jj),  sshn(ji+1,jj) + ht_0(ji+1,jj) )  & 
     996                  &    MAX(  ssh(ji,jj,Kmm) + ht_0(ji,jj),  ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) )  & 
    980997                  &                                                      > rn_wdmin1 + rn_wdmin2 
    981              ll_tmp2 = ( ABS( sshn(ji,jj)             -  sshn(ji+1,jj) ) > 1.E-12 ) .AND. (         & 
    982                   &    MAX(   sshn(ji,jj)             ,  sshn(ji+1,jj) ) >                & 
     998             ll_tmp2 = ( ABS( ssh(ji,jj,Kmm)             -  ssh(ji+1,jj,Kmm) ) > 1.E-12 ) .AND. (         & 
     999                  &    MAX(   ssh(ji,jj,Kmm)             ,  ssh(ji+1,jj,Kmm) ) >                & 
    9831000                  &    MAX(  -ht_0(ji,jj)             , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 
    9841001 
     
    9861003               zcpx(ji,jj) = 1.0_wp 
    9871004             ELSE IF(ll_tmp2) THEN 
    988                ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here 
    989                zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) & 
    990                            &    / (sshn(ji+1,jj) -  sshn(ji  ,jj)) ) 
     1005               ! no worries about  ssh(ji+1,jj,Kmm) -  ssh(ji  ,jj,Kmm) = 0, it won't happen ! here 
     1006               zcpx(ji,jj) = ABS( (ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) & 
     1007                           &    / (ssh(ji+1,jj,Kmm) -  ssh(ji  ,jj,Kmm)) ) 
    9911008               
    9921009                zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp) 
     
    9951012             END IF 
    9961013       
    997              ll_tmp1 = MIN(  sshn(ji,jj)              ,  sshn(ji,jj+1) ) >                & 
     1014             ll_tmp1 = MIN(  ssh(ji,jj,Kmm)              ,  ssh(ji,jj+1,Kmm) ) >                & 
    9981015                  &    MAX( -ht_0(ji,jj)              , -ht_0(ji,jj+1) ) .AND.            & 
    999                   &    MAX(  sshn(ji,jj) + ht_0(ji,jj),  sshn(ji,jj+1) + ht_0(ji,jj+1) )  & 
     1016                  &    MAX(  ssh(ji,jj,Kmm) + ht_0(ji,jj),  ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) )  & 
    10001017                  &                                                      > rn_wdmin1 + rn_wdmin2 
    1001              ll_tmp2 = ( ABS( sshn(ji,jj)             -  sshn(ji,jj+1) ) > 1.E-12 ) .AND. (      & 
    1002                   &    MAX(   sshn(ji,jj)             ,  sshn(ji,jj+1) ) >                & 
     1018             ll_tmp2 = ( ABS( ssh(ji,jj,Kmm)             -  ssh(ji,jj+1,Kmm) ) > 1.E-12 ) .AND. (      & 
     1019                  &    MAX(   ssh(ji,jj,Kmm)             ,  ssh(ji,jj+1,Kmm) ) >                & 
    10031020                  &    MAX(  -ht_0(ji,jj)             , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 
    10041021 
     
    10061023               zcpy(ji,jj) = 1.0_wp 
    10071024             ELSE IF(ll_tmp2) THEN 
    1008                ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here 
    1009                zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) & 
    1010                            &    / (sshn(ji,jj+1) - sshn(ji,jj  )) ) 
     1025               ! no worries about  ssh(ji,jj+1,Kmm) -  ssh(ji,jj  ,Kmm) = 0, it won't happen ! here 
     1026               zcpy(ji,jj) = ABS( (ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) & 
     1027                           &    / (ssh(ji,jj+1,Kmm) - ssh(ji,jj  ,Kmm)) ) 
    10111028                zcpy(ji,jj) = max(min( zcpy(ji,jj) , 1.0_wp),0.0_wp) 
    10121029 
     
    10311048          ELSEIF( jk < jpkm1 ) THEN 
    10321049             DO jkk = jk+1, jpk 
    1033                 zrhh(ji,jj,jkk) = interp1(gde3w_n(ji,jj,jkk  ), gde3w_n(ji,jj,jkk-1),   & 
    1034                    &                      gde3w_n(ji,jj,jkk-2), rhd    (ji,jj,jkk-1), rhd(ji,jj,jkk-2)) 
     1050                zrhh(ji,jj,jkk) = interp1(gde3w(ji,jj,jkk  ), gde3w(ji,jj,jkk-1),   & 
     1051                   &                      gde3w(ji,jj,jkk-2), rhd    (ji,jj,jkk-1), rhd(ji,jj,jkk-2)) 
    10351052             END DO 
    10361053          ENDIF 
     
    10411058      DO jj = 1, jpj 
    10421059         DO ji = 1, jpi 
    1043             zdept(ji,jj,1) = 0.5_wp * e3w_n(ji,jj,1) - sshn(ji,jj) * znad 
     1060            zdept(ji,jj,1) = 0.5_wp * e3w(ji,jj,1,Kmm) - ssh(ji,jj,Kmm) * znad 
    10441061         END DO 
    10451062      END DO 
     
    10481065         DO jj = 1, jpj 
    10491066            DO ji = 1, jpi 
    1050                zdept(ji,jj,jk) = zdept(ji,jj,jk-1) + e3w_n(ji,jj,jk) 
     1067               zdept(ji,jj,jk) = zdept(ji,jj,jk-1) + e3w(ji,jj,jk,Kmm) 
    10511068            END DO 
    10521069         END DO 
     
    10651082        DO ji = 2, jpi 
    10661083          zrhdt1 = zrhh(ji,jj,1) - interp3( zdept(ji,jj,1), asp(ji,jj,1), bsp(ji,jj,1),  & 
    1067              &                                              csp(ji,jj,1), dsp(ji,jj,1) ) * 0.25_wp * e3w_n(ji,jj,1) 
     1084             &                                              csp(ji,jj,1), dsp(ji,jj,1) ) * 0.25_wp * e3w(ji,jj,1,Kmm) 
    10681085 
    10691086          ! assuming linear profile across the top half surface layer 
    1070           zhpi(ji,jj,1) =  0.5_wp * e3w_n(ji,jj,1) * zrhdt1 
     1087          zhpi(ji,jj,1) =  0.5_wp * e3w(ji,jj,1,Kmm) * zrhdt1 
    10711088        END DO 
    10721089      END DO 
     
    10901107        DO ji = 2, jpim1 
    10911108!!gm BUG ?    if it is ssh at u- & v-point then it should be: 
    1092 !          zsshu_n(ji,jj) = (e1e2t(ji,jj) * sshn(ji,jj) + e1e2t(ji+1,jj) * sshn(ji+1,jj)) * & 
     1109!          zsshu_n(ji,jj) = (e1e2t(ji,jj) * ssh(ji,jj,Kmm) + e1e2t(ji+1,jj) * ssh(ji+1,jj,Kmm)) * & 
    10931110!                         & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp  
    1094 !          zsshv_n(ji,jj) = (e1e2t(ji,jj) * sshn(ji,jj) + e1e2t(ji,jj+1) * sshn(ji,jj+1)) * & 
     1111!          zsshv_n(ji,jj) = (e1e2t(ji,jj) * ssh(ji,jj,Kmm) + e1e2t(ji,jj+1) * ssh(ji,jj+1,Kmm)) * & 
    10951112!                         & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp  
    10961113!!gm not this: 
    1097           zsshu_n(ji,jj) = (e1e2u(ji,jj) * sshn(ji,jj) + e1e2u(ji+1, jj) * sshn(ji+1,jj)) * & 
     1114          zsshu_n(ji,jj) = (e1e2u(ji,jj) * ssh(ji,jj,Kmm) + e1e2u(ji+1, jj) * ssh(ji+1,jj,Kmm)) * & 
    10981115                         & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp  
    1099           zsshv_n(ji,jj) = (e1e2v(ji,jj) * sshn(ji,jj) + e1e2v(ji+1, jj) * sshn(ji,jj+1)) * & 
     1116          zsshv_n(ji,jj) = (e1e2v(ji,jj) * ssh(ji,jj,Kmm) + e1e2v(ji+1, jj) * ssh(ji,jj+1,Kmm)) * & 
    11001117                         & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp  
    11011118        END DO 
     
    11061123      DO jj = 2, jpjm1 
    11071124        DO ji = 2, jpim1 
    1108           zu(ji,jj,1) = - ( e3u_n(ji,jj,1) - zsshu_n(ji,jj) * znad)  
    1109           zv(ji,jj,1) = - ( e3v_n(ji,jj,1) - zsshv_n(ji,jj) * znad) 
     1125          zu(ji,jj,1) = - ( e3u(ji,jj,1,Kmm) - zsshu_n(ji,jj) * znad)  
     1126          zv(ji,jj,1) = - ( e3v(ji,jj,1,Kmm) - zsshv_n(ji,jj) * znad) 
    11101127        END DO 
    11111128      END DO 
     
    11141131        DO jj = 2, jpjm1 
    11151132          DO ji = 2, jpim1 
    1116             zu(ji,jj,jk) = zu(ji,jj,jk-1) - e3u_n(ji,jj,jk) 
    1117             zv(ji,jj,jk) = zv(ji,jj,jk-1) - e3v_n(ji,jj,jk) 
     1133            zu(ji,jj,jk) = zu(ji,jj,jk-1) - e3u(ji,jj,jk,Kmm) 
     1134            zv(ji,jj,jk) = zv(ji,jj,jk-1) - e3v(ji,jj,jk,Kmm) 
    11181135          END DO 
    11191136        END DO 
     
    11231140        DO jj = 2, jpjm1 
    11241141          DO ji = 2, jpim1 
    1125             zu(ji,jj,jk) = zu(ji,jj,jk) + 0.5_wp * e3u_n(ji,jj,jk) 
    1126             zv(ji,jj,jk) = zv(ji,jj,jk) + 0.5_wp * e3v_n(ji,jj,jk) 
     1142            zu(ji,jj,jk) = zu(ji,jj,jk) + 0.5_wp * e3u(ji,jj,jk,Kmm) 
     1143            zv(ji,jj,jk) = zv(ji,jj,jk) + 0.5_wp * e3v(ji,jj,jk,Kmm) 
    11271144          END DO 
    11281145        END DO 
     
    11761193               DO WHILE ( -zdept(jid,jj,jk1) < zuijk ) 
    11771194                 IF( jk1 == 1 ) THEN 
    1178                    zdeps = zdept(jid,jj,1) + MIN(zuijk, sshn(jid,jj)*znad) 
     1195                   zdeps = zdept(jid,jj,1) + MIN(zuijk, ssh(jid,jj,Kmm)*znad) 
    11791196                   zrhdt1 = zrhh(jid,jj,1) - interp3(zdept(jid,jj,1), asp(jid,jj,1), & 
    11801197                                                     bsp(jid,jj,1),   csp(jid,jj,1), & 
     
    11961213               IF( .NOT.ln_linssh ) THEN 
    11971214                 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * & 
    1198                     &    ( REAL(jis-jid, wp) * (zpwes + zpwed) + (sshn(ji+1,jj)-sshn(ji,jj)) ) 
     1215                    &    ( REAL(jis-jid, wp) * (zpwes + zpwed) + (ssh(ji+1,jj,Kmm)-ssh(ji,jj,Kmm)) ) 
    11991216                ELSE 
    12001217                 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed) 
     
    12041221                  zdpdx2 = zdpdx2 * zcpx(ji,jj) * wdrampu(ji,jj) 
    12051222               ENDIF 
    1206                ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk)  
     1223               puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk)  
    12071224            ENDIF 
    12081225 
     
    12341251               DO WHILE ( -zdept(ji,jjd,jk1) < zvijk ) 
    12351252                 IF( jk1 == 1 ) THEN 
    1236                    zdeps = zdept(ji,jjd,1) + MIN(zvijk, sshn(ji,jjd)*znad) 
     1253                   zdeps = zdept(ji,jjd,1) + MIN(zvijk, ssh(ji,jjd,Kmm)*znad) 
    12371254                   zrhdt1 = zrhh(ji,jjd,1) - interp3(zdept(ji,jjd,1), asp(ji,jjd,1), & 
    12381255                                                     bsp(ji,jjd,1),   csp(ji,jjd,1), & 
     
    12551272               IF( .NOT.ln_linssh ) THEN 
    12561273                  zdpdy2 = zcoef0 * r1_e2v(ji,jj) * & 
    1257                           ( REAL(jjs-jjd, wp) * (zpnss + zpnsd) + (sshn(ji,jj+1)-sshn(ji,jj)) ) 
     1274                          ( REAL(jjs-jjd, wp) * (zpnss + zpnsd) + (ssh(ji,jj+1,Kmm)-ssh(ji,jj,Kmm)) ) 
    12581275               ELSE 
    12591276                  zdpdy2 = zcoef0 * r1_e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd ) 
     
    12641281               ENDIF 
    12651282 
    1266                va(ji,jj,jk) = va(ji,jj,jk) + (zdpdy1 + zdpdy2) * vmask(ji,jj,jk) 
     1283               pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + (zdpdy1 + zdpdy2) * vmask(ji,jj,jk) 
    12671284            ENDIF 
    12681285               ! 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynkeg.F90

    r10893 r10928  
    187187      ENDIF 
    188188      ! 
    189       IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' keg  - Ua: ', mask1=umask,   & 
    190          &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
     189      IF(ln_ctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' keg  - Ua: ', mask1=umask,   & 
     190         &                       tab3d_2=pvv(:,:,:,Krhs), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    191191      ! 
    192192      IF( ln_timing )   CALL timing_stop('dyn_keg') 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynldf.F90

    r10874 r10928  
    4343CONTAINS 
    4444 
    45    SUBROUTINE dyn_ldf( kt ) 
     45   SUBROUTINE dyn_ldf( kt, Kbb, Kmm, puu, pvv, Krhs ) 
    4646      !!---------------------------------------------------------------------- 
    4747      !!                  ***  ROUTINE dyn_ldf  *** 
     
    4949      !! ** Purpose :   compute the lateral ocean dynamics physics. 
    5050      !!---------------------------------------------------------------------- 
    51       INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     51      INTEGER                             , INTENT( in )  ::  kt               ! ocean time-step index 
     52      INTEGER                             , INTENT( in )  ::  Kbb, Kmm, Krhs   ! ocean time level indices 
     53      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv         ! ocean velocities and RHS of momentum equation 
    5254      ! 
    5355      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdu, ztrdv 
     
    5860      IF( l_trddyn )   THEN                      ! temporary save of momentum trends 
    5961         ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) ) 
    60          ztrdu(:,:,:) = ua(:,:,:)  
    61          ztrdv(:,:,:) = va(:,:,:)  
     62         ztrdu(:,:,:) = puu(:,:,:,Krhs)  
     63         ztrdv(:,:,:) = pvv(:,:,:,Krhs)  
    6264      ENDIF 
    6365 
    6466      SELECT CASE ( nldf_dyn )                   ! compute lateral mixing trend and add it to the general trend 
    6567      ! 
    66       CASE ( np_lap   )    ;   CALL dyn_ldf_lap( kt, ub, vb, ua, va, 1 )      ! iso-level    laplacian 
    67       CASE ( np_lap_i )    ;   CALL dyn_ldf_iso( kt )                         ! rotated      laplacian 
    68       CASE ( np_blp   )    ;   CALL dyn_ldf_blp( kt, ub, vb, ua, va    )      ! iso-level bi-laplacian 
     68      CASE ( np_lap   )   
     69         CALL dyn_ldf_lap( kt, Kbb, Kmm, puu(:,:,:,Kbb), pvv(:,:,:,Kbb), puu(:,:,:,Krhs), pvv(:,:,:,Krhs), 1 ) ! iso-level    laplacian 
     70      CASE ( np_lap_i )  
     71         CALL dyn_ldf_iso( kt, Kbb, Kmm, puu, pvv, Krhs    )                                                   ! rotated      laplacian 
     72      CASE ( np_blp   )   
     73         CALL dyn_ldf_blp( kt, Kbb, Kmm, puu(:,:,:,Kbb), pvv(:,:,:,Kbb), puu(:,:,:,Krhs), pvv(:,:,:,Krhs)    ) ! iso-level bi-laplacian 
    6974      ! 
    7075      END SELECT 
    7176 
    7277      IF( l_trddyn ) THEN                        ! save the horizontal diffusive trends for further diagnostics 
    73          ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 
    74          ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 
     78         ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:) 
     79         ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:) 
    7580         CALL trd_dyn( ztrdu, ztrdv, jpdyn_ldf, kt ) 
    7681         DEALLOCATE ( ztrdu , ztrdv ) 
    7782      ENDIF 
    7883      !                                          ! print sum trends (used for debugging) 
    79       IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' ldf  - Ua: ', mask1=umask,   & 
    80          &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
     84      IF(ln_ctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' ldf  - Ua: ', mask1=umask,   & 
     85         &                       tab3d_2=pvv(:,:,:,Krhs), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    8186      ! 
    8287      IF( ln_timing )   CALL timing_stop('dyn_ldf') 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynldf_iso.F90

    r10425 r10928  
    6060 
    6161 
    62    SUBROUTINE dyn_ldf_iso( kt ) 
     62   SUBROUTINE dyn_ldf_iso( kt, Kbb, Kmm, puu, pvv, Krhs ) 
    6363      !!---------------------------------------------------------------------- 
    6464      !!                     ***  ROUTINE dyn_ldf_iso  *** 
     
    8181      !!      horizontal fluxes associated with the rotated lateral mixing: 
    8282      !!      u-component: 
    83       !!         ziut = ( ahmt + rn_ahm_b ) e2t * e3t / e1t  di[ ub ] 
    84       !!               -  ahmt              e2t * mi-1(uslp) dk[ mi(mk(ub)) ] 
    85       !!         zjuf = ( ahmf + rn_ahm_b ) e1f * e3f / e2f  dj[ ub ] 
    86       !!               -  ahmf              e1f * mi(vslp)   dk[ mj(mk(ub)) ] 
     83      !!         ziut = ( ahmt + rn_ahm_b ) e2t * e3t / e1t  di[ uu ] 
     84      !!               -  ahmt              e2t * mi-1(uslp) dk[ mi(mk(uu)) ] 
     85      !!         zjuf = ( ahmf + rn_ahm_b ) e1f * e3f / e2f  dj[ uu ] 
     86      !!               -  ahmf              e1f * mi(vslp)   dk[ mj(mk(uu)) ] 
    8787      !!      v-component: 
    88       !!         zivf = ( ahmf + rn_ahm_b ) e2t * e3t / e1t  di[ vb ] 
    89       !!               -  ahmf              e2t * mj(uslp)   dk[ mi(mk(vb)) ] 
    90       !!         zjvt = ( ahmt + rn_ahm_b ) e1f * e3f / e2f  dj[ ub ] 
    91       !!               -  ahmt              e1f * mj-1(vslp) dk[ mj(mk(vb)) ] 
     88      !!         zivf = ( ahmf + rn_ahm_b ) e2t * e3t / e1t  di[ vv ] 
     89      !!               -  ahmf              e2t * mj(uslp)   dk[ mi(mk(vv)) ] 
     90      !!         zjvt = ( ahmt + rn_ahm_b ) e1f * e3f / e2f  dj[ vv ] 
     91      !!               -  ahmt              e1f * mj-1(vslp) dk[ mj(mk(vv)) ] 
    9292      !!      take the horizontal divergence of the fluxes: 
    9393      !!         diffu = 1/(e1u*e2u*e3u) {  di  [ ziut ] + dj-1[ zjuf ]  } 
    9494      !!         diffv = 1/(e1v*e2v*e3v) {  di-1[ zivf ] + dj  [ zjvt ]  } 
    95       !!      Add this trend to the general trend (ua,va): 
    96       !!         ua = ua + diffu 
     95      !!      Add this trend to the general trend (uu(rhs),vv(rhs)): 
     96      !!         uu(rhs) = uu(rhs) + diffu 
    9797      !!      CAUTION: here the isopycnal part is with a coeff. of aht. This 
    9898      !!      should be modified for applications others than orca_r2 (!!bug) 
    9999      !! 
    100100      !! ** Action : 
    101       !!       -(ua,va) updated with the before geopotential harmonic mixing trend 
     101      !!       -(puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) updated with the before geopotential harmonic mixing trend 
    102102      !!       -(akzu,akzv) to accompt for the diagonal vertical component 
    103103      !!                    of the rotated operator in dynzdf module 
    104104      !!---------------------------------------------------------------------- 
    105       INTEGER, INTENT( in ) ::   kt   ! ocean time-step index 
     105      INTEGER                             , INTENT( in )  ::  kt               ! ocean time-step index 
     106      INTEGER                             , INTENT( in )  ::  Kbb, Kmm, Krhs   ! ocean time level indices 
     107      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv         ! ocean velocities and RHS of momentum equation 
    106108      ! 
    107109      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     
    128130            DO jj = 2, jpjm1 
    129131               DO ji = 2, jpim1 
    130                   uslp (ji,jj,jk) = - ( gdept_b(ji+1,jj,jk) - gdept_b(ji ,jj ,jk) ) * r1_e1u(ji,jj) * umask(ji,jj,jk) 
    131                   vslp (ji,jj,jk) = - ( gdept_b(ji,jj+1,jk) - gdept_b(ji ,jj ,jk) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk) 
    132                   wslpi(ji,jj,jk) = - ( gdepw_b(ji+1,jj,jk) - gdepw_b(ji-1,jj,jk) ) * r1_e1t(ji,jj) * tmask(ji,jj,jk) * 0.5 
    133                   wslpj(ji,jj,jk) = - ( gdepw_b(ji,jj+1,jk) - gdepw_b(ji,jj-1,jk) ) * r1_e2t(ji,jj) * tmask(ji,jj,jk) * 0.5 
     132                  uslp (ji,jj,jk) = - ( gdept(ji+1,jj,jk,Kbb) - gdept(ji ,jj ,jk,Kbb) ) * r1_e1u(ji,jj) * umask(ji,jj,jk) 
     133                  vslp (ji,jj,jk) = - ( gdept(ji,jj+1,jk,Kbb) - gdept(ji ,jj ,jk,Kbb) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk) 
     134                  wslpi(ji,jj,jk) = - ( gdepw(ji+1,jj,jk,Kbb) - gdepw(ji-1,jj,jk,Kbb) ) * r1_e1t(ji,jj) * tmask(ji,jj,jk) * 0.5 
     135                  wslpj(ji,jj,jk) = - ( gdepw(ji,jj+1,jk,Kbb) - gdepw(ji,jj-1,jk,Kbb) ) * r1_e2t(ji,jj) * tmask(ji,jj,jk) * 0.5 
    134136               END DO 
    135137            END DO 
     
    151153         !                             zdkv(jk=1)=zdkv(jk=2) 
    152154 
    153          zdk1u(:,:) = ( ub(:,:,jk) -ub(:,:,jk+1) ) * umask(:,:,jk+1) 
    154          zdk1v(:,:) = ( vb(:,:,jk) -vb(:,:,jk+1) ) * vmask(:,:,jk+1) 
     155         zdk1u(:,:) = ( puu(:,:,jk,Kbb) -puu(:,:,jk+1,Kbb) ) * umask(:,:,jk+1) 
     156         zdk1v(:,:) = ( pvv(:,:,jk,Kbb) -pvv(:,:,jk+1,Kbb) ) * vmask(:,:,jk+1) 
    155157 
    156158         IF( jk == 1 ) THEN 
     
    158160            zdkv(:,:) = zdk1v(:,:) 
    159161         ELSE 
    160             zdku(:,:) = ( ub(:,:,jk-1) - ub(:,:,jk) ) * umask(:,:,jk) 
    161             zdkv(:,:) = ( vb(:,:,jk-1) - vb(:,:,jk) ) * vmask(:,:,jk) 
     162            zdku(:,:) = ( puu(:,:,jk-1,Kbb) - puu(:,:,jk,Kbb) ) * umask(:,:,jk) 
     163            zdkv(:,:) = ( pvv(:,:,jk-1,Kbb) - pvv(:,:,jk,Kbb) ) * vmask(:,:,jk) 
    162164         ENDIF 
    163165 
     
    171173            DO jj = 2, jpjm1 
    172174               DO ji = fs_2, jpi   ! vector opt. 
    173                   zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj) * MIN( e3u_n(ji,jj,jk), e3u_n(ji-1,jj,jk) ) * r1_e1t(ji,jj) 
     175                  zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj) * MIN( e3u(ji,jj,jk,Kmm), e3u(ji-1,jj,jk,Kmm) ) * r1_e1t(ji,jj) 
    174176 
    175177                  zmskt = 1._wp / MAX(   umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)     & 
     
    178180                  zcof1 = - zaht_0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) 
    179181    
    180                   ziut(ji,jj) = (  zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) )    & 
     182                  ziut(ji,jj) = (  zabe1 * ( puu(ji,jj,jk,Kbb) - puu(ji-1,jj,jk,Kbb) )    & 
    181183                     &           + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)      & 
    182184                     &                      +zdk1u(ji,jj) + zdku (ji-1,jj) )  ) * tmask(ji,jj,jk) 
     
    186188            DO jj = 2, jpjm1 
    187189               DO ji = fs_2, jpi   ! vector opt. 
    188                   zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj) * e3t_n(ji,jj,jk) * r1_e1t(ji,jj) 
     190                  zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj) * e3t(ji,jj,jk,Kmm) * r1_e1t(ji,jj) 
    189191 
    190192                  zmskt = 1._wp / MAX(   umask(ji-1,jj,jk  ) + umask(ji,jj,jk+1)     & 
     
    193195                  zcof1 = - zaht_0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) 
    194196 
    195                   ziut(ji,jj) = (  zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) )   & 
     197                  ziut(ji,jj) = (  zabe1 * ( puu(ji,jj,jk,Kbb) - puu(ji-1,jj,jk,Kbb) )   & 
    196198                     &           + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)     & 
    197199                     &                      +zdk1u(ji,jj) + zdku (ji-1,jj) )  ) * tmask(ji,jj,jk) 
     
    203205         DO jj = 1, jpjm1 
    204206            DO ji = 1, fs_jpim1   ! vector opt. 
    205                zabe2 = ( ahmf(ji,jj,jk) + rn_ahm_b ) * e1f(ji,jj) * e3f_n(ji,jj,jk) * r1_e2f(ji,jj) 
     207               zabe2 = ( ahmf(ji,jj,jk) + rn_ahm_b ) * e1f(ji,jj) * e3f(ji,jj,jk) * r1_e2f(ji,jj) 
    206208 
    207209               zmskf = 1._wp / MAX(   umask(ji,jj+1,jk  )+umask(ji,jj,jk+1)     & 
     
    210212               zcof2 = - zaht_0 * e1f(ji,jj) * zmskf * 0.5  * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) ) 
    211213 
    212                zjuf(ji,jj) = (  zabe2 * ( ub(ji,jj+1,jk) - ub(ji,jj,jk) )   & 
     214               zjuf(ji,jj) = (  zabe2 * ( puu(ji,jj+1,jk,Kbb) - puu(ji,jj,jk,Kbb) )   & 
    213215                  &           + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj)     & 
    214216                  &                      +zdk1u(ji,jj+1) + zdku (ji,jj) )  ) * fmask(ji,jj,jk) 
     
    224226         DO jj = 2, jpjm1 
    225227            DO ji = 1, fs_jpim1   ! vector opt. 
    226                zabe1 = ( ahmf(ji,jj,jk) + rn_ahm_b ) * e2f(ji,jj) * e3f_n(ji,jj,jk) * r1_e1f(ji,jj) 
     228               zabe1 = ( ahmf(ji,jj,jk) + rn_ahm_b ) * e2f(ji,jj) * e3f(ji,jj,jk) * r1_e1f(ji,jj) 
    227229 
    228230               zmskf = 1._wp / MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)     & 
     
    231233               zcof1 = - zaht_0 * e2f(ji,jj) * zmskf * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) ) 
    232234 
    233                zivf(ji,jj) = (  zabe1 * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) )    & 
     235               zivf(ji,jj) = (  zabe1 * ( pvv(ji+1,jj,jk,Kbb) - pvv(ji,jj,jk,Kbb) )    & 
    234236                  &           + zcof1 * (  zdkv (ji,jj) + zdk1v(ji+1,jj)      & 
    235237                  &                      + zdk1v(ji,jj) + zdkv (ji+1,jj) )  ) * fmask(ji,jj,jk) 
     
    241243            DO jj = 2, jpj 
    242244               DO ji = 1, fs_jpim1   ! vector opt. 
    243                   zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj) * MIN( e3v_n(ji,jj,jk), e3v_n(ji,jj-1,jk) ) * r1_e2t(ji,jj) 
     245                  zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj) * MIN( e3v(ji,jj,jk,Kmm), e3v(ji,jj-1,jk,Kmm) ) * r1_e2t(ji,jj) 
    244246 
    245247                  zmskt = 1._wp / MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)     & 
     
    248250                  zcof2 = - zaht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) ) 
    249251 
    250                   zjvt(ji,jj) = (  zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) )    & 
     252                  zjvt(ji,jj) = (  zabe2 * ( pvv(ji,jj,jk,Kbb) - pvv(ji,jj-1,jk,Kbb) )    & 
    251253                     &           + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj)      & 
    252254                     &                      +zdk1v(ji,jj-1) + zdkv (ji,jj) )  ) * tmask(ji,jj,jk) 
     
    256258            DO jj = 2, jpj 
    257259               DO ji = 1, fs_jpim1   ! vector opt. 
    258                   zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj) * e3t_n(ji,jj,jk) * r1_e2t(ji,jj) 
     260                  zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj) * e3t(ji,jj,jk,Kmm) * r1_e2t(ji,jj) 
    259261 
    260262                  zmskt = 1./MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   & 
     
    263265                  zcof2 = - zaht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) ) 
    264266 
    265                   zjvt(ji,jj) = (  zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) )   & 
     267                  zjvt(ji,jj) = (  zabe2 * ( pvv(ji,jj,jk,Kbb) - pvv(ji,jj-1,jk,Kbb) )   & 
    266268                     &           + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj)     & 
    267269                     &                      +zdk1v(ji,jj-1) + zdkv (ji,jj) )  ) * tmask(ji,jj,jk) 
     
    275277         DO jj = 2, jpjm1 
    276278            DO ji = 2, jpim1          !!gm Question vectop possible??? !!bug 
    277                ua(ji,jj,jk) = ua(ji,jj,jk) + (  ziut(ji+1,jj) - ziut(ji,jj  )    & 
    278                   &                           + zjuf(ji  ,jj) - zjuf(ji,jj-1)  ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk) 
    279                va(ji,jj,jk) = va(ji,jj,jk) + (  zivf(ji,jj  ) - zivf(ji-1,jj)    & 
    280                   &                           + zjvt(ji,jj+1) - zjvt(ji,jj  )  ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk) 
     279               puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (  ziut(ji+1,jj) - ziut(ji,jj  )    & 
     280                  &                           + zjuf(ji  ,jj) - zjuf(ji,jj-1)  ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm) 
     281               pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + (  zivf(ji,jj  ) - zivf(ji-1,jj)    & 
     282                  &                           + zjvt(ji,jj+1) - zjvt(ji,jj  )  ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm) 
    281283            END DO 
    282284         END DO 
     
    286288 
    287289      ! print sum trends (used for debugging) 
    288       IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' ldfh - Ua: ', mask1=umask, & 
    289          &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
     290      IF(ln_ctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' ldfh - Ua: ', mask1=umask, & 
     291         &                       tab3d_2=pvv(:,:,:,Krhs), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    290292 
    291293 
     
    306308            DO ji = 2, jpi 
    307309               ! i-gradient of u at jj 
    308                zdiu (ji,jk) = tmask(ji,jj  ,jk) * ( ub(ji,jj  ,jk) - ub(ji-1,jj  ,jk) ) 
     310               zdiu (ji,jk) = tmask(ji,jj  ,jk) * ( puu(ji,jj  ,jk,Kbb) - puu(ji-1,jj  ,jk,Kbb) ) 
    309311               ! j-gradient of u and v at jj 
    310                zdju (ji,jk) = fmask(ji,jj  ,jk) * ( ub(ji,jj+1,jk) - ub(ji  ,jj  ,jk) ) 
    311                zdjv (ji,jk) = tmask(ji,jj  ,jk) * ( vb(ji,jj  ,jk) - vb(ji  ,jj-1,jk) ) 
     312               zdju (ji,jk) = fmask(ji,jj  ,jk) * ( puu(ji,jj+1,jk,Kbb) - puu(ji  ,jj  ,jk,Kbb) ) 
     313               zdjv (ji,jk) = tmask(ji,jj  ,jk) * ( pvv(ji,jj  ,jk,Kbb) - pvv(ji  ,jj-1,jk,Kbb) ) 
    312314               ! j-gradient of u and v at jj+1 
    313                zdj1u(ji,jk) = fmask(ji,jj-1,jk) * ( ub(ji,jj  ,jk) - ub(ji  ,jj-1,jk) ) 
    314                zdj1v(ji,jk) = tmask(ji,jj+1,jk) * ( vb(ji,jj+1,jk) - vb(ji  ,jj  ,jk) ) 
     315               zdj1u(ji,jk) = fmask(ji,jj-1,jk) * ( puu(ji,jj  ,jk,Kbb) - puu(ji  ,jj-1,jk,Kbb) ) 
     316               zdj1v(ji,jk) = tmask(ji,jj+1,jk) * ( pvv(ji,jj+1,jk,Kbb) - pvv(ji  ,jj  ,jk,Kbb) ) 
    315317            END DO 
    316318         END DO 
     
    318320            DO ji = 1, jpim1 
    319321               ! i-gradient of v at jj 
    320                zdiv (ji,jk) = fmask(ji,jj  ,jk) * ( vb(ji+1,jj,jk) - vb(ji  ,jj  ,jk) ) 
     322               zdiv (ji,jk) = fmask(ji,jj  ,jk) * ( pvv(ji+1,jj,jk,Kbb) - pvv(ji  ,jj  ,jk,Kbb) ) 
    321323            END DO 
    322324         END DO 
     
    391393         DO jk = 1, jpkm1 
    392394            DO ji = 2, jpim1 
    393                ua(ji,jj,jk) = ua(ji,jj,jk) + ( zfuw(ji,jk) - zfuw(ji,jk+1) ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk) 
    394                va(ji,jj,jk) = va(ji,jj,jk) + ( zfvw(ji,jk) - zfvw(ji,jk+1) ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk) 
     395               puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + ( zfuw(ji,jk) - zfuw(ji,jk+1) ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm) 
     396               pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + ( zfvw(ji,jk) - zfvw(ji,jk+1) ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm) 
    395397            END DO 
    396398         END DO 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynldf_lap_blp.F90

    r10874 r10928  
    3535CONTAINS 
    3636 
    37    SUBROUTINE dyn_ldf_lap( kt, pub, pvb, pua, pva, kpass ) 
     37   SUBROUTINE dyn_ldf_lap( kt, Kbb, Kmm, pu, pv, pu_rhs, pv_rhs, kpass ) 
    3838      !!---------------------------------------------------------------------- 
    3939      !!                     ***  ROUTINE dyn_ldf_lap  *** 
     
    4545      !!      writen as :   grad_h( ahmt div_h(U )) - curl_h( ahmf curl_z(U) )  
    4646      !! 
    47       !! ** Action : - pua, pva increased by the harmonic operator applied on pub, pvb. 
     47      !! ** Action : - pu_rhs, pv_rhs increased by the harmonic operator applied on pu, pv. 
    4848      !!---------------------------------------------------------------------- 
    4949      INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index 
     50      INTEGER                         , INTENT(in   ) ::   Kbb, Kmm   ! ocean time level indices 
    5051      INTEGER                         , INTENT(in   ) ::   kpass      ! =1/2 first or second passage 
    51       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pub, pvb   ! before velocity  [m/s] 
    52       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva   ! velocity trend   [m/s2] 
     52      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pu, pv     ! before velocity  [m/s] 
     53      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs   ! velocity trend   [m/s2] 
    5354      ! 
    5455      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     
    7677!!gm open question here : e3f  at before or now ?    probably now... 
    7778!!gm note that ahmf has already been multiplied by fmask 
    78                zcur(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) * e3f_n(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1)       & 
    79                   &     * (  e2v(ji  ,jj-1) * pvb(ji  ,jj-1,jk) - e2v(ji-1,jj-1) * pvb(ji-1,jj-1,jk)  & 
    80                   &        - e1u(ji-1,jj  ) * pub(ji-1,jj  ,jk) + e1u(ji-1,jj-1) * pub(ji-1,jj-1,jk)  ) 
     79               zcur(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) * e3f(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1)       & 
     80                  &     * (  e2v(ji  ,jj-1) * pv(ji  ,jj-1,jk) - e2v(ji-1,jj-1) * pv(ji-1,jj-1,jk)  & 
     81                  &        - e1u(ji-1,jj  ) * pu(ji-1,jj  ,jk) + e1u(ji-1,jj-1) * pu(ji-1,jj-1,jk)  ) 
    8182               !                                      ! ahm * div        (computed from 2 to jpi/jpj) 
    8283!!gm note that ahmt has already been multiplied by tmask 
    83                zdiv(ji,jj)     = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t_b(ji,jj,jk)                                         & 
    84                   &     * (  e2u(ji,jj)*e3u_b(ji,jj,jk) * pub(ji,jj,jk) - e2u(ji-1,jj)*e3u_b(ji-1,jj,jk) * pub(ji-1,jj,jk)  & 
    85                   &        + e1v(ji,jj)*e3v_b(ji,jj,jk) * pvb(ji,jj,jk) - e1v(ji,jj-1)*e3v_b(ji,jj-1,jk) * pvb(ji,jj-1,jk)  ) 
     84               zdiv(ji,jj)     = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kbb)                                         & 
     85                  &     * (  e2u(ji,jj)*e3u(ji,jj,jk,Kbb) * pu(ji,jj,jk) - e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kbb) * pu(ji-1,jj,jk)  & 
     86                  &        + e1v(ji,jj)*e3v(ji,jj,jk,Kbb) * pv(ji,jj,jk) - e1v(ji,jj-1)*e3v(ji,jj-1,jk,Kbb) * pv(ji,jj-1,jk)  ) 
    8687            END DO   
    8788         END DO   
     
    8990         DO jj = 2, jpjm1                             ! - curl( curl) + grad( div ) 
    9091            DO ji = fs_2, fs_jpim1   ! vector opt. 
    91                pua(ji,jj,jk) = pua(ji,jj,jk) + zsign * (                                                 & 
    92                   &              - ( zcur(ji  ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj) / e3u_n(ji,jj,jk)   & 
     92               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * (                                                 & 
     93                  &              - ( zcur(ji  ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj) / e3u(ji,jj,jk,Kmm)   & 
    9394                  &              + ( zdiv(ji+1,jj) - zdiv(ji,jj  ) ) * r1_e1u(ji,jj)                     ) 
    9495                  ! 
    95                pva(ji,jj,jk) = pva(ji,jj,jk) + zsign * (                                                 & 
    96                   &                ( zcur(ji,jj  ) - zcur(ji-1,jj) ) * r1_e1v(ji,jj) / e3v_n(ji,jj,jk)   & 
     96               pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zsign * (                                                 & 
     97                  &                ( zcur(ji,jj  ) - zcur(ji-1,jj) ) * r1_e1v(ji,jj) / e3v(ji,jj,jk,Kmm)   & 
    9798                  &              + ( zdiv(ji,jj+1) - zdiv(ji  ,jj) ) * r1_e2v(ji,jj)                     ) 
    9899            END DO 
     
    105106 
    106107 
    107    SUBROUTINE dyn_ldf_blp( kt, pub, pvb, pua, pva ) 
     108   SUBROUTINE dyn_ldf_blp( kt, Kbb, Kmm, pu, pv, pu_rhs, pv_rhs ) 
    108109      !!---------------------------------------------------------------------- 
    109110      !!                 ***  ROUTINE dyn_ldf_blp  *** 
     
    116117      !!      It is computed by two successive calls to dyn_ldf_lap routine 
    117118      !! 
    118       !! ** Action :   pta   updated with the before rotated bilaplacian diffusion 
     119      !! ** Action :   pt(:,:,:,:,Krhs)   updated with the before rotated bilaplacian diffusion 
    119120      !!---------------------------------------------------------------------- 
    120121      INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index 
    121       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pub, pvb   ! before velocity fields 
    122       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva   ! momentum trend 
     122      INTEGER                         , INTENT(in   ) ::   Kbb, Kmm   ! ocean time level indices 
     123      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pu, pv     ! before velocity fields 
     124      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs   ! momentum trend 
    123125      ! 
    124126      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zulap, zvlap   ! laplacian at u- and v-point 
     
    134136      zvlap(:,:,:) = 0._wp 
    135137      ! 
    136       CALL dyn_ldf_lap( kt, pub, pvb, zulap, zvlap, 1 )   ! rotated laplacian applied to ptb (output in zlap) 
     138      CALL dyn_ldf_lap( kt, Kbb, Kmm, pu, pv, zulap, zvlap, 1 )   ! rotated laplacian applied to pt (output in zlap,Kbb) 
    137139      ! 
    138140      CALL lbc_lnk_multi( 'dynldf_lap_blp', zulap, 'U', -1., zvlap, 'V', -1. )             ! Lateral boundary conditions 
    139141      ! 
    140       CALL dyn_ldf_lap( kt, zulap, zvlap, pua, pva, 2 )   ! rotated laplacian applied to zlap (output in pta) 
     142      CALL dyn_ldf_lap( kt, Kbb, Kmm, zulap, zvlap, pu_rhs, pv_rhs, 2 )   ! rotated laplacian applied to zlap (output in pt(:,:,:,:,Krhs)) 
    141143      ! 
    142144   END SUBROUTINE dyn_ldf_blp 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynspg.F90

    r10919 r10928  
    175175      ENDIF 
    176176      !                                      ! print mean trends (used for debugging) 
    177       IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' spg  - Ua: ', mask1=umask, & 
    178          &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
     177      IF(ln_ctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' spg  - Ua: ', mask1=umask, & 
     178         &                       tab3d_2=pvv(:,:,:,Krhs), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    179179      ! 
    180180      IF( ln_timing )   CALL timing_stop('dyn_spg') 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynvor.F90

    r10893 r10928  
    181181      ! 
    182182      !                       ! print sum trends (used for debugging) 
    183       IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' vor  - Ua: ', mask1=umask,               & 
    184          &                     tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
     183      IF(ln_ctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' vor  - Ua: ', mask1=umask,               & 
     184         &                     tab3d_2=pvv(:,:,:,Krhs), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    185185      ! 
    186186      IF( ln_timing )   CALL timing_stop('dyn_vor') 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynzad.F90

    r10893 r10928  
    116116      ENDIF 
    117117      !                             ! Control print 
    118       IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' zad  - Ua: ', mask1=umask,   & 
    119          &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
     118      IF(ln_ctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' zad  - Ua: ', mask1=umask,   & 
     119         &                       tab3d_2=pvv(:,:,:,Krhs), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    120120      ! 
    121121      IF( ln_timing )   CALL timing_stop('dyn_zad') 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynzdf.F90

    r10893 r10928  
    130130      IF( ln_drgimp .AND. ln_dynspg_ts ) THEN 
    131131         DO jk = 1, jpkm1        ! remove barotropic velocities 
    132             puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kaa) - ua_b(:,:) ) * umask(:,:,jk) 
    133             pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kaa) - va_b(:,:) ) * vmask(:,:,jk) 
     132            puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kaa) - uu_b(:,:,Kaa) ) * umask(:,:,jk) 
     133            pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kaa) - vv_b(:,:,Kaa) ) * vmask(:,:,jk) 
    134134         END DO 
    135135         DO jj = 2, jpjm1        ! Add bottom/top stress due to barotropic component only 
     
    139139               ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa) 
    140140               ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa) 
    141                puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * ua_b(ji,jj) / ze3ua 
    142                pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * va_b(ji,jj) / ze3va 
     141               puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * uu_b(ji,jj,Kaa) / ze3ua 
     142               pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * vv_b(ji,jj,Kaa) / ze3va 
    143143            END DO 
    144144         END DO 
     
    150150                  ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa) 
    151151                  ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa) 
    152                   puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * ua_b(ji,jj) / ze3ua 
    153                   pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * va_b(ji,jj) / ze3va 
     152                  puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * uu_b(ji,jj,Kaa) / ze3ua 
     153                  pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * vv_b(ji,jj,Kaa) / ze3va 
    154154               END DO 
    155155            END DO 
     
    494494      ENDIF 
    495495      !                                          ! print mean trends (used for debugging) 
    496       IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' zdf  - Ua: ', mask1=umask,               & 
    497          &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
     496      IF(ln_ctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Kaa), clinfo1=' zdf  - Ua: ', mask1=umask,               & 
     497         &                       tab3d_2=pvv(:,:,:,Kaa), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    498498         ! 
    499499      IF( ln_timing )   CALL timing_stop('dyn_zdf') 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/nemogcm.F90

    r10927 r10928  
    428428      !                             
    429429      IF( ln_diurnal_only ) THEN                   ! diurnal only: a subset of the initialisation routines 
    430          CALL  istate_init( Nbb, Nnn )                ! ocean initial state (Dynamics and tracers) 
     430         CALL  istate_init( Nbb, Nnn, Naa )           ! ocean initial state (Dynamics and tracers) 
    431431         CALL     sbc_init( Nbb, Nnn )                ! Forcings : surface module 
    432432         CALL tra_qsr_init                            ! penetrative solar radiation qsr 
     
    440440      ENDIF 
    441441       
    442                            CALL  istate_init( Nbb, Nnn ) ! ocean initial state (Dynamics and tracers) 
     442                           CALL  istate_init( Nbb, Nnn, Naa )    ! ocean initial state (Dynamics and tracers) 
    443443 
    444444      !                                      ! external forcing  
     
    464464 
    465465      !                                      ! Dynamics 
    466       IF( lk_c1d       )   CALL dyn_dmp_init      ! internal momentum damping 
    467                            CALL dyn_adv_init      ! advection (vector or flux form) 
    468                            CALL dyn_vor_init      ! vorticity term including Coriolis 
    469                            CALL dyn_ldf_init      ! lateral mixing 
    470                            CALL dyn_hpg_init      ! horizontal gradient of Hydrostatic pressure 
    471                            CALL dyn_spg_init      ! surface pressure gradient 
     466      IF( lk_c1d       )   CALL dyn_dmp_init         ! internal momentum damping 
     467                           CALL dyn_adv_init         ! advection (vector or flux form) 
     468                           CALL dyn_vor_init         ! vorticity term including Coriolis 
     469                           CALL dyn_ldf_init         ! lateral mixing 
     470                           CALL dyn_hpg_init( Nnn )  ! horizontal gradient of Hydrostatic pressure 
     471                           CALL dyn_spg_init         ! surface pressure gradient 
    472472 
    473473#if defined key_top 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/step.F90

    r10927 r10928  
    195195                         CALL dyn_adv( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! advection (VF or FF)   ==> RHS 
    196196                         CALL dyn_vor( kstp,      Nnn      , uu, vv, Nrhs )  ! vorticity              ==> RHS 
    197                          CALL dyn_ldf       ( kstp )  ! lateral mixing 
     197                         CALL dyn_ldf( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! lateral mixing 
    198198      IF( ln_zdfosm  )   CALL dyn_osm( kstp,      Nnn      , uu, vv, Nrhs )  ! OSMOSIS non-local velocity fluxes ==> RHS 
    199                          CALL dyn_hpg       ( kstp )  ! horizontal gradient of Hydrostatic pressure 
    200                          CALL dyn_spg       ( kstp, Nbb, Nnn, Nrhs, uu, vv, ssh, uu_b, vv_b, Naa )  ! surface pressure gradient 
     199                         CALL dyn_hpg( kstp,      Nnn,       uu, vv, Nrhs )  ! horizontal gradient of Hydrostatic pressure 
     200                         CALL dyn_spg( kstp, Nbb, Nnn, Nrhs, uu, vv, ssh, uu_b, vv_b, Naa )  ! surface pressure gradient 
    201201 
    202202                                                      ! With split-explicit free surface, since now transports have been updated and ssha as well 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OFF/nemogcm.F90

    r10922 r10928  
    315315      IF( ln_ctl       )   CALL prt_ctl_init        ! Print control 
    316316 
    317                            CALL  istate_init    ! ocean initial state (Dynamics and tracers) 
     317                           CALL  istate_init( Nnn, Naa )    ! ocean initial state (Dynamics and tracers) 
    318318 
    319319                           CALL     sbc_init( Nbb, Nnn )    ! Forcings : surface module 
     
    524524   END SUBROUTINE nemo_set_cfctl 
    525525 
    526    SUBROUTINE istate_init 
     526   SUBROUTINE istate_init( Kmm, Kaa ) 
    527527      !!---------------------------------------------------------------------- 
    528528      !!                   ***  ROUTINE istate_init  *** 
     
    530530      !! ** Purpose :   Initialization to zero of the dynamics and tracers. 
    531531      !!---------------------------------------------------------------------- 
     532      INTEGER, INTENT(in) ::   Kmm, Kaa  ! ocean time level indices 
    532533      ! 
    533534      !     now fields         !     after fields      ! 
    534       un   (:,:,:)   = 0._wp   ;   ua(:,:,:) = 0._wp   ! 
    535       vn   (:,:,:)   = 0._wp   ;   va(:,:,:) = 0._wp   ! 
    536       wn   (:,:,:)   = 0._wp   !                       ! 
     535      uu   (:,:,:,Kmm)   = 0._wp   ;   uu(:,:,:,Kaa) = 0._wp   ! 
     536      vv   (:,:,:,Kmm)   = 0._wp   ;   vv(:,:,:,Kaa) = 0._wp   ! 
     537      ww   (:,:,:)   = 0._wp   !                       ! 
    537538      hdivn(:,:,:)   = 0._wp   !                       ! 
    538       tsn  (:,:,:,:) = 0._wp   !                       ! 
     539      ts  (:,:,:,:,Kmm) = 0._wp   !                       ! 
    539540      ! 
    540541      rhd  (:,:,:) = 0.e0 
Note: See TracChangeset for help on using the changeset viewer.