- Timestamp:
- 2015-02-05T18:54:24+01:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90
r5053 r5064 116 116 CHARACTER (len=50) :: charout 117 117 REAL(wp) :: zt11, zt12, zt21, zt22, ztagnx, ztagny, delta ! 118 REAL(wp) :: za, zstms , zmask! local scalars119 REAL(wp) :: zc1, zc2, zc3 120 121 REAL(wp) :: dtevp ! time step for subcycling122 REAL(wp) :: dtotel, ecc2, ecci ! square of yield ellipse eccenticity123 REAL(wp) :: z0, zr, zcca, zccb ! temporary scalars124 REAL(wp) :: zu_ice2, zv_ice1 !125 REAL(wp) :: zddc, zdtc ! delta on corners and on centre126 REAL(wp) :: zdst ! shear at the center of the grid point127 REAL(wp) :: zdsshx, zdsshy ! term for the gradient of ocean surface128 REAL(wp) :: sigma1, sigma2 ! internal ice stress118 REAL(wp) :: za, zstms ! local scalars 119 REAL(wp) :: zc1, zc2, zc3 ! ice mass 120 121 REAL(wp) :: dtevp , z1_dtevp ! time step for subcycling 122 REAL(wp) :: dtotel, z1_dtotel, ecc2, ecci ! square of yield ellipse eccenticity 123 REAL(wp) :: z0, zr, zcca, zccb ! temporary scalars 124 REAL(wp) :: zu_ice2, zv_ice1 ! 125 REAL(wp) :: zddc, zdtc ! delta on corners and on centre 126 REAL(wp) :: zdst ! shear at the center of the grid point 127 REAL(wp) :: zdsshx, zdsshy ! term for the gradient of ocean surface 128 REAL(wp) :: sigma1, sigma2 ! internal ice stress 129 129 130 130 REAL(wp) :: zresm ! Maximal error on ice velocity … … 137 137 REAL(wp), POINTER, DIMENSION(:,:) :: zcorl1, zcorl2 ! coriolis parameter on U/V points 138 138 REAL(wp), POINTER, DIMENSION(:,:) :: za1ct , za2ct ! temporary arrays 139 REAL(wp), POINTER, DIMENSION(:,:) :: u_oce1, v_oce1! ocean u/v component on U points140 REAL(wp), POINTER, DIMENSION(:,:) :: u_oce2 , v_oce2! ocean u/v component on V points139 REAL(wp), POINTER, DIMENSION(:,:) :: v_oce1 ! ocean u/v component on U points 140 REAL(wp), POINTER, DIMENSION(:,:) :: u_oce2 ! ocean u/v component on V points 141 141 REAL(wp), POINTER, DIMENSION(:,:) :: u_ice2, v_ice1 ! ice u/v component on V/U point 142 142 REAL(wp), POINTER, DIMENSION(:,:) :: zf1 , zf2 ! arrays for internal stresses … … 156 156 157 157 CALL wrk_alloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct ) 158 CALL wrk_alloc( jpi,jpj, u_oce 1, u_oce2, u_ice2, v_oce1 , v_oce2, v_ice1 )158 CALL wrk_alloc( jpi,jpj, u_oce2, u_ice2, v_oce1 , v_ice1 ) 159 159 CALL wrk_alloc( jpi,jpj, zf1 , zu_ice, zf2 , zv_ice , zdt , zds ) 160 160 CALL wrk_alloc( jpi,jpj, zdt , zds , zs1 , zs2 , zs12 , zresr , zpice ) … … 228 228 ! zcorl2: Coriolis parameter on V-points 229 229 ! (ztagnx,ztagny): wind stress on U/V points 230 ! u_oce1: ocean u component on u points231 230 ! v_oce1: ocean v component on u points 232 231 ! u_oce2: ocean u component on v points 233 ! v_oce2: ocean v component on v points234 232 235 233 IF( nn_ice_embd == 2 ) THEN !== embedded sea ice: compute representative ice top surface ==! … … 266 264 267 265 ! Mass, coriolis coeff. and currents 268 zmass1(ji,jj) = ( zt12 *zc1 + zt11*zc2 ) / (zt11+zt12+zepsi)269 zmass2(ji,jj) = ( zt22 *zc1 + zt21*zc3 ) / (zt21+zt22+zepsi)270 zcorl1(ji,jj) = zmass1(ji,jj) * ( e1t(ji+1,jj) *fcor(ji,jj) + e1t(ji,jj)*fcor(ji+1,jj) ) &266 zmass1(ji,jj) = ( zt12 * zc1 + zt11 * zc2 ) / ( zt11 + zt12 + zepsi ) 267 zmass2(ji,jj) = ( zt22 * zc1 + zt21 * zc3 ) / ( zt21 + zt22 + zepsi ) 268 zcorl1(ji,jj) = zmass1(ji,jj) * ( e1t(ji+1,jj) * fcor(ji,jj) + e1t(ji,jj) * fcor(ji+1,jj) ) & 271 269 & / ( e1t(ji,jj) + e1t(ji+1,jj) + zepsi ) 272 zcorl2(ji,jj) = zmass2(ji,jj) * ( e2t(ji,jj+1) *fcor(ji,jj) + e2t(ji,jj)*fcor(ji,jj+1) ) &270 zcorl2(ji,jj) = zmass2(ji,jj) * ( e2t(ji,jj+1) * fcor(ji,jj) + e2t(ji,jj) * fcor(ji,jj+1) ) & 273 271 & / ( e2t(ji,jj+1) + e2t(ji,jj) + zepsi ) 274 272 ! 275 u_oce1(ji,jj) = u_oce(ji,jj)276 v_oce2(ji,jj) = v_oce(ji,jj)277 278 273 ! Ocean has no slip boundary condition 279 v_oce1(ji,jj) = 0.5 *( (v_oce(ji,jj)+v_oce(ji,jj-1))*e1t(ji,jj)&280 & +(v_oce(ji+1,jj)+v_oce(ji+1,jj-1))*e1t(ji+1,jj))&281 & /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)282 283 u_oce2(ji,jj) = 0.5 *((u_oce(ji,jj)+u_oce(ji-1,jj))*e2t(ji,jj)&284 & +(u_oce(ji,jj+1)+u_oce(ji-1,jj+1))*e2t(ji,jj+1))&285 & / (e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)274 v_oce1(ji,jj) = 0.5 * ( ( v_oce(ji ,jj) + v_oce(ji ,jj-1) ) * e1t(ji,jj) & 275 & + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji+1,jj) ) & 276 & / ( e1t(ji+1,jj) + e1t(ji,jj) ) * tmu(ji,jj) 277 278 u_oce2(ji,jj) = 0.5 * ( ( u_oce(ji,jj ) + u_oce(ji-1,jj ) ) * e2t(ji,jj) & 279 & + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj+1) ) & 280 & / ( e2t(ji,jj+1) + e2t(ji,jj) ) * tmv(ji,jj) 286 281 287 282 ! Wind stress at U,V-point … … 316 311 dtotel = dtevp / ( 2._wp * telast ) 317 312 #endif 313 z1_dtotel = 1._wp / ( 1._wp + dtotel ) 314 z1_dtevp = 1._wp / dtevp 318 315 !-ecc2: square of yield ellipse eccenticrity (reminder: must become a namelist parameter) 319 316 ecc2 = ecc * ecc … … 334 331 335 332 DO jj = k_j1+1, k_jpj-1 336 DO ji = fs_2, jpim1 !RB bug no vect opt due to tmi333 DO ji = fs_2, fs_jpim1 !RB bug no vect opt due to tmi 337 334 338 335 ! … … 357 354 divu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) & 358 355 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) & 359 & ) / area(ji,jj)356 & ) * r1_e12t(ji,jj) 360 357 361 358 zdt(ji,jj) = ( ( u_ice(ji,jj) / e2u(ji,jj) - u_ice(ji-1,jj) / e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) & 362 359 & - ( v_ice(ji,jj) / e1v(ji,jj) - v_ice(ji,jj-1) / e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) & 363 & ) / area(ji,jj)360 & ) * r1_e12t(ji,jj) 364 361 365 362 ! 366 363 zds(ji,jj) = ( ( u_ice(ji,jj+1) / e1u(ji,jj+1) - u_ice(ji,jj) / e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) & 367 364 & + ( v_ice(ji+1,jj) / e2v(ji+1,jj) - v_ice(ji,jj) / e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 368 & ) / ( e1f(ji,jj) * e2f(ji,jj)) * ( 2._wp - tmf(ji,jj) ) &365 & ) * r1_e12f(ji,jj) * ( 2._wp - tmf(ji,jj) ) & 369 366 & * tmi(ji,jj) * tmi(ji,jj+1) * tmi(ji+1,jj) * tmi(ji+1,jj+1) 370 367 371 368 372 v_ice1(ji,jj) = 0.5_wp * ( ( v_ice(ji ,jj) + v_ice(ji,jj-1) ) * e1t(ji+1,jj)&373 & + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji ,jj) ) &369 v_ice1(ji,jj) = 0.5_wp * ( ( v_ice(ji ,jj) + v_ice(ji ,jj-1) ) * e1t(ji+1,jj) & 370 & + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji ,jj) ) & 374 371 & / ( e1t(ji+1,jj) + e1t(ji,jj) ) * tmu(ji,jj) 375 372 376 u_ice2(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj) ) * e2t(ji,jj+1)&377 & + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj ) ) &373 u_ice2(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj+1) & 374 & + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj ) ) & 378 375 & / ( e2t(ji,jj+1) + e2t(ji,jj) ) * tmv(ji,jj) 379 380 END DO 381 END DO 382 CALL lbc_lnk( v_ice1, 'U', -1. ) ; CALL lbc_lnk( u_ice2, 'V', -1. ) ! lateral boundary cond. 376 END DO 377 END DO 378 CALL lbc_lnk( v_ice1 , 'U', -1. ) ; CALL lbc_lnk( u_ice2 , 'V', -1. ) ! lateral boundary cond. 383 379 384 380 DO jj = k_j1+1, k_jpj-1 … … 386 382 387 383 !- Calculate Delta at centre of grid cells 388 zdst = ( e2u(ji,jj) * v_ice1(ji,jj) - e2u(ji-1,jj ) * v_ice1(ji-1,jj ) &389 & + e1v(ji,jj) * u_ice2(ji,jj) - e1v(ji ,jj-1) * u_ice2(ji ,jj-1) &390 & ) / area(ji,jj)391 392 delta = SQRT( divu_i(ji,jj) * divu_i(ji,jj) + ( zdt(ji,jj) * zdt(ji,jj) + zdst * zdst) * usecc2 )384 zdst = ( e2u(ji,jj) * v_ice1(ji,jj) - e2u(ji-1,jj ) * v_ice1(ji-1,jj ) & 385 & + e1v(ji,jj) * u_ice2(ji,jj) - e1v(ji ,jj-1) * u_ice2(ji ,jj-1) & 386 & ) * r1_e12t(ji,jj) 387 388 delta = SQRT( divu_i(ji,jj)**2 + ( zdt(ji,jj)**2 + zdst**2 ) * usecc2 ) 393 389 delta_i(ji,jj) = delta + creepl 394 !-Calculate stress tensor components zs1 and zs2 395 !-at centre of grid cells (see section 3.5 of CICE user's guide). 396 zs1(ji,jj) = ( zs1(ji,jj) + dtotel * ( ( divu_i(ji,jj) / delta_i(ji,jj) - delta / delta_i(ji,jj) ) & 397 & * zpresh(ji,jj) ) ) / ( 1._wp + dtotel ) 398 zs2(ji,jj) = ( zs2(ji,jj) + dtotel * ( ecci * zdt(ji,jj) / delta_i(ji,jj) * zpresh(ji,jj) ) ) & 399 & / ( 1._wp + dtotel ) 400 401 END DO 402 END DO 403 404 CALL lbc_lnk( zs1(:,:), 'T', 1. ) 405 CALL lbc_lnk( zs2(:,:), 'T', 1. ) 406 407 DO jj = k_j1+1, k_jpj-1 408 DO ji = fs_2, fs_jpim1 390 409 391 !- Calculate Delta on corners 410 zddc = ( ( v_ice1(ji,jj+1)/e1u(ji,jj+1) & 411 & -v_ice1(ji,jj)/e1u(ji,jj) & 412 & )*e1f(ji,jj)*e1f(ji,jj) & 413 & +( u_ice2(ji+1,jj)/e2v(ji+1,jj) & 414 & -u_ice2(ji,jj)/e2v(ji,jj) & 415 & )*e2f(ji,jj)*e2f(ji,jj) & 416 & ) & 417 & / ( e1f(ji,jj) * e2f(ji,jj) ) 418 419 zdtc = (-( v_ice1(ji,jj+1)/e1u(ji,jj+1) & 420 & -v_ice1(ji,jj)/e1u(ji,jj) & 421 & )*e1f(ji,jj)*e1f(ji,jj) & 422 & +( u_ice2(ji+1,jj)/e2v(ji+1,jj) & 423 & -u_ice2(ji,jj)/e2v(ji,jj) & 424 & )*e2f(ji,jj)*e2f(ji,jj) & 425 & ) & 426 & / ( e1f(ji,jj) * e2f(ji,jj) ) 427 428 zddc = SQRT(zddc**2+(zdtc**2+zds(ji,jj)**2)*usecc2) + creepl 429 430 !-Calculate stress tensor component zs12 at corners (see section 3.5 of CICE user's guide). 431 zs12(ji,jj) = ( zs12(ji,jj) + dtotel * & 432 & ( ecci * zds(ji,jj) / ( 2._wp * zddc ) * zpreshc(ji,jj) ) ) & 433 & / ( 1.0 + dtotel ) 434 435 END DO ! ji 436 END DO ! jj 437 438 CALL lbc_lnk( zs12(:,:), 'F', 1. ) 392 zddc = ( ( v_ice1(ji,jj+1) / e1u(ji,jj+1) - v_ice1(ji,jj) / e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) & 393 & + ( u_ice2(ji+1,jj) / e2v(ji+1,jj) - u_ice2(ji,jj) / e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 394 & ) * r1_e12f(ji,jj) 395 396 zdtc = (- ( v_ice1(ji,jj+1) / e1u(ji,jj+1) - v_ice1(ji,jj) / e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) & 397 & + ( u_ice2(ji+1,jj) / e2v(ji+1,jj) - u_ice2(ji,jj) / e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 398 & ) * r1_e12f(ji,jj) 399 400 zddc = SQRT( zddc**2 + ( zdtc**2 + zds(ji,jj)**2 ) * usecc2 ) + creepl 401 402 !-Calculate stress tensor components zs1 and zs2 at centre of grid cells (see section 3.5 of CICE user's guide). 403 zs1(ji,jj) = ( zs1 (ji,jj) + dtotel * ( divu_i(ji,jj) - delta ) / delta_i(ji,jj) * zpresh(ji,jj) & 404 & ) * z1_dtotel 405 zs2(ji,jj) = ( zs2 (ji,jj) + dtotel * ecci * zdt(ji,jj) / delta_i(ji,jj) * zpresh(ji,jj) & 406 & ) * z1_dtotel 407 !-Calculate stress tensor component zs12 at corners 408 zs12(ji,jj) = ( zs12(ji,jj) + dtotel * ecci * zds(ji,jj) / ( 2._wp * zddc ) * zpreshc(ji,jj) & 409 & ) * z1_dtotel 410 411 END DO 412 END DO 413 CALL lbc_lnk( zs1 , 'T', 1. ) ; CALL lbc_lnk( zs2, 'T', 1. ) 414 CALL lbc_lnk( zs12, 'F', 1. ) 439 415 440 416 ! Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) … … 442 418 DO ji = fs_2, fs_jpim1 443 419 !- contribution of zs1, zs2 and zs12 to zf1 444 zf1(ji,jj) = 0.5 *( (zs1(ji+1,jj)-zs1(ji,jj))*e2u(ji,jj)&445 & +(zs2(ji+1,jj)*e2t(ji+1,jj)**2-zs2(ji,jj)*e2t(ji,jj)**2)/e2u(ji,jj)&446 & +2.0*(zs12(ji,jj)*e1f(ji,jj)**2-zs12(ji,jj-1)*e1f(ji,jj-1)**2)/e1u(ji,jj)&447 & ) / ( e1u(ji,jj)*e2u(ji,jj))420 zf1(ji,jj) = 0.5 * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj) & 421 & + ( zs2(ji+1,jj) * e2t(ji+1,jj)**2 - zs2(ji,jj) * e2t(ji,jj)**2 ) / e2u(ji,jj) & 422 & + 2.0 * ( zs12(ji,jj) * e1f(ji,jj)**2 - zs12(ji,jj-1) * e1f(ji,jj-1)**2 ) / e1u(ji,jj) & 423 & ) * r1_e12u(ji,jj) 448 424 ! contribution of zs1, zs2 and zs12 to zf2 449 zf2(ji,jj) = 0.5*( (zs1(ji,jj+1)-zs1(ji,jj))*e1v(ji,jj) & 450 & -(zs2(ji,jj+1)*e1t(ji,jj+1)**2 - zs2(ji,jj)*e1t(ji,jj)**2)/e1v(ji,jj) & 451 & + 2.0*(zs12(ji,jj)*e2f(ji,jj)**2 - & 452 zs12(ji-1,jj)*e2f(ji-1,jj)**2)/e2v(ji,jj) & 453 & ) / ( e1v(ji,jj)*e2v(ji,jj) ) 425 zf2(ji,jj) = 0.5 * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj) & 426 & - ( zs2(ji,jj+1) * e1t(ji,jj+1)**2 - zs2(ji,jj) * e1t(ji,jj)**2 ) / e1v(ji,jj) & 427 & + 2.0 * ( zs12(ji,jj) * e2f(ji,jj)**2 - zs12(ji-1,jj) * e2f(ji-1,jj)**2 ) / e2v(ji,jj) & 428 & ) * r1_e12v(ji,jj) 454 429 END DO 455 430 END DO … … 463 438 DO jj = k_j1+1, k_jpj-1 464 439 DO ji = fs_2, fs_jpim1 465 zmask = (1.0-MAX(0._wp,SIGN(1._wp,-zmass1(ji,jj))))*tmu(ji,jj)466 z0 = zmass1(ji,jj) /dtevp440 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) ) ) * tmu(ji,jj) 441 z0 = zmass1(ji,jj) * z1_dtevp 467 442 468 443 ! SB modif because ocean has no slip boundary condition 469 zv_ice1 = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji,jj) & 470 & +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji+1,jj)) & 471 & /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) 472 za = rhoco*SQRT((u_ice(ji,jj)-u_oce1(ji,jj))**2 + & 473 (zv_ice1-v_oce1(ji,jj))**2) * (1.0-zfrld1(ji,jj)) 474 zr = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + & 475 za*(u_oce1(ji,jj)) 476 zcca = z0+za 444 zv_ice1 = 0.5 * ( ( v_ice(ji ,jj) + v_ice(ji ,jj-1) ) * e1t(ji ,jj) & 445 & + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji+1,jj) ) & 446 & / ( e1t(ji+1,jj) + e1t(ji,jj) ) * tmu(ji,jj) 447 za = rhoco * SQRT( ( u_ice(ji,jj) - u_oce(ji,jj) )**2 + & 448 & ( zv_ice1 - v_oce1(ji,jj) )**2 ) * ( 1.0 - zfrld1(ji,jj) ) 449 zr = z0 * u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + za * u_oce(ji,jj) 450 zcca = z0 + za 477 451 zccb = zcorl1(ji,jj) 478 u_ice(ji,jj) = (zr+zccb*zv_ice1)/(zcca+zepsi)*zmask 479 452 u_ice(ji,jj) = ( zr + zccb * zv_ice1 ) / ( zcca + zepsi ) * rswitch 480 453 END DO 481 454 END DO … … 492 465 DO ji = fs_2, fs_jpim1 493 466 494 zmask = (1.0-MAX(0._wp,SIGN(1._wp,-zmass2(ji,jj))))*tmv(ji,jj)495 z0 = zmass2(ji,jj) /dtevp467 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) ) * tmv(ji,jj) 468 z0 = zmass2(ji,jj) * z1_dtevp 496 469 ! SB modif because ocean has no slip boundary condition 497 zu_ice2 = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj) & 498 & + (u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj+1)) & 499 & /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj) 500 za = rhoco*SQRT((zu_ice2-u_oce2(ji,jj))**2 + & 501 (v_ice(ji,jj)-v_oce2(ji,jj))**2)*(1.0-zfrld2(ji,jj)) 502 zr = z0*v_ice(ji,jj) + zf2(ji,jj) + & 503 za2ct(ji,jj) + za*(v_oce2(ji,jj)) 504 zcca = z0+za 470 zu_ice2 = 0.5 * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj) & 471 & + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj+1) ) & 472 & / ( e2t(ji,jj+1) + e2t(ji,jj) ) * tmv(ji,jj) 473 za = rhoco * SQRT( ( zu_ice2 - u_oce2(ji,jj) )**2 + & 474 & ( v_ice(ji,jj) - v_oce(ji,jj))**2 ) * ( 1.0 - zfrld2(ji,jj) ) 475 zr = z0 * v_ice(ji,jj) + zf2(ji,jj) + za2ct(ji,jj) + za * v_oce(ji,jj) 476 zcca = z0 + za 505 477 zccb = zcorl2(ji,jj) 506 v_ice(ji,jj) = (zr-zccb*zu_ice2)/(zcca+zepsi)*zmask 507 478 v_ice(ji,jj) = ( zr - zccb * zu_ice2 ) / ( zcca + zepsi ) * rswitch 508 479 END DO 509 480 END DO … … 520 491 DO jj = k_j1+1, k_jpj-1 521 492 DO ji = fs_2, fs_jpim1 522 zmask = (1.0-MAX(0._wp,SIGN(1._wp,-zmass2(ji,jj))))*tmv(ji,jj)523 z0 = zmass2(ji,jj) /dtevp493 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) ) * tmv(ji,jj) 494 z0 = zmass2(ji,jj) * z1_dtevp 524 495 ! SB modif because ocean has no slip boundary condition 525 zu_ice2 = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj) & 526 & +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj+1)) & 527 & /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj) 528 529 za = rhoco*SQRT((zu_ice2-u_oce2(ji,jj))**2 + & 530 (v_ice(ji,jj)-v_oce2(ji,jj))**2)*(1.0-zfrld2(ji,jj)) 531 zr = z0*v_ice(ji,jj) + zf2(ji,jj) + & 532 za2ct(ji,jj) + za*(v_oce2(ji,jj)) 533 zcca = z0+za 496 zu_ice2 = 0.5 * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj) & 497 & +( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj+1) ) & 498 & / ( e2t(ji,jj+1) + e2t(ji,jj) ) * tmv(ji,jj) 499 500 za = rhoco * SQRT( ( zu_ice2 - u_oce2(ji,jj) )**2 + & 501 & ( v_ice(ji,jj) - v_oce(ji,jj) )**2 ) * ( 1.0 - zfrld2(ji,jj) ) 502 zr = z0 * v_ice(ji,jj) + zf2(ji,jj) + za2ct(ji,jj) + za * v_oce(ji,jj) 503 zcca = z0 + za 534 504 zccb = zcorl2(ji,jj) 535 v_ice(ji,jj) = (zr-zccb*zu_ice2)/(zcca+zepsi)*zmask 536 505 v_ice(ji,jj) = ( zr - zccb * zu_ice2 ) / ( zcca + zepsi ) * rswitch 537 506 END DO 538 507 END DO … … 548 517 DO jj = k_j1+1, k_jpj-1 549 518 DO ji = fs_2, fs_jpim1 550 zmask = (1.0-MAX(0._wp,SIGN(1._wp,-zmass1(ji,jj))))*tmu(ji,jj) 551 z0 = zmass1(ji,jj)/dtevp 552 zv_ice1 = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji,jj) & 553 & +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji+1,jj)) & 554 & /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) 555 556 za = rhoco*SQRT((u_ice(ji,jj)-u_oce1(ji,jj))**2 + & 557 (zv_ice1-v_oce1(ji,jj))**2)*(1.0-zfrld1(ji,jj)) 558 zr = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + & 559 za*(u_oce1(ji,jj)) 560 zcca = z0+za 519 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) ) ) * tmu(ji,jj) 520 z0 = zmass1(ji,jj) * z1_dtevp 521 zv_ice1 = 0.5 * ( ( v_ice(ji ,jj) + v_ice(ji ,jj-1) ) * e1t(ji,jj) & 522 & + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji+1,jj) ) & 523 & / ( e1t(ji+1,jj) + e1t(ji,jj) ) * tmu(ji,jj) 524 525 za = rhoco * SQRT( ( u_ice(ji,jj) - u_oce(ji,jj) )**2 + & 526 & ( zv_ice1 - v_oce1(ji,jj) )**2 ) * ( 1.0 - zfrld1(ji,jj) ) 527 zr = z0 * u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + za * u_oce(ji,jj) 528 zcca = z0 + za 561 529 zccb = zcorl1(ji,jj) 562 u_ice(ji,jj) = ( zr+zccb*zv_ice1)/(zcca+zepsi)*zmask563 END DO ! ji564 END DO ! jj530 u_ice(ji,jj) = ( zr + zccb * zv_ice1 ) / ( zcca + zepsi ) * rswitch 531 END DO 532 END DO 565 533 566 534 CALL lbc_lnk( u_ice(:,:), 'U', -1. ) … … 577 545 !--- Convergence test. 578 546 DO jj = k_j1+1 , k_jpj-1 579 zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ) , & 580 ABS( v_ice(:,jj) - zv_ice(:,jj) ) ) 581 END DO 582 zresm = MAXVAL( zresr( 1:jpi , k_j1+1:k_jpj-1 ) ) 547 zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) ) 548 END DO 549 zresm = MAXVAL( zresr( 1:jpi, k_j1+1:k_jpj-1 ) ) 583 550 IF( lk_mpp ) CALL mpp_max( zresm ) ! max over the global domain 584 551 ENDIF … … 637 604 IF ( vt_i(ji,jj) <= zvmin ) THEN 638 605 639 divu_i(ji,jj) = ( e2u(ji,jj)*u_ice(ji,jj) & 640 & -e2u(ji-1,jj)*u_ice(ji-1,jj) & 641 & +e1v(ji,jj)*v_ice(ji,jj) & 642 & -e1v(ji,jj-1)*v_ice(ji,jj-1) & 643 & ) & 644 & / area(ji,jj) 645 646 zdt(ji,jj) = ( ( u_ice(ji,jj)/e2u(ji,jj) & 647 & -u_ice(ji-1,jj)/e2u(ji-1,jj) & 648 & )*e2t(ji,jj)*e2t(ji,jj) & 649 & -( v_ice(ji,jj)/e1v(ji,jj) & 650 & -v_ice(ji,jj-1)/e1v(ji,jj-1) & 651 & )*e1t(ji,jj)*e1t(ji,jj) & 652 & ) & 653 & / area(ji,jj) 606 divu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj ) * u_ice(ji-1,jj ) & 607 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji ,jj-1) * v_ice(ji ,jj-1) & 608 & ) * r1_e12t(ji,jj) 609 610 zdt(ji,jj) = ( ( u_ice(ji,jj) / e2u(ji,jj) - u_ice(ji-1,jj) / e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) & 611 & -( v_ice(ji,jj) / e1v(ji,jj) - v_ice(ji,jj-1) / e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) & 612 & ) * r1_e12t(ji,jj) 654 613 ! 655 614 ! SB modif because ocean has no slip boundary condition 656 zds(ji,jj) = ( ( u_ice(ji,jj+1) / e1u(ji,jj+1) & 657 & - u_ice(ji,jj) / e1u(ji,jj) ) & 658 & * e1f(ji,jj) * e1f(ji,jj) & 659 & + ( v_ice(ji+1,jj) / e2v(ji+1,jj) & 660 & - v_ice(ji,jj) / e2v(ji,jj) ) & 661 & * e2f(ji,jj) * e2f(ji,jj) ) & 662 & / ( e1f(ji,jj) * e2f(ji,jj) ) * ( 2.0 - tmf(ji,jj) ) & 663 & * tmi(ji,jj) * tmi(ji,jj+1) & 664 & * tmi(ji+1,jj) * tmi(ji+1,jj+1) 665 666 zdst = ( e2u( ji , jj ) * v_ice1(ji ,jj ) & 667 & - e2u( ji-1, jj ) * v_ice1(ji-1,jj ) & 668 & + e1v( ji , jj ) * u_ice2(ji ,jj ) & 669 & - e1v( ji , jj-1 ) * u_ice2(ji ,jj-1) ) / area(ji,jj) 670 671 delta = SQRT( divu_i(ji,jj)*divu_i(ji,jj) + ( zdt(ji,jj)*zdt(ji,jj) + zdst*zdst ) * usecc2 ) 615 zds(ji,jj) = ( ( u_ice(ji,jj+1) / e1u(ji,jj+1) - u_ice(ji,jj) / e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) & 616 & +( v_ice(ji+1,jj) / e2v(ji+1,jj) - v_ice(ji,jj) / e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 617 & ) * r1_e12f(ji,jj) * ( 2.0 - tmf(ji,jj) ) & 618 & * tmi(ji,jj) * tmi(ji,jj+1) * tmi(ji+1,jj) * tmi(ji+1,jj+1) 619 620 zdst = ( e2u(ji,jj) * v_ice1(ji,jj) - e2u(ji-1,jj ) * v_ice1(ji-1,jj ) & 621 & + e1v(ji,jj) * u_ice2(ji,jj) - e1v(ji ,jj-1) * u_ice2(ji ,jj-1) ) * r1_e12t(ji,jj) 622 623 delta = SQRT( divu_i(ji,jj)**2 + ( zdt(ji,jj)**2 + zdst**2 ) * usecc2 ) 672 624 delta_i(ji,jj) = delta + creepl 673 625 674 626 ENDIF 675 676 END DO !jj 677 END DO !ji 627 END DO 628 END DO 678 629 ! 679 630 !------------------------------------------------------------------------------! … … 684 635 DO jj = k_j1+1, k_jpj-1 685 636 DO ji = fs_2, fs_jpim1 686 zdst= ( e2u( ji , jj ) * v_ice1(ji,jj) & 687 & - e2u( ji-1, jj ) * v_ice1(ji-1,jj) & 688 & + e1v( ji , jj ) * u_ice2(ji,jj) & 689 & - e1v( ji , jj-1 ) * u_ice2(ji,jj-1) ) / area(ji,jj) 637 zdst = ( e2u(ji,jj) * v_ice1(ji,jj) - e2u( ji-1, jj ) * v_ice1(ji-1,jj) & 638 & + e1v(ji,jj) * u_ice2(ji,jj) - e1v( ji , jj-1 ) * u_ice2(ji,jj-1) ) * r1_e12t(ji,jj) 690 639 shear_i(ji,jj) = SQRT( zdt(ji,jj) * zdt(ji,jj) + zdst * zdst ) 691 640 END DO … … 741 690 ! 742 691 CALL wrk_dealloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct ) 743 CALL wrk_dealloc( jpi,jpj, u_oce 1, u_oce2, u_ice2, v_oce1 , v_oce2, v_ice1 )692 CALL wrk_dealloc( jpi,jpj, u_oce2, u_ice2, v_oce1 , v_ice1 ) 744 693 CALL wrk_dealloc( jpi,jpj, zf1 , zu_ice, zf2 , zv_ice , zdt , zds ) 745 694 CALL wrk_dealloc( jpi,jpj, zdt , zds , zs1 , zs2 , zs12 , zresr , zpice )
Note: See TracChangeset
for help on using the changeset viewer.
