Changeset 14037 for NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG/src/ICE/icethd_dh.F90
- Timestamp:
- 2020-12-03T12:20:38+01:00 (3 years ago)
- Location:
- NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/ICE/icethd_dh.F90 (modified) (18 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG
- Property svn:externals
-
old new 8 8 9 9 # SETTE 10 ^/utils/CI/sette @13292sette10 ^/utils/CI/sette_wave@13990 sette
-
- Property svn:externals
-
NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG/src/ICE/icethd_dh.F90
r13226 r14037 13 13 !!---------------------------------------------------------------------- 14 14 !! ice_thd_dh : vertical sea-ice growth and melt 15 !! ice_thd_snwblow : distribute snow fall between ice and ocean 16 !!---------------------------------------------------------------------- 15 !!---------------------------------------------------------------------- 17 16 USE dom_oce ! ocean space and time domain 18 17 USE phycst ! physical constants … … 20 19 USE ice1D ! sea-ice: thermodynamics variables 21 20 USE icethd_sal ! sea-ice: salinity profiles 21 USE icevar ! for CALL ice_var_snwblow 22 22 ! 23 23 USE in_out_manager ! I/O manager … … 29 29 30 30 PUBLIC ice_thd_dh ! called by ice_thd 31 PUBLIC ice_thd_snwblow ! called in sbcblk/sbccpl and here32 33 INTERFACE ice_thd_snwblow34 MODULE PROCEDURE ice_thd_snwblow_1d, ice_thd_snwblow_2d35 END INTERFACE36 31 37 32 !!---------------------------------------------------------------------- … … 59 54 !! - Bottom accretion/ablation 60 55 !! - Snow ice formation 56 !! 57 !! ** Note : h=max(0,h+dh) are often used to ensure positivity of h. 58 !! very small negative values can occur otherwise (e.g. -1.e-20) 61 59 !! 62 60 !! References : Bitz and Lipscomb, 1999, J. Geophys. Res. … … 84 82 REAL(wp) :: zfmdt ! exchange mass flux x time step (J/m2), >0 towards the ocean 85 83 86 REAL(wp), DIMENSION(jpij) :: zqprec ! energy of fallen snow (J.m-3)87 84 REAL(wp), DIMENSION(jpij) :: zq_top ! heat for surface ablation (J.m-2) 88 85 REAL(wp), DIMENSION(jpij) :: zq_bot ! heat for bottom ablation (J.m-2) … … 90 87 REAL(wp), DIMENSION(jpij) :: zf_tt ! Heat budget to determine melting or freezing(W.m-2) 91 88 REAL(wp), DIMENSION(jpij) :: zevap_rema ! remaining mass flux from sublimation (kg.m-2) 92 93 REAL(wp), DIMENSION(jpij) :: zdh_s_mel ! snow melt 94 REAL(wp), DIMENSION(jpij) :: zdh_s_pre ! snow precipitation 95 REAL(wp), DIMENSION(jpij) :: zdh_s_sub ! snow sublimation 96 97 REAL(wp), DIMENSION(jpij,nlay_s) :: zh_s ! snw layer thickness 98 REAL(wp), DIMENSION(jpij,nlay_i) :: zh_i ! ice layer thickness 99 REAL(wp), DIMENSION(jpij,nlay_i) :: zdeltah 100 INTEGER , DIMENSION(jpij,nlay_i) :: icount ! number of layers vanished by melting 101 89 REAL(wp), DIMENSION(jpij) :: zdeltah 102 90 REAL(wp), DIMENSION(jpij) :: zsnw ! distribution of snow after wind blowing 91 92 INTEGER , DIMENSION(jpij,nlay_i) :: icount ! number of layers vanishing by melting 93 REAL(wp), DIMENSION(jpij,0:nlay_i+1) :: zh_i ! ice layer thickness (m) 94 REAL(wp), DIMENSION(jpij,0:nlay_s ) :: zh_s ! snw layer thickness (m) 95 REAL(wp), DIMENSION(jpij,0:nlay_s ) :: ze_s ! snw layer enthalpy (J.m-3) 103 96 104 97 REAL(wp) :: zswitch_sal … … 113 106 END SELECT 114 107 115 ! initialize layer thicknesses and enthalpies 108 ! initialize ice layer thicknesses and enthalpies 109 eh_i_old(1:npti,0:nlay_i+1) = 0._wp 116 110 h_i_old (1:npti,0:nlay_i+1) = 0._wp 117 eh_i_old(1:npti,0:nlay_i+1) = 0._wp111 zh_i (1:npti,0:nlay_i+1) = 0._wp 118 112 DO jk = 1, nlay_i 119 113 DO ji = 1, npti 114 eh_i_old(ji,jk) = h_i_1d(ji) * r1_nlay_i * e_i_1d(ji,jk) 120 115 h_i_old (ji,jk) = h_i_1d(ji) * r1_nlay_i 121 eh_i_old(ji,jk) = e_i_1d(ji,jk) * h_i_old(ji,jk) 116 zh_i (ji,jk) = h_i_1d(ji) * r1_nlay_i 117 END DO 118 END DO 119 ! 120 ! initialize snw layer thicknesses and enthalpies 121 zh_s(1:npti,0) = 0._wp 122 ze_s(1:npti,0) = 0._wp 123 DO jk = 1, nlay_s 124 DO ji = 1, npti 125 zh_s(ji,jk) = h_s_1d(ji) * r1_nlay_s 126 ze_s(ji,jk) = e_s_1d(ji,jk) 122 127 END DO 123 128 END DO … … 144 149 ! 145 150 DO ji = 1, npti 146 zf_tt(ji) = qcn_ice_bot_1d(ji) + qsb_ice_bot_1d(ji) + fhld_1d(ji) 151 zf_tt(ji) = qcn_ice_bot_1d(ji) + qsb_ice_bot_1d(ji) + fhld_1d(ji) + qtr_ice_bot_1d(ji) * frq_m_1d(ji) 147 152 zq_bot(ji) = MAX( 0._wp, zf_tt(ji) * rDt_ice ) 148 END DO149 150 ! Ice and snow layer thicknesses151 !-------------------------------152 DO jk = 1, nlay_i153 DO ji = 1, npti154 zh_i(ji,jk) = h_i_1d(ji) * r1_nlay_i155 END DO156 END DO157 158 DO jk = 1, nlay_s159 DO ji = 1, npti160 zh_s(ji,jk) = h_s_1d(ji) * r1_nlay_s161 END DO162 153 END DO 163 154 … … 172 163 DO ji = 1, npti 173 164 IF( t_s_1d(ji,jk) > rt0 ) THEN 174 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], < 0175 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + rhos 165 hfx_res_1d (ji) = hfx_res_1d (ji) - ze_s(ji,jk) * zh_s(ji,jk) * a_i_1d(ji) * r1_Dt_ice ! heat flux to the ocean [W.m-2], < 0 166 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + rhos * zh_s(ji,jk) * a_i_1d(ji) * r1_Dt_ice ! mass flux 176 167 ! updates 177 dh_s_mlt(ji) = dh_s_mlt(ji) - zh_s(ji,jk)178 h_s_1d (ji) = h_s_1d(ji) - zh_s(ji,jk)168 dh_s_mlt(ji) = dh_s_mlt(ji) - zh_s(ji,jk) 169 h_s_1d (ji) = MAX( 0._wp, h_s_1d (ji) - zh_s(ji,jk) ) 179 170 zh_s (ji,jk) = 0._wp 180 e_s_1d (ji,jk) = 0._wp 181 t_s_1d (ji,jk) = rt0 171 ze_s (ji,jk) = 0._wp 182 172 END IF 183 173 END DO … … 186 176 ! Snow precipitation 187 177 !------------------- 188 CALL ice_thd_snwblow( 1.0_wp - at_i_1d(1:npti), zsnw(1:npti) ) ! snow distribution over ice after wind blowing 189 190 zdeltah(1:npti,:) = 0._wp 178 CALL ice_var_snwblow( 1._wp - at_i_1d(1:npti), zsnw(1:npti) ) ! snow distribution over ice after wind blowing 179 191 180 DO ji = 1, npti 192 181 IF( sprecip_1d(ji) > 0._wp ) THEN 182 zh_s(ji,0) = zsnw(ji) * sprecip_1d(ji) * rDt_ice * r1_rhos / at_i_1d(ji) ! thickness of precip 183 ze_s(ji,0) = MAX( 0._wp, - qprec_ice_1d(ji) ) ! enthalpy of the precip (>0, J.m-3) 193 184 ! 194 ! --- precipitation --- 195 zdh_s_pre (ji) = zsnw(ji) * sprecip_1d(ji) * rDt_ice * r1_rhos / at_i_1d(ji) ! thickness change 196 zqprec (ji) = - qprec_ice_1d(ji) ! enthalpy of the precip (>0, J.m-3) 185 hfx_spr_1d(ji) = hfx_spr_1d(ji) + ze_s(ji,0) * zh_s(ji,0) * a_i_1d(ji) * r1_Dt_ice ! heat flux from snow precip (>0, W.m-2) 186 wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhos * zh_s(ji,0) * a_i_1d(ji) * r1_Dt_ice ! mass flux, <0 197 187 ! 198 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)199 wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhos * a_i_1d(ji) * zdh_s_pre(ji) * r1_Dt_ice ! mass flux, <0200 201 ! --- melt of falling snow ---202 rswitch = MAX( 0._wp , SIGN( 1._wp , zqprec(ji) - epsi20 ) )203 zdeltah (ji,1) = - rswitch * zq_top(ji) / MAX( zqprec(ji) , epsi20 ) ! thickness change204 zdeltah (ji,1) = MAX( - zdh_s_pre(ji), zdeltah(ji,1) ) ! bound melting205 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)206 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), >0207 208 ! updates available heat + precipitations after melting209 dh_s_mlt (ji) = dh_s_mlt(ji) + zdeltah(ji,1)210 zq_top (ji) = MAX( 0._wp , zq_top (ji) + zdeltah(ji,1) * zqprec(ji) )211 zdh_s_pre(ji) = zdh_s_pre(ji) + zdeltah(ji,1)212 213 188 ! update thickness 214 h_s_1d(ji) = MAX( 0._wp , h_s_1d(ji) + zdh_s_pre(ji) ) 215 ! 216 ELSE 217 ! 218 zdh_s_pre(ji) = 0._wp 219 zqprec (ji) = 0._wp 220 ! 189 h_s_1d(ji) = h_s_1d(ji) + zh_s(ji,0) 221 190 ENDIF 222 END DO223 224 ! recalculate snow layers225 DO jk = 1, nlay_s226 DO ji = 1, npti227 zh_s(ji,jk) = h_s_1d(ji) * r1_nlay_s228 END DO229 191 END DO 230 192 231 193 ! Snow melting 232 194 ! ------------ 233 ! If heat still available (zq_top > 0), then melt more snow 234 zdeltah(1:npti,:) = 0._wp 235 zdh_s_mel(1:npti) = 0._wp 236 DO jk = 1, nlay_s 195 ! If heat still available (zq_top > 0) 196 ! then all snw precip has been melted and we need to melt more snow 197 DO jk = 0, nlay_s 237 198 DO ji = 1, npti 238 199 IF( zh_s(ji,jk) > 0._wp .AND. zq_top(ji) > 0._wp ) THEN 239 200 ! 240 rswitch = MAX( 0._wp, SIGN( 1._wp, e_s_1d(ji,jk) - epsi20 ) ) 241 zdeltah (ji,jk) = - rswitch * zq_top(ji) / MAX( e_s_1d(ji,jk), epsi20 ) ! thickness change 242 zdeltah (ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji,jk) ) ! bound melting 243 zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,jk) 244 245 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) 246 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) 201 rswitch = MAX( 0._wp , SIGN( 1._wp , ze_s(ji,jk) - epsi20 ) ) 202 zdum = - rswitch * zq_top(ji) / MAX( ze_s(ji,jk), epsi20 ) ! thickness change 203 zdum = MAX( zdum , - zh_s(ji,jk) ) ! bound melting 204 205 hfx_snw_1d (ji) = hfx_snw_1d (ji) - ze_s(ji,jk) * zdum * a_i_1d(ji) * r1_Dt_ice ! heat used to melt snow(W.m-2, >0) 206 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhos * zdum * a_i_1d(ji) * r1_Dt_ice ! snow melting only = water into the ocean 247 207 248 208 ! updates available heat + thickness 249 dh_s_mlt(ji) = dh_s_mlt(ji) + zdeltah(ji,jk) 250 zq_top (ji) = MAX( 0._wp , zq_top(ji) + zdeltah(ji,jk) * e_s_1d(ji,jk) ) 251 h_s_1d (ji) = MAX( 0._wp , h_s_1d(ji) + zdeltah(ji,jk) ) 252 zh_s (ji,jk) = MAX( 0._wp , zh_s(ji,jk) + zdeltah(ji,jk) ) 209 dh_s_mlt(ji) = dh_s_mlt(ji) + zdum 210 zq_top (ji) = MAX( 0._wp , zq_top (ji) + zdum * ze_s(ji,jk) ) 211 h_s_1d (ji) = MAX( 0._wp , h_s_1d (ji) + zdum ) 212 zh_s (ji,jk) = MAX( 0._wp , zh_s (ji,jk) + zdum ) 213 !!$ IF( zh_s(ji,jk) == 0._wp ) ze_s(ji,jk) = 0._wp 253 214 ! 254 215 ENDIF … … 260 221 ! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates 261 222 ! comment: not counted in mass/heat exchange in iceupdate.F90 since this is an exchange with atm. (not ocean) 262 zdeltah(1:npti,:) = 0._wp 223 zdeltah (1:npti) = 0._wp ! total snow thickness that sublimates, < 0 224 zevap_rema(1:npti) = 0._wp 263 225 DO ji = 1, npti 264 226 IF( evap_ice_1d(ji) > 0._wp ) THEN 227 zdeltah (ji) = MAX( - evap_ice_1d(ji) * r1_rhos * rDt_ice, - h_s_1d(ji) ) ! amount of snw that sublimates, < 0 228 zevap_rema(ji) = MAX( 0._wp, evap_ice_1d(ji) * rDt_ice + zdeltah(ji) * rhos ) ! remaining evap in kg.m-2 (used for ice sublimation later on) 229 ENDIF 230 END DO 231 232 DO jk = 0, nlay_s 233 DO ji = 1, npti 234 zdum = MAX( -zh_s(ji,jk), zdeltah(ji) ) ! snow layer thickness that sublimates, < 0 265 235 ! 266 zdh_s_sub (ji) = MAX( - h_s_1d(ji) , - evap_ice_1d(ji) * r1_rhos * rDt_ice ) 267 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) 268 zdeltah (ji,1) = MAX( zdh_s_sub(ji), - zdh_s_pre(ji) ) 269 270 hfx_sub_1d (ji) = hfx_sub_1d(ji) + & ! Heat flux by sublimation [W.m-2], < 0 (sublimate snow that had fallen, then pre-existing snow) 271 & ( zdeltah(ji,1) * zqprec(ji) + ( zdh_s_sub(ji) - zdeltah(ji,1) ) * e_s_1d(ji,1) ) & 272 & * a_i_1d(ji) * r1_Dt_ice 273 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 274 275 ! new snow thickness 276 h_s_1d(ji) = MAX( 0._wp , h_s_1d(ji) + zdh_s_sub(ji) ) 277 ! update precipitations after sublimation and correct sublimation 278 zdh_s_pre(ji) = zdh_s_pre(ji) + zdeltah(ji,1) 279 zdh_s_sub(ji) = zdh_s_sub(ji) - zdeltah(ji,1) 280 ! 281 ELSE 282 ! 283 zdh_s_sub (ji) = 0._wp 284 zevap_rema(ji) = 0._wp 285 ! 286 ENDIF 287 END DO 288 289 ! --- Update snow diags --- ! 290 DO ji = 1, npti 291 dh_s_tot(ji) = zdh_s_mel(ji) + zdh_s_pre(ji) + zdh_s_sub(ji) 292 END DO 293 294 ! Update temperature, energy 295 !--------------------------- 296 ! new temp and enthalpy of the snow (remaining snow precip + remaining pre-existing snow) 297 DO jk = 1, nlay_s 298 DO ji = 1,npti 299 rswitch = MAX( 0._wp , SIGN( 1._wp, h_s_1d(ji) - epsi20 ) ) 300 e_s_1d(ji,jk) = rswitch / MAX( h_s_1d(ji), epsi20 ) * & 301 & ( ( zdh_s_pre(ji) ) * zqprec(ji) + & 302 & ( h_s_1d(ji) - zdh_s_pre(ji) ) * rhos * ( rcpi * ( rt0 - t_s_1d(ji,jk) ) + rLfus ) ) 303 END DO 304 END DO 305 236 hfx_sub_1d (ji) = hfx_sub_1d (ji) + ze_s(ji,jk) * zdum * a_i_1d(ji) * r1_Dt_ice ! Heat flux of snw that sublimates [W.m-2], < 0 237 wfx_snw_sub_1d(ji) = wfx_snw_sub_1d(ji) - rhos * zdum * a_i_1d(ji) * r1_Dt_ice ! Mass flux by sublimation 238 239 ! update thickness 240 h_s_1d(ji) = MAX( 0._wp , h_s_1d(ji) + zdum ) 241 zh_s (ji,jk) = MAX( 0._wp , zh_s (ji,jk) + zdum ) 242 !!$ IF( zh_s(ji,jk) == 0._wp ) ze_s(ji,jk) = 0._wp 243 244 ! update sublimation left 245 zdeltah(ji) = MIN( zdeltah(ji) - zdum, 0._wp ) 246 END DO 247 END DO 248 249 ! 306 250 ! ! ============ ! 307 251 ! ! Ice ! … … 310 254 ! Surface ice melting 311 255 !-------------------- 312 zdeltah(1:npti,:) = 0._wp ! important313 256 DO jk = 1, nlay_i 314 257 DO ji = 1, npti … … 318 261 319 262 zEi = - e_i_1d(ji,jk) * r1_rhoi ! Specific enthalpy of layer k [J/kg, <0] 320 zdE = 0._wp ! Specific enthalpy difference (J/kg, <0) 321 ! set up at 0 since no energy is needed to melt water...(it is already melted) 322 zdeltah(ji,jk) = MIN( 0._wp , - zh_i(ji,jk) ) ! internal melting occurs when the internal temperature is above freezing 323 ! this should normally not happen, but sometimes, heat diffusion leads to this 324 zfmdt = - zdeltah(ji,jk) * rhoi ! Mass flux x time step > 0 325 326 dh_i_itm(ji) = dh_i_itm(ji) + zdeltah(ji,jk) ! Cumulate internal melting 327 328 zfmdt = - rhoi * zdeltah(ji,jk) ! Recompute mass flux [kg/m2, >0] 329 330 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 331 ! ice enthalpy zEi is "sent" to the ocean 332 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 333 ! using s_i_1d and not sz_i_1d(jk) is ok 334 wfx_res_1d(ji) = wfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice ! Mass flux 335 263 zdE = 0._wp ! Specific enthalpy difference (J/kg, <0) 264 ! set up at 0 since no energy is needed to melt water...(it is already melted) 265 zdum = MIN( 0._wp , - zh_i(ji,jk) ) ! internal melting occurs when the internal temperature is above freezing 266 ! this should normally not happen, but sometimes, heat diffusion leads to this 267 zfmdt = - zdum * rhoi ! Recompute mass flux [kg/m2, >0] 268 ! 269 dh_i_itm(ji) = dh_i_itm(ji) + zdum ! Cumulate internal melting 270 ! 271 hfx_res_1d(ji) = hfx_res_1d(ji) + zEi * zfmdt * a_i_1d(ji) * r1_Dt_ice ! Heat flux to the ocean [W.m-2], <0 272 ! ice enthalpy zEi is "sent" to the ocean 273 wfx_res_1d(ji) = wfx_res_1d(ji) - rhoi * zdum * a_i_1d(ji) * r1_Dt_ice ! Mass flux 274 sfx_res_1d(ji) = sfx_res_1d(ji) - rhoi * zdum * s_i_1d(ji) * a_i_1d(ji) * r1_Dt_ice ! Salt flux 275 ! using s_i_1d and not sz_i_1d(jk) is ok 336 276 ELSE !-- Surface melting 337 277 … … 342 282 zfmdt = - zq_top(ji) / zdE ! Mass flux to the ocean [kg/m2, >0] 343 283 344 zd eltah(ji,jk)= - zfmdt * r1_rhoi ! Melt of layer jk [m, <0]345 346 zd eltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) ) ! Melt of layer jk cannot exceed the layer thickness [m, <0]347 348 zq_top(ji) = MAX( 0._wp , zq_top(ji) - zdeltah(ji,jk)* rhoi * zdE ) ! update available heat349 350 dh_i_sum(ji) = dh_i_sum(ji) + zd eltah(ji,jk)! Cumulate surface melt351 352 zfmdt = - rhoi * zd eltah(ji,jk)! Recompute mass flux [kg/m2, >0]284 zdum = - zfmdt * r1_rhoi ! Melt of layer jk [m, <0] 285 286 zdum = MIN( 0._wp , MAX( zdum , - zh_i(ji,jk) ) ) ! Melt of layer jk cannot exceed the layer thickness [m, <0] 287 288 zq_top(ji) = MAX( 0._wp , zq_top(ji) - zdum * rhoi * zdE ) ! update available heat 289 290 dh_i_sum(ji) = dh_i_sum(ji) + zdum ! Cumulate surface melt 291 292 zfmdt = - rhoi * zdum ! Recompute mass flux [kg/m2, >0] 353 293 354 294 zQm = zfmdt * zEw ! Energy of the melt water sent to the ocean [J/m2, <0] 355 295 356 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 357 ! using s_i_1d and not sz_i_1d(jk) is ok) 358 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_Dt_ice ! Heat flux [W.m-2], < 0 359 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 360 ! 361 wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice ! Mass flux 362 296 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zEw * zfmdt * a_i_1d(ji) * r1_Dt_ice ! Heat flux [W.m-2], < 0 297 hfx_sum_1d(ji) = hfx_sum_1d(ji) - zdE * zfmdt * a_i_1d(ji) * r1_Dt_ice ! Heat flux used in this process [W.m-2], > 0 298 wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoi * zdum * a_i_1d(ji) * r1_Dt_ice ! Mass flux 299 sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoi * zdum * s_i_1d(ji) * a_i_1d(ji) * r1_Dt_ice ! Salt flux >0 300 ! using s_i_1d and not sz_i_1d(jk) is ok) 363 301 END IF 364 302 ! update thickness 303 zh_i(ji,jk) = MAX( 0._wp, zh_i(ji,jk) + zdum ) 304 h_i_1d(ji) = MAX( 0._wp, h_i_1d(ji) + zdum ) 305 ! 306 ! update heat content (J.m-2) and layer thickness 307 eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdum * e_i_1d(ji,jk) 308 h_i_old (ji,jk) = h_i_old (ji,jk) + zdum 309 ! 310 ! 365 311 ! Ice sublimation 366 312 ! --------------- 367 zdum = MAX( - ( zh_i(ji,jk) + zdeltah(ji,jk) ) , - zevap_rema(ji) * r1_rhoi ) 368 zdeltah (ji,jk) = zdeltah (ji,jk) + zdum 369 dh_i_sub(ji) = dh_i_sub(ji) + zdum 370 371 sfx_sub_1d(ji) = sfx_sub_1d(ji) - rhoi * a_i_1d(ji) * zdum * s_i_1d(ji) * r1_Dt_ice ! Salt flux >0 372 ! clem: flux is sent to the ocean for simplicity 373 ! but salt should remain in the ice except 374 ! if all ice is melted. => must be corrected 375 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 376 377 wfx_ice_sub_1d(ji) = wfx_ice_sub_1d(ji) - rhoi * a_i_1d(ji) * zdum * r1_Dt_ice ! Mass flux > 0 378 379 ! update remaining mass flux 380 zevap_rema(ji) = zevap_rema(ji) + zdum * rhoi 381 313 zdum = MAX( - zh_i(ji,jk) , - zevap_rema(ji) * r1_rhoi ) 314 ! 315 hfx_sub_1d(ji) = hfx_sub_1d(ji) + e_i_1d(ji,jk) * zdum * a_i_1d(ji) * r1_Dt_ice ! Heat flux [W.m-2], < 0 316 wfx_ice_sub_1d(ji) = wfx_ice_sub_1d(ji) - rhoi * zdum * a_i_1d(ji) * r1_Dt_ice ! Mass flux > 0 317 sfx_sub_1d(ji) = sfx_sub_1d(ji) - rhoi * zdum * s_i_1d(ji) * a_i_1d(ji) * r1_Dt_ice ! Salt flux >0 318 ! clem: flux is sent to the ocean for simplicity 319 ! but salt should remain in the ice except 320 ! if all ice is melted. => must be corrected 321 ! update remaining mass flux and thickness 322 zevap_rema(ji) = zevap_rema(ji) + zdum * rhoi 323 zh_i(ji,jk) = MAX( 0._wp, zh_i(ji,jk) + zdum ) 324 h_i_1d(ji) = MAX( 0._wp, h_i_1d(ji) + zdum ) 325 dh_i_sub(ji) = dh_i_sub(ji) + zdum 326 327 ! update heat content (J.m-2) and layer thickness 328 eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdum * e_i_1d(ji,jk) 329 h_i_old (ji,jk) = h_i_old (ji,jk) + zdum 330 382 331 ! record which layers have disappeared (for bottom melting) 383 332 ! => icount=0 : no layer has vanished 384 333 ! => icount=5 : 5 layers have vanished 385 rswitch = MAX( 0._wp , SIGN( 1._wp , - ( zh_i(ji,jk) + zdeltah(ji,jk)) ) )334 rswitch = MAX( 0._wp , SIGN( 1._wp , - zh_i(ji,jk) ) ) 386 335 icount(ji,jk) = NINT( rswitch ) 387 zh_i(ji,jk) = MAX( 0._wp , zh_i(ji,jk) + zdeltah(ji,jk) )388 336 389 ! update heat content (J.m-2) and layer thickness 390 eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk) 391 h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk) 392 END DO 393 END DO 394 395 ! update ice thickness 396 DO ji = 1, npti 397 h_i_1d(ji) = MAX( 0._wp , h_i_1d(ji) + dh_i_sum(ji) + dh_i_itm(ji) + dh_i_sub(ji) ) 398 END DO 399 337 END DO 338 END DO 339 400 340 ! remaining "potential" evap is sent to ocean 401 341 DO ji = 1, npti … … 435 375 & + zswi2 * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) ) , 0.5 ) 436 376 437 s_i_new(ji) = zswitch_sal * zfracs * sss_1d(ji) + ( 1. - zswitch_sal ) * s_i_1d(ji) ! New ice salinity438 439 ztmelts = - rTmlt * s_i_new(ji) ! New ice melting point (C)440 441 zt_i_new = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i)442 443 zEi = rcpi * ( zt_i_new - (ztmelts+rt0) ) & ! Specific enthalpy of forming ice (J/kg, <0)444 & - rLfus * ( 1.0 - ztmelts / ( MIN( zt_i_new - rt0, -epsi10 ) ) ) + rcp* ztmelts445 446 zEw = rcp * ( t_bo_1d(ji) - rt0 ) ! Specific enthalpy of seawater (J/kg, < 0)447 448 zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0)449 450 dh_i_bog(ji) = rDt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoi ) )377 s_i_new(ji) = zswitch_sal * zfracs * sss_1d(ji) + ( 1. - zswitch_sal ) * s_i_1d(ji) ! New ice salinity 378 379 ztmelts = - rTmlt * s_i_new(ji) ! New ice melting point (C) 380 381 zt_i_new = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i) 382 383 zEi = rcpi * ( zt_i_new - (ztmelts+rt0) ) & ! Specific enthalpy of forming ice (J/kg, <0) 384 & - rLfus * ( 1.0 - ztmelts / ( MIN( zt_i_new - rt0, -epsi10 ) ) ) + rcp * ztmelts 385 386 zEw = rcp * ( t_bo_1d(ji) - rt0 ) ! Specific enthalpy of seawater (J/kg, < 0) 387 388 zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0) 389 390 dh_i_bog(ji) = rDt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoi ) ) 451 391 452 392 END DO 453 393 ! Contribution to Energy and Salt Fluxes 454 zfmdt = - rhoi * dh_i_bog(ji)! Mass flux x time step (kg/m2, < 0)394 zfmdt = - rhoi * dh_i_bog(ji) ! Mass flux x time step (kg/m2, < 0) 455 395 456 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 457 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 458 459 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 460 461 wfx_bog_1d(ji) = wfx_bog_1d(ji) - rhoi * a_i_1d(ji) * dh_i_bog(ji) * r1_Dt_ice ! Mass flux, <0 396 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zEw * zfmdt * a_i_1d(ji) * r1_Dt_ice ! Heat flux to the ocean [W.m-2], >0 397 hfx_bog_1d(ji) = hfx_bog_1d(ji) - zdE * zfmdt * a_i_1d(ji) * r1_Dt_ice ! Heat flux used in this process [W.m-2], <0 398 wfx_bog_1d(ji) = wfx_bog_1d(ji) - rhoi * dh_i_bog(ji) * a_i_1d(ji) * r1_Dt_ice ! Mass flux, <0 399 sfx_bog_1d(ji) = sfx_bog_1d(ji) - rhoi * dh_i_bog(ji) * s_i_new(ji) * a_i_1d(ji) * r1_Dt_ice ! Salt flux, <0 400 401 ! update thickness 402 zh_i(ji,nlay_i+1) = zh_i(ji,nlay_i+1) + dh_i_bog(ji) 403 h_i_1d(ji) = h_i_1d(ji) + dh_i_bog(ji) 462 404 463 405 ! update heat content (J.m-2) and layer thickness … … 471 413 ! Ice Basal melt 472 414 !--------------- 473 zdeltah(1:npti,:) = 0._wp ! important474 415 DO jk = nlay_i, 1, -1 475 416 DO ji = 1, npti … … 480 421 IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN !-- Internal melting 481 422 482 zEi = - e_i_1d(ji,jk) * r1_rhoi ! Specific enthalpy of melting ice (J/kg, <0) 483 zdE = 0._wp ! Specific enthalpy difference (J/kg, <0) 484 ! set up at 0 since no energy is needed to melt water...(it is already melted) 485 zdeltah (ji,jk) = MIN( 0._wp , - zh_i(ji,jk) ) ! internal melting occurs when the internal temperature is above freezing 486 ! this should normally not happen, but sometimes, heat diffusion leads to this 487 488 dh_i_itm (ji) = dh_i_itm(ji) + zdeltah(ji,jk) 489 490 zfmdt = - zdeltah(ji,jk) * rhoi ! Mass flux x time step > 0 491 492 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 493 ! ice enthalpy zEi is "sent" to the ocean 494 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 495 ! using s_i_1d and not sz_i_1d(jk) is ok 496 wfx_res_1d(ji) = wfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice ! Mass flux 497 498 ! update heat content (J.m-2) and layer thickness 499 eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk) 500 h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk) 501 423 zEi = - e_i_1d(ji,jk) * r1_rhoi ! Specific enthalpy of melting ice (J/kg, <0) 424 zdE = 0._wp ! Specific enthalpy difference (J/kg, <0) 425 ! set up at 0 since no energy is needed to melt water...(it is already melted) 426 zdum = MIN( 0._wp , - zh_i(ji,jk) ) ! internal melting occurs when the internal temperature is above freezing 427 ! this should normally not happen, but sometimes, heat diffusion leads to this 428 dh_i_itm (ji) = dh_i_itm(ji) + zdum 429 ! 430 zfmdt = - zdum * rhoi ! Mass flux x time step > 0 431 ! 432 hfx_res_1d(ji) = hfx_res_1d(ji) + zEi * zfmdt * a_i_1d(ji) * r1_Dt_ice ! Heat flux to the ocean [W.m-2], <0 433 ! ice enthalpy zEi is "sent" to the ocean 434 wfx_res_1d(ji) = wfx_res_1d(ji) - rhoi * zdum * a_i_1d(ji) * r1_Dt_ice ! Mass flux 435 sfx_res_1d(ji) = sfx_res_1d(ji) - rhoi * zdum * s_i_1d(ji) * a_i_1d(ji) * r1_Dt_ice ! Salt flux 436 ! using s_i_1d and not sz_i_1d(jk) is ok 502 437 ELSE !-- Basal melting 503 438 504 zEi = - e_i_1d(ji,jk) * r1_rhoi! Specific enthalpy of melting ice (J/kg, <0)505 zEw = rcp * ztmelts! Specific enthalpy of meltwater (J/kg, <0)506 zdE = zEi - zEw! Specific enthalpy difference (J/kg, <0)507 508 zfmdt = - zq_bot(ji) / zdE! Mass flux x time step (kg/m2, >0)509 510 zd eltah(ji,jk) = - zfmdt * r1_rhoi! Gross thickness change511 512 zd eltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) ) ! bound thickness change439 zEi = - e_i_1d(ji,jk) * r1_rhoi ! Specific enthalpy of melting ice (J/kg, <0) 440 zEw = rcp * ztmelts ! Specific enthalpy of meltwater (J/kg, <0) 441 zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0) 442 443 zfmdt = - zq_bot(ji) / zdE ! Mass flux x time step (kg/m2, >0) 444 445 zdum = - zfmdt * r1_rhoi ! Gross thickness change 446 447 zdum = MIN( 0._wp , MAX( zdum, - zh_i(ji,jk) ) ) ! bound thickness change 513 448 514 zq_bot(ji) = MAX( 0._wp , zq_bot(ji) - zdeltah(ji,jk) * rhoi * zdE ) ! update available heat. MAX is necessary for roundup errors 515 516 dh_i_bom(ji) = dh_i_bom(ji) + zdeltah(ji,jk) ! Update basal melt 517 518 zfmdt = - zdeltah(ji,jk) * rhoi ! Mass flux x time step > 0 519 520 zQm = zfmdt * zEw ! Heat exchanged with ocean 521 522 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 523 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 524 525 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 526 ! using s_i_1d and not sz_i_1d(jk) is ok 527 528 wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice ! Mass flux 529 530 ! update heat content (J.m-2) and layer thickness 531 eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk) 532 h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk) 449 zq_bot(ji) = MAX( 0._wp , zq_bot(ji) - zdum * rhoi * zdE ) ! update available heat. MAX is necessary for roundup errors 450 451 dh_i_bom(ji) = dh_i_bom(ji) + zdum ! Update basal melt 452 453 zfmdt = - zdum * rhoi ! Mass flux x time step > 0 454 455 zQm = zfmdt * zEw ! Heat exchanged with ocean 456 457 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zEw * zfmdt * a_i_1d(ji) * r1_Dt_ice ! Heat flux to the ocean [W.m-2], <0 458 hfx_bom_1d(ji) = hfx_bom_1d(ji) - zdE * zfmdt * a_i_1d(ji) * r1_Dt_ice ! Heat used in this process [W.m-2], >0 459 wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoi * zdum * a_i_1d(ji) * r1_Dt_ice ! Mass flux 460 sfx_bom_1d(ji) = sfx_bom_1d(ji) - rhoi * zdum * s_i_1d(ji) * a_i_1d(ji) * r1_Dt_ice ! Salt flux 461 ! using s_i_1d and not sz_i_1d(jk) is ok 533 462 ENDIF 534 463 ! update thickness 464 zh_i(ji,jk) = MAX( 0._wp, zh_i(ji,jk) + zdum ) 465 h_i_1d(ji) = MAX( 0._wp, h_i_1d(ji) + zdum ) 466 ! 467 ! update heat content (J.m-2) and layer thickness 468 eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdum * e_i_1d(ji,jk) 469 h_i_old (ji,jk) = h_i_old (ji,jk) + zdum 535 470 ENDIF 536 471 END DO 537 472 END DO 538 473 539 ! Update temperature, energy 540 ! -------------------------- 541 DO ji = 1, npti 542 h_i_1d(ji) = MAX( 0._wp , h_i_1d(ji) + dh_i_bog(ji) + dh_i_bom(ji) ) 543 END DO 544 545 ! If heat still available then melt more snow 546 !------------------------------------------- 547 zdeltah(1:npti,:) = 0._wp ! important 548 DO ji = 1, npti 549 zq_rema (ji) = zq_top(ji) + zq_bot(ji) 550 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, - h_s_1d(ji) ) ) ! =1 if snow 551 rswitch = rswitch * MAX( 0._wp, SIGN( 1._wp, e_s_1d(ji,1) - epsi20 ) ) 552 zdeltah (ji,1) = - rswitch * zq_rema(ji) / MAX( e_s_1d(ji,1), epsi20 ) 553 zdeltah (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - h_s_1d(ji) ) ) ! bound melting 554 dh_s_tot(ji) = dh_s_tot(ji) + zdeltah(ji,1) 555 h_s_1d (ji) = h_s_1d (ji) + zdeltah(ji,1) 556 557 zq_rema(ji) = zq_rema(ji) + zdeltah(ji,1) * e_s_1d(ji,1) ! update available heat (J.m-2) 558 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) 559 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhos * a_i_1d(ji) * zdeltah(ji,1) * r1_Dt_ice ! Mass flux 560 dh_s_mlt(ji) = dh_s_mlt(ji) + zdeltah(ji,1) 561 ! 562 ! Remaining heat flux (W.m-2) is sent to the ocean heat budget 563 qt_oce_ai_1d(ji) = qt_oce_ai_1d(ji) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_Dt_ice 564 565 IF( ln_icectl .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji) 566 END DO 567 568 ! 474 ! Remove snow if ice has melted entirely 475 ! -------------------------------------- 476 DO jk = 0, nlay_s 477 DO ji = 1,npti 478 IF( h_i_1d(ji) == 0._wp ) THEN 479 ! mass & energy loss to the ocean 480 hfx_res_1d(ji) = hfx_res_1d(ji) - ze_s(ji,jk) * zh_s(ji,jk) * a_i_1d(ji) * r1_Dt_ice ! heat flux to the ocean [W.m-2], < 0 481 wfx_res_1d(ji) = wfx_res_1d(ji) + rhos * zh_s(ji,jk) * a_i_1d(ji) * r1_Dt_ice ! mass flux 482 483 ! update thickness and energy 484 h_s_1d(ji) = 0._wp 485 ze_s (ji,jk) = 0._wp 486 zh_s (ji,jk) = 0._wp 487 ENDIF 488 END DO 489 END DO 490 491 ! Snow load on ice 492 ! ----------------- 493 ! When snow load exceeds Archimede's limit and sst is positive, 494 ! snow-ice formation (next bloc) can lead to negative ice enthalpy. 495 ! Therefore we consider here that this excess of snow falls into the ocean 496 zdeltah(1:npti) = h_s_1d(1:npti) + h_i_1d(1:npti) * (rhoi-rho0) * r1_rhos 497 DO jk = 0, nlay_s 498 DO ji = 1, npti 499 IF( zdeltah(ji) > 0._wp .AND. sst_1d(ji) > 0._wp ) THEN 500 ! snow layer thickness that falls into the ocean 501 zdum = MIN( zdeltah(ji) , zh_s(ji,jk) ) 502 ! mass & energy loss to the ocean 503 hfx_res_1d(ji) = hfx_res_1d(ji) - ze_s(ji,jk) * zdum * a_i_1d(ji) * r1_Dt_ice ! heat flux to the ocean [W.m-2], < 0 504 wfx_res_1d(ji) = wfx_res_1d(ji) + rhos * zdum * a_i_1d(ji) * r1_Dt_ice ! mass flux 505 ! update thickness and energy 506 h_s_1d(ji) = MAX( 0._wp, h_s_1d(ji) - zdum ) 507 zh_s (ji,jk) = MAX( 0._wp, zh_s(ji,jk) - zdum ) 508 ! update snow thickness that still has to fall 509 zdeltah(ji) = MAX( 0._wp, zdeltah(ji) - zdum ) 510 ENDIF 511 END DO 512 END DO 513 569 514 ! Snow-Ice formation 570 515 ! ------------------ 571 ! When snow load exce sses Archimede's limit, snow-ice interface goes down under sea-level,572 ! flooding of seawater transforms snow into ice dh_snowice is positive for the ice516 ! When snow load exceeds Archimede's limit, snow-ice interface goes down under sea-level, 517 ! flooding of seawater transforms snow into ice. Thickness that is transformed is dh_snowice (positive for the ice) 573 518 z1_rho = 1._wp / ( rhos+rho0-rhoi ) 519 zdeltah(1:npti) = 0._wp 574 520 DO ji = 1, npti 575 521 ! 576 dh_snowice(ji) = MAX( 522 dh_snowice(ji) = MAX( 0._wp , ( rhos * h_s_1d(ji) + (rhoi-rho0) * h_i_1d(ji) ) * z1_rho ) 577 523 578 524 h_i_1d(ji) = h_i_1d(ji) + dh_snowice(ji) … … 584 530 zQm = zfmdt * zEw 585 531 586 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_Dt_ice ! Heat flux 587 588 sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_1d(ji) * a_i_1d(ji) * zfmdt * r1_Dt_ice ! Salt flux 532 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zEw * zfmdt * a_i_1d(ji) * r1_Dt_ice ! Heat flux 533 sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_1d(ji) * zfmdt * a_i_1d(ji) * r1_Dt_ice ! Salt flux 589 534 590 535 ! Case constant salinity in time: virtual salt flux to keep salinity constant 591 536 IF( nn_icesal /= 2 ) THEN 592 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 ocean593 & - s_i_1d(ji) * a_i_1d(ji) * dh_snowice(ji) * rhoi* r1_Dt_ice ! and get rn_icesal from the ocean537 sfx_bri_1d(ji) = sfx_bri_1d(ji) - sss_1d(ji) * zfmdt * a_i_1d(ji) * r1_Dt_ice & ! put back sss_m into the ocean 538 & - s_i_1d(ji) * dh_snowice(ji) * rhoi * a_i_1d(ji) * r1_Dt_ice ! and get rn_icesal from the ocean 594 539 ENDIF 595 540 596 541 ! Mass flux: All snow is thrown in the ocean, and seawater is taken to replace the volume 597 wfx_sni_1d(ji) = wfx_sni_1d(ji) - a_i_1d(ji) * dh_snowice(ji) * rhoi * r1_Dt_ice 598 wfx_snw_sni_1d(ji) = wfx_snw_sni_1d(ji) + a_i_1d(ji) * dh_snowice(ji) * rhos * r1_Dt_ice 542 wfx_sni_1d (ji) = wfx_sni_1d (ji) - dh_snowice(ji) * rhoi * a_i_1d(ji) * r1_Dt_ice 543 wfx_snw_sni_1d(ji) = wfx_snw_sni_1d(ji) + dh_snowice(ji) * rhos * a_i_1d(ji) * r1_Dt_ice 544 545 ! update thickness 546 zh_i(ji,0) = zh_i(ji,0) + dh_snowice(ji) 547 zdeltah(ji) = dh_snowice(ji) 599 548 600 549 ! update heat content (J.m-2) and layer thickness 601 eh_i_old(ji,0) = eh_i_old(ji,0) + dh_snowice(ji) * e_s_1d(ji,1) + zfmdt * zEw602 550 h_i_old (ji,0) = h_i_old (ji,0) + dh_snowice(ji) 603 604 END DO 605 606 ! 607 ! Update temperature, energy 608 ! -------------------------- 609 DO ji = 1, npti 610 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - h_i_1d(ji) ) ) 611 t_su_1d(ji) = rswitch * t_su_1d(ji) + ( 1._wp - rswitch ) * rt0 612 END DO 613 551 eh_i_old(ji,0) = eh_i_old(ji,0) + zfmdt * zEw ! 1st part (sea water enthalpy) 552 553 END DO 554 ! 555 DO jk = nlay_s, 0, -1 ! flooding of snow starts from the base 556 DO ji = 1, npti 557 zdum = MIN( zdeltah(ji), zh_s(ji,jk) ) ! amount of snw that floods, > 0 558 zh_s(ji,jk) = MAX( 0._wp, zh_s(ji,jk) - zdum ) ! remove some snow thickness 559 eh_i_old(ji,0) = eh_i_old(ji,0) + zdum * ze_s(ji,jk) ! 2nd part (snow enthalpy) 560 ! update dh_snowice 561 zdeltah(ji) = MAX( 0._wp, zdeltah(ji) - zdum ) 562 END DO 563 END DO 564 ! 565 ! 566 !!$ ! --- Update snow diags --- ! 567 !!$ !!clem: this is wrong. dh_s_tot is not used anyway 568 !!$ DO ji = 1, npti 569 !!$ dh_s_tot(ji) = dh_s_tot(ji) + dh_s_mlt(ji) + zdeltah(ji) + zdh_s_sub(ji) - dh_snowice(ji) 570 !!$ END DO 571 ! 572 ! 573 ! Remapping of snw enthalpy on a regular grid 574 !-------------------------------------------- 575 CALL snw_ent( zh_s, ze_s, e_s_1d ) 576 577 ! recalculate t_s_1d from e_s_1d 614 578 DO jk = 1, nlay_s 615 579 DO ji = 1,npti 616 ! where there is no ice or no snow 617 rswitch = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, - h_s_1d(ji) ) ) ) * ( 1._wp - MAX( 0._wp, SIGN(1._wp, - h_i_1d(ji) ) ) ) 618 ! mass & energy loss to the ocean 619 hfx_res_1d(ji) = hfx_res_1d(ji) + ( 1._wp - rswitch ) * & 620 & ( e_s_1d(ji,jk) * h_s_1d(ji) * r1_nlay_s * a_i_1d(ji) * r1_Dt_ice ) ! heat flux to the ocean [W.m-2], < 0 621 wfx_res_1d(ji) = wfx_res_1d(ji) + ( 1._wp - rswitch ) * & 622 & ( rhos * h_s_1d(ji) * r1_nlay_s * a_i_1d(ji) * r1_Dt_ice ) ! mass flux 623 ! update energy (mass is updated in the next loop) 624 e_s_1d(ji,jk) = rswitch * e_s_1d(ji,jk) 625 ! recalculate t_s_1d from e_s_1d 626 t_s_1d(ji,jk) = rt0 + rswitch * ( - e_s_1d(ji,jk) * r1_rhos * r1_rcpi + rLfus * r1_rcpi ) 627 END DO 628 END DO 580 IF( h_s_1d(ji) > 0._wp ) THEN 581 t_s_1d(ji,jk) = rt0 + ( - e_s_1d(ji,jk) * r1_rhos * r1_rcpi + rLfus * r1_rcpi ) 582 ELSE 583 t_s_1d(ji,jk) = rt0 584 ENDIF 585 END DO 586 END DO 587 588 ! Note: remapping of ice enthalpy is done in icethd.F90 629 589 630 590 ! --- ensure that a_i = 0 & h_s = 0 where h_i = 0 --- 631 591 WHERE( h_i_1d(1:npti) == 0._wp ) 632 a_i_1d(1:npti) = 0._wp 633 h_s_1d(1:npti) = 0._wp 592 a_i_1d (1:npti) = 0._wp 593 h_s_1d (1:npti) = 0._wp 594 t_su_1d(1:npti) = rt0 634 595 END WHERE 635 !596 636 597 END SUBROUTINE ice_thd_dh 637 598 638 639 !!-------------------------------------------------------------------------- 640 !! INTERFACE ice_thd_snwblow 641 !! 642 !! ** Purpose : Compute distribution of precip over the ice 643 !! 644 !! Snow accumulation in one thermodynamic time step 645 !! snowfall is partitionned between leads and ice. 646 !! If snow fall was uniform, a fraction (1-at_i) would fall into leads 647 !! but because of the winds, more snow falls on leads than on sea ice 648 !! and a greater fraction (1-at_i)^beta of the total mass of snow 649 !! (beta < 1) falls in leads. 650 !! In reality, beta depends on wind speed, 651 !! and should decrease with increasing wind speed but here, it is 652 !! considered as a constant. an average value is 0.66 653 !!-------------------------------------------------------------------------- 654 !!gm I think it can be usefull to set this as a FUNCTION, not a SUBROUTINE.... 655 SUBROUTINE ice_thd_snwblow_2d( pin, pout ) 656 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pin ! previous fraction lead ( 1. - a_i_b ) 657 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pout 658 pout = ( 1._wp - ( pin )**rn_blow_s ) 659 END SUBROUTINE ice_thd_snwblow_2d 660 661 SUBROUTINE ice_thd_snwblow_1d( pin, pout ) 662 REAL(wp), DIMENSION(:), INTENT(in ) :: pin 663 REAL(wp), DIMENSION(:), INTENT(inout) :: pout 664 pout = ( 1._wp - ( pin )**rn_blow_s ) 665 END SUBROUTINE ice_thd_snwblow_1d 666 599 SUBROUTINE snw_ent( ph_old, pe_old, pe_new ) 600 !!------------------------------------------------------------------- 601 !! *** ROUTINE snw_ent *** 602 !! 603 !! ** Purpose : 604 !! This routine computes new vertical grids in the snow, 605 !! and consistently redistributes temperatures. 606 !! Redistribution is made so as to ensure to energy conservation 607 !! 608 !! 609 !! ** Method : linear conservative remapping 610 !! 611 !! ** Steps : 1) cumulative integrals of old enthalpies/thicknesses 612 !! 2) linear remapping on the new layers 613 !! 614 !! ------------ cum0(0) ------------- cum1(0) 615 !! NEW ------------- 616 !! ------------ cum0(1) ==> ------------- 617 !! ... ------------- 618 !! ------------ ------------- 619 !! ------------ cum0(nlay_s+1) ------------- cum1(nlay_s) 620 !! 621 !! 622 !! References : Bitz & Lipscomb, JGR 99; Vancoppenolle et al., GRL, 2005 623 !!------------------------------------------------------------------- 624 REAL(wp), DIMENSION(jpij,0:nlay_s), INTENT(in ) :: ph_old ! old thicknesses (m) 625 REAL(wp), DIMENSION(jpij,0:nlay_s), INTENT(in ) :: pe_old ! old enthlapies (J.m-3) 626 REAL(wp), DIMENSION(jpij,1:nlay_s), INTENT(inout) :: pe_new ! new enthlapies (J.m-3, remapped) 627 ! 628 INTEGER :: ji ! dummy loop indices 629 INTEGER :: jk0, jk1 ! old/new layer indices 630 ! 631 REAL(wp), DIMENSION(jpij,0:nlay_s+1) :: zeh_cum0, zh_cum0 ! old cumulative enthlapies and layers interfaces 632 REAL(wp), DIMENSION(jpij,0:nlay_s) :: zeh_cum1, zh_cum1 ! new cumulative enthlapies and layers interfaces 633 REAL(wp), DIMENSION(jpij) :: zhnew ! new layers thicknesses 634 !!------------------------------------------------------------------- 635 636 !-------------------------------------------------------------------------- 637 ! 1) Cumulative integral of old enthalpy * thickness and layers interfaces 638 !-------------------------------------------------------------------------- 639 zeh_cum0(1:npti,0) = 0._wp 640 zh_cum0 (1:npti,0) = 0._wp 641 DO jk0 = 1, nlay_s+1 642 DO ji = 1, npti 643 zeh_cum0(ji,jk0) = zeh_cum0(ji,jk0-1) + pe_old(ji,jk0-1) * ph_old(ji,jk0-1) 644 zh_cum0 (ji,jk0) = zh_cum0 (ji,jk0-1) + ph_old(ji,jk0-1) 645 END DO 646 END DO 647 648 !------------------------------------ 649 ! 2) Interpolation on the new layers 650 !------------------------------------ 651 ! new layer thickesses 652 DO ji = 1, npti 653 zhnew(ji) = SUM( ph_old(ji,0:nlay_s) ) * r1_nlay_s 654 END DO 655 656 ! new layers interfaces 657 zh_cum1(1:npti,0) = 0._wp 658 DO jk1 = 1, nlay_s 659 DO ji = 1, npti 660 zh_cum1(ji,jk1) = zh_cum1(ji,jk1-1) + zhnew(ji) 661 END DO 662 END DO 663 664 zeh_cum1(1:npti,0:nlay_s) = 0._wp 665 ! new cumulative q*h => linear interpolation 666 DO jk0 = 1, nlay_s+1 667 DO jk1 = 1, nlay_s-1 668 DO ji = 1, npti 669 IF( zh_cum1(ji,jk1) <= zh_cum0(ji,jk0) .AND. zh_cum1(ji,jk1) > zh_cum0(ji,jk0-1) ) THEN 670 zeh_cum1(ji,jk1) = ( zeh_cum0(ji,jk0-1) * ( zh_cum0(ji,jk0) - zh_cum1(ji,jk1 ) ) + & 671 & zeh_cum0(ji,jk0 ) * ( zh_cum1(ji,jk1) - zh_cum0(ji,jk0-1) ) ) & 672 & / ( zh_cum0(ji,jk0) - zh_cum0(ji,jk0-1) ) 673 ENDIF 674 END DO 675 END DO 676 END DO 677 ! to ensure that total heat content is strictly conserved, set: 678 zeh_cum1(1:npti,nlay_s) = zeh_cum0(1:npti,nlay_s+1) 679 680 ! new enthalpies 681 DO jk1 = 1, nlay_s 682 DO ji = 1, npti 683 rswitch = MAX( 0._wp , SIGN( 1._wp , zhnew(ji) - epsi20 ) ) 684 pe_new(ji,jk1) = rswitch * ( zeh_cum1(ji,jk1) - zeh_cum1(ji,jk1-1) ) / MAX( zhnew(ji), epsi20 ) 685 END DO 686 END DO 687 688 END SUBROUTINE snw_ent 689 690 667 691 #else 668 692 !!----------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.
