Changeset 15741
- Timestamp:
- 2022-03-08T16:09:40+01:00 (2 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/UKMO/NEMO_4.0.4_ice_strength_EAP/src/ICE/icedyn_rhg_eap.F90
r15571 r15741 155 155 REAL(wp), DIMENSION(jpi,jpj) :: zyield11, zyield22, zyield12 ! yield surface tensor for history 156 156 REAL(wp), DIMENSION(jpi,jpj) :: zdelta, zp_delt ! delta and P/delta at T points 157 REAL(wp), DIMENSION(jpi,jpj) :: zten_i ! tension157 REAL(wp), DIMENSION(jpi,jpj) :: zten_i, zshear ! tension, shear 158 158 REAL(wp), DIMENSION(jpi,jpj) :: zbeta ! beta coef from Kimmritz 2017 159 159 ! … … 838 838 & ) * 0.25_wp * r1_e1e2t(ji,jj) 839 839 840 ! shear at T points 840 841 ! maximum shear rate at T points (includes tension, output only) 841 842 pshear_i(ji,jj) = SQRT( zdt2 + zds2 ) 843 844 ! shear at T-points 845 zshear(ji,jj) = SQRT( zds2 ) 842 846 843 847 ! divergence at T points … … 854 858 END DO 855 859 CALL lbc_lnk_multi( 'icedyn_rhg_eap', pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1., zten_i, 'T', 1., & 856 & zs1 , 'T', 1., zs2 , 'T', 1., zs12 , 'F', 1. )860 & zs1 , 'T', 1., zs2 , 'T', 1., zs12 , 'F', 1., zshear, 'T', 1. ) 857 861 858 862 ! --- Store the stress tensor for the next time step --- ! … … 901 905 zsig1 = zfac * ( pdivu_i(ji,jj) - pdelta_i(ji,jj) ) 902 906 zsig2 = zfac * z1_ecc2 * zten_i(ji,jj) 903 zsig12 = zfac * z1_ecc2 * pshear_i(ji,jj)907 zsig12 = zfac * z1_ecc2 * zshear(ji,jj) * 0.5_wp 904 908 905 909 ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008) 906 910 zsig_I (ji,jj) = zsig1 * 0.5_wp ! 1st stress invariant, aka average normal stress, aka negative pressure 907 zsig_II(ji,jj) = SQRT ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) ) ! 2nd '' '', aka maximum shear stress911 zsig_II(ji,jj) = SQRT( zsig2 * zsig2 * 0.25_wp + zsig12 * zsig12 ) ! 2nd '' '' , aka maximum shear stress 908 912 909 913 END DO … … 932 936 ! and **deformations** at current iterates 933 937 ! following Lemieux & Dupont (2020) 934 zfac = zp_delt(ji,jj)935 zsig1 = zfac * ( pdivu_i(ji,jj) - ( zdelta(ji,jj) + rn_creepl) )938 zfac = strength(ji,jj) / ( pdelta_i(ji,jj) + rn_creepl ) 939 zsig1 = zfac * ( pdivu_i(ji,jj) - zdelta(ji,jj) ) 936 940 zsig2 = zfac * z1_ecc2 * zten_i(ji,jj) 937 zsig12 = zfac * z1_ecc2 * pshear_i(ji,jj)941 zsig12 = zfac * z1_ecc2 * zshear(ji,jj) * 0.5_wp 938 942 939 943 ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008), T-point 940 944 zsig_I(ji,jj) = zsig1 * 0.5_wp ! 1st stress invariant, aka average normal stress, aka negative pressure 941 zsig_II(ji,jj) = SQRT ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) ) ! 2nd '' '', aka maximum shear stress945 zsig_II(ji,jj) = SQRT( zsig2 * zsig2 * 0.25_wp + zsig12 * zsig12 ) ! 2nd '' '' , aka maximum shear stress 942 946 943 947 ! Normalized principal stresses (used to display the ellipse) … … 948 952 END DO 949 953 ! 950 CALL iom_put( 'sig1_pnorm' , zsig1_p )951 CALL iom_put( 'sig2_pnorm' , zsig2_p)954 IF( iom_use('sig1_pnorm') ) CALL iom_put( 'sig1_pnorm' , zsig1_p(:,:) * zmsk00(:,:) ) 955 IF( iom_use('sig2_pnorm') ) CALL iom_put( 'sig2_pnorm' , zsig2_p(:,:) * zmsk00(:,:) ) 952 956 953 957 DEALLOCATE( zsig1_p , zsig2_p , zsig_I, zsig_II ) … … 1088 1092 1089 1093 ! time 1090 it = ( kt - 1) * kitermax + kiter1094 it = ( kt - nit000 ) * kitermax + kiter 1091 1095 1092 1096 ! convergence … … 1108 1112 istatus = NF90_PUT_VAR( ncvgid, nvarid, (/zresm/), (/it/), (/1/) ) 1109 1113 ! close file 1110 IF( kt == nitend ) istatus = NF90_CLOSE(ncvgid)1114 IF( kt == nitend - nn_fsbc + 1 .AND. kiter == kitermax ) istatus = NF90_CLOSE(ncvgid) 1111 1115 ENDIF 1112 1116
Note: See TracChangeset
for help on using the changeset viewer.