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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 9937 for NEMO – NEMO

Changeset 9937 for NEMO


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

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

Location:
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src
Files:
25 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icealb.F90

    r9604 r9937  
    148148               ! 
    149149               !                       !--- Snow-covered ice albedo (freezing, melting cases) 
    150                IF( pt_su(ji,jj,jl) < rt0_snow ) THEN 
     150               IF( pt_su(ji,jj,jl) < rt0 ) THEN 
    151151                  zalb_snw = rn_alb_sdry - ( rn_alb_sdry - zalb_ice ) * EXP( - ph_snw(ji,jj,jl) * z1_c3 ) 
    152152               ELSE 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icecor.F90

    r9923 r9937  
    8686      IF ( nn_icesal == 2 ) THEN    !  salinity must stay in bounds [Simin,Simax]        ! 
    8787      !                             !----------------------------------------------------- 
    88          zzc = rhoic * r1_Dt_ice 
     88         zzc = rhoi * r1_Dt_ice 
    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_Dt_ice 
    140140               !                 ! salt, volume 
    141                diag_sice(ji,jj) = SUM( sv_i(ji,jj,:) - sv_i_b(ji,jj,:) ) * rhoic * r1_Dt_ice 
    142                diag_vice(ji,jj) = SUM( v_i (ji,jj,:) - v_i_b (ji,jj,:) ) * rhoic * r1_Dt_ice 
    143                diag_vsnw(ji,jj) = SUM( v_s (ji,jj,:) - v_s_b (ji,jj,:) ) * rhosn * r1_Dt_ice 
     141               diag_sice(ji,jj) = SUM( sv_i(ji,jj,:) - sv_i_b(ji,jj,:) ) * rhoi * r1_Dt_ice 
     142               diag_vice(ji,jj) = SUM( v_i (ji,jj,:) - v_i_b (ji,jj,:) ) * rhoi * r1_Dt_ice 
     143               diag_vsnw(ji,jj) = SUM( v_s (ji,jj,:) - v_s_b (ji,jj,:) ) * rhos * r1_Dt_ice 
    144144            END DO 
    145145         END DO 
    … …  
    160160                  &                                   + SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) )  ) * r1_Dt_ice 
    161161               !                 ! salt, volume 
    162                diag_sice(ji,jj) = diag_sice(ji,jj) + SUM( sv_i(ji,jj,:) - sv_i_b(ji,jj,:) ) * rhoic * r1_Dt_ice 
    163                diag_vice(ji,jj) = diag_vice(ji,jj) + SUM( v_i (ji,jj,:) - v_i_b (ji,jj,:) ) * rhoic * r1_Dt_ice 
    164                diag_vsnw(ji,jj) = diag_vsnw(ji,jj) + SUM( v_s (ji,jj,:) - v_s_b (ji,jj,:) ) * rhosn * r1_Dt_ice 
     162               diag_sice(ji,jj) = diag_sice(ji,jj) + SUM( sv_i(ji,jj,:) - sv_i_b(ji,jj,:) ) * rhoi * r1_Dt_ice 
     163               diag_vice(ji,jj) = diag_vice(ji,jj) + SUM( v_i (ji,jj,:) - v_i_b (ji,jj,:) ) * rhoi * r1_Dt_ice 
     164               diag_vsnw(ji,jj) = diag_vsnw(ji,jj) + SUM( v_s (ji,jj,:) - v_s_b (ji,jj,:) ) * rhos * r1_Dt_ice 
    165165            END DO 
    166166         END DO 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icectl.F90

    r9923 r9937  
    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_Dt_ice - 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_Dt_ice + 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 ) 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icedia.F90

    r9923 r9937  
    110110      ! 3 -  Content variations ! 
    111111      ! ----------------------- ! 
    112       zdiff_vol = r1_rho0 * glob_sum( ( rhoic*vt_i(:,:) + rhosn*vt_s(:,:) - vol_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-9 ! freshwater trend (km3)  
    113       zdiff_sal = r1_rho0 * 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_rho0 * glob_sum( ( rhoi*vt_i(:,:) + rhos*vt_s(:,:) - vol_loc_ini(:,:) ) * e1e2t(:,:) ) * 1.e-9 ! freshwater trend (km3)  
     113      zdiff_sal = r1_rho0 * 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 
    … …  
    241241            frc_sal     = 0._wp                                                  
    242242            ! record initial ice volume, salt and temp 
    243             vol_loc_ini(:,:) = rhoic * vt_i(:,:) + rhosn * vt_s(:,:)  ! ice/snow volume (kg/m2) 
    244             tem_loc_ini(:,:) = et_i(:,:) + et_s(:,:)                  ! ice/snow heat content (J) 
    245             sal_loc_ini(:,:) = rhoic * SUM( sv_i(:,:,:), dim=3 )      ! ice salt content (pss*kg/m2) 
     243            vol_loc_ini(:,:) = rhoi * vt_i(:,:) + rhos * vt_s(:,:)   ! ice/snow volume (kg/m2) 
     244            tem_loc_ini(:,:) = et_i(:,:) + et_s(:,:)                 ! ice/snow heat content (J) 
     245            sal_loc_ini(:,:) = rhoi * SUM( sv_i(:,:,:), dim=3 )      ! ice salt content (pss*kg/m2) 
    246246         ENDIF 
    247247         ! 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icedyn_adv.F90

    r9923 r9937  
    103103      diag_trp_vi(:,:) = SUM(     v_i (:,:,:)          - v_i_b (:,:,:)                  , dim=3 ) * r1_Dt_ice 
    104104      diag_trp_vs(:,:) = SUM(     v_s (:,:,:)          - v_s_b (:,:,:)                  , dim=3 ) * r1_Dt_ice 
    105       IF( iom_use('icemtrp') )   CALL iom_put( "icemtrp" , diag_trp_vi * rhoic          )   ! ice mass transport 
    106       IF( iom_use('snwmtrp') )   CALL iom_put( "snwmtrp" , diag_trp_vs * rhosn          )   ! snw mass transport 
    107       IF( iom_use('salmtrp') )   CALL iom_put( "salmtrp" , diag_trp_sv * rhoic * 1.e-03 )   ! salt mass transport (kg/m2/s) 
     105      IF( iom_use('icemtrp') )   CALL iom_put( "icemtrp" , diag_trp_vi * rhoi           )   ! ice mass transport 
     106      IF( iom_use('snwmtrp') )   CALL iom_put( "snwmtrp" , diag_trp_vs * rhos           )   ! snw mass transport 
     107      IF( iom_use('salmtrp') )   CALL iom_put( "salmtrp" , diag_trp_sv * rhoi * 1.e-03 )   ! salt mass transport (kg/m2/s) 
    108108      IF( iom_use('dihctrp') )   CALL iom_put( "dihctrp" , -diag_trp_ei                 )   ! advected ice heat content (W/m2) 
    109109      IF( iom_use('dshctrp') )   CALL iom_put( "dshctrp" , -diag_trp_es                 )   ! advected snw heat content (W/m2) 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icedyn_rdgrft.F90

    r9923 r9937  
    543543               ! volume and enthalpy (J/m2, >0) of seawater trapped into ridges 
    544544               vsw = v_i_2d(ji,jl1) * afrdg * rn_porordg 
    545                ersw(ji) = -rhoic * vsw * rcp * sst_1d(ji)   ! clem: if sst>0, then ersw <0 (is that possible?) 
     545               ersw(ji) = -rhoi * vsw * rcp * sst_1d(ji)   ! clem: if sst>0, then ersw <0 (is that possible?) 
    546546 
    547547               ! volume etc of ridging / rafting ice and new ridges (vi, vs, sm, oi, es, ei) 
    … …  
    570570 
    571571               ! Ice-ocean exchanges associated with ice porosity 
    572                wfx_dyn_1d(ji) = wfx_dyn_1d(ji) - vsw              * rhoic * r1_Dt_ice   ! increase in ice volume due to seawater frozen in voids 
    573                sfx_dyn_1d(ji) = sfx_dyn_1d(ji) - vsw * sss_1d(ji) * rhoic * r1_Dt_ice 
    574                hfx_dyn_1d(ji) = hfx_dyn_1d(ji) + ersw(ji)                 * r1_Dt_ice   ! > 0 [W.m-2]  
     572               wfx_dyn_1d(ji) = wfx_dyn_1d(ji) - vsw              * rhoi * r1_Dt_ice   ! increase in ice volume due to seawater frozen in voids 
     573               sfx_dyn_1d(ji) = sfx_dyn_1d(ji) - vsw * sss_1d(ji) * rhoi * r1_Dt_ice 
     574               hfx_dyn_1d(ji) = hfx_dyn_1d(ji) + ersw(ji)                * r1_Dt_ice   ! > 0 [W.m-2]  
    575575 
    576576               ! Put the snow lost by ridging into the ocean 
    577577               !  Note that esrdg > 0; the ocean must cool to melt snow. If the ocean temp = Tf already, new ice must grow. 
    578                wfx_snw_dyn_1d(ji) = wfx_snw_dyn_1d(ji) + ( rhosn * vsrdg(ji) * ( 1._wp - rn_fsnwrdg )   &   ! fresh water source for ocean 
    579                   &                                      + rhosn * vsrft(ji) * ( 1._wp - rn_fsnwrft ) ) * r1_Dt_ice 
     578               wfx_snw_dyn_1d(ji) = wfx_snw_dyn_1d(ji) + ( rhos * vsrdg(ji) * ( 1._wp - rn_fsnwrdg )   &   ! fresh water source for ocean 
     579                  &                                      + rhos * vsrft(ji) * ( 1._wp - rn_fsnwrft ) ) * r1_Dt_ice 
    580580 
    581581               ! Put the melt pond water into the ocean 
    … …  
    583583               !       is no net mass flux between melt ponds and the ocean (see icethd_pnd.F90 for ex.) 
    584584               !IF ( ln_pnd_fwb ) THEN 
    585                !   wfx_pnd_1d(ji) = wfx_pnd_1d(ji) + ( rhofw * vprdg(ji) * ( 1._wp - rn_fpndrdg )   &        ! fresh water source for ocean 
    586                !      &                              + rhofw * vprft(ji) * ( 1._wp - rn_fpndrft ) ) * r1_Dt_ice 
     585               !   wfx_pnd_1d(ji) = wfx_pnd_1d(ji) + ( rhow * vprdg(ji) * ( 1._wp - rn_fpndrdg )   &        ! fresh water source for ocean 
     586               !      &                              + rhow * vprft(ji) * ( 1._wp - rn_fpndrft ) ) * r1_Dt_ice 
    587587               !ENDIF 
    588588 
    … …  
    590590               IF( nn_icesal /= 2 )  THEN 
    591591                  sirdg2(ji)     = sirdg2(ji)     - vsw * ( sss_1d(ji) - s_i_1d(ji) )        ! ridge salinity = s_i 
    592                   sfx_bri_1d(ji) = sfx_bri_1d(ji) + sss_1d(ji) * vsw * rhoic * r1_Dt_ice  &  ! put back sss_m into the ocean 
    593                      &                            - s_i_1d(ji) * vsw * rhoic * r1_Dt_ice     ! and get  s_i  from the ocean  
     592                  sfx_bri_1d(ji) = sfx_bri_1d(ji) + sss_1d(ji) * vsw * rhoi * r1_Dt_ice  &  ! put back sss_m into the ocean 
     593                     &                            - s_i_1d(ji) * vsw * rhoi * r1_Dt_ice     ! and get  s_i  from the ocean  
    594594               ENDIF 
    595595 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icedyn_rhg_evp.F90

    r9923 r9937  
    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) 
    … …  
    799799               zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * rswitch 
    800800                
    801                zdiag_xmtrp_ice(ji,jj) = rhoic * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component 
    802                zdiag_ymtrp_ice(ji,jj) = rhoic * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) !        ''           Y-   '' 
    803                 
    804                zdiag_xmtrp_snw(ji,jj) = rhosn * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component 
    805                zdiag_ymtrp_snw(ji,jj) = rhosn * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) !          ''          Y-   '' 
    806                 
    807                zdiag_xatrp(ji,jj)     = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) )         ! area transport,      X-component 
    808                zdiag_yatrp(ji,jj)     = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) )         !        ''            Y-   '' 
     801               zdiag_xmtrp_ice(ji,jj) = rhoi * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) )  ! ice mass transport, X-component 
     802               zdiag_ymtrp_ice(ji,jj) = rhoi * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) )  !        ''           Y-   '' 
     803                
     804               zdiag_xmtrp_snw(ji,jj) = rhos * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) )  ! snow mass transport, X-component 
     805               zdiag_ymtrp_snw(ji,jj) = rhos * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) )  !          ''          Y-   '' 
     806                
     807               zdiag_xatrp(ji,jj)     = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) )          ! area transport,      X-component 
     808               zdiag_yatrp(ji,jj)     = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) )          !        ''            Y-   '' 
    809809                
    810810            END DO 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/iceistate.F90

    r9923 r9937  
    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 - rho0 ) * h_i(ji,jj,jl) ) * r1_rho0 )  
     297                  zdh = MAX( 0._wp, ( rhos * h_s(ji,jj,jl) + ( rhoi - rho0 ) * h_i(ji,jj,jl) ) * r1_rho0 )  
    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 
    … …  
    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 )  )   & 
    344                         &             - rcp  * ( ztmelts - rt0 ) ) 
     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 )  )   & 
     344                        &             - rcp   * ( ztmelts - rt0 ) ) 
    345345                     ! 
    346346                     ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2 
    … …  
    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/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icethd.F90

    r9923 r9937  
    295295            ztmelts       = -tmut * sz_i_1d(ji,jk) 
    296296            ! Conversion q(S,T) -> T (second order equation) 
    297             zbbb          = ( rcp - cpic ) * ztmelts + e_i_1d(ji,jk) * r1_rhoic - lfus 
    298             zccc          = SQRT( MAX( zbbb * zbbb - 4._wp * cpic * lfus * ztmelts, 0._wp ) ) 
    299             t_i_1d(ji,jk) = rt0 - ( zbbb + zccc ) * 0.5_wp * r1_cpic 
     297            zbbb          = ( rcp - rcpi ) * ztmelts + e_i_1d(ji,jk) * r1_rhoi - rLfus 
     298            zccc          = SQRT( MAX( zbbb * zbbb - 4._wp * rcpi * rLfus * ztmelts, 0._wp ) ) 
     299            t_i_1d(ji,jk) = rt0 - ( zbbb + zccc ) * 0.5_wp * r1_cpi 
    300300             
    301301            ! mask temperature 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icethd_da.F90

    r9923 r9937  
    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_Dt_ice 
     139            sfx_lam_1d(ji) = sfx_lam_1d(ji) + rhoi *  h_i_1d(ji) * zda * s_i_1d(ji) * r1_Dt_ice 
    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_Dt_ice * ( rhoic * h_i_1d(ji) + rhosn * h_s_1d(ji) ) 
     146            wfx_lam_1d(ji) =  wfx_lam_1d(ji) + zda * r1_Dt_ice * ( rhoi * h_i_1d(ji) + rhos * h_s_1d(ji) ) 
    147147             
    148148            ! new concentration 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icethd_dh.F90

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

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icethd_do.F90

    r9923 r9937  
    140140         ! Physical constants 
    141141         zhicrit = 0.04                                          ! frazil ice thickness 
    142          ztwogp  = 2. * rho0 / ( grav * 0.3 * ( rho0 - rhoic ) ) ! reduced grav 
     142         ztwogp  = 2. * rho0 / ( grav * 0.3 * ( rho0 - rhoi ) ) ! reduced grav 
    143143         zsqcd   = 1.0 / SQRT( 1.3 * zcai )                      ! 1/SQRT(airdensity*drag) 
    144144         zgamafr = 0.03 
    … …  
    264264         DO ji = 1, npti 
    265265            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 ) 
     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_Dt_ice 
    294294            ! mass flux 
    295             wfx_opw_1d(ji) = wfx_opw_1d(ji) - zv_newice(ji) * rhoic * r1_Dt_ice 
     295            wfx_opw_1d(ji) = wfx_opw_1d(ji) - zv_newice(ji) * rhoi * r1_Dt_ice 
    296296            ! salt flux 
    297             sfx_opw_1d(ji) = sfx_opw_1d(ji) - zv_newice(ji) * rhoic * zs_newice(ji) * r1_Dt_ice 
     297            sfx_opw_1d(ji) = sfx_opw_1d(ji) - zv_newice(ji) * rhoi * zs_newice(ji) * r1_Dt_ice 
    298298         END DO 
    299299          

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icethd_pnd.F90

    r9923 r9937  
    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 
    136135      REAL(wp) ::   z1_zpnd_aspect   ! inverse pond aspect ratio 
    137136      REAL(wp) ::   zfac, zdum 
    … …  
    139138      INTEGER  ::   ji   ! loop indices 
    140139      !!------------------------------------------------------------------- 
    141       z1_rhofw       = 1._wp / rhofw  
    142140      z1_zpnd_aspect = 1._wp / zpnd_aspect 
    143141      z1_Tp          = 1._wp / zTp  
    … …  
    157155            ! 
    158156            ! 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) 
     157            zdv_mlt = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * r1_rhow * a_i_1d(ji) 
    160158            zfr_mlt = zrmin + ( zrmax - zrmin ) * a_i_1d(ji)  ! from CICE doc 
    161159            !zfr_mlt = zrmin + zrmax * a_i_1d(ji)             ! from Holland paper  
    … …  
    168166            ! melt pond mass flux (<0) 
    169167            IF( ln_pnd_fwb .AND. zdv_mlt > 0._wp ) THEN 
    170                zfac = zfr_mlt * zdv_mlt * rhofw * r1_Dt_ice 
     168               zfac = zfr_mlt * zdv_mlt * rhow * r1_Dt_ice 
    171169               wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac 
    172170               ! 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icethd_sal.F90

    r9923 r9937  
    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_Dt_ice 
     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_Dt_ice 
    101101            ENDIF 
    102102         END DO 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icethd_zdf_bl99.F90

    r9923 r9937  
    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!!gm faster coding 
     225                  ztcond_i(ji,jk) = rcnd_i + zbeta * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) )             & 
     226                     &                          / MIN( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) - 2._wp * rt0 , -epsi10 ) 
     227!!gm old 
     228!                  ztcond_i(ji,jk) = rcnd_i + zbeta * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) )   & 
     229!                     &                             / MIN( 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 , -epsi10 ) 
     230!!gm  
    226231               END DO 
    227232            END DO 
    … …  
    230235            ! 
    231236            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 ) 
     237               ztcond_i(ji,0)      = rcnd_i + 0.09_wp  *  sz_i_1d(ji,1)      / MIN( -epsi10, t_i_1d(ji,1) - rt0 )  & 
     238                  &                         - 0.011_wp * ( t_i_1d(ji,1) - rt0 ) 
     239               ztcond_i(ji,nlay_i) = rcnd_i + 0.09_wp  *  sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji)  - rt0 )  & 
     240                  &                         - 0.011_wp * ( t_bo_1d(ji) - rt0 ) 
    236241            END DO 
    237242            DO jk = 1, nlay_i-1 
    238243               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 ) 
     244!!gm faster coding 
     245                  zfac   = t_i_1d (ji,jk) + t_i_1d (ji,jk+1) - 2._wp * rt0 
     246                  ztcond_i(ji,jk) = rcnd_i + 0.09_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) / MIN( -epsi10, zfac )  & 
     247                     &                     - 0.0055_wp * zfac            !NB:  0.0055 = 1/2 * 0.011 
     248!!gm old 
     249!                  ztcond_i(ji,jk) = rcnd_i + 0.09_wp  *   0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) /        & 
     250!                     &                      MIN( -epsi10, 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 )  & 
     251!                     &                     - 0.011_wp * ( 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 ) 
     252!!gm  
    242253               END DO 
    243254            END DO 
    … …  
    299310         DO jk = 1, nlay_i 
    300311            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 )  
     312               zcpi = rcpi + zgamma * sz_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztiold(ji,jk) - rt0 ), epsi10 ) 
     313               zeta_i(ji,jk) = rdt_ice * r1_rhoi * z1_h_i(ji) / MAX( epsi10, zcpi )  
    303314            END DO 
    304315         END DO 
    … …  
    306317         DO jk = 1, nlay_s 
    307318            DO ji = 1, npti 
    308                zeta_s(ji,jk) = rdt_ice * r1_rhosn * r1_cpic * z1_h_s(ji) 
     319               zeta_s(ji,jk) = rdt_ice * r1_rhos * r1_cpi * z1_h_s(ji) 
    309320            END DO 
    310321         END DO 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/iceupdate.F90

    r9923 r9937  
    172172            snwice_mass_b(ji,jj) = snwice_mass(ji,jj)       ! save mass from the previous ice time step 
    173173            !                                               ! new mass per unit area 
    174             snwice_mass  (ji,jj) = tmask(ji,jj,1) * ( rhosn * vt_s(ji,jj) + rhoic * vt_i(ji,jj)  )  
     174            snwice_mass  (ji,jj) = tmask(ji,jj,1) * ( rhos * vt_s(ji,jj) + rhoi * vt_i(ji,jj)  )  
    175175            !                                               ! time evolution of snow+ice mass 
    176176            snwice_fmass (ji,jj) = ( snwice_mass(ji,jj) - snwice_mass_b(ji,jj) ) * r1_Dt_ice 
    … …  
    432432            ELSE                                     ! start from rest 
    433433               IF(lwp) WRITE(numout,*) '   ==>>   previous run without snow-ice mass output then set it' 
    434                snwice_mass  (:,:) = tmask(:,:,1) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:) ) 
     434               snwice_mass  (:,:) = tmask(:,:,1) * ( rhos * vt_s(:,:) + rhoi * vt_i(:,:) ) 
    435435               snwice_mass_b(:,:) = snwice_mass(:,:) 
    436436            ENDIF 
    437437         ELSE                                   !* Start from rest 
    438438            IF(lwp) WRITE(numout,*) '   ==>>   start from rest: set the snow-ice mass' 
    439             snwice_mass  (:,:) = tmask(:,:,1) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:) ) 
     439            snwice_mass  (:,:) = tmask(:,:,1) * ( rhos * vt_s(:,:) + rhoi * vt_i(:,:) ) 
    440440            snwice_mass_b(:,:) = snwice_mass(:,:) 
    441441         ENDIF 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icevar.F90

    r9923 r9937  
    228228                     ztmelts          = - sz_i(ji,jj,jk,jl) * tmut                                 ! 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_cpi , 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_cpi * ( -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_Dt_ice 
    501                wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_i (ji,jj,jl)   * rhoic * r1_Dt_ice 
    502                wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl)   * rhosn * r1_Dt_ice 
     500               sfx_res(ji,jj)  = sfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * sv_i(ji,jj,jl)   * rhoi * r1_Dt_ice 
     501               wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_i (ji,jj,jl)   * rhoi * r1_Dt_ice 
     502               wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl)   * rhos * r1_Dt_ice 
    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 - rho0 ) * zh_i(ji,jl) ) * r1_rho0 )  
     671               zdh = MAX( 0._wp, ( rhos * zh_s(ji,jl) + ( rhoi - rho0 ) * zh_i(ji,jl) ) * r1_rho0 )  
    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 
    … …  
    854854            ztmelts      = - tmut  * 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 
    856                                                                 !   (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 ) 
     856            !                                                   !   (sometimes zdf scheme produces abnormally high temperatures)    
     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/branches/2018/dev_r9838_ENHANCE04_MLF/src/ICE/icewri.F90

    r9923 r9937  
    8585      ! Standard outputs 
    8686      !----------------- 
    87       zrho1 = ( rho0 - rhoic ) * r1_rho0   ;   zrho2 = rhosn * r1_rho0 
     87      zrho1 = ( rho0 - rhoi ) * r1_rho0   ;   zrho2 = rhos * r1_rho0 
    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 
    … …  
    104104      IF( iom_use('snwvolu'  ) )   CALL iom_put( "snwvolu", vt_s  * zmsksn        )   ! snow volume 
    105105      IF( iom_use('icefrb') ) THEN 
     106!!gm remove the WHERE by using : 
     107!!         z2d(:,:) = MAX( zrho1 * hm_i(:,:) - zrho2 * hm_s(:,:) , 0._wp )                                          
     108!!gm end 
    106109         z2d(:,:) = ( zrho1 * hm_i(:,:) - zrho2 * hm_s(:,:) )                                          
    107110         WHERE( z2d < 0._wp )   z2d = 0._wp 
    … …  
    115118      ! salt 
    116119      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 
     120      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 
    118121 
    119122      ! heat 
    … …  
    164167      ! trends 
    165168      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) 
     169      IF( iom_use('dmidyn') )   CALL iom_put( "dmidyn", - wfx_dyn + rhoi * diag_trp_vi      )   ! Sea-ice mass change from dynamics(kg/m2/s) 
    167170      IF( iom_use('dmiopw') )   CALL iom_put( "dmiopw", - wfx_opw                           )   ! Sea-ice mass change through growth in open water 
    168171      IF( iom_use('dmibog') )   CALL iom_put( "dmibog", - wfx_bog                           )   ! Sea-ice mass change through basal growth 
    … …  
    174177      IF( iom_use('dmisub') )   CALL iom_put( "dmisub", - wfx_ice_sub                       )   ! Sea-ice mass change through sublimation 
    175178      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 
     179      IF( iom_use('dmsssi') )   CALL iom_put( "dmsssi",   wfx_sni*rhos*r1_rhoi              )   ! Snow mass change through snow-to-ice conversion 
    177180      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) 
     181      IF( iom_use('dmsdyn') )   CALL iom_put( "dmsdyn", - wfx_snw_dyn + rhos * diag_trp_vs )   ! Snow mass change through dynamics(kg/m2/s) 
    179182 
    180183      ! Global ice diagnostics 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/BDY/bdyice.F90

    r9923 r9937  
    124124 
    125125            ! Then, a) transfer the snow excess into the ice (different from icethd_dh) 
    126             zdh = MAX( 0._wp, ( rhosn * h_s(ji,jj,jl) + ( rhoic - rho0 ) * h_i(ji,jj,jl) ) * r1_rho0 ) 
     126            zdh = MAX( 0._wp, ( rhos * h_s(ji,jj,jl) + ( rhoi - rho0 ) * h_i(ji,jj,jl) ) * r1_rho0 ) 
    127127            ! Or, b) transfer all the snow into ice (if incoming ice is likely to melt as it comes into a warmer environment) 
    128             !zdh = MAX( 0._wp, h_s(ji,jj,jl) * rhosn / rhoic ) 
     128            !zdh = MAX( 0._wp, h_s(ji,jj,jl) * rhos / rhoi ) 
    129129 
    130130            ! recompute h_i, h_s 
    131131            h_i(ji,jj,jl) = MIN( hi_max(jl), h_i(ji,jj,jl) + zdh ) 
    132             h_s(ji,jj,jl) = MAX( 0._wp, h_s(ji,jj,jl) - zdh * rhoic / rhosn )  
    133  
    134          ENDDO 
     132            h_s(ji,jj,jl) = MAX( 0._wp, h_s(ji,jj,jl) - zdh * rhoi / rhos )  
     133 
     134         END DO 
    135135         CALL lbc_bdy_lnk( a_i(:,:,jl), 'T', 1., ib_bdy ) 
    136136         CALL lbc_bdy_lnk( h_i(:,:,jl), 'T', 1., ib_bdy ) 
    137137         CALL lbc_bdy_lnk( h_s(:,:,jl), 'T', 1., ib_bdy ) 
    138       ENDDO 
     138      END DO 
    139139      ! retrieve at_i 
    140140      at_i(:,:) = 0._wp 
    … …  
    212212            DO jk = 1, nlay_s 
    213213               ! Snow energy of melting 
    214                e_s(ji,jj,jk,jl) = rswitch * rhosn * ( cpic * ( rt0 - t_s(ji,jj,jk,jl) ) + lfus ) 
     214               e_s(ji,jj,jk,jl) = rswitch * rhos * ( rcpi * ( rt0 - t_s(ji,jj,jk,jl) ) + rLfus ) 
    215215               ! Multiply by volume, so that heat content in J/m2 
    216216               e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) * r1_nlay_s 
    … …  
    219219               ztmelts          = - tmut * sz_i(ji,jj,jk,jl) + rt0 !Melting temperature in K                   
    220220               ! heat content per unit volume 
    221                e_i(ji,jj,jk,jl) = rswitch * rhoic * & 
    222                   (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) & 
    223                   +   lfus    * ( 1.0 - (ztmelts-rt0) / MIN((t_i(ji,jj,jk,jl)-rt0),-epsi20) ) & 
    224                   - rcp      * ( ztmelts - rt0 ) ) 
     221               e_i(ji,jj,jk,jl) = rswitch * rhoi * & 
     222                  (   rcpi  * ( ztmelts - t_i(ji,jj,jk,jl) ) & 
     223                  +   rLfus * ( 1.0 - (ztmelts-rt0) / MIN((t_i(ji,jj,jk,jl)-rt0),-epsi20) ) & 
     224                  -   rcp   * ( ztmelts - rt0 ) ) 
    225225               ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2 
    226226               e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * a_i(ji,jj,jl) * h_i(ji,jj,jl) * r1_nlay_i 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DOM/phycst.F90

    r9923 r9937  
    3434   REAL(wp), PUBLIC ::   rhhmm =  60._wp        !: number of minutes in one hour 
    3535   REAL(wp), PUBLIC ::   rmmss =  60._wp        !: number of seconds in one minute 
    36    REAL(wp), PUBLIC ::   omega                  !: earth rotation parameter           [s-1] 
    37    REAL(wp), PUBLIC ::   ra    = 6371229._wp    !: earth radius                       [m] 
    38    REAL(wp), PUBLIC ::   grav  = 9.80665_wp     !: gravity                            [m/s2] 
     36   REAL(wp), PUBLIC ::   omega                  !: earth rotation parameter               [s-1] 
     37   REAL(wp), PUBLIC ::   ra    = 6371229._wp    !: earth radius                           [m] 
     38   REAL(wp), PUBLIC ::   grav  = 9.80665_wp     !: gravity                                [m/s2] 
    3939    
    40    REAL(wp), PUBLIC ::   rtt      = 273.16_wp        !: triple point of temperature   [Kelvin] 
    41    REAL(wp), PUBLIC ::   rt0      = 273.15_wp        !: freezing point of fresh water [Kelvin] 
    42    REAL(wp), PUBLIC ::   rt0_snow = 273.15_wp        !: melting point of snow         [Kelvin] 
    43 #if defined key_si3 
    44    REAL(wp), PUBLIC ::   rt0_ice  = 273.15_wp        !: melting point of ice          [Kelvin] 
    45 #else 
    46    REAL(wp), PUBLIC ::   rt0_ice  = 273.05_wp        !: melting point of ice          [Kelvin] 
    47 #endif 
    48    REAL(wp), PUBLIC ::   rho0                        !: volumic mass of reference     [kg/m3] 
    49    REAL(wp), PUBLIC ::   r1_rho0                     !: = 1. / rho0                   [m3/kg] 
    50    REAL(wp), PUBLIC ::   rcp                         !: ocean specific heat           [J/Kelvin] 
    51    REAL(wp), PUBLIC ::   r1_rcp                      !: = 1. / rcp                    [Kelvin/J] 
     40   REAL(wp), PUBLIC ::   rt0      = 273.15_wp        !: freezing point of fresh water     [Kelvin] 
     41   REAL(wp), PUBLIC ::   rho0                        !: volumic mass of reference         [kg/m3] 
     42   REAL(wp), PUBLIC ::   r1_rho0                     !: = 1. / rho0                       [m3/kg] 
     43   REAL(wp), PUBLIC ::   rcp                         !: ocean specific heat               [J/Kelvin] 
     44   REAL(wp), PUBLIC ::   r1_rcp                      !: = 1. / rcp                        [Kelvin/J] 
    5245   REAL(wp), PUBLIC ::   rho0_rcp                    !: = rho0 * rcp  
    5346   REAL(wp), PUBLIC ::   r1_rho0_rcp                 !: = 1. / ( rho0 * rcp ) 
    5447 
    55    REAL(wp), PUBLIC ::   rhosn    =  330._wp         !: volumic mass of snow          [kg/m3] 
    56    REAL(wp), PUBLIC ::   rhofw    = 1000._wp         !: volumic mass of freshwater in melt ponds [kg/m3] 
     48   REAL(wp), PUBLIC ::   rhoi     =  917._wp         !: sea ice density                   [kg/m3] 
     49   REAL(wp), PUBLIC ::   rhos     =  330._wp         !: snow    density                   [kg/m3] 
     50   REAL(wp), PUBLIC ::   rhow     = 1000._wp         !: water   density (in melt ponds)   [kg/m3] 
     51   REAL(wp), PUBLIC ::   rcnd_i   =    2.034396_wp   !: thermal conductivity of fresh ice [W/m/K] 
     52   REAL(wp), PUBLIC ::   rcpi     = 2067.0_wp        !: specific heat of fresh ice        [J/kg/K] 
     53   REAL(wp), PUBLIC ::   rLsub    =    2.834e+6_wp   !: pure ice latent heat of sublimation   [J/kg] 
     54   REAL(wp), PUBLIC ::   rLfus    =    0.334e+6_wp   !: latent heat of fusion of fresh ice    [J/kg] 
     55   REAL(wp), PUBLIC ::   tmut     =    0.054_wp      !: decrease of seawater meltpoint with salinity 
    5756   REAL(wp), PUBLIC ::   emic     =    0.97_wp       !: emissivity of snow or ice 
    58    REAL(wp), PUBLIC ::   sice     =    6.0_wp        !: salinity of ice               [psu] 
    59    REAL(wp), PUBLIC ::   soce     =   34.7_wp        !: salinity of sea               [psu] 
    60    REAL(wp), PUBLIC ::   cevap    =    2.5e+6_wp     !: latent heat of evaporation (water) 
    61    REAL(wp), PUBLIC ::   srgamma  =    0.9_wp        !: correction factor for solar radiation (Oberhuber, 1974) 
     57   REAL(wp), PUBLIC ::   sice     =    6.0_wp        !: salinity of ice                   [psu] 
     58   REAL(wp), PUBLIC ::   soce     =   34.7_wp        !: salinity of sea                   [psu] 
    6259   REAL(wp), PUBLIC ::   vkarmn   =    0.4_wp        !: von Karman constant 
    6360   REAL(wp), PUBLIC ::   stefan   =    5.67e-8_wp    !: Stefan-Boltzmann constant  
    6461 
    65 #if defined key_si3 || defined key_cice 
    66    REAL(wp), PUBLIC ::   rhoic    =  917._wp         !: volumic mass of sea ice                               [kg/m3] 
    67    REAL(wp), PUBLIC ::   rcdic    =    2.034396_wp   !: thermal conductivity of fresh ice                     [W/m/K] 
    68    REAL(wp), PUBLIC ::   cpic     = 2067.0_wp        !: specific heat of fresh ice                            [J/kg/K] 
    69    REAL(wp), PUBLIC ::   lsub     =    2.834e+6_wp   !: pure ice latent heat of sublimation                   [J/kg] 
    70    REAL(wp), PUBLIC ::   lfus     =    0.334e+6_wp   !: latent heat of fusion of fresh ice                    [J/kg] 
    71    REAL(wp), PUBLIC ::   tmut     =    0.054_wp      !: decrease of seawater meltpoint with salinity 
    72    REAL(wp), PUBLIC ::   xlsn                        !: = lfus*rhosn (volumetric latent heat fusion of snow)  [J/m3] 
    73 #else 
    74    REAL(wp), PUBLIC ::   rhoic    =  900._wp         !: volumic mass of sea ice                               [kg/m3] 
    75    REAL(wp), PUBLIC ::   rcdic    =    2.034396_wp   !: conductivity of the ice                               [W/m/K] 
    76    REAL(wp), PUBLIC ::   rcpic    =    1.8837e+6_wp  !: volumetric specific heat for ice                      [J/m3/K] 
    77    REAL(wp), PUBLIC ::   cpic                        !: = rcpic / rhoic  (specific heat for ice)              [J/Kg/K] 
    78    REAL(wp), PUBLIC ::   rcdsn    =    0.22_wp       !: conductivity of the snow                              [W/m/K] 
    79    REAL(wp), PUBLIC ::   rcpsn    =    6.9069e+5_wp  !: volumetric specific heat for snow                     [J/m3/K] 
    80    REAL(wp), PUBLIC ::   xlsn     =  110.121e+6_wp   !: volumetric latent heat fusion of snow                 [J/m3] 
    81    REAL(wp), PUBLIC ::   lfus                        !: = xlsn / rhosn   (latent heat of fusion of fresh ice) [J/Kg] 
    82    REAL(wp), PUBLIC ::   xlic     =  300.33e+6_wp    !: volumetric latent heat fusion of ice                  [J/m3] 
    83    REAL(wp), PUBLIC ::   xsn      =    2.8e+6_wp     !: volumetric latent heat of sublimation of snow         [J/m3] 
    84 #endif 
    85 #if defined key_cice 
    86    REAL(wp), PUBLIC ::   rcdsn    =    0.31_wp       !: thermal conductivity of snow                          [W/m/K] 
    87 #endif 
    88 #if defined key_si3 
    89    REAL(wp), PUBLIC ::   r1_rhoic                    !: 1 / rhoic 
    90    REAL(wp), PUBLIC ::   r1_rhosn                    !: 1 / rhosn 
    91    REAL(wp), PUBLIC ::   r1_cpic                     !: 1 / cpic 
    92 #endif 
     62   REAL(wp), PUBLIC ::   r1_rhoi                     !: 1 / rhoi 
     63   REAL(wp), PUBLIC ::   r1_rhos                     !: 1 / rhos 
     64   REAL(wp), PUBLIC ::   r1_rhow                     !: 1 / rhow 
     65   REAL(wp), PUBLIC ::   r1_cpi                      !: 1 / rcpi 
     66   REAL(wp), PUBLIC ::   r1_Lsub                     !: 1 / rLsub 
     67   REAL(wp), PUBLIC ::   r1_Lfus                     !: 1 / rLfus 
     68 
    9369   !!---------------------------------------------------------------------- 
    9470   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    10581      !! ** Purpose :   set and print the constants 
    10682      !!---------------------------------------------------------------------- 
    107  
     83      ! 
    10884      IF(lwp) WRITE(numout,*) 
    10985      IF(lwp) WRITE(numout,*) 'phy_cst : initialization of ocean parameters and constants' 
    11086      IF(lwp) WRITE(numout,*) '~~~~~~~' 
    11187 
    112       ! Define & print constants 
    113       ! ------------------------ 
    114       IF(lwp) WRITE(numout,*) 
    115       IF(lwp) WRITE(numout,*) '   Constants' 
    116  
    117       IF(lwp) WRITE(numout,*) 
    118       IF(lwp) WRITE(numout,*) '      mathematical constant                 rpi = ', rpi 
     88      !                 !==  Define derived constant  ==! 
    11989 
    12090      rsiyea = 365.25_wp * rday * 2._wp * rpi / 6.283076_wp 
    … …  
    12595      omega  = 2._wp * rpi / rsiday  
    12696#endif 
    127       IF(lwp) WRITE(numout,*) 
    128       IF(lwp) WRITE(numout,*) '      day                                rday   = ', rday,   ' s' 
    129       IF(lwp) WRITE(numout,*) '      sideral year                       rsiyea = ', rsiyea, ' s' 
    130       IF(lwp) WRITE(numout,*) '      sideral day                        rsiday = ', rsiday, ' s' 
    131       IF(lwp) WRITE(numout,*) '      omega                              omega  = ', omega,  ' s^-1' 
    132       IF(lwp) WRITE(numout,*) 
    133       IF(lwp) WRITE(numout,*) '      nb of months per year               raamo = ', raamo, ' months' 
    134       IF(lwp) WRITE(numout,*) '      nb of hours per day                 rjjhh = ', rjjhh, ' hours' 
    135       IF(lwp) WRITE(numout,*) '      nb of minutes per hour              rhhmm = ', rhhmm, ' mn' 
    136       IF(lwp) WRITE(numout,*) '      nb of seconds per minute            rmmss = ', rmmss, ' s' 
    137       IF(lwp) WRITE(numout,*) 
    138       IF(lwp) WRITE(numout,*) '      earth radius                         ra   = ', ra, ' m' 
    139       IF(lwp) WRITE(numout,*) '      gravity                              grav = ', grav , ' m/s^2' 
    140       IF(lwp) WRITE(numout,*) 
    141       IF(lwp) WRITE(numout,*) '      triple point of temperature      rtt      = ', rtt     , ' K' 
    142       IF(lwp) WRITE(numout,*) '      freezing point of water          rt0      = ', rt0     , ' K' 
    143       IF(lwp) WRITE(numout,*) '      melting point of snow            rt0_snow = ', rt0_snow, ' K' 
    144       IF(lwp) WRITE(numout,*) '      melting point of ice             rt0_ice  = ', rt0_ice , ' K' 
    145       IF(lwp) WRITE(numout,*) 
    146       IF(lwp) WRITE(numout,*) '   reference density and heat capacity now defined in eosbn2.f90' 
    147                
    148 #if defined key_si3 || defined key_cice 
    149       xlsn = lfus * rhosn        ! volumetric latent heat fusion of snow [J/m3] 
    150 #else 
    151       cpic = rcpic / rhoic       ! specific heat for ice   [J/Kg/K] 
    152       lfus = xlsn / rhosn        ! latent heat of fusion of fresh ice 
    153 #endif 
    154 #if defined key_si3 
    155       r1_rhoic = 1._wp / rhoic 
    156       r1_rhosn = 1._wp / rhosn 
    157       r1_cpic  = 1._wp / cpic 
    158 #endif 
    159       IF(lwp) THEN 
     97 
     98      r1_rhoi = 1._wp / rhoi 
     99      r1_rhos = 1._wp / rhos 
     100      r1_cpi  = 1._wp / rcpi 
     101      r1_Lsub = 1._wp / rLsub 
     102      r1_Lfus = 1._wp / rLfus 
     103 
     104      IF(lwp) THEN      !==  print constants  ==! 
    160105         WRITE(numout,*) 
    161 #if defined key_cice 
    162          WRITE(numout,*) '      thermal conductivity of the snow          = ', rcdsn   , ' J/s/m/K' 
    163 #endif 
    164          WRITE(numout,*) '      thermal conductivity of pure ice          = ', rcdic   , ' J/s/m/K' 
    165          WRITE(numout,*) '      fresh ice specific heat                   = ', cpic    , ' J/kg/K' 
    166          WRITE(numout,*) '      latent heat of fusion of fresh ice / snow = ', lfus    , ' J/kg' 
    167 #if defined key_si3 || defined key_cice 
    168          WRITE(numout,*) '      latent heat of subl.  of fresh ice / snow = ', lsub    , ' J/kg' 
    169 #else 
    170          WRITE(numout,*) '      density times specific heat for snow      = ', rcpsn   , ' J/m^3/K'  
    171          WRITE(numout,*) '      density times specific heat for ice       = ', rcpic   , ' J/m^3/K' 
    172          WRITE(numout,*) '      volumetric latent heat fusion of sea ice  = ', xlic    , ' J/m'  
    173          WRITE(numout,*) '      latent heat of sublimation of snow        = ', xsn     , ' J/kg'  
    174 #endif 
    175          WRITE(numout,*) '      volumetric latent heat fusion of snow     = ', xlsn    , ' J/m^3'  
    176          WRITE(numout,*) '      density of sea ice                        = ', rhoic   , ' kg/m^3' 
    177          WRITE(numout,*) '      density of snow                           = ', rhosn   , ' kg/m^3' 
    178          WRITE(numout,*) '      density of freshwater (in melt ponds)     = ', rhofw   , ' kg/m^3' 
    179          WRITE(numout,*) '      emissivity of snow or ice                 = ', emic   
    180          WRITE(numout,*) '      salinity of ice                           = ', sice    , ' psu' 
    181          WRITE(numout,*) '      salinity of sea                           = ', soce    , ' psu' 
    182          WRITE(numout,*) '      latent heat of evaporation (water)        = ', cevap   , ' J/m^3'  
    183          WRITE(numout,*) '      correction factor for solar radiation     = ', srgamma  
    184          WRITE(numout,*) '      von Karman constant                       = ', vkarmn  
    185          WRITE(numout,*) '      Stefan-Boltzmann constant                 = ', stefan  , ' J/s/m^2/K^4' 
     106         WRITE(numout,*) '   Constants' 
    186107         WRITE(numout,*) 
    187          WRITE(numout,*) '      conversion: degre ==> radian          rad = ', rad 
     108         WRITE(numout,*) '      mathematical constant              rpi    = ', rpi 
     109         WRITE(numout,*) '      conversion: degre ==> radian       rad    = ', rad 
    188110         WRITE(numout,*) 
    189          WRITE(numout,*) '      smallest real computer value       rsmall = ', rsmall 
     111         WRITE(numout,*) '      day in seconds                     rday   = ', rday  , ' s' 
     112         WRITE(numout,*) '      sideral year                       rsiyea = ', rsiyea, ' s' 
     113         WRITE(numout,*) '      sideral day                        rsiday = ', rsiday, ' s' 
     114         WRITE(numout,*) '      omega = 2 pi / rsiday              omega  = ', omega , ' s^-1' 
     115         WRITE(numout,*) '      earth radius                       ra     = ', ra    , ' m' 
     116         WRITE(numout,*) '      gravity                            grav   = ', grav  , ' m/s^2' 
     117         WRITE(numout,*) 
     118         WRITE(numout,*) '      nb of months per year              raamo  = ', raamo, ' months' 
     119         WRITE(numout,*) '      nb of hours per day                rjjhh  = ', rjjhh, ' hours' 
     120         WRITE(numout,*) '      nb of minutes per hour             rhhmm  = ', rhhmm, ' mn' 
     121         WRITE(numout,*) '      nb of seconds per minute           rmmss  = ', rmmss, ' s' 
     122         WRITE(numout,*) 
     123         WRITE(numout,*) '   reference ocean density and heat capacity now defined in eosbn2.f90' 
     124         WRITE(numout,*) 
     125         WRITE(numout,*) '      freezing point of freshwater                rt0    = ', rt0   , ' K' 
     126         WRITE(numout,*) '      sea ice density                             rhoi   = ', rhoi  , ' kg/m^3' 
     127         WRITE(numout,*) '      snow    density                             rhos   = ', rhos  , ' kg/m^3' 
     128         WRITE(numout,*) '      freshwater density (in melt ponds)          rhow   = ', rhow  , ' kg/m^3' 
     129         WRITE(numout,*) '      thermal conductivity of pure ice            rcnd_i = ', rcnd_i, ' J/s/m/K' 
     130         WRITE(numout,*) '      fresh ice specific heat                     rcpi   = ', rcpi  , ' J/kg/K' 
     131         WRITE(numout,*) '      latent heat of fusion of fresh ice / snow   rLfus  = ', rLfus , ' J/kg' 
     132         WRITE(numout,*) '      latent heat of subl.  of fresh ice / snow   rLsub  = ', rLsub , ' J/kg' 
     133         WRITE(numout,*) '      emissivity of snow or ice                   emic   = ', emic   
     134         WRITE(numout,*) '      salinity of ice                             sice   = ', sice  , ' psu' 
     135         WRITE(numout,*) '      salinity of sea                             soce   = ', soce  , ' psu' 
     136         WRITE(numout,*) '      von Karman constant                         vkarmn = ', vkarmn  
     137         WRITE(numout,*) '      Stefan-Boltzmann constant                   stefan = ', stefan, ' J/s/m^2/K^4' 
     138         WRITE(numout,*) 
     139         WRITE(numout,*) 
     140         WRITE(numout,*) '      smallest real computer value                rsmall = ', rsmall 
    190141      ENDIF 
    191  
     142      ! 
    192143   END SUBROUTINE phy_cst 
    193144 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/SBC/sbcblk.F90

    r9923 r9937  
    504504      ENDIF 
    505505 
    506       zqla(:,:) = L_vap(zst(:,:)) * zevap(:,:)     ! Latent Heat flux 
     506      zqla(:,:) = L_vap( zst(:,:) ) * zevap(:,:)     ! Latent Heat flux 
    507507 
    508508 
    … …  
    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_snow ) - rt0 ) * cpic 
     533         &     * ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi 
    534534      qns(:,:) = qns(:,:) * tmask(:,:,1) 
    535535      ! 
    … …  
    643643      !! ** Purpose : Compute the moist adiabatic lapse-rate. 
    644644      !!     => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate 
    645       !!     => http://www.geog.ucsb.edu/~joel/g266_s10/lecture_notes/chapt03/oh10_3_01/oh10_3_01.html 
    646645      !! 
    647646      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
    … …  
    652651      ! 
    653652      INTEGER  ::   ji, jj         ! dummy loop indices 
    654       REAL(wp) :: zrv, ziRT        ! local scalar 
     653      REAL(wp) ::   zrv, ziRT      ! local scalar 
     654      REAL(wp) ::   zLv = 2.5e+6_wp   ! latent heat of vaporisation  
    655655      !!---------------------------------------------------------------------------------- 
    656656      ! 
    … …  
    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. + zLv*zrv*ziRT ) / ( Cp_dry + zLv*zLv*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                    !   -      - 
    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 
     870      evap_ice (:,:,:) = rn_efac * qla_ice (:,:,:) * r1_Lsub    ! sublimation 
     871      devap_ice(:,:,:) = rn_efac * dqla_ice(:,:,:) * r1_Lsub    ! d(sublimation)/dT 
    873872      zevap    (:,:)   = rn_efac * ( emp(:,:) + tprecip(:,:) )  ! evaporation over ocean 
    874873 
    … …  
    884883         &          + ( tprecip(:,:) - sprecip(:,:) ) * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp  & ! liquid precip at Tair 
    885884         &          +   sprecip(:,:) * ( 1._wp - zsnw ) *                                        & ! solid precip at min(Tair,Tsnow) 
    886          &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus ) 
     885         &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus ) 
    887886      qemp_ice(:,:) =   sprecip(:,:) * zsnw *                                                    & ! solid precip (only) 
    888          &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus ) 
     887         &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus ) 
    889888 
    890889      ! --- total solar and non solar fluxes --- ! 
    … …  
    894893 
    895894      ! --- 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_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus ) 
     895      qprec_ice(:,:) = rhos * ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus ) 
    897896 
    898897      ! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) --- 
    899898      DO jl = 1, jpl 
    900          qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * cpic * tmask(:,:,1) ) 
     899         qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * rcpi * tmask(:,:,1) ) 
    901900         !                         ! But we do not have Tice => consider it at 0degC => evap=0  
    902901      END DO 
    … …  
    971970      CASE ( 1 , 2 ) 
    972971         ! 
    973          zfac  = 1._wp /  ( rn_cnd_s + rcdic ) 
     972         zfac  = 1._wp /  ( rn_cnd_s + rcnd_i ) 
    974973         zfac2 = EXP(1._wp) * 0.5_wp * zepsilon 
    975974         zfac3 = 2._wp / zepsilon 
    … …  
    978977            DO jj = 1 , jpj 
    979978               DO ji = 1, jpi 
    980                   zhe = ( rn_cnd_s * phi(ji,jj,jl) + rcdic * phs(ji,jj,jl) ) * zfac                             ! Effective thickness 
    981                   IF( zhe >=  zfac2 )   zgfac(ji,jj,jl) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( zhe * zfac3 ) ) ) ! Enhanced conduction factor 
     979                  zhe = ( rn_cnd_s * phi(ji,jj,jl) + rcnd_i * phs(ji,jj,jl) ) * zfac                             ! Effective thickness 
     980                  IF( zhe >=  zfac2 )   zgfac(ji,jj,jl) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( zhe * zfac3 ) ) )  ! Enhanced conduction factor 
    982981               END DO 
    983982            END DO 
    … …  
    990989      ! -------------------------------------------------------------! 
    991990      ! 
    992       zfac = rcdic * rn_cnd_s 
     991      zfac = rcnd_i * rn_cnd_s 
    993992      ! 
    994993      DO jl = 1, jpl 
    995994         DO jj = 1 , jpj 
    996995            DO ji = 1, jpi 
    997                !                     
    998                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) ) ) 
     996               !                                                                       ! Effective conductivity of the snow-ice system divided by thickness 
     997               zkeff_h = zfac * zgfac(ji,jj,jl) / ( rcnd_i * phs(ji,jj,jl) + rn_cnd_s * MAX( 0.01, phi(ji,jj,jl) ) ) 
    1000998               ztsu    = ptsu(ji,jj,jl)                                                ! Store current iteration temperature 
    1001999               ztsu0   = ptsu(ji,jj,jl)                                                ! Store initial surface temperature 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/SBC/sbccpl.F90

    r9923 r9937  
    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(:,:) 
    1827       ENDWHERE 
     1827      END WHERE 
    18281828      ! Heat content per unit mass of rain (J/kg) 
    18291829      zcptrain(:,:) = rcp * ( SUM( (tn_ice(:,:,:) - rt0) * a_i(:,:,:), dim=3 ) + sst_m(:,:) * ziceld(:,:) )  
    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) --- ! 
    … …  
    18751875      ! clem: this formulation is certainly wrong... but better than it was... 
    18761876      zqns_tot(:,:) = zqns_tot(:,:)                            &          ! zqns_tot update over free ocean with: 
    1877          &          - (  ziceld(:,:) * zsprecip(:,:) * lfus )  &          ! remove the latent heat flux of solid precip. melting 
     1877         &          - (  ziceld(:,:) * zsprecip(:,:) * rLfus )  &         ! remove the latent heat flux of solid precip. melting 
    18781878         &          - (  zemp_tot(:,:)                         &          ! remove the heat content of mass flux (assumed to be at SST) 
    18791879         &             - zemp_ice(:,:) ) * zcptn(:,:)  
    … …  
    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/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/SBC/sbcice_cice.F90

    r9923 r9937  
    1313   USE dom_oce         ! ocean space and time domain 
    1414   USE domvvl 
    15    USE phycst   , ONLY : rcp, rho0, r1_rho0, rhosn, rhoic 
     15   USE phycst   , ONLY : rcp, rho0, r1_rho0, 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 
    … …  
    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 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/SBC/sbcisf.F90

    r9923 r9937  
    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)        , SAVE ::   rcp_isf  = 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 ::   rho_isf  = 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] 
    5958 
    6059!: Variable used in fldread to read the forcing file (nn_isf == 4 .OR. nn_isf == 3) 
    … …  
    114113            ! compute fwf and heat flux 
    115114            IF( .NOT.l_isfcpl ) THEN    ;   CALL sbc_isf_cav (kt) 
    116             ELSE                        ;   qisf(:,:)  = fwfisf(:,:) * rlfusisf  ! heat        flux 
     115            ELSE                        ;   qisf(:,:)  = fwfisf(:,:) * rLfus    ! heat flux 
    117116            ENDIF 
    118117            ! 
    … …  
    127126               fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1)         ! fresh water flux from the isf (fwfisf <0 mean melting)  
    128127            ENDIF 
    129             qisf(:,:)   = fwfisf(:,:) * rlfusisf             ! heat flux 
     128            qisf(:,:)   = fwfisf(:,:) * rLfus                   ! heat flux 
    130129            stbl(:,:)   = soce 
    131130            ! 
    … …  
    137136               fwfisf(:,:) = -sf_fwfisf(1)%fnow(:,:,1)            ! fwf 
    138137            ENDIF 
    139             qisf(:,:)   = fwfisf(:,:) * rlfusisf               ! heat flux 
     138            qisf(:,:)   = fwfisf(:,:) * rLfus                     ! heat flux 
    140139            stbl(:,:)   = soce 
    141140            ! 
    … …  
    308307      qisf    (:,:)    = 0._wp   ;   fwfisf  (:,:) = 0._wp 
    309308      risf_tsc(:,:,:)  = 0._wp   ;   fwfisf_b(:,:) = 0._wp 
    310       ! 
    311       ! define isf tbl tickness, top and bottom indice 
    312       SELECT CASE ( nn_isf ) 
     309 
     310      SELECT CASE ( nn_isf )      ! define isf tbl tickness, top and bottom indice 
     311      ! 
    313312      CASE ( 1 )  
    314313         IF(lwp) WRITE(numout,*) 
    … …  
    455454                           & * r1_e1e2t(ji,jj) * tmask(ji,jj,jk) 
    456455              
    457                fwfisf(ji,jj) = qisf(ji,jj) / rlfusisf          !fresh water flux kg/(m2s)                   
     456               fwfisf(ji,jj) = qisf(ji,jj) * r1_Lfus                        ! fresh water flux kg/(m2s)                   
    458457               fwfisf(ji,jj) = fwfisf(ji,jj) * ( soce / stbl(ji,jj) ) 
    459458               !add to salinity trend 
    … …  
    527526               DO ji = 1, jpi 
    528527                  zhtflx(ji,jj) =   zgammat(ji,jj)*rcp*rho0*(ttbl(ji,jj)-zfrz(ji,jj)) 
    529                   zfwflx(ji,jj) = - zhtflx(ji,jj)/rlfusisf 
     528                  zfwflx(ji,jj) = - zhtflx(ji,jj) * r1_Lfus 
    530529               END DO 
    531530            END DO 
    … …  
    545544                  ! compute coeficient to solve the 2nd order equation 
    546545                  zeps1 = rcp*rho0*zgammat(ji,jj) 
    547                   zeps2 = rlfusisf*rho0*zgammas(ji,jj) 
    548                   zeps3 = rho_isf*rcpi*rkappa/MAX(risfdep(ji,jj),zeps) 
     546                  zeps2 = rLfus*rho0*zgammas(ji,jj) 
     547                  zeps3 = rho_isf*rcp_isf*rkappa/MAX(risfdep(ji,jj),zeps) 
    549548                  zeps4 = zlamb2+zlamb3*risfdep(ji,jj) 
    550549                  zeps6 = zeps4-ttbl(ji,jj) 

  • NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/SBC/sbcrnf.F90

    r9923 r9937  
    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_rho0 - rnf(:,:) * rlfusisf * r1_rho0_rcp 
     130               rnf_tsc(:,:,jp_tem) = ztfrz(:,:) * rnf(:,:) * r1_rho0 - rnf(:,:) * rLfus * r1_rho0_rcp 
    131131            END WHERE 
    132132         ELSE                                                        ! use SST as runoffs temperature 
Note: See TracChangeset for help on using the changeset viewer.