Changeset 5600 for branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90
- Timestamp:
- 2015-07-15T17:46:12+02:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90
r5034 r5600 102 102 !! and charge ellipse. 103 103 !! The user should make sure that the parameters 104 !! n evp, telast andcreepl maintain stress state104 !! nn_nevp, elastic time scale and rn_creepl maintain stress state 105 105 !! on the charge ellipse for plastic flow 106 106 !! e.g. in the Canadian Archipelago … … 108 108 !! References : Hunke and Dukowicz, JPO97 109 109 !! Bouillon et al., Ocean Modelling 2009 110 !! Vancoppenolle et al., Ocean Modelling 2008111 110 !!------------------------------------------------------------------- 112 111 INTEGER, INTENT(in) :: k_j1 ! southern j-index for ice computation … … 117 116 CHARACTER (len=50) :: charout 118 117 REAL(wp) :: zt11, zt12, zt21, zt22, ztagnx, ztagny, delta ! 119 REAL(wp) :: za, zstms , zmask! local scalars120 REAL(wp) :: zc1, zc2, zc3 121 122 REAL(wp) :: dtevp ! time step for subcycling123 REAL(wp) :: dtotel, ecc2, ecci ! square of yield ellipse eccenticity124 REAL(wp) :: z0, zr, zcca, zccb ! temporary scalars125 REAL(wp) :: zu_ice2, zv_ice1 !126 REAL(wp) :: zddc, zdtc ! delta on corners and on centre127 REAL(wp) :: zdst ! shear at the center of the grid point128 REAL(wp) :: zdsshx, zdsshy ! term for the gradient of ocean surface129 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 130 129 131 130 REAL(wp) :: zresm ! Maximal error on ice velocity 132 REAL(wp) :: zdummy ! dummy argument133 131 REAL(wp) :: zintb, zintn ! dummy argument 134 132 … … 139 137 REAL(wp), POINTER, DIMENSION(:,:) :: zcorl1, zcorl2 ! coriolis parameter on U/V points 140 138 REAL(wp), POINTER, DIMENSION(:,:) :: za1ct , za2ct ! temporary arrays 141 REAL(wp), POINTER, DIMENSION(:,:) :: u_oce1, v_oce1! ocean u/v component on U points142 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 143 141 REAL(wp), POINTER, DIMENSION(:,:) :: u_ice2, v_ice1 ! ice u/v component on V/U point 144 142 REAL(wp), POINTER, DIMENSION(:,:) :: zf1 , zf2 ! arrays for internal stresses 143 REAL(wp), POINTER, DIMENSION(:,:) :: zmask ! mask ocean grid points 145 144 146 145 REAL(wp), POINTER, DIMENSION(:,:) :: zdt ! tension at centre of grid cells … … 152 151 ! ocean surface (ssh_m) if ice is not embedded 153 152 ! ice top surface if ice is embedded 153 154 REAL(wp), PARAMETER :: zepsi = 1.0e-20_wp ! tolerance parameter 155 REAL(wp), PARAMETER :: zvmin = 1.0e-03_wp ! ice volume below which ice velocity equals ocean velocity 154 156 !!------------------------------------------------------------------- 155 157 156 158 CALL wrk_alloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct ) 157 CALL wrk_alloc( jpi,jpj, u_oce 1, u_oce2, u_ice2, v_oce1 , v_oce2, v_ice1)159 CALL wrk_alloc( jpi,jpj, u_oce2, u_ice2, v_oce1 , v_ice1 , zmask ) 158 160 CALL wrk_alloc( jpi,jpj, zf1 , zu_ice, zf2 , zv_ice , zdt , zds ) 159 161 CALL wrk_alloc( jpi,jpj, zdt , zds , zs1 , zs2 , zs12 , zresr , zpice ) … … 161 163 #if defined key_lim2 && ! defined key_lim2_vp 162 164 # if defined key_agrif 163 USE ice_2, vt_s => hsnm164 USE ice_2, vt_i => hicm165 USE ice_2, vt_s => hsnm 166 USE ice_2, vt_i => hicm 165 167 # else 166 vt_s => hsnm167 vt_i => hicm168 vt_s => hsnm 169 vt_i => hicm 168 170 # endif 169 at_i(:,:) = 1. - frld(:,:)171 at_i(:,:) = 1. - frld(:,:) 170 172 #endif 171 173 #if defined key_agrif && defined key_lim2 172 CALL agrif_rhg_lim2_load ! First interpolation of coarse values174 CALL agrif_rhg_lim2_load ! First interpolation of coarse values 173 175 #endif 174 176 ! … … 186 188 187 189 #if defined key_lim3 188 CALL lim_itd_me_icestrength( ridge_scheme_swi ) ! LIM-3: Ice strength on T-points 189 #endif 190 191 !CDIR NOVERRCHK 190 CALL lim_itd_me_icestrength( nn_icestr ) ! LIM-3: Ice strength on T-points 191 #endif 192 192 193 DO jj = k_j1 , k_jpj ! Ice mass and temp variables 193 !CDIR NOVERRCHK194 194 DO ji = 1 , jpi 195 195 #if defined key_lim3 196 zpresh(ji,jj) = tm s(ji,jj) * strength(ji,jj)196 zpresh(ji,jj) = tmask(ji,jj,1) * strength(ji,jj) 197 197 #endif 198 198 #if defined key_lim2 199 zpresh(ji,jj) = tm s(ji,jj) * pstar * vt_i(ji,jj) * EXP( -c_rhg * (1. - at_i(ji,jj) ) )200 #endif 201 ! tmi= 1 where there is ice or on land202 tmi(ji,jj) = 1._wp - ( 1._wp - MAX( 0._wp , SIGN ( 1._wp , vt_i(ji,jj) - epsd ) ) ) * 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) 203 203 END DO 204 204 END DO … … 206 206 ! Ice strength on grid cell corners (zpreshc) 207 207 ! needed for calculation of shear stress 208 !CDIR NOVERRCHK209 208 DO jj = k_j1+1, k_jpj-1 210 !CDIR NOVERRCHK211 209 DO ji = 2, jpim1 !RB caution no fs_ (ji+1,jj+1) 212 zstms = tms(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + & 213 & tms(ji,jj+1) * wght(ji+1,jj+1,1,2) + & 214 & tms(ji+1,jj) * wght(ji+1,jj+1,2,1) + & 215 & tms(ji,jj) * wght(ji+1,jj+1,1,1) 216 zpreshc(ji,jj) = ( zpresh(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + & 217 & zpresh(ji,jj+1) * wght(ji+1,jj+1,1,2) + & 218 & zpresh(ji+1,jj) * wght(ji+1,jj+1,2,1) + & 219 & zpresh(ji,jj) * wght(ji+1,jj+1,1,1) & 220 & ) / MAX( zstms, epsd ) 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) 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) + & 213 & zpresh(ji+1,jj) * wght(ji+1,jj+1,2,1) + zpresh(ji,jj) * wght(ji+1,jj+1,1,1) & 214 & ) / MAX( zstms, zepsi ) 221 215 END DO 222 216 END DO 223 224 217 CALL lbc_lnk( zpreshc(:,:), 'F', 1. ) 225 218 ! … … 236 229 ! zcorl2: Coriolis parameter on V-points 237 230 ! (ztagnx,ztagny): wind stress on U/V points 238 ! u_oce1: ocean u component on u points239 231 ! v_oce1: ocean v component on u points 240 232 ! u_oce2: ocean u component on v points 241 ! v_oce2: ocean v component on v points242 233 243 234 IF( nn_ice_embd == 2 ) THEN !== embedded sea ice: compute representative ice top surface ==! 244 245 246 235 ! 236 ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1} 237 ! = (1/nn_fsbc)^2 * {SUM[n], n=0,nn_fsbc-1} 247 238 zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp 248 249 250 239 ! 240 ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1} 241 ! = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1}) 251 242 zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp 252 243 ! 253 244 zpice(:,:) = ssh_m(:,:) + ( zintn * snwice_mass(:,:) + zintb * snwice_mass_b(:,:) ) * r1_rau0 254 245 ! 255 246 ELSE !== non-embedded sea ice: use ocean surface for slope calculation ==! 256 247 zpice(:,:) = ssh_m(:,:) … … 260 251 DO ji = fs_2, fs_jpim1 261 252 262 zc1 = tm s(ji ,jj) * ( rhosn * vt_s(ji ,jj ) + rhoic * vt_i(ji ,jj ) )263 zc2 = tm s(ji+1,jj) * ( rhosn * vt_s(ji+1,jj ) + rhoic * vt_i(ji+1,jj ) )264 zc3 = tm s(ji ,jj+1) * ( rhosn * vt_s(ji ,jj+1) + rhoic * vt_i(ji ,jj+1) )265 266 zt11 = tm s(ji ,jj) * e1t(ji ,jj)267 zt12 = tm s(ji+1,jj) * e1t(ji+1,jj)268 zt21 = tm s(ji,jj) * e2t(ji,jj )269 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) 270 261 271 262 ! Leads area. 272 zfrld1(ji,jj) = ( zt12 * ( 1.0 - at_i(ji,jj) ) + zt11 * ( 1.0 - at_i(ji+1,jj) ) ) / ( zt11 + zt12 + epsd)273 zfrld2(ji,jj) = ( zt22 * ( 1.0 - at_i(ji,jj) ) + zt21 * ( 1.0 - at_i(ji,jj+1) ) ) / ( zt21 + zt22 + epsd)263 zfrld1(ji,jj) = ( zt12 * ( 1.0 - at_i(ji,jj) ) + zt11 * ( 1.0 - at_i(ji+1,jj) ) ) / ( zt11 + zt12 + zepsi ) 264 zfrld2(ji,jj) = ( zt22 * ( 1.0 - at_i(ji,jj) ) + zt21 * ( 1.0 - at_i(ji,jj+1) ) ) / ( zt21 + zt22 + zepsi ) 274 265 275 266 ! Mass, coriolis coeff. and currents 276 zmass1(ji,jj) = ( zt12 *zc1 + zt11*zc2 ) / (zt11+zt12+epsd)277 zmass2(ji,jj) = ( zt22 *zc1 + zt21*zc3 ) / (zt21+zt22+epsd)278 zcorl1(ji,jj) = zmass1(ji,jj) * ( e1t(ji+1,jj) *fcor(ji,jj) + e1t(ji,jj)*fcor(ji+1,jj) ) &279 & / ( e1t(ji,jj) + e1t(ji+1,jj) + epsd)280 zcorl2(ji,jj) = zmass2(ji,jj) * ( e2t(ji,jj+1) *fcor(ji,jj) + e2t(ji,jj)*fcor(ji,jj+1) ) &281 & / ( e2t(ji,jj+1) + e2t(ji,jj) + epsd)267 zmass1(ji,jj) = ( zt12 * zc1 + zt11 * zc2 ) / ( zt11 + zt12 + zepsi ) 268 zmass2(ji,jj) = ( zt22 * zc1 + zt21 * zc3 ) / ( zt21 + zt22 + zepsi ) 269 zcorl1(ji,jj) = zmass1(ji,jj) * ( e1t(ji+1,jj) * fcor(ji,jj) + e1t(ji,jj) * fcor(ji+1,jj) ) & 270 & / ( e1t(ji,jj) + e1t(ji+1,jj) + zepsi ) 271 zcorl2(ji,jj) = zmass2(ji,jj) * ( e2t(ji,jj+1) * fcor(ji,jj) + e2t(ji,jj) * fcor(ji,jj+1) ) & 272 & / ( e2t(ji,jj+1) + e2t(ji,jj) + zepsi ) 282 273 ! 283 u_oce1(ji,jj) = u_oce(ji,jj)284 v_oce2(ji,jj) = v_oce(ji,jj)285 286 274 ! Ocean has no slip boundary condition 287 v_oce1(ji,jj) = 0.5 *( (v_oce(ji,jj)+v_oce(ji,jj-1))*e1t(ji,jj)&288 & +(v_oce(ji+1,jj)+v_oce(ji+1,jj-1))*e1t(ji+1,jj))&289 & /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)290 291 u_oce2(ji,jj) = 0.5 *((u_oce(ji,jj)+u_oce(ji-1,jj))*e2t(ji,jj)&292 & +(u_oce(ji,jj+1)+u_oce(ji-1,jj+1))*e2t(ji,jj+1))&293 & / (e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)275 v_oce1(ji,jj) = 0.5 * ( ( v_oce(ji ,jj) + v_oce(ji ,jj-1) ) * e1t(ji,jj) & 276 & + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji+1,jj) ) & 277 & / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 278 279 u_oce2(ji,jj) = 0.5 * ( ( u_oce(ji,jj ) + u_oce(ji-1,jj ) ) * e2t(ji,jj) & 280 & + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj+1) ) & 281 & / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1) 294 282 295 283 ! Wind stress at U,V-point … … 303 291 ! include it later 304 292 305 zdsshx = ( zpice(ji+1,jj) - zpice(ji,jj) ) /e1u(ji,jj)306 zdsshy = ( zpice(ji,jj+1) - zpice(ji,jj) ) /e2v(ji,jj)293 zdsshx = ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj) 294 zdsshy = ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj) 307 295 308 296 za1ct(ji,jj) = ztagnx - zmass1(ji,jj) * grav * zdsshx … … 318 306 ! 319 307 ! Time step for subcycling 320 dtevp = rdt_ice / nevp 308 dtevp = rdt_ice / nn_nevp 309 #if defined key_lim3 310 dtotel = dtevp / ( 2._wp * rn_relast * rdt_ice ) 311 #else 321 312 dtotel = dtevp / ( 2._wp * telast ) 322 313 #endif 314 z1_dtotel = 1._wp / ( 1._wp + dtotel ) 315 z1_dtevp = 1._wp / dtevp 323 316 !-ecc2: square of yield ellipse eccenticrity (reminder: must become a namelist parameter) 324 ecc2 = ecc *ecc317 ecc2 = rn_ecc * rn_ecc 325 318 ecci = 1. / ecc2 326 319 … … 331 324 332 325 ! !----------------------! 333 DO jter = 1 , n evp! loop over jter !326 DO jter = 1 , nn_nevp ! loop over jter ! 334 327 ! !----------------------! 335 328 DO jj = k_j1, k_jpj-1 … … 339 332 340 333 DO jj = k_j1+1, k_jpj-1 341 DO ji = fs_2, jpim1 !RB bug no vect opt due to tmi334 DO ji = fs_2, fs_jpim1 !RB bug no vect opt due to zmask 342 335 343 336 ! … … 360 353 ! 361 354 ! 362 divu_i(ji,jj) = ( e2u(ji,jj)*u_ice(ji,jj) & 363 & -e2u(ji-1,jj)*u_ice(ji-1,jj) & 364 & +e1v(ji,jj)*v_ice(ji,jj) & 365 & -e1v(ji,jj-1)*v_ice(ji,jj-1) & 366 & ) & 367 & / area(ji,jj) 368 369 zdt(ji,jj) = ( ( u_ice(ji,jj)/e2u(ji,jj) & 370 & -u_ice(ji-1,jj)/e2u(ji-1,jj) & 371 & )*e2t(ji,jj)*e2t(ji,jj) & 372 & -( v_ice(ji,jj)/e1v(ji,jj) & 373 & -v_ice(ji,jj-1)/e1v(ji,jj-1) & 374 & )*e1t(ji,jj)*e1t(ji,jj) & 375 & ) & 376 & / area(ji,jj) 355 divu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) & 356 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) & 357 & ) * r1_e12t(ji,jj) 358 359 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) & 360 & - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) & 361 & ) * r1_e12t(ji,jj) 377 362 378 363 ! 379 zds(ji,jj) = ( ( u_ice(ji,jj+1)/e1u(ji,jj+1) & 380 & -u_ice(ji,jj)/e1u(ji,jj) & 381 & )*e1f(ji,jj)*e1f(ji,jj) & 382 & +( v_ice(ji+1,jj)/e2v(ji+1,jj) & 383 & -v_ice(ji,jj)/e2v(ji,jj) & 384 & )*e2f(ji,jj)*e2f(ji,jj) & 385 & ) & 386 & / ( e1f(ji,jj) * e2f(ji,jj) ) * ( 2.0 - tmf(ji,jj) ) & 387 & * tmi(ji,jj) * tmi(ji,jj+1) & 388 & * tmi(ji+1,jj) * tmi(ji+1,jj+1) 389 390 391 v_ice1(ji,jj) = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji+1,jj) & 392 & +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj)) & 393 & /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) 394 395 u_ice2(ji,jj) = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj+1) & 396 & +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj)) & 397 & /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj) 398 399 END DO 400 END DO 401 CALL lbc_lnk( v_ice1, 'U', -1. ) ; CALL lbc_lnk( u_ice2, 'V', -1. ) ! lateral boundary cond. 402 403 !CDIR NOVERRCHK 364 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) & 365 & + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 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) 368 369 370 v_ice1(ji,jj) = 0.5_wp * ( ( v_ice(ji ,jj) + v_ice(ji ,jj-1) ) * e1t(ji+1,jj) & 371 & + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji ,jj) ) & 372 & / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 373 374 u_ice2(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj+1) & 375 & + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj ) ) & 376 & / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1) 377 END DO 378 END DO 379 380 CALL lbc_lnk_multi( v_ice1, 'U', -1., u_ice2, 'V', -1. ) ! lateral boundary cond. 381 404 382 DO jj = k_j1+1, k_jpj-1 405 !CDIR NOVERRCHK406 383 DO ji = fs_2, fs_jpim1 407 384 408 385 !- Calculate Delta at centre of grid cells 409 zdst = ( e2u(ji , jj) * v_ice1(ji ,jj) & 410 & - e2u(ji-1, jj) * v_ice1(ji-1,jj) & 411 & + e1v(ji, jj ) * u_ice2(ji,jj ) & 412 & - e1v(ji, jj-1) * u_ice2(ji,jj-1) & 413 & ) & 414 & / area(ji,jj) 415 416 delta = SQRT( divu_i(ji,jj)*divu_i(ji,jj) + ( zdt(ji,jj)*zdt(ji,jj) + zdst*zdst ) * usecc2 ) 417 delta_i(ji,jj) = delta + creepl 418 !-Calculate stress tensor components zs1 and zs2 419 !-at centre of grid cells (see section 3.5 of CICE user's guide). 420 zs1(ji,jj) = ( zs1(ji,jj) + dtotel * ( ( divu_i(ji,jj) / delta_i(ji,jj) - delta / delta_i(ji,jj) ) & 421 & * zpresh(ji,jj) ) ) / ( 1._wp + dtotel ) 422 zs2(ji,jj) = ( zs2(ji,jj) + dtotel * ( ecci * zdt(ji,jj) / delta_i(ji,jj) * zpresh(ji,jj) ) ) & 423 & / ( 1._wp + dtotel ) 424 425 END DO 426 END DO 427 428 CALL lbc_lnk( zs1(:,:), 'T', 1. ) 429 CALL lbc_lnk( zs2(:,:), 'T', 1. ) 430 431 !CDIR NOVERRCHK 432 DO jj = k_j1+1, k_jpj-1 433 !CDIR NOVERRCHK 434 DO ji = fs_2, fs_jpim1 386 zdst = ( e2u(ji,jj) * v_ice1(ji,jj) - e2u(ji-1,jj ) * v_ice1(ji-1,jj ) & 387 & + e1v(ji,jj) * u_ice2(ji,jj) - e1v(ji ,jj-1) * u_ice2(ji ,jj-1) & 388 & ) * r1_e12t(ji,jj) 389 390 delta = SQRT( divu_i(ji,jj)**2 + ( zdt(ji,jj)**2 + zdst**2 ) * usecc2 ) 391 delta_i(ji,jj) = delta + rn_creepl 392 435 393 !- Calculate Delta on corners 436 zddc = ( ( v_ice1(ji,jj+1)/e1u(ji,jj+1) & 437 & -v_ice1(ji,jj)/e1u(ji,jj) & 438 & )*e1f(ji,jj)*e1f(ji,jj) & 439 & +( u_ice2(ji+1,jj)/e2v(ji+1,jj) & 440 & -u_ice2(ji,jj)/e2v(ji,jj) & 441 & )*e2f(ji,jj)*e2f(ji,jj) & 442 & ) & 443 & / ( e1f(ji,jj) * e2f(ji,jj) ) 444 445 zdtc = (-( v_ice1(ji,jj+1)/e1u(ji,jj+1) & 446 & -v_ice1(ji,jj)/e1u(ji,jj) & 447 & )*e1f(ji,jj)*e1f(ji,jj) & 448 & +( u_ice2(ji+1,jj)/e2v(ji+1,jj) & 449 & -u_ice2(ji,jj)/e2v(ji,jj) & 450 & )*e2f(ji,jj)*e2f(ji,jj) & 451 & ) & 452 & / ( e1f(ji,jj) * e2f(ji,jj) ) 453 454 zddc = SQRT(zddc**2+(zdtc**2+zds(ji,jj)**2)*usecc2) + creepl 455 456 !-Calculate stress tensor component zs12 at corners (see section 3.5 of CICE user's guide). 457 zs12(ji,jj) = ( zs12(ji,jj) + dtotel * & 458 & ( ecci * zds(ji,jj) / ( 2._wp * zddc ) * zpreshc(ji,jj) ) ) & 459 & / ( 1.0 + dtotel ) 460 461 END DO ! ji 462 END DO ! jj 463 464 CALL lbc_lnk( zs12(:,:), 'F', 1. ) 465 394 zddc = ( ( v_ice1(ji,jj+1) * r1_e1u(ji,jj+1) - v_ice1(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) & 395 & + ( u_ice2(ji+1,jj) * r1_e2v(ji+1,jj) - u_ice2(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 396 & ) * r1_e12f(ji,jj) 397 398 zdtc = (- ( v_ice1(ji,jj+1) * r1_e1u(ji,jj+1) - v_ice1(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) & 399 & + ( u_ice2(ji+1,jj) * r1_e2v(ji+1,jj) - u_ice2(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 400 & ) * r1_e12f(ji,jj) 401 402 zddc = SQRT( zddc**2 + ( zdtc**2 + zds(ji,jj)**2 ) * usecc2 ) + rn_creepl 403 404 !-Calculate stress tensor components zs1 and zs2 at centre of grid cells (see section 3.5 of CICE user's guide). 405 zs1(ji,jj) = ( zs1 (ji,jj) + dtotel * ( divu_i(ji,jj) - delta ) / delta_i(ji,jj) * zpresh(ji,jj) & 406 & ) * z1_dtotel 407 zs2(ji,jj) = ( zs2 (ji,jj) + dtotel * ecci * zdt(ji,jj) / delta_i(ji,jj) * zpresh(ji,jj) & 408 & ) * z1_dtotel 409 !-Calculate stress tensor component zs12 at corners 410 zs12(ji,jj) = ( zs12(ji,jj) + dtotel * ecci * zds(ji,jj) / ( 2._wp * zddc ) * zpreshc(ji,jj) & 411 & ) * z1_dtotel 412 413 END DO 414 END DO 415 416 CALL lbc_lnk_multi( zs1 , 'T', 1., zs2, 'T', 1., zs12, 'F', 1. ) 417 466 418 ! Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) 467 419 DO jj = k_j1+1, k_jpj-1 468 420 DO ji = fs_2, fs_jpim1 469 421 !- contribution of zs1, zs2 and zs12 to zf1 470 zf1(ji,jj) = 0.5 *( (zs1(ji+1,jj)-zs1(ji,jj))*e2u(ji,jj)&471 & +(zs2(ji+1,jj)*e2t(ji+1,jj)**2-zs2(ji,jj)*e2t(ji,jj)**2)/e2u(ji,jj)&472 & +2.0*(zs12(ji,jj)*e1f(ji,jj)**2-zs12(ji,jj-1)*e1f(ji,jj-1)**2)/e1u(ji,jj)&473 & ) / ( e1u(ji,jj)*e2u(ji,jj))422 zf1(ji,jj) = 0.5 * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj) & 423 & + ( zs2(ji+1,jj) * e2t(ji+1,jj)**2 - zs2(ji,jj) * e2t(ji,jj)**2 ) * r1_e2u(ji,jj) & 424 & + 2.0 * ( zs12(ji,jj) * e1f(ji,jj)**2 - zs12(ji,jj-1) * e1f(ji,jj-1)**2 ) * r1_e1u(ji,jj) & 425 & ) * r1_e12u(ji,jj) 474 426 ! contribution of zs1, zs2 and zs12 to zf2 475 zf2(ji,jj) = 0.5*( (zs1(ji,jj+1)-zs1(ji,jj))*e1v(ji,jj) & 476 & -(zs2(ji,jj+1)*e1t(ji,jj+1)**2 - zs2(ji,jj)*e1t(ji,jj)**2)/e1v(ji,jj) & 477 & + 2.0*(zs12(ji,jj)*e2f(ji,jj)**2 - & 478 zs12(ji-1,jj)*e2f(ji-1,jj)**2)/e2v(ji,jj) & 479 & ) / ( e1v(ji,jj)*e2v(ji,jj) ) 427 zf2(ji,jj) = 0.5 * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj) & 428 & - ( zs2(ji,jj+1) * e1t(ji,jj+1)**2 - zs2(ji,jj) * e1t(ji,jj)**2 ) * r1_e1v(ji,jj) & 429 & + 2.0 * ( zs12(ji,jj) * e2f(ji,jj)**2 - zs12(ji-1,jj) * e2f(ji-1,jj)**2 ) * r1_e2v(ji,jj) & 430 & ) * r1_e12v(ji,jj) 480 431 END DO 481 432 END DO … … 487 438 IF (MOD(jter,2).eq.0) THEN 488 439 489 !CDIR NOVERRCHK490 440 DO jj = k_j1+1, k_jpj-1 491 !CDIR NOVERRCHK492 441 DO ji = fs_2, fs_jpim1 493 zmask = (1.0-MAX(0._wp,SIGN(1._wp,-zmass1(ji,jj))))*tmu(ji,jj)494 z0 = zmass1(ji,jj) /dtevp442 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) ) ) * umask(ji,jj,1) 443 z0 = zmass1(ji,jj) * z1_dtevp 495 444 496 445 ! SB modif because ocean has no slip boundary condition 497 zv_ice1 = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji,jj) & 498 & +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji+1,jj)) & 499 & /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) 500 za = rhoco*SQRT((u_ice(ji,jj)-u_oce1(ji,jj))**2 + & 501 (zv_ice1-v_oce1(ji,jj))**2) * (1.0-zfrld1(ji,jj)) 502 zr = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + & 503 za*(u_oce1(ji,jj)) 504 zcca = z0+za 446 zv_ice1 = 0.5 * ( ( v_ice(ji ,jj) + v_ice(ji ,jj-1) ) * e1t(ji ,jj) & 447 & + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji+1,jj) ) & 448 & / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 449 za = rhoco * SQRT( ( u_ice(ji,jj) - u_oce(ji,jj) )**2 + & 450 & ( zv_ice1 - v_oce1(ji,jj) )**2 ) * ( 1.0 - zfrld1(ji,jj) ) 451 zr = z0 * u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + za * u_oce(ji,jj) 452 zcca = z0 + za 505 453 zccb = zcorl1(ji,jj) 506 u_ice(ji,jj) = (zr+zccb*zv_ice1)/(zcca+epsd)*zmask 507 454 u_ice(ji,jj) = ( zr + zccb * zv_ice1 ) / ( zcca + zepsi ) * rswitch 508 455 END DO 509 456 END DO … … 511 458 CALL lbc_lnk( u_ice(:,:), 'U', -1. ) 512 459 #if defined key_agrif && defined key_lim2 513 CALL agrif_rhg_lim2( jter, n evp, 'U' )460 CALL agrif_rhg_lim2( jter, nn_nevp, 'U' ) 514 461 #endif 515 462 #if defined key_bdy … … 517 464 #endif 518 465 519 !CDIR NOVERRCHK520 466 DO jj = k_j1+1, k_jpj-1 521 !CDIR NOVERRCHK522 467 DO ji = fs_2, fs_jpim1 523 468 524 zmask = (1.0-MAX(0._wp,SIGN(1._wp,-zmass2(ji,jj))))*tmv(ji,jj)525 z0 = zmass2(ji,jj) /dtevp469 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) ) * vmask(ji,jj,1) 470 z0 = zmass2(ji,jj) * z1_dtevp 526 471 ! SB modif because ocean has no slip boundary condition 527 zu_ice2 = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj) & 528 & + (u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj+1)) & 529 & /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj) 530 za = rhoco*SQRT((zu_ice2-u_oce2(ji,jj))**2 + & 531 (v_ice(ji,jj)-v_oce2(ji,jj))**2)*(1.0-zfrld2(ji,jj)) 532 zr = z0*v_ice(ji,jj) + zf2(ji,jj) + & 533 za2ct(ji,jj) + za*(v_oce2(ji,jj)) 534 zcca = z0+za 472 zu_ice2 = 0.5 * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj) & 473 & + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj+1) ) & 474 & / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1) 475 za = rhoco * SQRT( ( zu_ice2 - u_oce2(ji,jj) )**2 + & 476 & ( v_ice(ji,jj) - v_oce(ji,jj))**2 ) * ( 1.0 - zfrld2(ji,jj) ) 477 zr = z0 * v_ice(ji,jj) + zf2(ji,jj) + za2ct(ji,jj) + za * v_oce(ji,jj) 478 zcca = z0 + za 535 479 zccb = zcorl2(ji,jj) 536 v_ice(ji,jj) = (zr-zccb*zu_ice2)/(zcca+epsd)*zmask 537 480 v_ice(ji,jj) = ( zr - zccb * zu_ice2 ) / ( zcca + zepsi ) * rswitch 538 481 END DO 539 482 END DO … … 541 484 CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 542 485 #if defined key_agrif && defined key_lim2 543 CALL agrif_rhg_lim2( jter, n evp, 'V' )486 CALL agrif_rhg_lim2( jter, nn_nevp, 'V' ) 544 487 #endif 545 488 #if defined key_bdy … … 548 491 549 492 ELSE 550 !CDIR NOVERRCHK551 493 DO jj = k_j1+1, k_jpj-1 552 !CDIR NOVERRCHK553 494 DO ji = fs_2, fs_jpim1 554 zmask = (1.0-MAX(0._wp,SIGN(1._wp,-zmass2(ji,jj))))*tmv(ji,jj)555 z0 = zmass2(ji,jj) /dtevp495 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) ) * vmask(ji,jj,1) 496 z0 = zmass2(ji,jj) * z1_dtevp 556 497 ! SB modif because ocean has no slip boundary condition 557 zu_ice2 = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj) & 558 & +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj+1)) & 559 & /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj) 560 561 za = rhoco*SQRT((zu_ice2-u_oce2(ji,jj))**2 + & 562 (v_ice(ji,jj)-v_oce2(ji,jj))**2)*(1.0-zfrld2(ji,jj)) 563 zr = z0*v_ice(ji,jj) + zf2(ji,jj) + & 564 za2ct(ji,jj) + za*(v_oce2(ji,jj)) 565 zcca = z0+za 498 zu_ice2 = 0.5 * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj) & 499 & +( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj+1) ) & 500 & / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1) 501 502 za = rhoco * SQRT( ( zu_ice2 - u_oce2(ji,jj) )**2 + & 503 & ( v_ice(ji,jj) - v_oce(ji,jj) )**2 ) * ( 1.0 - zfrld2(ji,jj) ) 504 zr = z0 * v_ice(ji,jj) + zf2(ji,jj) + za2ct(ji,jj) + za * v_oce(ji,jj) 505 zcca = z0 + za 566 506 zccb = zcorl2(ji,jj) 567 v_ice(ji,jj) = (zr-zccb*zu_ice2)/(zcca+epsd)*zmask 568 507 v_ice(ji,jj) = ( zr - zccb * zu_ice2 ) / ( zcca + zepsi ) * rswitch 569 508 END DO 570 509 END DO … … 572 511 CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 573 512 #if defined key_agrif && defined key_lim2 574 CALL agrif_rhg_lim2( jter, n evp, 'V' )513 CALL agrif_rhg_lim2( jter, nn_nevp, 'V' ) 575 514 #endif 576 515 #if defined key_bdy … … 578 517 #endif 579 518 580 !CDIR NOVERRCHK581 519 DO jj = k_j1+1, k_jpj-1 582 !CDIR NOVERRCHK583 520 DO ji = fs_2, fs_jpim1 584 zmask = (1.0-MAX(0._wp,SIGN(1._wp,-zmass1(ji,jj))))*tmu(ji,jj) 585 z0 = zmass1(ji,jj)/dtevp 586 zv_ice1 = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji,jj) & 587 & +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji+1,jj)) & 588 & /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) 589 590 za = rhoco*SQRT((u_ice(ji,jj)-u_oce1(ji,jj))**2 + & 591 (zv_ice1-v_oce1(ji,jj))**2)*(1.0-zfrld1(ji,jj)) 592 zr = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + & 593 za*(u_oce1(ji,jj)) 594 zcca = z0+za 521 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) ) ) * umask(ji,jj,1) 522 z0 = zmass1(ji,jj) * z1_dtevp 523 zv_ice1 = 0.5 * ( ( v_ice(ji ,jj) + v_ice(ji ,jj-1) ) * e1t(ji,jj) & 524 & + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji+1,jj) ) & 525 & / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 526 527 za = rhoco * SQRT( ( u_ice(ji,jj) - u_oce(ji,jj) )**2 + & 528 & ( zv_ice1 - v_oce1(ji,jj) )**2 ) * ( 1.0 - zfrld1(ji,jj) ) 529 zr = z0 * u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + za * u_oce(ji,jj) 530 zcca = z0 + za 595 531 zccb = zcorl1(ji,jj) 596 u_ice(ji,jj) = ( zr+zccb*zv_ice1)/(zcca+epsd)*zmask597 END DO ! ji598 END DO ! jj532 u_ice(ji,jj) = ( zr + zccb * zv_ice1 ) / ( zcca + zepsi ) * rswitch 533 END DO 534 END DO 599 535 600 536 CALL lbc_lnk( u_ice(:,:), 'U', -1. ) 601 537 #if defined key_agrif && defined key_lim2 602 CALL agrif_rhg_lim2( jter, n evp, 'U' )538 CALL agrif_rhg_lim2( jter, nn_nevp, 'U' ) 603 539 #endif 604 540 #if defined key_bdy … … 611 547 !--- Convergence test. 612 548 DO jj = k_j1+1 , k_jpj-1 613 zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ) , & 614 ABS( v_ice(:,jj) - zv_ice(:,jj) ) ) 615 END DO 616 zresm = MAXVAL( zresr( 1:jpi , k_j1+1:k_jpj-1 ) ) 549 zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) ) 550 END DO 551 zresm = MAXVAL( zresr( 1:jpi, k_j1+1:k_jpj-1 ) ) 617 552 IF( lk_mpp ) CALL mpp_max( zresm ) ! max over the global domain 618 553 ENDIF … … 625 560 ! 4) Prevent ice velocities when the ice is thin 626 561 !------------------------------------------------------------------------------! 627 ! If the ice thickness is below hminrhg (5cm) then ice velocity should equal the 628 ! ocean velocity, 629 ! This prevents high velocity when ice is thin 630 !CDIR NOVERRCHK 562 ! If the ice volume is below zvmin then ice velocity should equal the 563 ! ocean velocity. This prevents high velocity when ice is thin 631 564 DO jj = k_j1+1, k_jpj-1 632 !CDIR NOVERRCHK633 565 DO ji = fs_2, fs_jpim1 634 zdummy = vt_i(ji,jj) 635 IF ( zdummy .LE. hminrhg ) THEN 566 IF ( vt_i(ji,jj) <= zvmin ) THEN 636 567 u_ice(ji,jj) = u_oce(ji,jj) 637 568 v_ice(ji,jj) = v_oce(ji,jj) 638 ENDIF ! zdummy569 ENDIF 639 570 END DO 640 571 END DO 641 572 642 CALL lbc_lnk ( u_ice(:,:), 'U', -1. )643 CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 573 CALL lbc_lnk_multi( u_ice(:,:), 'U', -1., v_ice(:,:), 'V', -1. ) 574 644 575 #if defined key_agrif && defined key_lim2 645 CALL agrif_rhg_lim2( n evp ,nevp, 'U' )646 CALL agrif_rhg_lim2( n evp ,nevp, 'V' )576 CALL agrif_rhg_lim2( nn_nevp , nn_nevp, 'U' ) 577 CALL agrif_rhg_lim2( nn_nevp , nn_nevp, 'V' ) 647 578 #endif 648 579 #if defined key_bdy … … 653 584 DO jj = k_j1+1, k_jpj-1 654 585 DO ji = fs_2, fs_jpim1 655 zdummy = vt_i(ji,jj) 656 IF ( zdummy .LE. hminrhg ) THEN 657 v_ice1(ji,jj) = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji+1,jj) & 658 & +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj)) & 659 & /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) 660 661 u_ice2(ji,jj) = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj+1) & 662 & +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj)) & 663 & /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj) 664 ENDIF ! zdummy 586 IF ( vt_i(ji,jj) <= zvmin ) THEN 587 v_ice1(ji,jj) = 0.5_wp * ( ( v_ice(ji ,jj) + v_ice(ji, jj-1) ) * e1t(ji+1,jj) & 588 & + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji ,jj) ) & 589 & / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 590 591 u_ice2(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj+1) & 592 & + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj ) ) & 593 & / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1) 594 ENDIF 665 595 END DO 666 596 END DO 667 597 668 CALL lbc_lnk( u_ice2(:,:), 'V', -1. ) 669 CALL lbc_lnk( v_ice1(:,:), 'U', -1. ) 598 CALL lbc_lnk_multi( u_ice2(:,:), 'V', -1., v_ice1(:,:), 'U', -1. ) 670 599 671 600 ! Recompute delta, shear and div, inputs for mechanical redistribution 672 !CDIR NOVERRCHK673 601 DO jj = k_j1+1, k_jpj-1 674 !CDIR NOVERRCHK 675 DO ji = fs_2, jpim1 !RB bug no vect opt due to tmi 602 DO ji = fs_2, jpim1 !RB bug no vect opt due to zmask 676 603 !- divu_i(:,:), zdt(:,:): divergence and tension at centre 677 604 !- zds(:,:): shear on northeast corner of grid cells 678 zdummy = vt_i(ji,jj) 679 IF ( zdummy .LE. hminrhg ) THEN 680 681 divu_i(ji,jj) = ( e2u(ji,jj)*u_ice(ji,jj) & 682 & -e2u(ji-1,jj)*u_ice(ji-1,jj) & 683 & +e1v(ji,jj)*v_ice(ji,jj) & 684 & -e1v(ji,jj-1)*v_ice(ji,jj-1) & 685 & ) & 686 & / area(ji,jj) 687 688 zdt(ji,jj) = ( ( u_ice(ji,jj)/e2u(ji,jj) & 689 & -u_ice(ji-1,jj)/e2u(ji-1,jj) & 690 & )*e2t(ji,jj)*e2t(ji,jj) & 691 & -( v_ice(ji,jj)/e1v(ji,jj) & 692 & -v_ice(ji,jj-1)/e1v(ji,jj-1) & 693 & )*e1t(ji,jj)*e1t(ji,jj) & 694 & ) & 695 & / area(ji,jj) 605 IF ( vt_i(ji,jj) <= zvmin ) THEN 606 607 divu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj ) * u_ice(ji-1,jj ) & 608 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji ,jj-1) * v_ice(ji ,jj-1) & 609 & ) * r1_e12t(ji,jj) 610 611 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) & 612 & -( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) & 613 & ) * r1_e12t(ji,jj) 696 614 ! 697 615 ! SB modif because ocean has no slip boundary condition 698 zds(ji,jj) = ( ( u_ice(ji,jj+1) / e1u(ji,jj+1) & 699 & - u_ice(ji,jj) / e1u(ji,jj) ) & 700 & * e1f(ji,jj) * e1f(ji,jj) & 701 & + ( v_ice(ji+1,jj) / e2v(ji+1,jj) & 702 & - v_ice(ji,jj) / e2v(ji,jj) ) & 703 & * e2f(ji,jj) * e2f(ji,jj) ) & 704 & / ( e1f(ji,jj) * e2f(ji,jj) ) * ( 2.0 - tmf(ji,jj) ) & 705 & * tmi(ji,jj) * tmi(ji,jj+1) & 706 & * tmi(ji+1,jj) * tmi(ji+1,jj+1) 707 708 zdst = ( e2u( ji , jj ) * v_ice1(ji ,jj ) & 709 & - e2u( ji-1, jj ) * v_ice1(ji-1,jj ) & 710 & + e1v( ji , jj ) * u_ice2(ji ,jj ) & 711 & - e1v( ji , jj-1 ) * u_ice2(ji ,jj-1) ) / area(ji,jj) 712 713 delta = SQRT( divu_i(ji,jj)*divu_i(ji,jj) + ( zdt(ji,jj)*zdt(ji,jj) + zdst*zdst ) * usecc2 ) 714 delta_i(ji,jj) = delta + creepl 616 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) & 617 & +( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 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) 620 621 zdst = ( e2u(ji,jj) * v_ice1(ji,jj) - e2u(ji-1,jj ) * v_ice1(ji-1,jj ) & 622 & + e1v(ji,jj) * u_ice2(ji,jj) - e1v(ji ,jj-1) * u_ice2(ji ,jj-1) ) * r1_e12t(ji,jj) 623 624 delta = SQRT( divu_i(ji,jj)**2 + ( zdt(ji,jj)**2 + zdst**2 ) * usecc2 ) 625 delta_i(ji,jj) = delta + rn_creepl 715 626 716 ENDIF ! zdummy 717 718 END DO !jj 719 END DO !ji 627 ENDIF 628 END DO 629 END DO 720 630 ! 721 631 !------------------------------------------------------------------------------! 722 632 ! 5) Store stress tensor and its invariants 723 633 !------------------------------------------------------------------------------! 724 !725 634 ! * Invariants of the stress tensor are required for limitd_me 726 635 ! (accelerates convergence and improves stability) 727 636 DO jj = k_j1+1, k_jpj-1 728 637 DO ji = fs_2, fs_jpim1 729 ! begin TECLIM change 730 zdst= ( e2u( ji , jj ) * v_ice1(ji,jj) & 731 & - e2u( ji-1, jj ) * v_ice1(ji-1,jj) & 732 & + e1v( ji , jj ) * u_ice2(ji,jj) & 733 & - e1v( ji , jj-1 ) * u_ice2(ji,jj-1) ) / area(ji,jj) 638 zdst = ( e2u(ji,jj) * v_ice1(ji,jj) - e2u( ji-1, jj ) * v_ice1(ji-1,jj) & 639 & + e1v(ji,jj) * u_ice2(ji,jj) - e1v( ji , jj-1 ) * u_ice2(ji,jj-1) ) * r1_e12t(ji,jj) 734 640 shear_i(ji,jj) = SQRT( zdt(ji,jj) * zdt(ji,jj) + zdst * zdst ) 735 ! end TECLIM change736 641 END DO 737 642 END DO 738 643 739 644 ! Lateral boundary condition 740 CALL lbc_lnk( divu_i (:,:), 'T', 1. ) 741 CALL lbc_lnk( delta_i(:,:), 'T', 1. ) 742 ! CALL lbc_lnk( shear_i(:,:), 'F', 1. ) 743 CALL lbc_lnk( shear_i(:,:), 'T', 1. ) 645 CALL lbc_lnk_multi( divu_i (:,:), 'T', 1., delta_i(:,:), 'T', 1., shear_i(:,:), 'T', 1. ) 744 646 745 647 ! * Store the stress tensor for the next time step … … 772 674 DO jj = k_j1+1, k_jpj-1 773 675 DO ji = 2, jpim1 774 IF (zpresh(ji,jj) .GT.1.0) THEN676 IF (zpresh(ji,jj) > 1.0) THEN 775 677 sigma1 = ( zs1(ji,jj) + (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) ) 776 678 sigma2 = ( zs1(ji,jj) - (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) ) … … 786 688 ! 787 689 CALL wrk_dealloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct ) 788 CALL wrk_dealloc( jpi,jpj, u_oce 1, u_oce2, u_ice2, v_oce1 , v_oce2, v_ice1)690 CALL wrk_dealloc( jpi,jpj, u_oce2, u_ice2, v_oce1 , v_ice1 , zmask ) 789 691 CALL wrk_dealloc( jpi,jpj, zf1 , zu_ice, zf2 , zv_ice , zdt , zds ) 790 692 CALL wrk_dealloc( jpi,jpj, zdt , zds , zs1 , zs2 , zs12 , zresr , zpice )
Note: See TracChangeset
for help on using the changeset viewer.