Changeset 11412


Ignore:
Timestamp:
2019-08-06T17:57:38+02:00 (14 months ago)
Author:
acc
Message:

Branch dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps: commit changes to VORTEX/MY_SRC/domvvl.F90 that were omitted from previous commits. This still fails restartability (in common with all AGRIF tests on this branch) but will compile and run.

File:
1 edited

Legend:

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

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