Changeset 13646
- Timestamp:
- 2020-10-20T17:33:01+02:00 (4 years ago)
- Location:
- NEMO/releases/r4.0/r4.0-HEAD/src
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icedyn_rhg_evp.F90
r13589 r13646 137 137 ! 138 138 REAL(wp), DIMENSION(jpi,jpj) :: zdelta, zp_delt ! delta and P/delta at T points 139 REAL(wp), DIMENSION(jpi,jpj) :: zten_i ! tension 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, zfac, z1_strength 172 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zsig_I, zsig_II, zsig1_p, zsig2_p 171 173 !! --- SIMIP diags 172 174 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s) … … 764 766 zdt2 = zdt * zdt 765 767 768 zten_i(ji,jj) = zdt 769 766 770 ! shear**2 at T points (doc eq. A16) 767 771 zds2 = ( zds(ji,jj ) * zds(ji,jj ) * e1e2f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e1e2f(ji-1,jj ) & … … 778 782 779 783 ! delta at T points 780 z delta(ji,jj) = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 )781 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -z delta(ji,jj)) ) ! 0 if delta=0782 pdelta_i(ji,jj) = z delta(ji,jj) + rn_creepl * rswitch784 zfac = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) ! delta 785 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zfac ) ) ! 0 if delta=0 786 pdelta_i(ji,jj) = zfac + rn_creepl * rswitch ! delta+creepl 783 787 784 788 END DO 785 789 END DO 786 CALL lbc_lnk_multi( 'icedyn_rhg_evp', pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1. ) 790 CALL lbc_lnk_multi( 'icedyn_rhg_evp', pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1., zten_i, 'T', 1., & 791 & zs1 , 'T', 1., zs2 , 'T', 1., zs12 , 'F', 1. ) 787 792 788 793 ! --- Store the stress tensor for the next time step --- ! 789 CALL lbc_lnk_multi( 'icedyn_rhg_evp', zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. )790 794 pstress1_i (:,:) = zs1 (:,:) 791 795 pstress2_i (:,:) = zs2 (:,:) … … 816 820 IF( iom_use('icestr') ) CALL iom_put( 'icestr' , strength * zmsk00 ) ! strength 817 821 818 ! --- stress tensor--- !819 IF( iom_use(' isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') .OR. iom_use('normstr') .OR. iom_use('sheastr') ) THEN820 ! 821 ALLOCATE( zsig 1(jpi,jpj) , zsig2(jpi,jpj) , zsig3(jpi,jpj) )822 ! --- Stress tensor invariants (SIMIP diags) --- ! 823 IF( iom_use('normstr') .OR. iom_use('sheastr') ) THEN 824 ! 825 ALLOCATE( zsig_I(jpi,jpj) , zsig_II(jpi,jpj) ) 822 826 ! 823 DO jj = 2, jpjm1 824 DO ji = 2, jpim1 825 zdum1 = ( zmsk00(ji-1,jj) * pstress12_i(ji-1,jj) + zmsk00(ji ,jj-1) * pstress12_i(ji ,jj-1) + & ! stress12_i at T-point 826 & zmsk00(ji ,jj) * pstress12_i(ji ,jj) + zmsk00(ji-1,jj-1) * pstress12_i(ji-1,jj-1) ) & 827 & / MAX( 1._wp, zmsk00(ji-1,jj) + zmsk00(ji,jj-1) + zmsk00(ji,jj) + zmsk00(ji-1,jj-1) ) 828 829 zshear = SQRT( pstress2_i(ji,jj) * pstress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress 830 831 zdum2 = zmsk00(ji,jj) / MAX( 1._wp, strength(ji,jj) ) 832 833 !! zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002) 834 !! zsig2(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) - zshear ) ! principal stress (x-direction, see Hunke & Dukowicz 2002) 835 !! 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 836 !! ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11)) 837 zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) ) ! compressive stress, see Bouillon et al. 2015 838 zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear ) ! shear stress 839 zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) 840 END DO 841 END DO 842 CALL lbc_lnk_multi( 'icedyn_rhg_evp', zsig1, 'T', 1., zsig2, 'T', 1., zsig3, 'T', 1. ) 843 ! 844 CALL iom_put( 'isig1' , zsig1 ) 845 CALL iom_put( 'isig2' , zsig2 ) 846 CALL iom_put( 'isig3' , zsig3 ) 847 ! 848 ! Stress tensor invariants (normal and shear stress N/m) 849 IF( iom_use('normstr') ) CALL iom_put( 'normstr' , ( zs1(:,:) + zs2(:,:) ) * zmsk00(:,:) ) ! Normal stress 850 IF( iom_use('sheastr') ) CALL iom_put( 'sheastr' , SQRT( ( zs1(:,:) - zs2(:,:) )**2 + 4*zs12(:,:)**2 ) * zmsk00(:,:) ) ! Shear stress 851 852 DEALLOCATE( zsig1 , zsig2 , zsig3 ) 827 DO jj = 1, jpj 828 DO ji = 1, jpi 829 830 ! Ice stresses 831 ! sigma1, sigma2, sigma12 are some useful recombination of the stresses (Hunke and Dukowicz MWR 2002, Bouillon et al., OM2013) 832 ! These are NOT stress tensor components, neither stress invariants, neither stress principal components 833 ! I know, this can be confusing... 834 zfac = strength(ji,jj) / ( pdelta_i(ji,jj) + rn_creepl ) 835 zsig1 = zfac * ( pdivu_i(ji,jj) - pdelta_i(ji,jj) ) 836 zsig2 = zfac * z1_ecc2 * zten_i(ji,jj) 837 zsig12 = zfac * z1_ecc2 * pshear_i(ji,jj) 838 839 ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008) 840 zsig_I (ji,jj) = zsig1 * 0.5_wp ! 1st stress invariant, aka average normal stress, aka negative pressure 841 zsig_II(ji,jj) = SQRT ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) ) ! 2nd '' '', aka maximum shear stress 842 843 END DO 844 END DO 845 ! 846 ! Stress tensor invariants (normal and shear stress N/m) - SIMIP diags - definitions following Coon (1974) and Feltham (2008) 847 IF( iom_use('normstr') ) CALL iom_put( 'normstr', zsig_I (:,:) * zmsk00(:,:) ) ! Normal stress 848 IF( iom_use('sheastr') ) CALL iom_put( 'sheastr', zsig_II(:,:) * zmsk00(:,:) ) ! Maximum shear stress 849 850 DEALLOCATE ( zsig_I, zsig_II ) 851 853 852 ENDIF 853 854 ! --- Normalized stress tensor principal components --- ! 855 ! This are used to plot the normalized yield curve, see Lemieux & Dupont, 2020 856 ! Recommendation 1 : we use ice strength, not replacement pressure 857 ! Recommendation 2 : need to use deformations at PREVIOUS iterate for viscosities 858 !!$ IF( iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN 859 !!$ ! 860 !!$ ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) , zsig_I(jpi,jpj) , zsig_II(jpi,jpj) ) 861 !!$ ! 862 !!$ DO jj = 1, jpj 863 !!$ DO ji = 1, jpi 864 !!$ 865 !!$ ! Ice stresses computed with **viscosities** (delta, p/delta) at **previous** iterates 866 !!$ ! and **deformations** at current iterates 867 !!$ ! following Lemieux & Dupont (2020) 868 !!$ zfac = zp_delt(ji,jj) 869 !!$ zsig1 = zfac * ( pdivu_i(ji,jj) - ( zdelta(ji,jj) + rn_creepl ) ) 870 !!$ zsig2 = zfac * z1_ecc2 * zten_i(ji,jj) 871 !!$ zsig12 = zfac * z1_ecc2 * pshear_i(ji,jj) 872 !!$ 873 !!$ ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008), T-point 874 !!$ zsig_I(ji,jj) = zsig1 * 0.5_wp ! 1st stress invariant, aka average normal stress, aka negative pressure 875 !!$ zsig_II(ji,jj) = SQRT ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) ) ! 2nd '' '', aka maximum shear stress 876 !!$ 877 !!$ ! Normalized principal stresses (used to display the ellipse) 878 !!$ z1_strength = 1._wp / MAX( 1._wp, strength(ji,jj) ) 879 !!$ zsig1_p(ji,jj) = ( zsig_I(ji,jj) + zsig_II(ji,jj) ) * z1_strength 880 !!$ zsig2_p(ji,jj) = ( zsig_I(ji,jj) - zsig_II(ji,jj) ) * z1_strength 881 !!$ END DO 882 !!$ END DO 883 !!$ ! 884 !!$ CALL iom_put( 'sig1_pnorm' , zsig1_p ) 885 !!$ CALL iom_put( 'sig2_pnorm' , zsig2_p ) 886 !!$ 887 !!$ DEALLOCATE( zsig1_p , zsig2_p , zsig_I, zsig_II ) 888 !!$ 889 !!$ ENDIF 854 890 855 891 ! --- SIMIP --- ! -
NEMO/releases/r4.0/r4.0-HEAD/src/OCE/TRA/trabbl.F90
r11536 r13646 513 513 IF( tra_bbl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'tra_bbl_init : unable to allocate arrays' ) 514 514 ! 515 IF( nn_bbl_adv == 1 ) WRITE(numout,*) ' * Advective BBL using upper velocity' 516 IF( nn_bbl_adv == 2 ) WRITE(numout,*) ' * Advective BBL using velocity = F( delta rho)' 515 IF(lwp) THEN 516 IF( nn_bbl_adv == 1 ) WRITE(numout,*) ' * Advective BBL using upper velocity' 517 IF( nn_bbl_adv == 2 ) WRITE(numout,*) ' * Advective BBL using velocity = F( delta rho)' 518 ENDIF 517 519 ! 518 520 ! !* vertical index of "deep" bottom u- and v-points
Note: See TracChangeset
for help on using the changeset viewer.