Changeset 15374
- Timestamp:
- 2021-10-14T19:49:16+02:00 (3 years ago)
- File:
-
- 1 edited
-
NEMO/trunk/src/ICE/icethd.F90 (modified) (6 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/ICE/icethd.F90
r15334 r15374 92 92 ! 93 93 INTEGER :: ji, jj, jk, jl ! dummy loop indices 94 REAL(wp) :: zfric_u, zqld, zqfr, zqfr_neg, zqfr_pos95 REAL(wp), PARAMETER :: zfric_umin = 0._wp ! lower bound for the friction velocity (cice value=5.e-04)96 REAL(wp), PARAMETER :: zch = 0.0057_wp ! heat transfer coefficient97 REAL(wp), DIMENSION(jpi,jpj) :: zu_io, zv_io, zfric, zvel ! ice-ocean velocity (m/s) and frictional velocity (m2/s2)98 !99 ! for collection thickness100 INTEGER :: iter ! - -101 REAL(wp) :: zvfrx, zvgx, ztaux, zf ! - -102 REAL(wp) :: ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, ztwogp ! - -103 REAL(wp), PARAMETER :: zcai = 1.4e-3_wp ! ice-air drag (clem: should be dependent on coupling/forcing used)104 REAL(wp), PARAMETER :: zhicrit = 0.04_wp ! frazil ice thickness105 REAL(wp), PARAMETER :: zsqcd = 1.0_wp / SQRT( 1.3_wp * zcai ) ! 1/SQRT(airdensity*drag)106 REAL(wp), PARAMETER :: zgamafr = 0.03_wp107 94 !!------------------------------------------------------------------- 108 95 ! controls … … 123 110 ENDIF 124 111 125 !---------------------------------------------! 112 ! --- calculate (some) heat fluxes and frazil ice --- ! 113 ! out => sensible heat (qsb_ice_bot) & heat budget in the leads (fhld & qlead) 114 ! out => ice collection thickness (ht_i_new) and fraction of frazil (fraz_frac) 115 CALL thd_prep 116 117 !-------------------------------------------------------------------------------------------! 118 ! Thermodynamic computation (only on grid points covered by ice) => loop over ice categories 119 !-------------------------------------------------------------------------------------------! 120 DO jl = 1, jpl 121 122 ! select ice covered grid points 123 npti = 0 ; nptidx(:) = 0 124 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 125 IF ( a_i(ji,jj,jl) > epsi10 ) THEN 126 npti = npti + 1 127 nptidx(npti) = (jj - 1) * jpi + ji 128 ENDIF 129 END_2D 130 131 IF( npti > 0 ) THEN ! If there is no ice, do nothing. 132 ! 133 CALL ice_thd_1d2d( jl, 1 ) ! --- Move to 1D arrays --- ! 134 ! ! --- & Change units of e_i, e_s from J/m2 to J/m3 --- ! 135 ! 136 s_i_new (1:npti) = 0._wp ; dh_s_tot(1:npti) = 0._wp ! --- some init --- ! (important to have them here) 137 dh_i_sum (1:npti) = 0._wp ; dh_i_bom(1:npti) = 0._wp ; dh_i_itm (1:npti) = 0._wp 138 dh_i_sub (1:npti) = 0._wp ; dh_i_bog(1:npti) = 0._wp 139 dh_snowice(1:npti) = 0._wp ; dh_s_mlt(1:npti) = 0._wp 140 ! 141 CALL ice_thd_zdf ! --- Ice-Snow temperature --- ! 142 ! 143 IF( ln_icedH ) THEN ! --- Growing/Melting --- ! 144 CALL ice_thd_dh ! Ice-Snow thickness 145 CALL ice_thd_ent( e_i_1d(1:npti,:) ) ! Ice enthalpy remapping 146 ENDIF 147 CALL ice_thd_sal( ln_icedS ) ! --- Ice salinity --- ! 148 ! 149 CALL ice_thd_temp ! --- Temperature update --- ! 150 ! 151 IF( ln_icedH .AND. ln_virtual_itd ) & 152 & CALL ice_thd_mono ! --- Extra lateral melting if virtual_itd --- ! 153 ! 154 IF( ln_icedA ) CALL ice_thd_da ! --- Lateral melting --- ! 155 ! 156 CALL ice_thd_1d2d( jl, 2 ) ! --- Change units of e_i, e_s from J/m3 to J/m2 --- ! 157 ! ! --- & Move to 2D arrays --- ! 158 ENDIF 159 ! 160 END DO 161 ! 162 IF( ln_icediachk ) CALL ice_cons_hsm(1, 'icethd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 163 IF( ln_icediachk ) CALL ice_cons2D (1, 'icethd', diag_v, diag_s, diag_t, diag_fv, diag_fs, diag_ft) 164 ! 165 IF ( ln_pnd .AND. ln_icedH ) & 166 & CALL ice_thd_pnd ! --- Melt ponds --- ! 167 ! 168 IF( jpl > 1 ) CALL ice_itd_rem( kt ) ! --- Transport ice between thickness categories --- ! 169 ! 170 IF( ln_icedO ) CALL ice_thd_do ! --- Frazil ice growth in leads --- ! 171 ! 172 CALL ice_cor( kt , 2 ) ! --- Corrections --- ! 173 ! 174 oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * rDt_ice ! --- Ice natural aging incrementation 175 ! 176 DO_2D( 0, 0, 0, 0 ) ! --- Ice velocity corrections 177 IF( at_i(ji,jj) == 0._wp ) THEN ! if ice has melted 178 IF( at_i(ji+1,jj) == 0._wp ) u_ice(ji ,jj) = 0._wp ! right side 179 IF( at_i(ji-1,jj) == 0._wp ) u_ice(ji-1,jj) = 0._wp ! left side 180 IF( at_i(ji,jj+1) == 0._wp ) v_ice(ji,jj ) = 0._wp ! upper side 181 IF( at_i(ji,jj-1) == 0._wp ) v_ice(ji,jj-1) = 0._wp ! bottom side 182 ENDIF 183 END_2D 184 CALL lbc_lnk( 'icecor', u_ice, 'U', -1.0_wp, v_ice, 'V', -1.0_wp ) 185 ! 186 ! convergence tests 187 IF( ln_zdf_chkcvg ) THEN 188 CALL iom_put( 'tice_cvgerr', ztice_cvgerr ) ; DEALLOCATE( ztice_cvgerr ) 189 CALL iom_put( 'tice_cvgstp', ztice_cvgstp ) ; DEALLOCATE( ztice_cvgstp ) 190 ENDIF 191 ! 192 ! controls 193 IF( ln_icectl ) CALL ice_prt (kt, iiceprt, jiceprt, 1, ' - ice thermodyn. - ') ! prints 194 IF( sn_cfctl%l_prtctl ) & 195 & CALL ice_prt3D ('icethd') ! prints 196 IF( ln_timing ) CALL timing_stop('icethd') ! timing 197 ! 198 END SUBROUTINE ice_thd 199 200 201 SUBROUTINE ice_thd_temp 202 !!----------------------------------------------------------------------- 203 !! *** ROUTINE ice_thd_temp *** 204 !! 205 !! ** Purpose : Computes sea ice temperature (Kelvin) from enthalpy 206 !! 207 !! ** Method : Formula (Bitz and Lipscomb, 1999) 208 !!------------------------------------------------------------------- 209 INTEGER :: ji, jk ! dummy loop indices 210 REAL(wp) :: ztmelts, zbbb, zccc ! local scalar 211 !!------------------------------------------------------------------- 212 ! Recover ice temperature 213 DO jk = 1, nlay_i 214 DO ji = 1, npti 215 ztmelts = -rTmlt * sz_i_1d(ji,jk) 216 ! Conversion q(S,T) -> T (second order equation) 217 zbbb = ( rcp - rcpi ) * ztmelts + e_i_1d(ji,jk) * r1_rhoi - rLfus 218 zccc = SQRT( MAX( zbbb * zbbb - 4._wp * rcpi * rLfus * ztmelts, 0._wp ) ) 219 t_i_1d(ji,jk) = rt0 - ( zbbb + zccc ) * 0.5_wp * r1_rcpi 220 221 ! mask temperature 222 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - h_i_1d(ji) ) ) 223 t_i_1d(ji,jk) = rswitch * t_i_1d(ji,jk) + ( 1._wp - rswitch ) * rt0 224 END DO 225 END DO 226 ! 227 END SUBROUTINE ice_thd_temp 228 229 230 SUBROUTINE ice_thd_mono 231 !!----------------------------------------------------------------------- 232 !! *** ROUTINE ice_thd_mono *** 233 !! 234 !! ** Purpose : Lateral melting in case virtual_itd 235 !! ( dA = A/2h dh ) 236 !!----------------------------------------------------------------------- 237 INTEGER :: ji ! dummy loop indices 238 REAL(wp) :: zhi_bef ! ice thickness before thermo 239 REAL(wp) :: zdh_mel, zda_mel ! net melting 240 REAL(wp) :: zvi, zvs ! ice/snow volumes 241 !!----------------------------------------------------------------------- 242 ! 243 DO ji = 1, npti 244 zdh_mel = MIN( 0._wp, dh_i_itm(ji) + dh_i_sum(ji) + dh_i_bom(ji) + dh_snowice(ji) + dh_i_sub(ji) ) 245 IF( zdh_mel < 0._wp .AND. a_i_1d(ji) > 0._wp ) THEN 246 zvi = a_i_1d(ji) * h_i_1d(ji) 247 zvs = a_i_1d(ji) * h_s_1d(ji) 248 ! lateral melting = concentration change 249 zhi_bef = h_i_1d(ji) - zdh_mel 250 rswitch = MAX( 0._wp , SIGN( 1._wp , zhi_bef - epsi20 ) ) 251 zda_mel = rswitch * a_i_1d(ji) * zdh_mel / ( 2._wp * MAX( zhi_bef, epsi20 ) ) 252 a_i_1d(ji) = MAX( epsi20, a_i_1d(ji) + zda_mel ) 253 ! adjust thickness 254 h_i_1d(ji) = zvi / a_i_1d(ji) 255 h_s_1d(ji) = zvs / a_i_1d(ji) 256 ! retrieve total concentration 257 at_i_1d(ji) = a_i_1d(ji) 258 END IF 259 END DO 260 ! 261 END SUBROUTINE ice_thd_mono 262 263 SUBROUTINE thd_prep 264 !!----------------------------------------------------------------------- 265 !! *** ROUTINE thd_prep *** 266 !! 267 !! ** Purpose : prepare necessary fields for thermo calculations 268 !! 269 !! For the fluxes 270 !! ** Inputs : u_ice, v_ice, ssu_m, ssv_m, utau, vtau 271 !! frq_m, qsr_oce, qns_oce, qemp_oce, e3t_m, sst_m 272 !! ** Outputs : qsb_ice_bot, fhld, qlead 273 !! 274 !! For the collection thickness (frazil) 275 !! ** Inputs : u_ice, v_ice, utau_ice, vtau_ice 276 !! ** Ouputs : ht_i_new, fraz_frac 277 !!----------------------------------------------------------------------- 278 INTEGER :: ji, jj ! dummy loop indices 279 REAL(wp) :: zfric_u, zqld, zqfr, zqfr_neg, zqfr_pos, zu_io, zv_io, zu_iom1, zv_iom1 280 REAL(wp), PARAMETER :: zfric_umin = 0._wp ! lower bound for the friction velocity (cice value=5.e-04) 281 REAL(wp), PARAMETER :: zch = 0.0057_wp ! heat transfer coefficient 282 REAL(wp), DIMENSION(jpi,jpj) :: zfric, zvel ! ice-ocean velocity (m/s) and frictional velocity (m2/s2) 283 ! 284 ! for frazil ice 285 INTEGER :: iter 286 REAL(wp) :: zvfrx, zvgx, ztaux, zf, ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, ztwogp 287 REAL(wp), PARAMETER :: zcai = 1.4e-3_wp ! ice-air drag (clem: should be dependent on coupling/forcing used) 288 REAL(wp), PARAMETER :: zhicrit = 0.04_wp ! frazil ice thickness 289 REAL(wp), PARAMETER :: zsqcd = 1.0_wp / SQRT( 1.3_wp * zcai ) ! 1/SQRT(airdensity*drag) 290 REAL(wp), PARAMETER :: zgamafr = 0.03_wp 291 !!----------------------------------------------------------------------- 292 ! 126 293 ! computation of friction velocity at T points 127 !---------------------------------------------!128 294 IF( ln_icedyn ) THEN 129 zu_io(:,:) = u_ice(:,:) - ssu_m(:,:)130 zv_io(:,:) = v_ice(:,:) - ssv_m(:,:)131 295 DO_2D( 0, 0, 0, 0 ) 132 zfric(ji,jj) = rn_cio * ( 0.5_wp * & 133 & ( zu_io(ji,jj) * zu_io(ji,jj) + zu_io(ji-1,jj) * zu_io(ji-1,jj) & 134 & + zv_io(ji,jj) * zv_io(ji,jj) + zv_io(ji,jj-1) * zv_io(ji,jj-1) ) ) * tmask(ji,jj,1) 135 zvel(ji,jj) = 0.5_wp * SQRT( ( u_ice(ji-1,jj) + u_ice(ji,jj) ) * ( u_ice(ji-1,jj) + u_ice(ji,jj) ) + & 136 & ( v_ice(ji,jj-1) + v_ice(ji,jj) ) * ( v_ice(ji,jj-1) + v_ice(ji,jj) ) ) 296 zu_io = u_ice(ji ,jj ) - ssu_m(ji ,jj ) 297 zu_iom1 = u_ice(ji-1,jj ) - ssu_m(ji-1,jj ) 298 zv_io = v_ice(ji ,jj ) - ssv_m(ji ,jj ) 299 zv_iom1 = v_ice(ji ,jj-1) - ssv_m(ji ,jj-1) 300 ! 301 zfric(ji,jj) = rn_cio * ( 0.5_wp * ( zu_io*zu_io + zu_iom1*zu_iom1 + zv_io*zv_io + zv_iom1*zv_iom1 ) ) * tmask(ji,jj,1) 302 zvel (ji,jj) = 0.5_wp * SQRT( ( u_ice(ji-1,jj ) + u_ice(ji,jj) ) * ( u_ice(ji-1,jj ) + u_ice(ji,jj) ) + & 303 & ( v_ice(ji ,jj-1) + v_ice(ji,jj) ) * ( v_ice(ji ,jj-1) + v_ice(ji,jj) ) ) 137 304 END_2D 138 305 ELSE ! if no ice dynamics => transfer directly the atmospheric stress to the ocean … … 227 394 ENDIF 228 395 229 230 396 !---------------------------------------------------------! 231 397 ! Collection thickness of ice formed in leads and polynyas … … 251 417 IF ( qlead(ji,jj) < 0._wp ) THEN ! cooling 252 418 ! -- Wind stress -- ! 253 ztaux = ( utau_ice(ji-1,jj ) * umask(ji-1,jj ,1) & 254 & + utau_ice(ji ,jj ) * umask(ji ,jj ,1) ) * 0.5_wp 255 ztauy = ( vtau_ice(ji ,jj-1) * vmask(ji ,jj-1,1) & 256 & + vtau_ice(ji ,jj ) * vmask(ji ,jj ,1) ) * 0.5_wp 419 ztaux = ( utau_ice(ji-1,jj ) * umask(ji-1,jj ,1) + utau_ice(ji,jj) * umask(ji,jj,1) ) * 0.5_wp 420 ztauy = ( vtau_ice(ji ,jj-1) * vmask(ji ,jj-1,1) + vtau_ice(ji,jj) * vmask(ji,jj,1) ) * 0.5_wp 257 421 ! Square root of wind stress 258 ztenagm =SQRT( SQRT( ztaux * ztaux + ztauy * ztauy ) )422 ztenagm = SQRT( SQRT( ztaux * ztaux + ztauy * ztauy ) ) 259 423 260 424 ! -- Frazil ice velocity -- ! … … 268 432 269 433 ! -- Relative frazil/pack ice velocity -- ! 270 rswitch = MAX( 0._wp, SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) 271 zvrel2 = MAX( ( zvfrx - zvgx ) * ( zvfrx - zvgx ) & 272 & + ( zvfry - zvgy ) * ( zvfry - zvgy ) , 0.15_wp * 0.15_wp ) * rswitch 273 274 !!clem 434 rswitch = MAX( 0._wp, SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) 435 zvrel2 = MAX( (zvfrx - zvgx)*(zvfrx - zvgx) + (zvfry - zvgy)*(zvfry - zvgy), 0.15_wp*0.15_wp ) * rswitch 436 437 ! -- fraction of frazil ice -- ! 275 438 fraz_frac(ji,jj) = rswitch * ( TANH( rn_Cfraz * ( SQRT(zvrel2) - rn_vfraz ) ) + 1._wp ) * 0.5_wp * rn_maxfraz 276 !!clem277 439 278 440 ! -- new ice thickness (iterative loop) -- ! … … 301 463 302 464 ENDIF 303 304 !-------------------------------------------------------------------------------------------! 305 ! Thermodynamic computation (only on grid points covered by ice) => loop over ice categories 306 !-------------------------------------------------------------------------------------------! 307 DO jl = 1, jpl 308 309 ! select ice covered grid points 310 npti = 0 ; nptidx(:) = 0 311 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 312 IF ( a_i(ji,jj,jl) > epsi10 ) THEN 313 npti = npti + 1 314 nptidx(npti) = (jj - 1) * jpi + ji 315 ENDIF 316 END_2D 317 318 IF( npti > 0 ) THEN ! If there is no ice, do nothing. 319 ! 320 CALL ice_thd_1d2d( jl, 1 ) ! --- Move to 1D arrays --- ! 321 ! ! --- & Change units of e_i, e_s from J/m2 to J/m3 --- ! 322 ! 323 s_i_new (1:npti) = 0._wp ; dh_s_tot(1:npti) = 0._wp ! --- some init --- ! (important to have them here) 324 dh_i_sum (1:npti) = 0._wp ; dh_i_bom(1:npti) = 0._wp ; dh_i_itm (1:npti) = 0._wp 325 dh_i_sub (1:npti) = 0._wp ; dh_i_bog(1:npti) = 0._wp 326 dh_snowice(1:npti) = 0._wp ; dh_s_mlt(1:npti) = 0._wp 327 ! 328 CALL ice_thd_zdf ! --- Ice-Snow temperature --- ! 329 ! 330 IF( ln_icedH ) THEN ! --- Growing/Melting --- ! 331 CALL ice_thd_dh ! Ice-Snow thickness 332 CALL ice_thd_ent( e_i_1d(1:npti,:) ) ! Ice enthalpy remapping 333 ENDIF 334 CALL ice_thd_sal( ln_icedS ) ! --- Ice salinity --- ! 335 ! 336 CALL ice_thd_temp ! --- Temperature update --- ! 337 ! 338 IF( ln_icedH .AND. ln_virtual_itd ) & 339 & CALL ice_thd_mono ! --- Extra lateral melting if virtual_itd --- ! 340 ! 341 IF( ln_icedA ) CALL ice_thd_da ! --- Lateral melting --- ! 342 ! 343 CALL ice_thd_1d2d( jl, 2 ) ! --- Change units of e_i, e_s from J/m3 to J/m2 --- ! 344 ! ! --- & Move to 2D arrays --- ! 345 ENDIF 346 ! 347 END DO 348 ! 349 IF( ln_icediachk ) CALL ice_cons_hsm(1, 'icethd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 350 IF( ln_icediachk ) CALL ice_cons2D (1, 'icethd', diag_v, diag_s, diag_t, diag_fv, diag_fs, diag_ft) 351 ! 352 IF ( ln_pnd .AND. ln_icedH ) & 353 & CALL ice_thd_pnd ! --- Melt ponds --- ! 354 ! 355 IF( jpl > 1 ) CALL ice_itd_rem( kt ) ! --- Transport ice between thickness categories --- ! 356 ! 357 IF( ln_icedO ) CALL ice_thd_do ! --- Frazil ice growth in leads --- ! 358 ! 359 CALL ice_cor( kt , 2 ) ! --- Corrections --- ! 360 ! 361 oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * rDt_ice ! --- Ice natural aging incrementation 362 ! 363 DO_2D( 0, 0, 0, 0 ) ! --- Ice velocity corrections 364 IF( at_i(ji,jj) == 0._wp ) THEN ! if ice has melted 365 IF( at_i(ji+1,jj) == 0._wp ) u_ice(ji ,jj) = 0._wp ! right side 366 IF( at_i(ji-1,jj) == 0._wp ) u_ice(ji-1,jj) = 0._wp ! left side 367 IF( at_i(ji,jj+1) == 0._wp ) v_ice(ji,jj ) = 0._wp ! upper side 368 IF( at_i(ji,jj-1) == 0._wp ) v_ice(ji,jj-1) = 0._wp ! bottom side 369 ENDIF 370 END_2D 371 CALL lbc_lnk( 'icecor', u_ice, 'U', -1.0_wp, v_ice, 'V', -1.0_wp ) 372 ! 373 ! convergence tests 374 IF( ln_zdf_chkcvg ) THEN 375 CALL iom_put( 'tice_cvgerr', ztice_cvgerr ) ; DEALLOCATE( ztice_cvgerr ) 376 CALL iom_put( 'tice_cvgstp', ztice_cvgstp ) ; DEALLOCATE( ztice_cvgstp ) 377 ENDIF 378 ! 379 ! controls 380 IF( ln_icectl ) CALL ice_prt (kt, iiceprt, jiceprt, 1, ' - ice thermodyn. - ') ! prints 381 IF( sn_cfctl%l_prtctl ) & 382 & CALL ice_prt3D ('icethd') ! prints 383 IF( ln_timing ) CALL timing_stop('icethd') ! timing 384 ! 385 END SUBROUTINE ice_thd 386 387 388 SUBROUTINE ice_thd_temp 389 !!----------------------------------------------------------------------- 390 !! *** ROUTINE ice_thd_temp *** 391 !! 392 !! ** Purpose : Computes sea ice temperature (Kelvin) from enthalpy 393 !! 394 !! ** Method : Formula (Bitz and Lipscomb, 1999) 395 !!------------------------------------------------------------------- 396 INTEGER :: ji, jk ! dummy loop indices 397 REAL(wp) :: ztmelts, zbbb, zccc ! local scalar 398 !!------------------------------------------------------------------- 399 ! Recover ice temperature 400 DO jk = 1, nlay_i 401 DO ji = 1, npti 402 ztmelts = -rTmlt * sz_i_1d(ji,jk) 403 ! Conversion q(S,T) -> T (second order equation) 404 zbbb = ( rcp - rcpi ) * ztmelts + e_i_1d(ji,jk) * r1_rhoi - rLfus 405 zccc = SQRT( MAX( zbbb * zbbb - 4._wp * rcpi * rLfus * ztmelts, 0._wp ) ) 406 t_i_1d(ji,jk) = rt0 - ( zbbb + zccc ) * 0.5_wp * r1_rcpi 407 408 ! mask temperature 409 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - h_i_1d(ji) ) ) 410 t_i_1d(ji,jk) = rswitch * t_i_1d(ji,jk) + ( 1._wp - rswitch ) * rt0 411 END DO 412 END DO 413 ! 414 END SUBROUTINE ice_thd_temp 415 416 417 SUBROUTINE ice_thd_mono 418 !!----------------------------------------------------------------------- 419 !! *** ROUTINE ice_thd_mono *** 420 !! 421 !! ** Purpose : Lateral melting in case virtual_itd 422 !! ( dA = A/2h dh ) 423 !!----------------------------------------------------------------------- 424 INTEGER :: ji ! dummy loop indices 425 REAL(wp) :: zhi_bef ! ice thickness before thermo 426 REAL(wp) :: zdh_mel, zda_mel ! net melting 427 REAL(wp) :: zvi, zvs ! ice/snow volumes 428 !!----------------------------------------------------------------------- 429 ! 430 DO ji = 1, npti 431 zdh_mel = MIN( 0._wp, dh_i_itm(ji) + dh_i_sum(ji) + dh_i_bom(ji) + dh_snowice(ji) + dh_i_sub(ji) ) 432 IF( zdh_mel < 0._wp .AND. a_i_1d(ji) > 0._wp ) THEN 433 zvi = a_i_1d(ji) * h_i_1d(ji) 434 zvs = a_i_1d(ji) * h_s_1d(ji) 435 ! lateral melting = concentration change 436 zhi_bef = h_i_1d(ji) - zdh_mel 437 rswitch = MAX( 0._wp , SIGN( 1._wp , zhi_bef - epsi20 ) ) 438 zda_mel = rswitch * a_i_1d(ji) * zdh_mel / ( 2._wp * MAX( zhi_bef, epsi20 ) ) 439 a_i_1d(ji) = MAX( epsi20, a_i_1d(ji) + zda_mel ) 440 ! adjust thickness 441 h_i_1d(ji) = zvi / a_i_1d(ji) 442 h_s_1d(ji) = zvs / a_i_1d(ji) 443 ! retrieve total concentration 444 at_i_1d(ji) = a_i_1d(ji) 445 END IF 446 END DO 447 ! 448 END SUBROUTINE ice_thd_mono 449 465 466 END SUBROUTINE thd_prep 450 467 451 468 SUBROUTINE ice_thd_1d2d( kl, kn )
Note: See TracChangeset
for help on using the changeset viewer.
