- Timestamp:
- 2015-02-06T19:12:57+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
r5064 r5067 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 143 REAL(wp), POINTER, DIMENSION(:,:) :: zmask ! mask ocean grid points 143 144 144 145 REAL(wp), POINTER, DIMENSION(:,:) :: zdt ! tension at centre of grid cells … … 156 157 157 158 CALL wrk_alloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct ) 158 CALL wrk_alloc( jpi,jpj, u_oce2, u_ice2, v_oce1 , v_ice1 )159 CALL wrk_alloc( jpi,jpj, u_oce2, u_ice2, v_oce1 , v_ice1 , zmask ) 159 160 CALL wrk_alloc( jpi,jpj, zf1 , zu_ice, zf2 , zv_ice , zdt , zds ) 160 161 CALL wrk_alloc( jpi,jpj, zdt , zds , zs1 , zs2 , zs12 , zresr , zpice ) … … 187 188 188 189 #if defined key_lim3 189 CALL lim_itd_me_icestrength( ridge_scheme_swi) ! LIM-3: Ice strength on T-points190 CALL lim_itd_me_icestrength( nn_icestr ) ! LIM-3: Ice strength on T-points 190 191 #endif 191 192 … … 193 194 DO ji = 1 , jpi 194 195 #if defined key_lim3 195 zpresh(ji,jj) = tm s(ji,jj) * strength(ji,jj)196 zpresh(ji,jj) = tmask(ji,jj,1) * strength(ji,jj) 196 197 #endif 197 198 #if defined key_lim2 198 zpresh(ji,jj) = tm s(ji,jj) * pstar * vt_i(ji,jj) * EXP( -c_rhg * (1. - at_i(ji,jj) ) )199 #endif 200 ! tmi= 1 where there is ice or on land201 tmi(ji,jj) = 1._wp - ( 1._wp - MAX( 0._wp , SIGN ( 1._wp , vt_i(ji,jj) - zepsi ) ) ) * tms(ji,jj)199 zpresh(ji,jj) = tmask(ji,jj,1) * pstar * vt_i(ji,jj) * EXP( -c_rhg * (1. - at_i(ji,jj) ) ) 200 #endif 201 ! zmask = 1 where there is ice or on land 202 zmask(ji,jj) = 1._wp - ( 1._wp - MAX( 0._wp , SIGN ( 1._wp , vt_i(ji,jj) - zepsi ) ) ) * tmask(ji,jj,1) 202 203 END DO 203 204 END DO … … 207 208 DO jj = k_j1+1, k_jpj-1 208 209 DO ji = 2, jpim1 !RB caution no fs_ (ji+1,jj+1) 209 zstms = tm s(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + tms(ji,jj+1) * wght(ji+1,jj+1,1,2) + &210 & tm s(ji+1,jj) * wght(ji+1,jj+1,2,1) + tms(ji,jj) * wght(ji+1,jj+1,1,1)210 zstms = tmask(ji+1,jj+1,1) * wght(ji+1,jj+1,2,2) + tmask(ji,jj+1,1) * wght(ji+1,jj+1,1,2) + & 211 & tmask(ji+1,jj,1) * wght(ji+1,jj+1,2,1) + tmask(ji,jj,1) * wght(ji+1,jj+1,1,1) 211 212 zpreshc(ji,jj) = ( zpresh(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + zpresh(ji,jj+1) * wght(ji+1,jj+1,1,2) + & 212 213 & zpresh(ji+1,jj) * wght(ji+1,jj+1,2,1) + zpresh(ji,jj) * wght(ji+1,jj+1,1,1) & … … 250 251 DO ji = fs_2, fs_jpim1 251 252 252 zc1 = tm s(ji ,jj) * ( rhosn * vt_s(ji ,jj ) + rhoic * vt_i(ji ,jj ) )253 zc2 = tm s(ji+1,jj) * ( rhosn * vt_s(ji+1,jj ) + rhoic * vt_i(ji+1,jj ) )254 zc3 = tm s(ji ,jj+1) * ( rhosn * vt_s(ji ,jj+1) + rhoic * vt_i(ji ,jj+1) )255 256 zt11 = tm s(ji ,jj) * e1t(ji ,jj)257 zt12 = tm s(ji+1,jj) * e1t(ji+1,jj)258 zt21 = tm s(ji,jj) * e2t(ji,jj )259 zt22 = tm s(ji,jj+1) * e2t(ji,jj+1)253 zc1 = tmask(ji ,jj ,1) * ( rhosn * vt_s(ji ,jj ) + rhoic * vt_i(ji ,jj ) ) 254 zc2 = tmask(ji+1,jj ,1) * ( rhosn * vt_s(ji+1,jj ) + rhoic * vt_i(ji+1,jj ) ) 255 zc3 = tmask(ji ,jj+1,1) * ( rhosn * vt_s(ji ,jj+1) + rhoic * vt_i(ji ,jj+1) ) 256 257 zt11 = tmask(ji ,jj,1) * e1t(ji ,jj) 258 zt12 = tmask(ji+1,jj,1) * e1t(ji+1,jj) 259 zt21 = tmask(ji,jj ,1) * e2t(ji,jj ) 260 zt22 = tmask(ji,jj+1,1) * e2t(ji,jj+1) 260 261 261 262 ! Leads area. … … 274 275 v_oce1(ji,jj) = 0.5 * ( ( v_oce(ji ,jj) + v_oce(ji ,jj-1) ) * e1t(ji,jj) & 275 276 & + ( 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 & / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 277 278 278 279 u_oce2(ji,jj) = 0.5 * ( ( u_oce(ji,jj ) + u_oce(ji-1,jj ) ) * e2t(ji,jj) & 279 280 & + ( 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)281 & / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1) 281 282 282 283 ! Wind stress at U,V-point … … 305 306 ! 306 307 ! Time step for subcycling 307 dtevp = rdt_ice / n evp308 dtevp = rdt_ice / nn_nevp 308 309 #if defined key_lim3 309 dtotel = dtevp / ( 2._wp * r elast * rdt_ice )310 dtotel = dtevp / ( 2._wp * rn_relast * rdt_ice ) 310 311 #else 311 312 dtotel = dtevp / ( 2._wp * telast ) … … 314 315 z1_dtevp = 1._wp / dtevp 315 316 !-ecc2: square of yield ellipse eccenticrity (reminder: must become a namelist parameter) 316 ecc2 = ecc *ecc317 ecc2 = rn_ecc * rn_ecc 317 318 ecci = 1. / ecc2 318 319 … … 323 324 324 325 ! !----------------------! 325 DO jter = 1 , n evp! loop over jter !326 DO jter = 1 , nn_nevp ! loop over jter ! 326 327 ! !----------------------! 327 328 DO jj = k_j1, k_jpj-1 … … 331 332 332 333 DO jj = k_j1+1, k_jpj-1 333 DO ji = fs_2, fs_jpim1 !RB bug no vect opt due to tmi334 DO ji = fs_2, fs_jpim1 !RB bug no vect opt due to zmask 334 335 335 336 ! … … 363 364 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) & 364 365 & + ( v_ice(ji+1,jj) / e2v(ji+1,jj) - v_ice(ji,jj) / e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 365 & ) * r1_e12f(ji,jj) * ( 2._wp - tmf(ji,jj) ) &366 & * tmi(ji,jj) * tmi(ji,jj+1) * tmi(ji+1,jj) * tmi(ji+1,jj+1)366 & ) * r1_e12f(ji,jj) * ( 2._wp - fmask(ji,jj,1) ) & 367 & * zmask(ji,jj) * zmask(ji,jj+1) * zmask(ji+1,jj) * zmask(ji+1,jj+1) 367 368 368 369 369 370 v_ice1(ji,jj) = 0.5_wp * ( ( v_ice(ji ,jj) + v_ice(ji ,jj-1) ) * e1t(ji+1,jj) & 370 371 & + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji ,jj) ) & 371 & / ( e1t(ji+1,jj) + e1t(ji,jj) ) * tmu(ji,jj)372 & / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 372 373 373 374 u_ice2(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj+1) & 374 375 & + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj ) ) & 375 & / ( e2t(ji,jj+1) + e2t(ji,jj) ) * tmv(ji,jj)376 & / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1) 376 377 END DO 377 378 END DO … … 387 388 388 389 delta = SQRT( divu_i(ji,jj)**2 + ( zdt(ji,jj)**2 + zdst**2 ) * usecc2 ) 389 delta_i(ji,jj) = delta + creepl390 delta_i(ji,jj) = delta + rn_creepl 390 391 391 392 !- Calculate Delta on corners … … 398 399 & ) * r1_e12f(ji,jj) 399 400 400 zddc = SQRT( zddc**2 + ( zdtc**2 + zds(ji,jj)**2 ) * usecc2 ) + creepl401 zddc = SQRT( zddc**2 + ( zdtc**2 + zds(ji,jj)**2 ) * usecc2 ) + rn_creepl 401 402 402 403 !-Calculate stress tensor components zs1 and zs2 at centre of grid cells (see section 3.5 of CICE user's guide). … … 438 439 DO jj = k_j1+1, k_jpj-1 439 440 DO ji = fs_2, fs_jpim1 440 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) ) ) * tmu(ji,jj)441 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) ) ) * umask(ji,jj,1) 441 442 z0 = zmass1(ji,jj) * z1_dtevp 442 443 … … 444 445 zv_ice1 = 0.5 * ( ( v_ice(ji ,jj) + v_ice(ji ,jj-1) ) * e1t(ji ,jj) & 445 446 & + ( 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 & / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 447 448 za = rhoco * SQRT( ( u_ice(ji,jj) - u_oce(ji,jj) )**2 + & 448 449 & ( zv_ice1 - v_oce1(ji,jj) )**2 ) * ( 1.0 - zfrld1(ji,jj) ) … … 456 457 CALL lbc_lnk( u_ice(:,:), 'U', -1. ) 457 458 #if defined key_agrif && defined key_lim2 458 CALL agrif_rhg_lim2( jter, n evp, 'U' )459 CALL agrif_rhg_lim2( jter, nn_nevp, 'U' ) 459 460 #endif 460 461 #if defined key_bdy … … 465 466 DO ji = fs_2, fs_jpim1 466 467 467 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) ) * tmv(ji,jj)468 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) ) * vmask(ji,jj,1) 468 469 z0 = zmass2(ji,jj) * z1_dtevp 469 470 ! SB modif because ocean has no slip boundary condition 470 471 zu_ice2 = 0.5 * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj) & 471 472 & + ( 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 & / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1) 473 474 za = rhoco * SQRT( ( zu_ice2 - u_oce2(ji,jj) )**2 + & 474 475 & ( v_ice(ji,jj) - v_oce(ji,jj))**2 ) * ( 1.0 - zfrld2(ji,jj) ) … … 482 483 CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 483 484 #if defined key_agrif && defined key_lim2 484 CALL agrif_rhg_lim2( jter, n evp, 'V' )485 CALL agrif_rhg_lim2( jter, nn_nevp, 'V' ) 485 486 #endif 486 487 #if defined key_bdy … … 491 492 DO jj = k_j1+1, k_jpj-1 492 493 DO ji = fs_2, fs_jpim1 493 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) ) * tmv(ji,jj)494 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) ) * vmask(ji,jj,1) 494 495 z0 = zmass2(ji,jj) * z1_dtevp 495 496 ! SB modif because ocean has no slip boundary condition 496 497 zu_ice2 = 0.5 * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj) & 497 498 & +( 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 & / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1) 499 500 500 501 za = rhoco * SQRT( ( zu_ice2 - u_oce2(ji,jj) )**2 + & … … 509 510 CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 510 511 #if defined key_agrif && defined key_lim2 511 CALL agrif_rhg_lim2( jter, n evp, 'V' )512 CALL agrif_rhg_lim2( jter, nn_nevp, 'V' ) 512 513 #endif 513 514 #if defined key_bdy … … 517 518 DO jj = k_j1+1, k_jpj-1 518 519 DO ji = fs_2, fs_jpim1 519 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) ) ) * tmu(ji,jj)520 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) ) ) * umask(ji,jj,1) 520 521 z0 = zmass1(ji,jj) * z1_dtevp 521 522 zv_ice1 = 0.5 * ( ( v_ice(ji ,jj) + v_ice(ji ,jj-1) ) * e1t(ji,jj) & 522 523 & + ( 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 & / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 524 525 525 526 za = rhoco * SQRT( ( u_ice(ji,jj) - u_oce(ji,jj) )**2 + & … … 534 535 CALL lbc_lnk( u_ice(:,:), 'U', -1. ) 535 536 #if defined key_agrif && defined key_lim2 536 CALL agrif_rhg_lim2( jter, n evp, 'U' )537 CALL agrif_rhg_lim2( jter, nn_nevp, 'U' ) 537 538 #endif 538 539 #if defined key_bdy … … 572 573 CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 573 574 #if defined key_agrif && defined key_lim2 574 CALL agrif_rhg_lim2( n evp ,nevp, 'U' )575 CALL agrif_rhg_lim2( n evp ,nevp, 'V' )575 CALL agrif_rhg_lim2( nn_nevp , nn_nevp, 'U' ) 576 CALL agrif_rhg_lim2( nn_nevp , nn_nevp, 'V' ) 576 577 #endif 577 578 #if defined key_bdy … … 585 586 v_ice1(ji,jj) = 0.5_wp * ( ( v_ice(ji ,jj) + v_ice(ji, jj-1) ) * e1t(ji+1,jj) & 586 587 & + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji ,jj) ) & 587 & / ( e1t(ji+1,jj) + e1t(ji,jj) ) * tmu(ji,jj)588 & / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 588 589 589 590 u_ice2(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj+1) & 590 591 & + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj ) ) & 591 & / ( e2t(ji,jj+1) + e2t(ji,jj) ) * tmv(ji,jj)592 & / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1) 592 593 ENDIF 593 594 END DO … … 599 600 ! Recompute delta, shear and div, inputs for mechanical redistribution 600 601 DO jj = k_j1+1, k_jpj-1 601 DO ji = fs_2, jpim1 !RB bug no vect opt due to tmi602 DO ji = fs_2, jpim1 !RB bug no vect opt due to zmask 602 603 !- divu_i(:,:), zdt(:,:): divergence and tension at centre 603 604 !- zds(:,:): shear on northeast corner of grid cells … … 615 616 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 617 & +( 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)618 & ) * r1_e12f(ji,jj) * ( 2.0 - fmask(ji,jj,1) ) & 619 & * zmask(ji,jj) * zmask(ji,jj+1) * zmask(ji+1,jj) * zmask(ji+1,jj+1) 619 620 620 621 zdst = ( e2u(ji,jj) * v_ice1(ji,jj) - e2u(ji-1,jj ) * v_ice1(ji-1,jj ) & … … 622 623 623 624 delta = SQRT( divu_i(ji,jj)**2 + ( zdt(ji,jj)**2 + zdst**2 ) * usecc2 ) 624 delta_i(ji,jj) = delta + creepl625 delta_i(ji,jj) = delta + rn_creepl 625 626 626 627 ENDIF … … 690 691 ! 691 692 CALL wrk_dealloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct ) 692 CALL wrk_dealloc( jpi,jpj, u_oce2, u_ice2, v_oce1 , v_ice1 )693 CALL wrk_dealloc( jpi,jpj, u_oce2, u_ice2, v_oce1 , v_ice1 , zmask ) 693 694 CALL wrk_dealloc( jpi,jpj, zf1 , zu_ice, zf2 , zv_ice , zdt , zds ) 694 695 CALL wrk_dealloc( jpi,jpj, zdt , zds , zs1 , zs2 , zs12 , zresr , zpice )
Note: See TracChangeset
for help on using the changeset viewer.
