- Timestamp:
- 2020-09-15T12:49:18+02:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/temporary_r4_trunk/src/ICE/icedyn_rhg_evp.F90
r13467 r13469 182 182 ! for diagnostics and convergence tests 183 183 ALLOCATE( zmsk00(jpi,jpj), zmsk15(jpi,jpj) ) 184 DO jj = 1, jpj 185 DO ji = 1, jpi 186 zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice , 0 if no ice 187 zmsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less 188 END DO 189 END DO 184 DO_2D_11_11 185 zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice , 0 if no ice 186 zmsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less 187 END_2D 190 188 ! 191 189 !!gm for Clem: OPTIMIZATION: I think zfmask can be computed one for all at the initialization.... … … 194 192 !------------------------------------------------------------------------------! 195 193 ! ocean/land mask 196 DO jj = 1, jpjm1 197 DO ji = 1, jpim1 ! NO vector opt. 198 zfmask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1) 199 END DO 200 END DO 194 DO_2D_10_10 195 zfmask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1) 196 END_2D 201 197 CALL lbc_lnk( 'icedyn_rhg_evp', zfmask, 'F', 1._wp ) 202 198 203 199 ! Lateral boundary conditions on velocity (modify zfmask) 204 DO jj = 2, jpjm1 205 DO ji = fs_2, fs_jpim1 ! vector opt. 206 IF( zfmask(ji,jj) == 0._wp ) THEN 207 zfmask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(ji,jj,1), umask(ji,jj+1,1), & 208 & vmask(ji,jj,1), vmask(ji+1,jj,1) ) ) 209 ENDIF 210 END DO 211 END DO 200 DO_2D_00_00 201 IF( zfmask(ji,jj) == 0._wp ) THEN 202 zfmask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(ji,jj,1), umask(ji,jj+1,1), & 203 & vmask(ji,jj,1), vmask(ji+1,jj,1) ) ) 204 ENDIF 205 END_2D 212 206 DO jj = 2, jpjm1 213 207 IF( zfmask(1,jj) == 0._wp ) THEN … … 272 266 zsshdyn(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b) 273 267 274 DO jj = 2, jpjm1 275 DO ji = fs_2, fs_jpim1 276 277 ! ice fraction at U-V points 278 zaU(ji,jj) = 0.5_wp * ( at_i(ji,jj) * e1e2t(ji,jj) + at_i(ji+1,jj) * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1) 279 zaV(ji,jj) = 0.5_wp * ( at_i(ji,jj) * e1e2t(ji,jj) + at_i(ji,jj+1) * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1) 280 281 ! Ice/snow mass at U-V points 282 zm1 = ( rhos * vt_s(ji ,jj ) + rhoi * vt_i(ji ,jj ) ) 283 zm2 = ( rhos * vt_s(ji+1,jj ) + rhoi * vt_i(ji+1,jj ) ) 284 zm3 = ( rhos * vt_s(ji ,jj+1) + rhoi * vt_i(ji ,jj+1) ) 285 zmassU = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1) 286 zmassV = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1) 287 288 ! Ocean currents at U-V points 289 v_oceU(ji,jj) = 0.25_wp * ( v_oce(ji,jj) + v_oce(ji,jj-1) + v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * umask(ji,jj,1) 290 u_oceV(ji,jj) = 0.25_wp * ( u_oce(ji,jj) + u_oce(ji-1,jj) + u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * vmask(ji,jj,1) 291 292 ! Coriolis at T points (m*f) 293 zmf(ji,jj) = zm1 * ff_t(ji,jj) 294 295 ! dt/m at T points (for alpha and beta coefficients) 296 zdt_m(ji,jj) = zdtevp / MAX( zm1, zmmin ) 297 298 ! m/dt 299 zmU_t(ji,jj) = zmassU * z1_dtevp 300 zmV_t(ji,jj) = zmassV * z1_dtevp 301 302 ! Drag ice-atm. 303 ztaux_ai(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj) 304 ztauy_ai(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj) 305 306 ! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points 307 zspgU(ji,jj) = - zmassU * grav * ( zsshdyn(ji+1,jj) - zsshdyn(ji,jj) ) * r1_e1u(ji,jj) 308 zspgV(ji,jj) = - zmassV * grav * ( zsshdyn(ji,jj+1) - zsshdyn(ji,jj) ) * r1_e2v(ji,jj) 309 310 ! masks 311 zmsk00x(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) ) ! 0 if no ice 312 zmsk00y(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) ) ! 0 if no ice 313 314 ! switches 315 IF( zmassU <= zmmin .AND. zaU(ji,jj) <= zamin ) THEN ; zmsk01x(ji,jj) = 0._wp 316 ELSE ; zmsk01x(ji,jj) = 1._wp ; ENDIF 317 IF( zmassV <= zmmin .AND. zaV(ji,jj) <= zamin ) THEN ; zmsk01y(ji,jj) = 0._wp 318 ELSE ; zmsk01y(ji,jj) = 1._wp ; ENDIF 319 320 END DO 321 END DO 268 DO_2D_00_00 269 270 ! ice fraction at U-V points 271 zaU(ji,jj) = 0.5_wp * ( at_i(ji,jj) * e1e2t(ji,jj) + at_i(ji+1,jj) * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1) 272 zaV(ji,jj) = 0.5_wp * ( at_i(ji,jj) * e1e2t(ji,jj) + at_i(ji,jj+1) * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1) 273 274 ! Ice/snow mass at U-V points 275 zm1 = ( rhos * vt_s(ji ,jj ) + rhoi * vt_i(ji ,jj ) ) 276 zm2 = ( rhos * vt_s(ji+1,jj ) + rhoi * vt_i(ji+1,jj ) ) 277 zm3 = ( rhos * vt_s(ji ,jj+1) + rhoi * vt_i(ji ,jj+1) ) 278 zmassU = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1) 279 zmassV = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1) 280 281 ! Ocean currents at U-V points 282 v_oceU(ji,jj) = 0.25_wp * ( v_oce(ji,jj) + v_oce(ji,jj-1) + v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * umask(ji,jj,1) 283 u_oceV(ji,jj) = 0.25_wp * ( u_oce(ji,jj) + u_oce(ji-1,jj) + u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * vmask(ji,jj,1) 284 285 ! Coriolis at T points (m*f) 286 zmf(ji,jj) = zm1 * ff_t(ji,jj) 287 288 ! dt/m at T points (for alpha and beta coefficients) 289 zdt_m(ji,jj) = zdtevp / MAX( zm1, zmmin ) 290 291 ! m/dt 292 zmU_t(ji,jj) = zmassU * z1_dtevp 293 zmV_t(ji,jj) = zmassV * z1_dtevp 294 295 ! Drag ice-atm. 296 ztaux_ai(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj) 297 ztauy_ai(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj) 298 299 ! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points 300 zspgU(ji,jj) = - zmassU * grav * ( zsshdyn(ji+1,jj) - zsshdyn(ji,jj) ) * r1_e1u(ji,jj) 301 zspgV(ji,jj) = - zmassV * grav * ( zsshdyn(ji,jj+1) - zsshdyn(ji,jj) ) * r1_e2v(ji,jj) 302 303 ! masks 304 zmsk00x(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) ) ! 0 if no ice 305 zmsk00y(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) ) ! 0 if no ice 306 307 ! switches 308 IF( zmassU <= zmmin .AND. zaU(ji,jj) <= zamin ) THEN ; zmsk01x(ji,jj) = 0._wp 309 ELSE ; zmsk01x(ji,jj) = 1._wp ; ENDIF 310 IF( zmassV <= zmmin .AND. zaV(ji,jj) <= zamin ) THEN ; zmsk01y(ji,jj) = 0._wp 311 ELSE ; zmsk01y(ji,jj) = 1._wp ; ENDIF 312 313 END_2D 322 314 CALL lbc_lnk_multi( 'icedyn_rhg_evp', zmf, 'T', 1., zdt_m, 'T', 1. ) 323 315 ! … … 325 317 ! 326 318 IF( ln_landfast_L16 ) THEN !-- Lemieux 2016 327 DO jj = 2, jpjm1 328 DO ji = fs_2, fs_jpim1 329 ! ice thickness at U-V points 330 zvU = 0.5_wp * ( vt_i(ji,jj) * e1e2t(ji,jj) + vt_i(ji+1,jj) * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1) 331 zvV = 0.5_wp * ( vt_i(ji,jj) * e1e2t(ji,jj) + vt_i(ji,jj+1) * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1) 332 ! ice-bottom stress at U points 333 zvCr = zaU(ji,jj) * rn_lf_depfra * hu_n(ji,jj) 334 ztaux_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) ) 335 ! ice-bottom stress at V points 336 zvCr = zaV(ji,jj) * rn_lf_depfra * hv_n(ji,jj) 337 ztauy_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) ) 338 ! ice_bottom stress at T points 339 zvCr = at_i(ji,jj) * rn_lf_depfra * ht_n(ji,jj) 340 tau_icebfr(ji,jj) = - rn_lf_bfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) ) 341 END DO 342 END DO 319 DO_2D_00_00 320 ! ice thickness at U-V points 321 zvU = 0.5_wp * ( vt_i(ji,jj) * e1e2t(ji,jj) + vt_i(ji+1,jj) * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1) 322 zvV = 0.5_wp * ( vt_i(ji,jj) * e1e2t(ji,jj) + vt_i(ji,jj+1) * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1) 323 ! ice-bottom stress at U points 324 zvCr = zaU(ji,jj) * rn_lf_depfra * hu_n(ji,jj) 325 ztaux_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) ) 326 ! ice-bottom stress at V points 327 zvCr = zaV(ji,jj) * rn_lf_depfra * hv_n(ji,jj) 328 ztauy_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) ) 329 ! ice_bottom stress at T points 330 zvCr = at_i(ji,jj) * rn_lf_depfra * ht_n(ji,jj) 331 tau_icebfr(ji,jj) = - rn_lf_bfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) ) 332 END_2D 343 333 CALL lbc_lnk( 'icedyn_rhg_evp', tau_icebfr(:,:), 'T', 1. ) 344 334 ! 345 335 ELSE !-- no landfast 346 DO jj = 2, jpjm1 347 DO ji = fs_2, fs_jpim1 348 ztaux_base(ji,jj) = 0._wp 349 ztauy_base(ji,jj) = 0._wp 350 END DO 351 END DO 336 DO_2D_00_00 337 ztaux_base(ji,jj) = 0._wp 338 ztauy_base(ji,jj) = 0._wp 339 END_2D 352 340 ENDIF 353 341 … … 363 351 ! convergence test 364 352 IF( nn_rhg_chkcvg == 1 .OR. nn_rhg_chkcvg == 2 ) THEN 365 DO jj = 1, jpj 366 DO ji = 1, jpi 367 zu_ice(ji,jj) = u_ice(ji,jj) * umask(ji,jj,1) ! velocity at previous time step 368 zv_ice(ji,jj) = v_ice(ji,jj) * vmask(ji,jj,1) 369 END DO 370 END DO 353 DO_2D_11_11 354 zu_ice(ji,jj) = u_ice(ji,jj) * umask(ji,jj,1) ! velocity at previous time step 355 zv_ice(ji,jj) = v_ice(ji,jj) * vmask(ji,jj,1) 356 END_2D 371 357 ENDIF 372 358 373 359 ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- ! 374 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 375 DO ji = 1, jpim1 376 377 ! shear at F points 378 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) & 379 & + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 380 & ) * r1_e1e2f(ji,jj) * zfmask(ji,jj) 381 382 END DO 383 END DO 360 DO_2D_10_10 361 362 ! shear at F points 363 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) & 364 & + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 365 & ) * r1_e1e2f(ji,jj) * zfmask(ji,jj) 366 367 END_2D 384 368 CALL lbc_lnk( 'icedyn_rhg_evp', zds, 'F', 1. ) 385 369 386 DO jj = 2, jpj ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12 387 DO ji = 2, jpi ! no vector loop 388 389 ! shear**2 at T points (doc eq. A16) 390 zds2 = ( zds(ji,jj ) * zds(ji,jj ) * e1e2f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e1e2f(ji-1,jj ) & 391 & + zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1) & 392 & ) * 0.25_wp * r1_e1e2t(ji,jj) 393 394 ! divergence at T points 395 zdiv = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) & 396 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) & 397 & ) * r1_e1e2t(ji,jj) 398 zdiv2 = zdiv * zdiv 399 400 ! tension at T points 401 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) & 402 & - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) & 403 & ) * r1_e1e2t(ji,jj) 404 zdt2 = zdt * zdt 405 406 ! delta at T points 407 zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 ) 408 409 ! P/delta at T points 410 zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl ) 411 412 ! alpha for aEVP 413 ! gamma = 0.5*P/(delta+creepl) * (c*pi)**2/Area * dt/m 414 ! alpha = beta = sqrt(4*gamma) 415 IF( ln_aEVP ) THEN 416 zalph1 = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) ) 417 z1_alph1 = 1._wp / ( zalph1 + 1._wp ) 418 zalph2 = zalph1 419 z1_alph2 = z1_alph1 420 ! explicit: 421 ! z1_alph1 = 1._wp / zalph1 422 ! z1_alph2 = 1._wp / zalph1 423 ! zalph1 = zalph1 - 1._wp 424 ! zalph2 = zalph1 425 ENDIF 426 427 ! stress at T points (zkt/=0 if landfast) 428 zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv * (1._wp + zkt) - zdelta * (1._wp - zkt) ) ) * z1_alph1 429 zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2 430 431 END DO 432 END DO 370 DO_2D_01_01 371 372 ! shear**2 at T points (doc eq. A16) 373 zds2 = ( zds(ji,jj ) * zds(ji,jj ) * e1e2f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e1e2f(ji-1,jj ) & 374 & + zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1) & 375 & ) * 0.25_wp * r1_e1e2t(ji,jj) 376 377 ! divergence at T points 378 zdiv = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) & 379 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) & 380 & ) * r1_e1e2t(ji,jj) 381 zdiv2 = zdiv * zdiv 382 383 ! tension at T points 384 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) & 385 & - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) & 386 & ) * r1_e1e2t(ji,jj) 387 zdt2 = zdt * zdt 388 389 ! delta at T points 390 zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 ) 391 392 ! P/delta at T points 393 zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl ) 394 395 ! alpha for aEVP 396 ! gamma = 0.5*P/(delta+creepl) * (c*pi)**2/Area * dt/m 397 ! alpha = beta = sqrt(4*gamma) 398 IF( ln_aEVP ) THEN 399 zalph1 = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) ) 400 z1_alph1 = 1._wp / ( zalph1 + 1._wp ) 401 zalph2 = zalph1 402 z1_alph2 = z1_alph1 403 ! explicit: 404 ! z1_alph1 = 1._wp / zalph1 405 ! z1_alph2 = 1._wp / zalph1 406 ! zalph1 = zalph1 - 1._wp 407 ! zalph2 = zalph1 408 ENDIF 409 410 ! stress at T points (zkt/=0 if landfast) 411 zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv * (1._wp + zkt) - zdelta * (1._wp - zkt) ) ) * z1_alph1 412 zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2 413 414 END_2D 433 415 CALL lbc_lnk( 'icedyn_rhg_evp', zp_delt, 'T', 1. ) 434 416 435 417 ! Save beta at T-points for further computations 436 418 IF( ln_aEVP ) THEN 437 DO jj = 1, jpj 438 DO ji = 1, jpi 439 zbeta(ji,jj) = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) ) 440 END DO 441 END DO 419 DO_2D_11_11 420 zbeta(ji,jj) = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) ) 421 END_2D 442 422 ENDIF 443 423 444 DO jj = 1, jpjm1 445 DO ji = 1, jpim1 446 447 ! alpha for aEVP 448 IF( ln_aEVP ) THEN 449 zalph2 = MAX( zbeta(ji,jj), zbeta(ji+1,jj), zbeta(ji,jj+1), zbeta(ji+1,jj+1) ) 450 z1_alph2 = 1._wp / ( zalph2 + 1._wp ) 451 ! explicit: 452 ! z1_alph2 = 1._wp / zalph2 453 ! zalph2 = zalph2 - 1._wp 454 ENDIF 455 456 ! P/delta at F points 457 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) ) 458 459 ! stress at F points (zkt/=0 if landfast) 460 zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 * (1._wp + zkt) ) * 0.5_wp ) * z1_alph2 461 462 END DO 463 END DO 424 DO_2D_10_10 425 426 ! alpha for aEVP 427 IF( ln_aEVP ) THEN 428 zalph2 = MAX( zbeta(ji,jj), zbeta(ji+1,jj), zbeta(ji,jj+1), zbeta(ji+1,jj+1) ) 429 z1_alph2 = 1._wp / ( zalph2 + 1._wp ) 430 ! explicit: 431 ! z1_alph2 = 1._wp / zalph2 432 ! zalph2 = zalph2 - 1._wp 433 ENDIF 434 435 ! P/delta at F points 436 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) ) 437 438 ! stress at F points (zkt/=0 if landfast) 439 zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 * (1._wp + zkt) ) * 0.5_wp ) * z1_alph2 440 441 END_2D 464 442 465 443 ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- ! 466 DO jj = 2, jpjm1 467 DO ji = fs_2, fs_jpim1 468 ! !--- U points 469 zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj) & 470 & + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj) & 471 & ) * r1_e2u(ji,jj) & 472 & + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) & 473 & ) * 2._wp * r1_e1u(ji,jj) & 474 & ) * r1_e1e2u(ji,jj) 475 ! 476 ! !--- V points 477 zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj) & 478 & - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj) & 479 & ) * r1_e1v(ji,jj) & 480 & + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) & 481 & ) * 2._wp * r1_e2v(ji,jj) & 482 & ) * r1_e1e2v(ji,jj) 483 ! 484 ! !--- ice currents at U-V point 485 v_iceU(ji,jj) = 0.25_wp * ( v_ice(ji,jj) + v_ice(ji,jj-1) + v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * umask(ji,jj,1) 486 u_iceV(ji,jj) = 0.25_wp * ( u_ice(ji,jj) + u_ice(ji-1,jj) + u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * vmask(ji,jj,1) 487 ! 488 END DO 489 END DO 444 DO_2D_00_00 445 ! !--- U points 446 zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj) & 447 & + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj) & 448 & ) * r1_e2u(ji,jj) & 449 & + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) & 450 & ) * 2._wp * r1_e1u(ji,jj) & 451 & ) * r1_e1e2u(ji,jj) 452 ! 453 ! !--- V points 454 zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj) & 455 & - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj) & 456 & ) * r1_e1v(ji,jj) & 457 & + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) & 458 & ) * 2._wp * r1_e2v(ji,jj) & 459 & ) * r1_e1e2v(ji,jj) 460 ! 461 ! !--- ice currents at U-V point 462 v_iceU(ji,jj) = 0.25_wp * ( v_ice(ji,jj) + v_ice(ji,jj-1) + v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * umask(ji,jj,1) 463 u_iceV(ji,jj) = 0.25_wp * ( u_ice(ji,jj) + u_ice(ji-1,jj) + u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * vmask(ji,jj,1) 464 ! 465 END_2D 490 466 ! 491 467 ! --- Computation of ice velocity --- ! … … 494 470 IF( MOD(jter,2) == 0 ) THEN ! even iterations 495 471 ! 496 DO jj = 2, jpjm1 497 DO ji = fs_2, fs_jpim1 498 ! !--- tau_io/(v_oce - v_ice) 499 zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) ) & 500 & + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) ) 501 ! !--- Ocean-to-Ice stress 502 ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) ) 503 ! 504 ! !--- tau_bottom/v_ice 505 zvel = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) 506 zTauB = ztauy_base(ji,jj) / zvel 507 ! !--- OceanBottom-to-Ice stress 508 ztauy_bi(ji,jj) = zTauB * v_ice(ji,jj) 509 ! 510 ! !--- Coriolis at V-points (energy conserving formulation) 511 zCorV(ji,jj) = - 0.25_wp * r1_e2v(ji,jj) * & 512 & ( zmf(ji,jj ) * ( e2u(ji,jj ) * u_ice(ji,jj ) + e2u(ji-1,jj ) * u_ice(ji-1,jj ) ) & 513 & + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) ) 514 ! 515 ! !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 516 zRHS = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) 517 ! 518 ! !--- landfast switch => 0 = static friction : TauB > RHS & sign(TauB) /= sign(RHS) 519 ! 1 = sliding friction : TauB < RHS 520 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztauy_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) ) 521 ! 522 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 523 zbetav = MAX( zbeta(ji,jj), zbeta(ji,jj+1) ) 524 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbetav * v_ice(ji,jj) + v_ice_b(ji,jj) ) & ! previous velocity 525 & + zRHS + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 526 & ) / MAX( zepsi, zmV_t(ji,jj) * ( zbetav + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 527 & + ( 1._wp - rswitch ) * ( v_ice_b(ji,jj) & 528 & + v_ice (ji,jj) * MAX( 0._wp, zbetav - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 529 & ) / ( zbetav + 1._wp ) & 530 & ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 531 & ) * zmsk00y(ji,jj) 532 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 533 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) & ! previous velocity 534 & + zRHS + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 535 & ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 536 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 537 & ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 538 & ) * zmsk00y(ji,jj) 539 ENDIF 540 END DO 541 END DO 472 DO_2D_00_00 473 ! !--- tau_io/(v_oce - v_ice) 474 zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) ) & 475 & + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) ) 476 ! !--- Ocean-to-Ice stress 477 ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) ) 478 ! 479 ! !--- tau_bottom/v_ice 480 zvel = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) 481 zTauB = ztauy_base(ji,jj) / zvel 482 ! !--- OceanBottom-to-Ice stress 483 ztauy_bi(ji,jj) = zTauB * v_ice(ji,jj) 484 ! 485 ! !--- Coriolis at V-points (energy conserving formulation) 486 zCorV(ji,jj) = - 0.25_wp * r1_e2v(ji,jj) * & 487 & ( zmf(ji,jj ) * ( e2u(ji,jj ) * u_ice(ji,jj ) + e2u(ji-1,jj ) * u_ice(ji-1,jj ) ) & 488 & + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) ) 489 ! 490 ! !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 491 zRHS = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) 492 ! 493 ! !--- landfast switch => 0 = static friction : TauB > RHS & sign(TauB) /= sign(RHS) 494 ! 1 = sliding friction : TauB < RHS 495 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztauy_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) ) 496 ! 497 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 498 zbetav = MAX( zbeta(ji,jj), zbeta(ji,jj+1) ) 499 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbetav * v_ice(ji,jj) + v_ice_b(ji,jj) ) & ! previous velocity 500 & + zRHS + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 501 & ) / MAX( zepsi, zmV_t(ji,jj) * ( zbetav + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 502 & + ( 1._wp - rswitch ) * ( v_ice_b(ji,jj) & 503 & + v_ice (ji,jj) * MAX( 0._wp, zbetav - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 504 & ) / ( zbetav + 1._wp ) & 505 & ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 506 & ) * zmsk00y(ji,jj) 507 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 508 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) & ! previous velocity 509 & + zRHS + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 510 & ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 511 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 512 & ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 513 & ) * zmsk00y(ji,jj) 514 ENDIF 515 END_2D 542 516 CALL lbc_lnk( 'icedyn_rhg_evp', v_ice, 'V', -1. ) 543 517 ! … … 548 522 IF( ln_bdy ) CALL bdy_ice_dyn( 'V' ) 549 523 ! 550 DO jj = 2, jpjm1 551 DO ji = fs_2, fs_jpim1 552 ! !--- tau_io/(u_oce - u_ice) 553 zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) ) & 554 & + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) ) 555 ! !--- Ocean-to-Ice stress 556 ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) ) 557 ! 558 ! !--- tau_bottom/u_ice 559 zvel = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) 560 zTauB = ztaux_base(ji,jj) / zvel 561 ! !--- OceanBottom-to-Ice stress 562 ztaux_bi(ji,jj) = zTauB * u_ice(ji,jj) 563 ! 564 ! !--- Coriolis at U-points (energy conserving formulation) 565 zCorU(ji,jj) = 0.25_wp * r1_e1u(ji,jj) * & 566 & ( zmf(ji ,jj) * ( e1v(ji ,jj) * v_ice(ji ,jj) + e1v(ji ,jj-1) * v_ice(ji ,jj-1) ) & 567 & + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) ) 568 ! 569 ! !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 570 zRHS = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj) 571 ! 572 ! !--- landfast switch => 0 = static friction : TauB > RHS & sign(TauB) /= sign(RHS) 573 ! 1 = sliding friction : TauB < RHS 574 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztaux_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) ) 575 ! 576 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 577 zbetau = MAX( zbeta(ji,jj), zbeta(ji+1,jj) ) 578 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbetau * u_ice(ji,jj) + u_ice_b(ji,jj) ) & ! previous velocity 579 & + zRHS + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 580 & ) / MAX( zepsi, zmU_t(ji,jj) * ( zbetau + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 581 & + ( 1._wp - rswitch ) * ( u_ice_b(ji,jj) & 582 & + u_ice (ji,jj) * MAX( 0._wp, zbetau - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 583 & ) / ( zbetau + 1._wp ) & 584 & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 585 & ) * zmsk00x(ji,jj) 586 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 587 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) & ! previous velocity 588 & + zRHS + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 589 & ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 590 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 591 & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 592 & ) * zmsk00x(ji,jj) 593 ENDIF 594 END DO 595 END DO 524 DO_2D_00_00 525 ! !--- tau_io/(u_oce - u_ice) 526 zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) ) & 527 & + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) ) 528 ! !--- Ocean-to-Ice stress 529 ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) ) 530 ! 531 ! !--- tau_bottom/u_ice 532 zvel = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) 533 zTauB = ztaux_base(ji,jj) / zvel 534 ! !--- OceanBottom-to-Ice stress 535 ztaux_bi(ji,jj) = zTauB * u_ice(ji,jj) 536 ! 537 ! !--- Coriolis at U-points (energy conserving formulation) 538 zCorU(ji,jj) = 0.25_wp * r1_e1u(ji,jj) * & 539 & ( zmf(ji ,jj) * ( e1v(ji ,jj) * v_ice(ji ,jj) + e1v(ji ,jj-1) * v_ice(ji ,jj-1) ) & 540 & + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) ) 541 ! 542 ! !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 543 zRHS = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj) 544 ! 545 ! !--- landfast switch => 0 = static friction : TauB > RHS & sign(TauB) /= sign(RHS) 546 ! 1 = sliding friction : TauB < RHS 547 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztaux_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) ) 548 ! 549 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 550 zbetau = MAX( zbeta(ji,jj), zbeta(ji+1,jj) ) 551 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbetau * u_ice(ji,jj) + u_ice_b(ji,jj) ) & ! previous velocity 552 & + zRHS + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 553 & ) / MAX( zepsi, zmU_t(ji,jj) * ( zbetau + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 554 & + ( 1._wp - rswitch ) * ( u_ice_b(ji,jj) & 555 & + u_ice (ji,jj) * MAX( 0._wp, zbetau - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 556 & ) / ( zbetau + 1._wp ) & 557 & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 558 & ) * zmsk00x(ji,jj) 559 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 560 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) & ! previous velocity 561 & + zRHS + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 562 & ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 563 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 564 & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 565 & ) * zmsk00x(ji,jj) 566 ENDIF 567 END_2D 596 568 CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1. ) 597 569 ! … … 604 576 ELSE ! odd iterations 605 577 ! 606 DO jj = 2, jpjm1 607 DO ji = fs_2, fs_jpim1 608 ! !--- tau_io/(u_oce - u_ice) 609 zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) ) & 610 & + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) ) 611 ! !--- Ocean-to-Ice stress 612 ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) ) 613 ! 614 ! !--- tau_bottom/u_ice 615 zvel = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) 616 zTauB = ztaux_base(ji,jj) / zvel 617 ! !--- OceanBottom-to-Ice stress 618 ztaux_bi(ji,jj) = zTauB * u_ice(ji,jj) 619 ! 620 ! !--- Coriolis at U-points (energy conserving formulation) 621 zCorU(ji,jj) = 0.25_wp * r1_e1u(ji,jj) * & 622 & ( zmf(ji ,jj) * ( e1v(ji ,jj) * v_ice(ji ,jj) + e1v(ji ,jj-1) * v_ice(ji ,jj-1) ) & 623 & + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) ) 624 ! 625 ! !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 626 zRHS = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj) 627 ! 628 ! !--- landfast switch => 0 = static friction : TauB > RHS & sign(TauB) /= sign(RHS) 629 ! 1 = sliding friction : TauB < RHS 630 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztaux_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) ) 631 ! 632 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 633 zbetau = MAX( zbeta(ji,jj), zbeta(ji+1,jj) ) 634 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbetau * u_ice(ji,jj) + u_ice_b(ji,jj) ) & ! previous velocity 635 & + zRHS + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 636 & ) / MAX( zepsi, zmU_t(ji,jj) * ( zbetau + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 637 & + ( 1._wp - rswitch ) * ( u_ice_b(ji,jj) & 638 & + u_ice (ji,jj) * MAX( 0._wp, zbetau - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 639 & ) / ( zbetau + 1._wp ) & 640 & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 641 & ) * zmsk00x(ji,jj) 642 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 643 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) & ! previous velocity 644 & + zRHS + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 645 & ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 646 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 647 & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 648 & ) * zmsk00x(ji,jj) 649 ENDIF 650 END DO 651 END DO 578 DO_2D_00_00 579 ! !--- tau_io/(u_oce - u_ice) 580 zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) ) & 581 & + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) ) 582 ! !--- Ocean-to-Ice stress 583 ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) ) 584 ! 585 ! !--- tau_bottom/u_ice 586 zvel = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) 587 zTauB = ztaux_base(ji,jj) / zvel 588 ! !--- OceanBottom-to-Ice stress 589 ztaux_bi(ji,jj) = zTauB * u_ice(ji,jj) 590 ! 591 ! !--- Coriolis at U-points (energy conserving formulation) 592 zCorU(ji,jj) = 0.25_wp * r1_e1u(ji,jj) * & 593 & ( zmf(ji ,jj) * ( e1v(ji ,jj) * v_ice(ji ,jj) + e1v(ji ,jj-1) * v_ice(ji ,jj-1) ) & 594 & + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) ) 595 ! 596 ! !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 597 zRHS = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj) 598 ! 599 ! !--- landfast switch => 0 = static friction : TauB > RHS & sign(TauB) /= sign(RHS) 600 ! 1 = sliding friction : TauB < RHS 601 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztaux_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) ) 602 ! 603 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 604 zbetau = MAX( zbeta(ji,jj), zbeta(ji+1,jj) ) 605 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * ( zbetau * u_ice(ji,jj) + u_ice_b(ji,jj) ) & ! previous velocity 606 & + zRHS + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 607 & ) / MAX( zepsi, zmU_t(ji,jj) * ( zbetau + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 608 & + ( 1._wp - rswitch ) * ( u_ice_b(ji,jj) & 609 & + u_ice (ji,jj) * MAX( 0._wp, zbetau - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 610 & ) / ( zbetau + 1._wp ) & 611 & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 612 & ) * zmsk00x(ji,jj) 613 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 614 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) & ! previous velocity 615 & + zRHS + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 616 & ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 617 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 618 & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 619 & ) * zmsk00x(ji,jj) 620 ENDIF 621 END_2D 652 622 CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1. ) 653 623 ! … … 658 628 IF( ln_bdy ) CALL bdy_ice_dyn( 'U' ) 659 629 ! 660 DO jj = 2, jpjm1 661 DO ji = fs_2, fs_jpim1 662 ! !--- tau_io/(v_oce - v_ice) 663 zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) ) & 664 & + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) ) 665 ! !--- Ocean-to-Ice stress 666 ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) ) 667 ! 668 ! !--- tau_bottom/v_ice 669 zvel = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) 670 zTauB = ztauy_base(ji,jj) / zvel 671 ! !--- OceanBottom-to-Ice stress 672 ztauy_bi(ji,jj) = zTauB * v_ice(ji,jj) 673 ! 674 ! !--- Coriolis at v-points (energy conserving formulation) 675 zCorV(ji,jj) = - 0.25_wp * r1_e2v(ji,jj) * & 676 & ( zmf(ji,jj ) * ( e2u(ji,jj ) * u_ice(ji,jj ) + e2u(ji-1,jj ) * u_ice(ji-1,jj ) ) & 677 & + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) ) 678 ! 679 ! !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 680 zRHS = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) 681 ! 682 ! !--- landfast switch => 0 = static friction : TauB > RHS & sign(TauB) /= sign(RHS) 683 ! 1 = sliding friction : TauB < RHS 684 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztauy_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) ) 685 ! 686 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 687 zbetav = MAX( zbeta(ji,jj), zbeta(ji,jj+1) ) 688 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbetav * v_ice(ji,jj) + v_ice_b(ji,jj) ) & ! previous velocity 689 & + zRHS + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 690 & ) / MAX( zepsi, zmV_t(ji,jj) * ( zbetav + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 691 & + ( 1._wp - rswitch ) * ( v_ice_b(ji,jj) & 692 & + v_ice (ji,jj) * MAX( 0._wp, zbetav - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 693 & ) / ( zbetav + 1._wp ) & 694 & ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 695 & ) * zmsk00y(ji,jj) 696 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 697 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) & ! previous velocity 698 & + zRHS + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 699 & ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 700 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 701 & ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 702 & ) * zmsk00y(ji,jj) 703 ENDIF 704 END DO 705 END DO 630 DO_2D_00_00 631 ! !--- tau_io/(v_oce - v_ice) 632 zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) ) & 633 & + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) ) 634 ! !--- Ocean-to-Ice stress 635 ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) ) 636 ! 637 ! !--- tau_bottom/v_ice 638 zvel = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) 639 zTauB = ztauy_base(ji,jj) / zvel 640 ! !--- OceanBottom-to-Ice stress 641 ztauy_bi(ji,jj) = zTauB * v_ice(ji,jj) 642 ! 643 ! !--- Coriolis at v-points (energy conserving formulation) 644 zCorV(ji,jj) = - 0.25_wp * r1_e2v(ji,jj) * & 645 & ( zmf(ji,jj ) * ( e2u(ji,jj ) * u_ice(ji,jj ) + e2u(ji-1,jj ) * u_ice(ji-1,jj ) ) & 646 & + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) ) 647 ! 648 ! !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 649 zRHS = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) 650 ! 651 ! !--- landfast switch => 0 = static friction : TauB > RHS & sign(TauB) /= sign(RHS) 652 ! 1 = sliding friction : TauB < RHS 653 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztauy_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) ) 654 ! 655 IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 656 zbetav = MAX( zbeta(ji,jj), zbeta(ji,jj+1) ) 657 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbetav * v_ice(ji,jj) + v_ice_b(ji,jj) ) & ! previous velocity 658 & + zRHS + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 659 & ) / MAX( zepsi, zmV_t(ji,jj) * ( zbetav + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 660 & + ( 1._wp - rswitch ) * ( v_ice_b(ji,jj) & 661 & + v_ice (ji,jj) * MAX( 0._wp, zbetav - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 662 & ) / ( zbetav + 1._wp ) & 663 & ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 664 & ) * zmsk00y(ji,jj) 665 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 666 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) & ! previous velocity 667 & + zRHS + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 668 & ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 669 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 670 & ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 671 & ) * zmsk00y(ji,jj) 672 ENDIF 673 END_2D 706 674 CALL lbc_lnk( 'icedyn_rhg_evp', v_ice, 'V', -1. ) 707 675 ! … … 725 693 ! 4) Recompute delta, shear and div (inputs for mechanical redistribution) 726 694 !------------------------------------------------------------------------------! 727 DO jj = 1, jpjm1 728 DO ji = 1, jpim1 729 730 ! shear at F points 731 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) & 732 & + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 733 & ) * r1_e1e2f(ji,jj) * zfmask(ji,jj) 734 735 END DO 736 END DO 695 DO_2D_10_10 696 697 ! shear at F points 698 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) & 699 & + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 700 & ) * r1_e1e2f(ji,jj) * zfmask(ji,jj) 701 702 END_2D 737 703 738 DO jj = 2, jpjm1 739 DO ji = 2, jpim1 ! no vector loop 740 741 ! tension**2 at T points 742 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) & 743 & - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) & 744 & ) * r1_e1e2t(ji,jj) 745 zdt2 = zdt * zdt 746 747 ! shear**2 at T points (doc eq. A16) 748 zds2 = ( zds(ji,jj ) * zds(ji,jj ) * e1e2f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e1e2f(ji-1,jj ) & 749 & + zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1) & 750 & ) * 0.25_wp * r1_e1e2t(ji,jj) 751 752 ! shear at T points 753 pshear_i(ji,jj) = SQRT( zdt2 + zds2 ) 754 755 ! divergence at T points 756 pdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) & 757 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) & 758 & ) * r1_e1e2t(ji,jj) 759 760 ! delta at T points 761 zdelta = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) 762 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0 763 pdelta_i(ji,jj) = zdelta + rn_creepl * rswitch 764 765 END DO 766 END DO 704 DO_2D_00_00 705 706 ! tension**2 at T points 707 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) & 708 & - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) & 709 & ) * r1_e1e2t(ji,jj) 710 zdt2 = zdt * zdt 711 712 ! shear**2 at T points (doc eq. A16) 713 zds2 = ( zds(ji,jj ) * zds(ji,jj ) * e1e2f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e1e2f(ji-1,jj ) & 714 & + zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1) & 715 & ) * 0.25_wp * r1_e1e2t(ji,jj) 716 717 ! shear at T points 718 pshear_i(ji,jj) = SQRT( zdt2 + zds2 ) 719 720 ! divergence at T points 721 pdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) & 722 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) & 723 & ) * r1_e1e2t(ji,jj) 724 725 ! delta at T points 726 zdelta = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) 727 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0 728 pdelta_i(ji,jj) = zdelta + rn_creepl * rswitch 729 730 END_2D 767 731 CALL lbc_lnk_multi( 'icedyn_rhg_evp', pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1. ) 768 732 … … 802 766 ALLOCATE( zsig1(jpi,jpj) , zsig2(jpi,jpj) , zsig3(jpi,jpj) ) 803 767 ! 804 DO jj = 2, jpjm1 805 DO ji = 2, jpim1 806 zdum1 = ( zmsk00(ji-1,jj) * pstress12_i(ji-1,jj) + zmsk00(ji ,jj-1) * pstress12_i(ji ,jj-1) + & ! stress12_i at T-point 807 & zmsk00(ji ,jj) * pstress12_i(ji ,jj) + zmsk00(ji-1,jj-1) * pstress12_i(ji-1,jj-1) ) & 808 & / MAX( 1._wp, zmsk00(ji-1,jj) + zmsk00(ji,jj-1) + zmsk00(ji,jj) + zmsk00(ji-1,jj-1) ) 809 810 zshear = SQRT( pstress2_i(ji,jj) * pstress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress 811 812 zdum2 = zmsk00(ji,jj) / MAX( 1._wp, strength(ji,jj) ) 768 DO_2D_00_00 769 zdum1 = ( zmsk00(ji-1,jj) * pstress12_i(ji-1,jj) + zmsk00(ji ,jj-1) * pstress12_i(ji ,jj-1) + & ! stress12_i at T-point 770 & zmsk00(ji ,jj) * pstress12_i(ji ,jj) + zmsk00(ji-1,jj-1) * pstress12_i(ji-1,jj-1) ) & 771 & / MAX( 1._wp, zmsk00(ji-1,jj) + zmsk00(ji,jj-1) + zmsk00(ji,jj) + zmsk00(ji-1,jj-1) ) 772 773 zshear = SQRT( pstress2_i(ji,jj) * pstress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress 774 775 zdum2 = zmsk00(ji,jj) / MAX( 1._wp, strength(ji,jj) ) 813 776 814 777 !! zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002) … … 816 779 !! zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) ! quadratic relation linking compressive stress to shear stress 817 780 !! ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11)) 818 zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) ) ! compressive stress, see Bouillon et al. 2015 819 zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear ) ! shear stress 820 zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) 821 END DO 822 END DO 781 zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) ) ! compressive stress, see Bouillon et al. 2015 782 zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear ) ! shear stress 783 zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) 784 END_2D 823 785 CALL lbc_lnk_multi( 'icedyn_rhg_evp', zsig1, 'T', 1., zsig2, 'T', 1., zsig3, 'T', 1. ) 824 786 ! … … 855 817 & zdiag_xmtrp_snw(jpi,jpj) , zdiag_ymtrp_snw(jpi,jpj) , zdiag_xatrp(jpi,jpj) , zdiag_yatrp(jpi,jpj) ) 856 818 ! 857 DO jj = 2, jpjm1 858 DO ji = 2, jpim1 859 ! 2D ice mass, snow mass, area transport arrays (X, Y) 860 zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * zmsk00(ji,jj) 861 zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * zmsk00(ji,jj) 862 863 zdiag_xmtrp_ice(ji,jj) = rhoi * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component 864 zdiag_ymtrp_ice(ji,jj) = rhoi * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) ! '' Y- '' 865 866 zdiag_xmtrp_snw(ji,jj) = rhos * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component 867 zdiag_ymtrp_snw(ji,jj) = rhos * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) ! '' Y- '' 868 869 zdiag_xatrp(ji,jj) = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) ) ! area transport, X-component 870 zdiag_yatrp(ji,jj) = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) ) ! '' Y- '' 871 872 END DO 873 END DO 819 DO_2D_00_00 820 ! 2D ice mass, snow mass, area transport arrays (X, Y) 821 zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * zmsk00(ji,jj) 822 zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * zmsk00(ji,jj) 823 824 zdiag_xmtrp_ice(ji,jj) = rhoi * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component 825 zdiag_ymtrp_ice(ji,jj) = rhoi * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) ! '' Y- '' 826 827 zdiag_xmtrp_snw(ji,jj) = rhos * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component 828 zdiag_ymtrp_snw(ji,jj) = rhos * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) ! '' Y- '' 829 830 zdiag_xatrp(ji,jj) = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) ) ! area transport, X-component 831 zdiag_yatrp(ji,jj) = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) ) ! '' Y- '' 832 833 END_2D 874 834 875 835 CALL lbc_lnk_multi( 'icedyn_rhg_evp', zdiag_xmtrp_ice, 'U', -1., zdiag_ymtrp_ice, 'V', -1., & … … 957 917 zresm = 0._wp 958 918 ELSE 959 DO jj = 1, jpj 960 DO ji = 1, jpi 961 zres(ji,jj) = MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), & 962 & ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * zmsk15(ji,jj) 963 END DO 964 END DO 919 DO_2D_11_11 920 zres(ji,jj) = MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), & 921 & ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * zmsk15(ji,jj) 922 END_2D 965 923 zresm = MAXVAL( zres ) 966 924 CALL mpp_max( 'icedyn_rhg_evp', zresm ) ! max over the global domain
Note: See TracChangeset
for help on using the changeset viewer.
