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 9274 – NEMO

Changeset 9274


Ignore:
Timestamp:
2018-01-22T17:27:39+01:00 (6 years ago)
Author:
clem
Message:

cosmetics and progressing in preparing multi snow layers in thermo

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/LIM_SRC_3/icethd_dh.F90

    r9271 r9274  
    3131 
    3232   PUBLIC   ice_thd_dh        ! called by ice_thd 
    33    PUBLIC   ice_thd_snwblow   ! called in sbcblk/sbcclio/sbccpl and here 
     33   PUBLIC   ice_thd_snwblow   ! called in sbcblk/sbccpl and here 
    3434 
    3535   INTERFACE ice_thd_snwblow 
     
    4848      !!                ***  ROUTINE ice_thd_dh  *** 
    4949      !! 
    50       !! ** Purpose :   compute ice and snow thickness changes due to growing/melting 
     50      !! ** Purpose :   compute ice and snow thickness changes due to growth/melting 
    5151      !! 
    5252      !! ** Method  :   Ice/Snow surface melting arises from imbalance in surface fluxes 
    53       !!              Bottom accretion/ablation arises from flux budget 
    54       !!              Snow thickness can increase by precipitation and decrease by sublimation 
    55       !!              If snow load excesses Archmiede limit, snow-ice is formed by 
    56       !!              the flooding of sea-water in the snow 
     53      !!                Bottom accretion/ablation arises from flux budget 
     54      !!                Snow thickness can increase by precipitation and decrease by sublimation 
     55      !!                If snow load excesses Archmiede limit, snow-ice is formed by 
     56      !!                the flooding of sea-water in the snow 
    5757      !! 
    58       !!                 1) Compute available flux of heat for surface ablation 
    59       !!                 2) Compute snow and sea ice enthalpies 
    60       !!                 3) Surface ablation and sublimation 
    61       !!                 4) Bottom accretion/ablation 
    62       !!                 5) Case of Total ablation 
    63       !!                 6) Snow ice formation 
     58      !!                - Compute available flux of heat for surface ablation 
     59      !!                - Compute snow and sea ice enthalpies 
     60      !!                - Surface ablation and sublimation 
     61      !!                - Bottom accretion/ablation 
     62      !!                - Snow ice formation 
    6463      !! 
    6564      !! References : Bitz and Lipscomb, 1999, J. Geophys. Res. 
     
    126125      END DO 
    127126      ! 
    128       !------------------------------------------------------------------------------! 
    129       !  1) Calculate available heat for surface and bottom ablation                ! 
    130       !------------------------------------------------------------------------------! 
    131       ! 
    132       SELECT CASE ( nice_jules ) 
    133       ! 
    134       CASE ( np_jules_ACTIVE ) 
     127      ! ---------------------------------------------- ! 
     128      ! Available heat for surface and bottom ablation ! 
     129      ! ---------------------------------------------- ! 
     130      ! 
     131      SELECT CASE( nice_jules ) 
     132      ! 
     133      CASE( np_jules_ACTIVE ) 
    135134         ! 
    136135         DO ji = 1, npti 
    137             zq_su(ji)        =  MAX( 0._wp, qml_ice_1d(ji) * rdt_ice ) 
     136            zq_su(ji)      = MAX( 0._wp, qml_ice_1d(ji) * rdt_ice ) 
    138137         END DO 
    139138         ! 
    140       CASE ( np_jules_EMULE ) 
     139      CASE( np_jules_EMULE ) 
    141140         ! 
    142141         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 ) 
     142            zdum           = qns_ice_1d(ji) + qsr_ice_1d(ji) - qsr_ice_tr_1d(ji) - fc_su(ji) 
     143            qml_ice_1d(ji) = zdum * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) ) 
     144            zq_su(ji)      = MAX( 0._wp, qml_ice_1d(ji) * rdt_ice ) 
    146145         END DO 
    147146         ! 
    148       CASE ( np_jules_OFF )  
     147      CASE( np_jules_OFF )  
    149148         ! 
    150149         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 ) 
     150            zdum           = qns_ice_1d(ji) + qsr_ice_1d(ji) - qsr_ice_tr_1d(ji) - fc_su(ji)  
     151            zdum           = zdum * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) ) 
     152            zq_su(ji)      = MAX( 0._wp, zdum * rdt_ice ) 
    154153         END DO 
    155154         ! 
     
    157156      ! 
    158157      DO ji = 1, npti 
    159          zf_tt(ji)  = fc_bo_i(ji) + fhtur_1d(ji) + fhld_1d(ji)  
    160          zq_bo (ji) = MAX( 0._wp, zf_tt(ji) * rdt_ice ) 
    161       END DO 
    162  
    163       !--------------------------------------------! 
    164       !  2) Computing layer thicknesses            ! 
    165       !--------------------------------------------! 
     158         zf_tt(ji)         = fc_bo_i(ji) + fhtur_1d(ji) + fhld_1d(ji)  
     159         zq_bo(ji)         = MAX( 0._wp, zf_tt(ji) * rdt_ice ) 
     160      END DO 
     161 
     162      ! Ice and snow layer thicknesses 
     163      !------------------------------- 
    166164      DO jk = 1, nlay_i 
    167165         DO ji = 1, npti 
     
    176174      END DO 
    177175 
    178       !------------------------------------------------------------------------------! 
    179       ! If snow temperature is above freezing point, then snow melts  
    180       ! (should not happen but sometimes it does) 
    181       !------------------------------------------------------------------------------! 
     176      !                       ! ============ ! 
     177      !                       !     Snow     ! 
     178      !                       ! ============ ! 
     179      ! 
     180      ! Internal melting 
     181      ! ---------------- 
     182      ! IF snow temperature is above freezing point, THEN snow melts (should not happen but sometimes it does) 
    182183      DO jk = 1, nlay_s 
    183184         DO ji = 1, npti 
    184             IF( t_s_1d(ji,jk) > rt0 ) THEN !!! Internal melting 
    185                ! Contribution to heat flux to the ocean [W.m-2], < 0   
    186                hfx_res_1d(ji) = hfx_res_1d(ji) + e_s_1d(ji,jk) * zh_s(ji,jk) * a_i_1d(ji) * r1_rdtice 
    187                ! Contribution to mass flux 
    188                wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + rhosn * zh_s(ji,jk) * a_i_1d(ji) * r1_rdtice 
    189                dh_s_mlt(ji)       = dh_s_mlt(ji) - zh_s(ji,jk) 
     185            IF( t_s_1d(ji,jk) > rt0 ) THEN 
     186               hfx_res_1d    (ji) = hfx_res_1d    (ji) + e_s_1d(ji,jk) * zh_s(ji,jk) * a_i_1d(ji) * r1_rdtice   ! heat flux to the ocean [W.m-2], < 0 
     187               wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + rhosn         * zh_s(ji,jk) * a_i_1d(ji) * r1_rdtice   ! mass flux 
    190188               ! updates 
    191                h_s_1d(ji)    = h_s_1d(ji) - zh_s(ji,jk) 
    192                zh_s  (ji,jk) = 0._wp 
    193                e_s_1d(ji,jk) = 0._wp 
    194                t_s_1d(ji,jk) = rt0 
     189               dh_s_mlt(ji)    = dh_s_mlt(ji) - zh_s(ji,jk) 
     190               h_s_1d  (ji)    = h_s_1d(ji) - zh_s(ji,jk) 
     191               zh_s    (ji,jk) = 0._wp 
     192               e_s_1d  (ji,jk) = 0._wp 
     193               t_s_1d  (ji,jk) = rt0 
    195194            END IF 
    196195         END DO 
    197       END DO 
    198           
    199  
    200       !------------------------------------------------------------------------------| 
    201       !  3) Surface ablation and sublimation                                         | 
    202       !------------------------------------------------------------------------------| 
    203       ! 
    204       !------------------------- 
    205       ! 3.1 Snow precips / melt 
    206       !------------------------- 
    207       ! Snow accumulation in one thermodynamic time step 
    208       ! snowfall is partitionned between leads and ice 
    209       ! if snow fall was uniform, a fraction (1-at_i) would fall into leads 
    210       ! but because of the winds, more snow falls on leads than on sea ice 
    211       ! and a greater fraction (1-at_i)^beta of the total mass of snow  
    212       ! (beta < 1) falls in leads. 
    213       ! In reality, beta depends on wind speed,  
    214       ! and should decrease with increasing wind speed but here, it is  
    215       ! considered as a constant. an average value is 0.66 
    216       ! Martin Vancoppenolle, December 2006 
    217  
    218       CALL ice_thd_snwblow( 1. - at_i_1d(1:npti), zsnw(1:npti) ) ! snow distribution over ice after wind blowing 
     196      END DO          
     197 
     198      ! Snow precipitation 
     199      !------------------- 
     200      CALL ice_thd_snwblow( 1. - at_i_1d(1:npti), zsnw(1:npti) )   ! snow distribution over ice after wind blowing 
    219201 
    220202      zdeltah(1:npti,:) = 0._wp 
    221203      DO ji = 1, npti 
    222          !----------- 
    223          ! Snow fall 
    224          !----------- 
    225          ! thickness change 
    226          zdh_s_pre(ji) = zsnw(ji) * sprecip_1d(ji) * rdt_ice * r1_rhosn / at_i_1d(ji) 
    227          ! enthalpy of the precip (>0, J.m-3) 
    228          zqprec   (ji) = - qprec_ice_1d(ji)    
    229          IF( sprecip_1d(ji) == 0._wp ) zqprec(ji) = 0._wp 
    230          ! heat flux from snow precip (>0, W.m-2) 
    231          hfx_spr_1d(ji) = hfx_spr_1d(ji) + zdh_s_pre(ji) * a_i_1d(ji) * zqprec(ji) * r1_rdtice 
    232          ! mass flux, <0 
    233          wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_pre(ji) * r1_rdtice 
    234  
    235          !--------------------- 
    236          ! Melt of falling snow 
    237          !--------------------- 
    238          ! thickness change 
    239          rswitch        = MAX( 0._wp , SIGN( 1._wp , zqprec(ji) - epsi20 ) ) 
    240          zdeltah (ji,1) = - rswitch * zq_su(ji) / MAX( zqprec(ji) , epsi20 ) 
    241          zdeltah (ji,1) = MAX( - zdh_s_pre(ji), zdeltah(ji,1) ) ! bound melting  
    242          ! heat used to melt snow (W.m-2, >0) 
    243          hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * zqprec(ji) * r1_rdtice 
    244          ! snow melting only = water into the ocean (then without snow precip), >0 
    245          wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice 
    246          dh_s_mlt(ji)       = dh_s_mlt(ji) + zdeltah(ji,1) 
    247          ! updates available heat + precipitations after melting 
    248          zq_su     (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,1) * zqprec(ji) )       
    249          zdh_s_pre (ji) = zdh_s_pre(ji) + zdeltah(ji,1) 
    250  
    251          ! update thickness 
    252          h_s_1d(ji) = MAX( 0._wp , h_s_1d(ji) + zdh_s_pre(ji) ) 
     204         IF( sprecip_1d(ji) > 0._wp ) THEN 
     205            ! 
     206            ! --- precipitation --- 
     207            zdh_s_pre (ji) = zsnw(ji) * sprecip_1d(ji) * rdt_ice * r1_rhosn / at_i_1d(ji)   ! thickness change 
     208            zqprec    (ji) = - qprec_ice_1d(ji)                                             ! enthalpy of the precip (>0, J.m-3) 
     209            ! 
     210            hfx_spr_1d(ji) = hfx_spr_1d(ji) + zdh_s_pre(ji) * a_i_1d(ji) * zqprec(ji)    * r1_rdtice   ! heat flux from snow precip (>0, W.m-2) 
     211            wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhosn         * a_i_1d(ji) * zdh_s_pre(ji) * r1_rdtice   ! mass flux, <0 
     212             
     213            ! --- melt of falling snow --- 
     214            rswitch              = MAX( 0._wp , SIGN( 1._wp , zqprec(ji) - epsi20 ) ) 
     215            zdeltah       (ji,1) = - rswitch * zq_su(ji) / MAX( zqprec(ji) , epsi20 )   ! thickness change 
     216            zdeltah       (ji,1) = MAX( - zdh_s_pre(ji), zdeltah(ji,1) )                ! bound melting  
     217            hfx_snw_1d    (ji)   = hfx_snw_1d    (ji) - zdeltah(ji,1) * a_i_1d(ji) * zqprec(ji)    * r1_rdtice   ! heat used to melt snow (W.m-2, >0) 
     218            wfx_snw_sum_1d(ji)   = wfx_snw_sum_1d(ji) - rhosn         * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice   ! snow melting only = water into the ocean (then without snow precip), >0 
     219             
     220            ! updates available heat + precipitations after melting 
     221            dh_s_mlt (ji) = dh_s_mlt(ji) + zdeltah(ji,1) 
     222            zq_su    (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,1) * zqprec(ji) )       
     223            zdh_s_pre(ji) = zdh_s_pre(ji) + zdeltah(ji,1) 
     224             
     225            ! update thickness 
     226            h_s_1d(ji) = MAX( 0._wp , h_s_1d(ji) + zdh_s_pre(ji) ) 
     227            ! 
     228         ELSE 
     229            ! 
     230            zdh_s_pre(ji) = 0._wp 
     231            zqprec   (ji) = 0._wp 
     232            ! 
     233         ENDIF 
    253234      END DO 
    254235 
     
    260241      END DO 
    261242 
     243      ! Snow melting 
     244      ! ------------ 
    262245      ! If heat still available (zq_su > 0), then melt more snow 
    263246      zdeltah(1:npti,:) = 0._wp 
     
    265248      DO jk = 1, nlay_s 
    266249         DO ji = 1, npti 
    267             ! thickness change 
    268             rswitch          = 1._wp - MAX( 0._wp, SIGN( 1._wp, - zh_s(ji,jk) ) )  
    269             rswitch          = rswitch * ( MAX( 0._wp, SIGN( 1._wp, e_s_1d(ji,jk) - epsi20 ) ) )  
    270             zdeltah  (ji,jk) = - rswitch * zq_su(ji) / MAX( e_s_1d(ji,jk), epsi20 ) 
    271             zdeltah  (ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji,jk) ) ! bound melting 
    272             zdh_s_mel(ji)    = zdh_s_mel(ji) + zdeltah(ji,jk)     
    273             ! heat used to melt snow(W.m-2, >0) 
    274             hfx_snw_1d(ji)   = hfx_snw_1d(ji) - zdeltah(ji,jk) * a_i_1d(ji) * e_s_1d(ji,jk) * r1_rdtice  
    275             ! snow melting only = water into the ocean (then without snow precip) 
    276             wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 
    277             dh_s_mlt(ji)       = dh_s_mlt(ji) + zdeltah(ji,jk) 
    278             ! updates available heat + thickness 
    279             zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * e_s_1d(ji,jk) ) 
    280             h_s_1d(ji) = MAX( 0._wp , h_s_1d(ji) + zdeltah(ji,jk) ) 
    281             zh_s(ji,jk)= MAX( 0._wp , zh_s(ji,jk) + zdeltah(ji,jk) ) 
    282          END DO 
    283       END DO 
    284  
    285       !------------------------------ 
    286       ! 3.2 Sublimation (part1: snow)  
    287       !------------------------------ 
     250            IF( zh_s(ji,jk) > 0._wp .AND. zq_su(ji) > 0._wp ) THEN 
     251               ! 
     252               rswitch          = MAX( 0._wp, SIGN( 1._wp, e_s_1d(ji,jk) - epsi20 ) ) 
     253               zdeltah  (ji,jk) = - rswitch * zq_su(ji) / MAX( e_s_1d(ji,jk), epsi20 )   ! thickness change 
     254               zdeltah  (ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji,jk) )                  ! bound melting 
     255               zdh_s_mel(ji)    = zdh_s_mel(ji) + zdeltah(ji,jk) 
     256                
     257               hfx_snw_1d(ji)     = hfx_snw_1d(ji)     - zdeltah(ji,jk) * a_i_1d(ji) * e_s_1d (ji,jk) * r1_rdtice   ! heat used to melt snow(W.m-2, >0) 
     258               wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn          * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice   ! snow melting only = water into the ocean (then without snow precip) 
     259                
     260               ! updates available heat + thickness 
     261               dh_s_mlt(ji)    = dh_s_mlt(ji) + zdeltah(ji,jk) 
     262               zq_su   (ji)    = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * e_s_1d(ji,jk) ) 
     263               h_s_1d  (ji)    = MAX( 0._wp , h_s_1d(ji) + zdeltah(ji,jk) ) 
     264               zh_s    (ji,jk) = MAX( 0._wp , zh_s(ji,jk) + zdeltah(ji,jk) ) 
     265               ! 
     266            ENDIF 
     267         END DO 
     268      END DO 
     269 
     270      ! Snow sublimation  
     271      !----------------- 
    288272      ! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates 
    289273      !    comment: not counted in mass/heat exchange in iceupdate.F90 since this is an exchange with atm. (not ocean) 
    290274      zdeltah(1:npti,:) = 0._wp 
    291275      DO ji = 1, npti 
    292          zdh_s_sub(ji)  = MAX( - h_s_1d(ji) , - evap_ice_1d(ji) * r1_rhosn * rdt_ice ) 
    293          ! remaining evap in kg.m-2 (used for ice melting later on) 
    294          zevap_rema(ji)  = evap_ice_1d(ji) * rdt_ice + zdh_s_sub(ji) * rhosn 
    295          ! Heat flux by sublimation [W.m-2], < 0 (sublimate first snow that had fallen, then pre-existing snow) 
    296          zdeltah(ji,1)  = MAX( zdh_s_sub(ji), - zdh_s_pre(ji) ) 
    297          hfx_sub_1d(ji) = hfx_sub_1d(ji) + ( zdeltah(ji,1) * zqprec(ji) + ( zdh_s_sub(ji) - zdeltah(ji,1) ) * e_s_1d(ji,1)  & 
    298             &                              ) * a_i_1d(ji) * r1_rdtice 
    299          ! Mass flux by sublimation 
    300          wfx_snw_sub_1d(ji) = wfx_snw_sub_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_sub(ji) * r1_rdtice 
    301  
    302          ! new snow thickness 
    303          h_s_1d(ji)    =  MAX( 0._wp , h_s_1d(ji) + zdh_s_sub(ji) ) 
    304          ! update precipitations after sublimation and correct sublimation 
    305          zdh_s_pre(ji) = zdh_s_pre(ji) + zdeltah(ji,1) 
    306          zdh_s_sub(ji) = zdh_s_sub(ji) - zdeltah(ji,1) 
     276         IF( evap_ice_1d(ji) > 0._wp ) THEN 
     277            ! 
     278            zdh_s_sub (ji)   = MAX( - h_s_1d(ji) , - evap_ice_1d(ji) * r1_rhosn * rdt_ice ) 
     279            zevap_rema(ji)   = evap_ice_1d(ji) * rdt_ice + zdh_s_sub(ji) * rhosn   ! remaining evap in kg.m-2 (used for ice melting later on) 
     280            zdeltah   (ji,1) = MAX( zdh_s_sub(ji), - zdh_s_pre(ji) ) 
     281             
     282            hfx_sub_1d    (ji) = hfx_sub_1d(ji) + &   ! Heat flux by sublimation [W.m-2], < 0 (sublimate snow that had fallen, then pre-existing snow) 
     283               &                 ( zdeltah(ji,1) * zqprec(ji) + ( zdh_s_sub(ji) - zdeltah(ji,1) ) * e_s_1d(ji,1) )  & 
     284               &                 * a_i_1d(ji) * r1_rdtice 
     285            wfx_snw_sub_1d(ji) = wfx_snw_sub_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_sub(ji) * r1_rdtice   ! Mass flux by sublimation 
     286             
     287            ! new snow thickness 
     288            h_s_1d(ji)    =  MAX( 0._wp , h_s_1d(ji) + zdh_s_sub(ji) ) 
     289            ! update precipitations after sublimation and correct sublimation 
     290            zdh_s_pre(ji) = zdh_s_pre(ji) + zdeltah(ji,1) 
     291            zdh_s_sub(ji) = zdh_s_sub(ji) - zdeltah(ji,1) 
     292            ! 
     293         ELSE 
     294            ! 
     295            zdh_s_sub (ji) = 0._wp 
     296            zevap_rema(ji) = 0._wp 
     297            ! 
     298         ENDIF 
    307299      END DO 
    308300       
     
    312304      END DO 
    313305 
    314       !------------------------------------------- 
    315       ! 3.3 Update temperature, energy 
    316       !------------------------------------------- 
     306      ! Update temperature, energy 
     307      !--------------------------- 
    317308      ! new temp and enthalpy of the snow (remaining snow precip + remaining pre-existing snow) 
    318309      DO jk = 1, nlay_s 
    319310         DO ji = 1,npti 
    320311            rswitch       = MAX( 0._wp , SIGN( 1._wp, h_s_1d(ji) - epsi20 ) ) 
    321             e_s_1d(ji,jk) = rswitch / MAX( h_s_1d(ji), epsi20 ) *           & 
    322               &            ( ( zdh_s_pre(ji)               ) * zqprec(ji) +  & 
    323               &              ( h_s_1d(ji) - zdh_s_pre(ji) ) * rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) ) 
    324          END DO 
    325       END DO 
    326  
    327       !-------------------------- 
    328       ! 3.4 Surface ice ablation  
    329       !-------------------------- 
     312            e_s_1d(ji,jk) = rswitch / MAX( h_s_1d(ji), epsi20 ) *            & 
     313              &             ( ( zdh_s_pre(ji)              ) * zqprec(ji) +  & 
     314              &               ( h_s_1d(ji) - zdh_s_pre(ji) ) * rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) ) 
     315         END DO 
     316      END DO 
     317       
     318      !                       ! ============ ! 
     319      !                       !     Ice      ! 
     320      !                       ! ============ ! 
     321 
     322      ! Surface ice melting  
     323      !-------------------- 
    330324      zdeltah(1:npti,:) = 0._wp ! important 
    331325      DO jk = 1, nlay_i 
    332326         DO ji = 1, npti 
    333             ztmelts           = - tmut * sz_i_1d(ji,jk)          ! Melting point of layer k [C] 
    334              
    335             IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN !!! Internal melting 
     327            ztmelts = - tmut * sz_i_1d(ji,jk)   ! Melting point of layer k [C] 
     328             
     329            IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN   !-- Internal melting 
    336330 
    337331               zEi            = - e_i_1d(ji,jk) * r1_rhoic            ! Specific enthalpy of layer k [J/kg, <0]        
    338                zdE            = 0._wp                                 ! Specific enthalpy difference   (J/kg, <0) 
     332               zdE            = 0._wp                                 ! Specific enthalpy difference (J/kg, <0) 
    339333                                                                      ! set up at 0 since no energy is needed to melt water...(it is already melted) 
    340334               zdeltah(ji,jk) = MIN( 0._wp , - zh_i(ji,jk) )          ! internal melting occurs when the internal temperature is above freezing      
     
    346340               zfmdt          = - rhoic * zdeltah(ji,jk)              ! Recompute mass flux [kg/m2, >0] 
    347341 
    348                ! Contribution to heat flux to the ocean [W.m-2], <0 (ice enthalpy zEi is "sent" to the ocean)  
    349                hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice 
    350                 
    351                ! Contribution to salt flux (using s_i_1d and not sz_i_1d(jk) is ok) 
    352                sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice 
    353                 
    354                ! Contribution to mass flux 
    355                wfx_res_1d(ji) =  wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 
    356  
    357             ELSE                                !!! Surface melting 
     342               hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice                           ! Heat flux to the ocean [W.m-2], <0 
     343               !                                                                                                  ice enthalpy zEi is "sent" to the ocean 
     344               sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice   ! Salt flux 
     345               !                                                                                                  using s_i_1d and not sz_i_1d(jk) is ok 
     346               wfx_res_1d(ji) = wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice                ! Mass flux 
     347 
     348            ELSE                                        !-- Surface melting 
    358349                
    359350               zEi            = - e_i_1d(ji,jk) * r1_rhoic            ! Specific enthalpy of layer k [J/kg, <0] 
     
    375366               zQm            = zfmdt * zEw                           ! Energy of the melt water sent to the ocean [J/m2, <0] 
    376367                
    377                ! Contribution to salt flux >0 (using s_i_1d and not sz_i_1d(jk) is ok) 
    378                sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice 
    379                 
    380                ! Contribution to heat flux [W.m-2], < 0 
    381                hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 
    382                 
    383                ! Total heat flux used in this process [W.m-2], > 0   
    384                hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice 
    385                 
    386                ! Contribution to mass flux 
    387                wfx_sum_1d(ji) =  wfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 
     368               sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice   ! Salt flux >0 
     369               !                                                                                                  using s_i_1d and not sz_i_1d(jk) is ok) 
     370               hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice                           ! Heat flux [W.m-2], < 0 
     371               hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice                           ! Heat flux used in this process [W.m-2], > 0   
     372               !  
     373               wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice                ! Mass flux 
    388374                
    389375            END IF 
    390             ! ---------------------- 
    391             ! Sublimation part2: ice 
    392             ! ---------------------- 
    393             zdum      = MAX( - ( zh_i(ji,jk) + zdeltah(ji,jk) ) , - zevap_rema(ji) * r1_rhoic ) 
    394             zdeltah(ji,jk) = zdeltah(ji,jk) + zdum 
    395             dh_i_sub(ji)  = dh_i_sub(ji) + zdum 
    396             ! Salt flux > 0 (clem: flux is sent to the ocean for simplicity but salt should remain in the ice except if all ice is melted. 
    397             !                          It must be corrected at some point) 
    398             sfx_sub_1d(ji) = sfx_sub_1d(ji) - rhoic * a_i_1d(ji) * zdum * s_i_1d(ji) * r1_rdtice 
    399             ! Heat flux [W.m-2], < 0 
    400             hfx_sub_1d(ji) = hfx_sub_1d(ji) + zdum * e_i_1d(ji,jk) * a_i_1d(ji) * r1_rdtice 
    401             ! Mass flux > 0 
    402             wfx_ice_sub_1d(ji) =  wfx_ice_sub_1d(ji) - rhoic * a_i_1d(ji) * zdum * r1_rdtice 
     376             
     377            ! Ice sublimation 
     378            ! --------------- 
     379            zdum            = MAX( - ( zh_i(ji,jk) + zdeltah(ji,jk) ) , - zevap_rema(ji) * r1_rhoic ) 
     380            zdeltah (ji,jk) = zdeltah (ji,jk) + zdum 
     381            dh_i_sub(ji)    = dh_i_sub(ji)    + zdum 
     382             
     383            sfx_sub_1d(ji)     = sfx_sub_1d(ji) - rhoic * a_i_1d(ji) * zdum * s_i_1d(ji) * r1_rdtice ! Salt flux >0 
     384            !                                                                                          clem: flux is sent to the ocean for simplicity 
     385            !                                                                                                but salt should remain in the ice except 
     386            !                                                                                                if all ice is melted. => must be corrected 
     387            hfx_sub_1d(ji)     = hfx_sub_1d(ji) + zdum * e_i_1d(ji,jk) * a_i_1d(ji) * r1_rdtice      ! Heat flux [W.m-2], < 0 
     388 
     389            wfx_ice_sub_1d(ji) = wfx_ice_sub_1d(ji) - rhoic * a_i_1d(ji) * zdum * r1_rdtice          ! Mass flux > 0 
    403390 
    404391            ! update remaining mass flux 
     
    417404         END DO 
    418405      END DO 
     406       
    419407      ! update ice thickness 
    420408      DO ji = 1, npti 
     
    428416 
    429417 
    430       !------------------------------------------------------------------------------! 
    431       ! 4) Basal growth / melt                                                       ! 
    432       !------------------------------------------------------------------------------! 
    433       ! 
    434       !------------------ 
    435       ! 4.1 Basal growth  
     418      ! Ice Basal growth  
    436419      !------------------ 
    437420      ! Basal growth is driven by heat imbalance at the ice-ocean interface, 
     
    449432      IF( nn_icesal == 2 )   num_iter_max = 5  ! salinity varying in time 
    450433       
    451       ! Iterative procedure 
    452434      DO ji = 1, npti 
    453435         IF(  zf_tt(ji) < 0._wp  ) THEN 
    454             DO iter = 1, num_iter_max 
     436            DO iter = 1, num_iter_max   ! iterations 
    455437 
    456438               ! New bottom ice salinity (Cox & Weeks, JGR88 ) 
     
    458440               !--- zswi12 if 2.0e-8 < dh/dt < 3.6e-7  
    459441               !--- zswi2  if dh/dt > 3.6e-7 
    460                zgrr               = MIN( 1.0e-3, MAX ( dh_i_bott(ji) * r1_rdtice , epsi10 ) ) 
    461                zswi2              = MAX( 0._wp , SIGN( 1._wp , zgrr - 3.6e-7 ) ) 
    462                zswi12             = MAX( 0._wp , SIGN( 1._wp , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 ) 
    463                zswi1              = 1. - zswi2 * zswi12 
    464                zfracs             = MIN ( zswi1  * 0.12 + zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) )   & 
    465                   &               + zswi2  * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) )  , 0.5 ) 
    466  
    467                s_i_new(ji)        = zswitch_sal * zfracs * sss_1d(ji)  &  ! New ice salinity 
    468                                   + ( 1. - zswitch_sal ) * s_i_1d(ji)  
    469                ! New ice growth 
    470                ztmelts            = - tmut * s_i_new(ji)                 ! New ice melting point (C) 
    471  
    472                zt_i_new           = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i) 
    473                 
    474                zEi                = cpic * ( zt_i_new - (ztmelts+rt0) ) &     ! Specific enthalpy of forming ice (J/kg, <0)       
    475                   &               - lfus * ( 1.0 - ztmelts / ( zt_i_new - rt0 ) )   & 
    476                   &               + rcp  * ztmelts           
    477  
    478                zEw                = rcp  * ( t_bo_1d(ji) - rt0 )         ! Specific enthalpy of seawater (J/kg, < 0) 
    479  
    480                zdE                = zEi - zEw                           ! Specific enthalpy difference (J/kg, <0) 
    481  
    482                dh_i_bott(ji)      = rdt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoic ) ) 
     442               zgrr     = MIN( 1.0e-3, MAX ( dh_i_bott(ji) * r1_rdtice , epsi10 ) ) 
     443               zswi2    = MAX( 0._wp , SIGN( 1._wp , zgrr - 3.6e-7 ) ) 
     444               zswi12   = MAX( 0._wp , SIGN( 1._wp , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 ) 
     445               zswi1    = 1. - zswi2 * zswi12 
     446               zfracs   = MIN( zswi1  * 0.12 + zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) )   & 
     447                  &          + zswi2  * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) )  , 0.5 ) 
     448 
     449               s_i_new(ji)   = zswitch_sal * zfracs * sss_1d(ji) + ( 1. - zswitch_sal ) * s_i_1d(ji)  ! New ice salinity 
     450 
     451               ztmelts       = - tmut * s_i_new(ji)                                                   ! New ice melting point (C) 
     452 
     453               zt_i_new      = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i) 
     454                
     455               zEi           = cpic * ( zt_i_new - (ztmelts+rt0) ) &                                  ! Specific enthalpy of forming ice (J/kg, <0) 
     456                  &            - lfus * ( 1.0 - ztmelts / ( zt_i_new - rt0 ) ) + rcp  * ztmelts 
     457 
     458               zEw           = rcp  * ( t_bo_1d(ji) - rt0 )                                           ! Specific enthalpy of seawater (J/kg, < 0) 
     459 
     460               zdE           = zEi - zEw                                                              ! Specific enthalpy difference (J/kg, <0) 
     461 
     462               dh_i_bott(ji) = rdt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoic ) ) 
    483463                
    484464            END DO 
    485465            ! Contribution to Energy and Salt Fluxes                                     
    486             zfmdt          = - rhoic * dh_i_bott(ji)             ! Mass flux x time step (kg/m2, < 0) 
    487  
    488             ztmelts        = - tmut * s_i_new(ji)                ! New ice melting point (C) 
     466            zfmdt          = - rhoic * dh_i_bott(ji)                                                  ! Mass flux x time step (kg/m2, < 0) 
     467 
     468            ztmelts        = - tmut * s_i_new(ji)                                                     ! New ice melting point (C) 
    489469             
    490470            zt_i_new       = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i) 
    491471             
    492             zEi            = cpic * ( zt_i_new - (ztmelts+rt0) ) &     ! Specific enthalpy of forming ice (J/kg, <0)       
    493                &               - lfus * ( 1.0 - ztmelts / ( zt_i_new - rt0 ) )   & 
    494                &               + rcp  * ztmelts           
    495              
    496             zEw            = rcp  * ( t_bo_1d(ji) - rt0 )         ! Specific enthalpy of seawater (J/kg, < 0) 
    497              
    498             zdE            = zEi - zEw                           ! Specific enthalpy difference (J/kg, <0) 
    499              
    500             ! Contribution to heat flux to the ocean [W.m-2], >0   
    501             hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 
    502  
    503             ! Total heat flux used in this process [W.m-2], <0   
    504             hfx_bog_1d(ji) = hfx_bog_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice 
    505              
    506             ! Contribution to salt flux, <0 
    507             sfx_bog_1d(ji) = sfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bott(ji) * s_i_new(ji) * r1_rdtice 
    508  
    509             ! Contribution to mass flux, <0 
    510             wfx_bog_1d(ji) =  wfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bott(ji) * r1_rdtice 
     472            zEi            = cpic * ( zt_i_new - (ztmelts+rt0) ) &                                    ! Specific enthalpy of forming ice (J/kg, <0) 
     473               &           - lfus * ( 1.0 - ztmelts / ( zt_i_new - rt0 ) ) + rcp  * ztmelts           
     474             
     475            zEw            = rcp  * ( t_bo_1d(ji) - rt0 )                                             ! Specific enthalpy of seawater (J/kg, < 0) 
     476             
     477            zdE            = zEi - zEw                                                                ! Specific enthalpy difference (J/kg, <0) 
     478             
     479            hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice                           ! Heat flux to the ocean [W.m-2], >0 
     480            hfx_bog_1d(ji) = hfx_bog_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice                           ! Heat flux used in this process [W.m-2], <0 
     481             
     482            sfx_bog_1d(ji) = sfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bott(ji) * s_i_new(ji) * r1_rdtice   ! Salt flux, <0 
     483 
     484            wfx_bog_1d(ji) = wfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bott(ji) * r1_rdtice                 ! Mass flux, <0 
    511485 
    512486            ! update heat content (J.m-2) and layer thickness 
     
    518492      END DO 
    519493 
    520       !---------------- 
    521       ! 4.2 Basal melt 
    522       !---------------- 
     494      ! Ice Basal melt 
     495      !--------------- 
    523496      zdeltah(1:npti,:) = 0._wp ! important 
    524497      DO jk = nlay_i, 1, -1 
     
    528501               ztmelts = - tmut * sz_i_1d(ji,jk)  ! Melting point of layer jk (C) 
    529502 
    530                IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN !!! Internal melting 
     503               IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN   !-- Internal melting 
    531504 
    532505                  zEi               = - e_i_1d(ji,jk) * r1_rhoic    ! Specific enthalpy of melting ice (J/kg, <0) 
    533506                  zdE               = 0._wp                         ! Specific enthalpy difference   (J/kg, <0) 
    534                                                                     ! set up at 0 since no energy is needed to melt water...(it is already melted) 
     507                                                                    !    set up at 0 since no energy is needed to melt water...(it is already melted) 
    535508                  zdeltah   (ji,jk) = MIN( 0._wp , - zh_i(ji,jk) )  ! internal melting occurs when the internal temperature is above freezing      
    536509                                                                    ! this should normally not happen, but sometimes, heat diffusion leads to this 
     
    540513                  zfmdt             = - zdeltah(ji,jk) * rhoic      ! Mass flux x time step > 0 
    541514 
    542                   ! Contribution to heat flux to the ocean [W.m-2], <0 (ice enthalpy zEi is "sent" to the ocean)  
    543                   hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice 
    544  
    545                   ! Contribution to salt flux (using s_i_1d and not sz_i_1d(jk) is ok) 
    546                   sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice 
    547                                      
    548                   ! Contribution to mass flux 
    549                   wfx_res_1d(ji) =  wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 
     515                  hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice                           ! Heat flux to the ocean [W.m-2], <0 
     516                  !                                                                                                  ice enthalpy zEi is "sent" to the ocean 
     517                  sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice   ! Salt flux 
     518                  !                                                                                                  using s_i_1d and not sz_i_1d(jk) is ok 
     519                  wfx_res_1d(ji) = wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice                ! Mass flux 
    550520 
    551521                  ! update heat content (J.m-2) and layer thickness 
     
    553523                  h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk) 
    554524 
    555                ELSE                               !!! Basal melting 
    556  
    557                   zEi             = - e_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0) 
    558                   zEw             = rcp * ztmelts              ! Specific enthalpy of meltwater (J/kg, <0) 
    559                   zdE             = zEi - zEw                  ! Specific enthalpy difference   (J/kg, <0) 
    560  
    561                   zfmdt           = - zq_bo(ji) / zdE          ! Mass flux x time step (kg/m2, >0) 
    562  
    563                   zdeltah(ji,jk)  = - zfmdt * r1_rhoic         ! Gross thickness change 
    564  
    565                   zdeltah(ji,jk)  = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) ) ! bound thickness change 
     525               ELSE                                        !-- Basal melting 
     526 
     527                  zEi             = - e_i_1d(ji,jk) * r1_rhoic                                ! Specific enthalpy of melting ice (J/kg, <0) 
     528                  zEw             = rcp * ztmelts                                             ! Specific enthalpy of meltwater (J/kg, <0) 
     529                  zdE             = zEi - zEw                                                 ! Specific enthalpy difference   (J/kg, <0) 
     530 
     531                  zfmdt           = - zq_bo(ji) / zdE                                         ! Mass flux x time step (kg/m2, >0) 
     532 
     533                  zdeltah(ji,jk)  = - zfmdt * r1_rhoic                                        ! Gross thickness change 
     534 
     535                  zdeltah(ji,jk)  = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) )       ! bound thickness change 
    566536                   
    567                   zq_bo(ji)       = MAX( 0._wp , zq_bo(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat. MAX is necessary for roundup errors 
    568  
    569                   dh_i_bott(ji)   = dh_i_bott(ji) + zdeltah(ji,jk)    ! Update basal melt 
    570  
    571                   zfmdt           = - zdeltah(ji,jk) * rhoic          ! Mass flux x time step > 0 
    572  
    573                   zQm             = zfmdt * zEw         ! Heat exchanged with ocean 
    574  
    575                   ! Contribution to heat flux to the ocean [W.m-2], <0   
    576                   hfx_thd_1d(ji)  = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 
    577  
    578                   ! Contribution to salt flux (using s_i_1d and not sz_i_1d(jk) is ok) 
    579                   sfx_bom_1d(ji)  = sfx_bom_1d(ji) - rhoic *  a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice 
     537                  zq_bo(ji)       = MAX( 0._wp , zq_bo(ji) - zdeltah(ji,jk) * rhoic * zdE )   ! update available heat. MAX is necessary for roundup errors 
     538 
     539                  dh_i_bott(ji)   = dh_i_bott(ji) + zdeltah(ji,jk)                            ! Update basal melt 
     540 
     541                  zfmdt           = - zdeltah(ji,jk) * rhoic                                  ! Mass flux x time step > 0 
     542 
     543                  zQm             = zfmdt * zEw                                               ! Heat exchanged with ocean 
     544 
     545                  hfx_thd_1d(ji)  = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice                           ! Heat flux to the ocean [W.m-2], <0   
     546                  hfx_bom_1d(ji)  = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice                           ! Heat used in this process [W.m-2], >0   
     547 
     548                  sfx_bom_1d(ji)  = sfx_bom_1d(ji) - rhoic *  a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice  ! Salt flux 
     549                  !                                                                                                   using s_i_1d and not sz_i_1d(jk) is ok 
    580550                   
    581                   ! Total heat flux used in this process [W.m-2], >0   
    582                   hfx_bom_1d(ji)  = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice 
    583                    
    584                   ! Contribution to mass flux 
    585                   wfx_bom_1d(ji)  =  wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 
     551                  wfx_bom_1d(ji)  = wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice                ! Mass flux 
    586552 
    587553                  ! update heat content (J.m-2) and layer thickness 
     
    594560      END DO 
    595561 
    596       !------------------------------------------- 
    597562      ! Update temperature, energy 
    598       !------------------------------------------- 
    599       DO ji = 1, npti 
    600          h_i_1d(ji) =  MAX( 0._wp , h_i_1d(ji) + dh_i_bott(ji) ) 
     563      ! -------------------------- 
     564      DO ji = 1, npti 
     565         h_i_1d(ji) = MAX( 0._wp , h_i_1d(ji) + dh_i_bott(ji) ) 
    601566      END DO   
    602567 
    603       !------------------------------------------- 
    604       ! 5. What to do with remaining energy 
    605       !------------------------------------------- 
    606       ! If heat still available for melting and snow remains, then melt more snow 
     568      ! If heat still available then melt more snow 
    607569      !------------------------------------------- 
    608570      zdeltah(1:npti,:) = 0._wp ! important 
    609571      DO ji = 1, npti 
    610          zq_rema(ji)     = zq_su(ji) + zq_bo(ji)  
    611          rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, - h_s_1d(ji) ) )   ! =1 if snow 
    612          rswitch         = rswitch * MAX( 0._wp, SIGN( 1._wp, e_s_1d(ji,1) - epsi20 ) ) 
    613          zdeltah  (ji,1) = - rswitch * zq_rema(ji) / MAX( e_s_1d(ji,1), epsi20 ) 
    614          zdeltah  (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - h_s_1d(ji) ) ) ! bound melting 
    615          dh_s_tot (ji)   = dh_s_tot(ji) + zdeltah(ji,1) 
    616          h_s_1d   (ji)   = h_s_1d(ji)  + zdeltah(ji,1) 
     572         zq_rema (ji)   = zq_su(ji) + zq_bo(ji)  
     573         rswitch        = 1._wp - MAX( 0._wp, SIGN( 1._wp, - h_s_1d(ji) ) )   ! =1 if snow 
     574         rswitch        = rswitch * MAX( 0._wp, SIGN( 1._wp, e_s_1d(ji,1) - epsi20 ) ) 
     575         zdeltah (ji,1) = - rswitch * zq_rema(ji) / MAX( e_s_1d(ji,1), epsi20 ) 
     576         zdeltah (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - h_s_1d(ji) ) ) ! bound melting 
     577         dh_s_tot(ji)   = dh_s_tot(ji) + zdeltah(ji,1) 
     578         h_s_1d  (ji)   = h_s_1d  (ji) + zdeltah(ji,1) 
    617579         
    618          zq_rema(ji)     = zq_rema(ji) + zdeltah(ji,1) * e_s_1d(ji,1)                ! update available heat (J.m-2) 
    619          ! heat used to melt snow 
    620          hfx_snw_1d(ji)  = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * e_s_1d(ji,1) * r1_rdtice ! W.m-2 (>0) 
    621          ! Contribution to mass flux 
    622          wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice 
     580         zq_rema(ji)        = zq_rema(ji) + zdeltah(ji,1) * e_s_1d(ji,1)                               ! update available heat (J.m-2) 
     581         hfx_snw_1d(ji)     = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * e_s_1d(ji,1) * r1_rdtice   ! Heat used to melt snow, W.m-2 (>0) 
     582         wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice      ! Mass flux 
    623583         dh_s_mlt(ji)       = dh_s_mlt(ji) + zdeltah(ji,1) 
    624584         !     
    625585         ! Remaining heat flux (W.m-2) is sent to the ocean heat budget 
    626          hfx_out_1d(ji)  = hfx_out_1d(ji) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_rdtice 
     586         hfx_out_1d(ji) = hfx_out_1d(ji) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_rdtice 
    627587 
    628588         IF( ln_icectl .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji) 
     
    630590 
    631591      ! 
    632       !------------------------------------------------------------------------------| 
    633       !  6) Snow-Ice formation                                                       | 
    634       !------------------------------------------------------------------------------| 
     592      ! Snow-Ice formation 
     593      ! ------------------ 
    635594      ! When snow load excesses Archimede's limit, snow-ice interface goes down under sea-level,  
    636595      ! flooding of seawater transforms snow into ice dh_snowice is positive for the ice 
     
    648607         zQm            = zfmdt * zEw  
    649608          
    650          ! Contribution to heat flux 
    651          hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice  
    652  
    653          ! Contribution to salt flux 
    654          sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_1d(ji) * a_i_1d(ji) * zfmdt * r1_rdtice  
     609         hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice ! Heat flux 
     610 
     611         sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_1d(ji) * a_i_1d(ji) * zfmdt * r1_rdtice ! Salt flux 
    655612 
    656613         ! Case constant salinity in time: virtual salt flux to keep salinity constant 
     
    660617         ENDIF 
    661618 
    662          ! Contribution to mass flux 
    663          ! All snow is thrown in the ocean, and seawater is taken to replace the volume 
     619         ! Mass flux: All snow is thrown in the ocean, and seawater is taken to replace the volume 
    664620         wfx_sni_1d(ji)     = wfx_sni_1d(ji)     - a_i_1d(ji) * dh_snowice(ji) * rhoic * r1_rdtice 
    665621         wfx_snw_sni_1d(ji) = wfx_snw_sni_1d(ji) + a_i_1d(ji) * dh_snowice(ji) * rhosn * r1_rdtice 
     
    672628 
    673629      ! 
    674       !------------------------------------------- 
    675630      ! Update temperature, energy 
    676       !------------------------------------------- 
    677       DO ji = 1, npti 
    678          rswitch     =  1.0 - MAX( 0._wp , SIGN( 1._wp , - h_i_1d(ji) ) )  
    679          t_su_1d(ji) =  rswitch * t_su_1d(ji) + ( 1.0 - rswitch ) * rt0 
     631      ! -------------------------- 
     632      DO ji = 1, npti 
     633         rswitch     = 1._wp - MAX( 0._wp , SIGN( 1._wp , - h_i_1d(ji) ) )  
     634         t_su_1d(ji) = rswitch * t_su_1d(ji) + ( 1._wp - rswitch ) * rt0 
    680635      END DO 
    681636 
     
    698653   !!-------------------------------------------------------------------------- 
    699654   !! INTERFACE ice_thd_snwblow 
     655   !! 
    700656   !! ** Purpose :   Compute distribution of precip over the ice 
     657   !! 
     658   !!                Snow accumulation in one thermodynamic time step 
     659   !!                snowfall is partitionned between leads and ice. 
     660   !!                If snow fall was uniform, a fraction (1-at_i) would fall into leads 
     661   !!                but because of the winds, more snow falls on leads than on sea ice 
     662   !!                and a greater fraction (1-at_i)^beta of the total mass of snow  
     663   !!                (beta < 1) falls in leads. 
     664   !!                In reality, beta depends on wind speed,  
     665   !!                and should decrease with increasing wind speed but here, it is  
     666   !!                considered as a constant. an average value is 0.66 
    701667   !!-------------------------------------------------------------------------- 
    702668!!gm  I think it can be usefull to set this as a FUNCTION, not a SUBROUTINE.... 
Note: See TracChangeset for help on using the changeset viewer.