- Timestamp:
- 2020-09-29T11:22:04+02:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg_evp.F90
r13279 r13536 137 137 ! 138 138 REAL(wp), DIMENSION(jpi,jpj) :: zp_delt ! P/delta at T points 139 REAL(wp), DIMENSION(jpi,jpj) :: zdeltastar_t ! Delta* at T-points 139 140 REAL(wp), DIMENSION(jpi,jpj) :: zbeta ! beta coef from Kimmritz 2017 140 141 ! … … 168 169 REAL(wp), DIMENSION(jpi,jpj) :: zu_ice, zv_ice 169 170 !! --- diags 170 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zsig1, zsig2, zsig3 171 REAL(wp) :: zsig1, zsig2, zsig12 172 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zsig_I, zsig_II, zsig1_p, zsig2_p, zten_i 171 173 !! --- SIMIP diags 172 174 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s) … … 406 408 ! delta at T points 407 409 zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 ) 410 411 zdeltastar_t(ji,jj) = zdelta + rn_creepl ! store delta* at previous iterate (for ellipse computation) 408 412 409 413 ! P/delta at T points … … 695 699 & ) * zmsk00y(ji,jj) 696 700 ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 701 !---> MV : comment above is misleading as EVP is a purely explicit method 702 !---> MV : I found code below quite unreadable :-) 703 !---> MV : two maskings + landfast / no landfast options, I think this is a bit too much 704 !---> MV : I would suggest to split at least the landfast / masking steps if possible 705 !---> MV : Should also comment why 1% of ocean velocity is assumed - not sure it makes sense btw 706 !---> MV : I would say 100% of ocean current makes much more sense for free-drift 707 !---> MV : of very low concentration sea ice.. 708 !---> MV : but maybe it is a numerical trick... in brief it's n 697 709 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) & ! previous velocity 698 710 & + zRHS + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) … … 715 727 716 728 ! convergence test 717 IF( ln_rhg_chkcvg ) CALL rhg_cvg ( kt, jter, nn_nevp, u_ice, v_ice, zu_ice, zv_ice )729 IF( ln_rhg_chkcvg ) CALL rhg_cvg_evp( kt, jter, nn_nevp, u_ice, v_ice, zu_ice, zv_ice ) 718 730 ! 719 731 ! ! ==================== ! … … 745 757 zdt2 = zdt * zdt 746 758 759 zten_i(ji,jj) = zdt 760 747 761 ! shear**2 at T points (doc eq. A16) 748 762 zds2 = ( zds(ji,jj ) * zds(ji,jj ) * e1e2f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e1e2f(ji-1,jj ) & … … 797 811 IF( iom_use('icestr') ) CALL iom_put( 'icestr' , strength * zmsk00 ) ! strength 798 812 799 ! --- stress tensor--- !800 IF( iom_use(' isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') .OR. iom_use('normstr') .OR. iom_use('sheastr') ) THEN801 ! 802 ALLOCATE( zsig 1(jpi,jpj) , zsig2(jpi,jpj) , zsig3(jpi,jpj) )813 ! --- Stress tensor invariants (SIMIP diags) --- ! 814 IF( iom_use('normstr') .OR. iom_use('sheastr') ) THEN 815 ! 816 ALLOCATE( zsig_I(jpi,jpj) , zsig_II(jpi,jpj) ) 803 817 ! 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) ) 813 814 !! zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002) 815 !! zsig2(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) - zshear ) ! principal stress (x-direction, see Hunke & Dukowicz 2002) 816 !! 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 !! ! (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 823 CALL lbc_lnk_multi( 'icedyn_rhg_evp', zsig1, 'T', 1., zsig2, 'T', 1., zsig3, 'T', 1. ) 824 ! 825 CALL iom_put( 'isig1' , zsig1 ) 826 CALL iom_put( 'isig2' , zsig2 ) 827 CALL iom_put( 'isig3' , zsig3 ) 828 ! 829 ! Stress tensor invariants (normal and shear stress N/m) 830 IF( iom_use('normstr') ) CALL iom_put( 'normstr' , ( zs1(:,:) + zs2(:,:) ) * zmsk00(:,:) ) ! Normal stress 831 IF( iom_use('sheastr') ) CALL iom_put( 'sheastr' , SQRT( ( zs1(:,:) - zs2(:,:) )**2 + 4*zs12(:,:)**2 ) * zmsk00(:,:) ) ! Shear stress 832 833 DEALLOCATE( zsig1 , zsig2 , zsig3 ) 834 ENDIF 835 818 DO jj = 2, jpj - 1 819 DO ji = 2, jpi - 1 820 821 ! Ice stresses 822 ! sigma1, sigma2, sigma12 are some useful recombination of the stresses (Hunke and Dukowicz MWR 2002, Bouillon et al., OM2013) 823 ! These are NOT stress tensor components, neither stress invariants, neither stress principal components 824 ! I know, this can be confusing... 825 zfac = strength(ji,jj) / ( pdelta_i(ji,jj) + rn_creepl ) 826 zsig1 = zfac * ( pdivu_i(ji,jj) - pdelta_i(ji,jj) ) 827 zsig2 = zfac * z1_ecc2 * zten_i(ji,jj) 828 zsig12 = zfac * z1_ecc2 * pshear_i(ji,jj) 829 830 ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008) 831 zsig_I(ji,jj) = zsig1 * 0.5_wp ! 1st stress invariant, aka average normal stress, aka negative pressure 832 zsig_II(ji,jj) = - SQRT ( zsig2 * zsig2 * 0.25_wp + zsig12 ) ! 2nd '' '', aka maximum shear stress 833 834 END DO 835 END DO 836 837 CALL lbc_lnk_multi( 'icedyn_rhg_vp', zsig_I, 'T', 1., zsig_II, 'T', 1.) 838 839 ! 840 ! Stress tensor invariants (normal and shear stress N/m) - SIMIP diags - definitions following Coon (1974) and Feltham (2008) 841 IF( iom_use('normstr') ) CALL iom_put( 'normstr' , zsig_I(:,:) * zmsk00(:,:) ) ! Normal stress 842 IF( iom_use('sheastr') ) CALL iom_put( 'sheastr' , zsig_II(:,:) * zmsk00(:,:) ) ! Maximum shear stress 843 844 DEALLOCATE ( zsig_I, zsig_II ) 845 846 ENDIF 847 848 ! --- Normalized stress tensor principal components --- ! 849 ! This are used to plot the normalized yield curve, see Lemieux & Dupont, 2020 850 ! Recommendation 1 : we use ice strength, not replacement pressure 851 ! Recommendation 2 : need to use deformations at PREVIOUS iterate for viscosities 852 IF( iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN 853 ! 854 ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) ) 855 ! 856 DO jj = 2, jpj - 1 857 DO ji = 2, jpi - 1 858 859 ! Ice stresses computed with **viscosities** (delta, p/delta) at **previous** iterates 860 ! and **deformations** at current iterates 861 ! following Lemieux & Dupont (2020) 862 zfac = zp_delt(ji,jj) 863 zsig1 = zfac * ( pdivu_i(ji,jj) - zdeltastar_t(ji,jj) ) 864 zsig2 = zfac * z1_ecc2 * zten_i(ji,jj) 865 zsig12 = zfac * z1_ecc2 * pshear_i(ji,jj) 866 867 ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008), T-point 868 zsig_I(ji,jj) = zsig1 * 0.5_wp ! 1st stress invariant, aka average normal stress, aka negative pressure 869 zsig_II(ji,jj) = - SQRT ( zsig2 * zsig2 * 0.25_wp + zsig12 ) ! 2nd '' '', aka maximum shear stress 870 871 ! Normalized principal stresses (used to display the ellipse) 872 z1_strength = 1._wp / strength(ji,jj) 873 zsig1_p(ji,jj) = ( zsig_I(ji,jj) - zsig_II(ji,jj) ) * z1_strength 874 zsig2_p(ji,jj) = ( zsig_II(ji,jj) + zsig_II(ji,jj) ) * z1_strength 875 END DO 876 END DO 877 878 CALL lbc_lnk_multi( 'icedyn_rhg_vp', zsig1_p, 'T', 1., zsig2_p, 'T', 1.) 879 880 ! 881 CALL iom_put( 'sig1_pnorm' , zsig1_p ) 882 CALL iom_put( 'sig2_pnorm' , zsig2_p ) 883 884 DEALLOCATE( zsig1_p , zsig2_p ) 885 886 ENDIF 836 887 ! --- SIMIP --- ! 837 888 IF( iom_use('dssh_dx') .OR. iom_use('dssh_dy') .OR. & … … 907 958 908 959 909 SUBROUTINE rhg_cvg ( kt, kiter, kitermax, pu, pv, pub, pvb )960 SUBROUTINE rhg_cvg_evp( kt, kiter, kitermax, pu, pv, pub, pvb ) 910 961 !!---------------------------------------------------------------------- 911 !! *** ROUTINE rhg_cvg ***962 !! *** ROUTINE rhg_cvg_evp *** 912 963 !! 913 !! ** Purpose : check convergence of oce rheology964 !! ** Purpose : check convergence of evp ice rheology 914 965 !! 915 966 !! ** Method : create a file ice_cvg.nc containing the convergence of ice velocity … … 935 986 IF( lwp ) THEN 936 987 WRITE(numout,*) 937 WRITE(numout,*) 'rhg_cvg : ice rheology convergence control'938 WRITE(numout,*) '~~~~~~~ '988 WRITE(numout,*) 'rhg_cvg_evp : ice rheology convergence control' 989 WRITE(numout,*) '~~~~~~~~~~~' 939 990 ENDIF 940 991 ! … … 974 1025 ENDIF 975 1026 976 END SUBROUTINE rhg_cvg 1027 END SUBROUTINE rhg_cvg_evp 977 1028 978 1029
Note: See TracChangeset
for help on using the changeset viewer.
