Changeset 11771


Ignore:
Timestamp:
2019-10-23T15:04:25+02:00 (11 months ago)
Author:
acc
Message:

Branch 2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps. Updates to MY_SRC files in CANAL test case (previously untested on this branch).

Location:
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/tests/CANAL/MY_SRC
Files:
1 added
5 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/tests/CANAL/MY_SRC/diawri.F90

    r10425 r11771  
    5757   USE lib_mpp         ! MPP library 
    5858   USE timing          ! preformance summary 
    59    USE diurnal_bulk    ! diurnal warm layer 
    60    USE cool_skin       ! Cool skin 
     59   USE diu_bulk        ! diurnal warm layer 
     60   USE diu_coolskin    ! Cool skin 
    6161 
    6262   IMPLICIT NONE 
     
    9797 
    9898    
    99    SUBROUTINE dia_wri( kt ) 
     99   SUBROUTINE dia_wri( kt, Kmm ) 
    100100      !!--------------------------------------------------------------------- 
    101101      !!                  ***  ROUTINE dia_wri  *** 
     
    107107      !!---------------------------------------------------------------------- 
    108108      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index 
     109      INTEGER, INTENT( in ) ::   Kmm     ! ocean time level index 
    109110      !! 
    110111      INTEGER ::   ji, jj, jk       ! dummy loop indices 
     
    115116      REAL(wp), DIMENSION(jpi,jpj)     ::   z2d   ! 2D workspace 
    116117      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   z3d   ! 3D workspace 
    117       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   bu, bv   ! volume of u- and v-boxes 
    118       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   r1_bt    ! inverse of t-box volume 
    119118      !!---------------------------------------------------------------------- 
    120119      !  
     
    123122      ! Output the initial state and forcings 
    124123      IF( ninist == 1 ) THEN                        
    125          CALL dia_wri_state( 'output.init' ) 
     124         CALL dia_wri_state( Kmm, 'output.init' ) 
    126125         ninist = 0 
    127126      ENDIF 
     
    132131      CALL iom_put("e3v_0", e3v_0(:,:,:) ) 
    133132      ! 
    134       CALL iom_put( "e3t" , e3t_n(:,:,:) ) 
    135       CALL iom_put( "e3u" , e3u_n(:,:,:) ) 
    136       CALL iom_put( "e3v" , e3v_n(:,:,:) ) 
    137       CALL iom_put( "e3w" , e3w_n(:,:,:) ) 
     133      CALL iom_put( "e3t" , e3t(:,:,:,Kmm) ) 
     134      CALL iom_put( "e3u" , e3u(:,:,:,Kmm) ) 
     135      CALL iom_put( "e3v" , e3v(:,:,:,Kmm) ) 
     136      CALL iom_put( "e3w" , e3w(:,:,:,Kmm) ) 
    138137      IF( iom_use("e3tdef") )   & 
    139          CALL iom_put( "e3tdef"  , ( ( e3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 ) 
     138         CALL iom_put( "e3tdef"  , ( ( e3t(:,:,:,Kmm) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 ) 
    140139 
    141140      IF( ll_wd ) THEN 
    142          CALL iom_put( "ssh" , (sshn+ssh_ref)*tmask(:,:,1) )   ! sea surface height (brought back to the reference used for wetting and drying) 
     141         CALL iom_put( "ssh" , (ssh(:,:,Kmm)+ssh_ref)*tmask(:,:,1) )   ! sea surface height (brought back to the reference used for wetting and drying) 
    143142      ELSE 
    144          CALL iom_put( "ssh" , sshn )              ! sea surface height 
     143         CALL iom_put( "ssh" , ssh(:,:,Kmm) )              ! sea surface height 
    145144      ENDIF 
    146145 
    147146      IF( iom_use("wetdep") )   &                  ! wet depth 
    148          CALL iom_put( "wetdep" , ht_0(:,:) + sshn(:,:) ) 
     147         CALL iom_put( "wetdep" , ht_0(:,:) + ssh(:,:,Kmm) ) 
    149148       
    150       CALL iom_put( "toce", tsn(:,:,:,jp_tem) )    ! 3D temperature 
    151       CALL iom_put(  "sst", tsn(:,:,1,jp_tem) )    ! surface temperature 
     149      CALL iom_put( "toce", ts(:,:,:,jp_tem,Kmm) )    ! 3D temperature 
     150      CALL iom_put(  "sst", ts(:,:,1,jp_tem,Kmm) )    ! surface temperature 
    152151      IF ( iom_use("sbt") ) THEN 
    153152         DO jj = 1, jpj 
    154153            DO ji = 1, jpi 
    155154               ikbot = mbkt(ji,jj) 
    156                z2d(ji,jj) = tsn(ji,jj,ikbot,jp_tem) 
     155               z2d(ji,jj) = ts(ji,jj,ikbot,jp_tem,Kmm) 
    157156            END DO 
    158157         END DO 
     
    160159      ENDIF 
    161160       
    162       CALL iom_put( "soce", tsn(:,:,:,jp_sal) )    ! 3D salinity 
    163       CALL iom_put(  "sss", tsn(:,:,1,jp_sal) )    ! surface salinity 
     161      CALL iom_put( "soce", ts(:,:,:,jp_sal,Kmm) )    ! 3D salinity 
     162      CALL iom_put(  "sss", ts(:,:,1,jp_sal,Kmm) )    ! surface salinity 
    164163      IF ( iom_use("sbs") ) THEN 
    165164         DO jj = 1, jpj 
    166165            DO ji = 1, jpi 
    167166               ikbot = mbkt(ji,jj) 
    168                z2d(ji,jj) = tsn(ji,jj,ikbot,jp_sal) 
     167               z2d(ji,jj) = ts(ji,jj,ikbot,jp_sal,Kmm) 
    169168            END DO 
    170169         END DO 
     
    177176         DO jj = 2, jpjm1 
    178177            DO ji = fs_2, fs_jpim1   ! vector opt. 
    179                zztmp2 = (  ( rCdU_bot(ji+1,jj)+rCdU_bot(ji  ,jj) ) * un(ji  ,jj,mbku(ji  ,jj))  )**2   & 
    180                   &   + (  ( rCdU_bot(ji  ,jj)+rCdU_bot(ji-1,jj) ) * un(ji-1,jj,mbku(ji-1,jj))  )**2   & 
    181                   &   + (  ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj  ) ) * vn(ji,jj  ,mbkv(ji,jj  ))  )**2   & 
    182                   &   + (  ( rCdU_bot(ji,jj  )+rCdU_bot(ji,jj-1) ) * vn(ji,jj-1,mbkv(ji,jj-1))  )**2 
     178               zztmp2 = (  ( rCdU_bot(ji+1,jj)+rCdU_bot(ji  ,jj) ) * uu(ji  ,jj,mbku(ji  ,jj),Kmm)  )**2   & 
     179                  &   + (  ( rCdU_bot(ji  ,jj)+rCdU_bot(ji-1,jj) ) * uu(ji-1,jj,mbku(ji-1,jj),Kmm)  )**2   & 
     180                  &   + (  ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj  ) ) * vv(ji,jj  ,mbkv(ji,jj  ),Kmm)  )**2   & 
     181                  &   + (  ( rCdU_bot(ji,jj  )+rCdU_bot(ji,jj-1) ) * vv(ji,jj-1,mbkv(ji,jj-1),Kmm)  )**2 
    183182               z2d(ji,jj) = zztmp * SQRT( zztmp2 ) * tmask(ji,jj,1)  
    184183               ! 
     
    189188      ENDIF 
    190189          
    191       CALL iom_put( "uoce", un(:,:,:) )            ! 3D i-current 
    192       CALL iom_put(  "ssu", un(:,:,1) )            ! surface i-current 
     190      CALL iom_put( "uoce", uu(:,:,:,Kmm) )            ! 3D i-current 
     191      CALL iom_put(  "ssu", uu(:,:,1,Kmm) )            ! surface i-current 
    193192      IF ( iom_use("sbu") ) THEN 
    194193         DO jj = 1, jpj 
    195194            DO ji = 1, jpi 
    196195               ikbot = mbku(ji,jj) 
    197                z2d(ji,jj) = un(ji,jj,ikbot) 
     196               z2d(ji,jj) = uu(ji,jj,ikbot,Kmm) 
    198197            END DO 
    199198         END DO 
     
    201200      ENDIF 
    202201       
    203       CALL iom_put( "voce", vn(:,:,:) )            ! 3D j-current 
    204       CALL iom_put(  "ssv", vn(:,:,1) )            ! surface j-current 
     202      CALL iom_put( "voce", vv(:,:,:,Kmm) )            ! 3D j-current 
     203      CALL iom_put(  "ssv", vv(:,:,1,Kmm) )            ! surface j-current 
    205204      IF ( iom_use("sbv") ) THEN 
    206205         DO jj = 1, jpj 
    207206            DO ji = 1, jpi 
    208207               ikbot = mbkv(ji,jj) 
    209                z2d(ji,jj) = vn(ji,jj,ikbot) 
     208               z2d(ji,jj) = vv(ji,jj,ikbot,Kmm) 
    210209            END DO 
    211210         END DO 
     
    213212      ENDIF 
    214213 
    215       CALL iom_put( "woce", wn )                   ! vertical velocity 
     214      CALL iom_put( "woce", ww )                   ! vertical velocity 
    216215      IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN   ! vertical mass transport & its square value 
    217216         ! Caution: in the VVL case, it only correponds to the baroclinic mass transport. 
    218217         z2d(:,:) = rau0 * e1e2t(:,:) 
    219218         DO jk = 1, jpk 
    220             z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:) 
     219            z3d(:,:,jk) = ww(:,:,jk) * z2d(:,:) 
    221220         END DO 
    222221         CALL iom_put( "w_masstr" , z3d )   
     
    236235            DO jj = 2, jpjm1                                    ! sal gradient 
    237236               DO ji = fs_2, fs_jpim1   ! vector opt. 
    238                   zztmp  = tsn(ji,jj,jk,jp_sal) 
    239                   zztmpx = ( tsn(ji+1,jj,jk,jp_sal) - zztmp ) * r1_e1u(ji,jj) + ( zztmp - tsn(ji-1,jj  ,jk,jp_sal) ) * r1_e1u(ji-1,jj) 
    240                   zztmpy = ( tsn(ji,jj+1,jk,jp_sal) - zztmp ) * r1_e2v(ji,jj) + ( zztmp - tsn(ji  ,jj-1,jk,jp_sal) ) * r1_e2v(ji,jj-1) 
     237                  zztmp  = ts(ji,jj,jk,jp_sal,Kmm) 
     238                  zztmpx = ( ts(ji+1,jj,jk,jp_sal,Kmm) - zztmp ) * r1_e1u(ji,jj) + ( zztmp - ts(ji-1,jj  ,jk,jp_sal,Kmm) ) * r1_e1u(ji-1,jj) 
     239                  zztmpy = ( ts(ji,jj+1,jk,jp_sal,Kmm) - zztmp ) * r1_e2v(ji,jj) + ( zztmp - ts(ji  ,jj-1,jk,jp_sal,Kmm) ) * r1_e2v(ji,jj-1) 
    241240                  z3d(ji,jj,jk) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy )   & 
    242241                     &                 * umask(ji,jj,jk) * umask(ji-1,jj,jk) * vmask(ji,jj,jk) * umask(ji,jj-1,jk) 
     
    253252         DO jj = 2, jpjm1                                    ! sst gradient 
    254253            DO ji = fs_2, fs_jpim1   ! vector opt. 
    255                zztmp  = tsn(ji,jj,1,jp_tem) 
    256                zztmpx = ( tsn(ji+1,jj,1,jp_tem) - zztmp ) * r1_e1u(ji,jj) + ( zztmp - tsn(ji-1,jj  ,1,jp_tem) ) * r1_e1u(ji-1,jj) 
    257                zztmpy = ( tsn(ji,jj+1,1,jp_tem) - zztmp ) * r1_e2v(ji,jj) + ( zztmp - tsn(ji  ,jj-1,1,jp_tem) ) * r1_e2v(ji,jj-1) 
     254               zztmp  = ts(ji,jj,1,jp_tem,Kmm) 
     255               zztmpx = ( ts(ji+1,jj,1,jp_tem,Kmm) - zztmp ) * r1_e1u(ji,jj) + ( zztmp - ts(ji-1,jj  ,1,jp_tem,Kmm) ) * r1_e1u(ji-1,jj) 
     256               zztmpy = ( ts(ji,jj+1,1,jp_tem,Kmm) - zztmp ) * r1_e2v(ji,jj) + ( zztmp - ts(ji  ,jj-1,1,jp_tem,Kmm) ) * r1_e2v(ji,jj-1) 
    258257               z2d(ji,jj) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy )   & 
    259258                  &              * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * umask(ji,jj-1,1) 
     
    272271            DO jj = 1, jpj 
    273272               DO ji = 1, jpi 
    274                   z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk) 
     273                  z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm) * tmask(ji,jj,jk) 
    275274               END DO 
    276275            END DO 
     
    284283            DO jj = 1, jpj 
    285284               DO ji = 1, jpi 
    286                   z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk) 
     285                  z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * tmask(ji,jj,jk) 
    287286               END DO 
    288287            END DO 
     
    296295            DO jj = 1, jpj 
    297296               DO ji = 1, jpi 
    298                   z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) * tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk) 
     297                  z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * tmask(ji,jj,jk) 
    299298               END DO 
    300299            END DO 
     
    306305         z3d(:,:,jpk) = 0._wp  
    307306         DO jk = 1, jpkm1 
    308             DO jj = 2, jpj 
    309                DO ji = 2, jpi 
    310                   zztmpx = 0.5 * ( un(ji-1,jj  ,jk) + un(ji,jj,jk) ) 
    311                   zztmpy = 0.5 * ( vn(ji  ,jj-1,jk) + vn(ji,jj,jk) ) 
    312                   z3d(ji,jj,jk) = 0.5 * ( zztmpx*zztmpx + zztmpy*zztmpy ) 
     307            DO jj = 2, jpjm1 
     308               DO ji = fs_2, fs_jpim1   ! vector opt. 
     309                  zztmp  = 0.25_wp * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
     310                  z3d(ji,jj,jk) = zztmp * (  uu(ji-1,jj,jk,Kmm)**2 * e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm)   & 
     311                     &                     + uu(ji  ,jj,jk,Kmm)**2 * e2u(ji  ,jj) * e3u(ji  ,jj,jk,Kmm)   & 
     312                     &                     + vv(ji,jj-1,jk,Kmm)**2 * e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm)   & 
     313                     &                     + vv(ji,jj  ,jk,Kmm)**2 * e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm)   ) 
    313314               END DO 
    314315            END DO 
     
    326327            DO jj = 2, jpj 
    327328               DO ji = 2, jpi 
    328                   z3d(ji,jj,jk) = 0.25_wp * ( un(ji  ,jj,jk) * un(ji  ,jj,jk) * e1e2u(ji  ,jj) * e3u_n(ji  ,jj,jk)  & 
    329                      &                      + un(ji-1,jj,jk) * un(ji-1,jj,jk) * e1e2u(ji-1,jj) * e3u_n(ji-1,jj,jk)  & 
    330                      &                      + vn(ji,jj  ,jk) * vn(ji,jj  ,jk) * e1e2v(ji,jj  ) * e3v_n(ji,jj  ,jk)  & 
    331                      &                      + vn(ji,jj-1,jk) * vn(ji,jj-1,jk) * e1e2v(ji,jj-1) * e3v_n(ji,jj-1,jk)  )  & 
    332                      &                    * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) * tmask(ji,jj,jk) 
     329                  z3d(ji,jj,jk) = 0.25_wp * ( uu(ji  ,jj,jk,Kmm) * uu(ji  ,jj,jk,Kmm) * e1e2u(ji  ,jj) * e3u(ji  ,jj,jk,Kmm)  & 
     330                     &                      + uu(ji-1,jj,jk,Kmm) * uu(ji-1,jj,jk,Kmm) * e1e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm)  & 
     331                     &                      + vv(ji,jj  ,jk,Kmm) * vv(ji,jj  ,jk,Kmm) * e1e2v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm)  & 
     332                     &                      + vv(ji,jj-1,jk,Kmm) * vv(ji,jj-1,jk,Kmm) * e1e2v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm)  )  & 
     333                     &                    * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk) 
    333334               END DO 
    334335            END DO 
     
    342343            DO jj = 1, jpj 
    343344               DO ji = 1, jpi 
    344                   z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * z3d(ji,jj,jk) * tmask(ji,jj,jk) 
     345                  z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * z3d(ji,jj,jk) * tmask(ji,jj,jk) 
    345346               END DO 
    346347            END DO 
     
    350351      ENDIF 
    351352      ! 
    352       CALL iom_put( "hdiv", hdivn )                  ! Horizontal divergence 
     353      CALL iom_put( "hdiv", hdiv )                  ! Horizontal divergence 
    353354 
    354355      IF ( iom_use("relvor") .OR. iom_use("absvor") .OR. iom_use("potvor") ) THEN 
     
    358359            DO jj = 1, jpjm1 
    359360               DO ji = 1, fs_jpim1   ! vector opt. 
    360                   z3d(ji,jj,jk) = (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    & 
    361                      &              - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) * r1_e1e2f(ji,jj) 
     361                  z3d(ji,jj,jk) = (   e2v(ji+1,jj  ) * vv(ji+1,jj  ,jk,Kmm) - e2v(ji,jj) * vv(ji,jj,jk,Kmm)    & 
     362                     &              - e1u(ji  ,jj+1) * uu(ji  ,jj+1,jk,Kmm) + e1u(ji,jj) * uu(ji,jj,jk,Kmm)  ) * r1_e1e2f(ji,jj) 
    362363               END DO 
    363364            END DO 
     
    378379            DO jj = 1, jpjm1 
    379380               DO ji = 1, fs_jpim1   ! vector opt. 
    380                   ze3  = (  e3t_n(ji,jj+1,jk)*tmask(ji,jj+1,jk) + e3t_n(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   & 
    381                      &    + e3t_n(ji,jj  ,jk)*tmask(ji,jj  ,jk) + e3t_n(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk)  ) 
     381                  ze3  = (  e3t(ji,jj+1,jk,Kmm)*tmask(ji,jj+1,jk) + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   & 
     382                     &    + e3t(ji,jj  ,jk,Kmm)*tmask(ji,jj  ,jk) + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)  ) 
    382383                  IF( ze3 /= 0._wp ) THEN   ;   ze3 = 4._wp / ze3 
    383384                  ELSE                      ;   ze3 = 0._wp 
     
    397398         z2d(:,:) = 0.e0 
    398399         DO jk = 1, jpkm1 
    399             z3d(:,:,jk) = rau0 * un(:,:,jk) * e2u(:,:) * e3u_n(:,:,jk) * umask(:,:,jk) 
     400            z3d(:,:,jk) = rau0 * uu(:,:,jk,Kmm) * e2u(:,:) * e3u(:,:,jk,Kmm) * umask(:,:,jk) 
    400401            z2d(:,:) = z2d(:,:) + z3d(:,:,jk) 
    401402         END DO 
     
    409410            DO jj = 2, jpjm1 
    410411               DO ji = fs_2, fs_jpim1   ! vector opt. 
    411                   z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) ) 
     412                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_tem,Kmm) + ts(ji+1,jj,jk,jp_tem,Kmm) ) 
    412413               END DO 
    413414            END DO 
     
    422423            DO jj = 2, jpjm1 
    423424               DO ji = fs_2, fs_jpim1   ! vector opt. 
    424                   z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) ) 
     425                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_sal,Kmm) + ts(ji+1,jj,jk,jp_sal,Kmm) ) 
    425426               END DO 
    426427            END DO 
     
    434435         z3d(:,:,jpk) = 0.e0 
    435436         DO jk = 1, jpkm1 
    436             z3d(:,:,jk) = rau0 * vn(:,:,jk) * e1v(:,:) * e3v_n(:,:,jk) * vmask(:,:,jk) 
     437            z3d(:,:,jk) = rau0 * vv(:,:,jk,Kmm) * e1v(:,:) * e3v(:,:,jk,Kmm) * vmask(:,:,jk) 
    437438         END DO 
    438439         CALL iom_put( "v_masstr", z3d )              ! mass transport in j-direction 
     
    444445            DO jj = 2, jpjm1 
    445446               DO ji = fs_2, fs_jpim1   ! vector opt. 
    446                   z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji,jj+1,jk,jp_tem) ) 
     447                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_tem,Kmm) + ts(ji,jj+1,jk,jp_tem,Kmm) ) 
    447448               END DO 
    448449            END DO 
     
    457458            DO jj = 2, jpjm1 
    458459               DO ji = fs_2, fs_jpim1   ! vector opt. 
    459                   z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji,jj+1,jk,jp_sal) ) 
     460                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_sal,Kmm) + ts(ji,jj+1,jk,jp_sal,Kmm) ) 
    460461               END DO 
    461462            END DO 
     
    470471            DO jj = 2, jpjm1 
    471472               DO ji = fs_2, fs_jpim1   ! vector opt. 
    472                   z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) *  tsn(ji,jj,jk,jp_tem) 
     473                  z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) *  ts(ji,jj,jk,jp_tem,Kmm) 
    473474               END DO 
    474475            END DO 
     
    482483            DO jj = 2, jpjm1 
    483484               DO ji = fs_2, fs_jpim1   ! vector opt. 
    484                   z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) 
     485                  z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) 
    485486               END DO 
    486487            END DO 
     
    493494      ! 
    494495 
    495       IF (ln_diatmb)   CALL dia_tmb                   ! tmb values  
     496      IF (ln_diatmb)   CALL dia_tmb( Kmm )            ! tmb values  
    496497           
    497       IF (ln_dia25h)   CALL dia_25h( kt )             ! 25h averaging 
     498      IF (ln_dia25h)   CALL dia_25h( kt, Kmm )        ! 25h averaging 
    498499 
    499500      IF( ln_timing )   CALL timing_stop('dia_wri') 
     
    521522 
    522523    
    523    SUBROUTINE dia_wri( kt ) 
     524   SUBROUTINE dia_wri( kt, Kmm ) 
    524525      !!--------------------------------------------------------------------- 
    525526      !!                  ***  ROUTINE dia_wri  *** 
     
    534535      !!---------------------------------------------------------------------- 
    535536      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index 
     537      INTEGER, INTENT( in ) ::   Kmm  ! ocean time level index 
    536538      ! 
    537539      LOGICAL ::   ll_print = .FALSE.                        ! =T print and flush numout 
     
    551553      ! 
    552554      IF( ninist == 1 ) THEN     !==  Output the initial state and forcings  ==! 
    553          CALL dia_wri_state( 'output.init' ) 
     555         CALL dia_wri_state( Kmm, 'output.init' ) 
    554556         ninist = 0 
    555557      ENDIF 
     
    688690            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout ) 
    689691         IF(  .NOT.ln_linssh  ) THEN 
    690             CALL histdef( nid_T, "vovvle3t", "Level thickness"                    , "m"      ,&  ! e3t_n 
     692            CALL histdef( nid_T, "vovvle3t", "Level thickness"                    , "m"      ,&  ! e3t(:,:,:,Kmm) 
    691693            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout ) 
    692             CALL histdef( nid_T, "vovvldep", "T point depth"                      , "m"      ,&  ! e3t_n 
     694            CALL histdef( nid_T, "vovvldep", "T point depth"                      , "m"      ,&  ! e3t(:,:,:,Kmm) 
    693695            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout ) 
    694             CALL histdef( nid_T, "vovvldef", "Squared level deformation"          , "%^2"    ,&  ! e3t_n 
     696            CALL histdef( nid_T, "vovvldef", "Squared level deformation"          , "%^2"    ,&  ! e3t(:,:,:,Kmm) 
    695697            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout ) 
    696698         ENDIF 
     
    709711            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) 
    710712         IF(  ln_linssh  ) THEN 
    711             CALL histdef( nid_T, "sosst_cd", "Concentration/Dilution term on temperature"     &  ! emp * tsn(:,:,1,jp_tem) 
     713            CALL histdef( nid_T, "sosst_cd", "Concentration/Dilution term on temperature"     &  ! emp * ts(:,:,1,jp_tem,Kmm) 
    712714            &                                                                  , "KgC/m2/s",  &  ! sosst_cd 
    713715            &             jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) 
    714             CALL histdef( nid_T, "sosss_cd", "Concentration/Dilution term on salinity"        &  ! emp * tsn(:,:,1,jp_sal) 
     716            CALL histdef( nid_T, "sosss_cd", "Concentration/Dilution term on salinity"        &  ! emp * ts(:,:,1,jp_sal,Kmm) 
    715717            &                                                                  , "KgPSU/m2/s",&  ! sosss_cd 
    716718            &             jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) 
     
    797799 
    798800         !                                                                                      !!! nid_U : 3D 
    799          CALL histdef( nid_U, "vozocrtx", "Zonal Current"                      , "m/s"    ,   &  ! un 
     801         CALL histdef( nid_U, "vozocrtx", "Zonal Current"                      , "m/s"    ,   &  ! uu(:,:,:,Kmm) 
    800802            &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout ) 
    801803         IF( ln_wave .AND. ln_sdw) THEN 
     
    810812 
    811813         !                                                                                      !!! nid_V : 3D 
    812          CALL histdef( nid_V, "vomecrty", "Meridional Current"                 , "m/s"    ,   &  ! vn 
     814         CALL histdef( nid_V, "vomecrty", "Meridional Current"                 , "m/s"    ,   &  ! vv(:,:,:,Kmm) 
    813815            &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout ) 
    814816         IF( ln_wave .AND. ln_sdw) THEN 
     
    823825 
    824826         !                                                                                      !!! nid_W : 3D 
    825          CALL histdef( nid_W, "vovecrtz", "Vertical Velocity"                  , "m/s"    ,   &  ! wn 
     827         CALL histdef( nid_W, "vovecrtz", "Vertical Velocity"                  , "m/s"    ,   &  ! ww 
    826828            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout ) 
    827829         CALL histdef( nid_W, "votkeavt", "Vertical Eddy Diffusivity"          , "m2/s"   ,   &  ! avt 
     
    861863 
    862864      IF( .NOT.ln_linssh ) THEN 
    863          CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) * e3t_n(:,:,:) , ndim_T , ndex_T  )   ! heat content 
    864          CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) * e3t_n(:,:,:) , ndim_T , ndex_T  )   ! salt content 
    865          CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) * e3t_n(:,:,1) , ndim_hT, ndex_hT )   ! sea surface heat content 
    866          CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) * e3t_n(:,:,1) , ndim_hT, ndex_hT )   ! sea surface salinity content 
     865         CALL histwrite( nid_T, "votemper", it, ts(:,:,:,jp_tem,Kmm) * e3t(:,:,:,Kmm) , ndim_T , ndex_T  )   ! heat content 
     866         CALL histwrite( nid_T, "vosaline", it, ts(:,:,:,jp_sal,Kmm) * e3t(:,:,:,Kmm) , ndim_T , ndex_T  )   ! salt content 
     867         CALL histwrite( nid_T, "sosstsst", it, ts(:,:,1,jp_tem,Kmm) * e3t(:,:,1,Kmm) , ndim_hT, ndex_hT )   ! sea surface heat content 
     868         CALL histwrite( nid_T, "sosaline", it, ts(:,:,1,jp_sal,Kmm) * e3t(:,:,1,Kmm) , ndim_hT, ndex_hT )   ! sea surface salinity content 
    867869      ELSE 
    868          CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) , ndim_T , ndex_T  )   ! temperature 
    869          CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) , ndim_T , ndex_T  )   ! salinity 
    870          CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) , ndim_hT, ndex_hT )   ! sea surface temperature 
    871          CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) , ndim_hT, ndex_hT )   ! sea surface salinity 
     870         CALL histwrite( nid_T, "votemper", it, ts(:,:,:,jp_tem,Kmm) , ndim_T , ndex_T  )   ! temperature 
     871         CALL histwrite( nid_T, "vosaline", it, ts(:,:,:,jp_sal,Kmm) , ndim_T , ndex_T  )   ! salinity 
     872         CALL histwrite( nid_T, "sosstsst", it, ts(:,:,1,jp_tem,Kmm) , ndim_hT, ndex_hT )   ! sea surface temperature 
     873         CALL histwrite( nid_T, "sosaline", it, ts(:,:,1,jp_sal,Kmm) , ndim_hT, ndex_hT )   ! sea surface salinity 
    872874      ENDIF 
    873875      IF( .NOT.ln_linssh ) THEN 
    874          zw3d(:,:,:) = ( ( e3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 
    875          CALL histwrite( nid_T, "vovvle3t", it, e3t_n (:,:,:) , ndim_T , ndex_T  )   ! level thickness 
    876          CALL histwrite( nid_T, "vovvldep", it, gdept_n(:,:,:) , ndim_T , ndex_T  )   ! t-point depth 
     876         zw3d(:,:,:) = ( ( e3t(:,:,:,Kmm) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 
     877         CALL histwrite( nid_T, "vovvle3t", it, e3t (:,:,:,Kmm) , ndim_T , ndex_T  )   ! level thickness 
     878         CALL histwrite( nid_T, "vovvldep", it, gdept(:,:,:,Kmm) , ndim_T , ndex_T  )   ! t-point depth 
    877879         CALL histwrite( nid_T, "vovvldef", it, zw3d             , ndim_T , ndex_T  )   ! level thickness deformation 
    878880      ENDIF 
    879       CALL histwrite( nid_T, "sossheig", it, sshn          , ndim_hT, ndex_hT )   ! sea surface height 
     881      CALL histwrite( nid_T, "sossheig", it, ssh(:,:,Kmm)          , ndim_hT, ndex_hT )   ! sea surface height 
    880882      CALL histwrite( nid_T, "sowaflup", it, ( emp-rnf )   , ndim_hT, ndex_hT )   ! upward water flux 
    881883      CALL histwrite( nid_T, "sorunoff", it, rnf           , ndim_hT, ndex_hT )   ! river runoffs 
     
    884886                                                                                  ! in linear free surface case) 
    885887      IF( ln_linssh ) THEN 
    886          zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_tem) 
     888         zw2d(:,:) = emp (:,:) * ts(:,:,1,jp_tem,Kmm) 
    887889         CALL histwrite( nid_T, "sosst_cd", it, zw2d, ndim_hT, ndex_hT )          ! c/d term on sst 
    888          zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_sal) 
     890         zw2d(:,:) = emp (:,:) * ts(:,:,1,jp_sal,Kmm) 
    889891         CALL histwrite( nid_T, "sosss_cd", it, zw2d, ndim_hT, ndex_hT )          ! c/d term on sss 
    890892      ENDIF 
     
    922924         CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping 
    923925         CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping 
    924          IF( ln_ssr ) zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1) 
     926         IF( ln_ssr ) zw2d(:,:) = erp(:,:) * ts(:,:,1,jp_sal,Kmm) * tmask(:,:,1) 
    925927         CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping 
    926928      ENDIF 
     
    928930         CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping 
    929931         CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping 
    930          IF( ln_ssr ) zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1) 
     932         IF( ln_ssr ) zw2d(:,:) = erp(:,:) * ts(:,:,1,jp_sal,Kmm) * tmask(:,:,1) 
    931933         CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping 
    932934      ENDIF 
     
    941943#endif 
    942944 
    943       CALL histwrite( nid_U, "vozocrtx", it, un            , ndim_U , ndex_U )    ! i-current 
     945      CALL histwrite( nid_U, "vozocrtx", it, uu(:,:,:,Kmm)            , ndim_U , ndex_U )    ! i-current 
    944946      CALL histwrite( nid_U, "sozotaux", it, utau          , ndim_hU, ndex_hU )   ! i-wind stress 
    945947 
    946       CALL histwrite( nid_V, "vomecrty", it, vn            , ndim_V , ndex_V  )   ! j-current 
     948      CALL histwrite( nid_V, "vomecrty", it, vv(:,:,:,Kmm)            , ndim_V , ndex_V  )   ! j-current 
    947949      CALL histwrite( nid_V, "sometauy", it, vtau          , ndim_hV, ndex_hV )   ! j-wind stress 
    948950 
    949       CALL histwrite( nid_W, "vovecrtz", it, wn             , ndim_T, ndex_T )    ! vert. current 
     951      CALL histwrite( nid_W, "vovecrtz", it, ww             , ndim_T, ndex_T )    ! vert. current 
    950952      CALL histwrite( nid_W, "votkeavt", it, avt            , ndim_T, ndex_T )    ! T vert. eddy diff. coef. 
    951953      CALL histwrite( nid_W, "votkeavm", it, avm            , ndim_T, ndex_T )    ! T vert. eddy visc. coef. 
     
    974976#endif 
    975977 
    976    SUBROUTINE dia_wri_state( cdfile_name ) 
     978   SUBROUTINE dia_wri_state( Kmm, cdfile_name ) 
    977979      !!--------------------------------------------------------------------- 
    978980      !!                 ***  ROUTINE dia_wri_state  *** 
     
    987989      !!      File 'output.abort.nc' is created in case of abnormal job end 
    988990      !!---------------------------------------------------------------------- 
     991      INTEGER           , INTENT( in ) ::   Kmm              ! time level index 
    989992      CHARACTER (len=* ), INTENT( in ) ::   cdfile_name      ! name of the file created 
    990993      !! 
     
    10031006#endif 
    10041007 
    1005       CALL iom_rstput( 0, 0, inum, 'votemper', tsn(:,:,:,jp_tem) )    ! now temperature 
    1006       CALL iom_rstput( 0, 0, inum, 'vosaline', tsn(:,:,:,jp_sal) )    ! now salinity 
    1007       CALL iom_rstput( 0, 0, inum, 'sossheig', sshn              )    ! sea surface height 
    1008       CALL iom_rstput( 0, 0, inum, 'vozocrtx', un                )    ! now i-velocity 
    1009       CALL iom_rstput( 0, 0, inum, 'vomecrty', vn                )    ! now j-velocity 
    1010       CALL iom_rstput( 0, 0, inum, 'vovecrtz', wn                )    ! now k-velocity 
     1008      CALL iom_rstput( 0, 0, inum, 'votemper', ts(:,:,:,jp_tem,Kmm) )    ! now temperature 
     1009      CALL iom_rstput( 0, 0, inum, 'vosaline', ts(:,:,:,jp_sal,Kmm) )    ! now salinity 
     1010      CALL iom_rstput( 0, 0, inum, 'sossheig', ssh(:,:,Kmm)              )    ! sea surface height 
     1011      CALL iom_rstput( 0, 0, inum, 'vozocrtx', uu(:,:,:,Kmm)                )    ! now i-velocity 
     1012      CALL iom_rstput( 0, 0, inum, 'vomecrty', vv(:,:,:,Kmm)                )    ! now j-velocity 
     1013      CALL iom_rstput( 0, 0, inum, 'vovecrtz', ww                )    ! now k-velocity 
    10111014      IF( ALLOCATED(ahtu) ) THEN 
    10121015         CALL iom_rstput( 0, 0, inum,  'ahtu', ahtu              )    ! aht at u-point 
     
    10241027      CALL iom_rstput( 0, 0, inum, 'sometauy', vtau              )    ! j-wind stress 
    10251028      IF(  .NOT.ln_linssh  ) THEN              
    1026          CALL iom_rstput( 0, 0, inum, 'vovvldep', gdept_n        )    !  T-cell depth  
    1027          CALL iom_rstput( 0, 0, inum, 'vovvle3t', e3t_n          )    !  T-cell thickness   
     1029         CALL iom_rstput( 0, 0, inum, 'vovvldep', gdept(:,:,:,Kmm)        )    !  T-cell depth  
     1030         CALL iom_rstput( 0, 0, inum, 'vovvle3t', e3t(:,:,:,Kmm)          )    !  T-cell thickness   
    10281031      END IF 
    10291032      IF( ln_wave .AND. ln_sdw ) THEN 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/tests/CANAL/MY_SRC/domvvl.F90

    r10425 r11771  
    88   !!            3.3  !  2011-10  (M. Leclair) totally rewrote domvvl: vvl option includes z_star and z_tilde coordinates 
    99   !!            3.6  !  2014-11  (P. Mathiot) add ice shelf capability 
     10   !!            4.1  !  2019-08  (A. Coward, D. Storkey) rename dom_vvl_sf_swp -> dom_vvl_sf_update for new timestepping 
    1011   !!---------------------------------------------------------------------- 
    1112 
     
    1314   !!   dom_vvl_init     : define initial vertical scale factors, depths and column thickness 
    1415   !!   dom_vvl_sf_nxt   : Compute next vertical scale factors 
    15    !!   dom_vvl_sf_swp   : Swap vertical scale factors and update the vertical grid 
     16   !!   dom_vvl_sf_update   : Swap vertical scale factors and update the vertical grid 
    1617   !!   dom_vvl_interpol : Interpolate vertical scale factors from one grid point to another 
    1718   !!   dom_vvl_rst      : read/write restart file 
     
    3738   PUBLIC  dom_vvl_init       ! called by domain.F90 
    3839   PUBLIC  dom_vvl_sf_nxt     ! called by step.F90 
    39    PUBLIC  dom_vvl_sf_swp     ! called by step.F90 
     40   PUBLIC  dom_vvl_sf_update  ! called by step.F90 
    4041   PUBLIC  dom_vvl_interpol   ! called by dynnxt.F90 
    4142 
     
    9394 
    9495 
    95    SUBROUTINE dom_vvl_init 
     96   SUBROUTINE dom_vvl_init( Kbb, Kmm, Kaa ) 
    9697      !!---------------------------------------------------------------------- 
    9798      !!                ***  ROUTINE dom_vvl_init  *** 
     
    104105      !! 
    105106      !! ** 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 
     107      !!              - Regrid: e3[u/v](:,:,:,Kmm) 
     108      !!                        e3[u/v](:,:,:,Kmm)        
     109      !!                        e3w(:,:,:,Kmm)            
     110      !!                        e3[u/v]w_b 
     111      !!                        e3[u/v]w_n       
     112      !!                        gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm) and gde3w 
    112113      !!              - h(t/u/v)_0 
    113114      !!              - frq_rst_e3t and frq_rst_hdv 
     
    115116      !! Reference  : Leclair, M., and G. Madec, 2011, Ocean Modelling. 
    116117      !!---------------------------------------------------------------------- 
     118      INTEGER, INTENT(in) :: Kbb, Kmm, Kaa 
     119      ! 
    117120      INTEGER ::   ji, jj, jk 
    118121      INTEGER ::   ii0, ii1, ij0, ij1 
     
    130133      ! 
    131134      !                    ! 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 
     135      CALL dom_vvl_rst( nit000, Kbb, Kmm, 'READ' ) 
     136      e3t(:,:,jpk,Kaa) = e3t_0(:,:,jpk)  ! last level always inside the sea floor set one for all 
    134137      ! 
    135138      !                    !== Set of all other vertical scale factors  ==!  (now and before) 
    136139      !                                ! 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 
     140      CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3u(:,:,:,Kbb), 'U' )    ! from T to U 
     141      CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' ) 
     142      CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3v(:,:,:,Kbb), 'V' )    ! from T to V  
     143      CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' ) 
     144      CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' )    ! from U to F 
    142145      !                                ! 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' ) 
     146      CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W'  )  ! from T to W 
     147      CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3w (:,:,:,Kbb), 'W'  ) 
     148      CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' )  ! from U to UW 
     149      CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' ) 
     150      CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' )  ! from V to UW 
     151      CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' ) 
    149152 
    150153      ! 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(:,:,:) 
     154      e3t(:,:,:,Kaa) = e3t(:,:,:,Kmm) 
     155      e3u(:,:,:,Kaa) = e3u(:,:,:,Kmm) 
     156      e3v(:,:,:,Kaa) = e3v(:,:,:,Kmm) 
    154157      ! 
    155158      !                    !==  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 
     159      gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm)       ! reference to the ocean surface (used for MLD and light penetration) 
     160      gdepw(:,:,1,Kmm) = 0.0_wp 
     161      gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm)  ! reference to a common level z=0 for hpg 
     162      gdept(:,:,1,Kbb) = 0.5_wp * e3w(:,:,1,Kbb) 
     163      gdepw(:,:,1,Kbb) = 0.0_wp 
    161164      DO jk = 2, jpk                               ! vertical sum 
    162165         DO jj = 1,jpj 
     
    165168               !                             ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) 
    166169               !                             ! 0.5 where jk = mikt      
    167 !!gm ???????   BUG ?  gdept_n as well as gde3w_n  does not include the thickness of ISF ?? 
     170!!gm ???????   BUG ?  gdept(:,:,:,Kmm) as well as gde3w  does not include the thickness of ISF ?? 
    168171               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))  
     172               gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 
     173               gdept(ji,jj,jk,Kmm) =      zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm))  & 
     174                  &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm))  
     175               gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm) 
     176               gdepw(ji,jj,jk,Kbb) = gdepw(ji,jj,jk-1,Kbb) + e3t(ji,jj,jk-1,Kbb) 
     177               gdept(ji,jj,jk,Kbb) =      zcoef  * ( gdepw(ji,jj,jk  ,Kbb) + 0.5 * e3w(ji,jj,jk,Kbb))  & 
     178                  &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kbb) +       e3w(ji,jj,jk,Kbb))  
    176179            END DO 
    177180         END DO 
     
    179182      ! 
    180183      !                    !==  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) 
     184      ht(:,:) = e3t(:,:,1,Kmm) * tmask(:,:,1)   !!gm  BUG  :  this should be 1/2 * e3w(k=1) .... 
     185      hu(:,:,Kbb) = e3u(:,:,1,Kbb) * umask(:,:,1) 
     186      hu(:,:,Kmm) = e3u(:,:,1,Kmm) * umask(:,:,1) 
     187      hv(:,:,Kbb) = e3v(:,:,1,Kbb) * vmask(:,:,1) 
     188      hv(:,:,Kmm) = e3v(:,:,1,Kmm) * vmask(:,:,1) 
    186189      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) 
     190         ht(:,:) = ht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk) 
     191         hu(:,:,Kbb) = hu(:,:,Kbb) + e3u(:,:,jk,Kbb) * umask(:,:,jk) 
     192         hu(:,:,Kmm) = hu(:,:,Kmm) + e3u(:,:,jk,Kmm) * umask(:,:,jk) 
     193         hv(:,:,Kbb) = hv(:,:,Kbb) + e3v(:,:,jk,Kbb) * vmask(:,:,jk) 
     194         hv(:,:,Kmm) = hv(:,:,Kmm) + e3v(:,:,jk,Kmm) * vmask(:,:,jk) 
    192195      END DO 
    193196      ! 
    194197      !                    !==  inverse of water column thickness   ==!   (u- and v- points) 
    195       r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) )    ! _i mask due to ISF 
    196       r1_hu_n(:,:) = ssumask(:,:) / ( hu_n(:,:) + 1._wp - ssumask(:,:) ) 
    197       r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) ) 
    198       r1_hv_n(:,:) = ssvmask(:,:) / ( hv_n(:,:) + 1._wp - ssvmask(:,:) ) 
     198      r1_hu(:,:,Kbb) = ssumask(:,:) / ( hu(:,:,Kbb) + 1._wp - ssumask(:,:) )    ! _i mask due to ISF 
     199      r1_hu(:,:,Kmm) = ssumask(:,:) / ( hu(:,:,Kmm) + 1._wp - ssumask(:,:) ) 
     200      r1_hv(:,:,Kbb) = ssvmask(:,:) / ( hv(:,:,Kbb) + 1._wp - ssvmask(:,:) ) 
     201      r1_hv(:,:,Kmm) = ssvmask(:,:) / ( hv(:,:,Kmm) + 1._wp - ssvmask(:,:) ) 
    199202 
    200203      !                    !==   z_tilde coordinate case  ==!   (Restoring frequencies) 
     
    266269 
    267270 
    268    SUBROUTINE dom_vvl_sf_nxt( kt, kcall )  
     271   SUBROUTINE dom_vvl_sf_nxt( kt, Kbb, Kmm, Kaa, kcall )  
    269272      !!---------------------------------------------------------------------- 
    270273      !!                ***  ROUTINE dom_vvl_sf_nxt  *** 
     
    288291      !! Reference  : Leclair, M., and Madec, G. 2011, Ocean Modelling. 
    289292      !!---------------------------------------------------------------------- 
    290       INTEGER, INTENT( in )           ::   kt      ! time step 
    291       INTEGER, INTENT( in ), OPTIONAL ::   kcall   ! optional argument indicating call sequence 
     293      INTEGER, INTENT( in )           ::   kt             ! time step 
     294      INTEGER, INTENT( in )           ::   Kbb, Kmm, Kaa  ! time step 
     295      INTEGER, INTENT( in ), OPTIONAL ::   kcall          ! optional argument indicating call sequence 
    292296      ! 
    293297      INTEGER                ::   ji, jj, jk            ! dummy loop indices 
     
    321325      !                                                ! --------------------------------------------- ! 
    322326      ! 
    323       z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * ssmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) ) 
     327      z_scale(:,:) = ( ssh(:,:,Kaa) - ssh(:,:,Kbb) ) * ssmask(:,:) / ( ht_0(:,:) + ssh(:,:,Kmm) + 1. - ssmask(:,:) ) 
    324328      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) 
     329         ! formally this is the same as e3t(:,:,:,Kaa) = e3t_0*(1+ssha/ht_0) 
     330         e3t(:,:,jk,Kaa) = e3t(:,:,jk,Kbb) + e3t(:,:,jk,Kmm) * z_scale(:,:) * tmask(:,:,jk) 
    327331      END DO 
    328332      ! 
     
    337341         zht(:,:)   = 0._wp 
    338342         DO jk = 1, jpkm1 
    339             zhdiv(:,:) = zhdiv(:,:) + e3t_n(:,:,jk) * hdivn(:,:,jk) 
    340             zht  (:,:) = zht  (:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 
     343            zhdiv(:,:) = zhdiv(:,:) + e3t(:,:,jk,Kmm) * hdiv(:,:,jk) 
     344            zht  (:,:) = zht  (:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk) 
    341345         END DO 
    342346         zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) ) 
     
    348352               DO jk = 1, jpkm1 
    349353                  hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:)   & 
    350                      &          * ( hdiv_lf(:,:,jk) - e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) ) 
     354                     &          * ( hdiv_lf(:,:,jk) - e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) ) 
    351355               END DO 
    352356            ENDIF 
     
    361365         IF( ln_vvl_ztilde ) THEN     ! z_tilde case 
    362366            DO jk = 1, jpkm1 
    363                tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) ) 
     367               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) ) 
    364368            END DO 
    365369         ELSE                         ! layer case 
    366370            DO jk = 1, jpkm1 
    367                tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) -   e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk) 
     371               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) -   e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk) 
    368372            END DO 
    369373         ENDIF 
     
    476480            zht(:,:)  = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk) 
    477481         END DO 
    478          z_scale(:,:) =  - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) ) 
     482         z_scale(:,:) =  - zht(:,:) / ( ht_0(:,:) + ssh(:,:,Kmm) + 1. - ssmask(:,:) ) 
    479483         DO jk = 1, jpkm1 
    480             dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk) 
     484            dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t(:,:,jk,Kmm) * z_scale(:,:) * tmask(:,:,jk) 
    481485         END DO 
    482486 
     
    486490      !                                           ! ---baroclinic part--------- ! 
    487491         DO jk = 1, jpkm1 
    488             e3t_a(:,:,jk) = e3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk) 
     492            e3t(:,:,jk,Kaa) = e3t(:,:,jk,Kaa) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk) 
    489493         END DO 
    490494      ENDIF 
     
    501505         zht(:,:) = 0.0_wp 
    502506         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(:,:) ) ) 
     507            zht(:,:) = zht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk) 
     508         END DO 
     509         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kmm) - zht(:,:) ) ) 
    506510         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 
     511         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(e3t(:,:,:,Kmm)))) =', z_tmax 
    508512         ! 
    509513         zht(:,:) = 0.0_wp 
    510514         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(:,:) ) ) 
     515            zht(:,:) = zht(:,:) + e3t(:,:,jk,Kaa) * tmask(:,:,jk) 
     516         END DO 
     517         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kaa) - zht(:,:) ) ) 
    514518         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 
     519         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(e3t(:,:,:,Kaa)))) =', z_tmax 
    516520         ! 
    517521         zht(:,:) = 0.0_wp 
    518522         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(:,:) ) ) 
     523            zht(:,:) = zht(:,:) + e3t(:,:,jk,Kbb) * tmask(:,:,jk) 
     524         END DO 
     525         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kbb) - zht(:,:) ) ) 
    522526         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(:,:) ) ) 
     527         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(e3t(:,:,:,Kbb)))) =', z_tmax 
     528         ! 
     529         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssh(:,:,Kbb) ) ) 
    526530         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(:,:) ) ) 
     531         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kbb)))) =', z_tmax 
     532         ! 
     533         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssh(:,:,Kmm) ) ) 
    530534         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(:,:) ) ) 
     535         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kmm)))) =', z_tmax 
     536         ! 
     537         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssh(:,:,Kaa) ) ) 
    534538         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain 
    535          IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssha))) =', z_tmax 
     539         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kaa)))) =', z_tmax 
    536540      END IF 
    537541 
     
    540544      ! *********************************** ! 
    541545 
    542       CALL dom_vvl_interpol( e3t_a(:,:,:), e3u_a(:,:,:), 'U' ) 
    543       CALL dom_vvl_interpol( e3t_a(:,:,:), e3v_a(:,:,:), 'V' ) 
     546      CALL dom_vvl_interpol( e3t(:,:,:,Kaa), e3u(:,:,:,Kaa), 'U' ) 
     547      CALL dom_vvl_interpol( e3t(:,:,:,Kaa), e3v(:,:,:,Kaa), 'V' ) 
    544548 
    545549      ! *********************************** ! 
     
    547551      ! *********************************** ! 
    548552 
    549       hu_a(:,:) = e3u_a(:,:,1) * umask(:,:,1) 
    550       hv_a(:,:) = e3v_a(:,:,1) * vmask(:,:,1) 
     553      hu(:,:,Kaa) = e3u(:,:,1,Kaa) * umask(:,:,1) 
     554      hv(:,:,Kaa) = e3v(:,:,1,Kaa) * vmask(:,:,1) 
    551555      DO jk = 2, jpkm1 
    552          hu_a(:,:) = hu_a(:,:) + e3u_a(:,:,jk) * umask(:,:,jk) 
    553          hv_a(:,:) = hv_a(:,:) + e3v_a(:,:,jk) * vmask(:,:,jk) 
     556         hu(:,:,Kaa) = hu(:,:,Kaa) + e3u(:,:,jk,Kaa) * umask(:,:,jk) 
     557         hv(:,:,Kaa) = hv(:,:,Kaa) + e3v(:,:,jk,Kaa) * vmask(:,:,jk) 
    554558      END DO 
    555559      !                                        ! Inverse of the local depth 
    556560!!gm BUG ?  don't understand the use of umask_i here ..... 
    557       r1_hu_a(:,:) = ssumask(:,:) / ( hu_a(:,:) + 1._wp - ssumask(:,:) ) 
    558       r1_hv_a(:,:) = ssvmask(:,:) / ( hv_a(:,:) + 1._wp - ssvmask(:,:) ) 
     561      r1_hu(:,:,Kaa) = ssumask(:,:) / ( hu(:,:,Kaa) + 1._wp - ssumask(:,:) ) 
     562      r1_hv(:,:,Kaa) = ssvmask(:,:) / ( hv(:,:,Kaa) + 1._wp - ssvmask(:,:) ) 
    559563      ! 
    560564      IF( ln_timing )   CALL timing_stop('dom_vvl_sf_nxt') 
     
    563567 
    564568 
    565    SUBROUTINE dom_vvl_sf_swp( kt ) 
    566       !!---------------------------------------------------------------------- 
    567       !!                ***  ROUTINE dom_vvl_sf_swp  *** 
     569   SUBROUTINE dom_vvl_sf_update( kt, Kbb, Kmm, Kaa ) 
     570      !!---------------------------------------------------------------------- 
     571      !!                ***  ROUTINE dom_vvl_sf_update  *** 
    568572      !!                    
    569       !! ** Purpose :  compute time filter and swap of scale factors  
     573      !! ** Purpose :  for z tilde case: compute time filter and swap of scale factors  
    570574      !!               compute all depths and related variables for next time step 
    571575      !!               write outputs and restart file 
    572576      !! 
    573       !! ** Method  :  - swap of e3t with trick for volume/tracer conservation 
     577      !! ** Method  :  - swap of e3t with trick for volume/tracer conservation (ONLY FOR Z TILDE CASE) 
    574578      !!               - reconstruct scale factor at other grid points (interpolate) 
    575579      !!               - recompute depths and water height fields 
    576580      !! 
    577       !! ** Action  :  - e3t_(b/n), tilde_e3t_(b/n) and e3(u/v)_n ready for next time step 
     581      !! ** Action  :  - tilde_e3t_(b/n) ready for next time step 
    578582      !!               - Recompute: 
    579583      !!                    e3(u/v)_b        
    580       !!                    e3w_n            
     584      !!                    e3w(:,:,:,Kmm)            
    581585      !!                    e3(u/v)w_b       
    582586      !!                    e3(u/v)w_n       
    583       !!                    gdept_n, gdepw_n  and gde3w_n 
     587      !!                    gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm)  and gde3w 
    584588      !!                    h(u/v) and h(u/v)r 
    585589      !! 
     
    587591      !!              Leclair, M., and G. Madec, 2011, Ocean Modelling. 
    588592      !!---------------------------------------------------------------------- 
    589       INTEGER, INTENT( in ) ::   kt   ! time step 
     593      INTEGER, INTENT( in ) ::   kt              ! time step 
     594      INTEGER, INTENT( in ) ::   Kbb, Kmm, Kaa   ! time level indices 
    590595      ! 
    591596      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     
    595600      IF( ln_linssh )   RETURN      ! No calculation in linear free surface 
    596601      ! 
    597       IF( ln_timing )   CALL timing_start('dom_vvl_sf_swp') 
     602      IF( ln_timing )   CALL timing_start('dom_vvl_sf_update') 
    598603      ! 
    599604      IF( kt == nit000 )   THEN 
    600605         IF(lwp) WRITE(numout,*) 
    601          IF(lwp) WRITE(numout,*) 'dom_vvl_sf_swp : - time filter and swap of scale factors' 
    602          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~   - interpolate scale factors and compute depths for next time step' 
     606         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_update : - interpolate scale factors and compute depths for next time step' 
     607         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~' 
    603608      ENDIF 
    604609      ! 
     
    615620         tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:) 
    616621      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(:,:,:) 
     622      gdept(:,:,:,Kbb) = gdept(:,:,:,Kmm) 
     623      gdepw(:,:,:,Kbb) = gdepw(:,:,:,Kmm) 
     624 
     625      e3t(:,:,:,Kmm) = e3t(:,:,:,Kaa) 
     626      e3u(:,:,:,Kmm) = e3u(:,:,:,Kaa) 
     627      e3v(:,:,:,Kmm) = e3v(:,:,:,Kaa) 
    623628 
    624629      ! Compute all missing vertical scale factor and depths 
     
    626631      ! Horizontal scale factor interpolations 
    627632      ! -------------------------------------- 
    628       ! - ML - e3u_b and e3v_b are allready computed in dynnxt 
    629       ! - JC - hu_b, hv_b, hur_b, hvr_b also 
     633      ! - ML - e3u(:,:,:,Kbb) and e3v(:,:,:,Kbb) are already computed in dynnxt 
     634      ! - JC - hu(:,:,:,Kbb), hv(:,:,:,:,Kbb), hur_b, hvr_b also 
    630635       
    631       CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F'  ) 
     636      CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F'  ) 
    632637       
    633638      ! 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' ) 
     639      CALL dom_vvl_interpol( e3t(:,:,:,Kmm),  e3w(:,:,:,Kmm), 'W'  ) 
     640      CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' ) 
     641      CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' ) 
     642      CALL dom_vvl_interpol( e3t(:,:,:,Kbb),  e3w(:,:,:,Kbb), 'W'  ) 
     643      CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' ) 
     644      CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' ) 
    640645 
    641646      ! 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(:,:) 
     647      gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm) 
     648      gdepw(:,:,1,Kmm) = 0.0_wp 
     649      gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm) 
    645650      DO jk = 2, jpk 
    646651         DO jj = 1,jpj 
     
    649654                                                                 ! 1 for jk = mikt 
    650655               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) 
     656               gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 
     657               gdept(ji,jj,jk,Kmm) =    zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm) )  & 
     658                   &             + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm) )  
     659               gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm) 
    655660            END DO 
    656661         END DO 
     
    659664      ! Local depth and Inverse of the local depth of the water 
    660665      ! ------------------------------------------------------- 
    661       hu_n(:,:) = hu_a(:,:)   ;   r1_hu_n(:,:) = r1_hu_a(:,:) 
    662       hv_n(:,:) = hv_a(:,:)   ;   r1_hv_n(:,:) = r1_hv_a(:,:) 
    663       ! 
    664       ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1) 
     666      ! 
     667      ht(:,:) = e3t(:,:,1,Kmm) * tmask(:,:,1) 
    665668      DO jk = 2, jpkm1 
    666          ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 
     669         ht(:,:) = ht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk) 
    667670      END DO 
    668671 
    669672      ! write restart file 
    670673      ! ================== 
    671       IF( lrst_oce  )   CALL dom_vvl_rst( kt, 'WRITE' ) 
    672       ! 
    673       IF( ln_timing )   CALL timing_stop('dom_vvl_sf_swp') 
    674       ! 
    675    END SUBROUTINE dom_vvl_sf_swp 
     674      IF( lrst_oce  )   CALL dom_vvl_rst( kt, Kbb, Kmm, 'WRITE' ) 
     675      ! 
     676      IF( ln_timing )   CALL timing_stop('dom_vvl_sf_update') 
     677      ! 
     678   END SUBROUTINE dom_vvl_sf_update 
    676679 
    677680 
     
    783786 
    784787 
    785    SUBROUTINE dom_vvl_rst( kt, cdrw ) 
     788   SUBROUTINE dom_vvl_rst( kt, Kbb, Kmm, cdrw ) 
    786789      !!--------------------------------------------------------------------- 
    787790      !!                   ***  ROUTINE dom_vvl_rst  *** 
     
    795798      !!                they are set to 0. 
    796799      !!---------------------------------------------------------------------- 
    797       INTEGER         , INTENT(in) ::   kt     ! ocean time-step 
    798       CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag 
     800      INTEGER         , INTENT(in) ::   kt        ! ocean time-step 
     801      INTEGER         , INTENT(in) ::   Kbb, Kmm  ! ocean time level indices 
     802      CHARACTER(len=*), INTENT(in) ::   cdrw      ! "READ"/"WRITE" flag 
    799803      ! 
    800804      INTEGER ::   ji, jj, jk 
     
    806810         IF( ln_rstart ) THEN                   !* Read the restart file 
    807811            CALL rst_read_open                  !  open the restart file if necessary 
    808             CALL iom_get( numror, jpdom_autoglo, 'sshn'   , sshn, ldxios = lrxios    ) 
     812            CALL iom_get( numror, jpdom_autoglo, 'sshn'   , ssh(:,:,Kmm), ldxios = lrxios    ) 
    809813            ! 
    810814            id1 = iom_varid( numror, 'e3t_b', ldstop = .FALSE. ) 
     
    817821            !                             ! --------- ! 
    818822            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 ) 
     823               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios ) 
     824               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios ) 
    821825               ! 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' 
     826               IF(lwp) write(numout,*) 'dom_vvl_rst : e3t(:,:,:,Kbb) and e3t(:,:,:,Kmm) found in restart files' 
    823827               WHERE ( tmask(:,:,:) == 0.0_wp )  
    824                   e3t_n(:,:,:) = e3t_0(:,:,:) 
    825                   e3t_b(:,:,:) = e3t_0(:,:,:) 
     828                  e3t(:,:,:,Kmm) = e3t_0(:,:,:) 
     829                  e3t(:,:,:,Kbb) = e3t_0(:,:,:) 
    826830               END WHERE 
    827831               IF( neuler == 0 ) THEN 
    828                   e3t_b(:,:,:) = e3t_n(:,:,:) 
     832                  e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 
    829833               ENDIF 
    830834            ELSE IF( id1 > 0 ) THEN 
    831                IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart files' 
     835               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart files' 
    832836               IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.' 
    833837               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(:,:,:) 
     838               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios ) 
     839               e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb) 
    836840               neuler = 0 
    837841            ELSE IF( id2 > 0 ) THEN 
    838                IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_b not found in restart files' 
     842               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kbb) not found in restart files' 
    839843               IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.' 
    840844               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(:,:,:) 
     845               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios ) 
     846               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 
    843847               neuler = 0 
    844848            ELSE 
    845                IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart file' 
     849               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart file' 
    846850               IF(lwp) write(numout,*) 'Compute scale factor from sshn' 
    847851               IF(lwp) write(numout,*) 'neuler is forced to 0' 
    848852               DO jk = 1, jpk 
    849                   e3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) & 
     853                  e3t(:,:,jk,Kmm) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm) ) & 
    850854                      &                          / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   & 
    851855                      &          + e3t_0(:,:,jk)                               * (1._wp -tmask(:,:,jk)) 
    852856               END DO 
    853                e3t_b(:,:,:) = e3t_n(:,:,:) 
     857               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 
    854858               neuler = 0 
    855859            ENDIF 
     
    888892               IF( cn_cfg == 'wad' ) THEN 
    889893                  ! 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  (:,:,:) 
     894                  CALL usr_def_istate( gdept(:,:,:,Kbb), tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  ) 
     895                  ts  (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)       ! set now values from to before ones 
     896                  ssh (:,:,Kmm)     = ssh(:,:,Kbb) 
     897                  uu   (:,:,:,Kmm)   = uu  (:,:,:,Kbb) 
     898                  vv   (:,:,:,Kmm)   = vv  (:,:,:,Kbb) 
    895899               ELSE 
    896900                  ! if not test case 
    897                   sshn(:,:) = -ssh_ref 
    898                   sshb(:,:) = -ssh_ref 
     901                  ssh(:,:,Kmm) = -ssh_ref 
     902                  ssh(:,:,Kbb) = -ssh_ref 
    899903 
    900904                  DO jj = 1, jpj 
    901905                     DO ji = 1, jpi 
    902906                        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) ) 
     907                           ssh(ji,jj,Kbb) = rn_wdmin1 - (ht_0(ji,jj) ) 
     908                           ssh(ji,jj,Kmm) = rn_wdmin1 - (ht_0(ji,jj) ) 
    907909                        ENDIF 
    908910                     ENDDO 
     
    912914               ! Adjust vertical metrics for all wad 
    913915               DO jk = 1, jpk 
    914                   e3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:)  ) & 
     916                  e3t(:,:,jk,Kmm) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm)  ) & 
    915917                    &                            / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   & 
    916918                    &            + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) ) 
    917919               END DO 
    918                e3t_b(:,:,:) = e3t_n(:,:,:) 
     920               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 
    919921 
    920922               DO ji = 1, jpi 
     
    928930            ELSE 
    929931               ! 
    930                ! usr_def_istate called here only to get sshb, that is needed to initialize e3t_b and e3t_n 
    931                CALL usr_def_istate( gdept_0, tmask, tsb, ub, vb, sshb  )   
     932               ! usr_def_istate called here only to get sshb, that is needed to initialize e3t(Kbb) and e3t(Kmm) 
     933               CALL usr_def_istate( gdept_0, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  ) 
    932934               ! usr_def_istate will be called again in istate_init to initialize ts(bn), ssh(bn), u(bn) and v(bn) 
    933935               ! 
    934936               DO jk=1,jpk 
    935                   e3t_b(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshb(:,:) ) & 
    936                     &                            / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   & 
    937                     &            + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) )   ! make sure e3t_b != 0 on land points 
     937                  e3t(:,:,jk,Kmm) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kbb) ) & 
     938                    &                              / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   & 
     939                    &              + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) )   ! make sure e3t_b != 0 on land points 
    938940               END DO 
    939                e3t_n(:,:,:) = e3t_b(:,:,:) 
    940                sshn(:,:) = sshb(:,:)   ! needed later for gde3w 
    941 !!$                e3t_n(:,:,:)=e3t_0(:,:,:) 
    942 !!$                e3t_b(:,:,:)=e3t_0(:,:,:) 
     941               e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb) 
     942               ssh(:,:  ,Kmm) = ssh(:,:  ,Kbb)   ! needed later for gde3w 
    943943               ! 
    944944            END IF           ! end of ll_wd edits 
     
    958958         !                                           ! all cases ! 
    959959         !                                           ! --------- ! 
    960          CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t_b(:,:,:), ldxios = lwxios ) 
    961          CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t_n(:,:,:), ldxios = lwxios ) 
     960         CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lwxios ) 
     961         CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lwxios ) 
    962962         !                                           ! ----------------------- ! 
    963963         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/tests/CANAL/MY_SRC/stpctl.F90

    r10572 r11771  
    4242CONTAINS 
    4343 
    44    SUBROUTINE stp_ctl( kt, kindic ) 
     44   SUBROUTINE stp_ctl( kt, Kbb, Kmm, kindic ) 
    4545      !!---------------------------------------------------------------------- 
    4646      !!                    ***  ROUTINE stp_ctl  *** 
     
    6060      !!---------------------------------------------------------------------- 
    6161      INTEGER, INTENT(in   ) ::   kt       ! ocean time-step index 
     62      INTEGER, INTENT(in   ) ::   Kbb, Kmm      ! ocean time level index 
    6263      INTEGER, INTENT(inout) ::   kindic   ! error indicator 
    6364      !! 
     
    111112      !                                   !==  test of extrema  ==! 
    112113      IF( ll_wd ) THEN 
    113          zmax(1) = MAXVAL(  ABS( sshn(:,:) + ssh_ref*tmask(:,:,1) )  )        ! ssh max  
     114         zmax(1) = MAXVAL(  ABS( ssh(:,:,Kmm) + ssh_ref*tmask(:,:,1) )  )        ! ssh max  
    114115      ELSE 
    115          zmax(1) = MAXVAL(  ABS( sshn(:,:) )  )                               ! ssh max 
     116         zmax(1) = MAXVAL(  ABS( ssh(:,:,Kmm) )  )                               ! ssh max 
    116117      ENDIF 
    117       zmax(2) = MAXVAL(  ABS( un(:,:,:) )  )                                  ! velocity max (zonal only) 
    118       zmax(3) = MAXVAL( -tsn(:,:,:,jp_sal) , mask = tmask(:,:,:) == 1._wp )   ! minus salinity max 
    119       zmax(4) = MAXVAL(  tsn(:,:,:,jp_sal) , mask = tmask(:,:,:) == 1._wp )   !       salinity max 
    120       zmax(5) = MAXVAL( -tsn(:,:,:,jp_tem) , mask = tmask(:,:,:) == 1._wp )   ! minus temperature max 
    121       zmax(6) = MAXVAL(  tsn(:,:,:,jp_tem) , mask = tmask(:,:,:) == 1._wp )   !       temperature max 
     118      zmax(2) = MAXVAL(  ABS( uu(:,:,:,Kmm) )  )                                  ! velocity max (zonal only) 
     119      zmax(3) = MAXVAL( -ts(:,:,:,jp_sal,Kmm) , mask = tmask(:,:,:) == 1._wp )   ! minus salinity max 
     120      zmax(4) = MAXVAL(  ts(:,:,:,jp_sal,Kmm) , mask = tmask(:,:,:) == 1._wp )   !       salinity max 
     121      zmax(5) = MAXVAL( -ts(:,:,:,jp_tem,Kmm) , mask = tmask(:,:,:) == 1._wp )   ! minus temperature max 
     122      zmax(6) = MAXVAL(  ts(:,:,:,jp_tem,Kmm) , mask = tmask(:,:,:) == 1._wp )   !       temperature max 
    122123      zmax(7) = REAL( nstop , wp )                                            ! stop indicator 
    123124      IF( ln_zad_Aimp ) THEN 
     
    148149      !                                   !==  error handling  ==! 
    149150      IF( ( ln_ctl .OR. lsomeoce ) .AND. (   &             ! have use mpp_max (because ln_ctl=.T.) or contains some ocean points 
    150          &  zmax(1) >   50._wp .OR.   &                    ! too large sea surface height ( > 50 m ) 
    151          &  zmax(2) >   20._wp .OR.   &                    ! too large velocity ( > 20 m/s) 
    152 !!$         &  zmax(3) >=   0._wp .OR.   &                    ! negative or zero sea surface salinity 
    153 !!$         &  zmax(4) >= 100._wp .OR.   &                    ! too large sea surface salinity ( > 100 ) 
    154 !!$         &  zmax(4) <    0._wp .OR.   &                    ! too large sea surface salinity (keep this line for sea-ice) 
     151         &  zmax(1) >   20._wp .OR.   &                    ! too large sea surface height ( > 20 m ) 
     152         &  zmax(2) >   10._wp .OR.   &                    ! too large velocity ( > 10 m/s) 
     153!!$      &  zmax(3) >=   0._wp .OR.   &                    ! negative or zero sea surface salinity 
     154!!$      &  zmax(4) >= 100._wp .OR.   &                    ! too large sea surface salinity ( > 100 ) 
     155!!$      &  zmax(4) <    0._wp .OR.   &                    ! too large sea surface salinity (keep this line for sea-ice) 
    155156         &  ISNAN( zmax(1) + zmax(2) + zmax(3) ) ) ) THEN   ! NaN encounter in the tests 
    156157         IF( lk_mpp .AND. ln_ctl ) THEN 
    157             CALL mpp_maxloc( 'stpctl', ABS(sshn)        , ssmask(:,:)  , zzz, ih  ) 
    158             CALL mpp_maxloc( 'stpctl', ABS(un)          , umask (:,:,:), zzz, iu  ) 
    159             CALL mpp_minloc( 'stpctl', tsn(:,:,:,jp_sal), tmask (:,:,:), zzz, is1 ) 
    160             CALL mpp_maxloc( 'stpctl', tsn(:,:,:,jp_sal), tmask (:,:,:), zzz, is2 ) 
     158            CALL mpp_maxloc( 'stpctl', ABS(ssh(:,:,Kmm))        , ssmask(:,:)  , zzz, ih  ) 
     159            CALL mpp_maxloc( 'stpctl', ABS(uu(:,:,:,Kmm))          , umask (:,:,:), zzz, iu  ) 
     160            CALL mpp_minloc( 'stpctl', ts(:,:,:,jp_sal,Kmm), tmask (:,:,:), zzz, is1 ) 
     161            CALL mpp_maxloc( 'stpctl', ts(:,:,:,jp_sal,Kmm), tmask (:,:,:), zzz, is2 ) 
    161162         ELSE 
    162             ih(:)  = MAXLOC( ABS( sshn(:,:)   )                              ) + (/ nimpp - 1, njmpp - 1    /) 
    163             iu(:)  = MAXLOC( ABS( un  (:,:,:) )                              ) + (/ nimpp - 1, njmpp - 1, 0 /) 
    164             is1(:) = MINLOC( tsn(:,:,:,jp_sal), mask = tmask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /) 
    165             is2(:) = MAXLOC( tsn(:,:,:,jp_sal), mask = tmask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /) 
     163            ih(:)  = MAXLOC( ABS( ssh(:,:,Kmm)   )                              ) + (/ nimpp - 1, njmpp - 1    /) 
     164            iu(:)  = MAXLOC( ABS( uu  (:,:,:,Kmm) )                              ) + (/ nimpp - 1, njmpp - 1, 0 /) 
     165            is1(:) = MINLOC( ts(:,:,:,jp_sal,Kmm), mask = tmask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /) 
     166            is2(:) = MAXLOC( ts(:,:,:,jp_sal,Kmm), mask = tmask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /) 
    166167         ENDIF 
    167168          
    168          WRITE(ctmp1,*) ' stp_ctl: |ssh| > 50 m  or  |U| > 20 m/s  or  NaN encounter in the tests' 
     169         WRITE(ctmp1,*) ' stp_ctl: |ssh| > 20 m  or  |U| > 10 m/s  or  S <= 0  or  S >= 100  or  NaN encounter in the tests' 
    169170         WRITE(ctmp2,9100) kt,   zmax(1), ih(1) , ih(2) 
    170171         WRITE(ctmp3,9200) kt,   zmax(2), iu(1) , iu(2) , iu(3) 
     
    173174         WRITE(ctmp6,*) '      ===> output of last computed fields in output.abort.nc file' 
    174175          
    175          CALL dia_wri_state( 'output.abort' )     ! create an output.abort file 
     176         CALL dia_wri_state( Kmm, 'output.abort' )     ! create an output.abort file 
    176177          
    177178         IF( .NOT. ln_ctl ) THEN 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/tests/CANAL/MY_SRC/trazdf.F90

    r10572 r11771  
    4444CONTAINS 
    4545 
    46    SUBROUTINE tra_zdf( kt ) 
     46   SUBROUTINE tra_zdf( kt, Kbb, Kmm, Krhs, pts, Kaa ) 
    4747      !!---------------------------------------------------------------------- 
    4848      !!                  ***  ROUTINE tra_zdf  *** 
     
    5050      !! ** Purpose :   compute the vertical ocean tracer physics. 
    5151      !!--------------------------------------------------------------------- 
    52       INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     52      INTEGER                                  , INTENT(in)    :: kt                  ! ocean time-step index 
     53      INTEGER                                  , INTENT(in)    :: Kbb, Kmm, Krhs, Kaa ! time level indices 
     54      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts                 ! active tracers and RHS of tracer equation 
    5355      ! 
    5456      INTEGER  ::   jk   ! Dummy loop indices 
     
    7072      IF( l_trdtra )   THEN                  !* Save ta and sa trends 
    7173         ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) ) 
    72          ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 
    73          ztrds(:,:,:) = tsa(:,:,:,jp_sal) 
     74         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Kaa) 
     75         ztrds(:,:,:) = pts(:,:,:,jp_sal,Kaa) 
    7476      ENDIF 
    7577      ! 
    7678      !                                      !* compute lateral mixing trend and add it to the general trend 
    77       CALL tra_zdf_imp( kt, nit000, 'TRA', r2dt, tsb, tsa, jpts )  
     79      CALL tra_zdf_imp( kt, nit000, 'TRA', r2dt, Kbb, Kmm, Krhs, pts, Kaa, jpts )  
    7880 
    7981!!gm WHY here !   and I don't like that ! 
     
    8183      ! JMM avoid negative salinities near river outlet ! Ugly fix 
    8284      ! JMM : restore negative salinities to small salinities: 
    83 !!$      WHERE( tsa(:,:,:,jp_sal) < 0._wp )   tsa(:,:,:,jp_sal) = 0.1_wp 
     85!!$   WHERE( pts(:,:,:,jp_sal,Kaa) < 0._wp )   pts(:,:,:,jp_sal,Kaa) = 0.1_wp 
    8486!!gm 
    8587 
    8688      IF( l_trdtra )   THEN                      ! save the vertical diffusive trends for further diagnostics 
    8789         DO jk = 1, jpkm1 
    88             ztrdt(:,:,jk) = ( ( tsa(:,:,jk,jp_tem)*e3t_a(:,:,jk) - tsb(:,:,jk,jp_tem)*e3t_b(:,:,jk) ) & 
    89                &          / (e3t_n(:,:,jk)*r2dt) ) - ztrdt(:,:,jk) 
    90             ztrds(:,:,jk) = ( ( tsa(:,:,jk,jp_sal)*e3t_a(:,:,jk) - tsb(:,:,jk,jp_sal)*e3t_b(:,:,jk) ) & 
    91               &           / (e3t_n(:,:,jk)*r2dt) ) - ztrds(:,:,jk) 
     90            ztrdt(:,:,jk) = ( ( pts(:,:,jk,jp_tem,Kaa)*e3t(:,:,jk,Kaa) - pts(:,:,jk,jp_tem,Kbb)*e3t(:,:,jk,Kbb) ) & 
     91               &          / (e3t(:,:,jk,Kmm)*r2dt) ) - ztrdt(:,:,jk) 
     92            ztrds(:,:,jk) = ( ( pts(:,:,jk,jp_sal,Kaa)*e3t(:,:,jk,Kaa) - pts(:,:,jk,jp_sal,Kbb)*e3t(:,:,jk,Kbb) ) & 
     93              &           / (e3t(:,:,jk,Kmm)*r2dt) ) - ztrds(:,:,jk) 
    9294         END DO 
    9395!!gm this should be moved in trdtra.F90 and done on all trends 
    9496         CALL lbc_lnk_multi( 'trazdf', ztrdt, 'T', 1. , ztrds, 'T', 1. ) 
    9597!!gm 
    96          CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdf, ztrdt ) 
    97          CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdf, ztrds ) 
     98         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_zdf, ztrdt ) 
     99         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_zdf, ztrds ) 
    98100         DEALLOCATE( ztrdt , ztrds ) 
    99101      ENDIF 
    100102      !                                          ! print mean trends (used for debugging) 
    101       IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' zdf  - Ta: ', mask1=tmask,               & 
    102          &                       tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 
     103      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Kaa), clinfo1=' zdf  - Ta: ', mask1=tmask,               & 
     104         &                       tab3d_2=pts(:,:,:,jp_sal,Kaa), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 
    103105      ! 
    104106      IF( ln_timing )   CALL timing_stop('tra_zdf') 
     
    107109 
    108110  
    109    SUBROUTINE tra_zdf_imp( kt, kit000, cdtype, p2dt, ptb, pta, kjpt )  
     111   SUBROUTINE tra_zdf_imp( kt, kit000, cdtype, p2dt, Kbb, Kmm, Krhs, pt, Kaa, kjpt )  
    110112      !!---------------------------------------------------------------------- 
    111113      !!                  ***  ROUTINE tra_zdf_imp  *** 
     
    125127      !!      If iso-neutral mixing, add to avt the contribution due to lateral mixing. 
    126128      !! 
    127       !! ** Action  : - pta  becomes the after tracer 
    128       !!--------------------------------------------------------------------- 
    129       INTEGER                              , INTENT(in   ) ::   kt       ! ocean time-step index 
    130       INTEGER                              , INTENT(in   ) ::   kit000   ! first time step index 
    131       CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator) 
    132       INTEGER                              , INTENT(in   ) ::   kjpt     ! number of tracers 
    133       REAL(wp)                             , INTENT(in   ) ::   p2dt     ! tracer time-step 
    134       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb      ! before and now tracer fields 
    135       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta      ! in: tracer trend ; out: after tracer field 
     129      !! ** Action  : - pt(:,:,:,:,Kaa)  becomes the after tracer 
     130      !!--------------------------------------------------------------------- 
     131      INTEGER                                  , INTENT(in   ) ::   kt       ! ocean time-step index 
     132      INTEGER                                  , INTENT(in   ) ::   Kbb, Kmm, Krhs, Kaa  ! ocean time level indices 
     133      INTEGER                                  , INTENT(in   ) ::   kit000   ! first time step index 
     134      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator) 
     135      INTEGER                                  , INTENT(in   ) ::   kjpt     ! number of tracers 
     136      REAL(wp)                                 , INTENT(in   ) ::   p2dt     ! tracer time-step 
     137      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt       ! tracers and RHS of tracer equation 
    136138      ! 
    137139      INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices 
     
    181183                  DO jj = 2, jpjm1 
    182184                     DO ji = fs_2, fs_jpim1   ! vector opt. (ensure same order of calculation as below if wi=0.) 
    183                         zzwi = - p2dt * zwt(ji,jj,jk  ) / e3w_n(ji,jj,jk  ) 
    184                         zzws = - p2dt * zwt(ji,jj,jk+1) / e3w_n(ji,jj,jk+1) 
    185                         zwd(ji,jj,jk) = e3t_a(ji,jj,jk) - zzwi - zzws   & 
     185                        zzwi = - p2dt * zwt(ji,jj,jk  ) / e3w(ji,jj,jk  ,Kmm) 
     186                        zzws = - p2dt * zwt(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm) 
     187                        zwd(ji,jj,jk) = e3t(ji,jj,jk,Kaa) - zzwi - zzws   & 
    186188                           &                 + p2dt * ( MAX( wi(ji,jj,jk  ) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) ) 
    187189                        zwi(ji,jj,jk) = zzwi + p2dt *   MIN( wi(ji,jj,jk  ) , 0._wp ) 
     
    194196                  DO jj = 2, jpjm1 
    195197                     DO ji = fs_2, fs_jpim1   ! vector opt. 
    196                         zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk  ) / e3w_n(ji,jj,jk) 
    197                         zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / e3w_n(ji,jj,jk+1) 
    198                         zwd(ji,jj,jk) = e3t_a(ji,jj,jk) - zwi(ji,jj,jk) - zws(ji,jj,jk) 
     198                        zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk  ) / e3w(ji,jj,jk,Kmm) 
     199                        zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm) 
     200                        zwd(ji,jj,jk) = e3t(ji,jj,jk,Kaa) - zwi(ji,jj,jk) - zws(ji,jj,jk) 
    199201                    END DO 
    200202                  END DO 
     
    218220            !   The solution will be in the 4d array pta. 
    219221            !   The 3d array zwt is used as a work space array. 
    220             !   En route to the solution pta is used a to evaluate the rhs and then  
     222            !   En route to the solution pt(:,:,:,:,Kaa) is used a to evaluate the rhs and then  
    221223            !   used as a work space array: its value is modified. 
    222224            ! 
     
    238240         DO jj = 2, jpjm1           !* 2nd recurrence:    Zk = Yk - Ik / Tk-1  Zk-1 
    239241            DO ji = fs_2, fs_jpim1 
    240                pta(ji,jj,1,jn) = e3t_b(ji,jj,1) * ptb(ji,jj,1,jn) + p2dt * e3t_n(ji,jj,1) * pta(ji,jj,1,jn) 
     242               pt(ji,jj,1,jn,Kaa) = e3t(ji,jj,1,Kbb) * pt(ji,jj,1,jn,Kbb) + p2dt * e3t(ji,jj,1,Kmm) * pt(ji,jj,1,jn,Krhs) 
    241243            END DO 
    242244         END DO 
     
    244246            DO jj = 2, jpjm1 
    245247               DO ji = fs_2, fs_jpim1 
    246                   zrhs = e3t_b(ji,jj,jk) * ptb(ji,jj,jk,jn) + p2dt * e3t_n(ji,jj,jk) * pta(ji,jj,jk,jn)   ! zrhs=right hand side 
    247                   pta(ji,jj,jk,jn) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pta(ji,jj,jk-1,jn) 
     248                  zrhs = e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) + p2dt * e3t(ji,jj,jk,Kmm) * pt(ji,jj,jk,jn,Krhs)   ! zrhs=right hand side 
     249                  pt(ji,jj,jk,jn,Kaa) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pt(ji,jj,jk-1,jn,Kaa) 
    248250               END DO 
    249251            END DO 
     
    252254         DO jj = 2, jpjm1           !* 3d recurrence:    Xk = (Zk - Sk Xk+1 ) / Tk   (result is the after tracer) 
    253255            DO ji = fs_2, fs_jpim1 
    254                pta(ji,jj,jpkm1,jn) = pta(ji,jj,jpkm1,jn) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1) 
     256               pt(ji,jj,jpkm1,jn,Kaa) = pt(ji,jj,jpkm1,jn,Kaa) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1) 
    255257            END DO 
    256258         END DO 
     
    258260            DO jj = 2, jpjm1 
    259261               DO ji = fs_2, fs_jpim1 
    260                   pta(ji,jj,jk,jn) = ( pta(ji,jj,jk,jn) - zws(ji,jj,jk) * pta(ji,jj,jk+1,jn) )   & 
     262                  pt(ji,jj,jk,jn,Kaa) = ( pt(ji,jj,jk,jn,Kaa) - zws(ji,jj,jk) * pt(ji,jj,jk+1,jn,Kaa) )   & 
    261263                     &             / zwt(ji,jj,jk) * tmask(ji,jj,jk) 
    262264               END DO 
  • NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/tests/CANAL/MY_SRC/usrdef_sbc.F90

    r10976 r11771  
    4040CONTAINS 
    4141 
    42    SUBROUTINE usrdef_sbc_oce( kt, Kbb ) 
     42   SUBROUTINE usrdef_sbc_oce( kt, Kmm, Kbb ) 
    4343      !!--------------------------------------------------------------------- 
    4444      !!                    ***  ROUTINE usr_def_sbc  *** 
     
    5454      !! 
    5555      !!---------------------------------------------------------------------- 
    56       INTEGER, INTENT(in) ::   kt   ! ocean time step 
    57       INTEGER, INTENT(in) ::   Kbb  ! ocean time index 
     56      INTEGER, INTENT(in) ::   kt        ! ocean time step 
     57      INTEGER, INTENT(in) ::   Kbb, Kmm  ! ocean time index 
    5858      INTEGER  ::   ji, jj               ! dummy loop indices 
    5959      REAL(wp) :: zrhoair = 1.22     ! approximate air density [Kg/m3] 
     
    8888          
    8989         WHERE( ABS(gphit) <= rn_windszy/2. ) 
    90             zwndrel(:,:) = rn_u10 - rn_uofac * un(:,:,1) 
     90            zwndrel(:,:) = rn_u10 - rn_uofac * uu(:,:,1,Kmm) 
    9191         ELSEWHERE 
    92             zwndrel(:,:) =        - rn_uofac * un(:,:,1) 
     92            zwndrel(:,:) =        - rn_uofac * uu(:,:,1,Kmm) 
    9393         END WHERE 
    9494         utau(:,:) = zrhocd * zwndrel(:,:) * zwndrel(:,:) 
    9595 
    96          zwndrel(:,:) = - rn_uofac * vn(:,:,1) 
     96         zwndrel(:,:) = - rn_uofac * vv(:,:,1,Kmm) 
    9797         vtau(:,:) = zrhocd * zwndrel(:,:) * zwndrel(:,:) 
    9898 
Note: See TracChangeset for help on using the changeset viewer.