Changeset 9935


Ignore:
Timestamp:
2018-07-12T16:12:48+02:00 (2 years ago)
Author:
clem
Message:

changes to facilitate the future merge with Gurvan's crazy dev branch

Location:
NEMO/trunk/src
Files:
24 edited

Legend:

Unmodified
Added
Removed
  • NEMO/trunk/src/ICE/icecor.F90

    r9604 r9935  
    8686      IF ( nn_icesal == 2 ) THEN    !  salinity must stay in bounds [Simin,Simax]        ! 
    8787      !                             !----------------------------------------------------- 
    88          zzc = rhoic * r1_rdtice 
     88         zzc = rhoi * r1_rdtice 
    8989         DO jl = 1, jpl 
    9090            DO jj = 1, jpj  
     
    139139                  &                  + SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) )  ) * r1_rdtice 
    140140               !                 ! salt, volume 
    141                diag_sice(ji,jj) = SUM( sv_i(ji,jj,:) - sv_i_b(ji,jj,:) ) * rhoic * r1_rdtice 
    142                diag_vice(ji,jj) = SUM( v_i (ji,jj,:) - v_i_b (ji,jj,:) ) * rhoic * r1_rdtice 
    143                diag_vsnw(ji,jj) = SUM( v_s (ji,jj,:) - v_s_b (ji,jj,:) ) * rhosn * r1_rdtice 
     141               diag_sice(ji,jj) = SUM( sv_i(ji,jj,:) - sv_i_b(ji,jj,:) ) * rhoi * r1_rdtice 
     142               diag_vice(ji,jj) = SUM( v_i (ji,jj,:) - v_i_b (ji,jj,:) ) * rhoi * r1_rdtice 
     143               diag_vsnw(ji,jj) = SUM( v_s (ji,jj,:) - v_s_b (ji,jj,:) ) * rhos * r1_rdtice 
    144144            END DO 
    145145         END DO 
     
    160160                  &                                   + SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) )  ) * r1_rdtice 
    161161               !                 ! salt, volume 
    162                diag_sice(ji,jj) = diag_sice(ji,jj) + SUM( sv_i(ji,jj,:) - sv_i_b(ji,jj,:) ) * rhoic * r1_rdtice 
    163                diag_vice(ji,jj) = diag_vice(ji,jj) + SUM( v_i (ji,jj,:) - v_i_b (ji,jj,:) ) * rhoic * r1_rdtice 
    164                diag_vsnw(ji,jj) = diag_vsnw(ji,jj) + SUM( v_s (ji,jj,:) - v_s_b (ji,jj,:) ) * rhosn * r1_rdtice 
     162               diag_sice(ji,jj) = diag_sice(ji,jj) + SUM( sv_i(ji,jj,:) - sv_i_b(ji,jj,:) ) * rhoi * r1_rdtice 
     163               diag_vice(ji,jj) = diag_vice(ji,jj) + SUM( v_i (ji,jj,:) - v_i_b (ji,jj,:) ) * rhoi * r1_rdtice 
     164               diag_vsnw(ji,jj) = diag_vsnw(ji,jj) + SUM( v_s (ji,jj,:) - v_s_b (ji,jj,:) ) * rhos * r1_rdtice 
    165165            END DO 
    166166         END DO 
  • NEMO/trunk/src/ICE/icectl.F90

    r9913 r9935  
    9393            &                  ) *  e1e2t(:,:) ) * zconv 
    9494 
    95          pdiag_v = glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) * e1e2t * zconv ) 
    96  
    97          pdiag_s = glob_sum( SUM( sv_i * rhoic             , dim=3 ) * e1e2t * zconv ) 
     95         pdiag_v = glob_sum( SUM( v_i * rhoi + v_s * rhos, dim=3 ) * e1e2t * zconv ) 
     96 
     97         pdiag_s = glob_sum( SUM( sv_i * rhoi            , dim=3 ) * e1e2t * zconv ) 
    9898 
    9999         pdiag_t = glob_sum( (  SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 )     & 
     
    120120  
    121121         ! outputs 
    122          zv = ( ( glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) * e1e2t ) * zconv  & 
     122         zv = ( ( glob_sum( SUM( v_i * rhoi + v_s * rhos, dim=3 ) * e1e2t ) * zconv  & 
    123123            &     - pdiag_v ) * r1_rdtice - zfv ) * rday 
    124124 
    125          zs = ( ( glob_sum( SUM( sv_i * rhoic             , dim=3 ) * e1e2t ) * zconv  & 
     125         zs = ( ( glob_sum( SUM( sv_i * rhoi            , dim=3 ) * e1e2t ) * zconv  & 
    126126            &     - pdiag_s ) * r1_rdtice + zfs ) * rday 
    127127 
     
    131131 
    132132         ! zvtrp and zetrp must be close to 0 if the advection scheme is conservative 
    133          zvtrp = glob_sum( ( diag_trp_vi * rhoic + diag_trp_vs * rhosn ) * e1e2t  ) * zconv * rday  
    134          zetrp = glob_sum( ( diag_trp_ei         + diag_trp_es         ) * e1e2t  ) * zconv 
     133         zvtrp = glob_sum( ( diag_trp_vi * rhoi + diag_trp_vs * rhos ) * e1e2t  ) * zconv * rday  
     134         zetrp = glob_sum( ( diag_trp_ei        + diag_trp_es        ) * e1e2t  ) * zconv 
    135135 
    136136         zvmin = glob_min( v_i ) 
     
    382382            DO jj = 1, jpj 
    383383               DO ji = 1, jpi 
    384                   ztmelts    =  -tmut * sz_i(ji,jj,jk,jl) + rt0 
     384                  ztmelts    =  -rTmlt * sz_i(ji,jj,jk,jl) + rt0 
    385385                  IF( t_i(ji,jj,jk,jl) >= ztmelts  .AND.  v_i(ji,jj,jl) > 1.e-10   & 
    386386                     &                             .AND.  a_i(ji,jj,jl) > 0._wp   ) THEN 
  • NEMO/trunk/src/ICE/icedia.F90

    r9912 r9935  
    110110      ! 3 -  Content variations ! 
    111111      ! ----------------------- ! 
    112       zdiff_vol = r1_rau0 * glob_sum( ( rhoic*vt_i(:,:) + rhosn*vt_s(:,:) - vol_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-9   ! freshwater trend (km3)  
    113       zdiff_sal = r1_rau0 * glob_sum( ( rhoic* SUM( sv_i(:,:,:), dim=3 ) - sal_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-9   ! salt content trend (km3*pss) 
    114       zdiff_tem =           glob_sum( ( et_i(:,:) + et_s(:,:)             - tem_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-20  ! heat content trend (1.e20 J) 
     112      zdiff_vol = r1_rau0 * glob_sum( ( rhoi*vt_i(:,:) + rhos*vt_s(:,:) - vol_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-9   ! freshwater trend (km3)  
     113      zdiff_sal = r1_rau0 * glob_sum( ( rhoi* SUM( sv_i(:,:,:), dim=3 ) - sal_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-9   ! salt content trend (km3*pss) 
     114      zdiff_tem =           glob_sum( ( et_i(:,:) + et_s(:,:)           - tem_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-20  ! heat content trend (1.e20 J) 
    115115      !                               + SUM( qevap_ice * a_i_b, dim=3 )       !! clem: I think this term should not be there (but needs a check) 
    116116 
     
    246246            frc_sal     = 0._wp                                                  
    247247            ! record initial ice volume, salt and temp 
    248             vol_loc_ini(:,:) = rhoic * vt_i(:,:) + rhosn * vt_s(:,:)  ! ice/snow volume (kg/m2) 
    249             tem_loc_ini(:,:) = et_i(:,:) + et_s(:,:)                  ! ice/snow heat content (J) 
    250             sal_loc_ini(:,:) = rhoic * SUM( sv_i(:,:,:), dim=3 )      ! ice salt content (pss*kg/m2) 
     248            vol_loc_ini(:,:) = rhoi * vt_i(:,:) + rhos * vt_s(:,:)  ! ice/snow volume (kg/m2) 
     249            tem_loc_ini(:,:) = et_i(:,:) + et_s(:,:)                ! ice/snow heat content (J) 
     250            sal_loc_ini(:,:) = rhoi * SUM( sv_i(:,:,:), dim=3 )     ! ice salt content (pss*kg/m2) 
    251251         ENDIF 
    252252         ! 
  • NEMO/trunk/src/ICE/icedyn_adv.F90

    r9888 r9935  
    119119      diag_trp_vi(:,:) = SUM(     v_i (:,:,:)          - v_i_b (:,:,:)                  , dim=3 ) * r1_rdtice 
    120120      diag_trp_vs(:,:) = SUM(     v_s (:,:,:)          - v_s_b (:,:,:)                  , dim=3 ) * r1_rdtice 
    121       IF( iom_use('icemtrp') )   CALL iom_put( "icemtrp" , diag_trp_vi * rhoic          )   ! ice mass transport 
    122       IF( iom_use('snwmtrp') )   CALL iom_put( "snwmtrp" , diag_trp_vs * rhosn          )   ! snw mass transport 
    123       IF( iom_use('salmtrp') )   CALL iom_put( "salmtrp" , diag_trp_sv * rhoic * 1.e-03 )   ! salt mass transport (kg/m2/s) 
    124       IF( iom_use('dihctrp') )   CALL iom_put( "dihctrp" , -diag_trp_ei                 )   ! advected ice heat content (W/m2) 
    125       IF( iom_use('dshctrp') )   CALL iom_put( "dshctrp" , -diag_trp_es                 )   ! advected snw heat content (W/m2) 
     121      IF( iom_use('icemtrp') )   CALL iom_put( "icemtrp" , diag_trp_vi * rhoi          )   ! ice mass transport 
     122      IF( iom_use('snwmtrp') )   CALL iom_put( "snwmtrp" , diag_trp_vs * rhos          )   ! snw mass transport 
     123      IF( iom_use('salmtrp') )   CALL iom_put( "salmtrp" , diag_trp_sv * rhoi * 1.e-03 )   ! salt mass transport (kg/m2/s) 
     124      IF( iom_use('dihctrp') )   CALL iom_put( "dihctrp" , -diag_trp_ei                )   ! advected ice heat content (W/m2) 
     125      IF( iom_use('dshctrp') )   CALL iom_put( "dshctrp" , -diag_trp_es                )   ! advected snw heat content (W/m2) 
    126126 
    127127      ! controls 
  • NEMO/trunk/src/ICE/icedyn_rdgrft.F90

    r9880 r9935  
    547547               ! volume and enthalpy (J/m2, >0) of seawater trapped into ridges 
    548548               vsw = v_i_2d(ji,jl1) * afrdg * rn_porordg 
    549                ersw(ji) = -rhoic * vsw * rcp * sst_1d(ji)   ! clem: if sst>0, then ersw <0 (is that possible?) 
     549               ersw(ji) = -rhoi * vsw * rcp * sst_1d(ji)   ! clem: if sst>0, then ersw <0 (is that possible?) 
    550550 
    551551               ! volume etc of ridging / rafting ice and new ridges (vi, vs, sm, oi, es, ei) 
     
    574574 
    575575               ! Ice-ocean exchanges associated with ice porosity 
    576                wfx_dyn_1d(ji) = wfx_dyn_1d(ji) - vsw * rhoic * r1_rdtice   ! increase in ice volume due to seawater frozen in voids 
    577                sfx_dyn_1d(ji) = sfx_dyn_1d(ji) - vsw * sss_1d(ji) * rhoic * r1_rdtice 
     576               wfx_dyn_1d(ji) = wfx_dyn_1d(ji) - vsw * rhoi * r1_rdtice   ! increase in ice volume due to seawater frozen in voids 
     577               sfx_dyn_1d(ji) = sfx_dyn_1d(ji) - vsw * sss_1d(ji) * rhoi * r1_rdtice 
    578578               hfx_dyn_1d(ji) = hfx_dyn_1d(ji) + ersw(ji) * r1_rdtice          ! > 0 [W.m-2]  
    579579 
    580580               ! Put the snow lost by ridging into the ocean 
    581581               !  Note that esrdg > 0; the ocean must cool to melt snow. If the ocean temp = Tf already, new ice must grow. 
    582                wfx_snw_dyn_1d(ji) = wfx_snw_dyn_1d(ji) + ( rhosn * vsrdg(ji) * ( 1._wp - rn_fsnwrdg )   &   ! fresh water source for ocean 
    583                   &                                      + rhosn * vsrft(ji) * ( 1._wp - rn_fsnwrft ) ) * r1_rdtice 
     582               wfx_snw_dyn_1d(ji) = wfx_snw_dyn_1d(ji) + ( rhos * vsrdg(ji) * ( 1._wp - rn_fsnwrdg )   &   ! fresh water source for ocean 
     583                  &                                      + rhos * vsrft(ji) * ( 1._wp - rn_fsnwrft ) ) * r1_rdtice 
    584584 
    585585               ! Put the melt pond water into the ocean 
     
    587587               !       is no net mass flux between melt ponds and the ocean (see icethd_pnd.F90 for ex.) 
    588588               !IF ( ln_pnd_fwb ) THEN 
    589                !   wfx_pnd_1d(ji) = wfx_pnd_1d(ji) + ( rhofw * vprdg(ji) * ( 1._wp - rn_fpndrdg )   &        ! fresh water source for ocean 
    590                !      &                              + rhofw * vprft(ji) * ( 1._wp - rn_fpndrft ) ) * r1_rdtice 
     589               !   wfx_pnd_1d(ji) = wfx_pnd_1d(ji) + ( rhow * vprdg(ji) * ( 1._wp - rn_fpndrdg )   &        ! fresh water source for ocean 
     590               !      &                              + rhow * vprft(ji) * ( 1._wp - rn_fpndrft ) ) * r1_rdtice 
    591591               !ENDIF 
    592592 
     
    594594               IF( nn_icesal /= 2 )  THEN 
    595595                  sirdg2(ji)     = sirdg2(ji)     - vsw * ( sss_1d(ji) - s_i_1d(ji) )        ! ridge salinity = s_i 
    596                   sfx_bri_1d(ji) = sfx_bri_1d(ji) + sss_1d(ji) * vsw * rhoic * r1_rdtice  &  ! put back sss_m into the ocean 
    597                      &                            - s_i_1d(ji) * vsw * rhoic * r1_rdtice     ! and get  s_i  from the ocean  
     596                  sfx_bri_1d(ji) = sfx_bri_1d(ji) + sss_1d(ji) * vsw * rhoi * r1_rdtice  &  ! put back sss_m into the ocean 
     597                     &                            - s_i_1d(ji) * vsw * rhoi * r1_rdtice     ! and get  s_i  from the ocean  
    598598               ENDIF 
    599599 
  • NEMO/trunk/src/ICE/icedyn_rhg_evp.F90

    r9866 r9935  
    285285 
    286286            ! Ice/snow mass at U-V points 
    287             zm1 = ( rhosn * vt_s(ji  ,jj  ) + rhoic * vt_i(ji  ,jj  ) ) 
    288             zm2 = ( rhosn * vt_s(ji+1,jj  ) + rhoic * vt_i(ji+1,jj  ) ) 
    289             zm3 = ( rhosn * vt_s(ji  ,jj+1) + rhoic * vt_i(ji  ,jj+1) ) 
     287            zm1 = ( rhos * vt_s(ji  ,jj  ) + rhoi * vt_i(ji  ,jj  ) ) 
     288            zm2 = ( rhos * vt_s(ji+1,jj  ) + rhoi * vt_i(ji+1,jj  ) ) 
     289            zm3 = ( rhos * vt_s(ji  ,jj+1) + rhoi * vt_i(ji  ,jj+1) ) 
    290290            zmassU = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1) 
    291291            zmassV = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1) 
     
    795795               zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * rswitch 
    796796                
    797                zdiag_xmtrp_ice(ji,jj) = rhoic * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component 
    798                zdiag_ymtrp_ice(ji,jj) = rhoic * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) !        ''           Y-   '' 
    799                 
    800                zdiag_xmtrp_snw(ji,jj) = rhosn * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component 
    801                zdiag_ymtrp_snw(ji,jj) = rhosn * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) !          ''          Y-   '' 
    802                 
    803                zdiag_xatrp(ji,jj)     = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) )         ! area transport,      X-component 
    804                zdiag_yatrp(ji,jj)     = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) )         !        ''            Y-   '' 
     797               zdiag_xmtrp_ice(ji,jj) = rhoi * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component 
     798               zdiag_ymtrp_ice(ji,jj) = rhoi * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) !        ''           Y-   '' 
     799                
     800               zdiag_xmtrp_snw(ji,jj) = rhos * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component 
     801               zdiag_ymtrp_snw(ji,jj) = rhos * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) !          ''          Y-   '' 
     802                
     803               zdiag_xatrp(ji,jj)     = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) )        ! area transport,      X-component 
     804               zdiag_yatrp(ji,jj)     = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) )        !        ''            Y-   '' 
    805805                
    806806            END DO 
  • NEMO/trunk/src/ICE/iceistate.F90

    r9656 r9935  
    295295                  ! In case snow load is in excess that would lead to transformation from snow to ice 
    296296                  ! Then, transfer the snow excess into the ice (different from icethd_dh) 
    297                   zdh = MAX( 0._wp, ( rhosn * h_s(ji,jj,jl) + ( rhoic - rau0 ) * h_i(ji,jj,jl) ) * r1_rau0 )  
     297                  zdh = MAX( 0._wp, ( rhos * h_s(ji,jj,jl) + ( rhoi - rau0 ) * h_i(ji,jj,jl) ) * r1_rau0 )  
    298298                  ! recompute h_i, h_s avoiding out of bounds values 
    299299                  h_i(ji,jj,jl) = MIN( hi_max(jl), h_i(ji,jj,jl) + zdh ) 
    300                   h_s(ji,jj,jl) = MAX( 0._wp, h_s(ji,jj,jl) - zdh * rhoic * r1_rhosn ) 
     300                  h_s(ji,jj,jl) = MAX( 0._wp, h_s(ji,jj,jl) - zdh * rhoi * r1_rhos ) 
    301301                  ! 
    302302                  ! ice volume, salt content, age content 
     
    321321                     t_s(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0 
    322322                     ! Snow energy of melting 
    323                      e_s(ji,jj,jk,jl) = zswitch(ji,jj) * rhosn * ( cpic * ( rt0 - t_s(ji,jj,jk,jl) ) + lfus ) 
     323                     e_s(ji,jj,jk,jl) = zswitch(ji,jj) * rhos * ( rcpi * ( rt0 - t_s(ji,jj,jk,jl) ) + rLfus ) 
    324324                     ! 
    325325                     ! Mutliply by volume, and divide by number of layers to get heat content in J/m2 
     
    337337                     t_i (ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0  
    338338                     sz_i(ji,jj,jk,jl) = zswitch(ji,jj) * zsm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rn_simin 
    339                      ztmelts          = - tmut * sz_i(ji,jj,jk,jl) + rt0 !Melting temperature in K 
     339                     ztmelts          = - rTmlt * sz_i(ji,jj,jk,jl) + rt0 !Melting temperature in K 
    340340                     ! 
    341341                     ! heat content per unit volume 
    342                      e_i(ji,jj,jk,jl) = zswitch(ji,jj) * rhoic * (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) )         & 
    343                         &             + lfus * ( 1._wp - (ztmelts-rt0) / MIN( (t_i(ji,jj,jk,jl)-rt0) , -epsi20 )  )   & 
     342                     e_i(ji,jj,jk,jl) = zswitch(ji,jj) * rhoi * (   rcpi    * ( ztmelts - t_i(ji,jj,jk,jl) )           & 
     343                        &             + rLfus * ( 1._wp - (ztmelts-rt0) / MIN( (t_i(ji,jj,jk,jl)-rt0) , -epsi20 )  )   & 
    344344                        &             - rcp  * ( ztmelts - rt0 ) ) 
    345345                     ! 
     
    410410      ! 5) Snow-ice mass (case ice is fully embedded) 
    411411      !---------------------------------------------- 
    412       snwice_mass  (:,:) = tmask(:,:,1) * SUM( rhosn * v_s(:,:,:) + rhoic * v_i(:,:,:), dim=3  )   ! snow+ice mass 
     412      snwice_mass  (:,:) = tmask(:,:,1) * SUM( rhos * v_s(:,:,:) + rhoi * v_i(:,:,:), dim=3  )   ! snow+ice mass 
    413413      snwice_mass_b(:,:) = snwice_mass(:,:) 
    414414      ! 
  • NEMO/trunk/src/ICE/icethd.F90

    r9922 r9935  
    288288      DO jk = 1, nlay_i 
    289289         DO ji = 1, npti 
    290             ztmelts       = -tmut * sz_i_1d(ji,jk) 
     290            ztmelts       = -rTmlt * sz_i_1d(ji,jk) 
    291291            ! Conversion q(S,T) -> T (second order equation) 
    292             zbbb          = ( rcp - cpic ) * ztmelts + e_i_1d(ji,jk) * r1_rhoic - lfus 
    293             zccc          = SQRT( MAX( zbbb * zbbb - 4._wp * cpic * lfus * ztmelts, 0._wp ) ) 
    294             t_i_1d(ji,jk) = rt0 - ( zbbb + zccc ) * 0.5_wp * r1_cpic 
     292            zbbb          = ( rcp - rcpi ) * ztmelts + e_i_1d(ji,jk) * r1_rhoi - rLfus 
     293            zccc          = SQRT( MAX( zbbb * zbbb - 4._wp * rcpi * rLfus * ztmelts, 0._wp ) ) 
     294            t_i_1d(ji,jk) = rt0 - ( zbbb + zccc ) * 0.5_wp * r1_rcpi 
    295295             
    296296            ! mask temperature 
  • NEMO/trunk/src/ICE/icethd_da.F90

    r9604 r9935  
    137137             
    138138            ! Contribution to salt flux 
    139             sfx_lam_1d(ji) = sfx_lam_1d(ji) + rhoic *  h_i_1d(ji) * zda * s_i_1d(ji) * r1_rdtice 
     139            sfx_lam_1d(ji) = sfx_lam_1d(ji) + rhoi *  h_i_1d(ji) * zda * s_i_1d(ji) * r1_rdtice 
    140140             
    141141            ! Contribution to heat flux into the ocean [W.m-2], (<0)   
     
    144144             
    145145            ! Contribution to mass flux 
    146             wfx_lam_1d(ji) =  wfx_lam_1d(ji) + zda * r1_rdtice * ( rhoic * h_i_1d(ji) + rhosn * h_s_1d(ji) ) 
     146            wfx_lam_1d(ji) =  wfx_lam_1d(ji) + zda * r1_rdtice * ( rhoi * h_i_1d(ji) + rhos * h_s_1d(ji) ) 
    147147             
    148148            ! new concentration 
  • NEMO/trunk/src/ICE/icethd_dh.F90

    r9922 r9935  
    7676      REAL(wp) ::   zgrr         ! bottom growth rate 
    7777      REAL(wp) ::   zt_i_new     ! bottom formation temperature 
    78       REAL(wp) ::   z1_rho       ! 1/(rhosn+rau0-rhoic) 
     78      REAL(wp) ::   z1_rho       ! 1/(rhos+rau0-rhoi) 
    7979 
    8080      REAL(wp) ::   zQm          ! enthalpy exchanged with the ocean (J/m2), >0 towards the ocean 
     
    174174            IF( t_s_1d(ji,jk) > rt0 ) THEN 
    175175               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 
    176                wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + rhosn         * zh_s(ji,jk) * a_i_1d(ji) * r1_rdtice   ! mass flux 
     176               wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + rhos          * zh_s(ji,jk) * a_i_1d(ji) * r1_rdtice   ! mass flux 
    177177               ! updates 
    178178               dh_s_mlt(ji)    = dh_s_mlt(ji) - zh_s(ji,jk) 
     
    194194            ! 
    195195            ! --- precipitation --- 
    196             zdh_s_pre (ji) = zsnw(ji) * sprecip_1d(ji) * rdt_ice * r1_rhosn / at_i_1d(ji)   ! thickness change 
     196            zdh_s_pre (ji) = zsnw(ji) * sprecip_1d(ji) * rdt_ice * r1_rhos / at_i_1d(ji)   ! thickness change 
    197197            zqprec    (ji) = - qprec_ice_1d(ji)                                             ! enthalpy of the precip (>0, J.m-3) 
    198198            ! 
    199199            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) 
    200             wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhosn         * a_i_1d(ji) * zdh_s_pre(ji) * r1_rdtice   ! mass flux, <0 
     200            wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhos          * a_i_1d(ji) * zdh_s_pre(ji) * r1_rdtice   ! mass flux, <0 
    201201             
    202202            ! --- melt of falling snow --- 
     
    205205            zdeltah       (ji,1) = MAX( - zdh_s_pre(ji), zdeltah(ji,1) )                 ! bound melting  
    206206            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) 
    207             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 
     207            wfx_snw_sum_1d(ji)   = wfx_snw_sum_1d(ji) - rhos          * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice   ! snow melting only = water into the ocean (then without snow precip), >0 
    208208             
    209209            ! updates available heat + precipitations after melting 
     
    245245                
    246246               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) 
    247                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) 
     247               wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhos           * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice   ! snow melting only = water into the ocean (then without snow precip) 
    248248                
    249249               ! updates available heat + thickness 
     
    265265         IF( evap_ice_1d(ji) > 0._wp ) THEN 
    266266            ! 
    267             zdh_s_sub (ji)   = MAX( - h_s_1d(ji) , - evap_ice_1d(ji) * r1_rhosn * rdt_ice ) 
    268             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) 
     267            zdh_s_sub (ji)   = MAX( - h_s_1d(ji) , - evap_ice_1d(ji) * r1_rhos * rdt_ice ) 
     268            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) 
    269269            zdeltah   (ji,1) = MAX( zdh_s_sub(ji), - zdh_s_pre(ji) ) 
    270270             
     
    272272               &                 ( zdeltah(ji,1) * zqprec(ji) + ( zdh_s_sub(ji) - zdeltah(ji,1) ) * e_s_1d(ji,1) )  & 
    273273               &                 * a_i_1d(ji) * r1_rdtice 
    274             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 
     274            wfx_snw_sub_1d(ji) = wfx_snw_sub_1d(ji) - rhos * a_i_1d(ji) * zdh_s_sub(ji) * r1_rdtice   ! Mass flux by sublimation 
    275275             
    276276            ! new snow thickness 
     
    301301            e_s_1d(ji,jk) = rswitch / MAX( h_s_1d(ji), epsi20 ) *            & 
    302302              &             ( ( zdh_s_pre(ji)              ) * zqprec(ji) +  & 
    303               &               ( h_s_1d(ji) - zdh_s_pre(ji) ) * rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) ) 
     303              &               ( h_s_1d(ji) - zdh_s_pre(ji) ) * rhos * ( rcpi * ( rt0 - t_s_1d(ji,jk) ) + rLfus ) ) 
    304304         END DO 
    305305      END DO 
     
    314314      DO jk = 1, nlay_i 
    315315         DO ji = 1, npti 
    316             ztmelts = - tmut * sz_i_1d(ji,jk)   ! Melting point of layer k [C] 
     316            ztmelts = - rTmlt * sz_i_1d(ji,jk)   ! Melting point of layer k [C] 
    317317             
    318318            IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN   !-- Internal melting 
    319319 
    320                zEi            = - e_i_1d(ji,jk) * r1_rhoic            ! Specific enthalpy of layer k [J/kg, <0]        
     320               zEi            = - e_i_1d(ji,jk) * r1_rhoi             ! Specific enthalpy of layer k [J/kg, <0]        
    321321               zdE            = 0._wp                                 ! Specific enthalpy difference (J/kg, <0) 
    322322                                                                      ! set up at 0 since no energy is needed to melt water...(it is already melted) 
    323323               zdeltah(ji,jk) = MIN( 0._wp , - zh_i(ji,jk) )          ! internal melting occurs when the internal temperature is above freezing      
    324324                                                                      ! this should normally not happen, but sometimes, heat diffusion leads to this 
    325                zfmdt          = - zdeltah(ji,jk) * rhoic              ! Mass flux x time step > 0 
     325               zfmdt          = - zdeltah(ji,jk) * rhoi               ! Mass flux x time step > 0 
    326326                          
    327327               dh_i_itm(ji)   = dh_i_itm(ji) + zdeltah(ji,jk)         ! Cumulate internal melting 
    328328                
    329                zfmdt          = - rhoic * zdeltah(ji,jk)              ! Recompute mass flux [kg/m2, >0] 
     329               zfmdt          = - rhoi * zdeltah(ji,jk)               ! Recompute mass flux [kg/m2, >0] 
    330330 
    331331               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 
    332332               !                                                                                                  ice enthalpy zEi is "sent" to the ocean 
    333                sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice   ! Salt flux 
     333               sfx_res_1d(ji) = sfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice    ! Salt flux 
    334334               !                                                                                                  using s_i_1d and not sz_i_1d(jk) is ok 
    335                wfx_res_1d(ji) = wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice                ! Mass flux 
     335               wfx_res_1d(ji) = wfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice                 ! Mass flux 
    336336 
    337337            ELSE                                        !-- Surface melting 
    338338                
    339                zEi            = - e_i_1d(ji,jk) * r1_rhoic            ! Specific enthalpy of layer k [J/kg, <0] 
     339               zEi            = - e_i_1d(ji,jk) * r1_rhoi             ! Specific enthalpy of layer k [J/kg, <0] 
    340340               zEw            =    rcp * ztmelts                      ! Specific enthalpy of resulting meltwater [J/kg, <0] 
    341341               zdE            =    zEi - zEw                          ! Specific enthalpy difference < 0 
     
    343343               zfmdt          = - zq_top(ji) / zdE                    ! Mass flux to the ocean [kg/m2, >0] 
    344344                
    345                zdeltah(ji,jk) = - zfmdt * r1_rhoic                    ! Melt of layer jk [m, <0] 
     345               zdeltah(ji,jk) = - zfmdt * r1_rhoi                     ! Melt of layer jk [m, <0] 
    346346                
    347347               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] 
    348348                
    349                zq_top(ji)      = MAX( 0._wp , zq_top(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat 
     349               zq_top(ji)      = MAX( 0._wp , zq_top(ji) - zdeltah(ji,jk) * rhoi * zdE ) ! update available heat 
    350350                
    351351               dh_i_sum(ji)   = dh_i_sum(ji) + zdeltah(ji,jk)         ! Cumulate surface melt 
    352352                
    353                zfmdt          = - rhoic * zdeltah(ji,jk)              ! Recompute mass flux [kg/m2, >0] 
     353               zfmdt          = - rhoi * zdeltah(ji,jk)               ! Recompute mass flux [kg/m2, >0] 
    354354                
    355355               zQm            = zfmdt * zEw                           ! Energy of the melt water sent to the ocean [J/m2, <0] 
    356356                
    357                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 
     357               sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice    ! Salt flux >0 
    358358               !                                                                                                  using s_i_1d and not sz_i_1d(jk) is ok) 
    359359               hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice                           ! Heat flux [W.m-2], < 0 
    360360               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   
    361361               !  
    362                wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice                ! Mass flux 
     362               wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice                 ! Mass flux 
    363363                
    364364            END IF 
     
    366366            ! Ice sublimation 
    367367            ! --------------- 
    368             zdum            = MAX( - ( zh_i(ji,jk) + zdeltah(ji,jk) ) , - zevap_rema(ji) * r1_rhoic ) 
     368            zdum            = MAX( - ( zh_i(ji,jk) + zdeltah(ji,jk) ) , - zevap_rema(ji) * r1_rhoi ) 
    369369            zdeltah (ji,jk) = zdeltah (ji,jk) + zdum 
    370370            dh_i_sub(ji)    = dh_i_sub(ji)    + zdum 
    371371             
    372             sfx_sub_1d(ji)     = sfx_sub_1d(ji) - rhoic * a_i_1d(ji) * zdum * s_i_1d(ji) * r1_rdtice ! Salt flux >0 
     372            sfx_sub_1d(ji)     = sfx_sub_1d(ji) - rhoi * a_i_1d(ji) * zdum * s_i_1d(ji) * r1_rdtice ! Salt flux >0 
    373373            !                                                                                          clem: flux is sent to the ocean for simplicity 
    374374            !                                                                                                but salt should remain in the ice except 
     
    376376            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 
    377377 
    378             wfx_ice_sub_1d(ji) = wfx_ice_sub_1d(ji) - rhoic * a_i_1d(ji) * zdum * r1_rdtice          ! Mass flux > 0 
     378            wfx_ice_sub_1d(ji) = wfx_ice_sub_1d(ji) - rhoi * a_i_1d(ji) * zdum * r1_rdtice           ! Mass flux > 0 
    379379 
    380380            ! update remaining mass flux 
    381             zevap_rema(ji)  = zevap_rema(ji) + zdum * rhoic 
     381            zevap_rema(ji)  = zevap_rema(ji) + zdum * rhoi 
    382382             
    383383            ! record which layers have disappeared (for bottom melting)  
     
    438438               s_i_new(ji)   = zswitch_sal * zfracs * sss_1d(ji) + ( 1. - zswitch_sal ) * s_i_1d(ji)  ! New ice salinity 
    439439 
    440                ztmelts       = - tmut * s_i_new(ji)                                                   ! New ice melting point (C) 
     440               ztmelts       = - rTmlt * s_i_new(ji)                                                  ! New ice melting point (C) 
    441441 
    442442               zt_i_new      = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i) 
    443443                
    444                zEi           = cpic * ( zt_i_new - (ztmelts+rt0) ) &                                  ! Specific enthalpy of forming ice (J/kg, <0) 
    445                   &            - lfus * ( 1.0 - ztmelts / ( zt_i_new - rt0 ) ) + rcp  * ztmelts 
     444               zEi           = rcpi * ( zt_i_new - (ztmelts+rt0) ) &                                  ! Specific enthalpy of forming ice (J/kg, <0) 
     445                  &            - rLfus * ( 1.0 - ztmelts / ( zt_i_new - rt0 ) ) + rcp  * ztmelts 
    446446 
    447447               zEw           = rcp  * ( t_bo_1d(ji) - rt0 )                                           ! Specific enthalpy of seawater (J/kg, < 0) 
     
    449449               zdE           = zEi - zEw                                                              ! Specific enthalpy difference (J/kg, <0) 
    450450 
    451                dh_i_bog(ji)  = rdt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoic ) ) 
     451               dh_i_bog(ji)  = rdt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoi ) ) 
    452452                
    453453            END DO 
    454454            ! Contribution to Energy and Salt Fluxes                                     
    455             zfmdt          = - rhoic * dh_i_bog(ji)                                                   ! Mass flux x time step (kg/m2, < 0) 
     455            zfmdt          = - rhoi * dh_i_bog(ji)                                                   ! Mass flux x time step (kg/m2, < 0) 
    456456             
    457457            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 
    458458            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 
    459459             
    460             sfx_bog_1d(ji) = sfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bog(ji) * s_i_new(ji) * r1_rdtice    ! Salt flux, <0 
    461  
    462             wfx_bog_1d(ji) = wfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bog(ji) * r1_rdtice                  ! Mass flux, <0 
     460            sfx_bog_1d(ji) = sfx_bog_1d(ji) - rhoi * a_i_1d(ji) * dh_i_bog(ji) * s_i_new(ji) * r1_rdtice     ! Salt flux, <0 
     461 
     462            wfx_bog_1d(ji) = wfx_bog_1d(ji) - rhoi * a_i_1d(ji) * dh_i_bog(ji) * r1_rdtice                   ! Mass flux, <0 
    463463 
    464464            ! update heat content (J.m-2) and layer thickness 
    465             eh_i_old(ji,nlay_i+1) = eh_i_old(ji,nlay_i+1) + dh_i_bog(ji) * (-zEi * rhoic) 
     465            eh_i_old(ji,nlay_i+1) = eh_i_old(ji,nlay_i+1) + dh_i_bog(ji) * (-zEi * rhoi) 
    466466            h_i_old (ji,nlay_i+1) = h_i_old (ji,nlay_i+1) + dh_i_bog(ji) 
    467467 
     
    477477            IF(  zf_tt(ji)  >  0._wp  .AND. jk > icount(ji,jk) ) THEN   ! do not calculate where layer has already disappeared by surface melting  
    478478 
    479                ztmelts = - tmut * sz_i_1d(ji,jk)  ! Melting point of layer jk (C) 
     479               ztmelts = - rTmlt * sz_i_1d(ji,jk)  ! Melting point of layer jk (C) 
    480480 
    481481               IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN   !-- Internal melting 
    482482 
    483                   zEi               = - e_i_1d(ji,jk) * r1_rhoic    ! Specific enthalpy of melting ice (J/kg, <0) 
     483                  zEi               = - e_i_1d(ji,jk) * r1_rhoi     ! Specific enthalpy of melting ice (J/kg, <0) 
    484484                  zdE               = 0._wp                         ! Specific enthalpy difference   (J/kg, <0) 
    485485                                                                    !    set up at 0 since no energy is needed to melt water...(it is already melted) 
     
    489489                  dh_i_itm (ji)     = dh_i_itm(ji) + zdeltah(ji,jk) 
    490490 
    491                   zfmdt             = - zdeltah(ji,jk) * rhoic      ! Mass flux x time step > 0 
     491                  zfmdt             = - zdeltah(ji,jk) * rhoi      ! Mass flux x time step > 0 
    492492 
    493493                  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 
    494494                  !                                                                                                  ice enthalpy zEi is "sent" to the ocean 
    495                   sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice   ! Salt flux 
     495                  sfx_res_1d(ji) = sfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice    ! Salt flux 
    496496                  !                                                                                                  using s_i_1d and not sz_i_1d(jk) is ok 
    497                   wfx_res_1d(ji) = wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice                ! Mass flux 
     497                  wfx_res_1d(ji) = wfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice                 ! Mass flux 
    498498 
    499499                  ! update heat content (J.m-2) and layer thickness 
     
    503503               ELSE                                        !-- Basal melting 
    504504 
    505                   zEi             = - e_i_1d(ji,jk) * r1_rhoic                                ! Specific enthalpy of melting ice (J/kg, <0) 
     505                  zEi             = - e_i_1d(ji,jk) * r1_rhoi                                 ! Specific enthalpy of melting ice (J/kg, <0) 
    506506                  zEw             = rcp * ztmelts                                             ! Specific enthalpy of meltwater (J/kg, <0) 
    507507                  zdE             = zEi - zEw                                                 ! Specific enthalpy difference   (J/kg, <0) 
     
    509509                  zfmdt           = - zq_bot(ji) / zdE                                        ! Mass flux x time step (kg/m2, >0) 
    510510 
    511                   zdeltah(ji,jk)  = - zfmdt * r1_rhoic                                        ! Gross thickness change 
     511                  zdeltah(ji,jk)  = - zfmdt * r1_rhoi                                         ! Gross thickness change 
    512512 
    513513                  zdeltah(ji,jk)  = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) )       ! bound thickness change 
    514514                   
    515                   zq_bot(ji)      = MAX( 0._wp , zq_bot(ji) - zdeltah(ji,jk) * rhoic * zdE )  ! update available heat. MAX is necessary for roundup errors 
     515                  zq_bot(ji)      = MAX( 0._wp , zq_bot(ji) - zdeltah(ji,jk) * rhoi * zdE )   ! update available heat. MAX is necessary for roundup errors 
    516516 
    517517                  dh_i_bom(ji)    = dh_i_bom(ji) + zdeltah(ji,jk)                             ! Update basal melt 
    518518 
    519                   zfmdt           = - zdeltah(ji,jk) * rhoic                                  ! Mass flux x time step > 0 
     519                  zfmdt           = - zdeltah(ji,jk) * rhoi                                   ! Mass flux x time step > 0 
    520520 
    521521                  zQm             = zfmdt * zEw                                               ! Heat exchanged with ocean 
     
    524524                  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   
    525525 
    526                   sfx_bom_1d(ji)  = sfx_bom_1d(ji) - rhoic *  a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice  ! Salt flux 
     526                  sfx_bom_1d(ji)  = sfx_bom_1d(ji) - rhoi *  a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice   ! Salt flux 
    527527                  !                                                                                                   using s_i_1d and not sz_i_1d(jk) is ok 
    528528                   
    529                   wfx_bom_1d(ji)  = wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice                ! Mass flux 
     529                  wfx_bom_1d(ji)  = wfx_bom_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice                 ! Mass flux 
    530530 
    531531                  ! update heat content (J.m-2) and layer thickness 
     
    558558         zq_rema(ji)        = zq_rema(ji) + zdeltah(ji,1) * e_s_1d(ji,1)                               ! update available heat (J.m-2) 
    559559         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) 
    560          wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice      ! Mass flux 
     560         wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhos * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice       ! Mass flux 
    561561         dh_s_mlt(ji)       = dh_s_mlt(ji) + zdeltah(ji,1) 
    562562         !     
     
    572572      ! When snow load excesses Archimede's limit, snow-ice interface goes down under sea-level,  
    573573      ! flooding of seawater transforms snow into ice dh_snowice is positive for the ice 
    574       z1_rho = 1._wp / ( rhosn+rau0-rhoic ) 
     574      z1_rho = 1._wp / ( rhos+rau0-rhoi ) 
    575575      DO ji = 1, npti 
    576576         ! 
    577          dh_snowice(ji) = MAX(  0._wp , ( rhosn * h_s_1d(ji) + (rhoic-rau0) * h_i_1d(ji) ) * z1_rho ) 
     577         dh_snowice(ji) = MAX(  0._wp , ( rhos * h_s_1d(ji) + (rhoi-rau0) * h_i_1d(ji) ) * z1_rho ) 
    578578 
    579579         h_i_1d(ji)    = h_i_1d(ji) + dh_snowice(ji) 
     
    581581 
    582582         ! Contribution to energy flux to the ocean [J/m2], >0 (if sst<0) 
    583          zfmdt          = ( rhosn - rhoic ) * dh_snowice(ji)    ! <0 
     583         zfmdt          = ( rhos - rhoi ) * dh_snowice(ji)    ! <0 
    584584         zEw            = rcp * sst_1d(ji) 
    585585         zQm            = zfmdt * zEw  
     
    592592         IF( nn_icesal /= 2 )  THEN 
    593593            sfx_bri_1d(ji) = sfx_bri_1d(ji) - sss_1d (ji) * a_i_1d(ji) * zfmdt                  * r1_rdtice  & ! put back sss_m     into the ocean 
    594                &                            - s_i_1d(ji)  * a_i_1d(ji) * dh_snowice(ji) * rhoic * r1_rdtice    ! and get  rn_icesal from the ocean  
     594               &                            - s_i_1d(ji)  * a_i_1d(ji) * dh_snowice(ji) * rhoi * r1_rdtice     ! and get  rn_icesal from the ocean  
    595595         ENDIF 
    596596 
    597597         ! Mass flux: All snow is thrown in the ocean, and seawater is taken to replace the volume 
    598          wfx_sni_1d(ji)     = wfx_sni_1d(ji)     - a_i_1d(ji) * dh_snowice(ji) * rhoic * r1_rdtice 
    599          wfx_snw_sni_1d(ji) = wfx_snw_sni_1d(ji) + a_i_1d(ji) * dh_snowice(ji) * rhosn * r1_rdtice 
     598         wfx_sni_1d(ji)     = wfx_sni_1d(ji)     - a_i_1d(ji) * dh_snowice(ji) * rhoi * r1_rdtice 
     599         wfx_snw_sni_1d(ji) = wfx_snw_sni_1d(ji) + a_i_1d(ji) * dh_snowice(ji) * rhos * r1_rdtice 
    600600 
    601601         ! update heat content (J.m-2) and layer thickness 
     
    619619            e_s_1d(ji,jk) = rswitch * e_s_1d(ji,jk) 
    620620            ! recalculate t_s_1d from e_s_1d 
    621             t_s_1d(ji,jk) = rt0 + rswitch * ( - e_s_1d(ji,jk) * r1_rhosn * r1_cpic + lfus * r1_cpic ) 
     621            t_s_1d(ji,jk) = rt0 + rswitch * ( - e_s_1d(ji,jk) * r1_rhos * r1_rcpi + rLfus * r1_rcpi ) 
    622622         END DO 
    623623      END DO 
  • NEMO/trunk/src/ICE/icethd_do.F90

    r9604 r9935  
    140140         ! Physical constants 
    141141         zhicrit = 0.04                                          ! frazil ice thickness 
    142          ztwogp  = 2. * rau0 / ( grav * 0.3 * ( rau0 - rhoic ) ) ! reduced grav 
     142         ztwogp  = 2. * rau0 / ( grav * 0.3 * ( rau0 - rhoi ) ) ! reduced grav 
    143143         zsqcd   = 1.0 / SQRT( 1.3 * zcai )                      ! 1/SQRT(airdensity*drag) 
    144144         zgamafr = 0.03 
     
    263263         ! We assume that new ice is formed at the seawater freezing point 
    264264         DO ji = 1, npti 
    265             ztmelts       = - tmut * zs_newice(ji)                  ! Melting point (C) 
    266             ze_newice(ji) =   rhoic * (  cpic * ( ztmelts - ( t_bo_1d(ji) - rt0 ) )                     & 
    267                &                       + lfus * ( 1.0 - ztmelts / MIN( t_bo_1d(ji) - rt0, -epsi10 ) )   & 
    268                &                       - rcp  *         ztmelts ) 
     265            ztmelts       = - rTmlt * zs_newice(ji)                  ! Melting point (C) 
     266            ze_newice(ji) =   rhoi * (  rcpi * ( ztmelts - ( t_bo_1d(ji) - rt0 ) )                     & 
     267               &                      + rLfus * ( 1.0 - ztmelts / MIN( t_bo_1d(ji) - rt0, -epsi10 ) )   & 
     268               &                      - rcp   *         ztmelts ) 
    269269         END DO 
    270270 
     
    275275         DO ji = 1, npti 
    276276 
    277             zEi           = - ze_newice(ji) * r1_rhoic             ! specific enthalpy of forming ice [J/kg] 
     277            zEi           = - ze_newice(ji) * r1_rhoi              ! specific enthalpy of forming ice [J/kg] 
    278278 
    279279            zEw           = rcp * ( t_bo_1d(ji) - rt0 )            ! specific enthalpy of seawater at t_bo_1d [J/kg] 
     
    284284            zfmdt         = - qlead_1d(ji) / zdE                   ! Fm.dt [kg/m2] (<0)  
    285285                                                                   ! clem: we use qlead instead of zqld (icethd) because we suppose we are at the freezing point    
    286             zv_newice(ji) = - zfmdt * r1_rhoic 
     286            zv_newice(ji) = - zfmdt * r1_rhoi 
    287287 
    288288            zQm           = zfmdt * zEw                            ! heat to the ocean >0 associated with mass flux   
     
    293293            hfx_opw_1d(ji) = hfx_opw_1d(ji) - zfmdt * zdE * r1_rdtice 
    294294            ! mass flux 
    295             wfx_opw_1d(ji) = wfx_opw_1d(ji) - zv_newice(ji) * rhoic * r1_rdtice 
     295            wfx_opw_1d(ji) = wfx_opw_1d(ji) - zv_newice(ji) * rhoi * r1_rdtice 
    296296            ! salt flux 
    297             sfx_opw_1d(ji) = sfx_opw_1d(ji) - zv_newice(ji) * rhoic * zs_newice(ji) * r1_rdtice 
     297            sfx_opw_1d(ji) = sfx_opw_1d(ji) - zv_newice(ji) * rhoi * zs_newice(ji) * r1_rdtice 
    298298         END DO 
    299299          
  • NEMO/trunk/src/ICE/icethd_pnd.F90

    r9750 r9935  
    133133      REAL(wp) ::   zdv_mlt          ! available meltwater for melt ponding 
    134134      REAL(wp) ::   z1_Tp            ! inverse reference temperature 
    135       REAL(wp) ::   z1_rhofw         ! inverse freshwater density 
     135      REAL(wp) ::   z1_rhow          ! inverse freshwater density 
    136136      REAL(wp) ::   z1_zpnd_aspect   ! inverse pond aspect ratio 
    137137      REAL(wp) ::   zfac, zdum 
     
    139139      INTEGER  ::   ji   ! loop indices 
    140140      !!------------------------------------------------------------------- 
    141       z1_rhofw       = 1._wp / rhofw  
     141      z1_rhow        = 1._wp / rhow  
    142142      z1_zpnd_aspect = 1._wp / zpnd_aspect 
    143143      z1_Tp          = 1._wp / zTp  
     
    157157            ! 
    158158            ! available meltwater for melt ponding [m, >0] and fraction 
    159             zdv_mlt = -( dh_i_sum(ji)*rhoic + dh_s_mlt(ji)*rhosn ) * z1_rhofw * a_i_1d(ji) 
     159            zdv_mlt = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow * a_i_1d(ji) 
    160160            zfr_mlt = zrmin + ( zrmax - zrmin ) * a_i_1d(ji)  ! from CICE doc 
    161161            !zfr_mlt = zrmin + zrmax * a_i_1d(ji)             ! from Holland paper  
     
    168168            ! melt pond mass flux (<0) 
    169169            IF( ln_pnd_fwb .AND. zdv_mlt > 0._wp ) THEN 
    170                zfac = zfr_mlt * zdv_mlt * rhofw * r1_rdtice 
     170               zfac = zfr_mlt * zdv_mlt * rhow * r1_rdtice 
    171171               wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac 
    172172               ! 
  • NEMO/trunk/src/ICE/icethd_sal.F90

    r9750 r9935  
    7777            !--------------------------------------------------------- 
    7878            IF( h_i_1d(ji) > 0._wp ) THEN 
    79                zs_sni   = sss_1d(ji) * ( rhoic - rhosn ) * r1_rhoic                 ! Salinity of snow ice 
     79               zs_sni  = sss_1d(ji) * ( rhoi - rhos ) * r1_rhoi                     ! Salinity of snow ice 
    8080               zs_i_si = ( zs_sni      - s_i_1d(ji) ) * dh_snowice(ji) / h_i_1d(ji) ! snow-ice     
    8181               zs_i_bg = ( s_i_new(ji) - s_i_1d(ji) ) * dh_i_bog  (ji) / h_i_1d(ji) ! bottom growth 
     
    9898                
    9999               ! Salt flux 
    100                sfx_bri_1d(ji) = sfx_bri_1d(ji) - rhoic * a_i_1d(ji) * h_i_1d(ji) * ( zs_i_fl + zs_i_gd ) * r1_rdtice 
     100               sfx_bri_1d(ji) = sfx_bri_1d(ji) - rhoi * a_i_1d(ji) * h_i_1d(ji) * ( zs_i_fl + zs_i_gd ) * r1_rdtice 
    101101            ENDIF 
    102102         END DO 
  • NEMO/trunk/src/ICE/icethd_zdf_bl99.F90

    r9929 r9935  
    9393      REAL(wp) ::   zdti_bnd  =  1.e-4_wp     ! maximal authorized error on temperature  
    9494      REAL(wp) ::   zhs_min   =  0.01_wp      ! minimum snow thickness for conductivity calculation  
    95       REAL(wp) ::   ztmelt_i                  ! ice melting temperature 
     95      REAL(wp) ::   ztmelts                   ! ice melting temperature 
    9696      REAL(wp) ::   zdti_max                  ! current maximal error on temperature  
    9797      REAL(wp) ::   zcpi                      ! Ice specific heat 
     
    217217            ! 
    218218            DO ji = 1, npti 
    219                ztcond_i(ji,0)      = rcdic + zbeta * sz_i_1d(ji,1)      / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) 
    220                ztcond_i(ji,nlay_i) = rcdic + zbeta * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji)  - rt0 ) 
     219               ztcond_i(ji,0)      = rcnd_i + zbeta * sz_i_1d(ji,1)      / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) 
     220               ztcond_i(ji,nlay_i) = rcnd_i + zbeta * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji)  - rt0 ) 
    221221            END DO 
    222222            DO jk = 1, nlay_i-1 
    223223               DO ji = 1, npti 
    224                   ztcond_i(ji,jk) = rcdic + zbeta * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) /  & 
    225                      &                      MIN( -epsi10, 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 ) 
     224                  ztcond_i(ji,jk) = rcnd_i + zbeta * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) /  & 
     225                     &                       MIN( -epsi10, 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 ) 
    226226               END DO 
    227227            END DO 
     
    230230            ! 
    231231            DO ji = 1, npti 
    232                ztcond_i(ji,0)      = rcdic + 0.09_wp  *  sz_i_1d(ji,1)      / MIN( -epsi10, t_i_1d(ji,1) - rt0 )  & 
    233                   &                        - 0.011_wp * ( t_i_1d(ji,1) - rt0 ) 
    234                ztcond_i(ji,nlay_i) = rcdic + 0.09_wp  *  sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji)  - rt0 )  & 
    235                   &                        - 0.011_wp * ( t_bo_1d(ji) - rt0 ) 
     232               ztcond_i(ji,0)      = rcnd_i + 0.09_wp  *  sz_i_1d(ji,1)      / MIN( -epsi10, t_i_1d(ji,1) - rt0 )  & 
     233                  &                         - 0.011_wp * ( t_i_1d(ji,1) - rt0 ) 
     234               ztcond_i(ji,nlay_i) = rcnd_i + 0.09_wp  *  sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji)  - rt0 )  & 
     235                  &                         - 0.011_wp * ( t_bo_1d(ji) - rt0 ) 
    236236            END DO 
    237237            DO jk = 1, nlay_i-1 
    238238               DO ji = 1, npti 
    239                   ztcond_i(ji,jk) = rcdic + 0.09_wp  *   0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) /        & 
    240                      &                     MIN( -epsi10, 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 )  & 
    241                      &                    - 0.011_wp * ( 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 ) 
     239                  ztcond_i(ji,jk) = rcnd_i + 0.09_wp  *   0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) /        & 
     240                     &                      MIN( -epsi10, 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 )  & 
     241                     &                     - 0.011_wp * ( 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 ) 
    242242               END DO 
    243243            END DO 
     
    299299         DO jk = 1, nlay_i 
    300300            DO ji = 1, npti 
    301                zcpi = cpic + zgamma * sz_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztiold(ji,jk) - rt0 ), epsi10 ) 
    302                zeta_i(ji,jk) = rdt_ice * r1_rhoic * z1_h_i(ji) / MAX( epsi10, zcpi )  
     301               zcpi = rcpi + zgamma * sz_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztiold(ji,jk) - rt0 ), epsi10 ) 
     302               zeta_i(ji,jk) = rdt_ice * r1_rhoi * z1_h_i(ji) / MAX( epsi10, zcpi )  
    303303            END DO 
    304304         END DO 
     
    306306         DO jk = 1, nlay_s 
    307307            DO ji = 1, npti 
    308                zeta_s(ji,jk) = rdt_ice * r1_rhosn * r1_cpic * z1_h_s(ji) 
     308               zeta_s(ji,jk) = rdt_ice * r1_rhos * r1_rcpi * z1_h_s(ji) 
    309309            END DO 
    310310         END DO 
     
    539539            DO jk = 1, nlay_i 
    540540               DO ji = 1, npti 
    541                   ztmelt_i      = -tmut * sz_i_1d(ji,jk) + rt0  
    542                   t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), rt0 - 100._wp ) 
     541                  ztmelts       = -rTmlt * sz_i_1d(ji,jk) + rt0  
     542                  t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelts ), rt0 - 100._wp ) 
    543543                  zdti_max      =  MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) 
    544544               END DO 
     
    704704            DO jk = 1, nlay_i 
    705705               DO ji = 1, npti 
    706                   ztmelt_i      = -tmut * sz_i_1d(ji,jk) + rt0  
    707                   t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), rt0 - 100._wp ) 
     706                  ztmelts       = -rTmlt * sz_i_1d(ji,jk) + rt0  
     707                  t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelts ), rt0 - 100._wp ) 
    708708                  zdti_max      =  MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) 
    709709               END DO 
  • NEMO/trunk/src/ICE/iceupdate.F90

    r9922 r9935  
    171171            snwice_mass_b(ji,jj) = snwice_mass(ji,jj)       ! save mass from the previous ice time step 
    172172            !                                               ! new mass per unit area 
    173             snwice_mass  (ji,jj) = tmask(ji,jj,1) * ( rhosn * vt_s(ji,jj) + rhoic * vt_i(ji,jj)  )  
     173            snwice_mass  (ji,jj) = tmask(ji,jj,1) * ( rhos * vt_s(ji,jj) + rhoi * vt_i(ji,jj)  )  
    174174            !                                               ! time evolution of snow+ice mass 
    175175            snwice_fmass (ji,jj) = ( snwice_mass(ji,jj) - snwice_mass_b(ji,jj) ) * r1_rdtice 
     
    431431            ELSE                                     ! start from rest 
    432432               IF(lwp) WRITE(numout,*) '   ==>>   previous run without snow-ice mass output then set it' 
    433                snwice_mass  (:,:) = tmask(:,:,1) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:) ) 
     433               snwice_mass  (:,:) = tmask(:,:,1) * ( rhos * vt_s(:,:) + rhoi * vt_i(:,:) ) 
    434434               snwice_mass_b(:,:) = snwice_mass(:,:) 
    435435            ENDIF 
    436436         ELSE                                   !* Start from rest 
    437437            IF(lwp) WRITE(numout,*) '   ==>>   start from rest: set the snow-ice mass' 
    438             snwice_mass  (:,:) = tmask(:,:,1) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:) ) 
     438            snwice_mass  (:,:) = tmask(:,:,1) * ( rhos * vt_s(:,:) + rhoi * vt_i(:,:) ) 
    439439            snwice_mass_b(:,:) = snwice_mass(:,:) 
    440440         ENDIF 
  • NEMO/trunk/src/ICE/icevar.F90

    r9888 r9935  
    226226                     ! 
    227227                     ze_i             =   e_i (ji,jj,jk,jl) * z1_v_i(ji,jj,jl) * zlay_i               ! Energy of melting e(S,T) [J.m-3] 
    228                      ztmelts          = - sz_i(ji,jj,jk,jl) * tmut                                 ! Ice layer melt temperature [C] 
     228                     ztmelts          = - sz_i(ji,jj,jk,jl) * rTmlt                                 ! Ice layer melt temperature [C] 
    229229                     ! Conversion q(S,T) -> T (second order equation) 
    230                      zbbb             = ( rcp - cpic ) * ztmelts + ze_i * r1_rhoic - lfus 
    231                      zccc             = SQRT( MAX( zbbb * zbbb - 4._wp * cpic * lfus * ztmelts , 0._wp) ) 
    232                      t_i(ji,jj,jk,jl) = MAX( -100._wp , MIN( -( zbbb + zccc ) * 0.5_wp * r1_cpic , ztmelts ) ) + rt0   ! [K] with bounds: -100 < t_i < ztmelts 
     230                     zbbb             = ( rcp - rcpi ) * ztmelts + ze_i * r1_rhoi - rLfus 
     231                     zccc             = SQRT( MAX( zbbb * zbbb - 4._wp * rcpi * rLfus * ztmelts , 0._wp) ) 
     232                     t_i(ji,jj,jk,jl) = MAX( -100._wp , MIN( -( zbbb + zccc ) * 0.5_wp * r1_rcpi , ztmelts ) ) + rt0   ! [K] with bounds: -100 < t_i < ztmelts 
    233233                     ! 
    234234                  ELSE                                !--- no ice 
     
    247247         WHERE( v_s(:,:,:) > epsi20 )        !--- icy area 
    248248            t_s(:,:,jk,:) = rt0 + MAX( -100._wp ,  & 
    249                  &                MIN( r1_cpic * ( -r1_rhosn * ( e_s(:,:,jk,:) / v_s(:,:,:) * zlay_s ) + lfus ) , 0._wp ) ) 
     249                 &                MIN( r1_rcpi * ( -r1_rhos * ( e_s(:,:,jk,:) / v_s(:,:,:) * zlay_s ) + rLfus ) , 0._wp ) ) 
    250250         ELSEWHERE                           !--- no ice 
    251251            t_s(:,:,jk,:) = rt0 
     
    498498            DO ji = 1 , jpi 
    499499               ! update exchanges with ocean 
    500                sfx_res(ji,jj)  = sfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * sv_i(ji,jj,jl)   * rhoic * r1_rdtice 
    501                wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_i (ji,jj,jl)   * rhoic * r1_rdtice 
    502                wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl)   * rhosn * r1_rdtice 
     500               sfx_res(ji,jj)  = sfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * sv_i(ji,jj,jl)   * rhoi * r1_rdtice 
     501               wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_i (ji,jj,jl)   * rhoi * r1_rdtice 
     502               wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl)   * rhos * r1_rdtice 
    503503               ! 
    504504               !----------------------------------------------------------------- 
     
    669669               ! In case snow load is in excess that would lead to transformation from snow to ice 
    670670               ! Then, transfer the snow excess into the ice (different from icethd_dh) 
    671                zdh = MAX( 0._wp, ( rhosn * zh_s(ji,jl) + ( rhoic - rau0 ) * zh_i(ji,jl) ) * r1_rau0 )  
     671               zdh = MAX( 0._wp, ( rhos * zh_s(ji,jl) + ( rhoi - rau0 ) * zh_i(ji,jl) ) * r1_rau0 )  
    672672               ! recompute h_i, h_s avoiding out of bounds values 
    673673               zh_i(ji,jl) = MIN( hi_max(jl), zh_i(ji,jl) + zdh ) 
    674                zh_s(ji,jl) = MAX( 0._wp, zh_s(ji,jl) - zdh * rhoic * r1_rhosn ) 
     674               zh_s(ji,jl) = MAX( 0._wp, zh_s(ji,jl) - zdh * rhoi * r1_rhos ) 
    675675            ENDIF 
    676676         END DO 
     
    827827         DO jk = 1, nlay_i 
    828828            WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 )    
    829                bv_i(:,:,jl) = bv_i(:,:,jl) - tmut * sz_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 ) 
     829               bv_i(:,:,jl) = bv_i(:,:,jl) - rTmlt * sz_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 ) 
    830830            END WHERE 
    831831         END DO 
     
    852852      DO jk = 1, nlay_i             ! Sea ice energy of melting 
    853853         DO ji = 1, npti 
    854             ztmelts      = - tmut  * sz_i_1d(ji,jk) 
     854            ztmelts      = - rTmlt  * sz_i_1d(ji,jk) 
    855855            t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts + rt0 ) ! Force t_i_1d to be lower than melting point => likely conservation issue 
    856856                                                                !   (sometimes zdf scheme produces abnormally high temperatures)    
    857             e_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - ( t_i_1d(ji,jk) - rt0 ) )           & 
    858                &                    + lfus * ( 1._wp - ztmelts / ( t_i_1d(ji,jk) - rt0 ) )   & 
    859                &                    - rcp  *  ztmelts ) 
     857            e_i_1d(ji,jk) = rhoi * ( rcpi * ( ztmelts - ( t_i_1d(ji,jk) - rt0 ) )           & 
     858               &                   + rLfus * ( 1._wp - ztmelts / ( t_i_1d(ji,jk) - rt0 ) )   & 
     859               &                   - rcp   * ztmelts ) 
    860860         END DO 
    861861      END DO 
    862862      DO jk = 1, nlay_s             ! Snow energy of melting 
    863863         DO ji = 1, npti 
    864             e_s_1d(ji,jk) = rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) 
     864            e_s_1d(ji,jk) = rhos * ( rcpi * ( rt0 - t_s_1d(ji,jk) ) + rLfus ) 
    865865         END DO 
    866866      END DO 
  • NEMO/trunk/src/ICE/icewri.F90

    r9604 r9935  
    8585      ! Standard outputs 
    8686      !----------------- 
    87       zrho1 = ( rau0 - rhoic ) * r1_rau0; zrho2 = rhosn * r1_rau0 
     87      zrho1 = ( rau0 - rhoi ) * r1_rau0; zrho2 = rhos * r1_rau0 
    8888      ! masks 
    8989      IF( iom_use('icemask'  ) )   CALL iom_put( "icemask"  , zmsk00              )   ! ice mask 0% 
     
    9292      ! 
    9393      ! general fields 
    94       IF( iom_use('icemass'  ) )   CALL iom_put( "icemass", rhoic * vt_i * zmsk00 )   ! Ice mass per cell area  
    95       IF( iom_use('snwmass'  ) )   CALL iom_put( "snwmass", rhosn * vt_s * zmsksn )   ! Snow mass per cell area 
     94      IF( iom_use('icemass'  ) )   CALL iom_put( "icemass", rhoi * vt_i * zmsk00 )   ! Ice mass per cell area  
     95      IF( iom_use('snwmass'  ) )   CALL iom_put( "snwmass", rhos * vt_s * zmsksn )   ! Snow mass per cell area 
    9696      IF( iom_use('icepres'  ) )   CALL iom_put( "icepres", zmsk00                )   ! Ice presence (1 or 0)  
    9797      IF( iom_use('iceconc'  ) )   CALL iom_put( "iceconc", at_i  * zmsk00        )   ! ice concentration 
     
    115115      ! salt 
    116116      IF( iom_use('icesalt'  ) )   CALL iom_put( "icesalt", sm_i  * zmsk00        )   ! mean ice salinity 
    117       IF( iom_use('icesalm'  ) )   CALL iom_put( "icesalm", SUM( sv_i, DIM = 3 ) * rhoic * 1.0e-3 * zmsk00 )   ! Mass of salt in sea ice per cell area 
     117      IF( iom_use('icesalm'  ) )   CALL iom_put( "icesalm", SUM( sv_i, DIM = 3 ) * rhoi * 1.0e-3 * zmsk00 )   ! Mass of salt in sea ice per cell area 
    118118 
    119119      ! heat 
     
    164164      ! trends 
    165165      IF( iom_use('dmithd') )   CALL iom_put( "dmithd", - wfx_bog - wfx_bom - wfx_sum - wfx_sni - wfx_opw - wfx_lam - wfx_res ) ! Sea-ice mass change from thermodynamics 
    166       IF( iom_use('dmidyn') )   CALL iom_put( "dmidyn", - wfx_dyn + rhoic * diag_trp_vi     )   ! Sea-ice mass change from dynamics(kg/m2/s) 
     166      IF( iom_use('dmidyn') )   CALL iom_put( "dmidyn", - wfx_dyn + rhoi * diag_trp_vi      )   ! Sea-ice mass change from dynamics(kg/m2/s) 
    167167      IF( iom_use('dmiopw') )   CALL iom_put( "dmiopw", - wfx_opw                           )   ! Sea-ice mass change through growth in open water 
    168168      IF( iom_use('dmibog') )   CALL iom_put( "dmibog", - wfx_bog                           )   ! Sea-ice mass change through basal growth 
     
    174174      IF( iom_use('dmisub') )   CALL iom_put( "dmisub", - wfx_ice_sub                       )   ! Sea-ice mass change through sublimation 
    175175      IF( iom_use('dmsspr') )   CALL iom_put( "dmsspr", - wfx_spr                           )   ! Snow mass change through snow fall 
    176       IF( iom_use('dmsssi') )   CALL iom_put( "dmsssi",   wfx_sni*rhosn*r1_rhoic            )   ! Snow mass change through snow-to-ice conversion 
     176      IF( iom_use('dmsssi') )   CALL iom_put( "dmsssi",   wfx_sni*rhos*r1_rhoi              )   ! Snow mass change through snow-to-ice conversion 
    177177      IF( iom_use('dmsmel') )   CALL iom_put( "dmsmel", - wfx_snw_sum                       )   ! Snow mass change through melt 
    178       IF( iom_use('dmsdyn') )   CALL iom_put( "dmsdyn", - wfx_snw_dyn + rhosn * diag_trp_vs )   ! Snow mass change through dynamics(kg/m2/s) 
     178      IF( iom_use('dmsdyn') )   CALL iom_put( "dmsdyn", - wfx_snw_dyn + rhos * diag_trp_vs )   ! Snow mass change through dynamics(kg/m2/s) 
    179179 
    180180      ! Global ice diagnostics 
  • NEMO/trunk/src/OCE/BDY/bdyice.F90

    r9927 r9935  
    125125 
    126126            ! Then, a) transfer the snow excess into the ice (different from icethd_dh) 
    127             zdh = MAX( 0._wp, ( rhosn * h_s(ji,jj,jl) + ( rhoic - rau0 ) * h_i(ji,jj,jl) ) * r1_rau0 ) 
     127            zdh = MAX( 0._wp, ( rhos * h_s(ji,jj,jl) + ( rhoi - rau0 ) * h_i(ji,jj,jl) ) * r1_rau0 ) 
    128128            ! Or, b) transfer all the snow into ice (if incoming ice is likely to melt as it comes into a warmer environment) 
    129             !zdh = MAX( 0._wp, h_s(ji,jj,jl) * rhosn / rhoic ) 
     129            !zdh = MAX( 0._wp, h_s(ji,jj,jl) * rhos / rhoi ) 
    130130 
    131131            ! recompute h_i, h_s 
    132132            h_i(ji,jj,jl) = MIN( hi_max(jl), h_i(ji,jj,jl) + zdh ) 
    133             h_s(ji,jj,jl) = MAX( 0._wp, h_s(ji,jj,jl) - zdh * rhoic / rhosn )  
     133            h_s(ji,jj,jl) = MAX( 0._wp, h_s(ji,jj,jl) - zdh * rhoi / rhos )  
    134134 
    135135         ENDDO 
     
    198198               sv_i(ji,jj,jl) = MIN( s_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) ! salt content 
    199199               DO jk = 1, nlay_s 
    200                   e_s(ji,jj,jk,jl) = rhosn * ( cpic * ( rt0 - t_s(ji,jj,jk,jl) ) + lfus )   ! enthalpy in J/m3 
     200                  e_s(ji,jj,jk,jl) = rhos * ( rcpi * ( rt0 - t_s(ji,jj,jk,jl) ) + rLfus )   ! enthalpy in J/m3 
    201201                  e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) * r1_nlay_s           ! enthalpy in J/m2 
    202202               END DO                
    203203               DO jk = 1, nlay_i 
    204                   ztmelts          = - tmut  * sz_i(ji,jj,jk,jl)              ! Melting temperature in C 
     204                  ztmelts          = - rTmlt  * sz_i(ji,jj,jk,jl)             ! Melting temperature in C 
    205205                  t_i(ji,jj,jk,jl) = MIN( t_i(ji,jj,jk,jl), ztmelts + rt0 )   ! Force t_i to be lower than melting point => likely conservation issue 
    206206                  ! 
    207                   e_i(ji,jj,jk,jl) = rhoic * ( cpic * ( ztmelts - ( t_i(ji,jj,jk,jl) - rt0 ) )           &   ! enthalpy in J/m3 
    208                      &                       + lfus * ( 1._wp - ztmelts / ( t_i(ji,jj,jk,jl) - rt0 ) )   & 
    209                      &                       - rcp  *   ztmelts )                   
     207                  e_i(ji,jj,jk,jl) = rhoi * ( rcpi * ( ztmelts - ( t_i(ji,jj,jk,jl) - rt0 ) )           &   ! enthalpy in J/m3 
     208                     &                      + rLfus * ( 1._wp - ztmelts / ( t_i(ji,jj,jk,jl) - rt0 ) )   & 
     209                     &                      - rcp   *   ztmelts )                   
    210210                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i                            ! enthalpy in J/m2 
    211211               END DO 
  • NEMO/trunk/src/OCE/DOM/phycst.F90

    r9929 r9935  
    4646   REAL(wp), PUBLIC ::   r1_rau0_rcp                 !: = 1. / ( rau0 * rcp ) 
    4747 
    48 !clem: not sure these are needed for cice    
    49 #if defined key_cice 
    50    REAL(wp), PUBLIC ::   rt0_snow = 273.15_wp        !: melting point of snow         [Kelvin] 
    51    REAL(wp), PUBLIC ::   rt0_ice  = 273.05_wp        !: melting point of ice          [Kelvin] 
    52    REAL(wp), PUBLIC ::   emic     =    0.97_wp       !: emissivity of snow or ice 
    53    REAL(wp), PUBLIC ::   xlsn                        !: = lfus*rhosn (volumetric latent heat fusion of snow)  [J/m3] 
    54 #endif 
     48   REAL(wp), PUBLIC ::   emic     =    0.97_wp       !: emissivity of snow or ice (not used?) 
    5549 
    5650   REAL(wp), PUBLIC ::   sice     =    6.0_wp        !: salinity of ice (for pisces)          [psu] 
    5751   REAL(wp), PUBLIC ::   soce     =   34.7_wp        !: salinity of sea (for pisces and isf)  [psu] 
    58    REAL(wp), PUBLIC ::   cevap    =    2.5e+6_wp     !: latent heat of evaporation (water) 
    59    REAL(wp), PUBLIC ::   srgamma  =    0.9_wp        !: correction factor for solar radiation (Oberhuber, 1974) 
     52   REAL(wp), PUBLIC ::   rLevap   =    2.5e+6_wp     !: latent heat of evaporation (water) 
    6053   REAL(wp), PUBLIC ::   vkarmn   =    0.4_wp        !: von Karman constant 
    6154   REAL(wp), PUBLIC ::   stefan   =    5.67e-8_wp    !: Stefan-Boltzmann constant  
    6255 
    63    REAL(wp), PUBLIC ::   rhosn    =  330._wp         !: volumic mass of snow                                  [kg/m3] 
    64    REAL(wp), PUBLIC ::   rhoic    =  917._wp         !: volumic mass of sea ice                               [kg/m3] 
    65    REAL(wp), PUBLIC ::   rhofw    = 1000._wp         !: volumic mass of freshwater in melt ponds              [kg/m3] 
    66    REAL(wp), PUBLIC ::   rcdic    =    2.034396_wp   !: thermal conductivity of fresh ice                     [W/m/K] 
    67 #if defined key_cice 
    68    REAL(wp), PUBLIC ::   rcdsn    =    0.31_wp       !: thermal conductivity of snow                          [W/m/K] 
    69 #endif 
    70    REAL(wp), PUBLIC ::   cpic     = 2067.0_wp        !: specific heat of fresh ice                            [J/kg/K] 
    71    REAL(wp), PUBLIC ::   lsub     =    2.834e+6_wp   !: pure ice latent heat of sublimation                   [J/kg] 
    72    REAL(wp), PUBLIC ::   lfus     =    0.334e+6_wp   !: latent heat of fusion of fresh ice                    [J/kg] 
    73    REAL(wp), PUBLIC ::   tmut     =    0.054_wp      !: decrease of seawater meltpoint with salinity 
     56   REAL(wp), PUBLIC ::   rhos     =  330._wp         !: volumic mass of snow                                  [kg/m3] 
     57   REAL(wp), PUBLIC ::   rhoi     =  917._wp         !: volumic mass of sea ice                               [kg/m3] 
     58   REAL(wp), PUBLIC ::   rhow     = 1000._wp         !: volumic mass of freshwater in melt ponds              [kg/m3] 
     59   REAL(wp), PUBLIC ::   rcnd_i   =    2.034396_wp   !: thermal conductivity of fresh ice                     [W/m/K] 
     60   REAL(wp), PUBLIC ::   rcpi     = 2067.0_wp        !: specific heat of fresh ice                            [J/kg/K] 
     61   REAL(wp), PUBLIC ::   rLsub    =    2.834e+6_wp   !: pure ice latent heat of sublimation                   [J/kg] 
     62   REAL(wp), PUBLIC ::   rLfus    =    0.334e+6_wp   !: latent heat of fusion of fresh ice                    [J/kg] 
     63   REAL(wp), PUBLIC ::   rTmlt    =    0.054_wp      !: decrease of seawater meltpoint with salinity 
    7464 
    75    REAL(wp), PUBLIC ::   r1_rhoic                    !: 1 / rhoic 
    76    REAL(wp), PUBLIC ::   r1_rhosn                    !: 1 / rhosn 
    77    REAL(wp), PUBLIC ::   r1_cpic                     !: 1 / cpic 
     65   REAL(wp), PUBLIC ::   r1_rhoi                     !: 1 / rhoi 
     66   REAL(wp), PUBLIC ::   r1_rhos                     !: 1 / rhos 
     67   REAL(wp), PUBLIC ::   r1_rcpi                     !: 1 / rcpi 
    7868   !!---------------------------------------------------------------------- 
    7969   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    9989#endif 
    10090 
    101 #if defined key_cice 
    102       xlsn = lfus * rhosn 
    103 #endif 
    104  
    105       r1_rhoic = 1._wp / rhoic 
    106       r1_rhosn = 1._wp / rhosn 
    107       r1_cpic  = 1._wp / cpic 
     91      r1_rhoi = 1._wp / rhoi 
     92      r1_rhos = 1._wp / rhos 
     93      r1_rcpi = 1._wp / rcpi 
    10894 
    10995      IF(lwp) THEN 
     
    129115         WRITE(numout,*) '   reference density and heat capacity now defined in eosbn2.f90' 
    130116         WRITE(numout,*) 
    131 #if defined key_cice 
    132          WRITE(numout,*) '      thermal conductivity of the snow          = ', rcdsn   , ' J/s/m/K' 
    133 #endif 
    134          WRITE(numout,*) '      thermal conductivity of pure ice          = ', rcdic   , ' J/s/m/K' 
    135          WRITE(numout,*) '      fresh ice specific heat                   = ', cpic    , ' J/kg/K' 
    136          WRITE(numout,*) '      latent heat of fusion of fresh ice / snow = ', lfus    , ' J/kg' 
    137          WRITE(numout,*) '      latent heat of subl.  of fresh ice / snow = ', lsub    , ' J/kg' 
    138          WRITE(numout,*) '      density of sea ice                        = ', rhoic   , ' kg/m^3' 
    139          WRITE(numout,*) '      density of snow                           = ', rhosn   , ' kg/m^3' 
    140          WRITE(numout,*) '      density of freshwater (in melt ponds)     = ', rhofw   , ' kg/m^3' 
     117         WRITE(numout,*) '      thermal conductivity of pure ice          = ', rcnd_i  , ' J/s/m/K' 
     118         WRITE(numout,*) '      thermal conductivity of snow is defined in a namelist ' 
     119         WRITE(numout,*) '      fresh ice specific heat                   = ', rcpi    , ' J/kg/K' 
     120         WRITE(numout,*) '      latent heat of fusion of fresh ice / snow = ', rLfus   , ' J/kg' 
     121         WRITE(numout,*) '      latent heat of subl.  of fresh ice / snow = ', rLsub   , ' J/kg' 
     122         WRITE(numout,*) '      density of sea ice                        = ', rhoi    , ' kg/m^3' 
     123         WRITE(numout,*) '      density of snow                           = ', rhos    , ' kg/m^3' 
     124         WRITE(numout,*) '      density of freshwater (in melt ponds)     = ', rhow    , ' kg/m^3' 
    141125         WRITE(numout,*) '      salinity of ice (for pisces)              = ', sice    , ' psu' 
    142126         WRITE(numout,*) '      salinity of sea (for pisces and isf)      = ', soce    , ' psu' 
    143          WRITE(numout,*) '      latent heat of evaporation (water)        = ', cevap   , ' J/m^3'  
    144          WRITE(numout,*) '      correction factor for solar radiation     = ', srgamma  
     127         WRITE(numout,*) '      latent heat of evaporation (water)        = ', rLevap  , ' J/m^3'  
    145128         WRITE(numout,*) '      von Karman constant                       = ', vkarmn  
    146129         WRITE(numout,*) '      Stefan-Boltzmann constant                 = ', stefan  , ' J/s/m^2/K^4' 
  • NEMO/trunk/src/OCE/SBC/sbcblk.F90

    r9929 r9935  
    526526      ! 
    527527      qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                &   ! Downward Non Solar 
    528          &     - sf(jp_snow)%fnow(:,:,1) * rn_pfac * lfus                         &   ! remove latent melting heat for solid precip 
     528         &     - sf(jp_snow)%fnow(:,:,1) * rn_pfac * rLfus                        &   ! remove latent melting heat for solid precip 
    529529         &     - zevap(:,:) * pst(:,:) * rcp                                      &   ! remove evap heat content at SST 
    530530         &     + ( sf(jp_prec)%fnow(:,:,1) - sf(jp_snow)%fnow(:,:,1) ) * rn_pfac  &   ! add liquid precip heat content at Tair 
    531531         &     * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp                          & 
    532532         &     + sf(jp_snow)%fnow(:,:,1) * rn_pfac                                &   ! add solid  precip heat content at min(Tair,Tsnow) 
    533          &     * ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * cpic 
     533         &     * ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi 
    534534      qns(:,:) = qns(:,:) * tmask(:,:,1) 
    535535      ! 
     
    659659            zrv = pqa(ji,jj) / (1. - pqa(ji,jj)) 
    660660            ziRT = 1. / (R_dry*ptak(ji,jj))    ! 1/RT 
    661             gamma_moist(ji,jj) = grav * ( 1. + cevap*zrv*ziRT ) / ( Cp_dry + cevap*cevap*zrv*reps0*ziRT/ptak(ji,jj) ) 
     661            gamma_moist(ji,jj) = grav * ( 1. + rLevap*zrv*ziRT ) / ( Cp_dry + rLevap*rLevap*zrv*reps0*ziRT/ptak(ji,jj) ) 
    662662         END DO 
    663663      END DO 
     
    792792      REAL(wp) ::   zst3                     ! local variable 
    793793      REAL(wp) ::   zcoef_dqlw, zcoef_dqla   !   -      - 
    794       REAL(wp) ::   zztmp, z1_lsub           !   -      - 
     794      REAL(wp) ::   zztmp, z1_rLsub           !   -      - 
    795795      REAL(wp) ::   zfr1, zfr2               ! local variables 
    796796      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_st         ! inverse of surface temperature 
     
    868868 
    869869      ! --- evaporation --- ! 
    870       z1_lsub = 1._wp / Lsub 
    871       evap_ice (:,:,:) = rn_efac * qla_ice (:,:,:) * z1_lsub    ! sublimation 
    872       devap_ice(:,:,:) = rn_efac * dqla_ice(:,:,:) * z1_lsub    ! d(sublimation)/dT 
    873       zevap    (:,:)   = rn_efac * ( emp(:,:) + tprecip(:,:) )  ! evaporation over ocean 
     870      z1_rLsub = 1._wp / rLsub 
     871      evap_ice (:,:,:) = rn_efac * qla_ice (:,:,:) * z1_rLsub    ! sublimation 
     872      devap_ice(:,:,:) = rn_efac * dqla_ice(:,:,:) * z1_rLsub    ! d(sublimation)/dT 
     873      zevap    (:,:)   = rn_efac * ( emp(:,:) + tprecip(:,:) )   ! evaporation over ocean 
    874874 
    875875      ! --- evaporation minus precipitation --- ! 
     
    884884         &          + ( tprecip(:,:) - sprecip(:,:) ) * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp  & ! liquid precip at Tair 
    885885         &          +   sprecip(:,:) * ( 1._wp - zsnw ) *                                        & ! solid precip at min(Tair,Tsnow) 
    886          &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * cpic * tmask(:,:,1) - lfus ) 
     886         &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus ) 
    887887      qemp_ice(:,:) =   sprecip(:,:) * zsnw *                                                    & ! solid precip (only) 
    888          &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * cpic * tmask(:,:,1) - lfus ) 
     888         &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus ) 
    889889 
    890890      ! --- total solar and non solar fluxes --- ! 
     
    894894 
    895895      ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- ! 
    896       qprec_ice(:,:) = rhosn * ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * cpic * tmask(:,:,1) - lfus ) 
     896      qprec_ice(:,:) = rhos * ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus ) 
    897897 
    898898      ! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) --- 
    899899      DO jl = 1, jpl 
    900          qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * cpic * tmask(:,:,1) ) 
     900         qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * rcpi * tmask(:,:,1) ) 
    901901         !                         ! But we do not have Tice => consider it at 0degC => evap=0  
    902902      END DO 
     
    971971      CASE ( 1 , 2 ) 
    972972         ! 
    973          zfac  = 1._wp /  ( rn_cnd_s + rcdic ) 
     973         zfac  = 1._wp /  ( rn_cnd_s + rcnd_i ) 
    974974         zfac2 = EXP(1._wp) * 0.5_wp * zepsilon 
    975975         zfac3 = 2._wp / zepsilon 
     
    978978            DO jj = 1 , jpj 
    979979               DO ji = 1, jpi 
    980                   zhe = ( rn_cnd_s * phi(ji,jj,jl) + rcdic * phs(ji,jj,jl) ) * zfac                             ! Effective thickness 
     980                  zhe = ( rn_cnd_s * phi(ji,jj,jl) + rcnd_i * phs(ji,jj,jl) ) * zfac                            ! Effective thickness 
    981981                  IF( zhe >=  zfac2 )   zgfac(ji,jj,jl) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( zhe * zfac3 ) ) ) ! Enhanced conduction factor 
    982982               END DO 
     
    990990      ! -------------------------------------------------------------! 
    991991      ! 
    992       zfac = rcdic * rn_cnd_s 
     992      zfac = rcnd_i * rn_cnd_s 
    993993      ! 
    994994      DO jl = 1, jpl 
     
    997997               !                     
    998998               zkeff_h = zfac * zgfac(ji,jj,jl) / &                                    ! Effective conductivity of the snow-ice system divided by thickness 
    999                   &      ( rcdic * phs(ji,jj,jl) + rn_cnd_s * MAX( 0.01, phi(ji,jj,jl) ) ) 
     999                  &      ( rcnd_i * phs(ji,jj,jl) + rn_cnd_s * MAX( 0.01, phi(ji,jj,jl) ) ) 
    10001000               ztsu    = ptsu(ji,jj,jl)                                                ! Store current iteration temperature 
    10011001               ztsu0   = ptsu(ji,jj,jl)                                                ! Store initial surface temperature 
  • NEMO/trunk/src/OCE/SBC/sbccpl.F90

    r9921 r9935  
    14181418            zqns(:,:) =  zqns(:,:) - zemp(:,:) * sst_m(:,:) * rcp         ! remove heat content due to mass flux (assumed to be at SST) 
    14191419            IF( srcv(jpr_snow  )%laction ) THEN 
    1420                zqns(:,:) = zqns(:,:) - frcv(jpr_snow)%z3(:,:,1) * lfus    ! energy for melting solid precipitation over the free ocean 
     1420               zqns(:,:) = zqns(:,:) - frcv(jpr_snow)%z3(:,:,1) * rLfus   ! energy for melting solid precipitation over the free ocean 
    14211421            ENDIF 
    14221422         ENDIF 
    14231423         ! 
    1424          IF( srcv(jpr_icb)%laction )  zqns(:,:) = zqns(:,:) - frcv(jpr_icb)%z3(:,:,1) * lfus ! remove heat content associated to iceberg melting 
     1424         IF( srcv(jpr_icb)%laction )  zqns(:,:) = zqns(:,:) - frcv(jpr_icb)%z3(:,:,1) * rLfus ! remove heat content associated to iceberg melting 
    14251425         ! 
    14261426         IF( ln_mixcpl ) THEN   ;   qns(:,:) = qns(:,:) * xcplmask(:,:,0) + zqns(:,:) * zmsk(:,:) 
     
    18111811      !                                      
    18121812      ! --- calving (removed from qns_tot) --- ! 
    1813       IF( srcv(jpr_cal)%laction )   zqns_tot(:,:) = zqns_tot(:,:) - frcv(jpr_cal)%z3(:,:,1) * lfus  ! remove latent heat of calving 
    1814                                                                                                     ! we suppose it melts at 0deg, though it should be temp. of surrounding ocean 
     1813      IF( srcv(jpr_cal)%laction )   zqns_tot(:,:) = zqns_tot(:,:) - frcv(jpr_cal)%z3(:,:,1) * rLfus  ! remove latent heat of calving 
     1814                                                                                                     ! we suppose it melts at 0deg, though it should be temp. of surrounding ocean 
    18151815      ! --- iceberg (removed from qns_tot) --- ! 
    1816       IF( srcv(jpr_icb)%laction )   zqns_tot(:,:) = zqns_tot(:,:) - frcv(jpr_icb)%z3(:,:,1) * lfus  ! remove latent heat of iceberg melting 
     1816      IF( srcv(jpr_icb)%laction )   zqns_tot(:,:) = zqns_tot(:,:) - frcv(jpr_icb)%z3(:,:,1) * rLfus  ! remove latent heat of iceberg melting 
    18171817 
    18181818#if defined key_si3       
     
    18231823 
    18241824      ! Heat content per unit mass of snow (J/kg) 
    1825       WHERE( SUM( a_i, dim=3 ) > 1.e-10 )   ;   zcptsnw(:,:) = cpic * SUM( (tn_ice - rt0) * a_i, dim=3 ) / SUM( a_i, dim=3 ) 
     1825      WHERE( SUM( a_i, dim=3 ) > 1.e-10 )   ;   zcptsnw(:,:) = rcpi * SUM( (tn_ice - rt0) * a_i, dim=3 ) / SUM( a_i, dim=3 ) 
    18261826      ELSEWHERE                             ;   zcptsnw(:,:) = zcptn(:,:) 
    18271827      ENDWHERE 
     
    18301830 
    18311831      ! --- enthalpy of snow precip over ice in J/m3 (to be used in 1D-thermo) --- ! 
    1832       zqprec_ice(:,:) = rhosn * ( zcptsnw(:,:) - lfus ) 
     1832      zqprec_ice(:,:) = rhos * ( zcptsnw(:,:) - rLfus ) 
    18331833 
    18341834      ! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) --- ! 
    18351835      DO jl = 1, jpl 
    1836          zqevap_ice(:,:,jl) = 0._wp ! should be -evap * ( ( Tice - rt0 ) * cpic ) but atm. does not take it into account 
     1836         zqevap_ice(:,:,jl) = 0._wp ! should be -evap * ( ( Tice - rt0 ) * rcpi ) but atm. does not take it into account 
    18371837      END DO 
    18381838 
     
    18401840      zqemp_oce(:,:) = -  zevap_oce(:,:)                                      *   zcptn   (:,:)   &        ! evap 
    18411841         &             + ( ztprecip(:,:) - zsprecip(:,:) )                    *   zcptrain(:,:)   &        ! liquid precip 
    1842          &             +   zsprecip(:,:)                   * ( 1._wp - zsnw ) * ( zcptsnw (:,:) - lfus )   ! solid precip over ocean + snow melting 
    1843       zqemp_ice(:,:) =     zsprecip(:,:)                   * zsnw             * ( zcptsnw (:,:) - lfus )   ! solid precip over ice (qevap_ice=0 since atm. does not take it into account) 
     1842         &             +   zsprecip(:,:)                   * ( 1._wp - zsnw ) * ( zcptsnw (:,:) - rLfus )  ! solid precip over ocean + snow melting 
     1843      zqemp_ice(:,:) =     zsprecip(:,:)                   * zsnw             * ( zcptsnw (:,:) - rLfus )  ! solid precip over ice (qevap_ice=0 since atm. does not take it into account) 
    18441844!!    zqemp_ice(:,:) = -   frcv(jpr_ievp)%z3(:,:,1)        * picefr(:,:)      *   zcptsnw (:,:)   &        ! ice evap 
    1845 !!       &             +   zsprecip(:,:)                   * zsnw             * zqprec_ice(:,:) * r1_rhosn ! solid precip over ice 
     1845!!       &             +   zsprecip(:,:)                   * zsnw             * zqprec_ice(:,:) * r1_rhos  ! solid precip over ice 
    18461846       
    18471847      ! --- total non solar flux (including evap/precip) --- ! 
     
    18741874       
    18751875      ! clem: this formulation is certainly wrong... but better than it was... 
    1876       zqns_tot(:,:) = zqns_tot(:,:)                            &          ! zqns_tot update over free ocean with: 
    1877          &          - (  ziceld(:,:) * zsprecip(:,:) * lfus )  &          ! remove the latent heat flux of solid precip. melting 
    1878          &          - (  zemp_tot(:,:)                         &          ! remove the heat content of mass flux (assumed to be at SST) 
     1876      zqns_tot(:,:) = zqns_tot(:,:)                             &          ! zqns_tot update over free ocean with: 
     1877         &          - (  ziceld(:,:) * zsprecip(:,:) * rLfus )  &          ! remove the latent heat flux of solid precip. melting 
     1878         &          - (  zemp_tot(:,:)                          &          ! remove the heat content of mass flux (assumed to be at SST) 
    18791879         &             - zemp_ice(:,:) ) * zcptn(:,:)  
    18801880 
     
    18921892#endif 
    18931893      ! outputs 
    1894       IF ( srcv(jpr_cal)%laction       ) CALL iom_put('hflx_cal_cea'    , - frcv(jpr_cal)%z3(:,:,1) * lfus )                       ! latent heat from calving 
    1895       IF ( srcv(jpr_icb)%laction       ) CALL iom_put('hflx_icb_cea'    , - frcv(jpr_icb)%z3(:,:,1) * lfus )                       ! latent heat from icebergs melting 
     1894      IF ( srcv(jpr_cal)%laction       ) CALL iom_put('hflx_cal_cea'    , - frcv(jpr_cal)%z3(:,:,1) * rLfus )                      ! latent heat from calving 
     1895      IF ( srcv(jpr_icb)%laction       ) CALL iom_put('hflx_icb_cea'    , - frcv(jpr_icb)%z3(:,:,1) * rLfus )                      ! latent heat from icebergs melting 
    18961896      IF ( iom_use('hflx_rain_cea')    ) CALL iom_put('hflx_rain_cea'   , ( tprecip(:,:) - sprecip(:,:) ) * zcptrain(:,:) )        ! heat flux from rain (cell average) 
    18971897      IF ( iom_use('hflx_evap_cea')    ) CALL iom_put('hflx_evap_cea'   , ( frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) & 
    18981898           &                                                              * picefr(:,:) ) * zcptn(:,:) * tmask(:,:,1) )            ! heat flux from evap (cell average) 
    1899       IF ( iom_use('hflx_snow_cea')    ) CALL iom_put('hflx_snow_cea'   , sprecip(:,:) * ( zcptsnw(:,:) - Lfus )   )               ! heat flux from snow (cell average) 
    1900       IF ( iom_use('hflx_snow_ao_cea') ) CALL iom_put('hflx_snow_ao_cea', sprecip(:,:) * ( zcptsnw(:,:) - Lfus ) & 
     1899      IF ( iom_use('hflx_snow_cea')    ) CALL iom_put('hflx_snow_cea'   , sprecip(:,:) * ( zcptsnw(:,:) - rLfus )  )               ! heat flux from snow (cell average) 
     1900      IF ( iom_use('hflx_snow_ao_cea') ) CALL iom_put('hflx_snow_ao_cea', sprecip(:,:) * ( zcptsnw(:,:) - rLfus ) & 
    19011901           &                                                              * ( 1._wp - zsnw(:,:) )                  )               ! heat flux from snow (over ocean) 
    1902       IF ( iom_use('hflx_snow_ai_cea') ) CALL iom_put('hflx_snow_ai_cea', sprecip(:,:) * ( zcptsnw(:,:) - Lfus ) & 
     1902      IF ( iom_use('hflx_snow_ai_cea') ) CALL iom_put('hflx_snow_ai_cea', sprecip(:,:) * ( zcptsnw(:,:) - rLfus ) &  
    19031903           &                                                              *           zsnw(:,:)                    )               ! heat flux from snow (over ice) 
    19041904      ! note: hflx for runoff and iceshelf are done in sbcrnf and sbcisf resp. 
  • NEMO/trunk/src/OCE/SBC/sbcice_cice.F90

    r9927 r9935  
    1313   USE dom_oce         ! ocean space and time domain 
    1414   USE domvvl 
    15    USE phycst, only : rcp, rau0, r1_rau0, rhosn, rhoic 
     15   USE phycst, only : rcp, rau0, r1_rau0, rhos, rhoi 
    1616   USE in_out_manager  ! I/O manager 
    1717   USE iom, ONLY : iom_put,iom_use              ! I/O manager library !!Joakim edit 
     
    222222      CALL cice2nemo(vsno(:,:,:),ztmp1,'T', 1. ) 
    223223      CALL cice2nemo(vice(:,:,:),ztmp2,'T', 1. ) 
    224       snwice_mass  (:,:) = ( rhosn * ztmp1(:,:) + rhoic * ztmp2(:,:)  ) 
     224      snwice_mass  (:,:) = ( rhos * ztmp1(:,:) + rhoi * ztmp2(:,:)  ) 
    225225      snwice_mass_b(:,:) = snwice_mass(:,:) 
    226226 
     
    328328         ELSE 
    329329! emp_ice is set in sbc_cpl_ice_flx as sublimation-snow 
    330             qla_ice(:,:,1)= - ( emp_ice(:,:)+sprecip(:,:) ) * Lsub 
     330            qla_ice(:,:,1)= - ( emp_ice(:,:)+sprecip(:,:) ) * rLsub 
    331331! End of temporary code 
    332332            DO jj=1,jpj 
     
    644644      CALL cice2nemo(vsno(:,:,:),ztmp1,'T', 1. ) 
    645645      CALL cice2nemo(vice(:,:,:),ztmp2,'T', 1. ) 
    646       snwice_mass  (:,:) = ( rhosn * ztmp1(:,:) + rhoic * ztmp2(:,:)  ) 
     646      snwice_mass  (:,:) = ( rhos * ztmp1(:,:) + rhoi * ztmp2(:,:)  ) 
    647647      snwice_mass_b(:,:) = snwice_mass(:,:) 
    648648      snwice_fmass (:,:) = ( snwice_mass(:,:) - snwice_mass_b(:,:) ) / dt 
     
    801801      tprecip(:,:) = sf(jp_snow)%fnow(:,:,1)+sf(jp_rain)%fnow(:,:,1) 
    802802! May be better to do this conversion somewhere else 
    803       qla_ice(:,:,1) = -Lsub*sf(jp_sblm)%fnow(:,:,1) 
     803      qla_ice(:,:,1) = -rLsub*sf(jp_sblm)%fnow(:,:,1) 
    804804      topmelt(:,:,1) = sf(jp_top1)%fnow(:,:,1) 
    805805      topmelt(:,:,2) = sf(jp_top2)%fnow(:,:,1) 
  • NEMO/trunk/src/OCE/SBC/sbcisf.F90

    r9865 r9935  
    5252   LOGICAL, PUBLIC ::   l_isfcpl = .false.       !: isf recieved from oasis 
    5353 
    54    REAL(wp), PUBLIC, SAVE ::   rcpi     = 2000.0_wp     !: specific heat of ice shelf             [J/kg/K] 
     54   REAL(wp), PUBLIC, SAVE ::   rcpisf   = 2000.0_wp     !: specific heat of ice shelf             [J/kg/K] 
    5555   REAL(wp), PUBLIC, SAVE ::   rkappa   = 1.54e-6_wp    !: heat diffusivity through the ice-shelf [m2/s] 
    5656   REAL(wp), PUBLIC, SAVE ::   rhoisf   = 920.0_wp      !: volumic mass of ice shelf              [kg/m3] 
    5757   REAL(wp), PUBLIC, SAVE ::   tsurf    = -20.0_wp      !: air temperature on top of ice shelf    [C] 
    58    REAL(wp), PUBLIC, SAVE ::   rlfusisf = 0.334e6_wp    !: latent heat of fusion of ice shelf     [J/kg] 
     58   REAL(wp), PUBLIC, SAVE ::   rLfusisf = 0.334e6_wp    !: latent heat of fusion of ice shelf     [J/kg] 
    5959 
    6060!: Variable used in fldread to read the forcing file (nn_isf == 4 .OR. nn_isf == 3) 
     
    114114            ! compute fwf and heat flux 
    115115            IF( .NOT.l_isfcpl ) THEN    ;   CALL sbc_isf_cav (kt) 
    116             ELSE                        ;   qisf(:,:)  = fwfisf(:,:) * rlfusisf  ! heat        flux 
     116            ELSE                        ;   qisf(:,:)  = fwfisf(:,:) * rLfusisf  ! heat        flux 
    117117            ENDIF 
    118118            ! 
     
    127127               fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1)         ! fresh water flux from the isf (fwfisf <0 mean melting)  
    128128            ENDIF 
    129             qisf(:,:)   = fwfisf(:,:) * rlfusisf             ! heat flux 
     129            qisf(:,:)   = fwfisf(:,:) * rLfusisf             ! heat flux 
    130130            stbl(:,:)   = soce 
    131131            ! 
     
    137137               fwfisf(:,:) = -sf_fwfisf(1)%fnow(:,:,1)            ! fwf 
    138138            ENDIF 
    139             qisf(:,:)   = fwfisf(:,:) * rlfusisf               ! heat flux 
     139            qisf(:,:)   = fwfisf(:,:) * rLfusisf               ! heat flux 
    140140            stbl(:,:)   = soce 
    141141            ! 
     
    454454                           & * r1_e1e2t(ji,jj) * tmask(ji,jj,jk) 
    455455              
    456                fwfisf(ji,jj) = qisf(ji,jj) / rlfusisf          !fresh water flux kg/(m2s)                   
     456               fwfisf(ji,jj) = qisf(ji,jj) / rLfusisf          !fresh water flux kg/(m2s)                   
    457457               fwfisf(ji,jj) = fwfisf(ji,jj) * ( soce / stbl(ji,jj) ) 
    458458               !add to salinity trend 
     
    526526               DO ji = 1, jpi 
    527527                  zhtflx(ji,jj) =   zgammat(ji,jj)*rcp*rau0*(ttbl(ji,jj)-zfrz(ji,jj)) 
    528                   zfwflx(ji,jj) = - zhtflx(ji,jj)/rlfusisf 
     528                  zfwflx(ji,jj) = - zhtflx(ji,jj)/rLfusisf 
    529529               END DO 
    530530            END DO 
     
    544544                  ! compute coeficient to solve the 2nd order equation 
    545545                  zeps1 = rcp*rau0*zgammat(ji,jj) 
    546                   zeps2 = rlfusisf*rau0*zgammas(ji,jj) 
    547                   zeps3 = rhoisf*rcpi*rkappa/MAX(risfdep(ji,jj),zeps) 
     546                  zeps2 = rLfusisf*rau0*zgammas(ji,jj) 
     547                  zeps3 = rhoisf*rcpisf*rkappa/MAX(risfdep(ji,jj),zeps) 
    548548                  zeps4 = zlamb2+zlamb3*risfdep(ji,jj) 
    549549                  zeps6 = zeps4-ttbl(ji,jj) 
  • NEMO/trunk/src/OCE/SBC/sbcrnf.F90

    r9727 r9935  
    128128            END WHERE 
    129129            WHERE( sf_t_rnf(1)%fnow(:,:,1) == -222._wp )             ! where fwf comes from melting of ice shelves or iceberg 
    130                rnf_tsc(:,:,jp_tem) = ztfrz(:,:) * rnf(:,:) * r1_rau0 - rnf(:,:) * rlfusisf * r1_rau0_rcp 
     130               rnf_tsc(:,:,jp_tem) = ztfrz(:,:) * rnf(:,:) * r1_rau0 - rnf(:,:) * rLfusisf * r1_rau0_rcp 
    131131            END WHERE 
    132132         ELSE                                                        ! use SST as runoffs temperature 
Note: See TracChangeset for help on using the changeset viewer.