- Timestamp:
- 2017-09-01T15:49:35+02:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg_evp.F90
r8409 r8486 45 45 # include "vectopt_loop_substitute.h90" 46 46 !!---------------------------------------------------------------------- 47 !! NEMO/ LIM3 4.0 , UCL - NEMO Consortium (2011)47 !! NEMO/ICE 4.0 , NEMO Consortium (2017) 48 48 !! $Id: icerhg_evp.F90 8378 2017-07-26 13:55:59Z clem $ 49 49 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 386 386 DO jj = 2, jpjm1 387 387 DO ji = fs_2, fs_jpim1 388 389 ! U points 388 ! !--- U points 390 389 zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj) & 391 390 & + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj) & … … 394 393 & ) * 2._wp * r1_e1u(ji,jj) & 395 394 & ) * r1_e1e2u(ji,jj) 396 397 !V points395 ! 396 ! !--- V points 398 397 zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj) & 399 398 & - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj) & … … 402 401 & ) * 2._wp * r1_e2v(ji,jj) & 403 402 & ) * r1_e1e2v(ji,jj) 404 405 !u_ice at V point403 ! 404 ! !--- u_ice at V point 406 405 u_iceV(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj+1) & 407 406 & + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1) 408 409 !v_ice at U point407 ! 408 ! !--- v_ice at U point 410 409 v_iceU(ji,jj) = 0.5_wp * ( ( v_ice(ji ,jj) + v_ice(ji ,jj-1) ) * e1t(ji+1,jj) & 411 410 & + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1) 412 411 ! 413 412 END DO 414 413 END DO … … 417 416 ! Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta are chosen wisely and large nn_nevp 418 417 ! Bouillon et al. 2009 (eq 34-35) => stable 419 IF( MOD(jter,2) .EQ.0 ) THEN ! even iterations420 418 IF( MOD(jter,2) == 0 ) THEN ! even iterations 419 ! 421 420 DO jj = 2, jpjm1 422 421 DO ji = fs_2, fs_jpim1 423 424 ! tau_io/(v_oce - v_ice) 422 ! !--- tau_io/(v_oce - v_ice) 425 423 zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) ) & 426 424 & + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) ) 427 428 ! Ocean-to-Ice stress 425 ! !--- Ocean-to-Ice stress 429 426 ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) ) 430 431 !tau_bottom/v_ice427 ! 428 ! !--- tau_bottom/v_ice 432 429 zvel = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) ) 433 430 zTauB = - tau_icebfr(ji,jj) / zvel 434 435 !Coriolis at V-points (energy conserving formulation)431 ! 432 ! !--- Coriolis at V-points (energy conserving formulation) 436 433 zCory(ji,jj) = - 0.25_wp * r1_e2v(ji,jj) * & 437 434 & ( zmf(ji,jj ) * ( e2u(ji,jj ) * u_ice(ji,jj ) + e2u(ji-1,jj ) * u_ice(ji-1,jj ) ) & 438 435 & + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) ) 439 440 !Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io436 ! 437 ! !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 441 438 zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) 442 443 !landfast switch => 0 = static friction ; 1 = sliding friction439 ! 440 ! !--- landfast switch => 0 = static friction ; 1 = sliding friction 444 441 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 445 446 !ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)442 ! 443 ! !--- ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 447 444 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) & ! previous velocity 448 445 & + zTauE + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) … … 451 448 & ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce if mass < zmmin 452 449 & ) * zmaskV(ji,jj) 450 ! 453 451 ! Bouillon 2013 454 452 !!v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) ) & 455 453 !! & + zfV(ji,jj) + zCory(ji,jj) + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj) & 456 454 !! & ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj) 457 455 ! 458 456 END DO 459 457 END DO 460 458 CALL lbc_lnk( v_ice, 'V', -1. ) 461 459 ! 462 460 #if defined key_agrif 463 461 !! CALL agrif_interp_lim3( 'V', jter, nn_nevp ) … … 465 463 #endif 466 464 IF( ln_bdy ) CALL bdy_ice_dyn( 'V' ) 467 465 ! 468 466 DO jj = 2, jpjm1 469 467 DO ji = fs_2, fs_jpim1 … … 505 503 END DO 506 504 CALL lbc_lnk( u_ice, 'U', -1. ) 507 505 ! 508 506 #if defined key_agrif 509 507 !! CALL agrif_interp_lim3( 'U', jter, nn_nevp ) … … 511 509 #endif 512 510 IF( ln_bdy ) CALL bdy_ice_dyn( 'U' ) 513 511 ! 514 512 ELSE ! odd iterations 515 513 ! 516 514 DO jj = 2, jpjm1 517 515 DO ji = fs_2, fs_jpim1 … … 553 551 END DO 554 552 CALL lbc_lnk( u_ice, 'U', -1. ) 555 553 ! 556 554 #if defined key_agrif 557 555 !! CALL agrif_interp_lim3( 'U', jter, nn_nevp ) … … 559 557 #endif 560 558 IF( ln_bdy ) CALL bdy_ice_dyn( 'U' ) 561 559 ! 562 560 DO jj = 2, jpjm1 563 561 DO ji = fs_2, fs_jpim1 … … 599 597 END DO 600 598 CALL lbc_lnk( v_ice, 'V', -1. ) 601 599 ! 602 600 #if defined key_agrif 603 601 !! CALL agrif_interp_lim3( 'V', jter, nn_nevp ) … … 605 603 #endif 606 604 IF( ln_bdy ) CALL bdy_ice_dyn( 'V' ) 607 605 ! 608 606 ENDIF 609 607 … … 675 673 ! 5) SIMIP diagnostics 676 674 !------------------------------------------------------------------------------! 677 675 676 !!gm encapsulate with a flag (iom_use of the variable or better a flag defined one for all testing if one of the 677 !! diag is output. NB the diag_... are should only be ALLOCATED if the flag is true ! 678 678 679 DO jj = 2, jpjm1 679 680 DO ji = 2, jpim1 … … 714 715 715 716 CALL lbc_lnk_multi( diag_sig1 , 'T', 1., diag_sig2 , 'T', 1., & 716 &diag_dssh_dx, 'U', -1., diag_dssh_dy, 'V', -1., &717 &diag_corstrx, 'U', -1., diag_corstry, 'V', -1., &718 &diag_intstrx, 'U', -1., diag_intstry, 'V', -1. )717 & diag_dssh_dx, 'U', -1., diag_dssh_dy, 'V', -1., & 718 & diag_corstrx, 'U', -1., diag_corstry, 'V', -1., & 719 & diag_intstrx, 'U', -1., diag_intstry, 'V', -1. ) 719 720 720 721 CALL lbc_lnk_multi( diag_utau_oi, 'U', -1., diag_vtau_oi, 'V', -1. ) 721 722 722 CALL lbc_lnk_multi( diag_xmtrp_ice, 'U', -1., diag_xmtrp_snw, 'U', -1., &723 & diag_xatrp , 'U', -1., diag_ymtrp_ice, 'V', -1.,&724 & diag_ymtrp_snw, 'V', -1., diag_yatrp , 'V', -1.)723 CALL lbc_lnk_multi( diag_xmtrp_ice, 'U', -1., diag_xmtrp_snw, 'U', -1., & 724 & diag_xatrp , 'U', -1., diag_ymtrp_ice, 'V', -1., & 725 & diag_ymtrp_snw, 'V', -1., diag_yatrp , 'V', -1. ) 725 726 726 727 ! … … 734 735 CALL prt_ctl_info(charout) 735 736 CALL prt_ctl(tab2d_1=u_ice, clinfo1=' ice_rhg_evp : u_ice :', tab2d_2=v_ice, clinfo2=' v_ice :') 736 ENDIF 737 738 ! print charge ellipse 739 ! This can be desactivated once the user is sure that the stress state 740 ! lie on the charge ellipse. See Bouillon et al. 08 for more details 741 IF(ln_ctl) THEN 737 ! 738 ! print charge ellipse 739 ! This can be desactivated once the user is sure that the stress state 740 ! lie on the charge ellipse. See Bouillon et al. (2008) for more details 742 741 IF( MOD(kt_ice+nn_fsbc-1,nwrite) == 0 ) THEN 743 742 WRITE(charout,FMT="('ice_rhg_evp :', I4, I6, I1, I1, A10)") 1000, kt_ice, 0, 0, ' ch. ell. ' … … 760 759 END SUBROUTINE ice_rhg_evp 761 760 761 #else 762 !!---------------------------------------------------------------------- 763 !! Default option Empty module NO LIM-3 sea-ice model 764 !!---------------------------------------------------------------------- 762 765 #endif 763 766
Note: See TracChangeset
for help on using the changeset viewer.
