Changeset 10978 for NEMO


Ignore:
Timestamp:
2019-05-15T09:41:30+02:00 (19 months ago)
Author:
davestorkey
Message:

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : DOM and sshwzv.F90. (sshwzv.F90 will be rewritten somewhat when we rewrite the time-level swapping but I've done it anyway). Passes all non-AGRIF SETTE tests.

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

Legend:

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

    r10969 r10978  
    7979                         CALL zdf_phy( kstp )         ! vertical physics update (bfr, avt, avs, avm + MLD) 
    8080 
    81       IF(.NOT.ln_linssh )   CALL ssh_nxt       ( kstp, Nbb, Nnn )  ! after ssh (includes call to div_hor) 
    82       IF(.NOT.ln_linssh )   CALL dom_vvl_sf_nxt( kstp )  ! after vertical scale factors  
     81      IF(.NOT.ln_linssh )   CALL ssh_nxt       ( kstp, Nbb, Nnn, ssh, Naa )  ! after ssh (includes call to div_hor) 
     82      IF(.NOT.ln_linssh )   CALL dom_vvl_sf_nxt( kstp, Nbb, Nnn,      Naa )  ! after vertical scale factors  
    8383 
    84       IF(.NOT.ln_linssh )   CALL wzv           ( kstp )  ! now cross-level velocity  
     84      IF(.NOT.ln_linssh )   CALL wzv           ( kstp, Nbb, Nnn, ww,  Naa )  ! now cross-level velocity  
    8585      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    8686      ! diagnostics and outputs             (ua, va, ta, sa used as workspace) 
     
    125125                        CALL dyn_zdf    ( kstp )   ! vertical diffusion 
    126126                        CALL dyn_nxt    ( kstp )   ! lateral velocity at next time step 
    127       IF(.NOT.ln_linssh)CALL ssh_swp    ( kstp )   ! swap of sea surface height 
     127      IF(.NOT.ln_linssh)CALL ssh_swp    ( kstp, Nbb, Nnn, Naa )   ! swap of sea surface height 
    128128 
    129       IF(.NOT.ln_linssh)CALL dom_vvl_sf_swp( kstp )! swap of vertical scale factors 
     129      IF(.NOT.ln_linssh)CALL dom_vvl_sf_swp( kstp, Nbb, Nnn, Naa )! swap of vertical scale factors 
    130130 
    131131      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM/domain.F90

    r10425 r10978  
    5858CONTAINS 
    5959 
    60    SUBROUTINE dom_init(cdstr) 
     60   SUBROUTINE dom_init( Kbb, Kmm, Kaa, cdstr ) 
    6161      !!---------------------------------------------------------------------- 
    6262      !!                  ***  ROUTINE dom_init  *** 
     
    7373      !!              - 1D configuration, move Coriolis, u and v at T-point 
    7474      !!---------------------------------------------------------------------- 
     75      INTEGER          , INTENT(in) :: Kbb, Kmm, Kaa          ! ocean time level indices 
     76      CHARACTER (len=*), INTENT(in) :: cdstr                  ! model: NEMO or SAS. Determines core restart variables 
     77      ! 
    7578      INTEGER ::   ji, jj, jk, ik   ! dummy loop indices 
    7679      INTEGER ::   iconf = 0    ! local integers 
    7780      CHARACTER (len=64) ::   cform = "(A12, 3(A13, I7))"  
    78       CHARACTER (len=*), INTENT(IN) :: cdstr                  ! model: NEMO or SAS. Determines core restart variables 
    7981      INTEGER , DIMENSION(jpi,jpj) ::   ik_top , ik_bot       ! top and bottom ocean level 
    8082      REAL(wp), DIMENSION(jpi,jpj) ::   z1_hu_0, z1_hv_0 
     
    161163      ! 
    162164         !       before        !          now          !       after         ! 
    163             gdept_b = gdept_0  ;   gdept_n = gdept_0   !        ---          ! depth of grid-points 
    164             gdepw_b = gdepw_0  ;   gdepw_n = gdepw_0   !        ---          ! 
    165                                    gde3w_n = gde3w_0   !        ---          ! 
     165            gdept(:,:,:,Kbb) = gdept_0  ;   gdept(:,:,:,Kmm) = gdept_0   !        ---          ! depth of grid-points 
     166            gdepw(:,:,:,Kbb) = gdepw_0  ;   gdepw(:,:,:,Kmm) = gdepw_0   !        ---          ! 
     167                                   gde3w = gde3w_0   !        ---          ! 
    166168         !                                                                   
    167               e3t_b =   e3t_0  ;     e3t_n =   e3t_0   ;   e3t_a =  e3t_0    ! scale factors 
    168               e3u_b =   e3u_0  ;     e3u_n =   e3u_0   ;   e3u_a =  e3u_0    ! 
    169               e3v_b =   e3v_0  ;     e3v_n =   e3v_0   ;   e3v_a =  e3v_0    ! 
    170                                      e3f_n =   e3f_0   !        ---          ! 
    171               e3w_b =   e3w_0  ;     e3w_n =   e3w_0   !        ---          ! 
    172              e3uw_b =  e3uw_0  ;    e3uw_n =  e3uw_0   !        ---          ! 
    173              e3vw_b =  e3vw_0  ;    e3vw_n =  e3vw_0   !        ---          ! 
     169              e3t(:,:,:,Kbb) =   e3t_0  ;     e3t(:,:,:,Kmm) =   e3t_0   ;   e3t(:,:,:,Kaa) =  e3t_0    ! scale factors 
     170              e3u(:,:,:,Kbb) =   e3u_0  ;     e3u(:,:,:,Kmm) =   e3u_0   ;   e3u(:,:,:,Kaa) =  e3u_0    ! 
     171              e3v(:,:,:,Kbb) =   e3v_0  ;     e3v(:,:,:,Kmm) =   e3v_0   ;   e3v(:,:,:,Kaa) =  e3v_0    ! 
     172                                     e3f =   e3f_0   !        ---          ! 
     173              e3w(:,:,:,Kbb) =   e3w_0  ;     e3w(:,:,:,Kmm) =   e3w_0   !        ---          ! 
     174             e3uw(:,:,:,Kbb) =  e3uw_0  ;    e3uw(:,:,:,Kmm) =  e3uw_0   !        ---          ! 
     175             e3vw(:,:,:,Kbb) =  e3vw_0  ;    e3vw(:,:,:,Kmm) =  e3vw_0   !        ---          ! 
    174176         ! 
    175177         z1_hu_0(:,:) = ssumask(:,:) / ( hu_0(:,:) + 1._wp - ssumask(:,:) )     ! _i mask due to ISF 
     
    186188      ELSE                       != time varying : initialize before/now/after variables 
    187189         ! 
    188          IF( .NOT.l_offline )  CALL dom_vvl_init  
     190         IF( .NOT.l_offline )  CALL dom_vvl_init( Kbb, Kmm, Kaa ) 
    189191         ! 
    190192      ENDIF 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM/domvvl.F90

    r10799 r10978  
    9393 
    9494 
    95    SUBROUTINE dom_vvl_init 
     95   SUBROUTINE dom_vvl_init( Kbb, Kmm, Kaa ) 
    9696      !!---------------------------------------------------------------------- 
    9797      !!                ***  ROUTINE dom_vvl_init  *** 
     
    104104      !! 
    105105      !! ** Action  : - e3t_(n/b) and tilde_e3t_(n/b) 
    106       !!              - Regrid: e3(u/v)_n 
    107       !!                        e3(u/v)_b        
    108       !!                        e3w_n            
    109       !!                        e3(u/v)w_b       
    110       !!                        e3(u/v)w_n       
    111       !!                        gdept_n, gdepw_n and gde3w_n 
     106      !!              - Regrid: e3[u/v](:,:,:,Kmm) 
     107      !!                        e3[u/v](:,:,:,Kmm)        
     108      !!                        e3w(:,:,:,Kmm)            
     109      !!                        e3[u/v]w_b 
     110      !!                        e3[u/v]w_n       
     111      !!                        gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm) and gde3w 
    112112      !!              - h(t/u/v)_0 
    113113      !!              - frq_rst_e3t and frq_rst_hdv 
     
    115115      !! Reference  : Leclair, M., and G. Madec, 2011, Ocean Modelling. 
    116116      !!---------------------------------------------------------------------- 
     117      INTEGER, INTENT(in) :: Kbb, Kmm, Kaa 
     118      ! 
    117119      INTEGER ::   ji, jj, jk 
    118120      INTEGER ::   ii0, ii1, ij0, ij1 
     
    130132      ! 
    131133      !                    ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf 
    132       CALL dom_vvl_rst( nit000, 'READ' ) 
    133       e3t_a(:,:,jpk) = e3t_0(:,:,jpk)  ! last level always inside the sea floor set one for all 
     134      CALL dom_vvl_rst( nit000, Kbb, Kmm, 'READ' ) 
     135      e3t(:,:,jpk,Kaa) = e3t_0(:,:,jpk)  ! last level always inside the sea floor set one for all 
    134136      ! 
    135137      !                    !== Set of all other vertical scale factors  ==!  (now and before) 
    136138      !                                ! Horizontal interpolation of e3t 
    137       CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )    ! from T to U 
    138       CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' ) 
    139       CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )    ! from T to V  
    140       CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' ) 
    141       CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' )    ! from U to F 
     139      CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3u(:,:,:,Kbb), 'U' )    ! from T to U 
     140      CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' ) 
     141      CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3v(:,:,:,Kbb), 'V' )    ! from T to V  
     142      CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' ) 
     143      CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' )    ! from U to F 
    142144      !                                ! Vertical interpolation of e3t,u,v  
    143       CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W'  )  ! from T to W 
    144       CALL dom_vvl_interpol( e3t_b(:,:,:), e3w_b (:,:,:), 'W'  ) 
    145       CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )  ! from U to UW 
    146       CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' ) 
    147       CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )  ! from V to UW 
    148       CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' ) 
     145      CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W'  )  ! from T to W 
     146      CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3w (:,:,:,Kbb), 'W'  ) 
     147      CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' )  ! from U to UW 
     148      CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' ) 
     149      CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' )  ! from V to UW 
     150      CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' ) 
    149151 
    150152      ! We need to define e3[tuv]_a for AGRIF initialisation (should not be a problem for the restartability...) 
    151       e3t_a(:,:,:) = e3t_n(:,:,:) 
    152       e3u_a(:,:,:) = e3u_n(:,:,:) 
    153       e3v_a(:,:,:) = e3v_n(:,:,:) 
     153      e3t(:,:,:,Kaa) = e3t(:,:,:,Kmm) 
     154      e3u(:,:,:,Kaa) = e3u(:,:,:,Kmm) 
     155      e3v(:,:,:,Kaa) = e3v(:,:,:,Kmm) 
    154156      ! 
    155157      !                    !==  depth of t and w-point  ==!   (set the isf depth as it is in the initial timestep) 
    156       gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)       ! reference to the ocean surface (used for MLD and light penetration) 
    157       gdepw_n(:,:,1) = 0.0_wp 
    158       gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)  ! reference to a common level z=0 for hpg 
    159       gdept_b(:,:,1) = 0.5_wp * e3w_b(:,:,1) 
    160       gdepw_b(:,:,1) = 0.0_wp 
     158      gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm)       ! reference to the ocean surface (used for MLD and light penetration) 
     159      gdepw(:,:,1,Kmm) = 0.0_wp 
     160      gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm)  ! reference to a common level z=0 for hpg 
     161      gdept(:,:,1,Kbb) = 0.5_wp * e3w(:,:,1,Kbb) 
     162      gdepw(:,:,1,Kbb) = 0.0_wp 
    161163      DO jk = 2, jpk                               ! vertical sum 
    162164         DO jj = 1,jpj 
     
    165167               !                             ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) 
    166168               !                             ! 0.5 where jk = mikt      
    167 !!gm ???????   BUG ?  gdept_n as well as gde3w_n  does not include the thickness of ISF ?? 
     169!!gm ???????   BUG ?  gdept(:,:,:,Kmm) as well as gde3w  does not include the thickness of ISF ?? 
    168170               zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) ) 
    169                gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1) 
    170                gdept_n(ji,jj,jk) =      zcoef  * ( gdepw_n(ji,jj,jk  ) + 0.5 * e3w_n(ji,jj,jk))  & 
    171                   &                + (1-zcoef) * ( gdept_n(ji,jj,jk-1) +       e3w_n(ji,jj,jk))  
    172                gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj) 
    173                gdepw_b(ji,jj,jk) = gdepw_b(ji,jj,jk-1) + e3t_b(ji,jj,jk-1) 
    174                gdept_b(ji,jj,jk) =      zcoef  * ( gdepw_b(ji,jj,jk  ) + 0.5 * e3w_b(ji,jj,jk))  & 
    175                   &                + (1-zcoef) * ( gdept_b(ji,jj,jk-1) +       e3w_b(ji,jj,jk))  
     171               gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 
     172               gdept(ji,jj,jk,Kmm) =      zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm))  & 
     173                  &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm))  
     174               gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm) 
     175               gdepw(ji,jj,jk,Kbb) = gdepw(ji,jj,jk-1,Kbb) + e3t(ji,jj,jk-1,Kbb) 
     176               gdept(ji,jj,jk,Kbb) =      zcoef  * ( gdepw(ji,jj,jk  ,Kbb) + 0.5 * e3w(ji,jj,jk,Kbb))  & 
     177                  &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kbb) +       e3w(ji,jj,jk,Kbb))  
    176178            END DO 
    177179         END DO 
     
    179181      ! 
    180182      !                    !==  thickness of the water column  !!   (ocean portion only) 
    181       ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1)   !!gm  BUG  :  this should be 1/2 * e3w(k=1) .... 
    182       hu_b(:,:) = e3u_b(:,:,1) * umask(:,:,1) 
    183       hu_n(:,:) = e3u_n(:,:,1) * umask(:,:,1) 
    184       hv_b(:,:) = e3v_b(:,:,1) * vmask(:,:,1) 
    185       hv_n(:,:) = e3v_n(:,:,1) * vmask(:,:,1) 
     183      ht_n(:,:) = e3t(:,:,1,Kmm) * tmask(:,:,1)   !!gm  BUG  :  this should be 1/2 * e3w(k=1) .... 
     184      hu_b(:,:) = e3u(:,:,1,Kbb) * umask(:,:,1) 
     185      hu_n(:,:) = e3u(:,:,1,Kmm) * umask(:,:,1) 
     186      hv_b(:,:) = e3v(:,:,1,Kbb) * vmask(:,:,1) 
     187      hv_n(:,:) = e3v(:,:,1,Kmm) * vmask(:,:,1) 
    186188      DO jk = 2, jpkm1 
    187          ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 
    188          hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk) 
    189          hu_n(:,:) = hu_n(:,:) + e3u_n(:,:,jk) * umask(:,:,jk) 
    190          hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk) 
    191          hv_n(:,:) = hv_n(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk) 
     189         ht_n(:,:) = ht_n(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk) 
     190         hu_b(:,:) = hu_b(:,:) + e3u(:,:,jk,Kbb) * umask(:,:,jk) 
     191         hu_n(:,:) = hu_n(:,:) + e3u(:,:,jk,Kmm) * umask(:,:,jk) 
     192         hv_b(:,:) = hv_b(:,:) + e3v(:,:,jk,Kbb) * vmask(:,:,jk) 
     193         hv_n(:,:) = hv_n(:,:) + e3v(:,:,jk,Kmm) * vmask(:,:,jk) 
    192194      END DO 
    193195      ! 
     
    266268 
    267269 
    268    SUBROUTINE dom_vvl_sf_nxt( kt, kcall )  
     270   SUBROUTINE dom_vvl_sf_nxt( kt, Kbb, Kmm, Kaa, kcall )  
    269271      !!---------------------------------------------------------------------- 
    270272      !!                ***  ROUTINE dom_vvl_sf_nxt  *** 
     
    288290      !! Reference  : Leclair, M., and Madec, G. 2011, Ocean Modelling. 
    289291      !!---------------------------------------------------------------------- 
    290       INTEGER, INTENT( in )           ::   kt      ! time step 
    291       INTEGER, INTENT( in ), OPTIONAL ::   kcall   ! optional argument indicating call sequence 
     292      INTEGER, INTENT( in )           ::   kt             ! time step 
     293      INTEGER, INTENT( in )           ::   Kbb, Kmm, Kaa  ! time step 
     294      INTEGER, INTENT( in ), OPTIONAL ::   kcall          ! optional argument indicating call sequence 
    292295      ! 
    293296      INTEGER                ::   ji, jj, jk            ! dummy loop indices 
     
    321324      !                                                ! --------------------------------------------- ! 
    322325      ! 
    323       z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * ssmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) ) 
     326      z_scale(:,:) = ( ssh(:,:,Kaa) - ssh(:,:,Kbb) ) * ssmask(:,:) / ( ht_0(:,:) + ssh(:,:,Kmm) + 1. - ssmask(:,:) ) 
    324327      DO jk = 1, jpkm1 
    325          ! formally this is the same as e3t_a = e3t_0*(1+ssha/ht_0) 
    326          e3t_a(:,:,jk) = e3t_b(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk) 
     328         ! formally this is the same as e3t(:,:,:,Kaa) = e3t_0*(1+ssha/ht_0) 
     329         e3t(:,:,jk,Kaa) = e3t(:,:,jk,Kbb) + e3t(:,:,jk,Kmm) * z_scale(:,:) * tmask(:,:,jk) 
    327330      END DO 
    328331      ! 
     
    337340         zht(:,:)   = 0._wp 
    338341         DO jk = 1, jpkm1 
    339             zhdiv(:,:) = zhdiv(:,:) + e3t_n(:,:,jk) * hdivn(:,:,jk) 
    340             zht  (:,:) = zht  (:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 
     342            zhdiv(:,:) = zhdiv(:,:) + e3t(:,:,jk,Kmm) * hdiv(:,:,jk) 
     343            zht  (:,:) = zht  (:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk) 
    341344         END DO 
    342345         zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) ) 
     
    348351               DO jk = 1, jpkm1 
    349352                  hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:)   & 
    350                      &          * ( hdiv_lf(:,:,jk) - e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) ) 
     353                     &          * ( hdiv_lf(:,:,jk) - e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) ) 
    351354               END DO 
    352355            ENDIF 
     
    361364         IF( ln_vvl_ztilde ) THEN     ! z_tilde case 
    362365            DO jk = 1, jpkm1 
    363                tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) ) 
     366               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) ) 
    364367            END DO 
    365368         ELSE                         ! layer case 
    366369            DO jk = 1, jpkm1 
    367                tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) -   e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk) 
     370               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) -   e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk) 
    368371            END DO 
    369372         ENDIF 
     
    476479            zht(:,:)  = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk) 
    477480         END DO 
    478          z_scale(:,:) =  - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) ) 
     481         z_scale(:,:) =  - zht(:,:) / ( ht_0(:,:) + ssh(:,:,Kmm) + 1. - ssmask(:,:) ) 
    479482         DO jk = 1, jpkm1 
    480             dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk) 
     483            dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t(:,:,jk,Kmm) * z_scale(:,:) * tmask(:,:,jk) 
    481484         END DO 
    482485 
     
    486489      !                                           ! ---baroclinic part--------- ! 
    487490         DO jk = 1, jpkm1 
    488             e3t_a(:,:,jk) = e3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk) 
     491            e3t(:,:,jk,Kaa) = e3t(:,:,jk,Kaa) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk) 
    489492         END DO 
    490493      ENDIF 
     
    501504         zht(:,:) = 0.0_wp 
    502505         DO jk = 1, jpkm1 
    503             zht(:,:) = zht(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 
    504          END DO 
    505          z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshn(:,:) - zht(:,:) ) ) 
     506            zht(:,:) = zht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk) 
     507         END DO 
     508         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kmm) - zht(:,:) ) ) 
    506509         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain 
    507          IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(e3t_n))) =', z_tmax 
     510         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(e3t(:,:,:,Kmm)))) =', z_tmax 
    508511         ! 
    509512         zht(:,:) = 0.0_wp 
    510513         DO jk = 1, jpkm1 
    511             zht(:,:) = zht(:,:) + e3t_a(:,:,jk) * tmask(:,:,jk) 
    512          END DO 
    513          z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssha(:,:) - zht(:,:) ) ) 
     514            zht(:,:) = zht(:,:) + e3t(:,:,jk,Kaa) * tmask(:,:,jk) 
     515         END DO 
     516         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kaa) - zht(:,:) ) ) 
    514517         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain 
    515          IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(e3t_a))) =', z_tmax 
     518         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(e3t(:,:,:,Kaa)))) =', z_tmax 
    516519         ! 
    517520         zht(:,:) = 0.0_wp 
    518521         DO jk = 1, jpkm1 
    519             zht(:,:) = zht(:,:) + e3t_b(:,:,jk) * tmask(:,:,jk) 
    520          END DO 
    521          z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshb(:,:) - zht(:,:) ) ) 
     522            zht(:,:) = zht(:,:) + e3t(:,:,jk,Kbb) * tmask(:,:,jk) 
     523         END DO 
     524         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kbb) - zht(:,:) ) ) 
    522525         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain 
    523          IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(e3t_b))) =', z_tmax 
    524          ! 
    525          z_tmax = MAXVAL( tmask(:,:,1) *  ABS( sshb(:,:) ) ) 
     526         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(e3t(:,:,:,Kbb)))) =', z_tmax 
     527         ! 
     528         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssh(:,:,Kbb) ) ) 
    526529         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain 
    527          IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(sshb))) =', z_tmax 
    528          ! 
    529          z_tmax = MAXVAL( tmask(:,:,1) *  ABS( sshn(:,:) ) ) 
     530         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kbb)))) =', z_tmax 
     531         ! 
     532         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssh(:,:,Kmm) ) ) 
    530533         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain 
    531          IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(sshn))) =', z_tmax 
    532          ! 
    533          z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssha(:,:) ) ) 
     534         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kmm)))) =', z_tmax 
     535         ! 
     536         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssh(:,:,Kaa) ) ) 
    534537         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain 
    535          IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssha))) =', z_tmax 
     538         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kaa)))) =', z_tmax 
    536539      END IF 
    537540 
     
    540543      ! *********************************** ! 
    541544 
    542       CALL dom_vvl_interpol( e3t_a(:,:,:), e3u_a(:,:,:), 'U' ) 
    543       CALL dom_vvl_interpol( e3t_a(:,:,:), e3v_a(:,:,:), 'V' ) 
     545      CALL dom_vvl_interpol( e3t(:,:,:,Kaa), e3u(:,:,:,Kaa), 'U' ) 
     546      CALL dom_vvl_interpol( e3t(:,:,:,Kaa), e3v(:,:,:,Kaa), 'V' ) 
    544547 
    545548      ! *********************************** ! 
     
    547550      ! *********************************** ! 
    548551 
    549       hu_a(:,:) = e3u_a(:,:,1) * umask(:,:,1) 
    550       hv_a(:,:) = e3v_a(:,:,1) * vmask(:,:,1) 
     552      hu_a(:,:) = e3u(:,:,1,Kaa) * umask(:,:,1) 
     553      hv_a(:,:) = e3v(:,:,1,Kaa) * vmask(:,:,1) 
    551554      DO jk = 2, jpkm1 
    552          hu_a(:,:) = hu_a(:,:) + e3u_a(:,:,jk) * umask(:,:,jk) 
    553          hv_a(:,:) = hv_a(:,:) + e3v_a(:,:,jk) * vmask(:,:,jk) 
     555         hu_a(:,:) = hu_a(:,:) + e3u(:,:,jk,Kaa) * umask(:,:,jk) 
     556         hv_a(:,:) = hv_a(:,:) + e3v(:,:,jk,Kaa) * vmask(:,:,jk) 
    554557      END DO 
    555558      !                                        ! Inverse of the local depth 
     
    563566 
    564567 
    565    SUBROUTINE dom_vvl_sf_swp( kt ) 
     568   SUBROUTINE dom_vvl_sf_swp( kt, Kbb, Kmm, Kaa ) 
    566569      !!---------------------------------------------------------------------- 
    567570      !!                ***  ROUTINE dom_vvl_sf_swp  *** 
     
    578581      !!               - Recompute: 
    579582      !!                    e3(u/v)_b        
    580       !!                    e3w_n            
     583      !!                    e3w(:,:,:,Kmm)            
    581584      !!                    e3(u/v)w_b       
    582585      !!                    e3(u/v)w_n       
    583       !!                    gdept_n, gdepw_n  and gde3w_n 
     586      !!                    gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm)  and gde3w 
    584587      !!                    h(u/v) and h(u/v)r 
    585588      !! 
     
    587590      !!              Leclair, M., and G. Madec, 2011, Ocean Modelling. 
    588591      !!---------------------------------------------------------------------- 
    589       INTEGER, INTENT( in ) ::   kt   ! time step 
     592      INTEGER, INTENT( in ) ::   kt              ! time step 
     593      INTEGER, INTENT( in ) ::   Kbb, Kmm, Kaa   ! time level indices 
    590594      ! 
    591595      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     
    615619         tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:) 
    616620      ENDIF 
    617       gdept_b(:,:,:) = gdept_n(:,:,:) 
    618       gdepw_b(:,:,:) = gdepw_n(:,:,:) 
    619  
    620       e3t_n(:,:,:) = e3t_a(:,:,:) 
    621       e3u_n(:,:,:) = e3u_a(:,:,:) 
    622       e3v_n(:,:,:) = e3v_a(:,:,:) 
     621      gdept(:,:,:,Kbb) = gdept(:,:,:,Kmm) 
     622      gdepw(:,:,:,Kbb) = gdepw(:,:,:,Kmm) 
     623 
     624      e3t(:,:,:,Kmm) = e3t(:,:,:,Kaa) 
     625      e3u(:,:,:,Kmm) = e3u(:,:,:,Kaa) 
     626      e3v(:,:,:,Kmm) = e3v(:,:,:,Kaa) 
    623627 
    624628      ! Compute all missing vertical scale factor and depths 
     
    626630      ! Horizontal scale factor interpolations 
    627631      ! -------------------------------------- 
    628       ! - ML - e3u_b and e3v_b are allready computed in dynnxt 
     632      ! - ML - e3u(:,:,:,Kbb) and e3v(:,:,:,Kbb) are allready computed in dynnxt 
    629633      ! - JC - hu_b, hv_b, hur_b, hvr_b also 
    630634       
    631       CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F'  ) 
     635      CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F'  ) 
    632636       
    633637      ! Vertical scale factor interpolations 
    634       CALL dom_vvl_interpol( e3t_n(:,:,:),  e3w_n(:,:,:), 'W'  ) 
    635       CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' ) 
    636       CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' ) 
    637       CALL dom_vvl_interpol( e3t_b(:,:,:),  e3w_b(:,:,:), 'W'  ) 
    638       CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' ) 
    639       CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' ) 
     638      CALL dom_vvl_interpol( e3t(:,:,:,Kmm),  e3w(:,:,:,Kmm), 'W'  ) 
     639      CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' ) 
     640      CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' ) 
     641      CALL dom_vvl_interpol( e3t(:,:,:,Kbb),  e3w(:,:,:,Kbb), 'W'  ) 
     642      CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' ) 
     643      CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' ) 
    640644 
    641645      ! t- and w- points depth (set the isf depth as it is in the initial step) 
    642       gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1) 
    643       gdepw_n(:,:,1) = 0.0_wp 
    644       gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:) 
     646      gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm) 
     647      gdepw(:,:,1,Kmm) = 0.0_wp 
     648      gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm) 
    645649      DO jk = 2, jpk 
    646650         DO jj = 1,jpj 
     
    649653                                                                 ! 1 for jk = mikt 
    650654               zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 
    651                gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1) 
    652                gdept_n(ji,jj,jk) =    zcoef  * ( gdepw_n(ji,jj,jk  ) + 0.5 * e3w_n(ji,jj,jk) )  & 
    653                    &             + (1-zcoef) * ( gdept_n(ji,jj,jk-1) +       e3w_n(ji,jj,jk) )  
    654                gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj) 
     655               gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 
     656               gdept(ji,jj,jk,Kmm) =    zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm) )  & 
     657                   &             + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm) )  
     658               gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm) 
    655659            END DO 
    656660         END DO 
     
    662666      hv_n(:,:) = hv_a(:,:)   ;   r1_hv_n(:,:) = r1_hv_a(:,:) 
    663667      ! 
    664       ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1) 
     668      ht_n(:,:) = e3t(:,:,1,Kmm) * tmask(:,:,1) 
    665669      DO jk = 2, jpkm1 
    666          ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 
     670         ht_n(:,:) = ht_n(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk) 
    667671      END DO 
    668672 
    669673      ! write restart file 
    670674      ! ================== 
    671       IF( lrst_oce  )   CALL dom_vvl_rst( kt, 'WRITE' ) 
     675      IF( lrst_oce  )   CALL dom_vvl_rst( kt, Kbb, Kmm, 'WRITE' ) 
    672676      ! 
    673677      IF( ln_timing )   CALL timing_stop('dom_vvl_sf_swp') 
     
    783787 
    784788 
    785    SUBROUTINE dom_vvl_rst( kt, cdrw ) 
     789   SUBROUTINE dom_vvl_rst( kt, Kbb, Kmm, cdrw ) 
    786790      !!--------------------------------------------------------------------- 
    787791      !!                   ***  ROUTINE dom_vvl_rst  *** 
     
    795799      !!                they are set to 0. 
    796800      !!---------------------------------------------------------------------- 
    797       INTEGER         , INTENT(in) ::   kt     ! ocean time-step 
    798       CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag 
     801      INTEGER         , INTENT(in) ::   kt        ! ocean time-step 
     802      INTEGER         , INTENT(in) ::   Kbb, Kmm  ! ocean time level indices 
     803      CHARACTER(len=*), INTENT(in) ::   cdrw      ! "READ"/"WRITE" flag 
    799804      ! 
    800805      INTEGER ::   ji, jj, jk 
     
    806811         IF( ln_rstart ) THEN                   !* Read the restart file 
    807812            CALL rst_read_open                  !  open the restart file if necessary 
    808             CALL iom_get( numror, jpdom_autoglo, 'sshn'   , sshn, ldxios = lrxios    ) 
     813            CALL iom_get( numror, jpdom_autoglo, 'sshn'   , ssh(:,:,Kmm), ldxios = lrxios    ) 
    809814            ! 
    810815            id1 = iom_varid( numror, 'e3t_b', ldstop = .FALSE. ) 
     
    817822            !                             ! --------- ! 
    818823            IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist 
    819                CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:), ldxios = lrxios ) 
    820                CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:), ldxios = lrxios ) 
     824               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios ) 
     825               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios ) 
    821826               ! needed to restart if land processor not computed  
    822                IF(lwp) write(numout,*) 'dom_vvl_rst : e3t_b and e3t_n found in restart files' 
     827               IF(lwp) write(numout,*) 'dom_vvl_rst : e3t(:,:,:,Kbb) and e3t(:,:,:,Kmm) found in restart files' 
    823828               WHERE ( tmask(:,:,:) == 0.0_wp )  
    824                   e3t_n(:,:,:) = e3t_0(:,:,:) 
    825                   e3t_b(:,:,:) = e3t_0(:,:,:) 
     829                  e3t(:,:,:,Kmm) = e3t_0(:,:,:) 
     830                  e3t(:,:,:,Kbb) = e3t_0(:,:,:) 
    826831               END WHERE 
    827832               IF( neuler == 0 ) THEN 
    828                   e3t_b(:,:,:) = e3t_n(:,:,:) 
     833                  e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 
    829834               ENDIF 
    830835            ELSE IF( id1 > 0 ) THEN 
    831                IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart files' 
     836               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart files' 
    832837               IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.' 
    833838               IF(lwp) write(numout,*) 'neuler is forced to 0' 
    834                CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:), ldxios = lrxios ) 
    835                e3t_n(:,:,:) = e3t_b(:,:,:) 
     839               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios ) 
     840               e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb) 
    836841               neuler = 0 
    837842            ELSE IF( id2 > 0 ) THEN 
    838                IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_b not found in restart files' 
     843               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kbb) not found in restart files' 
    839844               IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.' 
    840845               IF(lwp) write(numout,*) 'neuler is forced to 0' 
    841                CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:), ldxios = lrxios ) 
    842                e3t_b(:,:,:) = e3t_n(:,:,:) 
     846               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios ) 
     847               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 
    843848               neuler = 0 
    844849            ELSE 
    845                IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart file' 
     850               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart file' 
    846851               IF(lwp) write(numout,*) 'Compute scale factor from sshn' 
    847852               IF(lwp) write(numout,*) 'neuler is forced to 0' 
    848853               DO jk = 1, jpk 
    849                   e3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) & 
     854                  e3t(:,:,jk,Kmm) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm) ) & 
    850855                      &                          / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   & 
    851856                      &          + e3t_0(:,:,jk)                               * (1._wp -tmask(:,:,jk)) 
    852857               END DO 
    853                e3t_b(:,:,:) = e3t_n(:,:,:) 
     858               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 
    854859               neuler = 0 
    855860            ENDIF 
     
    888893               IF( cn_cfg == 'wad' ) THEN 
    889894                  ! Wetting and drying test case 
    890                   CALL usr_def_istate( gdept_b, tmask, tsb, ub, vb, sshb  ) 
    891                   tsn  (:,:,:,:) = tsb (:,:,:,:)       ! set now values from to before ones 
    892                   sshn (:,:)     = sshb(:,:) 
    893                   un   (:,:,:)   = ub  (:,:,:) 
    894                   vn   (:,:,:)   = vb  (:,:,:) 
     895                  CALL usr_def_istate( gdept(:,:,:,Kbb), tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  ) 
     896                  ts  (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)       ! set now values from to before ones 
     897                  ssh (:,:,Kmm)     = ssh(:,:,Kbb) 
     898                  uu   (:,:,:,Kmm)   = uu  (:,:,:,Kbb) 
     899                  vv   (:,:,:,Kmm)   = vv  (:,:,:,Kbb) 
    895900               ELSE 
    896901                  ! if not test case 
    897                   sshn(:,:) = -ssh_ref 
    898                   sshb(:,:) = -ssh_ref 
     902                  ssh(:,:,Kmm) = -ssh_ref 
     903                  ssh(:,:,Kbb) = -ssh_ref 
    899904 
    900905                  DO jj = 1, jpj 
    901906                     DO ji = 1, jpi 
    902907                        IF( ht_0(ji,jj)-ssh_ref <  rn_wdmin1 ) THEN ! if total depth is less than min depth 
    903  
    904                            sshb(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) ) 
    905                            sshn(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) ) 
    906                            ssha(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) ) 
     908                           ssh(ji,jj,Kbb) = rn_wdmin1 - (ht_0(ji,jj) ) 
     909                           ssh(ji,jj,Kmm) = rn_wdmin1 - (ht_0(ji,jj) ) 
    907910                        ENDIF 
    908911                     ENDDO 
     
    912915               ! Adjust vertical metrics for all wad 
    913916               DO jk = 1, jpk 
    914                   e3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:)  ) & 
     917                  e3t(:,:,jk,Kmm) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm)  ) & 
    915918                    &                            / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   & 
    916919                    &            + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) ) 
    917920               END DO 
    918                e3t_b(:,:,:) = e3t_n(:,:,:) 
     921               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 
    919922 
    920923               DO ji = 1, jpi 
     
    930933               ! Just to read set ssh in fact, called latter once vertical grid 
    931934               ! is set up: 
    932 !               CALL usr_def_istate( gdept_0, tmask, tsb, ub, vb, sshb  ) 
     935!               CALL usr_def_istate( gdept_0, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  ) 
    933936!               ! 
    934937!               DO jk=1,jpk 
    935 !                  e3t_b(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshb(:,:) ) & 
     938!                  e3t(:,:,jk,Kbb) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kbb) ) & 
    936939!                     &            / ( ht_0(:,:) + 1._wp -ssmask(:,:) ) * tmask(:,:,jk) 
    937940!               END DO 
    938 !               e3t_n(:,:,:) = e3t_b(:,:,:) 
    939                 sshn(:,:)=0._wp 
    940                 e3t_n(:,:,:)=e3t_0(:,:,:) 
    941                 e3t_b(:,:,:)=e3t_0(:,:,:) 
     941!               e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb) 
     942                ssh(:,:,Kmm)=0._wp 
     943                e3t(:,:,:,Kmm)=e3t_0(:,:,:) 
     944                e3t(:,:,:,Kbb)=e3t_0(:,:,:) 
    942945               ! 
    943946            END IF           ! end of ll_wd edits 
     
    957960         !                                           ! all cases ! 
    958961         !                                           ! --------- ! 
    959          CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t_b(:,:,:), ldxios = lwxios ) 
    960          CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t_n(:,:,:), ldxios = lwxios ) 
     962         CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lwxios ) 
     963         CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lwxios ) 
    961964         !                                           ! ----------------------- ! 
    962965         IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN  ! z_tilde and layer cases ! 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM/iscplhsb.F90

    r10425 r10978  
    4040CONTAINS 
    4141 
    42    SUBROUTINE iscpl_cons(ptmask_b, psmask_b, pe3t_b, pts_flx, pvol_flx, prdt_iscpl) 
     42   SUBROUTINE iscpl_cons( Kbb, Kmm, ptmask_b, psmask_b, pe3t_b, pts_flx, pvol_flx, prdt_iscpl ) 
    4343      !!----------------------------------------------------------------------  
    4444      !!                   ***  ROUTINE iscpl_cons  *** 
     
    4848      !!                compute where the correction have to be applied 
    4949      !!  
    50       !! ** Method  :   compute tsn*e3t-tsb*e3tb and e3t-e3t_b 
    51       !!---------------------------------------------------------------------- 
     50      !! ** Method  :   compute tsn*e3tn-tsb*e3tb and e3tn-e3tb 
     51      !!---------------------------------------------------------------------- 
     52      INTEGER                     , INTENT(in ) :: Kbb, Kmm    !! time level indices 
    5253      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: ptmask_b    !! mask before 
    5354      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: pe3t_b      !! scale factor before 
     
    7879      z1_rdtiscpl = 1._wp / prdt_iscpl  
    7980 
    80       ! mask tsn and tsb  
    81       tsb(:,:,:,jp_tem) = tsb(:,:,:,jp_tem) * ptmask_b(:,:,:) 
    82       tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) *  tmask  (:,:,:) 
    83       tsb(:,:,:,jp_sal) = tsb(:,:,:,jp_sal) * ptmask_b(:,:,:) 
    84       tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) *  tmask  (:,:,:) 
     81      ! mask ts(:,:,:,:,Kmm) and ts(:,:,:,:,Kbb)  
     82      ts(:,:,:,jp_tem,Kbb) = ts(:,:,:,jp_tem,Kbb) * ptmask_b(:,:,:) 
     83      ts(:,:,:,jp_tem,Kmm) = ts(:,:,:,jp_tem,Kmm) *  tmask  (:,:,:) 
     84      ts(:,:,:,jp_sal,Kbb) = ts(:,:,:,jp_sal,Kbb) * ptmask_b(:,:,:) 
     85      ts(:,:,:,jp_sal,Kmm) = ts(:,:,:,jp_sal,Kmm) *  tmask  (:,:,:) 
    8586 
    8687      !============================================================================== 
     
    8990 
    9091      !  
    91       zdssh(:,:) = sshn(:,:) * ssmask(:,:) - sshb(:,:) * psmask_b(:,:) 
     92      zdssh(:,:) = ssh(:,:,Kmm) * ssmask(:,:) - ssh(:,:,Kbb) * psmask_b(:,:) 
    9293      IF (.NOT. ln_linssh ) zdssh = 0.0_wp ! already included in the levels by definition 
    9394       
     
    9899 
    99100                  ! volume differences 
    100                   zde3t = e3t_n(ji,jj,jk) * tmask(ji,jj,jk) - pe3t_b(ji,jj,jk) * ptmask_b(ji,jj,jk) 
     101                  zde3t = e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk) - pe3t_b(ji,jj,jk) * ptmask_b(ji,jj,jk) 
    101102 
    102103                  ! heat diff 
    103                   zdtem = tsn(ji,jj,jk,jp_tem) * e3t_n(ji,jj,jk) *  tmask  (ji,jj,jk)   & 
    104                         - tsb(ji,jj,jk,jp_tem) * pe3t_b (ji,jj,jk) * ptmask_b(ji,jj,jk) 
     104                  zdtem = ts(ji,jj,jk,jp_tem,Kmm) * e3t(ji,jj,jk,Kmm) *  tmask  (ji,jj,jk)   & 
     105                        - ts(ji,jj,jk,jp_tem,Kbb) * pe3t_b (ji,jj,jk) * ptmask_b(ji,jj,jk) 
    105106                  ! salt diff 
    106                   zdsal = tsn(ji,jj,jk,jp_sal) * e3t_n(ji,jj,jk) *  tmask  (ji,jj,jk)   & 
    107                         - tsb(ji,jj,jk,jp_sal) * pe3t_b (ji,jj,jk) * ptmask_b(ji,jj,jk) 
     107                  zdsal = ts(ji,jj,jk,jp_sal,Kmm) * e3t(ji,jj,jk,Kmm) *  tmask  (ji,jj,jk)   & 
     108                        - ts(ji,jj,jk,jp_sal,Kbb) * pe3t_b (ji,jj,jk) * ptmask_b(ji,jj,jk) 
    108109                
    109110                  ! shh changes 
     
    296297 
    297298 
    298    SUBROUTINE iscpl_div( phdivn ) 
     299   SUBROUTINE iscpl_div( Kmm, phdivn ) 
    299300      !!---------------------------------------------------------------------- 
    300301      !!                  ***  ROUTINE iscpl_div  *** 
     
    307308      !! ** Action  :   phdivn   increase by the iscpl correction term 
    308309      !!---------------------------------------------------------------------- 
     310      INTEGER                   , INTENT(in   ) ::   Kmm      ! time level index 
    309311      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   phdivn   ! horizontal divergence 
    310312      !! 
     
    315317         DO jj = 1, jpj 
    316318            DO ji = 1, jpi 
    317                phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + hdiv_iscpl(ji,jj,jk) / e3t_n(ji,jj,jk) 
     319               phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + hdiv_iscpl(ji,jj,jk) / e3t(ji,jj,jk,Kmm) 
    318320            END DO 
    319321         END DO 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM/iscplrst.F90

    r10425 r10978  
    4040CONTAINS 
    4141 
    42    SUBROUTINE iscpl_stp 
     42   SUBROUTINE iscpl_stp( Kbb, Kmm ) 
    4343      !!----------------------------------------------------------------------  
    4444      !!                   ***  ROUTINE iscpl_stp  *** 
    4545      !!  
    4646      !! ** Purpose : compute initialisation 
    47       !!              compute extrapolation of restart variable un, vn, tsn, sshn (wetting/drying)    
     47      !!              compute extrapolation of restart variable uu(Kmm), vv(Kmm), ts(Kmm), ssh(Kmm) (wetting/drying)    
    4848      !!              compute correction term if needed 
    4949      !!  
    5050      !!---------------------------------------------------------------------- 
     51      INTEGER, INTENT(in) :: Kbb, Kmm   ! time level indices 
     52      ! 
    5153      INTEGER  ::   inum0 
    5254      REAL(wp), DIMENSION(jpi,jpj)     ::   zsmask_b 
     
    6971      CALL iscpl_init()       ! read namelist 
    7072      !                       ! Extrapolation/interpolation of modify cell and new cells ... (maybe do it later after domvvl) 
    71       CALL iscpl_rst_interpol( ztmask_b, zumask_b, zvmask_b, zsmask_b, ze3t_b, ze3u_b, ze3v_b, zdepw_b ) 
     73      CALL iscpl_rst_interpol( Kbb, Kmm, ztmask_b, zumask_b, zvmask_b, zsmask_b, ze3t_b, ze3u_b, ze3v_b, zdepw_b ) 
    7274      ! 
    7375      IF ( ln_hsb ) THEN      ! compute correction if conservation needed 
    7476         IF( iscpl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'rst_iscpl : unable to allocate rst_iscpl arrays' ) 
    75          CALL iscpl_cons(ztmask_b, zsmask_b, ze3t_b, htsc_iscpl, hdiv_iscpl, rdt_iscpl) 
     77         CALL iscpl_cons( Kbb, Kmm, ztmask_b, zsmask_b, ze3t_b, htsc_iscpl, hdiv_iscpl, rdt_iscpl ) 
    7678      END IF 
    7779          
     
    9294      ! 
    9395      !                       ! set _b and _n variables equal 
    94       tsb (:,:,:,:) = tsn (:,:,:,:) 
    95       ub  (:,:,:)   = un  (:,:,:) 
    96       vb  (:,:,:)   = vn  (:,:,:) 
    97       sshb(:,:)     = sshn(:,:) 
     96      ts (:,:,:,:,Kbb) = ts (:,:,:,:,Kmm) 
     97      uu  (:,:,:,Kbb)   = uu  (:,:,:,Kmm) 
     98      vv  (:,:,:,Kbb)   = vv  (:,:,:,Kmm) 
     99      ssh(:,:,Kbb)     = ssh(:,:,Kmm) 
    98100      ! 
    99101      !                       ! set _b and _n vertical scale factor equal 
    100       e3t_b (:,:,:) = e3t_n (:,:,:) 
    101       e3u_b (:,:,:) = e3u_n (:,:,:) 
    102       e3v_b (:,:,:) = e3v_n (:,:,:) 
    103       ! 
    104       e3uw_b (:,:,:) = e3uw_n (:,:,:) 
    105       e3vw_b (:,:,:) = e3vw_n (:,:,:) 
    106       gdept_b(:,:,:) = gdept_n(:,:,:) 
    107       gdepw_b(:,:,:) = gdepw_n(:,:,:) 
     102      e3t (:,:,:,Kbb) = e3t (:,:,:,Kmm) 
     103      e3u (:,:,:,Kbb) = e3u (:,:,:,Kmm) 
     104      e3v (:,:,:,Kbb) = e3v (:,:,:,Kmm) 
     105      ! 
     106      e3uw (:,:,:,Kbb) = e3uw (:,:,:,Kmm) 
     107      e3vw (:,:,:,Kbb) = e3vw (:,:,:,Kmm) 
     108      gdept(:,:,:,Kbb) = gdept(:,:,:,Kmm) 
     109      gdepw(:,:,:,Kbb) = gdepw(:,:,:,Kmm) 
    108110      hu_b   (:,:)   = hu_n   (:,:) 
    109111      hv_b   (:,:)   = hv_n   (:,:) 
     
    114116 
    115117 
    116    SUBROUTINE iscpl_rst_interpol (ptmask_b, pumask_b, pvmask_b, psmask_b, pe3t_b, pe3u_b, pe3v_b, pdepw_b) 
     118   SUBROUTINE iscpl_rst_interpol ( Kbb, Kmm, ptmask_b, pumask_b, pvmask_b, psmask_b, pe3t_b, pe3u_b, pe3v_b, pdepw_b ) 
    117119      !!----------------------------------------------------------------------  
    118120      !!                   ***  ROUTINE iscpl_rst_interpol  *** 
    119121      !!  
    120       !! ** Purpose :   compute new tn, sn, un, vn and sshn in case of evolving geometry of ice shelves  
     122      !! ** Purpose :   compute new ts(Kmm), uu(Kmm), vv(Kmm) and ssh(Kmm) in case of evolving geometry of ice shelves  
    121123      !!                compute 2d fields of heat, salt and volume correction 
    122124      !!  
    123       !! ** Method  :   tn, sn : extrapolation from neigbourg cells 
    124       !!                un, vn : fill with 0 velocity and keep barotropic transport by modifing surface velocity or adjacent velocity 
     125      !! ** Method  :   ts(Kmm) : extrapolation from neigbourg cells 
     126      !!                uu(Kmm), vv(Kmm) : fill with 0 velocity and keep barotropic transport by modifing surface velocity or adjacent velocity 
    125127      !!---------------------------------------------------------------------- 
     128      INTEGER                     , INTENT(in ) :: Kbb, Kmm                        !! time level indices 
    126129      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: ptmask_b, pumask_b, pvmask_b    !! mask before 
    127130      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: pe3t_b  , pe3u_b  , pe3v_b      !! scale factor before 
     
    143146      ! 
    144147      !                 ! mask value to be sure 
    145       tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) * ptmask_b(:,:,:) 
    146       tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * ptmask_b(:,:,:) 
     148      ts(:,:,:,jp_tem,Kmm) = ts(:,:,:,jp_tem,Kmm) * ptmask_b(:,:,:) 
     149      ts(:,:,:,jp_sal,Kmm) = ts(:,:,:,jp_sal,Kmm) * ptmask_b(:,:,:) 
    147150      ! 
    148151      !                 ! compute wmask 
     
    155158      !     
    156159      !                 ! compute new ssh if we open a full water column (average of the closest neigbourgs)   
    157       sshb (:,:)=sshn(:,:) 
    158       zssh0(:,:)=sshn(:,:) 
     160      ssh (:,:,Kbb)=ssh(:,:,Kmm) 
     161      zssh0(:,:)=ssh(:,:,Kmm) 
    159162      zsmask0(:,:) = psmask_b(:,:) 
    160163      zsmask1(:,:) = psmask_b(:,:) 
     
    167170               summsk=(zsmask0(jip1,jj)+zsmask0(jim1,jj)+zsmask0(ji,jjp1)+zsmask0(ji,jjm1)) 
    168171               IF (zdsmask(ji,jj) == 1._wp .AND. summsk /= 0._wp) THEN 
    169                   sshn(ji,jj)=( zssh0(jip1,jj)*zsmask0(jip1,jj)     & 
     172                  ssh(ji,jj,Kmm)=( zssh0(jip1,jj)*zsmask0(jip1,jj)     & 
    170173                  &           + zssh0(jim1,jj)*zsmask0(jim1,jj)     & 
    171174                  &           + zssh0(ji,jjp1)*zsmask0(ji,jjp1)     & 
     
    175178            END DO 
    176179         END DO 
    177          CALL lbc_lnk_multi( 'iscplrst', sshn, 'T', 1., zsmask1, 'T', 1. ) 
    178          zssh0   = sshn 
     180         CALL lbc_lnk_multi( 'iscplrst', ssh(:,:,Kmm), 'T', 1., zsmask1, 'T', 1. ) 
     181         zssh0   = ssh(:,:,Kmm) 
    179182         zsmask0 = zsmask1 
    180183      END DO 
    181       sshn(:,:) = sshn(:,:) * ssmask(:,:) 
     184      ssh(:,:,Kmm) = ssh(:,:,Kmm) * ssmask(:,:) 
    182185 
    183186!============================================================================= 
     
    192195               DO ji=1,jpi 
    193196                  IF (tmask(ji,jj,1) == 0._wp .OR. ptmask_b(ji,jj,1) == 0._wp) THEN 
    194                      e3t_n(ji,jj,jk) = e3t_0(ji,jj,jk) * ( 1._wp + sshn(ji,jj) /   & 
     197                     e3t(ji,jj,jk,Kmm) = e3t_0(ji,jj,jk) * ( 1._wp + ssh(ji,jj,Kmm) /   & 
    195198                     &   ( ht_0(ji,jj) + 1._wp - ssmask(ji,jj) ) * tmask(ji,jj,jk) ) 
    196199                  ENDIF 
     
    199202         END DO 
    200203         ! 
    201          CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' ) 
    202          CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' ) 
    203          CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' ) 
     204         CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' ) 
     205         CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' ) 
     206         CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' ) 
    204207 
    205208         ! Vertical scale factor interpolations 
    206209         ! ------------------------------------ 
    207          CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W'  ) 
    208          CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' ) 
    209          CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' ) 
     210         CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W'  ) 
     211         CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' ) 
     212         CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' ) 
    210213          
    211214         ! t- and w- points depth 
    212215         ! ---------------------- 
    213          gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1) 
    214          gdepw_n(:,:,1) = 0.0_wp 
    215          gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:) 
     216         gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm) 
     217         gdepw(:,:,1,Kmm) = 0.0_wp 
     218         gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm) 
    216219         DO jj = 1,jpj 
    217220            DO ji = 1,jpi 
    218221               DO jk = 2,mikt(ji,jj)-1 
    219                   gdept_n(ji,jj,jk) = gdept_0(ji,jj,jk) 
    220                   gdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk) 
    221                   gde3w_n(ji,jj,jk) = gdept_0(ji,jj,jk) - sshn(ji,jj) 
     222                  gdept(ji,jj,jk,Kmm) = gdept_0(ji,jj,jk) 
     223                  gdepw(ji,jj,jk,Kmm) = gdepw_0(ji,jj,jk) 
     224                  gde3w(ji,jj,jk) = gdept_0(ji,jj,jk) - ssh(ji,jj,Kmm) 
    222225               END DO 
    223226               IF (mikt(ji,jj) > 1) THEN 
    224227                  jk = mikt(ji,jj) 
    225                   gdept_n(ji,jj,jk) = gdepw_0(ji,jj,jk) + 0.5_wp * e3w_n(ji,jj,jk) 
    226                   gdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk) 
    227                   gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk  ) - sshn   (ji,jj) 
     228                  gdept(ji,jj,jk,Kmm) = gdepw_0(ji,jj,jk) + 0.5_wp * e3w(ji,jj,jk,Kmm) 
     229                  gdepw(ji,jj,jk,Kmm) = gdepw_0(ji,jj,jk) 
     230                  gde3w(ji,jj,jk) = gdept(ji,jj,jk  ,Kmm) - ssh   (ji,jj,Kmm) 
    228231               END IF 
    229232               DO jk = mikt(ji,jj)+1, jpk 
    230                   gdept_n(ji,jj,jk) = gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk) 
    231                   gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1) 
    232                   gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk  ) - sshn (ji,jj) 
     233                  gdept(ji,jj,jk,Kmm) = gdept(ji,jj,jk-1,Kmm) + e3w(ji,jj,jk,Kmm) 
     234                  gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 
     235                  gde3w(ji,jj,jk) = gdept(ji,jj,jk  ,Kmm) - ssh (ji,jj,Kmm) 
    233236               END DO 
    234237            END DO 
     
    239242         ht_n(:,:) = 0._wp ; hu_n(:,:) = 0._wp ; hv_n(:,:) = 0._wp 
    240243         DO jk = 1, jpkm1 
    241             hu_n(:,:) = hu_n(:,:) + e3u_n(:,:,jk) * umask(:,:,jk) 
    242             hv_n(:,:) = hv_n(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk) 
    243             ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 
     244            hu_n(:,:) = hu_n(:,:) + e3u(:,:,jk,Kmm) * umask(:,:,jk) 
     245            hv_n(:,:) = hv_n(:,:) + e3v(:,:,jk,Kmm) * vmask(:,:,jk) 
     246            ht_n(:,:) = ht_n(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk) 
    244247         END DO 
    245248         !                                        ! Inverse of the local depth 
     
    252255! compute velocity 
    253256! compute velocity in order to conserve barotropic velocity (modification by poderation of the scale factor). 
    254       ub(:,:,:)=un(:,:,:) 
    255       vb(:,:,:)=vn(:,:,:) 
     257      uu(:,:,:,Kbb)=uu(:,:,:,Kmm) 
     258      vv(:,:,:,Kbb)=vv(:,:,:,Kmm) 
    256259      DO jk = 1,jpk 
    257260         DO jj = 1,jpj 
    258261            DO ji = 1,jpi 
    259                un(ji,jj,jk) = ub(ji,jj,jk)*pe3u_b(ji,jj,jk)*pumask_b(ji,jj,jk)/e3u_n(ji,jj,jk)*umask(ji,jj,jk); 
    260                vn(ji,jj,jk) = vb(ji,jj,jk)*pe3v_b(ji,jj,jk)*pvmask_b(ji,jj,jk)/e3v_n(ji,jj,jk)*vmask(ji,jj,jk); 
     262               uu(ji,jj,jk,Kmm) = uu(ji,jj,jk,Kbb)*pe3u_b(ji,jj,jk)*pumask_b(ji,jj,jk)/e3u(ji,jj,jk,Kmm)*umask(ji,jj,jk); 
     263               vv(ji,jj,jk,Kmm) = vv(ji,jj,jk,Kbb)*pe3v_b(ji,jj,jk)*pvmask_b(ji,jj,jk)/e3v(ji,jj,jk,Kmm)*vmask(ji,jj,jk); 
    261264            END DO 
    262265         END DO 
     
    265268! compute new velocity if we close a cell (check barotropic velocity and change velocity over the water column) 
    266269! compute barotropic velocity now and after  
    267       ztrp(:,:,:) = ub(:,:,:)*pe3u_b(:,:,:);  
     270      ztrp(:,:,:) = uu(:,:,:,Kbb)*pe3u_b(:,:,:);  
    268271      zbub(:,:)   = SUM(ztrp,DIM=3) 
    269       ztrp(:,:,:) = vb(:,:,:)*pe3v_b(:,:,:);  
     272      ztrp(:,:,:) = vv(:,:,:,Kbb)*pe3v_b(:,:,:);  
    270273      zbvb(:,:)   = SUM(ztrp,DIM=3) 
    271       ztrp(:,:,:) = un(:,:,:)*e3u_n(:,:,:);  
     274      ztrp(:,:,:) = uu(:,:,:,Kmm)*e3u(:,:,:,Kmm);  
    272275      zbun(:,:)   = SUM(ztrp,DIM=3) 
    273       ztrp(:,:,:) = vn(:,:,:)*e3v_n(:,:,:);  
     276      ztrp(:,:,:) = vv(:,:,:,Kmm)*e3v(:,:,:,Kmm);  
    274277      zbvn(:,:)   = SUM(ztrp,DIM=3) 
    275278 
     
    278281      zhv1=0.0_wp ; 
    279282      DO jk  = 1,jpk 
    280         zhu1(:,:) = zhu1(:,:) + e3u_n(:,:,jk) * umask(:,:,jk) 
    281         zhv1(:,:) = zhv1(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk) 
     283        zhu1(:,:) = zhu1(:,:) + e3u(:,:,jk,Kmm) * umask(:,:,jk) 
     284        zhv1(:,:) = zhv1(:,:) + e3v(:,:,jk,Kmm) * vmask(:,:,jk) 
    282285      END DO 
    283286       
     
    298301      ! update velocity 
    299302      DO jk  = 1,jpk 
    300          un(:,:,jk)=(un(:,:,jk) - zucorr(:,:))*umask(:,:,jk) 
    301          vn(:,:,jk)=(vn(:,:,jk) - zvcorr(:,:))*vmask(:,:,jk) 
     303         uu(:,:,jk,Kmm)=(uu(:,:,jk,Kmm) - zucorr(:,:))*umask(:,:,jk) 
     304         vv(:,:,jk,Kmm)=(vv(:,:,jk,Kmm) - zvcorr(:,:))*vmask(:,:,jk) 
    302305      END DO 
    303306 
     
    305308      ! compute temp and salt 
    306309      ! compute new tn and sn if we open a new cell 
    307       tsb (:,:,:,:) = tsn(:,:,:,:) 
    308       zts0(:,:,:,:) = tsn(:,:,:,:) 
     310      ts (:,:,:,:,Kbb) = ts(:,:,:,:,Kmm) 
     311      zts0(:,:,:,:) = ts(:,:,:,:,Kmm) 
    309312      ztmask1(:,:,:) = ptmask_b(:,:,:) 
    310313      ztmask0(:,:,:) = ptmask_b(:,:,:) 
     
    319322                      IF (zdmask(ji,jj) == 1._wp .AND. summsk /= 0._wp) THEN 
    320323                      !! horizontal basic extrapolation 
    321                          tsn(ji,jj,jk,1)=( zts0(jip1,jj  ,jk,1)*ztmask0(jip1,jj  ,jk) & 
     324                         ts(ji,jj,jk,1,Kmm)=( zts0(jip1,jj  ,jk,1)*ztmask0(jip1,jj  ,jk) & 
    322325                         &                +zts0(jim1,jj  ,jk,1)*ztmask0(jim1,jj  ,jk) & 
    323326                         &                +zts0(ji  ,jjp1,jk,1)*ztmask0(ji  ,jjp1,jk) & 
    324327                         &                +zts0(ji  ,jjm1,jk,1)*ztmask0(ji  ,jjm1,jk) ) / summsk 
    325                          tsn(ji,jj,jk,2)=( zts0(jip1,jj  ,jk,2)*ztmask0(jip1,jj  ,jk) & 
     328                         ts(ji,jj,jk,2,Kmm)=( zts0(jip1,jj  ,jk,2)*ztmask0(jip1,jj  ,jk) & 
    326329                         &                +zts0(jim1,jj  ,jk,2)*ztmask0(jim1,jj  ,jk) & 
    327330                         &                +zts0(ji  ,jjp1,jk,2)*ztmask0(ji  ,jjp1,jk) & 
     
    333336                         summsk=(ztmask0(ji,jj,jkm1)+ztmask0(ji,jj,jkp1)) 
    334337                         IF (zdmask(ji,jj) == 1._wp .AND. summsk /= 0._wp ) THEN 
    335                             tsn(ji,jj,jk,1)=( zts0(ji,jj,jkp1,1)*ztmask0(ji,jj,jkp1)     & 
     338                            ts(ji,jj,jk,1,Kmm)=( zts0(ji,jj,jkp1,1)*ztmask0(ji,jj,jkp1)     & 
    336339                            &                +zts0(ji,jj,jkm1,1)*ztmask0(ji,jj,jkm1))/summsk 
    337                             tsn(ji,jj,jk,2)=( zts0(ji,jj,jkp1,2)*ztmask0(ji,jj,jkp1)     & 
     340                            ts(ji,jj,jk,2,Kmm)=( zts0(ji,jj,jkp1,2)*ztmask0(ji,jj,jkp1)     & 
    338341                            &                +zts0(ji,jj,jkm1,2)*ztmask0(ji,jj,jkm1))/summsk 
    339342                            ztmask1(ji,jj,jk)=1._wp 
     
    344347          END DO 
    345348           
    346           CALL lbc_lnk_multi( 'iscplrst', tsn(:,:,:,jp_tem), 'T', 1., tsn(:,:,:,jp_sal), 'T', 1., ztmask1, 'T', 1.) 
     349          CALL lbc_lnk_multi( 'iscplrst', ts(:,:,:,jp_tem,Kmm), 'T', 1., ts(:,:,:,jp_sal,Kmm), 'T', 1., ztmask1, 'T', 1.) 
    347350 
    348351          ! update 
    349           zts0(:,:,:,:) = tsn(:,:,:,:) 
     352          zts0(:,:,:,:) = ts(:,:,:,:,Kmm) 
    350353          ztmask0 = ztmask1 
    351354 
    352355      END DO 
    353356 
    354       ! mask new tsn field 
    355       tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) * tmask(:,:,:) 
    356       tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * tmask(:,:,:) 
     357      ! mask new ts(:,:,:,:,Kmm) field 
     358      ts(:,:,:,jp_tem,Kmm) = ts(:,:,:,jp_tem,Kmm) * tmask(:,:,:) 
     359      ts(:,:,:,jp_sal,Kmm) = ts(:,:,:,jp_sal,Kmm) * tmask(:,:,:) 
    357360 
    358361      ! compute new T/S (interpolation) if vvl only for common wet cell in before and after wmask 
     
    365368                  &      (tmask(ji,jj,1)==0._wp .OR. ptmask_b(ji,jj,1)==0._wp) ) THEN 
    366369                     !compute weight 
    367                      zdzp1 = MAX(0._wp,gdepw_n(ji,jj,jk+1) - pdepw_b(ji,jj,jk+1)) 
    368                      zdz   =           gdepw_n(ji,jj,jk+1) - pdepw_b(ji,jj,jk  )  
    369                      zdzm1 = MAX(0._wp,pdepw_b(ji,jj,jk  ) - gdepw_n(ji,jj,jk  )) 
     370                     zdzp1 = MAX(0._wp,gdepw(ji,jj,jk+1,Kmm) - pdepw_b(ji,jj,jk+1)) 
     371                     zdz   =           gdepw(ji,jj,jk+1,Kmm) - pdepw_b(ji,jj,jk  )  
     372                     zdzm1 = MAX(0._wp,pdepw_b(ji,jj,jk  ) - gdepw(ji,jj,jk  ,Kmm)) 
    370373                     IF (zdz .LT. 0._wp) THEN  
    371374                        CALL ctl_stop( 'STOP', 'rst_iscpl : unable to compute the interpolation' ) 
    372375                     END IF 
    373                      tsn(ji,jj,jk,jp_tem) = ( zdzp1*tsb(ji,jj,jk+1,jp_tem)      & 
    374                         &                   + zdz  *tsb(ji,jj,jk  ,jp_tem)      & 
    375                         &                   + zdzm1*tsb(ji,jj,jk-1,jp_tem)      )/e3t_n(ji,jj,jk) 
    376                      tsn(ji,jj,jk,jp_sal) = ( zdzp1*tsb(ji,jj,jk+1,jp_sal)      & 
    377                         &                   + zdz  *tsb(ji,jj,jk  ,jp_sal)      & 
    378                         &                   + zdzm1*tsb(ji,jj,jk-1,jp_sal)      )/e3t_n(ji,jj,jk) 
     376                     ts(ji,jj,jk,jp_tem,Kmm) = ( zdzp1*ts(ji,jj,jk+1,jp_tem,Kbb)      & 
     377                        &                   + zdz  *ts(ji,jj,jk  ,jp_tem,Kbb)      & 
     378                        &                   + zdzm1*ts(ji,jj,jk-1,jp_tem,Kbb)      )/e3t(ji,jj,jk,Kmm) 
     379                     ts(ji,jj,jk,jp_sal,Kmm) = ( zdzp1*ts(ji,jj,jk+1,jp_sal,Kbb)      & 
     380                        &                   + zdz  *ts(ji,jj,jk  ,jp_sal,Kbb)      & 
     381                        &                   + zdzm1*ts(ji,jj,jk-1,jp_sal,Kbb)      )/e3t(ji,jj,jk,Kmm) 
    379382                  END IF 
    380383               END DO 
     
    386389      ! ----------------------------------------------------------------------------------------- 
    387390      ! case we open a cell but no neigbour cells available to get an estimate of T and S 
    388       WHERE (tmask(:,:,:) == 1._wp .AND. tsn(:,:,:,2) == 0._wp)  
    389          tsn(:,:,:,2) = -99._wp  ! Special value for closed pool (checking purpose in output.init) 
     391      WHERE (tmask(:,:,:) == 1._wp .AND. ts(:,:,:,2,Kmm) == 0._wp)  
     392         ts(:,:,:,2,Kmm) = -99._wp  ! Special value for closed pool (checking purpose in output.init) 
    390393         tmask(:,:,:) = 0._wp    ! set mask to 0 to run 
    391394         umask(:,:,:) = 0._wp 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM/istate.F90

    r10969 r10978  
    8888         !                                    ! ------------------- 
    8989         CALL rst_read( Kbb, Kmm )            ! Read the restart file 
    90          IF (ln_iscpl)       CALL iscpl_stp   ! extrapolate restart to wet and dry 
     90         IF (ln_iscpl)       CALL iscpl_stp( Kbb, Kmm )   ! extrapolate restart to wet and dry 
    9191         CALL day_init                        ! model calendar (using both namelist and restart infos) 
    9292         ! 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/divhor.F90

    r10969 r10978  
    103103      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( hdiv ) !==  ice sheet  ==!   (update hdiv field) 
     105      IF( ln_iscpl .AND. ln_hsb )   CALL iscpl_div( Kmm, hdiv ) !==  ice sheet  ==!   (update hdiv field) 
    106106      ! 
    107107      CALL lbc_lnk( 'divhor', hdiv, 'T', 1. )   !   (no sign change) 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/sshwzv.F90

    r10969 r10978  
    5454CONTAINS 
    5555 
    56    SUBROUTINE ssh_nxt( kt, Kbb, Kmm ) 
     56   SUBROUTINE ssh_nxt( kt, Kbb, Kmm, pssh, Kaa ) 
    5757      !!---------------------------------------------------------------------- 
    5858      !!                ***  ROUTINE ssh_nxt  *** 
    5959      !!                    
    60       !! ** Purpose :   compute the after ssh (ssha) 
     60      !! ** Purpose :   compute the after ssh (ssh(Kaa)) 
    6161      !! 
    6262      !! ** Method  : - Using the incompressibility hypothesis, the ssh increment 
     
    6464      !!      by the time step. 
    6565      !! 
    66       !! ** action  :   ssha, after sea surface height 
     66      !! ** action  :   ssh(:,:,Kaa), after sea surface height 
    6767      !! 
    6868      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling. 
    6969      !!---------------------------------------------------------------------- 
    70       INTEGER, INTENT(in) ::   kt        ! time step 
    71       INTEGER, INTENT(in) ::   Kbb, Kmm  ! time level index 
     70      INTEGER                         , INTENT(in   ) ::   kt             ! time step 
     71      INTEGER                         , INTENT(in   ) ::   Kbb, Kmm, Kaa  ! time level index 
     72      REAL(wp), DIMENSION(jpi,jpj,jpt), INTENT(inout) ::   pssh           ! sea-surface height 
    7273      !  
    7374      INTEGER  ::   jk            ! dummy loop indice 
     
    9293      !                                           !------------------------------! 
    9394      IF(ln_wd_il) THEN 
    94          CALL wad_lmt(sshb, zcoef * (emp_b(:,:) + emp(:,:)), z2dt) 
     95         CALL wad_lmt(pssh(:,:,Kbb), zcoef * (emp_b(:,:) + emp(:,:)), z2dt) 
    9596      ENDIF 
    9697 
     
    99100      zhdiv(:,:) = 0._wp 
    100101      DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports 
    101         zhdiv(:,:) = zhdiv(:,:) + e3t_n(:,:,jk) * hdivn(:,:,jk) 
     102        zhdiv(:,:) = zhdiv(:,:) + e3t(:,:,jk,Kmm) * hdiv(:,:,jk) 
    102103      END DO 
    103104      !                                                ! Sea surface elevation time stepping 
     
    105106      ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations. 
    106107      !  
    107       ssha(:,:) = (  sshb(:,:) - z2dt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * ssmask(:,:) 
     108      pssh(:,:,Kaa) = (  pssh(:,:,Kbb) - z2dt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * ssmask(:,:) 
    108109      ! 
    109110#if defined key_agrif 
     
    113114      IF ( .NOT.ln_dynspg_ts ) THEN 
    114115         IF( ln_bdy ) THEN 
    115             CALL lbc_lnk( 'sshwzv', ssha, 'T', 1. )    ! Not sure that's necessary 
    116             CALL bdy_ssh( ssha )             ! Duplicate sea level across open boundaries 
     116            CALL lbc_lnk( 'sshwzv', pssh(:,:,Kaa), 'T', 1. )    ! Not sure that's necessary 
     117            CALL bdy_ssh( pssh(:,:,Kaa) )             ! Duplicate sea level across open boundaries 
    117118         ENDIF 
    118119      ENDIF 
     
    121122      !                                           !------------------------------! 
    122123      ! 
    123       IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha  - : ', mask1=tmask ) 
     124      IF(ln_ctl)   CALL prt_ctl( tab2d_1=pssh(:,:,Kaa), clinfo1=' pssh(:,:,Kaa)  - : ', mask1=tmask ) 
    124125      ! 
    125126      IF( ln_timing )   CALL timing_stop('ssh_nxt') 
     
    128129 
    129130    
    130    SUBROUTINE wzv( kt ) 
     131   SUBROUTINE wzv( kt, Kbb, Kmm, pww, Kaa ) 
    131132      !!---------------------------------------------------------------------- 
    132133      !!                ***  ROUTINE wzv  *** 
     
    139140      !!        The boundary conditions are w=0 at the bottom (no flux) and. 
    140141      !! 
    141       !! ** action  :   wn      : now vertical velocity 
     142      !! ** action  :   pww      : now vertical velocity 
    142143      !! 
    143144      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling. 
    144145      !!---------------------------------------------------------------------- 
    145       INTEGER, INTENT(in) ::   kt   ! time step 
     146      INTEGER                         , INTENT(in)    ::   kt             ! time step 
     147      INTEGER                         , INTENT(in)    ::   Kbb, Kmm, Kaa  ! time level indices 
     148      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pww            ! now vertical velocity 
    146149      ! 
    147150      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     
    157160         IF(lwp) WRITE(numout,*) '~~~~~ ' 
    158161         ! 
    159          wn(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all) 
     162         pww(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all) 
    160163      ENDIF 
    161164      !                                           !------------------------------! 
     
    179182         CALL lbc_lnk('sshwzv', zhdiv, 'T', 1.)  ! - ML - Perhaps not necessary: not used for horizontal "connexions" 
    180183         !                             ! Is it problematic to have a wrong vertical velocity in boundary cells? 
    181          !                             ! Same question holds for hdivn. Perhaps just for security 
     184         !                             ! Same question holds for hdiv. Perhaps just for security 
    182185         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence 
    183186            ! computation of w 
    184             wn(:,:,jk) = wn(:,:,jk+1) - (  e3t_n(:,:,jk) * hdivn(:,:,jk) + zhdiv(:,:,jk)    & 
    185                &                         + z1_2dt * ( e3t_a(:,:,jk) - e3t_b(:,:,jk) )     ) * tmask(:,:,jk) 
    186          END DO 
    187          !          IF( ln_vvl_layer ) wn(:,:,:) = 0.e0 
     187            pww(:,:,jk) = pww(:,:,jk+1) - (  e3t(:,:,jk,Kmm) * hdiv(:,:,jk) + zhdiv(:,:,jk)    & 
     188               &                         + z1_2dt * ( e3t(:,:,jk,Kaa) - e3t(:,:,jk,Kbb) )     ) * tmask(:,:,jk) 
     189         END DO 
     190         !          IF( ln_vvl_layer ) pww(:,:,:) = 0.e0 
    188191         DEALLOCATE( zhdiv )  
    189192      ELSE   ! z_star and linear free surface cases 
    190193         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence 
    191194            ! computation of w 
    192             wn(:,:,jk) = wn(:,:,jk+1) - (  e3t_n(:,:,jk) * hdivn(:,:,jk)                 & 
    193                &                         + z1_2dt * ( e3t_a(:,:,jk) - e3t_b(:,:,jk) )  ) * tmask(:,:,jk) 
     195            pww(:,:,jk) = pww(:,:,jk+1) - (  e3t(:,:,jk,Kmm) * hdiv(:,:,jk)                 & 
     196               &                         + z1_2dt * ( e3t(:,:,jk,Kaa) - e3t(:,:,jk,Kbb) )  ) * tmask(:,:,jk) 
    194197         END DO 
    195198      ENDIF 
     
    197200      IF( ln_bdy ) THEN 
    198201         DO jk = 1, jpkm1 
    199             wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:) 
     202            pww(:,:,jk) = pww(:,:,jk) * bdytmask(:,:) 
    200203         END DO 
    201204      ENDIF 
     
    203206#if defined key_agrif  
    204207      IF( .NOT. AGRIF_Root() ) THEN  
    205          IF ((nbondi ==  1).OR.(nbondi == 2)) wn(nlci-1 , :     ,:) = 0.e0      ! east  
    206          IF ((nbondi == -1).OR.(nbondi == 2)) wn(2      , :     ,:) = 0.e0      ! west  
    207          IF ((nbondj ==  1).OR.(nbondj == 2)) wn(:      ,nlcj-1 ,:) = 0.e0      ! north  
    208          IF ((nbondj == -1).OR.(nbondj == 2)) wn(:      ,2      ,:) = 0.e0      ! south  
     208         IF ((nbondi ==  1).OR.(nbondi == 2)) pww(nlci-1 , :     ,:) = 0.e0      ! east  
     209         IF ((nbondi == -1).OR.(nbondi == 2)) pww(2      , :     ,:) = 0.e0      ! west  
     210         IF ((nbondj ==  1).OR.(nbondj == 2)) pww(:      ,nlcj-1 ,:) = 0.e0      ! north  
     211         IF ((nbondj == -1).OR.(nbondj == 2)) pww(:      ,2      ,:) = 0.e0      ! south  
    209212      ENDIF  
    210213#endif  
     
    215218 
    216219 
    217    SUBROUTINE ssh_swp( kt ) 
     220   SUBROUTINE ssh_swp( kt, Kbb, Kmm, Kaa ) 
    218221      !!---------------------------------------------------------------------- 
    219222      !!                    ***  ROUTINE ssh_nxt  *** 
     
    221224      !! ** Purpose :   achieve the sea surface  height time stepping by  
    222225      !!              applying Asselin time filter and swapping the arrays 
    223       !!              ssha  already computed in ssh_nxt   
     226      !!              ssh(:,:,Kaa)  already computed in ssh_nxt   
    224227      !! 
    225228      !! ** Method  : - apply Asselin time fiter to now ssh (excluding the forcing 
    226229      !!              from the filter, see Leclair and Madec 2010) and swap : 
    227       !!                sshn = ssha + atfp * ( sshb -2 sshn + ssha ) 
     230      !!                ssh(:,:,Kmm) = ssh(:,:,Kaa) + atfp * ( ssh(:,:,Kbb) -2 ssh(:,:,Kmm) + ssh(:,:,Kaa) ) 
    228231      !!                            - atfp * rdt * ( emp_b - emp ) / rau0 
    229       !!                sshn = ssha 
    230       !! 
    231       !! ** action  : - sshb, sshn   : before & now sea surface height 
     232      !!                ssh(:,:,Kmm) = ssh(:,:,Kaa) 
     233      !! 
     234      !! ** action  : - ssh(:,:,Kbb), ssh(:,:,Kmm)   : before & now sea surface height 
    232235      !!                               ready for the next time step 
    233236      !! 
    234237      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling. 
    235238      !!---------------------------------------------------------------------- 
    236       INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     239      INTEGER, INTENT(in) ::   kt              ! ocean time-step index 
     240      INTEGER, INTENT(in) ::   Kbb, Kmm, Kaa   ! ocean time-step index 
    237241      ! 
    238242      REAL(wp) ::   zcoef   ! local scalar 
     
    248252      !              !==  Euler time-stepping: no filter, just swap  ==! 
    249253      IF ( neuler == 0 .AND. kt == nit000 ) THEN 
    250          sshn(:,:) = ssha(:,:)                              ! now    <-- after  (before already = now) 
     254         ssh(:,:,Kmm) = ssh(:,:,Kaa)                              ! now    <-- after  (before already = now) 
    251255         ! 
    252256      ELSE           !==  Leap-Frog time-stepping: Asselin filter + swap  ==! 
    253257         !                                                  ! before <-- now filtered 
    254          sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) ) 
     258         ssh(:,:,Kbb) = ssh(:,:,Kmm) + atfp * ( ssh(:,:,Kbb) - 2 * ssh(:,:,Kmm) + ssh(:,:,Kaa) ) 
    255259         IF( .NOT.ln_linssh ) THEN                          ! before <-- with forcing removed 
    256260            zcoef = atfp * rdt * r1_rau0 
    257             sshb(:,:) = sshb(:,:) - zcoef * (     emp_b(:,:) - emp   (:,:)   & 
     261            ssh(:,:,Kbb) = ssh(:,:,Kbb) - zcoef * (     emp_b(:,:) - emp   (:,:)   & 
    258262               &                             -    rnf_b(:,:) + rnf   (:,:)   & 
    259263               &                             + fwfisf_b(:,:) - fwfisf(:,:)   ) * ssmask(:,:) 
    260264         ENDIF 
    261          sshn(:,:) = ssha(:,:)                              ! now <-- after 
    262       ENDIF 
    263       ! 
    264       IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb  - : ', mask1=tmask ) 
     265         ssh(:,:,Kmm) = ssh(:,:,Kaa)                              ! now <-- after 
     266      ENDIF 
     267      ! 
     268      IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssh(:,:,Kbb), clinfo1=' ssh(:,:,Kbb)  - : ', mask1=tmask ) 
    265269      ! 
    266270      IF( ln_timing )   CALL timing_stop('ssh_swp') 
     
    268272   END SUBROUTINE ssh_swp 
    269273 
    270    SUBROUTINE wAimp( kt ) 
     274   SUBROUTINE wAimp( kt, Kmm ) 
    271275      !!---------------------------------------------------------------------- 
    272276      !!                ***  ROUTINE wAimp  *** 
     
    277281      !! ** Method  : -  
    278282      !! 
    279       !! ** action  :   wn      : now vertical velocity (to be handled explicitly) 
     283      !! ** action  :   ww      : now vertical velocity (to be handled explicitly) 
    280284      !!            :   wi      : now vertical velocity (for implicit treatment) 
    281285      !! 
     
    283287      !!---------------------------------------------------------------------- 
    284288      INTEGER, INTENT(in) ::   kt   ! time step 
     289      INTEGER, INTENT(in) ::   Kmm  ! time level index 
    285290      ! 
    286291      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     
    305310         DO jj = 2, jpjm1 
    306311            DO ji = 2, fs_jpim1   ! vector opt. 
    307                z1_e3w = 1._wp / e3w_n(ji,jj,jk) 
    308                Cu_adv(ji,jj,jk) = r2dt * ( ( MAX( wn(ji,jj,jk) , 0._wp ) - MIN( wn(ji,jj,jk+1) , 0._wp ) )    & 
    309                   &                      + ( MAX( e2u(ji  ,jj)*e3uw_n(ji  ,jj,jk)*un(ji  ,jj,jk), 0._wp ) -   & 
    310                   &                          MIN( e2u(ji-1,jj)*e3uw_n(ji-1,jj,jk)*un(ji-1,jj,jk), 0._wp ) )   & 
     312               z1_e3w = 1._wp / e3w(ji,jj,jk,Kmm) 
     313               Cu_adv(ji,jj,jk) = r2dt * ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) )    & 
     314                  &                      + ( MAX( e2u(ji  ,jj)*e3uw(ji  ,jj,jk,Kmm)*uu(ji  ,jj,jk,Kmm), 0._wp ) -   & 
     315                  &                          MIN( e2u(ji-1,jj)*e3uw(ji-1,jj,jk,Kmm)*uu(ji-1,jj,jk,Kmm), 0._wp ) )   & 
    311316                  &                        * r1_e1e2t(ji,jj)                                                  & 
    312                   &                      + ( MAX( e1v(ji,jj  )*e3vw_n(ji,jj  ,jk)*vn(ji,jj  ,jk), 0._wp ) -   & 
    313                   &                          MIN( e1v(ji,jj-1)*e3vw_n(ji,jj-1,jk)*vn(ji,jj-1,jk), 0._wp ) )   & 
     317                  &                      + ( MAX( e1v(ji,jj  )*e3vw(ji,jj  ,jk,Kmm)*vv(ji,jj  ,jk,Kmm), 0._wp ) -   & 
     318                  &                          MIN( e1v(ji,jj-1)*e3vw(ji,jj-1,jk,Kmm)*vv(ji,jj-1,jk,Kmm), 0._wp ) )   & 
    314319                  &                        * r1_e1e2t(ji,jj)                                                  & 
    315320                  &                      ) * z1_e3w 
     
    338343                  zcff = MIN(1._wp, zcff) 
    339344                  ! 
    340                   wi(ji,jj,jk) =           zcff   * wn(ji,jj,jk) 
    341                   wn(ji,jj,jk) = ( 1._wp - zcff ) * wn(ji,jj,jk) 
     345                  wi(ji,jj,jk) =           zcff   * ww(ji,jj,jk) 
     346                  ww(ji,jj,jk) = ( 1._wp - zcff ) * ww(ji,jj,jk) 
    342347                  ! 
    343348                  Cu_adv(ji,jj,jk) = zcff               ! Reuse array to output coefficient 
     
    351356      CALL iom_put("wimp",wi)  
    352357      CALL iom_put("wi_cff",Cu_adv) 
    353       CALL iom_put("wexp",wn) 
     358      CALL iom_put("wexp",ww) 
    354359      ! 
    355360      IF( ln_timing )   CALL timing_stop('wAimp') 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/nemogcm.F90

    r10970 r10978  
    420420      IF( lk_c1d       )   CALL     c1d_init        ! 1D column configuration 
    421421                           CALL     wad_init        ! Wetting and drying options 
    422                            CALL     dom_init("OPA") ! Domain 
    423       IF( ln_crs       )   CALL     crs_init( Nnn ) ! coarsened grid: domain initialization  
     422                           CALL     dom_init( Nbb, Nnn, Naa, "OPA") ! Domain 
     423      IF( ln_crs       )   CALL     crs_init(      Nnn )      ! coarsened grid: domain initialization  
    424424      IF( ln_ctl       )   CALL prt_ctl_init        ! Print control 
    425425       
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/step.F90

    r10970 r10978  
    165165      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    166166 
    167                             CALL ssh_nxt       ( kstp, Nbb, Nnn )  ! after ssh (includes call to div_hor) 
    168       IF( .NOT.ln_linssh )  CALL dom_vvl_sf_nxt( kstp )  ! after vertical scale factors  
    169                             CALL wzv           ( kstp )  ! now cross-level velocity  
    170       IF( ln_zad_Aimp )     CALL wAimp         ( kstp )  ! Adaptive-implicit vertical advection partitioning 
     167                            CALL ssh_nxt       ( kstp, Nbb, Nnn, ssh, Naa )    ! after ssh (includes call to div_hor) 
     168      IF( .NOT.ln_linssh )  CALL dom_vvl_sf_nxt( kstp, Nbb, Nnn,      Naa )    ! after vertical scale factors  
     169                            CALL wzv           ( kstp, Nbb, Nnn, ww,  Naa )    ! now cross-level velocity  
     170      IF( ln_zad_Aimp )     CALL wAimp         ( kstp,      Nnn          )  ! Adaptive-implicit vertical advection partitioning 
    171171                            CALL eos    ( tsn, rhd, rhop, gdept_n(:,:,:) )  ! now in situ density for hpg computation 
    172172                             
     
    202202                                                      ! With split-explicit free surface, since now transports have been updated and ssha as well 
    203203      IF( ln_dynspg_ts ) THEN                         ! vertical scale factors and vertical velocity need to be updated 
    204                             CALL div_hor    ( kstp, Nbb, Nnn )    ! Horizontal divergence  (2nd call in time-split case) 
    205          IF(.NOT.ln_linssh) CALL dom_vvl_sf_nxt( kstp, kcall=2 )  ! after vertical scale factors (update depth average component) 
    206                             CALL wzv        ( kstp )              ! now cross-level velocity  
    207          IF( ln_zad_Aimp )  CALL wAimp      ( kstp )  ! Adaptive-implicit vertical advection partitioning 
     204                            CALL div_hor       ( kstp, Nbb, Nnn )    ! Horizontal divergence  (2nd call in time-split case) 
     205         IF(.NOT.ln_linssh) CALL dom_vvl_sf_nxt( kstp, Nbb, Nnn, Naa, kcall=2 )  ! after vertical scale factors (update depth average component) 
     206                            CALL wzv        ( kstp, Nbb, Nnn, ww, Naa )  ! now cross-level velocity  
     207         IF( ln_zad_Aimp )  CALL wAimp      ( kstp,      Nnn )           ! Adaptive-implicit vertical advection partitioning 
    208208      ENDIF 
    209209       
     
    283283                         CALL tra_nxt       ( kstp, Nbb, Nnn, Nrhs, Naa )  ! finalize (bcs) tracer fields at next time step and swap 
    284284                         CALL dyn_nxt       ( kstp, Nbb, Nnn, Naa  )  ! finalize (bcs) velocities at next time step and swap (always called after tra_nxt) 
    285                          CALL ssh_swp       ( kstp )  ! swap of sea surface height 
    286       IF(.NOT.ln_linssh) CALL dom_vvl_sf_swp( kstp )  ! swap of vertical scale factors 
     285                         CALL ssh_swp       ( kstp, Nbb, Nnn, Naa )  ! swap of sea surface height 
     286      IF(.NOT.ln_linssh) CALL dom_vvl_sf_swp( kstp, Nbb, Nnn, Naa )  ! swap of vertical scale factors 
    287287      ! 
    288288      IF( ln_diahsb  )   CALL dia_hsb       ( kstp, Nbb, Nnn )  ! - ML - global conservation diagnostics 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OFF/nemogcm.F90

    r10928 r10978  
    312312                           CALL     eos_init        ! Equation of state 
    313313      IF( lk_c1d       )   CALL     c1d_init        ! 1D column configuration 
    314                            CALL     dom_init("OPA") ! Domain 
     314                           CALL     dom_init( Nbb, Nnn, Naa, "OPA") ! Domain 
    315315      IF( ln_ctl       )   CALL prt_ctl_init        ! Print control 
    316316 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/SAO/nemogcm.F90

    r10874 r10978  
    246246                           CALL phy_cst            ! Physical constants 
    247247                           CALL eos_init           ! Equation of state 
    248                            CALL dom_init('SAO')    ! Domain 
     248                           CALL dom_init( Nbb, Nnn, Naa, 'SAO')    ! Domain 
    249249 
    250250 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/SAS/nemogcm.F90

    r10927 r10978  
    369369                           CALL phy_cst         ! Physical constants 
    370370                           CALL eos_init        ! Equation of seawater 
    371                            CALL dom_init('SAS') ! Domain 
     371                           CALL dom_init( Nbb, Nnn, Naa, 'SAS') ! Domain 
    372372      IF( ln_ctl      )    CALL prt_ctl_init    ! Print control 
    373373       
Note: See TracChangeset for help on using the changeset viewer.