New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 10001 – NEMO

Changeset 10001


Ignore:
Timestamp:
2018-07-26T09:50:51+02:00 (5 years ago)
Author:
gm
Message:

#1911 (ENHANCE-04): RK3 branch - step I.1 and I.2 (see wiki page)

Location:
NEMO/branches/2018/dev_r9838_ENHANCE04_RK3
Files:
27 edited
1 copied

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/cfgs/AGRIF_DEMO/EXPREF/1_namelist_cfg

    r9939 r10001  
    2323   nn_it000    =       1   !  first time step 
    2424   nn_itend    =      75   !  last  time step (std 5475) 
    25    nn_istate   =       0   !  output the initial state (1) or not (0) 
     25   ln_istate   = .false.   !  output the initial state 
    2626/ 
    2727!----------------------------------------------------------------------- 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/cfgs/AGRIF_DEMO/EXPREF/2_namelist_cfg

    r9939 r10001  
    2323   nn_it000    =    1      !  first time step 
    2424   nn_itend    =  300      !  last  time step 
    25    nn_istate   =      0    !  output the initial state (1) or not (0) 
     25   ln_istate   = .false.   !  output the initial state 
    2626   ln_clobber  = .true.    !  clobber (overwrite) an existing file 
    2727/ 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/cfgs/AGRIF_DEMO/EXPREF/3_namelist_cfg

    r9939 r10001  
    2323   nn_it000    =    1      !  first time step 
    2424   nn_itend    =  900      !  last  time step 
    25    nn_istate   =      0    !  output the initial state (1) or not (0) 
     25   ln_istate   = .false.   !  output the initial state 
    2626   ln_clobber  = .true.    !  clobber (overwrite) an existing file 
    2727/ 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/cfgs/AGRIF_DEMO/EXPREF/namelist_cfg

    r9939 r10001  
    2323   nn_it000    =       1   !  first time step 
    2424   nn_itend    =      75   !  last  time step (std 5475) 
    25    nn_istate   =       0   !  output the initial state (1) or not (0) 
     25   ln_istate   = .false.   !  output the initial state 
    2626/ 
    2727!----------------------------------------------------------------------- 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/cfgs/GYRE_BFM/EXPREF/namelist_cfg

    r9939 r10001  
    2626   nn_stock    =    4320   !  frequency of creation of a restart file (modulo referenced to 1) 
    2727   nn_write    =      60   !  frequency of write in the output file   (modulo referenced to nn_it000) 
    28    nn_istate   =       0   !  output the initial state (1) or not (0) 
     28   ln_istate   = .false.   !  output the initial state 
    2929/ 
    3030!----------------------------------------------------------------------- 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/cfgs/GYRE_PISCES/EXPREF/namelist_cfg

    r9939 r10001  
    157157&namdyn_adv    !   formulation of the momentum advection                (default: NO selection) 
    158158!----------------------------------------------------------------------- 
    159    ln_dynadv_vec = .true.  !  vector form - 2nd centered scheme 
    160      nn_dynkeg     = 0        ! grad(KE) scheme: =0   C2  ;  =1   Hollingsworth correction 
     159   ln_dynadv_cen2= .true.  !  flux form - 2nd order centered scheme 
     160   ln_dynadv_ubs = .false. !  flux form - 3rd order UBS      scheme 
    161161/ 
    162162!----------------------------------------------------------------------- 
    163163&namdyn_vor    !   Vorticity / Coriolis scheme                          (default: NO selection) 
    164164!----------------------------------------------------------------------- 
    165    ln_dynvor_ene = .true.  !  energy conserving scheme 
     165   ln_dynvor_enT = .true.  !  energy conserving scheme (T-point) 
    166166/ 
    167167!----------------------------------------------------------------------- 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_cfg

    r9939 r10001  
    2323   nn_it000    =       1   !  first time step 
    2424   nn_itend    =    5475   !  last  time step (std 5475) 
    25    nn_istate   =       0   !  output the initial state (1) or not (0) 
     25   ln_istate   = .false.   !  output the initial state 
    2626/ 
    2727!----------------------------------------------------------------------- 
    … …  
    295295&namdyn_adv    !   formulation of the momentum advection                (default: NO selection) 
    296296!----------------------------------------------------------------------- 
    297    ln_dynadv_vec = .true.  !  vector form - 2nd centered scheme 
    298      nn_dynkeg     = 0        ! grad(KE) scheme: =0   C2  ;  =1   Hollingsworth correction 
     297   ln_dynadv_cen2= .true.  !  flux form - 2nd order centered scheme 
     298   ln_dynadv_ubs = .false. !  flux form - 3rd order UBS      scheme 
    299299/ 
    300300!----------------------------------------------------------------------- 
    301301&namdyn_vor    !   Vorticity / Coriolis scheme                          (default: NO selection) 
    302302!----------------------------------------------------------------------- 
    303    ln_dynvor_een = .true.  !  energy & enstrophy scheme 
    304       nn_een_e3f = 0          ! =0   e3f = mean masked e3t divided by 4 
     303   ln_dynvor_enT = .true.  !  energy conserving scheme (T-point) 
    305304/ 
    306305!----------------------------------------------------------------------- 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/cfgs/SHARED/namelist_ref

    r9939 r10001  
    5252      cn_ocerst_outdir= "."         !  directory in which to write output ocean restarts 
    5353   ln_iscpl    = .false.   !  cavity evolution forcing or coupling to ice sheet model 
    54    nn_istate   =       0   !  output the initial state (1) or not (0) 
     54   ln_istate   = .false.   !  output the initial state 
    5555   ln_rst_list = .false.   !  output restarts at list of times using nn_stocklist (T) or at set frequency with nn_stock (F) 
    5656   nn_stock    =    5475   !  frequency of creation of a restart file (modulo referenced to 1) 
    … …  
    6767&namdom        !   time and space domain 
    6868!----------------------------------------------------------------------- 
     69   ln_RK3      = .false.   !  3rd order Runge-Kutta time stepping scheme 
     70   ln_MLF      = .false.   !  Modified Leap-Frog    time stepping scheme 
     71      rn_atfp     =    0.1    !  asselin time filter parameter 
     72   ! 
     73   rn_Dt       = 5760.     !  ocean time step (dynamics and T-S) 
     74   ! 
    6975   ln_linssh   = .false.   !  =T  linear free surface  ==>>  model level are fixed in time 
     76   ! 
    7077   rn_isfhmin  =    1.00   !  treshold [m] to discriminate grounding ice from floating ice 
    71    ! 
    72    rn_Dt       = 5760.     !  time step for the dynamics and tracer 
    73    rn_atfp     =    0.1    !  asselin time filter parameter 
    7478   ! 
    7579   ln_crs      = .false.   !  Logical switch for coarsening module      (T => fill namcrs) 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/ICE/iceistate.F90

    r9939 r10001  
    9494      REAL(wp) ::   zarg, zV, zconv, zdv, zfac 
    9595      INTEGER , DIMENSION(4)           ::   itest 
    96       REAL(wp), DIMENSION(jpi,jpj)     ::   z2d 
    9796      REAL(wp), DIMENSION(jpi,jpj)     ::   zswitch    ! ice indicator 
    9897      REAL(wp), DIMENSION(jpi,jpj)     ::   zht_i_ini, zat_i_ini, zvt_i_ini            !data from namelist or nc file 
    9998      REAL(wp), DIMENSION(jpi,jpj)     ::   zts_u_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file 
    10099      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zh_i_ini , za_i_ini                        !data by cattegories to fill 
     100      REAL(wp), DIMENSION(jpi,jpj)     ::   z_ssh_h0, zsshu, zsshv, zsshf 
    101101      !-------------------------------------------------------------------- 
    102102 
    … …  
    413413      snwice_mass_b(:,:) = snwice_mass(:,:) 
    414414      ! 
    415       IF( ln_ice_embd ) THEN            ! embedded sea-ice: deplete the initial ssh below sea-ice area 
     415      IF( ln_ice_embd ) THEN           ! embedded sea-ice: deplete the initial ssh below sea-ice area 
    416416         ! 
    417417         sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rho0 
    418418         sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rho0 
    419419         ! 
    420          IF( .NOT.ln_linssh ) THEN 
    421             ! 
    422             WHERE( ht_0(:,:) > 0 )   ;   z2d(:,:) = 1._wp + sshn(:,:)*tmask(:,:,1) / ht_0(:,:) 
    423             ELSEWHERE                ;   z2d(:,:) = 1._wp   ;   END WHERE 
    424             ! 
    425             DO jk = 1,jpkm1                     ! adjust initial vertical scale factors                 
    426                e3t_n(:,:,jk) = e3t_0(:,:,jk) * z2d(:,:) 
    427                e3t_b(:,:,jk) = e3t_n(:,:,jk) 
    428                e3t_a(:,:,jk) = e3t_n(:,:,jk) 
    429             END DO 
    430             ! 
    431             ! Reconstruction of all vertical scale factors at now and before time-steps 
    432             ! ========================================================================= 
    433             ! Horizontal scale factor interpolations 
    434             ! -------------------------------------- 
    435             CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' ) 
    436             CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' ) 
    437             CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' ) 
    438             CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' ) 
    439             CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' ) 
    440             ! Vertical scale factor interpolations 
    441             ! ------------------------------------ 
    442             CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W'  ) 
    443             CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' ) 
    444             CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' ) 
    445             CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' ) 
    446             CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' ) 
    447             ! t- and w- points depth 
    448             ! ---------------------- 
    449             !!gm not sure of that.... 
    450             gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1) 
    451             gdepw_n(:,:,1) = 0.0_wp 
    452             gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:) 
    453             DO jk = 2, jpk 
    454                gdept_n(:,:,jk) = gdept_n(:,:,jk-1) + e3w_n(:,:,jk  ) 
    455                gdepw_n(:,:,jk) = gdepw_n(:,:,jk-1) + e3t_n(:,:,jk-1) 
    456                gde3w_n(:,:,jk) = gdept_n(:,:,jk  ) - sshn (:,:) 
    457             END DO 
     420         IF( .NOT.ln_linssh ) THEN     ! modified the now and before vertical mesh and scale factors  
     421            ! 
     422            !                             !* BEFORE fields :  
     423            CALL ssh2e3_before               ! set:      hu , hv , r1_hu, r1_hv  
     424            !                                    !  e3t, e3w, e3u, e3uw, e3v, e3vw 
     425            ! 
     426            !                             !* NOW fields :  
     427            CALL ssh2e3_now                  ! set: ht , hu , hv , r1_hu, r1_hv 
     428            !                                !      e3t, e3w, e3u, e3uw, e3v, e3vw, e3f 
     429            !                                !      gdept_n, gdepw_n, gde3w_n 
    458430         ENDIF 
    459431      ENDIF 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/ASM/asminc.F90

    r9939 r10001  
    5858    LOGICAL, PUBLIC, PARAMETER :: lk_asminc = .FALSE.  !: No assimilation increments 
    5959#endif 
    60    LOGICAL, PUBLIC :: ln_bkgwri     !: No output of the background state fields 
    61    LOGICAL, PUBLIC :: ln_asmiau     !: No applying forcing with an assimilation increment 
    62    LOGICAL, PUBLIC :: ln_asmdin     !: No direct initialization 
    63    LOGICAL, PUBLIC :: ln_trainc     !: No tracer (T and S) assimilation increments 
    64    LOGICAL, PUBLIC :: ln_dyninc     !: No dynamics (u and v) assimilation increments 
    65    LOGICAL, PUBLIC :: ln_sshinc     !: No sea surface height assimilation increment 
    66    LOGICAL, PUBLIC :: ln_seaiceinc  !: No sea ice concentration increment 
    67    LOGICAL, PUBLIC :: ln_salfix     !: Apply minimum salinity check 
    68    LOGICAL, PUBLIC :: ln_temnofreeze = .FALSE. !: Don't allow the temperature to drop below freezing 
    69    INTEGER, PUBLIC :: nn_divdmp     !: Apply divergence damping filter nn_divdmp times 
    70  
     60   LOGICAL, PUBLIC ::   ln_bkgwri     !: No output of the background state fields 
     61   LOGICAL, PUBLIC ::   ln_asmiau     !: No applying forcing with an assimilation increment 
     62   LOGICAL, PUBLIC ::   ln_asmdin     !: No direct initialization 
     63   LOGICAL, PUBLIC ::   ln_trainc     !: No tracer (T and S) assimilation increments 
     64   LOGICAL, PUBLIC ::   ln_dyninc     !: No dynamics (u and v) assimilation increments 
     65   LOGICAL, PUBLIC ::   ln_sshinc     !: No sea surface height assimilation increment 
     66   LOGICAL, PUBLIC ::   ln_seaiceinc  !: No sea ice concentration increment 
     67   LOGICAL, PUBLIC ::   ln_salfix     !: Apply minimum salinity check 
     68   LOGICAL, PUBLIC ::   ln_temnofreeze = .FALSE. !: Don't allow the temperature to drop below freezing 
     69   INTEGER, PUBLIC ::   nn_divdmp     !: Apply divergence damping filter nn_divdmp times 
     70 
     71   LOGICAL, PUBLIC ::   l_dynasm      !: flag derived from ln_asm... flags for calling  
     72   LOGICAL, PUBLIC ::   l_traasm      !: dyn_asm_inc and tra_asm_inc routines in step.F90 
     73    
    7174   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   t_bkg   , s_bkg      !: Background temperature and salinity 
    7275   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   u_bkg   , v_bkg      !: Background u- & v- velocity components 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DIA/diawri.F90

    r9939 r10001  
    118118      !  
    119119      ! Output the initial state and forcings 
    120       IF( ninist == 1 ) THEN                        
     120      IF( ln_istate ) THEN                        
    121121         CALL dia_wri_state( 'output.init', kt ) 
    122          ninist = 0 
     122         ln_istate = .FALSE. 
    123123      ENDIF 
    124124 
    … …  
    444444      IF( ln_timing )   CALL timing_start('dia_wri') 
    445445      ! 
    446       IF( ninist == 1 ) THEN     !==  Output the initial state and forcings  ==! 
     446      IF( ln_istate ) THEN     !==  Output the initial state and forcings  ==! 
    447447         CALL dia_wri_state( 'output.init', kt ) 
    448          ninist = 0 
     448         ln_istate = .FALSE. 
    449449      ENDIF 
    450450      ! 
    … …  
    878878      !! 
    879879      !! ** Method  :   NetCDF files using ioipsl 
    880       !!      File 'output.init.nc'  is created if ninist = 1 (namelist) 
     880      !!      File 'output.init.nc'  is created if ln_istate = T (namelist) 
    881881      !!      File 'output.abort.nc' is created in case of abnormal job end 
    882882      !!---------------------------------------------------------------------- 
    … …  
    10291029      CALL histclo( id_i ) 
    10301030#if ! defined key_iomput 
    1031       IF( ninist /= 1  ) THEN 
     1031      IF( ln_istate  ) THEN 
    10321032         CALL histclo( nid_T ) 
    10331033         CALL histclo( nid_U ) 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/dom_oce.F90

    r9939 r10001  
    3030   !! ---------------------------- 
    3131   !                                   !!* Namelist namdom : time & space domain * 
     32   LOGICAL , PUBLIC ::   ln_RK3         !: 2rd order Runge-Kutta scheme  
     33   LOGICAL , PUBLIC ::   ln_MLF         !: Modified Leap-Frog    scheme (Leclair & Madec OM2010) 
     34   REAL(wp), PUBLIC ::      rn_atfp        !: asselin time filter parameter 
     35   REAL(wp), PUBLIC ::   rn_Dt          !: time step for the dynamics and tracer (for both RK3 & MLF) 
    3236   LOGICAL , PUBLIC ::   ln_linssh      !: =T  linear free surface ==>> model level are fixed in time 
    3337   LOGICAL , PUBLIC ::   ln_meshmask    !: =T  create a mesh-mask file (mesh_mask.nc) 
    3438   REAL(wp), PUBLIC ::   rn_isfhmin     !: threshold to discriminate grounded ice to floating ice 
    35    REAL(wp), PUBLIC ::   rn_dt          !: time step for the dynamics and tracer 
    36    REAL(wp), PUBLIC ::   rn_atfp        !: asselin time filter parameter 
    3739   LOGICAL , PUBLIC ::   ln_1st_euler   !: =0 start with forward time step or not (=1) 
    3840   LOGICAL , PUBLIC ::   ln_iscpl       !: coupling with ice sheet 
    … …  
    143145   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hu_0  ,    hu_b ,    hu_n ,    hu_a   !: u-depth              [m] 
    144146   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hv_0  ,    hv_b ,    hv_n ,    hv_a   !: v-depth              [m] 
    145    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::           r1_hu_b , r1_hu_n , r1_hu_a   !: inverse of u-depth [1/m] 
    146    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::           r1_hv_b , r1_hv_n , r1_hv_a   !: inverse of v-depth [1/m] 
     147   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hf_0                                  !: v-depth              [m] 
     148   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_ht_0                                 !: inverse of u-depth [1/m] 
     149   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_hu_0 , r1_hu_b , r1_hu_n , r1_hu_a   !: inverse of u-depth [1/m] 
     150   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_hv_0 , r1_hv_b , r1_hv_n , r1_hv_a   !: inverse of v-depth [1/m] 
     151   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_hf_0                                 !: inverse of v-depth [1/m] 
    147152 
    148153   INTEGER, PUBLIC ::   nla10              !: deepest    W level Above  ~10m (nlb10 - 1) 
    149154   INTEGER, PUBLIC ::   nlb10              !: shallowest W level Bellow ~10m (nla10 + 1)  
    150155 
    151    !! 1D reference  vertical coordinate 
    152    !! =-----------------====------ 
     156   !! 1D reference vertical coordinate 
     157   !! ==------------------------------ 
    153158   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)   ::   gdept_1d, gdepw_1d !: reference depth of t- and w-points (m) 
    154159   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)   ::   e3t_1d  , e3w_1d   !: reference vertical scale factors at T- and W-pts (m) 
    … …  
    170175   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   risfdep                 !: Iceshelf draft                    (ISF) 
    171176 
    172    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ssmask, ssumask, ssvmask             !: surface mask at T-,U-, V- and F-pts 
     177   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ssmask, ssumask, ssvmask, ssfmask    !: surface mask at T-,U-, V- and F-pts 
    173178   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: tmask, umask, vmask, fmask   !: land/ocean mask at T-, U-, V- and F-pts 
    174179   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: wmask, wumask, wvmask        !: land/ocean mask at WT-, WU- and WV-pts 
    … …  
    265270         &      e3uw_n(jpi,jpj,jpk) , e3vw_n(jpi,jpj,jpk) ,     STAT=ierr(5) )                        
    266271         ! 
    267       ALLOCATE( ht_0(jpi,jpj) , hu_0(jpi,jpj) , hv_0(jpi,jpj) ,                                           & 
    268          &                      hu_b(jpi,jpj) , hv_b(jpi,jpj) , r1_hu_b(jpi,jpj) , r1_hv_b(jpi,jpj) ,     & 
    269          &      ht_n(jpi,jpj) , hu_n(jpi,jpj) , hv_n(jpi,jpj) , r1_hu_n(jpi,jpj) , r1_hv_n(jpi,jpj) ,     & 
    270          &                      hu_a(jpi,jpj) , hv_a(jpi,jpj) , r1_hu_a(jpi,jpj) , r1_hv_a(jpi,jpj) , STAT=ierr(6)  ) 
    271          ! 
     272      ALLOCATE( ht_0(jpi,jpj)    , hu_0(jpi,jpj)   , hv_0(jpi,jpj)    , hf_0(jpi,jpj)    ,                     & 
     273         &                         hu_b(jpi,jpj)   , hv_b(jpi,jpj)    , r1_hu_b(jpi,jpj) , r1_hv_b(jpi,jpj) ,  & 
     274         &      ht_n(jpi,jpj)    , hu_n(jpi,jpj)   , hv_n(jpi,jpj)    , r1_hu_n(jpi,jpj) , r1_hv_n(jpi,jpj) ,  & 
     275         &                         hu_a(jpi,jpj)   , hv_a(jpi,jpj)    , r1_hu_a(jpi,jpj) , r1_hv_a(jpi,jpj) ,  & 
     276         &      r1_ht_0(jpi,jpj) , r1_hu_0(jpi,jpj), r1_hv_0(jpi,jpj) , r1_hf_0(jpi,jpj) , STAT=ierr(6)  ) 
    272277         ! 
    273278      ALLOCATE( gdept_1d(jpk) , gdepw_1d(jpk) , e3t_1d(jpk) , e3w_1d(jpk) , STAT=ierr(7) ) 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/domain.F90

    r9939 r10001  
    7878      CHARACTER (len=*), INTENT(IN) :: cdstr                  ! model: NEMO or SAS. Determines core restart variables 
    7979      INTEGER , DIMENSION(jpi,jpj) ::   ik_top , ik_bot       ! top and bottom ocean level 
    80       REAL(wp), DIMENSION(jpi,jpj) ::   z1_hu_0, z1_hv_0 
    8180      !!---------------------------------------------------------------------- 
    8281      ! 
    … …  
    138137      CALL dom_msk( ik_top, ik_bot )   ! Masks 
    139138      IF( ln_closea )   CALL dom_clo   ! ln_closea=T : closed seas included in the simulation 
    140                                        ! Read in masks to define closed seas and lakes  
     139      !                                ! Read in masks to define closed seas and lakes  
    141140      ! 
    142141      DO jj = 1, jpj                   ! depth of the iceshelves 
    143142         DO ji = 1, jpi 
    144143            ik = mikt(ji,jj) 
    145             risfdep(ji,jj) = gdepw_0(ji,jj,ik) 
     144            risfdep(ji,jj) = gdepw_0(ji,jj,ik)        !!gm  RENAME it as h_isf(:,:)  better no? 
    146145         END DO 
    147146      END DO 
    148147      ! 
    149       ht_0(:,:) = 0._wp  ! Reference ocean thickness 
    150       hu_0(:,:) = 0._wp 
     148      ht_0(:,:) = 0._wp  ! Reference water column thickness  (wet point only, i.e. without ISF thickness 
     149      hu_0(:,:) = 0._wp  !           ------------                                  and zero over land    )     
    151150      hv_0(:,:) = 0._wp 
     151      hf_0(:,:) = 0._wp 
    152152      DO jk = 1, jpk 
    153153         ht_0(:,:) = ht_0(:,:) + e3t_0(:,:,jk) * tmask(:,:,jk) 
    154154         hu_0(:,:) = hu_0(:,:) + e3u_0(:,:,jk) * umask(:,:,jk) 
    155155         hv_0(:,:) = hv_0(:,:) + e3v_0(:,:,jk) * vmask(:,:,jk) 
     156         hf_0(:,:) = hf_0(:,:) + e3f_0(:,:,jk) * fmask(:,:,jk) 
    156157      END DO 
    157158      ! 
     159      r1_ht_0(:,:) =  ssmask(:,:) / ( ht_0(:,:) + 1._wp -  ssmask(:,:) )    ! =0 on lands 
     160      r1_hu_0(:,:) = ssumask(:,:) / ( hu_0(:,:) + 1._wp - ssumask(:,:) )    ! not in ISF areas 
     161      r1_hv_0(:,:) = ssvmask(:,:) / ( hv_0(:,:) + 1._wp - ssvmask(:,:) ) 
     162      r1_hf_0(:,:) = ssfmask(:,:) / ( hf_0(:,:) + 1._wp - ssfmask(:,:) ) 
     163      ! 
    158164      !           !==  time varying part of coordinate system  ==! 
    159165      ! 
    160166      IF( ln_linssh ) THEN       != Fix in time : set to the reference one for all 
    161       ! 
    162          !       before        !          now          !       after         ! 
    163             gdept_b = gdept_0  ;   gdept_n = gdept_0   !        ---          ! depth of grid-points 
    164             gdepw_b = gdepw_0  ;   gdepw_n = gdepw_0   !        ---          ! 
    165                                    gde3w_n = gde3w_0   !        ---          ! 
    166          !                                                                   
    167               e3t_b =   e3t_0  ;     e3t_n =   e3t_0   ;   e3t_a =  e3t_0    ! scale factors 
    168               e3u_b =   e3u_0  ;     e3u_n =   e3u_0   ;   e3u_a =  e3u_0    ! 
    169               e3v_b =   e3v_0  ;     e3v_n =   e3v_0   ;   e3v_a =  e3v_0    ! 
    170                                      e3f_n =   e3f_0   !        ---          ! 
    171               e3w_b =   e3w_0  ;     e3w_n =   e3w_0   !        ---          ! 
    172              e3uw_b =  e3uw_0  ;    e3uw_n =  e3uw_0   !        ---          ! 
    173              e3vw_b =  e3vw_0  ;    e3vw_n =  e3vw_0   !        ---          ! 
    174          ! 
    175          z1_hu_0(:,:) = ssumask(:,:) / ( hu_0(:,:) + 1._wp - ssumask(:,:) )     ! _i mask due to ISF 
    176          z1_hv_0(:,:) = ssvmask(:,:) / ( hv_0(:,:) + 1._wp - ssvmask(:,:) ) 
    177          ! 
    178          !        before       !          now          !       after         ! 
    179                                       ht_n =    ht_0   !                     ! water column thickness 
    180                hu_b =    hu_0  ;      hu_n =    hu_0   ;    hu_a =    hu_0   !  
    181                hv_b =    hv_0  ;      hv_n =    hv_0   ;    hv_a =    hv_0   ! 
    182             r1_hu_b = z1_hu_0  ;   r1_hu_n = z1_hu_0   ; r1_hu_a = z1_hu_0   ! inverse of water column thickness 
    183             r1_hv_b = z1_hv_0  ;   r1_hv_n = z1_hv_0   ; r1_hv_a = z1_hv_0   ! 
     167         ! 
     168         !    before         !          now          !       after         ! 
     169         gdept_b = gdept_0   ;   gdept_n = gdept_0   !        ---          !    depth  
     170         gdepw_b = gdepw_0   ;   gdepw_n = gdepw_0   !        ---          !     of 
     171                                 gde3w_n = gde3w_0   !        ---          ! grid-points 
     172         !                   !                       !                     !                                               
     173          e3t_b  =   e3t_0   ;     e3t_n =   e3t_0   ;   e3t_a =  e3t_0    !  scale  
     174          e3u_b  =   e3u_0   ;     e3u_n =   e3u_0   ;   e3u_a =  e3u_0    ! factors 
     175          e3v_b  =   e3v_0   ;     e3v_n =   e3v_0   ;   e3v_a =  e3v_0    ! 
     176                                   e3f_n =   e3f_0   !        ---          ! 
     177          e3w_b  =   e3w_0   ;     e3w_n =   e3w_0   !        ---          ! 
     178         e3uw_b  =  e3uw_0   ;    e3uw_n =  e3uw_0   !        ---          ! 
     179         e3vw_b  =  e3vw_0   ;    e3vw_n =  e3vw_0   !        ---          ! 
     180         !                   ! 
     181                                    ht_n =    ht_0   !                     ! water column 
     182            hu_b =    hu_0   ;      hu_n =    hu_0   ;    hu_a =    hu_0   !  thickness 
     183            hv_b =    hv_0   ;      hv_n =    hv_0   ;    hv_a =    hv_0   ! 
     184         r1_hu_b = r1_hu_0   ;   r1_hu_n = r1_hu_0   ; r1_hu_a = r1_hu_0   ! 1 / water column 
     185         r1_hv_b = r1_hv_0   ;   r1_hv_n = r1_hv_0   ; r1_hv_a = r1_hv_0   !       thickness 
    184186         ! 
    185187         ! 
    … …  
    192194      IF( lk_c1d         )   CALL cor_c1d       ! 1D configuration: Coriolis set at T-point 
    193195      ! 
    194       IF( ln_meshmask .AND. .NOT.ln_iscpl )                        CALL dom_wri     ! Create a domain file 
     196      IF( ln_meshmask .AND. .NOT.ln_iscpl                      )   CALL dom_wri     ! Create a domain file 
    195197      IF( ln_meshmask .AND.      ln_iscpl .AND. .NOT.ln_rstart )   CALL dom_wri     ! Create a domain file 
    196198      IF(                                       .NOT.ln_rstart )   CALL dom_ctl     ! Domain control 
    197199      ! 
    198       IF( ln_write_cfg )   CALL cfg_write         ! create the configuration file 
     200      IF( ln_write_cfg )   CALL cfg_write       ! create the configuration file 
    199201      ! 
    200202      IF(lwp) THEN 
    … …  
    204206         WRITE(numout,*)  
    205207      ENDIF 
    206       ! 
     208       ! 
    207209   END SUBROUTINE dom_init 
    208210 
    … …  
    290292      NAMELIST/namrun/ cn_ocerst_indir, cn_ocerst_outdir, nn_stocklist, ln_rst_list,                   & 
    291293         &             nn_no   , cn_exp   , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl   ,     & 
    292          &             nn_it000, nn_itend , nn_date0    , nn_time0     , nn_leapy  , nn_istate   ,     & 
     294         &             nn_it000, nn_itend , nn_date0    , nn_time0     , nn_leapy  , ln_istate   ,     & 
    293295         &             nn_stock, nn_write , ln_mskland  , ln_clobber   , nn_chunksz, ln_1st_euler,     & 
    294296         &             ln_cfmeta, ln_iscpl, ln_xios_read, nn_wxios 
    295       NAMELIST/namdom/ ln_linssh, rn_isfhmin, rn_Dt, rn_atfp, ln_crs, ln_meshmask 
     297      NAMELIST/namdom/ ln_RK3, ln_MLF, rn_Dt, rn_atfp, ln_linssh, rn_isfhmin, ln_crs, ln_meshmask 
    296298#if defined key_netcdf4 
    297299      NAMELIST/namnc4/ nn_nchunks_i, nn_nchunks_j, nn_nchunks_k, ln_nc4zip 
    … …  
    305307      ENDIF 
    306308      ! 
     309      REWIND( numnam_ref )              ! Namelist namdom in reference namelist : space & time domain (bathymetry, mesh, timestep) 
     310      READ  ( numnam_ref, namdom, IOSTAT = ios, ERR = 903) 
     311903   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdom in reference namelist', lwp ) 
     312      REWIND( numnam_cfg )              ! Namelist namdom in configuration namelist : space & time domain (bathymetry, mesh, timestep) 
     313      READ  ( numnam_cfg, namdom, IOSTAT = ios, ERR = 904 ) 
     314904   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdom in configuration namelist', lwp ) 
     315      IF(lwm) WRITE( numond, namdom ) 
     316      ! 
     317      IF(lwp) THEN 
     318         WRITE(numout,*) 
     319         WRITE(numout,*) '   Namelist : namdom   ---   space & time domain' 
     320         WRITE(numout,*) '      3rd order Runge-Kutta scheme            ln_RK3      = ', ln_RK3 
     321         WRITE(numout,*) '      Modified Leap-Frog    scheme            ln_MLF      = ', ln_MLF 
     322         WRITE(numout,*) '         asselin time filter parameter           rn_atfp  = ', rn_atfp 
     323         WRITE(numout,*) '      ocean time step (RK3 & MLF)             rn_Dt       = ', rn_Dt 
     324         WRITE(numout,*) '      linear free surface (=T)                ln_linssh   = ', ln_linssh 
     325         WRITE(numout,*) '      create mesh/mask file                   ln_meshmask = ', ln_meshmask 
     326         WRITE(numout,*) '      treshold to open the isf cavity         rn_isfhmin  = ', rn_isfhmin, ' [m]' 
     327         WRITE(numout,*) '      online coarsening of dynamical fields   ln_crs      = ', ln_crs 
     328         WRITE(numout,*)   
     329      ENDIF 
     330      ! 
     331      IF ( (      ln_RK3 .AND. .NOT.ln_MLF ) .OR.   & 
     332         & ( .NOT.ln_RK3 .AND.      ln_MLF ) ) THEN 
     333         IF( ln_RK3 .AND. lwp )   WRITE(numout,*)'   ==>>>   a RK3 time-stepping scheme is used' 
     334         IF( ln_MLF .AND. lwp )   WRITE(numout,*)'   ==>>>   a MLF time-stepping scheme is used' 
     335      ELSE 
     336         CALL ctl_stop( 'dom_nam: select one and only one time stepping scheme') 
     337      ENDIF 
    307338      ! 
    308339      REWIND( numnam_ref )              ! Namelist namrun in reference namelist : Parameters of the run 
    … …  
    315346      ! 
    316347      IF(lwp) THEN                  ! control print 
     348         WRITE(numout,*) 
    317349         WRITE(numout,*) '   Namelist : namrun   ---   run parameters' 
    318350         WRITE(numout,*) '      Assimilation cycle              nn_no           = ', nn_no 
    … …  
    330362         WRITE(numout,*) '      initial time of day in hhmm     nn_time0        = ', nn_time0 
    331363         WRITE(numout,*) '      leap year calendar (0/1)        nn_leapy        = ', nn_leapy 
    332          WRITE(numout,*) '      initial state output            nn_istate       = ', nn_istate 
     364         WRITE(numout,*) '      output the initial state        ln_istate       = ', ln_istate 
    333365         IF( ln_rst_list ) THEN 
    334366            WRITE(numout,*) '      list of restart dump times      nn_stocklist    =', nn_stocklist 
    … …  
    357389      ndate0 = nn_date0 
    358390      nleapy = nn_leapy 
    359       ninist = nn_istate 
    360391      nstock = nn_stock 
    361392      nstocklist = nn_stocklist 
    362393      nwrite = nn_write 
    363       IF( ln_rstart ) THEN       ! restart : set 1st time-step scheme (LF or forward)  
    364          l_1st_euler = ln_1st_euler 
    365       ELSE                       ! start from rest : always an Euler scheme for the 1st time-step 
    366          IF(lwp) WRITE(numout,*)   
    367          IF(lwp) WRITE(numout,*)'   ==>>>   Start from rest (ln_rstart=F)' 
    368          IF(lwp) WRITE(numout,*)'           an Euler initial time step is used '    
    369          l_1st_euler = .TRUE. 
     394       
     395      IF( ln_MLF ) THEN       ! Leap-Frog only 
     396         IF( ln_rstart ) THEN       ! restart : set 1st time-step scheme (LF or forward)  
     397            l_1st_euler = ln_1st_euler 
     398         ELSE                       ! start from rest : always an Euler scheme for the 1st time-step 
     399            IF(lwp) WRITE(numout,*)   
     400            IF(lwp) WRITE(numout,*)'   ==>>>   Start from rest (ln_rstart=F)' 
     401            IF(lwp) WRITE(numout,*)'           an Euler initial time step is used '    
     402            l_1st_euler = .TRUE. 
     403         ENDIF 
    370404      ENDIF 
    371405      !                             ! control of output frequency 
    … …  
    399433      ENDIF 
    400434#endif 
    401  
    402       REWIND( numnam_ref )              ! Namelist namdom in reference namelist : space & time domain (bathymetry, mesh, timestep) 
    403       READ  ( numnam_ref, namdom, IOSTAT = ios, ERR = 903) 
    404 903   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdom in reference namelist', lwp ) 
    405       REWIND( numnam_cfg )              ! Namelist namdom in configuration namelist : space & time domain (bathymetry, mesh, timestep) 
    406       READ  ( numnam_cfg, namdom, IOSTAT = ios, ERR = 904 ) 
    407 904   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdom in configuration namelist', lwp ) 
    408       IF(lwm) WRITE( numond, namdom ) 
    409       ! 
    410       IF(lwp) THEN 
    411          WRITE(numout,*) 
    412          WRITE(numout,*) '   Namelist : namdom   ---   space & time domain' 
    413          WRITE(numout,*) '      linear free surface (=T)                ln_linssh   = ', ln_linssh 
    414          WRITE(numout,*) '      create mesh/mask file                   ln_meshmask = ', ln_meshmask 
    415          WRITE(numout,*) '      treshold to open the isf cavity         rn_isfhmin  = ', rn_isfhmin, ' [m]' 
    416          WRITE(numout,*) '      ocean time step                         rn_dt       = ', rn_dt 
    417          WRITE(numout,*) '      asselin time filter parameter           rn_atfp     = ', rn_atfp 
    418          WRITE(numout,*) '      online coarsening of dynamical fields   ln_crs      = ', ln_crs 
    419       ENDIF 
    420       ! 
    421       IF( TRIM(Agrif_CFixed()) == '0' ) THEN 
     435      ! 
     436      IF( TRIM(Agrif_CFixed()) == '0' ) THEN    !set output file type for XIOS based on NEMO namelist  
    422437         lrxios = ln_xios_read.AND.ln_rstart 
    423 !set output file type for XIOS based on NEMO namelist  
    424438         IF (nn_wxios > 0) lwxios = .TRUE.  
    425439         nxioso = nn_wxios 
    … …  
    553567            CALL iom_getatt( inum, 'cn_cfg', cd_cfg )  ! returns   !  if not found 
    554568            CALL iom_getatt( inum, 'nn_cfg', kk_cfg )  ! returns -999 if not found 
    555             IF( TRIM(cd_cfg) == '!') cd_cfg = 'UNKNOWN' 
    556             IF( kk_cfg == -999     ) kk_cfg = -9999999 
     569            IF( TRIM(cd_cfg) == '!')   cd_cfg = 'UNKNOWN' 
     570            IF( kk_cfg == -999     )   kk_cfg = -9999999 
    557571         ENDIF 
    558572         ! 
    … …  
    701715      ! 
    702716   END SUBROUTINE cfg_write 
    703  
     717    
    704718   !!====================================================================== 
    705719END MODULE domain 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/dommsk.F90

    r9657 r10001  
    202202      ssumask(:,:) = MAXVAL( umask(:,:,:), DIM=3 ) 
    203203      ssvmask(:,:) = MAXVAL( vmask(:,:,:), DIM=3 ) 
     204      ssfmask(:,:) = MAXVAL( fmask(:,:,:), DIM=3 ) 
    204205 
    205206 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/domvvl.F90

    r9939 r10001  
    44   !! Ocean :  
    55   !!====================================================================== 
    6    !! History :  2.0  !  2006-06  (B. Levier, L. Marie)  original code 
    7    !!            3.1  !  2009-02  (G. Madec, M. Leclair, R. Benshila)  pure z* coordinate 
    8    !!            3.3  !  2011-10  (M. Leclair) totally rewrote domvvl: vvl option includes z_star and z_tilde coordinates 
    9    !!            3.6  !  2014-11  (P. Mathiot) add ice shelf capability 
     6   !! History :  5.0  !  2018-07  (G. Madec)  Flux Form with Kinetic Energy conservation 
     7   !!                             ==>>>  here z* and s* only (no z-tilde)  
     8    
     9   ! 1- remove   z-tilde          ==>>>  pure z-star (or s-star) 
     10   ! 2- remove   dom_vvl_interpol   
     11    
    1012   !!---------------------------------------------------------------------- 
    1113 
    … …  
    1416   !!   dom_vvl_sf_nxt   : Compute next vertical scale factors 
    1517   !!   dom_vvl_sf_swp   : Swap vertical scale factors and update the vertical grid 
    16    !!   dom_vvl_interpol : Interpolate vertical scale factors from one grid point to another 
    1718   !!   dom_vvl_rst      : read/write restart file 
    1819   !!   dom_vvl_ctl      : Check the vvl options 
    … …  
    3839   PUBLIC  dom_vvl_sf_nxt     ! called by step.F90 
    3940   PUBLIC  dom_vvl_sf_swp     ! called by step.F90 
    40    PUBLIC  dom_vvl_interpol   ! called by dynnxt.F90 
     41   PUBLIC  ssh2e3_before      ! ... 
     42   PUBLIC  ssh2e3_now         ! ... 
    4143 
    4244   !                                                      !!* Namelist nam_vvl 
    … …  
    6163   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_hdv       ! retoring period for low freq. divergence 
    6264 
     65 
     66!!gm add 
     67!!gm  
     68 
     69 
     70 
    6371   !! * Substitutions 
    6472#  include "vectopt_loop_substitute.h90" 
    … …  
    7482      !!                ***  FUNCTION dom_vvl_alloc  *** 
    7583      !!---------------------------------------------------------------------- 
     84      ! 
    7685      IF( ln_vvl_zstar )   dom_vvl_alloc = 0 
    7786      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 
    … …  
    115124      !! Reference  : Leclair, M., and G. Madec, 2011, Ocean Modelling. 
    116125      !!---------------------------------------------------------------------- 
    117       INTEGER ::   ji, jj, jk 
    118       INTEGER ::   ii0, ii1, ij0, ij1 
    119       REAL(wp)::   zcoef, z1_Dt 
    120126      !!---------------------------------------------------------------------- 
    121127      ! 
    … …  
    129135      IF( dom_vvl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' ) 
    130136      ! 
    131       !                    ! Read or initialize e3t_(b/n), te3t_(b/n) and hdiv_lf 
     137      !                    ! Read or initialize e3t_(b/n), ssh(b/n) 
    132138      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       ! 
    135       !                    !== Set of all other vertical scale factors  ==!  (now and before) 
    136       !                                ! 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 
    142       !                                ! 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' ) 
    149  
    150       ! 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       ! 
    155       !                    !==  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 
    161       DO jk = 2, jpk                               ! vertical sum 
    162          DO jj = 1,jpj 
    163             DO ji = 1,jpi 
    164                !    zcoef = tmask - wmask    ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 
    165                !                             ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) 
    166                !                             ! 0.5 where jk = mikt      
    167 !!gm ???????   BUG ?  gdept_n as well as gde3w_n  does not include the thickness of ISF ?? 
    168                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))  
    176             END DO 
    177          END DO 
    178       END DO 
    179       ! 
    180       !                    !==  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) 
    186       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) 
    192       END DO 
    193       ! 
    194       !                    !==  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(:,:) ) 
    199  
    200       !                    !==   z_tilde coordinate case  ==!   (Restoring frequencies) 
    201       IF( ln_vvl_ztilde ) THEN 
    202 !!gm : idea: add here a READ in a file of custumized restoring frequency 
    203          !                                   ! Values in days provided via the namelist 
    204          !                                   ! use rsmall to avoid possible division by zero errors with faulty settings 
    205          frq_rst_e3t(:,:) = 2._wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400.0_wp ) 
    206          frq_rst_hdv(:,:) = 2._wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp ) 
    207          ! 
    208          IF( ln_vvl_ztilde_as_zstar ) THEN   ! z-star emulation using z-tile 
    209             frq_rst_e3t(:,:) = 0._wp               !Ignore namelist settings 
    210             frq_rst_hdv(:,:) = 1._wp / rn_Dt 
    211          ENDIF 
    212          IF ( ln_vvl_zstar_at_eqtor ) THEN   ! use z-star in vicinity of the Equator 
    213             z1_Dt = 1._wp / rn_Dt 
    214             DO jj = 1, jpj 
    215                DO ji = 1, jpi 
    216 !!gm  case |gphi| >= 6 degrees is useless   initialized just above by default 
    217                   IF( ABS(gphit(ji,jj)) >= 6.) THEN 
    218                      ! values outside the equatorial band and transition zone (ztilde) 
    219                      frq_rst_e3t(ji,jj) =  2._wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400._wp ) 
    220                      frq_rst_hdv(ji,jj) =  2._wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400._wp ) 
    221                   ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN    ! Equator strip ==> z-star 
    222                      ! values inside the equatorial band (ztilde as zstar) 
    223                      frq_rst_e3t(ji,jj) =  0._wp 
    224                      frq_rst_hdv(ji,jj) =  z1_Dt 
    225                   ELSE                                      ! transition band (2.5 to 6 degrees N/S) 
    226                      !                                      ! (linearly transition from z-tilde to z-star) 
    227                      frq_rst_e3t(ji,jj) = 0._wp + ( frq_rst_e3t(ji,jj) - 0._wp  ) * 0.5_wp                             & 
    228                         &                       * (  1._wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) * 180._wp / 3.5_wp )  ) 
    229                      frq_rst_hdv(ji,jj) = z1_Dt + (  frq_rst_hdv(ji,jj) - z1_Dt ) * 0.5_wp                             & 
    230                         &                       * (  1._wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) * 180._wp / 3.5_wp )  ) 
    231                   ENDIF 
    232                END DO 
    233             END DO 
    234             IF( cn_cfg == "orca" .AND. nn_cfg == 3 ) THEN   ! ORCA2: Suppress ztilde in the Foxe Basin for ORCA2 
    235                ii0 = 103   ;   ii1 = 111        
    236                ij0 = 128   ;   ij1 = 135   ;    
    237                frq_rst_e3t( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  0._wp 
    238                frq_rst_hdv( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  z1_Dt 
    239             ENDIF 
    240          ENDIF 
    241       ENDIF 
    242       ! 
    243       IF(lwxios) THEN 
    244 ! define variables in restart file when writing with XIOS 
     139      ! 
     140      !                    !== Set of all other vertical mesh fields  ==!  (now and before)       
     141      ! 
     142      !                          !* BEFORE fields :  
     143      CALL ssh2e3_before               ! set:      hu , hv , r1_hu, r1_hv  
     144      !                                    !  e3t, e3w, e3u, e3uw, e3v, e3vw 
     145      ! 
     146      !                                ! set one for all last level to the e3._0 value 
     147      e3t_b(:,:,jpk) = e3t_0(:,:,jpk)  ;   e3u_b(:,:,jpk) =  e3w_0(:,:,jpk)  ;   e3v_b(:,:,jpk) =  e3v_0(:,:,jpk) 
     148      e3w_b(:,:,jpk) = e3w_0(:,:,jpk)  ;  e3uw_b(:,:,jpk) = e3uw_0(:,:,jpk)  ;  e3vw_b(:,:,jpk) = e3vw_0(:,:,jpk) 
     149      ! 
     150      !                          !* NOW fields :  
     151      CALL ssh2e3_now                  ! set: ht , hu , hv , r1_hu, r1_hv 
     152      !                                !      e3t, e3w, e3u, e3uw, e3v, e3vw, e3f 
     153      !                                !      gdept_n, gdepw_n, gde3w_n 
     154      ! 
     155      !                                ! set one for all last level to the e3._0 value 
     156      e3t_n(:,:,jpk) = e3t_0(:,:,jpk)  ;   e3u_n(:,:,jpk) =  e3w_0(:,:,jpk)  ;   e3v_n(:,:,jpk) =  e3v_0(:,:,jpk) 
     157      e3w_n(:,:,jpk) = e3w_0(:,:,jpk)  ;  e3uw_n(:,:,jpk) = e3uw_0(:,:,jpk)  ;  e3vw_n(:,:,jpk) = e3vw_0(:,:,jpk) 
     158      e3f_n(:,:,jpk) = e3f_0(:,:,jpk) 
     159      ! 
     160      !                          !* AFTER fields : (last level for OPA, 3D required for AGRIF initialisation) 
     161      e3t_a(:,:,:) = e3t_n(:,:,:)   ;   e3u_a(:,:,:) = e3u_n(:,:,:)   ;   e3v_a(:,:,:) = e3v_n(:,:,:) 
     162       
     163!!gm            
     164!!gm        ===>>>>   below:  some issues to think about !!! 
     165!!gm   
     166!!gm   fmask definition checked (0 or 1 nothing else) 
     167!!        in z-tilde or ALE    e3._0  should be the time varying fields step forward with an euler scheme 
     168!!gm   e3w_b  & gdept_b are not used   except its update in agrif   
     169!!     and used to compute before slope of surface in dynldf_iso ==>>>  remove it !!!! 
     170!!         NB: in triads on TRA, gdept_n  is used !!!!   BUG? 
     171!!gm   e3f_n  almost not used  ===>>>>  verify whether it can be removed or not... 
     172!!       verify the use of wumask & wvmask   mau be replaced by a multiplication by umask(k)*umask(k+1)  ??? 
     173!! 
     174!!gm  ISF case to be checked by Pierre Mathiot 
     175!! 
     176!!gm  setting of e3._a  for agrif....  TO BE CHECKED 
     177 
     178      ! 
     179      IF(lwxios) THEN      ! define variables in restart file when writing with XIOS 
    245180         CALL iom_set_rstw_var_active('e3t_b') 
    246181         CALL iom_set_rstw_var_active('e3t_n') 
    247          !                                           ! ----------------------- ! 
    248          IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN  ! z_tilde and layer cases ! 
    249             !                                        ! ----------------------- ! 
    250             CALL iom_set_rstw_var_active('tilde_e3t_b') 
    251             CALL iom_set_rstw_var_active('tilde_e3t_n') 
    252          END IF 
    253          !                                           ! -------------!     
    254          IF( ln_vvl_ztilde ) THEN                    ! z_tilde case ! 
    255             !                                        ! ------------ ! 
    256             CALL iom_set_rstw_var_active('hdiv_lf') 
    257          ENDIF 
    258          ! 
    259182      ENDIF 
    260183      ! 
    … …  
    287210      INTEGER, INTENT( in ), OPTIONAL ::   kcall   ! optional argument indicating call sequence 
    288211      ! 
    289       INTEGER                ::   ji, jj, jk            ! dummy loop indices 
    290       INTEGER , DIMENSION(3) ::   ijk_max, ijk_min      ! temporary integers 
    291       REAL(wp)               ::   z2dt, z_tmin, z_tmax  ! local scalars 
    292       LOGICAL                ::   ll_do_bclinic         ! local logical 
    293       REAL(wp), DIMENSION(jpi,jpj)     ::   zht, z_scale, zwu, zwv, zhdiv 
    294       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ze3t 
     212      INTEGER ::   ji, jj, jk   ! dummy loop indices 
     213      REAL(wp), DIMENSION(jpi,jpj) ::   zssht_h, zsshu_h, zsshv_h 
    295214      !!---------------------------------------------------------------------- 
    296215      ! 
    … …  
    305224      ENDIF 
    306225 
    307       ll_do_bclinic = .TRUE. 
    308       IF( PRESENT(kcall) ) THEN 
    309          IF( kcall == 2 .AND. ln_vvl_ztilde )   ll_do_bclinic = .FALSE. 
    310       ENDIF 
    311  
    312       ! ******************************* ! 
    313       ! After acale factors at t-points ! 
    314       ! ******************************* ! 
    315       !                                                ! --------------------------------------------- ! 
    316       !                                                ! z_star coordinate and barotropic z-tilde part ! 
    317       !                                                ! --------------------------------------------- ! 
    318       ! 
    319       z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * ssmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) ) 
    320       DO jk = 1, jpkm1 
    321          ! formally this is the same as e3t_a = e3t_0*(1+ssha/ht_0) 
    322          e3t_a(:,:,jk) = e3t_b(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk) 
    323       END DO 
    324       ! 
    325       IF( ln_vvl_ztilde .OR. ln_vvl_layer .AND. ll_do_bclinic ) THEN   ! z_tilde or layer coordinate ! 
    326          !                                                            ! ------baroclinic part------ ! 
    327          ! I - initialization 
    328          ! ================== 
    329  
    330          ! 1 - barotropic divergence 
    331          ! ------------------------- 
    332          zhdiv(:,:) = 0._wp 
    333          zht(:,:)   = 0._wp 
    334          DO jk = 1, jpkm1 
    335             zhdiv(:,:) = zhdiv(:,:) + e3t_n(:,:,jk) * hdivn(:,:,jk) 
    336             zht  (:,:) = zht  (:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 
    337          END DO 
    338          zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) ) 
    339  
    340          ! 2 - Low frequency baroclinic horizontal divergence  (z-tilde case only) 
    341          ! -------------------------------------------------- 
    342          IF( ln_vvl_ztilde ) THEN 
    343             IF( kt > nit000 ) THEN 
    344                DO jk = 1, jpkm1 
    345                   hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rn_Dt * frq_rst_hdv(:,:)   & 
    346                      &          * ( hdiv_lf(:,:,jk) - e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) ) 
    347                END DO 
    348             ENDIF 
    349          ENDIF 
    350  
    351          ! II - after z_tilde increments of vertical scale factors 
    352          ! ======================================================= 
    353          te3t_a(:,:,:) = 0._wp  ! te3t_a used to store tendency terms 
    354  
    355          ! 1 - High frequency divergence term 
    356          ! ---------------------------------- 
    357          IF( ln_vvl_ztilde ) THEN     ! z_tilde case 
    358             DO jk = 1, jpkm1 
    359                te3t_a(:,:,jk) = te3t_a(:,:,jk) - ( e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) ) 
    360             END DO 
    361          ELSE                         ! layer case 
    362             DO jk = 1, jpkm1 
    363                te3t_a(:,:,jk) = te3t_a(:,:,jk) -   e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk) 
    364             END DO 
    365          ENDIF 
    366  
    367          ! 2 - Restoring term (z-tilde case only) 
    368          ! ------------------ 
    369          IF( ln_vvl_ztilde ) THEN 
    370             DO jk = 1, jpk 
    371                te3t_a(:,:,jk) = te3t_a(:,:,jk) - frq_rst_e3t(:,:) * te3t_b(:,:,jk) 
    372             END DO 
    373          ENDIF 
    374  
    375          ! 3 - Thickness diffusion term 
    376          ! ---------------------------- 
    377          zwu(:,:) = 0._wp 
    378          zwv(:,:) = 0._wp 
    379          DO jk = 1, jpkm1        ! a - first derivative: diffusive fluxes 
    380             DO jj = 1, jpjm1 
    381                DO ji = 1, fs_jpim1   ! vector opt. 
    382                   un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj)           & 
    383                      &            * ( te3t_b(ji,jj,jk) - te3t_b(ji+1,jj  ,jk) ) 
    384                   vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj)           &  
    385                      &            * ( te3t_b(ji,jj,jk) - te3t_b(ji  ,jj+1,jk) ) 
    386                   zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk) 
    387                   zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk) 
    388                END DO 
    389             END DO 
    390          END DO 
    391          DO jj = 1, jpj          ! b - correction for last oceanic u-v points 
    392             DO ji = 1, jpi 
    393                un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj) 
    394                vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj) 
    395             END DO 
    396          END DO 
    397          DO jk = 1, jpkm1        ! c - second derivative: divergence of diffusive fluxes 
    398             DO jj = 2, jpjm1 
    399                DO ji = fs_2, fs_jpim1   ! vector opt. 
    400                   te3t_a(ji,jj,jk) = te3t_a(ji,jj,jk) + (   un_td(ji-1,jj  ,jk) - un_td(ji,jj,jk)    & 
    401                      &                                          +     vn_td(ji  ,jj-1,jk) - vn_td(ji,jj,jk)    & 
    402                      &                                            ) * r1_e1e2t(ji,jj) 
    403                END DO 
    404             END DO 
    405          END DO 
    406          !                       ! d - thickness diffusion transport: boundary conditions 
    407          !                             (stored for tracer advction and continuity equation) 
    408          CALL lbc_lnk_multi( un_td , 'U' , -1._wp, vn_td , 'V' , -1._wp) 
    409  
    410          ! 4 - Time stepping of baroclinic scale factors 
    411          ! --------------------------------------------- 
    412          ! Leapfrog time stepping 
    413          ! ~~~~~~~~~~~~~~~~~~~~~~ 
    414          CALL lbc_lnk( te3t_a(:,:,:), 'T', 1._wp ) 
    415          te3t_a(:,:,:) = te3t_b(:,:,:) + z2dt * tmask(:,:,:) * te3t_a(:,:,:) 
    416  
    417          ! Maximum deformation control 
    418          ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
    419          ze3t(:,:,jpk) = 0._wp 
    420          DO jk = 1, jpkm1 
    421             ze3t(:,:,jk) = te3t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:) 
    422          END DO 
    423          z_tmax = MAXVAL( ze3t(:,:,:) ) 
    424          IF( lk_mpp )   CALL mpp_max( z_tmax )                 ! max over the global domain 
    425          z_tmin = MINVAL( ze3t(:,:,:) ) 
    426          IF( lk_mpp )   CALL mpp_min( z_tmin )                 ! min over the global domain 
    427          ! - ML - test: for the moment, stop simulation for too large e3_t variations 
    428          IF( ( z_tmax >  rn_zdef_max ) .OR. ( z_tmin < - rn_zdef_max ) ) THEN 
    429             IF( lk_mpp ) THEN 
    430                CALL mpp_maxloc( ze3t, tmask, z_tmax, ijk_max(1), ijk_max(2), ijk_max(3) ) 
    431                CALL mpp_minloc( ze3t, tmask, z_tmin, ijk_min(1), ijk_min(2), ijk_min(3) ) 
    432             ELSE 
    433                ijk_max = MAXLOC( ze3t(:,:,:) ) 
    434                ijk_max(1) = ijk_max(1) + nimpp - 1 
    435                ijk_max(2) = ijk_max(2) + njmpp - 1 
    436                ijk_min = MINLOC( ze3t(:,:,:) ) 
    437                ijk_min(1) = ijk_min(1) + nimpp - 1 
    438                ijk_min(2) = ijk_min(2) + njmpp - 1 
    439             ENDIF 
    440             IF (lwp) THEN 
    441                WRITE(numout, *) 'MAX( te3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmax 
    442                WRITE(numout, *) 'at i, j, k=', ijk_max 
    443                WRITE(numout, *) 'MIN( te3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmin 
    444                WRITE(numout, *) 'at i, j, k=', ijk_min             
    445                CALL ctl_warn('MAX( ABS( te3t_a(:,:,:) ) / e3t_0(:,:,:) ) too high') 
    446             ENDIF 
    447          ENDIF 
    448          ! - ML - end test 
    449          ! - ML - Imposing these limits will cause a baroclinicity error which is corrected for below 
    450          te3t_a(:,:,:) = MIN( te3t_a(:,:,:),   rn_zdef_max * e3t_0(:,:,:) ) 
    451          te3t_a(:,:,:) = MAX( te3t_a(:,:,:), - rn_zdef_max * e3t_0(:,:,:) ) 
    452  
    453          ! 
    454          ! "tilda" change in the after scale factor 
    455          ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
    456          DO jk = 1, jpkm1 
    457             dte3t_a(:,:,jk) = te3t_a(:,:,jk) - te3t_b(:,:,jk) 
    458          END DO 
    459          ! III - Barotropic repartition of the sea surface height over the baroclinic profile 
    460          ! ================================================================================== 
    461          ! add ( ssh increment + "baroclinicity error" ) proportionly to e3t(n) 
    462          ! - ML - baroclinicity error should be better treated in the future 
    463          !        i.e. locally and not spread over the water column. 
    464          !        (keep in mind that the idea is to reduce Eulerian velocity as much as possible) 
    465          zht(:,:) = 0._wp 
    466          DO jk = 1, jpkm1 
    467             zht(:,:)  = zht(:,:) + te3t_a(:,:,jk) * tmask(:,:,jk) 
    468          END DO 
    469          z_scale(:,:) =  - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) ) 
    470          DO jk = 1, jpkm1 
    471             dte3t_a(:,:,jk) = dte3t_a(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk) 
    472          END DO 
    473  
    474       ENDIF 
    475  
    476       IF( ln_vvl_ztilde .OR. ln_vvl_layer )  THEN   ! z_tilde or layer coordinate ! 
    477       !                                           ! ---baroclinic part--------- ! 
    478          DO jk = 1, jpkm1 
    479             e3t_a(:,:,jk) = e3t_a(:,:,jk) + dte3t_a(:,:,jk) * tmask(:,:,jk) 
    480          END DO 
    481       ENDIF 
    482  
    483       IF( ln_vvl_dbg .AND. .NOT. ll_do_bclinic ) THEN   ! - ML - test: control prints for debuging 
    484          ! 
    485          IF( lwp ) WRITE(numout, *) 'kt =', kt 
    486          IF ( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 
    487             z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( zht(:,:) ) ) 
    488             IF( lk_mpp ) CALL mpp_max( z_tmax )                             ! max over the global domain 
    489             IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(SUM(te3t_a))) =', z_tmax 
    490          END IF 
    491          ! 
    492          zht(:,:) = 0.0_wp 
    493          DO jk = 1, jpkm1 
    494             zht(:,:) = zht(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 
    495          END DO 
    496          z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshn(:,:) - zht(:,:) ) ) 
    497          IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain 
    498          IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(e3t_n))) =', z_tmax 
    499          ! 
    500          zht(:,:) = 0.0_wp 
    501          DO jk = 1, jpkm1 
    502             zht(:,:) = zht(:,:) + e3t_a(:,:,jk) * tmask(:,:,jk) 
    503          END DO 
    504          z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssha(:,:) - zht(:,:) ) ) 
    505          IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain 
    506          IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(e3t_a))) =', z_tmax 
    507          ! 
    508          zht(:,:) = 0.0_wp 
    509          DO jk = 1, jpkm1 
    510             zht(:,:) = zht(:,:) + e3t_b(:,:,jk) * tmask(:,:,jk) 
    511          END DO 
    512          z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshb(:,:) - zht(:,:) ) ) 
    513          IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain 
    514          IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(e3t_b))) =', z_tmax 
    515          ! 
    516          z_tmax = MAXVAL( tmask(:,:,1) *  ABS( sshb(:,:) ) ) 
    517          IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain 
    518          IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(sshb))) =', z_tmax 
    519          ! 
    520          z_tmax = MAXVAL( tmask(:,:,1) *  ABS( sshn(:,:) ) ) 
    521          IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain 
    522          IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(sshn))) =', z_tmax 
    523          ! 
    524          z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssha(:,:) ) ) 
    525          IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain 
    526          IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssha))) =', z_tmax 
    527       END IF 
    528  
    529       ! *********************************** ! 
    530       ! After scale factors at u- v- points ! 
    531       ! *********************************** ! 
    532  
    533       CALL dom_vvl_interpol( e3t_a(:,:,:), e3u_a(:,:,:), 'U' ) 
    534       CALL dom_vvl_interpol( e3t_a(:,:,:), e3v_a(:,:,:), 'V' ) 
    535  
    536       ! *********************************** ! 
    537       ! After depths at u- v points         ! 
    538       ! *********************************** ! 
    539  
    540       hu_a(:,:) = e3u_a(:,:,1) * umask(:,:,1) 
    541       hv_a(:,:) = e3v_a(:,:,1) * vmask(:,:,1) 
    542       DO jk = 2, jpkm1 
    543          hu_a(:,:) = hu_a(:,:) + e3u_a(:,:,jk) * umask(:,:,jk) 
    544          hv_a(:,:) = hv_a(:,:) + e3v_a(:,:,jk) * vmask(:,:,jk) 
    545       END DO 
    546       !                                        ! Inverse of the local depth 
    547 !!gm BUG ?  don't understand the use of umask_i here ..... 
     226      !                             ! ------------------! 
     227      !                             ! z_star coordinate !     (and barotropic z-tilde part) 
     228      !                             ! ------------------! 
     229      ! 
     230      !                                   !==  after ssh  ==!  (u- and v-points) 
     231      DO jj = 2, jpjm1   ;   DO ji = 2, jpim1 
     232         zsshu_h(ji,jj) = 0.5_wp * ( ssha(ji,jj) + ssha(ji+1,jj) ) * ssumask(ji,jj) 
     233         zsshv_h(ji,jj) = 0.5_wp * ( ssha(ji,jj) + ssha(ji,jj+1) ) * ssvmask(ji,jj) 
     234      END DO             ;   END DO       
     235      CALL lbc_lnk_multi( zsshu_h(:,:), 'U', 1._wp , zsshv_h(:,:), 'V', 1._wp ) 
     236      ! 
     237      !                                   !==  after depths and its inverse  ==!  
     238         hu_a(:,:) = hu_0(:,:) + zsshu_h(:,:) 
     239         hv_a(:,:) = hv_0(:,:) + zsshv_h(:,:) 
    548240      r1_hu_a(:,:) = ssumask(:,:) / ( hu_a(:,:) + 1._wp - ssumask(:,:) ) 
    549241      r1_hv_a(:,:) = ssvmask(:,:) / ( hv_a(:,:) + 1._wp - ssvmask(:,:) ) 
     242      ! 
     243      !                                   !==  after scale factors  ==!  (e3t , e3u , e3v) 
     244      zssht_h(:,:) = ssha   (:,:) * r1_ht_0(:,:)           ! t-point 
     245      zsshu_h(:,:) = zsshu_h(:,:) * r1_hu_0(:,:)           ! u-point 
     246      zsshv_h(:,:) = zsshv_h(:,:) * r1_hv_0(:,:)           ! v-point 
     247      DO jk = 1, jpkm1 
     248         e3t_a(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * tmask(:,:,jk) ) 
     249         e3u_a(:,:,jk) = e3u_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) * umask(:,:,jk) ) 
     250         e3v_a(:,:,jk) = e3v_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * vmask(:,:,jk) ) 
     251      END DO 
    550252      ! 
    551253      IF( ln_timing )   CALL timing_stop('dom_vvl_sf_nxt') 
    … …  
    581283      ! 
    582284      INTEGER  ::   ji, jj, jk    ! dummy loop indices 
    583       REAL(wp) ::   zcoef, ze3f   ! local scalar 
     285      REAL(wp), DIMENSION(jpi,jpj) ::   zssht_h, zsshu_h, zsshv_h, zsshf_h 
    584286      !!---------------------------------------------------------------------- 
    585287      ! 
    … …  
    594296      ENDIF 
    595297      ! 
    596       ! Time filter and swap of scale factors 
    597       ! ===================================== 
    598       ! - ML - e3(t/u/v)_b are allready computed in dynnxt. 
    599       IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 
    600          IF( l_1st_euler ) THEN 
    601             te3t_n(:,:,:) = te3t_a(:,:,:) 
    602          ELSE 
    603             DO jk = 1, jpk 
    604                DO jj = 1, jpj 
    605                   DO ji = 1, jpi 
    606                      ze3f = te3t_n(ji,jj,jk)   & 
    607                         & + rn_atfp * ( te3t_b(ji,jj,jk) - 2.0_wp * te3t_n(ji,jj,jk) + te3t_a(ji,jj,jk) ) 
    608                      te3t_b(ji,jj,jk) = ze3f 
    609                      te3t_n(ji,jj,jk) = te3t_a(ji,jj,jk) 
    610                   END DO 
    611                END DO 
    612             END DO 
    613          ENDIF 
    614       ENDIF 
    615       gdept_b(:,:,:) = gdept_n(:,:,:) 
    616       gdepw_b(:,:,:) = gdepw_n(:,:,:) 
    617  
    618       e3t_n(:,:,:) = e3t_a(:,:,:) 
    619       e3u_n(:,:,:) = e3u_a(:,:,:) 
    620       e3v_n(:,:,:) = e3v_a(:,:,:) 
    621  
    622       ! Compute all missing vertical scale factor and depths 
    623       ! ==================================================== 
    624       ! Horizontal scale factor interpolations 
    625       ! -------------------------------------- 
     298      ! Swap and Compute all missing vertical scale factor and depths 
     299      ! ============================================================= 
    626300      ! - ML - e3u_b and e3v_b are allready computed in dynnxt 
    627301      ! - JC - hu_b, hv_b, hur_b, hvr_b also 
    628        
    629       CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F'  ) 
    630        
    631       ! Vertical scale factor interpolations 
    632       CALL dom_vvl_interpol( e3t_n(:,:,:),  e3w_n(:,:,:), 'W'  ) 
    633       CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' ) 
    634       CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' ) 
    635       CALL dom_vvl_interpol( e3t_b(:,:,:),  e3w_b(:,:,:), 'W'  ) 
    636       CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' ) 
    637       CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' ) 
    638  
    639       ! t- and w- points depth (set the isf depth as it is in the initial step) 
    640       gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1) 
    641       gdepw_n(:,:,1) = 0.0_wp 
    642       gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:) 
    643       DO jk = 2, jpk 
    644          DO jj = 1,jpj 
    645             DO ji = 1,jpi 
    646               !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 
    647                                                                  ! 1 for jk = mikt 
    648                zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 
    649                gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1) 
    650                gdept_n(ji,jj,jk) =    zcoef  * ( gdepw_n(ji,jj,jk  ) + 0.5 * e3w_n(ji,jj,jk) )  & 
    651                    &             + (1-zcoef) * ( gdept_n(ji,jj,jk-1) +       e3w_n(ji,jj,jk) )  
    652                gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj) 
    653             END DO 
    654          END DO 
    655       END DO 
    656  
     302      ! 
     303      ! - GM - to be updated :   e3f_n,  e3w_n , e3uw_n , e3vw_n  
     304      !                                  e3w_b , e3uw_b , e3vw_b 
     305      !                          gdept_n , gdepw_n , gde3w_n 
     306      !                          ht_n 
     307      !        to be swap    :   hu_a , hv_a  , r1_hu_a , r1_hv_a 
     308      ! 
    657309      ! Local depth and Inverse of the local depth of the water 
    658310      ! ------------------------------------------------------- 
     311      !                       ! swap of depth and scale factors 
     312      !                       ! =============================== 
     313      DO jk = 1, jpkm1                          ! swap n--> a   
     314         gdept_b(:,:,jk) = gdept_n(:,:,jk)         ! depth at t and w 
     315         gdepw_b(:,:,jk) = gdepw_n(:,:,jk) 
     316         e3t_n  (:,:,jk) = e3t_a  (:,:,jk)         ! e3t, e3u, e3v 
     317         e3u_n  (:,:,jk) = e3u_a  (:,:,jk) 
     318         e3v_n  (:,:,jk) = e3v_a  (:,:,jk) 
     319      END DO 
     320      ht_n(:,:) = ht_0(:,:) + sshn(:,:)            ! ocean thickness 
     321      ! 
    659322      hu_n(:,:) = hu_a(:,:)   ;   r1_hu_n(:,:) = r1_hu_a(:,:) 
    660323      hv_n(:,:) = hv_a(:,:)   ;   r1_hv_n(:,:) = r1_hv_a(:,:) 
    661324      ! 
    662       ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1) 
    663       DO jk = 2, jpkm1 
    664          ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 
     325      !                    !==  before :  
     326      !                                            !* ssh at u- and v-points) 
     327      DO jj = 2, jpjm1   ;   DO ji = 2, jpim1 
     328         zsshu_h(ji,jj) = 0.5_wp  * ( sshb(ji  ,jj) + sshb(ji+1,jj  ) ) * ssumask(ji,jj) 
     329         zsshv_h(ji,jj) = 0.5_wp  * ( sshb(ji  ,jj) + sshb(ji  ,jj+1) ) * ssvmask(ji,jj) 
     330      END DO             ;   END DO       
     331      CALL lbc_lnk_multi( zsshu_h(:,:),'U', 1._wp , zsshv_h(:,:),'V', 1._wp ) 
     332      ! 
     333      !                                            !*  e3w_b , e3uw_b , e3vw_b 
     334      zssht_h(:,:) = sshb (:,:) * r1_ht_0(:,:)           ! w-point 
     335      zsshu_h(:,:) = zsshu_h(:,:) * r1_hu_0(:,:)           ! uw-point 
     336      zsshv_h(:,:) = zsshv_h(:,:) * r1_hv_0(:,:)           ! vw-point 
     337      DO jk = 1, jpkm1 
     338          e3w_b(:,:,jk) =  e3w_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * tmask(:,:,jk) ) 
     339         e3uw_b(:,:,jk) = e3uw_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) * umask(:,:,jk) ) 
     340         e3vw_b(:,:,jk) = e3vw_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * vmask(:,:,jk) ) 
    665341      END DO 
    666  
     342      !       
     343      !                    !==  now    :  
     344      !                                            !* ssh at u- and v-points) 
     345      DO jj = 1, jpjm1   ;   DO ji = 1, jpim1            ! start from 1 for f-point 
     346         zsshu_h(ji,jj) = 0.50_wp * ( sshn(ji  ,jj) + sshn(ji+1,jj  ) ) * ssumask(ji,jj) 
     347         zsshv_h(ji,jj) = 0.50_wp * ( sshn(ji  ,jj) + sshn(ji  ,jj+1) ) * ssvmask(ji,jj) 
     348         zsshf_h(ji,jj) = 0.25_wp * ( sshn(ji  ,jj) + sshn(ji  ,jj+1)   &  
     349            &                       + sshn(ji+1,jj) + sshn(ji+1,jj+1) ) * ssfmask(ji,jj) 
     350      END DO             ;   END DO       
     351      CALL lbc_lnk_multi( zsshu_h(:,:),'U', 1._wp , zsshv_h(:,:),'V', 1._wp , zsshf_h(:,:),'F', 1._wp )       
     352      ! 
     353      !                                            !* e3w_n , e3uw_n , e3vw_n, e3f_n  
     354      zssht_h(:,:) = sshn   (:,:) * r1_ht_0(:,:)           ! t- & w-point 
     355      zsshu_h(:,:) = zsshu_h(:,:) * r1_hu_0(:,:)           ! uw-point 
     356      zsshv_h(:,:) = zsshv_h(:,:) * r1_hv_0(:,:)           ! vw-point 
     357      zsshf_h(:,:) = zsshf_h(:,:) * r1_hf_0(:,:)           ! f-point 
     358      DO jk = 1, jpkm1 
     359          e3w_n(:,:,jk) =  e3w_0(:,:,jk) * ( 1._wp + zssht_h(:,:) *  tmask(:,:,jk) ) 
     360         e3uw_n(:,:,jk) = e3uw_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) * wumask(:,:,jk) ) 
     361         e3vw_n(:,:,jk) = e3vw_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * wvmask(:,:,jk) ) 
     362          e3f_n(:,:,jk) =  e3f_0(:,:,jk) * ( 1._wp + zsshf_h(:,:) *  fmask(:,:,jk) ) 
     363      END DO       
     364      !  
     365      zssht_h(:,:) = 1._wp + sshn (:,:) * r1_ht_0(:,:)   ! t-point 
     366      ! 
     367      IF( ln_isfcav ) THEN    ! ISF cavities : ssh scaling not applied over the iceshelf thickness  
     368         DO jk = 1, jpkm1 
     369            gdept_n(:,:,jk) = ( gdept_0(:,:,jk) - risfdep(:,:) ) * zssht_h(:,:) + risfdep(:,:) 
     370            gdepw_n(:,:,jk) = ( gdepw_0(:,:,jk) - risfdep(:,:) ) * zssht_h(:,:) + risfdep(:,:) 
     371            gde3w_n(:,:,jk) =   gdept_n(:,:,jk) - sshn   (:,:) 
     372         END DO 
     373      ELSE                    ! no ISF cavities  
     374         DO jk = 1, jpkm1 
     375            gdept_n(:,:,jk) = gdept_0(:,:,jk) * zssht_h(:,:) 
     376            gdepw_n(:,:,jk) = gdepw_0(:,:,jk) * zssht_h(:,:) 
     377            gde3w_n(:,:,jk) = gdept_n(:,:,jk) - sshn(:,:) 
     378         END DO 
     379      ENDIF 
     380      ! 
    667381      ! write restart file 
    668382      ! ================== 
    … …  
    672386      ! 
    673387   END SUBROUTINE dom_vvl_sf_swp 
    674  
    675  
    676    SUBROUTINE dom_vvl_interpol( pe3_in, pe3_out, pout ) 
    677       !!--------------------------------------------------------------------- 
    678       !!                  ***  ROUTINE dom_vvl__interpol  *** 
    679       !! 
    680       !! ** Purpose :   interpolate scale factors from one grid point to another 
    681       !! 
    682       !! ** Method  :   e3_out = e3_0 + interpolation(e3_in - e3_0) 
    683       !!                - horizontal interpolation: grid cell surface averaging 
    684       !!                - vertical interpolation: simple averaging 
    685       !!---------------------------------------------------------------------- 
    686       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::  pe3_in    ! input e3 to be interpolated 
    687       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::  pe3_out   ! output interpolated e3 
    688       CHARACTER(LEN=*)                , INTENT(in   ) ::  pout      ! grid point of out scale factors 
    689       !                                                             !   =  'U', 'V', 'W, 'F', 'UW' or 'VW' 
    690       ! 
    691       INTEGER ::   ji, jj, jk                                       ! dummy loop indices 
    692       REAL(wp) ::  zlnwd                                            ! =1./0. when ln_wd_il = T/F 
    693       !!---------------------------------------------------------------------- 
    694       ! 
    695       IF(ln_wd_il) THEN 
    696         zlnwd = 1.0_wp 
    697       ELSE 
    698         zlnwd = 0.0_wp 
    699       END IF 
    700       ! 
    701       SELECT CASE ( pout )    !==  type of interpolation  ==! 
    702          ! 
    703       CASE( 'U' )                   !* from T- to U-point : hor. surface weighted mean 
    704          DO jk = 1, jpk 
    705             DO jj = 1, jpjm1 
    706                DO ji = 1, fs_jpim1   ! vector opt. 
    707                   pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj)   & 
    708                      &                       * (   e1e2t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     & 
    709                      &                           + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) ) 
    710                END DO 
    711             END DO 
    712          END DO 
    713          CALL lbc_lnk( pe3_out(:,:,:), 'U', 1._wp ) 
    714          pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:) 
    715          ! 
    716       CASE( 'V' )                   !* from T- to V-point : hor. surface weighted mean 
    717          DO jk = 1, jpk 
    718             DO jj = 1, jpjm1 
    719                DO ji = 1, fs_jpim1   ! vector opt. 
    720                   pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk)  * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj)   & 
    721                      &                       * (   e1e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     & 
    722                      &                           + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) ) 
    723                END DO 
    724             END DO 
    725          END DO 
    726          CALL lbc_lnk( pe3_out(:,:,:), 'V', 1._wp ) 
    727          pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:) 
    728          ! 
    729       CASE( 'F' )                   !* from U-point to F-point : hor. surface weighted mean 
    730          DO jk = 1, jpk 
    731             DO jj = 1, jpjm1 
    732                DO ji = 1, fs_jpim1   ! vector opt. 
    733                   pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) & 
    734                      &                       *    r1_e1e2f(ji,jj)                                                  & 
    735                      &                       * (   e1e2u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     & 
    736                      &                           + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) ) 
    737                END DO 
    738             END DO 
    739          END DO 
    740          CALL lbc_lnk( pe3_out(:,:,:), 'F', 1._wp ) 
    741          pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:) 
    742          ! 
    743       CASE( 'W' )                   !* from T- to W-point : vertical simple mean 
    744          ! 
    745          pe3_out(:,:,1) = e3w_0(:,:,1) + pe3_in(:,:,1) - e3t_0(:,:,1) 
    746          ! - ML - The use of mask in this formulea enables the special treatment of the last w-point without indirect adressing 
    747 !!gm BUG? use here wmask in case of ISF ?  to be checked 
    748          DO jk = 2, jpk 
    749             pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )   & 
    750                &                            * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) )                               & 
    751                &                            +            0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )     & 
    752                &                            * ( pe3_in(:,:,jk  ) - e3t_0(:,:,jk  ) ) 
    753          END DO 
    754          ! 
    755       CASE( 'UW' )                  !* from U- to UW-point : vertical simple mean 
    756          ! 
    757          pe3_out(:,:,1) = e3uw_0(:,:,1) + pe3_in(:,:,1) - e3u_0(:,:,1) 
    758          ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing 
    759 !!gm BUG? use here wumask in case of ISF ?  to be checked 
    760          DO jk = 2, jpk 
    761             pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )  & 
    762                &                             * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) )                              & 
    763                &                             +            0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )    & 
    764                &                             * ( pe3_in(:,:,jk  ) - e3u_0(:,:,jk  ) ) 
    765          END DO 
    766          ! 
    767       CASE( 'VW' )                  !* from V- to VW-point : vertical simple mean 
    768          ! 
    769          pe3_out(:,:,1) = e3vw_0(:,:,1) + pe3_in(:,:,1) - e3v_0(:,:,1) 
    770          ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing 
    771 !!gm BUG? use here wvmask in case of ISF ?  to be checked 
    772          DO jk = 2, jpk 
    773             pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )  & 
    774                &                             * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) )                              & 
    775                &                             +            0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )    & 
    776                &                             * ( pe3_in(:,:,jk  ) - e3v_0(:,:,jk  ) ) 
    777          END DO 
    778       END SELECT 
    779       ! 
    780    END SUBROUTINE dom_vvl_interpol 
    781388 
    782389 
    … …  
    804411         IF( ln_rstart ) THEN                   !* Read the restart file 
    805412            CALL rst_read_open                  !  open the restart file if necessary 
    806             CALL iom_get( numror, jpdom_autoglo, 'sshn'   , sshn, ldxios = lrxios    ) 
    807413            ! 
    808             id1 = iom_varid( numror, 'e3t_b'      , ldstop = .FALSE. ) 
    809             id2 = iom_varid( numror, 'e3t_n'      , ldstop = .FALSE. ) 
    810             id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. ) 
    811             id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. ) 
    812             id5 = iom_varid( numror, 'hdiv_lf'    , ldstop = .FALSE. ) 
    813             !                             ! --------- ! 
    814             !                             ! all cases ! 
    815             !                             ! --------- ! 
     414            id1 = iom_varid( numror, 'sshb'      , ldstop = .FALSE. ) 
     415            id2 = iom_varid( numror, 'sshn'      , ldstop = .FALSE. ) 
     416            ! 
    816417            IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist 
    817                CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:), ldxios = lrxios ) 
    818                CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:), ldxios = lrxios ) 
    819                ! needed to restart if land processor not computed  
    820                IF(lwp) write(numout,*) 'dom_vvl_rst : e3t_b and e3t_n found in restart files' 
    821                WHERE ( tmask(:,:,:) == 0.0_wp )  
    822                   e3t_n(:,:,:) = e3t_0(:,:,:) 
    823                   e3t_b(:,:,:) = e3t_0(:,:,:) 
    824                END WHERE 
     418               IF(lwp) write(numout,*) 'dom_vvl_rst : both sshb and sshn found in restart files' 
     419               ! 
     420!!gm  Question: use jpdom_data above to read data over jpi x jpj    (like is dom_hgr_read and dom_zgr_read) 
     421!!              so that it will work with land processor suppression 
     422!               CALL iom_get( numror, jpdom_autoglo, 'sshn'   , sshn, ldxios = lrxios    ) 
     423!               CALL iom_get( numror, jpdom_autoglo, 'sshb'   , sshb, ldxios = lrxios    ) 
     424!!gm  
     425               CALL iom_get( numror, jpdom_data, 'sshn'   , sshn, ldxios = lrxios    ) 
     426               CALL iom_get( numror, jpdom_data, 'sshb'   , sshb, ldxios = lrxios    ) 
     427!!gm end 
    825428               IF( l_1st_euler ) THEN 
    826                   e3t_b(:,:,:) = e3t_n(:,:,:) 
     429                  sshb(:,:) = sshn(:,:) 
    827430               ENDIF 
    828431            ELSE IF( id1 > 0 ) THEN 
    829                IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart files' 
    830                IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.' 
    831                IF(lwp) write(numout,*) 'l_1st_euler is forced to true' 
    832                CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:), ldxios = lrxios ) 
    833                e3t_n(:,:,:) = e3t_b(:,:,:) 
     432               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : sshn not found in restart files' 
     433               IF(lwp) write(numout,*) '   set sshn = sshb  and force l_1st_euler = true' 
     434!!gm               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:), ldxios = lrxios ) 
     435               CALL iom_get( numror, jpdom_data, 'sshb', sshb, ldxios = lrxios ) 
     436               sshn(:,:) = sshb(:,:) 
    834437               l_1st_euler = .TRUE. 
    835438            ELSE IF( id2 > 0 ) THEN 
    836                IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_b not found in restart files' 
    837                IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.' 
    838                IF(lwp) write(numout,*) 'l_1st_euler is forced to true' 
    839                CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:), ldxios = lrxios ) 
    840                e3t_b(:,:,:) = e3t_n(:,:,:) 
     439               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : sshb not found in restart files' 
     440               IF(lwp) write(numout,*) 'set sshb = sshn  and force l_1st_euler = true' 
     441               CALL iom_get( numror, jpdom_data, 'sshn', sshb, ldxios = lrxios ) 
     442               sshb(:,:) = sshn(:,:) 
    841443               l_1st_euler = .TRUE. 
    842444            ELSE 
    843                IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart file' 
    844                IF(lwp) write(numout,*) 'Compute scale factor from sshn' 
    845                IF(lwp) write(numout,*) 'l_1st_euler is forced to true' 
    846                DO jk = 1, jpk 
    847                   e3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) & 
    848                       &                          / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   & 
    849                       &          + e3t_0(:,:,jk)                               * (1._wp -tmask(:,:,jk)) 
    850                END DO 
    851                e3t_b(:,:,:) = e3t_n(:,:,:) 
     445               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : sshb and sshn not found in restart file' 
     446               IF(lwp) write(numout,*) 'set sshb = sshn = 0  and force l_1st_euler = true' 
     447               sshb(:,:) = 0._wp 
     448               sshn(:,:) = 0._wp 
    852449               l_1st_euler = .TRUE. 
    853450            ENDIF 
    854             !                             ! ----------- ! 
    855             IF( ln_vvl_zstar ) THEN       ! z_star case ! 
    856                !                          ! ----------- ! 
    857                IF( MIN( id3, id4 ) > 0 ) THEN 
    858                   CALL ctl_stop( 'dom_vvl_rst: z_star cannot restart from a z_tilde or layer run' ) 
    859                ENDIF 
    860                !                          ! ----------------------- ! 
    861             ELSE                          ! z_tilde and layer cases ! 
    862                !                          ! ----------------------- ! 
    863                IF( MIN( id3, id4 ) > 0 ) THEN  ! all required arrays exist 
    864                   CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_b', te3t_b(:,:,:), ldxios = lrxios ) 
    865                   CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_n', te3t_n(:,:,:), ldxios = lrxios ) 
    866                ELSE                            ! one at least array is missing 
    867                   te3t_b(:,:,:) = 0.0_wp 
    868                   te3t_n(:,:,:) = 0.0_wp 
    869                ENDIF 
    870                !                          ! ------------ ! 
    871                IF( ln_vvl_ztilde ) THEN   ! z_tilde case ! 
    872                   !                       ! ------------ ! 
    873                   IF( id5 > 0 ) THEN  ! required array exists 
    874                      CALL iom_get( numror, jpdom_autoglo, 'hdiv_lf', hdiv_lf(:,:,:), ldxios = lrxios ) 
    875                   ELSE                ! array is missing 
    876                      hdiv_lf(:,:,:) = 0.0_wp 
    877                   ENDIF 
    878                ENDIF 
    879             ENDIF 
    880             ! 
    881451         ELSE                                   !* Initialize at "rest" 
    882452            ! 
    883  
    884453            IF( ll_wd ) THEN   ! MJB ll_wd edits start here - these are essential  
    885454               ! 
    886                IF( cn_cfg == 'wad' ) THEN 
    887                   ! Wetting and drying test case 
     455               IF( cn_cfg == 'wad' ) THEN             ! Wetting and drying test case 
    888456                  CALL usr_def_istate( gdept_b, tmask, tsb, ub, vb, sshb  ) 
    889                   tsn  (:,:,:,:) = tsb (:,:,:,:)       ! set now values from to before ones 
     457                  tsn  (:,:,:,:) = tsb (:,:,:,:)            ! set now values from to before ones 
    890458                  sshn (:,:)     = sshb(:,:) 
    891459                  un   (:,:,:)   = ub  (:,:,:) 
    892460                  vn   (:,:,:)   = vb  (:,:,:) 
    893                ELSE 
    894                   ! if not test case 
     461               ELSE                                   ! Not the test case 
    895462                  sshn(:,:) = -ssh_ref 
    896463                  sshb(:,:) = -ssh_ref 
    897  
     464                  ! 
    898465                  DO jj = 1, jpj 
    899466                     DO ji = 1, jpi 
    900                         IF( ht_0(ji,jj)-ssh_ref <  rn_wdmin1 ) THEN ! if total depth is less than min depth 
    901  
     467                        IF( ht_0(ji,jj)-ssh_ref <  rn_wdmin1 ) THEN  ! if total depth is less than min depth 
    902468                           sshb(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) ) 
    903469                           sshn(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) ) 
    904470                           ssha(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) ) 
    905471                        ENDIF 
    906                      ENDDO 
    907                   ENDDO 
    908                ENDIF !If test case else 
    909  
    910                ! Adjust vertical metrics for all wad 
    911                DO jk = 1, jpk 
    912                   e3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:)  ) & 
    913                     &                            / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   & 
    914                     &            + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) ) 
    915                END DO 
    916                e3t_b(:,:,:) = e3t_n(:,:,:) 
    917  
     472                     END DO 
     473                  END DO 
     474               ENDIF     ! If test case else 
     475               ! 
    918476               DO ji = 1, jpi 
    919477                  DO jj = 1, jpj 
    920                      IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN 
     478                     IF ( ht_0(ji,jj) /= 0._wp .AND. NINT( ssmask(ji,jj) ) == 1 ) THEN 
    921479                       CALL ctl_stop( 'dom_vvl_rst: ht_0 must be positive at potentially wet points' ) 
    922480                     ENDIF 
    … …  
    926484            ELSE 
    927485               ! 
    928                ! Just to read set ssh in fact, called latter once vertical grid 
    929                ! is set up: 
     486               ! Just to read set ssh in fact, called latter once vertical grid is set up: 
    930487!               CALL usr_def_istate( gdept_0, tmask, tsb, ub, vb, sshb  ) 
    931 !               ! 
    932 !               DO jk=1,jpk 
    933 !                  e3t_b(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshb(:,:) ) & 
    934 !                     &            / ( ht_0(:,:) + 1._wp -ssmask(:,:) ) * tmask(:,:,jk) 
    935 !               END DO 
    936 !               e3t_n(:,:,:) = e3t_b(:,:,:) 
    937                 sshn(:,:)=0._wp 
    938                 e3t_n(:,:,:)=e3t_0(:,:,:) 
    939                 e3t_b(:,:,:)=e3t_0(:,:,:) 
     488               sshn(:,:) = 0._wp 
     489               sshb(:,:) = 0._wp 
    940490               ! 
    941             END IF           ! end of ll_wd edits 
    942  
    943             IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN 
    944                te3t_b(:,:,:) = 0._wp 
    945                te3t_n(:,:,:) = 0._wp 
    946                IF( ln_vvl_ztilde ) hdiv_lf(:,:,:) = 0._wp 
    947491            END IF 
     492            ! 
    948493         ENDIF 
    949494         ! 
    950495      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file 
    951496         !                                   ! =================== 
    952          IF(lwp) WRITE(numout,*) '---- dom_vvl_rst ----' 
    953          IF( lwxios ) CALL iom_swap(      cwxios_context          ) 
    954          !                                           ! --------- ! 
    955          !                                           ! all cases ! 
    956          !                                           ! --------- ! 
    957          CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t_b(:,:,:), ldxios = lwxios ) 
    958          CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t_n(:,:,:), ldxios = lwxios ) 
    959          !                                           ! ----------------------- ! 
    960          IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN  ! z_tilde and layer cases ! 
    961             !                                        ! ----------------------- ! 
    962             CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_b', te3t_b(:,:,:), ldxios = lwxios) 
    963             CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_n', te3t_n(:,:,:), ldxios = lwxios) 
    964          END IF 
    965          !                                           ! -------------!     
    966          IF( ln_vvl_ztilde ) THEN                    ! z_tilde case ! 
    967             !                                        ! ------------ ! 
    968             CALL iom_rstput( kt, nitrst, numrow, 'hdiv_lf', hdiv_lf(:,:,:), ldxios = lwxios) 
    969          ENDIF 
    970          ! 
    971          IF( lwxios ) CALL iom_swap(      cxios_context          ) 
     497 
     498!!gm      DO NOTHING,   sshb and sshn  are written in restart.F90 
     499 
    972500      ENDIF 
    973501      ! 
    … …  
    1024552      ENDIF 
    1025553      ! 
     554 
     555!!gm  
     556      IF ( ln_vvl_ztilde .OR. ln_vvl_ztilde_as_zstar )   CALL ctl_stop( 'z-tilde NOT available in this branch' ) 
     557!!gm 
     558 
     559      ! 
    1026560      ioptio = 0                      ! Parameter control 
    1027561      IF( ln_vvl_ztilde_as_zstar )   ln_vvl_ztilde = .true. 
    … …  
    1047581   END SUBROUTINE dom_vvl_ctl 
    1048582 
     583    
     584   SUBROUTINE ssh2e3_now 
     585      !!---------------------------------------------------------------------- 
     586      !!                  ***  ROUTINE ssh2e3_now  *** 
     587      !!---------------------------------------------------------------------- 
     588      INTEGER ::   ji, jj, jk 
     589      REAL(wp), DIMENSION(jpi,jpj) ::   zssht_h, zsshu_h, zsshv_h, zsshf_h 
     590      !!---------------------------------------------------------------------- 
     591      ! 
     592      !                             !==  ssh at u- and v-points  ==! 
     593      ! 
     594      DO jj = 1, jpjm1                    ! start from 1 due to f-point 
     595         DO ji = 1, jpim1 
     596            zsshu_h(ji,jj) = 0.50_wp * ( sshn(ji  ,jj) + sshn(ji+1,jj  ) ) * ssumask(ji,jj) 
     597            zsshv_h(ji,jj) = 0.50_wp * ( sshn(ji  ,jj) + sshn(ji  ,jj+1) ) * ssvmask(ji,jj) 
     598            zsshf_h(ji,jj) = 0.25_wp * ( sshn(ji  ,jj) + sshn(ji  ,jj+1)   &  
     599               &                       + sshn(ji+1,jj) + sshn(ji+1,jj+1) ) * ssfmask(ji,jj) 
     600         END DO 
     601      END DO       
     602      CALL lbc_lnk_multi( zsshu_h(:,:),'U', 1._wp , zsshv_h(:,:),'V', 1._wp , zsshf_h(:,:),'F', 1._wp ) 
     603      ! 
     604      !                             !==  ht, hu and hv  == !   (and their inverse) 
     605      ! 
     606      ht_n   (:,:) = ht_0(:,:) +  sshn  (:,:) 
     607      hu_n   (:,:) = hu_0(:,:) + zsshu_h(:,:) 
     608      hv_n   (:,:) = hv_0(:,:) + zsshv_h(:,:) 
     609      r1_hu_n(:,:) = ssumask(:,:) / ( hu_n(:,:) + 1._wp - ssumask(:,:) )    ! ss mask mask due to ISF 
     610      r1_hv_n(:,:) = ssvmask(:,:) / ( hv_n(:,:) + 1._wp - ssvmask(:,:) ) 
     611      !       
     612      !                             !==  ssh / h  factor at t-, u- ,v- & f-points  ==! 
     613      ! 
     614      zssht_h(:,:) =  sshn  (:,:) * r1_ht_0(:,:) 
     615      zsshu_h(:,:) = zsshu_h(:,:) * r1_hu_0(:,:) 
     616      zsshv_h(:,:) = zsshv_h(:,:) * r1_hv_0(:,:) 
     617      zsshf_h(:,:) = zsshf_h(:,:) * r1_hf_0(:,:) 
     618      ! 
     619      !                             !==  e3t, e3w  ,  e3u, e3uw ,  e3v, e3vw  , and e3f  ==! 
     620      !       
     621      DO jk = 1, jpkm1 
     622          e3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * tmask(:,:,jk) ) 
     623          e3w_n(:,:,jk) =  e3w_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * tmask(:,:,jk) ) 
     624          ! 
     625          e3u_n(:,:,jk) =  e3u_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) *  umask(:,:,jk) ) 
     626         e3uw_n(:,:,jk) = e3uw_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) * wumask(:,:,jk) ) 
     627         ! 
     628          e3v_n(:,:,jk) =  e3v_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) *  vmask(:,:,jk) ) 
     629         e3vw_n(:,:,jk) = e3vw_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * wvmask(:,:,jk) ) 
     630          ! 
     631          e3f_n(:,:,jk) =  e3f_0(:,:,jk) * ( 1._wp + zsshf_h(:,:) * fmask(:,:,jk) ) 
     632      END DO 
     633      !       
     634      !                             !== depth of t- and w-points  ==! 
     635      ! 
     636      zssht_h(:,:) = 1._wp + zssht_h(:,:)     ! = 1 + sshn / ht_0 
     637      ! 
     638      IF( ln_isfcav ) THEN    ! ISF cavities : ssh scaling not applied over the iceshelf thickness  
     639         DO jk = 1, jpkm1 
     640            gdept_n(:,:,jk) = ( gdept_0(:,:,jk) - risfdep(:,:) ) * zssht_h(:,:) + risfdep(:,:) 
     641            gdepw_n(:,:,jk) = ( gdepw_0(:,:,jk) - risfdep(:,:) ) * zssht_h(:,:) + risfdep(:,:) 
     642            gde3w_n(:,:,jk) =   gdept_n(:,:,jk) - sshn(:,:) 
     643         END DO 
     644      ELSE                    ! no ISF cavities 
     645         DO jk = 1, jpkm1 
     646            gdept_n(:,:,jk) = gdept_0(:,:,jk) * zssht_h(:,:) 
     647            gdepw_n(:,:,jk) = gdepw_0(:,:,jk) * zssht_h(:,:) 
     648            gde3w_n(:,:,jk) = gdept_n(:,:,jk) - sshn(:,:) 
     649         END DO 
     650      ENDIF 
     651      ! 
     652   END SUBROUTINE ssh2e3_now 
     653    
     654    
     655   SUBROUTINE ssh2e3_before 
     656      !!---------------------------------------------------------------------- 
     657      !!                  ***  ROUTINE ssh2e3_before  *** 
     658      !!---------------------------------------------------------------------- 
     659      INTEGER ::   ji, jj, jk 
     660      REAL(wp), DIMENSION(jpi,jpj) ::   zssht_h, zsshu_h, zsshv_h 
     661      !!---------------------------------------------------------------------- 
     662      ! 
     663      !                             !==  ssh at u- and v-points  ==! 
     664      DO jj = 2, jpjm1 
     665         DO ji = 2, jpim1 
     666            zsshu_h(ji,jj) = 0.5_wp  * ( sshb(ji  ,jj) + sshb(ji+1,jj  ) ) * ssumask(ji,jj) 
     667            zsshv_h(ji,jj) = 0.5_wp  * ( sshb(ji  ,jj) + sshb(ji  ,jj+1) ) * ssvmask(ji,jj) 
     668         END DO 
     669      END DO       
     670      CALL lbc_lnk_multi( zsshu_h(:,:),'U', 1._wp , zsshv_h(:,:),'V', 1._wp ) 
     671      ! 
     672      !                             !==  ht, hu and hv  == !   (and their inverse) 
     673         hu_b(:,:) = hu_0(:,:) + zsshu_h(:,:) 
     674         hv_b(:,:) = hv_0(:,:) + zsshv_h(:,:) 
     675      r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) ) 
     676      r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) ) 
     677      ! 
     678      !       
     679      !                             !==  ssh / h  factor at t-, u- ,v- & f-points  ==! 
     680      zssht_h(:,:) = sshb (:,:) * r1_ht_0(:,:) 
     681      zsshu_h   (:,:) = zsshu_h(:,:) * r1_hu_0(:,:) 
     682      zsshv_h   (:,:) = zsshv_h(:,:) * r1_hv_0(:,:) 
     683      ! 
     684      !                             !==  e3t, e3w  ,  e3u, e3uw , and  e3v, e3vw  ==! 
     685      DO jk = 1, jpkm1 
     686          e3t_b(:,:,jk) =  e3t_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * tmask(:,:,jk) ) 
     687          e3w_b(:,:,jk) =  e3w_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * tmask(:,:,jk) ) 
     688          ! 
     689          e3u_b(:,:,jk) =  e3u_0(:,:,jk) * ( 1._wp + zsshu_h  (:,:) *  umask(:,:,jk) ) 
     690         e3uw_b(:,:,jk) = e3uw_0(:,:,jk) * ( 1._wp + zsshu_h  (:,:) * wumask(:,:,jk) ) 
     691          ! 
     692          e3v_b(:,:,jk) =  e3v_0(:,:,jk) * ( 1._wp + zsshv_h  (:,:) *  vmask(:,:,jk) ) 
     693         e3vw_b(:,:,jk) = e3vw_0(:,:,jk) * ( 1._wp + zsshv_h  (:,:) * wvmask(:,:,jk) ) 
     694      END DO 
     695      !    
     696   END SUBROUTINE ssh2e3_before 
     697    
    1049698   !!====================================================================== 
    1050699END MODULE domvvl 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/iscplrst.F90

    r9939 r10001  
    1414   USE dom_oce         ! ocean space and time domain 
    1515   USE domwri          ! ocean space and time domain 
    16    USE domvvl   , ONLY : dom_vvl_interpol 
    1716   USE phycst          ! physical constants 
    1817   USE sbc_oce         ! surface boundary condition variables 
    … …  
    137136      REAL(wp), DIMENSION(jpi,jpj)          :: zdmask , zsmask0, zucorr, zbub, zbun, zssh0, zhu1, zde3t 
    138137      REAL(wp), DIMENSION(jpi,jpj)          :: zdsmask, zsmask1, zvcorr, zbvb, zbvn, zssh1, zhv1 
    139       REAL(wp), DIMENSION(jpi,jpj,jpk)      :: ztmask0, zwmaskn, ztrp 
     138      REAL(wp), DIMENSION(jpi,jpj,jpk)      :: ztmask0, zwmaskn 
    140139      REAL(wp), DIMENSION(jpi,jpj,jpk)      :: ztmask1, zwmaskb, ztmp3d 
    141140      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts) :: zts0 
     141      REAL(wp), DIMENSION(jpi,jpj) ::   z_ssh_h0, zsshu, zsshv, zsshf 
    142142      !!---------------------------------------------------------------------- 
    143143      ! 
    … …  
    163163         DO jj = 2,jpj-1 
    164164            DO ji = fs_2, fs_jpim1   ! vector opt. 
    165                jip1=ji+1; jim1=ji-1; 
    166                jjp1=jj+1; jjm1=jj-1; 
    167                summsk=(zsmask0(jip1,jj)+zsmask0(jim1,jj)+zsmask0(ji,jjp1)+zsmask0(ji,jjm1)) 
     165               summsk = zsmask0(ji+1,jj)+zsmask0(ji-1,jj)+zsmask0(ji,jj+1)+zsmask0(ji,jj-1) 
    168166               IF (zdsmask(ji,jj) == 1._wp .AND. summsk /= 0._wp) THEN 
    169                   sshn(ji,jj)=( zssh0(jip1,jj)*zsmask0(jip1,jj)     & 
    170                   &           + zssh0(jim1,jj)*zsmask0(jim1,jj)     & 
    171                   &           + zssh0(ji,jjp1)*zsmask0(ji,jjp1)     & 
    172                   &           + zssh0(ji,jjm1)*zsmask0(ji,jjm1))/summsk 
    173                   zsmask1(ji,jj)=1._wp 
     167                  sshn(ji,jj)=( zssh0(ji+1,jj)*zsmask0(ji+1,jj)     & 
     168                  &           + zssh0(ji-1,jj)*zsmask0(ji-1,jj)     & 
     169                  &           + zssh0(ji,jj+1)*zsmask0(ji,jj+1)     & 
     170                  &           + zssh0(ji,jj-1)*zsmask0(ji,jj+1)) / summsk 
     171                  zsmask1(ji,jj) = 1._wp 
    174172               ENDIF 
    175173            END DO 
    … …  
    185183      IF ( .NOT.ln_linssh ) THEN 
    186184      ! Reconstruction of all vertical scale factors at now time steps 
    187       ! ============================================================================= 
     185      ! ====================================================================== 
     186       
     187!!gm Question : bug ???? 
     188      ! at this stage  is ht_0, r1_ht0 already computed ???? 
     189      ! I have the feeling that this should be done in a different manner... 
     190      ! that is iscpl_rst_interpol  should be call directly in restart !!! 
     191       
     192      ! in my changes here I assume that ht_0 AND  r1_ht_0  have been updated  
     193      ! Note that the former calculation were using ht_0  so if it as not been updated ===>>> BUG 
     194!!gm 
     195       
     196       
    188197      ! Horizontal scale factor interpolations 
    189198      ! -------------------------------------- 
    190          DO jk = 1,jpk 
    191             DO jj=1,jpj 
    192                DO ji=1,jpi 
    193                   IF (tmask(ji,jj,1) == 0._wp .OR. ptmask_b(ji,jj,1) == 0._wp) THEN 
    194                      e3t_n(ji,jj,jk) = e3t_0(ji,jj,jk) * ( 1._wp + sshn(ji,jj) /   & 
    195                      &   ( ht_0(ji,jj) + 1._wp - ssmask(ji,jj) ) * tmask(ji,jj,jk) ) 
    196                   ENDIF 
    197                END DO 
    198             END DO 
    199          END DO 
    200          ! 
    201          CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' ) 
    202          CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' ) 
    203          CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' ) 
    204  
    205          ! Vertical scale factor interpolations 
    206          ! ------------------------------------ 
    207          CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W'  ) 
    208          CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' ) 
    209          CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' ) 
     199         DO jj = 1, jpj 
     200            DO ji = 1, jpi 
     201               IF ( tmask(ji,jj,1) == 0._wp .OR. ptmask_b(ji,jj,1) == 0._wp ) THEN 
     202                  DO jk = 1, jpk 
     203                     e3t_n(ji,jj,jk) = e3t_0(ji,jj,jk) * ( 1._wp + sshn(ji,jj) * r1_ht_0(ji,jj) * tmask(ji,jj,jk) ) 
     204                  END DO 
     205               ENDIF 
     206            END DO 
     207         END DO 
     208         ! 
     209!!gm  Note that if this routine is called in dom_vvl_init then all the lines below are uselss !!! 
     210!!        they are a duplication of dom_vvl_init lines 
     211 
     212         !                                   !==  now fields  ==! 
     213         ! 
     214         !                                            !* ssh at u- and v-points) 
     215         DO jj = 1, jpjm1   ;   DO ji = 1, jpim1            ! start from 1 due to f-point 
     216            zsshu(ji,jj) = 0.5_wp  * ( sshn(ji  ,jj) + sshn(ji+1,jj  ) ) * ssumask(ji,jj) 
     217            zsshv(ji,jj) = 0.5_wp  * ( sshn(ji  ,jj) + sshn(ji  ,jj+1) ) * ssvmask(ji,jj) 
     218            zsshf(ji,jj) = 0.25_wp * ( sshn(ji  ,jj) + sshn(ji  ,jj+1)   &  
     219               &                     + sshn(ji+1,jj) + sshn(ji+1,jj+1) ) * ssfmask(ji,jj) 
     220         END DO             ;   END DO       
     221         CALL lbc_lnk_multi( zsshu(:,:),'U', 1._wp , zsshu(:,:),'V', 1._wp , zsshf(:,:),'F', 1._wp ) 
     222         ! 
     223         !                                            !* hu and hv (and their inverse)  
     224         ht_n   (:,:) = ht_0(:,:) +  sshn(:,:) 
     225         hu_n   (:,:) = hu_0(:,:) + zsshu(:,:) 
     226         hv_n   (:,:) = hv_0(:,:) + zsshv(:,:) 
     227         r1_hu_n(:,:) = ssumask(:,:) / ( hu_n(:,:) + 1._wp - ssumask(:,:) )    ! ss mask mask due to ISF 
     228         r1_hv_n(:,:) = ssvmask(:,:) / ( hv_n(:,:) + 1._wp - ssvmask(:,:) ) 
     229         ! 
     230         !                                            !* e3u, e3uw  and  e3v, e3vw 
     231         z_ssh_h0(:,:) = sshn (:,:) * r1_ht_0(:,:)           ! t-point 
     232         DO jk = 1, jpkm1 
     233            e3w_n(:,:,jk) = e3w_0(:,:,jk) * ( 1._wp + z_ssh_h0(:,:) * tmask(:,:,jk) ) 
     234         END DO 
     235         z_ssh_h0(:,:) = zsshu(:,:) * r1_hu_0(:,:)           ! u-point 
     236         DO jk = 1, jpkm1 
     237            e3u_n (:,:,jk) = e3u_0 (:,:,jk) * ( 1._wp + z_ssh_h0(:,:) *  umask(:,:,jk) ) 
     238            e3uw_n(:,:,jk) = e3uw_0(:,:,jk) * ( 1._wp + z_ssh_h0(:,:) * wumask(:,:,jk) ) 
     239         END DO 
     240         z_ssh_h0(:,:) = zsshv(:,:) * r1_hv_0(:,:)           ! v-point 
     241         DO jk = 1, jpkm1 
     242            e3v_n (:,:,jk) = e3v_0 (:,:,jk) * ( 1._wp + z_ssh_h0(:,:) *  vmask(:,:,jk) ) 
     243            e3vw_n(:,:,jk) = e3vw_0(:,:,jk) * ( 1._wp + z_ssh_h0(:,:) * wvmask(:,:,jk) ) 
     244         END DO 
     245         z_ssh_h0(:,:) = zsshf(:,:) * r1_hf_0(:,:)           ! f-point 
     246         DO jk = 1, jpkm1 
     247            e3f_n(:,:,jk) = e3f_0(:,:,jk) * ( 1._wp + z_ssh_h0(:,:) * fmask(:,:,jk) ) 
     248         END DO 
     249 
     250         z_ssh_h0(:,:) = 1._wp + sshn(:,:) * r1_ht_0(:,:)    ! t-point 
     251         ! 
     252         IF( ln_isfcav ) THEN    ! iceshelf cavities : ssh scaling not applied over the iceshelf thickness  
     253            DO jk = 1, jpkm1 
     254               gdept_n(:,:,jk) = ( gdept_0(:,:,jk) - risfdep(:,:) ) * z_ssh_h0(:,:) + risfdep(:,:) 
     255               gdepw_n(:,:,jk) = ( gdepw_0(:,:,jk) - risfdep(:,:) ) * z_ssh_h0(:,:) + risfdep(:,:) 
     256               gde3w_n(:,:,jk) = gdept_n(:,:,jk) - sshn(:,:) 
     257            END DO 
     258         ELSE 
     259            DO jk = 1, jpkm1 
     260               gdept_n(:,:,jk) = gdept_0(:,:,jk) * z_ssh_h0(:,:) 
     261               gdepw_n(:,:,jk) = gdepw_0(:,:,jk) * z_ssh_h0(:,:) 
     262               gde3w_n(:,:,jk) = gdept_n(:,:,jk) - sshn(:,:) 
     263            END DO 
     264         ENDIF 
    210265          
    211          ! t- and w- points depth 
    212          ! ---------------------- 
    213          gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1) 
    214          gdepw_n(:,:,1) = 0.0_wp 
    215          gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:) 
    216          DO jj = 1,jpj 
    217             DO ji = 1,jpi 
    218                DO jk = 2,mikt(ji,jj)-1 
    219                   gdept_n(ji,jj,jk) = gdept_0(ji,jj,jk) 
    220                   gdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk) 
    221                   gde3w_n(ji,jj,jk) = gdept_0(ji,jj,jk) - sshn(ji,jj) 
    222                END DO 
    223                IF (mikt(ji,jj) > 1) THEN 
    224                   jk = mikt(ji,jj) 
    225                   gdept_n(ji,jj,jk) = gdepw_0(ji,jj,jk) + 0.5_wp * e3w_n(ji,jj,jk) 
    226                   gdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk) 
    227                   gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk  ) - sshn   (ji,jj) 
    228                END IF 
    229                DO jk = mikt(ji,jj)+1, jpk 
    230                   gdept_n(ji,jj,jk) = gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk) 
    231                   gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1) 
    232                   gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk  ) - sshn (ji,jj) 
    233                END DO 
    234             END DO 
    235          END DO 
    236  
    237       ! t-, u- and v- water column thickness 
    238       ! ------------------------------------ 
    239          ht_n(:,:) = 0._wp ; hu_n(:,:) = 0._wp ; hv_n(:,:) = 0._wp 
    240          DO jk = 1, jpkm1 
    241             hu_n(:,:) = hu_n(:,:) + e3u_n(:,:,jk) * umask(:,:,jk) 
    242             hv_n(:,:) = hv_n(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk) 
    243             ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 
    244          END DO 
    245          !                                        ! Inverse of the local depth 
    246          r1_hu_n(:,:) = 1._wp / ( hu_n(:,:) + 1._wp - ssumask(:,:) ) * ssumask(:,:) 
    247          r1_hv_n(:,:) = 1._wp / ( hv_n(:,:) + 1._wp - ssvmask(:,:) ) * ssvmask(:,:) 
    248  
    249       END IF 
     266      ENDIF 
    250267 
    251268!============================================================================= 
    252269! compute velocity 
    253270! compute velocity in order to conserve barotropic velocity (modification by poderation of the scale factor). 
    254       ub(:,:,:)=un(:,:,:) 
    255       vb(:,:,:)=vn(:,:,:) 
    256       DO jk = 1,jpk 
    257          DO jj = 1,jpj 
    258             DO ji = 1,jpi 
     271      ub(:,:,:) = un(:,:,:) 
     272      vb(:,:,:) = vn(:,:,:) 
     273      DO jk = 1, jpk 
     274         DO jj = 1, jpj 
     275            DO ji = 1, jpi 
    259276               un(ji,jj,jk) = ub(ji,jj,jk)*pe3u_b(ji,jj,jk)*pumask_b(ji,jj,jk)/e3u_n(ji,jj,jk)*umask(ji,jj,jk); 
    260277               vn(ji,jj,jk) = vb(ji,jj,jk)*pe3v_b(ji,jj,jk)*pvmask_b(ji,jj,jk)/e3v_n(ji,jj,jk)*vmask(ji,jj,jk); 
    … …  
    265282! compute new velocity if we close a cell (check barotropic velocity and change velocity over the water column) 
    266283! compute barotropic velocity now and after  
    267       ztrp(:,:,:) = ub(:,:,:)*pe3u_b(:,:,:);  
    268       zbub(:,:)   = SUM(ztrp,DIM=3) 
    269       ztrp(:,:,:) = vb(:,:,:)*pe3v_b(:,:,:);  
    270       zbvb(:,:)   = SUM(ztrp,DIM=3) 
    271       ztrp(:,:,:) = un(:,:,:)*e3u_n(:,:,:);  
    272       zbun(:,:)   = SUM(ztrp,DIM=3) 
    273       ztrp(:,:,:) = vn(:,:,:)*e3v_n(:,:,:);  
    274       zbvn(:,:)   = SUM(ztrp,DIM=3) 
     284      zbub(:,:)   = SUM( ub(:,:,:)*pe3u_b(:,:,:) , DIM=3 ) 
     285      zbvb(:,:)   = SUM( vb(:,:,:)*pe3v_b(:,:,:) , DIM=3 ) 
     286      zbun(:,:)   = SUM( un(:,:,:)* e3u_n(:,:,:) , DIM=3 ) 
     287      zbvn(:,:)   = SUM( vn(:,:,:)* e3v_n(:,:,:) , DIM=3 ) 
    275288 
    276289      ! new water column 
    … …  
    285298      zucorr = 0._wp 
    286299      zvcorr = 0._wp 
    287       DO jj = 1,jpj 
    288          DO ji = 1,jpi 
     300      DO jj = 1, jpj 
     301         DO ji = 1, jpi 
    289302            IF (zbun(ji,jj) /= zbub(ji,jj) .AND. zhu1(ji,jj) /= 0._wp ) THEN 
    290303               zucorr(ji,jj) = (zbun(ji,jj) - zbub(ji,jj))/zhu1(ji,jj) 
    … …  
    297310       
    298311      ! update velocity 
    299       DO jk  = 1,jpk 
    300          un(:,:,jk)=(un(:,:,jk) - zucorr(:,:))*umask(:,:,jk) 
    301          vn(:,:,jk)=(vn(:,:,jk) - zvcorr(:,:))*vmask(:,:,jk) 
     312      DO jk  = 1, jpk 
     313         un(:,:,jk) = ( un(:,:,jk) - zucorr(:,:) ) * umask(:,:,jk) 
     314         vn(:,:,jk) = ( vn(:,:,jk) - zvcorr(:,:) ) * vmask(:,:,jk) 
    302315      END DO 
    303316 
    … …  
    305318      ! compute temp and salt 
    306319      ! compute new tn and sn if we open a new cell 
    307       tsb (:,:,:,:) = tsn(:,:,:,:) 
    308       zts0(:,:,:,:) = tsn(:,:,:,:) 
     320      tsb  (:,:,:,:) = tsn   (:,:,:,:) 
     321      zts0 (:,:,:,:) = tsn   (:,:,:,:) 
    309322      ztmask1(:,:,:) = ptmask_b(:,:,:) 
    310323      ztmask0(:,:,:) = ptmask_b(:,:,:) 
    311324      DO iz = 1,nn_drown ! resolution dependent (OK for ISOMIP+ case) 
    312           DO jk = 1,jpk-1 
    313               zdmask=tmask(:,:,jk)-ztmask0(:,:,jk); 
    314               DO jj = 2,jpj-1 
    315                  DO ji = fs_2,fs_jpim1 
    316                       jip1=ji+1; jim1=ji-1; 
    317                       jjp1=jj+1; jjm1=jj-1; 
    318                       summsk= (ztmask0(jip1,jj  ,jk)+ztmask0(jim1,jj  ,jk)+ztmask0(ji  ,jjp1,jk)+ztmask0(ji  ,jjm1,jk)) 
    319                       IF (zdmask(ji,jj) == 1._wp .AND. summsk /= 0._wp) THEN 
    320                       !! horizontal basic extrapolation 
    321                          tsn(ji,jj,jk,1)=( zts0(jip1,jj  ,jk,1)*ztmask0(jip1,jj  ,jk) & 
    322                          &                +zts0(jim1,jj  ,jk,1)*ztmask0(jim1,jj  ,jk) & 
    323                          &                +zts0(ji  ,jjp1,jk,1)*ztmask0(ji  ,jjp1,jk) & 
    324                          &                +zts0(ji  ,jjm1,jk,1)*ztmask0(ji  ,jjm1,jk) ) / summsk 
    325                          tsn(ji,jj,jk,2)=( zts0(jip1,jj  ,jk,2)*ztmask0(jip1,jj  ,jk) & 
    326                          &                +zts0(jim1,jj  ,jk,2)*ztmask0(jim1,jj  ,jk) & 
    327                          &                +zts0(ji  ,jjp1,jk,2)*ztmask0(ji  ,jjp1,jk) & 
    328                          &                +zts0(ji  ,jjm1,jk,2)*ztmask0(ji  ,jjm1,jk) ) / summsk 
    329                          ztmask1(ji,jj,jk)=1 
    330                       ELSEIF (zdmask(ji,jj) == 1._wp .AND. summsk == 0._wp) THEN 
    331                       !! vertical extrapolation if horizontal extrapolation failed 
    332                          jkm1=max(1,jk-1) ; jkp1=min(jpk,jk+1) 
    333                          summsk=(ztmask0(ji,jj,jkm1)+ztmask0(ji,jj,jkp1)) 
    334                          IF (zdmask(ji,jj) == 1._wp .AND. summsk /= 0._wp ) THEN 
    335                             tsn(ji,jj,jk,1)=( zts0(ji,jj,jkp1,1)*ztmask0(ji,jj,jkp1)     & 
    336                             &                +zts0(ji,jj,jkm1,1)*ztmask0(ji,jj,jkm1))/summsk 
    337                             tsn(ji,jj,jk,2)=( zts0(ji,jj,jkp1,2)*ztmask0(ji,jj,jkp1)     & 
    338                             &                +zts0(ji,jj,jkm1,2)*ztmask0(ji,jj,jkm1))/summsk 
    339                             ztmask1(ji,jj,jk)=1._wp 
    340                          END IF 
    341                       END IF 
    342                   END DO 
    343               END DO 
    344           END DO 
    345            
    346           CALL lbc_lnk_multi( tsn(:,:,:,jp_tem), 'T', 1., tsn(:,:,:,jp_sal), 'T', 1., ztmask1, 'T', 1.) 
    347  
    348           ! update 
    349           zts0(:,:,:,:) = tsn(:,:,:,:) 
    350           ztmask0 = ztmask1 
    351  
     325         DO jk = 1,jpk-1 
     326            zdmask=tmask(:,:,jk)-ztmask0(:,:,jk); 
     327            DO jj = 2,jpj-1 
     328               DO ji = fs_2,fs_jpim1 
     329                  jip1=ji+1; jim1=ji-1; 
     330                  jjp1=jj+1; jjm1=jj-1; 
     331                  summsk= (ztmask0(jip1,jj  ,jk)+ztmask0(jim1,jj  ,jk)+ztmask0(ji  ,jjp1,jk)+ztmask0(ji  ,jjm1,jk)) 
     332                  IF (zdmask(ji,jj) == 1._wp .AND. summsk /= 0._wp) THEN 
     333                     !! horizontal basic extrapolation 
     334                     tsn(ji,jj,jk,1)=( zts0(jip1,jj  ,jk,1)*ztmask0(jip1,jj  ,jk) & 
     335                         &            +zts0(jim1,jj  ,jk,1)*ztmask0(jim1,jj  ,jk) & 
     336                         &            +zts0(ji  ,jjp1,jk,1)*ztmask0(ji  ,jjp1,jk) & 
     337                         &            +zts0(ji  ,jjm1,jk,1)*ztmask0(ji  ,jjm1,jk) ) / summsk 
     338                     tsn(ji,jj,jk,2)=( zts0(jip1,jj  ,jk,2)*ztmask0(jip1,jj  ,jk) & 
     339                         &            +zts0(jim1,jj  ,jk,2)*ztmask0(jim1,jj  ,jk) & 
     340                         &            +zts0(ji  ,jjp1,jk,2)*ztmask0(ji  ,jjp1,jk) & 
     341                         &            +zts0(ji  ,jjm1,jk,2)*ztmask0(ji  ,jjm1,jk) ) / summsk 
     342                     ztmask1(ji,jj,jk)=1 
     343                  ELSEIF (zdmask(ji,jj) == 1._wp .AND. summsk == 0._wp) THEN 
     344                     !! vertical extrapolation if horizontal extrapolation failed 
     345                     jkm1=max(1,jk-1) ; jkp1=min(jpk,jk+1) 
     346                     summsk=(ztmask0(ji,jj,jkm1)+ztmask0(ji,jj,jkp1)) 
     347                     IF (zdmask(ji,jj) == 1._wp .AND. summsk /= 0._wp ) THEN 
     348                        tsn(ji,jj,jk,1)=( zts0(ji,jj,jkp1,1)*ztmask0(ji,jj,jkp1)     & 
     349                            &            +zts0(ji,jj,jkm1,1)*ztmask0(ji,jj,jkm1)) / summsk 
     350                        tsn(ji,jj,jk,2)=( zts0(ji,jj,jkp1,2)*ztmask0(ji,jj,jkp1)     & 
     351                            &            +zts0(ji,jj,jkm1,2)*ztmask0(ji,jj,jkm1)) / summsk 
     352                        ztmask1(ji,jj,jk)=1._wp 
     353                     ENDIF 
     354                  ENDIF 
     355               END DO 
     356            END DO 
     357         END DO 
     358         ! 
     359         CALL lbc_lnk_multi( tsn(:,:,:,jp_tem), 'T', 1., tsn(:,:,:,jp_sal), 'T', 1., ztmask1, 'T', 1.) 
     360         ! 
     361         ! update 
     362         zts0(:,:,:,:) = tsn(:,:,:,:) 
     363         ztmask0 = ztmask1 
     364         ! 
    352365      END DO 
    353366 
    … …  
    395408      ! set mbkt and mikt to 1 in thiese location 
    396409      WHERE (SUM(tmask,dim=3) == 0) 
    397          mbkt(:,:)=1 ; mbku(:,:)=1 ; mbkv(:,:)=1 
    398          mikt(:,:)=1 ; miku(:,:)=1 ; mikv(:,:)=1 
     410         mbkt(:,:) = 1   ;   mbku(:,:)=1   ;   mbkv(:,:) = 1 
     411         mikt(:,:) = 1   ;   miku(:,:)=1   ;   mikv(:,:) = 1 
    399412      END WHERE 
    400413      ! ------------------------------------------------------------------------------------------- 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DYN/dynadv.F90

    r9598 r10001  
    121121      ENDIF 
    122122 
     123!!gm  
     124      IF ( ln_dynadv_vec )   CALL ctl_stop( 'Vector form not available in this branch' ) 
     125!!gm 
     126 
    123127      ioptio = 0                      ! parameter control and set n_dynadv 
    124128      IF( ln_dynadv_OFF  ) THEN   ;   ioptio = ioptio + 1   ;   n_dynadv = np_LIN_dyn   ;   ENDIF 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DYN/dynnxt.F90

    r9939 r10001  
    2020   !!            3.6  !  2014-04  (G. Madec) add the diagnostic of the time filter trends 
    2121   !!            3.7  !  2015-11  (J. Chanut) Free surface simplification 
     22   !!            4.0  !  2018-07  (G. Madec)   1- z-star (s-star) only   
     23   !!                                          2- remove   dom_vvl_interpol   
    2224   !!------------------------------------------------------------------------- 
    2325   
    … …  
    3436   USE dynspg_ts      ! surface pressure gradient: split-explicit scheme 
    3537   USE domvvl         ! variable volume 
    36    USE bdy_oce   , ONLY: ln_bdy 
     38   USE bdy_oce , ONLY : ln_bdy 
    3739   USE bdydta         ! ocean open boundary conditions 
    3840   USE bdydyn         ! ocean open boundary conditions 
    … …  
    9294      !!               un,vn   now horizontal velocity of next time-step 
    9395      !!---------------------------------------------------------------------- 
    94       INTEGER, INTENT( in ) ::   kt   ! ocean time-step index 
     96      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
    9597      ! 
    9698      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    … …  
    100102      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zue, zve 
    101103      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ze3u_f, ze3v_f, zua, zva  
     104      REAL(wp), DIMENSION(jpi,jpj) ::   z_ssh_h0, zsshu, zsshv 
    102105      !!---------------------------------------------------------------------- 
    103106      ! 
    … …  
    212215            !    => time filter + conservation correction (only at the first level) 
    213216            zcoef = rn_atfp * rn_Dt * r1_rho0 
    214  
     217            ! 
    215218            e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1) 
    216  
     219            ! 
    217220            IF ( ln_rnf ) THEN 
    218221               IF( ln_rnf_depth ) THEN 
    … …  
    228231                  END DO 
    229232               ELSE 
    230                   e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef *  ( -rnf_b(:,:) + rnf(:,:))*tmask(:,:,1) 
     233                  e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef *  ( -rnf_b(:,:) + rnf(:,:) )*tmask(:,:,1) 
    231234               ENDIF 
    232235            ENDIF 
    233  
     236            ! 
    234237            IF ( ln_isf ) THEN   ! if ice shelf melting 
    235238               DO jk = 1, jpkm1 ! Deal with isf separetely, as can be through depth too 
    … …  
    245248            ENDIF 
    246249            ! 
     250            ! 
    247251            IF( ln_dynadv_vec ) THEN      ! Asselin filter applied on velocity 
    248                ! Before filtered scale factor at (u/v)-points 
    249                CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' ) 
    250                CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' ) 
     252               !                                            !* ssh at u- and v-points) 
     253               DO jj = 2, jpjm1   ;   DO ji = 2, jpim1 
     254                  zsshu(ji,jj) = 0.5_wp  * ( sshb(ji  ,jj) + sshb(ji+1,jj  ) ) * ssumask(ji,jj) 
     255                  zsshv(ji,jj) = 0.5_wp  * ( sshb(ji  ,jj) + sshb(ji  ,jj+1) ) * ssvmask(ji,jj) 
     256               END DO             ;   END DO       
     257               CALL lbc_lnk_multi( zsshu(:,:),'U', 1._wp , zsshu(:,:),'V', 1._wp ) 
     258               ! 
     259               ! 
     260               !                                            !* e3u and e3v  
     261               z_ssh_h0(:,:) = zsshu(:,:) * r1_hu_0(:,:)           ! u-point 
     262               DO jk = 1, jpkm1 
     263                  e3u_b (:,:,jk) = e3u_0 (:,:,jk) * ( 1._wp + z_ssh_h0(:,:) * umask(:,:,jk) ) 
     264               END DO 
     265               z_ssh_h0(:,:) = zsshv(:,:) * r1_hv_0(:,:)           ! v-point 
     266               DO jk = 1, jpkm1 
     267                  e3v_b (:,:,jk) = e3v_0 (:,:,jk) * ( 1._wp + z_ssh_h0(:,:) * vmask(:,:,jk) ) 
     268               END DO 
     269               ! 
    251270               DO jk = 1, jpkm1 
    252271                  DO jj = 1, jpj 
    … …  
    266285               ! 
    267286               ALLOCATE( ze3u_f(jpi,jpj,jpk) , ze3v_f(jpi,jpj,jpk) ) 
    268                ! Before filtered scale factor at (u/v)-points stored in ze3u_f, ze3v_f 
    269                CALL dom_vvl_interpol( e3t_b(:,:,:), ze3u_f, 'U' ) 
    270                CALL dom_vvl_interpol( e3t_b(:,:,:), ze3v_f, 'V' ) 
     287               ! 
     288               !                                            !* ssh at u- and v-points) 
     289               DO jj = 2, jpjm1   ;   DO ji = 2, jpim1 
     290                  zsshu(ji,jj) = 0.5_wp  * ( sshb(ji  ,jj) + sshb(ji+1,jj  ) ) * ssumask(ji,jj) 
     291                  zsshv(ji,jj) = 0.5_wp  * ( sshb(ji  ,jj) + sshb(ji  ,jj+1) ) * ssvmask(ji,jj) 
     292               END DO             ;   END DO       
     293               CALL lbc_lnk_multi( zsshu(:,:),'U', 1._wp , zsshu(:,:),'V', 1._wp ) 
     294               ! 
     295               ! 
     296               !                                            !* e3u and e3v  
     297               z_ssh_h0(:,:) = zsshu(:,:) * r1_hu_0(:,:)           ! u-point 
     298               DO jk = 1, jpkm1 
     299                  ze3u_f(:,:,jk) = e3u_0 (:,:,jk) * ( 1._wp + z_ssh_h0(:,:) * umask(:,:,jk) ) 
     300               END DO 
     301               z_ssh_h0(:,:) = zsshv(:,:) * r1_hv_0(:,:)           ! v-point 
     302               DO jk = 1, jpkm1 
     303                  ze3u_f(:,:,jk) = e3v_0 (:,:,jk) * ( 1._wp + z_ssh_h0(:,:) * vmask(:,:,jk) ) 
     304               END DO 
     305               ! 
    271306               DO jk = 1, jpkm1 
    272307                  DO jj = 1, jpj 
    … …  
    319354      ! 
    320355      IF(.NOT.ln_linssh ) THEN 
    321          hu_b(:,:) = e3u_b(:,:,1) * umask(:,:,1) 
    322          hv_b(:,:) = e3v_b(:,:,1) * vmask(:,:,1) 
    323          DO jk = 2, jpkm1 
    324             hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk) 
    325             hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk) 
    326          END DO 
    327          r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) ) 
     356         DO jj = 2, jpjm1   ;   DO ji = 2, jpim1 
     357            zsshu(ji,jj) = 0.5_wp  * ( sshb(ji  ,jj) + sshb(ji+1,jj  ) ) * ssumask(ji,jj) 
     358            zsshv(ji,jj) = 0.5_wp  * ( sshb(ji  ,jj) + sshb(ji  ,jj+1) ) * ssvmask(ji,jj) 
     359         END DO             ;   END DO       
     360         CALL lbc_lnk_multi( zsshu(:,:),'U', 1._wp , zsshu(:,:),'V', 1._wp ) 
     361         ! 
     362         hu_b   (:,:) = hu_0(:,:) + zsshu(:,:) 
     363         hv_b   (:,:) = hv_0(:,:) + zsshv(:,:) 
     364         r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) )    ! _i mask due to ISF 
    328365         r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) ) 
    329366      ENDIF 
    330367      ! 
     368      !                                   !  
    331369      un_b(:,:) = e3u_a(:,:,1) * un(:,:,1) * umask(:,:,1) 
    332370      ub_b(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1) 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DYN/dynzdf.F90

    r9939 r10001  
    458458         !                                ! Note also that now e3 scale factors are used as after one are not computed ! 
    459459         ! 
    460          CALL wrk_alloc(jpi,jpj,   z2d ) 
     460         ALLOCATE( z2d(jpi,jpj) ) 
    461461         z2d(:,:) = 0._wp 
    462462         DO jk = 1, jpkm1 
    … …  
    498498         CALL lbc_lnk( z2d,'T', 1. ) 
    499499         CALL iom_put( 'dispkevfo', z2d ) 
     500         DEALLOCATE( z2d ) 
    500501      ENDIF 
    501502      ! 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/IOM/in_out_manager.F90

    r9598 r10001  
    3535   INTEGER       ::   nn_time0         !: initial time of day in hhmm 
    3636   INTEGER       ::   nn_leapy         !: Leap year calendar flag (0/1 or 30) 
    37    INTEGER       ::   nn_istate        !: initial state output flag (0/1) 
     37   LOGICAL       ::   ln_istate        !: initial state output flag 
    3838   INTEGER       ::   nn_write         !: model standard output frequency 
    3939   INTEGER       ::   nn_stock         !: restart file frequency 
    … …  
    7979   INTEGER       ::   ndate0                      !: initial calendar date aammjj 
    8080   INTEGER       ::   nleapy                      !: Leap year calendar flag (0/1 or 30) 
    81    INTEGER       ::   ninist                      !: initial state output flag (0/1) 
    8281   INTEGER       ::   nwrite                      !: model standard output frequency 
    8382   INTEGER       ::   nstock                      !: restart file frequency 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/SBC/sbcice_cice.F90

    r9939 r10001  
    205205         DO jl=1,ncat 
    206206            CALL cice2nemo(aicen(:,:,jl,:),a_i(:,:,jl), 'T', 1. ) 
    207          ENDDO 
     207         END DO 
    208208      ENDIF 
    209209 
    … …  
    214214            fr_iu(ji,jj)=0.5*(fr_i(ji,jj)+fr_i(ji+1,jj))*umask(ji,jj,1) 
    215215            fr_iv(ji,jj)=0.5*(fr_i(ji,jj)+fr_i(ji,jj+1))*vmask(ji,jj,1) 
    216          ENDDO 
    217       ENDDO 
     216         END DO 
     217      END DO 
    218218 
    219219      CALL lbc_lnk_multi( fr_iu , 'U', 1.,  fr_iv , 'V', 1. ) 
    … …  
    235235               ! 
    236236               DO jk = 1,jpkm1                     ! adjust initial vertical scale factors 
    237                   e3t_n(:,:,jk) = e3t_0(:,:,jk)*( 1._wp + sshn(:,:)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) ) 
    238                   e3t_b(:,:,jk) = e3t_0(:,:,jk)*( 1._wp + sshb(:,:)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) ) 
    239                ENDDO 
     237                  e3t_n(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + sshn(:,:) * r1_ht_0(:,:) ) 
     238                  e3t_b(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + sshb(:,:) * r1_ht_0(:,:) ) 
     239               END DO 
    240240               e3t_a(:,:,:) = e3t_b(:,:,:) 
    241241               ! Reconstruction of all vertical scale factors at now and before time-steps 
    … …  
    307307               ztmp(ji,jj) = 0.5 * (  fr_iu(ji,jj) * utau(ji,jj)      & 
    308308                                    + fr_iu(ji,jj+1) * utau(ji,jj+1) ) * fmask(ji,jj,1) 
    309             ENDDO 
    310          ENDDO 
     309            END DO 
     310         END DO 
    311311         CALL nemo2cice(ztmp,strax,'F', -1. ) 
    312312 
    … …  
    317317               ztmp(ji,jj) = 0.5 * (  fr_iv(ji,jj) * vtau(ji,jj)      & 
    318318                                    + fr_iv(ji+1,jj) * vtau(ji+1,jj) ) * fmask(ji,jj,1) 
    319             ENDDO 
    320          ENDDO 
     319            END DO 
     320         END DO 
    321321         CALL nemo2cice(ztmp,stray,'F', -1. ) 
    322322 
    … …  
    325325            DO jl=1,ncat 
    326326               ztmpn(:,:,jl)=qla_ice(:,:,1)*a_i(:,:,jl) 
    327             ENDDO 
     327            END DO 
    328328         ELSE 
    329329! emp_ice is set in sbc_cpl_ice_flx as sublimation-snow 
    … …  
    332332            DO jj=1,jpj 
    333333               DO ji=1,jpi 
    334                   IF (fr_i(ji,jj).eq.0.0) THEN 
     334                  IF (fr_i(ji,jj) == 0._wp ) THEN 
    335335                     DO jl=1,ncat 
    336                         ztmpn(ji,jj,jl)=0.0 
    337                      ENDDO 
     336                        ztmpn(ji,jj,jl) = 0._wp 
     337                     END DO 
    338338                     ! This will then be conserved in CICE 
    339339                     ztmpn(ji,jj,1)=qla_ice(ji,jj,1) 
    … …  
    341341                     DO jl=1,ncat 
    342342                        ztmpn(ji,jj,jl)=qla_ice(ji,jj,1)*a_i(ji,jj,jl)/fr_i(ji,jj) 
    343                      ENDDO 
     343                     END DO 
    344344                  ENDIF 
    345                ENDDO 
    346             ENDDO 
     345               END DO 
     346            END DO 
    347347         ENDIF 
    348348         DO jl=1,ncat 
    … …  
    366366            ENDIF 
    367367            CALL nemo2cice(ztmp,fsurfn_f(:,:,jl,:),'T', 1. ) 
    368          ENDDO 
     368         END DO 
    369369 
    370370      ELSE IF (ksbc == jp_blk) THEN 
    … …  
    437437         DO ji=1,jpi 
    438438            ztmp(ji,jj)=0.5*(ssu_m(ji,jj)+ssu_m(ji,jj+1))*fmask(ji,jj,1) 
    439          ENDDO 
    440       ENDDO 
     439         END DO 
     440      END DO 
    441441      CALL nemo2cice(ztmp,uocn,'F', -1. ) 
    442442 
    … …  
    445445         DO ji=1,jpim1 
    446446            ztmp(ji,jj)=0.5*(ssv_m(ji,jj)+ssv_m(ji+1,jj))*fmask(ji,jj,1) 
    447          ENDDO 
    448       ENDDO 
     447         END DO 
     448      END DO 
    449449      CALL nemo2cice(ztmp,vocn,'F', -1. ) 
    450450 
    … …  
    511511         DO ji=2,jpim1 
    512512            ss_iou(ji,jj) = 0.5 * ( ztmp1(ji,jj-1) + ztmp1(ji,jj) ) * umask(ji,jj,1) 
    513          ENDDO 
    514       ENDDO 
     513         END DO 
     514      END DO 
    515515      CALL lbc_lnk( ss_iou , 'U', -1. ) 
    516516 
    … …  
    523523         DO ji=2,jpim1 
    524524            ss_iov(ji,jj) = 0.5 * ( ztmp1(ji-1,jj) + ztmp1(ji,jj) ) * vmask(ji,jj,1) 
    525          ENDDO 
    526       ENDDO 
     525         END DO 
     526      END DO 
    527527      CALL lbc_lnk( ss_iov , 'V', -1. ) 
    528528 
    … …  
    609609         DO ji=1,jpi 
    610610            nfrzmlt(ji,jj)=MAX(nfrzmlt(ji,jj),0.0) 
    611          ENDDO 
    612       ENDDO 
     611         END DO 
     612      END DO 
    613613 
    614614#if defined key_cice4 
    … …  
    627627         DO jl=1,ncat 
    628628            CALL cice2nemo(aicen(:,:,jl,:),a_i(:,:,jl), 'T', 1. ) 
    629          ENDDO 
     629         END DO 
    630630      ENDIF 
    631631 
    … …  
    636636            fr_iu(ji,jj)=0.5*(fr_i(ji,jj)+fr_i(ji+1,jj))*umask(ji,jj,1) 
    637637            fr_iv(ji,jj)=0.5*(fr_i(ji,jj)+fr_i(ji,jj+1))*vmask(ji,jj,1) 
    638          ENDDO 
    639       ENDDO 
     638         END DO 
     639      END DO 
    640640 
    641641      CALL lbc_lnk_multi( fr_iu , 'U', 1., fr_iv , 'V', 1. ) 
    … …  
    872872         DO ji=2,nx_block-1 
    873873            pc(ji,jj,1)=pn(ji-1+ji_off,jj-1+jj_off) 
    874          ENDDO 
    875       ENDDO 
     874         END DO 
     875      END DO 
    876876 
    877877#else 
    … …  
    898898               DO ji=nldit(jn),nleit(jn) 
    899899                  png2(ji+nimppt(jn)-1,jj+njmppt(jn)-1)=png(ji,jj,jn) 
    900                ENDDO 
    901             ENDDO 
    902          ENDDO 
     900               END DO 
     901            END DO 
     902         END DO 
    903903         DO jj=1,ny_global 
    904904            DO ji=1,nx_global 
    905905               pcg(ji,jj)=png2(ji+ji_off,jj+jj_off) 
    906             ENDDO 
    907          ENDDO 
     906            END DO 
     907         END DO 
    908908      ENDIF 
    909909 
    … …  
    999999         DO ji=1,jpim1 
    10001000            pn(ji,jj)=pc(ji+1-ji_off,jj+1-jj_off,1) 
    1001          ENDDO 
    1002       ENDDO 
     1001         END DO 
     1002      END DO 
    10031003 
    10041004#else 
    … …  
    10211021               DO ji=nldit(jn),nleit(jn) 
    10221022                  png(ji,jj,jn)=pcg(ji+nimppt(jn)-1-ji_off,jj+njmppt(jn)-1-jj_off) 
    1023                ENDDO 
    1024             ENDDO 
    1025          ENDDO 
     1023               END DO 
     1024            END DO 
     1025         END DO 
    10261026      ENDIF 
    10271027 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/TRA/trazdf.F90

    r9939 r10001  
    175175               DO jj = 2, jpjm1 
    176176                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     177!!gm BUG   here e3w_a  should be used !!!!!   but then should be added in the system  
    177178                     zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk  ) / e3w_n(ji,jj,jk  ) 
    178179                     zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / e3w_n(ji,jj,jk+1) 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/nemogcm.F90

    r9939 r10001  
    3636   !!---------------------------------------------------------------------- 
    3737   !!   nemo_gcm      : solve ocean dynamics, tracer, biogeochemistry and/or sea-ice 
     38   !!   nemo_MLF      : Modified Leap-Frog time stepping loop 
     39   !!   nemo_RK3      : 3rd order Runge-Kutta time stepping loop 
    3840   !!   nemo_init     : initialization of the NEMO system 
    3941   !!   nemo_ctl      : initialisation of the contol print 
    … …  
    6062   USE diacfl         ! CFL diagnostics               (dia_cfl_init routine) 
    6163   USE step           ! NEMO time-stepping                 (stp     routine) 
     64   USE stpRK3         ! NEMO time-stepping                 (stp     routine) 
    6265   USE icbini         ! handle bergs, initialisation 
    6366   USE icbstp         ! handle bergs, calving, themodynamics and transport 
    … …  
    123126      !!              Madec, 2008, internal report, IPSL. 
    124127      !!---------------------------------------------------------------------- 
    125       INTEGER ::   istp   ! time step index 
    126       !!---------------------------------------------------------------------- 
    127128      ! 
    128129#if defined key_agrif 
    … …  
    152153      !                            !-----------------------! 
    153154      ! 
    154       !                                               !== set the model time-step  ==! 
     155      IF( ln_MLF )   CALL nemo_MLF        ! Modified Leap-Frog 
     156      ! 
     157      IF( ln_RK3 )   CALL nemo_RK3        ! 3rd order Runge-Kutta 
     158      ! 
     159      ! 
     160      IF( ln_diaobs   )   CALL dia_obs_wri 
     161      ! 
     162      IF( ln_icebergs )   CALL icb_end( nitend ) 
     163 
     164      !                            !------------------------! 
     165      !                            !==  finalize the run  ==! 
     166      !                            !------------------------! 
     167      IF(lwp) WRITE(numout,cform_aaa)        ! Flag AAAAAAA 
     168      ! 
     169      IF( nstop /= 0 .AND. lwp ) THEN        ! error print 
     170         WRITE(numout,cform_err) 
     171         WRITE(numout,*) '   ==>>>   nemo_gcm: a total of ', nstop, ' errors have been found' 
     172         WRITE(numout,*) 
     173      ENDIF 
     174      ! 
     175      IF( ln_timing )   CALL timing_finalize 
     176      ! 
     177      CALL nemo_closefile 
     178      ! 
     179#if defined key_iomput 
     180                                    CALL xios_finalize  ! end mpp communications with xios 
     181      IF( lk_oasis     )            CALL cpl_finalize   ! end coupling and mpp communications with OASIS 
     182#else 
     183      IF    ( lk_oasis ) THEN   ;   CALL cpl_finalize   ! end coupling and mpp communications with OASIS 
     184      ELSEIF( lk_mpp   ) THEN   ;   CALL mppstop        ! end mpp communications 
     185      ENDIF 
     186#endif 
     187      ! 
     188   END SUBROUTINE nemo_gcm 
     189 
     190 
     191   SUBROUTINE nemo_MLF 
     192      !!---------------------------------------------------------------------- 
     193      !!                     ***  ROUTINE nemo_MLF  *** 
     194      !! 
     195      !! ** Purpose :   Modified Leap-Frog time stepping loop 
     196      !!---------------------------------------------------------------------- 
     197      INTEGER ::   istp   ! time step index 
     198      !!---------------------------------------------------------------------- 
     199      !                                               !== set the MLF time-step  ==! 
    155200      ! 
    156201      IF( l_1st_euler ) THEN   ;   rDt =         rn_Dt   ;   l_1st_euler = .TRUE.    ! start or restart with Euler 1st time-step 
    … …  
    212257#endif 
    213258      ! 
    214       IF( ln_diaobs   )   CALL dia_obs_wri 
    215       ! 
    216       IF( ln_icebergs )   CALL icb_end( nitend ) 
    217  
    218       !                            !------------------------! 
    219       !                            !==  finalize the run  ==! 
    220       !                            !------------------------! 
    221       IF(lwp) WRITE(numout,cform_aaa)        ! Flag AAAAAAA 
    222       ! 
    223       IF( nstop /= 0 .AND. lwp ) THEN        ! error print 
    224          WRITE(numout,cform_err) 
    225          WRITE(numout,*) '   ==>>>   nemo_gcm: a total of ', nstop, ' errors have been found' 
    226          WRITE(numout,*) 
    227       ENDIF 
    228       ! 
    229       IF( ln_timing )   CALL timing_finalize 
    230       ! 
    231       CALL nemo_closefile 
    232       ! 
    233 #if defined key_iomput 
    234                                     CALL xios_finalize  ! end mpp communications with xios 
    235       IF( lk_oasis     )            CALL cpl_finalize   ! end coupling and mpp communications with OASIS 
     259   END SUBROUTINE nemo_MLF 
     260 
     261 
     262   SUBROUTINE nemo_RK3 
     263      !!---------------------------------------------------------------------- 
     264      !!                     ***  ROUTINE nemo_RK3  *** 
     265      !! 
     266      !! ** Purpose :   3rd order Runge-Kutta time stepping loop 
     267      !!---------------------------------------------------------------------- 
     268      INTEGER ::   istp   ! time step index 
     269      !!---------------------------------------------------------------------- 
     270      ! 
     271      !                                         !== set the model time-step  ==! 
     272      rDt = rn_Dt 
     273      r1_Dt = 1._wp / rDt 
     274      istp = nit000 
     275      ! 
     276#if defined key_agrif 
     277      !                                         !==  AGRIF time-stepping  ==! 
     278      CALL Agrif_Regrid() 
     279      ! 
     280      CALL Agrif_step_child_adj(Agrif_Update_All)        ! Recursive update from highest to lowest nested levels 
     281      ! 
     282      DO WHILE( istp <= nitend .AND. nstop == 0 ) 
     283         CALL stp_RK3 
     284         istp = istp + 1 
     285      END DO 
     286      ! 
     287      IF( .NOT. Agrif_Root() ) THEN 
     288         CALL Agrif_ParentGrid_To_ChildGrid() 
     289         IF( ln_diaobs )   CALL dia_obs_wri 
     290         IF( ln_timing )   CALL timing_finalize 
     291         CALL Agrif_ChildGrid_To_ParentGrid() 
     292      ENDIF 
     293      ! 
    236294#else 
    237       IF    ( lk_oasis ) THEN   ;   CALL cpl_finalize   ! end coupling and mpp communications with OASIS 
    238       ELSEIF( lk_mpp   ) THEN   ;   CALL mppstop        ! end mpp communications 
    239       ENDIF 
    240 #endif 
    241       ! 
    242    END SUBROUTINE nemo_gcm 
     295      !                                         !==  Standard time-stepping  ==! 
     296      ! 
     297      DO WHILE( istp <= nitend .AND. nstop == 0 ) 
     298         CALL stp_RK3( istp )  
     299         istp = istp + 1 
     300      END DO 
     301      ! 
     302#endif 
     303      ! 
     304   END SUBROUTINE nemo_RK3 
    243305 
    244306 
    … …  
    459521 
    460522      !                                      ! Assimilation increments 
    461       IF( lk_asminc    )   CALL asm_inc_init    ! Initialize assimilation increments 
     523      IF( lk_asminc    ) THEN 
     524                           CALL asm_inc_init    ! Initialize assimilation increments 
     525                           l_dynasm = ln_asmiau .AND. ln_dyninc 
     526                           l_traasm = ln_asmiau .AND. ln_trainc 
     527      ELSE 
     528                           l_dynasm = .FALSE. 
     529                           l_traasm = .FALSE. 
     530      ENDIF 
    462531      ! 
    463532      IF(lwp) WRITE(numout,cform_aaa)           ! Flag AAAAAAA 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/oce.F90

    r9939 r10001  
    3434   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ub_b   ,  un_b  ,  ua_b  !: Barotropic velocities at u-point [m/s] 
    3535   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   vb_b   ,  vn_b  ,  va_b  !: Barotropic velocities at v-point [m/s] 
    36    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sshb   ,  sshn  ,  ssha  !: sea surface height at t-point [m] 
     36   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sshb   ,  sshn  , ssha   !: sea surface height at t-point [m] 
     37 
     38!!gm??   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ssh  , sshu   , sshv   !: sea surface height at t-, u- v-points [m] 
     39 
    3740 
    3841   !! Arrays at barotropic time step:                   ! befbefore! before !  now   ! after  ! 
    … …  
    9093         &      rhd  (jpi,jpj,jpk)      , rhop (jpi,jpj,jpk)                              , STAT=ierr(1) ) 
    9194         ! 
    92       ALLOCATE( sshb(jpi,jpj)     , sshn(jpi,jpj)  , ssha(jpi,jpj)   ,     & 
     95      ALLOCATE( sshb  (jpi,jpj)  , sshn   (jpi,jpj) , ssha(jpi,jpj)   ,     & 
    9396         &      ub_b(jpi,jpj)     , un_b(jpi,jpj)   , ua_b(jpi,jpj)   ,     & 
    9497         &      vb_b(jpi,jpj)     , vn_b(jpi,jpj)   , va_b(jpi,jpj)   ,     & 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/par_oce.F90

    r9787 r10001  
    1212   IMPLICIT NONE 
    1313   PUBLIC 
     14 
     15 
     16   !! 
     17   INTEGER, PUBLIC ::   Nnn, Np1, Nm1   ! =now, before, after 
     18 
     19   INTEGER, PUBLIC ::   Nb, Nn, Na      ! before, now, after  index 
     20 
     21 
    1422 
    1523   !!---------------------------------------------------------------------- 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/step.F90

    r9939 r10001  
    179179                         va(:,:,:) = 0._wp 
    180180 
    181       IF(  lk_asminc .AND. ln_asmiau .AND. ln_dyninc )   & 
    182                &         CALL dyn_asm_inc   ( kstp )  ! apply dynamics assimilation increment 
     181      IF( l_dynasm   )   CALL dyn_asm_inc   ( kstp )  ! apply dynamics assimilation increment 
    183182      IF( ln_bdy     )   CALL bdy_dyn3d_dmp ( kstp )  ! bdy damping trends 
    184183#if defined key_agrif 
    … …  
    232231                         tsa(:,:,:,:) = 0._wp         ! set tracer trends to zero 
    233232 
    234       IF(  lk_asminc .AND. ln_asmiau .AND. & 
    235          & ln_trainc )   CALL tra_asm_inc   ( kstp )  ! apply tracer assimilation increment 
     233      IF(  l_traasm  )   CALL tra_asm_inc   ( kstp )  ! apply tracer assimilation increment 
    236234                         CALL tra_sbc       ( kstp )  ! surface boundary condition 
    237235      IF( ln_traqsr  )   CALL tra_qsr       ( kstp )  ! penetrative solar radiation qsr 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/step_oce.F90

    r9780 r10001  
    108108#if defined key_top 
    109109   USE trcstp           ! passive tracer time-stepping      (trc_stp routine) 
     110   USE trcadv           ! passive tracer time-stepping      (trc_stp routine) 
    110111#endif 
    111112   !!---------------------------------------------------------------------- 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/stpRK3.F90

    r9939 r10001  
    1 MODULE step 
     1MODULE stpRK3 
    22   !!====================================================================== 
    3    !!                       ***  MODULE step  *** 
    4    !! Time-stepping   : manager of the ocean, tracer and ice time stepping 
     3   !!                       ***  MODULE stpRK3  *** 
     4   !! RK3 Time-stepping: manager of the ocean, tracer and ice time stepping 
    55   !!====================================================================== 
    6    !! History :  OPA  !  1991-03  (G. Madec)  Original code 
    7    !!             -   !  1991-11  (G. Madec) 
    8    !!             -   !  1992-06  (M. Imbard)  add a first output record 
    9    !!             -   !  1996-04  (G. Madec)  introduction of dynspg 
    10    !!             -   !  1996-04  (M.A. Foujols)  introduction of passive tracer 
    11    !!            8.0  !  1997-06  (G. Madec)  new architecture of call 
    12    !!            8.2  !  1997-06  (G. Madec, M. Imbard, G. Roullet)  free surface 
    13    !!             -   !  1999-02  (G. Madec, N. Grima)  hpg implicit 
    14    !!             -   !  2000-07  (J-M Molines, M. Imbard)  Open Bondary Conditions 
    15    !!   NEMO     1.0  !  2002-06  (G. Madec)  free form, suppress macro-tasking 
    16    !!             -   !  2004-08  (C. Talandier) New trends organization 
    17    !!             -   !  2005-01  (C. Ethe) Add the KPP closure scheme 
    18    !!             -   !  2005-11  (G. Madec)  Reorganisation of tra and dyn calls 
    19    !!             -   !  2006-01  (L. Debreu, C. Mazauric)  Agrif implementation 
    20    !!             -   !  2006-07  (S. Masson)  restart using iom 
    21    !!            3.2  !  2009-02  (G. Madec, R. Benshila)  reintroduicing z*-coordinate 
    22    !!             -   !  2009-06  (S. Masson, G. Madec)  TKE restart compatible with key_cpl 
    23    !!            3.3  !  2010-05  (K. Mogensen, A. Weaver, M. Martin, D. Lea) Assimilation interface 
    24    !!             -   !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase + merge TRC-TRA 
    25    !!            3.4  !  2011-04  (G. Madec, C. Ethe) Merge of dtatem and dtasal 
    26    !!            3.6  !  2012-07  (J. Simeon, G. Madec. C. Ethe)  Online coarsening of outputs 
    27    !!            3.6  !  2014-04  (F. Roquet, G. Madec) New equations of state 
    28    !!            3.6  !  2014-10  (E. Clementi, P. Oddo) Add Qiao vertical mixing in case of waves 
    29    !!            3.7  !  2014-10  (G. Madec)  LDF simplication  
    30    !!             -   !  2014-12  (G. Madec) remove KPP scheme 
    31    !!             -   !  2015-11  (J. Chanut) free surface simplification (remove filtered free surface) 
    32    !!            4.0  !  2017-05  (G. Madec)  introduction of the vertical physics manager (zdfphy) 
    33    !!---------------------------------------------------------------------- 
    34  
    35    !!---------------------------------------------------------------------- 
    36    !!   stp           : NEMO system time-stepping 
     6   !! History :  5.0  !  2018-07  (G. Madec)  original code 
     7   !!---------------------------------------------------------------------- 
     8 
     9   !!---------------------------------------------------------------------- 
     10   !!   stp_RK3       : NEMO system RK3 time-stepping scheme 
     11   !!   stp_RK3_init  : initialize the RK3 scheme 
    3712   !!---------------------------------------------------------------------- 
    3813   USE step_oce       ! time stepping definition modules 
    … …  
    4318   PRIVATE 
    4419 
    45    PUBLIC   stp   ! called by nemogcm.F90 
    46  
    47    !!---------------------------------------------------------------------- 
    48    !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     20   PUBLIC   stp_RK3   ! called by nemogcm.F90 
     21    
     22   LOGICAL ::   l_1st_stg = .TRUE.     ! 1st stage only flag 
     23   LOGICAL ::   l_2nd_stg = .TRUE.     ! 2nd stage only flag 
     24   LOGICAL ::   l_3rd_stg = .TRUE.     ! 3rd stage only flag 
     25    
     26   !!---------------------------------------------------------------------- 
     27   !! NEMO/OCE 5.0 , NEMO Consortium (2018) 
    4928   !! $Id$ 
    5029   !! Software governed by the CeCILL licence     (./LICENSE) 
    … …  
    5332 
    5433#if defined key_agrif 
    55    RECURSIVE SUBROUTINE stp( ) 
     34   RECURSIVE SUBROUTINE stp_RK3( ) 
    5635      INTEGER             ::   kstp   ! ocean time-step index 
    5736#else 
    58    SUBROUTINE stp( kstp ) 
     37   SUBROUTINE stp_RK3( kstp ) 
    5938      INTEGER, INTENT(in) ::   kstp   ! ocean time-step index 
    6039#endif 
    … …  
    7554      !!              -8- Outputs and diagnostics 
    7655      !!---------------------------------------------------------------------- 
    77       INTEGER ::   ji, jj, jk   ! dummy loop indice 
    78       INTEGER ::   indic        ! error indicator if < 0 
     56      INTEGER ::   ji, jj, jk, jstg   ! dummy loop indice 
     57      INTEGER ::   indic              ! error indicator if < 0 
    7958!!gm kcall can be removed, I guess 
    8059      INTEGER ::   kcall        ! optional integer argument (dom_vvl_sf_nxt) 
    … …  
    9372      ! 
    9473      IF( ln_timing )   CALL timing_start('stp') 
    95       ! 
     74 
     75!!!======================!!! 
     76!!!   First STAGE only   !!! 
     77!!!======================!!! 
     78 
    9679      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    9780      ! update I/O and calendar  
    … …  
    11497      IF( ln_bdy     )   CALL bdy_dta ( kstp, time_offset=+1 )   ! update dynamic & tracer data at open boundaries 
    11598                         CALL sbc     ( kstp )                   ! Sea Boundary Condition (including sea-ice) 
     99      IF ( ln_diurnal )  CALL stp_diurnal( kstp )                ! diagnose cool skin 
     100      ! 
    116101 
    117102      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    118103      ! Update stochastic parameters and random T/S fluctuations 
    119104      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    120       IF( ln_sto_eos ) CALL sto_par( kstp )          ! Stochastic parameters 
    121       IF( ln_sto_eos ) CALL sto_pts( tsn  )          ! Random T/S fluctuations 
     105      IF( ln_sto_eos )   CALL sto_par( kstp )          ! Stochastic parameters 
     106      IF( ln_sto_eos )   CALL sto_pts( tsn  )          ! Random T/S fluctuations 
    122107 
    123108      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    … …  
    155140      IF( l_ldfdyn_time                    )   CALL ldf_dyn( kstp )       ! eddy viscosity coeff.  
    156141 
    157       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    158       !  Ocean dynamics : hdiv, ssh, e3, u, v, w 
    159       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    160  
    161                             CALL ssh_nxt       ( kstp )  ! after ssh (includes call to div_hor) 
    162       IF(.NOT.ln_linssh )   CALL dom_vvl_sf_nxt( kstp )  ! after vertical scale factors  
    163                             CALL wzv           ( kstp )  ! now cross-level velocity  
    164                             CALL eos    ( tsn, rhd, rhop, gdept_n(:,:,:) )  ! now in situ density for hpg computation 
     142 
     143!!!======================!!! 
     144!!!   Loop over stages   !!! 
     145!!!======================!!! 
     146 
     147      DO jstg = 1, 3 
     148 
     149         SELECT CASE( jstg ) 
     150         CASE( 1 )   ;   rDt = rn_Dt / 3._wp 
     151         CASE( 2 )   ;   rDt = rn_Dt / 2._wp 
     152         CASE( 3 )   ;   rDt = rn_Dt 
     153         END SELECT 
     154          
     155         !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     156         !  Ocean dynamics : hdiv, ssh, e3, u, v, w 
     157         !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     158 
     159                               CALL ssh_nxt       ( kstp )  ! after ssh (includes call to div_hor) 
     160         IF(.NOT.ln_linssh )   CALL dom_vvl_sf_nxt( kstp )  ! after vertical scale factors  
     161                               CALL wzv           ( kstp )  ! now cross-level velocity  
     162                               CALL eos    ( tsb, rhd, rhop, gdept_n(:,:,:) )  ! now in situ density for hpg computation 
    165163                             
    166 !!jc: fs simplification 
    167 !!jc: lines below are useless if ln_linssh=F. Keep them here (which maintains a bug if ln_linssh=T and ln_zps=T, cf ticket #1636)  
    168 !!                                         but ensures reproductible results 
    169 !!                                         with previous versions using split-explicit free surface           
    170             IF( ln_zps .AND. .NOT. ln_isfcav )                               & 
    171                &            CALL zps_hde    ( kstp, jpts, tsn, gtsu, gtsv,   &  ! Partial steps: before horizontal gradient 
    172                &                                          rhd, gru , grv     )  ! of t, s, rd at the last ocean level 
    173             IF( ln_zps .AND.       ln_isfcav )                                          & 
    174                &            CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv, gtui, gtvi,  &  ! Partial steps for top cell (ISF) 
    175                &                                          rhd, gru , grv , grui, grvi   )  ! of t, s, rd at the first ocean level 
    176 !!jc: fs simplification 
    177                              
    178                          ua(:,:,:) = 0._wp            ! set dynamics trends to zero 
    179                          va(:,:,:) = 0._wp 
    180  
    181       IF(  lk_asminc .AND. ln_asmiau .AND. ln_dyninc )   & 
    182                &         CALL dyn_asm_inc   ( kstp )  ! apply dynamics assimilation increment 
    183       IF( ln_bdy     )   CALL bdy_dyn3d_dmp ( kstp )  ! bdy damping trends 
    184 #if defined key_agrif 
    185       IF(.NOT. Agrif_Root())  &  
    186                &         CALL Agrif_Sponge_dyn        ! momentum sponge 
    187 #endif 
    188                          CALL dyn_adv       ( kstp )  ! advection (vector or flux form) 
    189                          CALL dyn_vor       ( kstp )  ! vorticity term including Coriolis 
    190                          CALL dyn_ldf       ( kstp )  ! lateral mixing 
    191       IF( ln_zdfosm  )   CALL dyn_osm       ( kstp )  ! OSMOSIS non-local velocity fluxes 
    192                          CALL dyn_hpg       ( kstp )  ! horizontal gradient of Hydrostatic pressure 
    193                          CALL dyn_spg       ( kstp )  ! surface pressure gradient 
    194  
    195                                                       ! With split-explicit free surface, since now transports have been updated and ssha as well 
    196       IF( ln_dynspg_ts ) THEN                         ! vertical scale factors and vertical velocity need to be updated 
    197                             CALL div_hor    ( kstp )              ! Horizontal divergence  (2nd call in time-split case) 
    198          IF(.NOT.ln_linssh) CALL dom_vvl_sf_nxt( kstp, kcall=2 )  ! after vertical scale factors (update depth average component) 
    199                             CALL wzv        ( kstp )              ! now cross-level velocity  
    200       ENDIF 
    201        
    202                          CALL dyn_zdf       ( kstp )  ! vertical diffusion 
    203  
    204       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    205       ! cool skin 
    206       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<       
    207       IF ( ln_diurnal )  CALL stp_diurnal( kstp ) 
    208        
    209       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    210       ! diagnostics and outputs 
    211       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    212       IF( lk_floats  )   CALL flo_stp ( kstp )        ! drifting Floats 
    213       IF( ln_diacfl  )   CALL dia_cfl ( kstp )        ! Courant number diagnostics 
    214       IF( lk_diahth  )   CALL dia_hth ( kstp )        ! Thermocline depth (20 degres isotherm depth) 
    215       IF( lk_diadct  )   CALL dia_dct ( kstp )        ! Transports 
    216                          CALL dia_ar5 ( kstp )        ! ar5 diag 
    217       IF( lk_diaharm )   CALL dia_harm( kstp )        ! Tidal harmonic analysis 
    218                          CALL dia_wri ( kstp )        ! ocean model: outputs 
    219       ! 
    220       IF( ln_crs     )   CALL crs_fld       ( kstp )  ! ocean model: online field coarsening & output 
    221        
     164         !!jc: fs simplification 
     165         !!jc: lines below are useless if ln_linssh=F. Keep them here (which maintains a bug if ln_linssh=T and ln_zps=T, cf ticket #1636)  
     166         !!                                         but ensures reproductible results 
     167         !!                                         with previous versions using split-explicit free surface           
     168         IF( ln_zps .AND. .NOT. ln_isfcav )                               & 
     169            &            CALL zps_hde    ( kstp, jpts, tsn, gtsu, gtsv,   &  ! Partial steps: before horizontal gradient 
     170            &                                          rhd, gru , grv     )  ! of t, s, rd at the last ocean level 
     171         IF( ln_zps .AND.       ln_isfcav )                                          & 
     172            &            CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv, gtui, gtvi,  &  ! Partial steps for top cell (ISF) 
     173            &                                          rhd, gru , grv , grui, grvi   )  ! of t, s, rd at the first ocean level 
     174         !!jc: fs simplification 
     175 
     176                        ua (:,:,:)  = 0._wp          ! set the RHS of dyn Eq. to zero 
     177                        va (:,:,:)  = 0._wp 
     178                        tsa(:,:,:,:) = 0._wp         ! set tracer trends to zero 
     179         ! 
     180         !                          ! ================ !     
     181         IF ( jstg <= 3 ) THEN      !   stages 1 & 2   :   ADV, COR, HPG and SPG trends only 
     182            !                       ! ================ ! 
     183            ! 
     184            !                             !==  dynamics  ==! 
     185                        CALL dyn_adv( kstp )    ! advection (vector or flux form) 
     186                        CALL dyn_vor( kstp )    ! vorticity term including Coriolis 
     187                        CALL dyn_hpg( kstp )    ! horizontal gradient of Hydrostatic pressure 
     188                        CALL dyn_spg( kstp )    ! surface pressure gradient 
     189                                                ! With split-explicit free surface, since now transports have been updated and ssha as well 
     190            IF( ln_dynspg_ts ) THEN             ! vertical scale factors and vertical velocity need to be updated 
     191                        CALL div_hor( kstp )          ! Horizontal divergence  (2nd call in time-split case) 
     192               IF(.NOT.ln_linssh) CALL dom_vvl_sf_nxt( kstp, kcall=2 )  ! after vertical scale factors (update depth average component) 
     193                        CALL wzv    ( kstp )    ! now cross-level velocity  
     194         ENDIF 
     195!!gm  to be added here :  time stepping ==>>>   un & vn 
     196 
     197 
     198            !                             !==  tracers  ==! 
    222199#if defined key_top 
    223       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    224       ! Passive Tracer Model 
    225       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     200                        CALL trc_adv( kstp )    ! horizontal & vertical advection 
     201#endif 
     202                        CALL tra_adv( kstp )    ! horizontal & vertical advection 
     203 
     204!!gm  to be added here :  time stepping ==>>>   tsn 
     205                         
     206            ! 
     207            !                       ! ================ ! 
     208         ELSE                       !     stage 3      :   add all dynamical trends  
     209            !                       ! ================ ! 
     210            ! 
     211            !                             !==  dynamics  ==! 
     212                        CALL dyn_adv( kstp )  ! advection (vector or flux form) 
     213                        CALL dyn_vor( kstp )  ! vorticity term including Coriolis 
     214                        CALL dyn_hpg( kstp )  ! horizontal gradient of Hydrostatic pressure 
     215            ! 
     216            IF(l_dynasm)   CALL dyn_asm_inc   ( kstp )  ! apply dynamics assimilation increment 
     217            IF( ln_bdy )   CALL bdy_dyn3d_dmp ( kstp )  ! bdy damping trends 
     218#if defined key_agrif 
     219            IF(.NOT. Agrif_Root())   CALL Agrif_Sponge_dyn        ! momentum sponge 
     220#endif 
     221            IF( ln_zdfosm )   CALL dyn_osm( kstp )  ! OSMOSIS non-local velocity fluxes 
     222                        CALL dyn_ldf( kstp )  ! lateral mixing 
     223                        CALL dyn_spg( kstp )    ! surface pressure gradient 
     224                                                ! With split-explicit free surface, since now transports have been updated and ssha as well 
     225            IF( ln_dynspg_ts ) THEN             ! vertical scale factors and vertical velocity need to be updated 
     226                        CALL div_hor( kstp )          ! Horizontal divergence  (2nd call in time-split case) 
     227               IF(.NOT.ln_linssh) CALL dom_vvl_sf_nxt( kstp, kcall=2 )  ! after vertical scale factors (update depth average component) 
     228                        CALL wzv    ( kstp )    ! now cross-level velocity  
     229            ENDIF 
     230                        CALL dyn_zdf( kstp )    ! vertical diffusion  & time-stepping 
     231 
     232            !                             !==  tracers  ==! 
     233            IF( ln_crs )   CALL crs_fld( kstp )   ! ocean model: online field coarsening & output 
     234#if defined key_top 
    226235                         CALL trc_stp       ( kstp )  ! time-stepping 
    227236#endif 
    228237 
    229       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    230       ! Active tracers                               
    231       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    232238                         tsa(:,:,:,:) = 0._wp         ! set tracer trends to zero 
    233239 
    234       IF(  lk_asminc .AND. ln_asmiau .AND. & 
    235          & ln_trainc )   CALL tra_asm_inc   ( kstp )  ! apply tracer assimilation increment 
     240                         CALL tra_adv       ( kstp )  ! horizontal & vertical advection 
     241      IF(  l_traasm )   CALL tra_asm_inc   ( kstp )  ! apply tracer assimilation increment 
    236242                         CALL tra_sbc       ( kstp )  ! surface boundary condition 
    237243      IF( ln_traqsr  )   CALL tra_qsr       ( kstp )  ! penetrative solar radiation qsr 
    … …  
    241247      IF( ln_bdy     )   CALL bdy_tra_dmp   ( kstp )  ! bdy damping trends 
    242248#if defined key_agrif 
    243       IF(.NOT. Agrif_Root())  &  
    244                &         CALL Agrif_Sponge_tra        ! tracers sponge 
    245 #endif 
    246                          CALL tra_adv       ( kstp )  ! horizontal & vertical advection 
     249      IF(.NOT. Agrif_Root())   CALL Agrif_Sponge_tra        ! tracers sponge 
     250#endif 
    247251      IF( ln_zdfosm  )   CALL tra_osm       ( kstp )  ! OSMOSIS non-local tracer fluxes 
    248252      IF( lrst_oce .AND. ln_zdfosm ) & 
    … …  
    255259                         CALL tra_zdf       ( kstp )  ! vertical mixing and after tracer fields 
    256260      IF( ln_zdfnpc  )   CALL tra_npc       ( kstp )  ! update after fields by non-penetrative convection 
     261 
     262          
     263         ENDIF 
     264 
     265 
     266 
     267 
    257268 
    258269      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    … …  
    278289      IF(.NOT.ln_linssh) CALL dom_vvl_sf_swp( kstp )  ! swap of vertical scale factors 
    279290      ! 
     291 
     292 
     293 
     294 
     295!!!==========================!!! 
     296!!!   end Loop over stages   !!! 
     297!!!==========================!!! 
     298 
     299      END DO 
     300 
     301 
    280302      IF( ln_diahsb  )   CALL dia_hsb       ( kstp )  ! - ML - global conservation diagnostics 
    281303 
    … …  
    296318                         IF( Agrif_NbStepint() == 0 ) CALL Agrif_update_all( ) ! Update all components 
    297319#endif 
    298       IF( ln_diaobs  )   CALL dia_obs      ( kstp )      ! obs-minus-model (assimilation) diagnostics (call after dynamics update) 
     320 
     321      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     322      ! diagnostics and outputs 
     323      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     324      IF( ln_diaobs  )   CALL dia_obs ( kstp )        ! obs-minus-model (assimilation) diagnostics (call after dynamics update) 
     325      IF( lk_floats  )   CALL flo_stp ( kstp )        ! drifting Floats 
     326      IF( ln_diacfl  )   CALL dia_cfl ( kstp )        ! Courant number diagnostics 
     327      IF( lk_diahth  )   CALL dia_hth ( kstp )        ! Thermocline depth (20 degres isotherm depth) 
     328      IF( lk_diadct  )   CALL dia_dct ( kstp )        ! Transports 
     329                         CALL dia_ar5 ( kstp )        ! ar5 diag 
     330      IF( lk_diaharm )   CALL dia_harm( kstp )        ! Tidal harmonic analysis 
     331                         CALL dia_wri ( kstp )        ! ocean model: outputs 
    299332 
    300333      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    … …  
    323356#endif 
    324357      ! 
    325       IF( l_1st_euler ) THEN 
    326          rDt    = 2._wp * rn_Dt                             ! recover Leap-frog time-step 
    327          r1_Dt = 1._wp / rDt 
    328          l_1st_euler = .FALSE. 
    329       ENDIF 
    330       ! 
    331358      IF( ln_timing )   CALL timing_stop('stp') 
    332359      ! 
    333    END SUBROUTINE stp 
     360   END SUBROUTINE stp_RK3 
     361 
     362 
     363   SUBROUTINE stp_RK3_init 
     364      !!---------------------------------------------------------------------- 
     365      !!                     ***  ROUTINE stp_RK3_init  *** 
     366      !! 
     367      !! ** Purpose :   RK3 time stepping initialization 
     368      !! 
     369      !! ** Method  :  
     370      !!---------------------------------------------------------------------- 
     371       
     372       
     373   END SUBROUTINE stp_RK3_init 
    334374    
    335375   !!====================================================================== 
    336 END MODULE step 
     376END MODULE stpRK3 
Note: See TracChangeset for help on using the changeset viewer.