Changeset 8726


Ignore:
Timestamp:
2017-11-17T10:22:38+01:00 (3 years ago)
Author:
vancop
Message:

Jules coupler functionalities

Location:
branches/2017/dev_r8183_ICEMODEL/NEMOGCM
Files:
11 edited

Legend:

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

    r8680 r8726  
    122122                                      !     = 2  Redistribute a single flux over categories 
    123123                                      !          ==> coupled mode only 
     124   nice_jules       = 1               !  Jules coupling (0=OFF, 1=EMULATED, 2=ACTIVE) 
    124125/ 
    125126!------------------------------------------------------------------------------ 
     
    140141                                      !     Obs: 0.1-0.5 (Lecomte et al, JAMES 2013) 
    141142   rn_kappa_i       =   1.0           !  radiation attenuation coefficient in sea ice [1/m] 
    142    ln_dqns_i        = .true.          !  change the surface non-solar flux with surface temperature (T) or not (F) 
    143143/ 
    144144!------------------------------------------------------------------------------ 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice.F90

    r8678 r8726  
    196196   !                                      !   = 2  Redistribute a single flux over categories 
    197197 
     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) 
     203 
    198204   !                                     !!** ice-salinity namelist (namthd_sal) ** 
    199205   INTEGER , PUBLIC ::   nn_icesal        !: salinity configuration used in the model 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice1d.F90

    r8691 r8726  
    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),  & 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceforcing.F90

    r8597 r8726  
    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_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90

    r8598 r8726  
    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_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90

    r8626 r8726  
    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_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_dh.F90

    r8626 r8726  
    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_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_zdf.F90

    r8585 r8726  
    3434   LOGICAL  ::   ln_cndi_U64      ! thermal conductivity: Untersteiner (1964) 
    3535   LOGICAL  ::   ln_cndi_P07      ! thermal conductivity: Pringle et al (2007) 
    36    REAL(wp) ::   rn_cnd_s         ! thermal conductivity of the snow [W/m/K] 
    3736   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  
     37   REAL(wp), PUBLIC ::   rn_cnd_s ! thermal conductivity of the snow [W/m/K] 
     38 
     39   INTEGER  ::   nice_zdf         ! Choice of the type of vertical heat diffusion formulation 
     40   !                                           ! associated indices: 
     41   INTEGER, PARAMETER ::   np_BL99    = 1      ! Bitz and Lipscomb (1999) 
     42 
     43   INTEGER , PARAMETER ::   np_zdf_jules_OFF    = 0   !   compute all temperatures from qsr and qns 
     44   INTEGER , PARAMETER ::   np_zdf_jules_SND    = 1   !   compute conductive heat flux and surface temperature from qsr and qns 
     45   INTEGER , PARAMETER ::   np_zdf_jules_RCV    = 2   !   compute snow and inner ice temperatures from qcnd 
     46    
    4047   !!---------------------------------------------------------------------- 
    4148   !! NEMO/ICE 4.0 , NEMO Consortium (2017) 
     
    4552CONTAINS 
    4653 
    47    SUBROUTINE ice_thd_zdf 
     54      SUBROUTINE ice_thd_zdf 
     55       
     56      !! 
    4857      !!------------------------------------------------------------------- 
    4958      !!                ***  ROUTINE ice_thd_zdf  *** 
    5059      !! ** Purpose : 
     60      !!           This chooses between the appropriate routine for the  
     61      !!           computation of vertical diffusion 
     62      !! 
     63      !!------------------------------------------------------------------- 
     64      !! 
     65    
     66      SELECT CASE ( nice_zdf )      ! Choose the vertical heat diffusion solver 
     67       
     68                                       !-------------       
     69         CASE( np_BL99 )               ! BL99 solver 
     70                                       !------------- 
     71          
     72            IF ( nice_jules == np_jules_OFF ) THEN   ! No Jules coupler      
     73 
     74               CALL ice_thd_zdf_BL99 ( np_zdf_jules_OFF ) 
     75                                        
     76            ENDIF 
     77             
     78            IF ( nice_jules == np_jules_EMULE ) THEN   ! Jules coupler is emulated 
     79             
     80               CALL ice_thd_zdf_BL99 ( np_zdf_jules_SND ) 
     81               CALL ice_thd_zdf_BL99 ( np_zdf_jules_RCV ) 
     82                                                    
     83            ENDIF 
     84 
     85            IF ( nice_jules == np_jules_ACTIVE ) THEN   ! Jules coupler is emulated 
     86             
     87               CALL ice_thd_zdf_BL99 ( np_zdf_jules_RCV ) 
     88                                                    
     89            ENDIF 
     90             
     91      END SELECT 
     92       
     93   END SUBROUTINE ice_thd_zdf 
     94    
     95    
     96    
     97   SUBROUTINE ice_thd_zdf_BL99(k_jules) 
     98      !!------------------------------------------------------------------- 
     99      !!                ***  ROUTINE ice_thd_zdf_BL99  *** 
     100      !! ** Purpose : 
    51101      !!           This routine determines the time evolution of snow and sea-ice  
    52       !!           temperature profiles. 
     102      !!           temperature profiles, using the original Bitz and Lipscomb (1999) algorithm 
     103      !! 
    53104      !! ** Method  : 
    54105      !!           This is done by solving the heat equation diffusion with 
     
    83134      !!           total ice/snow thickness : h_i_1d, h_s_1d 
    84135      !!------------------------------------------------------------------- 
     136      INTEGER, INTENT(in) ::   k_jules     ! Jules coupling (0=OFF, 1=RECEIVE, 2=SEND) 
     137       
    85138      INTEGER ::   ji, jk         ! spatial loop index 
    86139      INTEGER ::   jm             ! current reference number of equation 
     
    101154      REAL(wp) ::   zdti_bnd  =  1.e-4_wp     ! maximal authorized error on temperature  
    102155      REAL(wp) ::   ztmelt_i                  ! ice melting temperature 
    103       REAL(wp) ::   z1_hsu 
    104156      REAL(wp) ::   zdti_max                  ! current maximal error on temperature  
    105157      REAL(wp) ::   zcpi                      ! Ice specific heat 
     
    111163      REAL(wp), DIMENSION(jpij)     ::   zh_i, z1_h_i ! ice layer thickness 
    112164      REAL(wp), DIMENSION(jpij)     ::   zh_s, z1_h_s ! snow layer thickness 
    113       REAL(wp), DIMENSION(jpij)     ::   zfsw        ! solar radiation absorbed at the surface 
    114165      REAL(wp), DIMENSION(jpij)     ::   zqns_ice_b  ! solar radiation absorbed at the surface 
    115166      REAL(wp), DIMENSION(jpij)     ::   zfnet       ! surface flux function 
    116167      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        
     168       
     169      REAL(wp), DIMENSION(jpij       )     ::   ztsuold     ! Old surface temperature in the ice 
    119170      REAL(wp), DIMENSION(jpij,nlay_i)     ::   ztiold      ! Old temperature in the ice 
    120171      REAL(wp), DIMENSION(jpij,nlay_s)     ::   ztsold      ! Old temperature in the snow 
     
    136187      REAL(wp), DIMENSION(jpij)            ::   zq_ini      ! diag errors on heat 
    137188      REAL(wp), DIMENSION(jpij)            ::   zghe        ! G(he), th. conduct enhancement factor, mono-cat 
     189 
     190      REAL(wp) ::   zfr1, zfr2, zfrqsr_tr_i   ! dummy factor 
    138191       
    139192      ! Mono-category 
     
    166219      ELSEWHERE                         ;   z1_h_s(1:npti) = 0._wp 
    167220      END WHERE 
    168       ! 
    169       ! temperatures 
    170       ztsub  (1:npti)   = t_su_1d(1:npti)  ! temperature at the previous iteration 
     221 
     222      ! Store initial temperatures and non solar heat fluxes 
     223      IF ( k_jules == np_zdf_jules_OFF .OR. k_jules == np_zdf_jules_SND  ) THEN ! OFF or SND mode 
     224 
     225         ztsub  (1:npti)       =   t_su_1d(1:npti)                          ! surface temperature at iteration n-1 
     226         ztsuold(1:npti)       =   t_su_1d(1:npti)                          ! surface temperature initial value 
     227         zdqns_ice_b(1:npti)   =   dqns_ice_1d(1:npti)                      ! derivative of incoming nonsolar flux  
     228         zqns_ice_b (1:npti)   =   qns_ice_1d(1:npti)                       ! store previous qns_ice_1d value 
     229 
     230         t_su_1d(1:npti)       =   MIN( t_su_1d(1:npti), rt0 - ztsu_err )   ! required to leave the choice between melting or not 
     231 
     232      ENDIF    
     233 
    171234      ztsold (1:npti,:) = t_s_1d(1:npti,:)   ! Old snow temperature 
    172235      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       ! 
     236 
    175237      !------------- 
    176238      ! 2) Radiation 
    177239      !------------- 
    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  
    200240      ! --- Transmission/absorption of solar radiation in the ice --- ! 
    201       zradtr_s(1:npti,0) = zftrice(1:npti) 
     241!     zfr1 = ( 0.18 * ( 1.0 - 0.81 ) + 0.35 * 0.81 )            !   standard value 
     242!     zfr2 = ( 0.82 * ( 1.0 - 0.81 ) + 0.65 * 0.81 )            !   zfr2 such that zfr1 + zfr2 to equal 1 
     243 
     244!     DO ji = 1, npti 
     245 
     246!        zfac = MAX( 0._wp , 1._wp - ( h_i_1d(ji) * 10._wp ) )      
     247 
     248!              zfrqsr_tr_i = zfr1 + zfac * zfr2                         !   below 10 cm, linearly increase zfrqsr_tr_i until 1 at zero thickness 
     249!              IF ( h_s_1d(ji) >= 0.0_wp )   zfrqsr_tr_i = 0._wp     !   snow fully opaque 
     250 
     251!              qsr_ice_tr_1d(ji) = zfrqsr_tr_i * qsr_ice_1d(ji)   !   transmitted solar radiation  
     252 
     253!              zfsw(ji) = qsr_ice_1d(ji) - qsr_ice_tr_1d(ji) 
     254!              zftrice(ji) = qsr_ice_tr_1d(ji) 
     255!              i0(ji) = zfrqsr_tr_i 
     256 
     257!     END DO 
     258 
     259      zradtr_s(1:npti,0) = qsr_ice_tr_1d(1:npti) 
    202260      DO jk = 1, nlay_s 
    203261         DO ji = 1, npti 
     
    209267      END DO 
    210268          
    211       zradtr_i(1:npti,0) = zradtr_s(1:npti,nlay_s) * isnow(1:npti) + zftrice(1:npti) * ( 1._wp - isnow(1:npti) ) 
     269      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) ) 
    212270      DO jk = 1, nlay_i  
    213271         DO ji = 1, npti 
     
    330388            END DO 
    331389         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 
     390          
     391         ! 
     392         !----------------------------------------! 
     393         !                                        ! 
     394         !       JULES IF (OFF or SND MODE)       ! 
     395         !                                        ! 
     396         !----------------------------------------! 
     397         ! 
     398          
     399         IF ( k_jules == np_zdf_jules_OFF .OR. k_jules == np_zdf_jules_SND  ) THEN ! OFF or SND mode 
     400          
     401            ! 
     402            ! In OFF mode the original BL99 temperature computation is used 
     403            ! (with qsr_ice, qns_ice and dqns_ice as inputs) 
     404            ! 
     405            ! In SND mode, the computation is required to compute the conduction fluxes 
     406            ! 
     407          
     408            !---------------------------- 
     409            ! 6) surface flux computation 
     410            !---------------------------- 
     411 
     412            DO ji = 1, npti 
     413            ! update of the non solar flux according to the update in T_su 
    339414               qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) ) 
    340415            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  
     416 
     417            DO ji = 1, npti 
     418               zfnet(ji) = qsr_ice_1d(ji) - qsr_ice_tr_1d(ji) + qns_ice_1d(ji) ! net heat flux = net solar - transmitted solar + non solar 
     419            END DO 
     420            ! 
     421            !---------------------------- 
     422            ! 7) tridiagonal system terms 
     423            !---------------------------- 
     424            !!layer denotes the number of the layer in the snow or in the ice 
     425            !!jm denotes the reference number of the equation in the tridiagonal 
     426            !!system, terms of tridiagonal system are indexed as following : 
     427            !!1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one 
     428 
     429            !!ice interior terms (top equation has the same form as the others) 
     430            ztrid   (1:npti,:,:) = 0._wp 
     431            zindterm(1:npti,:)   = 0._wp 
     432            zindtbis(1:npti,:)   = 0._wp 
     433            zdiagbis(1:npti,:)   = 0._wp 
     434 
     435            DO jm = nlay_s + 2, nlay_s + nlay_i  
    362436            DO ji = 1, npti 
    363437               jk = jm - nlay_s - 1 
     
    367441               zindterm(ji,jm) = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk) 
    368442            END DO 
    369          ENDDO 
    370  
    371          jm =  nlay_s + nlay_i + 1 
    372          DO ji = 1, npti 
     443            ENDDO 
     444 
     445            jm =  nlay_s + nlay_i + 1 
     446            DO ji = 1, npti 
    373447            !!ice bottom term 
    374448            ztrid(ji,jm,1)  = - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1)    
     
    377451            zindterm(ji,jm) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) *  & 
    378452               &            ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 *  t_bo_1d(ji) )  
    379          ENDDO 
    380  
    381  
    382          DO ji = 1, npti 
     453            ENDDO 
     454 
     455 
     456            DO ji = 1, npti 
    383457            !                               !---------------------! 
    384             IF ( h_s_1d(ji) > 0.0 ) THEN   !  snow-covered cells ! 
     458            IF ( h_s_1d(ji) > 0.0 ) THEN    !  snow-covered cells ! 
    385459               !                            !---------------------! 
    386460               ! snow interior terms (bottom equation has the same form as the others) 
    387461               DO jm = 3, nlay_s + 1 
    388                   jk = jm - 1 
    389                   ztrid(ji,jm,1)  = - zeta_s(ji,jk) * zkappa_s(ji,jk-1) 
    390                   ztrid(ji,jm,2)  = 1.0 + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) ) 
    391                   ztrid(ji,jm,3)  = - zeta_s(ji,jk)*zkappa_s(ji,jk) 
    392                   zindterm(ji,jm) = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk) 
     462               jk = jm - 1 
     463               ztrid(ji,jm,1)  = - zeta_s(ji,jk) * zkappa_s(ji,jk-1) 
     464               ztrid(ji,jm,2)  = 1.0 + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) ) 
     465               ztrid(ji,jm,3)  = - zeta_s(ji,jk)*zkappa_s(ji,jk) 
     466               zindterm(ji,jm) = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk) 
    393467               END DO 
    394468 
    395469               ! case of only one layer in the ice (ice equation is altered) 
    396470               IF ( nlay_i == 1 ) THEN 
    397                   ztrid(ji,nlay_s+2,3)  = 0.0 
    398                   zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1) * t_bo_1d(ji)  
     471               ztrid(ji,nlay_s+2,3)  = 0.0 
     472               zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1) * t_bo_1d(ji)  
    399473               ENDIF 
    400474 
    401475               IF ( t_su_1d(ji) < rt0 ) THEN   !--  case 1 : no surface melting 
    402476 
    403                   jm_min(ji) = 1 
    404                   jm_max(ji) = nlay_i + nlay_s + 1 
    405  
    406                   ! surface equation 
    407                   ztrid(ji,1,1)  = 0.0 
    408                   ztrid(ji,1,2)  = zdqns_ice_b(ji) - zg1s * zkappa_s(ji,0) 
    409                   ztrid(ji,1,3)  = zg1s * zkappa_s(ji,0) 
    410                   zindterm(ji,1) = zdqns_ice_b(ji) * t_su_1d(ji) - zfnet(ji) 
    411  
    412                   ! first layer of snow equation 
    413                   ztrid(ji,2,1)  = - zkappa_s(ji,0) * zg1s * zeta_s(ji,1) 
    414                   ztrid(ji,2,2)  = 1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s ) 
    415                   ztrid(ji,2,3)  = - zeta_s(ji,1)* zkappa_s(ji,1) 
    416                   zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * zradab_s(ji,1) 
     477               jm_min(ji) = 1 
     478               jm_max(ji) = nlay_i + nlay_s + 1 
     479 
     480               ! surface equation 
     481               ztrid(ji,1,1)  = 0.0 
     482               ztrid(ji,1,2)  = zdqns_ice_b(ji) - zg1s * zkappa_s(ji,0) 
     483               ztrid(ji,1,3)  = zg1s * zkappa_s(ji,0) 
     484               zindterm(ji,1) = zdqns_ice_b(ji) * t_su_1d(ji) - zfnet(ji) 
     485 
     486               ! first layer of snow equation 
     487               ztrid(ji,2,1)  = - zkappa_s(ji,0) * zg1s * zeta_s(ji,1) 
     488               ztrid(ji,2,2)  = 1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s ) 
     489               ztrid(ji,2,3)  = - zeta_s(ji,1)* zkappa_s(ji,1) 
     490               zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * zradab_s(ji,1) 
    417491 
    418492               ELSE                            !--  case 2 : surface is melting 
    419                   ! 
    420                   jm_min(ji) = 2 
    421                   jm_max(ji) = nlay_i + nlay_s + 1 
    422  
    423                   ! first layer of snow equation 
    424                   ztrid(ji,2,1)  = 0.0 
    425                   ztrid(ji,2,2)  = 1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s ) 
    426                   ztrid(ji,2,3)  = - zeta_s(ji,1)*zkappa_s(ji,1)  
    427                   zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) *  & 
    428                      &           ( zradab_s(ji,1) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) )  
     493               ! 
     494               jm_min(ji) = 2 
     495               jm_max(ji) = nlay_i + nlay_s + 1 
     496 
     497               ! first layer of snow equation 
     498               ztrid(ji,2,1)  = 0.0 
     499               ztrid(ji,2,2)  = 1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s ) 
     500               ztrid(ji,2,3)  = - zeta_s(ji,1)*zkappa_s(ji,1)  
     501               zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) *  & 
     502               &           ( zradab_s(ji,1) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) )  
    429503               ENDIF 
    430504            !                               !---------------------! 
     
    433507               ! 
    434508               IF ( t_su_1d(ji) < rt0 ) THEN   !--  case 1 : no surface melting 
    435                   ! 
    436                   jm_min(ji) = nlay_s + 1 
    437                   jm_max(ji) = nlay_i + nlay_s + 1 
    438  
    439                   ! surface equation    
    440                   ztrid(ji,jm_min(ji),1)  = 0.0 
    441                   ztrid(ji,jm_min(ji),2)  = zdqns_ice_b(ji) - zkappa_i(ji,0)*zg1     
    442                   ztrid(ji,jm_min(ji),3)  = zkappa_i(ji,0)*zg1 
    443                   zindterm(ji,jm_min(ji)) = zdqns_ice_b(ji)*t_su_1d(ji) - zfnet(ji) 
    444  
    445                   ! first layer of ice equation 
    446                   ztrid(ji,jm_min(ji)+1,1)  = - zkappa_i(ji,0) * zg1 * zeta_i(ji,1) 
    447                   ztrid(ji,jm_min(ji)+1,2)  = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 ) 
    448                   ztrid(ji,jm_min(ji)+1,3)  = - zeta_i(ji,1) * zkappa_i(ji,1)   
    449                   zindterm(ji,jm_min(ji)+1) = ztiold(ji,1) + zeta_i(ji,1) * zradab_i(ji,1)   
    450  
    451                   ! case of only one layer in the ice (surface & ice equations are altered) 
    452                   IF ( nlay_i == 1 ) THEN 
    453                      ztrid(ji,jm_min(ji),1)    = 0.0 
    454                      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 
    456                      ztrid(ji,jm_min(ji)+1,1)  = -zkappa_i(ji,0) * 2.0 * zeta_i(ji,1) 
    457                      ztrid(ji,jm_min(ji)+1,2)  = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) ) 
    458                      ztrid(ji,jm_min(ji)+1,3)  = 0.0 
    459                      zindterm(ji,jm_min(ji)+1) = ztiold(ji,1) + zeta_i(ji,1) *  & 
    460                         &                      ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) ) 
    461                   ENDIF 
     509               ! 
     510               jm_min(ji) = nlay_s + 1 
     511               jm_max(ji) = nlay_i + nlay_s + 1 
     512 
     513               ! surface equation    
     514               ztrid(ji,jm_min(ji),1)  = 0.0 
     515               ztrid(ji,jm_min(ji),2)  = zdqns_ice_b(ji) - zkappa_i(ji,0)*zg1     
     516               ztrid(ji,jm_min(ji),3)  = zkappa_i(ji,0)*zg1 
     517               zindterm(ji,jm_min(ji)) = zdqns_ice_b(ji)*t_su_1d(ji) - zfnet(ji) 
     518 
     519               ! first layer of ice equation 
     520               ztrid(ji,jm_min(ji)+1,1)  = - zkappa_i(ji,0) * zg1 * zeta_i(ji,1) 
     521               ztrid(ji,jm_min(ji)+1,2)  = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 ) 
     522               ztrid(ji,jm_min(ji)+1,3)  = - zeta_i(ji,1) * zkappa_i(ji,1)   
     523               zindterm(ji,jm_min(ji)+1) = ztiold(ji,1) + zeta_i(ji,1) * zradab_i(ji,1)   
     524 
     525               ! case of only one layer in the ice (surface & ice equations are altered) 
     526               IF ( nlay_i == 1 ) THEN 
     527               ztrid(ji,jm_min(ji),1)    = 0.0 
     528               ztrid(ji,jm_min(ji),2)    = zdqns_ice_b(ji) - zkappa_i(ji,0) * 2.0 
     529               ztrid(ji,jm_min(ji),3)    = zkappa_i(ji,0) * 2.0 
     530               ztrid(ji,jm_min(ji)+1,1)  = -zkappa_i(ji,0) * 2.0 * zeta_i(ji,1) 
     531               ztrid(ji,jm_min(ji)+1,2)  = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) ) 
     532               ztrid(ji,jm_min(ji)+1,3)  = 0.0 
     533               zindterm(ji,jm_min(ji)+1) = ztiold(ji,1) + zeta_i(ji,1) *  & 
     534               &                      ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) ) 
     535               ENDIF 
    462536 
    463537               ELSE                            !--  case 2 : surface is melting 
    464538 
    465                   jm_min(ji)    =  nlay_s + 2 
    466                   jm_max(ji)    =  nlay_i + nlay_s + 1 
    467  
    468                   ! first layer of ice equation 
    469                   ztrid(ji,jm_min(ji),1)  = 0.0 
    470                   ztrid(ji,jm_min(ji),2)  = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 )   
    471                   ztrid(ji,jm_min(ji),3)  = - zeta_i(ji,1) * zkappa_i(ji,1) 
    472                   zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) *  & 
    473                      &                    ( zradab_i(ji,1) + zkappa_i(ji,0) * zg1 * t_su_1d(ji) )  
    474  
    475                   ! case of only one layer in the ice (surface & ice equations are altered) 
    476                   IF ( nlay_i == 1 ) THEN 
    477                      ztrid(ji,jm_min(ji),1)  = 0.0 
    478                      ztrid(ji,jm_min(ji),2)  = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) ) 
    479                      ztrid(ji,jm_min(ji),3)  = 0.0 
    480                      zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) )  & 
    481                         &                    + t_su_1d(ji) * zeta_i(ji,1) * zkappa_i(ji,0) * 2.0 
    482                   ENDIF 
     539               jm_min(ji)    =  nlay_s + 2 
     540               jm_max(ji)    =  nlay_i + nlay_s + 1 
     541 
     542               ! first layer of ice equation 
     543               ztrid(ji,jm_min(ji),1)  = 0.0 
     544               ztrid(ji,jm_min(ji),2)  = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 )   
     545               ztrid(ji,jm_min(ji),3)  = - zeta_i(ji,1) * zkappa_i(ji,1) 
     546               zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) *  & 
     547               &                    ( zradab_i(ji,1) + zkappa_i(ji,0) * zg1 * t_su_1d(ji) )  
     548 
     549               ! case of only one layer in the ice (surface & ice equations are altered) 
     550               IF ( nlay_i == 1 ) THEN 
     551               ztrid(ji,jm_min(ji),1)  = 0.0 
     552               ztrid(ji,jm_min(ji),2)  = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) ) 
     553               ztrid(ji,jm_min(ji),3)  = 0.0 
     554               zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) )  & 
     555               &                    + t_su_1d(ji) * zeta_i(ji,1) * zkappa_i(ji,0) * 2.0 
     556               ENDIF 
    483557 
    484558               ENDIF 
     
    488562            zdiagbis(ji,jm_min(ji)) = ztrid(ji,jm_min(ji),2) 
    489563            ! 
    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 
     564            END DO 
     565            ! 
     566            !------------------------------ 
     567            ! 8) tridiagonal system solving 
     568            !------------------------------ 
     569            ! Solve the tridiagonal system with Gauss elimination method. 
     570            ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984 
     571            jm_maxt = 0 
     572            jm_mint = nlay_i+5 
     573            DO ji = 1, npti 
    500574            jm_mint = MIN(jm_min(ji),jm_mint) 
    501575            jm_maxt = MAX(jm_max(ji),jm_maxt) 
    502          END DO 
    503  
    504          DO jk = jm_mint+1, jm_maxt 
     576            END DO 
     577 
     578            DO jk = jm_mint+1, jm_maxt 
    505579            DO ji = 1, npti 
    506580               jm = min(max(jm_min(ji)+1,jk),jm_max(ji)) 
     
    508582               zindtbis(ji,jm) = zindterm(ji,jm) - ztrid(ji,jm,1) * zindtbis(ji,jm-1) / zdiagbis(ji,jm-1) 
    509583            END DO 
    510          END DO 
    511  
    512          DO ji = 1, npti 
     584            END DO 
     585 
     586            DO ji = 1, npti 
    513587            ! ice temperatures 
    514588            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 
     589            END DO 
     590 
     591            DO jm = nlay_i + nlay_s, nlay_s + 2, -1 
    518592            DO ji = 1, npti 
    519593               jk    =  jm - nlay_s - 1 
    520594               t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm) 
    521595            END DO 
    522          END DO 
    523  
    524          DO ji = 1, npti 
     596            END DO 
     597 
     598            DO ji = 1, npti 
    525599            ! snow temperatures       
    526600            IF( h_s_1d(ji) > 0._wp ) THEN 
    527601               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) 
     602               &                / zdiagbis(ji,nlay_s+1) 
    529603            ENDIF 
    530604            ! surface temperature 
     
    532606            IF( t_su_1d(ji) < rt0 ) THEN 
    533607               t_su_1d(ji) = ( zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) *  & 
    534                   &          ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji)) 
    535             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 
     608               &          ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji)) 
     609            ENDIF 
     610            END DO 
     611            ! 
     612            !-------------------------------------------------------------- 
     613            ! 9) Has the scheme converged ?, end of the iterative procedure 
     614            !-------------------------------------------------------------- 
     615            ! check that nowhere it has started to melt 
     616            ! zdti_max is a measure of error, it has to be under zdti_bnd 
     617            zdti_max = 0._wp 
     618            DO ji = 1, npti 
    545619            t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , rt0 - 100._wp ) 
    546620            zdti_max = MAX( zdti_max, ABS( t_su_1d(ji) - ztsub(ji) ) )      
    547          END DO 
    548  
    549          DO jk  =  1, nlay_s 
     621            END DO 
     622 
     623            DO jk  =  1, nlay_s 
    550624            DO ji = 1, npti 
    551625               t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp ) 
    552626               zdti_max = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) ) 
    553627            END DO 
    554          END DO 
    555  
    556          DO jk  =  1, nlay_i 
     628            END DO 
     629 
     630            DO jk  =  1, nlay_i 
    557631            DO ji = 1, npti 
    558632               ztmelt_i      = -tmut * sz_i_1d(ji,jk) + rt0  
     
    560634               zdti_max      =  MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) 
    561635            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  
     636            END DO 
     637 
     638            ! Compute spatial maximum over all errors 
     639            ! note that this could be optimized substantially by iterating only the non-converging points 
     640            IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice ) 
     641         ! 
     642         !----------------------------------------! 
     643         !                                        ! 
     644         !       JULES IF (RCV MODE)              ! 
     645         !                                        ! 
     646         !----------------------------------------! 
     647         ! 
     648 
     649         ELSE IF ( k_jules == np_zdf_jules_RCV  ) THEN ! RCV mode 
     650          
     651            ! 
     652            ! In RCV mode, we use a modified BL99 solver  
     653            ! with conduction flux (qcn_ice) as forcing term 
     654            ! 
     655            !---------------------------- 
     656            ! 7) tridiagonal system terms 
     657            !---------------------------- 
     658            !!layer denotes the number of the layer in the snow or in the ice 
     659            !!jm denotes the reference number of the equation in the tridiagonal 
     660            !!system, terms of tridiagonal system are indexed as following : 
     661            !!1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one 
     662 
     663            !!ice interior terms (top equation has the same form as the others) 
     664            ztrid   (1:npti,:,:) = 0._wp 
     665            zindterm(1:npti,:)   = 0._wp 
     666            zindtbis(1:npti,:)   = 0._wp 
     667            zdiagbis(1:npti,:)   = 0._wp 
     668 
     669            DO jm = nlay_s + 2, nlay_s + nlay_i  
     670            DO ji = 1, npti 
     671               jk = jm - nlay_s - 1 
     672               ztrid(ji,jm,1)  = - zeta_i(ji,jk) * zkappa_i(ji,jk-1) 
     673               ztrid(ji,jm,2)  = 1.0 + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) ) 
     674               ztrid(ji,jm,3)  = - zeta_i(ji,jk) * zkappa_i(ji,jk) 
     675               zindterm(ji,jm) = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk) 
     676            END DO 
     677            ENDDO 
     678 
     679            jm =  nlay_s + nlay_i + 1 
     680            DO ji = 1, npti 
     681            !!ice bottom term 
     682            ztrid(ji,jm,1)  = - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1)    
     683            ztrid(ji,jm,2)  = 1.0 + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i) * zg1 + zkappa_i(ji,nlay_i-1) ) 
     684            ztrid(ji,jm,3)  = 0.0 
     685            zindterm(ji,jm) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) *  & 
     686               &            ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 *  t_bo_1d(ji) )  
     687            ENDDO 
     688 
     689 
     690            DO ji = 1, npti 
     691            !                              !---------------------! 
     692            IF ( h_s_1d(ji) > 0.0 ) THEN   !  snow-covered cells ! 
     693               !                           !---------------------! 
     694               ! snow interior terms (bottom equation has the same form as the others) 
     695               DO jm = 3, nlay_s + 1 
     696               jk = jm - 1 
     697               ztrid(ji,jm,1)  = - zeta_s(ji,jk) * zkappa_s(ji,jk-1) 
     698               ztrid(ji,jm,2)  = 1.0 + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) ) 
     699               ztrid(ji,jm,3)  = - zeta_s(ji,jk)*zkappa_s(ji,jk) 
     700               zindterm(ji,jm) = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk) 
     701               END DO 
     702 
     703               ! case of only one layer in the ice (ice equation is altered) 
     704               IF ( nlay_i == 1 ) THEN 
     705               ztrid(ji,nlay_s+2,3)  = 0.0 
     706               zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1) * t_bo_1d(ji)  
     707               ENDIF 
     708 
     709               jm_min(ji) = 2 
     710               jm_max(ji) = nlay_i + nlay_s + 1 
     711 
     712               ! first layer of snow equation 
     713               ztrid(ji,2,1)  = 0.0 
     714               ztrid(ji,2,2)  = 1.0 + zeta_s(ji,1) * zkappa_s(ji,1) 
     715               ztrid(ji,2,3)  = - zeta_s(ji,1)*zkappa_s(ji,1)  
     716               zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) *  & 
     717               &           ( zradab_s(ji,1) + qcn_ice_1d(ji) )  
     718 
     719            !                               !---------------------! 
     720            ELSE                            ! cells without snow  ! 
     721            !                               !---------------------! 
     722 
     723               jm_min(ji)    =  nlay_s + 2 
     724               jm_max(ji)    =  nlay_i + nlay_s + 1 
     725 
     726               ! first layer of ice equation 
     727               ztrid(ji,jm_min(ji),1)  = 0.0 
     728               ztrid(ji,jm_min(ji),2)  = 1.0 + zeta_i(ji,1) * zkappa_i(ji,1)   
     729               ztrid(ji,jm_min(ji),3)  = - zeta_i(ji,1) * zkappa_i(ji,1) 
     730               zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) *  & 
     731               &                    ( zradab_i(ji,1) + qcn_ice_1d(ji) ) 
     732 
     733               ! case of only one layer in the ice (surface & ice equations are altered) 
     734               IF ( nlay_i == 1 ) THEN 
     735               ztrid(ji,jm_min(ji),1)  = 0.0 
     736               ztrid(ji,jm_min(ji),2)  = 1.0 + zeta_i(ji,1) * zkappa_i(ji,1) 
     737               ztrid(ji,jm_min(ji),3)  = 0.0 
     738               zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji)    & 
     739               &                    + qcn_ice_1d(ji) ) 
     740 
     741               ENDIF 
     742 
     743            ENDIF 
     744            ! 
     745            zindtbis(ji,jm_min(ji)) = zindterm(ji,jm_min(ji)) 
     746            zdiagbis(ji,jm_min(ji)) = ztrid(ji,jm_min(ji),2) 
     747            ! 
     748            END DO 
     749            ! 
     750            !------------------------------ 
     751            ! 8) tridiagonal system solving 
     752            !------------------------------ 
     753            ! Solve the tridiagonal system with Gauss elimination method. 
     754            ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984 
     755            jm_maxt = 0 
     756            jm_mint = nlay_i+5 
     757            DO ji = 1, npti 
     758            jm_mint = MIN(jm_min(ji),jm_mint) 
     759            jm_maxt = MAX(jm_max(ji),jm_maxt) 
     760            END DO 
     761 
     762            DO jk = jm_mint+1, jm_maxt 
     763            DO ji = 1, npti 
     764               jm = min(max(jm_min(ji)+1,jk),jm_max(ji)) 
     765               zdiagbis(ji,jm) = ztrid(ji,jm,2)  - ztrid(ji,jm,1) * ztrid(ji,jm-1,3)  / zdiagbis(ji,jm-1) 
     766               zindtbis(ji,jm) = zindterm(ji,jm) - ztrid(ji,jm,1) * zindtbis(ji,jm-1) / zdiagbis(ji,jm-1) 
     767            END DO 
     768            END DO 
     769 
     770            DO ji = 1, npti 
     771            ! ice temperatures 
     772            t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji)) 
     773            END DO 
     774 
     775            DO jm = nlay_i + nlay_s, nlay_s + 2, -1 
     776            DO ji = 1, npti 
     777               jk    =  jm - nlay_s - 1 
     778               t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm) 
     779            END DO 
     780            END DO 
     781 
     782            DO ji = 1, npti 
     783            ! snow temperatures       
     784            IF( h_s_1d(ji) > 0._wp ) THEN 
     785               t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) )  & 
     786               &                / zdiagbis(ji,nlay_s+1) 
     787            ENDIF 
     788            END DO 
     789            ! 
     790            !-------------------------------------------------------------- 
     791            ! 9) Has the scheme converged ?, end of the iterative procedure 
     792            !-------------------------------------------------------------- 
     793            ! check that nowhere it has started to melt 
     794            ! zdti_max is a measure of error, it has to be under zdti_bnd 
     795            zdti_max = 0._wp 
     796 
     797            DO jk  =  1, nlay_s 
     798            DO ji = 1, npti 
     799               t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp ) 
     800               zdti_max = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) ) 
     801            END DO 
     802            END DO 
     803 
     804            DO jk  =  1, nlay_i 
     805            DO ji = 1, npti 
     806               ztmelt_i      = -tmut * sz_i_1d(ji,jk) + rt0  
     807               t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), rt0 - 100._wp ) 
     808               zdti_max      =  MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) 
     809            END DO 
     810            END DO 
     811 
     812            ! Compute spatial maximum over all errors 
     813            ! note that this could be optimized substantially by iterating only the non-converging points 
     814            IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice ) 
     815 
     816         ENDIF ! k_jules 
     817          
    568818      END DO  ! End of the do while iterative procedure 
    569  
     819       
    570820      IF( ln_icectl .AND. lwp ) THEN 
    571821         WRITE(numout,*) ' zdti_max : ', zdti_max 
    572822         WRITE(numout,*) ' iconv    : ', iconv 
    573823      ENDIF 
     824       
    574825      ! 
    575826      !----------------------------- 
    576827      ! 10) Fluxes at the interfaces 
    577828      !----------------------------- 
     829      ! 
     830      ! --- update conduction fluxes 
     831      ! 
    578832      DO ji = 1, npti 
    579          !                                ! surface ice conduction flux 
     833      !                                ! surface ice conduction flux 
    580834         fc_su(ji)   =  -           isnow(ji)   * zkappa_s(ji,0) * zg1s * (t_s_1d(ji,1) - t_su_1d(ji))   & 
    581835            &           - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1  * (t_i_1d(ji,1) - t_su_1d(ji)) 
    582          !                                ! bottom ice conduction flux 
     836      !                                ! bottom ice conduction flux 
    583837         fc_bo_i(ji) =  - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_1d(ji) - t_i_1d(ji,nlay_i)) ) 
    584838      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 
     839       
     840      ! 
     841      ! --- Diagnose the heat loss due to changing non-solar / conduction flux --- ! 
     842      ! 
     843      DO ji = 1, npti 
     844            IF ( k_jules == np_zdf_jules_OFF .OR. k_jules == np_zdf_jules_SND  ) THEN  
     845                   ! OFF or SND mode 
     846               hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qns_ice_1d(ji)  - zqns_ice_b(ji) ) * a_i_1d(ji)  
     847            ELSE   ! RCV mode 
     848               hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( fc_su(ji) - qcn_ice_1d(ji) ) * a_i_1d(ji)  
     849            ENDIF 
     850      END DO 
     851       
     852      ! 
     853      ! --- Diagnose the heat loss due to non-fully converged temperature solution (should not be above 10-4 W-m2) --- ! 
     854      ! 
     855 
     856      IF ( ( k_jules == np_zdf_jules_OFF ) .OR. ( k_jules == np_zdf_jules_RCV ) ) THEN ! OFF 
     857       
     858         CALL ice_thd_enmelt        
     859 
     860         !     zhfx_err = correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation 
    591861         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)  
     862            zdq = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * h_i_1d(ji) * r1_nlay_i +  & 
     863               &                   SUM( e_s_1d(ji,1:nlay_s) ) * h_s_1d(ji) * r1_nlay_s ) 
     864                
     865            IF ( ( k_jules == np_zdf_jules_OFF ) ) THEN 
     866                
     867                  IF( t_su_1d(ji) < rt0 ) THEN  ! case T_su < 0degC 
     868                     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)  
     869                  ELSE                          ! case T_su = 0degC 
     870                     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) 
     871                  ENDIF 
     872             
     873            ELSE ! RCV CASE 
     874             
     875               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) 
     876             
     877            ENDIF 
     878             
     879            ! total heat sink to be sent to the ocean 
     880            hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err 
     881 
     882            ! hfx_dif = Heat flux diagnostic of sensible heat used to warm/cool ice in W.m-2    
     883            hfx_dif_1d(ji) = hfx_dif_1d(ji) - zdq * r1_rdtice * a_i_1d(ji) 
     884    
     885         END DO  
     886            
     887         ! 
     888         ! --- SIMIP diagnostics 
     889         ! 
     890 
     891         DO ji = 1, npti 
     892            !--- Conduction fluxes (positive downwards) 
     893            diag_fc_bo_1d(ji) = diag_fc_bo_1d(ji) + fc_bo_i(ji) * a_i_1d(ji) / at_i_1d(ji) 
     894            diag_fc_su_1d(ji) = diag_fc_su_1d(ji) + fc_su(ji)   * a_i_1d(ji) / at_i_1d(ji) 
     895    
     896            !--- Snow-ice interfacial temperature (diagnostic SIMIP) 
     897            zfac = rn_cnd_s * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji) 
     898            IF( zh_s(ji) >= 1.e-3 .AND. zfac > epsi10 ) THEN 
     899               t_si_1d(ji) = ( rn_cnd_s       * zh_i(ji) * t_s_1d(ji,1) +   & 
     900                  &            ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) / zfac 
     901            ELSE 
     902               t_si_1d(ji) = t_su_1d(ji) 
     903            ENDIF 
    593904         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 
     905      
     906      ENDIF 
     907      ! 
     908      !--------------------------------------------------------------------------------------- 
     909      ! 11) Jules coupling: reset inner snow and ice temperatures, update conduction fluxes 
     910      !--------------------------------------------------------------------------------------- 
     911      ! 
     912      IF ( k_jules == np_zdf_jules_SND ) THEN   ! --- Jules coupling in "SND" mode 
     913       
     914         ! Restore temperatures to their initial values 
     915         t_s_1d(1:npti,:)        = ztsold (1:npti,:) 
     916         t_i_1d(1:npti,:)        = ztiold (1:npti,:) 
     917         qcn_ice_1d(1:npti)      = fc_su(1:npti) 
     918          
     919      ENDIF 
     920       
     921   END SUBROUTINE ice_thd_zdf_BL99 
     922 
    632923 
    633924 
     
    663954 
    664955 
     956 
    665957   SUBROUTINE ice_thd_zdf_init 
    666958      !!----------------------------------------------------------------------- 
     
    675967      !! ** input   :   Namelist namthd_zdf 
    676968      !!------------------------------------------------------------------- 
    677       INTEGER  ::   ios   ! Local integer output status for namelist read 
     969      INTEGER  ::   ios, ioptio   ! Local integer output status for namelist read 
    678970      !! 
    679       NAMELIST/namthd_zdf/ ln_zdf_BL99, ln_cndi_U64, ln_cndi_P07, rn_cnd_s, rn_kappa_i, ln_dqns_i 
     971      NAMELIST/namthd_zdf/ ln_zdf_BL99, ln_cndi_U64, ln_cndi_P07, rn_cnd_s, rn_kappa_i 
    680972      !!------------------------------------------------------------------- 
    681973      ! 
     
    694986         WRITE(numout,*) '~~~~~~~~~~~~~~~~' 
    695987         WRITE(numout,*) '   Namelist namthd_zdf:' 
    696          WRITE(numout,*) '      Diffusion follows a Bitz and Lipscomb (1999)            ln_zdf_BL99  = ', ln_zdf_BL99 
     988         WRITE(numout,*) '      Bitz and Lipscomb (1999) formulation                    ln_zdf_BL99  = ', ln_zdf_BL99 
    697989         WRITE(numout,*) '      thermal conductivity in the ice (Untersteiner 1964)     ln_cndi_U64  = ', ln_cndi_U64 
    698990         WRITE(numout,*) '      thermal conductivity in the ice (Pringle et al 2007)    ln_cndi_P07  = ', ln_cndi_P07 
    699991         WRITE(numout,*) '      thermal conductivity in the snow                        rn_cnd_s     = ', rn_cnd_s 
    700992         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 
    702993     ENDIF 
     994 
    703995      ! 
    704996      IF ( ( ln_cndi_U64 .AND. ln_cndi_P07 ) .OR. ( .NOT.ln_cndi_U64 .AND. .NOT.ln_cndi_P07 ) ) THEN 
    705997         CALL ctl_stop( 'ice_thd_zdf_init: choose one and only one formulation for thermal conduction (ln_cndi_U64 or ln_cndi_P07)' ) 
    706998      ENDIF 
     999 
     1000      !                             !== set the choice of ice vertical thermodynamic formulation ==! 
     1001      ioptio = 0  
     1002      !      !--- BL99 thermo dynamics                               (linear liquidus + constant thermal properties) 
     1003      IF( ln_zdf_BL99 ) THEN   ;   ioptio = ioptio + 1   ;   nice_zdf = np_BL99   ;   ENDIF 
     1004      IF( ioptio /= 1 )    CALL ctl_stop( 'ice_thd_init: one and only one ice thermo option has to be defined ' ) 
    7071005      ! 
    7081006   END SUBROUTINE ice_thd_zdf_init 
     
    7111009   !!---------------------------------------------------------------------- 
    7121010   !!   Default option       Dummy Module             No ESIM sea-ice model 
    713    !!---------------------------------------------------------------------- 
     1011   !!--------------------------------------------------------------------- 
    7141012#endif 
    7151013 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_ice.F90

    r8563 r8726  
    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) )       
     
    168170   REAL            , 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 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90

    r8589 r8726  
    1717   !!                                          ==> 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 
     
    4243   USE ice     , ONLY :   u_ice, v_ice, jpl, a_i_b, at_i_b, tm_su 
    4344   USE icethd_dh      ! for CALL ice_thd_snwblow 
     45   USE icethd_zdf, ONLY:   rn_cnd_s  ! for blk_ice_qcn 
    4446#endif 
    4547   USE sbcblk_algo_ncar     ! => turb_ncar     : NCAR - CORE (Large & Yeager, 2009)  
     
    6466   PUBLIC   blk_ice_tau   ! routine called in icestp module 
    6567   PUBLIC   blk_ice_flx   ! routine called in icestp module 
    66 #endif 
     68   PUBLIC   blk_ice_qcn   ! routine called in icestp module 
     69#endif  
    6770 
    6871!!Lolo: should ultimately be moved in the module with all physical constants ? 
     
    697700 
    698701 
    699    SUBROUTINE blk_ice_flx( ptsu, palb ) 
     702   SUBROUTINE blk_ice_flx( ptsu, phs, phi, palb ) 
    700703      !!--------------------------------------------------------------------- 
    701704      !!                     ***  ROUTINE blk_ice_flx  *** 
     
    710713      !!--------------------------------------------------------------------- 
    711714      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   ptsu   ! sea ice surface temperature 
     715      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phs    ! snow thickness 
     716      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phi    ! ice thickness 
    712717      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   palb   ! ice albedo (all skies) 
    713718      !! 
     
    716721      REAL(wp) ::   zcoef_dqlw, zcoef_dqla   !   -      - 
    717722      REAL(wp) ::   zztmp, z1_lsub           !   -      - 
     723      REAL(wp) ::   zfrqsr_tr_i              ! sea ice shortwave fraction transmitted below through the ice 
     724      REAL(wp) ::   zfr1, zfr2, zfac         ! local variables 
    718725      REAL(wp), DIMENSION(:,:,:), POINTER ::   z_qlw         ! long wave heat flux over ice 
    719726      REAL(wp), DIMENSION(:,:,:), POINTER ::   z_qsb         ! sensible  heat flux over ice 
     
    722729      REAL(wp), DIMENSION(:,:)  , POINTER ::   zevap, zsnw   ! evaporation and snw distribution after wind blowing (LIM3) 
    723730      REAL(wp), DIMENSION(:,:)  , POINTER ::   zrhoa 
     731 
    724732      !!--------------------------------------------------------------------- 
    725733      ! 
     
    830838      CALL wrk_dealloc( jpi,jpj,   zevap, zsnw ) 
    831839 
    832       !-------------------------------------------------------------------- 
    833       ! FRACTIONs of net shortwave radiation which is not absorbed in the 
    834       ! thin surface layer and penetrates inside the ice cover 
    835       ! ( Maykut and Untersteiner, 1971 ; Ebert and Curry, 1993 ) 
    836       ! 
    837       fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice ) 
    838       fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) 
    839       ! 
    840       ! 
     840      ! --- absorbed and transmitted shortwave radiation (W/m2) --- ! 
     841      ! 
     842      ! former coding was 
     843      !     fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice ) 
     844      !     fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) 
     845 
     846      ! --- surface transmission parameter (i0, Grenfell Maykut 77) --- ! 
     847      zfr1 = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )            !   standard value 
     848      zfr2 = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )            !   zfr2 such that zfr1 + zfr2 to equal 1 
     849 
     850      qsr_ice_tr(:,:,:) = 0._wp 
     851 
     852      DO jl = 1, jpl 
     853         DO jj = 1, jpj 
     854            DO ji = 1, jpi 
     855 
     856               zfac = MAX( 0._wp , 1._wp - phi(ji,jj,jl) * 10._wp )     !   linear weighting factor: =0 for phi=0, =1 for phi = 0.1 
     857               zfrqsr_tr_i = zfr1 + zfac * zfr2                         !   below 10 cm, linearly increase zfrqsr_tr_i until 1 at zero thickness 
     858 
     859               IF ( phs(ji,jj,jl) <= 0.0_wp ) THEN  
     860                  zfrqsr_tr_i  = zfr1            + zfac * zfr2 
     861               ELSE  
     862                  zfrqsr_tr_i  = 0._wp                                  !   snow fully opaque 
     863               ENDIF 
     864 
     865               qsr_ice_tr(ji,jj,jl) = zfrqsr_tr_i * qsr_ice(ji,jj,jl)   !   transmitted solar radiation  
     866 
     867            END DO 
     868         END DO 
     869      END DO 
     870                                                                                       
     871 
    841872      IF(ln_ctl) THEN 
    842873         CALL prt_ctl(tab3d_1=qla_ice , clinfo1=' blk_ice: qla_ice  : ', tab3d_2=z_qsb   , clinfo2=' z_qsb    : ', kdim=jpl) 
     
    854885 
    855886   END SUBROUTINE blk_ice_flx 
     887    
     888    
     889 
     890   SUBROUTINE blk_ice_qcn( k_monocat, ptsu, ptb, phs, phi ) 
     891    
     892      !!--------------------------------------------------------------------- 
     893      !!                     ***  ROUTINE blk_ice_qcn  *** 
     894      !! 
     895      !! ** Purpose :   Compute surface temperature and snow/ice conduction flux 
     896      !!                to force sea ice / snow thermodynamics 
     897      !!                in the case JULES coupler is emulated 
     898      !!                 
     899      !! ** Method  :   compute surface energy balance assuming neglecting heat storage 
     900      !!                following the 0-layer Semtner (1976) approach 
     901      !! 
     902      !! ** Outputs : - ptsu    : sea-ice / snow surface temperature (K) 
     903      !!              - qcn_ice : surface inner conduction flux (W/m2) 
     904      !! 
     905      !!--------------------------------------------------------------------- 
     906      !! 
     907      INTEGER, INTENT(in)                        ::   k_monocat                 ! single-category option 
     908 
     909      REAL(wp), DIMENSION(:,:,:), INTENT(inout)  ::   ptsu                      ! sea ice / snow surface temperature 
     910 
     911      REAL(wp), DIMENSION(:,:)  , INTENT(in)     ::   ptb                       ! sea ice base temperature 
     912      REAL(wp), DIMENSION(:,:,:), INTENT(in)     ::   phs                       ! snow thickness 
     913      REAL(wp), DIMENSION(:,:,:), INTENT(in)     ::   phi                       ! sea ice thickness 
     914       
     915      INTEGER                                    ::   ji, jj, jl                ! dummy loop indices 
     916      INTEGER                                    ::   iter                      ! 
     917      REAL(wp)                                   ::   zfac, zfac2, zfac3        ! dummy factors 
     918      REAL(wp)                                   ::   zkeff_h, ztsu             ! 
     919      REAL(wp)                                   ::   zqc, zqnet                ! 
     920      REAL(wp)                                   ::   zhe, zqa0                 ! 
     921            
     922      INTEGER , PARAMETER                        ::   nit = 10                  ! number of iterations 
     923      REAL(wp), PARAMETER                        ::   zepsilon = 0.1_wp         ! characteristic thickness for enhanced conduction 
     924 
     925      REAL(wp), DIMENSION(:,:,:), POINTER        ::   zgfac                    ! enhanced conduction factor 
     926       
     927      !!--------------------------------------------------------------------- 
     928 
     929      IF( nn_timing == 1 )  CALL timing_start('blk_ice_qcn') 
     930      ! 
     931      CALL wrk_alloc( jpi,jpj,jpl,   zgfac ) 
     932 
     933      ! -------------------------------------! 
     934      !      I   Enhanced conduction factor  ! 
     935      ! -------------------------------------! 
     936      ! 
     937      ! Emulates the enhancement of conduction by unresolved thin ice (k_monocat = 1/3) 
     938      ! Fichefet and Morales Maqueda, JGR 1997 
     939      ! 
     940      zgfac(:,:,:) = 1._wp 
     941       
     942      SELECT CASE ( k_monocat )  
     943       
     944         CASE ( 1 , 3 ) 
     945       
     946            zfac    =   1._wp /  ( rn_cnd_s + rcdic ) 
     947            zfac2   =   EXP(1._wp) * 0.5_wp * zepsilon 
     948            zfac3   =   2._wp / zepsilon 
     949       
     950            DO jl = 1, jpl                 
     951               DO jj = 1 , jpj 
     952                  DO ji = 1, jpi 
     953                     !                                ! Effective thickness 
     954                     zhe               =   ( rn_cnd_s * phi(ji,jj,jl) + rcdic * phs(ji,jj,jl) ) * zfac 
     955                     !                                ! Enhanced conduction factor 
     956                     IF( zhe >=  zfac2 )  & 
     957                        zgfac(ji,jj,jl)   =   MIN( 2._wp, ( 0.5_wp + 0.5 * LOG( zhe * zfac3 ) ) ) 
     958                  END DO 
     959               END DO 
     960            END DO 
     961       
     962      END SELECT 
     963       
     964      ! -------------------------------------------------------------! 
     965      !      II   Surface temperature and conduction flux            ! 
     966      ! -------------------------------------------------------------! 
     967 
     968      zfac   =   rcdic * rn_cnd_s  
     969                                      ! ========================== ! 
     970      DO jl = 1, jpl                  !  Loop over ice categories  ! 
     971         !                            ! ========================== ! 
     972         DO jj = 1 , jpj 
     973            DO ji = 1, jpi 
     974               !                      ! Effective conductivity of the snow-ice system divided by thickness 
     975               zkeff_h   =   zfac * zgfac(ji,jj,jl) / ( rcdic * phs(ji,jj,jl) + rn_cnd_s * phi(ji,jj,jl) ) 
     976                                      ! Store initial surface temperature 
     977               ztsu      =   ptsu(ji,jj,jl) 
     978                                      ! Net initial atmospheric heat flux 
     979               zqa0      =   qsr_ice(ji,jj,jl) - qsr_ice_tr(ji,jj,jl) + qns_ice(ji,jj,jl) 
     980                
     981               DO iter = 1, nit       ! --- Iteration loop 
     982                
     983                   !                  ! Conduction heat flux through snow-ice system (>0 downwards) 
     984                   zqc   =   zkeff_h * ( ztsu - ptb(ji,jj) ) 
     985                   !                  ! Surface energy budget 
     986                   zqnet =   zqa0 + dqns_ice(ji,jj,jl) * ( ztsu - ptsu(ji,jj,jl) ) - zqc 
     987                   !                  ! Temperature update 
     988                   ztsu  =   ztsu - zqnet / ( dqns_ice(ji,jj,jl) - zkeff_h ) 
     989                
     990               END DO 
     991                
     992               ptsu(ji,jj,jl)    = MIN( rt0, ztsu ) 
     993                
     994               qcn_ice(ji,jj,jl) = zkeff_h * ( ptsu(ji,jj,jl) - ptb(ji,jj) ) 
     995                
     996            END DO 
     997         END DO 
     998          
     999      END DO  
     1000       
     1001      CALL wrk_dealloc( jpi,jpj,jpl,   zgfac ) 
     1002       
     1003      IF( nn_timing == 1 )  CALL timing_stop('blk_ice_qcn') 
     1004 
     1005   END SUBROUTINE blk_ice_qcn 
    8561006    
    8571007#endif 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90

    r8592 r8726  
    15241524    
    15251525 
    1526    SUBROUTINE sbc_cpl_ice_flx( picefr, palbi, psst, pist ) 
     1526   SUBROUTINE sbc_cpl_ice_flx( picefr, palbi, psst, pist, phs, phi ) 
    15271527      !!---------------------------------------------------------------------- 
    15281528      !!             ***  ROUTINE sbc_cpl_ice_flx  *** 
     
    15791579      REAL(wp), INTENT(in), DIMENSION(:,:  ), OPTIONAL ::   psst       ! sea surface temperature     [Celsius] 
    15801580      REAL(wp), INTENT(in), DIMENSION(:,:,:), OPTIONAL ::   pist       ! ice surface temperature     [Kelvin] 
    1581       ! 
    1582       INTEGER ::   jl         ! dummy loop index 
     1581      REAL(wp), INTENT(in), DIMENSION(:,:,:), OPTIONAL ::   phs        ! snow depth                  [m] 
     1582      REAL(wp), INTENT(in), DIMENSION(:,:,:), OPTIONAL ::   phi        ! ice thickness               [m] 
     1583      ! 
     1584      INTEGER ::   ji,jj,jl         ! dummy loop index 
    15831585      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zcptn, zcptrain, zcptsnw, ziceld, zmsk, zsnw 
    15841586      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap_oce, zevap_ice, zdevap_ice 
    15851587      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice 
    1586       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice 
     1588      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice, zfrqsr_tr_i 
    15871589      !!---------------------------------------------------------------------- 
    15881590      ! 
     
    15921594      CALL wrk_alloc( jpi,jpj,     zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap_oce, zevap_ice, zdevap_ice ) 
    15931595      CALL wrk_alloc( jpi,jpj,     zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice ) 
    1594       CALL wrk_alloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice ) 
     1596      CALL wrk_alloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice, zfrqsr_tr_i ) 
    15951597 
    15961598      IF( ln_mixcpl )   zmsk(:,:) = 1. - xcplmask(:,:,0) 
     
    19451947      END SELECT 
    19461948 
    1947       ! Surface transimission parameter io (Maykut Untersteiner , 1971 ; Ebert and Curry, 1993 ) 
    1948       ! Used for LIM3 
    1949       ! Coupled case: since cloud cover is not received from atmosphere  
    1950       !               ===> used prescribed cloud fraction representative for polar oceans in summer (0.81) 
    1951       fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice ) 
    1952       fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) 
     1949      ! --- Transmitted shortwave radiation (W/m2) --- ! 
     1950       
     1951      IF ( nice_jules == 0 ) THEN 
     1952                
     1953         zfrqsr_tr_i(:,:,:) = 0._wp   !   surface transmission parameter 
     1954      
     1955         ! former coding was     
     1956         ! Coupled case: since cloud cover is not received from atmosphere  
     1957         !               ===> used prescribed cloud fraction representative for polar oceans in summer (0.81) 
     1958         !     fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice ) 
     1959         !     fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) 
     1960          
     1961         ! to retrieve that coding, we needed to access h_i & h_s from here 
     1962         ! we could even retrieve cloud fraction from the coupler 
     1963                
     1964         DO jl = 1, jpl 
     1965            DO jj = 1 , jpj 
     1966               DO ji = 1, jpi 
     1967             
     1968                  !--- surface transmission parameter (Grenfell Maykut 77) --- ! 
     1969                  zfrqsr_tr_i(ji,jj,jl) = 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice  
     1970                
     1971                  ! --- influence of snow and thin ice --- ! 
     1972                  IF ( phs(ji,jj,jl) >= 0.0_wp ) zfrqsr_tr_i(ji,jj,jl) = 0._wp   !   snow fully opaque 
     1973                  IF ( phi(ji,jj,jl) <= 0.1_wp ) zfrqsr_tr_i(ji,jj,jl) = 1._wp   !   thin ice transmits all solar radiation 
     1974               END DO 
     1975            END DO 
     1976         END DO 
     1977          
     1978         qsr_ice_tr(:,:,:) =   zfrqsr_tr_i(:,:,:) * qsr_ice(:,:,:)               !   transmitted solar radiation  
     1979                
     1980      ENDIF 
     1981       
     1982      IF ( nice_jules == 2 ) THEN 
     1983       
     1984         ! here we must receive the qsr_ice_tr array from the coupler 
     1985         ! for now just assume zero 
     1986          
     1987         qsr_ice_tr(:,:,:) = 0.0_wp 
     1988       
     1989      ENDIF 
     1990 
     1991 
    19531992 
    19541993      CALL wrk_dealloc( jpi,jpj,     zcptn, zcptrain, zcptsnw, ziceld, zmsk, zsnw ) 
    19551994      CALL wrk_dealloc( jpi,jpj,     zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap_oce, zevap_ice, zdevap_ice ) 
    19561995      CALL wrk_dealloc( jpi,jpj,     zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice ) 
    1957       CALL wrk_dealloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice ) 
     1996      CALL wrk_dealloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice, zfrqsr_tr_i ) 
    19581997      ! 
    19591998      IF( nn_timing == 1 )   CALL timing_stop('sbc_cpl_ice_flx') 
Note: See TracChangeset for help on using the changeset viewer.