Changeset 8559
- Timestamp:
- 2017-09-22T16:55:24+02:00 (7 years ago)
- Location:
- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3
- Files:
-
- 5 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice1D.F90
r8554 r8559 100 100 ! ! to reintegrate longwave flux inside the ice thermodynamics 101 101 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: i0 !: fraction of radiation transmitted to the ice 102 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: hicol_1d !: Ice collection thickness accumulated in leads103 102 104 103 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: t_su_1d !: <==> the 2D t_su … … 197 196 & sfx_bri_1d (jpij) , sfx_bog_1d (jpij) , sfx_bom_1d (jpij) , sfx_sum_1d (jpij), & 198 197 & sfx_sni_1d (jpij) , sfx_opw_1d (jpij) , sfx_res_1d (jpij) , sfx_sub_1d (jpij), & 199 & sfx_lam_1d (jpij) , sfx_dyn_1d(jpij) , hicol_1d (jpij), STAT=ierr(ii) )198 & sfx_lam_1d (jpij) , sfx_dyn_1d(jpij) , STAT=ierr(ii) ) 200 199 ! 201 200 ii = ii + 1 -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_dh.F90
r8534 r8559 293 293 DO jk = 1, nlay_i 294 294 DO ji = 1, nidx 295 ztmelts = - tmut * s_i_1d(ji,jk) + rt0 ! Melting point of layer k [K]296 297 IF( t_i_1d(ji,jk) >= ztmelts) THEN !!! Internal melting295 ztmelts = - tmut * s_i_1d(ji,jk) ! Melting point of layer k [C] 296 297 IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN !!! Internal melting 298 298 299 299 zEi = - e_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of layer k [J/kg, <0] … … 320 320 321 321 zEi = - e_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of layer k [J/kg, <0] 322 zEw = rcp * ( ztmelts - rt0 )! Specific enthalpy of resulting meltwater [J/kg, <0]322 zEw = rcp * ztmelts ! Specific enthalpy of resulting meltwater [J/kg, <0] 323 323 zdE = zEi - zEw ! Specific enthalpy difference < 0 324 324 … … 430 430 + ( 1. - zswitch_sal ) * sm_i_1d(ji) 431 431 ! New ice growth 432 ztmelts = - tmut * s_i_new(ji) + rt0 ! New ice melting point (K)432 ztmelts = - tmut * s_i_new(ji) ! New ice melting point (C) 433 433 434 434 zt_i_new = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i) 435 435 436 zEi = cpic * ( zt_i_new - ztmelts) & ! Specific enthalpy of forming ice (J/kg, <0)437 & - lfus * ( 1.0 - ( ztmelts - rt0 )/ ( zt_i_new - rt0 ) ) &438 & + rcp * ( ztmelts-rt0 )436 zEi = cpic * ( zt_i_new - (ztmelts+rt0) ) & ! Specific enthalpy of forming ice (J/kg, <0) 437 & - lfus * ( 1.0 - ztmelts / ( zt_i_new - rt0 ) ) & 438 & + rcp * ztmelts 439 439 440 440 zEw = rcp * ( t_bo_1d(ji) - rt0 ) ! Specific enthalpy of seawater (J/kg, < 0) … … 450 450 zfmdt = - rhoic * dh_i_bott(ji) ! Mass flux x time step (kg/m2, < 0) 451 451 452 ztmelts = - tmut * s_i_new(ji) + rt0 ! New ice melting point (K)452 ztmelts = - tmut * s_i_new(ji) ! New ice melting point (C) 453 453 454 454 zt_i_new = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i) 455 455 456 zEi = cpic * ( zt_i_new - ztmelts) & ! Specific enthalpy of forming ice (J/kg, <0)457 & - lfus * ( 1.0 - ( ztmelts - rt0 )/ ( zt_i_new - rt0 ) ) &458 & + rcp * ( ztmelts-rt0 )456 zEi = cpic * ( zt_i_new - (ztmelts+rt0) ) & ! Specific enthalpy of forming ice (J/kg, <0) 457 & - lfus * ( 1.0 - ztmelts / ( zt_i_new - rt0 ) ) & 458 & + rcp * ztmelts 459 459 460 460 zEw = rcp * ( t_bo_1d(ji) - rt0 ) ! Specific enthalpy of seawater (J/kg, < 0) … … 490 490 IF( zf_tt(ji) > 0._wp .AND. jk > icount(ji,jk) ) THEN ! do not calculate where layer has already disappeared by surface melting 491 491 492 ztmelts = - tmut * s_i_1d(ji,jk) + rt0 ! Melting point of layer jk (K)493 494 IF( t_i_1d(ji,jk) >= ztmelts) THEN !!! Internal melting492 ztmelts = - tmut * s_i_1d(ji,jk) ! Melting point of layer jk (C) 493 494 IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN !!! Internal melting 495 495 496 496 zEi = - e_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0) … … 520 520 521 521 zEi = - e_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0) 522 zEw = rcp * ( ztmelts - rt0 )! Specific enthalpy of meltwater (J/kg, <0)522 zEw = rcp * ztmelts ! Specific enthalpy of meltwater (J/kg, <0) 523 523 zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0) 524 524 -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_do.F90
r8534 r8559 80 80 INTEGER :: ji,jj,jk,jl ! dummy loop indices 81 81 INTEGER :: iter ! - - 82 REAL(wp) :: ztmelts, z dv, zfrazb, zweight, zde ! local scalars82 REAL(wp) :: ztmelts, zfrazb, zweight, zde ! local scalars 83 83 REAL(wp) :: zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf ! - - 84 84 REAL(wp) :: ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, zsqcd , zhicrit ! - - … … 102 102 REAL(wp), DIMENSION(jpij) :: zdv_res ! residual volume in case of excessive heat budget 103 103 REAL(wp), DIMENSION(jpij) :: zda_res ! residual area in case of excessive heat budget 104 REAL(wp), DIMENSION(jpij) :: zat_i_1d ! total ice fraction105 104 REAL(wp), DIMENSION(jpij) :: zv_frazb ! accretion of frazil ice at the ice bottom 106 105 REAL(wp), DIMENSION(jpij) :: zvrel_1d ! relative ice / frazil velocity (1D vector) … … 108 107 REAL(wp), DIMENSION(jpij,jpl) :: zv_b ! old volume of ice in category jl 109 108 REAL(wp), DIMENSION(jpij,jpl) :: za_b ! old area of ice in category jl 110 REAL(wp), DIMENSION(jpij,jpl) :: za_i_1d ! 1-D version of a_i 111 REAL(wp), DIMENSION(jpij,jpl) :: zv_i_1d ! 1-D version of v_i 112 REAL(wp), DIMENSION(jpij,jpl) :: zsmv_i_1d ! 1-D version of smv_i 113 114 REAL(wp), DIMENSION(jpij,nlay_i,jpl) :: ze_i_1d !: 1-D version of e_i 109 110 REAL(wp), DIMENSION(jpij,nlay_i,jpl) :: ze_i_2d !: 1-D version of e_i 115 111 116 112 REAL(wp), DIMENSION(jpi,jpj) :: zvrel ! relative ice / frazil velocity … … 243 239 IF ( nidx > 0 ) THEN 244 240 245 CALL tab_2d_1d( nidx, idxice(1:nidx), zat_i_1d (1:nidx) , at_i ) 241 CALL tab_2d_1d( nidx, idxice(1:nidx), at_i_1d (1:nidx) , at_i ) 242 CALL tab_3d_2d( nidx, idxice(1:nidx), a_i_2d (1:nidx,1:jpl), a_i (:,:,:) ) 243 CALL tab_3d_2d( nidx, idxice(1:nidx), v_i_2d (1:nidx,1:jpl), v_i (:,:,:) ) 244 CALL tab_3d_2d( nidx, idxice(1:nidx), smv_i_2d(1:nidx,1:jpl), smv_i(:,:,:) ) 246 245 DO jl = 1, jpl 247 CALL tab_2d_1d( nidx, idxice(1:nidx), za_i_1d (1:nidx,jl), a_i (:,:,jl) )248 CALL tab_2d_1d( nidx, idxice(1:nidx), zv_i_1d (1:nidx,jl), v_i (:,:,jl) )249 CALL tab_2d_1d( nidx, idxice(1:nidx), zsmv_i_1d(1:nidx,jl), smv_i(:,:,jl) )250 246 DO jk = 1, nlay_i 251 CALL tab_2d_1d( nidx, idxice(1:nidx), ze_i_1d(1:nidx,jk,jl), e_i(:,:,jk,jl) ) 252 END DO 253 END DO 254 255 CALL tab_2d_1d( nidx, idxice(1:nidx), qlead_1d (1:nidx) , qlead ) 256 CALL tab_2d_1d( nidx, idxice(1:nidx), t_bo_1d (1:nidx) , t_bo ) 257 CALL tab_2d_1d( nidx, idxice(1:nidx), sfx_opw_1d(1:nidx) , sfx_opw ) 258 CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_opw_1d(1:nidx) , wfx_opw ) 259 CALL tab_2d_1d( nidx, idxice(1:nidx), hicol_1d (1:nidx) , hicol ) 260 CALL tab_2d_1d( nidx, idxice(1:nidx), zvrel_1d (1:nidx) , zvrel ) 261 262 CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_thd_1d(1:nidx) , hfx_thd ) 263 CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_opw_1d(1:nidx) , hfx_opw ) 264 CALL tab_2d_1d( nidx, idxice(1:nidx), rn_amax_1d(1:nidx) , rn_amax_2d ) 265 CALL tab_2d_1d( nidx, idxice(1:nidx), sss_1d (1:nidx) , sss_m ) 247 CALL tab_2d_1d( nidx, idxice(1:nidx), ze_i_2d(1:nidx,jk,jl), e_i(:,:,jk,jl) ) 248 END DO 249 END DO 250 CALL tab_2d_1d( nidx, idxice(1:nidx), qlead_1d (1:nidx) , qlead ) 251 CALL tab_2d_1d( nidx, idxice(1:nidx), t_bo_1d (1:nidx) , t_bo ) 252 CALL tab_2d_1d( nidx, idxice(1:nidx), sfx_opw_1d(1:nidx) , sfx_opw ) 253 CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_opw_1d(1:nidx) , wfx_opw ) 254 CALL tab_2d_1d( nidx, idxice(1:nidx), zh_newice (1:nidx) , hicol ) 255 CALL tab_2d_1d( nidx, idxice(1:nidx), zvrel_1d (1:nidx) , zvrel ) 256 257 CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_thd_1d(1:nidx) , hfx_thd ) 258 CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_opw_1d(1:nidx) , hfx_opw ) 259 CALL tab_2d_1d( nidx, idxice(1:nidx), rn_amax_1d(1:nidx) , rn_amax_2d ) 260 CALL tab_2d_1d( nidx, idxice(1:nidx), sss_1d (1:nidx) , sss_m ) 266 261 267 262 !------------------------------------------------------------------------------| … … 269 264 !------------------------------------------------------------------------------| 270 265 DO jl = 1, jpl 271 DO jk = 1, nlay_i 272 DO ji = 1, nidx 273 IF( zv_i_1d(ji,jl) > 0._wp ) ze_i_1d(ji,jk,jl) = ze_i_1d(ji,jk,jl) / zv_i_1d(ji,jl) * REAL( nlay_i ) 274 END DO 266 DO jk = 1, nlay_i 267 WHERE( v_i_2d(1:nidx,jl) > 0._wp ) 268 ze_i_2d(1:nidx,jk,jl) = ze_i_2d(1:nidx,jk,jl) / v_i_2d(1:nidx,jl) * REAL( nlay_i ) 269 ELSEWHERE 270 ze_i_2d(1:nidx,jk,jl) = 0._wp 271 END WHERE 275 272 END DO 276 273 END DO … … 282 279 ! Keep old ice areas and volume in memory 283 280 !----------------------------------------- 284 zv_b(1:nidx,:) = zv_i_1d(1:nidx,:) 285 za_b(1:nidx,:) = za_i_1d(1:nidx,:) 286 287 !---------------------- 288 ! Thickness of new ice 289 !---------------------- 290 zh_newice(1:nidx) = hicol_1d(1:nidx) 281 zv_b(1:nidx,:) = v_i_2d(1:nidx,:) 282 za_b(1:nidx,:) = a_i_2d(1:nidx,:) 291 283 292 284 !---------------------- … … 309 301 ! We assume that new ice is formed at the seawater freezing point 310 302 DO ji = 1, nidx 311 ztmelts = - tmut * zs_newice(ji) + rt0 ! Melting point (K)312 ze_newice(ji) = rhoic * ( cpic * ( ztmelts - t_bo_1d(ji) )&313 & + lfus * ( 1.0 - ( ztmelts - rt0 )/ MIN( t_bo_1d(ji) - rt0, -epsi10 ) ) &314 & - rcp * ( ztmelts - rt0 ))303 ztmelts = - tmut * zs_newice(ji) ! Melting point (C) 304 ze_newice(ji) = rhoic * ( cpic * ( ztmelts - ( t_bo_1d(ji) - rt0 ) ) & 305 & + lfus * ( 1.0 - ztmelts / MIN( t_bo_1d(ji) - rt0, -epsi10 ) ) & 306 & - rcp * ztmelts ) 315 307 END DO 316 308 … … 318 310 ! Age of new ice 319 311 !---------------- 320 DO ji = 1, nidx 321 zo_newice(ji) = 0._wp 322 END DO 312 zo_newice(1:nidx) = 0._wp 323 313 324 314 !------------------- … … 350 340 END DO 351 341 352 zv_frazb( :) = 0._wp342 zv_frazb(1:nidx) = 0._wp 353 343 IF( ln_frazil ) THEN 354 344 ! A fraction zfrazb of frazil ice is accreted at the ice bottom 355 345 DO ji = 1, nidx 356 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp , - zat_i_1d(ji) ) )346 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp , - at_i_1d(ji) ) ) 357 347 zfrazb = rswitch * ( TANH( rn_Cfraz * ( zvrel_1d(ji) - rn_vfraz ) ) + 1.0 ) * 0.5 * rn_maxfraz 358 348 zv_frazb(ji) = zfrazb * zv_newice(ji) … … 378 368 ! we keep the excessive volume in memory and attribute it later to bottom accretion 379 369 DO ji = 1, nidx 380 IF ( za_newice(ji) > ( rn_amax_1d(ji) - zat_i_1d(ji) ) ) THEN381 zda_res(ji) = za_newice(ji) - ( rn_amax_1d(ji) - zat_i_1d(ji) )370 IF ( za_newice(ji) > ( rn_amax_1d(ji) - at_i_1d(ji) ) ) THEN 371 zda_res(ji) = za_newice(ji) - ( rn_amax_1d(ji) - at_i_1d(ji) ) 382 372 zdv_res(ji) = zda_res (ji) * zh_newice(ji) 383 373 za_newice(ji) = za_newice(ji) - zda_res (ji) … … 390 380 391 381 ! find which category to fill 392 zat_i_1d(:) = 0._wp393 382 DO jl = 1, jpl 394 383 DO ji = 1, nidx 395 384 IF( zh_newice(ji) > hi_max(jl-1) .AND. zh_newice(ji) <= hi_max(jl) ) THEN 396 za_i_1d (ji,jl) = za_i_1d(ji,jl) + za_newice(ji)397 zv_i_1d (ji,jl) = zv_i_1d(ji,jl) + zv_newice(ji)398 jcat (ji)= jl385 a_i_2d(ji,jl) = a_i_2d(ji,jl) + za_newice(ji) 386 v_i_2d(ji,jl) = v_i_2d(ji,jl) + zv_newice(ji) 387 jcat(ji) = jl 399 388 ENDIF 400 zat_i_1d(ji) = zat_i_1d(ji) + za_i_1d (ji,jl)401 402 END DO389 END DO 390 END DO 391 at_i_1d(1:nidx) = SUM( a_i_2d(1:nidx,:), dim=2 ) 403 392 404 393 ! Heat content … … 411 400 DO ji = 1, nidx 412 401 jl = jcat(ji) 413 rswitch = MAX( 0._wp, SIGN( 1._wp , zv_i_1d(ji,jl) - epsi20 ) )414 ze_i_ 1d(ji,jk,jl) = zswinew(ji) * ze_newice(ji) + &415 & ( 1.0 - zswinew(ji) ) * ( ze_newice(ji) * zv_newice(ji) + ze_i_ 1d(ji,jk,jl) * zv_b(ji,jl) ) &416 & * rswitch / MAX( zv_i_1d(ji,jl), epsi20 )402 rswitch = MAX( 0._wp, SIGN( 1._wp , v_i_2d(ji,jl) - epsi20 ) ) 403 ze_i_2d(ji,jk,jl) = zswinew(ji) * ze_newice(ji) + & 404 & ( 1.0 - zswinew(ji) ) * ( ze_newice(ji) * zv_newice(ji) + ze_i_2d(ji,jk,jl) * zv_b(ji,jl) ) & 405 & * rswitch / MAX( v_i_2d(ji,jl), epsi20 ) 417 406 END DO 418 407 END DO … … 428 417 DO jk = 1, nlay_i 429 418 DO ji = 1, nidx 430 h_i_old (ji,jk) = zv_i_1d(ji,jl) * r1_nlay_i431 eh_i_old(ji,jk) = ze_i_ 1d(ji,jk,jl) * h_i_old(ji,jk)419 h_i_old (ji,jk) = v_i_2d(ji,jl) * r1_nlay_i 420 eh_i_old(ji,jk) = ze_i_2d(ji,jk,jl) * h_i_old(ji,jk) 432 421 END DO 433 422 END DO … … 435 424 ! new volumes including lateral/bottom accretion + residual 436 425 DO ji = 1, nidx 437 rswitch = MAX( 0._wp, SIGN( 1._wp , zat_i_1d(ji) - epsi20 ) )438 zv_newfra = rswitch * ( zdv_res(ji) + zv_frazb(ji) ) * za_i_1d(ji,jl) / MAX( zat_i_1d(ji) , epsi20 )439 za_i_1d(ji,jl) = rswitch * za_i_1d(ji,jl)440 zv_i_1d(ji,jl) = zv_i_1d(ji,jl) + zv_newfra426 rswitch = MAX( 0._wp, SIGN( 1._wp , at_i_1d(ji) - epsi20 ) ) 427 zv_newfra = rswitch * ( zdv_res(ji) + zv_frazb(ji) ) * a_i_2d(ji,jl) / MAX( at_i_1d(ji) , epsi20 ) 428 a_i_2d(ji,jl) = rswitch * a_i_2d(ji,jl) 429 v_i_2d(ji,jl) = v_i_2d(ji,jl) + zv_newfra 441 430 ! for remapping 442 431 h_i_old (ji,nlay_i+1) = zv_newfra … … 444 433 ENDDO 445 434 ! --- Ice enthalpy remapping --- ! 446 CALL ice_thd_ent( ze_i_ 1d(1:nidx,:,jl) )435 CALL ice_thd_ent( ze_i_2d(1:nidx,:,jl) ) 447 436 ENDDO 448 437 … … 452 441 DO jl = 1, jpl 453 442 DO ji = 1, nidx 454 zdv = zv_i_1d(ji,jl) - zv_b(ji,jl) 455 zsmv_i_1d(ji,jl) = zsmv_i_1d(ji,jl) + zdv * zs_newice(ji) 443 smv_i_2d(ji,jl) = smv_i_2d(ji,jl) + zs_newice(ji) * ( v_i_2d(ji,jl) - zv_b(ji,jl) ) 456 444 END DO 457 445 END DO … … 462 450 DO jl = 1, jpl 463 451 DO jk = 1, nlay_i 464 DO ji = 1, nidx 465 ze_i_1d(ji,jk,jl) = ze_i_1d(ji,jk,jl) * zv_i_1d(ji,jl) * r1_nlay_i 466 END DO 452 ze_i_2d(1:nidx,jk,jl) = ze_i_2d(1:nidx,jk,jl) * v_i_2d(1:nidx,jl) * r1_nlay_i 467 453 END DO 468 454 END DO … … 470 456 ! 7) Change 2D vectors to 1D vectors 471 457 !------------------------------------------------------------------------------! 472 DO jl = 1, jpl473 CALL tab_1d_2d( nidx, idxice(1:nidx), za_i_1d (1:nidx,jl), a_i (:,:,jl))474 CALL tab_1d_2d( nidx, idxice(1:nidx), zv_i_1d (1:nidx,jl), v_i (:,:,jl))475 CALL tab_1d_2d( nidx, idxice(1:nidx), zsmv_i_1d(1:nidx,jl), smv_i (:,:,jl) )458 CALL tab_2d_3d( nidx, idxice(1:nidx), a_i_2d (1:nidx,1:jpl), a_i (:,:,:) ) 459 CALL tab_2d_3d( nidx, idxice(1:nidx), v_i_2d (1:nidx,1:jpl), v_i (:,:,:) ) 460 CALL tab_2d_3d( nidx, idxice(1:nidx), smv_i_2d(1:nidx,1:jpl), smv_i(:,:,:) ) 461 DO jl = 1, jpl 476 462 DO jk = 1, nlay_i 477 CALL tab_1d_2d( nidx, idxice(1:nidx), ze_i_ 1d(1:nidx,jk,jl), e_i(:,:,jk,jl))478 END DO 479 END DO 480 CALL tab_1d_2d( nidx, idxice(1:nidx), sfx_opw_1d(1:nidx), sfx_opw 481 CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_opw_1d(1:nidx), wfx_opw 482 CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_thd_1d(1:nidx), hfx_thd 483 CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_opw_1d(1:nidx), hfx_opw 463 CALL tab_1d_2d( nidx, idxice(1:nidx), ze_i_2d(1:nidx,jk,jl), e_i(:,:,jk,jl) ) 464 END DO 465 END DO 466 CALL tab_1d_2d( nidx, idxice(1:nidx), sfx_opw_1d(1:nidx), sfx_opw ) 467 CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_opw_1d(1:nidx), wfx_opw ) 468 CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_thd_1d(1:nidx), hfx_thd ) 469 CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_opw_1d(1:nidx), hfx_opw ) 484 470 ! 485 471 ENDIF ! nidx > 0 -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_sal.F90
r8534 r8559 59 59 REAL(wp) :: iflush, igravdr ! local scalars 60 60 REAL(wp) :: zs_sni, zsm_i_gd, zsm_i_fl, zsm_i_si, zsm_i_bg ! local scalars 61 REAL(wp) :: z1_time_gd, z1_time_fl 61 62 !!--------------------------------------------------------------------- 62 63 … … 66 67 CASE( 2 ) ! time varying salinity with linear profile ! 67 68 ! !---------------------------------------------! 69 z1_time_gd = 1._wp / rn_time_gd * rdt_ice 70 z1_time_fl = 1._wp / rn_time_fl * rdt_ice 71 ! 68 72 DO ji = 1, nidx 69 73 … … 71 75 ! Update ice salinity from snow-ice and bottom growth 72 76 !--------------------------------------------------------- 73 zs_sni = sss_1d(ji) * ( rhoic - rhosn ) * r1_rhoic ! Salinity of snow ice74 rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi20 ) )75 zsm_i_si = ( zs_sni - sm_i_1d(ji) ) * dh_snowice(ji) / MAX( ht_i_1d(ji), epsi20 ) * rswitch! snow-ice76 zsm_i_bg = ( s_i_new(ji) - sm_i_1d(ji) ) * MAX( 0._wp, dh_i_bott(ji) ) / MAX( ht_i_1d(ji), epsi20 ) * rswitch! bottom growth77 78 ! Update salinity (nb: salt flux already included in icethd_dh)79 sm_i_1d(ji) = sm_i_1d(ji) + zsm_i_bg + zsm_i_si77 IF( ht_i_1d(ji) > 0._wp ) THEN 78 zs_sni = sss_1d(ji) * ( rhoic - rhosn ) * r1_rhoic ! Salinity of snow ice 79 zsm_i_si = ( zs_sni - sm_i_1d(ji) ) * dh_snowice(ji) / ht_i_1d(ji) ! snow-ice 80 zsm_i_bg = ( s_i_new(ji) - sm_i_1d(ji) ) * MAX( 0._wp, dh_i_bott(ji) ) / ht_i_1d(ji) ! bottom growth 81 ! Update salinity (nb: salt flux already included in icethd_dh) 82 sm_i_1d(ji) = sm_i_1d(ji) + zsm_i_bg + zsm_i_si 83 ENDIF 80 84 81 85 IF( ld_sal ) THEN … … 85 89 iflush = MAX( 0._wp , SIGN( 1._wp , t_su_1d(ji) - rt0 ) ) ! =1 if summer 86 90 igravdr = MAX( 0._wp , SIGN( 1._wp , t_bo_1d(ji) - t_su_1d(ji) ) ) ! =1 if t_su < t_bo 87 zsm_i_gd = - igravdr * MAX( sm_i_1d(ji) - rn_sal_gd , 0._wp ) / rn_time_gd * rdt_ice ! gravity drainage 88 zsm_i_fl = - iflush * MAX( sm_i_1d(ji) - rn_sal_fl , 0._wp ) / rn_time_fl * rdt_ice ! flushing 91 92 zsm_i_gd = - igravdr * MAX( sm_i_1d(ji) - rn_sal_gd , 0._wp ) * z1_time_gd ! gravity drainage 93 zsm_i_fl = - iflush * MAX( sm_i_1d(ji) - rn_sal_fl , 0._wp ) * z1_time_fl ! flushing 89 94 90 95 ! Update salinity -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90
r8534 r8559 42 42 !! ice_var_eqv2glo : transform from VEQV to VGLO 43 43 !! ice_var_salprof : salinity profile in the ice 44 !! ice_var_bv : brine volume45 44 !! ice_var_salprof1d : salinity profile in the ice 1D 46 45 !! ice_var_zapsmall : remove very small area and volume 47 46 !! ice_var_itd : convert 1-cat to multiple cat 47 !! ice_var_bv : brine volume 48 48 !!---------------------------------------------------------------------- 49 49 USE dom_oce ! ocean space and time domain … … 64 64 PUBLIC ice_var_eqv2glo 65 65 PUBLIC ice_var_salprof 66 PUBLIC ice_var_bv67 66 PUBLIC ice_var_salprof1d 68 67 PUBLIC ice_var_zapsmall 69 68 PUBLIC ice_var_itd 69 PUBLIC ice_var_bv 70 70 71 71 !!---------------------------------------------------------------------- … … 170 170 END WHERE 171 171 ! 172 WHERE( v_i(:,:,:) > epsi20 ) ; z1_v_i(:,:,:) = 1._wp / v_i(:,:,:) 173 ELSEWHERE ; z1_v_i(:,:,:) = 0._wp 174 END WHERE 175 ! 172 176 ht_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:) !--- ice thickness 173 177 … … 177 181 ht_i (:,:,jpl) = zhmax 178 182 a_i (:,:,jpl) = v_i(:,:,jpl) * z1_zhmax 179 z1_a_i(:,:,jpl) = zhmax /v_i(:,:,jpl) ! NB: v_i always /=0 as ht_i > hi_max183 z1_a_i(:,:,jpl) = zhmax * z1_v_i(:,:,jpl) ! NB: v_i always /=0 as ht_i > hi_max 180 184 END WHERE 181 185 … … 185 189 186 190 IF( nn_icesal == 2 ) THEN !--- salinity (with a minimum value imposed everywhere) 187 WHERE( v_i(:,:,:) > epsi20 ) ; sm_i(:,:,:) = MAX( rn_simin , smv_i(:,:,:) / v_i(:,:,:) )191 WHERE( v_i(:,:,:) > epsi20 ) ; sm_i(:,:,:) = MAX( rn_simin , MIN( rn_simax, smv_i(:,:,:) * z1_v_i(:,:,:) ) ) 188 192 ELSEWHERE ; sm_i(:,:,:) = rn_simin 189 193 END WHERE … … 202 206 IF ( v_i(ji,jj,jl) > epsi20 ) THEN !--- icy area 203 207 ! 204 ze_i = e_i(ji,jj,jk,jl) /v_i(ji,jj,jl) * zlay_i ! Energy of melting e(S,T) [J.m-3]208 ze_i = e_i(ji,jj,jk,jl) * z1_v_i(ji,jj,jl) * zlay_i ! Energy of melting e(S,T) [J.m-3] 205 209 ztmelts = - s_i(ji,jj,jk,jl) * tmut ! Ice layer melt temperature [C] 206 210 ! Conversion q(S,T) -> T (second order equation) … … 295 299 ALLOCATE( z_slope_s(jpi,jpj,jpl) , zalpha(jpi,jpj,jpl) ) 296 300 ! 297 DO jk = 1, nlay_i 298 s_i(:,:,jk,:) = sm_i(:,:,:) 301 DO jl = 1, jpl 302 DO jk = 1, nlay_i 303 s_i(:,:,jk,jl) = sm_i(:,:,jl) 304 END DO 299 305 END DO 300 306 ! ! Slope of the linear profile … … 353 359 END SUBROUTINE ice_var_salprof 354 360 355 356 SUBROUTINE ice_var_bv357 !!-------------------------------------------------------------------358 !! *** ROUTINE ice_var_bv ***359 !!360 !! ** Purpose : computes mean brine volume (%) in sea ice361 !!362 !! ** Method : e = - 0.054 * S (ppt) / T (C)363 !!364 !! References : Vancoppenolle et al., JGR, 2007365 !!-------------------------------------------------------------------366 INTEGER :: ji, jj, jk, jl ! dummy loop indices367 !!-------------------------------------------------------------------368 !369 !!gm I prefere to use WHERE / ELSEWHERE to set it to zero only where needed <<<=== to be done370 !! instead of setting everything to zero as just below371 bv_i (:,:,:) = 0._wp372 DO jl = 1, jpl373 DO jk = 1, nlay_i374 WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 )375 bv_i(:,:,jl) = bv_i(:,:,jl) - tmut * s_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 )376 END WHERE377 END DO378 END DO379 WHERE( vt_i(:,:) > epsi20 ) ; bvm_i(:,:) = SUM( bv_i(:,:,:) * v_i(:,:,:) , dim=3 ) / vt_i(:,:)380 ELSEWHERE ; bvm_i(:,:) = 0._wp381 END WHERE382 !383 END SUBROUTINE ice_var_bv384 385 386 361 SUBROUTINE ice_var_salprof1d 387 362 !!------------------------------------------------------------------- … … 393 368 INTEGER :: ji, jk ! dummy loop indices 394 369 REAL(wp) :: zargtemp, zsal, z1_dS ! local scalars 395 REAL(wp) :: z alpha, zs, zs0 ! - -396 ! 397 REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_slope_s !370 REAL(wp) :: zs, zs0 ! - - 371 ! 372 REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_slope_s, zalpha ! 398 373 REAL(wp), PARAMETER :: zsi0 = 3.5_wp 399 374 REAL(wp), PARAMETER :: zsi1 = 4.5_wp … … 405 380 CASE( 1 ) ! constant salinity in time and space ! 406 381 ! !---------------------------------------! 407 s_i_1d( :,:) = rn_icesal382 s_i_1d(1:nidx,:) = rn_icesal 408 383 ! 409 384 ! !---------------------------------------------! … … 411 386 ! !---------------------------------------------! 412 387 ! 413 ALLOCATE( z_slope_s(jpij) )414 ! 415 DO ji = 1, nidx ! Slope of the linear profile416 rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi20 ))417 z_slope_s(ji) = rswitch * 2._wp * sm_i_1d(ji) / MAX( epsi20 , ht_i_1d(ji) )418 END DO419 388 ALLOCATE( z_slope_s(jpij), zalpha(jpij) ) 389 ! 390 ! ! Slope of the linear profile 391 WHERE( ht_i_1d(1:nidx) > epsi20 ) ; z_slope_s(1:nidx) = 2._wp * sm_i_1d(1:nidx) / ht_i_1d(1:nidx) 392 ELSEWHERE ; z_slope_s(1:nidx) = 0._wp 393 END WHERE 394 420 395 z1_dS = 1._wp / ( zsi1 - zsi0 ) 396 DO ji = 1, nidx 397 zalpha(ji) = MAX( 0._wp , MIN( ( zsi1 - sm_i_1d(ji) ) * z1_dS , 1._wp ) ) 398 ! ! force a constant profile when SSS too low (Baltic Sea) 399 IF( 2._wp * sm_i_1d(ji) >= sss_1d(ji) ) zalpha(ji) = 0._wp 400 END DO 401 ! 402 ! Computation of the profile 421 403 DO jk = 1, nlay_i 422 404 DO ji = 1, nidx 423 zalpha = MAX( 0._wp , MIN( ( zsi1 - sm_i_1d(ji) ) * z1_dS , 1._wp ) ) 424 IF( 2._wp * sm_i_1d(ji) >= sss_1d(ji) ) zalpha = 0._wp ! cst profile when SSS too low (Baltic Sea) 425 ! 426 zs0 = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_1d(ji) * r1_nlay_i ! linear profile with 0 surface value 427 zs = zalpha * zs0 + ( 1._wp - zalpha ) * sm_i_1d(ji) ! weighting the profile 428 s_i_1d(ji,jk) = MIN( rn_simax , MAX( zs , rn_simin ) ) ! bounding salinity 429 END DO 430 END DO 431 ! 432 DEALLOCATE( z_slope_s ) 405 ! ! linear profile with 0 surface value 406 zs0 = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_1d(ji) * r1_nlay_i 407 zs = zalpha(ji) * zs0 + ( 1._wp - zalpha(ji) ) * sm_i_1d(ji) 408 s_i_1d(ji,jk) = MIN( rn_simax , MAX( zs , rn_simin ) ) 409 END DO 410 END DO 411 ! 412 DEALLOCATE( z_slope_s, zalpha ) 433 413 434 414 ! !-------------------------------------------! … … 436 416 ! !-------------------------------------------! (mean = 2.30) 437 417 ! 438 sm_i_1d( :) = 2.30_wp418 sm_i_1d(1:nidx) = 2.30_wp 439 419 ! 440 420 !!gm cf remark in ice_var_salprof routine, CASE( 3 ) … … 459 439 !!------------------------------------------------------------------- 460 440 INTEGER :: ji, jj, jl, jk ! dummy loop indices 461 REAL(wp) :: zsal, zvi, zvs, zei, zes, zvp441 REAL(wp), DIMENSION(jpi,jpj) :: zswitch 462 442 !!------------------------------------------------------------------- 463 443 ! … … 467 447 ! Zap ice energy and use ocean heat to melt ice 468 448 !----------------------------------------------------------------- 449 WHERE( a_i(:,:,jl) > epsi20 ) ; ht_i(:,:,jl) = v_i(:,:,jl) / a_i(:,:,jl) 450 ELSEWHERE ; ht_i(:,:,jl) = 0._wp 451 END WHERE 452 ! 453 WHERE( a_i(:,:,jl) < epsi10 .OR. v_i(:,:,jl) < epsi10 .OR. ht_i(:,:,jl) < epsi10 ) ; zswitch(:,:) = 0._wp 454 ELSEWHERE ; zswitch(:,:) = 1._wp 455 END WHERE 456 469 457 DO jk = 1, nlay_i 470 458 DO jj = 1 , jpj 471 459 DO ji = 1 , jpi 472 !!gm I think we can do better/faster : to be discussed473 rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) )474 rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi10 ) ) * rswitch475 rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) * rswitch &476 & / MAX( a_i(ji,jj,jl), epsi10 ) - epsi10 ) ) * rswitch477 zei = e_i(ji,jj,jk,jl)478 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * rswitch479 t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * rswitch + rt0 * ( 1._wp - rswitch )480 460 ! update exchanges with ocean 481 hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_i(ji,jj,jk,jl) - zei ) * r1_rdtice ! W.m-2 <0 461 hfx_res(ji,jj) = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_i(ji,jj,jk,jl) * r1_rdtice ! W.m-2 <0 462 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * zswitch(ji,jj) 463 t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) ) 482 464 END DO 483 465 END DO … … 486 468 DO jj = 1 , jpj 487 469 DO ji = 1 , jpi 488 rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) ) 489 rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi10 ) ) * rswitch 490 rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) * rswitch & 491 & / MAX( a_i(ji,jj,jl), epsi10 ) - epsi10 ) ) * rswitch 492 zsal = smv_i(ji,jj, jl) 493 zvi = v_i (ji,jj, jl) 494 zvs = v_s (ji,jj, jl) 495 zes = e_s (ji,jj,1,jl) 496 IF ( ( nn_pnd_scheme > 0 ) .AND. ln_pnd_fw ) zvp = v_ip (ji,jj ,jl) 470 ! update exchanges with ocean 471 sfx_res(ji,jj) = sfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * smv_i(ji,jj,jl) * rhoic * r1_rdtice 472 wfx_res(ji,jj) = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_i (ji,jj,jl) * rhoic * r1_rdtice 473 wfx_res(ji,jj) = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl) * rhosn * r1_rdtice 474 hfx_res(ji,jj) = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_s (ji,jj,1,jl) * r1_rdtice ! W.m-2 <0 497 475 !----------------------------------------------------------------- 498 476 ! Zap snow energy 499 477 !----------------------------------------------------------------- 500 t_s(ji,jj,1,jl) = t_s(ji,jj,1,jl) * rswitch + rt0 * ( 1._wp - rswitch)501 e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * rswitch478 t_s(ji,jj,1,jl) = t_s(ji,jj,1,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) ) 479 e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * zswitch(ji,jj) 502 480 503 481 !----------------------------------------------------------------- 504 482 ! zap ice and snow volume, add water and salt to ocean 505 483 !----------------------------------------------------------------- 506 ato_i(ji,jj) = a_i (ji,jj,jl) * ( 1._wp - rswitch ) + ato_i(ji,jj) 507 a_i (ji,jj,jl) = a_i (ji,jj,jl) * rswitch 508 v_i (ji,jj,jl) = v_i (ji,jj,jl) * rswitch 509 v_s (ji,jj,jl) = v_s (ji,jj,jl) * rswitch 510 t_su (ji,jj,jl) = t_su (ji,jj,jl) * rswitch + t_bo(ji,jj) * ( 1._wp - rswitch ) 511 oa_i (ji,jj,jl) = oa_i (ji,jj,jl) * rswitch 512 smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * rswitch 513 514 ht_i (ji,jj,jl) = ht_i (ji,jj,jl) * rswitch 515 ht_s (ji,jj,jl) = ht_s (ji,jj,jl) * rswitch 516 517 ! MV MP 2016 518 IF ( ln_pnd ) THEN 519 a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * rswitch 520 v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * rswitch 521 IF ( ln_pnd_fw ) wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_ip(ji,jj,jl) - zvp ) * rhofw * r1_rdtice 522 ENDIF 523 ! END MV MP 2016 524 525 ! update exchanges with ocean 526 sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice 527 wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice 528 wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice 529 hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * r1_rdtice ! W.m-2 <0 530 END DO 531 END DO 484 ato_i(ji,jj) = a_i (ji,jj,jl) * ( 1._wp - zswitch(ji,jj) ) + ato_i(ji,jj) 485 a_i (ji,jj,jl) = a_i (ji,jj,jl) * zswitch(ji,jj) 486 v_i (ji,jj,jl) = v_i (ji,jj,jl) * zswitch(ji,jj) 487 v_s (ji,jj,jl) = v_s (ji,jj,jl) * zswitch(ji,jj) 488 t_su (ji,jj,jl) = t_su (ji,jj,jl) * zswitch(ji,jj) + t_bo(ji,jj) * ( 1._wp - zswitch(ji,jj) ) 489 oa_i (ji,jj,jl) = oa_i (ji,jj,jl) * zswitch(ji,jj) 490 smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * zswitch(ji,jj) 491 492 ht_i (ji,jj,jl) = ht_i (ji,jj,jl) * zswitch(ji,jj) 493 ht_s (ji,jj,jl) = ht_s (ji,jj,jl) * zswitch(ji,jj) 494 495 END DO 496 END DO 497 498 IF( ln_pnd ) THEN 499 DO jj = 1 , jpj 500 DO ji = 1 , jpi 501 IF( ln_pnd_fw ) & 502 & wfx_res(ji,jj) = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_ip(ji,jj,jl) * rhofw * r1_rdtice 503 a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * zswitch(ji,jj) 504 v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj) 505 END DO 506 END DO 507 ENDIF 508 532 509 END DO 533 510 … … 688 665 END SUBROUTINE ice_var_itd 689 666 667 668 SUBROUTINE ice_var_bv 669 !!------------------------------------------------------------------- 670 !! *** ROUTINE ice_var_bv *** 671 !! 672 !! ** Purpose : computes mean brine volume (%) in sea ice 673 !! 674 !! ** Method : e = - 0.054 * S (ppt) / T (C) 675 !! 676 !! References : Vancoppenolle et al., JGR, 2007 677 !!------------------------------------------------------------------- 678 INTEGER :: ji, jj, jk, jl ! dummy loop indices 679 !!------------------------------------------------------------------- 680 ! 681 !!gm I prefere to use WHERE / ELSEWHERE to set it to zero only where needed <<<=== to be done 682 !! instead of setting everything to zero as just below 683 bv_i (:,:,:) = 0._wp 684 DO jl = 1, jpl 685 DO jk = 1, nlay_i 686 WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 ) 687 bv_i(:,:,jl) = bv_i(:,:,jl) - tmut * s_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 ) 688 END WHERE 689 END DO 690 END DO 691 WHERE( vt_i(:,:) > epsi20 ) ; bvm_i(:,:) = SUM( bv_i(:,:,:) * v_i(:,:,:) , dim=3 ) / vt_i(:,:) 692 ELSEWHERE ; bvm_i(:,:) = 0._wp 693 END WHERE 694 ! 695 END SUBROUTINE ice_var_bv 696 697 690 698 #else 691 699 !!----------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.