Changeset 9939


Ignore:
Timestamp:
2018-07-13T09:28:50+02:00 (2 years ago)
Author:
gm
Message:

#1911 (ENHANCE-04): RK3 branche phased with MLF@9937 branche

Location:
NEMO/branches/2018/dev_r9838_ENHANCE04_RK3
Files:
1 deleted
154 edited
2 copied

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/cfgs/AGRIF_DEMO/EXPREF/1_namelist_cfg

    r9773 r9939  
    3030   ln_linssh   = .false.   !  =T  linear free surface  ==>>  model level are fixed in time 
    3131   ! 
    32    rn_rdt      = 5760.     !  time step for the dynamics and tracer 
     32   rn_Dt       = 5760.     !  time step for the dynamics and tracer 
    3333/ 
    3434!----------------------------------------------------------------------- 
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/cfgs/AGRIF_DEMO/EXPREF/2_namelist_cfg

    r9773 r9939  
    3131   ln_linssh   = .false.   !  =T  linear free surface  ==>>  model level are fixed in time 
    3232   ! 
    33    rn_rdt      = 1440.     !  time step for the dynamics (and tracer if nn_acc=0) 
     33   rn_Dt       = 1440.     !  time step for the dynamics (and tracer if nn_acc=0) 
    3434/ 
    3535!----------------------------------------------------------------------- 
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/cfgs/AGRIF_DEMO/EXPREF/3_namelist_cfg

    r9773 r9939  
    3131   ln_linssh   = .false.   !  =T  linear free surface  ==>>  model level are fixed in time 
    3232   ! 
    33    rn_rdt      = 480.     !  time step for the dynamics (and tracer if nn_acc=0) 
     33   rn_Dt       = 480.     !  time step for the dynamics (and tracer if nn_acc=0) 
    3434/ 
    3535!----------------------------------------------------------------------- 
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/cfgs/AGRIF_DEMO/EXPREF/namelist_cfg

    r9771 r9939  
    3030   ln_linssh   = .false.   !  =T  linear free surface  ==>>  model level are fixed in time 
    3131   ! 
    32    rn_rdt      = 5760.     !  time step for the dynamics and tracer 
     32   rn_Dt       = 5760.     !  time step for the dynamics and tracer 
    3333/ 
    3434!----------------------------------------------------------------------- 
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/cfgs/AMM12/EXPREF/namelist_cfg

    r9742 r9939  
    3333&namdom        !   time and space domain 
    3434!----------------------------------------------------------------------- 
    35    rn_rdt      =   600.    !  time step for the dynamics (and tracer if nn_acc=0) 
     35   rn_Dt       =   600.    !  time step for the dynamics (and tracer if nn_acc=0) 
    3636/ 
    3737!----------------------------------------------------------------------- 
     
    301301   ln_dynspg_ts = .true.   ! split-explicit free surface 
    302302   ln_bt_auto   = .false.  ! Number of sub-step defined from: 
    303    nn_baro      = 30       ! =F : the number of sub-step in rn_rdt seconds 
     303   nn_e         = 30       ! =F : the number of sub-step in rn_Dt seconds 
    304304/ 
    305305!----------------------------------------------------------------------- 
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/cfgs/C1D_PAPA/EXPREF/namelist_cfg

    r9799 r9939  
    4949&namdom        !   time and space domain 
    5050!----------------------------------------------------------------------- 
    51    rn_rdt      =  360.     !  time step for the dynamics and tracer 
     51   rn_Dt       =  360.     !  time step for the dynamics and tracer 
    5252/ 
    5353!----------------------------------------------------------------------- 
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/cfgs/GYRE_BFM/EXPREF/namelist_cfg

    r9560 r9939  
    4545   ln_linssh   = .true.    !  =T  linear free surface  ==>>  model level are fixed in time 
    4646   ! 
    47    rn_rdt      = 7200.     !  time step for the dynamics 
     47   rn_Dt       = 7200.     !  time step for the dynamics 
    4848/ 
    4949 
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/cfgs/GYRE_PISCES/EXPREF/namelist_cfg

    r9742 r9939  
    4545   ln_linssh   = .true.    !  =T  linear free surface  ==>>  model level are fixed in time 
    4646   ! 
    47    rn_rdt      = 7200.     !  time step for the dynamics 
     47   rn_Dt       = 7200.     !  time step for the dynamics 
    4848/ 
    4949!!====================================================================== 
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/cfgs/GYRE_PISCES/cpp_GYRE_PISCES.fcm

    r9139 r9939  
    1 bld::tool::fppkeys   key_top key_mpp_mpi key_iomput 
     1bld::tool::fppkeys   key_top     key_mpp_mpi key_iomput    key_nosignedzero 
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_cfg

    r9838 r9939  
    2828&namdom        !   time and space domain 
    2929!----------------------------------------------------------------------- 
    30    rn_rdt      = 5760.     !  time step for the dynamics and tracer 
     30   rn_Dt      = 5760.     !  time step for the dynamics and tracer 
    3131/ 
    3232!----------------------------------------------------------------------- 
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/cfgs/ORCA2_OFF_PISCES/EXPREF/namelist_cfg

    r9751 r9939  
    3434   ln_linssh   = .true.   !  =T  linear free surface  ==>>  model level are fixed in time 
    3535   ! 
    36    rn_rdt      = 21600.     !  time step for the dynamics and tracer 
     36   rn_Dt       = 21600.     !  time step for the dynamics and tracer 
    3737/ 
    3838!----------------------------------------------------------------------- 
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/cfgs/ORCA2_OFF_TRC/EXPREF/namelist_cfg

    r9742 r9939  
    3434   ln_linssh   = .true.   !  =T  linear free surface  ==>>  model level are fixed in time 
    3535   ! 
    36    rn_rdt      = 21600.     !  time step for the dynamics and tracer 
     36   rn_Dt       = 21600.     !  time step for the dynamics and tracer 
    3737/ 
    3838!----------------------------------------------------------------------- 
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/cfgs/SHARED/namelist_ref

    r9838 r9939  
    4242   nn_leapy    =       0   !  Leap year calendar (1) or not (0) 
    4343   ln_rstart   = .false.   !  start from rest (F) or from a restart file (T) 
    44       nn_euler    =    1      !  = 0 : start with forward time step if ln_rstart=T 
    45       nn_rstctl   =    0      !  restart control ==> activated only if ln_rstart=T 
     44      ln_1st_euler = .false.  !  =T force a start with forward time step (ln_rstart=T) 
     45      nn_rstctl    =    0     !  restart control ==> activated only if ln_rstart=T 
    4646      !                          !    = 0 nn_date0 read in namelist ; nn_it000 : read in namelist 
    4747      !                          !    = 1 nn_date0 read in namelist ; nn_it000 : check consistancy between namelist and restart 
     
    7070   rn_isfhmin  =    1.00   !  treshold [m] to discriminate grounding ice from floating ice 
    7171   ! 
    72    rn_rdt      = 5760.     !  time step for the dynamics and tracer 
     72   rn_Dt       = 5760.     !  time step for the dynamics and tracer 
    7373   rn_atfp     =    0.1    !  asselin time filter parameter 
    7474   ! 
     
    562562!----------------------------------------------------------------------- 
    563563   ln_tide     = .false.      ! Activate tides 
    564       ln_tide_pot   = .true.                !  use tidal potential forcing 
    565          ln_scal_load  = .false.               ! Use scalar approximation for 
    566             rn_scal_load = 0.094               !     load potential 
    567          ln_read_load  = .false.               ! Or read load potential from file 
    568             cn_tide_load = 'tide_LOAD_grid_T.nc'  ! filename for load potential 
     564      ln_tide_pot   = .true.        !  use tidal potential forcing 
     565         ln_scal_load  = .false.       ! Use scalar approximation for 
     566            rn_load    = 0.094         !     load potential 
     567         ln_read_load  = .false.    ! read load potential from file 
     568            cn_tide_load = 'tide_LOAD_grid_T.nc'  ! load potential filename  
    569569            !       
    570570      ln_tide_ramp  = .false.               !  Use linear ramp for tides at startup 
    571          rdttideramp   =    0.                 !  ramp duration in days 
     571         rn_ramp    =    0.                 !  ramp duration in days 
    572572      clname(1)     = 'DUMMY'               !  name of constituent - all tidal components must be set in namelist_cfg 
    573573/ 
     
    888888   ln_dynspg_exp  = .false.   ! explicit free surface 
    889889   ln_dynspg_ts   = .false.   ! split-explicit free surface 
    890       ln_bt_fw      = .true.     ! Forward integration of barotropic Eqs. 
     890      ln_bt_fw      = .true.     ! Forward integration of external mode Eqs. 
    891891      ln_bt_av      = .true.     ! Time filtering of barotropic variables 
    892892         nn_bt_flt     = 1          ! Time filter choice  = 0 None 
     
    894894         !                          !                     = 2 Boxcar over 2*nn_baro  "    " 
    895895      ln_bt_auto    = .true.     ! Number of sub-step defined from: 
    896          rn_bt_cmax   =  0.8        ! =T : the Maximum Courant Number allowed 
    897          nn_baro      = 30          ! =F : the number of sub-step in rn_rdt seconds 
     896         rn_bt_cmax    = 0.8         ! =T : the Maximum Courant Number allowed 
     897         nn_e          = 30          ! =F : the number of sub-step in rn_Dt seconds 
    898898      rn_bt_alpha   = 0.         ! Temporal diffusion parameter (if ln_bt_av=F) 
    899899/ 
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/cfgs/SPITZ12/EXPREF/namelist_cfg

    r9793 r9939  
    2727&namdom        !   time and space domain 
    2828!----------------------------------------------------------------------- 
    29    rn_rdt      =  720.     !  time step for the dynamics and tracer 
     29   rn_Dt       =  720.     !  time step for the dynamics and tracer 
    3030   ! 
    3131/ 
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/ICE/ice.F90

    r9604 r9939  
    188188   !                                     !!** some other parameters  
    189189   INTEGER , PUBLIC ::   kt_ice           !: iteration number 
    190    REAL(wp), PUBLIC ::   rdt_ice          !: ice time step 
    191    REAL(wp), PUBLIC ::   r1_rdtice        !: = 1. / rdt_ice 
     190   REAL(wp), PUBLIC ::   rDt_ice          !: ice time step 
     191   REAL(wp), PUBLIC ::   r1_Dt_ice        !: = 1. / rdt_ice 
    192192   REAL(wp), PUBLIC ::   r1_nlay_i        !: 1 / nlay_i 
    193193   REAL(wp), PUBLIC ::   r1_nlay_s        !: 1 / nlay_s  
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/ICE/icealb.F90

    r9604 r9939  
    148148               ! 
    149149               !                       !--- Snow-covered ice albedo (freezing, melting cases) 
    150                IF( pt_su(ji,jj,jl) < rt0_snow ) THEN 
     150               IF( pt_su(ji,jj,jl) < rt0 ) THEN 
    151151                  zalb_snw = rn_alb_sdry - ( rn_alb_sdry - zalb_ice ) * EXP( - ph_snw(ji,jj,jl) * z1_c3 ) 
    152152               ELSE 
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/ICE/icecor.F90

    r9604 r9939  
    8686      IF ( nn_icesal == 2 ) THEN    !  salinity must stay in bounds [Simin,Simax]        ! 
    8787      !                             !----------------------------------------------------- 
    88          zzc = rhoic * r1_rdtice 
     88         zzc = rhoi * r1_Dt_ice 
    8989         DO jl = 1, jpl 
    9090            DO jj = 1, jpj  
     
    137137               !                 ! heat content variation (W.m-2) 
    138138               diag_heat(ji,jj) = - (  SUM( e_i(ji,jj,1:nlay_i,:) - e_i_b(ji,jj,1:nlay_i,:) )    &  
    139                   &                  + SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) )  ) * r1_rdtice 
     139                  &                  + SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) )  ) * r1_Dt_ice 
    140140               !                 ! salt, volume 
    141                diag_sice(ji,jj) = SUM( sv_i(ji,jj,:) - sv_i_b(ji,jj,:) ) * rhoic * r1_rdtice 
    142                diag_vice(ji,jj) = SUM( v_i (ji,jj,:) - v_i_b (ji,jj,:) ) * rhoic * r1_rdtice 
    143                diag_vsnw(ji,jj) = SUM( v_s (ji,jj,:) - v_s_b (ji,jj,:) ) * rhosn * r1_rdtice 
     141               diag_sice(ji,jj) = SUM( sv_i(ji,jj,:) - sv_i_b(ji,jj,:) ) * rhoi * r1_Dt_ice 
     142               diag_vice(ji,jj) = SUM( v_i (ji,jj,:) - v_i_b (ji,jj,:) ) * rhoi * r1_Dt_ice 
     143               diag_vsnw(ji,jj) = SUM( v_s (ji,jj,:) - v_s_b (ji,jj,:) ) * rhos * r1_Dt_ice 
    144144            END DO 
    145145         END DO 
    146146         !                       ! concentration tendency (dynamics) 
    147          zafx   (:,:) = SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_rdtice  
     147         zafx   (:,:) = SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_Dt_ice  
    148148         afx_tot(:,:) = zafx(:,:) 
    149149         IF( iom_use('afxdyn') )   CALL iom_put( 'afxdyn' , zafx(:,:) ) 
     
    158158               !                 ! heat content variation (W.m-2) 
    159159               diag_heat(ji,jj) = diag_heat(ji,jj) - (  SUM( e_i(ji,jj,1:nlay_i,:) - e_i_b(ji,jj,1:nlay_i,:) )    &  
    160                   &                                   + SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) )  ) * r1_rdtice 
     160                  &                                   + SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) )  ) * r1_Dt_ice 
    161161               !                 ! salt, volume 
    162                diag_sice(ji,jj) = diag_sice(ji,jj) + SUM( sv_i(ji,jj,:) - sv_i_b(ji,jj,:) ) * rhoic * r1_rdtice 
    163                diag_vice(ji,jj) = diag_vice(ji,jj) + SUM( v_i (ji,jj,:) - v_i_b (ji,jj,:) ) * rhoic * r1_rdtice 
    164                diag_vsnw(ji,jj) = diag_vsnw(ji,jj) + SUM( v_s (ji,jj,:) - v_s_b (ji,jj,:) ) * rhosn * r1_rdtice 
     162               diag_sice(ji,jj) = diag_sice(ji,jj) + SUM( sv_i(ji,jj,:) - sv_i_b(ji,jj,:) ) * rhoi * r1_Dt_ice 
     163               diag_vice(ji,jj) = diag_vice(ji,jj) + SUM( v_i (ji,jj,:) - v_i_b (ji,jj,:) ) * rhoi * r1_Dt_ice 
     164               diag_vsnw(ji,jj) = diag_vsnw(ji,jj) + SUM( v_s (ji,jj,:) - v_s_b (ji,jj,:) ) * rhos * r1_Dt_ice 
    165165            END DO 
    166166         END DO 
    167167         !                       ! concentration tendency (total + thermo) 
    168          zafx   (:,:) = SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_rdtice 
     168         zafx   (:,:) = SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_Dt_ice 
    169169         afx_tot(:,:) = afx_tot(:,:) + zafx(:,:) 
    170170         IF( iom_use('afxthd') )   CALL iom_put( 'afxthd' , zafx(:,:) ) 
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/ICE/icectl.F90

    r9604 r9939  
    9393            &                  ) *  e1e2t(:,:) ) * zconv 
    9494 
    95          pdiag_v = glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) * e1e2t * zconv ) 
    96  
    97          pdiag_s = glob_sum( SUM( sv_i * rhoic             , dim=3 ) * e1e2t * zconv ) 
     95         pdiag_v = glob_sum( SUM( v_i  * rhoi + v_s * rhos, dim=3 ) * e1e2t * zconv ) 
     96 
     97         pdiag_s = glob_sum( SUM( sv_i * rhoi             , dim=3 ) * e1e2t * zconv ) 
    9898 
    9999         pdiag_t = glob_sum( (  SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 )     & 
     
    120120  
    121121         ! outputs 
    122          zv = ( ( glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) * e1e2t ) * zconv  & 
    123             &     - pdiag_v ) * r1_rdtice - zfv ) * rday 
    124  
    125          zs = ( ( glob_sum( SUM( sv_i * rhoic             , dim=3 ) * e1e2t ) * zconv  & 
    126             &     - pdiag_s ) * r1_rdtice + zfs ) * rday 
     122         zv = ( ( glob_sum( SUM( v_i * rhoi + v_s * rhos, dim=3 ) * e1e2t ) * zconv  & 
     123            &     - pdiag_v ) * r1_Dt_ice - zfv ) * rday 
     124 
     125         zs = ( ( glob_sum( SUM( sv_i * rhoi             , dim=3 ) * e1e2t ) * zconv  & 
     126            &     - pdiag_s ) * r1_Dt_ice + zfs ) * rday 
    127127 
    128128         zt = ( glob_sum( (  SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 )   & 
    129129            &              + SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) ) * e1e2t ) * zconv   & 
    130             &   - pdiag_t ) * r1_rdtice + zft 
     130            &   - pdiag_t ) * r1_Dt_ice + zft 
    131131 
    132132         ! zvtrp and zetrp must be close to 0 if the advection scheme is conservative 
    133          zvtrp = glob_sum( ( diag_trp_vi * rhoic + diag_trp_vs * rhosn ) * e1e2t  ) * zconv * rday  
    134          zetrp = glob_sum( ( diag_trp_ei         + diag_trp_es         ) * e1e2t  ) * zconv 
     133         zvtrp = glob_sum( ( diag_trp_vi * rhoi + diag_trp_vs * rhos ) * e1e2t  ) * zconv * rday  
     134         zetrp = glob_sum( ( diag_trp_ei        + diag_trp_es        ) * e1e2t  ) * zconv 
    135135 
    136136         zvmin = glob_min( v_i ) 
     
    580580               WRITE(numout,*) ' hfx_res      : ', hfx_res(ji,jj) 
    581581               WRITE(numout,*) ' fhtur        : ', fhtur(ji,jj)  
    582                WRITE(numout,*) ' qlead        : ', qlead(ji,jj) * r1_rdtice 
     582               WRITE(numout,*) ' qlead        : ', qlead(ji,jj) * r1_Dt_ice 
    583583               WRITE(numout,*) 
    584584               WRITE(numout,*) ' - Salt fluxes at bottom interface ***' 
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/ICE/icedia.F90

    r9604 r9939  
    9595      ! 2 - Trends due to forcing  ! 
    9696      ! ---------------------------! 
    97       z_frc_volbot = r1_rau0 * glob_sum( - ( wfx_ice(:,:) + wfx_snw(:,:) + wfx_err_sub(:,:) ) * e1e2t(:,:) ) * 1.e-9  ! freshwater flux ice/snow-ocean  
    98       z_frc_voltop = r1_rau0 * glob_sum( - ( wfx_sub(:,:) + wfx_spr(:,:) ) * e1e2t(:,:) ) * 1.e-9                     ! freshwater flux ice/snow-atm 
    99       z_frc_sal    = r1_rau0 * glob_sum(   - sfx(:,:) * e1e2t(:,:) ) * 1.e-9                                          ! salt fluxes ice/snow-ocean 
     97      z_frc_volbot = r1_rho0 * glob_sum( - ( wfx_ice(:,:) + wfx_snw(:,:) + wfx_err_sub(:,:) ) * e1e2t(:,:) ) * 1.e-9  ! freshwater flux ice/snow-ocean  
     98      z_frc_voltop = r1_rho0 * glob_sum( - ( wfx_sub(:,:) + wfx_spr(:,:) ) * e1e2t(:,:) ) * 1.e-9                     ! freshwater flux ice/snow-atm 
     99      z_frc_sal    = r1_rho0 * glob_sum(   - sfx(:,:) * e1e2t(:,:) ) * 1.e-9                                          ! salt fluxes ice/snow-ocean 
    100100      z_frc_tembot =           glob_sum( hfx_out(:,:) * e1e2t(:,:) ) * 1.e-20                                         ! heat on top of ocean (and below ice) 
    101101      z_frc_temtop =           glob_sum( hfx_in (:,:) * e1e2t(:,:) ) * 1.e-20                                         ! heat on top of ice-coean 
     
    110110      ! 3 -  Content variations ! 
    111111      ! ----------------------- ! 
    112       zdiff_vol = r1_rau0 * glob_sum( ( rhoic*vt_i(:,:) + rhosn*vt_s(:,:) - vol_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-9 ! freshwater trend (km3)  
    113       zdiff_sal = r1_rau0 * glob_sum( ( rhoic* SUM( sv_i(:,:,:), dim=3 ) - sal_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-9 ! salt content trend (km3*pss) 
    114       zdiff_tem =           glob_sum( ( et_i(:,:) + et_s(:,:)             - tem_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-20 ! heat content trend (1.e20 J) 
     112      zdiff_vol = r1_rho0 * glob_sum( ( rhoi*vt_i(:,:) + rhos*vt_s(:,:) - vol_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-9 ! freshwater trend (km3)  
     113      zdiff_sal = r1_rho0 * glob_sum( ( rhoi* SUM( sv_i(:,:,:), dim=3 ) - sal_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-9 ! salt content trend (km3*pss) 
     114      zdiff_tem =           glob_sum( ( et_i(:,:) + et_s(:,:)           - tem_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-20 ! heat content trend (1.e20 J) 
    115115      !                               + SUM( qevap_ice * a_i_b, dim=3 )       !! clem: I think this term should not be there (but needs a check) 
    116116 
     
    125125      ! 5 - Diagnostics writing ! 
    126126      ! ----------------------- ! 
    127 !!gm I don't understand the division by the ocean surface (i.e. glob_sum( e1e2t(:,:) ) * 1.e-20 * kt*rdt ) 
    128 !!   and its multiplication bu kt ! is it really what we want ? what is this quantity ? 
    129 !!   IF it is really what we want, compute it at kt=nit000, not 3 time by time-step ! 
    130 !!   kt*rdt  : you mean rdtice ? 
    131 !!gm 
    132127      ! 
    133128      IF( iom_use('ibgvolume')    )   CALL iom_put( 'ibgvolume' , zdiff_vol     )   ! ice/snow volume  drift            (km3 equivalent ocean water)          
     
    135130      IF( iom_use('ibgheatco')    )   CALL iom_put( 'ibgheatco' , zdiff_tem     )   ! ice/snow heat content drift       (1.e20 J) 
    136131      IF( iom_use('ibgheatfx')    )   CALL iom_put( 'ibgheatfx' ,               &   ! ice/snow heat flux drift          (W/m2) 
    137          &                                                     zdiff_tem /glob_sum( e1e2t(:,:) * 1.e-20 * kt*rdt ) ) 
     132         &                                                     zdiff_tem /glob_sum( e1e2t(:,:) * 1.e-20 * kt*rn_Dt ) ) 
    138133 
    139134      IF( iom_use('ibgfrcvoltop') )   CALL iom_put( 'ibgfrcvoltop' , frc_voltop )   ! vol  forcing ice/snw-atm          (km3 equivalent ocean water)  
     
    143138      IF( iom_use('ibgfrctembot') )   CALL iom_put( 'ibgfrctembot' , frc_tembot )   ! heat on top of ocean(below ice)   (1.e20 J)    
    144139      IF( iom_use('ibgfrchfxtop') )   CALL iom_put( 'ibgfrchfxtop' ,            &   ! heat on top of ice/snw/ocean      (W/m2)  
    145          &                                                          frc_temtop / glob_sum( e1e2t(:,:) ) * 1.e-20 * kt*rdt  ) 
     140         &                                                          frc_temtop / glob_sum( e1e2t(:,:) ) * 1.e-20 * kt*rn_Dt  ) 
    146141      IF( iom_use('ibgfrchfxbot') )   CALL iom_put( 'ibgfrchfxbot' ,            &   ! heat on top of ocean(below ice)   (W/m2)  
    147          &                                                          frc_tembot / glob_sum( e1e2t(:,:) ) * 1.e-20 * kt*rdt  ) 
     142         &                                                          frc_tembot / glob_sum( e1e2t(:,:) ) * 1.e-20 * kt*rn_Dt  ) 
    148143 
    149144      IF( iom_use('ibgvol_tot' )  )   CALL iom_put( 'ibgvol_tot'  , zbg_ivol     )   ! ice volume                       (km3) 
     
    246241            frc_sal     = 0._wp                                                  
    247242            ! record initial ice volume, salt and temp 
    248             vol_loc_ini(:,:) = rhoic * vt_i(:,:) + rhosn * vt_s(:,:)  ! ice/snow volume (kg/m2) 
    249             tem_loc_ini(:,:) = et_i(:,:) + et_s(:,:)                  ! ice/snow heat content (J) 
    250             sal_loc_ini(:,:) = rhoic * SUM( sv_i(:,:,:), dim=3 )      ! ice salt content (pss*kg/m2) 
     243            vol_loc_ini(:,:) = rhoi * vt_i(:,:) + rhos * vt_s(:,:)   ! ice/snow volume (kg/m2) 
     244            tem_loc_ini(:,:) = et_i(:,:) + et_s(:,:)                 ! ice/snow heat content (J) 
     245            sal_loc_ini(:,:) = rhoi * SUM( sv_i(:,:,:), dim=3 )      ! ice salt content (pss*kg/m2) 
    251246         ENDIF 
    252247         ! 
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/ICE/icedyn_adv.F90

    r9604 r9939  
    9898      ! diagnostics 
    9999      !------------ 
    100       diag_trp_ei(:,:) = SUM(SUM( e_i (:,:,1:nlay_i,:) - e_i_b (:,:,1:nlay_i,:), dim=4 ), dim=3 ) * r1_rdtice 
    101       diag_trp_es(:,:) = SUM(SUM( e_s (:,:,1:nlay_s,:) - e_s_b (:,:,1:nlay_s,:), dim=4 ), dim=3 ) * r1_rdtice 
    102       diag_trp_sv(:,:) = SUM(     sv_i(:,:,:)          - sv_i_b(:,:,:)                  , dim=3 ) * r1_rdtice 
    103       diag_trp_vi(:,:) = SUM(     v_i (:,:,:)          - v_i_b (:,:,:)                  , dim=3 ) * r1_rdtice 
    104       diag_trp_vs(:,:) = SUM(     v_s (:,:,:)          - v_s_b (:,:,:)                  , dim=3 ) * r1_rdtice 
    105       IF( iom_use('icemtrp') )   CALL iom_put( "icemtrp" , diag_trp_vi * rhoic          )   ! ice mass transport 
    106       IF( iom_use('snwmtrp') )   CALL iom_put( "snwmtrp" , diag_trp_vs * rhosn          )   ! snw mass transport 
    107       IF( iom_use('salmtrp') )   CALL iom_put( "salmtrp" , diag_trp_sv * rhoic * 1.e-03 )   ! salt mass transport (kg/m2/s) 
     100      diag_trp_ei(:,:) = SUM(SUM( e_i (:,:,1:nlay_i,:) - e_i_b (:,:,1:nlay_i,:), dim=4 ), dim=3 ) * r1_Dt_ice 
     101      diag_trp_es(:,:) = SUM(SUM( e_s (:,:,1:nlay_s,:) - e_s_b (:,:,1:nlay_s,:), dim=4 ), dim=3 ) * r1_Dt_ice 
     102      diag_trp_sv(:,:) = SUM(     sv_i(:,:,:)          - sv_i_b(:,:,:)                  , dim=3 ) * r1_Dt_ice 
     103      diag_trp_vi(:,:) = SUM(     v_i (:,:,:)          - v_i_b (:,:,:)                  , dim=3 ) * r1_Dt_ice 
     104      diag_trp_vs(:,:) = SUM(     v_s (:,:,:)          - v_s_b (:,:,:)                  , dim=3 ) * r1_Dt_ice 
     105      IF( iom_use('icemtrp') )   CALL iom_put( "icemtrp" , diag_trp_vi * rhoi           )   ! ice mass transport 
     106      IF( iom_use('snwmtrp') )   CALL iom_put( "snwmtrp" , diag_trp_vs * rhos           )   ! snw mass transport 
     107      IF( iom_use('salmtrp') )   CALL iom_put( "salmtrp" , diag_trp_sv * rhoi * 1.e-03 )   ! salt mass transport (kg/m2/s) 
    108108      IF( iom_use('dihctrp') )   CALL iom_put( "dihctrp" , -diag_trp_ei                 )   ! advected ice heat content (W/m2) 
    109109      IF( iom_use('dshctrp') )   CALL iom_put( "dshctrp" , -diag_trp_es                 )   ! advected snw heat content (W/m2) 
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/ICE/icedyn_rdgrft.F90

    r9604 r9939  
    189189            ! divergence given by the advection scheme 
    190190            !   (which may not be equal to divu as computed from the velocity field) 
    191             zdivu_adv(ji) = ( 1._wp - ato_i_1d(ji) - SUM( a_i_2d(ji,:) ) ) * r1_rdtice 
     191            zdivu_adv(ji) = ( 1._wp - ato_i_1d(ji) - SUM( a_i_2d(ji,:) ) ) * r1_Dt_ice 
    192192            ! 
    193193            IF( zdivu_adv(ji) < 0._wp )   closing_net(ji) = MAX( closing_net(ji), -zdivu_adv(ji) )   ! make sure the closing rate is large enough 
     
    255255               ELSE 
    256256                  iterate_ridging  = 1 
    257                   zdivu_adv  (ji) = zfac * r1_rdtice 
     257                  zdivu_adv  (ji) = zfac * r1_Dt_ice 
    258258                  closing_net(ji) = MAX( 0._wp, -zdivu_adv(ji) ) 
    259259                  opning     (ji) = MAX( 0._wp,  zdivu_adv(ji) ) 
     
    460460            zfac = apartf(ji,jl) * closing_gross(ji) * rdt_ice 
    461461            IF( zfac > pa_i(ji,jl) ) THEN 
    462                closing_gross(ji) = pa_i(ji,jl) / apartf(ji,jl) * r1_rdtice 
     462               closing_gross(ji) = pa_i(ji,jl) / apartf(ji,jl) * r1_Dt_ice 
    463463            ENDIF 
    464464         END DO 
     
    472472         zfac = pato_i(ji) + ( opning(ji) - apartf(ji,0) * closing_gross(ji) ) * rdt_ice 
    473473         IF( zfac < 0._wp ) THEN           ! would lead to negative ato_i 
    474             opning(ji) = apartf(ji,0) * closing_gross(ji) - pato_i(ji) * r1_rdtice  
     474            opning(ji) = apartf(ji,0) * closing_gross(ji) - pato_i(ji) * r1_Dt_ice  
    475475         ELSEIF( zfac > zasum(ji) ) THEN   ! would lead to ato_i > asum 
    476             opning(ji) = apartf(ji,0) * closing_gross(ji) + ( zasum(ji) - pato_i(ji) ) * r1_rdtice  
     476            opning(ji) = apartf(ji,0) * closing_gross(ji) + ( zasum(ji) - pato_i(ji) ) * r1_Dt_ice  
    477477         ENDIF 
    478478      END DO 
     
    543543               ! volume and enthalpy (J/m2, >0) of seawater trapped into ridges 
    544544               vsw = v_i_2d(ji,jl1) * afrdg * rn_porordg 
    545                ersw(ji) = -rhoic * vsw * rcp * sst_1d(ji)   ! clem: if sst>0, then ersw <0 (is that possible?) 
     545               ersw(ji) = -rhoi * vsw * rcp * sst_1d(ji)   ! clem: if sst>0, then ersw <0 (is that possible?) 
    546546 
    547547               ! volume etc of ridging / rafting ice and new ridges (vi, vs, sm, oi, es, ei) 
     
    570570 
    571571               ! Ice-ocean exchanges associated with ice porosity 
    572                wfx_dyn_1d(ji) = wfx_dyn_1d(ji) - vsw * rhoic * r1_rdtice   ! increase in ice volume due to seawater frozen in voids 
    573                sfx_dyn_1d(ji) = sfx_dyn_1d(ji) - vsw * sss_1d(ji) * rhoic * r1_rdtice 
    574                hfx_dyn_1d(ji) = hfx_dyn_1d(ji) + ersw(ji) * r1_rdtice          ! > 0 [W.m-2]  
     572               wfx_dyn_1d(ji) = wfx_dyn_1d(ji) - vsw              * rhoi * r1_Dt_ice   ! increase in ice volume due to seawater frozen in voids 
     573               sfx_dyn_1d(ji) = sfx_dyn_1d(ji) - vsw * sss_1d(ji) * rhoi * r1_Dt_ice 
     574               hfx_dyn_1d(ji) = hfx_dyn_1d(ji) + ersw(ji)                * r1_Dt_ice   ! > 0 [W.m-2]  
    575575 
    576576               ! Put the snow lost by ridging into the ocean 
    577577               !  Note that esrdg > 0; the ocean must cool to melt snow. If the ocean temp = Tf already, new ice must grow. 
    578                wfx_snw_dyn_1d(ji) = wfx_snw_dyn_1d(ji) + ( rhosn * vsrdg(ji) * ( 1._wp - rn_fsnwrdg )   &   ! fresh water source for ocean 
    579                   &                                      + rhosn * vsrft(ji) * ( 1._wp - rn_fsnwrft ) ) * r1_rdtice 
     578               wfx_snw_dyn_1d(ji) = wfx_snw_dyn_1d(ji) + ( rhos * vsrdg(ji) * ( 1._wp - rn_fsnwrdg )   &   ! fresh water source for ocean 
     579                  &                                      + rhos * vsrft(ji) * ( 1._wp - rn_fsnwrft ) ) * r1_Dt_ice 
    580580 
    581581               ! Put the melt pond water into the ocean 
     
    583583               !       is no net mass flux between melt ponds and the ocean (see icethd_pnd.F90 for ex.) 
    584584               !IF ( ln_pnd_fwb ) THEN 
    585                !   wfx_pnd_1d(ji) = wfx_pnd_1d(ji) + ( rhofw * vprdg(ji) * ( 1._wp - rn_fpndrdg )   &        ! fresh water source for ocean 
    586                !      &                              + rhofw * vprft(ji) * ( 1._wp - rn_fpndrft ) ) * r1_rdtice 
     585               !   wfx_pnd_1d(ji) = wfx_pnd_1d(ji) + ( rhow * vprdg(ji) * ( 1._wp - rn_fpndrdg )   &        ! fresh water source for ocean 
     586               !      &                              + rhow * vprft(ji) * ( 1._wp - rn_fpndrft ) ) * r1_Dt_ice 
    587587               !ENDIF 
    588588 
     
    590590               IF( nn_icesal /= 2 )  THEN 
    591591                  sirdg2(ji)     = sirdg2(ji)     - vsw * ( sss_1d(ji) - s_i_1d(ji) )        ! ridge salinity = s_i 
    592                   sfx_bri_1d(ji) = sfx_bri_1d(ji) + sss_1d(ji) * vsw * rhoic * r1_rdtice  &  ! put back sss_m into the ocean 
    593                      &                            - s_i_1d(ji) * vsw * rhoic * r1_rdtice     ! and get  s_i  from the ocean  
     592                  sfx_bri_1d(ji) = sfx_bri_1d(ji) + sss_1d(ji) * vsw * rhoi * r1_Dt_ice  &  ! put back sss_m into the ocean 
     593                     &                            - s_i_1d(ji) * vsw * rhoi * r1_Dt_ice     ! and get  s_i  from the ocean  
    594594               ENDIF 
    595595 
     
    621621                  ! Put the snow lost by ridging into the ocean 
    622622                  hfx_dyn_1d(ji) = hfx_dyn_1d(ji) + ( - esrdg(ji,jk) * ( 1._wp - rn_fsnwrdg )   &                 ! heat sink for ocean (<0, W.m-2) 
    623                      &                                - esrft(ji,jk) * ( 1._wp - rn_fsnwrft ) ) * r1_rdtice 
     623                     &                                - esrft(ji,jk) * ( 1._wp - rn_fsnwrft ) ) * r1_Dt_ice 
    624624                  ! 
    625625                  ! Remove energy of new ridge to each category jl1 
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/ICE/icedyn_rhg_evp.F90

    r9660 r9939  
    114114      INTEGER ::   jter         ! local integers 
    115115      ! 
    116       REAL(wp) ::   zrhoco                                              ! rau0 * rn_cio 
     116      REAL(wp) ::   zrhoco                                              ! rho0 * rn_cio 
    117117      REAL(wp) ::   zdtevp, z1_dtevp                                    ! time step for subcycling 
    118118      REAL(wp) ::   ecc2, z1_ecc2                                       ! square of yield ellipse eccenticity 
     
    221221      ! 1) define some variables and initialize arrays 
    222222      !------------------------------------------------------------------------------! 
    223       zrhoco = rau0 * rn_cio  
     223      zrhoco = rho0 * rn_cio  
    224224 
    225225      ! ecc2: square of yield ellipse eccenticrity 
     
    271271         zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp 
    272272         ! 
    273          zpice(:,:) = ssh_m(:,:) + ( zintn * snwice_mass(:,:) + zintb * snwice_mass_b(:,:) ) * r1_rau0 
     273         zpice(:,:) = ssh_m(:,:) + ( zintn * snwice_mass(:,:) + zintb * snwice_mass_b(:,:) ) * r1_rho0 
    274274         ! 
    275275      ELSE                                    !== non-embedded sea ice: use ocean surface for slope calculation ==! 
     
    285285 
    286286            ! Ice/snow mass at U-V points 
    287             zm1 = ( rhosn * vt_s(ji  ,jj  ) + rhoic * vt_i(ji  ,jj  ) ) 
    288             zm2 = ( rhosn * vt_s(ji+1,jj  ) + rhoic * vt_i(ji+1,jj  ) ) 
    289             zm3 = ( rhosn * vt_s(ji  ,jj+1) + rhoic * vt_i(ji  ,jj+1) ) 
     287            zm1 = ( rhos * vt_s(ji  ,jj  ) + rhoi * vt_i(ji  ,jj  ) ) 
     288            zm2 = ( rhos * vt_s(ji+1,jj  ) + rhoi * vt_i(ji+1,jj  ) ) 
     289            zm3 = ( rhos * vt_s(ji  ,jj+1) + rhoi * vt_i(ji  ,jj+1) ) 
    290290            zmassU = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1) 
    291291            zmassV = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1) 
     
    799799               zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * rswitch 
    800800                
    801                zdiag_xmtrp_ice(ji,jj) = rhoic * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component 
    802                zdiag_ymtrp_ice(ji,jj) = rhoic * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) !        ''           Y-   '' 
    803                 
    804                zdiag_xmtrp_snw(ji,jj) = rhosn * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component 
    805                zdiag_ymtrp_snw(ji,jj) = rhosn * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) !          ''          Y-   '' 
    806                 
    807                zdiag_xatrp(ji,jj)     = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) )         ! area transport,      X-component 
    808                zdiag_yatrp(ji,jj)     = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) )         !        ''            Y-   '' 
     801               zdiag_xmtrp_ice(ji,jj) = rhoi * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) )  ! ice mass transport, X-component 
     802               zdiag_ymtrp_ice(ji,jj) = rhoi * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) )  !        ''           Y-   '' 
     803                
     804               zdiag_xmtrp_snw(ji,jj) = rhos * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) )  ! snow mass transport, X-component 
     805               zdiag_ymtrp_snw(ji,jj) = rhos * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) )  !          ''          Y-   '' 
     806                
     807               zdiag_xatrp(ji,jj)     = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) )          ! area transport,      X-component 
     808               zdiag_yatrp(ji,jj)     = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) )          !        ''            Y-   '' 
    809809                
    810810            END DO 
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/ICE/iceistate.F90

    r9656 r9939  
    295295                  ! In case snow load is in excess that would lead to transformation from snow to ice 
    296296                  ! Then, transfer the snow excess into the ice (different from icethd_dh) 
    297                   zdh = MAX( 0._wp, ( rhosn * h_s(ji,jj,jl) + ( rhoic - rau0 ) * h_i(ji,jj,jl) ) * r1_rau0 )  
     297                  zdh = MAX( 0._wp, ( rhos * h_s(ji,jj,jl) + ( rhoi - rho0 ) * h_i(ji,jj,jl) ) * r1_rho0 )  
    298298                  ! recompute h_i, h_s avoiding out of bounds values 
    299299                  h_i(ji,jj,jl) = MIN( hi_max(jl), h_i(ji,jj,jl) + zdh ) 
    300                   h_s(ji,jj,jl) = MAX( 0._wp, h_s(ji,jj,jl) - zdh * rhoic * r1_rhosn ) 
     300                  h_s(ji,jj,jl) = MAX( 0._wp, h_s(ji,jj,jl) - zdh * rhoi * r1_rhos ) 
    301301                  ! 
    302302                  ! ice volume, salt content, age content 
     
    321321                     t_s(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0 
    322322                     ! Snow energy of melting 
    323                      e_s(ji,jj,jk,jl) = zswitch(ji,jj) * rhosn * ( cpic * ( rt0 - t_s(ji,jj,jk,jl) ) + lfus ) 
     323                     e_s(ji,jj,jk,jl) = zswitch(ji,jj) * rhos * ( rcpi * ( rt0 - t_s(ji,jj,jk,jl) ) + rLfus ) 
    324324                     ! 
    325325                     ! Mutliply by volume, and divide by number of layers to get heat content in J/m2 
     
    340340                     ! 
    341341                     ! heat content per unit volume 
    342                      e_i(ji,jj,jk,jl) = zswitch(ji,jj) * rhoic * (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) )         & 
    343                         &             + lfus * ( 1._wp - (ztmelts-rt0) / MIN( (t_i(ji,jj,jk,jl)-rt0) , -epsi20 )  )   & 
    344                         &             - rcp  * ( ztmelts - rt0 ) ) 
     342                     e_i(ji,jj,jk,jl) = zswitch(ji,jj) * rhoi * (   rcpi    * ( ztmelts - t_i(ji,jj,jk,jl) )           & 
     343                        &             + rLfus * ( 1._wp - (ztmelts-rt0) / MIN( (t_i(ji,jj,jk,jl)-rt0) , -epsi20 )  )   & 
     344                        &             - rcp   * ( ztmelts - rt0 ) ) 
    345345                     ! 
    346346                     ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2 
     
    410410      ! 5) Snow-ice mass (case ice is fully embedded) 
    411411      !---------------------------------------------- 
    412       snwice_mass  (:,:) = tmask(:,:,1) * SUM( rhosn * v_s(:,:,:) + rhoic * v_i(:,:,:), dim=3  )   ! snow+ice mass 
     412      snwice_mass  (:,:) = tmask(:,:,1) * SUM( rhos * v_s(:,:,:) + rhoi * v_i(:,:,:), dim=3  )   ! snow+ice mass 
    413413      snwice_mass_b(:,:) = snwice_mass(:,:) 
    414414      ! 
    415415      IF( ln_ice_embd ) THEN            ! embedded sea-ice: deplete the initial ssh below sea-ice area 
    416416         ! 
    417          sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0 
    418          sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0 
     417         sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rho0 
     418         sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rho0 
    419419         ! 
    420420         IF( .NOT.ln_linssh ) THEN 
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/ICE/icestp.F90

    r9725 r9939  
    341341      IF( ln_bdy .AND. ln_icediachk )   CALL ctl_warn('par_init: online conservation check does not work with BDY') 
    342342      ! 
    343       rdt_ice   = REAL(nn_fsbc) * rdt          !--- sea-ice timestep and its inverse 
    344       r1_rdtice = 1._wp / rdt_ice 
     343      rdt_ice   = REAL(nn_fsbc) * rn_Dt        !--- sea-ice timestep and its inverse 
     344      r1_Dt_ice = 1._wp / rdt_ice 
    345345      IF(lwp) WRITE(numout,*) 
    346       IF(lwp) WRITE(numout,*) '      ice timestep rdt_ice = nn_fsbc*rdt = ', rdt_ice 
     346      IF(lwp) WRITE(numout,*) '      ice timestep rdt_ice = nn_fsbc*rn_Dt = ', rdt_ice, ' [s]' 
    347347      ! 
    348348      r1_nlay_i = 1._wp / REAL( nlay_i, wp )   !--- inverse of nlay_i and nlay_s 
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/ICE/icethd.F90

    r9750 r9939  
    120120         DO jj = 2, jpjm1 
    121121            DO ji = fs_2, fs_jpim1 
    122                zfric(ji,jj) = r1_rau0 * SQRT( 0.5_wp *  & 
     122               zfric(ji,jj) = r1_rho0 * SQRT( 0.5_wp *  & 
    123123                  &                         (  utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj)   & 
    124124                  &                          + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1) ) ) * tmask(ji,jj,1) 
     
    150150            ! --- Energy needed to bring ocean surface layer until its freezing (<0, J.m-2) --- ! 
    151151            ! includes supercooling potential energy (>0) or "above-freezing" energy (<0) 
    152             zqfr = tmask(ji,jj,1) * rau0 * rcp * e3t_m(ji,jj) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ) 
     152            zqfr = tmask(ji,jj,1) * rho0_rcp * e3t_m(ji,jj) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ) 
    153153 
    154154            ! --- Above-freezing sensible heat content (J/m2 grid) 
    155             zqfr_neg = tmask(ji,jj,1) * rau0 * rcp * e3t_m(ji,jj) * MIN( ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ), 0._wp ) 
     155            zqfr_neg = tmask(ji,jj,1) * rho0_rcp * e3t_m(ji,jj) * MIN( ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ), 0._wp ) 
    156156 
    157157            ! --- Sensible ocean-to-ice heat flux (W/m2) 
    158158            zfric_u      = MAX( SQRT( zfric(ji,jj) ), zfric_umin )  
    159             fhtur(ji,jj) = rswitch * rau0 * rcp * zch  * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) ! W.m-2 
    160  
    161             fhtur(ji,jj) = rswitch * MIN( fhtur(ji,jj), - zqfr_neg * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) ) 
     159            fhtur(ji,jj) = rswitch * rho0_rcp * zch  * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) ! W.m-2 
     160 
     161            fhtur(ji,jj) = rswitch * MIN( fhtur(ji,jj), - zqfr_neg * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) ) 
    162162            ! upper bound for fhtur: the heat retrieved from the ocean must be smaller than the heat necessary to reach  
    163163            !                        the freezing point, so that we do not have SST < T_freeze 
     
    169169            ! If there is ice and leads are warming, then transfer energy from the lead budget and use it for bottom melting  
    170170            IF( zqld > 0._wp ) THEN 
    171                fhld (ji,jj) = rswitch * zqld * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) ! divided by at_i since this is (re)multiplied by a_i in icethd_dh.F90 
     171               fhld (ji,jj) = rswitch * zqld * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) ! divided by at_i since this is (re)multiplied by a_i in icethd_dh.F90 
    172172               qlead(ji,jj) = 0._wp 
    173173            ELSE 
     
    197197      !     Third  step in iceupdate.F90  :  heat from ice-ocean mass exchange (zf_mass) + solar 
    198198      hfx_out(:,:) = ( 1._wp - at_i_b(:,:) ) * qns_oce(:,:) + qemp_oce(:,:)  &  ! Non solar heat flux received by the ocean                
    199          &           - qlead(:,:) * r1_rdtice                                &  ! heat flux taken from the ocean where there is open water ice formation 
     199         &           - qlead(:,:) * r1_Dt_ice                                &  ! heat flux taken from the ocean where there is open water ice formation 
    200200         &           - at_i (:,:) * fhtur(:,:)                               &  ! heat flux taken by turbulence 
    201201         &           - at_i (:,:) *  fhld(:,:)                                  ! heat flux taken during bottom growth/melt  
     
    295295            ztmelts       = -tmut * sz_i_1d(ji,jk) 
    296296            ! Conversion q(S,T) -> T (second order equation) 
    297             zbbb          = ( rcp - cpic ) * ztmelts + e_i_1d(ji,jk) * r1_rhoic - lfus 
    298             zccc          = SQRT( MAX( zbbb * zbbb - 4._wp * cpic * lfus * ztmelts, 0._wp ) ) 
    299             t_i_1d(ji,jk) = rt0 - ( zbbb + zccc ) * 0.5_wp * r1_cpic 
     297            zbbb          = ( rcp - rcpi ) * ztmelts + e_i_1d(ji,jk) * r1_rhoi - rLfus 
     298            zccc          = SQRT( MAX( zbbb * zbbb - 4._wp * rcpi * rLfus * ztmelts, 0._wp ) ) 
     299            t_i_1d(ji,jk) = rt0 - ( zbbb + zccc ) * 0.5_wp * r1_cpi 
    300300             
    301301            ! mask temperature 
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/ICE/icethd_da.F90

    r9604 r9939  
    137137             
    138138            ! Contribution to salt flux 
    139             sfx_lam_1d(ji) = sfx_lam_1d(ji) + rhoic *  h_i_1d(ji) * zda * s_i_1d(ji) * r1_rdtice 
     139            sfx_lam_1d(ji) = sfx_lam_1d(ji) + rhoi *  h_i_1d(ji) * zda * s_i_1d(ji) * r1_Dt_ice 
    140140             
    141141            ! Contribution to heat flux into the ocean [W.m-2], (<0)   
    142             hfx_thd_1d(ji) = hfx_thd_1d(ji) - zda * r1_rdtice * ( h_i_1d(ji) * r1_nlay_i * SUM( e_i_1d(ji,1:nlay_i) )  & 
     142            hfx_thd_1d(ji) = hfx_thd_1d(ji) - zda * r1_Dt_ice * ( h_i_1d(ji) * r1_nlay_i * SUM( e_i_1d(ji,1:nlay_i) )  & 
    143143                                                                + h_s_1d(ji) * r1_nlay_s * SUM( e_s_1d(ji,1:nlay_s) ) )  
    144144             
    145145            ! Contribution to mass flux 
    146             wfx_lam_1d(ji) =  wfx_lam_1d(ji) + zda * r1_rdtice * ( rhoic * h_i_1d(ji) + rhosn * h_s_1d(ji) ) 
     146            wfx_lam_1d(ji) =  wfx_lam_1d(ji) + zda * r1_Dt_ice * ( rhoi * h_i_1d(ji) + rhos * h_s_1d(ji) ) 
    147147             
    148148            ! new concentration 
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/ICE/icethd_dh.F90

    r9767 r9939  
    7676      REAL(wp) ::   zgrr         ! bottom growth rate 
    7777      REAL(wp) ::   zt_i_new     ! bottom formation temperature 
    78       REAL(wp) ::   z1_rho       ! 1/(rhosn+rau0-rhoic) 
     78      REAL(wp) ::   z1_rho       ! 1/(rhos+rho0-rhoi) 
    7979 
    8080      REAL(wp) ::   zQm          ! enthalpy exchanged with the ocean (J/m2), >0 towards the ocean 
     
    181181         DO ji = 1, npti 
    182182            IF( t_s_1d(ji,jk) > rt0 ) THEN 
    183                hfx_res_1d    (ji) = hfx_res_1d    (ji) + e_s_1d(ji,jk) * zh_s(ji,jk) * a_i_1d(ji) * r1_rdtice   ! heat flux to the ocean [W.m-2], < 0 
    184                wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + rhosn         * zh_s(ji,jk) * a_i_1d(ji) * r1_rdtice   ! mass flux 
     183               hfx_res_1d    (ji) = hfx_res_1d    (ji) + e_s_1d(ji,jk) * zh_s(ji,jk) * a_i_1d(ji) * r1_Dt_ice   ! heat flux to the ocean [W.m-2], < 0 
     184               wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + rhos          * zh_s(ji,jk) * a_i_1d(ji) * r1_Dt_ice   ! mass flux 
    185185               ! updates 
    186186               dh_s_mlt(ji)    = dh_s_mlt(ji) - zh_s(ji,jk) 
     
    202202            ! 
    203203            ! --- precipitation --- 
    204             zdh_s_pre (ji) = zsnw(ji) * sprecip_1d(ji) * rdt_ice * r1_rhosn / at_i_1d(ji)   ! thickness change 
     204            zdh_s_pre (ji) = zsnw(ji) * sprecip_1d(ji) * rdt_ice * r1_rhos / at_i_1d(ji)    ! thickness change 
    205205            zqprec    (ji) = - qprec_ice_1d(ji)                                             ! enthalpy of the precip (>0, J.m-3) 
    206206            ! 
    207             hfx_spr_1d(ji) = hfx_spr_1d(ji) + zdh_s_pre(ji) * a_i_1d(ji) * zqprec(ji)    * r1_rdtice   ! heat flux from snow precip (>0, W.m-2) 
    208             wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhosn         * a_i_1d(ji) * zdh_s_pre(ji) * r1_rdtice   ! mass flux, <0 
     207            hfx_spr_1d(ji) = hfx_spr_1d(ji) + zdh_s_pre(ji) * a_i_1d(ji) * zqprec(ji)    * r1_Dt_ice   ! heat flux from snow precip (>0, W.m-2) 
     208            wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhos          * a_i_1d(ji) * zdh_s_pre(ji) * r1_Dt_ice   ! mass flux, <0 
    209209             
    210210            ! --- melt of falling snow --- 
     
    212212            zdeltah       (ji,1) = - rswitch * zq_su(ji) / MAX( zqprec(ji) , epsi20 )   ! thickness change 
    213213            zdeltah       (ji,1) = MAX( - zdh_s_pre(ji), zdeltah(ji,1) )                ! bound melting  
    214             hfx_snw_1d    (ji)   = hfx_snw_1d    (ji) - zdeltah(ji,1) * a_i_1d(ji) * zqprec(ji)    * r1_rdtice   ! heat used to melt snow (W.m-2, >0) 
    215             wfx_snw_sum_1d(ji)   = wfx_snw_sum_1d(ji) - rhosn         * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice   ! snow melting only = water into the ocean (then without snow precip), >0 
     214            hfx_snw_1d    (ji)   = hfx_snw_1d    (ji) - zdeltah(ji,1) * a_i_1d(ji) * zqprec(ji)    * r1_Dt_ice   ! heat used to melt snow (W.m-2, >0) 
     215            wfx_snw_sum_1d(ji)   = wfx_snw_sum_1d(ji) - rhos          * a_i_1d(ji) * zdeltah(ji,1) * r1_Dt_ice   ! snow melting only = water into the ocean (then without snow precip), >0 
    216216             
    217217            ! updates available heat + precipitations after melting 
     
    252252               zdh_s_mel(ji)    = zdh_s_mel(ji) + zdeltah(ji,jk) 
    253253                
    254                hfx_snw_1d(ji)     = hfx_snw_1d(ji)     - zdeltah(ji,jk) * a_i_1d(ji) * e_s_1d (ji,jk) * r1_rdtice   ! heat used to melt snow(W.m-2, >0) 
    255                wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn          * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice   ! snow melting only = water into the ocean (then without snow precip) 
     254               hfx_snw_1d(ji)     = hfx_snw_1d(ji)     - zdeltah(ji,jk) * a_i_1d(ji) * e_s_1d (ji,jk) * r1_Dt_ice   ! heat used to melt snow(W.m-2, >0) 
     255               wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhos           * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice   ! snow melting only = water into the ocean (then without snow precip) 
    256256                
    257257               ! updates available heat + thickness 
     
    273273         IF( evap_ice_1d(ji) > 0._wp ) THEN 
    274274            ! 
    275             zdh_s_sub (ji)   = MAX( - h_s_1d(ji) , - evap_ice_1d(ji) * r1_rhosn * rdt_ice ) 
    276             zevap_rema(ji)   = evap_ice_1d(ji) * rdt_ice + zdh_s_sub(ji) * rhosn   ! remaining evap in kg.m-2 (used for ice melting later on) 
     275            zdh_s_sub (ji)   = MAX( - h_s_1d(ji) , - evap_ice_1d(ji) * r1_rhos  * rdt_ice ) 
     276            zevap_rema(ji)   = evap_ice_1d(ji) * rdt_ice + zdh_s_sub(ji) * rhos   ! remaining evap in kg.m-2 (used for ice melting later on) 
    277277            zdeltah   (ji,1) = MAX( zdh_s_sub(ji), - zdh_s_pre(ji) ) 
    278278             
    279279            hfx_sub_1d    (ji) = hfx_sub_1d(ji) + &   ! Heat flux by sublimation [W.m-2], < 0 (sublimate snow that had fallen, then pre-existing snow) 
    280280               &                 ( zdeltah(ji,1) * zqprec(ji) + ( zdh_s_sub(ji) - zdeltah(ji,1) ) * e_s_1d(ji,1) )  & 
    281                &                 * a_i_1d(ji) * r1_rdtice 
    282             wfx_snw_sub_1d(ji) = wfx_snw_sub_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_sub(ji) * r1_rdtice   ! Mass flux by sublimation 
     281               &                 * a_i_1d(ji) * r1_Dt_ice 
     282            wfx_snw_sub_1d(ji) = wfx_snw_sub_1d(ji) - rhos * a_i_1d(ji) * zdh_s_sub(ji) * r1_Dt_ice   ! Mass flux by sublimation 
    283283             
    284284            ! new snow thickness 
     
    309309            e_s_1d(ji,jk) = rswitch / MAX( h_s_1d(ji), epsi20 ) *            & 
    310310              &             ( ( zdh_s_pre(ji)              ) * zqprec(ji) +  & 
    311               &               ( h_s_1d(ji) - zdh_s_pre(ji) ) * rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) ) 
     311              &               ( h_s_1d(ji) - zdh_s_pre(ji) ) * rhos * ( rcpi * ( rt0 - t_s_1d(ji,jk) ) + rLfus ) ) 
    312312         END DO 
    313313      END DO 
     
    326326            IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN   !-- Internal melting 
    327327 
    328                zEi            = - e_i_1d(ji,jk) * r1_rhoic            ! Specific enthalpy of layer k [J/kg, <0]        
     328               zEi            = - e_i_1d(ji,jk) * r1_rhoi             ! Specific enthalpy of layer k [J/kg, <0]        
    329329               zdE            = 0._wp                                 ! Specific enthalpy difference (J/kg, <0) 
    330330                                                                      ! set up at 0 since no energy is needed to melt water...(it is already melted) 
    331331               zdeltah(ji,jk) = MIN( 0._wp , - zh_i(ji,jk) )          ! internal melting occurs when the internal temperature is above freezing      
    332332                                                                      ! this should normally not happen, but sometimes, heat diffusion leads to this 
    333                zfmdt          = - zdeltah(ji,jk) * rhoic              ! Mass flux x time step > 0 
     333               zfmdt          = - zdeltah(ji,jk) * rhoi               ! Mass flux x time step > 0 
    334334                          
    335335               dh_i_itm(ji)   = dh_i_itm(ji) + zdeltah(ji,jk)         ! Cumulate internal melting 
    336336                
    337                zfmdt          = - rhoic * zdeltah(ji,jk)              ! Recompute mass flux [kg/m2, >0] 
    338  
    339                hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice                           ! Heat flux to the ocean [W.m-2], <0 
     337               zfmdt          = - rhoi * zdeltah(ji,jk)               ! Recompute mass flux [kg/m2, >0] 
     338 
     339               hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_Dt_ice                           ! Heat flux to the ocean [W.m-2], <0 
    340340               !                                                                                                  ice enthalpy zEi is "sent" to the ocean 
    341                sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice   ! Salt flux 
     341               sfx_res_1d(ji) = sfx_res_1d(ji) - rhoi  * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_Dt_ice   ! Salt flux 
    342342               !                                                                                                  using s_i_1d and not sz_i_1d(jk) is ok 
    343                wfx_res_1d(ji) = wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice                ! Mass flux 
     343               wfx_res_1d(ji) = wfx_res_1d(ji) - rhoi  * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice                ! Mass flux 
    344344 
    345345            ELSE                                        !-- Surface melting 
    346346                
    347                zEi            = - e_i_1d(ji,jk) * r1_rhoic            ! Specific enthalpy of layer k [J/kg, <0] 
     347               zEi            = - e_i_1d(ji,jk) * r1_rhoi             ! Specific enthalpy of layer k [J/kg, <0] 
    348348               zEw            =    rcp * ztmelts                      ! Specific enthalpy of resulting meltwater [J/kg, <0] 
    349349               zdE            =    zEi - zEw                          ! Specific enthalpy difference < 0 
     
    351351               zfmdt          = - zq_su(ji) / zdE                     ! Mass flux to the ocean [kg/m2, >0] 
    352352                
    353                zdeltah(ji,jk) = - zfmdt * r1_rhoic                    ! Melt of layer jk [m, <0] 
     353               zdeltah(ji,jk) = - zfmdt * r1_rhoi                     ! Melt of layer jk [m, <0] 
    354354                
    355355               zdeltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk) , - zh_i(ji,jk) ) )    ! Melt of layer jk cannot exceed the layer thickness [m, <0] 
    356356                
    357                zq_su(ji)      = MAX( 0._wp , zq_su(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat 
     357               zq_su(ji)      = MAX( 0._wp , zq_su(ji) - zdeltah(ji,jk) * rhoi * zdE ) ! update available heat 
    358358                
    359359               dh_i_sum(ji)   = dh_i_sum(ji) + zdeltah(ji,jk)         ! Cumulate surface melt 
    360360                
    361                zfmdt          = - rhoic * zdeltah(ji,jk)              ! Recompute mass flux [kg/m2, >0] 
     361               zfmdt          = - rhoi * zdeltah(ji,jk)               ! Recompute mass flux [kg/m2, >0] 
    362362                
    363363               zQm            = zfmdt * zEw                           ! Energy of the melt water sent to the ocean [J/m2, <0] 
    364364                
    365                sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice   ! Salt flux >0 
     365               sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_Dt_ice   ! Salt flux >0 
    366366               !                                                                                                  using s_i_1d and not sz_i_1d(jk) is ok) 
    367                hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice                           ! Heat flux [W.m-2], < 0 
    368                hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice                           ! Heat flux used in this process [W.m-2], > 0   
     367               hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_Dt_ice                           ! Heat flux [W.m-2], < 0 
     368               hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_Dt_ice                           ! Heat flux used in this process [W.m-2], > 0   
    369369               !  
    370                wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice                ! Mass flux 
     370               wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice                ! Mass flux 
    371371                
    372372            END IF 
     
    374374            ! Ice sublimation 
    375375            ! --------------- 
    376             zdum            = MAX( - ( zh_i(ji,jk) + zdeltah(ji,jk) ) , - zevap_rema(ji) * r1_rhoic ) 
     376            zdum            = MAX( - ( zh_i(ji,jk) + zdeltah(ji,jk) ) , - zevap_rema(ji) * r1_rhoi ) 
    377377            zdeltah (ji,jk) = zdeltah (ji,jk) + zdum 
    378378            dh_i_sub(ji)    = dh_i_sub(ji)    + zdum 
    379379             
    380             sfx_sub_1d(ji)     = sfx_sub_1d(ji) - rhoic * a_i_1d(ji) * zdum * s_i_1d(ji) * r1_rdtice ! Salt flux >0 
     380            sfx_sub_1d(ji)     = sfx_sub_1d(ji) - rhoi * a_i_1d(ji) * zdum * s_i_1d(ji) * r1_Dt_ice ! Salt flux >0 
    381381            !                                                                                          clem: flux is sent to the ocean for simplicity 
    382382            !                                                                                                but salt should remain in the ice except 
    383383            !                                                                                                if all ice is melted. => must be corrected 
    384             hfx_sub_1d(ji)     = hfx_sub_1d(ji) + zdum * e_i_1d(ji,jk) * a_i_1d(ji) * r1_rdtice      ! Heat flux [W.m-2], < 0 
    385  
    386             wfx_ice_sub_1d(ji) = wfx_ice_sub_1d(ji) - rhoic * a_i_1d(ji) * zdum * r1_rdtice          ! Mass flux > 0 
     384            hfx_sub_1d(ji)     = hfx_sub_1d(ji) + zdum * e_i_1d(ji,jk) * a_i_1d(ji) * r1_Dt_ice      ! Heat flux [W.m-2], < 0 
     385 
     386            wfx_ice_sub_1d(ji) = wfx_ice_sub_1d(ji) - rhoi * a_i_1d(ji) * zdum * r1_Dt_ice          ! Mass flux > 0 
    387387 
    388388            ! update remaining mass flux 
    389             zevap_rema(ji)  = zevap_rema(ji) + zdum * rhoic 
     389            zevap_rema(ji)  = zevap_rema(ji) + zdum * rhoi 
    390390             
    391391            ! record which layers have disappeared (for bottom melting)  
     
    409409      ! remaining "potential" evap is sent to ocean 
    410410      DO ji = 1, npti 
    411          wfx_err_sub_1d(ji) = wfx_err_sub_1d(ji) - zevap_rema(ji) * a_i_1d(ji) * r1_rdtice  ! <=0 (net evap for the ocean in kg.m-2.s-1) 
     411         wfx_err_sub_1d(ji) = wfx_err_sub_1d(ji) - zevap_rema(ji) * a_i_1d(ji) * r1_Dt_ice  ! <=0 (net evap for the ocean in kg.m-2.s-1) 
    412412      END DO 
    413413 
     
    437437               !--- zswi12 if 2.0e-8 < dh/dt < 3.6e-7  
    438438               !--- zswi2  if dh/dt > 3.6e-7 
    439                zgrr     = MIN( 1.0e-3, MAX ( dh_i_bog(ji) * r1_rdtice , epsi10 ) ) 
     439               zgrr     = MIN( 1.0e-3, MAX ( dh_i_bog(ji) * r1_Dt_ice , epsi10 ) ) 
    440440               zswi2    = MAX( 0._wp , SIGN( 1._wp , zgrr - 3.6e-7 ) ) 
    441441               zswi12   = MAX( 0._wp , SIGN( 1._wp , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 ) 
     
    450450               zt_i_new      = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i) 
    451451                
    452                zEi           = cpic * ( zt_i_new - (ztmelts+rt0) ) &                                  ! Specific enthalpy of forming ice (J/kg, <0) 
    453                   &            - lfus * ( 1.0 - ztmelts / ( zt_i_new - rt0 ) ) + rcp * ztmelts 
     452               zEi           = rcpi * ( zt_i_new - (ztmelts+rt0) ) &                                  ! Specific enthalpy of forming ice (J/kg, <0) 
     453                  &          - rLfus * ( 1.0 - ztmelts / ( zt_i_new - rt0 ) ) + rcp * ztmelts 
    454454 
    455455               zEw           = rcp  * ( t_bo_1d(ji) - rt0 )                                           ! Specific enthalpy of seawater (J/kg, < 0) 
     
    457457               zdE           = zEi - zEw                                                              ! Specific enthalpy difference (J/kg, <0) 
    458458 
    459                dh_i_bog(ji)  = rdt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoic ) ) 
     459               dh_i_bog(ji)  = rdt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoi ) ) 
    460460                
    461461            END DO 
    462462            ! Contribution to Energy and Salt Fluxes                                     
    463             zfmdt          = - rhoic * dh_i_bog(ji)                                                   ! Mass flux x time step (kg/m2, < 0) 
    464              
    465             hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice                           ! Heat flux to the ocean [W.m-2], >0 
    466             hfx_bog_1d(ji) = hfx_bog_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice                           ! Heat flux used in this process [W.m-2], <0 
    467              
    468             sfx_bog_1d(ji) = sfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bog(ji) * s_i_new(ji) * r1_rdtice    ! Salt flux, <0 
    469  
    470             wfx_bog_1d(ji) = wfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bog(ji) * r1_rdtice                  ! Mass flux, <0 
     463            zfmdt          = - rhoi * dh_i_bog(ji)                                                   ! Mass flux x time step (kg/m2, < 0) 
     464             
     465            hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_Dt_ice                           ! Heat flux to the ocean [W.m-2], >0 
     466            hfx_bog_1d(ji) = hfx_bog_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_Dt_ice                           ! Heat flux used in this process [W.m-2], <0 
     467             
     468            sfx_bog_1d(ji) = sfx_bog_1d(ji) - rhoi * a_i_1d(ji) * dh_i_bog(ji) * s_i_new(ji) * r1_Dt_ice    ! Salt flux, <0 
     469 
     470            wfx_bog_1d(ji) = wfx_bog_1d(ji) - rhoi * a_i_1d(ji) * dh_i_bog(ji) * r1_Dt_ice                  ! Mass flux, <0 
    471471 
    472472            ! update heat content (J.m-2) and layer thickness 
    473             eh_i_old(ji,nlay_i+1) = eh_i_old(ji,nlay_i+1) + dh_i_bog(ji) * (-zEi * rhoic) 
     473            eh_i_old(ji,nlay_i+1) = eh_i_old(ji,nlay_i+1) + dh_i_bog(ji) * (-zEi * rhoi) 
    474474            h_i_old (ji,nlay_i+1) = h_i_old (ji,nlay_i+1) + dh_i_bog(ji) 
    475475 
     
    489489               IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN   !-- Internal melting 
    490490 
    491                   zEi               = - e_i_1d(ji,jk) * r1_rhoic    ! Specific enthalpy of melting ice (J/kg, <0) 
     491                  zEi               = - e_i_1d(ji,jk) * r1_rhoi     ! Specific enthalpy of melting ice (J/kg, <0) 
    492492                  zdE               = 0._wp                         ! Specific enthalpy difference   (J/kg, <0) 
    493493                                                                    !    set up at 0 since no energy is needed to melt water...(it is already melted) 
     
    497497                  dh_i_itm (ji)     = dh_i_itm(ji) + zdeltah(ji,jk) 
    498498 
    499                   zfmdt             = - zdeltah(ji,jk) * rhoic      ! Mass flux x time step > 0 
    500  
    501                   hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice                           ! Heat flux to the ocean [W.m-2], <0 
     499                  zfmdt             = - zdeltah(ji,jk) * rhoi       ! Mass flux x time step > 0 
     500 
     501                  hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_Dt_ice                           ! Heat flux to the ocean [W.m-2], <0 
    502502                  !                                                                                                  ice enthalpy zEi is "sent" to the ocean 
    503                   sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice   ! Salt flux 
     503                  sfx_res_1d(ji) = sfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_Dt_ice   ! Salt flux 
    504504                  !                                                                                                  using s_i_1d and not sz_i_1d(jk) is ok 
    505                   wfx_res_1d(ji) = wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice                ! Mass flux 
     505                  wfx_res_1d(ji) = wfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice                ! Mass flux 
    506506 
    507507                  ! update heat content (J.m-2) and layer thickness 
     
    511511               ELSE                                        !-- Basal melting 
    512512 
    513                   zEi             = - e_i_1d(ji,jk) * r1_rhoic                                ! Specific enthalpy of melting ice (J/kg, <0) 
     513                  zEi             = - e_i_1d(ji,jk) * r1_rhoi                                 ! Specific enthalpy of melting ice (J/kg, <0) 
    514514                  zEw             = rcp * ztmelts                                             ! Specific enthalpy of meltwater (J/kg, <0) 
    515515                  zdE             = zEi - zEw                                                 ! Specific enthalpy difference   (J/kg, <0) 
     
    517517                  zfmdt           = - zq_bo(ji) / zdE                                         ! Mass flux x time step (kg/m2, >0) 
    518518 
    519                   zdeltah(ji,jk)  = - zfmdt * r1_rhoic                                        ! Gross thickness change 
     519                  zdeltah(ji,jk)  = - zfmdt * r1_rhoi                                         ! Gross thickness change 
    520520 
    521521                  zdeltah(ji,jk)  = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) )       ! bound thickness change 
    522522                   
    523                   zq_bo(ji)       = MAX( 0._wp , zq_bo(ji) - zdeltah(ji,jk) * rhoic * zdE )   ! update available heat. MAX is necessary for roundup errors 
    524  
    525                   dh_i_bom(ji)    = dh_i_bom(ji) + zdeltah(ji,jk)                            ! Update basal melt 
    526  
    527                   zfmdt           = - zdeltah(ji,jk) * rhoic                                  ! Mass flux x time step > 0 
     523                  zq_bo(ji)       = MAX( 0._wp , zq_bo(ji) - zdeltah(ji,jk) * rhoi * zdE )    ! update available heat. MAX is necessary for roundup errors 
     524 
     525                  dh_i_bom(ji)    = dh_i_bom(ji) + zdeltah(ji,jk)                             ! Update basal melt 
     526 
     527                  zfmdt           = - zdeltah(ji,jk) * rhoi                                   ! Mass flux x time step > 0 
    528528 
    529529                  zQm             = zfmdt * zEw                                               ! Heat exchanged with ocean 
    530530 
    531                   hfx_thd_1d(ji)  = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice                           ! Heat flux to the ocean [W.m-2], <0   
    532                   hfx_bom_1d(ji)  = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice                           ! Heat used in this process [W.m-2], >0   
    533  
    534                   sfx_bom_1d(ji)  = sfx_bom_1d(ji) - rhoic *  a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice  ! Salt flux 
     531                  hfx_thd_1d(ji)  = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_Dt_ice                           ! Heat flux to the ocean [W.m-2], <0   
     532                  hfx_bom_1d(ji)  = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_Dt_ice                           ! Heat used in this process [W.m-2], >0   
     533 
     534                  sfx_bom_1d(ji)  = sfx_bom_1d(ji) - rhoi *  a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_Dt_ice  ! Salt flux 
    535535                  !                                                                                                   using s_i_1d and not sz_i_1d(jk) is ok 
    536536                   
    537                   wfx_bom_1d(ji)  = wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice                ! Mass flux 
     537                  wfx_bom_1d(ji)  = wfx_bom_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice                ! Mass flux 
    538538 
    539539                  ! update heat content (J.m-2) and layer thickness 
     
    565565         
    566566         zq_rema(ji)        = zq_rema(ji) + zdeltah(ji,1) * e_s_1d(ji,1)                               ! update available heat (J.m-2) 
    567          hfx_snw_1d(ji)     = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * e_s_1d(ji,1) * r1_rdtice   ! Heat used to melt snow, W.m-2 (>0) 
    568          wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice      ! Mass flux 
     567         hfx_snw_1d(ji)     = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * e_s_1d(ji,1) * r1_Dt_ice   ! Heat used to melt snow, W.m-2 (>0) 
     568         wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhos * a_i_1d(ji) * zdeltah(ji,1) * r1_Dt_ice      ! Mass flux 
    569569         dh_s_mlt(ji)       = dh_s_mlt(ji) + zdeltah(ji,1) 
    570570         !     
    571571         ! Remaining heat flux (W.m-2) is sent to the ocean heat budget 
    572          hfx_out_1d(ji) = hfx_out_1d(ji) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_rdtice 
     572         hfx_out_1d(ji) = hfx_out_1d(ji) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_Dt_ice 
    573573 
    574574         IF( ln_icectl .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji) 
     
    580580      ! When snow load excesses Archimede's limit, snow-ice interface goes down under sea-level,  
    581581      ! flooding of seawater transforms snow into ice dh_snowice is positive for the ice 
    582       z1_rho = 1._wp / ( rhosn+rau0-rhoic ) 
     582      z1_rho = 1._wp / ( rhos + rho0 - rhoi ) 
    583583      DO ji = 1, npti 
    584584         ! 
    585          dh_snowice(ji) = MAX(  0._wp , ( rhosn * h_s_1d(ji) + (rhoic-rau0) * h_i_1d(ji) ) * z1_rho ) 
     585         dh_snowice(ji) = MAX(  0._wp , ( rhos * h_s_1d(ji) + (rhoi-rho0) * h_i_1d(ji) ) * z1_rho ) 
    586586 
    587587         h_i_1d(ji)    = h_i_1d(ji) + dh_snowice(ji) 
     
    589589 
    590590         ! Contribution to energy flux to the ocean [J/m2], >0 (if sst<0) 
    591          zfmdt          = ( rhosn - rhoic ) * dh_snowice(ji)    ! <0 
     591         zfmdt          = ( rhos - rhoi ) * dh_snowice(ji)    ! <0 
    592592         zEw            = rcp * sst_1d(ji) 
    593593         zQm            = zfmdt * zEw  
    594594          
    595          hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice ! Heat flux 
    596  
    597          sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_1d(ji) * a_i_1d(ji) * zfmdt * r1_rdtice ! Salt flux 
     595         hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_Dt_ice ! Heat flux 
     596 
     597         sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_1d(ji) * a_i_1d(ji) * zfmdt * r1_Dt_ice ! Salt flux 
    598598 
    599599         ! Case constant salinity in time: virtual salt flux to keep salinity constant 
    600600         IF( nn_icesal /= 2 )  THEN 
    601             sfx_bri_1d(ji) = sfx_bri_1d(ji) - sss_1d (ji) * a_i_1d(ji) * zfmdt                  * r1_rdtice  & ! put back sss_m     into the ocean 
    602                &                            - s_i_1d(ji)  * a_i_1d(ji) * dh_snowice(ji) * rhoic * r1_rdtice    ! and get  rn_icesal from the ocean  
     601            sfx_bri_1d(ji) = sfx_bri_1d(ji) - sss_1d (ji) * a_i_1d(ji) * zfmdt                 * r1_Dt_ice  & ! put back sss_m     into the ocean 
     602               &                            - s_i_1d(ji)  * a_i_1d(ji) * dh_snowice(ji) * rhoi * r1_Dt_ice    ! and get  rn_icesal from the ocean  
    603603         ENDIF 
    604604 
    605605         ! Mass flux: All snow is thrown in the ocean, and seawater is taken to replace the volume 
    606          wfx_sni_1d(ji)     = wfx_sni_1d(ji)     - a_i_1d(ji) * dh_snowice(ji) * rhoic * r1_rdtice 
    607          wfx_snw_sni_1d(ji) = wfx_snw_sni_1d(ji) + a_i_1d(ji) * dh_snowice(ji) * rhosn * r1_rdtice 
     606         wfx_sni_1d(ji)     = wfx_sni_1d(ji)     - a_i_1d(ji) * dh_snowice(ji) * rhoi * r1_Dt_ice 
     607         wfx_snw_sni_1d(ji) = wfx_snw_sni_1d(ji) + a_i_1d(ji) * dh_snowice(ji) * rhos * r1_Dt_ice 
    608608 
    609609         ! update heat content (J.m-2) and layer thickness 
     
    627627            e_s_1d(ji,jk) = rswitch * e_s_1d(ji,jk) 
    628628            ! recalculate t_s_1d from e_s_1d 
    629             t_s_1d(ji,jk) = rt0 + rswitch * ( - e_s_1d(ji,jk) * r1_rhosn * r1_cpic + lfus * r1_cpic ) 
     629            t_s_1d(ji,jk) = rt0 + rswitch * ( - e_s_1d(ji,jk) * r1_rhos * r1_cpi + rLfus * r1_cpi ) 
    630630         END DO 
    631631      END DO 
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/ICE/icethd_do.F90

    r9604 r9939  
    140140         ! Physical constants 
    141141         zhicrit = 0.04                                          ! frazil ice thickness 
    142          ztwogp  = 2. * rau0 / ( grav * 0.3 * ( rau0 - rhoic ) ) ! reduced grav 
     142         ztwogp  = 2. * rho0 / ( grav * 0.3 * ( rho0 - rhoi ) ) ! reduced grav 
    143143         zsqcd   = 1.0 / SQRT( 1.3 * zcai )                      ! 1/SQRT(airdensity*drag) 
    144144         zgamafr = 0.03 
     
    264264         DO ji = 1, npti 
    265265            ztmelts       = - tmut * zs_newice(ji)                  ! Melting point (C) 
    266             ze_newice(ji) =   rhoic * (  cpic * ( ztmelts - ( t_bo_1d(ji) - rt0 ) )                     & 
    267                &                       + lfus * ( 1.0 - ztmelts / MIN( t_bo_1d(ji) - rt0, -epsi10 ) )   & 
    268                &                       - rcp  *         ztmelts ) 
     266            ze_newice(ji) =   rhoi * (  rcpi * ( ztmelts - ( t_bo_1d(ji) - rt0 ) )                     & 
     267               &                      + rLfus * ( 1.0 - ztmelts / MIN( t_bo_1d(ji) - rt0, -epsi10 ) )   & 
     268               &                      - rcp   *         ztmelts ) 
    269269         END DO 
    270270 
     
    275275         DO ji = 1, npti 
    276276 
    277             zEi           = - ze_newice(ji) * r1_rhoic             ! specific enthalpy of forming ice [J/kg] 
     277            zEi           = - ze_newice(ji) * r1_rhoi              ! specific enthalpy of forming ice [J/kg] 
    278278 
    279279            zEw           = rcp * ( t_bo_1d(ji) - rt0 )            ! specific enthalpy of seawater at t_bo_1d [J/kg] 
     
    284284            zfmdt         = - qlead_1d(ji) / zdE                   ! Fm.dt [kg/m2] (<0)  
    285285                                                                   ! clem: we use qlead instead of zqld (icethd) because we suppose we are at the freezing point    
    286             zv_newice(ji) = - zfmdt * r1_rhoic 
     286            zv_newice(ji) = - zfmdt * r1_rhoi 
    287287 
    288288            zQm           = zfmdt * zEw                            ! heat to the ocean >0 associated with mass flux   
    289289 
    290290            ! Contribution to heat flux to the ocean [W.m-2], >0   
    291             hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * zEw * r1_rdtice 
     291            hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * zEw * r1_Dt_ice 
    292292            ! Total heat flux used in this process [W.m-2]   
    293             hfx_opw_1d(ji) = hfx_opw_1d(ji) - zfmdt * zdE * r1_rdtice 
     293            hfx_opw_1d(ji) = hfx_opw_1d(ji) - zfmdt * zdE * r1_Dt_ice 
    294294            ! mass flux 
    295             wfx_opw_1d(ji) = wfx_opw_1d(ji) - zv_newice(ji) * rhoic * r1_rdtice 
     295            wfx_opw_1d(ji) = wfx_opw_1d(ji) - zv_newice(ji) * rhoi * r1_Dt_ice 
    296296            ! salt flux 
    297             sfx_opw_1d(ji) = sfx_opw_1d(ji) - zv_newice(ji) * rhoic * zs_newice(ji) * r1_rdtice 
     297            sfx_opw_1d(ji) = sfx_opw_1d(ji) - zv_newice(ji) * rhoi * zs_newice(ji) * r1_Dt_ice 
    298298         END DO 
    299299          
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/ICE/icethd_ent.F90

    r9604 r9939  
    129129      ! then we should not (* a_i) again but not important since this is just to check that remap error is ~0 
    130130      DO ji = 1, npti 
    131          hfx_err_rem_1d(ji) = hfx_err_rem_1d(ji) + a_i_1d(ji) * r1_rdtice *  & 
     131         hfx_err_rem_1d(ji) = hfx_err_rem_1d(ji) + a_i_1d(ji) * r1_Dt_ice *  & 
    132132            &               ( SUM( qnew(ji,1:nlay_i) ) * zhnew(ji) - SUM( eh_i_old(ji,0:nlay_i+1) ) )  
    133133      END DO 
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/ICE/icethd_pnd.F90

    r9750 r9939  
    133133      REAL(wp) ::   zdv_mlt          ! available meltwater for melt ponding 
    134134      REAL(wp) ::   z1_Tp            ! inverse reference temperature 
    135       REAL(wp) ::   z1_rhofw         ! inverse freshwater density 
    136135      REAL(wp) ::   z1_zpnd_aspect   ! inverse pond aspect ratio 
    137136      REAL(wp) ::   zfac, zdum 
     
    139138      INTEGER  ::   ji   ! loop indices 
    140139      !!------------------------------------------------------------------- 
    141       z1_rhofw       = 1._wp / rhofw  
    142140      z1_zpnd_aspect = 1._wp / zpnd_aspect 
    143141      z1_Tp          = 1._wp / zTp  
     
    157155            ! 
    158156            ! available meltwater for melt ponding [m, >0] and fraction 
    159             zdv_mlt = -( dh_i_sum(ji)*rhoic + dh_s_mlt(ji)*rhosn ) * z1_rhofw * a_i_1d(ji) 
     157            zdv_mlt = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * r1_rhow * a_i_1d(ji) 
    160158            zfr_mlt = zrmin + ( zrmax - zrmin ) * a_i_1d(ji)  ! from CICE doc 
    161159            !zfr_mlt = zrmin + zrmax * a_i_1d(ji)             ! from Holland paper  
     
    168166            ! melt pond mass flux (<0) 
    169167            IF( ln_pnd_fwb .AND. zdv_mlt > 0._wp ) THEN 
    170                zfac = zfr_mlt * zdv_mlt * rhofw * r1_rdtice 
     168               zfac = zfr_mlt * zdv_mlt * rhow * r1_Dt_ice 
    171169               wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac 
    172170               ! 
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/ICE/icethd_sal.F90

    r9750 r9939  
    7777            !--------------------------------------------------------- 
    7878            IF( h_i_1d(ji) > 0._wp ) THEN 
    79                zs_sni   = sss_1d(ji) * ( rhoic - rhosn ) * r1_rhoic                 ! Salinity of snow ice 
     79               zs_sni   = sss_1d(ji) * ( rhoi - rhos ) * r1_rhoi                    ! Salinity of snow ice 
    8080               zs_i_si = ( zs_sni      - s_i_1d(ji) ) * dh_snowice(ji) / h_i_1d(ji) ! snow-ice     
    8181               zs_i_bg = ( s_i_new(ji) - s_i_1d(ji) ) * dh_i_bog  (ji) / h_i_1d(ji) ! bottom growth 
     
    9898                
    9999               ! Salt flux 
    100                sfx_bri_1d(ji) = sfx_bri_1d(ji) - rhoic * a_i_1d(ji) * h_i_1d(ji) * ( zs_i_fl + zs_i_gd ) * r1_rdtice 
     100               sfx_bri_1d(ji) = sfx_bri_1d(ji) - rhoi * a_i_1d(ji) * h_i_1d(ji) * ( zs_i_fl + zs_i_gd ) * r1_Dt_ice 
    101101            ENDIF 
    102102         END DO 
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/ICE/icethd_zdf_bl99.F90

    r9656 r9939  
    217217            ! 
    218218            DO ji = 1, npti 
    219                ztcond_i(ji,0)      = rcdic + zbeta * sz_i_1d(ji,1)      / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) 
    220                ztcond_i(ji,nlay_i) = rcdic + zbeta * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji)  - rt0 ) 
     219               ztcond_i(ji,0)      = rcnd_i + zbeta * sz_i_1d(ji,1)      / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) 
     220               ztcond_i(ji,nlay_i) = rcnd_i + zbeta * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji)  - rt0 ) 
    221221            END DO 
    222222            DO jk = 1, nlay_i-1 
    223223               DO ji = 1, npti 
    224                   ztcond_i(ji,jk) = rcdic + zbeta * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) /  & 
    225                      &                      MIN( -epsi10, 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 ) 
     224!!gm faster coding 
     225                  ztcond_i(ji,jk) = rcnd_i + zbeta * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) )             & 
     226                     &                          / MIN( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) - 2._wp * rt0 , -epsi10 ) 
     227!!gm old 
     228!                  ztcond_i(ji,jk) = rcnd_i + zbeta * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) )   & 
     229!                     &                             / MIN( 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 , -epsi10 ) 
     230!!gm  
    226231               END DO 
    227232            END DO 
     
    230235            ! 
    231236            DO ji = 1, npti 
    232                ztcond_i(ji,0)      = rcdic + 0.09_wp  *  sz_i_1d(ji,1)      / MIN( -epsi10, t_i_1d(ji,1) - rt0 )  & 
    233                   &                        - 0.011_wp * ( t_i_1d(ji,1) - rt0 ) 
    234                ztcond_i(ji,nlay_i) = rcdic + 0.09_wp  *  sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji)  - rt0 )  & 
    235                   &                        - 0.011_wp * ( t_bo_1d(ji) - rt0 ) 
     237               ztcond_i(ji,0)      = rcnd_i + 0.09_wp  *  sz_i_1d(ji,1)      / MIN( -epsi10, t_i_1d(ji,1) - rt0 )  & 
     238                  &                         - 0.011_wp * ( t_i_1d(ji,1) - rt0 ) 
     239               ztcond_i(ji,nlay_i) = rcnd_i + 0.09_wp  *  sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji)  - rt0 )  & 
     240                  &                         - 0.011_wp * ( t_bo_1d(ji) - rt0 ) 
    236241            END DO 
    237242            DO jk = 1, nlay_i-1 
    238243               DO ji = 1, npti 
    239                   ztcond_i(ji,jk) = rcdic + 0.09_wp  *   0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) /        & 
    240                      &                     MIN( -epsi10, 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 )  & 
    241                      &                    - 0.011_wp * ( 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 ) 
     244!!gm faster coding 
     245                  zfac   = t_i_1d (ji,jk) + t_i_1d (ji,jk+1) - 2._wp * rt0 
     246                  ztcond_i(ji,jk) = rcnd_i + 0.09_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) / MIN( -epsi10, zfac )  & 
     247                     &                     - 0.0055_wp * zfac            !NB:  0.0055 = 1/2 * 0.011 
     248!!gm old 
     249!                  ztcond_i(ji,jk) = rcnd_i + 0.09_wp  *   0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) /        & 
     250!                     &                      MIN( -epsi10, 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 )  & 
     251!                     &                     - 0.011_wp * ( 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 ) 
     252!!gm  
    242253               END DO 
    243254            END DO 
     
    299310         DO jk = 1, nlay_i 
    300311            DO ji = 1, npti 
    301                zcpi = cpic + zgamma * sz_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztiold(ji,jk) - rt0 ), epsi10 ) 
    302                zeta_i(ji,jk) = rdt_ice * r1_rhoic * z1_h_i(ji) / MAX( epsi10, zcpi )  
     312               zcpi = rcpi + zgamma * sz_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztiold(ji,jk) - rt0 ), epsi10 ) 
     313               zeta_i(ji,jk) = rdt_ice * r1_rhoi * z1_h_i(ji) / MAX( epsi10, zcpi )  
    303314            END DO 
    304315         END DO 
     
    306317         DO jk = 1, nlay_s 
    307318            DO ji = 1, npti 
    308                zeta_s(ji,jk) = rdt_ice * r1_rhosn * r1_cpic * z1_h_s(ji) 
     319               zeta_s(ji,jk) = rdt_ice * r1_rhos * r1_cpi * z1_h_s(ji) 
    309320            END DO 
    310321         END DO 
     
    770781                
    771782               IF( t_su_1d(ji) < rt0 ) THEN  ! case T_su < 0degC 
    772                   zhfx_err = ( qns_ice_1d(ji) + qsr_ice_1d(ji)    - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice )*a_i_1d(ji) 
     783                  zhfx_err = ( qns_ice_1d(ji) + qsr_ice_1d(ji)    - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_Dt_ice )*a_i_1d(ji) 
    773784               ELSE                          ! case T_su = 0degC 
    774                   zhfx_err = ( fc_su(ji)      + qsr_ice_tr_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice )*a_i_1d(ji) 
     785                  zhfx_err = ( fc_su(ji)      + qsr_ice_tr_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_Dt_ice )*a_i_1d(ji) 
    775786               ENDIF 
    776787                
    777788            ELSEIF( k_jules == np_jules_ACTIVE ) THEN 
    778789             
    779                zhfx_err = ( fc_su(ji) + qsr_ice_tr_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji) 
     790               zhfx_err = ( fc_su(ji) + qsr_ice_tr_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_Dt_ice ) * a_i_1d(ji) 
    780791             
    781792            ENDIF 
     
    785796            ! 
    786797            ! hfx_dif = Heat flux diagnostic of sensible heat used to warm/cool ice in W.m-2    
    787             hfx_dif_1d(ji) = hfx_dif_1d(ji) - zdq * r1_rdtice * a_i_1d(ji) 
     798            hfx_dif_1d(ji) = hfx_dif_1d(ji) - zdq * r1_Dt_ice * a_i_1d(ji) 
    788799            ! 
    789800         END DO 
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/ICE/iceupdate.F90

    r9784 r9939  
    172172            snwice_mass_b(ji,jj) = snwice_mass(ji,jj)       ! save mass from the previous ice time step 
    173173            !                                               ! new mass per unit area 
    174             snwice_mass  (ji,jj) = tmask(ji,jj,1) * ( rhosn * vt_s(ji,jj) + rhoic * vt_i(ji,jj)  )  
     174            snwice_mass  (ji,jj) = tmask(ji,jj,1) * ( rhos * vt_s(ji,jj) + rhoi * vt_i(ji,jj)  )  
    175175            !                                               ! time evolution of snow+ice mass 
    176             snwice_fmass (ji,jj) = ( snwice_mass(ji,jj) - snwice_mass_b(ji,jj) ) * r1_rdtice 
     176            snwice_fmass (ji,jj) = ( snwice_mass(ji,jj) - snwice_mass_b(ji,jj) ) * r1_Dt_ice 
    177177             
    178178         END DO 
     
    336336      ENDIF 
    337337 
    338       zrhoco = rau0 * rn_cio 
     338      zrhoco = rho0 * rn_cio 
    339339      ! 
    340340      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN     !==  Ice time-step only  ==!   (i.e. surface module time-step) 
     
    432432            ELSE                                     ! start from rest 
    433433               IF(lwp) WRITE(numout,*) '   ==>>   previous run without snow-ice mass output then set it' 
    434                snwice_mass  (:,:) = tmask(:,:,1) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:) ) 
     434               snwice_mass  (:,:) = tmask(:,:,1) * ( rhos * vt_s(:,:) + rhoi * vt_i(:,:) ) 
    435435               snwice_mass_b(:,:) = snwice_mass(:,:) 
    436436            ENDIF 
    437437         ELSE                                   !* Start from rest 
    438438            IF(lwp) WRITE(numout,*) '   ==>>   start from rest: set the snow-ice mass' 
    439             snwice_mass  (:,:) = tmask(:,:,1) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:) ) 
     439            snwice_mass  (:,:) = tmask(:,:,1) * ( rhos * vt_s(:,:) + rhoi * vt_i(:,:) ) 
    440440            snwice_mass_b(:,:) = snwice_mass(:,:) 
    441441         ENDIF 
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/ICE/icevar.F90

    r9725 r9939  
    228228                     ztmelts          = - sz_i(ji,jj,jk,jl) * tmut                                 ! Ice layer melt temperature [C] 
    229229                     ! Conversion q(S,T) -> T (second order equation) 
    230                      zbbb             = ( rcp - cpic ) * ztmelts + ze_i * r1_rhoic - lfus 
    231                      zccc             = SQRT( MAX( zbbb * zbbb - 4._wp * cpic * lfus * ztmelts , 0._wp) ) 
    232                      t_i(ji,jj,jk,jl) = MAX( -100._wp , MIN( -( zbbb + zccc ) * 0.5_wp * r1_cpic , ztmelts ) ) + rt0   ! [K] with bounds: -100 < t_i < ztmelts 
     230                     zbbb             = ( rcp - rcpi ) * ztmelts + ze_i * r1_rhoi - rLfus 
     231                     zccc             = SQRT( MAX( zbbb * zbbb - 4._wp * rcpi * rLfus * ztmelts , 0._wp) ) 
     232                     t_i(ji,jj,jk,jl) = MAX( -100._wp , MIN( -( zbbb + zccc ) * 0.5_wp * r1_cpi , ztmelts ) ) + rt0   ! [K] with bounds: -100 < t_i < ztmelts 
    233233                     ! 
    234234                  ELSE                                !--- no ice 
     
    247247         WHERE( v_s(:,:,:) > epsi20 )        !--- icy area 
    248248            t_s(:,:,jk,:) = rt0 + MAX( -100._wp ,  & 
    249                  &                MIN( r1_cpic * ( -r1_rhosn * ( e_s(:,:,jk,:) / v_s(:,:,:) * zlay_s ) + lfus ) , 0._wp ) ) 
     249                 &                MIN( r1_cpi * ( -r1_rhos * ( e_s(:,:,jk,:) / v_s(:,:,:) * zlay_s ) + rLfus ) , 0._wp ) ) 
    250250         ELSEWHERE                           !--- no ice 
    251251            t_s(:,:,jk,:) = rt0 
     
    477477               DO ji = 1 , jpi 
    478478                  ! update exchanges with ocean 
    479                   hfx_res(ji,jj)   = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_i(ji,jj,jk,jl) * r1_rdtice ! W.m-2 <0 
     479                  hfx_res(ji,jj)   = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_i(ji,jj,jk,jl) * r1_Dt_ice ! W.m-2 <0 
    480480                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * zswitch(ji,jj) 
    481481                  t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) ) 
     
    488488               DO ji = 1 , jpi 
    489489                  ! update exchanges with ocean 
    490                   hfx_res(ji,jj)   = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_s(ji,jj,jk,jl) * r1_rdtice ! W.m-2 <0 
     490                  hfx_res(ji,jj)   = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_s(ji,jj,jk,jl) * r1_Dt_ice ! W.m-2 <0 
    491491                  e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * zswitch(ji,jj) 
    492492                  t_s(ji,jj,jk,jl) = t_s(ji,jj,jk,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) ) 
     
    498498            DO ji = 1 , jpi 
    499499               ! update exchanges with ocean 
    500                sfx_res(ji,jj)  = sfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * sv_i(ji,jj,jl)   * rhoic * r1_rdtice 
    501                wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_i (ji,jj,jl)   * rhoic * r1_rdtice 
    502                wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl)   * rhosn * r1_rdtice 
     500               sfx_res(ji,jj)  = sfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * sv_i(ji,jj,jl)   * rhoi * r1_Dt_ice 
     501               wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_i (ji,jj,jl)   * rhoi * r1_Dt_ice 
     502               wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl)   * rhos * r1_Dt_ice 
    503503               ! 
    504504               !----------------------------------------------------------------- 
     
    669669               ! In case snow load is in excess that would lead to transformation from snow to ice 
    670670               ! Then, transfer the snow excess into the ice (different from icethd_dh) 
    671                zdh = MAX( 0._wp, ( rhosn * zh_s(ji,jl) + ( rhoic - rau0 ) * zh_i(ji,jl) ) * r1_rau0 )  
     671               zdh = MAX( 0._wp, ( rhos * zh_s(ji,jl) + ( rhoi - rho0 ) * zh_i(ji,jl) ) * r1_rho0 )  
    672672               ! recompute h_i, h_s avoiding out of bounds values 
    673673               zh_i(ji,jl) = MIN( hi_max(jl), zh_i(ji,jl) + zdh ) 
    674                zh_s(ji,jl) = MAX( 0._wp, zh_s(ji,jl) - zdh * rhoic * r1_rhosn ) 
     674               zh_s(ji,jl) = MAX( 0._wp, zh_s(ji,jl) - zdh * rhoi * r1_rhos ) 
    675675            ENDIF 
    676676         END DO 
     
    854854            ztmelts      = - tmut  * sz_i_1d(ji,jk) 
    855855            t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts + rt0 ) ! Force t_i_1d to be lower than melting point 
    856                                                                 !   (sometimes zdf scheme produces abnormally high temperatures)    
    857             e_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - ( t_i_1d(ji,jk) - rt0 ) )           & 
    858                &                    + lfus * ( 1._wp - ztmelts / ( t_i_1d(ji,jk) - rt0 ) )   & 
    859                &                    - rcp  *   ztmelts ) 
     856            !                                                   !   (sometimes zdf scheme produces abnormally high temperatures)    
     857            e_i_1d(ji,jk) = rhoi * ( rcpi * ( ztmelts - ( t_i_1d(ji,jk) - rt0 ) )           & 
     858               &                   + rLfus * ( 1._wp - ztmelts / ( t_i_1d(ji,jk) - rt0 ) )   & 
     859               &                   - rcp   *   ztmelts ) 
    860860         END DO 
    861861      END DO 
    862862      DO jk = 1, nlay_s             ! Snow energy of melting 
    863863         DO ji = 1, npti 
    864             e_s_1d(ji,jk) = rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) 
     864            e_s_1d(ji,jk) = rhos * ( rcpi * ( rt0 - t_s_1d(ji,jk) ) + rLfus ) 
    865865         END DO 
    866866      END DO 
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/ICE/icewri.F90

    r9604 r9939  
    8585      ! Standard outputs 
    8686      !----------------- 
    87       zrho1 = ( rau0 - rhoic ) * r1_rau0; zrho2 = rhosn * r1_rau0 
     87      zrho1 = ( rho0 - rhoi ) * r1_rho0   ;   zrho2 = rhos * r1_rho0 
    8888      ! masks 
    8989      IF( iom_use('icemask'  ) )   CALL iom_put( "icemask"  , zmsk00              )   ! ice mask 0% 
     
    9292      ! 
    9393      ! general fields 
    94       IF( iom_use('icemass'  ) )   CALL iom_put( "icemass", rhoic * vt_i * zmsk00 )   ! Ice mass per cell area  
    95       IF( iom_use('snwmass'  ) )   CALL iom_put( "snwmass", rhosn * vt_s * zmsksn )   ! Snow mass per cell area 
     94      IF( iom_use('icemass'  ) )   CALL iom_put( "icemass", rhoi * vt_i * zmsk00 )   ! Ice mass per cell area  
     95      IF( iom_use('snwmass'  ) )   CALL iom_put( "snwmass", rhos * vt_s * zmsksn )   ! Snow mass per cell area 
    9696      IF( iom_use('icepres'  ) )   CALL iom_put( "icepres", zmsk00                )   ! Ice presence (1 or 0)  
    9797      IF( iom_use('iceconc'  ) )   CALL iom_put( "iceconc", at_i  * zmsk00        )   ! ice concentration 
     
    104104      IF( iom_use('snwvolu'  ) )   CALL iom_put( "snwvolu", vt_s  * zmsksn        )   ! snow volume 
    105105      IF( iom_use('icefrb') ) THEN 
     106!!gm remove the WHERE by using : 
     107!!         z2d(:,:) = MAX( zrho1 * hm_i(:,:) - zrho2 * hm_s(:,:) , 0._wp )                                          
     108!!gm end 
    106109         z2d(:,:) = ( zrho1 * hm_i(:,:) - zrho2 * hm_s(:,:) )                                          
    107110         WHERE( z2d < 0._wp )   z2d = 0._wp 
     
    115118      ! salt 
    116119      IF( iom_use('icesalt'  ) )   CALL iom_put( "icesalt", sm_i  * zmsk00        )   ! mean ice salinity 
    117       IF( iom_use('icesalm'  ) )   CALL iom_put( "icesalm", SUM( sv_i, DIM = 3 ) * rhoic * 1.0e-3 * zmsk00 )   ! Mass of salt in sea ice per cell area 
     120      IF( iom_use('icesalm'  ) )   CALL iom_put( "icesalm", SUM( sv_i, DIM = 3 ) * rhoi * 1.0e-3 * zmsk00 )   ! Mass of salt in sea ice per cell area 
    118121 
    119122      ! heat 
     
    164167      ! trends 
    165168      IF( iom_use('dmithd') )   CALL iom_put( "dmithd", - wfx_bog - wfx_bom - wfx_sum - wfx_sni - wfx_opw - wfx_lam - wfx_res ) ! Sea-ice mass change from thermodynamics 
    166       IF( iom_use('dmidyn') )   CALL iom_put( "dmidyn", - wfx_dyn + rhoic * diag_trp_vi     )   ! Sea-ice mass change from dynamics(kg/m2/s) 
     169      IF( iom_use('dmidyn') )   CALL iom_put( "dmidyn", - wfx_dyn + rhoi * diag_trp_vi      )   ! Sea-ice mass change from dynamics(kg/m2/s) 
    167170      IF( iom_use('dmiopw') )   CALL iom_put( "dmiopw", - wfx_opw                           )   ! Sea-ice mass change through growth in open water 
    168171      IF( iom_use('dmibog') )   CALL iom_put( "dmibog", - wfx_bog                           )   ! Sea-ice mass change through basal growth 
     
    174177      IF( iom_use('dmisub') )   CALL iom_put( "dmisub", - wfx_ice_sub                       )   ! Sea-ice mass change through sublimation 
    175178      IF( iom_use('dmsspr') )   CALL iom_put( "dmsspr", - wfx_spr                           )   ! Snow mass change through snow fall 
    176       IF( iom_use('dmsssi') )   CALL iom_put( "dmsssi",   wfx_sni*rhosn*r1_rhoic            )   ! Snow mass change through snow-to-ice conversion 
     179      IF( iom_use('dmsssi') )   CALL iom_put( "dmsssi",   wfx_sni*rhos*r1_rhoi              )   ! Snow mass change through snow-to-ice conversion 
    177180      IF( iom_use('dmsmel') )   CALL iom_put( "dmsmel", - wfx_snw_sum                       )   ! Snow mass change through melt 
    178       IF( iom_use('dmsdyn') )   CALL iom_put( "dmsdyn", - wfx_snw_dyn + rhosn * diag_trp_vs )   ! Snow mass change through dynamics(kg/m2/s) 
     181      IF( iom_use('dmsdyn') )   CALL iom_put( "dmsdyn", - wfx_snw_dyn + rhos * diag_trp_vs )   ! Snow mass change through dynamics(kg/m2/s) 
    179182 
    180183      ! Global ice diagnostics 
     
    250253      CALL histvert( kid, "ncatice", "Ice Categories","", jpl, jcat, nz_i, "up") 
    251254 
    252       CALL histdef( kid, "sithic", "Ice thickness"          , "m"      , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    253       CALL histdef( kid, "siconc", "Ice concentration"      , "%"      , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    254       CALL histdef( kid, "sitemp", "Ice temperature"        , "C"      , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    255       CALL histdef( kid, "sivelu", "i-Ice speed "           , "m/s"    , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    256       CALL histdef( kid, "sivelv", "j-Ice speed "           , "m/s"    , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )  
    257       CALL histdef( kid, "sistru", "i-Wind stress over ice" , "Pa"     , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    258       CALL histdef( kid, "sistrv", "j-Wind stress over ice" , "Pa"     , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )  
    259       CALL histdef( kid, "sisflx", "Solar flx over ocean"   , "W/m2"   , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )  
    260       CALL histdef( kid, "sinflx", "NonSolar flx over ocean", "W/m2"   , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    261       CALL histdef( kid, "snwpre", "Snow precipitation"     , "kg/m2/s", jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )  
    262       CALL histdef( kid, "sisali", "Ice salinity"           , "PSU"    , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )  
    263       CALL histdef( kid, "sivolu", "Ice volume"             , "m"      , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )  
    264       CALL histdef( kid, "sidive", "Ice divergence"         , "10-8s-1", jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )  
    265       CALL histdef( kid, "si_amp", "Melt pond fraction"     , "%"      , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    266       CALL histdef( kid, "si_vmp", "Melt pond volume"       ,  "m"     , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    267       ! 
    268       CALL histdef( kid, "sithicat", "Ice thickness"        , "m"      , jpi,jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt ) 
    269       CALL histdef( kid, "siconcat", "Ice concentration"    , "%"      , jpi,jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt ) 
    270       CALL histdef( kid, "sisalcat", "Ice salinity"         , ""       , jpi,jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt ) 
    271       CALL histdef( kid, "snthicat", "Snw thickness"        , "m"      , jpi,jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt ) 
     255      CALL histdef( kid, "sithic", "Ice thickness"          , "m"      , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rn_Dt, rn_Dt ) 
     256      CALL histdef( kid, "siconc", "Ice concentration"      , "%"      , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rn_Dt, rn_Dt ) 
     257      CALL histdef( kid, "sitemp", "Ice temperature"        , "C"      , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rn_Dt, rn_Dt ) 
     258      CALL histdef( kid, "sivelu", "i-Ice speed "           , "m/s"    , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rn_Dt, rn_Dt ) 
     259      CALL histdef( kid, "sivelv", "j-Ice speed "           , "m/s"    , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rn_Dt, rn_Dt )  
     260      CALL histdef( kid, "sistru", "i-Wind stress over ice" , "Pa"     , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rn_Dt, rn_Dt ) 
     261      CALL histdef( kid, "sistrv", "j-Wind stress over ice" , "Pa"     , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rn_Dt, rn_Dt )  
     262      CALL histdef( kid, "sisflx", "Solar flx over ocean"   , "W/m2"   , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rn_Dt, rn_Dt )  
     263      CALL histdef( kid, "sinflx", "NonSolar flx over ocean", "W/m2"   , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rn_Dt, rn_Dt ) 
     264      CALL histdef( kid, "snwpre", "Snow precipitation"     , "kg/m2/s", jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rn_Dt, rn_Dt )  
     265      CALL histdef( kid, "sisali", "Ice salinity"           , "PSU"    , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rn_Dt, rn_Dt )  
     266      CALL histdef( kid, "sivolu", "Ice volume"             , "m"      , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rn_Dt, rn_Dt )  
     267      CALL histdef( kid, "sidive", "Ice divergence"         , "10-8s-1", jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rn_Dt, rn_Dt )  
     268      CALL histdef( kid, "si_amp", "Melt pond fraction"     , "%"      , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rn_Dt, rn_Dt ) 
     269      CALL histdef( kid, "si_vmp", "Melt pond volume"       ,  "m"     , jpi,jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rn_Dt, rn_Dt ) 
     270      ! 
     271      CALL histdef( kid, "sithicat", "Ice thickness"        , "m" , jpi,jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rn_Dt, rn_Dt ) 
     272      CALL histdef( kid, "siconcat", "Ice concentration"    , "%" , jpi,jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rn_Dt, rn_Dt ) 
     273      CALL histdef( kid, "sisalcat", "Ice salinity"         , ""  , jpi,jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rn_Dt, rn_Dt ) 
     274      CALL histdef( kid, "snthicat", "Snw thickness"        , "m" , jpi,jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rn_Dt, rn_Dt ) 
    272275 
    273276      CALL histend( kid, snc4set )   ! end of the file definition 
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/NST/agrif_oce_update.F90

    r9780 r9939  
    1212   !!            3.6  !  2014-09  (R. Benshila)  
    1313   !!---------------------------------------------------------------------- 
     14 
     15   !!---------------------------------------------------------------------- 
     16   !!   Agrif_Update_Tra   : T-S agrif update 
     17   !!   Agrif_Update_Dyn   : dynamics agrif update 
     18   !!   Agrif_Update_ssh   : sea surface height update 
     19   !!   Agrif_Update_Tke   :  
     20   !!   Agrif_Update_vvl   :  
     21   !!   dom_vvl_update_UVF :  
     22   !!   updateTS           :  
     23   !!   updateu            : 
     24   !!   correct_u_bdy      : 
     25   !!   updatev            : 
     26   !!   correct_v_bdy      : 
     27   !!   updateu2d          : 
     28   !!   updatev2d          : 
     29   !!   updateSSH          : 
     30   !!   updateub2b         : 
     31   !!   reflux_sshu        : 
     32   !!   updatevb2b         : 
     33   !!   reflux_sshv        : 
     34   !!   update_scales      : 
     35   !!   updateEN           : 
     36   !!   updateAVT          : 
     37   !!   updateAVM          : 
     38   !!   updatee3t          : 
     39   !!---------------------------------------------------------------------- 
     40 
    1441#if defined key_agrif  
    1542   !!---------------------------------------------------------------------- 
    1643   !!   'key_agrif'                                              AGRIF zoom 
    1744   !!---------------------------------------------------------------------- 
    18    USE par_oce 
    19    USE oce 
    20    USE dom_oce 
     45   USE par_oce        ! ocean parameter 
     46   USE oce            ! ocean variables 
     47   USE dom_oce        ! ocean domain 
    2148   USE zdf_oce        ! vertical physics: ocean variables  
    22    USE agrif_oce 
     49   USE agrif_oce      !  
    2350   ! 
    2451   USE in_out_manager ! I/O manager 
     
    6794      ! 
    6895   END SUBROUTINE Agrif_Update_Tra 
     96 
    6997 
    7098   SUBROUTINE Agrif_Update_Dyn( ) 
     
    125153   END SUBROUTINE Agrif_Update_Dyn 
    126154 
     155 
    127156   SUBROUTINE Agrif_Update_ssh( ) 
    128       !!--------------------------------------------- 
    129       !!   *** ROUTINE Agrif_Update_ssh *** 
    130       !!--------------------------------------------- 
     157      !!---------------------------------------------------------------------- 
     158      !!                   *** ROUTINE Agrif_Update_ssh *** 
     159      !!---------------------------------------------------------------------- 
    131160      !  
    132161      IF (Agrif_Root()) RETURN 
     
    163192 
    164193   SUBROUTINE Agrif_Update_Tke( ) 
    165       !!--------------------------------------------- 
    166       !!   *** ROUTINE Agrif_Update_Tke *** 
    167       !!--------------------------------------------- 
    168       !! 
     194      !!---------------------------------------------------------------------- 
     195      !!                   *** ROUTINE Agrif_Update_Tke *** 
     196      !!---------------------------------------------------------------------- 
    169197      !  
    170198      IF (Agrif_Root()) RETURN 
    171199      !        
    172200#  if defined TWO_WAY 
    173  
     201      ! 
    174202      Agrif_UseSpecialValueInUpdate = .TRUE. 
    175203      Agrif_SpecialValueFineGrid = 0. 
    176  
     204      ! 
    177205      CALL Agrif_Update_Variable( en_id, locupdate=(/0,0/), procname=updateEN  ) 
    178206      CALL Agrif_Update_Variable(avt_id, locupdate=(/0,0/), procname=updateAVT ) 
    179207      CALL Agrif_Update_Variable(avm_id, locupdate=(/0,0/), procname=updateAVM ) 
    180  
     208      ! 
    181209      Agrif_UseSpecialValueInUpdate = .FALSE. 
    182  
     210      ! 
    183211#  endif 
    184        
     212      ! 
    185213   END SUBROUTINE Agrif_Update_Tke 
    186214 
    187215 
    188216   SUBROUTINE Agrif_Update_vvl( ) 
    189       !!--------------------------------------------- 
    190       !!   *** ROUTINE Agrif_Update_vvl *** 
    191       !!--------------------------------------------- 
    192       ! 
    193       IF (Agrif_Root()) RETURN 
     217      !!---------------------------------------------------------------------- 
     218      !!                   *** ROUTINE Agrif_Update_vvl *** 
     219      !!---------------------------------------------------------------------- 
     220      ! 
     221      IF ( Agrif_Root() )  RETURN 
    194222      ! 
    195223#if defined TWO_WAY   
     
    214242   END SUBROUTINE Agrif_Update_vvl 
    215243 
     244 
    216245   SUBROUTINE dom_vvl_update_UVF 
    217       !!--------------------------------------------- 
    218       !!       *** ROUTINE dom_vvl_update_UVF *** 
    219       !!--------------------------------------------- 
    220       !! 
    221       INTEGER :: jk 
    222       REAL(wp):: zcoef 
    223       !!--------------------------------------------- 
    224  
    225       IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Finalize e3 on grid Number', & 
    226                   & Agrif_Fixed(), 'Step', Agrif_Nb_Step() 
    227  
    228       ! Save "old" scale factor (prior update) for subsequent asselin correction 
    229       ! of prognostic variables 
     246      !!---------------------------------------------------------------------- 
     247      !!                   *** ROUTINE dom_vvl_update_UVF *** 
     248      !!---------------------------------------------------------------------- 
     249      INTEGER ::   jk      ! dummy loop index 
     250      REAL(wp)::   zcoef   ! local scalar 
     251      !!---------------------------------------------------------------------- 
     252      ! 
     253      IF (lwp.AND.lk_agrif_debug)   Write(*,*) 'Finalize e3 on grid Number', & 
     254         &                                      Agrif_Fixed(), 'Step', Agrif_Nb_Step() 
     255 
     256      ! Save "old" scale factor (prior update) for subsequent asselin correction of prognostic variables 
    230257      ! ----------------------- 
    231       ! 
    232258      e3u_a(:,:,:) = e3u_n(:,:,:) 
    233259      e3v_a(:,:,:) = e3v_n(:,:,:) 
     
    239265      ! 1) NOW fields 
    240266      !-------------- 
    241        
    242          ! Vertical scale factor interpolations 
    243          ! ------------------------------------ 
    244       CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:) ,  'U' ) 
    245       CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:) ,  'V' ) 
    246       CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:) ,  'F' ) 
    247  
     267      !                       ! Vertical scale factor interpolations 
     268      CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n (:,:,:) , 'U' ) 
     269      CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n (:,:,:) , 'V' ) 
     270      CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n (:,:,:) , 'F' ) 
    248271      CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' ) 
    249272      CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' ) 
    250  
    251          ! Update total depths: 
    252          ! -------------------- 
     273      ! 
     274      !                       ! Update total depths 
    253275      hu_n(:,:) = 0._wp                        ! Ocean depth at U-points 
    254276      hv_n(:,:) = 0._wp                        ! Ocean depth at V-points 
     
    264286      ! 2) BEFORE fields: 
    265287      !------------------ 
    266       IF (.NOT.(lk_agrif_fstep.AND.(neuler==0) )) THEN 
    267          ! 
    268          ! Vertical scale factor interpolations 
    269          ! ------------------------------------ 
    270          CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:),  'U'  ) 
    271          CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:),  'V'  ) 
    272  
     288      IF (.NOT.( lk_agrif_fstep .AND. l_1st_euler ) ) THEN 
     289!!gm      IF (.NOT.(lk_agrif_fstep.AND.(neuler==0) )) THEN 
     290         !                    ! Vertical scale factor interpolations 
     291         CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b (:,:,:), 'U'  ) 
     292         CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b (:,:,:), 'V'  ) 
    273293         CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' ) 
    274294         CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' ) 
    275  
    276          ! Update total depths: 
    277          ! -------------------- 
     295         ! 
     296         !                    ! Update total depths: 
    278297         hu_b(:,:) = 0._wp                     ! Ocean depth at U-points 
    279298         hv_b(:,:) = 0._wp                     ! Ocean depth at V-points 
     
    289308   END SUBROUTINE dom_vvl_update_UVF 
    290309 
    291 #if defined key_vertical 
     310# if defined key_vertical 
    292311 
    293312   SUBROUTINE updateTS( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 
    294313      !!---------------------------------------------------------------------- 
    295       !!           *** ROUTINE updateT *** 
    296       !!--------------------------------------------- 
     314      !!                   *** ROUTINE updateT *** 
     315      !!---------------------------------------------------------------------- 
    297316      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2 
    298317      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 
     
    306325      REAL(wp) :: zrho_xy, h_diff 
    307326      REAL(wp) :: tabin(k1:k2,n1:n2) 
    308       !!--------------------------------------------- 
     327      !!---------------------------------------------------------------------- 
    309328      ! 
    310329      IF (before) THEN 
    311330         AGRIF_SpecialValue = -999._wp 
    312331         zrho_xy = Agrif_rhox() * Agrif_rhoy()  
    313          DO jn = n1,n2-1 
    314             DO jk=k1,k2 
    315                DO jj=j1,j2 
    316                   DO ji=i1,i2 
     332         DO jn = n1, n2-1 
     333            DO jk = k1, k2 
     334               DO jj = j1, j2 
     335                  DO ji = i1, i2 
    317336                     tabres(ji,jj,jk,jn) = (tsn(ji,jj,jk,jn) * e3t_n(ji,jj,jk) ) & 
    318337                                           * tmask(ji,jj,jk) + (tmask(ji,jj,jk)-1)*999._wp 
     
    321340            END DO 
    322341         END DO 
    323          DO jk=k1,k2 
    324             DO jj=j1,j2 
    325                DO ji=i1,i2 
     342         DO jk = k1, k2 
     343            DO jj = j1, j2 
     344               DO ji = i1, i2 
    326345                  tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk) & 
    327346                                           + (tmask(ji,jj,jk)-1)*999._wp 
     
    332351         tabres_child(:,:,:,:) = 0. 
    333352         AGRIF_SpecialValue = 0._wp 
    334          DO jj=j1,j2 
    335             DO ji=i1,i2 
     353         DO jj = j1 , j2 
     354            DO ji = i1, i2 
    336355               N_in = 0 
    337                DO jk=k1,k2 !k2 = jpk of child grid 
    338                   IF (tabres(ji,jj,jk,n2) == 0  ) EXIT 
     356               DO jk = k1, k2 !k2 = jpk of child grid 
     357                  IF ( tabres(ji,jj,jk,n2) == 0  )  EXIT 
    339358                  N_in = N_in + 1 
    340359                  tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2) 
    341                   h_in(N_in) = tabres(ji,jj,jk,n2) 
    342                ENDDO 
     360                  h_in (N_in) = tabres(ji,jj,jk,n2) 
     361               END DO 
    343362               N_out = 0 
    344                DO jk=1,jpk ! jpk of parent grid 
    345                   IF (tmask(ji,jj,jk) < -900) EXIT ! TODO: Will not work with ISF 
     363               DO jk = 1, jpk ! jpk of parent grid 
     364                  IF (tmask(ji,jj,jk) < -900)   EXIT ! TODO: Will not work with ISF 
    346365                  N_out = N_out + 1 
    347366                  h_out(N_out) = e3t_n(ji,jj,jk)  
    348                ENDDO 
     367               END DO 
    349368               IF (N_in > 0) THEN !Remove this? 
    350369                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
     
    355374                     STOP 
    356375                  ENDIF 
    357                   DO jn=n1,n2-1 
     376                  DO jn = n1, n2-1 
    358377                     CALL reconstructandremap(tabin(1:N_in,jn),h_in(1:N_in),tabres_child(ji,jj,1:N_out,jn),h_out(1:N_out),N_in,N_out) 
    359                   ENDDO 
     378                  END DO 
    360379               ENDIF 
    361             ENDDO 
    362          ENDDO 
    363  
    364          IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 
    365             ! Add asselin part 
    366             DO jn = n1,n2-1 
    367                DO jk=1,jpk 
    368                   DO jj=j1,j2 
    369                      DO ji=i1,i2 
    370                         IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN 
    371                            tsb(ji,jj,jk,jn) = tsb(ji,jj,jk,jn) &  
    372                                  & + atfp * ( tabres_child(ji,jj,jk,jn) & 
    373                                  &          - tsn(ji,jj,jk,jn) ) * tmask(ji,jj,jk) 
     380            END DO 
     381         END DO 
     382 
     383         IF ( .NOT.( lk_agrif_fstep .AND. l_1st_euler ) ) THEN       ! Add asselin part 
     384 
     385!!gm         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 
     386            DO jn = n1, n2-1 
     387               DO jk = 1, jpk 
     388                  DO jj = j1, j2 
     389                     DO ji = i1, i2 
     390                        IF( tabres_child(ji,jj,jk,jn) /= 0. ) THEN 
     391                           tsb(ji,jj,jk,jn) = tsb(ji,jj,jk,jn)   &  
     392                              &             + rn_atfp * ( tabres_child(ji,jj,jk,jn) - tsn(ji,jj,jk,jn) ) * tmask(ji,jj,jk) 
    374393                        ENDIF 
    375                      ENDDO 
    376                   ENDDO 
    377                ENDDO 
    378             ENDDO 
    379          ENDIF 
    380          DO jn = n1,n2-1 
    381             DO jk=1,jpk 
    382                DO jj=j1,j2 
    383                   DO ji=i1,i2 
    384                      IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN  
     394                     END DO 
     395                  END DO 
     396               END DO 
     397            END DO 
     398         ENDIF 
     399         DO jn = n1, n2-1 
     400            DO jk = 1, jpk 
     401               DO jj = j1, j2 
     402                  DO ji = i1, i2 
     403                     IF( tabres_child(ji,jj,jk,jn) /= 0. ) THEN  
    385404                        tsn(ji,jj,jk,jn) = tabres_child(ji,jj,jk,jn) * tmask(ji,jj,jk) 
    386405                     END IF 
     
    396415 
    397416   SUBROUTINE updateTS( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 
    398       !!--------------------------------------------- 
    399       !!           *** ROUTINE updateT *** 
    400       !!--------------------------------------------- 
     417      !!---------------------------------------------------------------------- 
     418      !!                   *** ROUTINE ROUTINE updateT *** 
     419      !!---------------------------------------------------------------------- 
    401420      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2 
    402421      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 
    403422      LOGICAL, INTENT(in) :: before 
    404       !! 
     423      ! 
    405424      INTEGER :: ji,jj,jk,jn 
    406425      REAL(wp) :: ztb, ztnu, ztno 
    407       !!--------------------------------------------- 
     426      !!---------------------------------------------------------------------- 
    408427      ! 
    409428      IF (before) THEN 
     
    425444            tabres(i1:i2,j1:j2,k1:k2,jn) =  tabres(i1:i2,j1:j2,k1:k2,jn) * e3t_0(i1:i2,j1:j2,k1:k2) & 
    426445                                         & * tmask(i1:i2,j1:j2,k1:k2) 
    427          ENDDO 
     446         END DO 
    428447!< jc tmp 
    429          IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 
     448         IF ( .NOT.( lk_agrif_fstep .AND. l_1st_euler ) ) THEN 
     449!!gm         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 
    430450            ! Add asselin part 
    431451            DO jn = 1,jpts 
     
    437457                           ztnu = tabres(ji,jj,jk,jn) 
    438458                           ztno = tsn(ji,jj,jk,jn) * e3t_a(ji,jj,jk) 
    439                            tsb(ji,jj,jk,jn) = ( ztb + atfp * ( ztnu - ztno) )  &  
    440                                      &        * tmask(ji,jj,jk) / e3t_b(ji,jj,jk) 
     459                           tsb(ji,jj,jk,jn) = ( ztb + rn_rn_atfp * ( ztnu - ztno) ) / e3t_b(ji,jj,jk) * tmask(ji,jj,jk) 
    441460                        ENDIF 
    442461                     END DO 
     
    457476         END DO 
    458477         ! 
    459          IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 
     478         IF ( l_1st_euler .AND. Agrif_Nb_Step() == 0 ) THEN 
     479!!gm         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 
    460480            tsb(i1:i2,j1:j2,k1:k2,1:jpts)  = tsn(i1:i2,j1:j2,k1:k2,1:jpts) 
    461481         ENDIF 
     
    470490 
    471491   SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 
    472       !!--------------------------------------------- 
    473       !!           *** ROUTINE updateu *** 
    474       !!--------------------------------------------- 
     492      !!---------------------------------------------------------------------- 
     493      !!                   *** ROUTINE updateu *** 
     494      !!---------------------------------------------------------------------- 
    475495      INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2 
    476496      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 
     
    487507      REAL(wp) :: tabin(k1:k2) 
    488508! VERTICAL REFINEMENT END 
    489       !!--------------------------------------------- 
     509      !!---------------------------------------------------------------------- 
    490510      !  
    491511      IF( before ) THEN 
     
    515535                  tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2) 
    516536                  h_in(N_in) = tabres(ji,jj,jk,2)/e2u(ji,jj) 
    517                ENDDO 
     537               END DO 
    518538               N_out = 0 
    519539               DO jk=1,jpk 
     
    521541                  N_out = N_out + 1 
    522542                  h_out(N_out) = e3u_n(ji,jj,jk) 
    523                ENDDO 
     543               END DO 
    524544               IF (N_in * N_out > 0) THEN 
    525545                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
     
    538558                           EXIT 
    539559                        ENDIF 
    540                      ENDDO 
     560                     END DO 
    541561                  ENDIF 
    542562                  CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out) 
    543563                  tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e2u(ji,jj)*h_out(N_out)) 
    544564               ENDIF 
    545             ENDDO 
    546          ENDDO 
    547  
    548          DO jk=1,jpk 
    549             DO jj=j1,j2 
    550                DO ji=i1,i2 
    551                   IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 
    552                      ub(ji,jj,jk) = ub(ji,jj,jk) &  
    553                            & + atfp * ( tabres_child(ji,jj,jk) - un(ji,jj,jk) ) * umask(ji,jj,jk) 
     565            END DO 
     566         END DO 
     567 
     568         DO jk = 1, jpk 
     569            DO jj = j1, j2 
     570               DO ji = i1, i2 
     571                  IF ( .NOT.( lk_agrif_fstep .AND. l_1st_euler) ) THEN    ! Add asselin part 
     572!!gm                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN   ! Add asselin part 
     573                     ub(ji,jj,jk) = ub(ji,jj,jk) + rn_atfp * ( tabres_child(ji,jj,jk) - un(ji,jj,jk) ) * umask(ji,jj,jk) 
    554574                  ENDIF 
    555                   ! 
    556575                  un(ji,jj,jk) = tabres_child(ji,jj,jk) * umask(ji,jj,jk) 
    557576               END DO 
     
    565584 
    566585   SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 
    567       !!--------------------------------------------- 
    568       !!           *** ROUTINE updateu *** 
    569       !!--------------------------------------------- 
     586      !!---------------------------------------------------------------------- 
     587      !!                   *** ROUTINE updateu *** 
     588      !!---------------------------------------------------------------------- 
    570589      INTEGER                               , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2 
    571590      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 
     
    574593      INTEGER  :: ji, jj, jk 
    575594      REAL(wp) :: zrhoy, zub, zunu, zuno 
    576       !!--------------------------------------------- 
     595      !!---------------------------------------------------------------------- 
    577596      !  
    578597      IF( before ) THEN 
     
    587606                  tabres(ji,jj,jk,1) = tabres(ji,jj,jk,1) * r1_e2u(ji,jj)  
    588607                  ! 
    589                   IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 
     608                  IF ( .NOT.( lk_agrif_fstep .AND. l_1st_euler ) ) THEN    ! Add asselin part 
     609!!gm                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN   ! Add asselin part 
    590610                     zub  = ub(ji,jj,jk) * e3u_b(ji,jj,jk)  ! fse3t_b prior update should be used 
    591611                     zuno = un(ji,jj,jk) * e3u_a(ji,jj,jk) 
    592612                     zunu = tabres(ji,jj,jk,1) 
    593                      ub(ji,jj,jk) = ( zub + atfp * ( zunu - zuno) ) &       
    594                                     & * umask(ji,jj,jk) / e3u_b(ji,jj,jk) 
     613                     ub(ji,jj,jk) = ( zub + rn_atfp * ( zunu - zuno) ) / e3u_b(ji,jj,jk) * umask(ji,jj,jk) 
    595614                  ENDIF 
    596615                  ! 
    597                   un(ji,jj,jk) = tabres(ji,jj,jk,1) * umask(ji,jj,jk) / e3u_n(ji,jj,jk) 
    598                END DO 
    599             END DO 
    600          END DO 
    601          ! 
    602          IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 
     616                  un(ji,jj,jk) = tabres(ji,jj,jk,1) / e3u_n(ji,jj,jk) * umask(ji,jj,jk) 
     617               END DO 
     618            END DO 
     619         END DO 
     620         ! 
     621         IF ( l_1st_euler .AND. Agrif_Nb_Step() == 0 ) THEN 
     622!!gm         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 
    603623            ub(i1:i2,j1:j2,k1:k2)  = un(i1:i2,j1:j2,k1:k2) 
    604624         ENDIF 
     
    611631 
    612632   SUBROUTINE correct_u_bdy( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir ) 
    613       !!--------------------------------------------- 
    614       !!           *** ROUTINE correct_u_bdy *** 
    615       !!--------------------------------------------- 
     633      !!---------------------------------------------------------------------- 
     634      !!                   *** ROUTINE correct_u_bdy *** 
     635      !!---------------------------------------------------------------------- 
    616636      INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2 
    617637      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 
    618638      LOGICAL                                     , INTENT(in   ) :: before 
    619       INTEGER                                     , INTENT(in)    :: nb, ndir 
     639      INTEGER                                     , INTENT(in   ) :: nb, ndir 
    620640      !! 
    621641      LOGICAL :: western_side, eastern_side  
    622       ! 
    623       INTEGER  :: jj, jk 
    624       REAL(wp) :: zcor 
    625       !!--------------------------------------------- 
     642      INTEGER ::   jj, jk 
     643      REAL(wp)::   zcor 
     644      !!---------------------------------------------------------------------- 
    626645      !  
    627646      IF( .NOT.before ) THEN 
     
    657676 
    658677   SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 
    659       !!--------------------------------------------- 
    660       !!           *** ROUTINE updatev *** 
    661       !!--------------------------------------------- 
     678      !!---------------------------------------------------------------------- 
     679      !!                   *** ROUTINE updatev *** 
     680      !!---------------------------------------------------------------------- 
    662681      INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2 
    663682      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 
     
    674693      REAL(wp) :: tabin(k1:k2) 
    675694! VERTICAL REFINEMENT END 
    676       !!---------------------------------------------       
     695      !!----------------------------------------------------------------------       
    677696      ! 
    678697      IF( before ) THEN 
     
    700719                  tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2) 
    701720                  h_in(N_in) = tabres(ji,jj,jk,2)/e1v(ji,jj) 
    702                ENDDO 
     721               END DO 
    703722               N_out = 0 
    704723               DO jk=1,jpk 
     
    706725                  N_out = N_out + 1 
    707726                  h_out(N_out) = e3v_n(ji,jj,jk) 
    708                ENDDO 
     727               END DO 
    709728               IF (N_in * N_out > 0) THEN 
    710729                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
     
    723742                           EXIT 
    724743                        ENDIF 
    725                      ENDDO 
     744                     END DO 
    726745                  ENDIF 
    727746                  CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out) 
    728747                  tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e1v(ji,jj)*h_out(N_out)) 
    729748               ENDIF 
    730             ENDDO 
    731          ENDDO 
     749            END DO 
     750         END DO 
    732751 
    733752         DO jk=1,jpk 
     
    735754               DO ji=i1,i2 
    736755                  ! 
    737                   IF( .NOT.(lk_agrif_fstep.AND.(neuler==0)) ) THEN ! Add asselin part 
    738                      vb(ji,jj,jk) = vb(ji,jj,jk) &  
    739                            & + atfp * ( tabres_child(ji,jj,jk) - vn(ji,jj,jk) ) * vmask(ji,jj,jk) 
     756                  IF( .NOT.( lk_agrif_fstep .AND. l_1st_euler ) ) THEN ! Add asselin part 
     757!!gm                  IF( .NOT.(lk_agrif_fstep.AND.(neuler==0)) ) THEN  ! Add asselin part 
     758                     vb(ji,jj,jk) = vb(ji,jj,jk) + rn_atfp * ( tabres_child(ji,jj,jk) - vn(ji,jj,jk) ) * vmask(ji,jj,jk) 
    740759                  ENDIF 
    741760                  ! 
     
    751770 
    752771   SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before) 
    753       !!--------------------------------------------- 
    754       !!           *** ROUTINE updatev *** 
    755       !!--------------------------------------------- 
     772      !!---------------------------------------------------------------------- 
     773      !!                   *** ROUTINE updatev *** 
     774      !!---------------------------------------------------------------------- 
    756775      INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2 
    757776      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 
     
    760779      INTEGER  :: ji, jj, jk 
    761780      REAL(wp) :: zrhox, zvb, zvnu, zvno 
    762       !!---------------------------------------------       
     781      !!----------------------------------------------------------------------       
    763782      ! 
    764783      IF (before) THEN 
     
    777796                  tabres(ji,jj,jk,1) = tabres(ji,jj,jk,1) * r1_e1v(ji,jj) 
    778797                  ! 
    779                   IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 
     798                  IF ( .NOT.( lk_agrif_fstep .AND. l_1st_euler ) ) THEN    ! Add asselin part 
     799!!gm                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN   ! Add asselin part 
    780800                     zvb  = vb(ji,jj,jk) * e3v_b(ji,jj,jk) ! fse3t_b prior update should be used 
    781801                     zvno = vn(ji,jj,jk) * e3v_a(ji,jj,jk) 
    782802                     zvnu = tabres(ji,jj,jk,1) 
    783                      vb(ji,jj,jk) = ( zvb + atfp * ( zvnu - zvno) ) &       
    784                                     & * vmask(ji,jj,jk) / e3v_b(ji,jj,jk) 
     803                     vb(ji,jj,jk) = ( zvb + rn_atfp * ( zvnu - zvno) ) / e3v_b(ji,jj,jk) * vmask(ji,jj,jk) 
    785804                  ENDIF 
    786805                  ! 
    787                   vn(ji,jj,jk) = tabres(ji,jj,jk,1) * vmask(ji,jj,jk) / e3v_n(ji,jj,jk) 
    788                END DO 
    789             END DO 
    790          END DO 
    791          ! 
    792          IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 
     806                  vn(ji,jj,jk) = tabres(ji,jj,jk,1) / e3v_n(ji,jj,jk) * vmask(ji,jj,jk) 
     807               END DO 
     808            END DO 
     809         END DO 
     810         ! 
     811         IF ( l_1st_euler .AND. Agrif_Nb_Step() == 0 ) THEN 
     812!!gm         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 
    793813            vb(i1:i2,j1:j2,k1:k2)  = vn(i1:i2,j1:j2,k1:k2) 
    794814         ENDIF 
     
    801821 
    802822   SUBROUTINE correct_v_bdy( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir ) 
    803       !!--------------------------------------------- 
    804       !!           *** ROUTINE correct_u_bdy *** 
    805       !!--------------------------------------------- 
     823      !!---------------------------------------------------------------------- 
     824      !!                   *** ROUTINE correct_v_bdy *** 
     825      !!---------------------------------------------------------------------- 
    806826      INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2 
    807827      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 
     
    813833      INTEGER  :: ji, jk 
    814834      REAL(wp) :: zcor 
    815       !!--------------------------------------------- 
     835      !!---------------------------------------------------------------------- 
    816836      !  
    817837      IF( .NOT.before ) THEN 
     
    847867   SUBROUTINE updateu2d( tabres, i1, i2, j1, j2, before ) 
    848868      !!---------------------------------------------------------------------- 
    849       !!                      *** ROUTINE updateu2d *** 
     869      !!                   *** ROUTINE updateu2d *** 
    850870      !!---------------------------------------------------------------------- 
    851871      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2 
    852872      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   tabres 
    853873      LOGICAL                         , INTENT(in   ) ::   before 
    854       !!  
     874      ! 
    855875      INTEGER  :: ji, jj, jk 
    856876      REAL(wp) :: zrhoy 
    857877      REAL(wp) :: zcorr 
    858       !!--------------------------------------------- 
     878      !!---------------------------------------------------------------------- 
    859879      ! 
    860880      IF( before ) THEN 
     
    883903               ! Update barotropic velocities: 
    884904               IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN 
    885                   IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 
     905                  IF ( .NOT.( lk_agrif_fstep .AND. l_1st_euler ) ) THEN    ! Add asselin part 
     906!!gm                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 
    886907                     zcorr = (tabres(ji,jj) - un_b(ji,jj) * hu_a(ji,jj)) * r1_hu_b(ji,jj) 
    887                      ub_b(ji,jj) = ub_b(ji,jj) + atfp * zcorr * umask(ji,jj,1) 
     908                     ub_b(ji,jj) = ub_b(ji,jj) + rn_atfp * zcorr * umask(ji,jj,1) 
    888909                  END IF 
    889910               ENDIF     
     
    904925         END DO 
    905926         ! 
    906          IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 
     927         IF ( l_1st_euler .AND. Agrif_Nb_Step() == 0 ) THEN 
     928!!gm         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 
    907929            ub_b(i1:i2,j1:j2)  = un_b(i1:i2,j1:j2) 
    908930         ENDIF 
     
    948970               ! 
    949971               ! Update barotropic velocities: 
    950                IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN 
    951                   IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 
     972               IF ( .NOT.ln_dynspg_ts .OR. ( ln_dynspg_ts .AND. .NOT.ln_bt_fw ) ) THEN 
     973                  IF ( .NOT.( lk_agrif_fstep. AND. l_1st_euler ) ) THEN    ! Add asselin part 
     974!!gm                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 
    952975                     zcorr = (tabres(ji,jj) - vn_b(ji,jj) * hv_a(ji,jj)) * r1_hv_b(ji,jj) 
    953                      vb_b(ji,jj) = vb_b(ji,jj) + atfp * zcorr * vmask(ji,jj,1) 
     976                     vb_b(ji,jj) = vb_b(ji,jj) + rn_atfp * zcorr * vmask(ji,jj,1) 
    954977                  END IF 
    955978               ENDIF               
     
    970993         END DO 
    971994         ! 
    972          IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 
     995         IF ( l_1st_euler .AND. Agrif_Nb_Step() == 0 ) THEN 
     996!!gm         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 
    973997            vb_b(i1:i2,j1:j2)  = vn_b(i1:i2,j1:j2) 
    974998         ENDIF 
     
    9861010      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   tabres 
    9871011      LOGICAL                         , INTENT(in   ) ::   before 
    988       !! 
     1012      ! 
    9891013      INTEGER :: ji, jj 
    9901014      !!---------------------------------------------------------------------- 
     
    9971021         END DO 
    9981022      ELSE 
    999          IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 
     1023         IF ( .NOT.( lk_agrif_fstep .AND. l_1st_euler ) ) THEN 
     1024!!gm         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 
    10001025            DO jj=j1,j2 
    10011026               DO ji=i1,i2 
    1002                   sshb(ji,jj) =   sshb(ji,jj) & 
    1003                         & + atfp * ( tabres(ji,jj) - sshn(ji,jj) ) * tmask(ji,jj,1) 
     1027                  sshb(ji,jj) = sshb(ji,jj) + rn_atfp * ( tabres(ji,jj) - sshn(ji,jj) ) * tmask(ji,jj,1) 
    10041028               END DO 
    10051029            END DO 
     
    10121036         END DO 
    10131037         ! 
    1014          IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 
     1038         IF ( l_1st_euler .AND. Agrif_Nb_Step() == 0 ) THEN 
     1039!!gm         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 
    10151040            sshb(i1:i2,j1:j2)  = sshn(i1:i2,j1:j2) 
    10161041         ENDIF 
    10171042         ! 
    1018  
    10191043      ENDIF 
    10201044      ! 
     
    10621086   END SUBROUTINE updateub2b 
    10631087 
     1088 
    10641089   SUBROUTINE reflux_sshu( tabres, i1, i2, j1, j2, before, nb, ndir ) 
    1065       !!--------------------------------------------- 
    1066       !!          *** ROUTINE reflux_sshu *** 
    1067       !!--------------------------------------------- 
    1068       INTEGER, INTENT(in) :: i1, i2, j1, j2 
    1069       REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres 
    1070       LOGICAL, INTENT(in) :: before 
    1071       INTEGER, INTENT(in) :: nb, ndir 
    1072       !! 
    1073       LOGICAL :: western_side, eastern_side  
    1074       INTEGER :: ji, jj 
    1075       REAL(wp) :: zrhoy, za1, zcor 
    1076       !!--------------------------------------------- 
     1090      !!---------------------------------------------------------------------- 
     1091      !!                   *** ROUTINE reflux_sshu *** 
     1092      !!---------------------------------------------------------------------- 
     1093      INTEGER                         , INTENT(in   ) ::  i1, i2, j1, j2 
     1094      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   tabres 
     1095      LOGICAL                         , INTENT(in   ) ::  before 
     1096      INTEGER                         , INTENT(in   ) ::  nb, ndir 
     1097      ! 
     1098      LOGICAL ::   western_side, eastern_side  
     1099      INTEGER ::   ji, jj 
     1100      REAL(wp)::  zrhoy, za1, zcor 
     1101      !!---------------------------------------------------------------------- 
    10771102      ! 
    10781103      IF (before) THEN 
     
    10911116         eastern_side  = (nb == 1).AND.(ndir == 2) 
    10921117         ! 
    1093          IF (western_side) THEN 
     1118         IF ( western_side ) THEN 
    10941119            DO jj=j1,j2 
    1095                zcor = rdt * r1_e1e2t(i1  ,jj) * e2u(i1,jj) * (ub2_b(i1,jj)-tabres(i1,jj))  
     1120               zcor = rn_Dt * r1_e1e2t(i1  ,jj) * e2u(i1,jj) * (ub2_b(i1,jj)-tabres(i1,jj))  
    10961121               sshn(i1  ,jj) = sshn(i1  ,jj) + zcor 
    1097                IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(i1  ,jj) = sshb(i1  ,jj) + atfp * zcor 
    1098             END DO 
    1099          ENDIF 
    1100          IF (eastern_side) THEN 
     1122               IF ( .NOT.( lk_agrif_fstep .AND. l_1st_euler ) )   sshb(i1  ,jj) = sshb(i1  ,jj) + rn_atfp * zcor 
     1123!!gm               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(i1  ,jj) = sshb(i1  ,jj) + rn_atfp * zcor 
     1124            END DO 
     1125         ENDIF 
     1126         IF ( eastern_side ) THEN 
    11011127            DO jj=j1,j2 
    1102                zcor = - rdt * r1_e1e2t(i2+1,jj) * e2u(i2,jj) * (ub2_b(i2,jj)-tabres(i2,jj)) 
     1128               zcor = - rn_Dt * r1_e1e2t(i2+1,jj) * e2u(i2,jj) * (ub2_b(i2,jj)-tabres(i2,jj)) 
    11031129               sshn(i2+1,jj) = sshn(i2+1,jj) + zcor 
    1104                IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(i2+1,jj) = sshb(i2+1,jj) + atfp * zcor 
     1130               IF (.NOT.( lk_agrif_fstep .AND. l_1st_euler ) )   sshb(i2+1,jj) = sshb(i2+1,jj) + rn_atfp * zcor 
     1131!!gm               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(i2+1,jj) = sshb(i2+1,jj) + rn_atfp * zcor 
    11051132            END DO 
    11061133         ENDIF 
     
    11101137   END SUBROUTINE reflux_sshu 
    11111138 
     1139 
    11121140   SUBROUTINE updatevb2b( tabres, i1, i2, j1, j2, before ) 
    11131141      !!---------------------------------------------------------------------- 
    1114       !!                      *** ROUTINE updatevb2b *** 
     1142      !!                    *** ROUTINE updatevb2b *** 
    11151143      !!---------------------------------------------------------------------- 
    11161144      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2 
    11171145      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   tabres 
    11181146      LOGICAL                         , INTENT(in   ) ::   before 
    1119       !! 
     1147      ! 
    11201148      INTEGER :: ji, jj 
    11211149      REAL(wp) :: zrhox, za1, zcor 
    1122       !!--------------------------------------------- 
     1150      !!--------------------------------------------------------------------- 
    11231151      ! 
    11241152      IF( before ) THEN 
     
    11501178   END SUBROUTINE updatevb2b 
    11511179 
     1180 
    11521181   SUBROUTINE reflux_sshv( tabres, i1, i2, j1, j2, before, nb, ndir ) 
    1153       !!--------------------------------------------- 
    1154       !!          *** ROUTINE reflux_sshv *** 
    1155       !!--------------------------------------------- 
    1156       INTEGER, INTENT(in) :: i1, i2, j1, j2 
    1157       REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres 
    1158       LOGICAL, INTENT(in) :: before 
    1159       INTEGER, INTENT(in) :: nb, ndir 
     1182      !!---------------------------------------------------------------------- 
     1183      !!                   *** ROUTINE reflux_sshv *** 
     1184      !!---------------------------------------------------------------------- 
     1185      INTEGER                         , INTENT(in   ) ::  i1, i2, j1, j2 
     1186      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   tabres 
     1187      LOGICAL                         , INTENT(in   ) ::  before 
     1188      INTEGER                         , INTENT(in   ) ::  nb, ndir 
    11601189      !! 
    11611190      LOGICAL :: southern_side, northern_side  
    11621191      INTEGER :: ji, jj 
    11631192      REAL(wp) :: zrhox, za1, zcor 
    1164       !!--------------------------------------------- 
     1193      !!---------------------------------------------------------------------- 
    11651194      ! 
    11661195      IF (before) THEN 
     
    11811210         IF (southern_side) THEN 
    11821211            DO ji=i1,i2 
    1183                zcor = rdt * r1_e1e2t(ji,j1  ) * e1v(ji,j1  ) * (vb2_b(ji,j1)-tabres(ji,j1)) 
     1212               zcor = rn_Dt * r1_e1e2t(ji,j1  ) * e1v(ji,j1  ) * (vb2_b(ji,j1)-tabres(ji,j1)) 
    11841213               sshn(ji,j1  ) = sshn(ji,j1  ) + zcor 
    1185                IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(ji,j1  ) = sshb(ji,j1) + atfp * zcor 
     1214               IF ( .NOT.( lk_agrif_fstep .AND. l_euler ) )   sshb(ji,j1  ) = sshb(ji,j1) + rn_atfp * zcor 
     1215!!gm               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(ji,j1  ) = sshb(ji,j1) + rn_atfp * zcor 
    11861216            END DO 
    11871217         ENDIF 
    11881218         IF (northern_side) THEN                
    11891219            DO ji=i1,i2 
    1190                zcor = - rdt * r1_e1e2t(ji,j2+1) * e1v(ji,j2  ) * (vb2_b(ji,j2)-tabres(ji,j2)) 
     1220               zcor = - rn_Dt * r1_e1e2t(ji,j2+1) * e1v(ji,j2  ) * (vb2_b(ji,j2)-tabres(ji,j2)) 
    11911221               sshn(ji,j2+1) = sshn(ji,j2+1) + zcor 
    1192                IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(ji,j2+1) = sshb(ji,j2+1) + atfp * zcor 
     1222               IF ( .NOT.( lk_agrif_fstep .AND. l_1st_euler ) )   sshb(ji,j2+1) = sshb(ji,j2+1) + rn_atfp * zcor 
     1223!!gm               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(ji,j2+1) = sshb(ji,j2+1) + rn_atfp * zcor 
    11931224            END DO 
    11941225         ENDIF 
     
    11981229   END SUBROUTINE reflux_sshv 
    11991230 
     1231 
    12001232   SUBROUTINE update_scales( tabres, i1, i2, j1, j2, k1, k2, n1,n2, before ) 
     1233      !!---------------------------------------------------------------------- 
     1234      !!                      *** ROUTINE updateT *** 
    12011235      ! 
    12021236      ! ====>>>>>>>>>>    currently not used 
    12031237      ! 
    1204       !!---------------------------------------------------------------------- 
    1205       !!                      *** ROUTINE updateT *** 
    12061238      !!---------------------------------------------------------------------- 
    12071239      INTEGER                                    , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2, n1, n2 
     
    12841316 
    12851317   SUBROUTINE updateAVM( ptab, i1, i2, j1, j2, k1, k2, before ) 
    1286       !!--------------------------------------------- 
    1287       !!           *** ROUTINE updateavm *** 
     1318      !!---------------------------------------------------------------------- 
     1319      !!                      *** ROUTINE updateavm *** 
    12881320      !!---------------------------------------------------------------------- 
    12891321      INTEGER                               , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2 
     
    12981330   END SUBROUTINE updateAVM 
    12991331 
     1332 
    13001333   SUBROUTINE updatee3t(ptab_dum, i1, i2, j1, j2, k1, k2, before ) 
    1301       !!--------------------------------------------- 
    1302       !!           *** ROUTINE updatee3t *** 
    1303       !!--------------------------------------------- 
     1334      !!---------------------------------------------------------------------- 
     1335      !!                   *** ROUTINE updatee3t *** 
     1336      !!---------------------------------------------------------------------- 
    13041337      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: ptab_dum 
    13051338      INTEGER, INTENT(in) :: i1, i2, j1, j2, k1, k2 
     
    13131346      IF (.NOT.before) THEN 
    13141347         ! 
    1315          ALLOCATE(ptab(i1:i2,j1:j2,1:jpk))  
     1348         ALLOCATE( ptab(i1:i2,j1:j2,1:jpk) )  
    13161349         ! 
    13171350         ! Update e3t from ssh (z* case only) 
     
    13351368!         hdivn(i1:i2,j1:j2,1:jpkm1)   = e3t_b(i1:i2,j1:j2,1:jpkm1) 
    13361369 
    1337          IF (.NOT.(lk_agrif_fstep.AND.(neuler==0) )) THEN 
     1370         IF ( .NOT.( lk_agrif_fstep .AND. l_1st_euler==0 ) ) THEN 
     1371!!gm         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0) )) THEN 
    13381372            DO jk = 1, jpkm1 
    13391373               DO jj=j1,j2 
    13401374                  DO ji=i1,i2 
    1341                      e3t_b(ji,jj,jk) =  e3t_b(ji,jj,jk) & 
    1342                            & + atfp * ( ptab(ji,jj,jk) - e3t_n(ji,jj,jk) ) 
     1375                     e3t_b(ji,jj,jk) =  e3t_b(ji,jj,jk) + rn_atfp * ( ptab(ji,jj,jk) - e3t_n(ji,jj,jk) ) 
    13431376                  END DO 
    13441377               END DO 
     
    13981431         END DO 
    13991432         ! 
    1400          IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 
     1433         IF ( l_1st_euler .AND. Agrif_Nb_Step()==0 ) THEN 
     1434!!gm         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 
    14011435            e3t_b (i1:i2,j1:j2,1:jpk)  = e3t_n (i1:i2,j1:j2,1:jpk) 
    14021436            e3w_b (i1:i2,j1:j2,1:jpk)  = e3w_n (i1:i2,j1:j2,1:jpk) 
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/NST/agrif_top_update.F90

    r9598 r9939  
    109109                  tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2) 
    110110                  h_in(N_in) = tabres(ji,jj,jk,n2) 
    111                ENDDO 
     111               END DO 
    112112               N_out = 0 
    113113               DO jk=1,jpk ! jpk of parent grid 
     
    115115                  N_out = N_out + 1 
    116116                  h_out(N_out) = e3t_n(ji,jj,jk) !Parent grid scale factors. Could multiply by e1e2t here instead of division above 
    117                ENDDO 
     117               END DO 
    118118               IF (N_in > 0) THEN !Remove this? 
    119119                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
     
    126126                  DO jn=1,jptra 
    127127                     CALL reconstructandremap(tabin(1:N_in,jn),h_in(1:N_in),tabres_child(ji,jj,1:N_out,jn),h_out(1:N_out),N_in,N_out) 
    128                   ENDDO 
     128                  END DO 
    129129               ENDIF 
    130             ENDDO 
    131          ENDDO 
    132  
    133          IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 
     130            END DO 
     131         END DO 
     132 
     133         IF ( .NOT.( lk_agrif_fstep .AND. l_1st_euler ) ) THEN 
     134!!gm         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 
    134135            ! Add asselin part 
    135136            DO jn = 1,jptra 
     
    139140                        IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN 
    140141                           trb(ji,jj,jk,jn) = tsb(ji,jj,jk,jn) &  
    141                                  & + atfp * ( tabres_child(ji,jj,jk,jn) & 
    142                                  &          - trn(ji,jj,jk,jn) ) * tmask(ji,jj,jk) 
     142                                 & + rn_atfp * ( tabres_child(ji,jj,jk,jn) - trn(ji,jj,jk,jn) ) * tmask(ji,jj,jk) 
    143143                        ENDIF 
    144                      ENDDO 
    145                   ENDDO 
    146                ENDDO 
    147             ENDDO 
     144                     END DO 
     145                  END DO 
     146               END DO 
     147            END DO 
    148148         ENDIF 
    149149         DO jn = 1,jptra 
     
    195195            tabres(i1:i2,j1:j2,k1:k2,jn) =  tabres(i1:i2,j1:j2,k1:k2,jn) * e3t_0(i1:i2,j1:j2,k1:k2) & 
    196196                                         & * tmask(i1:i2,j1:j2,k1:k2) 
    197          ENDDO 
     197         END DO 
    198198!< jc tmp 
    199          IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 
     199         IF (.NOT.( lk_agrif_fstep .AND. l_1st_euler ) ) THEN 
     200!!gm         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 
    200201            ! Add asselin part 
    201202            DO jn = n1,n2 
     
    207208                           ztnu = tabres(ji,jj,jk,jn) 
    208209                           ztno = trn(ji,jj,jk,jn) * e3t_a(ji,jj,jk) 
    209                            trb(ji,jj,jk,jn) = ( ztb + atfp * ( ztnu - ztno) )  &  
    210                                      &        * tmask(ji,jj,jk) / e3t_b(ji,jj,jk) 
     210                           trb(ji,jj,jk,jn) = ( ztb + rn_atfp * ( ztnu - ztno) ) / e3t_b(ji,jj,jk) * tmask(ji,jj,jk) 
    211211                        ENDIF 
    212                      ENDDO 
    213                   ENDDO 
    214                ENDDO 
    215             ENDDO 
     212                     END DO 
     213                  END DO 
     214               END DO 
     215            END DO 
    216216         ENDIF 
    217217         DO jn = n1,n2 
     
    227227         END DO 
    228228         ! 
    229          IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 
     229         IF ( l_1st_euler .AND. Agrif_Nb_Step() == 0 ) THEN 
     230!!gm         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 
    230231            trb(i1:i2,j1:j2,k1:k2,n1:n2)  = trn(i1:i2,j1:j2,k1:k2,n1:n2) 
    231232         ENDIF 
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/NST/agrif_user.F90

    r9788 r9939  
    217217 
    218218      ! Check time steps            
    219       IF( NINT(Agrif_Rhot()) * NINT(rdt) .NE. Agrif_Parent(rdt) ) THEN 
    220          WRITE(cl_check1,*)  NINT(Agrif_Parent(rdt)) 
    221          WRITE(cl_check2,*)  NINT(rdt) 
    222          WRITE(cl_check3,*)  NINT(Agrif_Parent(rdt)/Agrif_Rhot()) 
     219      IF( NINT(Agrif_Rhot()) * NINT(rn_Dt) /= Agrif_Parent(rn_Dt) ) THEN 
     220         WRITE(cl_check1,*)  NINT(Agrif_Parent(rn_Dt)) 
     221         WRITE(cl_check2,*)  NINT(rn_Dt) 
     222         WRITE(cl_check3,*)  NINT(Agrif_Parent(rn_Dt)/Agrif_Rhot()) 
    223223         CALL ctl_stop( 'Incompatible time step between ocean grids',   & 
    224224               &               'parent grid value : '//cl_check1    ,   &  
     
    229229      ! Check run length 
    230230      IF( Agrif_IRhot() * (Agrif_Parent(nitend)- & 
    231             Agrif_Parent(nit000)+1) .NE. (nitend-nit000+1) ) THEN 
     231            Agrif_Parent(nit000)+1) /= (nitend-nit000+1) ) THEN 
    232232         WRITE(cl_check1,*)  (Agrif_Parent(nit000)-1)*Agrif_IRhot() + 1 
    233233         WRITE(cl_check2,*)   Agrif_Parent(nitend)   *Agrif_IRhot() 
     
    601601   IF( check_namelist ) THEN 
    602602      ! Check time steps 
    603       IF( NINT(Agrif_Rhot()) * NINT(rdt) .NE. Agrif_Parent(rdt) ) THEN 
    604          WRITE(cl_check1,*)  Agrif_Parent(rdt) 
    605          WRITE(cl_check2,*)  rdt 
    606          WRITE(cl_check3,*)  rdt*Agrif_Rhot() 
     603      IF( NINT(Agrif_Rhot()) * NINT(rn_Dt) .NE. Agrif_Parent(rn_Dt) ) THEN 
     604         WRITE(cl_check1,*)  Agrif_Parent(rn_Dt) 
     605         WRITE(cl_check2,*)  rn_Dt 
     606         WRITE(cl_check3,*)  rn_Dt*Agrif_Rhot() 
    607607         CALL ctl_stop( 'incompatible time step between grids',   & 
    608608               &               'parent grid value : '//cl_check1    ,   &  
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/ASM/asminc.F90

    r9656 r9939  
    491491      ENDIF 
    492492      ! 
    493       IF(lwp) WRITE(numout,*) '   ==>>>   Euler time step switch is ', neuler 
     493      IF(lwp) WRITE(numout,*) '   ==>>>   Euler time step switch is ', ln_1st_euler 
    494494      ! 
    495495      IF( lk_asminc ) THEN                            !==  data assimilation  ==! 
     
    536536            ! 
    537537            it = kt - nit000 + 1 
    538             zincwgt = wgtiau(it) / rdt   ! IAU weight for the current time step 
     538            zincwgt = wgtiau(it) / rn_Dt   ! IAU weight for the current time step 
    539539            ! 
    540540            IF(lwp) THEN 
     
    579579         IF ( kt == nitdin_r ) THEN 
    580580            ! 
    581             neuler = 0  ! Force Euler forward step 
     581            l_1st_euler = .TRUE.  ! Force Euler forward step 
    582582            ! 
    583583            ! Initialize the now fields with the background + increment 
     
    651651            ! 
    652652            it = kt - nit000 + 1 
    653             zincwgt = wgtiau(it) / rdt   ! IAU weight for the current time step 
     653            zincwgt = wgtiau(it) / rn_Dt   ! IAU weight for the current time step 
    654654            ! 
    655655            IF(lwp) THEN 
     
    677677         IF ( kt == nitdin_r ) THEN 
    678678            ! 
    679             neuler = 0                    ! Force Euler forward step 
     679            l_1st_euler = .TRUE.          ! Force Euler forward step 
    680680            ! 
    681681            ! Initialize the now fields with the background + increment 
     
    721721            ! 
    722722            it = kt - nit000 + 1 
    723             zincwgt = wgtiau(it) / rdt   ! IAU weight for the current time step 
     723            zincwgt = wgtiau(it) / rn_Dt   ! IAU weight for the current time step 
    724724            ! 
    725725            IF(lwp) THEN 
     
    752752         IF ( kt == nitdin_r ) THEN 
    753753            ! 
    754             neuler = 0                                   ! Force Euler forward step 
     754            l_1st_euler = .TRUE.                         ! Force Euler forward step 
    755755            ! 
    756756            sshn(:,:) = ssh_bkg(:,:) + ssh_bkginc(:,:)   ! Initialize the now fields the background + increment 
     
    758758            sshb(:,:) = sshn(:,:)                        ! Update before fields 
    759759            e3t_b(:,:,:) = e3t_n(:,:,:) 
    760 !!gm why not e3u_b, e3v_b, gdept_b ???? 
     760             
     761!!gm BUG :   missing the update of all other scale factors (e3u e3v e3w  etc... _n and _b)  
     762!!           see dom_vvl_init  
    761763            ! 
    762764            DEALLOCATE( ssh_bkg    ) 
     
    839841            it = kt - nit000 + 1 
    840842            zincwgt = wgtiau(it)      ! IAU weight for the current time step  
    841             ! note this is not a tendency so should not be divided by rdt (as with the tracer and other increments) 
     843            ! note this is not a tendency so should not be divided by rn_Dt (as with the tracer and other increments) 
    842844            ! 
    843845            IF(lwp) THEN 
     
    874876#if defined key_cice && defined key_asminc 
    875877            ! Sea-ice : CICE case. Pass ice increment tendency into CICE 
    876             ndaice_da(:,:) = seaice_bkginc(:,:) * zincwgt / rdt 
     878            ndaice_da(:,:) = seaice_bkginc(:,:) * zincwgt / rn_Dt 
    877879#endif 
    878880            ! 
     
    894896         IF ( kt == nitdin_r ) THEN 
    895897            ! 
    896             neuler = 0                    ! Force Euler forward step 
     898            l_1st_euler = .TRUE.             ! Force Euler forward step 
    897899            ! 
    898900            ! Sea-ice : SI3 case 
     
    924926#if defined key_cice && defined key_asminc 
    925927            ! Sea-ice : CICE case. Pass ice increment tendency into CICE 
    926            ndaice_da(:,:) = seaice_bkginc(:,:) / rdt 
     928           ndaice_da(:,:) = seaice_bkginc(:,:) / rn_Dt 
    927929#endif 
    928930            IF ( .NOT. PRESENT(kindic) ) THEN 
     
    957959!           ! fwf : ice formation and melting 
    958960! 
    959 !                 zfons = ( -nfresh_da(ji,jj)*soce + nfsalt_da(ji,jj) )*rdt 
     961!                 zfons = ( -nfresh_da(ji,jj)*soce + nfsalt_da(ji,jj) ) * rn_Dt 
    960962! 
    961963!           ! change salinity down to mixed layer depth 
     
    10061008!      !!                                                     ! E-P (kg m-2 s-2) 
    10071009!      !            emp(ji,jj) = emp(ji,jj) + zpmess          ! E-P (kg m-2 s-2) 
    1008 !               ENDDO !ji 
    1009 !             ENDDO !jj! 
     1010!               END DO !ji 
     1011!             END DO !jj! 
    10101012! 
    10111013!            ENDIF !ln_seaicebal 
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/BDY/bdyice.F90

    r9657 r9939  
    124124 
    125125            ! Then, a) transfer the snow excess into the ice (different from icethd_dh) 
    126             zdh = MAX( 0._wp, ( rhosn * h_s(ji,jj,jl) + ( rhoic - rau0 ) * h_i(ji,jj,jl) ) * r1_rau0 ) 
     126            zdh = MAX( 0._wp, ( rhos * h_s(ji,jj,jl) + ( rhoi - rho0 ) * h_i(ji,jj,jl) ) * r1_rho0 ) 
    127127            ! Or, b) transfer all the snow into ice (if incoming ice is likely to melt as it comes into a warmer environment) 
    128             !zdh = MAX( 0._wp, h_s(ji,jj,jl) * rhosn / rhoic ) 
     128            !zdh = MAX( 0._wp, h_s(ji,jj,jl) * rhos / rhoi ) 
    129129 
    130130            ! recompute h_i, h_s 
    131131            h_i(ji,jj,jl) = MIN( hi_max(jl), h_i(ji,jj,jl) + zdh ) 
    132             h_s(ji,jj,jl) = MAX( 0._wp, h_s(ji,jj,jl) - zdh * rhoic / rhosn )  
    133  
    134          ENDDO 
     132            h_s(ji,jj,jl) = MAX( 0._wp, h_s(ji,jj,jl) - zdh * rhoi / rhos )  
     133 
     134         END DO 
    135135         CALL lbc_bdy_lnk( a_i(:,:,jl), 'T', 1., ib_bdy ) 
    136136         CALL lbc_bdy_lnk( h_i(:,:,jl), 'T', 1., ib_bdy ) 
    137137         CALL lbc_bdy_lnk( h_s(:,:,jl), 'T', 1., ib_bdy ) 
    138       ENDDO 
     138      END DO 
    139139      ! retrieve at_i 
    140140      at_i(:,:) = 0._wp 
     
    212212            DO jk = 1, nlay_s 
    213213               ! Snow energy of melting 
    214                e_s(ji,jj,jk,jl) = rswitch * rhosn * ( cpic * ( rt0 - t_s(ji,jj,jk,jl) ) + lfus ) 
     214               e_s(ji,jj,jk,jl) = rswitch * rhos * ( rcpi * ( rt0 - t_s(ji,jj,jk,jl) ) + rLfus ) 
    215215               ! Multiply by volume, so that heat content in J/m2 
    216216               e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) * r1_nlay_s 
     
    219219               ztmelts          = - tmut * sz_i(ji,jj,jk,jl) + rt0 !Melting temperature in K                   
    220220               ! heat content per unit volume 
    221                e_i(ji,jj,jk,jl) = rswitch * rhoic * & 
    222                   (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) & 
    223                   +   lfus    * ( 1.0 - (ztmelts-rt0) / MIN((t_i(ji,jj,jk,jl)-rt0),-epsi20) ) & 
    224                   - rcp      * ( ztmelts - rt0 ) ) 
     221               e_i(ji,jj,jk,jl) = rswitch * rhoi * & 
     222                  (   rcpi  * ( ztmelts - t_i(ji,jj,jk,jl) ) & 
     223                  +   rLfus * ( 1.0 - (ztmelts-rt0) / MIN((t_i(ji,jj,jk,jl)-rt0),-epsi20) ) & 
     224                  -   rcp   * ( ztmelts - rt0 ) ) 
    225225               ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2 
    226226               e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * a_i(ji,jj,jl) * h_i(ji,jj,jl) * r1_nlay_i 
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/BDY/bdylib.F90

    r9598 r9939  
    44   !! Unstructured Open Boundary Cond. :  Library module of generic boundary algorithms. 
    55   !!====================================================================== 
    6    !! History :  3.6  !  2013     (D. Storkey) original code 
    7    !!            4.0  !  2014     (T. Lovato) Generalize OBC structure 
     6   !! History :  3.6  !  2013     (D. Storkey)  original code 
     7   !!            4.0  !  2014     (T. Lovato)  Generalize OBC structure 
    88   !!---------------------------------------------------------------------- 
     9    
    910   !!---------------------------------------------------------------------- 
    10    !!   bdy_orlanski_2d 
    11    !!   bdy_orlanski_3d 
     11   !!  bdy_frs        : Apply the Flow Relaxation Scheme (tracers) 
     12   !!  bdy_spe        : Apply a specified value (tracers) 
     13   !!  bdy_orl        : Apply Orlanski radiation (tracers) 
     14   !!  bdy_orlanski_2d:   2D      -        -         - 
     15   !!  bdy_orlanski_3d:   3D      -        -         - 
     16   !!  bdy_nmn        : Duplicate the value at open boundaries (zero gradient) 
    1217   !!---------------------------------------------------------------------- 
    1318   USE oce            ! ocean dynamics and tracers  
     
    2227   PRIVATE 
    2328 
    24    PUBLIC   bdy_frs, bdy_spe, bdy_nmn, bdy_orl 
    25    PUBLIC   bdy_orlanski_2d 
    26    PUBLIC   bdy_orlanski_3d 
     29   PUBLIC   bdy_frs, bdy_spe, bdy_nmn 
     30   PUBLIC   bdy_orl, bdy_orlanski_2d, bdy_orlanski_3d 
    2731 
    2832   !!---------------------------------------------------------------------- 
     
    230234         ! Note no rdt factor in expression for zdt because it cancels in the expressions for  
    231235         ! zrx and zry. 
    232          zdt = phia(iibm1,ijbm1) - phib(iibm1,ijbm1) 
     236         zdt =     phia(iibm1,ijbm1) - phib(iibm1,ijbm1) 
    233237         zdx = ( ( phia(iibm1,ijbm1) - phia(iibm2,ijbm2) ) / zex2 ) * zmask_x  
    234238         zdy_1 = ( ( phib(iibm1   ,ijbm1   ) - phib(iibm1jm1,ijbm1jm1) ) / zey1 ) * zmask_y1     
     
    247251         zout = sign( 1., zrx ) 
    248252         zout = 0.5*( zout + abs(zout) ) 
    249          zwgt = 2.*rdt*( (1.-zout) * idx%nbd(jb,igrd) + zout * idx%nbdout(jb,igrd) ) 
     253         zwgt = rDt*( (1.-zout) * idx%nbd(jb,igrd) + zout * idx%nbdout(jb,igrd) ) 
    250254         ! only apply radiation on outflow points  
    251255         if( ll_npo ) then     !! NPO version !! 
     
    385389            ! Centred derivative is calculated as average of "left" and "right" derivatives for  
    386390            ! this reason.  
    387             zdt = phia(iibm1,ijbm1,jk) - phib(iibm1,ijbm1,jk) 
     391            zdt =     phia(iibm1,ijbm1,jk) - phib(iibm1,ijbm1,jk) 
    388392            zdx = ( ( phia(iibm1,ijbm1,jk) - phia(iibm2,ijbm2,jk) ) / zex2 ) * zmask_x                   
    389393            zdy_1 = ( ( phib(iibm1   ,ijbm1   ,jk) - phib(iibm1jm1,ijbm1jm1,jk) ) / zey1 ) * zmask_y1   
     
    402406!!$            zrx = min(zrx,2.0_wp) 
    403407            zout = sign( 1., zrx ) 
    404             zout = 0.5*( zout + abs(zout) ) 
    405             zwgt = 2.*rdt*( (1.-zout) * idx%nbd(jb,igrd) + zout * idx%nbdout(jb,igrd) ) 
     408            zout = 0.5 * ( zout + abs(zout) ) 
     409            zwgt = rDt * ( (1.-zout) * idx%nbd(jb,igrd) + zout * idx%nbdout(jb,igrd) ) 
    406410            ! only apply radiation on outflow points  
    407411            if( ll_npo ) then     !! NPO version !! 
     
    426430      ! 
    427431   END SUBROUTINE bdy_orlanski_3d 
     432 
    428433 
    429434   SUBROUTINE bdy_nmn( idx, igrd, phia ) 
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/BDY/bdytides.F90

    r9598 r9939  
    295295      !!---------------------------------------------------------------------- 
    296296      ! 
    297       ilen0(1) =  SIZE(td%ssh(:,1,1)) 
    298       ilen0(2) =  SIZE(td%u(:,1,1)) 
    299       ilen0(3) =  SIZE(td%v(:,1,1)) 
     297      ilen0(1) =  SIZE( td%ssh(:,1,1) ) 
     298      ilen0(2) =  SIZE( td%u  (:,1,1) ) 
     299      ilen0(3) =  SIZE( td%v  (:,1,1) ) 
    300300 
    301301      zflag=1 
    302302      IF ( PRESENT(jit) ) THEN 
    303         IF ( jit /= 1 ) zflag=0 
     303        IF ( jit /= 1 )   zflag=0 
    304304      ENDIF 
    305305 
    306       IF ( (nsec_day == NINT(0.5_wp * rdt) .OR. kt==nit000) .AND. zflag==1 ) THEN 
     306      IF ( ( nsec_day == NINT( 0.5_wp * rn_Dt )  .OR.  kt == nit000 ) .AND. zflag==1 ) THEN 
    307307        ! 
    308         kt_tide = kt - (nsec_day - 0.5_wp * rdt)/rdt 
     308        kt_tide = kt - (nsec_day - 0.5_wp * rn_Dt) / rn_Dt 
    309309        ! 
    310310        IF(lwp) THEN 
    311311           WRITE(numout,*) 
    312            WRITE(numout,*) 'bdytide_update : (re)Initialization of the tidal bdy forcing at kt=',kt 
     312           WRITE(numout,*) 'bdytide_update : (re)Initialization of the tidal bdy forcing at kt=', kt 
    313313           WRITE(numout,*) '~~~~~~~~~~~~~~ ' 
    314314        ENDIF 
     
    325325          
    326326      IF( PRESENT(jit) ) THEN   
    327          z_arg = ((kt-kt_tide) * rdt + (jit+0.5_wp*(time_add-1)) * rdt / REAL(nn_baro,wp) ) 
     327         z_arg = ((kt-kt_tide) * rn_Dt + (jit+0.5_wp*(time_add-1)) * rn_Dt / REAL(nn_e,wp) ) 
    328328      ELSE                               
    329          z_arg = ((kt-kt_tide)+time_add) * rdt 
     329         z_arg = ((kt-kt_tide)+time_add) * rn_Dt 
    330330      ENDIF 
    331331 
    332332      ! Linear ramp on tidal component at open boundaries  
    333333      zramp = 1._wp 
    334       IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg + (kt_tide-nit000)*rdt)/(rdttideramp*rday),0._wp),1._wp) 
     334      IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg + (kt_tide-nit000)*rn_Dt)/(rn_ramp*rday),0._wp),1._wp) 
    335335 
    336336      DO itide = 1, nb_harmo 
     
    392392      ! Absolute time from model initialization:    
    393393      IF( PRESENT(kit) ) THEN   
    394          z_arg = ( kt + (kit+time_add-1) / REAL(nn_baro,wp) ) * rdt 
     394         z_arg = ( kt + (kit+time_add-1) / REAL(nn_e,wp) ) * rn_Dt 
    395395      ELSE                               
    396          z_arg = ( kt + time_add ) * rdt 
     396         z_arg = ( kt + time_add ) * rn_Dt 
    397397      ENDIF 
    398398 
    399399      ! Linear ramp on tidal component at open boundaries  
    400400      zramp = 1. 
    401       IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg - nit000*rdt)/(rdttideramp*rday),0.),1.) 
     401      IF ( ln_tide_ramp )   zramp = MIN(  MAX( 0. , (z_arg - nit000*rn_Dt)/(rn_ramp*rday) ) , 1.  ) 
    402402 
    403403      DO ib_bdy = 1,nb_bdy 
     
    414414            ! We refresh nodal factors every day below 
    415415            ! This should be done somewhere else 
    416             IF ( ( nsec_day == NINT(0.5_wp * rdt) .OR. kt==nit000 ) .AND. lk_first_btstp ) THEN 
    417                ! 
    418                kt_tide = kt - (nsec_day - 0.5_wp * rdt)/rdt 
     416            IF ( ( nsec_day == NINT(0.5_wp * rn_Dt) .OR. kt==nit000 ) .AND. lk_first_btstp ) THEN 
     417               ! 
     418               kt_tide = kt - (nsec_day - 0.5_wp * rn_Dt) / rn_Dt 
    419419               ! 
    420420               IF(lwp) THEN 
     
    428428               ! 
    429429            ENDIF 
    430             zoff = -kt_tide * rdt ! time offset relative to nodal factor computation time 
     430            zoff = -kt_tide * rn_Dt    ! time offset relative to nodal factor computation time 
    431431            ! 
    432432            ! If time splitting, initialize arrays from slow varying open boundary data: 
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/BDY/bdyvol.F90

    r9598 r9939  
    8484      ! ----------------------------------------------------------------------- 
    8585!!gm replace these lines : 
    86       z_cflxemp = SUM ( ( emp(:,:) - rnf(:,:) + fwfisf(:,:) ) * bdytmask(:,:) * e1e2t(:,:) ) / rau0 
     86      z_cflxemp = SUM ( ( emp(:,:) - rnf(:,:) + fwfisf(:,:) ) * bdytmask(:,:) * e1e2t(:,:) ) / rho0 
    8787      IF( lk_mpp )   CALL mpp_sum( z_cflxemp )     ! sum over the global domain 
    8888!!gm   by : 
    89 !!gm      z_cflxemp = glob_sum(  ( emp(:,:)-rnf(:,:)+fwfisf(:,:) ) * bdytmask(:,:) * e1e2t(:,:)  ) / rau0 
    90 !!gm 
     89!!gm      z_cflxemp = glob_sum(  ( emp(:,:)-rnf(:,:)+fwfisf(:,:) ) * bdytmask(:,:) * e1e2t(:,:)  ) / rho0 
     90!!gm   ??? 
    9191 
    9292      ! Transport through the unstructured open boundary 
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DIA/dia25h.F90

    r9598 r9939  
    139139      ! ----------------- 
    140140      ! Define frequency of summing to create 25 h mean 
    141       IF( MOD( 3600,INT(rdt) ) == 0 ) THEN 
    142          i_steps = 3600/INT(rdt) 
     141      IF( MOD( 3600 , INT(rn_Dt) ) == 0 ) THEN 
     142         i_steps = 3600 / INT( rn_Dt ) 
    143143      ELSE 
    144          CALL ctl_stop('STOP', 'dia_wri_tide: timestep must give MOD(3600,rdt) = 0 otherwise no hourly values are possible') 
     144         CALL ctl_stop('STOP', 'dia_wri_tide: timestep must give MOD(3600,rn_Dt) = 0 otherwise no hourly values are possible') 
    145145      ENDIF 
    146146 
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DIA/diaar5.F90

    r9598 r9939  
    161161       
    162162         !                                         ! ocean bottom pressure 
    163          zztmp = rau0 * grav * 1.e-4_wp               ! recover pressure from pressure anomaly and cover to dbar = 1.e4 Pa 
     163         zztmp = rho0 * grav * 1.e-4_wp               ! recover pressure from pressure anomaly and cover to dbar = 1.e4 Pa 
    164164         zbotpres(:,:) = zztmp * ( zbotpres(:,:) + sshn(:,:) + thick0(:,:) ) 
    165165         CALL iom_put( 'botpres', zbotpres ) 
     
    198198         END IF 
    199199         ! 
    200          zmass = rau0 * ( zarho + zvol )                 ! total mass of liquid seawater 
     200         zmass = rho0 * ( zarho + zvol )                 ! total mass of liquid seawater 
    201201         ztemp = ztemp / zvol                            ! potential temperature in liquid seawater 
    202202         zsal  = zsal  / zvol                            ! Salinity of liquid seawater 
     
    239239               DO ji = 1, jpi 
    240240                  DO jj = 1, jpj 
    241                      zpe(ji,jj) = zpe(ji,jj) + avt(ji, jj, jk) * MIN(0._wp,rn2(ji, jj, jk)) * rau0 * e3w_n(ji, jj, jk) 
     241                     zpe(ji,jj) = zpe(ji,jj) + avt(ji, jj, jk) * MIN(0._wp,rn2(ji, jj, jk)) * rho0 * e3w_n(ji, jj, jk) 
    242242                  END DO 
    243243               END DO 
     
    287287       CALL lbc_lnk( z2d, 'U', -1. ) 
    288288       IF( cptr == 'adv' ) THEN 
    289           IF( ktra == jp_tem ) CALL iom_put( "uadv_heattr" , rau0_rcp * z2d )  ! advective heat transport in i-direction 
    290           IF( ktra == jp_sal ) CALL iom_put( "uadv_salttr" , rau0     * z2d )  ! advective salt transport in i-direction 
     289          IF( ktra == jp_tem ) CALL iom_put( "uadv_heattr" , rho0_rcp * z2d )  ! advective heat transport in i-direction 
     290          IF( ktra == jp_sal ) CALL iom_put( "uadv_salttr" , rho0     * z2d )  ! advective salt transport in i-direction 
    291291       ENDIF 
    292292       IF( cptr == 'ldf' ) THEN 
    293           IF( ktra == jp_tem ) CALL iom_put( "udiff_heattr" , rau0_rcp * z2d ) ! diffusive heat transport in i-direction 
    294           IF( ktra == jp_sal ) CALL iom_put( "udiff_salttr" , rau0     * z2d ) ! diffusive salt transport in i-direction 
     293          IF( ktra == jp_tem ) CALL iom_put( "udiff_heattr" , rho0_rcp * z2d ) ! diffusive heat transport in i-direction 
     294          IF( ktra == jp_sal ) CALL iom_put( "udiff_salttr" , rho0     * z2d ) ! diffusive salt transport in i-direction 
    295295       ENDIF 
    296296       ! 
     
    305305       CALL lbc_lnk( z2d, 'V', -1. ) 
    306306       IF( cptr == 'adv' ) THEN 
    307           IF( ktra == jp_tem ) CALL iom_put( "vadv_heattr" , rau0_rcp * z2d )  ! advective heat transport in j-direction 
    308           IF( ktra == jp_sal ) CALL iom_put( "vadv_salttr" , rau0     * z2d )  ! advective salt transport in j-direction 
     307          IF( ktra == jp_tem ) CALL iom_put( "vadv_heattr" , rho0_rcp * z2d )  ! advective heat transport in j-direction 
     308          IF( ktra == jp_sal ) CALL iom_put( "vadv_salttr" , rho0     * z2d )  ! advective salt transport in j-direction 
    309309       ENDIF 
    310310       IF( cptr == 'ldf' ) THEN 
    311           IF( ktra == jp_tem ) CALL iom_put( "vdiff_heattr" , rau0_rcp * z2d ) ! diffusive heat transport in j-direction 
    312           IF( ktra == jp_sal ) CALL iom_put( "vdiff_salttr" , rau0     * z2d ) ! diffusive salt transport in j-direction 
     311          IF( ktra == jp_tem ) CALL iom_put( "vdiff_heattr" , rho0_rcp * z2d ) ! diffusive heat transport in j-direction 
     312          IF( ktra == jp_sal ) CALL iom_put( "vdiff_salttr" , rho0     * z2d ) ! diffusive salt transport in j-direction 
    313313       ENDIF 
    314314           
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DIA/diacfl.F90

    r9598 r9939  
    5555      ! 
    5656      INTEGER :: ji, jj, jk   ! dummy loop indices 
    57       REAL(wp)::   z2dt, zCu_max, zCv_max, zCw_max       ! local scalars 
     57      REAL(wp)::   zCu_max, zCv_max, zCw_max       ! local scalars 
    5858      INTEGER , DIMENSION(3)           ::   iloc_u , iloc_v , iloc_w , iloc   ! workspace 
    5959!!gm this does not work      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zCu_cfl, zCv_cfl, zCw_cfl         ! workspace 
     
    6262      IF( ln_timing )   CALL timing_start('dia_cfl') 
    6363      ! 
    64       !                       ! setup timestep multiplier to account for initial Eulerian timestep 
    65       IF( neuler == 0 .AND. kt == nit000 ) THEN   ;    z2dt = rdt 
    66       ELSE                                        ;    z2dt = rdt * 2._wp 
    67       ENDIF 
    68       ! 
    6964      !                 
    7065      DO jk = 1, jpk       ! calculate Courant numbers 
    7166         DO jj = 1, jpj 
    7267            DO ji = 1, fs_jpim1   ! vector opt. 
    73                zCu_cfl(ji,jj,jk) = ABS( un(ji,jj,jk) ) * z2dt / e1u  (ji,jj)      ! for i-direction 
    74                zCv_cfl(ji,jj,jk) = ABS( vn(ji,jj,jk) ) * z2dt / e2v  (ji,jj)      ! for j-direction 
    75                zCw_cfl(ji,jj,jk) = ABS( wn(ji,jj,jk) ) * z2dt / e3w_n(ji,jj,jk)   ! for k-direction 
     68               zCu_cfl(ji,jj,jk) = ABS( un(ji,jj,jk) ) * rDt / e1u  (ji,jj)      ! for i-direction 
     69               zCv_cfl(ji,jj,jk) = ABS( vn(ji,jj,jk) ) * rDt / e2v  (ji,jj)      ! for j-direction 
     70               zCw_cfl(ji,jj,jk) = ABS( wn(ji,jj,jk) ) * rDt / e3w_n(ji,jj,jk)   ! for k-direction 
    7671            END DO 
    7772         END DO          
     
    120115         WRITE(numcfl,*) '******************************************' 
    121116         WRITE(numcfl,FMT='(3x,a12,6x,f7.4,1x,i4,1x,i4,1x,i4)') 'Run Max Cu', rCu_max, nCu_loc(1), nCu_loc(2), nCu_loc(3) 
    122          WRITE(numcfl,FMT='(3x,a8,11x,f15.1)') ' => dt/C', z2dt/rCu_max 
     117         WRITE(numcfl,FMT='(3x,a8,11x,f15.1)') ' => dt/C', rDt/rCu_max 
    123118         WRITE(numcfl,*) '******************************************' 
    124119         WRITE(numcfl,FMT='(3x,a12,6x,f7.4,1x,i4,1x,i4,1x,i4)') 'Run Max Cv', rCv_max, nCv_loc(1), nCv_loc(2), nCv_loc(3) 
    125          WRITE(numcfl,FMT='(3x,a8,11x,f15.1)') ' => dt/C', z2dt/rCv_max 
     120         WRITE(numcfl,FMT='(3x,a8,11x,f15.1)') ' => dt/C', rDt/rCv_max 
    126121         WRITE(numcfl,*) '******************************************' 
    127122         WRITE(numcfl,FMT='(3x,a12,6x,f7.4,1x,i4,1x,i4,1x,i4)') 'Run Max Cw', rCw_max, nCw_loc(1), nCw_loc(2), nCw_loc(3) 
    128          WRITE(numcfl,FMT='(3x,a8,11x,f15.1)') ' => dt/C', z2dt/rCw_max 
     123         WRITE(numcfl,FMT='(3x,a8,11x,f15.1)') ' => dt/C', rDt/rCw_max 
    129124         CLOSE( numcfl )  
    130125         ! 
     
    133128         WRITE(numout,*) 'dia_cfl : Maximum Courant number information for the run ' 
    134129         WRITE(numout,*) '~~~~~~~' 
    135          WRITE(numout,*) '   Max Cu = ', rCu_max, ' at (i,j,k) = (',nCu_loc(1),nCu_loc(2),nCu_loc(3),') => dt/C = ', z2dt/rCu_max 
    136          WRITE(numout,*) '   Max Cv = ', rCv_max, ' at (i,j,k) = (',nCv_loc(1),nCv_loc(2),nCv_loc(3),') => dt/C = ', z2dt/rCv_max 
    137          WRITE(numout,*) '   Max Cw = ', rCw_max, ' at (i,j,k) = (',nCw_loc(1),nCw_loc(2),nCw_loc(3),') => dt/C = ', z2dt/rCw_max 
     130         WRITE(numout,*) '   Max Cu = ', rCu_max, ' at (i,j,k) = (',nCu_loc(1),nCu_loc(2),nCu_loc(3),') => dt/C = ', rDt/rCu_max 
     131         WRITE(numout,*) '   Max Cv = ', rCv_max, ' at (i,j,k) = (',nCv_loc(1),nCv_loc(2),nCv_loc(3),') => dt/C = ', rDt/rCv_max 
     132         WRITE(numout,*) '   Max Cw = ', rCw_max, ' at (i,j,k) = (',nCw_loc(1),nCw_loc(2),nCw_loc(3),') => dt/C = ', rDt/rCw_max 
    138133      ENDIF 
    139134      ! 
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DIA/diadct.F90

    r9598 r9939  
    679679                  zsn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) )  
    680680                  zrhop = interp(k%I,k%J,jk,'V',rhop)  
    681                   zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0)  
     681                  zrhoi = interp(k%I,k%J,jk,'V',rhd*rho0+rho0)  
    682682                  zsshn =  0.5*( sshn(k%I,k%J) + sshn(k%I,k%J+1)    ) * vmask(k%I,k%J,1)  
    683683               CASE(2,3)  
     
    685685                  zsn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) )  
    686686                  zrhop = interp(k%I,k%J,jk,'U',rhop)  
    687                   zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0)  
     687                  zrhoi = interp(k%I,k%J,jk,'U',rhd*rho0+rho0)  
    688688                  zsshn =  0.5*( sshn(k%I,k%J) + sshn(k%I+1,k%J)    ) * umask(k%I,k%J,1)   
    689689               END SELECT  
     
    851851                 zsn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) )  
    852852                 zrhop = interp(k%I,k%J,jk,'V',rhop)  
    853                  zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0)  
     853                 zrhoi = interp(k%I,k%J,jk,'V',rhd*rho0+rho0)  
    854854 
    855855              CASE(2,3)  
     
    857857                 zsn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) )  
    858858                 zrhop = interp(k%I,k%J,jk,'U',rhop)  
    859                  zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0)  
     859                 zrhoi = interp(k%I,k%J,jk,'U',rhd*rho0+rho0)  
    860860                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I+1,k%J)    ) * umask(k%I,k%J,1)   
    861861              END SELECT  
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DIA/diaharm.F90

    r9598 r9939  
    181181      IF( kt >= nit000_han .AND. kt <= nitend_han .AND. MOD(kt,nstep_han) == 0 ) THEN 
    182182         ! 
    183          ztime = (kt-nit000+1) * rdt  
     183         ztime = ( kt - nit000+1 ) * rn_Dt  
    184184         ! 
    185185         nhc = 0 
     
    231231      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 
    232232 
    233       ztime_ini = nit000_han*rdt                 ! Initial time in seconds at the beginning of analysis 
    234       ztime_end = nitend_han*rdt                 ! Final time in seconds at the end of analysis 
     233      ztime_ini = nit000_han*rn_Dt                 ! Initial time in seconds at the beginning of analysis 
     234      ztime_end = nitend_han*rn_Dt                 ! Final time in seconds at the end of analysis 
    235235      nhan = (nitend_han-nit000_han+1)/nstep_han ! Number of dumps used for analysis 
    236236 
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DIA/diahsb.F90

    r9598 r9939  
    9191      ! 1 - Trends due to forcing ! 
    9292      ! ------------------------- ! 
    93       z_frc_trd_v = r1_rau0 * glob_sum( - ( emp(:,:) - rnf(:,:) + fwfisf(:,:) ) * surf(:,:) )   ! volume fluxes 
     93      z_frc_trd_v = r1_rho0 * glob_sum( - ( emp(:,:) - rnf(:,:) + fwfisf(:,:) ) * surf(:,:) )   ! volume fluxes 
    9494      z_frc_trd_t =           glob_sum( sbc_tsc(:,:,jp_tem) * surf(:,:) )                       ! heat fluxes 
    9595      z_frc_trd_s =           glob_sum( sbc_tsc(:,:,jp_sal) * surf(:,:) )                       ! salt fluxes 
     
    100100      IF( ln_isf    )   z_frc_trd_t = z_frc_trd_t + glob_sum( risf_tsc(:,:,jp_tem) * surf(:,:) ) 
    101101      !                    ! Add penetrative solar radiation 
    102       IF( ln_traqsr )   z_frc_trd_t = z_frc_trd_t + r1_rau0_rcp * glob_sum( qsr     (:,:) * surf(:,:) ) 
     102      IF( ln_traqsr )   z_frc_trd_t = z_frc_trd_t + r1_rho0_rcp * glob_sum( qsr     (:,:) * surf(:,:) ) 
    103103      !                    ! Add geothermal heat flux 
    104104      IF( ln_trabbc )   z_frc_trd_t = z_frc_trd_t +               glob_sum( qgh_trd0(:,:) * surf(:,:) ) 
     
    120120      ENDIF 
    121121 
    122       frc_v = frc_v + z_frc_trd_v * rdt 
    123       frc_t = frc_t + z_frc_trd_t * rdt 
    124       frc_s = frc_s + z_frc_trd_s * rdt 
     122      frc_v = frc_v + z_frc_trd_v * rn_Dt 
     123      frc_t = frc_t + z_frc_trd_t * rn_Dt 
     124      frc_s = frc_s + z_frc_trd_s * rn_Dt 
    125125      !                                          ! Advection flux through fixed surface (z=0) 
    126126      IF( ln_linssh ) THEN 
    127          frc_wn_t = frc_wn_t + z_wn_trd_t * rdt 
    128          frc_wn_s = frc_wn_s + z_wn_trd_s * rdt 
     127         frc_wn_t = frc_wn_t + z_wn_trd_t * rn_Dt 
     128         frc_wn_s = frc_wn_s + z_wn_trd_s * rn_Dt 
    129129      ENDIF 
    130130 
     
    196196 
    197197      CALL iom_put(   'bgfrcvol' , frc_v    * 1.e-9    )              ! vol - surface forcing (km3)  
    198       CALL iom_put(   'bgfrctem' , frc_t    * rau0 * rcp * 1.e-20 )   ! hc  - surface forcing (1.e20 J)  
    199       CALL iom_put(   'bgfrchfx' , frc_t    * rau0 * rcp /  &         ! hc  - surface forcing (W/m2)  
    200          &                       ( surf_tot * kt * rdt )        ) 
     198      CALL iom_put(   'bgfrctem' , frc_t    * rho0_rcp * 1.e-20 )     ! hc  - surface forcing (1.e20 J)  
     199      CALL iom_put(   'bgfrchfx' , frc_t    * rho0_rcp /  &           ! hc  - surface forcing (W/m2)  
     200         &                         ( surf_tot * kt * rn_Dt )    ) 
    201201      CALL iom_put(   'bgfrcsal' , frc_s    * 1.e-9    )              ! sc  - surface forcing (psu*km3)  
    202202 
     
    204204         CALL iom_put( 'bgtemper' , zdiff_hc / zvol_tot )              ! Temperature drift     (C)  
    205205         CALL iom_put( 'bgsaline' , zdiff_sc / zvol_tot )              ! Salinity    drift     (PSU) 
    206          CALL iom_put( 'bgheatco' , zdiff_hc * 1.e-20 * rau0 * rcp )   ! Heat content drift    (1.e20 J)  
    207          CALL iom_put( 'bgheatfx' , zdiff_hc * rau0 * rcp /  &         ! Heat flux drift       (W/m2)  
    208             &                       ( surf_tot * kt * rdt )        ) 
     206         CALL iom_put( 'bgheatco' , zdiff_hc * 1.e-20 * rho0_rcp )     ! Heat content drift    (1.e20 J)  
     207         CALL iom_put( 'bgheatfx' , zdiff_hc * rho0_rcp /  &           ! Heat flux drift       (W/m2)  
     208            &                       ( surf_tot * kt * rn_Dt )    ) 
    209209         CALL iom_put( 'bgsaltco' , zdiff_sc * 1.e-9    )              ! Salt content drift    (psu*km3) 
    210210         CALL iom_put( 'bgvolssh' , zdiff_v1 * 1.e-9    )              ! volume ssh drift      (km3)   
     
    224224         CALL iom_put( 'bgtemper' , zdiff_hc1 / zvol_tot)              ! Heat content drift    (C)  
    225225         CALL iom_put( 'bgsaline' , zdiff_sc1 / zvol_tot)              ! Salt content drift    (PSU) 
    226          CALL iom_put( 'bgheatco' , zdiff_hc1 * 1.e-20 * rau0 * rcp )  ! Heat content drift    (1.e20 J)  
    227          CALL iom_put( 'bgheatfx' , zdiff_hc1 * rau0 * rcp /  &        ! Heat flux drift       (W/m2)  
    228             &                       ( surf_tot * kt * rdt )         ) 
     226         CALL iom_put( 'bgheatco' , zdiff_hc1 * 1.e-20 * rho0_rcp )    ! Heat content drift    (1.e20 J)  
     227         CALL iom_put( 'bgheatfx' , zdiff_hc1 * rho0_rcp /  &          ! Heat flux drift       (W/m2)  
     228            &                       ( surf_tot * kt * rn_Dt )     ) 
    229229         CALL iom_put( 'bgsaltco' , zdiff_sc1 * 1.e-9    )             ! Salt content drift    (psu*km3) 
    230230         CALL iom_put( 'bgvolssh' , zdiff_v1 * 1.e-9    )              ! volume ssh drift      (km3)   
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DIA/diahth.F90

    r9598 r9939  
    8989      REAL(wp)                         ::   zrho1 = 0.01_wp       ! density     criterion for mixed layer depth 
    9090      REAL(wp)                         ::   ztem2 = 0.2_wp        ! temperature criterion for mixed layer depth 
    91       REAL(wp)                         ::   zthick_0, zcoef       ! temporary scalars 
    92       REAL(wp)                         ::   zztmp, zzdep          ! temporary scalars inside do loop 
    93       REAL(wp)                         ::   zu, zv, zw, zut, zvt  ! temporary workspace 
     91      REAL(wp)                         ::   zthick_0              ! local scalars 
     92      REAL(wp)                         ::   zztmp, zzdep          ! local scalars inside do loop 
     93      REAL(wp)                         ::   zu, zv, zw, zut, zvt  ! local workspace 
    9494      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zabs2      ! MLD: abs( tn - tn(10m) ) = ztem2  
    9595      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ztm2       ! Top of thermocline: tn = tn(10m) - ztem2      
     
    328328      END DO 
    329329      ! from temperature to heat contain 
    330       zcoef = rau0 * rcp 
    331       htc3(:,:) = zcoef * htc3(:,:) 
     330      htc3(:,:) = rho0_rcp * htc3(:,:) 
    332331      CALL iom_put( "hc300", htc3 )      ! first 300m heat content 
    333332      ! 
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DIA/dianam.F90

    r9598 r9939  
    7171      ENDIF 
    7272 
    73       IF( llfsec .OR. kfreq < 0 ) THEN   ;   inbsec = kfreq                       ! output frequency already in seconds 
    74       ELSE                               ;   inbsec = kfreq * NINT( rdt )   ! from time-step to seconds 
     73      IF( llfsec .OR. kfreq < 0 ) THEN   ;   inbsec = kfreq                    ! output frequency already in seconds 
     74      ELSE                               ;   inbsec = kfreq * NINT( rn_Dt )    ! from time-step to seconds 
    7575      ENDIF 
    7676      iddss = NINT( rday          )                                         ! number of seconds in 1 day 
     
    116116      ! date of the beginning and the end of the run 
    117117 
    118       zdrun = rdt / rday * REAL( nitend - nit000, wp )                ! length of the run in days 
    119       zjul  = fjulday - rdt / rday 
     118      zdrun = rn_Dt / rday * REAL( nitend - nit000, wp )              ! length of the run in days 
     119      zjul  = fjulday - rn_Dt / rday 
    120120      CALL ju2ymds( zjul        , iyear1, imonth1, iday1, zsec1 )           ! year/month/day of the beginning of run 
    121121      CALL ju2ymds( zjul + zdrun, iyear2, imonth2, iday2, zsec2 )           ! year/month/day of the end       of run 
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DIA/diaptr.F90

    r9598 r9939  
    5252 
    5353   REAL(wp) ::   rc_sv    = 1.e-6_wp   ! conversion from m3/s to Sverdrup 
    54    REAL(wp) ::   rc_pwatt = 1.e-15_wp  ! conversion from W    to PW (further x rau0 x Cp) 
     54   REAL(wp) ::   rc_pwatt = 1.e-15_wp  ! conversion from W    to PW (further x rho0 x Cp) 
    5555   REAL(wp) ::   rc_ggram = 1.e-6_wp   ! conversion from g    to Pg 
    5656 
     
    424424         IF( dia_ptr_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dia_ptr_init : unable to allocate arrays' ) 
    425425 
    426          rc_pwatt = rc_pwatt * rau0_rcp          ! conversion from K.s-1 to PetaWatt 
     426         rc_pwatt = rc_pwatt * rho0_rcp          ! conversion from K.s-1 to PetaWatt 
    427427 
    428428         IF( lk_mpp )   CALL mpp_ini_znl( numout )     ! Define MPI communicator for zonal sum 
     
    448448         ! Initialise arrays to zero because diatpr is called before they are first calculated 
    449449         ! Note that this means diagnostics will not be exactly correct when model run is restarted. 
    450          htr_adv(:,:) = 0._wp  ;  str_adv(:,:) =  0._wp  
    451          htr_ldf(:,:) = 0._wp  ;  str_ldf(:,:) =  0._wp  
    452          htr_eiv(:,:) = 0._wp  ;  str_eiv(:,:) =  0._wp  
     450         htr_adv(:,:) = 0._wp  ;   str_adv(:,:) =  0._wp  
     451         htr_ldf(:,:) = 0._wp  ;   str_ldf(:,:) =  0._wp  
     452         htr_eiv(:,:) = 0._wp  ;   str_eiv(:,:) =  0._wp  
    453453         htr_ove(:,:) = 0._wp  ;   str_ove(:,:) =  0._wp 
    454454         htr_btr(:,:) = 0._wp  ;   str_btr(:,:) =  0._wp 
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DIA/diawri.F90

    r9652 r9939  
    169169 
    170170      IF ( iom_use("taubot") ) THEN                ! bottom stress 
    171          zztmp = rau0 * 0.25 
     171         zztmp = rho0 * 0.25 
    172172         z2d(:,:) = 0._wp 
    173173         DO jj = 2, jpjm1 
     
    212212      IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN   ! vertical mass transport & its square value 
    213213         ! Caution: in the VVL case, it only correponds to the baroclinic mass transport. 
    214          z2d(:,:) = rau0 * e1e2t(:,:) 
     214         z2d(:,:) = rho0 * e1e2t(:,:) 
    215215         DO jk = 1, jpk 
    216216            z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:) 
     
    253253            END DO 
    254254         END DO 
    255          CALL iom_put( "heatc", rau0_rcp * z2d )   ! vertically integrated heat content (J/m2) 
     255         CALL iom_put( "heatc", rho0_rcp * z2d )   ! vertically integrated heat content (J/m2) 
    256256      ENDIF 
    257257 
     
    265265            END DO 
    266266         END DO 
    267          CALL iom_put( "saltc", rau0 * z2d )          ! vertically integrated salt content (PSU*kg/m2) 
     267         CALL iom_put( "saltc", rho0 * z2d )          ! vertically integrated salt content (PSU*kg/m2) 
    268268      ENDIF 
    269269      ! 
     
    291291         z2d(:,:) = 0.e0 
    292292         DO jk = 1, jpkm1 
    293             z3d(:,:,jk) = rau0 * un(:,:,jk) * e2u(:,:) * e3u_n(:,:,jk) * umask(:,:,jk) 
     293            z3d(:,:,jk) = rho0 * un(:,:,jk) * e2u(:,:) * e3u_n(:,:,jk) * umask(:,:,jk) 
    294294            z2d(:,:) = z2d(:,:) + z3d(:,:,jk) 
    295295         END DO 
     
    328328         z3d(:,:,jpk) = 0.e0 
    329329         DO jk = 1, jpkm1 
    330             z3d(:,:,jk) = rau0 * vn(:,:,jk) * e1v(:,:) * e3v_n(:,:,jk) * vmask(:,:,jk) 
     330            z3d(:,:,jk) = rho0 * vn(:,:,jk) * e1v(:,:) * e3v_n(:,:,jk) * vmask(:,:,jk) 
    331331         END DO 
    332332         CALL iom_put( "v_masstr", z3d )              ! mass transport in j-direction 
     
    369369         END DO 
    370370         CALL lbc_lnk( z2d, 'T', -1. ) 
    371          CALL iom_put( "tosmint", rau0 * z2d )        ! Vertical integral of temperature 
     371         CALL iom_put( "tosmint", rho0 * z2d )        ! Vertical integral of temperature 
    372372      ENDIF 
    373373      IF( iom_use("somint") ) THEN 
     
    381381         END DO 
    382382         CALL lbc_lnk( z2d, 'T', -1. ) 
    383          CALL iom_put( "somint", rau0 * z2d )         ! Vertical integral of salinity 
     383         CALL iom_put( "somint", rho0 * z2d )         ! Vertical integral of salinity 
    384384      ENDIF 
    385385 
     
    458458      clop = "x"         ! no use of the mask value (require less cpu time and otherwise the model crashes) 
    459459#if defined key_diainstant 
    460       zsto = nwrite * rdt 
     460      zsto = nwrite * rn_Dt 
    461461      clop = "inst("//TRIM(clop)//")" 
    462462#else 
    463       zsto=rdt 
     463      zsto = rn_Dt 
    464464      clop = "ave("//TRIM(clop)//")" 
    465465#endif 
    466       zout = nwrite * rdt 
    467       zmax = ( nitend - nit000 + 1 ) * rdt 
     466      zout = nwrite * rn_Dt 
     467      zmax = ( nitend - nit000 + 1 ) * rn_Dt 
    468468 
    469469      ! Define indices of the horizontal output zoom and vertical limit storage 
     
    485485 
    486486         ! Compute julian date from starting date of the run 
    487          CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian ) 
     487         CALL ymds2ju( nyear, nmonth, nday, rn_Dt, zjulian ) 
    488488         zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment 
    489489         IF(lwp)WRITE(numout,*) 
     
    507507         CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit 
    508508            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       & 
    509             &          nit000-1, zjulian, rdt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set ) 
     509            &          nit000-1, zjulian, rn_Dt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set ) 
    510510         CALL histvert( nid_T, "deptht", "Vertical T levels",      &  ! Vertical grid: gdept 
    511511            &           "m", ipk, gdept_1d, nz_T, "down" ) 
     
    543543         CALL histbeg( clhstnam, jpi, glamu, jpj, gphiu,           &  ! Horizontal grid: glamu and gphiu 
    544544            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       & 
    545             &          nit000-1, zjulian, rdt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set ) 
     545            &          nit000-1, zjulian, rn_Dt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set ) 
    546546         CALL histvert( nid_U, "depthu", "Vertical U levels",      &  ! Vertical grid: gdept 
    547547            &           "m", ipk, gdept_1d, nz_U, "down" ) 
     
    556556         CALL histbeg( clhstnam, jpi, glamv, jpj, gphiv,           &  ! Horizontal grid: glamv and gphiv 
    557557            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       & 
    558             &          nit000-1, zjulian, rdt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set ) 
     558            &          nit000-1, zjulian, rn_Dt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set ) 
    559559         CALL histvert( nid_V, "depthv", "Vertical V levels",      &  ! Vertical grid : gdept 
    560560            &          "m", ipk, gdept_1d, nz_V, "down" ) 
     
    569569         CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit 
    570570            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       & 
    571             &          nit000-1, zjulian, rdt, nh_W, nid_W, domain_id=nidom, snc4chunks=snc4set ) 
     571            &          nit000-1, zjulian, rn_Dt, nh_W, nid_W, domain_id=nidom, snc4chunks=snc4set ) 
    572572         CALL histvert( nid_W, "depthw", "Vertical W levels",      &  ! Vertical grid: gdepw 
    573573            &          "m", ipk, gdepw_1d, nz_W, "down" ) 
     
    897897      clname = cdfile_name 
    898898      IF( .NOT. Agrif_Root() ) clname = TRIM(Agrif_CFixed())//'_'//TRIM(clname) 
    899       zsto = rdt 
     899      zsto = rn_Dt 
    900900      clop = "inst(x)"           ! no use of the mask value (require less cpu time) 
    901       zout = rdt 
    902       zmax = ( nitend - nit000 + 1 ) * rdt 
     901      zout = rn_Dt 
     902      zmax = ( nitend - nit000 + 1 ) * rn_Dt 
    903903 
    904904      IF(lwp) WRITE(numout,*) 
     
    912912 
    913913      ! Compute julian date from starting date of the run 
    914       CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )         ! time axis  
     914      CALL ymds2ju( nyear, nmonth, nday, rn_Dt, zjulian )         ! time axis  
    915915      zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment 
    916916      CALL histbeg( clname, jpi, glamt, jpj, gphit,   & 
    917           1, jpi, 1, jpj, nit000-1, zjulian, rdt, nh_i, id_i, domain_id=nidom, snc4chunks=snc4set ) ! Horizontal grid : glamt and gphit 
     917          1, jpi, 1, jpj, nit000-1, zjulian, rn_Dt, nh_i, id_i, domain_id=nidom, snc4chunks=snc4set ) ! Horizontal grid : glamt and gphit 
    918918      CALL histvert( id_i, "deptht", "Vertical T levels",   &    ! Vertical grid : gdept 
    919919          "m", jpk, gdept_1d, nz_i, "down") 
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DIU/cool_skin.F90

    r9598 r9939  
    6868 
    6969 
    70    SUBROUTINE diurnal_sst_coolskin_step(psqflux, pstauflux, psrho, rdt) 
     70   SUBROUTINE diurnal_sst_coolskin_step(psqflux, pstauflux, psrho, pdt) 
    7171      !!---------------------------------------------------------------------- 
    7272      !! *** ROUTINE diurnal_sst_takaya_step *** 
     
    8282      REAL(wp), INTENT(IN), DIMENSION(jpi,jpj) :: pstauflux   ! Wind stress (kg/ m s^2) 
    8383      REAL(wp), INTENT(IN), DIMENSION(jpi,jpj) :: psrho       ! Water density (kg/m^3) 
    84       REAL(wp), INTENT(IN) :: rdt                             ! Time-step 
     84      REAL(wp), INTENT(IN)                     :: pdt         ! Time-step (s) 
    8585      
    8686      ! Local variables 
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DIU/diurnal_bulk.F90

    r9168 r9939  
    7878 
    7979 
    80    SUBROUTINE diurnal_sst_takaya_step(kt, psolflux, pqflux, ptauflux, prho, p_rdt,   & 
     80   SUBROUTINE diurnal_sst_takaya_step(kt, psolflux, pqflux, ptauflux, prho, pdt,   & 
    8181            &                  pla, pthick, pcoolthick, pmu, & 
    8282            &                  p_fvel_bkginc, p_hflux_bkginc) 
     
    9898      REAL(wp), DIMENSION(jpi,jpj)          , INTENT(in) ::   ptauflux       ! wind stress  (kg/ m s^2) 
    9999      REAL(wp), DIMENSION(jpi,jpj)          , INTENT(in) ::   prho           ! water density  (kg/m^3) 
    100       REAL(wp)                              , INTENT(in) ::   p_rdt          ! time-step 
     100      REAL(wp)                              , INTENT(in) ::   pdt            ! time-step (s) 
    101101      REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in) ::   pLa            ! Langmuir number 
    102102      REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in) ::   pthick         ! warm layer thickness (m) 
     
    167167       
    168168      ! Increment the temperature using the implicit solution 
    169       x_dsst(:,:) = t_imp( x_dsst(:,:), p_rdt, z_abflux(:,:), z_fvel(:,:),   & 
     169      x_dsst(:,:) = t_imp( x_dsst(:,:), pdt, z_abflux(:,:), z_fvel(:,:),   & 
    170170         &                       z_fla(:,:), zmu(:,:), zthick(:,:), prho(:,:) ) 
    171171      ! 
     
    173173 
    174174    
    175    FUNCTION t_imp(p_dsst, p_rdt, p_abflux, p_fvel, & 
     175   FUNCTION t_imp(p_dsst, pdt, p_abflux, p_fvel, & 
    176176                          p_fla, pmu, pthick, prho ) 
    177177                           
     
    182182      ! Dummy variables 
    183183      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: p_dsst     ! Delta SST 
    184       REAL(wp), INTENT(IN)                     :: p_rdt      ! Time-step 
     184      REAL(wp), INTENT(IN)                     :: pdt        ! Time-step (s) 
    185185      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: p_abflux   ! Heat forcing 
    186186      REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: p_fvel     ! Friction velocity 
     
    257257            &      ( pthick(ji,jj) * z_stabfunc ) )      
    258258           
    259             t_imp(ji,jj) = ( p_dsst(ji,jj) + p_rdt * z_term1 ) / & 
    260                            ( 1._wp - p_rdt * z_term2 ) 
     259            t_imp(ji,jj) = ( p_dsst(ji,jj) + pdt * z_term1 ) / & 
     260                           ( 1._wp - pdt * z_term2 ) 
    261261 
    262262         END DO 
    263263      END DO 
    264264       
    265       END FUNCTION t_imp 
    266  
     265   END FUNCTION t_imp 
     266 
     267   !!====================================================================== 
    267268END MODULE diurnal_bulk 
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DIU/step_diu.F90

    r9598 r9939  
    55   !!====================================================================== 
    66   !! History :  3.7  ! 2015-11  (J. While)  Original code 
     7   !!---------------------------------------------------------------------- 
    78 
    89   USE diurnal_bulk    ! diurnal SST bulk routines  (diurnal_sst_takaya routine)  
     
    2728   !! Software governed by the CeCILL licence     (./LICENSE) 
    2829   !!---------------------------------------------------------------------- 
    29  
    3030   CONTAINS 
    3131 
    3232   SUBROUTINE stp_diurnal( kstp )  
    33       INTEGER, INTENT(in) ::   kstp   ! ocean time-step index  
    3433      !!----------------------------------------------------------------------  
    3534      !!                     ***  ROUTINE stp_diurnal  ***  
     
    4645      !!              -8- Outputs and diagnostics  
    4746      !!----------------------------------------------------------------------  
     47      INTEGER, INTENT(in) ::   kstp   ! ocean time-step index  
     48      ! 
    4849      INTEGER ::   jk       ! dummy loop indices 
    4950      INTEGER ::   indic    ! error indicator if < 0  
     
    5152      !! ---------------------------------------------------------------------  
    5253       
    53       IF(ln_diurnal_only) THEN 
     54      IF( ln_diurnal_only ) THEN 
    5455         indic = 0                                 ! reset to no error condition  
    5556         IF( kstp /= nit000 )   CALL day( kstp )   ! Calendar (day was already called at nit000 in day_init)  
     
    6061         ENDIF 
    6162        
    62             CALL sbc    ( kstp )                      ! Sea Boundary Conditions  
     63         CALL sbc( kstp )                          ! Sea Surface Boundary Conditions  
    6364      ENDIF 
    6465      
    65       ! Cool skin 
    6666      IF( .NOT.ln_diurnal )   CALL ctl_stop( "stp_diurnal: ln_diurnal not set" ) 
    6767          
    6868      IF( .NOT. ln_blk    )   CALL ctl_stop( "stp_diurnal: diurnal flux processing only implemented for bulk forcing" )  
    6969 
    70       CALL diurnal_sst_coolskin_step( qns, taum, rhop(:,:,1), rdt) 
     70      !                                            ! Cool skin 
     71      CALL diurnal_sst_coolskin_step( qns, taum, rhop(:,:,1), rn_Dt ) 
    7172 
    72       CALL iom_put( "sst_wl"   , x_dsst               )    ! warm layer (write out before update below). 
    73       CALL iom_put( "sst_cs"   , x_csdsst             )    ! cool skin 
     73      CALL iom_put( "sst_wl", x_dsst   )                 ! warm layer (write out before update below). 
     74      CALL iom_put( "sst_cs", x_csdsst )                 ! cool skin 
    7475 
    75       ! Diurnal warm layer model        
    76       CALL diurnal_sst_takaya_step( kstp, &  
    77       &    qsr, qns, taum, rhop(:,:,1), rdt)  
     76      !                                            ! Diurnal warm layer model        
     77      CALL diurnal_sst_takaya_step( kstp, qsr, qns, taum, rhop(:,:,1), rn_Dt )  
    7878 
    7979      IF( ln_diurnal_only ) THEN 
    80          IF( ln_diaobs )         CALL dia_obs( kstp )         ! obs-minus-model (assimilation) diagnostics (call after dynamics update) 
     80         IF( ln_diaobs )   CALL dia_obs( kstp )    ! obs-minus-model (assimilation) diagnostics (call after dynamics update) 
    8181      
    8282         !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>  
     
    8484         !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<  
    8585         IF( kstp == nit000   )   CALL iom_close( numror )     ! close input  ocean restart file  
    86          IF( lrst_oce         )   CALL rst_write    ( kstp )   ! write output ocean restart file 
     86         IF( lrst_oce         )   CALL rst_write( kstp   )     ! write output ocean restart file 
    8787      
    8888         IF( ln_timing .AND.  kstp == nit000  )   CALL timing_reset  
     
    9191   END SUBROUTINE stp_diurnal   
    9292    
     93   !!====================================================================== 
    9394END MODULE step_diu 
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/daymod.F90

    r9598 r9939  
    2020   !!                    ------------------------------- 
    2121   !!   sbcmod assume that the time step is dividing the number of second of  
    22    !!   in a day, i.e. ===> MOD( rday, rdt ) == 0  
     22   !!   in a day, i.e. ===> MOD( rday, rn_Dt ) == 0  
    2323   !!   except when user defined forcing is used (see sbcmod.F90) 
    2424   !!---------------------------------------------------------------------- 
     
    7272      ! 
    7373      ! max number of seconds between each restart 
    74       IF( REAL( nitend - nit000 + 1 ) * rdt > REAL( HUGE( nsec1jan000 ) ) ) THEN 
     74      IF( REAL( nitend - nit000 + 1 ) * rn_Dt > REAL( HUGE( nsec1jan000 ) ) ) THEN 
    7575         CALL ctl_stop( 'The number of seconds between each restart exceeds the integer 4 max value: 2^31-1. ',   & 
    7676            &           'You must do a restart at higher frequency (or remove this stop and recompile the code in I8)' ) 
    7777      ENDIF 
    78       nsecd   = NINT( rday       ) 
    79       nsecd05 = NINT( 0.5 * rday ) 
    80       ndt     = NINT(       rdt ) 
    81       ndt05   = NINT( 0.5 * rdt ) 
     78      nsecd   = NINT( rday         ) 
     79      nsecd05 = NINT( 0.5 * rday   ) 
     80      ndt     = NINT(       rn_Dt ) 
     81      ndt05   = NINT( 0.5 * rn_Dt ) 
    8282 
    8383      IF( .NOT. l_offline )   CALL day_rst( nit000, 'READ' ) 
     
    239239      nsec_week  = nsec_week  + ndt 
    240240      nsec_day   = nsec_day   + ndt 
    241       adatrj  = adatrj  + rdt / rday 
    242       fjulday = fjulday + rdt / rday 
     241      adatrj  = adatrj  + rn_Dt / rday 
     242      fjulday = fjulday + rn_Dt / rday 
    243243      IF( ABS(fjulday - REAL(NINT(fjulday),wp)) < zprec )   fjulday = REAL(NINT(fjulday),wp)   ! avoid truncation error 
    244244      IF( ABS(adatrj  - REAL(NINT(adatrj ),wp)) < zprec )   adatrj  = REAL(NINT(adatrj ),wp)   ! avoid truncation error 
     
    309309      !!       In both those options, the  exact duration of the experiment 
    310310      !!       since the beginning (cumulated duration of all previous restart runs) 
    311       !!       is not stored in the restart and is assumed to be (nit000-1)*rdt. 
     311      !!       is not stored in the restart and is assumed to be (nit000-1)*rn_Dt. 
    312312      !!       This is valid is the time step has remained constant. 
    313313      !! 
     
    378378               nminute = ( nn_time0 - nhour * 100 ) 
    379379               IF( nhour*3600+nminute*60-ndt05 .lt. 0 )  ndastp=ndastp-1      ! Start hour is specified in the namelist (default 0) 
    380                adatrj = ( REAL( nit000-1, wp ) * rdt ) / rday 
     380               adatrj = ( REAL( nit000-1, wp ) * rn_Dt ) / rday 
    381381               ! note this is wrong if time step has changed during run 
    382382            ENDIF 
     
    387387       nminute = ( nn_time0 - nhour * 100 ) 
    388388            IF( nhour*3600+nminute*60-ndt05 .lt. 0 )  ndastp=ndastp-1      ! Start hour is specified in the namelist (default 0) 
    389             adatrj = ( REAL( nit000-1, wp ) * rdt ) / rday 
     389            adatrj = ( REAL( nit000-1, wp ) * rn_Dt ) / rday 
    390390         ENDIF 
    391391         IF( ABS(adatrj  - REAL(NINT(adatrj),wp)) < 0.1 / rday )   adatrj = REAL(NINT(adatrj),wp)   ! avoid truncation error 
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/dom_oce.F90

    r9667 r9939  
    3333   LOGICAL , PUBLIC ::   ln_meshmask    !: =T  create a mesh-mask file (mesh_mask.nc) 
    3434   REAL(wp), PUBLIC ::   rn_isfhmin     !: threshold to discriminate grounded ice to floating ice 
    35    REAL(wp), PUBLIC ::   rn_rdt         !: time step for the dynamics and tracer 
     35   REAL(wp), PUBLIC ::   rn_dt          !: time step for the dynamics and tracer 
    3636   REAL(wp), PUBLIC ::   rn_atfp        !: asselin time filter parameter 
    37    INTEGER , PUBLIC ::   nn_euler       !: =0 start with forward time step or not (=1) 
     37   LOGICAL , PUBLIC ::   ln_1st_euler   !: =0 start with forward time step or not (=1) 
    3838   LOGICAL , PUBLIC ::   ln_iscpl       !: coupling with ice sheet 
    3939   LOGICAL , PUBLIC ::   ln_crs         !: Apply grid coarsening to dynamical model output or online passive tracers 
     
    5050   LOGICAL,  PUBLIC :: ln_bt_auto       !: Set number of barotropic iterations automatically 
    5151   INTEGER,  PUBLIC :: nn_bt_flt        !: Filter choice 
    52    INTEGER,  PUBLIC :: nn_baro          !: Number of barotropic iterations during one baroclinic step (rdt) 
     52   INTEGER,  PUBLIC :: nn_e             !: Number of external mode sub-step used at each ocean time-step 
    5353   REAL(wp), PUBLIC :: rn_bt_cmax       !: Maximum allowed courant number (used if ln_bt_auto=T) 
    5454   REAL(wp), PUBLIC :: rn_bt_alpha      !: Time stepping diffusion parameter 
    5555 
    56  
    57    !                                   !! old non-DOCTOR names still used in the model 
    58    REAL(wp), PUBLIC ::   atfp           !: asselin time filter parameter 
    59    REAL(wp), PUBLIC ::   rdt            !: time step for the dynamics and tracer 
    60  
    6156   !                                   !!! associated variables 
    62    INTEGER , PUBLIC ::   neuler         !: restart euler forward option (0=Euler) 
    63    REAL(wp), PUBLIC ::   r2dt           !: = 2*rdt except at nit000 (=rdt) if neuler=0 
     57   LOGICAL , PUBLIC ::   l_1st_euler    !: Euler 1st time-step flag (=T if ln_restart=F or ln_1st_euler=T) 
     58   REAL(wp), PUBLIC ::   rDt, r1_Dt     !: MLF: = 2*rn_Dt and 1/(2*rn_Dt) except if l_1st_euler=T where half the value is used 
     59   !                                    !  RK3: = rn_Dt 
    6460 
    6561   !!---------------------------------------------------------------------- 
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/domain.F90

    r9598 r9939  
    288288      INTEGER  ::   ios   ! Local integer 
    289289      ! 
    290       NAMELIST/namrun/ cn_ocerst_indir, cn_ocerst_outdir, nn_stocklist, ln_rst_list,                 & 
    291          &             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 ,     & 
    293          &             nn_stock, nn_write , ln_mskland  , ln_clobber   , nn_chunksz, nn_euler  ,     & 
     290      NAMELIST/namrun/ cn_ocerst_indir, cn_ocerst_outdir, nn_stocklist, ln_rst_list,                   & 
     291         &             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   ,     & 
     293         &             nn_stock, nn_write , ln_mskland  , ln_clobber   , nn_chunksz, ln_1st_euler,     & 
    294294         &             ln_cfmeta, ln_iscpl, ln_xios_read, nn_wxios 
    295       NAMELIST/namdom/ ln_linssh, rn_isfhmin, rn_rdt, rn_atfp, ln_crs, ln_meshmask 
     295      NAMELIST/namdom/ ln_linssh, rn_isfhmin, rn_Dt, rn_atfp, ln_crs, ln_meshmask 
    296296#if defined key_netcdf4 
    297297      NAMELIST/namnc4/ nn_nchunks_i, nn_nchunks_j, nn_nchunks_k, ln_nc4zip 
     
    323323         WRITE(numout,*) '      restart output directory        cn_ocerst_outdir= ', TRIM( cn_ocerst_outdir ) 
    324324         WRITE(numout,*) '      restart logical                 ln_rstart       = ', ln_rstart 
    325          WRITE(numout,*) '      start with forward time step    nn_euler        = ', nn_euler 
     325         WRITE(numout,*) '      start with forward time step    ln_1st_euler    = ', ln_1st_euler 
    326326         WRITE(numout,*) '      control of time step            nn_rstctl       = ', nn_rstctl 
    327327         WRITE(numout,*) '      number of the first time step   nn_it000        = ', nn_it000 
     
    361361      nstocklist = nn_stocklist 
    362362      nwrite = nn_write 
    363       neuler = nn_euler 
    364       IF( neuler == 1 .AND. .NOT. ln_rstart ) THEN 
     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 
    365366         IF(lwp) WRITE(numout,*)   
    366367         IF(lwp) WRITE(numout,*)'   ==>>>   Start from rest (ln_rstart=F)' 
    367          IF(lwp) WRITE(numout,*)'           an Euler initial time step is used : nn_euler is forced to 0 '    
    368          neuler = 0 
     368         IF(lwp) WRITE(numout,*)'           an Euler initial time step is used '    
     369         l_1st_euler = .TRUE. 
    369370      ENDIF 
    370371      !                             ! control of output frequency 
     
    374375         nstock = nitend 
    375376      ENDIF 
    376       IF ( nwrite == 0 ) THEN 
     377      IF( nwrite == 0 ) THEN 
    377378         WRITE(ctmp1,*) 'nwrite = ', nwrite, ' it is forced to ', nitend 
    378379         CALL ctl_warn( ctmp1 ) 
     
    413414         WRITE(numout,*) '      create mesh/mask file                   ln_meshmask = ', ln_meshmask 
    414415         WRITE(numout,*) '      treshold to open the isf cavity         rn_isfhmin  = ', rn_isfhmin, ' [m]' 
    415          WRITE(numout,*) '      ocean time step                         rn_rdt      = ', rn_rdt 
     416         WRITE(numout,*) '      ocean time step                         rn_dt       = ', rn_dt 
    416417         WRITE(numout,*) '      asselin time filter parameter           rn_atfp     = ', rn_atfp 
    417418         WRITE(numout,*) '      online coarsening of dynamical fields   ln_crs      = ', ln_crs 
    418419      ENDIF 
    419420      ! 
    420       !          ! conversion DOCTOR names into model names (this should disappear soon) 
    421       atfp = rn_atfp 
    422       rdt  = rn_rdt 
    423  
    424421      IF( TRIM(Agrif_CFixed()) == '0' ) THEN 
    425422         lrxios = ln_xios_read.AND.ln_rstart 
  • NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/domvvl.F90

    r9598 r9939  
    5454   LOGICAL , PUBLIC :: ln_vvl_dbg = .FALSE.                ! debug control prints 
    5555 
    56    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td                ! thickness diffusion transport 
    57    REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdiv_lf                     ! low frequency part of hz divergence 
    58    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_b, tilde_e3t_n    ! baroclinic scale factors 
    59    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_a, dtilde_e3t_a   ! baroclinic scale factors 
    60    REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_e3t                 ! retoring period for scale factors 
    61    REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_hdv                 ! retoring period for low freq. divergence 
     56   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td      ! thickness diffusion transport 
     57   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdiv_lf           ! low frequency part of hz divergence 
     58   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: te3t_b, te3t_n    ! baroclinic scale factors 
     59   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: te3t_a, dte3t_a   ! baroclinic scale factors 
     60   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_e3t       ! retoring period for scale factors 
     61   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_hdv       ! retoring period for low freq. divergence 
    6262 
    6363   !! * Substitutions 
     
    7676      IF( ln_vvl_zstar )   dom_vvl_alloc = 0 
    7777      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 
    78          ALLOCATE( tilde_e3t_b(jpi,jpj,jpk)  , tilde_e3t_n(jpi,jpj,jpk) , tilde_e3t_a(jpi,jpj,jpk) ,   & 
    79             &      dtilde_e3t_a(jpi,jpj,jpk) , un_td  (jpi,jpj,jpk)     , vn_td  (jpi,jpj,jpk)     ,   & 
     78         ALLOCATE( te3t_b(jpi,jpj,jpk)  , te3t_n(jpi,jpj,jpk) , te3t_a(jpi,jpj,jpk) ,   & 
     79            &      dte3t_a(jpi,jpj,jpk) , un_td  (jpi,jpj,jpk)     , vn_td  (jpi,jpj,jpk)     ,   & 
    8080            &      STAT = dom_vvl_alloc        ) 
    8181         IF( lk_mpp             )   CALL mpp_sum ( dom_vvl_alloc ) 
     
    103103      !!               - interpolate scale factors 
    104104      !! 
    105       !! ** Action  : - e3t_(n/b) and tilde_e3t_(n/b) 
     105      !! ** Action  : - e3t_(n/b) and te3t_(n/b) 
    106106      !!              - Regrid: e3(u/v)_n 
    107107      !!                        e3(u/v)_b        
     
    117117      INTEGER ::   ji, jj, jk 
    118118      INTEGER ::   ii0, ii1, ij0, ij1 
    119       REAL(wp)::   zcoef 
     119      REAL(wp)::   zcoef, z1_Dt 
    120120      !!---------------------------------------------------------------------- 
    121121      ! 
     
    129129      IF( dom_vvl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' ) 
    130130      ! 
    131       !                    ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf 
     131      !                    ! Read or initialize e3t_(b/n), te3t_(b/n) and hdiv_lf 
    132132      CALL dom_vvl_rst( nit000, 'READ' ) 
    133133      e3t_a(:,:,jpk) = e3t_0(:,:,jpk)  ! last level always inside the sea floor set one for all 
     
    208208         IF( ln_vvl_ztilde_as_zstar ) THEN   ! z-star emulation using z-tile 
    209209            frq_rst_e3t(:,:) = 0._wp               !Ignore namelist settings 
    210             frq_rst_hdv(:,:) = 1._wp / rdt 
     210            frq_rst_hdv(:,:) = 1._wp / rn_Dt 
    211211         ENDIF 
    212212         IF ( ln_vvl_zstar_at_eqtor ) THEN   ! use z-star in vicinity of the Equator 
     213            z1_Dt = 1._wp / rn_Dt 
    213214            DO jj = 1, jpj 
    214215               DO ji = 1, jpi 
     
    216217                  IF( ABS(gphit(ji,jj)) >= 6.) THEN 
    217218                     ! values outside the equatorial band and transition zone (ztilde) 
    218                      frq_rst_e3t(ji,jj) =  2.0_wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400.e0_wp ) 
    219                      frq_rst_hdv(ji,jj) =  2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp ) 
     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 ) 
    220221                  ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN    ! Equator strip ==> z-star 
    221222                     ! values inside the equatorial band (ztilde as zstar) 
    222                      frq_rst_e3t(ji,jj) =  0.0_wp 
    223                      frq_rst_hdv(ji,jj) =  1.0_wp / rdt 
     223                     frq_rst_e3t(ji,jj) =  0._wp 
     224                     frq_rst_hdv(ji,jj) =  z1_Dt 
    224225                  ELSE                                      ! transition band (2.5 to 6 degrees N/S) 
    225226                     !                                      ! (linearly transition from z-tilde to z-star) 
    226                      frq_rst_e3t(ji,jj) = 0.0_wp + (frq_rst_e3t(ji,jj)-0.0_wp)*0.5_wp   & 
    227                         &            * (  1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  & 
    228                         &                                          * 180._wp / 3.5_wp ) ) 
    229                      frq_rst_hdv(ji,jj) = (1.0_wp / rdt)                                & 
    230                         &            + (  frq_rst_hdv(ji,jj)-(1.e0_wp / rdt) )*0.5_wp   & 
    231                         &            * (  1._wp  - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  & 
    232                         &                                          * 180._wp / 3.5_wp ) ) 
     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 )  ) 
    233231                  ENDIF 
    234232               END DO 
     
    237235               ii0 = 103   ;   ii1 = 111        
    238236               ij0 = 128   ;   ij1 = 135   ;    
    239                frq_rst_e3t( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  0.0_wp 
    240                frq_rst_hdv( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  1.e0_wp / rdt 
     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 
    241239            ENDIF 
    242240         ENDIF 
     
    280278      !! 
    281279      !! ** Action  :  - hdiv_lf    : restoring towards full baroclinic divergence in z_tilde case 
    282       !!               - tilde_e3t_a: after increment of vertical scale factor  
     280      !!               - te3t_a: after increment of vertical scale factor  
    283281      !!                              in z_tilde case 
    284282      !!               - e3(t/u/v)_a 
     
    345343            IF( kt > nit000 ) THEN 
    346344               DO jk = 1, jpkm1 
    347                   hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:)   & 
     345                  hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rn_Dt * frq_rst_hdv(:,:)   & 
    348346                     &          * ( hdiv_lf(:,:,jk) - e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) ) 
    349347               END DO 
     
    353351         ! II - after z_tilde increments of vertical scale factors 
    354352         ! ======================================================= 
    355          tilde_e3t_a(:,:,:) = 0._wp  ! tilde_e3t_a used to store tendency terms 
     353         te3t_a(:,:,:) = 0._wp  ! te3t_a used to store tendency terms 
    356354 
    357355         ! 1 - High frequency divergence term 
     
    359357         IF( ln_vvl_ztilde ) THEN     ! z_tilde case 
    360358            DO jk = 1, jpkm1 
    361                tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) ) 
     359               te3t_a(:,:,jk) = te3t_a(:,:,jk) - ( e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) ) 
    362360            END DO 
    363361         ELSE                         ! layer case 
    364362            DO jk = 1, jpkm1 
    365                tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) -   e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk) 
     363               te3t_a(:,:,jk) = te3t_a(:,:,jk) -   e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk) 
    366364            END DO 
    367365         ENDIF 
     
    371369         IF( ln_vvl_ztilde ) THEN 
    372370            DO jk = 1, jpk 
    373                tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - frq_rst_e3t(:,:) * tilde_e3t_b(:,:,jk) 
     371               te3t_a(:,:,jk) = te3t_a(:,:,jk) - frq_rst_e3t(:,:) * te3t_b(:,:,jk) 
    374372            END DO 
    375373         ENDIF 
     
    383381               DO ji = 1, fs_jpim1   ! vector opt. 
    384382                  un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj)           & 
    385                      &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) ) 
     383                     &            * ( te3t_b(ji,jj,jk) - te3t_b(ji+1,jj  ,jk) ) 
    386384                  vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj)           &  
    387                      &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) ) 
     385                     &            * ( te3t_b(ji,jj,jk) - te3t_b(ji  ,jj+1,jk) ) 
    388386                  zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk) 
    389387                  zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk) 
     
    400398            DO jj = 2, jpjm1 
    401399               DO ji = fs_2, fs_jpim1   ! vector opt. 
    402                   tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + (   un_td(ji-1,jj  ,jk) - un_td(ji,jj,jk)    & 
     400                  te3t_a(ji,jj,jk) = te3t_a(ji,jj,jk) + (   un_td(ji-1,jj  ,jk) - un_td(ji,jj,jk)    & 
    403401                     &                                          +     vn_td(ji  ,jj-1,jk) - vn_td(ji,jj,jk)    & 
    404402                     &                                            ) * r1_e1e2t(ji,jj) 
     
    414412         ! Leapfrog time stepping 
    415413         ! ~~~~~~~~~~~~~~~~~~~~~~ 
    416          IF( neuler == 0 .AND. kt == nit000 ) THEN 
    417             z2dt =  rdt 
    418          ELSE 
    419             z2dt = 2.0_wp * rdt 
    420          ENDIF 
    421          CALL lbc_lnk( tilde_e3t_a(:,:,:), 'T', 1._wp ) 
    422          tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + z2dt * tmask(:,:,:) * tilde_e3t_a(:,:,:) 
     414         CALL lbc_lnk( te3t_a(:,:,:), 'T', 1._wp ) 
     415         te3t_a(:,:,:) = te3t_b(:,:,:) + z2dt * tmask(:,:,:) * te3t_a(:,:,:) 
    423416 
    424417         ! Maximum deformation control 
     
    426419         ze3t(:,:,jpk) = 0._wp 
    427420         DO jk = 1, jpkm1 
    428