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 9937 for NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icethd_dh.F90 – NEMO

Ignore:
Timestamp:
2018-07-12T17:55:41+02:00 (6 years ago)
Author:
gm
Message:

#1911 (ENHANCE-04): step I.2 (end): clean sea ice related physical constant in dev_r9838_ENHANCE04_MLF

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icethd_dh.F90

    r9923 r9937  
    7676      REAL(wp) ::   zgrr         ! bottom growth rate 
    7777      REAL(wp) ::   zt_i_new     ! bottom formation temperature 
    78       REAL(wp) ::   z1_rho       ! 1/(rhosn+rho0-rhoic) 
     78      REAL(wp) ::   z1_rho       ! 1/(rhos+rho0-rhoi) 
    7979 
    8080      REAL(wp) ::   zQm          ! enthalpy exchanged with the ocean (J/m2), >0 towards the ocean 
     
    182182            IF( t_s_1d(ji,jk) > rt0 ) THEN 
    183183               hfx_res_1d    (ji) = hfx_res_1d    (ji) + e_s_1d(ji,jk) * zh_s(ji,jk) * a_i_1d(ji) * r1_Dt_ice   ! heat flux to the ocean [W.m-2], < 0 
    184                wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + rhosn         * zh_s(ji,jk) * a_i_1d(ji) * r1_Dt_ice   ! mass flux 
     184               wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + rhos          * zh_s(ji,jk) * a_i_1d(ji) * r1_Dt_ice   ! mass flux 
    185185               ! updates 
    186186               dh_s_mlt(ji)    = dh_s_mlt(ji) - zh_s(ji,jk) 
     
    202202            ! 
    203203            ! --- precipitation --- 
    204             zdh_s_pre (ji) = zsnw(ji) * sprecip_1d(ji) * rdt_ice * r1_rhosn / at_i_1d(ji)   ! thickness change 
     204            zdh_s_pre (ji) = zsnw(ji) * sprecip_1d(ji) * rdt_ice * r1_rhos / at_i_1d(ji)    ! thickness change 
    205205            zqprec    (ji) = - qprec_ice_1d(ji)                                             ! enthalpy of the precip (>0, J.m-3) 
    206206            ! 
    207207            hfx_spr_1d(ji) = hfx_spr_1d(ji) + zdh_s_pre(ji) * a_i_1d(ji) * zqprec(ji)    * r1_Dt_ice   ! heat flux from snow precip (>0, W.m-2) 
    208             wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhosn         * a_i_1d(ji) * zdh_s_pre(ji) * r1_Dt_ice   ! mass flux, <0 
     208            wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhos          * a_i_1d(ji) * zdh_s_pre(ji) * r1_Dt_ice   ! mass flux, <0 
    209209             
    210210            ! --- melt of falling snow --- 
     
    213213            zdeltah       (ji,1) = MAX( - zdh_s_pre(ji), zdeltah(ji,1) )                ! bound melting  
    214214            hfx_snw_1d    (ji)   = hfx_snw_1d    (ji) - zdeltah(ji,1) * a_i_1d(ji) * zqprec(ji)    * r1_Dt_ice   ! heat used to melt snow (W.m-2, >0) 
    215             wfx_snw_sum_1d(ji)   = wfx_snw_sum_1d(ji) - rhosn         * a_i_1d(ji) * zdeltah(ji,1) * r1_Dt_ice   ! snow melting only = water into the ocean (then without snow precip), >0 
     215            wfx_snw_sum_1d(ji)   = wfx_snw_sum_1d(ji) - rhos          * a_i_1d(ji) * zdeltah(ji,1) * r1_Dt_ice   ! snow melting only = water into the ocean (then without snow precip), >0 
    216216             
    217217            ! updates available heat + precipitations after melting 
     
    253253                
    254254               hfx_snw_1d(ji)     = hfx_snw_1d(ji)     - zdeltah(ji,jk) * a_i_1d(ji) * e_s_1d (ji,jk) * r1_Dt_ice   ! heat used to melt snow(W.m-2, >0) 
    255                wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn          * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice   ! snow melting only = water into the ocean (then without snow precip) 
     255               wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhos           * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice   ! snow melting only = water into the ocean (then without snow precip) 
    256256                
    257257               ! updates available heat + thickness 
     
    273273         IF( evap_ice_1d(ji) > 0._wp ) THEN 
    274274            ! 
    275             zdh_s_sub (ji)   = MAX( - h_s_1d(ji) , - evap_ice_1d(ji) * r1_rhosn * rdt_ice ) 
    276             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) 
     275            zdh_s_sub (ji)   = MAX( - h_s_1d(ji) , - evap_ice_1d(ji) * r1_rhos  * rdt_ice ) 
     276            zevap_rema(ji)   = evap_ice_1d(ji) * rdt_ice + zdh_s_sub(ji) * rhos   ! remaining evap in kg.m-2 (used for ice melting later on) 
    277277            zdeltah   (ji,1) = MAX( zdh_s_sub(ji), - zdh_s_pre(ji) ) 
    278278             
     
    280280               &                 ( zdeltah(ji,1) * zqprec(ji) + ( zdh_s_sub(ji) - zdeltah(ji,1) ) * e_s_1d(ji,1) )  & 
    281281               &                 * a_i_1d(ji) * r1_Dt_ice 
    282             wfx_snw_sub_1d(ji) = wfx_snw_sub_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_sub(ji) * r1_Dt_ice   ! Mass flux by sublimation 
     282            wfx_snw_sub_1d(ji) = wfx_snw_sub_1d(ji) - rhos * a_i_1d(ji) * zdh_s_sub(ji) * r1_Dt_ice   ! Mass flux by sublimation 
    283283             
    284284            ! new snow thickness 
     
    309309            e_s_1d(ji,jk) = rswitch / MAX( h_s_1d(ji), epsi20 ) *            & 
    310310              &             ( ( zdh_s_pre(ji)              ) * zqprec(ji) +  & 
    311               &               ( h_s_1d(ji) - zdh_s_pre(ji) ) * rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) ) 
     311              &               ( h_s_1d(ji) - zdh_s_pre(ji) ) * rhos * ( rcpi * ( rt0 - t_s_1d(ji,jk) ) + rLfus ) ) 
    312312         END DO 
    313313      END DO 
     
    326326            IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN   !-- Internal melting 
    327327 
    328                zEi            = - e_i_1d(ji,jk) * r1_rhoic            ! Specific enthalpy of layer k [J/kg, <0]        
     328               zEi            = - e_i_1d(ji,jk) * r1_rhoi             ! Specific enthalpy of layer k [J/kg, <0]        
    329329               zdE            = 0._wp                                 ! Specific enthalpy difference (J/kg, <0) 
    330330                                                                      ! set up at 0 since no energy is needed to melt water...(it is already melted) 
    331331               zdeltah(ji,jk) = MIN( 0._wp , - zh_i(ji,jk) )          ! internal melting occurs when the internal temperature is above freezing      
    332332                                                                      ! this should normally not happen, but sometimes, heat diffusion leads to this 
    333                zfmdt          = - zdeltah(ji,jk) * rhoic              ! Mass flux x time step > 0 
     333               zfmdt          = - zdeltah(ji,jk) * rhoi               ! Mass flux x time step > 0 
    334334                          
    335335               dh_i_itm(ji)   = dh_i_itm(ji) + zdeltah(ji,jk)         ! Cumulate internal melting 
    336336                
    337                zfmdt          = - rhoic * zdeltah(ji,jk)              ! Recompute mass flux [kg/m2, >0] 
     337               zfmdt          = - rhoi * zdeltah(ji,jk)               ! Recompute mass flux [kg/m2, >0] 
    338338 
    339339               hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_Dt_ice                           ! Heat flux to the ocean [W.m-2], <0 
    340340               !                                                                                                  ice enthalpy zEi is "sent" to the ocean 
    341                sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_Dt_ice   ! Salt flux 
     341               sfx_res_1d(ji) = sfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_Dt_ice   ! Salt flux 
    342342               !                                                                                                  using s_i_1d and not sz_i_1d(jk) is ok 
    343                wfx_res_1d(ji) = wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice                ! Mass flux 
     343               wfx_res_1d(ji) = wfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice                ! Mass flux 
    344344 
    345345            ELSE                                        !-- Surface melting 
    346346                
    347                zEi            = - e_i_1d(ji,jk) * r1_rhoic            ! Specific enthalpy of layer k [J/kg, <0] 
     347               zEi            = - e_i_1d(ji,jk) * r1_rhoi             ! Specific enthalpy of layer k [J/kg, <0] 
    348348               zEw            =    rcp * ztmelts                      ! Specific enthalpy of resulting meltwater [J/kg, <0] 
    349349               zdE            =    zEi - zEw                          ! Specific enthalpy difference < 0 
     
    351351               zfmdt          = - zq_su(ji) / zdE                     ! Mass flux to the ocean [kg/m2, >0] 
    352352                
    353                zdeltah(ji,jk) = - zfmdt * r1_rhoic                    ! Melt of layer jk [m, <0] 
     353               zdeltah(ji,jk) = - zfmdt * r1_rhoi                     ! Melt of layer jk [m, <0] 
    354354                
    355355               zdeltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk) , - zh_i(ji,jk) ) )    ! Melt of layer jk cannot exceed the layer thickness [m, <0] 
    356356                
    357                zq_su(ji)      = MAX( 0._wp , zq_su(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat 
     357               zq_su(ji)      = MAX( 0._wp , zq_su(ji) - zdeltah(ji,jk) * rhoi * zdE ) ! update available heat 
    358358                
    359359               dh_i_sum(ji)   = dh_i_sum(ji) + zdeltah(ji,jk)         ! Cumulate surface melt 
    360360                
    361                zfmdt          = - rhoic * zdeltah(ji,jk)              ! Recompute mass flux [kg/m2, >0] 
     361               zfmdt          = - rhoi * zdeltah(ji,jk)               ! Recompute mass flux [kg/m2, >0] 
    362362                
    363363               zQm            = zfmdt * zEw                           ! Energy of the melt water sent to the ocean [J/m2, <0] 
    364364                
    365                sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_Dt_ice   ! Salt flux >0 
     365               sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_Dt_ice   ! Salt flux >0 
    366366               !                                                                                                  using s_i_1d and not sz_i_1d(jk) is ok) 
    367367               hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_Dt_ice                           ! Heat flux [W.m-2], < 0 
    368368               hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_Dt_ice                           ! Heat flux used in this process [W.m-2], > 0   
    369369               !  
    370                wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice                ! Mass flux 
     370               wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice                ! Mass flux 
    371371                
    372372            END IF 
     
    374374            ! Ice sublimation 
    375375            ! --------------- 
    376             zdum            = MAX( - ( zh_i(ji,jk) + zdeltah(ji,jk) ) , - zevap_rema(ji) * r1_rhoic ) 
     376            zdum            = MAX( - ( zh_i(ji,jk) + zdeltah(ji,jk) ) , - zevap_rema(ji) * r1_rhoi ) 
    377377            zdeltah (ji,jk) = zdeltah (ji,jk) + zdum 
    378378            dh_i_sub(ji)    = dh_i_sub(ji)    + zdum 
    379379             
    380             sfx_sub_1d(ji)     = sfx_sub_1d(ji) - rhoic * a_i_1d(ji) * zdum * s_i_1d(ji) * r1_Dt_ice ! Salt flux >0 
     380            sfx_sub_1d(ji)     = sfx_sub_1d(ji) - rhoi * a_i_1d(ji) * zdum * s_i_1d(ji) * r1_Dt_ice ! Salt flux >0 
    381381            !                                                                                          clem: flux is sent to the ocean for simplicity 
    382382            !                                                                                                but salt should remain in the ice except 
     
    384384            hfx_sub_1d(ji)     = hfx_sub_1d(ji) + zdum * e_i_1d(ji,jk) * a_i_1d(ji) * r1_Dt_ice      ! Heat flux [W.m-2], < 0 
    385385 
    386             wfx_ice_sub_1d(ji) = wfx_ice_sub_1d(ji) - rhoic * a_i_1d(ji) * zdum * r1_Dt_ice          ! Mass flux > 0 
     386            wfx_ice_sub_1d(ji) = wfx_ice_sub_1d(ji) - rhoi * a_i_1d(ji) * zdum * r1_Dt_ice          ! Mass flux > 0 
    387387 
    388388            ! update remaining mass flux 
    389             zevap_rema(ji)  = zevap_rema(ji) + zdum * rhoic 
     389            zevap_rema(ji)  = zevap_rema(ji) + zdum * rhoi 
    390390             
    391391            ! record which layers have disappeared (for bottom melting)  
     
    450450               zt_i_new      = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i) 
    451451                
    452                zEi           = cpic * ( zt_i_new - (ztmelts+rt0) ) &                                  ! Specific enthalpy of forming ice (J/kg, <0) 
    453                   &            - lfus * ( 1.0 - ztmelts / ( zt_i_new - rt0 ) ) + rcp * ztmelts 
     452               zEi           = rcpi * ( zt_i_new - (ztmelts+rt0) ) &                                  ! Specific enthalpy of forming ice (J/kg, <0) 
     453                  &          - rLfus * ( 1.0 - ztmelts / ( zt_i_new - rt0 ) ) + rcp * ztmelts 
    454454 
    455455               zEw           = rcp  * ( t_bo_1d(ji) - rt0 )                                           ! Specific enthalpy of seawater (J/kg, < 0) 
     
    457457               zdE           = zEi - zEw                                                              ! Specific enthalpy difference (J/kg, <0) 
    458458 
    459                dh_i_bog(ji)  = rdt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoic ) ) 
     459               dh_i_bog(ji)  = rdt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoi ) ) 
    460460                
    461461            END DO 
    462462            ! Contribution to Energy and Salt Fluxes                                     
    463             zfmdt          = - rhoic * dh_i_bog(ji)                                                   ! Mass flux x time step (kg/m2, < 0) 
     463            zfmdt          = - rhoi * dh_i_bog(ji)                                                   ! Mass flux x time step (kg/m2, < 0) 
    464464             
    465465            hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_Dt_ice                           ! Heat flux to the ocean [W.m-2], >0 
    466466            hfx_bog_1d(ji) = hfx_bog_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_Dt_ice                           ! Heat flux used in this process [W.m-2], <0 
    467467             
    468             sfx_bog_1d(ji) = sfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bog(ji) * s_i_new(ji) * r1_Dt_ice    ! Salt flux, <0 
    469  
    470             wfx_bog_1d(ji) = wfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bog(ji) * r1_Dt_ice                  ! Mass flux, <0 
     468            sfx_bog_1d(ji) = sfx_bog_1d(ji) - rhoi * a_i_1d(ji) * dh_i_bog(ji) * s_i_new(ji) * r1_Dt_ice    ! Salt flux, <0 
     469 
     470            wfx_bog_1d(ji) = wfx_bog_1d(ji) - rhoi * a_i_1d(ji) * dh_i_bog(ji) * r1_Dt_ice                  ! Mass flux, <0 
    471471 
    472472            ! update heat content (J.m-2) and layer thickness 
    473             eh_i_old(ji,nlay_i+1) = eh_i_old(ji,nlay_i+1) + dh_i_bog(ji) * (-zEi * rhoic) 
     473            eh_i_old(ji,nlay_i+1) = eh_i_old(ji,nlay_i+1) + dh_i_bog(ji) * (-zEi * rhoi) 
    474474            h_i_old (ji,nlay_i+1) = h_i_old (ji,nlay_i+1) + dh_i_bog(ji) 
    475475 
     
    489489               IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN   !-- Internal melting 
    490490 
    491                   zEi               = - e_i_1d(ji,jk) * r1_rhoic    ! Specific enthalpy of melting ice (J/kg, <0) 
     491                  zEi               = - e_i_1d(ji,jk) * r1_rhoi     ! Specific enthalpy of melting ice (J/kg, <0) 
    492492                  zdE               = 0._wp                         ! Specific enthalpy difference   (J/kg, <0) 
    493493                                                                    !    set up at 0 since no energy is needed to melt water...(it is already melted) 
     
    497497                  dh_i_itm (ji)     = dh_i_itm(ji) + zdeltah(ji,jk) 
    498498 
    499                   zfmdt             = - zdeltah(ji,jk) * rhoic      ! Mass flux x time step > 0 
     499                  zfmdt             = - zdeltah(ji,jk) * rhoi       ! Mass flux x time step > 0 
    500500 
    501501                  hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_Dt_ice                           ! Heat flux to the ocean [W.m-2], <0 
    502502                  !                                                                                                  ice enthalpy zEi is "sent" to the ocean 
    503                   sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_Dt_ice   ! Salt flux 
     503                  sfx_res_1d(ji) = sfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_Dt_ice   ! Salt flux 
    504504                  !                                                                                                  using s_i_1d and not sz_i_1d(jk) is ok 
    505                   wfx_res_1d(ji) = wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice                ! Mass flux 
     505                  wfx_res_1d(ji) = wfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice                ! Mass flux 
    506506 
    507507                  ! update heat content (J.m-2) and layer thickness 
     
    511511               ELSE                                        !-- Basal melting 
    512512 
    513                   zEi             = - e_i_1d(ji,jk) * r1_rhoic                                ! Specific enthalpy of melting ice (J/kg, <0) 
     513                  zEi             = - e_i_1d(ji,jk) * r1_rhoi                                 ! Specific enthalpy of melting ice (J/kg, <0) 
    514514                  zEw             = rcp * ztmelts                                             ! Specific enthalpy of meltwater (J/kg, <0) 
    515515                  zdE             = zEi - zEw                                                 ! Specific enthalpy difference   (J/kg, <0) 
     
    517517                  zfmdt           = - zq_bo(ji) / zdE                                         ! Mass flux x time step (kg/m2, >0) 
    518518 
    519                   zdeltah(ji,jk)  = - zfmdt * r1_rhoic                                        ! Gross thickness change 
     519                  zdeltah(ji,jk)  = - zfmdt * r1_rhoi                                         ! Gross thickness change 
    520520 
    521521                  zdeltah(ji,jk)  = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) )       ! bound thickness change 
    522522                   
    523                   zq_bo(ji)       = MAX( 0._wp , zq_bo(ji) - zdeltah(ji,jk) * rhoic * zdE )   ! update available heat. MAX is necessary for roundup errors 
    524  
    525                   dh_i_bom(ji)    = dh_i_bom(ji) + zdeltah(ji,jk)                            ! Update basal melt 
    526  
    527                   zfmdt           = - zdeltah(ji,jk) * rhoic                                  ! Mass flux x time step > 0 
     523                  zq_bo(ji)       = MAX( 0._wp , zq_bo(ji) - zdeltah(ji,jk) * rhoi * zdE )    ! update available heat. MAX is necessary for roundup errors 
     524 
     525                  dh_i_bom(ji)    = dh_i_bom(ji) + zdeltah(ji,jk)                             ! Update basal melt 
     526 
     527                  zfmdt           = - zdeltah(ji,jk) * rhoi                                   ! Mass flux x time step > 0 
    528528 
    529529                  zQm             = zfmdt * zEw                                               ! Heat exchanged with ocean 
     
    532532                  hfx_bom_1d(ji)  = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_Dt_ice                           ! Heat used in this process [W.m-2], >0   
    533533 
    534                   sfx_bom_1d(ji)  = sfx_bom_1d(ji) - rhoic *  a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_Dt_ice  ! Salt flux 
     534                  sfx_bom_1d(ji)  = sfx_bom_1d(ji) - rhoi *  a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_Dt_ice  ! Salt flux 
    535535                  !                                                                                                   using s_i_1d and not sz_i_1d(jk) is ok 
    536536                   
    537                   wfx_bom_1d(ji)  = wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice                ! Mass flux 
     537                  wfx_bom_1d(ji)  = wfx_bom_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice                ! Mass flux 
    538538 
    539539                  ! update heat content (J.m-2) and layer thickness 
     
    566566         zq_rema(ji)        = zq_rema(ji) + zdeltah(ji,1) * e_s_1d(ji,1)                               ! update available heat (J.m-2) 
    567567         hfx_snw_1d(ji)     = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * e_s_1d(ji,1) * r1_Dt_ice   ! Heat used to melt snow, W.m-2 (>0) 
    568          wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_Dt_ice      ! Mass flux 
     568         wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhos * a_i_1d(ji) * zdeltah(ji,1) * r1_Dt_ice      ! Mass flux 
    569569         dh_s_mlt(ji)       = dh_s_mlt(ji) + zdeltah(ji,1) 
    570570         !     
     
    580580      ! When snow load excesses Archimede's limit, snow-ice interface goes down under sea-level,  
    581581      ! flooding of seawater transforms snow into ice dh_snowice is positive for the ice 
    582       z1_rho = 1._wp / ( rhosn+rho0-rhoic ) 
     582      z1_rho = 1._wp / ( rhos + rho0 - rhoi ) 
    583583      DO ji = 1, npti 
    584584         ! 
    585          dh_snowice(ji) = MAX(  0._wp , ( rhosn * h_s_1d(ji) + (rhoic-rho0) * h_i_1d(ji) ) * z1_rho ) 
     585         dh_snowice(ji) = MAX(  0._wp , ( rhos * h_s_1d(ji) + (rhoi-rho0) * h_i_1d(ji) ) * z1_rho ) 
    586586 
    587587         h_i_1d(ji)    = h_i_1d(ji) + dh_snowice(ji) 
     
    589589 
    590590         ! Contribution to energy flux to the ocean [J/m2], >0 (if sst<0) 
    591          zfmdt          = ( rhosn - rhoic ) * dh_snowice(ji)    ! <0 
     591         zfmdt          = ( rhos - rhoi ) * dh_snowice(ji)    ! <0 
    592592         zEw            = rcp * sst_1d(ji) 
    593593         zQm            = zfmdt * zEw  
     
    599599         ! Case constant salinity in time: virtual salt flux to keep salinity constant 
    600600         IF( nn_icesal /= 2 )  THEN 
    601             sfx_bri_1d(ji) = sfx_bri_1d(ji) - sss_1d (ji) * a_i_1d(ji) * zfmdt                  * r1_Dt_ice  & ! put back sss_m     into the ocean 
    602                &                            - s_i_1d(ji)  * a_i_1d(ji) * dh_snowice(ji) * rhoic * r1_Dt_ice    ! and get  rn_icesal from the ocean  
     601            sfx_bri_1d(ji) = sfx_bri_1d(ji) - sss_1d (ji) * a_i_1d(ji) * zfmdt                 * r1_Dt_ice  & ! put back sss_m     into the ocean 
     602               &                            - s_i_1d(ji)  * a_i_1d(ji) * dh_snowice(ji) * rhoi * r1_Dt_ice    ! and get  rn_icesal from the ocean  
    603603         ENDIF 
    604604 
    605605         ! Mass flux: All snow is thrown in the ocean, and seawater is taken to replace the volume 
    606          wfx_sni_1d(ji)     = wfx_sni_1d(ji)     - a_i_1d(ji) * dh_snowice(ji) * rhoic * r1_Dt_ice 
    607          wfx_snw_sni_1d(ji) = wfx_snw_sni_1d(ji) + a_i_1d(ji) * dh_snowice(ji) * rhosn * r1_Dt_ice 
     606         wfx_sni_1d(ji)     = wfx_sni_1d(ji)     - a_i_1d(ji) * dh_snowice(ji) * rhoi * r1_Dt_ice 
     607         wfx_snw_sni_1d(ji) = wfx_snw_sni_1d(ji) + a_i_1d(ji) * dh_snowice(ji) * rhos * r1_Dt_ice 
    608608 
    609609         ! update heat content (J.m-2) and layer thickness 
     
    627627            e_s_1d(ji,jk) = rswitch * e_s_1d(ji,jk) 
    628628            ! recalculate t_s_1d from e_s_1d 
    629             t_s_1d(ji,jk) = rt0 + rswitch * ( - e_s_1d(ji,jk) * r1_rhosn * r1_cpic + lfus * r1_cpic ) 
     629            t_s_1d(ji,jk) = rt0 + rswitch * ( - e_s_1d(ji,jk) * r1_rhos * r1_cpi + rLfus * r1_cpi ) 
    630630         END DO 
    631631      END DO 
Note: See TracChangeset for help on using the changeset viewer.