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 4506 for branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90 – NEMO

Ignore:
Timestamp:
2014-02-21T15:20:39+01:00 (10 years ago)
Author:
vancop
Message:

[Heat conservation in LIM3, part HC1 (LIM_SRC_3_HC17)]

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90

    r4332 r4506  
    7474      INTEGER  ::   ji , jk        ! dummy loop indices 
    7575      INTEGER  ::   ii, ij         ! 2D corresponding indices to ji 
     76      INTEGER  ::   i_use          ! debugging flag 
    7677      INTEGER  ::   isnow          ! switch for presence (1) or absence (0) of snow 
    7778      INTEGER  ::   isnowic        ! snow ice formation not 
     
    8586      REAL(wp) ::   zfracs       ! fractionation coefficient for bottom salt entrapment 
    8687      REAL(wp) ::   zcoeff       ! dummy argument for snowfall partitioning over ice and leads 
    87       REAL(wp) ::   zsm_snowice  ! snow-ice salinity 
     88      REAL(wp) ::   zs_snic  ! snow-ice salinity 
    8889      REAL(wp) ::   zswi1        ! switch for computation of bottom salinity 
    8990      REAL(wp) ::   zswi12       ! switch for computation of bottom salinity 
    9091      REAL(wp) ::   zswi2        ! switch for computation of bottom salinity 
    9192      REAL(wp) ::   zgrr         ! bottom growth rate 
    92       REAL(wp) ::   ztform       ! bottom formation temperature 
    93       ! 
     93      REAL(wp) ::   zt_i_new     ! bottom formation temperature 
     94 
     95      REAL(wp) ::   zQm          ! enthalpy exchanged with the ocean (J/m2), >0 towards the ocean 
     96      REAL(wp) ::   zEi          ! specific enthalpy of sea ice (J/kg) 
     97      REAL(wp) ::   zEw          ! specific enthalpy of exchanged water (J/kg) 
     98      REAL(wp) ::   zdE          ! specific enthalpy difference (J/kg) 
     99      REAL(wp) ::   zfmdt        ! exchange mass flux x time step (J/m2), >0 towards the ocean 
     100      REAL(wp) ::   zsstK        ! SST in Kelvin 
     101 
    94102      REAL(wp), POINTER, DIMENSION(:) ::   zh_i        ! ice layer thickness 
    95103      REAL(wp), POINTER, DIMENSION(:) ::   zh_s        ! snow layer thickness 
     
    119127 
    120128      ! mass and salt flux (clem) 
    121       REAL(wp) :: zdvres, zdvsur, zdvbot 
     129      REAL(wp) :: zdvres, zdvsur, zdvbot, zswitch_sal 
    122130      REAL(wp), POINTER, DIMENSION(:) ::   zviold, zvsold   ! old ice volume... 
    123131 
     
    129137      REAL(wp), POINTER, DIMENSION(:) ::   zfbase, zdq_i 
    130138      !!------------------------------------------------------------------ 
     139 
     140      ! Discriminate between varying salinity (num_sal=2) and prescribed cases (other values) 
     141      SELECT CASE( num_sal )                       ! varying salinity or not 
     142         CASE( 1, 3, 4 ) ;   zswitch_sal = 0       ! prescribed salinity profile 
     143         CASE( 2 )       ;   zswitch_sal = 1       ! varying salinity profile 
     144      END SELECT 
    131145 
    132146      CALL wrk_alloc( jpij, zh_i, zh_s, ztfs, zhsold, zqprec, zqfont_su, zqfont_bo, z_f_surf, zhgnew, zfmass_i ) 
     
    222236      DO ji = kideb, kiut 
    223237         ! tatm_ice is now in K 
    224          zqprec   (ji)   =  rhosn * ( cpic * ( rtt - tatm_ice_1d(ji) ) + lfus )   
     238         zqprec   (ji)   =  rhosn * ( cpic * ( rtt - MIN( tatm_ice_1d(ji), rt0_snow) ) + lfus )   
    225239         zqfont_su(ji)   =  z_f_surf(ji) * rdt_ice 
    226240         zdeltah  (ji,1) =  MIN( 0.e0 , - zqfont_su(ji) / MAX( zqprec(ji) , epsi13 ) ) 
     
    272286      DO jk = 1, nlay_i 
    273287         DO ji = kideb, kiut  
    274             !                                                    ! melt of layer jk 
    275             zdeltah  (ji,jk) = - zqfont_su(ji) / q_i_b(ji,jk) 
    276             !                                                    ! recompute heat available 
    277             zqfont_su(ji   ) = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * q_i_b(ji,jk)  
    278             !                                                    ! melt of layer jk cannot be higher than its thickness 
    279             zdeltah  (ji,jk) = MAX( zdeltah(ji,jk) , - zh_i(ji) ) 
    280             !                                                    ! update surface melt 
    281             dh_i_surf(ji   ) = dh_i_surf(ji) + zdeltah(ji,jk)  
    282             !                                                    ! for energy conservation 
    283             zdq_i    (ji   ) = zdq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk) * r1_rdtice 
    284             ! 
    285             ! clem 
    286             sfx_thd_1d(ji) = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji)    & 
     288            zEi            = - q_i_b(ji,jk) / rhoic                ! Specific enthalpy of layer k [J/kg, <0] 
     289 
     290            ztmelts        = - tmut * s_i_b(ji,jk) + rtt           ! Melting point of layer k [K] 
     291 
     292            zEw            =    rcp * ( ztmelts - rt0 )            ! Specific enthalpy of resulting meltwater [J/kg, <0] 
     293 
     294            zdE            =    zEi - zEw                          ! Specific enthalpy difference < 0 
     295 
     296            zfmdt          = - zqfont_su(ji) / zdE                 ! Mass flux to the ocean [kg/m2, >0] 
     297 
     298            zdeltah(ji,jk) = - zfmdt / rhoic                       ! Melt of layer jk [m, <0] 
     299 
     300            zqfont_su(ji)  = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * rhoic * ( - zdE ) 
     301                                                                   ! Energy remaining in case of melting of the full layer [J/m2, >0] 
     302            zdeltah(ji,jk) = MAX( zdeltah(ji,jk) , - zh_i(ji) )    ! Melt of layer jk cannot exceed the layer thickness [m, <0] 
     303 
     304            dh_i_surf(ji)  = dh_i_surf(ji) + zdeltah(ji,jk)        ! Cumulate surface melt 
     305 
     306            zfmdt          = - rhoic*zdeltah(ji,jk)                ! Recompute mass flux [kg/m2, >0] 
     307 
     308            zQm            = MAX ( zfmdt, 0._wp ) * zEw            ! Energy of the melt water sent to the ocean [J/m2, <0] 
     309 
     310            ! Contribution to salt flux  
     311            sfx_thd_1d(ji)   = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji)    & 
    287312               &                              * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic / rdt_ice 
     313            ! Contribution to heat flux 
     314            rdq_ice_1d(ji) = rdq_ice_1d(ji) + zQm * a_i_b(ji) 
     315             
     316            ! Conservation test 
     317            zdq_i(ji)      = zdq_i(ji) + zdeltah(ji,jk) * rhoic * zdE * r1_rdtice  ! ? still don't know 
    288318         END DO 
    289319      END DO 
     
    364394      !------------------------------------------------------------------------------! 
    365395      ! 
    366       ! Ice basal growth / melt is given by the ratio of heat budget over basal 
    367       ! ice heat content.  Basal heat budget is given by the difference between 
    368       ! the inner conductive flux  (fc_bo_i), from the open water heat flux  
     396      !------------------ 
     397      ! 4.1 Basal growth  
     398      !------------------ 
     399      ! Basal growth is driven by heat imbalance at the ice-ocean interface, 
     400      ! between the inner conductive flux  (fc_bo_i), from the open water heat flux  
    369401      ! (qlbbqb) and the turbulent ocean flux (fbif).  
    370402      ! fc_bo_i is positive downwards. fbif and qlbbq are positive to the ice  
    371403 
    372       !----------------------------------------------------- 
    373       ! 4.1 Basal growth - (a) salinity not varying in time  
    374       !----------------------------------------------------- 
    375       IF(  num_sal /= 2  ) THEN   ! ice salinity constant in time 
     404      ! If salinity varies in time, an iterative procedure is required, because 
     405      ! the involved quantities are inter-dependent. 
     406      ! Basal growth (dh_i_bott) depends upon new ice specific enthalpy (zEi), 
     407      ! which depends on forming ice salinity (s_i_new), which depends on dh/dt (dh_i_bott) 
     408      ! -> need for an iterative procedure, which converges quickly 
     409 
     410      IF ( num_sal == 2 ) THEN 
     411         num_iter_max = 5 
     412      ELSE 
     413         num_iter_max = 1 
     414      ENDIF 
     415 
     416      ! Initialize dh_i_bott 
     417      dh_i_bott(:) = 0.0e0 
     418 
     419      ! Iterative procedure 
     420      DO iter = 1, num_iter_max 
    376421         DO ji = kideb, kiut 
    377             IF(  ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) < 0._wp  ) THEN 
    378                s_i_new(ji)         =  sm_i_b(ji) 
    379                ! Melting point in K 
    380                ztmelts             =  - tmut * s_i_new(ji) + rtt  
    381                ! New ice heat content (Bitz and Lipscomb, 1999) 
    382                ztform              =  t_i_b(ji,nlay_i)  ! t_bo_b crashes in the 
    383                ! Baltic 
    384                q_i_b(ji,nlay_i+1)  = rhoic * (  cpic * ( ztmelts - ztform )                                & 
    385                   &                           + lfus * (  1.0 - ( ztmelts - rtt ) / ( ztform - rtt )  )    & 
    386                   &                           - rcp  * ( ztmelts - rtt )                                 ) 
    387                ! Basal growth rate = - F*dt / q 
    388                dh_i_bott(ji)       =  - rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1)  
    389                sfx_thd_1d(ji) = sfx_thd_1d(ji) - s_i_new(ji) * a_i_b(ji) * dh_i_bott(ji) * rhoic * r1_rdtice 
    390             ENDIF 
    391          END DO 
    392       ENDIF 
    393  
    394       !------------------------------------------------- 
    395       ! 4.1 Basal growth - (b) salinity varying in time  
    396       !------------------------------------------------- 
    397       IF(  num_sal == 2  ) THEN 
    398          ! the growth rate (dh_i_bott) is function of the new ice heat content (q_i_b(nlay_i+1)).  
    399          ! q_i_b depends on the new ice salinity (snewice).  
    400          ! snewice depends on dh_i_bott ; it converges quickly, so, no problem 
    401          ! See Vancoppenolle et al., OM08 for more info on this 
    402  
    403          ! Initial value (tested 1D, can be anything between 1 and 20) 
    404          num_iter_max = 4 
    405          s_i_new(:)   = 4.0 
    406  
    407          ! Iterative procedure 
    408          DO iter = 1, num_iter_max 
    409             DO ji = kideb, kiut 
    410                IF(  fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) < 0.e0  ) THEN 
    411                   ii = MOD( npb(ji) - 1, jpi ) + 1 
    412                   ij = ( npb(ji) - 1 ) / jpi + 1 
    413                   ! Melting point in K 
    414                   ztmelts             =   - tmut * s_i_new(ji) + rtt  
    415                   ! New ice heat content (Bitz and Lipscomb, 1999) 
    416                   q_i_b(ji,nlay_i+1)  =  rhoic * (  cpic * ( ztmelts - t_bo_b(ji) )                             & 
    417                      &                            + lfus * ( 1.0 - ( ztmelts - rtt ) / ( t_bo_b(ji) - rtt ) )   & 
    418                      &                            - rcp * ( ztmelts-rtt )                                     ) 
    419                   ! Bottom growth rate = - F*dt / q 
    420                   dh_i_bott(ji) =  - rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 
    421                   ! New ice salinity ( Cox and Weeks, JGR, 1988 ) 
    422                   ! zswi2  (1) if dh_i_bott/rdt .GT. 3.6e-7 
    423                   ! zswi12 (1) if dh_i_bott/rdt .LT. 3.6e-7 and .GT. 2.0e-8 
    424                   ! zswi1  (1) if dh_i_bott/rdt .LT. 2.0e-8 
    425                   zgrr   = MIN( 1.0e-3, MAX ( dh_i_bott(ji) * r1_rdtice , epsi13 ) ) 
    426                   zswi2  = MAX( zzero , SIGN( zone , zgrr - 3.6e-7 ) )  
    427                   zswi12 = MAX( zzero , SIGN( zone , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 ) 
    428                   zswi1  = 1. - zswi2 * zswi12  
    429                   zfracs = zswi1  * 0.12 + zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) )   & 
    430                      &                   + zswi2  * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) )  
    431                   zfracs = MIN( 0.5 , zfracs ) 
    432                   s_i_new(ji) = zfracs * sss_m(ii,ij) 
    433                ENDIF ! fc_bo_i 
    434             END DO ! ji 
    435          END DO ! iter 
    436  
    437          ! Final values 
    438          DO ji = kideb, kiut 
    439             IF( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .LT. 0.0  ) THEN 
    440                ! New ice salinity must not exceed 20 psu 
    441                s_i_new(ji) = MIN( s_i_new(ji), s_i_max ) 
    442                ! Metling point in K 
    443                ztmelts     =   - tmut * s_i_new(ji) + rtt  
    444                ! New ice heat content (Bitz and Lipscomb, 1999) 
    445                q_i_b(ji,nlay_i+1)  =  rhoic * (  cpic * ( ztmelts - t_bo_b(ji) )                              & 
    446                   &                            + lfus * ( 1.0 - ( ztmelts - rtt ) / ( t_bo_b(ji) - rtt ) )    & 
    447                   &                            - rcp * ( ztmelts - rtt )                                    ) 
    448                ! Basal growth rate = - F*dt / q 
    449                dh_i_bott(ji)       =  - rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 
    450                ! Salinity update 
    451                ! entrapment during bottom growth 
    452                sfx_thd_1d(ji) = sfx_thd_1d(ji) - s_i_new(ji) * a_i_b(ji) * dh_i_bott(ji) * rhoic * r1_rdtice 
    453             ENDIF ! heat budget 
    454          END DO 
    455       ENDIF 
     422 
     423            IF(  fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) < 0.e0  ) THEN 
     424 
     425               ! New bottom ice salinity (Cox & Weeks, JGR88 ) 
     426               !--- zswi1  if dh/dt < 2.0e-8 
     427               !--- zswi12 if 2.0e-8 < dh/dt < 3.6e-7  
     428               !--- zswi2  if dh/dt > 3.6e-7 
     429               zgrr               = MIN( 1.0e-3, MAX ( dh_i_bott(ji) * r1_rdtice , epsi13 ) ) 
     430               zswi2              = MAX( zzero , SIGN( zone , zgrr - 3.6e-7 ) ) 
     431               zswi12             = MAX( zzero , SIGN( zone , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 ) 
     432               zswi1              = 1. - zswi2 * zswi12 
     433               zfracs             = MIN ( zswi1  * 0.12 + zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) )   & 
     434                  &               + zswi2  * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) )  , 0.5 ) 
     435 
     436               ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 
     437 
     438               s_i_new(ji)        = zswitch_sal * zfracs * sss_m(ii,ij)  &  ! New ice salinity 
     439                                  + ( 1. - zswitch_sal ) * sm_i_b(ji)  
     440               ! New ice growth 
     441               ztmelts            = - tmut * s_i_new(ji) + rtt          ! New ice melting point (K) 
     442 
     443               zt_i_new           = zswitch_sal * t_bo_b(ji) + ( 1. - zswitch_sal) * t_i_b(ji, nlay_i) 
     444                
     445               zEi                = cpic * ( zt_i_new - ztmelts ) &     ! Specific enthalpy of forming ice (J/kg, <0)       
     446                  &               - lfus * ( 1.0 - ( ztmelts - rtt ) / ( zt_i_new - rtt ) )   & 
     447                  &               + rcp  * ( ztmelts-rtt )           
     448 
     449               zEw                = rcp  * ( t_bo_b(ji) - rt0 )         ! Specific enthalpy of seawater (J/kg, < 0) 
     450 
     451               zdE                = zEi - zEw                           ! Specific enthalpy difference (J/kg, <0) 
     452 
     453               dh_i_bott(ji)      = rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / ( zdE * rhoic ) 
     454 
     455               !!! not sure we still need the next line... useful to keep this in memory ? check limthd_ent... 
     456               q_i_b(ji,nlay_i+1) = -zEi * rhoic                        ! New ice energy of melting (J/m3, >0) 
     457 
     458            ENDIF ! fc_bo_i 
     459         END DO ! ji 
     460      END DO ! iter 
     461 
     462      ! Contribution to Energy and Salt Fluxes 
     463      DO ji = kideb, kiut 
     464 
     465         zEw            = rcp  * ( t_bo_b(ji) - rt0 )                   ! Specific enthalpy of seawater (J/kg, < 0) 
     466 
     467         zfmdt          = - rhoic * dh_i_bott(ji)                       ! Mass flux x time step (kg/m2, < 0) 
     468 
     469         ! Contribution to energy flux to the ocean (J/m2) 
     470         rdq_ice_1d(ji) = rdq_ice_1d(ji) + zEw         * a_i_b(ji) * zfmdt 
     471 
     472         ! Contribution to salt flux  () 
     473         sfx_thd_1d(ji) = sfx_thd_1d(ji) + s_i_new(ji) * a_i_b(ji) * zfmdt * r1_rdtice 
     474 
     475      END DO 
    456476 
    457477      !---------------- 
     
    465485         ! heat convergence at the surface > 0 
    466486         IF(  ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) >= 0._wp  ) THEN 
    467             s_i_new(ji)   =  s_i_b(ji,nlay_i) 
     487            s_i_new(ji) = s_i_b(ji,nlay_i) ! is this line still useful ? 
    468488            zqfont_bo(ji) =  rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) 
    469489            zfbase(ji)    =  zqfont_bo(ji) * r1_rdtice     ! heat conservation test 
     
    476496         DO ji = kideb, kiut 
    477497            IF(  fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji)  >=  0._wp  ) THEN 
    478                ztmelts = - tmut * s_i_b(ji,jk) + rtt  
    479                IF( t_i_b(ji,jk) >= ztmelts ) THEN   !!gm : a comment is needed 
    480                   zdeltah   (ji,jk) = - zh_i(ji) 
    481                   dh_i_bott (ji   ) = dh_i_bott(ji) + zdeltah(ji,jk) 
    482                   zinnermelt(ji   ) = 1._wp 
    483                ELSE                                  ! normal ablation 
    484                   zdeltah  (ji,jk) = - zqfont_bo(ji) / q_i_b(ji,jk) 
    485                   zqfont_bo(ji   ) = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * q_i_b(ji,jk) 
    486                   zdeltah  (ji,jk) = MAX(zdeltah(ji,jk), - zh_i(ji) ) 
    487                   dh_i_bott(ji   ) = dh_i_bott(ji) + zdeltah(ji,jk) 
    488                   zdq_i    (ji   ) = zdq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk) * r1_rdtice 
     498 
     499               ztmelts = - tmut * s_i_b(ji,jk) + rtt  ! Melting point of layer jk (K) 
     500 
     501               IF( t_i_b(ji,jk) >= ztmelts ) THEN !!! Internal melting 
     502 
     503                  zdeltah   (ji,jk) = - zh_i(ji)  ! internal melting occurs when the internal temperature is above freezing      
     504                                                  ! this should normally not happen, but sometimes, heat diffusion leads to this 
     505 
     506                  dh_i_bott (ji)    = dh_i_bott(ji) + zdeltah(ji,jk) 
     507 
     508                  zinnermelt(ji)    = 1._wp 
     509 
     510                  zQm               = 0.          ! Not sure which specific enthalpy we should use here (MV HC 2014) 
     511                                                  ! If that happens, heat is probably not well counted 
     512                                                  ! Put zero by default, but more bug analysis should be done to investigate this case 
     513 
     514               ELSE                               !!! Basal melting 
     515 
     516                  zEi               = - q_i_b(ji,jk) / rhoic ! Specific enthalpy of melting ice (J/kg, <0) 
     517 
     518                  zEw               = rcp * ( ztmelts - rtt )! Specific enthalpy of meltwater (J/kg, <0) 
     519 
     520                  zdE               = zEi - zEw              ! Specific enthalpy difference   (J/kg, <0) 
     521 
     522                  zfmdt             = - zqfont_bo(ji) / zdE  ! Mass flux x time step (kg/m2, >0) 
     523 
     524                  zdeltah(ji,jk)    = - zfmdt / rhoic        ! Gross thickness change 
     525 
     526                  zqfont_bo(ji)     = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * rhoic * ( - zdE ) ! Update heat available 
     527 
     528                  zdeltah(ji,jk)    = MAX( zdeltah(ji,jk), - zh_i(ji) ) ! Update thickness change 
     529 
     530                  dh_i_bott(ji)     = dh_i_bott(ji) + zdeltah(ji,jk)    ! Update basal melt 
     531 
     532                  zfmdt             = - zdeltah(ji,jk) * rhoic          ! Mass flux x time step 
     533 
     534                  zQm               = MAX ( zfmdt , 0.0 ) * zEw         ! Heat exchanged with ocean 
     535 
     536                  zdq_i(ji)         = zdq_i(ji) + zdeltah(ji,jk) * rhoic * zdE * r1_rdtice ! for heat conservation  
     537 
    489538               ENDIF 
    490                ! clem: contribution to salt flux 
     539 
     540               ! Contribution to salt flux 
    491541               sfx_thd_1d(ji) = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji)    & 
    492542                    &                              * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic * r1_rdtice 
     543 
     544               ! Contribution to energy flux to the ocean (J/m2) 
     545               rdq_ice_1d(ji) = rdq_ice_1d(ji) + zQm * a_i_b(ji) 
     546 
    493547            ENDIF 
    494548         END DO ! ji 
     
    571625         dvbbq_1d (ji) = a_i_b(ji) * dh_i_bott(ji) 
    572626 
    573         ! remaining heat 
     627         ! remaining heat 
    574628         zfdt_final(ji) = ( 1.0 - zihgnew ) * ( zqfont_su(ji) +  zqfont_bo(ji) )  
    575629 
     
    661715      DO ji = kideb, kiut 
    662716         ht_i_b(ji) = zhgnew(ji) 
    663       END DO  ! ji 
     717      END DO ! ji  
     718       
    664719      ! 
    665720      !------------------------------------------------------------------------------| 
     
    681736         !clem rdm_snw_1d(ji) = rdm_snw_1d(ji) + a_i_b(ji) * ( zhnnew     - ht_s_b(ji) ) * rhosn  
    682737 
    683          !        Equivalent salt flux (1) Snow-ice formation component 
    684          !        ----------------------------------------------------- 
    685          ii = MOD( npb(ji) - 1, jpi ) + 1 
    686          ij =    ( npb(ji) - 1 ) / jpi + 1 
    687  
    688          IF( num_sal == 2 ) THEN   ;   zsm_snowice = sss_m(ii,ij) * ( rhoic - rhosn ) / rhoic 
    689          ELSE                      ;   zsm_snowice = sm_i_b(ji)    
    690          ENDIF 
     738         ! Salinity of snow ice 
     739         ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 
     740 
     741         zs_snic = zswitch_sal          * sss_m(ii,ij) * ( rhoic - rhosn ) / rhoic     & 
     742                 + ( 1. - zswitch_sal ) * sm_i_b(ji) 
     743 
    691744         ! entrapment during snow ice formation 
    692745         ! clem: new salinity difference stored (to be used in limthd_ent.F90) 
     
    694747            i_ice_switch = MAX( 0._wp , SIGN( 1._wp , zhgnew(ji) - epsi10 ) ) 
    695748            ! salinity dif due to snow-ice formation 
    696             dsm_i_si_1d(ji) = ( zsm_snowice - sm_i_b(ji) ) * dh_snowice(ji) / MAX( zhgnew(ji), epsi10 ) * i_ice_switch      
     749            dsm_i_si_1d(ji) = ( zs_snic - sm_i_b(ji) ) * dh_snowice(ji) / MAX( zhgnew(ji), epsi10 ) * i_ice_switch      
    697750            ! salinity dif due to bottom growth  
    698751            IF (  fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji)  < 0._wp ) THEN 
     
    709762 
    710763         ! diagnostic ( snow ice growth ) 
    711          ii = MOD( npb(ji) - 1, jpi ) + 1 
    712          ij =    ( npb(ji) - 1 ) / jpi + 1 
     764         ii = MOD( npb(ji) - 1, jpi ) + 1 !MV HC 2014 useless 
     765         ij =    ( npb(ji) - 1 ) / jpi + 1 ! MV HC 2014 useless 
    713766         diag_sni_gr(ii,ij)  = diag_sni_gr(ii,ij) + dh_snowice(ji)*a_i_b(ji) * r1_rdtice 
    714          ! 
    715          ! salt flux 
    716          sfx_thd_1d(ji) = sfx_thd_1d(ji) - zsm_snowice * a_i_b(ji) * dh_snowice(ji) * rhoic * r1_rdtice 
    717          !-------------------------------- 
    718          ! Update mass fluxes (clem) 
    719          !-------------------------------- 
     767 
     768         ! Contribution to energy flux to the ocean [J/m2], >0 
     769         zfmdt          = ( rhosn - rhoic ) * MAX( dh_snowice(ji), 0.0 )    ! <0 
     770         zsstK          = sst_m(ii,ij) + rt0                                 
     771         zEw            = rcp * ( zsstK - rt0 ) 
     772         zQm            = zfmdt * zEw  
     773         rdq_ice_1d(ji) = rdq_ice_1d(ji) + zQm * a_i_b(ji) 
     774 
     775         ! Contribution to salt flux 
     776         sfx_thd_1d(ji) = sfx_thd_1d(ji) + sss_m(ii,ij) * a_i_b(ji) * zfmdt * r1_rdtice  
     777           
     778         ! Contribution to mass fluxes 
    720779         rdm_ice_1d(ji) = rdm_ice_1d(ji) + ( a_i_b(ji) * ht_i_b(ji) - zviold(ji) ) * rhoic  
    721780         rdm_snw_1d(ji) = rdm_snw_1d(ji) + ( a_i_b(ji) * ht_s_b(ji) - zvsold(ji) ) * rhosn  
Note: See TracChangeset for help on using the changeset viewer.