Changeset 8813


Ignore:
Timestamp:
2017-11-24T17:56:51+01:00 (3 years ago)
Author:
gm
Message:

#1911 (ENHANCE-09): PART I.3 - phasing with updated branch dev_r8183_ICEMODEL revision 8787

Location:
branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM
Files:
20 edited

Legend:

Unmodified
Added
Removed
  • branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref

    r8637 r8813  
    4141&namitd         !   Ice discretization 
    4242!------------------------------------------------------------------------------ 
    43    rn_himean        =   2.0           !  expected domain-average ice thickness (m) 
     43   ln_cat_hfn       = .true.          !  ice categories are defined by a function following rn_himean**(-0.05) 
     44      rn_himean     =   2.0           !  expected domain-average ice thickness (m) 
     45   ln_cat_usr       = .false.         !  ice categories are defined by rn_catbnd below (m) 
     46      rn_catbnd     =   0.,0.45,1.1,2.1,3.7,6.0   
    4447   rn_himin         =   0.1           !  minimum ice thickness (m) used in remapping 
    4548/ 
     
    9093!------------------------------------------------------------------------------ 
    9194   ln_rhg_EVP       = .true.          !  EVP rheology 
    92       rn_creepl     =   1.0e-12       !     creep limit [1/s] 
     95      ln_aEVP       = .false.         !     adaptive rheology (Kimmritz et al. 2016 & 2017) 
     96      rn_creepl     =   2.0e-9        !     creep limit [1/s] 
    9397      rn_ecc        =   2.0           !     eccentricity of the elliptical yield curve           
    9498      nn_nevp       = 120             !     number of EVP subcycles                              
     
    118122                                      !     = 2  Redistribute a single flux over categories 
    119123                                      !          ==> coupled mode only 
     124   nice_jules       = 1               !  Jules coupling (0=OFF, 1=EMULATED, 2=ACTIVE) 
    120125/ 
    121126!------------------------------------------------------------------------------ 
     
    136141                                      !     Obs: 0.1-0.5 (Lecomte et al, JAMES 2013) 
    137142   rn_kappa_i       =   1.0           !  radiation attenuation coefficient in sea ice [1/m] 
    138    ln_dqns_i        = .true.          !  change the surface non-solar flux with surface temperature (T) or not (F) 
    139143/ 
    140144!------------------------------------------------------------------------------ 
  • branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/LIM_SRC_3/ice.F90

    r8637 r8813  
    178178   ! 
    179179   !                                     !!** ice-rheology namelist (namrhg) ** 
     180   LOGICAL , PUBLIC ::   ln_aEVP          !: using adaptive EVP (T or F)  
    180181   REAL(wp), PUBLIC ::   rn_creepl        !: creep limit : has to be under 1.0e-9 
    181182   REAL(wp), PUBLIC ::   rn_ecc           !: eccentricity of the elliptical yield curve 
     
    194195   !                                      !                                   using T-ice and albedo sensitivity 
    195196   !                                      !   = 2  Redistribute a single flux over categories 
     197 
     198   INTEGER , PUBLIC            ::   nice_jules            !: Choice of jules coupling  
     199   !                                                      !  Associated indices 
     200   INTEGER , PUBLIC, PARAMETER ::   np_jules_OFF    = 0   !: no Jules coupling (ice thermodynamics forced via qsr and qns) 
     201   INTEGER , PUBLIC, PARAMETER ::   np_jules_EMULE  = 1   !: emulated Jules coupling via icethd_zdf.F90 (BL99) (1st round compute qcn and qsr_tr, 2nd round use it) 
     202   INTEGER , PUBLIC, PARAMETER ::   np_jules_ACTIVE = 2   !: active Jules coupling                      (SM0L) (compute qcn and qsr_tr via sbcblk.F90 or sbccpl.F90) 
    196203 
    197204   !                                     !!** ice-salinity namelist (namthd_sal) ** 
  • branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/LIM_SRC_3/ice1d.F90

    r8692 r8813  
    3333   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   ftr_ice_1d    
    3434   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qsr_ice_1d   
    35    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fr1_i0_1d    
    36    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fr2_i0_1d    
    3735   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qns_ice_1d   
    3836   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   t_bo_1d      
    3937   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   rn_amax_1d 
     38    
     39   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qml_ice_1d     !: heat available for snow / ice surface melting [W/m2]  
     40   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qcn_ice_1d     !: heat available for snow / ice surface sublimation¬ [W/m2]  
     41   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qsr_ice_tr_1d  !: solar flux transmitted below the ice surface¬ [W/m2]  
    4042 
    4143   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_sum_1d 
     
    9698   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   evap_ice_1d   !: <==> the 2D  evap_ice 
    9799   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qprec_ice_1d  !: <==> the 2D  qprec_ice 
    98    !                                                     ! to reintegrate longwave flux inside the ice thermodynamics 
    99    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   i0            !: fraction of radiation transmitted to the ice 
    100100 
    101101   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   t_su_1d       !: <==> the 2D  t_su 
     
    178178      ALLOCATE( nptidx   (jpij) ,   & 
    179179         &      qlead_1d (jpij) , ftr_ice_1d(jpij) , qsr_ice_1d (jpij) ,   & 
    180          &      fr1_i0_1d(jpij) , fr2_i0_1d (jpij) , qns_ice_1d(jpij)  ,   & 
     180         &      qns_ice_1d(jpij)  ,   & 
     181         &      qml_ice_1d(jpij), qcn_ice_1d(jpij) , qsr_ice_tr_1d(jpij) , & 
    181182         &      t_bo_1d   (jpij) ,                                         & 
    182183         &      hfx_sum_1d(jpij) , hfx_bom_1d(jpij) ,hfx_bog_1d(jpij) ,    &  
     
    194195         &      wfx_snw_sub_1d(jpij), wfx_snw_dyn_1d(jpij), wfx_ice_sub_1d(jpij), wfx_err_sub_1d(jpij) ,   & 
    195196         &      wfx_lam_1d(jpij)  , wfx_dyn_1d(jpij), wfx_pnd_1d(jpij),dqns_ice_1d(jpij) , evap_ice_1d (jpij),   & 
    196          &      qprec_ice_1d(jpij), i0         (jpij) ,                                         &   
     197         &      qprec_ice_1d(jpij),                                                             &   
    197198         &      sfx_bri_1d (jpij) , sfx_bog_1d (jpij) , sfx_bom_1d (jpij) , sfx_sum_1d (jpij),  & 
    198199         &      sfx_sni_1d (jpij) , sfx_opw_1d (jpij) , sfx_res_1d (jpij) , sfx_sub_1d (jpij),  & 
     
    200201      ! 
    201202      ii = ii + 1 
    202       ALLOCATE( t_su_1d  (jpij) , t_si_1d   (jpij) , a_i_1d  (jpij) , a_ib_1d(jpij) ,                & 
    203          &      h_i_1d  (jpij) , h_ib_1d  (jpij) , h_s_1d (jpij) , fc_su  (jpij) , fc_bo_i(jpij) ,   &     
    204          &      dh_s_tot (jpij) , dh_i_surf (jpij) , dh_i_sub(jpij) ,                                &     
     203      ALLOCATE( t_su_1d  (jpij) , t_si_1d   (jpij) , a_i_1d  (jpij) , a_ib_1d(jpij) ,                  & 
     204         &      h_i_1d  (jpij) , h_ib_1d  (jpij) , h_s_1d (jpij) , fc_su  (jpij) , fc_bo_i(jpij) ,  &     
     205         &      dh_s_tot (jpij) , dh_i_surf (jpij) , dh_i_sub(jpij) ,                                  &     
    205206         &      dh_i_bott(jpij) , dh_s_mlt(jpij), dh_snowice(jpij) , s_i_1d (jpij) , s_i_new(jpij) , & 
    206207         &      a_ip_1d  (jpij) , v_ip_1d   (jpij) , v_i_1d  (jpij) , v_s_1d (jpij) , & 
     
    210211      ii = ii + 1 
    211212      ALLOCATE( t_s_1d  (jpij,nlay_s)     , t_i_1d (jpij,nlay_i)     , sz_i_1d(jpij,nlay_i) ,  &             
    212          &      e_i_1d  (jpij,nlay_i)     , e_s_1d (jpij,nlay_s)     ,                         & 
     213         &      e_i_1d  (jpij,nlay_i)     , e_s_1d (jpij,nlay_s)     ,                        & 
    213214         &      eh_i_old(jpij,0:nlay_i+1) , h_i_old(jpij,0:nlay_i+1) , STAT=ierr(ii) ) 
    214215      ! 
     
    220221      ! 
    221222      ii = ii + 1 
    222       ALLOCATE( a_i_2d (jpij,jpl) , a_ib_2d(jpij,jpl) , h_i_2d (jpij,jpl) , h_ib_2d(jpij,jpl) ,   & 
    223          &      v_i_2d (jpij,jpl) , v_s_2d (jpij,jpl) , oa_i_2d(jpij,jpl) , sv_i_2d(jpij,jpl) ,   & 
    224          &      a_ip_2d(jpij,jpl) , v_ip_2d(jpij,jpl) , t_su_2d(jpij,jpl) ,                   STAT=ierr(ii) ) 
     223      ALLOCATE( a_i_2d(jpij,jpl) , a_ib_2d(jpij,jpl) , h_i_2d(jpij,jpl) , h_ib_2d(jpij,jpl) , & 
     224         &      v_i_2d(jpij,jpl) ,v_s_2d(jpij,jpl) ,oa_i_2d(jpij,jpl) ,sv_i_2d(jpij,jpl) ,  & 
     225         &      a_ip_2d(jpij,jpl) ,v_ip_2d(jpij,jpl) ,t_su_2d(jpij,jpl) ,  & 
     226         &      STAT=ierr(ii) ) 
    225227 
    226228      ice1D_alloc = MAXVAL( ierr(:) ) 
  • branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/LIM_SRC_3/icedyn_rhg.F90

    r8586 r8813  
    7878      CASE( np_rhgEVP )                ! Elasto-Viscous-Plastic ! 
    7979         !                             !------------------------! 
    80          CALL ice_dyn_rhg_evp( kt, stress1_i, stress2_i, stress12_i, u_ice, v_ice, shear_i, divu_i, delta_i ) 
     80         CALL ice_dyn_rhg_evp( kt, stress1_i, stress2_i, stress12_i, shear_i, divu_i, delta_i ) 
    8181         !          
    8282      END SELECT 
     
    107107      INTEGER ::   ios, ioptio   ! Local integer output status for namelist read 
    108108      !! 
    109       NAMELIST/namdyn_rhg/  ln_rhg_EVP, rn_creepl, rn_ecc , nn_nevp, rn_relast 
     109      NAMELIST/namdyn_rhg/  ln_rhg_EVP, ln_aEVP, rn_creepl, rn_ecc , nn_nevp, rn_relast 
    110110      !!------------------------------------------------------------------- 
    111111      ! 
     
    125125         WRITE(numout,*) '   Namelist namdyn_rhg:' 
    126126         WRITE(numout,*) '      rheology EVP (icedyn_rhg_evp)                        ln_rhg_EVP = ', ln_rhg_EVP 
     127         WRITE(numout,*) '         use adaptive EVP (aEVP)                           ln_aEVP    = ', ln_aEVP 
    127128         WRITE(numout,*) '         creep limit                                       rn_creepl  = ', rn_creepl 
    128129         WRITE(numout,*) '         eccentricity of the elliptical yield curve        rn_ecc     = ', rn_ecc 
     
    137138      IF( ioptio /= 1 )   CALL ctl_stop( 'ice_dyn_rhg_init: choose one and only one ice rheology' ) 
    138139      ! 
    139       IF( ln_rhg_EVP )   CALL rhg_evp_rst( 'READ' )  !* read or initialize all required files 
     140      IF( ln_rhg_EVP  )   CALL rhg_evp_rst( 'READ' )  !* read or initialize all required files 
    140141      ! 
    141142   END SUBROUTINE ice_dyn_rhg_init 
  • branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/LIM_SRC_3/icedyn_rhg_evp.F90

    r8637 r8813  
    1616   !!   'key_lim3'                                       ESIM sea-ice model 
    1717   !!---------------------------------------------------------------------- 
    18    !!   ice_dyn_rhg_evp       : computes ice velocities from EVP rheology 
     18   !!   ice_dyn_rhg_evp : computes ice velocities from EVP rheology 
     19   !!   rhg_evp_rst     : read/write EVP fields in ice restart 
    1920   !!---------------------------------------------------------------------- 
    2021   USE phycst         ! Physical constant 
     
    2425   USE ice            ! sea-ice: ice variables 
    2526   USE icedyn_rdgrft  ! sea-ice: ice strength 
     27   USE bdy_oce , ONLY : ln_bdy  
     28   USE bdyice  
     29#if defined key_agrif 
     30   USE agrif_lim3_interp 
     31#endif 
    2632   ! 
    2733   USE in_out_manager ! I/O manager 
     
    3137   USE lbclnk         ! lateral boundary conditions (or mpp links) 
    3238   USE prtctl         ! Print control 
    33    ! 
    34 #if defined key_agrif 
    35    USE agrif_lim3_interp 
    36 #endif 
    37    USE bdy_oce   , ONLY: ln_bdy  
    38    USE bdyice  
    3939 
    4040   IMPLICIT NONE 
     
    5353CONTAINS 
    5454 
    55    SUBROUTINE ice_dyn_rhg_evp( kt, pstress1_i, pstress2_i, pstress12_i, u_ice, v_ice, pshear_i, pdivu_i, pdelta_i ) 
     55   SUBROUTINE ice_dyn_rhg_evp( kt, pstress1_i, pstress2_i, pstress12_i, pshear_i, pdivu_i, pdelta_i ) 
    5656      !!------------------------------------------------------------------- 
    5757      !!                 ***  SUBROUTINE ice_dyn_rhg_evp  *** 
    58       !!                          EVP-C-grid 
     58      !!                             EVP-C-grid 
    5959      !! 
    6060      !! ** purpose : determines sea ice drift from wind stress, ice-ocean 
     
    9595      !!              5) Diagnostics including charge ellipse 
    9696      !! 
    97       !! ** Notes   : There is the possibility to use mEVP from Bouillon 2013 
    98       !!              (by uncommenting some lines in part 3 and changing alpha and beta parameters) 
    99       !!              but this solution appears very unstable (see Kimmritz et al 2016) 
     97      !! ** Notes   : There is the possibility to use aEVP from the nice work of Kimmritz et al. (2016 & 2017) 
     98      !!              by setting up ln_aEVP=T (i.e. changing alpha and beta parameters). 
     99      !!              This is an upgraded version of mEVP from Bouillon et al. 2013 
     100      !!              (i.e. more stable and better convergence) 
    100101      !! 
    101102      !! References : Hunke and Dukowicz, JPO97 
    102103      !!              Bouillon et al., Ocean Modelling 2009 
    103104      !!              Bouillon et al., Ocean Modelling 2013 
     105      !!              Kimmritz et al., Ocean Modelling 2016 & 2017 
    104106      !!------------------------------------------------------------------- 
    105       INTEGER, INTENT(in) ::   kt      ! time step 
    106       REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pstress1_i, pstress2_i, pstress12_i 
    107       REAL(wp), DIMENSION(:,:), INTENT(inout) ::   u_ice, v_ice  
    108       REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pshear_i, pdivu_i, pdelta_i  
     107      INTEGER                 , INTENT(in   ) ::   kt                                    ! time step 
     108      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pstress1_i, pstress2_i, pstress12_i   ! 
     109      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pshear_i  , pdivu_i   , pdelta_i      ! 
    109110      !! 
    110111      INTEGER ::   ji, jj       ! dummy loop indices 
    111112      INTEGER ::   jter         ! local integers 
    112  
     113      ! 
    113114      REAL(wp) ::   zrhoco                                                   ! rau0 * rn_cio 
    114115      REAL(wp) ::   zdtevp, z1_dtevp                                         ! time step for subcycling 
    115116      REAL(wp) ::   ecc2, z1_ecc2                                            ! square of yield ellipse eccenticity 
    116       REAL(wp) ::   zbeta, zalph1, z1_alph1, zalph2, z1_alph2                ! alpha and beta from Bouillon 2009 and 2013 
     117      REAL(wp) ::   zalph1, z1_alph1, zalph2, z1_alph2                       ! alpha coef from Bouillon 2009 or Kimmritz 2017 
    117118      REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV                            ! ice/snow mass 
    118119      REAL(wp) ::   zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2            ! temporary scalars 
    119120      REAL(wp) ::   zTauO, zTauB, zTauE, zvel                                ! temporary scalars 
    120  
     121      ! 
    121122      REAL(wp) ::   zresm                                                    ! Maximal error on ice velocity 
    122123      REAL(wp) ::   zintb, zintn                                             ! dummy argument 
    123124      REAL(wp) ::   zfac_x, zfac_y 
    124125      REAL(wp) ::   zshear, zdum1, zdum2 
    125        
     126      ! 
    126127      REAL(wp), DIMENSION(jpi,jpj) ::   z1_e1t0, z1_e2t0                ! scale factors 
    127128      REAL(wp), DIMENSION(jpi,jpj) ::   zp_delt                         ! P/delta at T points 
    128       ! 
     129      REAL(wp), DIMENSION(jpi,jpj) ::   zbeta                           ! beta coef from Kimmritz 2017 
     130      ! 
     131      REAL(wp), DIMENSION(jpi,jpj) ::   zdt_m                           ! (dt / ice-snow_mass) on T points 
    129132      REAL(wp), DIMENSION(jpi,jpj) ::   zaU   , zaV                     ! ice fraction on U/V points 
    130       REAL(wp), DIMENSION(jpi,jpj) ::   zmU_t, zmV_t                    ! ice/snow mass/dt on U/V points 
     133      REAL(wp), DIMENSION(jpi,jpj) ::   zmU_t, zmV_t                    ! (ice-snow_mass / dt) on U/V points 
    131134      REAL(wp), DIMENSION(jpi,jpj) ::   zmf                             ! coriolis parameter at T points 
    132135      REAL(wp), DIMENSION(jpi,jpj) ::   zTauU_ia , ztauV_ia             ! ice-atm. stress at U-V points 
     
    134137      REAL(wp), DIMENSION(jpi,jpj) ::   v_oceU, u_oceV, v_iceU, u_iceV  ! ocean/ice u/v component on V/U points                            
    135138      REAL(wp), DIMENSION(jpi,jpj) ::   zfU   , zfV                     ! internal stresses 
    136        
     139      ! 
    137140      REAL(wp), DIMENSION(jpi,jpj) ::   zds                             ! shear 
    138141      REAL(wp), DIMENSION(jpi,jpj) ::   zs1, zs2, zs12                  ! stress tensor components 
    139142      REAL(wp), DIMENSION(jpi,jpj) ::   zu_ice, zv_ice, zresr           ! check convergence 
    140143      REAL(wp), DIMENSION(jpi,jpj) ::   zpice                           ! array used for the calculation of ice surface slope: 
    141                                                                              !   ocean surface (ssh_m) if ice is not embedded 
    142                                                                              !   ice top surface if ice is embedded    
     144      !                                                                      !   ocean surface (ssh_m) if ice is not embedded 
     145      !                                                                       !   ice top surface if ice is embedded    
    143146      REAL(wp), DIMENSION(jpi,jpj) ::   zCorx, zCory                    ! Coriolis stress array 
    144147      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_oi, ztauy_oi              ! Ocean-to-ice stress array 
    145  
     148      ! 
    146149      REAL(wp), DIMENSION(jpi,jpj) ::   zswitchU, zswitchV              ! dummy arrays 
    147150      REAL(wp), DIMENSION(jpi,jpj) ::   zmaskU, zmaskV                  ! mask for ice presence 
     
    179182#endif 
    180183      ! 
     184!!gm for Clem:  OPTIMIZATION:  I think zfmask can be computed one for all at the initialization.... 
    181185      !------------------------------------------------------------------------------! 
    182186      ! 0) mask at F points for the ice 
     
    231235 
    232236      ! alpha parameters (Bouillon 2009) 
    233       zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp 
    234       zalph2 = zalph1 * z1_ecc2 
    235  
    236       ! alpha and beta parameters (Bouillon 2013) 
    237       !!zalph1 = 40. 
    238       !!zalph2 = 40. 
    239       !!zbeta  = 3000. 
    240       !!zbeta = REAL( nn_nevp )   ! close to classical EVP of Hunke (2001) 
    241  
    242       z1_alph1 = 1._wp / ( zalph1 + 1._wp ) 
    243       z1_alph2 = 1._wp / ( zalph2 + 1._wp ) 
    244  
     237      IF( .NOT. ln_aEVP ) THEN 
     238         zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp 
     239         zalph2 = zalph1 * z1_ecc2 
     240 
     241         z1_alph1 = 1._wp / ( zalph1 + 1._wp ) 
     242         z1_alph2 = 1._wp / ( zalph2 + 1._wp ) 
     243      ENDIF 
     244          
    245245      ! Initialise stress tensor  
    246246      zs1 (:,:) = pstress1_i (:,:)  
     
    266266      IF( ln_ice_embd ) THEN             !== embedded sea ice: compute representative ice top surface ==! 
    267267         !                                             
    268          ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1} 
    269          !                                               = (1/nn_fsbc)^2 * {SUM[n], n=0,nn_fsbc-1} 
     268         ! average interpolation coeff as used in dynspg = (1/nn_fsbc)   * {SUM[n/nn_fsbc], n=0,nn_fsbc-1} 
     269         !                                               = (1/nn_fsbc)^2 * {SUM[n]        , n=0,nn_fsbc-1} 
    270270         zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp      
    271271         ! 
    272          ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1} 
     272         ! average interpolation coeff as used in dynspg = (1/nn_fsbc)   *    {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1} 
    273273         !                                               = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1}) 
    274274         zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp 
     
    304304            zmf(ji,jj)      = zm1 * ff_t(ji,jj) 
    305305 
     306            ! dt/m at T points (for alpha and beta coefficients) 
     307            zdt_m(ji,jj)    = zdtevp / MAX( zm1, zmmin ) 
     308             
    306309            ! m/dt 
    307310            zmU_t(ji,jj)    = zmassU * z1_dtevp 
    308311            zmV_t(ji,jj)    = zmassV * z1_dtevp 
    309  
     312             
    310313            ! Drag ice-atm. 
    311314            zTauU_ia(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj) 
     
    326329         END DO 
    327330      END DO 
    328       CALL lbc_lnk( zmf, 'T', 1. ) 
     331      CALL lbc_lnk_multi( zmf, 'T', 1., zdt_m, 'T', 1. ) 
    329332      ! 
    330333      !------------------------------------------------------------------------------! 
     
    380383               ! P/delta at T points 
    381384               zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl ) 
     385 
     386               ! alpha & beta for aEVP 
     387               !   gamma = 0.5*P/(delta+creepl) * (c*pi)**2/Area * dt/m 
     388               !   alpha = beta = sqrt(4*gamma) 
     389               IF( ln_aEVP ) THEN 
     390                  zalph1   = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) ) 
     391                  z1_alph1 = 1._wp / ( zalph1 + 1._wp ) 
     392                  zalph2   = zalph1 
     393                  z1_alph2 = z1_alph1 
     394               ENDIF 
    382395                
    383396               ! stress at T points 
     
    392405            DO ji = 1, jpim1 
    393406 
     407               ! alpha & beta for aEVP 
     408               IF( ln_aEVP ) THEN 
     409                  zalph2   = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) ) 
     410                  z1_alph2 = 1._wp / ( zalph2 + 1._wp ) 
     411                  zbeta(ji,jj) = zalph2 
     412               ENDIF 
     413                
    394414               ! P/delta at F points 
    395415               zp_delf = 0.25_wp * ( zp_delt(ji,jj) + zp_delt(ji+1,jj) + zp_delt(ji,jj+1) + zp_delt(ji+1,jj+1) ) 
     
    411431                  &                    ) * 2._wp * r1_e1u(ji,jj)                                                              & 
    412432                  &                  ) * r1_e1e2u(ji,jj) 
    413                   ! 
    414                   !                !--- V points 
     433               ! 
     434               !                !--- V points 
    415435               zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)                                             & 
    416436                  &                  - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj)    & 
     
    419439                  &                    ) * 2._wp * r1_e2v(ji,jj)                                                              & 
    420440                  &                  ) * r1_e1e2v(ji,jj) 
    421                   ! 
    422                   !                !--- u_ice at V point 
     441               ! 
     442               !                !--- u_ice at V point 
    423443               u_iceV(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     & 
    424444                  &                     + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1) 
    425                   ! 
    426                   !                !--- v_ice at U point 
     445               ! 
     446               !                !--- v_ice at U point 
    427447               v_iceU(ji,jj) = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji+1,jj)     & 
    428448                  &                     + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1) 
    429                ! 
    430449            END DO 
    431450         END DO 
    432451         ! 
    433452         ! --- Computation of ice velocity --- ! 
    434          !  Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta are chosen wisely and large nn_nevp 
     453         !  Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta vary as in Kimmritz 2016 & 2017 
    435454         !  Bouillon et al. 2009 (eq 34-35) => stable 
    436455         IF( MOD(jter,2) == 0 ) THEN ! even iterations 
     
    438457            DO jj = 2, jpjm1 
    439458               DO ji = fs_2, fs_jpim1 
    440                      !                 !--- tau_io/(v_oce - v_ice) 
     459                  !                 !--- tau_io/(v_oce - v_ice) 
    441460                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  & 
    442461                     &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) ) 
    443                      !                 !--- Ocean-to-Ice stress 
     462                  !                 !--- Ocean-to-Ice stress 
    444463                  ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) ) 
    445                      ! 
    446                      !                 !--- tau_bottom/v_ice 
     464                  ! 
     465                  !                 !--- tau_bottom/v_ice 
    447466                  zvel  = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) ) 
    448467                  zTauB = - tau_icebfr(ji,jj) / zvel 
    449                      ! 
    450                      !                 !--- Coriolis at V-points (energy conserving formulation) 
     468                  ! 
     469                  !                 !--- Coriolis at V-points (energy conserving formulation) 
    451470                  zCory(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  & 
    452471                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  & 
    453472                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) ) 
    454                      ! 
    455                      !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
     473                  ! 
     474                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
    456475                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) 
    457                      ! 
    458                      !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction 
     476                  ! 
     477                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction 
    459478                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
    460                      ! 
    461                      !                 !--- ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 
     479                  ! 
     480                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 
     481                  v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )         & ! previous velocity 
     482                     &                                  + zTauE + zTauO * v_ice(ji,jj)                                            & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     483                     &                                  ) / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
     484                     &                 + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0 
     485                     &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                               & ! v_ice = v_oce if mass < zmmin 
     486                     &           ) * zmaskV(ji,jj) 
     487                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
    462488                  v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                             &  ! previous velocity 
    463489                     &                                     + zTauE + zTauO * v_ice(ji,jj)                             &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     
    466492                     &              ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin 
    467493                     &            ) * zmaskV(ji,jj) 
    468                      ! 
    469                   ! Bouillon 2013 
    470                   !!v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) )                  & 
    471                   !!   &           + zfV(ji,jj) + zCory(ji,jj) + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj)  & 
    472                   !!   &           ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj) 
    473                   ! 
     494                  ENDIF 
    474495               END DO 
    475496            END DO 
     
    483504            ! 
    484505            DO jj = 2, jpjm1 
    485                DO ji = fs_2, fs_jpim1 
    486                                 
    487                   ! tau_io/(u_oce - u_ice) 
     506               DO ji = fs_2, fs_jpim1           
     507                  !                 !--- tau_io/(u_oce - u_ice) 
    488508                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  & 
    489509                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) ) 
    490  
    491                   ! Ocean-to-Ice stress 
     510                  !                 !--- Ocean-to-Ice stress 
    492511                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) ) 
    493  
    494                   ! tau_bottom/u_ice 
     512                  ! 
     513                  !                 !--- tau_bottom/u_ice 
    495514                  zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) ) 
    496515                  zTauB = - tau_icebfr(ji,jj) / zvel 
    497  
    498                   ! Coriolis at U-points (energy conserving formulation) 
     516                  ! 
     517                  !                 !--- Coriolis at U-points (energy conserving formulation) 
    499518                  zCorx(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  & 
    500519                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  & 
    501520                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) ) 
    502                    
    503                   ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
     521                  ! 
     522                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
    504523                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj) 
    505  
    506                   ! landfast switch => 0 = static friction ; 1 = sliding friction 
     524                  ! 
     525                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction 
    507526                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
    508  
    509                   ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 
     527                  ! 
     528                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 
     529                  u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )         & ! previous velocity 
     530                     &                                     + zTauE + zTauO * u_ice(ji,jj)                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     531                     &                                  ) / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
     532                     &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )              & ! static friction => slow decrease to v=0 
     533                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                              & ! v_ice = v_oce if mass < zmmin  
     534                     &            ) * zmaskU(ji,jj) 
     535                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
    510536                  u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                             &  ! previous velocity 
    511537                     &                                     + zTauE + zTauO * u_ice(ji,jj)                             &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     
    514540                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin  
    515541                     &            ) * zmaskU(ji,jj) 
    516                   ! Bouillon 2013 
    517                   !!u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) )                  & 
    518                   !!   &           + zfU(ji,jj) + zCorx(ji,jj) + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj)  & 
    519                   !!   &           ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchU(ji,jj) 
     542                  ENDIF 
    520543               END DO 
    521544            END DO 
     
    532555            DO jj = 2, jpjm1 
    533556               DO ji = fs_2, fs_jpim1 
    534  
    535                   ! tau_io/(u_oce - u_ice) 
     557                  !                 !--- tau_io/(u_oce - u_ice) 
    536558                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  & 
    537559                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) ) 
    538  
    539                   ! Ocean-to-Ice stress 
     560                  !                 !--- Ocean-to-Ice stress 
    540561                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) ) 
    541  
    542                   ! tau_bottom/u_ice 
     562                  ! 
     563                  !                 !--- tau_bottom/u_ice 
    543564                  zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) ) 
    544565                  zTauB = - tau_icebfr(ji,jj) / zvel 
    545  
    546                   ! Coriolis at U-points (energy conserving formulation) 
     566                  ! 
     567                  !                 !--- Coriolis at U-points (energy conserving formulation) 
    547568                  zCorx(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  & 
    548569                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  & 
    549570                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) ) 
    550  
    551                   ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
     571                  ! 
     572                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
    552573                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj) 
    553  
    554                   ! landfast switch => 0 = static friction ; 1 = sliding friction 
     574                  ! 
     575                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction 
    555576                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
    556  
    557                   ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 
     577                  ! 
     578                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 
     579                  u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )         & ! previous velocity 
     580                     &                                     + zTauE + zTauO * u_ice(ji,jj)                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     581                     &                                  ) / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
     582                     &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )              & ! static friction => slow decrease to v=0 
     583                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                              & ! v_ice = v_oce if mass < zmmin  
     584                     &            ) * zmaskU(ji,jj) 
     585                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
    558586                  u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                             &  ! previous velocity 
    559587                     &                                     + zTauE + zTauO * u_ice(ji,jj)                             &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    560588                     &                                     ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )             &  ! m/dt + tau_io(only ice part) + landfast 
    561589                     &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0 
    562                      &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin 
     590                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin  
    563591                     &            ) * zmaskU(ji,jj) 
    564                   ! Bouillon 2013 
    565                   !!u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) )                  & 
    566                   !!   &           + zfU(ji,jj) + zCorx(ji,jj) + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj)  & 
    567                   !!   &           ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchU(ji,jj) 
     592                  ENDIF 
    568593               END DO 
    569594            END DO 
     
    578603            DO jj = 2, jpjm1 
    579604               DO ji = fs_2, fs_jpim1 
    580           
    581                   ! tau_io/(v_oce - v_ice) 
     605                  !                 !--- tau_io/(v_oce - v_ice) 
    582606                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  & 
    583607                     &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) ) 
    584  
    585                   ! Ocean-to-Ice stress 
     608                  !                 !--- Ocean-to-Ice stress 
    586609                  ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) ) 
    587  
    588                   ! tau_bottom/v_ice 
     610                  ! 
     611                  !                 !--- tau_bottom/v_ice 
    589612                  zvel  = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) ) 
    590613                  ztauB = - tau_icebfr(ji,jj) / zvel 
    591                    
    592                   ! Coriolis at V-points (energy conserving formulation) 
     614                  ! 
     615                  !                 !--- Coriolis at v-points (energy conserving formulation) 
    593616                  zCory(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  & 
    594617                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  & 
    595618                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) ) 
    596  
    597                   ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
     619                  ! 
     620                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
    598621                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) 
    599  
    600                   ! landfast switch => 0 = static friction (tau_icebfr > zTauE); 1 = sliding friction 
     622                  ! 
     623                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction 
    601624                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
    602                    
    603                   ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 
     625                  ! 
     626                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 
     627                  v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )         & ! previous velocity 
     628                     &                                  + zTauE + zTauO * v_ice(ji,jj)                                            & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     629                     &                                  ) / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
     630                     &                 + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0 
     631                     &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                               & ! v_ice = v_oce if mass < zmmin 
     632                     &           ) * zmaskV(ji,jj) 
     633                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
    604634                  v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                             &  ! previous velocity 
    605635                     &                                     + zTauE + zTauO * v_ice(ji,jj)                             &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     
    608638                     &              ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin 
    609639                     &            ) * zmaskV(ji,jj) 
    610                   ! Bouillon 2013 
    611                   !!v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) )                  & 
    612                   !!   &           + zfV(ji,jj) + zCory(ji,jj) + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj)  & 
    613                   !!   &           ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj) 
     640                  ENDIF 
    614641               END DO 
    615642            END DO 
     
    820847   END SUBROUTINE ice_dyn_rhg_evp 
    821848 
     849 
    822850   SUBROUTINE rhg_evp_rst( cdrw, kt ) 
    823851      !!--------------------------------------------------------------------- 
  • branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/LIM_SRC_3/iceforcing.F90

    r8637 r8813  
    9999      !! 
    100100      !! ** Action  : It provides the following fields used in sea ice model: 
    101       !!                fr1_i0  , fr2_i0                         = 1sr & 2nd fraction of qsr penetration in ice  [%] 
    102101      !!                emp_oce , emp_ice                        = E-P over ocean and sea ice                    [Kg/m2/s] 
    103102      !!                sprecip                                  = solid precipitation                           [Kg/m2/s] 
     
    142141         ! 
    143142      CASE( jp_blk )                !--- bulk formulation 
    144                                 CALL blk_ice_flx( t_su, alb_ice )    !  
    145          IF( ln_mixcpl      )   CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su ) 
     143                                CALL blk_ice_flx( t_su, h_s, h_i, alb_ice )    !  
     144         IF( ln_mixcpl      )   CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su, phs=h_s, phi=h_i ) 
    146145         IF( nn_iceflx /= 2 )   CALL ice_flx_dist( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_iceflx ) 
    147146         ! 
    148147      CASE ( jp_purecpl )           !--- coupled formulation 
    149                                 CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su ) 
     148                                CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su, phs=h_s, phi=h_i ) 
    150149         IF( nn_iceflx == 2 )   CALL ice_flx_dist( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_iceflx ) 
    151150         ! 
     
    268267      INTEGER ::   ios, ioptio   ! Local integer output status for namelist read 
    269268      !! 
    270       NAMELIST/namforcing/ rn_cio, rn_blow_s, nn_iceflx 
     269      NAMELIST/namforcing/ rn_cio, rn_blow_s, nn_iceflx, nice_jules 
    271270      !!------------------------------------------------------------------- 
    272271      ! 
     
    285284         WRITE(numout,*) '~~~~~~~~~~~~~~~' 
    286285         WRITE(numout,*) '   Namelist namforcing:' 
    287          WRITE(numout,*) '      drag coefficient for oceanic stress              rn_cio    = ', rn_cio 
    288          WRITE(numout,*) '      coefficient for ice-lead partition of snowfall   rn_blow_s = ', rn_blow_s 
    289          WRITE(numout,*) '      Multicategory heat flux formulation              nn_iceflx = ', nn_iceflx 
     286         WRITE(numout,*) '      drag coefficient for oceanic stress              rn_cio     = ', rn_cio 
     287         WRITE(numout,*) '      coefficient for ice-lead partition of snowfall   rn_blow_s  = ', rn_blow_s 
     288         WRITE(numout,*) '      Multicategory heat flux formulation              nn_iceflx  = ', nn_iceflx 
     289         WRITE(numout,*) '      Jules coupling (0=no, 1=emulated, 2=active)      nice_jules = ', nice_jules 
    290290      ENDIF 
    291291      ! 
  • branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/LIM_SRC_3/iceistate.F90

    r8637 r8813  
    4848   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   si  ! structure of input fields (file informations, fields read) 
    4949   ! 
    50    ! ** namelist (namini) ** 
     50   !                             !! ** namelist (namini) ** 
    5151   LOGICAL  ::   ln_iceini        ! initialization or not 
    5252   LOGICAL  ::   ln_iceini_file   ! Ice initialization state from 2D netcdf file 
     
    9191      !!              where there is no ice (clem: I do not know why, is it mandatory?)  
    9292      !!-------------------------------------------------------------------- 
    93       INTEGER  ::   ji, jj, jk, jl             ! dummy loop indices 
     93      INTEGER  ::   ji, jj, jk, jl         ! dummy loop indices 
     94      INTEGER  ::   i_hemis, i_fill, jl0   ! local integers 
    9495      REAL(wp) ::   ztmelts, zdh 
    95       INTEGER  ::   i_hemis, i_fill, jl0 
    9696      REAL(wp) ::   zarg, zV, zconv, zdv, zfac 
    9797      INTEGER , DIMENSION(4)           ::   itest 
     
    100100      REAL(wp), DIMENSION(jpi,jpj)     ::   zht_i_ini, zat_i_ini, zvt_i_ini            !data from namelist or nc file 
    101101      REAL(wp), DIMENSION(jpi,jpj)     ::   zts_u_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file 
    102       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zh_i_ini, za_i_ini                         !data by cattegories to fill 
     102      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zh_i_ini , za_i_ini                        !data by cattegories to fill 
    103103      !-------------------------------------------------------------------- 
    104104 
     
    116116         tn_ice(:,:,jl) = rt0 * tmask(:,:,1) 
    117117      END DO 
    118  
    119       ! init basal temperature (considered at freezing point) 
     118      ! 
     119      ! init basal temperature (considered at freezing point)   [Kelvin] 
    120120      CALL eos_fzp( sss_m(:,:), t_bo(:,:) ) 
    121121      t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1)  
     
    142142            zvt_i_ini(:,:) = zht_i_ini(:,:) * zat_i_ini(:,:) 
    143143            ! 
    144          !                             !---------------! 
     144            !                          !---------------! 
    145145         ELSE                          ! Read namelist ! 
    146146            !                          !---------------! 
    147  
    148            ! no ice if sst <= t-freez + ttest 
    149             WHERE( ( sst_m(:,:) - (t_bo(:,:) - rt0) ) * tmask(:,:,1) >= rn_thres_sst ) ; zswitch(:,:) = 0._wp  
    150             ELSEWHERE                                                                  ; zswitch(:,:) = tmask(:,:,1) 
     147            ! no ice if sst <= t-freez + ttest 
     148            WHERE( ( sst_m(:,:) - (t_bo(:,:) - rt0) ) * tmask(:,:,1) >= rn_thres_sst )   ;   zswitch(:,:) = 0._wp  
     149            ELSEWHERE                                                                    ;   zswitch(:,:) = tmask(:,:,1) 
    151150            END WHERE 
    152  
     151            ! 
    153152            ! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array 
    154153            WHERE( ff_t(:,:) >= 0._wp ) 
     
    242241                        ! 
    243242                     ENDIF 
    244  
     243                     ! 
    245244                     ! Compatibility tests 
    246245                     zconv = ABS( zat_i_ini(ji,jj) - SUM( za_i_ini(ji,jj,1:jpl) ) )           ! Test 1: area conservation 
    247246                     IF ( zconv < epsi06 ) itest(1) = 1 
    248                       
     247                     ! 
    249248                     zconv = ABS(       zat_i_ini(ji,jj)       * zht_i_ini(ji,jj)   &         ! Test 2: volume conservation 
    250249                        &        - SUM( za_i_ini (ji,jj,1:jpl) * zh_i_ini (ji,jj,1:jpl) ) ) 
    251250                     IF ( zconv < epsi06 ) itest(2) = 1 
    252                       
     251                     ! 
    253252                     IF ( zh_i_ini(ji,jj,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1           ! Test 3: thickness of the last category is in-bounds ? 
    254                       
     253                     ! 
    255254                     itest(4) = 1 
    256255                     DO jl = 1, i_fill 
     
    260259                  END DO                                                        ! end iteration on categories 
    261260                  !                                                             !---------------------------- 
    262                   ! 
    263261                  IF( lwp .AND. SUM(itest) /= 4 ) THEN  
    264262                     WRITE(numout,*) 
     
    270268                     WRITE(numout,*) ' zht_i_ini : ', zht_i_ini(ji,jj) 
    271269                  ENDIF 
    272                 
     270                  ! 
    273271               ENDIF 
    274272               ! 
     
    279277         ! 4) Fill in sea ice arrays 
    280278         !--------------------------------------------------------------------- 
    281  
     279         ! 
    282280         ! Ice concentration, thickness and volume, ice salinity, ice age, surface temperature 
    283281         DO jl = 1, jpl ! loop over categories 
     
    289287                  o_i(ji,jj,jl)  = 0._wp                                                     ! age (0 day) 
    290288                  t_su(ji,jj,jl) = zswitch(ji,jj) * zts_u_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0 ! surf temp 
    291  
     289                  ! 
    292290                  IF( zht_i_ini(ji,jj) > 0._wp )THEN 
    293291                    h_s(ji,jj,jl)= h_i(ji,jj,jl) * ( zht_s_ini(ji,jj) / zht_i_ini(ji,jj) )  ! snow depth 
     
    295293                    h_s(ji,jj,jl)= 0._wp 
    296294                  ENDIF 
    297  
     295                  ! 
    298296                  ! This case below should not be used if (h_s/h_i) is ok in namelist 
    299297                  ! In case snow load is in excess that would lead to transformation from snow to ice 
     
    303301                  h_i(ji,jj,jl) = MIN( hi_max(jl), h_i(ji,jj,jl) + zdh ) 
    304302                  h_s(ji,jj,jl) = MAX( 0._wp, h_s(ji,jj,jl) - zdh * rhoic * r1_rhosn ) 
    305  
     303                  ! 
    306304                  ! ice volume, salt content, age content 
    307305                  v_i (ji,jj,jl) = h_i(ji,jj,jl) * a_i(ji,jj,jl)              ! ice volume 
     
    312310            END DO 
    313311         END DO 
    314  
    315          ! for constant salinity in time 
    316          IF( nn_icesal /= 2 )  THEN 
     312         ! 
     313         IF( nn_icesal /= 2 )  THEN         ! for constant salinity in time 
    317314            CALL ice_var_salprof 
    318315            sv_i = s_i * v_i 
    319316         ENDIF 
    320              
     317         !   
    321318         ! Snow temperature and heat content 
    322319         DO jk = 1, nlay_s 
     
    327324                     ! Snow energy of melting 
    328325                     e_s(ji,jj,jk,jl) = zswitch(ji,jj) * rhosn * ( cpic * ( rt0 - t_s(ji,jj,jk,jl) ) + lfus ) 
    329  
     326                     ! 
    330327                     ! Mutliply by volume, and divide by number of layers to get heat content in J/m2 
    331328                     e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) * r1_nlay_s 
     
    334331            END DO 
    335332         END DO 
    336  
     333         ! 
    337334         ! Ice salinity, temperature and heat content 
    338335         DO jk = 1, nlay_i 
     
    343340                     sz_i(ji,jj,jk,jl) = zswitch(ji,jj) * zsm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rn_simin 
    344341                     ztmelts          = - tmut * sz_i(ji,jj,jk,jl) + rt0 !Melting temperature in K 
    345  
     342                     ! 
    346343                     ! heat content per unit volume 
    347                      e_i(ji,jj,jk,jl) = zswitch(ji,jj) * rhoic * (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) & 
    348                         +   lfus    * ( 1._wp - (ztmelts-rt0) / MIN((t_i(ji,jj,jk,jl)-rt0),-epsi20) ) & 
    349                         -   rcp     * ( ztmelts - rt0 ) ) 
    350  
     344                     e_i(ji,jj,jk,jl) = zswitch(ji,jj) * rhoic * (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) )         & 
     345                        &             + lfus * ( 1._wp - (ztmelts-rt0) / MIN( (t_i(ji,jj,jk,jl)-rt0) , -epsi20 )  )  & 
     346                        &             - rcp  * ( ztmelts - rt0 ) ) 
     347                     ! 
    351348                     ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2 
    352349                     e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i 
     
    355352            END DO 
    356353         END DO 
    357  
     354         ! 
    358355         tn_ice (:,:,:) = t_su (:,:,:) 
    359  
     356         ! 
    360357         ! Melt pond volume and fraction 
    361358         IF ( ln_pnd_CST .OR. ln_pnd_H12 ) THEN   ;   zfac = 1._wp 
     
    368365         a_ip(:,:,:) = a_ip_frac(:,:,:) * a_i (:,:,:)  
    369366         v_ip(:,:,:) = h_ip     (:,:,:) * a_ip(:,:,:) 
    370  
     367         ! 
    371368      ELSE ! if ln_iceini=false 
    372369         a_i  (:,:,:) = 0._wp 
     
    379376         s_i  (:,:,:) = 0._wp 
    380377         o_i  (:,:,:) = 0._wp 
    381  
     378         ! 
    382379         e_i(:,:,:,:) = 0._wp 
    383380         e_s(:,:,:,:) = 0._wp 
    384  
     381         ! 
    385382         DO jl = 1, jpl 
    386383            DO jk = 1, nlay_i 
     
    391388            END DO 
    392389         END DO 
    393  
     390         ! 
    394391         a_ip(:,:,:)      = 0._wp 
    395392         v_ip(:,:,:)      = 0._wp 
    396393         a_ip_frac(:,:,:) = 0._wp 
    397394         h_ip     (:,:,:) = 0._wp 
    398  
     395         ! 
    399396      ENDIF ! ln_iceini 
    400        
     397      ! 
    401398      at_i (:,:) = 0.0_wp 
    402399      DO jl = 1, jpl 
     
    405402      ! 
    406403      ! --- set ice velocities --- ! 
    407       u_ice (:,:)     = 0._wp 
    408       v_ice (:,:)     = 0._wp 
     404      u_ice (:,:) = 0._wp 
     405      v_ice (:,:) = 0._wp 
    409406      ! 
    410407      !---------------------------------------------- 
     
    415412      ! 
    416413      IF( ln_ice_embd ) THEN            ! embedded sea-ice: deplete the initial ssh below sea-ice area 
    417  
     414         ! 
    418415         sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0 
    419416         sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0 
    420  
     417         ! 
    421418         IF( .NOT.ln_linssh ) THEN 
    422              
     419            ! 
    423420            WHERE( ht_0(:,:) > 0 )   ;   z2d(:,:) = 1._wp + sshn(:,:)*tmask(:,:,1) / ht_0(:,:) 
    424421            ELSEWHERE                ;   z2d(:,:) = 1._wp   ;   END WHERE 
    425           
     422            ! 
    426423            DO jk = 1,jpkm1                     ! adjust initial vertical scale factors                 
    427424               e3t_n(:,:,jk) = e3t_0(:,:,jk) * z2d(:,:) 
     
    429426               e3t_a(:,:,jk) = e3t_n(:,:,jk) 
    430427            END DO 
    431              
     428            ! 
    432429            ! Reconstruction of all vertical scale factors at now and before time-steps 
    433430            ! ========================================================================= 
     
    478475!!      CALL dia_wri_state( 'output.init', nit000 ) 
    479476!!!       
    480  
     477      ! 
    481478   END SUBROUTINE ice_istate 
     479 
    482480 
    483481   SUBROUTINE ice_istate_init 
     
    485483      !!                   ***  ROUTINE ice_istate_init  *** 
    486484      !!         
    487       !! ** Purpose : Definition of initial state of the ice  
    488       !! 
    489       !! ** Method : Read the namini namelist and check the parameter  
    490       !!       values called at the first timestep (nit000) 
    491       !! 
    492       !! ** input :  
    493       !!        Namelist namini 
    494       !! 
    495       !! history : 
    496       !!  8.5  ! 03-08 (C. Ethe) original code  
    497       !!  8.5  ! 07-11 (M. Vancoppenolle) rewritten initialization 
     485      !! ** Purpose :   Definition of initial state of the ice  
     486      !! 
     487      !! ** Method  :   Read the namini namelist and check the parameter  
     488      !!              values called at the first timestep (nit000) 
     489      !! 
     490      !! ** input   :  Namelist namini 
     491      !! 
    498492      !!----------------------------------------------------------------------------- 
    499       ! 
    500       INTEGER :: ios,ierr,inum_ice                 ! Local integer output status for namelist read 
    501       INTEGER :: ji,jj 
    502       INTEGER :: ifpr, ierror 
     493      INTEGER ::   ji, jj 
     494      INTEGER ::   ios, ierr, inum_ice   ! Local integer output status for namelist read 
     495      INTEGER ::   ifpr, ierror 
    503496      ! 
    504497      CHARACTER(len=256) ::  cn_dir          ! Root directory for location of ice files 
     
    515508      READ  ( numnam_ice_ref, namini, IOSTAT = ios, ERR = 901) 
    516509901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namini in reference namelist', lwp ) 
    517  
     510      ! 
    518511      REWIND( numnam_ice_cfg )              ! Namelist namini in configuration namelist : Ice initial state 
    519512      READ  ( numnam_ice_cfg, namini, IOSTAT = ios, ERR = 902 ) 
    520513902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namini in configuration namelist', lwp ) 
    521514      IF(lwm) WRITE ( numoni, namini ) 
    522  
     515      ! 
    523516      slf_i(jp_hti) = sn_hti  ;  slf_i(jp_hts) = sn_hts 
    524517      slf_i(jp_ati) = sn_ati  ;  slf_i(jp_tsu) = sn_tsu 
     
    545538         WRITE(numout,*) '      initial  ice/snw temp  in the south                    rn_tmi_ini_s    = ', rn_tmi_ini_s 
    546539      ENDIF 
    547  
     540      ! 
    548541      IF( ln_iceini_file ) THEN                      ! Ice initialization using input file 
    549542         ! 
     
    553546            CALL ctl_stop( 'Ice_ini in iceistate: unable to allocate si structure' )   ;   RETURN 
    554547         ENDIF 
    555  
     548         ! 
    556549         DO ifpr = 1, jpfldi 
    557550            ALLOCATE( si(ifpr)%fnow(jpi,jpj,1) ) 
    558551            ALLOCATE( si(ifpr)%fdta(jpi,jpj,1,2) ) 
    559552         END DO 
    560  
     553         ! 
    561554         ! fill si with slf_i and control print 
    562555         CALL fld_fill( si, slf_i, cn_dir, 'ice_istate', 'ice istate ini', 'numnam_ice' ) 
    563  
     556         ! 
    564557         CALL fld_read( nit000, 1, si )                ! input fields provided at the current time-step 
    565  
     558         ! 
    566559      ENDIF 
    567  
     560      ! 
    568561   END SUBROUTINE ice_istate_init 
    569562 
  • branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/LIM_SRC_3/iceitd.F90

    r8637 r8813  
    1313   !!   'key_lim3'                                       ESIM sea-ice model 
    1414   !!---------------------------------------------------------------------- 
    15    !!   ice_itd_init   : read ice thicknesses mean and min from namelist 
    16    !!   ice_itd_rem    : redistribute ice thicknesses after thermo growth and melt 
    17    !!   ice_itd_reb    : rebin ice thicknesses into bounded categories 
     15   !!   ice_itd_rem   : redistribute ice thicknesses after thermo growth and melt 
     16   !!   itd_glinear   : build g(h) satisfying area and volume constraints 
     17   !!   itd_shiftice  : shift ice across category boundaries, conserving everything 
     18   !!   ice_itd_reb   : rebin ice thicknesses into bounded categories 
     19   !!   ice_itd_init  : read ice thicknesses mean and min from namelist 
    1820   !!---------------------------------------------------------------------- 
    1921   USE dom_oce        ! ocean domain 
     
    3638   PUBLIC   ice_itd_reb   ! called in icecor 
    3739 
    38    ! ** namelist (namitd) ** 
    39    REAL(wp) ::   rn_himean   ! mean thickness of the domain 
    40  
     40   INTEGER            ::   nice_catbnd     ! choice of the type of ice category function 
     41   !                                       ! associated indices: 
     42   INTEGER, PARAMETER ::   np_cathfn = 1   ! categories defined by a function 
     43   INTEGER, PARAMETER ::   np_catusr = 2   ! categories defined by the user 
     44   ! 
     45   !                                           !! ** namelist (namitd) ** 
     46   LOGICAL                    ::   ln_cat_hfn   ! ice categories are defined by function like rn_himean**(-0.05) 
     47   REAL(wp)                   ::   rn_himean    ! mean thickness of the domain 
     48   LOGICAL                    ::   ln_cat_usr   ! ice categories are defined by rn_catbnd 
     49   REAL(wp), DIMENSION(0:100) ::   rn_catbnd    ! ice categories bounds 
     50   ! 
    4151   !!---------------------------------------------------------------------- 
    4252   !! NEMO/ICE 4.0 , NEMO Consortium (2017) 
     
    6575      REAL(wp) ::   zx3         
    6676      REAL(wp) ::   zslope          ! used to compute local thermodynamic "speeds" 
    67  
     77      ! 
    6878      INTEGER , DIMENSION(jpij)       ::   iptidx          ! compute remapping or not 
    6979      INTEGER , DIMENSION(jpij,jpl-1) ::   jdonor          ! donor category index 
     
    97107      !----------------------------------------------------------------------------------------------- 
    98108      IF( npti > 0 ) THEN 
    99  
     109         ! 
    100110         zdhice(:,:) = 0._wp 
    101111         zhbnew(:,:) = 0._wp 
    102  
     112         ! 
    103113         CALL tab_3d_2d( npti, nptidx(1:npti), h_i_2d (1:npti,1:jpl), h_i   ) 
    104114         CALL tab_3d_2d( npti, nptidx(1:npti), h_ib_2d(1:npti,1:jpl), h_i_b ) 
    105115         CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d  (1:npti,1:jpl), a_i    ) 
    106116         CALL tab_3d_2d( npti, nptidx(1:npti), a_ib_2d (1:npti,1:jpl), a_i_b  ) 
    107  
     117         ! 
    108118         DO jl = 1, jpl 
    109119            ! Compute thickness change in each ice category 
     
    112122            END DO 
    113123         END DO 
    114           
     124         ! 
    115125         ! --- New boundaries for category 1:jpl-1 --- ! 
    116126         DO jl = 1, jpl - 1 
     
    135145               ! 1) hn(t+1)+espi < Hn* < hn+1(t+1)-epsi                
    136146               !    Note: hn(t+1) must not be too close to either HR or HL otherwise a division by nearly 0 is possible  
    137                !          in ice_itd_glinear in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice) 
     147               !          in itd_glinear in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice) 
    138148               IF( a_i_2d(ji,jl  ) > epsi10 .AND. h_i_2d(ji,jl  ) > ( zhbnew(ji,jl) - epsi10 ) )   nptidx(ji) = 0 
    139149               IF( a_i_2d(ji,jl+1) > epsi10 .AND. h_i_2d(ji,jl+1) < ( zhbnew(ji,jl) + epsi10 ) )   nptidx(ji) = 0 
    140                 
     150               ! 
    141151               ! 2) Hn-1 < Hn* < Hn+1   
    142152               IF( zhbnew(ji,jl) < hi_max(jl-1) )   nptidx(ji) = 0 
    143153               IF( zhbnew(ji,jl) > hi_max(jl+1) )   nptidx(ji) = 0 
    144                 
     154               ! 
    145155            END DO 
    146156         END DO 
     
    153163               zhbnew(ji,jpl) = hi_max(jpl)   
    154164            ENDIF 
    155              
     165            ! 
    156166            ! --- 1 additional condition for remapping (1st category) --- ! 
    157167            ! H0+epsi < h1(t) < H1-epsi  
    158168            !    h1(t) must not be too close to either HR or HL otherwise a division by nearly 0 is possible  
    159             !    in ice_itd_glinear in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice) 
     169            !    in itd_glinear in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice) 
    160170            IF( h_ib_2d(ji,1) < ( hi_max(0) + epsi10 ) )   nptidx(ji) = 0 
    161171            IF( h_ib_2d(ji,1) > ( hi_max(1) - epsi10 ) )   nptidx(ji) = 0 
     
    165175         !  3) Identify cells where remapping 
    166176         !----------------------------------------------------------------------------------------------- 
    167          ipti = 0 ; iptidx(:) = 0 
     177         ipti = 0   ;  iptidx(:) = 0 
    168178         DO ji = 1, npti 
    169179            IF( nptidx(ji) /= 0 ) THEN 
     
    197207               !   
    198208               ! --- g(h) for category 1 --- ! 
    199                CALL ice_itd_glinear( zhb0(1:npti)  , zhb1(1:npti)  , h_ib_1d(1:npti)  , a_i_1d(1:npti)  ,  &   ! in 
    200                   &                  g0  (1:npti,1), g1  (1:npti,1), hL      (1:npti,1), hR    (1:npti,1)   )   ! out 
     209               CALL itd_glinear( zhb0(1:npti)  , zhb1(1:npti)  , h_ib_1d(1:npti)  , a_i_1d(1:npti)  ,  &   ! in 
     210                  &              g0  (1:npti,1), g1  (1:npti,1), hL      (1:npti,1), hR    (1:npti,1)   )   ! out 
    201211                  ! 
    202212               ! Area lost due to melting of thin ice 
     
    240250            ! 
    241251            ! --- g(h) for each thickness category --- !   
    242             CALL ice_itd_glinear( zhbnew(1:npti,jl-1), zhbnew(1:npti,jl), h_i_1d(1:npti)   , a_i_1d(1:npti)   ,  &   ! in 
    243                &                  g0    (1:npti,jl  ), g1    (1:npti,jl), hL     (1:npti,jl), hR    (1:npti,jl)   )   ! out 
     252            CALL itd_glinear( zhbnew(1:npti,jl-1), zhbnew(1:npti,jl), h_i_1d(1:npti)   , a_i_1d(1:npti)   ,  &   ! in 
     253               &              g0    (1:npti,jl  ), g1    (1:npti,jl), hL     (1:npti,jl), hR    (1:npti,jl)   )   ! out 
    244254            ! 
    245255         END DO 
     
    281291         ! 6) Shift ice between categories 
    282292         !---------------------------------------------------------------------------------------------- 
    283          CALL ice_itd_shiftice ( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) ) 
     293         CALL itd_shiftice ( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) ) 
    284294          
    285295         !---------------------------------------------------------------------------------------------- 
     
    289299         CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d  (1:npti), a_i(:,:,1)    ) 
    290300         CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d (1:npti), a_ip(:,:,1)   ) 
    291           
     301         ! 
    292302         DO ji = 1, npti 
    293303            IF ( a_i_1d(ji) > epsi10 .AND. h_i_1d(ji) < rn_himin ) THEN 
     
    309319 
    310320 
    311    SUBROUTINE ice_itd_glinear( HbL, Hbr, phice, paice, pg0, pg1, phL, phR ) 
    312       !!------------------------------------------------------------------ 
    313       !!                ***  ROUTINE ice_itd_glinear *** 
     321   SUBROUTINE itd_glinear( HbL, Hbr, phice, paice, pg0, pg1, phL, phR ) 
     322      !!------------------------------------------------------------------ 
     323      !!                ***  ROUTINE itd_glinear *** 
    314324      !! 
    315325      !! ** Purpose :   build g(h) satisfying area and volume constraints (Eq. 6 and 9) 
     
    367377      END DO 
    368378      ! 
    369    END SUBROUTINE ice_itd_glinear 
    370  
    371  
    372    SUBROUTINE ice_itd_shiftice( kdonor, pdaice, pdvice ) 
    373       !!------------------------------------------------------------------ 
    374       !!                ***  ROUTINE ice_itd_shiftice *** 
     379   END SUBROUTINE itd_glinear 
     380 
     381 
     382   SUBROUTINE itd_shiftice( kdonor, pdaice, pdvice ) 
     383      !!------------------------------------------------------------------ 
     384      !!                ***  ROUTINE itd_shiftice *** 
    375385      !! 
    376386      !! ** Purpose :   shift ice across category boundaries, conserving everything 
     
    486496            END DO 
    487497         END DO 
    488  
     498         ! 
    489499         DO jk = 1, nlay_i         !--- Ice heat content 
    490500            DO ji = 1, npti 
     
    529539      CALL tab_2d_3d( npti, nptidx(1:npti), t_su_2d(1:npti,1:jpl), t_su ) 
    530540      ! 
    531    END SUBROUTINE ice_itd_shiftice 
     541   END SUBROUTINE itd_shiftice 
    532542    
    533543 
     
    550560      ! 
    551561      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_itd_reb: rebining ice thickness distribution'  
    552  
     562      ! 
    553563      jdonor(:,:) = 0 
    554564      zdaice(:,:) = 0._wp 
     
    587597         ! 
    588598         IF( npti > 0 ) THEN 
    589             CALL ice_itd_shiftice( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) )  ! Shift jl=>jl+1 
     599            CALL itd_shiftice( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) )  ! Shift jl=>jl+1 
    590600            ! Reset shift parameters 
    591601            jdonor(1:npti,jl) = 0 
     
    618628         ! 
    619629         IF( npti > 0 ) THEN 
    620             CALL ice_itd_shiftice( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) )  ! Shift jl+1=>jl 
     630            CALL itd_shiftice( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) )  ! Shift jl+1=>jl 
    621631            ! Reset shift parameters 
    622632            jdonor(1:npti,jl) = 0 
     
    629639   END SUBROUTINE ice_itd_reb 
    630640 
     641 
    631642   SUBROUTINE ice_itd_init 
    632643      !!------------------------------------------------------------------ 
     
    637648      !! ** input   :   Namelist namitd 
    638649      !!------------------------------------------------------------------- 
    639       INTEGER  ::   jl    ! dummy loop index 
    640       INTEGER  ::   ios   ! Local integer output status for namelist read 
     650      INTEGER  ::   jl            ! dummy loop index 
     651      INTEGER  ::   ios, ioptio   ! Local integer output status for namelist read 
    641652      REAL(wp) ::   zhmax, znum, zden, zalpha   !   -      - 
    642       !! 
    643       NAMELIST/namitd/ rn_himean, rn_himin 
     653      ! 
     654      NAMELIST/namitd/ ln_cat_hfn, rn_himean, ln_cat_usr, rn_catbnd, rn_himin 
    644655      !!------------------------------------------------------------------ 
    645656      ! 
     
    647658      READ  ( numnam_ice_ref, namitd, IOSTAT = ios, ERR = 901) 
    648659901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namitd in reference namelist', lwp ) 
    649  
     660      ! 
    650661      REWIND( numnam_ice_cfg )      ! Namelist namitd in configuration namelist : Parameters for ice 
    651662      READ  ( numnam_ice_cfg, namitd, IOSTAT = ios, ERR = 902 ) 
     
    658669         WRITE(numout,*) '~~~~~~~~~~~~' 
    659670         WRITE(numout,*) '   Namelist namitd: ' 
    660          WRITE(numout,*) '      mean ice thickness in the domain               rn_himean = ', rn_himean 
    661          WRITE(numout,*) '      minimum ice thickness                          rn_himin  = ', rn_himin  
     671         WRITE(numout,*) '      Ice categories are defined by a function of rn_himean**(-0.05)    ln_cat_hfn = ', ln_cat_hfn 
     672         WRITE(numout,*) '         mean ice thickness in the domain                               rn_himean  = ', rn_himean 
     673         WRITE(numout,*) '      Ice categories are defined by rn_catbnd                           ln_cat_usr = ', ln_cat_usr 
     674         WRITE(numout,*) '         ice categories boundaries (m)                                  rn_catbnd  = ', rn_catbnd  
     675         WRITE(numout,*) '      minimum ice thickness                                             rn_himin   = ', rn_himin  
    662676      ENDIF 
    663677      ! 
     
    665679      !  Thickness categories boundaries  ! 
    666680      !-----------------------------------! 
    667       ! 
    668       zalpha = 0.05_wp              ! max of each category (from h^(-alpha) function) 
    669       zhmax  = 3._wp * rn_himean 
    670       DO jl = 1, jpl 
    671          znum = jpl * ( zhmax+1 )**zalpha 
    672          zden = REAL( jpl-jl , wp ) * ( zhmax + 1._wp )**zalpha + REAL( jl , wp ) 
    673          hi_max(jl) = ( znum / zden )**(1./zalpha) - 1 
    674       END DO 
     681      !                             !== set the choice of ice categories ==! 
     682      ioptio = 0  
     683      IF( ln_cat_hfn ) THEN   ;   ioptio = ioptio + 1   ;   nice_catbnd = np_cathfn    ;   ENDIF 
     684      IF( ln_cat_usr ) THEN   ;   ioptio = ioptio + 1   ;   nice_catbnd = np_catusr    ;   ENDIF 
     685      IF( ioptio /= 1 )   CALL ctl_stop( 'ice_itd_init: choose one and only one ice categories boundaries' ) 
     686      ! 
     687      SELECT CASE( nice_catbnd ) 
     688      !                                !------------------------! 
     689      CASE( np_cathfn )                ! h^(-alpha) function 
     690         !                             !------------------------! 
     691         zalpha = 0.05_wp 
     692         zhmax  = 3._wp * rn_himean 
     693         DO jl = 1, jpl 
     694            znum = jpl * ( zhmax+1 )**zalpha 
     695            zden = REAL( jpl-jl , wp ) * ( zhmax + 1._wp )**zalpha + REAL( jl , wp ) 
     696            hi_max(jl) = ( znum / zden )**(1./zalpha) - 1 
     697         END DO 
     698         !                             !------------------------! 
     699      CASE( np_catusr )                ! user defined 
     700         !                             !------------------------! 
     701         DO jl = 0, jpl 
     702            hi_max(jl) = rn_catbnd(jl) 
     703         END DO 
     704         ! 
     705      END SELECT 
    675706      ! 
    676707      DO jl = 1, jpl                ! mean thickness by category 
  • branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90

    r8637 r8813  
    153153         !------------------------------------------------------! 
    154154         ! It provides the following fields used in sea ice model: 
    155          !    fr1_i0  , fr2_i0     = 1sr & 2nd fraction of qsr penetration in ice  [%] 
    156155         !    emp_oce , emp_ice    = E-P over ocean and sea ice                    [Kg/m2/s] 
    157156         !    sprecip              = solid precipitation                           [Kg/m2/s] 
     
    320319         IF(lwp) WRITE(numout,*) '   nn_monocat forced to 0 as jpl>1, i.e. multi-category case is chosen' 
    321320      ENDIF 
    322       IF ( jpl == 1 .AND. nn_monocat == 0 ) THEN 
    323          CALL ctl_stop( 'STOP', 'par_init : if jpl=1 then nn_monocat should be between 1 and 4' ) 
    324       ENDIF 
     321!     IF ( jpl == 1 .AND. nn_monocat == 0 ) THEN 
     322!        CALL ctl_stop( 'STOP', 'par_init : if jpl=1 then nn_monocat should be between 1 and 4' ) 
     323!     ENDIF 
    325324      ! 
    326325      IF( ln_bdy .AND. ln_icediachk )   CALL ctl_warn('par_init: online conservation check does not work with BDY') 
  • branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90

    r8637 r8813  
    2626   USE sbc_oce , ONLY : sss_m, sst_m, e3t_m, utau, vtau, ssu_m, ssv_m, frq_m, qns_tot, qsr_tot, sprecip, ln_cpl 
    2727   USE sbc_ice , ONLY : qsr_oce, qns_oce, qemp_oce, qsr_ice, qns_ice, dqns_ice, evap_ice, qprec_ice, qevap_ice, & 
    28       &                 fr1_i0, fr2_i0 
     28      &                 qml_ice, qcn_ice, qsr_ice_tr 
    2929   USE ice1D          ! sea-ice: thermodynamics variables 
    3030   USE icethd_zdf     ! sea-ice: vertical heat diffusion 
     
    383383         CALL tab_2d_1d( npti, nptidx(1:npti), qprec_ice_1d(1:npti), qprec_ice        ) 
    384384         CALL tab_2d_1d( npti, nptidx(1:npti), qsr_ice_1d  (1:npti), qsr_ice (:,:,kl) ) 
    385          CALL tab_2d_1d( npti, nptidx(1:npti), fr1_i0_1d   (1:npti), fr1_i0           ) 
    386          CALL tab_2d_1d( npti, nptidx(1:npti), fr2_i0_1d   (1:npti), fr2_i0           ) 
    387385         CALL tab_2d_1d( npti, nptidx(1:npti), qns_ice_1d  (1:npti), qns_ice (:,:,kl) ) 
    388386         CALL tab_2d_1d( npti, nptidx(1:npti), ftr_ice_1d  (1:npti), ftr_ice (:,:,kl) ) 
     
    393391         CALL tab_2d_1d( npti, nptidx(1:npti), fhtur_1d    (1:npti), fhtur            ) 
    394392         CALL tab_2d_1d( npti, nptidx(1:npti), fhld_1d     (1:npti), fhld             ) 
     393          
     394         CALL tab_2d_1d( npti, nptidx(1:npti), qml_ice_1d   (1:npti), qml_ice      (:,:,kl)  ) 
     395         CALL tab_2d_1d( npti, nptidx(1:npti), qcn_ice_1d   (1:npti), qcn_ice      (:,:,kl) ) 
     396         CALL tab_2d_1d( npti, nptidx(1:npti), qsr_ice_tr_1d(1:npti), qsr_ice_tr   (:,:,kl) ) 
    395397         ! 
    396398         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_sni_1d(1:npti), wfx_snw_sni   ) 
     
    498500         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_spr_1d (1:npti), wfx_spr        ) 
    499501         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_lam_1d (1:npti), wfx_lam        ) 
    500          CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd        ) 
     502         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd        ) 
    501503         ! 
    502504         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_bog_1d (1:npti), sfx_bog        ) 
  • branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/LIM_SRC_3/icethd_dh.F90

    r8637 r8813  
    130130      !------------------------------------------------------------------------------! 
    131131      ! 
    132       DO ji = 1, npti 
    133          zdum       = qns_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji)  
     132      SELECT CASE ( nice_jules ) 
     133      ! 
     134      CASE ( np_jules_ACTIVE ) 
     135         ! 
     136         DO ji = 1, npti 
     137            zq_su(ji)        =   MAX( 0._wp, qml_ice_1d(ji) * rdt_ice ) 
     138         END DO 
     139         ! 
     140      CASE ( np_jules_EMULE ) 
     141         ! 
     142         DO ji = 1, npti 
     143            zdum             =   ( qns_ice_1d(ji) + qsr_ice_1d(ji) - qsr_ice_tr_1d(ji) - fc_su(ji) )  
     144            qml_ice_1d(ji)   =   zdum * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) ) 
     145            zq_su(ji)        =   MAX( 0._wp, qml_ice_1d(ji) * rdt_ice ) 
     146         END DO 
     147         ! 
     148      CASE ( np_jules_OFF )  
     149         ! 
     150         DO ji = 1, npti 
     151            zdum             =   ( qns_ice_1d(ji) + qsr_ice_1d(ji) - qsr_ice_tr_1d(ji) - fc_su(ji) )  
     152            zdum             =   zdum * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) ) 
     153            zq_su(ji)        =   MAX( 0._wp, zdum             * rdt_ice ) 
     154         END DO 
     155         ! 
     156      END SELECT 
     157      ! 
     158      DO ji = 1, npti 
    134159         zf_tt(ji)  = fc_bo_i(ji) + fhtur_1d(ji) + fhld_1d(ji)  
    135  
    136          zq_su (ji) = MAX( 0._wp, zdum      * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) ) 
    137160         zq_bo (ji) = MAX( 0._wp, zf_tt(ji) * rdt_ice ) 
    138161      END DO 
  • branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/LIM_SRC_3/icethd_zdf.F90

    r8637 r8813  
    1515   !!   'key_lim3'                                       ESIM sea-ice model 
    1616   !!---------------------------------------------------------------------- 
     17   !!  ice_thd_zdf      : select the appropriate routine for vertical heat diffusion calculation 
     18   !!  ice_thd_zdf_BL99 :  
     19   !!  ice_thd_enmelt   : 
     20   !!  ice_thd_zdf_init : 
     21   !!---------------------------------------------------------------------- 
    1722   USE dom_oce        ! ocean space and time domain 
    1823   USE phycst         ! physical constants (ocean directory)  
     
    3439   LOGICAL  ::   ln_cndi_U64      ! thermal conductivity: Untersteiner (1964) 
    3540   LOGICAL  ::   ln_cndi_P07      ! thermal conductivity: Pringle et al (2007) 
    36    REAL(wp) ::   rn_cnd_s         ! thermal conductivity of the snow [W/m/K] 
    3741   REAL(wp) ::   rn_kappa_i       ! coef. for the extinction of radiation Grenfell et al. (2006) [1/m] 
    38    LOGICAL  ::   ln_dqns_i        ! change non-solar surface flux with changing surface temperature (T) or not (F) 
    39  
     42   REAL(wp), PUBLIC ::   rn_cnd_s ! thermal conductivity of the snow [W/m/K] 
     43 
     44   INTEGER  ::   nice_zdf         ! Choice of the type of vertical heat diffusion formulation 
     45   !                                           ! associated indices: 
     46   INTEGER, PARAMETER ::   np_BL99    = 1      ! Bitz and Lipscomb (1999) 
     47 
     48   INTEGER, PARAMETER ::   np_zdf_jules_OFF    = 0   !   compute all temperatures from qsr and qns 
     49   INTEGER, PARAMETER ::   np_zdf_jules_SND    = 1   !   compute conductive heat flux and surface temperature from qsr and qns 
     50   INTEGER, PARAMETER ::   np_zdf_jules_RCV    = 2   !   compute snow and inner ice temperatures from qcnd 
     51    
    4052   !!---------------------------------------------------------------------- 
    4153   !! NEMO/ICE 4.0 , NEMO Consortium (2017) 
     
    4860      !!------------------------------------------------------------------- 
    4961      !!                ***  ROUTINE ice_thd_zdf  *** 
    50       !! ** Purpose : 
    51       !!           This routine determines the time evolution of snow and sea-ice  
    52       !!           temperature profiles. 
    53       !! ** Method  : 
    54       !!           This is done by solving the heat equation diffusion with 
    55       !!           a Neumann boundary condition at the surface and a Dirichlet one 
    56       !!           at the bottom. Solar radiation is partially absorbed into the ice. 
    57       !!           The specific heat and thermal conductivities depend on ice salinity 
    58       !!           and temperature to take into account brine pocket melting. The  
    59       !!           numerical 
    60       !!           scheme is an iterative Crank-Nicolson on a non-uniform multilayer grid  
    61       !!           in the ice and snow system. 
     62      !! 
     63      !! ** Purpose :   select the appropriate routine for the computation 
     64      !!              of vertical diffusion 
     65      !!------------------------------------------------------------------- 
     66    
     67      SELECT CASE ( nice_zdf )      ! Choose the vertical heat diffusion solver 
     68      ! 
     69      !                             !-------------       
     70      CASE( np_BL99 )               ! BL99 solver 
     71         !                          !------------- 
     72         SELECT CASE( nice_jules ) 
     73         CASE( np_jules_OFF    )   ;   CALL ice_thd_zdf_BL99 ( np_zdf_jules_OFF )   ! No Jules coupler      
     74         CASE( np_jules_EMULE  )   ;   CALL ice_thd_zdf_BL99 ( np_zdf_jules_SND )   ! Jules coupler is emulated (send/recieve) 
     75                                       CALL ice_thd_zdf_BL99 ( np_zdf_jules_RCV ) 
     76         CASE( np_jules_ACTIVE )   ;   CALL ice_thd_zdf_BL99 ( np_zdf_jules_RCV )   ! Jules coupler is active (receive only) 
     77         END SELECT 
     78      END SELECT 
     79       
     80   END SUBROUTINE ice_thd_zdf 
     81    
     82    
     83   SUBROUTINE ice_thd_zdf_BL99( k_jules ) 
     84      !!------------------------------------------------------------------- 
     85      !!                ***  ROUTINE ice_thd_zdf_BL99  *** 
     86      !! 
     87      !! ** Purpose :   computes the time evolution of snow and sea-ice temperature 
     88      !!              profiles, using the original Bitz and Lipscomb (1999) algorithm 
     89      !! 
     90      !! ** Method  :   solves the heat equation diffusion with a Neumann boundary 
     91      !!              condition at the surface and a Dirichlet one at the bottom.  
     92      !!              Solar radiation is partially absorbed into the ice. 
     93      !!              The specific heat and thermal conductivities depend on ice  
     94      !!              salinity and temperature to take into account brine pocket    
     95      !!              melting. The numerical scheme is an iterative Crank-Nicolson 
     96      !!              on a non-uniform multilayer grid in the ice and snow system. 
    6297      !! 
    6398      !!           The successive steps of this routine are 
     
    83118      !!           total ice/snow thickness : h_i_1d, h_s_1d 
    84119      !!------------------------------------------------------------------- 
     120      INTEGER, INTENT(in) ::   k_jules     ! Jules coupling (0=OFF, 1=RECEIVE, 2=SEND) 
     121      ! 
    85122      INTEGER ::   ji, jk         ! spatial loop index 
    86123      INTEGER ::   jm             ! current reference number of equation 
     
    88125      INTEGER ::   iconv          ! number of iterations in iterative procedure 
    89126      INTEGER ::   iconv_max = 50 ! max number of iterations in iterative procedure 
    90        
     127      ! 
    91128      INTEGER, DIMENSION(jpij) ::   jm_min    ! reference number of top equation 
    92129      INTEGER, DIMENSION(jpij) ::   jm_max    ! reference number of bottom equation 
    93        
     130      ! 
    94131      REAL(wp) ::   zg1s      =  2._wp        ! for the tridiagonal system 
    95132      REAL(wp) ::   zg1       =  2._wp        ! 
     
    101138      REAL(wp) ::   zdti_bnd  =  1.e-4_wp     ! maximal authorized error on temperature  
    102139      REAL(wp) ::   ztmelt_i                  ! ice melting temperature 
    103       REAL(wp) ::   z1_hsu 
    104140      REAL(wp) ::   zdti_max                  ! current maximal error on temperature  
    105141      REAL(wp) ::   zcpi                      ! Ice specific heat 
    106142      REAL(wp) ::   zhfx_err, zdq             ! diag errors on heat 
    107143      REAL(wp) ::   zfac                      ! dummy factor 
    108        
     144      ! 
    109145      REAL(wp), DIMENSION(jpij)     ::   isnow       ! switch for presence (1) or absence (0) of snow 
    110146      REAL(wp), DIMENSION(jpij)     ::   ztsub       ! surface temperature at previous iteration 
    111147      REAL(wp), DIMENSION(jpij)     ::   zh_i, z1_h_i ! ice layer thickness 
    112148      REAL(wp), DIMENSION(jpij)     ::   zh_s, z1_h_s ! snow layer thickness 
    113       REAL(wp), DIMENSION(jpij)     ::   zfsw        ! solar radiation absorbed at the surface 
    114149      REAL(wp), DIMENSION(jpij)     ::   zqns_ice_b  ! solar radiation absorbed at the surface 
    115150      REAL(wp), DIMENSION(jpij)     ::   zfnet       ! surface flux function 
    116151      REAL(wp), DIMENSION(jpij)     ::   zdqns_ice_b ! derivative of the surface flux function 
    117       REAL(wp), DIMENSION(jpij)     ::   zftrice     ! solar radiation transmitted through the ice 
    118        
     152      ! 
     153      REAL(wp), DIMENSION(jpij       )     ::   ztsuold     ! Old surface temperature in the ice 
    119154      REAL(wp), DIMENSION(jpij,nlay_i)     ::   ztiold      ! Old temperature in the ice 
    120155      REAL(wp), DIMENSION(jpij,nlay_s)     ::   ztsold      ! Old temperature in the snow 
     
    136171      REAL(wp), DIMENSION(jpij)            ::   zq_ini      ! diag errors on heat 
    137172      REAL(wp), DIMENSION(jpij)            ::   zghe        ! G(he), th. conduct enhancement factor, mono-cat 
    138        
     173      ! 
     174      REAL(wp) ::   zfr1, zfr2, zfrqsr_tr_i   ! dummy factor 
     175      ! 
    139176      ! Mono-category 
    140177      REAL(wp) :: zepsilon      ! determines thres. above which computation of G(h) is done 
     
    162199      ELSEWHERE                         ;   z1_h_i(1:npti) = 0._wp 
    163200      END WHERE 
    164  
     201      ! 
    165202      WHERE( zh_s(1:npti) >= epsi10 )   ;   z1_h_s(1:npti) = 1._wp / zh_s(1:npti) 
    166203      ELSEWHERE                         ;   z1_h_s(1:npti) = 0._wp 
    167204      END WHERE 
    168205      ! 
    169       ! temperatures 
    170       ztsub  (1:npti)   = t_su_1d(1:npti)  ! temperature at the previous iteration 
     206      ! Store initial temperatures and non solar heat fluxes 
     207      IF ( k_jules == np_zdf_jules_OFF .OR. k_jules == np_zdf_jules_SND  ) THEN ! OFF or SND mode 
     208         ! 
     209         ztsub  (1:npti)       =   t_su_1d(1:npti)                          ! surface temperature at iteration n-1 
     210         ztsuold(1:npti)       =   t_su_1d(1:npti)                          ! surface temperature initial value 
     211         zdqns_ice_b(1:npti)   =   dqns_ice_1d(1:npti)                      ! derivative of incoming nonsolar flux  
     212         zqns_ice_b (1:npti)   =   qns_ice_1d(1:npti)                       ! store previous qns_ice_1d value 
     213         ! 
     214         t_su_1d(1:npti)       =   MIN( t_su_1d(1:npti), rt0 - ztsu_err )   ! required to leave the choice between melting or not 
     215         ! 
     216      ENDIF    
     217      ! 
    171218      ztsold (1:npti,:) = t_s_1d(1:npti,:)   ! Old snow temperature 
    172219      ztiold (1:npti,:) = t_i_1d(1:npti,:)   ! Old ice temperature 
    173       t_su_1d(1:npti)   = MIN( t_su_1d(1:npti), rt0 - ztsu_err )   ! necessary 
    174       ! 
     220 
    175221      !------------- 
    176222      ! 2) Radiation 
    177223      !------------- 
    178       z1_hsu = 1._wp / 0.1_wp ! threshold for the computation of i0 
    179       DO ji = 1, npti 
    180          ! --- Computation of i0 --- ! 
    181          ! i0 describes the fraction of solar radiation which does not contribute 
    182          ! to the surface energy budget but rather penetrates inside the ice. 
    183          ! We assume that no radiation is transmitted through the snow 
    184          ! If there is no no snow 
    185          ! zfsw    = (1-i0).qsr_ice   is absorbed at the surface  
    186          ! zftrice = io.qsr_ice       is below the surface  
    187          ! ftr_ice = io.qsr_ice.exp(-k(h_i)) transmitted below the ice  
    188          ! fr1_i0_1d = i0 for a thin ice cover, fr1_i0_2d = i0 for a thick ice cover 
    189          zfac = MAX( 0._wp , 1._wp - ( h_i_1d(ji) * z1_hsu ) )      
    190          i0(ji) = ( 1._wp - isnow(ji) ) * ( fr1_i0_1d(ji) + zfac * fr2_i0_1d(ji) ) 
    191  
    192          ! --- Solar radiation absorbed / transmitted at the surface --- ! 
    193          !     Derivative of the non solar flux 
    194          zfsw   (ji)     =  qsr_ice_1d(ji) * ( 1 - i0(ji) )   ! Shortwave radiation absorbed at surface 
    195          zftrice(ji)     =  qsr_ice_1d(ji) *       i0(ji)     ! Solar radiation transmitted below the surface layer 
    196          zdqns_ice_b(ji) = dqns_ice_1d(ji)                    ! derivative of incoming nonsolar flux  
    197          zqns_ice_b (ji) =  qns_ice_1d(ji)                    ! store previous qns_ice_1d value 
    198       END DO 
    199  
    200224      ! --- Transmission/absorption of solar radiation in the ice --- ! 
    201       zradtr_s(1:npti,0) = zftrice(1:npti) 
     225!     zfr1 = ( 0.18 * ( 1.0 - 0.81 ) + 0.35 * 0.81 )            !   standard value 
     226!     zfr2 = ( 0.82 * ( 1.0 - 0.81 ) + 0.65 * 0.81 )            !   zfr2 such that zfr1 + zfr2 to equal 1 
     227! 
     228!     DO ji = 1, npti 
     229! 
     230!        zfac = MAX( 0._wp , 1._wp - ( h_i_1d(ji) * 10._wp ) )      
     231! 
     232!              zfrqsr_tr_i = zfr1 + zfac * zfr2                         !   below 10 cm, linearly increase zfrqsr_tr_i until 1 at zero thickness 
     233!              IF ( h_s_1d(ji) >= 0.0_wp )   zfrqsr_tr_i = 0._wp     !   snow fully opaque 
     234! 
     235!              qsr_ice_tr_1d(ji) = zfrqsr_tr_i * qsr_ice_1d(ji)   !   transmitted solar radiation  
     236! 
     237!              zfsw(ji) = qsr_ice_1d(ji) - qsr_ice_tr_1d(ji) 
     238!              zftrice(ji) = qsr_ice_tr_1d(ji) 
     239!              i0(ji) = zfrqsr_tr_i 
     240! 
     241!     END DO 
     242 
     243      zradtr_s(1:npti,0) = qsr_ice_tr_1d(1:npti) 
    202244      DO jk = 1, nlay_s 
    203245         DO ji = 1, npti 
     
    208250         END DO 
    209251      END DO 
    210           
    211       zradtr_i(1:npti,0) = zradtr_s(1:npti,nlay_s) * isnow(1:npti) + zftrice(1:npti) * ( 1._wp - isnow(1:npti) ) 
     252      ! 
     253      zradtr_i(1:npti,0) = zradtr_s(1:npti,nlay_s) * isnow(1:npti) + qsr_ice_tr_1d(1:npti) * ( 1._wp - isnow(1:npti) ) 
    212254      DO jk = 1, nlay_i  
    213255         DO ji = 1, npti 
     
    218260         END DO 
    219261      END DO 
    220  
     262      ! 
    221263      ftr_ice_1d(1:npti) = zradtr_i(1:npti,nlay_i)   ! record radiation transmitted below the ice 
    222264      ! 
     
    273315         ! 
    274316         SELECT CASE ( nn_monocat ) 
    275  
     317         ! 
    276318         CASE ( 1 , 3 ) 
    277  
     319            ! 
    278320            zepsilon = 0.1_wp 
    279321            DO ji = 1, npti 
     
    284326               ENDIF 
    285327            END DO 
    286  
     328            ! 
    287329         END SELECT 
    288330         ! 
     
    330372            END DO 
    331373         END DO 
    332          ! 
    333          !---------------------------- 
    334          ! 6) surface flux computation 
    335          !---------------------------- 
    336          IF ( ln_dqns_i ) THEN  
    337             DO ji = 1, npti 
    338                ! update of the non solar flux according to the update in T_su 
     374          
     375         ! 
     376         !----------------------------------------! 
     377         !                                        ! 
     378         !       JULES IF (OFF or SND MODE)       ! 
     379         !                                        ! 
     380         !----------------------------------------! 
     381         ! 
     382          
     383         IF ( k_jules == np_zdf_jules_OFF .OR. k_jules == np_zdf_jules_SND  ) THEN ! OFF or SND mode 
     384          
     385            ! 
     386            ! In OFF mode the original BL99 temperature computation is used 
     387            ! (with qsr_ice, qns_ice and dqns_ice as inputs) 
     388            ! 
     389            ! In SND mode, the computation is required to compute the conduction fluxes 
     390            ! 
     391          
     392            !---------------------------- 
     393            ! 6) surface flux computation 
     394            !---------------------------- 
     395 
     396            DO ji = 1, npti 
     397            ! update of the non solar flux according to the update in T_su 
    339398               qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) ) 
    340399            END DO 
    341          ENDIF 
    342  
    343          DO ji = 1, npti 
    344             zfnet(ji) = zfsw(ji) + qns_ice_1d(ji) ! incoming = net absorbed solar radiation + non solar total flux (LWup, LWdw, SH, LH) 
    345          END DO 
    346          ! 
    347          !---------------------------- 
    348          ! 7) tridiagonal system terms 
    349          !---------------------------- 
    350          !!layer denotes the number of the layer in the snow or in the ice 
    351          !!jm denotes the reference number of the equation in the tridiagonal 
    352          !!system, terms of tridiagonal system are indexed as following : 
    353          !!1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one 
    354  
    355          !!ice interior terms (top equation has the same form as the others) 
    356          ztrid   (1:npti,:,:) = 0._wp 
    357          zindterm(1:npti,:)   = 0._wp 
    358          zindtbis(1:npti,:)   = 0._wp 
    359          zdiagbis(1:npti,:)   = 0._wp 
    360  
    361          DO jm = nlay_s + 2, nlay_s + nlay_i  
     400 
     401            DO ji = 1, npti 
     402               zfnet(ji) = qsr_ice_1d(ji) - qsr_ice_tr_1d(ji) + qns_ice_1d(ji) ! net heat flux = net solar - transmitted solar + non solar 
     403            END DO 
     404            ! 
     405            !---------------------------- 
     406            ! 7) tridiagonal system terms 
     407            !---------------------------- 
     408            !!layer denotes the number of the layer in the snow or in the ice 
     409            !!jm denotes the reference number of the equation in the tridiagonal 
     410            !!system, terms of tridiagonal system are indexed as following : 
     411            !!1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one 
     412 
     413            !!ice interior terms (top equation has the same form as the others) 
     414            ztrid   (1:npti,:,:) = 0._wp 
     415            zindterm(1:npti,:)   = 0._wp 
     416            zindtbis(1:npti,:)   = 0._wp 
     417            zdiagbis(1:npti,:)   = 0._wp 
     418 
     419            DO jm = nlay_s + 2, nlay_s + nlay_i  
    362420            DO ji = 1, npti 
    363421               jk = jm - nlay_s - 1 
     
    367425               zindterm(ji,jm) = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk) 
    368426            END DO 
    369          ENDDO 
    370  
    371          jm =  nlay_s + nlay_i + 1 
    372          DO ji = 1, npti 
     427            END DO 
     428 
     429            jm =  nlay_s + nlay_i + 1 
     430            DO ji = 1, npti 
    373431            !!ice bottom term 
    374432            ztrid(ji,jm,1)  = - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1)    
     
    377435            zindterm(ji,jm) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) *  & 
    378436               &            ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 *  t_bo_1d(ji) )  
    379          ENDDO 
    380  
    381  
    382          DO ji = 1, npti 
     437            END DO 
     438 
     439 
     440            DO ji = 1, npti 
    383441            !                               !---------------------! 
    384             IF ( h_s_1d(ji) > 0.0 ) THEN   !  snow-covered cells ! 
     442            IF ( h_s_1d(ji) > 0.0 ) THEN    !  snow-covered cells ! 
    385443               !                            !---------------------! 
    386444               ! snow interior terms (bottom equation has the same form as the others) 
     
    428486                     &           ( zradab_s(ji,1) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) )  
    429487               ENDIF 
    430             !                               !---------------------! 
     488               !                            !---------------------! 
    431489            ELSE                            ! cells without snow  ! 
    432490               !                            !---------------------! 
     
    453511                     ztrid(ji,jm_min(ji),1)    = 0.0 
    454512                     ztrid(ji,jm_min(ji),2)    = zdqns_ice_b(ji) - zkappa_i(ji,0) * 2.0 
    455                      ztrid(ji,jm_min(ji),3)    = zkappa_i(ji,0) * 2.0 
     513                     ztrid(ji,jm_min(ji),3)    =  zkappa_i(ji,0) * 2.0 
    456514                     ztrid(ji,jm_min(ji)+1,1)  = -zkappa_i(ji,0) * 2.0 * zeta_i(ji,1) 
    457515                     ztrid(ji,jm_min(ji)+1,2)  = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) ) 
     
    488546            zdiagbis(ji,jm_min(ji)) = ztrid(ji,jm_min(ji),2) 
    489547            ! 
    490          END DO 
    491          ! 
    492          !------------------------------ 
    493          ! 8) tridiagonal system solving 
    494          !------------------------------ 
    495          ! Solve the tridiagonal system with Gauss elimination method. 
    496          ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984 
    497          jm_maxt = 0 
    498          jm_mint = nlay_i+5 
    499          DO ji = 1, npti 
    500             jm_mint = MIN(jm_min(ji),jm_mint) 
    501             jm_maxt = MAX(jm_max(ji),jm_maxt) 
    502          END DO 
    503  
    504          DO jk = jm_mint+1, jm_maxt 
     548            END DO 
     549            ! 
     550            !------------------------------ 
     551            ! 8) tridiagonal system solving 
     552            !------------------------------ 
     553            ! Solve the tridiagonal system with Gauss elimination method. 
     554            ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984 
     555            jm_maxt = 0 
     556            jm_mint = nlay_i+5 
     557            DO ji = 1, npti 
     558               jm_mint = MIN(jm_min(ji),jm_mint) 
     559               jm_maxt = MAX(jm_max(ji),jm_maxt) 
     560            END DO 
     561 
     562            DO jk = jm_mint+1, jm_maxt 
    505563            DO ji = 1, npti 
    506564               jm = min(max(jm_min(ji)+1,jk),jm_max(ji)) 
     
    508566               zindtbis(ji,jm) = zindterm(ji,jm) - ztrid(ji,jm,1) * zindtbis(ji,jm-1) / zdiagbis(ji,jm-1) 
    509567            END DO 
    510          END DO 
    511  
    512          DO ji = 1, npti 
     568            END DO 
     569 
     570            DO ji = 1, npti 
    513571            ! ice temperatures 
    514572            t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji)) 
    515          END DO 
    516  
    517          DO jm = nlay_i + nlay_s, nlay_s + 2, -1 
     573            END DO 
     574 
     575            DO jm = nlay_i + nlay_s, nlay_s + 2, -1 
    518576            DO ji = 1, npti 
    519577               jk    =  jm - nlay_s - 1 
    520578               t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm) 
    521579            END DO 
    522          END DO 
    523  
    524          DO ji = 1, npti 
     580            END DO 
     581 
     582            DO ji = 1, npti 
    525583            ! snow temperatures       
    526584            IF( h_s_1d(ji) > 0._wp ) THEN 
    527585               t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) )  & 
    528                   &                / zdiagbis(ji,nlay_s+1) 
     586               &                / zdiagbis(ji,nlay_s+1) 
    529587            ENDIF 
    530588            ! surface temperature 
     
    534592                  &          ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji)) 
    535593            ENDIF 
    536          END DO 
    537          ! 
    538          !-------------------------------------------------------------- 
    539          ! 9) Has the scheme converged ?, end of the iterative procedure 
    540          !-------------------------------------------------------------- 
    541          ! check that nowhere it has started to melt 
    542          ! zdti_max is a measure of error, it has to be under zdti_bnd 
    543          zdti_max = 0._wp 
    544          DO ji = 1, npti 
     594            END DO 
     595            ! 
     596            !-------------------------------------------------------------- 
     597            ! 9) Has the scheme converged ?, end of the iterative procedure 
     598            !-------------------------------------------------------------- 
     599            ! check that nowhere it has started to melt 
     600            ! zdti_max is a measure of error, it has to be under zdti_bnd 
     601            zdti_max = 0._wp 
     602            DO ji = 1, npti 
    545603            t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , rt0 - 100._wp ) 
    546604            zdti_max = MAX( zdti_max, ABS( t_su_1d(ji) - ztsub(ji) ) )      
    547          END DO 
    548  
    549          DO jk  =  1, nlay_s 
     605            END DO 
     606 
     607            DO jk  =  1, nlay_s 
    550608            DO ji = 1, npti 
    551609               t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp ) 
    552610               zdti_max = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) ) 
    553611            END DO 
    554          END DO 
    555  
    556          DO jk  =  1, nlay_i 
     612            END DO 
     613 
     614            DO jk  =  1, nlay_i 
    557615            DO ji = 1, npti 
    558616               ztmelt_i      = -tmut * sz_i_1d(ji,jk) + rt0  
     
    560618               zdti_max      =  MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) 
    561619            END DO 
    562          END DO 
    563  
    564          ! Compute spatial maximum over all errors 
    565          ! note that this could be optimized substantially by iterating only the non-converging points 
    566          IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice ) 
    567  
     620            END DO 
     621 
     622            ! Compute spatial maximum over all errors 
     623            ! note that this could be optimized substantially by iterating only the non-converging points 
     624            IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice ) 
     625         ! 
     626         !----------------------------------------! 
     627         !                                        ! 
     628         !       JULES IF (RCV MODE)              ! 
     629         !                                        ! 
     630         !----------------------------------------! 
     631         ! 
     632 
     633         ELSE IF ( k_jules == np_zdf_jules_RCV  ) THEN ! RCV mode 
     634          
     635            ! 
     636            ! In RCV mode, we use a modified BL99 solver  
     637            ! with conduction flux (qcn_ice) as forcing term 
     638            ! 
     639            !---------------------------- 
     640            ! 7) tridiagonal system terms 
     641            !---------------------------- 
     642            !!layer denotes the number of the layer in the snow or in the ice 
     643            !!jm denotes the reference number of the equation in the tridiagonal 
     644            !!system, terms of tridiagonal system are indexed as following : 
     645            !!1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one 
     646 
     647            !!ice interior terms (top equation has the same form as the others) 
     648            ztrid   (1:npti,:,:) = 0._wp 
     649            zindterm(1:npti,:)   = 0._wp 
     650            zindtbis(1:npti,:)   = 0._wp 
     651            zdiagbis(1:npti,:)   = 0._wp 
     652 
     653            DO jm = nlay_s + 2, nlay_s + nlay_i  
     654            DO ji = 1, npti 
     655               jk = jm - nlay_s - 1 
     656               ztrid(ji,jm,1)  = - zeta_i(ji,jk) * zkappa_i(ji,jk-1) 
     657               ztrid(ji,jm,2)  = 1.0 + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) ) 
     658               ztrid(ji,jm,3)  = - zeta_i(ji,jk) * zkappa_i(ji,jk) 
     659               zindterm(ji,jm) = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk) 
     660            END DO 
     661            ENDDO 
     662 
     663            jm =  nlay_s + nlay_i + 1 
     664            DO ji = 1, npti 
     665            !!ice bottom term 
     666            ztrid(ji,jm,1)  = - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1)    
     667            ztrid(ji,jm,2)  = 1.0 + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i) * zg1 + zkappa_i(ji,nlay_i-1) ) 
     668            ztrid(ji,jm,3)  = 0.0 
     669            zindterm(ji,jm) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) *  & 
     670               &            ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 *  t_bo_1d(ji) )  
     671            ENDDO 
     672 
     673 
     674            DO ji = 1, npti 
     675            !                              !---------------------! 
     676            IF ( h_s_1d(ji) > 0.0 ) THEN   !  snow-covered cells ! 
     677               !                           !---------------------! 
     678               ! snow interior terms (bottom equation has the same form as the others) 
     679               DO jm = 3, nlay_s + 1 
     680               jk = jm - 1 
     681               ztrid(ji,jm,1)  = - zeta_s(ji,jk) * zkappa_s(ji,jk-1) 
     682               ztrid(ji,jm,2)  = 1.0 + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) ) 
     683               ztrid(ji,jm,3)  = - zeta_s(ji,jk)*zkappa_s(ji,jk) 
     684               zindterm(ji,jm) = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk) 
     685               END DO 
     686 
     687               ! case of only one layer in the ice (ice equation is altered) 
     688               IF ( nlay_i == 1 ) THEN 
     689               ztrid(ji,nlay_s+2,3)  = 0.0 
     690               zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1) * t_bo_1d(ji)  
     691               ENDIF 
     692 
     693               jm_min(ji) = 2 
     694               jm_max(ji) = nlay_i + nlay_s + 1 
     695 
     696               ! first layer of snow equation 
     697               ztrid(ji,2,1)  = 0.0 
     698               ztrid(ji,2,2)  = 1.0 + zeta_s(ji,1) * zkappa_s(ji,1) 
     699               ztrid(ji,2,3)  = - zeta_s(ji,1)*zkappa_s(ji,1)  
     700               zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) *  & 
     701               &           ( zradab_s(ji,1) + qcn_ice_1d(ji) )  
     702 
     703            !                               !---------------------! 
     704            ELSE                            ! cells without snow  ! 
     705            !                               !---------------------! 
     706 
     707               jm_min(ji)    =  nlay_s + 2 
     708               jm_max(ji)    =  nlay_i + nlay_s + 1 
     709 
     710               ! first layer of ice equation 
     711               ztrid(ji,jm_min(ji),1)  = 0.0 
     712               ztrid(ji,jm_min(ji),2)  = 1.0 + zeta_i(ji,1) * zkappa_i(ji,1)   
     713               ztrid(ji,jm_min(ji),3)  = - zeta_i(ji,1) * zkappa_i(ji,1) 
     714               zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) *  & 
     715               &                    ( zradab_i(ji,1) + qcn_ice_1d(ji) ) 
     716 
     717               ! case of only one layer in the ice (surface & ice equations are altered) 
     718               IF ( nlay_i == 1 ) THEN 
     719               ztrid(ji,jm_min(ji),1)  = 0.0 
     720               ztrid(ji,jm_min(ji),2)  = 1.0 + zeta_i(ji,1) * zkappa_i(ji,1) 
     721               ztrid(ji,jm_min(ji),3)  = 0.0 
     722               zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji)    & 
     723               &                    + qcn_ice_1d(ji) ) 
     724 
     725               ENDIF 
     726 
     727            ENDIF 
     728            ! 
     729            zindtbis(ji,jm_min(ji)) = zindterm(ji,jm_min(ji)) 
     730            zdiagbis(ji,jm_min(ji)) = ztrid(ji,jm_min(ji),2) 
     731            ! 
     732            END DO 
     733            ! 
     734            !------------------------------ 
     735            ! 8) tridiagonal system solving 
     736            !------------------------------ 
     737            ! Solve the tridiagonal system with Gauss elimination method. 
     738            ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984 
     739            jm_maxt = 0 
     740            jm_mint = nlay_i+5 
     741            DO ji = 1, npti 
     742            jm_mint = MIN(jm_min(ji),jm_mint) 
     743            jm_maxt = MAX(jm_max(ji),jm_maxt) 
     744            END DO 
     745 
     746            DO jk = jm_mint+1, jm_maxt 
     747            DO ji = 1, npti 
     748               jm = min(max(jm_min(ji)+1,jk),jm_max(ji)) 
     749               zdiagbis(ji,jm) = ztrid(ji,jm,2)  - ztrid(ji,jm,1) * ztrid(ji,jm-1,3)  / zdiagbis(ji,jm-1) 
     750               zindtbis(ji,jm) = zindterm(ji,jm) - ztrid(ji,jm,1) * zindtbis(ji,jm-1) / zdiagbis(ji,jm-1) 
     751            END DO 
     752            END DO 
     753 
     754            DO ji = 1, npti 
     755            ! ice temperatures 
     756            t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji)) 
     757            END DO 
     758 
     759            DO jm = nlay_i + nlay_s, nlay_s + 2, -1 
     760            DO ji = 1, npti 
     761               jk    =  jm - nlay_s - 1 
     762               t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm) 
     763            END DO 
     764            END DO 
     765 
     766            DO ji = 1, npti 
     767            ! snow temperatures       
     768            IF( h_s_1d(ji) > 0._wp ) THEN 
     769               t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) )  & 
     770               &                / zdiagbis(ji,nlay_s+1) 
     771            ENDIF 
     772            END DO 
     773            ! 
     774            !-------------------------------------------------------------- 
     775            ! 9) Has the scheme converged ?, end of the iterative procedure 
     776            !-------------------------------------------------------------- 
     777            ! check that nowhere it has started to melt 
     778            ! zdti_max is a measure of error, it has to be under zdti_bnd 
     779            zdti_max = 0._wp 
     780 
     781            DO jk  =  1, nlay_s 
     782            DO ji = 1, npti 
     783               t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp ) 
     784               zdti_max = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) ) 
     785            END DO 
     786            END DO 
     787 
     788            DO jk  =  1, nlay_i 
     789            DO ji = 1, npti 
     790               ztmelt_i      = -tmut * sz_i_1d(ji,jk) + rt0  
     791               t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), rt0 - 100._wp ) 
     792               zdti_max      =  MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) 
     793            END DO 
     794            END DO 
     795 
     796            ! Compute spatial maximum over all errors 
     797            ! note that this could be optimized substantially by iterating only the non-converging points 
     798            IF( lk_mpp )   CALL mpp_max( zdti_max, kcom=ncomm_ice ) 
     799 
     800         ENDIF ! k_jules 
     801          
    568802      END DO  ! End of the do while iterative procedure 
    569  
     803       
    570804      IF( ln_icectl .AND. lwp ) THEN 
    571805         WRITE(numout,*) ' zdti_max : ', zdti_max 
    572806         WRITE(numout,*) ' iconv    : ', iconv 
    573807      ENDIF 
     808       
    574809      ! 
    575810      !----------------------------- 
    576811      ! 10) Fluxes at the interfaces 
    577812      !----------------------------- 
     813      ! 
     814      ! --- update conduction fluxes 
     815      ! 
    578816      DO ji = 1, npti 
    579          !                                ! surface ice conduction flux 
     817      !                                ! surface ice conduction flux 
    580818         fc_su(ji)   =  -           isnow(ji)   * zkappa_s(ji,0) * zg1s * (t_s_1d(ji,1) - t_su_1d(ji))   & 
    581819            &           - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1  * (t_i_1d(ji,1) - t_su_1d(ji)) 
    582          !                                ! bottom ice conduction flux 
     820      !                                ! bottom ice conduction flux 
    583821         fc_bo_i(ji) =  - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_1d(ji) - t_i_1d(ji,nlay_i)) ) 
    584822      END DO 
    585  
    586       ! --- computes sea ice energy of melting compulsory for icethd_dh --- ! 
    587       CALL ice_thd_enmelt 
    588  
    589       ! --- diagnose the change in non-solar flux due to surface temperature change --- ! 
    590       IF ( ln_dqns_i ) THEN 
     823       
     824      ! 
     825      ! --- Diagnose the heat loss due to changing non-solar / conduction flux --- ! 
     826      ! 
     827      IF ( k_jules == np_zdf_jules_OFF .OR. k_jules == np_zdf_jules_SND  ) THEN  ! OFF or SND mode 
    591828         DO ji = 1, npti 
    592             hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qns_ice_1d(ji)  - zqns_ice_b(ji) ) * a_i_1d(ji)  
    593          END DO 
    594       END IF 
    595  
    596       ! --- diag conservation imbalance on heat diffusion - PART 2 --- ! 
    597       !     hfx_dif = Heat flux used to warm/cool ice in W.m-2 
    598       !     zhfx_err = correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation 
    599       DO ji = 1, npti 
    600          zdq = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * h_i_1d(ji) * r1_nlay_i +  & 
    601             &                   SUM( e_s_1d(ji,1:nlay_s) ) * h_s_1d(ji) * r1_nlay_s ) 
    602  
    603          IF( t_su_1d(ji) < rt0 ) THEN  ! case T_su < 0degC 
    604             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)  
    605          ELSE                          ! case T_su = 0degC 
    606             zhfx_err = ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji) 
    607          ENDIF 
    608          hfx_dif_1d(ji) = hfx_dif_1d(ji) - zdq * r1_rdtice * a_i_1d(ji) 
    609  
    610          ! total heat that is sent to the ocean (i.e. not used in the heat diffusion equation) 
    611          hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err 
    612  
    613       END DO    
    614  
    615       ! --- Diagnostics SIMIP --- ! 
    616       DO ji = 1, npti 
    617          !--- Conduction fluxes (positive downwards) 
    618          diag_fc_bo_1d(ji) = diag_fc_bo_1d(ji) + fc_bo_i(ji) * a_i_1d(ji) / at_i_1d(ji) 
    619          diag_fc_su_1d(ji) = diag_fc_su_1d(ji) + fc_su(ji)   * a_i_1d(ji) / at_i_1d(ji) 
    620  
    621          !--- Snow-ice interfacial temperature (diagnostic SIMIP) 
    622          zfac = rn_cnd_s * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji) 
    623          IF( zh_s(ji) >= 1.e-3 .AND. zfac > epsi10 ) THEN 
    624             t_si_1d(ji) = ( rn_cnd_s       * zh_i(ji) * t_s_1d(ji,1) +   & 
    625                &            ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) / zfac 
    626          ELSE 
    627             t_si_1d(ji) = t_su_1d(ji) 
    628          ENDIF 
    629       END DO 
    630       ! 
    631    END SUBROUTINE ice_thd_zdf 
     829            hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qns_ice_1d(ji) - zqns_ice_b(ji) ) * a_i_1d(ji)  
     830         END DO 
     831      ELSE        ! RCV mode 
     832         DO ji = 1, npti 
     833            hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( fc_su(ji)      - qcn_ice_1d(ji) ) * a_i_1d(ji)  
     834         END DO 
     835      ENDIF 
     836       
     837      ! 
     838      ! --- Diagnose the heat loss due to non-fully converged temperature solution (should not be above 10-4 W-m2) --- ! 
     839      ! 
     840 
     841      IF ( ( k_jules == np_zdf_jules_OFF ) .OR. ( k_jules == np_zdf_jules_RCV ) ) THEN ! OFF 
     842       
     843         CALL ice_thd_enmelt        
     844 
     845         !     zhfx_err = correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation 
     846         DO ji = 1, npti 
     847            zdq = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * h_i_1d(ji) * r1_nlay_i +  & 
     848               &                   SUM( e_s_1d(ji,1:nlay_s) ) * h_s_1d(ji) * r1_nlay_s ) 
     849                
     850            IF ( ( k_jules == np_zdf_jules_OFF ) ) THEN 
     851                
     852                  IF( t_su_1d(ji) < rt0 ) THEN  ! case T_su < 0degC 
     853                     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)  
     854                  ELSE                          ! case T_su = 0degC 
     855                     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) 
     856                  ENDIF 
     857             
     858            ELSE ! RCV CASE 
     859             
     860               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) 
     861             
     862            ENDIF 
     863            ! 
     864            ! total heat sink to be sent to the ocean 
     865            hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err 
     866            ! 
     867            ! hfx_dif = Heat flux diagnostic of sensible heat used to warm/cool ice in W.m-2    
     868            hfx_dif_1d(ji) = hfx_dif_1d(ji) - zdq * r1_rdtice * a_i_1d(ji) 
     869            ! 
     870         END DO 
     871         ! 
     872         ! --- SIMIP diagnostics 
     873         ! 
     874         DO ji = 1, npti 
     875            !--- Conduction fluxes (positive downwards) 
     876            diag_fc_bo_1d(ji) = diag_fc_bo_1d(ji) + fc_bo_i(ji) * a_i_1d(ji) / at_i_1d(ji) 
     877            diag_fc_su_1d(ji) = diag_fc_su_1d(ji) + fc_su(ji)   * a_i_1d(ji) / at_i_1d(ji) 
     878    
     879            !--- Snow-ice interfacial temperature (diagnostic SIMIP) 
     880            zfac = rn_cnd_s * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji) 
     881            IF( zh_s(ji) >= 1.e-3 .AND. zfac > epsi10 ) THEN 
     882               t_si_1d(ji) = ( rn_cnd_s       * zh_i(ji) * t_s_1d(ji,1) +   & 
     883                  &            ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) / zfac 
     884            ELSE 
     885               t_si_1d(ji) = t_su_1d(ji) 
     886            ENDIF 
     887         END DO 
     888         ! 
     889      ENDIF 
     890      ! 
     891      !--------------------------------------------------------------------------------------- 
     892      ! 11) Jules coupling: reset inner snow and ice temperatures, update conduction fluxes 
     893      !--------------------------------------------------------------------------------------- 
     894      ! 
     895      IF ( k_jules == np_zdf_jules_SND ) THEN   ! --- Jules coupling in "SND" mode 
     896         ! 
     897         ! Restore temperatures to their initial values 
     898         t_s_1d(1:npti,:)        = ztsold (1:npti,:) 
     899         t_i_1d(1:npti,:)        = ztiold (1:npti,:) 
     900         qcn_ice_1d(1:npti)      = fc_su(1:npti) 
     901         !   
     902      ENDIF 
     903      ! 
     904   END SUBROUTINE ice_thd_zdf_BL99 
    632905 
    633906 
     
    675948      !! ** input   :   Namelist namthd_zdf 
    676949      !!------------------------------------------------------------------- 
    677       INTEGER  ::   ios   ! Local integer output status for namelist read 
     950      INTEGER  ::   ios, ioptio   ! Local integer 
    678951      !! 
    679       NAMELIST/namthd_zdf/ ln_zdf_BL99, ln_cndi_U64, ln_cndi_P07, rn_cnd_s, rn_kappa_i, ln_dqns_i 
     952      NAMELIST/namthd_zdf/ ln_zdf_BL99, ln_cndi_U64, ln_cndi_P07, rn_cnd_s, rn_kappa_i 
    680953      !!------------------------------------------------------------------- 
    681954      ! 
     
    694967         WRITE(numout,*) '~~~~~~~~~~~~~~~~' 
    695968         WRITE(numout,*) '   Namelist namthd_zdf:' 
    696          WRITE(numout,*) '      Diffusion follows a Bitz and Lipscomb (1999)            ln_zdf_BL99  = ', ln_zdf_BL99 
     969         WRITE(numout,*) '      Bitz and Lipscomb (1999) formulation                    ln_zdf_BL99  = ', ln_zdf_BL99 
    697970         WRITE(numout,*) '      thermal conductivity in the ice (Untersteiner 1964)     ln_cndi_U64  = ', ln_cndi_U64 
    698971         WRITE(numout,*) '      thermal conductivity in the ice (Pringle et al 2007)    ln_cndi_P07  = ', ln_cndi_P07 
    699972         WRITE(numout,*) '      thermal conductivity in the snow                        rn_cnd_s     = ', rn_cnd_s 
    700973         WRITE(numout,*) '      extinction radiation parameter in sea ice               rn_kappa_i   = ', rn_kappa_i 
    701          WRITE(numout,*) '      change the surface non-solar flux with Tsu or not       ln_dqns_i    = ', ln_dqns_i 
    702974     ENDIF 
     975 
    703976      ! 
    704977      IF ( ( ln_cndi_U64 .AND. ln_cndi_P07 ) .OR. ( .NOT.ln_cndi_U64 .AND. .NOT.ln_cndi_P07 ) ) THEN 
    705978         CALL ctl_stop( 'ice_thd_zdf_init: choose one and only one formulation for thermal conduction (ln_cndi_U64 or ln_cndi_P07)' ) 
    706979      ENDIF 
     980 
     981      !                             !== set the choice of ice vertical thermodynamic formulation ==! 
     982      ioptio = 0  
     983      !      !--- BL99 thermo dynamics                               (linear liquidus + constant thermal properties) 
     984      IF( ln_zdf_BL99 ) THEN   ;   ioptio = ioptio + 1   ;   nice_zdf = np_BL99   ;   ENDIF 
     985      IF( ioptio /= 1 )    CALL ctl_stop( 'ice_thd_init: one and only one ice thermo option has to be defined ' ) 
    707986      ! 
    708987   END SUBROUTINE ice_thd_zdf_init 
  • branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90

    r8637 r8813  
    4444   !!   ice_var_salprof1d : salinity profile in the ice 1D 
    4545   !!   ice_var_zapsmall  : remove very small area and volume 
    46    !!   ice_var_itd       : convert 1-cat to multiple cat 
     46   !!   ice_var_itd       : convert 1-cat to jpl-cat 
     47   !!   ice_var_itd2      : convert N-cat to jpl-cat 
    4748   !!   ice_var_bv        : brine volume 
    4849   !!---------------------------------------------------------------------- 
     
    6768   PUBLIC   ice_var_zapsmall 
    6869   PUBLIC   ice_var_itd 
     70   PUBLIC   ice_var_itd2 
    6971   PUBLIC   ice_var_bv            
    7072 
     
    9698      et_s(:,:)  = SUM( SUM( e_s(:,:,:,:), dim=4 ), dim=3 ) 
    9799      et_i(:,:)  = SUM( SUM( e_i(:,:,:,:), dim=4 ), dim=3 ) 
    98  
     100      ! 
    99101      at_ip(:,:) = SUM( a_ip(:,:,:), dim=3 ) ! melt ponds 
    100102      vt_ip(:,:) = SUM( v_ip(:,:,:), dim=3 ) 
    101  
    102       ato_i(:,:) = 1._wp - at_i(:,:)    ! open water fraction   
     103      ! 
     104      ato_i(:,:) = 1._wp - at_i(:,:)         ! open water fraction   
    103105 
    104106      IF( kn > 1 ) THEN 
     
    238240         END WHERE 
    239241      END DO 
    240  
     242      ! 
    241243      ! integrated values  
    242244      vt_i (:,:) = SUM( v_i, dim=3 ) 
     
    322324            END DO 
    323325         END DO 
    324  
     326         ! 
    325327         ! Computation of the profile 
    326328         DO jl = 1, jpl 
     
    362364   END SUBROUTINE ice_var_salprof 
    363365 
     366 
    364367   SUBROUTINE ice_var_salprof1d 
    365368      !!------------------------------------------------------------------- 
     
    457460         ELSEWHERE                                                                           ;   zswitch(:,:) = 1._wp 
    458461         END WHERE 
    459           
     462         ! 
    460463         DO jk = 1, nlay_i 
    461464            DO jj = 1 , jpj 
     
    468471            END DO 
    469472         END DO 
    470  
     473         ! 
    471474         DO jj = 1 , jpj 
    472475            DO ji = 1 , jpi 
     
    484487               t_s(ji,jj,1,jl) = t_s(ji,jj,1,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) ) 
    485488               e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * zswitch(ji,jj) 
    486  
     489               ! 
    487490               !----------------------------------------------------------------- 
    488491               ! zap ice and snow volume, add water and salt to ocean 
     
    495498               oa_i (ji,jj,jl) = oa_i(ji,jj,jl) * zswitch(ji,jj) 
    496499               sv_i (ji,jj,jl) = sv_i(ji,jj,jl) * zswitch(ji,jj) 
    497  
     500               ! 
    498501               h_i (ji,jj,jl) = h_i (ji,jj,jl) * zswitch(ji,jj) 
    499502               h_s (ji,jj,jl) = h_s (ji,jj,jl) * zswitch(ji,jj) 
    500  
     503               ! 
    501504               a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * zswitch(ji,jj) 
    502505               v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj) 
    503  
    504             END DO 
    505          END DO 
    506          
     506               ! 
     507            END DO 
     508         END DO 
     509         ! 
    507510      END DO  
    508511 
     
    548551      !!------------------------------------------------------------------- 
    549552      INTEGER  :: ji, jk, jl             ! dummy loop indices 
    550       INTEGER  :: ijpij, i_fill, jl0   
     553      INTEGER  :: idim, i_fill, jl0   
    551554      REAL(wp) :: zarg, zV, zconv, zdh, zdv 
    552555      REAL(wp), DIMENSION(:),   INTENT(in)    ::   zhti, zhts, zati    ! input ice/snow variables 
     
    561564      ! then we check whether the distribution fullfills 
    562565      ! volume and area conservation, positivity and ice categories bounds 
    563       ijpij = SIZE( zhti , 1 ) 
    564       zh_i(1:ijpij,1:jpl) = 0._wp 
    565       zh_s(1:ijpij,1:jpl) = 0._wp 
    566       za_i (1:ijpij,1:jpl) = 0._wp 
    567  
    568       DO ji = 1, ijpij 
    569           
     566      idim = SIZE( zhti , 1 ) 
     567      zh_i(1:idim,1:jpl) = 0._wp 
     568      zh_s(1:idim,1:jpl) = 0._wp 
     569      za_i(1:idim,1:jpl) = 0._wp 
     570      ! 
     571      DO ji = 1, idim 
     572         ! 
    570573         IF( zhti(ji) > 0._wp ) THEN 
    571  
     574            ! 
    572575            ! find which category (jl0) the input ice thickness falls into 
    573576            jl0 = jpl 
     
    578581               ENDIF 
    579582            END DO 
    580  
     583            ! 
    581584            itest(:) = 0 
    582585            i_fill   = jpl + 1                                            !------------------------------------ 
     
    586589               ! 
    587590               zh_i(ji,1:jpl) = 0._wp 
    588                za_i (ji,1:jpl) = 0._wp 
    589                itest(:)        = 0       
    590                 
     591               za_i(ji,1:jpl) = 0._wp 
     592               itest(:)       = 0       
     593               ! 
    591594               IF ( i_fill == 1 ) THEN      !-- case very thin ice: fill only category 1 
    592595                  zh_i(ji,1) = zhti(ji) 
     
    597600                     zh_i(ji,jl) = hi_mean(jl) 
    598601                  END DO 
    599                    
     602                  ! 
    600603                  ! concentration 
    601604                  za_i(ji,jl0) = zati(ji) / SQRT(REAL(jpl)) 
     
    606609                     ENDIF 
    607610                  END DO 
    608                    
     611                  ! 
    609612                  ! last category 
    610613                  za_i(ji,i_fill) = zati(ji) - SUM( za_i(ji,1:i_fill-1) ) 
    611614                  zV = SUM( za_i(ji,1:i_fill-1) * zh_i(ji,1:i_fill-1) ) 
    612615                  zh_i(ji,i_fill) = ( zhti(ji) * zati(ji) - zV ) / MAX( za_i(ji,i_fill), epsi10 )  
    613                    
     616                  ! 
    614617                  ! clem: correction if concentration of upper cat is greater than lower cat 
    615618                  !       (it should be a gaussian around jl0 but sometimes it is not) 
     
    622625                           za_i (ji,1:jl-1) = za_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * zhti(ji), epsi10 ) 
    623626                        END IF 
    624                      ENDDO 
     627                     END DO 
    625628                  ENDIF 
    626                 
     629                  ! 
    627630               ENDIF 
    628              
     631               ! 
    629632               ! Compatibility tests 
    630633               zconv = ABS( zati(ji) - SUM( za_i(ji,1:jpl) ) )  
    631                IF ( zconv < epsi06 ) itest(1) = 1                                        ! Test 1: area conservation 
    632              
     634               IF ( zconv < epsi06 )   itest(1) = 1                                        ! Test 1: area conservation 
     635               ! 
    633636               zconv = ABS( zhti(ji)*zati(ji) - SUM( za_i(ji,1:jpl)*zh_i(ji,1:jpl) ) ) 
    634                IF ( zconv < epsi06 ) itest(2) = 1                                        ! Test 2: volume conservation 
    635                 
    636                IF ( zh_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1                  ! Test 3: thickness of the last category is in-bounds ? 
    637                 
     637               IF ( zconv < epsi06 )   itest(2) = 1                                        ! Test 2: volume conservation 
     638               ! 
     639               IF ( zh_i(ji,i_fill) >= hi_max(i_fill-1) )   itest(3) = 1                  ! Test 3: thickness of the last category is in-bounds ? 
     640               ! 
    638641               itest(4) = 1 
    639642               DO jl = 1, i_fill 
     
    642645               !                                         !---------------------------- 
    643646            END DO                                       ! end iteration on categories 
    644                !                                         !---------------------------- 
     647            !                                            !---------------------------- 
    645648         ENDIF 
    646649      END DO 
     
    648651      ! Add Snow in each category where za_i is not 0 
    649652      DO jl = 1, jpl 
    650          DO ji = 1, ijpij 
     653         DO ji = 1, idim 
    651654            IF( za_i(ji,jl) > 0._wp ) THEN 
    652655               zh_s(ji,jl) = zh_i(ji,jl) * ( zhts(ji) / zhti(ji) ) 
     
    661664      END DO 
    662665      ! 
    663     END SUBROUTINE ice_var_itd 
    664  
    665  
    666     SUBROUTINE ice_var_bv 
     666   END SUBROUTINE ice_var_itd 
     667 
     668 
     669   SUBROUTINE ice_var_itd2( zhti, zhts, zati, zh_i, zh_s, za_i ) 
     670      !!------------------------------------------------------------------- 
     671      !!                ***  ROUTINE ice_var_itd2   *** 
     672      !! 
     673      !! ** Purpose :  converting N-cat ice to jpl ice categories 
     674      !! 
     675      !!                  ice thickness distribution follows a gaussian law 
     676      !!               around the concentration of the most likely ice thickness 
     677      !!                           (similar as iceistate.F90) 
     678      !! 
     679      !! ** Method:   Iterative procedure 
     680      !!                 
     681      !!               1) Fill ice cat that correspond to input thicknesses 
     682      !!                  Find the lowest(jlmin) and highest(jlmax) cat that are filled 
     683      !! 
     684      !!               2) Expand the filling to the cat jlmin-1 and jlmax+1 
     685      !!                   by removing 25% ice area from jlmin and jlmax (resp.)  
     686      !!               
     687      !!               3) Expand the filling to the empty cat between jlmin and jlmax  
     688      !!                   by a) removing 25% ice area from the lower cat (ascendant loop jlmin=>jlmax) 
     689      !!                      b) removing 25% ice area from the higher cat (descendant loop jlmax=>jlmin) 
     690      !! 
     691      !! ** Arguments : zhti: N-cat ice thickness 
     692      !!                zhts: N-cat snow depth 
     693      !!                zati: N-cat ice concentration 
     694      !! 
     695      !! ** Output    : jpl-cat  
     696      !! 
     697      !!  (Example of application: BDY forcings when inputs have N-cat /= jpl)   
     698      !!------------------------------------------------------------------- 
     699      INTEGER  ::   ji, jl, jl1, jl2             ! dummy loop indices 
     700      INTEGER  ::   idim, icat   
     701      INTEGER, PARAMETER ::   ztrans = 0.25_wp 
     702      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   zhti, zhts, zati    ! input ice/snow variables 
     703      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zh_i, zh_s, za_i    ! output ice/snow variables 
     704      INTEGER , DIMENSION(:,:), ALLOCATABLE   ::   jlfil, jlfil2 
     705      INTEGER , DIMENSION(:)  , ALLOCATABLE   ::   jlmax, jlmin 
     706      !!------------------------------------------------------------------- 
     707      ! 
     708      idim = SIZE( zhti, 1 ) 
     709      icat = SIZE( zhti, 2 ) 
     710      ! 
     711      ALLOCATE( jlfil(idim,jpl), jlfil2(idim,jpl) )       ! allocate arrays 
     712      ALLOCATE( jlmin(idim), jlmax(idim) ) 
     713 
     714      ! --- initialize output fields to 0 --- ! 
     715      zh_i(1:idim,1:jpl) = 0._wp 
     716      zh_s(1:idim,1:jpl) = 0._wp 
     717      za_i(1:idim,1:jpl) = 0._wp 
     718      ! 
     719      ! --- fill the categories --- ! 
     720      !     find where cat-input = cat-output and fill cat-output fields   
     721      jlmax(:) = 0 
     722      jlmin(:) = 999 
     723      jlfil(:,:) = 0 
     724      DO jl1 = 1, jpl 
     725         DO jl2 = 1, icat 
     726            DO ji = 1, idim 
     727               IF( hi_max(jl1-1) <= zhti(ji,jl2) .AND. hi_max(jl1) > zhti(ji,jl2) ) THEN 
     728                  ! fill the right category 
     729                  zh_i(ji,jl1) = zhti(ji,jl2) 
     730                  zh_s(ji,jl1) = zhts(ji,jl2) 
     731                  za_i(ji,jl1) = zati(ji,jl2) 
     732                  ! record categories that are filled 
     733                  jlmax(ji) = MAX( jlmax(ji), jl1 ) 
     734                  jlmin(ji) = MIN( jlmin(ji), jl1 ) 
     735                  jlfil(ji,jl1) = jl1 
     736               ENDIF 
     737            END DO 
     738         END DO 
     739      END DO 
     740      ! 
     741      ! --- fill the gaps between categories --- !   
     742      !     transfer from categories filled at the previous step to the empty ones in between 
     743      DO ji = 1, idim 
     744         jl1 = jlmin(ji) 
     745         jl2 = jlmax(ji) 
     746         IF( jl1 > 1 ) THEN 
     747            ! fill the lower cat (jl1-1) 
     748            za_i(ji,jl1-1) = ztrans * za_i(ji,jl1) 
     749            zh_i(ji,jl1-1) = hi_mean(jl1-1) 
     750            ! remove from cat jl1 
     751            za_i(ji,jl1  ) = ( 1._wp - ztrans ) * za_i(ji,jl1) 
     752         ENDIF 
     753         IF( jl2 < jpl ) THEN 
     754            ! fill the upper cat (jl2+1) 
     755            za_i(ji,jl2+1) = ztrans * za_i(ji,jl2) 
     756            zh_i(ji,jl2+1) = hi_mean(jl2+1) 
     757            ! remove from cat jl2 
     758            za_i(ji,jl2  ) = ( 1._wp - ztrans ) * za_i(ji,jl2) 
     759         ENDIF 
     760      END DO 
     761      ! 
     762      jlfil2(:,:) = jlfil(:,:)  
     763      ! fill categories from low to high 
     764      DO jl = 2, jpl-1 
     765         DO ji = 1, idim 
     766            IF( jlfil(ji,jl-1) /= 0 .AND. jlfil(ji,jl) == 0 ) THEN 
     767               ! fill high 
     768               za_i(ji,jl) = ztrans * za_i(ji,jl-1) 
     769               zh_i(ji,jl) = hi_mean(jl) 
     770               jlfil(ji,jl) = jl 
     771               ! remove low 
     772               za_i(ji,jl-1) = ( 1._wp - ztrans ) * za_i(ji,jl-1) 
     773            ENDIF 
     774         END DO 
     775      END DO 
     776      ! 
     777      ! fill categories from high to low 
     778      DO jl = jpl-1, 2, -1 
     779         DO ji = 1, idim 
     780            IF( jlfil2(ji,jl+1) /= 0 .AND. jlfil2(ji,jl) == 0 ) THEN 
     781               ! fill low 
     782               za_i(ji,jl) = za_i(ji,jl) + ztrans * za_i(ji,jl+1) 
     783               zh_i(ji,jl) = hi_mean(jl)  
     784               jlfil2(ji,jl) = jl 
     785               ! remove high 
     786               za_i(ji,jl+1) = ( 1._wp - ztrans ) * za_i(ji,jl+1) 
     787            ENDIF 
     788         END DO 
     789      END DO 
     790      ! 
     791      DEALLOCATE( jlfil, jlfil2 )      ! deallocate arrays 
     792      DEALLOCATE( jlmin, jlmax ) 
     793      ! 
     794   END SUBROUTINE ice_var_itd2 
     795 
     796 
     797   SUBROUTINE ice_var_bv 
    667798      !!------------------------------------------------------------------- 
    668799      !!                ***  ROUTINE ice_var_bv *** 
  • branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90

    r8586 r8813  
    5151 
    5252#if defined key_lim3 
    53    LOGICAL :: ll_bdylim3                  ! determine whether ice input is 1cat (F) or Xcat (T) type 
    54    INTEGER :: jfld_hti, jfld_hts, jfld_ai ! indices of ice thickness, snow thickness and concentration in bf structure 
     53   INTEGER ::   nice_cat                      ! number of categories in the input file 
     54   INTEGER ::   jfld_hti, jfld_hts, jfld_ai  ! indices of ice thickness, snow thickness and concentration in bf structure 
    5555#endif 
    5656 
    5757   !!---------------------------------------------------------------------- 
    58    !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
     58   !! NEMO/OPA 4.0 , NEMO Consortium (2017) 
    5959   !! $Id$  
    6060   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
     
    8080      !                                               ! etc. 
    8181      ! 
    82       INTEGER ::  ib_bdy, jfld, jstart, jend, ib, ii, ij, ik, igrd, jl  ! local indices 
     82      INTEGER ::  jbdy, jfld, jstart, jend, ib, jl  ! dummy loop indices 
     83      INTEGER ::  ii, ij, ik, igrd                  ! local integers 
    8384      INTEGER,          DIMENSION(jpbgrd) ::   ilen1  
    8485      INTEGER, POINTER, DIMENSION(:)      ::   nblen, nblenrim  ! short cuts 
     
    9596         !----------------------------- 
    9697          
    97          DO ib_bdy = 1, nb_bdy 
     98         DO jbdy = 1, nb_bdy 
    9899            ! 
    99             nblen => idx_bdy(ib_bdy)%nblen 
    100             nblenrim => idx_bdy(ib_bdy)%nblenrim 
    101             dta => dta_bdy(ib_bdy) 
    102  
    103             IF( nn_dyn2d_dta(ib_bdy) == 0 ) THEN  
     100            nblen    => idx_bdy(jbdy)%nblen 
     101            nblenrim => idx_bdy(jbdy)%nblenrim 
     102            dta      => dta_bdy(jbdy) 
     103            ! 
     104            IF( nn_dyn2d_dta(jbdy) == 0 ) THEN  
    104105               ilen1(:) = nblen(:) 
    105106               IF( dta%ll_ssh ) THEN  
    106107                  igrd = 1 
    107108                  DO ib = 1, ilen1(igrd) 
    108                      ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    109                      ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
    110                      dta_bdy(ib_bdy)%ssh(ib) = sshn(ii,ij) * tmask(ii,ij,1)          
     109                     ii = idx_bdy(jbdy)%nbi(ib,igrd) 
     110                     ij = idx_bdy(jbdy)%nbj(ib,igrd) 
     111                     dta_bdy(jbdy)%ssh(ib) = sshn(ii,ij) * tmask(ii,ij,1)          
    111112                  END DO  
    112                END IF 
     113               ENDIF 
    113114               IF( dta%ll_u2d ) THEN  
    114115                  igrd = 2 
    115116                  DO ib = 1, ilen1(igrd) 
    116                      ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    117                      ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
    118                      dta_bdy(ib_bdy)%u2d(ib) = un_b(ii,ij) * umask(ii,ij,1)          
     117                     ii = idx_bdy(jbdy)%nbi(ib,igrd) 
     118                     ij = idx_bdy(jbdy)%nbj(ib,igrd) 
     119                     dta_bdy(jbdy)%u2d(ib) = un_b(ii,ij) * umask(ii,ij,1)          
    119120                  END DO  
    120                END IF 
     121               ENDIF 
    121122               IF( dta%ll_v2d ) THEN  
    122123                  igrd = 3 
    123124                  DO ib = 1, ilen1(igrd) 
    124                      ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    125                      ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
    126                      dta_bdy(ib_bdy)%v2d(ib) = vn_b(ii,ij) * vmask(ii,ij,1)          
     125                     ii = idx_bdy(jbdy)%nbi(ib,igrd) 
     126                     ij = idx_bdy(jbdy)%nbj(ib,igrd) 
     127                     dta_bdy(jbdy)%v2d(ib) = vn_b(ii,ij) * vmask(ii,ij,1)          
    127128                  END DO  
    128                END IF 
    129             ENDIF 
    130  
    131             IF( nn_dyn3d_dta(ib_bdy) == 0 ) THEN  
     129               ENDIF 
     130            ENDIF 
     131            ! 
     132            IF( nn_dyn3d_dta(jbdy) == 0 ) THEN  
    132133               ilen1(:) = nblen(:) 
    133134               IF( dta%ll_u3d ) THEN  
     
    135136                  DO ib = 1, ilen1(igrd) 
    136137                     DO ik = 1, jpkm1 
    137                         ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    138                         ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
    139                         dta_bdy(ib_bdy)%u3d(ib,ik) =  ( un(ii,ij,ik) - un_b(ii,ij) ) * umask(ii,ij,ik)          
     138                        ii = idx_bdy(jbdy)%nbi(ib,igrd) 
     139                        ij = idx_bdy(jbdy)%nbj(ib,igrd) 
     140                        dta_bdy(jbdy)%u3d(ib,ik) =  ( un(ii,ij,ik) - un_b(ii,ij) ) * umask(ii,ij,ik)          
    140141                     END DO 
    141142                  END DO  
    142                END IF 
     143               ENDIF 
    143144               IF( dta%ll_v3d ) THEN  
    144145                  igrd = 3  
    145146                  DO ib = 1, ilen1(igrd) 
    146147                     DO ik = 1, jpkm1 
    147                         ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    148                         ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
    149                         dta_bdy(ib_bdy)%v3d(ib,ik) =  ( vn(ii,ij,ik) - vn_b(ii,ij) ) * vmask(ii,ij,ik)          
     148                        ii = idx_bdy(jbdy)%nbi(ib,igrd) 
     149                        ij = idx_bdy(jbdy)%nbj(ib,igrd) 
     150                        dta_bdy(jbdy)%v3d(ib,ik) =  ( vn(ii,ij,ik) - vn_b(ii,ij) ) * vmask(ii,ij,ik)          
    150151                        END DO 
    151152                  END DO  
    152                END IF 
    153             ENDIF 
    154  
    155             IF( nn_tra_dta(ib_bdy) == 0 ) THEN  
     153               ENDIF 
     154            ENDIF 
     155 
     156            IF( nn_tra_dta(jbdy) == 0 ) THEN  
    156157               ilen1(:) = nblen(:) 
    157158               IF( dta%ll_tem ) THEN 
     
    159160                  DO ib = 1, ilen1(igrd) 
    160161                     DO ik = 1, jpkm1 
    161                         ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    162                         ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
    163                         dta_bdy(ib_bdy)%tem(ib,ik) = tsn(ii,ij,ik,jp_tem) * tmask(ii,ij,ik)          
     162                        ii = idx_bdy(jbdy)%nbi(ib,igrd) 
     163                        ij = idx_bdy(jbdy)%nbj(ib,igrd) 
     164                        dta_bdy(jbdy)%tem(ib,ik) = tsn(ii,ij,ik,jp_tem) * tmask(ii,ij,ik)          
    164165                     END DO 
    165166                  END DO  
    166                END IF 
     167               ENDIF 
    167168               IF( dta%ll_sal ) THEN 
    168169                  igrd = 1  
    169170                  DO ib = 1, ilen1(igrd) 
    170171                     DO ik = 1, jpkm1 
    171                         ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    172                         ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
    173                         dta_bdy(ib_bdy)%sal(ib,ik) = tsn(ii,ij,ik,jp_sal) * tmask(ii,ij,ik)          
     172                        ii = idx_bdy(jbdy)%nbi(ib,igrd) 
     173                        ij = idx_bdy(jbdy)%nbj(ib,igrd) 
     174                        dta_bdy(jbdy)%sal(ib,ik) = tsn(ii,ij,ik,jp_sal) * tmask(ii,ij,ik)          
    174175                     END DO 
    175176                  END DO  
    176                END IF 
    177             ENDIF 
    178  
    179 #if defined key_lim3 
    180             IF( nn_ice_lim_dta(ib_bdy) == 0 ) THEN  
     177               ENDIF 
     178            ENDIF 
     179 
     180#if defined key_lim3 
     181            IF( nn_ice_lim_dta(jbdy) == 0 ) THEN    ! set ice to initial values 
    181182               ilen1(:) = nblen(:) 
    182183               IF( dta%ll_a_i ) THEN 
     
    184185                  DO jl = 1, jpl 
    185186                     DO ib = 1, ilen1(igrd) 
    186                         ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    187                         ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
    188                         dta_bdy(ib_bdy)%a_i (ib,jl) =  a_i(ii,ij,jl) * tmask(ii,ij,1)  
     187                        ii = idx_bdy(jbdy)%nbi(ib,igrd) 
     188                        ij = idx_bdy(jbdy)%nbj(ib,igrd) 
     189                        dta_bdy(jbdy)%a_i (ib,jl) =  a_i(ii,ij,jl) * tmask(ii,ij,1)  
    189190                     END DO 
    190191                  END DO 
     
    194195                  DO jl = 1, jpl 
    195196                     DO ib = 1, ilen1(igrd) 
    196                         ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    197                         ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
    198                         dta_bdy(ib_bdy)%h_i (ib,jl) =  h_i(ii,ij,jl) * tmask(ii,ij,1)  
     197                        ii = idx_bdy(jbdy)%nbi(ib,igrd) 
     198                        ij = idx_bdy(jbdy)%nbj(ib,igrd) 
     199                        dta_bdy(jbdy)%h_i (ib,jl) =  h_i(ii,ij,jl) * tmask(ii,ij,1)  
    199200                     END DO 
    200201                  END DO 
     
    204205                  DO jl = 1, jpl 
    205206                     DO ib = 1, ilen1(igrd) 
    206                         ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    207                         ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
    208                         dta_bdy(ib_bdy)%h_s (ib,jl) =  h_s(ii,ij,jl) * tmask(ii,ij,1)  
     207                        ii = idx_bdy(jbdy)%nbi(ib,igrd) 
     208                        ij = idx_bdy(jbdy)%nbj(ib,igrd) 
     209                        dta_bdy(jbdy)%h_s (ib,jl) =  h_s(ii,ij,jl) * tmask(ii,ij,1)  
    209210                     END DO 
    210211                  END DO 
     
    212213            ENDIF 
    213214#endif 
    214          END DO ! ib_bdy 
     215         END DO ! jbdy 
    215216         ! 
    216217      ENDIF ! kt == nit000 
     
    220221      
    221222      jstart = 1 
    222       DO ib_bdy = 1, nb_bdy    
    223          dta => dta_bdy(ib_bdy) 
    224          IF( nn_dta(ib_bdy) == 1 ) THEN ! skip this bit if no external data required 
     223      DO jbdy = 1, nb_bdy    
     224         dta => dta_bdy(jbdy) 
     225         IF( nn_dta(jbdy) == 1 ) THEN ! skip this bit if no external data required 
    225226       
    226227            IF( PRESENT(jit) ) THEN 
    227228               ! Update barotropic boundary conditions only 
    228229               ! jit is optional argument for fld_read and bdytide_update 
    229                IF( cn_dyn2d(ib_bdy) /= 'none' ) THEN 
    230                   IF( nn_dyn2d_dta(ib_bdy) == 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays 
     230               IF( cn_dyn2d(jbdy) /= 'none' ) THEN 
     231                  IF( nn_dyn2d_dta(jbdy) == 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays 
    231232                     IF( dta%ll_ssh ) dta%ssh(:) = 0._wp 
    232233                     IF( dta%ll_u2d ) dta%u2d(:) = 0._wp 
    233234                     IF( dta%ll_u3d ) dta%v2d(:) = 0._wp 
    234235                  ENDIF 
    235                   IF (cn_tra(ib_bdy) /= 'runoff') THEN 
    236                      IF( nn_dyn2d_dta(ib_bdy) == 1 .OR. nn_dyn2d_dta(ib_bdy) == 3 ) THEN 
     236                  IF (cn_tra(jbdy) /= 'runoff') THEN 
     237                     IF( nn_dyn2d_dta(jbdy) == 1 .OR. nn_dyn2d_dta(jbdy) == 3 ) THEN 
    237238 
    238239                        jend = jstart + dta%nread(2) - 1 
    239                         IF( ln_full_vel_array(ib_bdy) ) THEN 
     240                        IF( ln_full_vel_array(jbdy) ) THEN 
    240241                           CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend),  & 
    241                                      & kit=jit, kt_offset=time_offset , jpk_bdy=nb_jpk_bdy, fvl=ln_full_vel_array(ib_bdy)  ) 
     242                                     & kit=jit, kt_offset=time_offset , jpk_bdy=nb_jpk_bdy, fvl=ln_full_vel_array(jbdy)  ) 
    242243                        ELSE 
    243244                           CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend),  & 
     
    246247 
    247248                        ! If full velocities in boundary data then extract barotropic velocities from 3D fields 
    248                         IF( ln_full_vel_array(ib_bdy) .AND.                                             & 
    249                           &    ( nn_dyn2d_dta(ib_bdy) == 1 .OR. nn_dyn2d_dta(ib_bdy) == 3 .OR.  & 
    250                           &      nn_dyn3d_dta(ib_bdy) == 1 ) )THEN 
     249                        IF( ln_full_vel_array(jbdy) .AND.                                             & 
     250                          &    ( nn_dyn2d_dta(jbdy) == 1 .OR. nn_dyn2d_dta(jbdy) == 3 .OR.  & 
     251                          &      nn_dyn3d_dta(jbdy) == 1 ) )THEN 
    251252 
    252253                           igrd = 2                      ! zonal velocity 
    253254                           dta%u2d(:) = 0._wp 
    254                            DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 
    255                               ii   = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    256                               ij   = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     255                           DO ib = 1, idx_bdy(jbdy)%nblen(igrd) 
     256                              ii   = idx_bdy(jbdy)%nbi(ib,igrd) 
     257                              ij   = idx_bdy(jbdy)%nbj(ib,igrd) 
    257258                              DO ik = 1, jpkm1 
    258259                                 dta%u2d(ib) = dta%u2d(ib) & 
     
    263264                           igrd = 3                      ! meridional velocity 
    264265                           dta%v2d(:) = 0._wp 
    265                            DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 
    266                               ii   = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    267                               ij   = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     266                           DO ib = 1, idx_bdy(jbdy)%nblen(igrd) 
     267                              ii   = idx_bdy(jbdy)%nbi(ib,igrd) 
     268                              ij   = idx_bdy(jbdy)%nbj(ib,igrd) 
    268269                              DO ik = 1, jpkm1 
    269270                                 dta%v2d(ib) = dta%v2d(ib) & 
     
    274275                        ENDIF                     
    275276                     ENDIF 
    276                      IF( nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN ! update tidal harmonic forcing 
    277                         CALL bdytide_update( kt=kt, idx=idx_bdy(ib_bdy), dta=dta, td=tides(ib_bdy),   &  
     277                     IF( nn_dyn2d_dta(jbdy) .ge. 2 ) THEN ! update tidal harmonic forcing 
     278                        CALL bdytide_update( kt=kt, idx=idx_bdy(jbdy), dta=dta, td=tides(jbdy),   &  
    278279                          &                 jit=jit, time_offset=time_offset ) 
    279280                     ENDIF 
     
    281282               ENDIF 
    282283            ELSE 
    283                IF (cn_tra(ib_bdy) == 'runoff') then      ! runoff condition 
    284                   jend = nb_bdy_fld(ib_bdy) 
     284               IF (cn_tra(jbdy) == 'runoff') then      ! runoff condition 
     285                  jend = nb_bdy_fld(jbdy) 
    285286                  CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend),  & 
    286287                               & map=nbmap_ptr(jstart:jend), kt_offset=time_offset ) 
    287288                  ! 
    288289                  igrd = 2                      ! zonal velocity 
    289                   DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 
    290                      ii   = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    291                      ij   = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     290                  DO ib = 1, idx_bdy(jbdy)%nblen(igrd) 
     291                     ii   = idx_bdy(jbdy)%nbi(ib,igrd) 
     292                     ij   = idx_bdy(jbdy)%nbj(ib,igrd) 
    292293                     dta%u2d(ib) = dta%u2d(ib) / ( e2u(ii,ij) * hu_0(ii,ij) ) 
    293294                  END DO 
    294295                  ! 
    295296                  igrd = 3                      ! meridional velocity 
    296                   DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 
    297                      ii   = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    298                      ij   = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     297                  DO ib = 1, idx_bdy(jbdy)%nblen(igrd) 
     298                     ii   = idx_bdy(jbdy)%nbi(ib,igrd) 
     299                     ij   = idx_bdy(jbdy)%nbj(ib,igrd) 
    299300                     dta%v2d(ib) = dta%v2d(ib) / ( e1v(ii,ij) * hv_0(ii,ij) ) 
    300301                  END DO 
    301302               ELSE 
    302                   IF( nn_dyn2d_dta(ib_bdy) == 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays 
     303                  IF( nn_dyn2d_dta(jbdy) == 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays 
    303304                     IF( dta%ll_ssh ) dta%ssh(:) = 0._wp 
    304305                     IF( dta%ll_u2d ) dta%u2d(:) = 0._wp 
     
    308309                     jend = jstart + dta%nread(1) - 1 
    309310                     CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), & 
    310                                   & map=nbmap_ptr(jstart:jend), kt_offset=time_offset, jpk_bdy=nb_jpk_bdy, fvl=ln_full_vel_array(ib_bdy) ) 
     311                                  & map=nbmap_ptr(jstart:jend), kt_offset=time_offset, jpk_bdy=nb_jpk_bdy, fvl=ln_full_vel_array(jbdy) ) 
    311312                  ENDIF 
    312313                  ! If full velocities in boundary data then split into barotropic and baroclinic data 
    313                   IF( ln_full_vel_array(ib_bdy) .and.                                             & 
    314                     & ( nn_dyn2d_dta(ib_bdy) == 1 .OR. nn_dyn2d_dta(ib_bdy) == 3 .OR. & 
    315                     &   nn_dyn3d_dta(ib_bdy) == 1 ) ) THEN 
     314                  IF( ln_full_vel_array(jbdy) .and.                                             & 
     315                    & ( nn_dyn2d_dta(jbdy) == 1 .OR. nn_dyn2d_dta(jbdy) == 3 .OR. & 
     316                    &   nn_dyn3d_dta(jbdy) == 1 ) ) THEN 
    316317                     igrd = 2                      ! zonal velocity 
    317318                     dta%u2d(:) = 0._wp 
    318                      DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 
    319                         ii   = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    320                         ij   = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     319                     DO ib = 1, idx_bdy(jbdy)%nblen(igrd) 
     320                        ii   = idx_bdy(jbdy)%nbi(ib,igrd) 
     321                        ij   = idx_bdy(jbdy)%nbj(ib,igrd) 
    321322                        DO ik = 1, jpkm1 
    322323                           dta%u2d(ib) = dta%u2d(ib) & 
     
    330331                     igrd = 3                      ! meridional velocity 
    331332                     dta%v2d(:) = 0._wp 
    332                      DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 
    333                         ii   = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    334                         ij   = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     333                     DO ib = 1, idx_bdy(jbdy)%nblen(igrd) 
     334                        ii   = idx_bdy(jbdy)%nbi(ib,igrd) 
     335                        ij   = idx_bdy(jbdy)%nbj(ib,igrd) 
    335336                        DO ik = 1, jpkm1 
    336337                           dta%v2d(ib) = dta%v2d(ib) & 
     
    346347               ENDIF 
    347348#if defined key_lim3 
    348                IF( .NOT. ll_bdylim3 .AND. cn_ice_lim(ib_bdy) /= 'none' .AND. nn_ice_lim_dta(ib_bdy) == 1 ) THEN ! bdy ice input (case input is 1cat) 
    349                 CALL ice_var_itd ( bf(jfld_hti)%fnow(:,1,1), bf(jfld_hts)%fnow(:,1,1), bf(jfld_ai)%fnow(:,1,1), & 
    350                                   & dta_bdy(ib_bdy)%h_i,     dta_bdy(ib_bdy)%h_s,     dta_bdy(ib_bdy)%a_i     ) 
     349               IF( cn_ice_lim(jbdy) /= 'none' .AND. nn_ice_lim_dta(jbdy) == 1 ) THEN 
     350                  IF( nice_cat == 1 ) THEN ! case input cat = 1 
     351                     CALL ice_var_itd ( bf(jfld_hti)%fnow(:,1,1), bf(jfld_hts)%fnow(:,1,1), bf(jfld_ai)%fnow(:,1,1), & 
     352                        &               dta_bdy(jbdy)%h_i     , dta_bdy(jbdy)%h_s     , dta_bdy(jbdy)%a_i    ) 
     353                  ELSEIF( nice_cat /= 1 .AND. nice_cat /= jpl ) THEN ! case input cat /=1 and /=jpl 
     354                     CALL ice_var_itd2( bf(jfld_hti)%fnow(:,1,:), bf(jfld_hts)%fnow(:,1,:), bf(jfld_ai)%fnow(:,1,:), & 
     355                        &               dta_bdy(jbdy)%h_i     , dta_bdy(jbdy)%h_s     , dta_bdy(jbdy)%a_i    ) 
     356                  ENDIF 
    351357               ENDIF 
    352358#endif 
    353359            ENDIF 
    354360            jstart = jstart + dta%nread(1) 
    355          END IF ! nn_dta(ib_bdy) = 1 
    356       END DO  ! ib_bdy 
     361         ENDIF    ! nn_dta(jbdy) = 1 
     362      END DO  ! jbdy 
    357363 
    358364      IF ( ln_tide ) THEN 
    359365         IF (ln_dynspg_ts) THEN      ! Fill temporary arrays with slow-varying bdy data                            
    360             DO ib_bdy = 1, nb_bdy    ! Tidal component added in ts loop 
    361                IF ( nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN 
    362                   nblen => idx_bdy(ib_bdy)%nblen 
    363                   nblenrim => idx_bdy(ib_bdy)%nblenrim 
    364                   IF( cn_dyn2d(ib_bdy) == 'frs' ) THEN; ilen1(:)=nblen(:) ; ELSE ; ilen1(:)=nblenrim(:) ; ENDIF  
    365                   IF ( dta_bdy(ib_bdy)%ll_ssh ) dta_bdy_s(ib_bdy)%ssh(1:ilen1(1)) = dta_bdy(ib_bdy)%ssh(1:ilen1(1)) 
    366                   IF ( dta_bdy(ib_bdy)%ll_u2d ) dta_bdy_s(ib_bdy)%u2d(1:ilen1(2)) = dta_bdy(ib_bdy)%u2d(1:ilen1(2)) 
    367                   IF ( dta_bdy(ib_bdy)%ll_v2d ) dta_bdy_s(ib_bdy)%v2d(1:ilen1(3)) = dta_bdy(ib_bdy)%v2d(1:ilen1(3)) 
     366            DO jbdy = 1, nb_bdy    ! Tidal component added in ts loop 
     367               IF ( nn_dyn2d_dta(jbdy) .ge. 2 ) THEN 
     368                  nblen => idx_bdy(jbdy)%nblen 
     369                  nblenrim => idx_bdy(jbdy)%nblenrim 
     370                  IF( cn_dyn2d(jbdy) == 'frs' ) THEN; ilen1(:)=nblen(:) ; ELSE ; ilen1(:)=nblenrim(:) ; ENDIF  
     371                  IF ( dta_bdy(jbdy)%ll_ssh ) dta_bdy_s(jbdy)%ssh(1:ilen1(1)) = dta_bdy(jbdy)%ssh(1:ilen1(1)) 
     372                  IF ( dta_bdy(jbdy)%ll_u2d ) dta_bdy_s(jbdy)%u2d(1:ilen1(2)) = dta_bdy(jbdy)%u2d(1:ilen1(2)) 
     373                  IF ( dta_bdy(jbdy)%ll_v2d ) dta_bdy_s(jbdy)%v2d(1:ilen1(3)) = dta_bdy(jbdy)%v2d(1:ilen1(3)) 
    368374               ENDIF 
    369375            END DO 
     
    375381 
    376382      IF ( ln_apr_obc ) THEN 
    377          DO ib_bdy = 1, nb_bdy 
    378             IF (cn_tra(ib_bdy) /= 'runoff')THEN 
     383         DO jbdy = 1, nb_bdy 
     384            IF (cn_tra(jbdy) /= 'runoff')THEN 
    379385               igrd = 1                      ! meridional velocity 
    380                DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
    381                   ii   = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    382                   ij   = idx_bdy(ib_bdy)%nbj(ib,igrd) 
    383                   dta_bdy(ib_bdy)%ssh(ib) = dta_bdy(ib_bdy)%ssh(ib) + ssh_ib(ii,ij) 
     386               DO ib = 1, idx_bdy(jbdy)%nblenrim(igrd) 
     387                  ii   = idx_bdy(jbdy)%nbi(ib,igrd) 
     388                  ij   = idx_bdy(jbdy)%nbj(ib,igrd) 
     389                  dta_bdy(jbdy)%ssh(ib) = dta_bdy(jbdy)%ssh(ib) + ssh_ib(ii,ij) 
    384390               END DO 
    385391            ENDIF 
     
    402408      !!                 
    403409      !!---------------------------------------------------------------------- 
    404       INTEGER ::   ib_bdy, jfld, jstart, jend, ierror, ios     ! Local integers 
     410      INTEGER ::   jbdy, jfld, jstart, jend, ierror, ios     ! Local integers 
    405411      ! 
    406412      CHARACTER(len=100)                     ::   cn_dir        ! Root directory for location of data files 
     
    408414      CHARACTER(len = 256)::   clname                           ! temporary file name 
    409415      LOGICAL                                ::   ln_full_vel   ! =T => full velocities in 3D boundary data 
    410                                                                 ! =F => baroclinic velocities in 3D boundary data 
     416      !                                                         ! =F => baroclinic velocities in 3D boundary data 
    411417      INTEGER                                ::   ilen_global   ! Max length required for global bdy dta arrays 
    412418      INTEGER, ALLOCATABLE, DIMENSION(:)     ::   ilen1, ilen3  ! size of 1st and 3rd dimensions of local arrays 
     
    416422      TYPE(OBC_DATA), POINTER                ::   dta           ! short cut 
    417423#if defined key_lim3 
    418       INTEGER               ::   zndims   ! number of dimensions in an array (i.e. 3 = wo ice cat; 4 = w ice cat) 
     424      INTEGER               ::   kndims   ! number of dimensions in an array (i.e. 3 = wo ice cat; 4 = w ice cat) 
     425      INTEGER, DIMENSION(4) ::   kdimsz   ! size   of dimensions 
    419426      INTEGER               ::   inum,id1 ! local integer 
    420427#endif 
     
    440447 
    441448      ! Set nn_dta 
    442       DO ib_bdy = 1, nb_bdy 
    443          nn_dta(ib_bdy) = MAX(  nn_dyn2d_dta(ib_bdy)       & 
    444                                ,nn_dyn3d_dta(ib_bdy)       & 
    445                                ,nn_tra_dta(ib_bdy)         & 
    446 #if defined key_lim3 
    447                               ,nn_ice_lim_dta(ib_bdy)    & 
     449      DO jbdy = 1, nb_bdy 
     450         nn_dta(jbdy) = MAX(   nn_dyn2d_dta  (jbdy)    & 
     451            &                , nn_dyn3d_dta  (jbdy)    & 
     452            &                , nn_tra_dta    (jbdy)    & 
     453#if defined key_lim3 
     454            &                , nn_ice_lim_dta(jbdy)    & 
    448455#endif 
    449456                              ) 
    450          IF(nn_dta(ib_bdy) > 1) nn_dta(ib_bdy) = 1 
     457         IF(nn_dta(jbdy) > 1)   nn_dta(jbdy) = 1 
    451458      END DO 
    452459 
     
    455462      ALLOCATE( nb_bdy_fld(nb_bdy) ) 
    456463      nb_bdy_fld(:) = 0 
    457       DO ib_bdy = 1, nb_bdy          
    458          IF( cn_dyn2d(ib_bdy) /= 'none' .and. ( nn_dyn2d_dta(ib_bdy) == 1 .or. nn_dyn2d_dta(ib_bdy) == 3 ) ) THEN 
    459             nb_bdy_fld(ib_bdy) = nb_bdy_fld(ib_bdy) + 3 
    460          ENDIF 
    461          IF( cn_dyn3d(ib_bdy) /= 'none' .and. nn_dyn3d_dta(ib_bdy) == 1 ) THEN 
    462             nb_bdy_fld(ib_bdy) = nb_bdy_fld(ib_bdy) + 2 
    463          ENDIF 
    464          IF( cn_tra(ib_bdy) /= 'none' .and. nn_tra_dta(ib_bdy) == 1  ) THEN 
    465             nb_bdy_fld(ib_bdy) = nb_bdy_fld(ib_bdy) + 2 
    466          ENDIF 
    467 #if defined key_lim3 
    468          IF( cn_ice_lim(ib_bdy) /= 'none' .and. nn_ice_lim_dta(ib_bdy) == 1  ) THEN 
    469             nb_bdy_fld(ib_bdy) = nb_bdy_fld(ib_bdy) + 3 
     464      DO jbdy = 1, nb_bdy          
     465         IF( cn_dyn2d(jbdy) /= 'none' .AND. ( nn_dyn2d_dta(jbdy) == 1 .or. nn_dyn2d_dta(jbdy) == 3 ) ) THEN 
     466            nb_bdy_fld(jbdy) = nb_bdy_fld(jbdy) + 3 
     467         ENDIF 
     468         IF( cn_dyn3d(jbdy) /= 'none' .AND. nn_dyn3d_dta(jbdy) == 1 ) THEN 
     469            nb_bdy_fld(jbdy) = nb_bdy_fld(jbdy) + 2 
     470         ENDIF 
     471         IF( cn_tra(jbdy) /= 'none' .AND. nn_tra_dta(jbdy) == 1  ) THEN 
     472            nb_bdy_fld(jbdy) = nb_bdy_fld(jbdy) + 2 
     473         ENDIF 
     474#if defined key_lim3 
     475         IF( cn_ice_lim(jbdy) /= 'none' .AND. nn_ice_lim_dta(jbdy) == 1  ) THEN 
     476            nb_bdy_fld(jbdy) = nb_bdy_fld(jbdy) + 3 
    470477         ENDIF 
    471478#endif                
    472          IF(lwp) WRITE(numout,*) 'Maximum number of files to open =',nb_bdy_fld(ib_bdy) 
     479         IF(lwp) WRITE(numout,*) 'Maximum number of files to open =', nb_bdy_fld(jbdy) 
    473480      END DO             
    474481 
     
    496503      REWIND(numnam_cfg) 
    497504      jfld = 0  
    498       DO ib_bdy = 1, nb_bdy          
    499          IF( nn_dta(ib_bdy) == 1 ) THEN 
     505      DO jbdy = 1, nb_bdy          
     506         IF( nn_dta(jbdy) == 1 ) THEN 
    500507            READ  ( numnam_ref, nambdy_dta, IOSTAT = ios, ERR = 901) 
    501508901         IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy_dta in reference namelist', lwp ) 
     
    505512            IF(lwm) WRITE( numond, nambdy_dta ) 
    506513 
    507             cn_dir_array(ib_bdy) = cn_dir 
    508             ln_full_vel_array(ib_bdy) = ln_full_vel 
    509  
    510             nblen => idx_bdy(ib_bdy)%nblen 
    511             nblenrim => idx_bdy(ib_bdy)%nblenrim 
    512             dta => dta_bdy(ib_bdy) 
     514            cn_dir_array(jbdy) = cn_dir 
     515            ln_full_vel_array(jbdy) = ln_full_vel 
     516 
     517            nblen => idx_bdy(jbdy)%nblen 
     518            nblenrim => idx_bdy(jbdy)%nblenrim 
     519            dta => dta_bdy(jbdy) 
    513520            dta%nread(2) = 0 
    514521 
    515522            ! Only read in necessary fields for this set. 
    516523            ! Important that barotropic variables come first. 
    517             IF( nn_dyn2d_dta(ib_bdy) == 1 .or. nn_dyn2d_dta(ib_bdy) == 3 ) THEN  
     524            IF( nn_dyn2d_dta(jbdy) == 1 .or. nn_dyn2d_dta(jbdy) == 3 ) THEN  
    518525 
    519526               IF( dta%ll_ssh ) THEN  
     
    521528                  jfld = jfld + 1 
    522529                  blf_i(jfld) = bn_ssh 
    523                   ibdy(jfld) = ib_bdy 
     530                  ibdy(jfld) = jbdy 
    524531                  igrid(jfld) = 1 
    525532                  ilen1(jfld) = nblen(igrid(jfld)) 
     
    528535               ENDIF 
    529536 
    530                IF( dta%ll_u2d .and. .not. ln_full_vel_array(ib_bdy) ) THEN 
     537               IF( dta%ll_u2d .and. .not. ln_full_vel_array(jbdy) ) THEN 
    531538                  if(lwp) write(numout,*) '++++++ reading in u2d field' 
    532539                  jfld = jfld + 1 
    533540                  blf_i(jfld) = bn_u2d 
    534                   ibdy(jfld) = ib_bdy 
     541                  ibdy(jfld) = jbdy 
    535542                  igrid(jfld) = 2 
    536543                  ilen1(jfld) = nblen(igrid(jfld)) 
     
    539546               ENDIF 
    540547 
    541                IF( dta%ll_v2d .and. .not. ln_full_vel_array(ib_bdy) ) THEN 
     548               IF( dta%ll_v2d .and. .not. ln_full_vel_array(jbdy) ) THEN 
    542549                  if(lwp) write(numout,*) '++++++ reading in v2d field' 
    543550                  jfld = jfld + 1 
    544551                  blf_i(jfld) = bn_v2d 
    545                   ibdy(jfld) = ib_bdy 
     552                  ibdy(jfld) = jbdy 
    546553                  igrid(jfld) = 3 
    547554                  ilen1(jfld) = nblen(igrid(jfld)) 
     
    554561            ! read 3D velocities if baroclinic velocities require OR if 
    555562            ! barotropic velocities required and ln_full_vel set to .true. 
    556             IF( nn_dyn3d_dta(ib_bdy) == 1 .OR. & 
    557            &  ( ln_full_vel_array(ib_bdy) .AND. ( nn_dyn2d_dta(ib_bdy) == 1 .or. nn_dyn2d_dta(ib_bdy) == 3 ) ) ) THEN 
    558  
    559                IF( dta%ll_u3d .OR. ( ln_full_vel_array(ib_bdy) .and. dta%ll_u2d ) ) THEN 
     563            IF( nn_dyn3d_dta(jbdy) == 1 .OR. & 
     564           &  ( ln_full_vel_array(jbdy) .AND. ( nn_dyn2d_dta(jbdy) == 1 .OR. nn_dyn2d_dta(jbdy) == 3 ) ) ) THEN 
     565 
     566               IF( dta%ll_u3d .OR. ( ln_full_vel_array(jbdy) .and. dta%ll_u2d ) ) THEN 
    560567                  if(lwp) write(numout,*) '++++++ reading in u3d field' 
    561568                  jfld = jfld + 1 
    562569                  blf_i(jfld) = bn_u3d 
    563                   ibdy(jfld) = ib_bdy 
     570                  ibdy(jfld) = jbdy 
    564571                  igrid(jfld) = 2 
    565572                  ilen1(jfld) = nblen(igrid(jfld)) 
    566573                  ilen3(jfld) = jpk 
    567                   IF( ln_full_vel_array(ib_bdy) .and. dta%ll_u2d ) dta%nread(2) = dta%nread(2) + 1 
    568                ENDIF 
    569  
    570                IF( dta%ll_v3d .OR. ( ln_full_vel_array(ib_bdy) .and. dta%ll_v2d ) ) THEN 
     574                  IF( ln_full_vel_array(jbdy) .and. dta%ll_u2d ) dta%nread(2) = dta%nread(2) + 1 
     575               ENDIF 
     576 
     577               IF( dta%ll_v3d .OR. ( ln_full_vel_array(jbdy) .and. dta%ll_v2d ) ) THEN 
    571578                  if(lwp) write(numout,*) '++++++ reading in v3d field' 
    572579                  jfld = jfld + 1 
    573580                  blf_i(jfld) = bn_v3d 
    574                   ibdy(jfld) = ib_bdy 
     581                  ibdy(jfld) = jbdy 
    575582                  igrid(jfld) = 3 
    576583                  ilen1(jfld) = nblen(igrid(jfld)) 
    577584                  ilen3(jfld) = jpk 
    578                   IF( ln_full_vel_array(ib_bdy) .and. dta%ll_v2d ) dta%nread(2) = dta%nread(2) + 1 
     585                  IF( ln_full_vel_array(jbdy) .and. dta%ll_v2d ) dta%nread(2) = dta%nread(2) + 1 
    579586               ENDIF 
    580587 
     
    582589 
    583590            ! temperature and salinity 
    584             IF( nn_tra_dta(ib_bdy) == 1 ) THEN 
     591            IF( nn_tra_dta(jbdy) == 1 ) THEN 
    585592 
    586593               IF( dta%ll_tem ) THEN 
     
    588595                  jfld = jfld + 1 
    589596                  blf_i(jfld) = bn_tem 
    590                   ibdy(jfld) = ib_bdy 
     597                  ibdy(jfld) = jbdy 
    591598                  igrid(jfld) = 1 
    592599                  ilen1(jfld) = nblen(igrid(jfld)) 
     
    598605                  jfld = jfld + 1 
    599606                  blf_i(jfld) = bn_sal 
    600                   ibdy(jfld) = ib_bdy 
     607                  ibdy(jfld) = jbdy 
    601608                  igrid(jfld) = 1 
    602609                  ilen1(jfld) = nblen(igrid(jfld)) 
     
    608615#if defined key_lim3 
    609616            ! sea ice 
    610             IF( nn_ice_lim_dta(ib_bdy) == 1 ) THEN 
     617            IF( nn_ice_lim_dta(jbdy) == 1 ) THEN 
    611618               ! Test for types of ice input (1cat or Xcat)  
    612619               ! Build file name to find dimensions  
     
    622629               ! 
    623630               CALL iom_open  ( clname, inum ) 
    624                id1 = iom_varid( inum, bn_a_i%clvar, kndims=zndims, ldstop = .FALSE. ) 
     631               id1 = iom_varid( inum, bn_a_i%clvar, kdimsz=kdimsz, kndims=kndims, ldstop = .FALSE. ) 
    625632               CALL iom_close ( inum ) 
    626633 
    627                 IF ( zndims == 4 ) THEN 
    628                  ll_bdylim3 = .TRUE.   ! Xcat input 
     634                IF ( kndims == 4 ) THEN 
     635                 nice_cat = kdimsz(4)   ! Xcat input 
    629636               ELSE 
    630                  ll_bdylim3 = .FALSE.  ! 1cat input       
     637                 nice_cat = 1           ! 1cat input       
    631638               ENDIF 
    632639               ! End test 
     
    635642                  jfld = jfld + 1 
    636643                  blf_i(jfld) = bn_a_i 
    637                   ibdy(jfld) = ib_bdy 
     644                  ibdy(jfld)  = jbdy 
    638645                  igrid(jfld) = 1 
    639646                  ilen1(jfld) = nblen(igrid(jfld)) 
    640                   IF ( ll_bdylim3 ) THEN ; ilen3(jfld)=jpl ; ELSE ; ilen3(jfld)=1 ; ENDIF 
     647                  ilen3(jfld) = nice_cat 
    641648               ENDIF 
    642649 
     
    644651                  jfld = jfld + 1 
    645652                  blf_i(jfld) = bn_h_i 
    646                   ibdy(jfld) = ib_bdy 
     653                  ibdy(jfld)  = jbdy 
    647654                  igrid(jfld) = 1 
    648655                  ilen1(jfld) = nblen(igrid(jfld)) 
    649                   IF ( ll_bdylim3 ) THEN ; ilen3(jfld)=jpl ; ELSE ; ilen3(jfld)=1 ; ENDIF 
     656                  ilen3(jfld) = nice_cat 
    650657               ENDIF 
    651658 
    652659               IF( dta%ll_h_s ) THEN 
    653660                  jfld = jfld + 1 
    654                    blf_i(jfld) = bn_h_s 
    655                   ibdy(jfld) = ib_bdy 
     661                  blf_i(jfld) = bn_h_s 
     662                  ibdy(jfld)  = jbdy 
    656663                  igrid(jfld) = 1 
    657664                  ilen1(jfld) = nblen(igrid(jfld)) 
    658                   IF ( ll_bdylim3 ) THEN ; ilen3(jfld)=jpl ; ELSE ; ilen3(jfld)=1 ; ENDIF 
     665                  ilen3(jfld) = nice_cat 
    659666               ENDIF 
    660667 
     
    663670            ! Recalculate field counts 
    664671            !------------------------- 
    665             IF( ib_bdy == 1 ) THEN  
     672            IF( jbdy == 1 ) THEN  
    666673               nb_bdy_fld_sum = 0 
    667                nb_bdy_fld(ib_bdy) = jfld 
     674               nb_bdy_fld(jbdy) = jfld 
    668675               nb_bdy_fld_sum     = jfld               
    669676            ELSE 
    670                nb_bdy_fld(ib_bdy) = jfld - nb_bdy_fld_sum 
    671                nb_bdy_fld_sum = nb_bdy_fld_sum + nb_bdy_fld(ib_bdy) 
    672             ENDIF 
    673  
    674             dta%nread(1) = nb_bdy_fld(ib_bdy) 
     677               nb_bdy_fld(jbdy) = jfld - nb_bdy_fld_sum 
     678               nb_bdy_fld_sum = nb_bdy_fld_sum + nb_bdy_fld(jbdy) 
     679            ENDIF 
     680 
     681            dta%nread(1) = nb_bdy_fld(jbdy) 
    675682 
    676683         ENDIF ! nn_dta == 1 
    677       ENDDO ! ib_bdy 
     684      ENDDO ! jbdy 
    678685 
    679686      DO jfld = 1, nb_bdy_fld_sum 
     
    687694      !------------------------------------- 
    688695      jstart = 1 
    689       DO ib_bdy = 1, nb_bdy 
    690          jend = jstart - 1 + nb_bdy_fld(ib_bdy)  
    691          CALL fld_fill( bf(jstart:jend), blf_i(jstart:jend), cn_dir_array(ib_bdy), 'bdy_dta',   & 
     696      DO jbdy = 1, nb_bdy 
     697         jend = jstart - 1 + nb_bdy_fld(jbdy)  
     698         CALL fld_fill( bf(jstart:jend), blf_i(jstart:jend), cn_dir_array(jbdy), 'bdy_dta',   & 
    692699         &              'open boundary conditions', 'nambdy_dta' ) 
    693700         jstart = jend + 1 
     
    700707 
    701708      jfld = 0 
    702       DO ib_bdy=1, nb_bdy 
    703  
    704          nblen => idx_bdy(ib_bdy)%nblen 
    705          dta => dta_bdy(ib_bdy) 
     709      DO jbdy=1, nb_bdy 
     710 
     711         nblen => idx_bdy(jbdy)%nblen 
     712         dta => dta_bdy(jbdy) 
    706713 
    707714         if(lwp) then 
     
    715722         endif 
    716723 
    717          IF ( nn_dyn2d_dta(ib_bdy) == 0 .or. nn_dyn2d_dta(ib_bdy) == 2 ) THEN 
     724         IF ( nn_dyn2d_dta(jbdy) == 0 .or. nn_dyn2d_dta(jbdy) == 2 ) THEN 
    718725            if(lwp) write(numout,*) '++++++ dta%ssh/u2d/u3d allocated space' 
    719726            IF( dta%ll_ssh ) ALLOCATE( dta%ssh(nblen(1)) ) 
     
    721728            IF( dta%ll_v2d ) ALLOCATE( dta%v2d(nblen(3)) ) 
    722729         ENDIF 
    723          IF ( nn_dyn2d_dta(ib_bdy) == 1 .or. nn_dyn2d_dta(ib_bdy) == 3 ) THEN 
     730         IF ( nn_dyn2d_dta(jbdy) == 1 .or. nn_dyn2d_dta(jbdy) == 3 ) THEN 
    724731            IF( dta%ll_ssh ) THEN 
    725732               if(lwp) write(numout,*) '++++++ dta%ssh pointing to fnow' 
     
    728735            ENDIF 
    729736            IF ( dta%ll_u2d ) THEN 
    730                IF ( ln_full_vel_array(ib_bdy) ) THEN 
     737               IF ( ln_full_vel_array(jbdy) ) THEN 
    731738                  if(lwp) write(numout,*) '++++++ dta%u2d allocated space' 
    732739                  ALLOCATE( dta%u2d(nblen(2)) ) 
     
    738745            ENDIF 
    739746            IF ( dta%ll_v2d ) THEN 
    740                IF ( ln_full_vel_array(ib_bdy) ) THEN 
     747               IF ( ln_full_vel_array(jbdy) ) THEN 
    741748                  if(lwp) write(numout,*) '++++++ dta%v2d allocated space' 
    742749                  ALLOCATE( dta%v2d(nblen(3)) ) 
     
    749756         ENDIF 
    750757 
    751          IF ( nn_dyn3d_dta(ib_bdy) == 0 ) THEN 
     758         IF ( nn_dyn3d_dta(jbdy) == 0 ) THEN 
    752759            if(lwp) write(numout,*) '++++++ dta%u3d/v3d allocated space' 
    753             IF( dta%ll_u3d ) ALLOCATE( dta_bdy(ib_bdy)%u3d(nblen(2),jpk) ) 
    754             IF( dta%ll_v3d ) ALLOCATE( dta_bdy(ib_bdy)%v3d(nblen(3),jpk) ) 
    755          ENDIF 
    756          IF ( nn_dyn3d_dta(ib_bdy) == 1 .or. & 
    757            &  ( ln_full_vel_array(ib_bdy) .and. ( nn_dyn2d_dta(ib_bdy) == 1 .or. nn_dyn2d_dta(ib_bdy) == 3 ) ) ) THEN 
    758             IF ( dta%ll_u3d .or. ( ln_full_vel_array(ib_bdy) .and. dta%ll_u2d ) ) THEN 
     760            IF( dta%ll_u3d ) ALLOCATE( dta_bdy(jbdy)%u3d(nblen(2),jpk) ) 
     761            IF( dta%ll_v3d ) ALLOCATE( dta_bdy(jbdy)%v3d(nblen(3),jpk) ) 
     762         ENDIF 
     763         IF ( nn_dyn3d_dta(jbdy) == 1 .or. & 
     764           &  ( ln_full_vel_array(jbdy) .and. ( nn_dyn2d_dta(jbdy) == 1 .or. nn_dyn2d_dta(jbdy) == 3 ) ) ) THEN 
     765            IF ( dta%ll_u3d .or. ( ln_full_vel_array(jbdy) .and. dta%ll_u2d ) ) THEN 
    759766               if(lwp) write(numout,*) '++++++ dta%u3d pointing to fnow' 
    760767               jfld = jfld + 1 
    761                dta_bdy(ib_bdy)%u3d => bf(jfld)%fnow(:,1,:) 
    762             ENDIF 
    763             IF ( dta%ll_v3d .or. ( ln_full_vel_array(ib_bdy) .and. dta%ll_v2d ) ) THEN 
     768               dta_bdy(jbdy)%u3d => bf(jfld)%fnow(:,1,:) 
     769            ENDIF 
     770            IF ( dta%ll_v3d .or. ( ln_full_vel_array(jbdy) .and. dta%ll_v2d ) ) THEN 
    764771               if(lwp) write(numout,*) '++++++ dta%v3d pointing to fnow' 
    765772               jfld = jfld + 1 
    766                dta_bdy(ib_bdy)%v3d => bf(jfld)%fnow(:,1,:) 
    767             ENDIF 
    768          ENDIF 
    769  
    770          IF( nn_tra_dta(ib_bdy) == 0 ) THEN 
     773               dta_bdy(jbdy)%v3d => bf(jfld)%fnow(:,1,:) 
     774            ENDIF 
     775         ENDIF 
     776 
     777         IF( nn_tra_dta(jbdy) == 0 ) THEN 
    771778            if(lwp) write(numout,*) '++++++ dta%tem/sal allocated space' 
    772             IF( dta%ll_tem ) ALLOCATE( dta_bdy(ib_bdy)%tem(nblen(1),jpk) ) 
    773             IF( dta%ll_sal ) ALLOCATE( dta_bdy(ib_bdy)%sal(nblen(1),jpk) ) 
     779            IF( dta%ll_tem ) ALLOCATE( dta_bdy(jbdy)%tem(nblen(1),jpk) ) 
     780            IF( dta%ll_sal ) ALLOCATE( dta_bdy(jbdy)%sal(nblen(1),jpk) ) 
    774781         ELSE 
    775782            IF( dta%ll_tem ) THEN 
    776783               if(lwp) write(numout,*) '++++++ dta%tem pointing to fnow' 
    777784               jfld = jfld + 1 
    778                dta_bdy(ib_bdy)%tem => bf(jfld)%fnow(:,1,:) 
     785               dta_bdy(jbdy)%tem => bf(jfld)%fnow(:,1,:) 
    779786            ENDIF 
    780787            IF( dta%ll_sal ) THEN  
    781788               if(lwp) write(numout,*) '++++++ dta%sal pointing to fnow' 
    782789               jfld = jfld + 1 
    783                dta_bdy(ib_bdy)%sal => bf(jfld)%fnow(:,1,:) 
    784             ENDIF 
    785          ENDIF 
    786  
    787 #if defined key_lim3 
    788          IF (cn_ice_lim(ib_bdy) /= 'none') THEN 
    789             IF( nn_ice_lim_dta(ib_bdy) == 0 ) THEN 
    790                ALLOCATE( dta_bdy(ib_bdy)%a_i(nblen(1),jpl) ) 
    791                ALLOCATE( dta_bdy(ib_bdy)%h_i(nblen(1),jpl) ) 
    792                ALLOCATE( dta_bdy(ib_bdy)%h_s(nblen(1),jpl) ) 
     790               dta_bdy(jbdy)%sal => bf(jfld)%fnow(:,1,:) 
     791            ENDIF 
     792         ENDIF 
     793 
     794#if defined key_lim3 
     795         IF (cn_ice_lim(jbdy) /= 'none') THEN 
     796            IF( nn_ice_lim_dta(jbdy) == 0 ) THEN 
     797               ALLOCATE( dta_bdy(jbdy)%a_i(nblen(1),jpl) ) 
     798               ALLOCATE( dta_bdy(jbdy)%h_i(nblen(1),jpl) ) 
     799               ALLOCATE( dta_bdy(jbdy)%h_s(nblen(1),jpl) ) 
    793800            ELSE 
    794                IF ( ll_bdylim3 ) THEN ! case input is Xcat 
    795                   jfld = jfld + 1 
    796                   dta_bdy(ib_bdy)%a_i => bf(jfld)%fnow(:,1,:) 
    797                   jfld = jfld + 1 
    798                   dta_bdy(ib_bdy)%h_i => bf(jfld)%fnow(:,1,:) 
    799                   jfld = jfld + 1 
    800                   dta_bdy(ib_bdy)%h_s => bf(jfld)%fnow(:,1,:) 
    801                ELSE ! case input is 1cat 
     801               IF ( nice_cat == jpl ) THEN ! case input cat = jpl 
     802                  jfld = jfld + 1 
     803                  dta_bdy(jbdy)%a_i => bf(jfld)%fnow(:,1,:) 
     804                  jfld = jfld + 1 
     805                  dta_bdy(jbdy)%h_i => bf(jfld)%fnow(:,1,:) 
     806                  jfld = jfld + 1 
     807                  dta_bdy(jbdy)%h_s => bf(jfld)%fnow(:,1,:) 
     808               ELSE                        ! case input cat = 1 OR (/=1 and /=jpl) 
    802809                  jfld_ai  = jfld + 1 
    803810                  jfld_hti = jfld + 2 
    804811                  jfld_hts = jfld + 3 
    805812                  jfld     = jfld + 3 
    806                   ALLOCATE( dta_bdy(ib_bdy)%a_i(nblen(1),jpl) ) 
    807                   ALLOCATE( dta_bdy(ib_bdy)%h_i(nblen(1),jpl) ) 
    808                   ALLOCATE( dta_bdy(ib_bdy)%h_s(nblen(1),jpl) ) 
    809                   dta_bdy(ib_bdy)%a_i(:,:) = 0._wp 
    810                   dta_bdy(ib_bdy)%h_i(:,:) = 0._wp 
    811                   dta_bdy(ib_bdy)%h_s(:,:) = 0._wp 
     813                  ALLOCATE( dta_bdy(jbdy)%a_i(nblen(1),jpl) ) 
     814                  ALLOCATE( dta_bdy(jbdy)%h_i(nblen(1),jpl) ) 
     815                  ALLOCATE( dta_bdy(jbdy)%h_s(nblen(1),jpl) ) 
     816                  dta_bdy(jbdy)%a_i(:,:) = 0._wp 
     817                  dta_bdy(jbdy)%h_i(:,:) = 0._wp 
     818                  dta_bdy(jbdy)%h_s(:,:) = 0._wp 
    812819               ENDIF 
    813820 
     
    816823#endif 
    817824         ! 
    818       END DO ! ib_bdy  
     825      END DO ! jbdy  
    819826      ! 
    820827      IF( nn_timing == 1 ) CALL timing_stop('bdy_dta_init') 
  • branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_ice.F90

    r8586 r8813  
    4848   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   alb_ice        !: ice albedo                                       [-] 
    4949 
     50   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qml_ice        !: heat available for snow / ice surface melting     [W/m2]  
     51   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qcn_ice        !: heat conduction flux in the layer below surface   [W/m2]  
     52   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qsr_ice_tr     !: solar flux transmitted below the ice surface      [W/m2] 
     53 
    5054   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   utau_ice       !: atmos-ice u-stress. VP: I-pt ; EVP: U,V-pts   [N/m2] 
    5155   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   vtau_ice       !: atmos-ice v-stress. VP: I-pt ; EVP: U,V-pts   [N/m2] 
    52    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   fr1_i0         !: Solar surface transmission parameter, thick ice  [-] 
    53    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   fr2_i0         !: Solar surface transmission parameter, thin ice   [-] 
    5456   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   emp_ice        !: sublimation - precip over sea ice          [kg/m2/s] 
    5557 
     
    123125      ALLOCATE( qns_ice (jpi,jpj,jpl) , qsr_ice (jpi,jpj,jpl) ,     & 
    124126         &      qla_ice (jpi,jpj,jpl) , dqla_ice(jpi,jpj,jpl) ,     & 
    125          &      dqns_ice(jpi,jpj,jpl) , tn_ice  (jpi,jpj,jpl) , alb_ice (jpi,jpj,jpl) ,   & 
     127         &      dqns_ice(jpi,jpj,jpl) , tn_ice  (jpi,jpj,jpl)  , alb_ice (jpi,jpj,jpl) ,   & 
     128         &      qml_ice(jpi,jpj,jpl)  , qcn_ice(jpi,jpj,jpl)   , qsr_ice_tr(jpi,jpj,jpl),  & 
    126129         &      utau_ice(jpi,jpj)     , vtau_ice(jpi,jpj)     , wndm_ice(jpi,jpj)     ,   & 
    127          &      fr1_i0  (jpi,jpj)     , fr2_i0  (jpi,jpj)     ,     & 
    128130         &      evap_ice(jpi,jpj,jpl) , devap_ice(jpi,jpj,jpl) , qprec_ice(jpi,jpj) ,   & 
    129131         &      qemp_ice(jpi,jpj)     , qevap_ice(jpi,jpj,jpl) , qemp_oce (jpi,jpj) ,   & 
     
    139141                a_i(jpi,jpj,ncat)     , topmelt(jpi,jpj,ncat) , botmelt(jpi,jpj,ncat) , & 
    140142                STAT= ierr(2) ) 
    141       IF( ln_cpl )   ALLOCATE( u_ice(jpi,jpj)        , fr1_i0(jpi,jpj)       , tn_ice (jpi,jpj,1)    , & 
    142          &                     v_ice(jpi,jpj)        , fr2_i0(jpi,jpj)       , alb_ice(jpi,jpj,1)    , & 
     143      IF( ln_cpl )   ALLOCATE( u_ice(jpi,jpj)        , tn_ice (jpi,jpj,1)    , & 
     144         &                     v_ice(jpi,jpj)        , alb_ice(jpi,jpj,1)    , & 
    143145         &                     emp_ice(jpi,jpj)      , qns_ice(jpi,jpj,1)    , dqns_ice(jpi,jpj,1)   , & 
    144146         &                     STAT= ierr(3) )       
     
    161163   PRIVATE 
    162164 
    163    PUBLIC sbc_ice_alloc 
     165   PUBLIC   sbc_ice_alloc   ! 
    164166 
    165167   LOGICAL         , PUBLIC, PARAMETER ::   lk_lim3    = .FALSE.  !: no LIM-3 ice model 
     
    168170   REAL(wp)        , PUBLIC, PARAMETER ::   cldf_ice = 0.81       !: cloud fraction over sea ice, summer CLIO value   [-] 
    169171   INTEGER         , PUBLIC, PARAMETER ::   jpl = 1  
    170    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   u_ice, v_ice,fr1_i0,fr2_i0          ! jpi, jpj 
     172   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   u_ice, v_ice                        ! jpi, jpj 
    171173   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   tn_ice, alb_ice, qns_ice, dqns_ice  ! (jpi,jpj,jpl) 
    172174   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   a_i 
     
    176178   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   topmelt, botmelt 
    177179   ! 
    178    !! arrays relating to embedding ice in the ocean. These arrays need to be declared  
    179    !! even if no ice model is required. In the no ice model or traditional levitating  
    180    !! ice cases they contain only zeros 
     180   !! arrays related to embedding ice in the ocean.  
     181   !! These arrays need to be declared even if no ice model is required.  
     182   !! In the no ice model or traditional levitating ice cases they contain only zeros 
    181183   !! --------------------- 
    182184   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   snwice_mass        !: mass of snow and ice at current  ice time step   [Kg/m2] 
  • branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90

    r8637 r8813  
    1515   !!            3.7  !  2014-06  (L. Brodeau)  simplification and optimization of CORE bulk 
    1616   !!            4.0  !  2016-06  (L. Brodeau)  sbcblk_core becomes sbcblk and is not restricted to the CORE algorithm anymore 
    17    !!                                          ==> based on AeroBulk (http://aerobulk.sourceforge.net/) 
     17   !!                 !                        ==> based on AeroBulk (http://aerobulk.sourceforge.net/) 
    1818   !!            4.0  !  2016-10  (G. Madec)  introduce a sbc_blk_init routine 
     19   !!            4.0  !  2016-10  (M. Vancoppenolle)  Introduce Jules emulator (M. Vancoppenolle)  
    1920   !!---------------------------------------------------------------------- 
    2021 
     
    2324   !!   sbc_blk       : bulk formulation as ocean surface boundary condition 
    2425   !!   blk_oce       : computes momentum, heat and freshwater fluxes over ocean 
    25    !!   blk_ice       : computes momentum, heat and freshwater fluxes over sea ice 
    2626   !!   rho_air       : density of (moist) air (depends on T_air, q_air and SLP 
    2727   !!   cp_air        : specific heat of (moist) air (depends spec. hum. q_air) 
    2828   !!   q_sat         : saturation humidity as a function of SLP and temperature 
    2929   !!   L_vap         : latent heat of vaporization of water as a function of temperature 
     30   !!             sea-ice case only :  
     31   !!   blk_ice_tau   : provide the air-ice stress 
     32   !!   blk_ice_flx   : provide the heat and mass fluxes at air-ice interface 
     33   !!   blk_ice_qcn   : provide ice surface temperature and snow/ice conduction flux (emulating JULES coupler) 
     34   !!   Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag 
     35   !!   Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag  
    3036   !!---------------------------------------------------------------------- 
    3137   USE oce            ! ocean dynamics and tracers 
     
    4248   USE ice     , ONLY :   u_ice, v_ice, jpl, a_i_b, at_i_b, tm_su 
    4349   USE icethd_dh      ! for CALL ice_thd_snwblow 
     50   USE icethd_zdf,ONLY:   rn_cnd_s  ! for blk_ice_qcn 
    4451#endif 
    4552   USE sbcblk_algo_ncar     ! => turb_ncar     : NCAR - CORE (Large & Yeager, 2009)  
     
    6370   PUBLIC   blk_ice_tau   ! routine called in icestp module 
    6471   PUBLIC   blk_ice_flx   ! routine called in icestp module 
    65 #endif 
     72   PUBLIC   blk_ice_qcn   ! routine called in icestp module 
     73#endif  
    6674 
    6775!!Lolo: should ultimately be moved in the module with all physical constants ? 
     
    111119   LOGICAL  ::   ln_Cd_L15 = .FALSE. !  Modify the drag ice-atm depending on ice concentration (from Lupkes et al. JGR2015) 
    112120   ! 
    113 !!cr   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Cd_oce   ! air-ocean drag (clem) 
    114121   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Cd_atm                    ! transfer coefficient for momentum      (tau) 
    115122   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Ch_atm                    ! transfer coefficient for sensible heat (Q_sens) 
     
    139146      !!             ***  ROUTINE sbc_blk_alloc *** 
    140147      !!------------------------------------------------------------------- 
    141 !!cr      ALLOCATE( Cd_oce(jpi,jpj) , STAT=sbc_blk_alloc ) 
    142148      ALLOCATE( Cd_atm (jpi,jpj), Ch_atm (jpi,jpj), Ce_atm (jpi,jpj), t_zu(jpi,jpj), q_zu(jpi,jpj), & 
    143149         &      cdn_oce(jpi,jpj), chn_oce(jpi,jpj), cen_oce(jpi,jpj), STAT=sbc_blk_alloc ) 
     
    147153   END FUNCTION sbc_blk_alloc 
    148154 
     155 
    149156   SUBROUTINE sbc_blk_init 
    150157      !!--------------------------------------------------------------------- 
     
    154161      !! 
    155162      !! ** Method  :  
    156       !! 
    157       !!      C A U T I O N : never mask the surface stress fields 
    158       !!                      the stress is assumed to be in the (i,j) mesh referential 
    159       !! 
    160       !! ** Action  :    
    161163      !! 
    162164      !!---------------------------------------------------------------------- 
     
    285287      !! 
    286288      !! ** Purpose :   provide at each time step the surface ocean fluxes 
    287       !!      (momentum, heat, freshwater and runoff) 
     289      !!              (momentum, heat, freshwater and runoff) 
    288290      !! 
    289291      !! ** Method  : (1) READ each fluxes in NetCDF files: 
     
    566568   END SUBROUTINE blk_oce 
    567569 
     570 
     571 
     572   FUNCTION rho_air( ptak, pqa, pslp ) 
     573      !!------------------------------------------------------------------------------- 
     574      !!                           ***  FUNCTION rho_air  *** 
     575      !! 
     576      !! ** Purpose : compute density of (moist) air using the eq. of state of the atmosphere 
     577      !! 
     578      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)  
     579      !!------------------------------------------------------------------------------- 
     580      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak      ! air temperature             [K] 
     581      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa       ! air specific humidity   [kg/kg] 
     582      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pslp      ! pressure in                [Pa] 
     583      REAL(wp), DIMENSION(jpi,jpj)             ::   rho_air   ! density of moist air   [kg/m^3] 
     584      !!------------------------------------------------------------------------------- 
     585      ! 
     586      rho_air = pslp / (  R_dry*ptak * ( 1._wp + rctv0*pqa )  ) 
     587      ! 
     588   END FUNCTION rho_air 
     589 
     590 
     591   FUNCTION cp_air( pqa ) 
     592      !!------------------------------------------------------------------------------- 
     593      !!                           ***  FUNCTION cp_air  *** 
     594      !! 
     595      !! ** Purpose : Compute specific heat (Cp) of moist air 
     596      !! 
     597      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
     598      !!------------------------------------------------------------------------------- 
     599      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa      ! air specific humidity         [kg/kg] 
     600      REAL(wp), DIMENSION(jpi,jpj)             ::   cp_air   ! specific heat of moist air   [J/K/kg] 
     601      !!------------------------------------------------------------------------------- 
     602      ! 
     603      Cp_air = Cp_dry + Cp_vap * pqa 
     604      ! 
     605   END FUNCTION cp_air 
     606 
     607 
     608   FUNCTION q_sat( ptak, pslp ) 
     609      !!---------------------------------------------------------------------------------- 
     610      !!                           ***  FUNCTION q_sat  *** 
     611      !! 
     612      !! ** Purpose : Specific humidity at saturation in [kg/kg] 
     613      !!              Based on accurate estimate of "e_sat" 
     614      !!              aka saturation water vapor (Goff, 1957) 
     615      !! 
     616      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
     617      !!---------------------------------------------------------------------------------- 
     618      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak    ! air temperature                       [K] 
     619      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pslp    ! sea level atmospheric pressure       [Pa] 
     620      REAL(wp), DIMENSION(jpi,jpj)             ::   q_sat   ! Specific humidity at saturation   [kg/kg] 
     621      ! 
     622      INTEGER  ::   ji, jj         ! dummy loop indices 
     623      REAL(wp) ::   ze_sat, ztmp   ! local scalar 
     624      !!---------------------------------------------------------------------------------- 
     625      ! 
     626      DO jj = 1, jpj 
     627         DO ji = 1, jpi 
     628            ! 
     629            ztmp = rt0 / ptak(ji,jj) 
     630            ! 
     631            ! Vapour pressure at saturation [hPa] : WMO, (Goff, 1957) 
     632            ze_sat = 10.**( 10.79574*(1. - ztmp) - 5.028*LOG10(ptak(ji,jj)/rt0)        & 
     633               &    + 1.50475*10.**(-4)*(1. - 10.**(-8.2969*(ptak(ji,jj)/rt0 - 1.)) )  & 
     634               &    + 0.42873*10.**(-3)*(10.**(4.76955*(1. - ztmp)) - 1.) + 0.78614  ) 
     635               ! 
     636            q_sat(ji,jj) = reps0 * ze_sat/( 0.01_wp*pslp(ji,jj) - (1._wp - reps0)*ze_sat )   ! 0.01 because SLP is in [Pa] 
     637            ! 
     638         END DO 
     639      END DO 
     640      ! 
     641   END FUNCTION q_sat 
     642 
     643 
     644   FUNCTION gamma_moist( ptak, pqa ) 
     645      !!---------------------------------------------------------------------------------- 
     646      !!                           ***  FUNCTION gamma_moist  *** 
     647      !! 
     648      !! ** Purpose : Compute the moist adiabatic lapse-rate. 
     649      !!     => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate 
     650      !!     => http://www.geog.ucsb.edu/~joel/g266_s10/lecture_notes/chapt03/oh10_3_01/oh10_3_01.html 
     651      !! 
     652      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
     653      !!---------------------------------------------------------------------------------- 
     654      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak          ! air temperature       [K] 
     655      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa           ! specific humidity [kg/kg] 
     656      REAL(wp), DIMENSION(jpi,jpj)             ::   gamma_moist   ! moist adiabatic lapse-rate 
     657      ! 
     658      INTEGER  ::   ji, jj         ! dummy loop indices 
     659      REAL(wp) :: zrv, ziRT        ! local scalar 
     660      !!---------------------------------------------------------------------------------- 
     661      ! 
     662      DO jj = 1, jpj 
     663         DO ji = 1, jpi 
     664            zrv = pqa(ji,jj) / (1. - pqa(ji,jj)) 
     665            ziRT = 1. / (R_dry*ptak(ji,jj))    ! 1/RT 
     666            gamma_moist(ji,jj) = grav * ( 1. + cevap*zrv*ziRT ) / ( Cp_dry + cevap*cevap*zrv*reps0*ziRT/ptak(ji,jj) ) 
     667         END DO 
     668      END DO 
     669      ! 
     670   END FUNCTION gamma_moist 
     671 
     672 
     673   FUNCTION L_vap( psst ) 
     674      !!--------------------------------------------------------------------------------- 
     675      !!                           ***  FUNCTION L_vap  *** 
     676      !! 
     677      !! ** Purpose : Compute the latent heat of vaporization of water from temperature 
     678      !! 
     679      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
     680      !!---------------------------------------------------------------------------------- 
     681      REAL(wp), DIMENSION(jpi,jpj)             ::   L_vap   ! latent heat of vaporization   [J/kg] 
     682      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   psst   ! water temperature                [K] 
     683      !!---------------------------------------------------------------------------------- 
     684      ! 
     685      L_vap = (  2.501 - 0.00237 * ( psst(:,:) - rt0)  ) * 1.e6 
     686      ! 
     687   END FUNCTION L_vap 
     688 
    568689#if defined key_lim3 
     690   !!---------------------------------------------------------------------- 
     691   !!   'key_lim3'                                       ESIM sea-ice model 
     692   !!---------------------------------------------------------------------- 
     693   !!   blk_ice_tau : provide the air-ice stress 
     694   !!   blk_ice_flx : provide the heat and mass fluxes at air-ice interface 
     695   !!   blk_ice_qcn : provide ice surface temperature and snow/ice conduction flux (emulating JULES coupler) 
     696   !!   Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag 
     697   !!   Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag  
     698   !!---------------------------------------------------------------------- 
    569699 
    570700   SUBROUTINE blk_ice_tau 
     
    688818 
    689819 
    690    SUBROUTINE blk_ice_flx( ptsu, palb ) 
     820   SUBROUTINE blk_ice_flx( ptsu, phs, phi, palb ) 
    691821      !!--------------------------------------------------------------------- 
    692822      !!                     ***  ROUTINE blk_ice_flx  *** 
    693823      !! 
    694       !! ** Purpose :   provide the surface boundary condition over sea-ice 
     824      !! ** Purpose :   provide the heat and mass fluxes at air-ice interface 
    695825      !! 
    696826      !! ** Method  :   compute heat and freshwater exchanged 
     
    701831      !!--------------------------------------------------------------------- 
    702832      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   ptsu   ! sea ice surface temperature 
     833      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phs    ! snow thickness 
     834      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phi    ! ice thickness 
    703835      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   palb   ! ice albedo (all skies) 
    704836      !! 
     
    707839      REAL(wp) ::   zcoef_dqlw, zcoef_dqla   !   -      - 
    708840      REAL(wp) ::   zztmp, z1_lsub           !   -      - 
     841      REAL(wp) ::   zfrqsr_tr_i              ! sea ice shortwave fraction transmitted below through the ice 
     842      REAL(wp) ::   zfr1, zfr2, zfac         ! local variables 
    709843      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_qlw         ! long wave heat flux over ice 
    710844      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_qsb         ! sensible  heat flux over ice 
     
    717851      IF( ln_timing )   CALL timing_start('blk_ice_flx') 
    718852      ! 
    719       ! 
    720       ! local scalars 
    721       zcoef_dqlw   = 4.0 * 0.95 * Stef 
    722       zcoef_dqla   = -Ls * 11637800. * (-5897.8) 
     853      zcoef_dqlw = 4.0 * 0.95 * Stef             ! local scalars 
     854      zcoef_dqla = -Ls * 11637800. * (-5897.8) 
    723855      ! 
    724856      zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 
     
    815947      END DO 
    816948 
    817       !-------------------------------------------------------------------- 
    818       ! FRACTIONs of net shortwave radiation which is not absorbed in the 
    819       ! thin surface layer and penetrates inside the ice cover 
    820       ! ( Maykut and Untersteiner, 1971 ; Ebert and Curry, 1993 ) 
    821       ! 
    822       fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice ) 
    823       fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) 
     949      ! --- absorbed and transmitted shortwave radiation (W/m2) --- ! 
     950      ! 
     951      ! former coding was 
     952      !     fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice ) 
     953      !     fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) 
     954 
     955      ! --- surface transmission parameter (i0, Grenfell Maykut 77) --- ! 
     956      zfr1 = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )            !   standard value 
     957      zfr2 = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )            !   zfr2 such that zfr1 + zfr2 to equal 1 
     958 
     959      qsr_ice_tr(:,:,:) = 0._wp 
     960 
     961      DO jl = 1, jpl 
     962         DO jj = 1, jpj 
     963            DO ji = 1, jpi 
     964               ! 
     965               zfac = MAX( 0._wp , 1._wp - phi(ji,jj,jl) * 10._wp )     !   linear weighting factor: =0 for phi=0, =1 for phi = 0.1 
     966               zfrqsr_tr_i = zfr1 + zfac * zfr2                         !   below 10 cm, linearly increase zfrqsr_tr_i until 1 at zero thickness 
     967               ! 
     968               IF ( phs(ji,jj,jl) <= 0._wp ) THEN   ;    zfrqsr_tr_i  = zfr1 + zfac * zfr2 
     969               ELSE                                 ;    zfrqsr_tr_i  = 0._wp                                  !   snow fully opaque 
     970               ENDIF 
     971               ! 
     972               qsr_ice_tr(ji,jj,jl) = zfrqsr_tr_i * qsr_ice(ji,jj,jl)   !   transmitted solar radiation  
     973               ! 
     974            END DO 
     975         END DO 
     976      END DO 
    824977      ! 
    825978      IF(ln_ctl) THEN 
     
    836989   END SUBROUTINE blk_ice_flx 
    837990    
    838 #endif 
    839  
    840    FUNCTION rho_air( ptak, pqa, pslp ) 
    841       !!------------------------------------------------------------------------------- 
    842       !!                           ***  FUNCTION rho_air  *** 
    843       !! 
    844       !! ** Purpose : compute density of (moist) air using the eq. of state of the atmosphere 
    845       !! 
    846       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)  
    847       !!------------------------------------------------------------------------------- 
    848       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak      ! air temperature             [K] 
    849       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa       ! air specific humidity   [kg/kg] 
    850       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pslp      ! pressure in                [Pa] 
    851       REAL(wp), DIMENSION(jpi,jpj)             ::   rho_air   ! density of moist air   [kg/m^3] 
    852       !!------------------------------------------------------------------------------- 
    853       ! 
    854       rho_air = pslp / (  R_dry*ptak * ( 1._wp + rctv0*pqa )  ) 
    855       ! 
    856    END FUNCTION rho_air 
    857  
    858  
    859    FUNCTION cp_air( pqa ) 
    860       !!------------------------------------------------------------------------------- 
    861       !!                           ***  FUNCTION cp_air  *** 
    862       !! 
    863       !! ** Purpose : Compute specific heat (Cp) of moist air 
    864       !! 
    865       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
    866       !!------------------------------------------------------------------------------- 
    867       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa      ! air specific humidity         [kg/kg] 
    868       REAL(wp), DIMENSION(jpi,jpj)             ::   cp_air   ! specific heat of moist air   [J/K/kg] 
    869       !!------------------------------------------------------------------------------- 
    870       ! 
    871       Cp_air = Cp_dry + Cp_vap * pqa 
    872       ! 
    873    END FUNCTION cp_air 
    874  
    875  
    876    FUNCTION q_sat( ptak, pslp ) 
    877       !!---------------------------------------------------------------------------------- 
    878       !!                           ***  FUNCTION q_sat  *** 
    879       !! 
    880       !! ** Purpose : Specific humidity at saturation in [kg/kg] 
    881       !!              Based on accurate estimate of "e_sat" 
    882       !!              aka saturation water vapor (Goff, 1957) 
    883       !! 
    884       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
    885       !!---------------------------------------------------------------------------------- 
    886       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak    ! air temperature                       [K] 
    887       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pslp    ! sea level atmospheric pressure       [Pa] 
    888       REAL(wp), DIMENSION(jpi,jpj)             ::   q_sat   ! Specific humidity at saturation   [kg/kg] 
    889       ! 
    890       INTEGER  ::   ji, jj         ! dummy loop indices 
    891       REAL(wp) ::   ze_sat, ztmp   ! local scalar 
    892       !!---------------------------------------------------------------------------------- 
    893       ! 
    894       DO jj = 1, jpj 
    895          DO ji = 1, jpi 
    896             ! 
    897             ztmp = rt0 / ptak(ji,jj) 
    898             ! 
    899             ! Vapour pressure at saturation [hPa] : WMO, (Goff, 1957) 
    900             ze_sat = 10.**( 10.79574*(1. - ztmp) - 5.028*LOG10(ptak(ji,jj)/rt0)        & 
    901                &    + 1.50475*10.**(-4)*(1. - 10.**(-8.2969*(ptak(ji,jj)/rt0 - 1.)) )  & 
    902                &    + 0.42873*10.**(-3)*(10.**(4.76955*(1. - ztmp)) - 1.) + 0.78614  ) 
     991 
     992   SUBROUTINE blk_ice_qcn( k_monocat, ptsu, ptb, phs, phi ) 
     993      !!--------------------------------------------------------------------- 
     994      !!                     ***  ROUTINE blk_ice_qcn  *** 
     995      !! 
     996      !! ** Purpose :   Compute surface temperature and snow/ice conduction flux 
     997      !!                to force sea ice / snow thermodynamics 
     998      !!                in the case JULES coupler is emulated 
     999      !!                 
     1000      !! ** Method  :   compute surface energy balance assuming neglecting heat storage 
     1001      !!                following the 0-layer Semtner (1976) approach 
     1002      !! 
     1003      !! ** Outputs : - ptsu    : sea-ice / snow surface temperature (K) 
     1004      !!              - qcn_ice : surface inner conduction flux (W/m2) 
     1005      !! 
     1006      !!--------------------------------------------------------------------- 
     1007      INTEGER                   , INTENT(in   ) ::   k_monocat   ! single-category option 
     1008      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptsu        ! sea ice / snow surface temperature 
     1009      REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   ptb         ! sea ice base temperature 
     1010      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   phs         ! snow thickness 
     1011      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   phi         ! sea ice thickness 
     1012      ! 
     1013      INTEGER , PARAMETER ::   nit = 10                  ! number of iterations 
     1014      REAL(wp), PARAMETER ::   zepsilon = 0.1_wp         ! characteristic thickness for enhanced conduction 
     1015      ! 
     1016      INTEGER  ::   ji, jj, jl           ! dummy loop indices 
     1017      INTEGER  ::   iter                 ! local integer 
     1018      REAL(wp) ::   zfac, zfac2, zfac3   ! local scalars 
     1019      REAL(wp) ::   zkeff_h, ztsu        ! 
     1020      REAL(wp) ::   zqc, zqnet           ! 
     1021      REAL(wp) ::   zhe, zqa0            ! 
     1022      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zgfac   ! enhanced conduction factor 
     1023      !!--------------------------------------------------------------------- 
     1024 
     1025      IF( nn_timing == 1 )  CALL timing_start('blk_ice_qcn') 
     1026       
     1027      ! -------------------------------------! 
     1028      !      I   Enhanced conduction factor  ! 
     1029      ! -------------------------------------! 
     1030      ! Emulates the enhancement of conduction by unresolved thin ice (k_monocat = 1/3) 
     1031      ! Fichefet and Morales Maqueda, JGR 1997 
     1032      ! 
     1033      zgfac(:,:,:) = 1._wp 
     1034       
     1035      SELECT CASE ( k_monocat )  
     1036      ! 
     1037      CASE ( 1 , 3 ) 
     1038         ! 
     1039         zfac    =   1._wp /  ( rn_cnd_s + rcdic ) 
     1040         zfac2   =   EXP(1._wp) * 0.5_wp * zepsilon 
     1041         zfac3   =   2._wp / zepsilon 
     1042         !    
     1043         DO jl = 1, jpl                 
     1044            DO jj = 1 , jpj 
     1045               DO ji = 1, jpi 
     1046                  !                                ! Effective thickness 
     1047                  zhe               =   ( rn_cnd_s * phi(ji,jj,jl) + rcdic * phs(ji,jj,jl) ) * zfac 
     1048                  !                                ! Enhanced conduction factor 
     1049                  IF( zhe >=  zfac2 )  & 
     1050                     zgfac(ji,jj,jl)   =   MIN( 2._wp, ( 0.5_wp + 0.5 * LOG( zhe * zfac3 ) ) ) 
     1051               END DO 
     1052            END DO 
     1053         END DO 
     1054         !       
     1055      END SELECT 
     1056       
     1057      ! -------------------------------------------------------------! 
     1058      !      II   Surface temperature and conduction flux            ! 
     1059      ! -------------------------------------------------------------! 
     1060      ! 
     1061      zfac   =   rcdic * rn_cnd_s  
     1062      !                             ! ========================== ! 
     1063      DO jl = 1, jpl                !  Loop over ice categories  ! 
     1064         !                          ! ========================== ! 
     1065         DO jj = 1 , jpj 
     1066            DO ji = 1, jpi 
     1067               !                    ! Effective conductivity of the snow-ice system divided by thickness 
     1068               zkeff_h =   zfac * zgfac(ji,jj,jl) / ( rcdic * phs(ji,jj,jl) + rn_cnd_s * phi(ji,jj,jl) ) 
     1069               !                    ! Store initial surface temperature 
     1070               ztsu    =   ptsu(ji,jj,jl) 
     1071               !                    ! Net initial atmospheric heat flux 
     1072               zqa0    =   qsr_ice(ji,jj,jl) - qsr_ice_tr(ji,jj,jl) + qns_ice(ji,jj,jl) 
    9031073               ! 
    904             q_sat(ji,jj) = reps0 * ze_sat/( 0.01_wp*pslp(ji,jj) - (1._wp - reps0)*ze_sat )   ! 0.01 because SLP is in [Pa] 
    905             ! 
    906          END DO 
    907       END DO 
    908       ! 
    909    END FUNCTION q_sat 
    910  
    911  
    912    FUNCTION gamma_moist( ptak, pqa ) 
    913       !!---------------------------------------------------------------------------------- 
    914       !!                           ***  FUNCTION gamma_moist  *** 
    915       !! 
    916       !! ** Purpose : Compute the moist adiabatic lapse-rate. 
    917       !!     => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate 
    918       !!     => http://www.geog.ucsb.edu/~joel/g266_s10/lecture_notes/chapt03/oh10_3_01/oh10_3_01.html 
    919       !! 
    920       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
    921       !!---------------------------------------------------------------------------------- 
    922       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak          ! air temperature       [K] 
    923       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa           ! specific humidity [kg/kg] 
    924       REAL(wp), DIMENSION(jpi,jpj)             ::   gamma_moist   ! moist adiabatic lapse-rate 
    925       ! 
    926       INTEGER  ::   ji, jj         ! dummy loop indices 
    927       REAL(wp) :: zrv, ziRT        ! local scalar 
    928       !!---------------------------------------------------------------------------------- 
    929       ! 
    930       DO jj = 1, jpj 
    931          DO ji = 1, jpi 
    932             zrv = pqa(ji,jj) / (1. - pqa(ji,jj)) 
    933             ziRT = 1. / (R_dry*ptak(ji,jj))    ! 1/RT 
    934             gamma_moist(ji,jj) = grav * ( 1. + cevap*zrv*ziRT ) / ( Cp_dry + cevap*cevap*zrv*reps0*ziRT/ptak(ji,jj) ) 
    935          END DO 
    936       END DO 
    937       ! 
    938    END FUNCTION gamma_moist 
    939  
    940  
    941    FUNCTION L_vap( psst ) 
    942       !!--------------------------------------------------------------------------------- 
    943       !!                           ***  FUNCTION L_vap  *** 
    944       !! 
    945       !! ** Purpose : Compute the latent heat of vaporization of water from temperature 
    946       !! 
    947       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
    948       !!---------------------------------------------------------------------------------- 
    949       REAL(wp), DIMENSION(jpi,jpj)             ::   L_vap   ! latent heat of vaporization   [J/kg] 
    950       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   psst   ! water temperature                [K] 
    951       !!---------------------------------------------------------------------------------- 
    952       ! 
    953       L_vap = (  2.501 - 0.00237 * ( psst(:,:) - rt0)  ) * 1.e6 
    954       ! 
    955    END FUNCTION L_vap 
    956  
    957 #if defined key_lim3 
     1074               DO iter = 1, nit     ! --- Iteration loop 
     1075                   !                      ! Conduction heat flux through snow-ice system (>0 downwards) 
     1076                   zqc   =   zkeff_h * ( ztsu - ptb(ji,jj) ) 
     1077                   !                      ! Surface energy budget 
     1078                   zqnet =   zqa0 + dqns_ice(ji,jj,jl) * ( ztsu - ptsu(ji,jj,jl) ) - zqc 
     1079                   !                      ! Temperature update 
     1080                   ztsu  =   ztsu - zqnet / ( dqns_ice(ji,jj,jl) - zkeff_h ) 
     1081               END DO 
     1082               ! 
     1083               ptsu   (ji,jj,jl) = MIN( rt0, ztsu ) 
     1084               qcn_ice(ji,jj,jl) = zkeff_h * ( ptsu(ji,jj,jl) - ptb(ji,jj) ) 
     1085               ! 
     1086            END DO 
     1087         END DO 
     1088         ! 
     1089      END DO  
     1090      !       
     1091      IF( nn_timing == 1 )  CALL timing_stop('blk_ice_qcn') 
     1092      ! 
     1093   END SUBROUTINE blk_ice_qcn 
     1094    
    9581095 
    9591096   SUBROUTINE Cdn10_Lupkes2012( Cd ) 
     
    9611098      !!                      ***  ROUTINE  Cdn10_Lupkes2012  *** 
    9621099      !! 
    963       !! ** Purpose :    Recompute the ice-atm drag at 10m height to make 
    964       !!                 it dependent on edges at leads, melt ponds and flows. 
     1100      !! ** Purpose :    Recompute the neutral air-ice drag referenced at 10m  
     1101      !!                 to make it dependent on edges at leads, melt ponds and flows. 
    9651102      !!                 After some approximations, this can be resumed to a dependency 
    9661103      !!                 on ice concentration. 
     
    10121149      !!                      ***  ROUTINE  Cdn10_Lupkes2015  *** 
    10131150      !! 
    1014       !! ** pUrpose :    1lternative turbulent transfert coefficients formulation 
     1151      !! ** pUrpose :    Alternative turbulent transfert coefficients formulation 
    10151152      !!                 between sea-ice and atmosphere with distinct momentum  
    10161153      !!                 and heat coefficients depending on sea-ice concentration  
  • branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_algo_coare.F90

    r8637 r8813  
    4949   PUBLIC ::   TURB_COARE   ! called by sbcblk.F90 
    5050 
    51    !                                           !! COARE own values for given constants: 
    52    REAL(wp), PARAMETER ::   zi0     = 600.      ! scale height of the atmospheric boundary layer...1 
    53    REAL(wp), PARAMETER ::   Beta0   =   1.25    ! gustiness parameter 
    54    REAL(wp), PARAMETER ::   rctv0   =   0.608   ! constant to obtain virtual temperature... 
     51   !                                              !! COARE own values for given constants: 
     52   REAL(wp), PARAMETER ::   zi0     = 600._wp      ! scale height of the atmospheric boundary layer... 
     53   REAL(wp), PARAMETER ::   Beta0   =   1.250_wp   ! gustiness parameter 
     54   REAL(wp), PARAMETER ::   rctv0   =   0.608_wp   ! constant to obtain virtual temperature... 
    5555 
    5656   !!---------------------------------------------------------------------- 
  • branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90

    r8637 r8813  
    971971      !!                        emp          upward mass flux [evap. - precip. (- runoffs) (- calving)] (ocean only case) 
    972972      !!---------------------------------------------------------------------- 
    973       USE zdf_oce,  ONLY : ln_zdfswm 
    974  
    975       IMPLICIT NONE 
    976  
    977       INTEGER, INTENT(in)           ::   kt          ! ocean model time step index 
    978       INTEGER, INTENT(in)           ::   k_fsbc      ! frequency of sbc (-> ice model) computation  
    979       INTEGER, INTENT(in)           ::   k_ice       ! ice management in the sbc (=0/1/2/3) 
     973      USE zdf_oce,  ONLY :   ln_zdfswm 
     974      ! 
     975      INTEGER, INTENT(in) ::   kt          ! ocean model time step index 
     976      INTEGER, INTENT(in) ::   k_fsbc      ! frequency of sbc (-> ice model) computation  
     977      INTEGER, INTENT(in) ::   k_ice       ! ice management in the sbc (=0/1/2/3) 
    980978      !! 
    981979      LOGICAL  ::   llnewtx, llnewtau      ! update wind stress components and module?? 
     
    11701168         ! Calculate the 3D Stokes drift both in coupled and not fully uncoupled mode 
    11711169         IF( srcv(jpr_sdrftx)%laction .OR. srcv(jpr_sdrfty)%laction .OR. srcv(jpr_wper)%laction & 
    1172                                                                     .OR. srcv(jpr_hsig)%laction ) THEN 
     1170            &                                                       .OR. srcv(jpr_hsig)%laction ) THEN 
    11731171            CALL sbc_stokes() 
    11741172         ENDIF 
     
    15251523    
    15261524 
    1527    SUBROUTINE sbc_cpl_ice_flx( picefr, palbi, psst, pist ) 
     1525   SUBROUTINE sbc_cpl_ice_flx( picefr, palbi, psst, pist, phs, phi ) 
    15281526      !!---------------------------------------------------------------------- 
    15291527      !!             ***  ROUTINE sbc_cpl_ice_flx  *** 
     
    15761574      !!---------------------------------------------------------------------- 
    15771575      REAL(wp), INTENT(in), DIMENSION(:,:)             ::   picefr     ! ice fraction                [0 to 1] 
    1578       ! optional arguments, used only in 'mixed oce-ice' case 
     1576      !                                                !!           ! optional arguments, used only in 'mixed oce-ice' case 
    15791577      REAL(wp), INTENT(in), DIMENSION(:,:,:), OPTIONAL ::   palbi      ! all skies ice albedo  
    15801578      REAL(wp), INTENT(in), DIMENSION(:,:  ), OPTIONAL ::   psst       ! sea surface temperature     [Celsius] 
    15811579      REAL(wp), INTENT(in), DIMENSION(:,:,:), OPTIONAL ::   pist       ! ice surface temperature     [Kelvin] 
    1582       ! 
    1583       INTEGER ::   jl         ! dummy loop index 
    1584       REAL(wp), POINTER, DIMENSION(:,:  ) ::   zcptn, zcptrain, zcptsnw, ziceld, zmsk, zsnw 
    1585       REAL(wp), POINTER, DIMENSION(:,:  ) ::   zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap_oce, zevap_ice, zdevap_ice 
    1586       REAL(wp), POINTER, DIMENSION(:,:  ) ::   zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice 
    1587       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice 
     1580      REAL(wp), INTENT(in), DIMENSION(:,:,:), OPTIONAL ::   phs        ! snow depth                  [m] 
     1581      REAL(wp), INTENT(in), DIMENSION(:,:,:), OPTIONAL ::   phi        ! ice thickness               [m] 
     1582      ! 
     1583      INTEGER  ::   ji, jj, jl   ! dummy loop index 
     1584      REAL(wp) ::   ztri         ! local scalar 
     1585      REAL(wp), DIMENSION(jpi,jpj)     ::   zcptn, zcptrain, zcptsnw, ziceld, zmsk, zsnw 
     1586      REAL(wp), DIMENSION(jpi,jpj)     ::   zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip  , zevap_oce, zevap_ice, zdevap_ice 
     1587      REAL(wp), DIMENSION(jpi,jpj)     ::   zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice 
     1588      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice    !!gm , zfrqsr_tr_i 
    15881589      !!---------------------------------------------------------------------- 
    15891590      ! 
    15901591      IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_ice_flx') 
    15911592      ! 
    1592       CALL wrk_alloc( jpi,jpj,     zcptn, zcptrain, zcptsnw, ziceld, zmsk, zsnw ) 
    1593       CALL wrk_alloc( jpi,jpj,     zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap_oce, zevap_ice, zdevap_ice ) 
    1594       CALL wrk_alloc( jpi,jpj,     zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice ) 
    1595       CALL wrk_alloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice ) 
    1596  
    15971593      IF( ln_mixcpl )   zmsk(:,:) = 1. - xcplmask(:,:,0) 
    1598       ziceld(:,:) = 1. - picefr(:,:) 
    1599       zcptn(:,:) = rcp * sst_m(:,:) 
     1594      ziceld(:,:) = 1._wp - picefr(:,:) 
     1595      zcptn (:,:) = rcp * sst_m(:,:) 
    16001596      ! 
    16011597      !                                                      ! ========================= ! 
     
    16221618#if defined key_lim3 
    16231619      ! zsnw = snow fraction over ice after wind blowing (=picefr if no blowing) 
    1624       zsnw(:,:) = 0._wp  ;  CALL ice_thd_snwblow( ziceld, zsnw ) 
     1620      zsnw(:,:) = 0._wp   ;   CALL ice_thd_snwblow( ziceld, zsnw ) 
    16251621       
    16261622      ! --- evaporation minus precipitation corrected (because of wind blowing on snow) --- ! 
     
    16591655         sprecip(:,:) = sprecip(:,:) * xcplmask(:,:,0) + zsprecip(:,:) * zmsk(:,:) 
    16601656         tprecip(:,:) = tprecip(:,:) * xcplmask(:,:,0) + ztprecip(:,:) * zmsk(:,:) 
    1661          DO jl=1,jpl 
     1657         DO jl = 1, jpl 
    16621658            evap_ice (:,:,jl) = evap_ice (:,:,jl) * xcplmask(:,:,0) + zevap_ice (:,:) * zmsk(:,:) 
    16631659            devap_ice(:,:,jl) = devap_ice(:,:,jl) * xcplmask(:,:,0) + zdevap_ice(:,:) * zmsk(:,:) 
    1664          ENDDO 
     1660         END DO 
    16651661      ELSE 
    1666          emp_tot(:,:) =         zemp_tot(:,:) 
    1667          emp_ice(:,:) =         zemp_ice(:,:) 
    1668          emp_oce(:,:) =         zemp_oce(:,:)      
    1669          sprecip(:,:) =         zsprecip(:,:) 
    1670          tprecip(:,:) =         ztprecip(:,:) 
    1671          DO jl=1,jpl 
     1662         emp_tot(:,:) = zemp_tot(:,:) 
     1663         emp_ice(:,:) = zemp_ice(:,:) 
     1664         emp_oce(:,:) = zemp_oce(:,:)      
     1665         sprecip(:,:) = zsprecip(:,:) 
     1666         tprecip(:,:) = ztprecip(:,:) 
     1667         DO jl = 1, jpl 
    16721668            evap_ice (:,:,jl) = zevap_ice (:,:) 
    16731669            devap_ice(:,:,jl) = zdevap_ice(:,:) 
    1674          ENDDO 
     1670         END DO 
    16751671      ENDIF 
    16761672 
     
    16911687        fwfisf(:,:) = - frcv(jpr_isf)%z3(:,:,1) 
    16921688      ENDIF 
    1693  
     1689      ! 
    16941690      IF( ln_mixcpl ) THEN 
    16951691         emp_tot(:,:) = emp_tot(:,:) * xcplmask(:,:,0) + zemp_tot(:,:) * zmsk(:,:) 
     
    17031699         tprecip(:,:) =                                  ztprecip(:,:) 
    17041700      ENDIF 
    1705  
     1701      ! 
    17061702#endif 
     1703 
    17071704      ! outputs 
    17081705!!      IF( srcv(jpr_rnf)%laction )   CALL iom_put( 'runoffs' , rnf(:,:) * tmask(:,:,1)                                 )  ! runoff 
     
    17301727            zqns_ice(:,:,1:jpl) = frcv(jpr_qnsice)%z3(:,:,1:jpl) 
    17311728         ELSE 
    1732             DO jl=1,jpl 
     1729            DO jl = 1, jpl 
    17331730               zqns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,1) ! Set all category values equal 
    1734             ENDDO 
     1731            END DO 
    17351732         ENDIF 
    17361733      CASE( 'oce and ice' )      ! the total flux is computed from ocean and ice fluxes 
     
    17431740         ELSE 
    17441741            qns_tot(:,:) = qns_tot(:,:) + picefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1) 
    1745             DO jl=1,jpl 
     1742            DO jl = 1, jpl 
    17461743               zqns_tot(:,:   ) = zqns_tot(:,:) + picefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1) 
    17471744               zqns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,1) 
    1748             ENDDO 
     1745            END DO 
    17491746         ENDIF 
    17501747      CASE( 'mixed oce-ice' )    ! the ice flux is cumputed from the total flux, the SST and ice informations 
     
    17661763      !         note: ziceld cannot be = 0 since we limit the ice concentration to amax 
    17671764      zqns_oce = 0._wp 
    1768       WHERE( ziceld /= 0._wp )  zqns_oce(:,:) = ( zqns_tot(:,:) - SUM( a_i * zqns_ice, dim=3 ) ) / ziceld(:,:) 
     1765      WHERE( ziceld /= 0._wp )   zqns_oce(:,:) = ( zqns_tot(:,:) - SUM( a_i * zqns_ice, dim=3 ) ) / ziceld(:,:) 
    17691766 
    17701767      ! Heat content per unit mass of snow (J/kg) 
     
    18591856         ELSE 
    18601857            ! Set all category values equal for the moment 
    1861             DO jl=1,jpl 
     1858            DO jl = 1, jpl 
    18621859               zqsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,1) 
    1863             ENDDO 
     1860            END DO 
    18641861         ENDIF 
    18651862         zqsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1) 
     
    18681865         zqsr_tot(:,:  ) =  ziceld(:,:) * frcv(jpr_qsroce)%z3(:,:,1) 
    18691866         IF ( TRIM(sn_rcv_qsr%clcat) == 'yes' ) THEN 
    1870             DO jl=1,jpl 
     1867            DO jl = 1, jpl 
    18711868               zqsr_tot(:,:   ) = zqsr_tot(:,:) + a_i(:,:,jl) * frcv(jpr_qsrice)%z3(:,:,jl)    
    18721869               zqsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,jl) 
    1873             ENDDO 
     1870            END DO 
    18741871         ELSE 
    18751872            qsr_tot(:,:   ) = qsr_tot(:,:) + picefr(:,:) * frcv(jpr_qsrice)%z3(:,:,1) 
    1876             DO jl=1,jpl 
     1873            DO jl = 1, jpl 
    18771874               zqsr_tot(:,:   ) = zqsr_tot(:,:) + picefr(:,:) * frcv(jpr_qsrice)%z3(:,:,1) 
    18781875               zqsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,1) 
    1879             ENDDO 
     1876            END DO 
    18801877         ENDIF 
    18811878      CASE( 'mixed oce-ice' ) 
     
    18901887      IF( ln_dm2dc .AND. ln_cpl ) THEN   ! modify qsr to include the diurnal cycle 
    18911888         zqsr_tot(:,:  ) = sbc_dcy( zqsr_tot(:,:  ) ) 
    1892          DO jl=1,jpl 
     1889         DO jl = 1, jpl 
    18931890            zqsr_ice(:,:,jl) = sbc_dcy( zqsr_ice(:,:,jl) ) 
    1894          ENDDO 
     1891         END DO 
    18951892      ENDIF 
    18961893 
     
    19081905         qsr_tot(:,:) = qsr(:,:) * ziceld(:,:) + SUM( qsr_ice(:,:,:) * a_i(:,:,:), dim=3 )   ! total flux from blk 
    19091906         qsr_tot(:,:) = qsr_tot(:,:) * xcplmask(:,:,0) +  zqsr_tot(:,:)* zmsk(:,:) 
    1910          DO jl=1,jpl 
     1907         DO jl = 1, jpl 
    19111908            qsr_ice(:,:,jl) = qsr_ice(:,:,jl) * xcplmask(:,:,0) +  zqsr_ice(:,:,jl)* zmsk(:,:) 
    1912          ENDDO 
     1909         END DO 
    19131910      ELSE 
    19141911         qsr_tot(:,:  ) = zqsr_tot(:,:  ) 
     
    19461943      END SELECT 
    19471944 
    1948       ! Surface transimission parameter io (Maykut Untersteiner , 1971 ; Ebert and Curry, 1993 ) 
    1949       ! Used for LIM3 
    1950       ! Coupled case: since cloud cover is not received from atmosphere  
    1951       !               ===> used prescribed cloud fraction representative for polar oceans in summer (0.81) 
    1952       fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice ) 
    1953       fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) 
    1954  
    1955       CALL wrk_dealloc( jpi,jpj,     zcptn, zcptrain, zcptsnw, ziceld, zmsk, zsnw ) 
    1956       CALL wrk_dealloc( jpi,jpj,     zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap_oce, zevap_ice, zdevap_ice ) 
    1957       CALL wrk_dealloc( jpi,jpj,     zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice ) 
    1958       CALL wrk_dealloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice ) 
     1945      !                                                      ! ========================= ! 
     1946      !                                                      !      Transmitted Qsr      !   [W/m2] 
     1947      !                                                      ! ========================= ! 
     1948      SELECT CASE( nice_jules ) 
     1949      CASE( np_jules_OFF    )       !==  No Jules coupler  ==! 
     1950         ! 
     1951!!gm         ! former coding was     
     1952!!gm         ! Coupled case: since cloud cover is not received from atmosphere  
     1953!!gm         !               ===> used prescribed cloud fraction representative for polar oceans in summer (0.81) 
     1954!!gm         !     fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice ) 
     1955!!gm         !     fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) 
     1956!!gm          
     1957!!gm         ! to retrieve that coding, we needed to access h_i & h_s from here 
     1958!!gm         ! we could even retrieve cloud fraction from the coupler 
     1959!!gm         ! 
     1960!!gm         zfrqsr_tr_i(:,:,:) = 0._wp   !   surface transmission parameter 
     1961!!gm         ! 
     1962!!gm         DO jl = 1, jpl 
     1963!!gm            DO jj = 1 , jpj 
     1964!!gm               DO ji = 1, jpi 
     1965!!gm                  !              !--- surface transmission parameter (Grenfell Maykut 77) --- ! 
     1966!!gm                  zfrqsr_tr_i(ji,jj,jl) = 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice  
     1967!!gm                  ! 
     1968!!gm                  !              ! --- influence of snow and thin ice --- ! 
     1969!!gm                  IF ( phs(ji,jj,jl) >= 0.0_wp )   zfrqsr_tr_i(ji,jj,jl) = 0._wp   !   snow fully opaque 
     1970!!gm                  IF ( phi(ji,jj,jl) <= 0.1_wp )   zfrqsr_tr_i(ji,jj,jl) = 1._wp   !   thin ice transmits all solar radiation 
     1971!!gm               END DO 
     1972!!gm            END DO 
     1973!!gm         END DO 
     1974!!gm         ! 
     1975!!gm         qsr_ice_tr(:,:,:) =   zfrqsr_tr_i(:,:,:) * qsr_ice(:,:,:)               !   transmitted solar radiation  
     1976!!gm         ! 
     1977!!gm better coding of the above calculation: 
     1978         ! 
     1979         !                    ! ===> used prescribed cloud fraction representative for polar oceans in summer (0.81) 
     1980         ztri = 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice    ! surface transmission parameter (Grenfell Maykut 77) 
     1981         ! 
     1982         qsr_ice_tr(:,:,:) = ztri * qsr_ice(:,:,:) 
     1983         WHERE( phs(:,:,:) >= 0.0_wp )   qsr_ice_tr(:,:,:) = 0._wp            ! snow fully opaque 
     1984    &n