Changeset 6763
- Timestamp:
- 2016-06-30T17:17:35+02:00 (9 years ago)
- Location:
- branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM
- Files:
-
- 8 edited
-
CONFIG/SHARED/namelist_ice_lim3_ref (modified) (1 diff)
-
CONFIG/SHARED/namelist_ref (modified) (1 diff)
-
NEMO/LIM_SRC_3/ice.F90 (modified) (2 diffs)
-
NEMO/LIM_SRC_3/limrhg.F90 (modified) (17 diffs)
-
NEMO/LIM_SRC_3/limsbc.F90 (modified) (1 diff)
-
NEMO/LIM_SRC_3/limwri.F90 (modified) (1 diff)
-
NEMO/OPA_SRC/SBC/sbcblk_core.F90 (modified) (18 diffs)
-
NEMO/OPA_SRC/SBC/sbcice_lim.F90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref
r6747 r6763 90 90 rn_gamma = 0.2 ! (ln_landfast=T) fraction of ocean depth that ice must reach to initiate landfast 91 91 ! recommended range: [0.1 ; 0.25] 92 rn_icebfr = 1.e-05 ! (ln_landfast=T) maximum bottom stress per unit area of contact 92 rn_icebfr = 5.6e-06 ! (ln_landfast=T) maximum bottom stress per unit area of contact 93 ! a very large value ensures ice velocity=0 even with a small contact area 93 94 ! recommended range: ?? 94 95 / -
branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/CONFIG/SHARED/namelist_ref
r6319 r6763 336 336 rn_vfac = 0. ! multiplicative factor for ocean/ice velocity 337 337 ! in the calculation of the wind stress (0.=absolute winds or 1.=relative winds) 338 ln_Cd_L12 = .false. ! Modify the drag ice-atm and oce-atm depending on ice concentration 339 ! This parameterization is from Lupkes et al. (JGR 2012) 338 340 / 339 341 !----------------------------------------------------------------------- -
branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
r6746 r6763 422 422 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: e_i_b !: ice temperatures 423 423 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: u_ice_b, v_ice_b !: ice velocity 424 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: at_i_b !: ice concentration (total) 424 425 425 426 !!-------------------------------------------------------------------------- … … 528 529 & oa_i_b (jpi,jpj,jpl) , STAT=ierr(ii) ) 529 530 ii = ii + 1 530 ALLOCATE( u_ice_b(jpi,jpj) , v_ice_b(jpi,jpj) , STAT=ierr(ii) )531 ALLOCATE( u_ice_b(jpi,jpj) , v_ice_b(jpi,jpj) , at_i_b(jpi,jpj) , STAT=ierr(ii) ) 531 532 532 533 ! * Ice thickness distribution variables -
branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90
r6746 r6763 116 116 REAL(wp) :: zbeta, zalph1, z1_alph1, zalph2, z1_alph2 ! alpha and beta from Bouillon 2009 and 2013 117 117 REAL(wp) :: zm1, zm2, zm3 ! ice/snow mass 118 REAL(wp) :: delta, zp_delf, zds2118 REAL(wp) :: zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2 119 119 REAL(wp) :: zTauO, zTauA, zTauB, ZCor, zmt, zu_ice2, zv_ice1, zvel ! temporary scalars 120 120 … … 134 134 REAL(wp), POINTER, DIMENSION(:,:) :: zf1 , zf2 ! internal stresses 135 135 136 REAL(wp), POINTER, DIMENSION(:,:) :: zd t, zds ! tension andshear136 REAL(wp), POINTER, DIMENSION(:,:) :: zds ! shear 137 137 REAL(wp), POINTER, DIMENSION(:,:) :: zs1, zs2, zs12 ! stress tensor components 138 138 REAL(wp), POINTER, DIMENSION(:,:) :: zu_ice, zv_ice, zresr ! check convergence … … 149 149 CALL wrk_alloc( jpi,jpj, za1, za2, zmass1, zmass2, zcor1, zcor2 ) 150 150 CALL wrk_alloc( jpi,jpj, zspgu, zspgv, v_oce1, u_oce2, v_ice1, u_ice2, zf1, zf2 ) 151 CALL wrk_alloc( jpi,jpj, zds, z dt, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice )151 CALL wrk_alloc( jpi,jpj, zds, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice ) 152 152 CALL wrk_alloc( jpi,jpj, zswitch1, zswitch2 ) 153 153 … … 257 257 END DO 258 258 END DO 259 260 259 ! 261 260 !------------------------------------------------------------------------------! … … 266 265 DO jter = 1 , nn_nevp ! loop over jter ! 267 266 ! !----------------------! 268 ! Convergence test 269 IF(ln_ctl) THEN 267 IF(ln_ctl) THEN ! Convergence test 270 268 DO jj = 1, jpjm1 271 269 zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step … … 275 273 276 274 ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- ! 277 DO jj = 2, jpjm1 278 DO ji = fs_2, fs_jpim1 279 280 ! divergence at T points 281 divu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) & 282 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) & 283 & ) * r1_e12t(ji,jj) 284 285 ! tension at T points 286 zdt(ji,jj) = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) & 287 & - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) & 288 & ) * r1_e12t(ji,jj) 275 DO jj = 1, jpjm1 ! loops start at 1 since there is no boundary condition (lbc_lnk) at i=1 and j=1 for F points 276 DO ji = 1, fs_jpim1 289 277 290 278 ! shear at F points … … 292 280 & + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 293 281 & ) * r1_e12f(ji,jj) 282 283 END DO 284 END DO 285 CALL lbc_lnk( zds, 'F', 1. ) 286 287 DO jj = 2, jpjm1 288 DO ji = 2, jpim1 ! no vector loop 289 290 ! shear**2 at T points (doc eq. A16) 291 zds2 = ( zds(ji,jj ) * zds(ji,jj ) * e12f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e12f(ji-1,jj ) & 292 & + zds(ji,jj-1) * zds(ji,jj-1) * e12f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e12f(ji-1,jj-1) & 293 & ) * 0.25_wp * r1_e12t(ji,jj) 294 295 ! divergence at T points 296 zdiv = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) & 297 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) & 298 & ) * r1_e12t(ji,jj) 299 zdiv2 = zdiv * zdiv 300 301 ! tension at T points 302 zdt = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) & 303 & - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) & 304 & ) * r1_e12t(ji,jj) 305 zdt2 = zdt * zdt 306 307 ! delta at T points 308 zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * usecc2 ) 309 310 ! P/delta at T points 311 zp_delt(ji,jj) = zpresh(ji,jj) / ( zdelta + rn_creepl ) 312 313 ! stress at T points 314 zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv - zdelta ) ) * z1_alph1 315 zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 ) ) * z1_alph2 316 317 END DO 318 END DO 319 CALL lbc_lnk( zp_delt, 'T', 1. ) 320 321 DO jj = 1, jpjm1 322 DO ji = 1, jpim1 323 324 ! P/delta at F points 325 zp_delf = 0.25_wp * ( zp_delt(ji,jj) + zp_delt(ji+1,jj) + zp_delt(ji,jj+1) + zp_delt(ji+1,jj+1) ) 326 327 ! stress at F points 328 zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 ) * 0.5_wp ) * z1_alph2 329 330 END DO 331 END DO 332 CALL lbc_lnk_multi( zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. ) 333 334 335 ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- ! 336 DO jj = 2, jpjm1 337 DO ji = fs_2, fs_jpim1 338 339 ! U points 340 zf1(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj) & 341 & + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj) & 342 & ) * r1_e2u(ji,jj) & 343 & + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) & 344 & ) * 2._wp * r1_e1u(ji,jj) & 345 & ) * r1_e12u(ji,jj) 346 347 ! V points 348 zf2(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj) & 349 & - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj) & 350 & ) * r1_e1v(ji,jj) & 351 & + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) & 352 & ) * 2._wp * r1_e2v(ji,jj) & 353 & ) * r1_e12v(ji,jj) 294 354 295 355 ! u_ice at V point … … 303 363 END DO 304 364 END DO 305 CALL lbc_lnk( zds, 'F', 1. )306 307 DO jj = 2, jpjm1308 DO ji = 2, jpim1 ! no vector loop309 310 ! shear**2 at T points (doc eq. A16)311 zds2 = ( zds(ji,jj ) * zds(ji,jj ) * e12f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e12f(ji-1,jj ) &312 & + zds(ji,jj-1) * zds(ji,jj-1) * e12f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e12f(ji-1,jj-1) &313 & ) * 0.25_wp * r1_e12t(ji,jj)314 315 ! delta at T points316 delta = SQRT( divu_i(ji,jj) * divu_i(ji,jj) + ( zdt(ji,jj) * zdt(ji,jj) + zds2 ) * usecc2 )317 delta_i(ji,jj) = delta + rn_creepl318 319 ! P/delta at T points320 zp_delt(ji,jj) = zpresh(ji,jj) / delta_i(ji,jj)321 END DO322 END DO323 CALL lbc_lnk( zp_delt, 'T', 1. )324 325 ! --- Stress tensor --- !326 DO jj = 2, jpjm1327 DO ji = 2, jpim1 ! no vector loop328 329 ! P/delta at F points330 zp_delf = 0.25_wp * ( zp_delt(ji,jj) + zp_delt(ji+1,jj) + zp_delt(ji,jj+1) + zp_delt(ji+1,jj+1) )331 332 ! stress tensor at T points333 zs1(ji,jj) = ( zs1 (ji,jj) * zalph1 + zp_delt(ji,jj) * ( divu_i(ji,jj) - (delta_i(ji,jj)-rn_creepl) ) ) * z1_alph1334 zs2(ji,jj) = ( zs2 (ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt(ji,jj) * z1_ecc2 ) ) * z1_alph2335 ! F points336 zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 ) * 0.5_wp ) * z1_alph2337 END DO338 END DO339 CALL lbc_lnk_multi( zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. )340 341 ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- !342 DO jj = 2, jpjm1343 DO ji = fs_2, fs_jpim1344 ! U points345 zf1(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj) &346 & + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj) &347 & ) * r1_e2u(ji,jj) &348 & + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) &349 & ) * 2._wp * r1_e1u(ji,jj) &350 & ) * r1_e12u(ji,jj)351 352 ! V points353 zf2(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj) &354 & - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj) &355 & ) * r1_e1v(ji,jj) &356 & + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) &357 & ) * 2._wp * r1_e2v(ji,jj) &358 & ) * r1_e12v(ji,jj)359 END DO360 END DO361 365 ! 362 366 ! --- Computation of ice velocity --- ! … … 368 372 DO ji = fs_2, fs_jpim1 369 373 370 zvel = SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) ) & 371 & + ( u_ice2(ji,jj) - u_oce2(ji,jj) ) * ( u_ice2(ji,jj) - u_oce2(ji,jj) ) ) 372 374 zvel = SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_ice2(ji,jj) * u_ice2(ji,jj) ) 375 373 376 zTauA = za2(ji,jj) * vtau_ice(ji,jj) 374 zTauO = za2(ji,jj) * rhoco * zvel 375 ztauB = - tau_icebfr(ji,jj) / zvel 377 zTauO = za2(ji,jj) * rhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) ) & 378 & + ( u_ice2(ji,jj) - u_oce2(ji,jj) ) * ( u_ice2(ji,jj) - u_oce2(ji,jj) ) ) 379 ztauB = - tau_icebfr(ji,jj) / MAX( zvel, zepsi ) 376 380 zCor = zcor2(ji,jj) * u_ice2(ji,jj) 377 381 zmt = zmass2(ji,jj) * z1_dtevp … … 403 407 & + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1) 404 408 405 zvel = SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) ) & 406 & + ( v_ice1(ji,jj) - v_oce1(ji,jj) ) * ( v_ice1(ji,jj) - v_oce1(ji,jj) ) ) 407 409 zvel = SQRT( v_ice1(ji,jj) * v_ice1(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) 410 408 411 zTauA = za1(ji,jj) * utau_ice(ji,jj) 409 zTauO = za1(ji,jj) * rhoco * zvel 410 ztauB = - tau_icebfr(ji,jj) / zvel 412 zTauO = za1(ji,jj) * rhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) ) & 413 & + ( v_ice1(ji,jj) - v_oce1(ji,jj) ) * ( v_ice1(ji,jj) - v_oce1(ji,jj) ) ) 414 ztauB = - tau_icebfr(ji,jj) / MAX( zvel, zepsi ) 411 415 zCor = zcor1(ji,jj) * zv_ice1 412 416 zmt = zmass1(ji,jj) * z1_dtevp … … 435 439 DO ji = fs_2, fs_jpim1 436 440 437 zvel = SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) ) & 438 & + ( v_ice1(ji,jj) - v_oce1(ji,jj) ) * ( v_ice1(ji,jj) - v_oce1(ji,jj) ) ) 441 zvel = SQRT( v_ice1(ji,jj) * v_ice1(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) 439 442 440 443 zTauA = za1(ji,jj) * utau_ice(ji,jj) 441 zTauO = za1(ji,jj) * rhoco * zvel 442 ztauB = - tau_icebfr(ji,jj) / zvel 444 zTauO = za1(ji,jj) * rhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) ) & 445 & + ( v_ice1(ji,jj) - v_oce1(ji,jj) ) * ( v_ice1(ji,jj) - v_oce1(ji,jj) ) ) 446 ztauB = - tau_icebfr(ji,jj) / MAX( zvel, zepsi ) 443 447 zCor = zcor1(ji,jj) * v_ice1(ji,jj) 444 448 zmt = zmass1(ji,jj) * z1_dtevp … … 469 473 & + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1) 470 474 471 zvel = SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) ) & 472 & + ( u_ice2(ji,jj) - u_oce2(ji,jj) ) * ( u_ice2(ji,jj) - u_oce2(ji,jj) ) ) 475 zvel = SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_ice2(ji,jj) * u_ice2(ji,jj) ) 473 476 474 477 zTauA = za2(ji,jj) * vtau_ice(ji,jj) 475 zTauO = za2(ji,jj) * rhoco * zvel 476 ztauB = - tau_icebfr(ji,jj) / zvel 478 zTauO = za2(ji,jj) * rhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) ) & 479 & + ( u_ice2(ji,jj) - u_oce2(ji,jj) ) * ( u_ice2(ji,jj) - u_oce2(ji,jj) ) ) 480 ztauB = - tau_icebfr(ji,jj) / MAX( zvel, zepsi ) 477 481 zCor = zcor2(ji,jj) * zu_ice2 478 482 zmt = zmass2(ji,jj) * z1_dtevp … … 485 489 !! & + zf2(ji,jj) + zCor + zTauA + zTauO * v_oce(ji,jj) + zspgv(ji,jj) & 486 490 !! & ) / MAX( zmt * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitch2(ji,jj) 487 488 491 END DO 489 492 END DO … … 499 502 ENDIF 500 503 501 ! Convergence test 502 IF(ln_ctl) THEN 504 IF(ln_ctl) THEN ! Convergence test 503 505 DO jj = 2 , jpjm1 504 506 zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) ) … … 539 541 ! 5) Recompute delta, shear and div (inputs for mechanical redistribution) 540 542 !------------------------------------------------------------------------------! 541 542 ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- ! 543 DO jj = 1, jpjm1 544 DO ji = 1, fs_jpim1 545 546 ! shear at F points 547 zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) & 548 & + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 549 & ) * r1_e12f(ji,jj) 550 551 END DO 552 END DO 553 CALL lbc_lnk( zds, 'F', 1. ) 554 543 555 DO jj = 2, jpjm1 544 DO ji = fs_2, fs_jpim1 545 556 DO ji = 2, jpim1 ! no vector loop 557 558 ! tension**2 at T points 559 zdt = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) & 560 & - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) & 561 & ) * r1_e12t(ji,jj) 562 zdt2 = zdt * zdt 563 564 ! shear**2 at T points (doc eq. A16) 565 zds2 = ( zds(ji,jj ) * zds(ji,jj ) * e12f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e12f(ji-1,jj ) & 566 & + zds(ji,jj-1) * zds(ji,jj-1) * e12f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e12f(ji-1,jj-1) & 567 & ) * 0.25_wp * r1_e12t(ji,jj) 568 569 ! shear at T points 570 shear_i(ji,jj) = SQRT( zdt2 + zds2 ) 571 546 572 ! divergence at T points 547 573 divu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) & … … 549 575 & ) * r1_e12t(ji,jj) 550 576 551 ! tension at T points 552 zdt(ji,jj) = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) & 553 & - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) & 554 & ) * r1_e12t(ji,jj) 555 556 ! shear at F points 557 zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) & 558 & + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 559 & ) * r1_e12f(ji,jj) 560 577 ! delta at T points 578 zdelta = SQRT( divu_i(ji,jj) * divu_i(ji,jj) + ( zdt2 + zds2 ) * usecc2 ) 579 delta_i(ji,jj) = zdelta + rn_creepl 580 561 581 END DO 562 582 END DO 563 CALL lbc_lnk_multi( divu_i, 'T', 1., zds, 'F', 1. ) 564 565 566 DO jj = 2, jpjm1 567 DO ji = 2, jpim1 ! no vector loop 568 569 ! shear**2 at T points (doc eq. A16) 570 zds2 = ( zds(ji,jj ) * zds(ji,jj ) * e12f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e12f(ji-1,jj ) & 571 & + zds(ji,jj-1) * zds(ji,jj-1) * e12f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e12f(ji-1,jj-1) & 572 & ) * 0.25_wp * r1_e12t(ji,jj) 573 574 ! delta at T points 575 delta = SQRT( divu_i(ji,jj)**2 + ( zdt(ji,jj)**2 + zds2 ) * usecc2 ) 576 delta_i(ji,jj) = delta + rn_creepl 577 578 shear_i(ji,jj) = SQRT( zdt(ji,jj) * zdt(ji,jj) + zds2 ) 579 END DO 580 END DO 581 CALL lbc_lnk_multi( delta_i, 'T', 1., shear_i, 'T', 1. ) 583 CALL lbc_lnk_multi( shear_i, 'T', 1., divu_i, 'T', 1., delta_i, 'T', 1. ) 582 584 583 585 ! --- Store the stress tensor for the next time step --- ! … … 626 628 CALL wrk_dealloc( jpi,jpj, za1, za2, zmass1, zmass2, zcor1, zcor2 ) 627 629 CALL wrk_dealloc( jpi,jpj, zspgu, zspgv, v_oce1, u_oce2, v_ice1, u_ice2, zf1, zf2 ) 628 CALL wrk_dealloc( jpi,jpj, zds, z dt, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice )630 CALL wrk_dealloc( jpi,jpj, zds, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice ) 629 631 CALL wrk_dealloc( jpi,jpj, zswitch1, zswitch2 ) 630 632 -
branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90
r6515 r6763 137 137 138 138 zalb(:,:) = 0._wp 139 WHERE ( SUM( a_i_b, dim=3 )<= epsi06 ) ; zalb(:,:) = 0.066_wp140 ELSEWHERE ; zalb(:,:) = SUM( alb_ice * a_i_b, dim=3 ) / SUM( a_i_b, dim=3 )139 WHERE ( at_i_b <= epsi06 ) ; zalb(:,:) = 0.066_wp 140 ELSEWHERE ; zalb(:,:) = SUM( alb_ice * a_i_b, dim=3 ) / at_i_b 141 141 END WHERE 142 142 IF( iom_use('alb_ice' ) ) CALL iom_put( "alb_ice" , zalb(:,:) ) ! ice albedo output 143 143 144 zalb(:,:) = SUM( alb_ice * a_i_b, dim=3 ) + 0.066_wp * ( 1._wp - SUM( a_i_b, dim=3 ))144 zalb(:,:) = SUM( alb_ice * a_i_b, dim=3 ) + 0.066_wp * ( 1._wp - at_i_b ) 145 145 IF( iom_use('albedo' ) ) CALL iom_put( "albedo" , zalb(:,:) ) ! ice albedo output 146 146 -
branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90
r6746 r6763 185 185 CALL iom_put ('hfxdif' , hfx_dif(:,:) ) ! 186 186 CALL iom_put ('hfxopw' , hfx_opw(:,:) ) ! 187 CALL iom_put ('hfxtur' , fhtur(:,:) * SUM(a_i_b(:,:,:), dim=3) ) ! turbulent heat flux at ice base187 CALL iom_put ('hfxtur' , fhtur(:,:) * at_i_b(:,:) ) ! turbulent heat flux at ice base 188 188 CALL iom_put ('hfxdhc' , diag_heat(:,:) ) ! Heat content variation in snow and ice 189 189 CALL iom_put ('hfxspr' , hfx_spr(:,:) ) ! Heat content of snow precip -
branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90
r6399 r6763 16 16 !! 3.4 ! 2011-11 (C. Harris) Fill arrays required by CICE 17 17 !! 3.7 ! 2014-06 (L. Brodeau) simplification and optimization of CORE bulk 18 !! 4.0 ! 2016-06 (C. Rousset) Add new param of drags with sea-ice (Lupkes at al 2012) 18 19 !!---------------------------------------------------------------------- 19 20 … … 45 46 USE lib_fortran ! to use key_nosignedzero 46 47 #if defined key_lim3 47 USE ice, ONLY : u_ice, v_ice, jpl, pfrld, a_i_b 48 USE ice, ONLY : u_ice, v_ice, jpl, pfrld, a_i_b, at_i_b 48 49 USE limthd_dh ! for CALL lim_thd_snwblow 49 50 #elif defined key_lim2 … … 60 61 PUBLIC blk_ice_core_flx ! routine called in sbc_ice_lim module 61 62 #endif 62 PUBLIC turb_core_2z ! routine calle sin sbcblk_mfs module63 PUBLIC turb_core_2z ! routine called in sbcblk_mfs module 63 64 64 65 INTEGER , PARAMETER :: jpfld = 9 ! maximum number of files to read … … 76 77 77 78 ! !!! CORE bulk parameters 78 REAL(wp), PARAMETER :: rhoa = 1.22 ! air density79 REAL(wp), PARAMETER :: cpa = 1000.5 ! specific heat of air80 REAL(wp), PARAMETER :: Lv = 2.5e6 ! latent heat of vaporization81 REAL(wp), PARAMETER :: Ls = 2.839e6 ! latent heat of sublimation82 REAL(wp), PARAMETER :: Stef = 5.67e-8 ! Stefan Boltzmann constant83 REAL(wp), PARAMETER :: C ice = 1.4e-3 ! iovi 1.63e-3! transfer coefficient over ice84 REAL(wp), PARAMETER :: albo = 0.066 ! ocean albedo assumed to be constant79 REAL(wp), PARAMETER :: rhoa = 1.22 ! air density 80 REAL(wp), PARAMETER :: cpa = 1000.5 ! specific heat of air 81 REAL(wp), PARAMETER :: Lv = 2.5e6 ! latent heat of vaporization 82 REAL(wp), PARAMETER :: Ls = 2.839e6 ! latent heat of sublimation 83 REAL(wp), PARAMETER :: Stef = 5.67e-8 ! Stefan Boltzmann constant 84 REAL(wp), PARAMETER :: Cd_ice = 1.4e-3 ! transfer coefficient over ice 85 REAL(wp), PARAMETER :: albo = 0.066 ! ocean albedo assumed to be constant 85 86 86 87 ! !!* Namelist namsbc_core : CORE bulk parameters … … 91 92 REAL(wp) :: rn_zqt ! z(q,t) : height of humidity and temperature measurements 92 93 REAL(wp) :: rn_zu ! z(u) : height of wind measurements 93 94 LOGICAL :: ln_Cd_L12 = .FALSE. ! Modify the drag ice-atm and oce-atm depending on ice concentration (from Lupkes et al. JGR2012) 95 96 ! 97 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: Cd_oce ! air-ocean drag (clem) 98 94 99 !! * Substitutions 95 100 # include "domzgr_substitute.h90" … … 102 107 CONTAINS 103 108 109 INTEGER FUNCTION sbc_blk_core_alloc() 110 !!------------------------------------------------------------------- 111 !! *** ROUTINE sbc_blk_core_alloc (clem) *** 112 !!------------------------------------------------------------------- 113 ALLOCATE( Cd_oce(jpi,jpj) , STAT=sbc_blk_core_alloc ) 114 ! 115 IF( lk_mpp ) CALL mpp_sum( sbc_blk_core_alloc ) 116 IF( sbc_blk_core_alloc /= 0 ) CALL ctl_warn('sbc_blk_core_alloc: failed to allocate arrays') 117 END FUNCTION sbc_blk_core_alloc 118 119 104 120 SUBROUTINE sbc_blk_core( kt ) 105 121 !!--------------------------------------------------------------------- … … 151 167 & sn_wndi, sn_wndj, sn_humi , sn_qsr , & 152 168 & sn_qlw , sn_tair, sn_prec , sn_snow, & 153 & sn_tdif, rn_zqt, rn_zu 169 & sn_tdif, rn_zqt, rn_zu , ln_Cd_L12 154 170 !!--------------------------------------------------------------------- 155 171 ! … … 157 173 IF( kt == nit000 ) THEN ! First call kt=nit000 ! 158 174 ! ! ====================== ! 175 ! 176 ! ! allocate sbc_blk_core array (clem) 177 IF( sbc_blk_core_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_blk_core : unable to allocate standard arrays' ) 159 178 ! 160 179 REWIND( numnam_ref ) ! Namelist namsbc_core in reference namelist : CORE bulk parameters … … 319 338 & Cd, Ch, Ce, zt_zu, zq_zu ) 320 339 340 ! Make ocean-atm. drag dependent on ice concentration (see Lupkes et al. 2012) (clem) 341 #if defined key_lim3 342 IF( ln_Cd_L12 ) THEN 343 344 Cd_oce(:,:) = Cd(:,:) ! record value of pure ocean-atm. drag 345 346 CALL Cdn10_Lupkes2012( Cd ) ! calculate new drag from Lupkes(2012) equations 347 348 ENDIF 349 #endif 350 321 351 ! ... tau module, i and j component 322 352 DO jj = 1, jpj … … 437 467 !!--------------------------------------------------------------------- 438 468 INTEGER :: ji, jj ! dummy loop indices 439 REAL(wp) :: zcoef_wnorm, zcoef_wnorm2440 469 REAL(wp) :: zwnorm_f, zwndi_f , zwndj_f ! relative wind module and components at F-point 441 470 REAL(wp) :: zwndi_t , zwndj_t ! relative wind components at T-point 471 REAL(wp), DIMENSION(:,:), POINTER :: Cd ! transfer coefficient for momentum (tau) 442 472 !!--------------------------------------------------------------------- 443 473 ! 444 474 IF( nn_timing == 1 ) CALL timing_start('blk_ice_core_tau') 445 475 ! 446 ! local scalars ( place there for vector optimisation purposes) 447 zcoef_wnorm = rhoa * Cice 448 zcoef_wnorm2 = rhoa * Cice * 0.5 476 CALL wrk_alloc( jpi,jpj, Cd ) 477 478 Cd(:,:) = Cd_ice 479 480 ! Make ice-atm. drag dependent on ice concentration (see Lupkes et al. 2012) (clem) 481 #if defined key_lim3 482 IF( ln_Cd_L12 ) THEN 483 484 CALL Cdn10_Lupkes2012( Cd ) ! calculate new drag from Lupkes(2012) equations 485 486 ENDIF 487 #endif 449 488 450 489 !!gm brutal.... … … 467 506 zwndj_f = 0.25 * ( sf(jp_wndj)%fnow(ji-1,jj ,1) + sf(jp_wndj)%fnow(ji ,jj ,1) & 468 507 & + sf(jp_wndj)%fnow(ji-1,jj-1,1) + sf(jp_wndj)%fnow(ji ,jj-1,1) ) - rn_vfac * v_ice(ji,jj) 469 zwnorm_f = zcoef_wnorm* SQRT( zwndi_f * zwndi_f + zwndj_f * zwndj_f )508 zwnorm_f = rhoa * Cd(ji,jj) * SQRT( zwndi_f * zwndi_f + zwndj_f * zwndj_f ) 470 509 ! ... ice stress at I-point 471 510 utau_ice(ji,jj) = zwnorm_f * zwndi_f … … 476 515 zwndj_t = sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.25 * ( v_ice(ji,jj+1) + v_ice(ji+1,jj+1) & 477 516 & + v_ice(ji,jj ) + v_ice(ji+1,jj ) ) 478 wndm_ice(ji,jj) 517 wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1) 479 518 END DO 480 519 END DO … … 493 532 DO jj = 2, jpjm1 494 533 DO ji = fs_2, fs_jpim1 ! vect. opt. 495 utau_ice(ji,jj) = zcoef_wnorm2* ( wndm_ice(ji+1,jj ) + wndm_ice(ji,jj) ) &534 utau_ice(ji,jj) = rhoa * Cd(ji,jj) * 0.5_wp * ( wndm_ice(ji+1,jj ) + wndm_ice(ji,jj) ) & 496 535 & * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) - rn_vfac * u_ice(ji,jj) ) 497 vtau_ice(ji,jj) = zcoef_wnorm2* ( wndm_ice(ji,jj+1 ) + wndm_ice(ji,jj) ) &536 vtau_ice(ji,jj) = rhoa * Cd(ji,jj) * 0.5_wp * ( wndm_ice(ji,jj+1 ) + wndm_ice(ji,jj) ) & 498 537 & * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) - rn_vfac * v_ice(ji,jj) ) 499 538 END DO … … 510 549 ENDIF 511 550 551 CALL wrk_dealloc( jpi,jpj, Cd ) 552 512 553 IF( nn_timing == 1 ) CALL timing_stop('blk_ice_core_tau') 513 554 … … 548 589 ! local scalars ( place there for vector optimisation purposes) 549 590 zcoef_dqlw = 4.0 * 0.95 * Stef 550 zcoef_dqla = -Ls * C ice * 11637800. * (-5897.8)551 zcoef_dqsb = rhoa * cpa * C ice591 zcoef_dqla = -Ls * Cd_ice * 11637800. * (-5897.8) 592 zcoef_dqsb = rhoa * cpa * Cd_ice 552 593 553 594 zztmp = 1. / ( 1. - albo ) … … 575 616 ! ... turbulent heat fluxes 576 617 ! Sensible Heat 577 z_qsb(ji,jj,jl) = rhoa * cpa * C ice * wndm_ice(ji,jj) * ( ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1) )618 z_qsb(ji,jj,jl) = rhoa * cpa * Cd_ice * wndm_ice(ji,jj) * ( ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1) ) 578 619 ! Latent Heat 579 qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa * Ls * C ice * wndm_ice(ji,jj) &620 qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa * Ls * Cd_ice * wndm_ice(ji,jj) & 580 621 & * ( 11637800. * EXP( -5897.8 / ptsu(ji,jj,jl) ) / rhoa - sf(jp_humi)%fnow(ji,jj,1) ) ) 581 622 ! Latent heat sensitivity for ice (Dqla/Dt) … … 820 861 ! 821 862 END DO 822 863 823 864 CALL wrk_dealloc( jpi,jpj, U_zu, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd ) 824 865 CALL wrk_dealloc( jpi,jpj, zeta_u, stab ) … … 903 944 END FUNCTION psi_h 904 945 946 947 SUBROUTINE Cdn10_Lupkes2012( Cd ) 948 !!---------------------------------------------------------------------- 949 !! *** ROUTINE Cdn10_Lupkes2012 *** 950 !! 951 !! ** Purpose : Recompute the ice-atm and ocean-atm drags at 10m height to make 952 !! them dependent on edges at leads, melt ponds and flows. 953 !! After some approximations, this can be resumed to a dependency 954 !! on ice concentration. 955 !! 956 !! ** Method : The parameterization is taken from Lupkes et al. (2012) eq.(50) 957 !! with the highest level of approximation: level4, eq.(59) 958 !! The drag can be re-written as follows: 959 !! 960 !! Cd = Cdw * (1-A) + Cdi * A + Ce * (1-A)**(nu+1/(10*beta)) * A**mu 961 !! 962 !! Ce = 2.23e-3 , as suggested by Lupkes (eq. 59) 963 !! nu = mu = beta = 1 , as suggested by Lupkes (eq. 59) 964 !! A is the concentration of ice minus melt ponds (if any) 965 !! 966 !! This new drag has a parabolic shape (as a function of A) starting at 967 !! Cdw(say 1.5e-3) for A=0, reaching 1.97e-3 for A~0.5 968 !! and going down to Cdi(say 1.4e-3) for A=1 969 !! 970 !! It is theoretically applicable to all ice conditions (not only MIZ) 971 !! => see Lupkes et al (2013) 972 !! 973 !! ** References : Lupkes et al. JGR 2012 (theory) 974 !! Lupkes et al. GRL 2013 (application to GCM) 975 !! 976 !!---------------------------------------------------------------------- 977 REAL(wp), DIMENSION(:,:), INTENT(inout) :: Cd 978 REAL(wp), PARAMETER :: zCe = 2.23e-03_wp 979 REAL(wp), PARAMETER :: znu = 1._wp 980 REAL(wp), PARAMETER :: zmu = 1._wp 981 REAL(wp), PARAMETER :: zbeta = 1._wp 982 REAL(wp) :: zcoef 983 !!---------------------------------------------------------------------- 984 zcoef = znu + 1._wp / ( 10._wp * zbeta ) 985 986 Cd(:,:) = Cd_oce(:,:) * ( 1._wp - at_i_b(:,:) ) + & ! pure ocean drag 987 & Cd_ice * at_i_b(:,:) + & ! pure ice drag 988 & zCe * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**zmu ! change due to sea-ice morphology 989 990 END SUBROUTINE Cdn10_Lupkes2012 991 905 992 !!====================================================================== 906 993 END MODULE sbcblk_core -
branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
r6746 r6763 617 617 u_ice_b(:,:) = u_ice(:,:) 618 618 v_ice_b(:,:) = v_ice(:,:) 619 at_i_b (:,:) = SUM( a_i_b(:,:,:), dim=3 ) 619 620 620 621 END SUBROUTINE sbc_lim_bef
Note: See TracChangeset
for help on using the changeset viewer.
