Changeset 13959 for NEMO/branches
- Timestamp:
- 2020-12-01T23:47:44+01:00 (3 years ago)
- Location:
- NEMO/branches/2020/SI3_martin_ponds/src/ICE
- Files:
-
- 3 edited
-
icethd.F90 (modified) (1 diff)
-
icethd_dh.F90 (modified) (15 diffs)
-
icethd_zdf_bl99.F90 (modified) (6 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/SI3_martin_ponds/src/ICE/icethd.F90
r13957 r13959 166 166 qsb_ice_bot(ji,jj) = rswitch * MIN( qsb_ice_bot(ji,jj), - zqfr_neg * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) ) 167 167 168 ! If conditions are always supercooled (such as at the mouth of ice-shelves), then ice grows continuously 169 ! ==> stop ice formation by artificially setting up the turbulent fluxes to 0 when volume > 20m (arbitrary) 170 IF( ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ) > 0._wp .AND. vt_i(ji,jj) >= 20._wp ) THEN 171 zqfr = 0._wp 172 zqfr_pos = 0._wp 173 qsb_ice_bot(ji,jj) = 0._wp 174 ENDIF 175 ! 168 176 ! --- Energy Budget of the leads (qlead, J.m-2) --- ! 169 177 ! qlead is the energy received from the atm. in the leads. -
NEMO/branches/2020/SI3_martin_ponds/src/ICE/icethd_dh.F90
r13643 r13959 55 55 !! - Snow ice formation 56 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) 59 !! 57 60 !! References : Bitz and Lipscomb, 1999, J. Geophys. Res. 58 61 !! Fichefet T. and M. Maqueda 1997, J. Geophys. Res., 102(C6), 12609-12646 … … 79 82 REAL(wp) :: zfmdt ! exchange mass flux x time step (J/m2), >0 towards the ocean 80 83 81 REAL(wp), DIMENSION(jpij) :: zqprec ! energy of fallen snow (J.m-3)82 84 REAL(wp), DIMENSION(jpij) :: zq_top ! heat for surface ablation (J.m-2) 83 85 REAL(wp), DIMENSION(jpij) :: zq_bot ! heat for bottom ablation (J.m-2) … … 85 87 REAL(wp), DIMENSION(jpij) :: zf_tt ! Heat budget to determine melting or freezing(W.m-2) 86 88 REAL(wp), DIMENSION(jpij) :: zevap_rema ! remaining mass flux from sublimation (kg.m-2) 87 88 REAL(wp), DIMENSION(jpij) :: zdh_s_mel ! snow melt 89 REAL(wp), DIMENSION(jpij) :: zdh_s_pre ! snow precipitation 90 REAL(wp), DIMENSION(jpij) :: zdh_s_sub ! snow sublimation 91 92 REAL(wp), DIMENSION(jpij,nlay_s) :: zh_s ! snw layer thickness 93 REAL(wp), DIMENSION(jpij,nlay_i) :: zh_i ! ice layer thickness 94 REAL(wp), DIMENSION(jpij,nlay_i) :: zdeltah 95 INTEGER , DIMENSION(jpij,nlay_i) :: icount ! number of layers vanished by melting 96 89 REAL(wp), DIMENSION(jpij) :: zdeltah 97 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) 98 96 99 97 REAL(wp) :: zswitch_sal … … 108 106 END SELECT 109 107 110 ! initialize layer thicknesses and enthalpies 108 ! initialize ice layer thicknesses and enthalpies 109 eh_i_old(1:npti,0:nlay_i+1) = 0._wp 111 110 h_i_old (1:npti,0:nlay_i+1) = 0._wp 112 eh_i_old(1:npti,0:nlay_i+1) = 0._wp111 zh_i (1:npti,0:nlay_i+1) = 0._wp 113 112 DO jk = 1, nlay_i 114 113 DO ji = 1, npti 114 eh_i_old(ji,jk) = h_i_1d(ji) * r1_nlay_i * e_i_1d(ji,jk) 115 115 h_i_old (ji,jk) = h_i_1d(ji) * r1_nlay_i 116 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) 117 127 END DO 118 128 END DO … … 141 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) 142 152 zq_bot(ji) = MAX( 0._wp, zf_tt(ji) * rDt_ice ) 143 END DO144 145 ! Ice and snow layer thicknesses146 !-------------------------------147 DO jk = 1, nlay_i148 DO ji = 1, npti149 zh_i(ji,jk) = h_i_1d(ji) * r1_nlay_i150 END DO151 END DO152 153 DO jk = 1, nlay_s154 DO ji = 1, npti155 zh_s(ji,jk) = h_s_1d(ji) * r1_nlay_s156 END DO157 153 END DO 158 154 … … 167 163 DO ji = 1, npti 168 164 IF( t_s_1d(ji,jk) > rt0 ) THEN 169 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], < 0170 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 171 167 ! updates 172 dh_s_mlt(ji) = dh_s_mlt(ji) - zh_s(ji,jk)173 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) ) 174 170 zh_s (ji,jk) = 0._wp 175 e_s_1d (ji,jk) = 0._wp 176 t_s_1d (ji,jk) = rt0 171 ze_s (ji,jk) = 0._wp 177 172 END IF 178 173 END DO … … 181 176 ! Snow precipitation 182 177 !------------------- 183 CALL ice_var_snwblow( 1.0_wp - at_i_1d(1:npti), zsnw(1:npti) ) ! snow distribution over ice after wind blowing 184 185 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 186 180 DO ji = 1, npti 187 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) 188 184 ! 189 ! --- precipitation --- 190 zdh_s_pre (ji) = zsnw(ji) * sprecip_1d(ji) * rDt_ice * r1_rhos / at_i_1d(ji) ! thickness change 191 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 192 187 ! 193 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)194 wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhos * a_i_1d(ji) * zdh_s_pre(ji) * r1_Dt_ice ! mass flux, <0195 196 ! --- melt of falling snow ---197 rswitch = MAX( 0._wp , SIGN( 1._wp , zqprec(ji) - epsi20 ) )198 zdeltah (ji,1) = - rswitch * zq_top(ji) / MAX( zqprec(ji) , epsi20 ) ! thickness change199 zdeltah (ji,1) = MAX( - zdh_s_pre(ji), zdeltah(ji,1) ) ! bound melting200 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)201 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), >0202 203 ! updates available heat + precipitations after melting204 dh_s_mlt (ji) = dh_s_mlt(ji) + zdeltah(ji,1)205 zq_top (ji) = MAX( 0._wp , zq_top (ji) + zdeltah(ji,1) * zqprec(ji) )206 zdh_s_pre(ji) = zdh_s_pre(ji) + zdeltah(ji,1)207 208 188 ! update thickness 209 h_s_1d(ji) = MAX( 0._wp , h_s_1d(ji) + zdh_s_pre(ji) ) 210 ! 211 ELSE 212 ! 213 zdh_s_pre(ji) = 0._wp 214 zqprec (ji) = 0._wp 215 ! 189 h_s_1d(ji) = h_s_1d(ji) + zh_s(ji,0) 216 190 ENDIF 217 END DO218 219 ! recalculate snow layers220 DO jk = 1, nlay_s221 DO ji = 1, npti222 zh_s(ji,jk) = h_s_1d(ji) * r1_nlay_s223 END DO224 191 END DO 225 192 226 193 ! Snow melting 227 194 ! ------------ 228 ! If heat still available (zq_top > 0), then melt more snow 229 zdeltah(1:npti,:) = 0._wp 230 zdh_s_mel(1:npti) = 0._wp 231 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 232 198 DO ji = 1, npti 233 199 IF( zh_s(ji,jk) > 0._wp .AND. zq_top(ji) > 0._wp ) THEN 234 200 ! 235 rswitch = MAX( 0._wp, SIGN( 1._wp, e_s_1d(ji,jk) - epsi20 ) ) 236 zdeltah (ji,jk) = - rswitch * zq_top(ji) / MAX( e_s_1d(ji,jk), epsi20 ) ! thickness change 237 zdeltah (ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji,jk) ) ! bound melting 238 zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,jk) 239 240 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) 241 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 242 207 243 208 ! updates available heat + thickness 244 dh_s_mlt(ji) = dh_s_mlt(ji) + zdeltah(ji,jk) 245 zq_top (ji) = MAX( 0._wp , zq_top(ji) + zdeltah(ji,jk) * e_s_1d(ji,jk) ) 246 h_s_1d (ji) = MAX( 0._wp , h_s_1d(ji) + zdeltah(ji,jk) ) 247 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 248 214 ! 249 215 ENDIF … … 255 221 ! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates 256 222 ! comment: not counted in mass/heat exchange in iceupdate.F90 since this is an exchange with atm. (not ocean) 257 zdeltah(1:npti,:) = 0._wp 223 zdeltah (1:npti) = 0._wp ! total snow thickness that sublimates, < 0 224 zevap_rema(1:npti) = 0._wp 258 225 DO ji = 1, npti 259 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 260 235 ! 261 zdh_s_sub (ji) = MAX( - h_s_1d(ji) , - evap_ice_1d(ji) * r1_rhos * rDt_ice ) 262 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) 263 zdeltah (ji,1) = MAX( zdh_s_sub(ji), - zdh_s_pre(ji) ) 264 265 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) 266 & ( zdeltah(ji,1) * zqprec(ji) + ( zdh_s_sub(ji) - zdeltah(ji,1) ) * e_s_1d(ji,1) ) & 267 & * a_i_1d(ji) * r1_Dt_ice 268 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 269 270 ! new snow thickness 271 h_s_1d(ji) = MAX( 0._wp , h_s_1d(ji) + zdh_s_sub(ji) ) 272 ! update precipitations after sublimation and correct sublimation 273 zdh_s_pre(ji) = zdh_s_pre(ji) + zdeltah(ji,1) 274 zdh_s_sub(ji) = zdh_s_sub(ji) - zdeltah(ji,1) 275 ! 276 ELSE 277 ! 278 zdh_s_sub (ji) = 0._wp 279 zevap_rema(ji) = 0._wp 280 ! 281 ENDIF 282 END DO 283 284 ! --- Update snow diags --- ! 285 DO ji = 1, npti 286 dh_s_tot(ji) = zdh_s_mel(ji) + zdh_s_pre(ji) + zdh_s_sub(ji) 287 END DO 288 289 ! Update temperature, energy 290 !--------------------------- 291 ! new temp and enthalpy of the snow (remaining snow precip + remaining pre-existing snow) 292 DO jk = 1, nlay_s 293 DO ji = 1,npti 294 rswitch = MAX( 0._wp , SIGN( 1._wp, h_s_1d(ji) - epsi20 ) ) 295 e_s_1d(ji,jk) = rswitch / MAX( h_s_1d(ji), epsi20 ) * & 296 & ( ( zdh_s_pre(ji) ) * zqprec(ji) + & 297 & ( h_s_1d(ji) - zdh_s_pre(ji) ) * rhos * ( rcpi * ( rt0 - t_s_1d(ji,jk) ) + rLfus ) ) 298 END DO 299 END DO 300 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 ! 301 250 ! ! ============ ! 302 251 ! ! Ice ! … … 305 254 ! Surface ice melting 306 255 !-------------------- 307 zdeltah(1:npti,:) = 0._wp ! important308 256 DO jk = 1, nlay_i 309 257 DO ji = 1, npti … … 313 261 314 262 zEi = - e_i_1d(ji,jk) * r1_rhoi ! Specific enthalpy of layer k [J/kg, <0] 315 zdE = 0._wp ! Specific enthalpy difference (J/kg, <0) 316 ! set up at 0 since no energy is needed to melt water...(it is already melted) 317 zdeltah(ji,jk) = MIN( 0._wp , - zh_i(ji,jk) ) ! internal melting occurs when the internal temperature is above freezing 318 ! this should normally not happen, but sometimes, heat diffusion leads to this 319 zfmdt = - zdeltah(ji,jk) * rhoi ! Mass flux x time step > 0 320 321 dh_i_itm(ji) = dh_i_itm(ji) + zdeltah(ji,jk) ! Cumulate internal melting 322 323 zfmdt = - rhoi * zdeltah(ji,jk) ! Recompute mass flux [kg/m2, >0] 324 325 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 326 ! ice enthalpy zEi is "sent" to the ocean 327 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 328 ! using s_i_1d and not sz_i_1d(jk) is ok 329 wfx_res_1d(ji) = wfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice ! Mass flux 330 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 331 276 ELSE !-- Surface melting 332 277 … … 337 282 zfmdt = - zq_top(ji) / zdE ! Mass flux to the ocean [kg/m2, >0] 338 283 339 zd eltah(ji,jk)= - zfmdt * r1_rhoi ! Melt of layer jk [m, <0]340 341 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]342 343 zq_top(ji) = MAX( 0._wp , zq_top(ji) - zdeltah(ji,jk)* rhoi * zdE ) ! update available heat344 345 dh_i_sum(ji) = dh_i_sum(ji) + zd eltah(ji,jk)! Cumulate surface melt346 347 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] 348 293 349 294 zQm = zfmdt * zEw ! Energy of the melt water sent to the ocean [J/m2, <0] 350 295 351 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 352 ! using s_i_1d and not sz_i_1d(jk) is ok) 353 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_Dt_ice ! Heat flux [W.m-2], < 0 354 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 355 ! 356 wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice ! Mass flux 357 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) 358 301 END IF 359 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 ! 360 311 ! Ice sublimation 361 312 ! --------------- 362 zdum = MAX( - ( zh_i(ji,jk) + zdeltah(ji,jk) ) , - zevap_rema(ji) * r1_rhoi ) 363 zdeltah (ji,jk) = zdeltah (ji,jk) + zdum 364 dh_i_sub(ji) = dh_i_sub(ji) + zdum 365 366 sfx_sub_1d(ji) = sfx_sub_1d(ji) - rhoi * a_i_1d(ji) * zdum * s_i_1d(ji) * r1_Dt_ice ! Salt flux >0 367 ! clem: flux is sent to the ocean for simplicity 368 ! but salt should remain in the ice except 369 ! if all ice is melted. => must be corrected 370 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 371 372 wfx_ice_sub_1d(ji) = wfx_ice_sub_1d(ji) - rhoi * a_i_1d(ji) * zdum * r1_Dt_ice ! Mass flux > 0 373 374 ! update remaining mass flux 375 zevap_rema(ji) = zevap_rema(ji) + zdum * rhoi 376 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 377 331 ! record which layers have disappeared (for bottom melting) 378 332 ! => icount=0 : no layer has vanished 379 333 ! => icount=5 : 5 layers have vanished 380 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) ) ) 381 335 icount(ji,jk) = NINT( rswitch ) 382 zh_i(ji,jk) = MAX( 0._wp , zh_i(ji,jk) + zdeltah(ji,jk) )383 336 384 ! update heat content (J.m-2) and layer thickness 385 eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk) 386 h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk) 387 END DO 388 END DO 389 390 ! update ice thickness 391 DO ji = 1, npti 392 h_i_1d(ji) = MAX( 0._wp , h_i_1d(ji) + dh_i_sum(ji) + dh_i_itm(ji) + dh_i_sub(ji) ) 393 END DO 394 337 END DO 338 END DO 339 395 340 ! remaining "potential" evap is sent to ocean 396 341 DO ji = 1, npti … … 430 375 & + zswi2 * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) ) , 0.5 ) 431 376 432 s_i_new(ji) = zswitch_sal * zfracs * sss_1d(ji) + ( 1. - zswitch_sal ) * s_i_1d(ji) ! New ice salinity433 434 ztmelts = - rTmlt * s_i_new(ji) ! New ice melting point (C)435 436 zt_i_new = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i)437 438 zEi = rcpi * ( zt_i_new - (ztmelts+rt0) ) & ! Specific enthalpy of forming ice (J/kg, <0)439 & - rLfus * ( 1.0 - ztmelts / ( MIN( zt_i_new - rt0, -epsi10 ) ) ) + rcp* ztmelts440 441 zEw = rcp * ( t_bo_1d(ji) - rt0 ) ! Specific enthalpy of seawater (J/kg, < 0)442 443 zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0)444 445 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 ) ) 446 391 447 392 END DO 448 393 ! Contribution to Energy and Salt Fluxes 449 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) 450 395 451 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 452 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 453 454 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 455 456 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) 457 404 458 405 ! update heat content (J.m-2) and layer thickness … … 466 413 ! Ice Basal melt 467 414 !--------------- 468 zdeltah(1:npti,:) = 0._wp ! important469 415 DO jk = nlay_i, 1, -1 470 416 DO ji = 1, npti … … 475 421 IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN !-- Internal melting 476 422 477 zEi = - e_i_1d(ji,jk) * r1_rhoi ! Specific enthalpy of melting ice (J/kg, <0) 478 zdE = 0._wp ! Specific enthalpy difference (J/kg, <0) 479 ! set up at 0 since no energy is needed to melt water...(it is already melted) 480 zdeltah (ji,jk) = MIN( 0._wp , - zh_i(ji,jk) ) ! internal melting occurs when the internal temperature is above freezing 481 ! this should normally not happen, but sometimes, heat diffusion leads to this 482 483 dh_i_itm (ji) = dh_i_itm(ji) + zdeltah(ji,jk) 484 485 zfmdt = - zdeltah(ji,jk) * rhoi ! Mass flux x time step > 0 486 487 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 488 ! ice enthalpy zEi is "sent" to the ocean 489 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 490 ! using s_i_1d and not sz_i_1d(jk) is ok 491 wfx_res_1d(ji) = wfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice ! Mass flux 492 493 ! update heat content (J.m-2) and layer thickness 494 eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk) 495 h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk) 496 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 497 437 ELSE !-- Basal melting 498 438 499 zEi = - e_i_1d(ji,jk) * r1_rhoi! Specific enthalpy of melting ice (J/kg, <0)500 zEw = rcp * ztmelts! Specific enthalpy of meltwater (J/kg, <0)501 zdE = zEi - zEw! Specific enthalpy difference (J/kg, <0)502 503 zfmdt = - zq_bot(ji) / zdE! Mass flux x time step (kg/m2, >0)504 505 zd eltah(ji,jk) = - zfmdt * r1_rhoi! Gross thickness change506 507 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 508 448 509 zq_bot(ji) = MAX( 0._wp , zq_bot(ji) - zdeltah(ji,jk) * rhoi * zdE ) ! update available heat. MAX is necessary for roundup errors 510 511 dh_i_bom(ji) = dh_i_bom(ji) + zdeltah(ji,jk) ! Update basal melt 512 513 zfmdt = - zdeltah(ji,jk) * rhoi ! Mass flux x time step > 0 514 515 zQm = zfmdt * zEw ! Heat exchanged with ocean 516 517 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 518 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 519 520 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 521 ! using s_i_1d and not sz_i_1d(jk) is ok 522 523 wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice ! Mass flux 524 525 ! update heat content (J.m-2) and layer thickness 526 eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk) 527 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 528 462 ENDIF 529 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 530 470 ENDIF 531 471 END DO 532 472 END DO 533 473 534 ! Update temperature, energy 535 ! -------------------------- 536 DO ji = 1, npti 537 h_i_1d(ji) = MAX( 0._wp , h_i_1d(ji) + dh_i_bog(ji) + dh_i_bom(ji) ) 538 END DO 539 540 ! If heat still available then melt more snow 541 !------------------------------------------- 542 zdeltah(1:npti,:) = 0._wp ! important 543 DO ji = 1, npti 544 zq_rema (ji) = zq_top(ji) + zq_bot(ji) 545 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, - h_s_1d(ji) ) ) ! =1 if snow 546 rswitch = rswitch * MAX( 0._wp, SIGN( 1._wp, e_s_1d(ji,1) - epsi20 ) ) 547 zdeltah (ji,1) = - rswitch * zq_rema(ji) / MAX( e_s_1d(ji,1), epsi20 ) 548 zdeltah (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - h_s_1d(ji) ) ) ! bound melting 549 dh_s_tot(ji) = dh_s_tot(ji) + zdeltah(ji,1) 550 h_s_1d (ji) = h_s_1d (ji) + zdeltah(ji,1) 551 552 zq_rema(ji) = zq_rema(ji) + zdeltah(ji,1) * e_s_1d(ji,1) ! update available heat (J.m-2) 553 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) 554 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhos * a_i_1d(ji) * zdeltah(ji,1) * r1_Dt_ice ! Mass flux 555 dh_s_mlt(ji) = dh_s_mlt(ji) + zdeltah(ji,1) 556 ! 557 ! Remaining heat flux (W.m-2) is sent to the ocean heat budget 558 !!hfx_res_1d(ji) = hfx_res_1d(ji) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_Dt_ice 559 560 IF( ln_icectl .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji) 561 END DO 562 563 ! 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 564 514 ! Snow-Ice formation 565 515 ! ------------------ 566 ! When snow load exce sses Archimede's limit, snow-ice interface goes down under sea-level,567 ! 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) 568 518 z1_rho = 1._wp / ( rhos+rho0-rhoi ) 519 zdeltah(1:npti) = 0._wp 569 520 DO ji = 1, npti 570 521 ! 571 dh_snowice(ji) = MAX( 522 dh_snowice(ji) = MAX( 0._wp , ( rhos * h_s_1d(ji) + (rhoi-rho0) * h_i_1d(ji) ) * z1_rho ) 572 523 573 524 h_i_1d(ji) = h_i_1d(ji) + dh_snowice(ji) … … 579 530 zQm = zfmdt * zEw 580 531 581 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_Dt_ice ! Heat flux 582 583 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 584 534 585 535 ! Case constant salinity in time: virtual salt flux to keep salinity constant 586 536 IF( nn_icesal /= 2 ) THEN 587 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 ocean588 & - 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 589 539 ENDIF 590 540 591 541 ! Mass flux: All snow is thrown in the ocean, and seawater is taken to replace the volume 592 wfx_sni_1d(ji) = wfx_sni_1d(ji) - a_i_1d(ji) * dh_snowice(ji) * rhoi * r1_Dt_ice 593 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) 594 548 595 549 ! update heat content (J.m-2) and layer thickness 596 eh_i_old(ji,0) = eh_i_old(ji,0) + dh_snowice(ji) * e_s_1d(ji,1) + zfmdt * zEw597 550 h_i_old (ji,0) = h_i_old (ji,0) + dh_snowice(ji) 598 599 END DO 600 601 ! 602 ! Update temperature, energy 603 ! -------------------------- 604 DO ji = 1, npti 605 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - h_i_1d(ji) ) ) 606 t_su_1d(ji) = rswitch * t_su_1d(ji) + ( 1._wp - rswitch ) * rt0 607 END DO 608 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 609 578 DO jk = 1, nlay_s 610 579 DO ji = 1,npti 611 ! where there is no ice or no snow 612 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) ) ) ) 613 ! mass & energy loss to the ocean 614 hfx_res_1d(ji) = hfx_res_1d(ji) + ( 1._wp - rswitch ) * & 615 & ( 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 616 wfx_res_1d(ji) = wfx_res_1d(ji) + ( 1._wp - rswitch ) * & 617 & ( rhos * h_s_1d(ji) * r1_nlay_s * a_i_1d(ji) * r1_Dt_ice ) ! mass flux 618 ! update energy (mass is updated in the next loop) 619 e_s_1d(ji,jk) = rswitch * e_s_1d(ji,jk) 620 ! recalculate t_s_1d from e_s_1d 621 t_s_1d(ji,jk) = rt0 + rswitch * ( - e_s_1d(ji,jk) * r1_rhos * r1_rcpi + rLfus * r1_rcpi ) 622 END DO 623 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 624 589 625 590 ! --- ensure that a_i = 0 & h_s = 0 where h_i = 0 --- 626 591 WHERE( h_i_1d(1:npti) == 0._wp ) 627 a_i_1d(1:npti) = 0._wp 628 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 629 595 END WHERE 630 !596 631 597 END SUBROUTINE ice_thd_dh 632 598 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 633 691 #else 634 692 !!---------------------------------------------------------------------- -
NEMO/branches/2020/SI3_martin_ponds/src/ICE/icethd_zdf_bl99.F90
r13472 r13959 109 109 REAL(wp), DIMENSION(jpij) :: zdqns_ice_b ! derivative of the surface flux function 110 110 ! 111 REAL(wp), DIMENSION(jpij ) 112 REAL(wp), DIMENSION(jpij,nlay_i) 113 REAL(wp), DIMENSION(jpij,nlay_s) 114 REAL(wp), DIMENSION(jpij,nlay_i) 115 REAL(wp), DIMENSION(jpij,nlay_s) 116 REAL(wp), DIMENSION(jpij,0:nlay_i) 117 REAL(wp), DIMENSION(jpij,0:nlay_i) 118 REAL(wp), DIMENSION(jpij,0:nlay_i) 119 REAL(wp), DIMENSION(jpij,0:nlay_i) 120 REAL(wp), DIMENSION(jpij,0:nlay_i) 121 REAL(wp), DIMENSION(jpij,0:nlay_i) 122 REAL(wp), DIMENSION(jpij,0:nlay_s) 123 REAL(wp), DIMENSION(jpij,0:nlay_s) 124 REAL(wp), DIMENSION(jpij,0:nlay_s) 125 REAL(wp), DIMENSION(jpij,0:nlay_s) 126 REAL(wp), DIMENSION(jpij) 127 REAL(wp), DIMENSION(jpij ,nlay_i+3) :: zindterm ! 'Ind'ependent term128 REAL(wp), DIMENSION(jpij ,nlay_i+3) :: zindtbis ! Temporary 'ind'ependent term129 REAL(wp), DIMENSION(jpij ,nlay_i+3) :: zdiagbis ! Temporary 'dia'gonal term130 REAL(wp), DIMENSION(jpij ,nlay_i+3,3) :: ztrid ! Tridiagonal system terms131 REAL(wp), DIMENSION(jpij) :: zq_ini ! diag errors on heat132 REAL(wp), DIMENSION(jpij ) :: zghe ! G(he), th. conduct enhancement factor, mono-cat133 REAL(wp), DIMENSION(jpij ) :: za_s_fra ! ice fraction covered by snow134 REAL(wp), DIMENSION(jpij ) :: isnow ! snow presence (1) or not (0)135 REAL(wp), DIMENSION(jpij ) :: isnow_comb ! snow presence for met-office111 REAL(wp), DIMENSION(jpij ) :: ztsuold ! Old surface temperature in the ice 112 REAL(wp), DIMENSION(jpij,nlay_i) :: ztiold ! Old temperature in the ice 113 REAL(wp), DIMENSION(jpij,nlay_s) :: ztsold ! Old temperature in the snow 114 REAL(wp), DIMENSION(jpij,nlay_i) :: ztib ! Temporary temperature in the ice to check the convergence 115 REAL(wp), DIMENSION(jpij,nlay_s) :: ztsb ! Temporary temperature in the snow to check the convergence 116 REAL(wp), DIMENSION(jpij,0:nlay_i) :: ztcond_i ! Ice thermal conductivity 117 REAL(wp), DIMENSION(jpij,0:nlay_i) :: ztcond_i_cp ! copy 118 REAL(wp), DIMENSION(jpij,0:nlay_i) :: zradtr_i ! Radiation transmitted through the ice 119 REAL(wp), DIMENSION(jpij,0:nlay_i) :: zradab_i ! Radiation absorbed in the ice 120 REAL(wp), DIMENSION(jpij,0:nlay_i) :: zkappa_i ! Kappa factor in the ice 121 REAL(wp), DIMENSION(jpij,0:nlay_i) :: zeta_i ! Eta factor in the ice 122 REAL(wp), DIMENSION(jpij,0:nlay_s) :: zradtr_s ! Radiation transmited through the snow 123 REAL(wp), DIMENSION(jpij,0:nlay_s) :: zradab_s ! Radiation absorbed in the snow 124 REAL(wp), DIMENSION(jpij,0:nlay_s) :: zkappa_s ! Kappa factor in the snow 125 REAL(wp), DIMENSION(jpij,0:nlay_s) :: zeta_s ! Eta factor in the snow 126 REAL(wp), DIMENSION(jpij) :: zkappa_comb ! Combined snow and ice surface conductivity 127 REAL(wp), DIMENSION(jpij) :: zq_ini ! diag errors on heat 128 REAL(wp), DIMENSION(jpij) :: zghe ! G(he), th. conduct enhancement factor, mono-cat 129 REAL(wp), DIMENSION(jpij) :: za_s_fra ! ice fraction covered by snow 130 REAL(wp), DIMENSION(jpij) :: isnow ! snow presence (1) or not (0) 131 REAL(wp), DIMENSION(jpij) :: isnow_comb ! snow presence for met-office 132 REAL(wp), DIMENSION(jpij,nlay_i+nlay_s+1) :: zindterm ! 'Ind'ependent term 133 REAL(wp), DIMENSION(jpij,nlay_i+nlay_s+1) :: zindtbis ! Temporary 'ind'ependent term 134 REAL(wp), DIMENSION(jpij,nlay_i+nlay_s+1) :: zdiagbis ! Temporary 'dia'gonal term 135 REAL(wp), DIMENSION(jpij,nlay_i+nlay_s+1,3) :: ztrid ! Tridiagonal system terms 136 136 ! 137 137 ! Mono-category … … 533 533 ! Solve the tridiagonal system with Gauss elimination method. 534 534 ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984 535 jm_maxt = 0 536 jm_mint = nlay_i+5 537 DO ji = 1, npti 538 jm_mint = MIN(jm_min(ji),jm_mint) 539 jm_maxt = MAX(jm_max(ji),jm_maxt) 540 END DO 541 542 DO jk = jm_mint+1, jm_maxt 543 DO ji = 1, npti 544 jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji)) 535 !!$ jm_maxt = 0 536 !!$ jm_mint = nlay_i+5 537 !!$ DO ji = 1, npti 538 !!$ jm_mint = MIN(jm_min(ji),jm_mint) 539 !!$ jm_maxt = MAX(jm_max(ji),jm_maxt) 540 !!$ END DO 541 !!$ !!clem SNWLAY => check why LIM1D does not get this loop. Is nlay_i+5 correct? 542 !!$ 543 !!$ DO jk = jm_mint+1, jm_maxt 544 !!$ DO ji = 1, npti 545 !!$ jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji)) 546 !!$ zdiagbis(ji,jm) = ztrid (ji,jm,2) - ztrid(ji,jm,1) * ztrid (ji,jm-1,3) / zdiagbis(ji,jm-1) 547 !!$ zindtbis(ji,jm) = zindterm(ji,jm ) - ztrid(ji,jm,1) * zindtbis(ji,jm-1 ) / zdiagbis(ji,jm-1) 548 !!$ END DO 549 !!$ END DO 550 ! clem: maybe one should find a way to reverse this loop for mpi performance 551 DO ji = 1, npti 552 jm_mint = jm_min(ji) 553 jm_maxt = jm_max(ji) 554 DO jm = jm_mint+1, jm_maxt 545 555 zdiagbis(ji,jm) = ztrid (ji,jm,2) - ztrid(ji,jm,1) * ztrid (ji,jm-1,3) / zdiagbis(ji,jm-1) 546 556 zindtbis(ji,jm) = zindterm(ji,jm ) - ztrid(ji,jm,1) * zindtbis(ji,jm-1 ) / zdiagbis(ji,jm-1) … … 564 574 END DO 565 575 576 ! snow temperatures 566 577 DO ji = 1, npti 567 578 ! Variables used after iterations 568 579 ! Value must be frozen after convergence for MPP independance reason 569 IF ( .NOT. l_T_converged(ji) ) THEN 570 ! snow temperatures 571 IF( h_s_1d(ji) > 0._wp ) THEN 572 t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1) 573 ENDIF 574 ! surface temperature 580 IF ( .NOT. l_T_converged(ji) .AND. h_s_1d(ji) > 0._wp ) & 581 & t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1) 582 END DO 583 !!clem SNWLAY 584 DO jm = nlay_s, 2, -1 585 DO ji = 1, npti 586 jk = jm - 1 587 IF ( .NOT. l_T_converged(ji) .AND. h_s_1d(ji) > 0._wp ) & 588 & t_s_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_s_1d(ji,jk+1) ) / zdiagbis(ji,jm) 589 END DO 590 END DO 591 592 ! surface temperature 593 DO ji = 1, npti 594 IF( .NOT. l_T_converged(ji) ) THEN 575 595 ztsub(ji) = t_su_1d(ji) 576 596 IF( t_su_1d(ji) < rt0 ) THEN 577 t_su_1d(ji) = ( 578 & ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) *t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji))597 t_su_1d(ji) = ( zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) * & 598 & ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji)) 579 599 ENDIF 580 600 ENDIF 581 601 END DO 582 !clem: in order to have several layers of snow, there is a missing loop here for t_s_1d(1:nlay_s-1)583 602 ! 584 603 !-------------------------------------------------------------- … … 727 746 ! Solve the tridiagonal system with Gauss elimination method. 728 747 ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984 729 jm_maxt = 0 730 jm_mint = nlay_i+5 731 DO ji = 1, npti 732 jm_mint = MIN(jm_min(ji),jm_mint) 733 jm_maxt = MAX(jm_max(ji),jm_maxt) 734 END DO 735 736 DO jk = jm_mint+1, jm_maxt 737 DO ji = 1, npti 738 jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji)) 748 !!$ jm_maxt = 0 749 !!$ jm_mint = nlay_i+5 750 !!$ DO ji = 1, npti 751 !!$ jm_mint = MIN(jm_min(ji),jm_mint) 752 !!$ jm_maxt = MAX(jm_max(ji),jm_maxt) 753 !!$ END DO 754 !!$ 755 !!$ DO jk = jm_mint+1, jm_maxt 756 !!$ DO ji = 1, npti 757 !!$ jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji)) 758 !!$ zdiagbis(ji,jm) = ztrid (ji,jm,2) - ztrid(ji,jm,1) * ztrid (ji,jm-1,3) / zdiagbis(ji,jm-1) 759 !!$ zindtbis(ji,jm) = zindterm(ji,jm) - ztrid(ji,jm,1) * zindtbis(ji,jm-1) / zdiagbis(ji,jm-1) 760 !!$ END DO 761 !!$ END DO 762 ! clem: maybe one should find a way to reverse this loop for mpi performance 763 DO ji = 1, npti 764 jm_mint = jm_min(ji) 765 jm_maxt = jm_max(ji) 766 DO jm = jm_mint+1, jm_maxt 739 767 zdiagbis(ji,jm) = ztrid (ji,jm,2) - ztrid(ji,jm,1) * ztrid (ji,jm-1,3) / zdiagbis(ji,jm-1) 740 zindtbis(ji,jm) = zindterm(ji,jm ) - ztrid(ji,jm,1) * zindtbis(ji,jm-1)/ zdiagbis(ji,jm-1)768 zindtbis(ji,jm) = zindterm(ji,jm ) - ztrid(ji,jm,1) * zindtbis(ji,jm-1 ) / zdiagbis(ji,jm-1) 741 769 END DO 742 770 END DO 743 771 744 772 ! ice temperatures 745 773 DO ji = 1, npti … … 761 789 ! snow temperatures 762 790 DO ji = 1, npti 763 ! Variable used after iterations791 ! Variables used after iterations 764 792 ! Value must be frozen after convergence for MPP independance reason 765 IF ( .NOT. l_T_converged(ji) ) THEN 766 IF( h_s_1d(ji) > 0._wp ) THEN 767 t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1) 768 ENDIF 769 ENDIF 770 END DO 771 !clem: in order to have several layers of snow, there is a missing loop here for t_s_1d(1:nlay_s-1) 793 IF ( .NOT. l_T_converged(ji) .AND. h_s_1d(ji) > 0._wp ) & 794 & t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1) 795 END DO 796 !!clem SNWLAY 797 DO jm = nlay_s, 2, -1 798 DO ji = 1, npti 799 jk = jm - 1 800 IF ( .NOT. l_T_converged(ji) .AND. h_s_1d(ji) > 0._wp ) & 801 & t_s_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_s_1d(ji,jk+1) ) / zdiagbis(ji,jm) 802 END DO 803 END DO 772 804 ! 773 805 !-------------------------------------------------------------- … … 923 955 !--- Snow-ice interfacial temperature (diagnostic SIMIP) 924 956 IF( h_s_1d(ji) >= zhs_ssl ) THEN 925 t_si_1d(ji) = ( rn_cnd_s * h_i_1d(ji) * r1_nlay_i * t_s_1d(ji, 1) &926 & + ztcond_i(ji,1) * h_s_1d(ji) * r1_nlay_s * t_i_1d(ji,1) ) &957 t_si_1d(ji) = ( rn_cnd_s * h_i_1d(ji) * r1_nlay_i * t_s_1d(ji,nlay_s) & 958 & + ztcond_i(ji,1) * h_s_1d(ji) * r1_nlay_s * t_i_1d(ji,1) ) & 927 959 & / ( rn_cnd_s * h_i_1d(ji) * r1_nlay_i & 928 960 & + ztcond_i(ji,1) * h_s_1d(ji) * r1_nlay_s )
Note: See TracChangeset
for help on using the changeset viewer.
