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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 8752 for branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3/icethd_zdf.F90 – NEMO

Ignore:
Timestamp:
2017-11-20T13:54:32+01:00 (6 years ago)
Author:
dancopsey
Message:

Merged in main ICEMODEL branch (branches/2017/dev_r8183_ICEMODEL) from revision 8587 to 8726.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3/icethd_zdf.F90

    r8738 r8752  
    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 
Note: See TracChangeset for help on using the changeset viewer.