Changeset 9274
- Timestamp:
- 2018-01-22T17:27:39+01:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_merge_2017/NEMOGCM/NEMO/LIM_SRC_3/icethd_dh.F90
r9271 r9274 31 31 32 32 PUBLIC ice_thd_dh ! called by ice_thd 33 PUBLIC ice_thd_snwblow ! called in sbcblk/sbcc lio/sbccpl and here33 PUBLIC ice_thd_snwblow ! called in sbcblk/sbccpl and here 34 34 35 35 INTERFACE ice_thd_snwblow … … 48 48 !! *** ROUTINE ice_thd_dh *** 49 49 !! 50 !! ** Purpose : compute ice and snow thickness changes due to grow ing/melting50 !! ** Purpose : compute ice and snow thickness changes due to growth/melting 51 51 !! 52 52 !! ** Method : Ice/Snow surface melting arises from imbalance in surface fluxes 53 !! Bottom accretion/ablation arises from flux budget54 !! Snow thickness can increase by precipitation and decrease by sublimation55 !! If snow load excesses Archmiede limit, snow-ice is formed by56 !! the flooding of sea-water in the snow53 !! Bottom accretion/ablation arises from flux budget 54 !! Snow thickness can increase by precipitation and decrease by sublimation 55 !! If snow load excesses Archmiede limit, snow-ice is formed by 56 !! the flooding of sea-water in the snow 57 57 !! 58 !! 1) Compute available flux of heat for surface ablation 59 !! 2) Compute snow and sea ice enthalpies 60 !! 3) Surface ablation and sublimation 61 !! 4) Bottom accretion/ablation 62 !! 5) Case of Total ablation 63 !! 6) Snow ice formation 58 !! - Compute available flux of heat for surface ablation 59 !! - Compute snow and sea ice enthalpies 60 !! - Surface ablation and sublimation 61 !! - Bottom accretion/ablation 62 !! - Snow ice formation 64 63 !! 65 64 !! References : Bitz and Lipscomb, 1999, J. Geophys. Res. … … 126 125 END DO 127 126 ! 128 ! ------------------------------------------------------------------------------!129 ! 1) Calculate available heat for surface and bottom ablation!130 ! ------------------------------------------------------------------------------!131 ! 132 SELECT CASE 133 ! 134 CASE 127 ! ---------------------------------------------- ! 128 ! Available heat for surface and bottom ablation ! 129 ! ---------------------------------------------- ! 130 ! 131 SELECT CASE( nice_jules ) 132 ! 133 CASE( np_jules_ACTIVE ) 135 134 ! 136 135 DO ji = 1, npti 137 zq_su(ji) =MAX( 0._wp, qml_ice_1d(ji) * rdt_ice )136 zq_su(ji) = MAX( 0._wp, qml_ice_1d(ji) * rdt_ice ) 138 137 END DO 139 138 ! 140 CASE 139 CASE( np_jules_EMULE ) 141 140 ! 142 141 DO ji = 1, npti 143 zdum = ( qns_ice_1d(ji) + qsr_ice_1d(ji) - qsr_ice_tr_1d(ji) - fc_su(ji) )144 qml_ice_1d(ji) =zdum * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) )145 zq_su(ji) =MAX( 0._wp, qml_ice_1d(ji) * rdt_ice )142 zdum = qns_ice_1d(ji) + qsr_ice_1d(ji) - qsr_ice_tr_1d(ji) - fc_su(ji) 143 qml_ice_1d(ji) = zdum * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) ) 144 zq_su(ji) = MAX( 0._wp, qml_ice_1d(ji) * rdt_ice ) 146 145 END DO 147 146 ! 148 CASE 147 CASE( np_jules_OFF ) 149 148 ! 150 149 DO ji = 1, npti 151 zdum = ( qns_ice_1d(ji) + qsr_ice_1d(ji) - qsr_ice_tr_1d(ji) - fc_su(ji))152 zdum =zdum * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) )153 zq_su(ji) = MAX( 0._wp, zdum* rdt_ice )150 zdum = qns_ice_1d(ji) + qsr_ice_1d(ji) - qsr_ice_tr_1d(ji) - fc_su(ji) 151 zdum = zdum * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) ) 152 zq_su(ji) = MAX( 0._wp, zdum * rdt_ice ) 154 153 END DO 155 154 ! … … 157 156 ! 158 157 DO ji = 1, npti 159 zf_tt(ji) = fc_bo_i(ji) + fhtur_1d(ji) + fhld_1d(ji) 160 zq_bo (ji) = MAX( 0._wp, zf_tt(ji) * rdt_ice ) 161 END DO 162 163 !--------------------------------------------! 164 ! 2) Computing layer thicknesses ! 165 !--------------------------------------------! 158 zf_tt(ji) = fc_bo_i(ji) + fhtur_1d(ji) + fhld_1d(ji) 159 zq_bo(ji) = MAX( 0._wp, zf_tt(ji) * rdt_ice ) 160 END DO 161 162 ! Ice and snow layer thicknesses 163 !------------------------------- 166 164 DO jk = 1, nlay_i 167 165 DO ji = 1, npti … … 176 174 END DO 177 175 178 !------------------------------------------------------------------------------! 179 ! If snow temperature is above freezing point, then snow melts 180 ! (should not happen but sometimes it does) 181 !------------------------------------------------------------------------------! 176 ! ! ============ ! 177 ! ! Snow ! 178 ! ! ============ ! 179 ! 180 ! Internal melting 181 ! ---------------- 182 ! IF snow temperature is above freezing point, THEN snow melts (should not happen but sometimes it does) 182 183 DO jk = 1, nlay_s 183 184 DO ji = 1, npti 184 IF( t_s_1d(ji,jk) > rt0 ) THEN !!! Internal melting 185 ! Contribution to heat flux to the ocean [W.m-2], < 0 186 hfx_res_1d(ji) = hfx_res_1d(ji) + e_s_1d(ji,jk) * zh_s(ji,jk) * a_i_1d(ji) * r1_rdtice 187 ! Contribution to mass flux 188 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + rhosn * zh_s(ji,jk) * a_i_1d(ji) * r1_rdtice 189 dh_s_mlt(ji) = dh_s_mlt(ji) - zh_s(ji,jk) 185 IF( t_s_1d(ji,jk) > rt0 ) THEN 186 hfx_res_1d (ji) = hfx_res_1d (ji) + e_s_1d(ji,jk) * zh_s(ji,jk) * a_i_1d(ji) * r1_rdtice ! heat flux to the ocean [W.m-2], < 0 187 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + rhosn * zh_s(ji,jk) * a_i_1d(ji) * r1_rdtice ! mass flux 190 188 ! updates 191 h_s_1d(ji) = h_s_1d(ji) - zh_s(ji,jk) 192 zh_s (ji,jk) = 0._wp 193 e_s_1d(ji,jk) = 0._wp 194 t_s_1d(ji,jk) = rt0 189 dh_s_mlt(ji) = dh_s_mlt(ji) - zh_s(ji,jk) 190 h_s_1d (ji) = h_s_1d(ji) - zh_s(ji,jk) 191 zh_s (ji,jk) = 0._wp 192 e_s_1d (ji,jk) = 0._wp 193 t_s_1d (ji,jk) = rt0 195 194 END IF 196 195 END DO 197 END DO 198 199 200 !------------------------------------------------------------------------------| 201 ! 3) Surface ablation and sublimation | 202 !------------------------------------------------------------------------------| 203 ! 204 !------------------------- 205 ! 3.1 Snow precips / melt 206 !------------------------- 207 ! Snow accumulation in one thermodynamic time step 208 ! snowfall is partitionned between leads and ice 209 ! if snow fall was uniform, a fraction (1-at_i) would fall into leads 210 ! but because of the winds, more snow falls on leads than on sea ice 211 ! and a greater fraction (1-at_i)^beta of the total mass of snow 212 ! (beta < 1) falls in leads. 213 ! In reality, beta depends on wind speed, 214 ! and should decrease with increasing wind speed but here, it is 215 ! considered as a constant. an average value is 0.66 216 ! Martin Vancoppenolle, December 2006 217 218 CALL ice_thd_snwblow( 1. - at_i_1d(1:npti), zsnw(1:npti) ) ! snow distribution over ice after wind blowing 196 END DO 197 198 ! Snow precipitation 199 !------------------- 200 CALL ice_thd_snwblow( 1. - at_i_1d(1:npti), zsnw(1:npti) ) ! snow distribution over ice after wind blowing 219 201 220 202 zdeltah(1:npti,:) = 0._wp 221 203 DO ji = 1, npti 222 !----------- 223 ! Snow fall 224 !----------- 225 ! thickness change 226 zdh_s_pre(ji) = zsnw(ji) * sprecip_1d(ji) * rdt_ice * r1_rhosn / at_i_1d(ji) 227 ! enthalpy of the precip (>0, J.m-3) 228 zqprec (ji) = - qprec_ice_1d(ji) 229 IF( sprecip_1d(ji) == 0._wp ) zqprec(ji) = 0._wp 230 ! heat flux from snow precip (>0, W.m-2) 231 hfx_spr_1d(ji) = hfx_spr_1d(ji) + zdh_s_pre(ji) * a_i_1d(ji) * zqprec(ji) * r1_rdtice 232 ! mass flux, <0 233 wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_pre(ji) * r1_rdtice 234 235 !--------------------- 236 ! Melt of falling snow 237 !--------------------- 238 ! thickness change 239 rswitch = MAX( 0._wp , SIGN( 1._wp , zqprec(ji) - epsi20 ) ) 240 zdeltah (ji,1) = - rswitch * zq_su(ji) / MAX( zqprec(ji) , epsi20 ) 241 zdeltah (ji,1) = MAX( - zdh_s_pre(ji), zdeltah(ji,1) ) ! bound melting 242 ! heat used to melt snow (W.m-2, >0) 243 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * zqprec(ji) * r1_rdtice 244 ! snow melting only = water into the ocean (then without snow precip), >0 245 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice 246 dh_s_mlt(ji) = dh_s_mlt(ji) + zdeltah(ji,1) 247 ! updates available heat + precipitations after melting 248 zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,1) * zqprec(ji) ) 249 zdh_s_pre (ji) = zdh_s_pre(ji) + zdeltah(ji,1) 250 251 ! update thickness 252 h_s_1d(ji) = MAX( 0._wp , h_s_1d(ji) + zdh_s_pre(ji) ) 204 IF( sprecip_1d(ji) > 0._wp ) THEN 205 ! 206 ! --- precipitation --- 207 zdh_s_pre (ji) = zsnw(ji) * sprecip_1d(ji) * rdt_ice * r1_rhosn / at_i_1d(ji) ! thickness change 208 zqprec (ji) = - qprec_ice_1d(ji) ! enthalpy of the precip (>0, J.m-3) 209 ! 210 hfx_spr_1d(ji) = hfx_spr_1d(ji) + zdh_s_pre(ji) * a_i_1d(ji) * zqprec(ji) * r1_rdtice ! heat flux from snow precip (>0, W.m-2) 211 wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_pre(ji) * r1_rdtice ! mass flux, <0 212 213 ! --- melt of falling snow --- 214 rswitch = MAX( 0._wp , SIGN( 1._wp , zqprec(ji) - epsi20 ) ) 215 zdeltah (ji,1) = - rswitch * zq_su(ji) / MAX( zqprec(ji) , epsi20 ) ! thickness change 216 zdeltah (ji,1) = MAX( - zdh_s_pre(ji), zdeltah(ji,1) ) ! bound melting 217 hfx_snw_1d (ji) = hfx_snw_1d (ji) - zdeltah(ji,1) * a_i_1d(ji) * zqprec(ji) * r1_rdtice ! heat used to melt snow (W.m-2, >0) 218 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice ! snow melting only = water into the ocean (then without snow precip), >0 219 220 ! updates available heat + precipitations after melting 221 dh_s_mlt (ji) = dh_s_mlt(ji) + zdeltah(ji,1) 222 zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,1) * zqprec(ji) ) 223 zdh_s_pre(ji) = zdh_s_pre(ji) + zdeltah(ji,1) 224 225 ! update thickness 226 h_s_1d(ji) = MAX( 0._wp , h_s_1d(ji) + zdh_s_pre(ji) ) 227 ! 228 ELSE 229 ! 230 zdh_s_pre(ji) = 0._wp 231 zqprec (ji) = 0._wp 232 ! 233 ENDIF 253 234 END DO 254 235 … … 260 241 END DO 261 242 243 ! Snow melting 244 ! ------------ 262 245 ! If heat still available (zq_su > 0), then melt more snow 263 246 zdeltah(1:npti,:) = 0._wp … … 265 248 DO jk = 1, nlay_s 266 249 DO ji = 1, npti 267 ! thickness change 268 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, - zh_s(ji,jk) ) ) 269 rswitch = rswitch * ( MAX( 0._wp, SIGN( 1._wp, e_s_1d(ji,jk) - epsi20 ) ) ) 270 zdeltah (ji,jk) = - rswitch * zq_su(ji) / MAX( e_s_1d(ji,jk), epsi20 ) 271 zdeltah (ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji,jk) ) ! bound melting 272 zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,jk) 273 ! heat used to melt snow(W.m-2, >0) 274 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,jk) * a_i_1d(ji) * e_s_1d(ji,jk) * r1_rdtice 275 ! snow melting only = water into the ocean (then without snow precip) 276 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 277 dh_s_mlt(ji) = dh_s_mlt(ji) + zdeltah(ji,jk) 278 ! updates available heat + thickness 279 zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * e_s_1d(ji,jk) ) 280 h_s_1d(ji) = MAX( 0._wp , h_s_1d(ji) + zdeltah(ji,jk) ) 281 zh_s(ji,jk)= MAX( 0._wp , zh_s(ji,jk) + zdeltah(ji,jk) ) 282 END DO 283 END DO 284 285 !------------------------------ 286 ! 3.2 Sublimation (part1: snow) 287 !------------------------------ 250 IF( zh_s(ji,jk) > 0._wp .AND. zq_su(ji) > 0._wp ) THEN 251 ! 252 rswitch = MAX( 0._wp, SIGN( 1._wp, e_s_1d(ji,jk) - epsi20 ) ) 253 zdeltah (ji,jk) = - rswitch * zq_su(ji) / MAX( e_s_1d(ji,jk), epsi20 ) ! thickness change 254 zdeltah (ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji,jk) ) ! bound melting 255 zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,jk) 256 257 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,jk) * a_i_1d(ji) * e_s_1d (ji,jk) * r1_rdtice ! heat used to melt snow(W.m-2, >0) 258 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice ! snow melting only = water into the ocean (then without snow precip) 259 260 ! updates available heat + thickness 261 dh_s_mlt(ji) = dh_s_mlt(ji) + zdeltah(ji,jk) 262 zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * e_s_1d(ji,jk) ) 263 h_s_1d (ji) = MAX( 0._wp , h_s_1d(ji) + zdeltah(ji,jk) ) 264 zh_s (ji,jk) = MAX( 0._wp , zh_s(ji,jk) + zdeltah(ji,jk) ) 265 ! 266 ENDIF 267 END DO 268 END DO 269 270 ! Snow sublimation 271 !----------------- 288 272 ! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates 289 273 ! comment: not counted in mass/heat exchange in iceupdate.F90 since this is an exchange with atm. (not ocean) 290 274 zdeltah(1:npti,:) = 0._wp 291 275 DO ji = 1, npti 292 zdh_s_sub(ji) = MAX( - h_s_1d(ji) , - evap_ice_1d(ji) * r1_rhosn * rdt_ice ) 293 ! remaining evap in kg.m-2 (used for ice melting later on) 294 zevap_rema(ji) = evap_ice_1d(ji) * rdt_ice + zdh_s_sub(ji) * rhosn 295 ! Heat flux by sublimation [W.m-2], < 0 (sublimate first snow that had fallen, then pre-existing snow) 296 zdeltah(ji,1) = MAX( zdh_s_sub(ji), - zdh_s_pre(ji) ) 297 hfx_sub_1d(ji) = hfx_sub_1d(ji) + ( zdeltah(ji,1) * zqprec(ji) + ( zdh_s_sub(ji) - zdeltah(ji,1) ) * e_s_1d(ji,1) & 298 & ) * a_i_1d(ji) * r1_rdtice 299 ! Mass flux by sublimation 300 wfx_snw_sub_1d(ji) = wfx_snw_sub_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_sub(ji) * r1_rdtice 301 302 ! new snow thickness 303 h_s_1d(ji) = MAX( 0._wp , h_s_1d(ji) + zdh_s_sub(ji) ) 304 ! update precipitations after sublimation and correct sublimation 305 zdh_s_pre(ji) = zdh_s_pre(ji) + zdeltah(ji,1) 306 zdh_s_sub(ji) = zdh_s_sub(ji) - zdeltah(ji,1) 276 IF( evap_ice_1d(ji) > 0._wp ) THEN 277 ! 278 zdh_s_sub (ji) = MAX( - h_s_1d(ji) , - evap_ice_1d(ji) * r1_rhosn * rdt_ice ) 279 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) 280 zdeltah (ji,1) = MAX( zdh_s_sub(ji), - zdh_s_pre(ji) ) 281 282 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) 283 & ( zdeltah(ji,1) * zqprec(ji) + ( zdh_s_sub(ji) - zdeltah(ji,1) ) * e_s_1d(ji,1) ) & 284 & * a_i_1d(ji) * r1_rdtice 285 wfx_snw_sub_1d(ji) = wfx_snw_sub_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_sub(ji) * r1_rdtice ! Mass flux by sublimation 286 287 ! new snow thickness 288 h_s_1d(ji) = MAX( 0._wp , h_s_1d(ji) + zdh_s_sub(ji) ) 289 ! update precipitations after sublimation and correct sublimation 290 zdh_s_pre(ji) = zdh_s_pre(ji) + zdeltah(ji,1) 291 zdh_s_sub(ji) = zdh_s_sub(ji) - zdeltah(ji,1) 292 ! 293 ELSE 294 ! 295 zdh_s_sub (ji) = 0._wp 296 zevap_rema(ji) = 0._wp 297 ! 298 ENDIF 307 299 END DO 308 300 … … 312 304 END DO 313 305 314 !------------------------------------------- 315 ! 3.3 Update temperature, energy 316 !------------------------------------------- 306 ! Update temperature, energy 307 !--------------------------- 317 308 ! new temp and enthalpy of the snow (remaining snow precip + remaining pre-existing snow) 318 309 DO jk = 1, nlay_s 319 310 DO ji = 1,npti 320 311 rswitch = MAX( 0._wp , SIGN( 1._wp, h_s_1d(ji) - epsi20 ) ) 321 e_s_1d(ji,jk) = rswitch / MAX( h_s_1d(ji), epsi20 ) * & 322 & ( ( zdh_s_pre(ji) ) * zqprec(ji) + & 323 & ( h_s_1d(ji) - zdh_s_pre(ji) ) * rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) ) 324 END DO 325 END DO 326 327 !-------------------------- 328 ! 3.4 Surface ice ablation 329 !-------------------------- 312 e_s_1d(ji,jk) = rswitch / MAX( h_s_1d(ji), epsi20 ) * & 313 & ( ( zdh_s_pre(ji) ) * zqprec(ji) + & 314 & ( h_s_1d(ji) - zdh_s_pre(ji) ) * rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) ) 315 END DO 316 END DO 317 318 ! ! ============ ! 319 ! ! Ice ! 320 ! ! ============ ! 321 322 ! Surface ice melting 323 !-------------------- 330 324 zdeltah(1:npti,:) = 0._wp ! important 331 325 DO jk = 1, nlay_i 332 326 DO ji = 1, npti 333 ztmelts = - tmut * sz_i_1d(ji,jk)! Melting point of layer k [C]334 335 IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN !!!Internal melting327 ztmelts = - tmut * sz_i_1d(ji,jk) ! Melting point of layer k [C] 328 329 IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN !-- Internal melting 336 330 337 331 zEi = - e_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of layer k [J/kg, <0] 338 zdE = 0._wp ! Specific enthalpy difference 332 zdE = 0._wp ! Specific enthalpy difference (J/kg, <0) 339 333 ! set up at 0 since no energy is needed to melt water...(it is already melted) 340 334 zdeltah(ji,jk) = MIN( 0._wp , - zh_i(ji,jk) ) ! internal melting occurs when the internal temperature is above freezing … … 346 340 zfmdt = - rhoic * zdeltah(ji,jk) ! Recompute mass flux [kg/m2, >0] 347 341 348 ! Contribution to heat flux to the ocean [W.m-2], <0 (ice enthalpy zEi is "sent" to the ocean) 349 hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice 350 351 ! Contribution to salt flux (using s_i_1d and not sz_i_1d(jk) is ok) 352 sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice 353 354 ! Contribution to mass flux 355 wfx_res_1d(ji) = wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 356 357 ELSE !!! Surface melting 342 hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice ! Heat flux to the ocean [W.m-2], <0 343 ! ice enthalpy zEi is "sent" to the ocean 344 sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice ! Salt flux 345 ! using s_i_1d and not sz_i_1d(jk) is ok 346 wfx_res_1d(ji) = wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice ! Mass flux 347 348 ELSE !-- Surface melting 358 349 359 350 zEi = - e_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of layer k [J/kg, <0] … … 375 366 zQm = zfmdt * zEw ! Energy of the melt water sent to the ocean [J/m2, <0] 376 367 377 ! Contribution to salt flux >0 (using s_i_1d and not sz_i_1d(jk) is ok) 378 sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice 379 380 ! Contribution to heat flux [W.m-2], < 0 381 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 382 383 ! Total heat flux used in this process [W.m-2], > 0 384 hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice 385 386 ! Contribution to mass flux 387 wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 368 sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice ! Salt flux >0 369 ! using s_i_1d and not sz_i_1d(jk) is ok) 370 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice ! Heat flux [W.m-2], < 0 371 hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice ! Heat flux used in this process [W.m-2], > 0 372 ! 373 wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice ! Mass flux 388 374 389 375 END IF 390 ! ---------------------- 391 ! Sublimation part2: ice 392 ! ---------------------- 393 zdum = MAX( - ( zh_i(ji,jk) + zdeltah(ji,jk) ) , - zevap_rema(ji) * r1_rhoic ) 394 zdeltah(ji,jk) = zdeltah(ji,jk) + zdum 395 dh_i_sub(ji) = dh_i_sub(ji) + zdum 396 ! Salt flux > 0 (clem: flux is sent to the ocean for simplicity but salt should remain in the ice except if all ice is melted. 397 ! It must be corrected at some point) 398 sfx_sub_1d(ji) = sfx_sub_1d(ji) - rhoic * a_i_1d(ji) * zdum * s_i_1d(ji) * r1_rdtice 399 ! Heat flux [W.m-2], < 0 400 hfx_sub_1d(ji) = hfx_sub_1d(ji) + zdum * e_i_1d(ji,jk) * a_i_1d(ji) * r1_rdtice 401 ! Mass flux > 0 402 wfx_ice_sub_1d(ji) = wfx_ice_sub_1d(ji) - rhoic * a_i_1d(ji) * zdum * r1_rdtice 376 377 ! Ice sublimation 378 ! --------------- 379 zdum = MAX( - ( zh_i(ji,jk) + zdeltah(ji,jk) ) , - zevap_rema(ji) * r1_rhoic ) 380 zdeltah (ji,jk) = zdeltah (ji,jk) + zdum 381 dh_i_sub(ji) = dh_i_sub(ji) + zdum 382 383 sfx_sub_1d(ji) = sfx_sub_1d(ji) - rhoic * a_i_1d(ji) * zdum * s_i_1d(ji) * r1_rdtice ! Salt flux >0 384 ! clem: flux is sent to the ocean for simplicity 385 ! but salt should remain in the ice except 386 ! if all ice is melted. => must be corrected 387 hfx_sub_1d(ji) = hfx_sub_1d(ji) + zdum * e_i_1d(ji,jk) * a_i_1d(ji) * r1_rdtice ! Heat flux [W.m-2], < 0 388 389 wfx_ice_sub_1d(ji) = wfx_ice_sub_1d(ji) - rhoic * a_i_1d(ji) * zdum * r1_rdtice ! Mass flux > 0 403 390 404 391 ! update remaining mass flux … … 417 404 END DO 418 405 END DO 406 419 407 ! update ice thickness 420 408 DO ji = 1, npti … … 428 416 429 417 430 !------------------------------------------------------------------------------! 431 ! 4) Basal growth / melt ! 432 !------------------------------------------------------------------------------! 433 ! 434 !------------------ 435 ! 4.1 Basal growth 418 ! Ice Basal growth 436 419 !------------------ 437 420 ! Basal growth is driven by heat imbalance at the ice-ocean interface, … … 449 432 IF( nn_icesal == 2 ) num_iter_max = 5 ! salinity varying in time 450 433 451 ! Iterative procedure452 434 DO ji = 1, npti 453 435 IF( zf_tt(ji) < 0._wp ) THEN 454 DO iter = 1, num_iter_max 436 DO iter = 1, num_iter_max ! iterations 455 437 456 438 ! New bottom ice salinity (Cox & Weeks, JGR88 ) … … 458 440 !--- zswi12 if 2.0e-8 < dh/dt < 3.6e-7 459 441 !--- zswi2 if dh/dt > 3.6e-7 460 zgrr = MIN( 1.0e-3, MAX ( dh_i_bott(ji) * r1_rdtice , epsi10 ) ) 461 zswi2 = MAX( 0._wp , SIGN( 1._wp , zgrr - 3.6e-7 ) ) 462 zswi12 = MAX( 0._wp , SIGN( 1._wp , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 ) 463 zswi1 = 1. - zswi2 * zswi12 464 zfracs = MIN ( zswi1 * 0.12 + zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) ) & 465 & + zswi2 * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) ) , 0.5 ) 466 467 s_i_new(ji) = zswitch_sal * zfracs * sss_1d(ji) & ! New ice salinity 468 + ( 1. - zswitch_sal ) * s_i_1d(ji) 469 ! New ice growth 470 ztmelts = - tmut * s_i_new(ji) ! New ice melting point (C) 471 472 zt_i_new = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i) 473 474 zEi = cpic * ( zt_i_new - (ztmelts+rt0) ) & ! Specific enthalpy of forming ice (J/kg, <0) 475 & - lfus * ( 1.0 - ztmelts / ( zt_i_new - rt0 ) ) & 476 & + rcp * ztmelts 477 478 zEw = rcp * ( t_bo_1d(ji) - rt0 ) ! Specific enthalpy of seawater (J/kg, < 0) 479 480 zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0) 481 482 dh_i_bott(ji) = rdt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoic ) ) 442 zgrr = MIN( 1.0e-3, MAX ( dh_i_bott(ji) * r1_rdtice , epsi10 ) ) 443 zswi2 = MAX( 0._wp , SIGN( 1._wp , zgrr - 3.6e-7 ) ) 444 zswi12 = MAX( 0._wp , SIGN( 1._wp , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 ) 445 zswi1 = 1. - zswi2 * zswi12 446 zfracs = MIN( zswi1 * 0.12 + zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) ) & 447 & + zswi2 * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) ) , 0.5 ) 448 449 s_i_new(ji) = zswitch_sal * zfracs * sss_1d(ji) + ( 1. - zswitch_sal ) * s_i_1d(ji) ! New ice salinity 450 451 ztmelts = - tmut * s_i_new(ji) ! New ice melting point (C) 452 453 zt_i_new = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i) 454 455 zEi = cpic * ( zt_i_new - (ztmelts+rt0) ) & ! Specific enthalpy of forming ice (J/kg, <0) 456 & - lfus * ( 1.0 - ztmelts / ( zt_i_new - rt0 ) ) + rcp * ztmelts 457 458 zEw = rcp * ( t_bo_1d(ji) - rt0 ) ! Specific enthalpy of seawater (J/kg, < 0) 459 460 zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0) 461 462 dh_i_bott(ji) = rdt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoic ) ) 483 463 484 464 END DO 485 465 ! Contribution to Energy and Salt Fluxes 486 zfmdt = - rhoic * dh_i_bott(ji) ! Mass flux x time step (kg/m2, < 0)487 488 ztmelts = - tmut * s_i_new(ji) ! New ice melting point (C)466 zfmdt = - rhoic * dh_i_bott(ji) ! Mass flux x time step (kg/m2, < 0) 467 468 ztmelts = - tmut * s_i_new(ji) ! New ice melting point (C) 489 469 490 470 zt_i_new = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i) 491 471 492 zEi = cpic * ( zt_i_new - (ztmelts+rt0) ) & ! Specific enthalpy of forming ice (J/kg, <0) 493 & - lfus * ( 1.0 - ztmelts / ( zt_i_new - rt0 ) ) & 494 & + rcp * ztmelts 495 496 zEw = rcp * ( t_bo_1d(ji) - rt0 ) ! Specific enthalpy of seawater (J/kg, < 0) 497 498 zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0) 499 500 ! Contribution to heat flux to the ocean [W.m-2], >0 501 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 502 503 ! Total heat flux used in this process [W.m-2], <0 504 hfx_bog_1d(ji) = hfx_bog_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice 505 506 ! Contribution to salt flux, <0 507 sfx_bog_1d(ji) = sfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bott(ji) * s_i_new(ji) * r1_rdtice 508 509 ! Contribution to mass flux, <0 510 wfx_bog_1d(ji) = wfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bott(ji) * r1_rdtice 472 zEi = cpic * ( zt_i_new - (ztmelts+rt0) ) & ! Specific enthalpy of forming ice (J/kg, <0) 473 & - lfus * ( 1.0 - ztmelts / ( zt_i_new - rt0 ) ) + rcp * ztmelts 474 475 zEw = rcp * ( t_bo_1d(ji) - rt0 ) ! Specific enthalpy of seawater (J/kg, < 0) 476 477 zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0) 478 479 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice ! Heat flux to the ocean [W.m-2], >0 480 hfx_bog_1d(ji) = hfx_bog_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice ! Heat flux used in this process [W.m-2], <0 481 482 sfx_bog_1d(ji) = sfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bott(ji) * s_i_new(ji) * r1_rdtice ! Salt flux, <0 483 484 wfx_bog_1d(ji) = wfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bott(ji) * r1_rdtice ! Mass flux, <0 511 485 512 486 ! update heat content (J.m-2) and layer thickness … … 518 492 END DO 519 493 520 !---------------- 521 ! 4.2 Basal melt 522 !---------------- 494 ! Ice Basal melt 495 !--------------- 523 496 zdeltah(1:npti,:) = 0._wp ! important 524 497 DO jk = nlay_i, 1, -1 … … 528 501 ztmelts = - tmut * sz_i_1d(ji,jk) ! Melting point of layer jk (C) 529 502 530 IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN !!!Internal melting503 IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN !-- Internal melting 531 504 532 505 zEi = - e_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0) 533 506 zdE = 0._wp ! Specific enthalpy difference (J/kg, <0) 534 ! set up at 0 since no energy is needed to melt water...(it is already melted)507 ! set up at 0 since no energy is needed to melt water...(it is already melted) 535 508 zdeltah (ji,jk) = MIN( 0._wp , - zh_i(ji,jk) ) ! internal melting occurs when the internal temperature is above freezing 536 509 ! this should normally not happen, but sometimes, heat diffusion leads to this … … 540 513 zfmdt = - zdeltah(ji,jk) * rhoic ! Mass flux x time step > 0 541 514 542 ! Contribution to heat flux to the ocean [W.m-2], <0 (ice enthalpy zEi is "sent" to the ocean) 543 hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice 544 545 ! Contribution to salt flux (using s_i_1d and not sz_i_1d(jk) is ok) 546 sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice 547 548 ! Contribution to mass flux 549 wfx_res_1d(ji) = wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 515 hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice ! Heat flux to the ocean [W.m-2], <0 516 ! ice enthalpy zEi is "sent" to the ocean 517 sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice ! Salt flux 518 ! using s_i_1d and not sz_i_1d(jk) is ok 519 wfx_res_1d(ji) = wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice ! Mass flux 550 520 551 521 ! update heat content (J.m-2) and layer thickness … … 553 523 h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk) 554 524 555 ELSE !!!Basal melting556 557 zEi = - e_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0)558 zEw = rcp * ztmelts ! Specific enthalpy of meltwater (J/kg, <0)559 zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0)560 561 zfmdt = - zq_bo(ji) / zdE ! Mass flux x time step (kg/m2, >0)562 563 zdeltah(ji,jk) = - zfmdt * r1_rhoic ! Gross thickness change564 565 zdeltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) ) ! bound thickness change525 ELSE !-- Basal melting 526 527 zEi = - e_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0) 528 zEw = rcp * ztmelts ! Specific enthalpy of meltwater (J/kg, <0) 529 zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0) 530 531 zfmdt = - zq_bo(ji) / zdE ! Mass flux x time step (kg/m2, >0) 532 533 zdeltah(ji,jk) = - zfmdt * r1_rhoic ! Gross thickness change 534 535 zdeltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) ) ! bound thickness change 566 536 567 zq_bo(ji) = MAX( 0._wp , zq_bo(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat. MAX is necessary for roundup errors568 569 dh_i_bott(ji) = dh_i_bott(ji) + zdeltah(ji,jk) ! Update basal melt570 571 zfmdt = - zdeltah(ji,jk) * rhoic ! Mass flux x time step > 0572 573 zQm = zfmdt * zEw ! Heat exchanged with ocean574 575 ! Contribution to heat flux to the ocean [W.m-2], <0576 hfx_ thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice577 578 ! Contribution to salt flux (using s_i_1d and not sz_i_1d(jk) is ok)579 sfx_bom_1d(ji) = sfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice537 zq_bo(ji) = MAX( 0._wp , zq_bo(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat. MAX is necessary for roundup errors 538 539 dh_i_bott(ji) = dh_i_bott(ji) + zdeltah(ji,jk) ! Update basal melt 540 541 zfmdt = - zdeltah(ji,jk) * rhoic ! Mass flux x time step > 0 542 543 zQm = zfmdt * zEw ! Heat exchanged with ocean 544 545 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice ! Heat flux to the ocean [W.m-2], <0 546 hfx_bom_1d(ji) = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice ! Heat used in this process [W.m-2], >0 547 548 sfx_bom_1d(ji) = sfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice ! Salt flux 549 ! using s_i_1d and not sz_i_1d(jk) is ok 580 550 581 ! Total heat flux used in this process [W.m-2], >0 582 hfx_bom_1d(ji) = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice 583 584 ! Contribution to mass flux 585 wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 551 wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice ! Mass flux 586 552 587 553 ! update heat content (J.m-2) and layer thickness … … 594 560 END DO 595 561 596 !-------------------------------------------597 562 ! Update temperature, energy 598 ! -------------------------------------------599 DO ji = 1, npti 600 h_i_1d(ji) = 563 ! -------------------------- 564 DO ji = 1, npti 565 h_i_1d(ji) = MAX( 0._wp , h_i_1d(ji) + dh_i_bott(ji) ) 601 566 END DO 602 567 603 !------------------------------------------- 604 ! 5. What to do with remaining energy 605 !------------------------------------------- 606 ! If heat still available for melting and snow remains, then melt more snow 568 ! If heat still available then melt more snow 607 569 !------------------------------------------- 608 570 zdeltah(1:npti,:) = 0._wp ! important 609 571 DO ji = 1, npti 610 zq_rema (ji)= zq_su(ji) + zq_bo(ji)611 rswitch 612 rswitch 613 zdeltah 614 zdeltah 615 dh_s_tot (ji) = dh_s_tot(ji)+ zdeltah(ji,1)616 h_s_1d (ji) = h_s_1d(ji)+ zdeltah(ji,1)572 zq_rema (ji) = zq_su(ji) + zq_bo(ji) 573 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, - h_s_1d(ji) ) ) ! =1 if snow 574 rswitch = rswitch * MAX( 0._wp, SIGN( 1._wp, e_s_1d(ji,1) - epsi20 ) ) 575 zdeltah (ji,1) = - rswitch * zq_rema(ji) / MAX( e_s_1d(ji,1), epsi20 ) 576 zdeltah (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - h_s_1d(ji) ) ) ! bound melting 577 dh_s_tot(ji) = dh_s_tot(ji) + zdeltah(ji,1) 578 h_s_1d (ji) = h_s_1d (ji) + zdeltah(ji,1) 617 579 618 zq_rema(ji) = zq_rema(ji) + zdeltah(ji,1) * e_s_1d(ji,1) ! update available heat (J.m-2) 619 ! heat used to melt snow 620 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * e_s_1d(ji,1) * r1_rdtice ! W.m-2 (>0) 621 ! Contribution to mass flux 622 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice 580 zq_rema(ji) = zq_rema(ji) + zdeltah(ji,1) * e_s_1d(ji,1) ! update available heat (J.m-2) 581 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * e_s_1d(ji,1) * r1_rdtice ! Heat used to melt snow, W.m-2 (>0) 582 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice ! Mass flux 623 583 dh_s_mlt(ji) = dh_s_mlt(ji) + zdeltah(ji,1) 624 584 ! 625 585 ! Remaining heat flux (W.m-2) is sent to the ocean heat budget 626 hfx_out_1d(ji) 586 hfx_out_1d(ji) = hfx_out_1d(ji) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_rdtice 627 587 628 588 IF( ln_icectl .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji) … … 630 590 631 591 ! 632 !------------------------------------------------------------------------------| 633 ! 6) Snow-Ice formation | 634 !------------------------------------------------------------------------------| 592 ! Snow-Ice formation 593 ! ------------------ 635 594 ! When snow load excesses Archimede's limit, snow-ice interface goes down under sea-level, 636 595 ! flooding of seawater transforms snow into ice dh_snowice is positive for the ice … … 648 607 zQm = zfmdt * zEw 649 608 650 ! Contribution to heat flux 651 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 652 653 ! Contribution to salt flux 654 sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_1d(ji) * a_i_1d(ji) * zfmdt * r1_rdtice 609 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice ! Heat flux 610 611 sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_1d(ji) * a_i_1d(ji) * zfmdt * r1_rdtice ! Salt flux 655 612 656 613 ! Case constant salinity in time: virtual salt flux to keep salinity constant … … 660 617 ENDIF 661 618 662 ! Contribution to mass flux 663 ! All snow is thrown in the ocean, and seawater is taken to replace the volume 619 ! Mass flux: All snow is thrown in the ocean, and seawater is taken to replace the volume 664 620 wfx_sni_1d(ji) = wfx_sni_1d(ji) - a_i_1d(ji) * dh_snowice(ji) * rhoic * r1_rdtice 665 621 wfx_snw_sni_1d(ji) = wfx_snw_sni_1d(ji) + a_i_1d(ji) * dh_snowice(ji) * rhosn * r1_rdtice … … 672 628 673 629 ! 674 !-------------------------------------------675 630 ! Update temperature, energy 676 ! -------------------------------------------677 DO ji = 1, npti 678 rswitch = 1.0- MAX( 0._wp , SIGN( 1._wp , - h_i_1d(ji) ) )679 t_su_1d(ji) = rswitch * t_su_1d(ji) + ( 1.0- rswitch ) * rt0631 ! -------------------------- 632 DO ji = 1, npti 633 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - h_i_1d(ji) ) ) 634 t_su_1d(ji) = rswitch * t_su_1d(ji) + ( 1._wp - rswitch ) * rt0 680 635 END DO 681 636 … … 698 653 !!-------------------------------------------------------------------------- 699 654 !! INTERFACE ice_thd_snwblow 655 !! 700 656 !! ** Purpose : Compute distribution of precip over the ice 657 !! 658 !! Snow accumulation in one thermodynamic time step 659 !! snowfall is partitionned between leads and ice. 660 !! If snow fall was uniform, a fraction (1-at_i) would fall into leads 661 !! but because of the winds, more snow falls on leads than on sea ice 662 !! and a greater fraction (1-at_i)^beta of the total mass of snow 663 !! (beta < 1) falls in leads. 664 !! In reality, beta depends on wind speed, 665 !! and should decrease with increasing wind speed but here, it is 666 !! considered as a constant. an average value is 0.66 701 667 !!-------------------------------------------------------------------------- 702 668 !!gm I think it can be usefull to set this as a FUNCTION, not a SUBROUTINE....
Note: See TracChangeset
for help on using the changeset viewer.