Changeset 14998 for NEMO/branches
- Timestamp:
- 2021-06-16T10:06:19+02:00 (3 years ago)
- Location:
- NEMO/branches/UKMO/NEMO_4.0.4_EAP_rheology
- Files:
-
- 20 added
- 6 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/UKMO/NEMO_4.0.4_EAP_rheology/cfgs/SHARED/field_def_nemo-ice.xml
r13648 r14998 81 81 <field id="sig1_pnorm" long_name="P-normalized 1st principal stress component" unit="" /> 82 82 <field id="sig2_pnorm" long_name="P-normalized 2nd principal stress component" unit="" /> 83 <field id="icedlt" long_name="delta" standard_name="delta" unit="" /> 83 84 <field id="normstr" long_name="Average normal stress in sea ice" standard_name="average_normal_stress" unit="N/m" /> 84 85 <field id="sheastr" long_name="Maximum shear stress in sea ice" standard_name="maximum_shear_stress" unit="N/m" /> … … 86 87 <field id="icediv" long_name="Divergence of the sea-ice velocity field" standard_name="divergence_of_sea_ice_velocity" unit="s-1" /> 87 88 <field id="iceshe" long_name="Maximum shear of sea-ice velocity field" standard_name="maximum_shear_of_sea_ice_velocity" unit="s-1" /> 89 <field id="aniso" long_name="anisotropy of sea ice floe orientation (0.5 - 1)" standard_name="anisotropy" unit="" /> 90 <field id="yield11" long_name="yield surface tensor component 11" standard_name="yield11" unit="N/m" /> 91 <field id="yield22" long_name="yield surface tensor component 22" standard_name="yield22" unit="N/m" /> 92 <field id="yield12" long_name="yield surface tensor component 12" standard_name="yield12" unit="N/m" /> 88 93 <field id="beta_evp" long_name="Relaxation parameter of ice rheology (beta)" standard_name="relaxation_parameter_of_ice_rheology" unit="" /> 94 89 95 90 96 <!-- surface heat fluxes --> … … 167 173 <!-- diags --> 168 174 <field id="hfxdhc" long_name="Heat content variation in snow and ice (neg = ice cooling)" unit="W/m2" /> 175 <field id="icedrift_mass" long_name="Ice mass drift" unit="kg/m2/s" /> 176 <field id="icedrift_salt" long_name="Ice salt drift" unit="g/m2/s" /> 177 <field id="icedrift_heat" long_name="Ice heat drift" unit="W/m2" /> 169 178 170 179 <!-- sbcssm variables --> … … 401 410 <field field_ref="normstr" name="normstr" /> 402 411 <field field_ref="sheastr" name="sheastr" /> 412 <field field_ref="sig1_pnorm" name="sig1_pnorm"/> 413 <field field_ref="sig2_pnorm" name="sig2_pnorm"/> 414 <field field_ref="icedlt" name="sidelta" /> 403 415 404 416 <!-- heat fluxes --> -
NEMO/branches/UKMO/NEMO_4.0.4_EAP_rheology/src/ICE/ice.F90
r14075 r14998 150 150 ! 151 151 ! !!** ice-rheology namelist (namdyn_rhg) ** 152 LOGICAL , PUBLIC :: ln_rhg_EVP ! EVP rheology switch, used for rdgrft and rheology 153 LOGICAL , PUBLIC :: ln_rhg_EAP ! EAP rheology switch, used for rdgrft and rheology 152 154 LOGICAL , PUBLIC :: ln_aEVP !: using adaptive EVP (T or F) 153 155 REAL(wp), PUBLIC :: rn_creepl !: creep limit : has to be under 1.0e-9 … … 246 248 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: divu_i !: Divergence of the velocity field [s-1] 247 249 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: shear_i !: Shear of the velocity field [s-1] 250 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: aniso_11, aniso_12 !: structure tensor elements 251 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rdg_conv 248 252 ! 249 253 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: t_bo !: Sea-Ice bottom temperature [Kelvin] … … 436 440 ALLOCATE( u_oce (jpi,jpj) , v_oce (jpi,jpj) , ht_i_new (jpi,jpj) , strength(jpi,jpj) , & 437 441 & stress1_i(jpi,jpj) , stress2_i(jpi,jpj) , stress12_i(jpi,jpj) , & 438 & delta_i (jpi,jpj) , divu_i (jpi,jpj) , shear_i (jpi,jpj) , STAT=ierr(ii) ) 442 & delta_i (jpi,jpj) , divu_i (jpi,jpj) , shear_i (jpi,jpj) , & 443 & aniso_11 (jpi,jpj) , aniso_12 (jpi,jpj) , rdg_conv (jpi,jpj) , STAT=ierr(ii) ) 439 444 440 445 ii = ii + 1 -
NEMO/branches/UKMO/NEMO_4.0.4_EAP_rheology/src/ICE/icedyn_rdgrft.F90
r14075 r14998 23 23 USE icevar ! sea-ice: operations 24 24 USE icectl ! sea-ice: control prints 25 !USE icedyn_rhg ! sea-ice: rheology choice 25 26 ! 26 27 USE in_out_manager ! I/O manager … … 138 139 INTEGER , DIMENSION(jpij) :: iptidx ! compute ridge/raft or not 139 140 REAL(wp), DIMENSION(jpij) :: zdivu, zdelt ! 1D divu_i & delta_i 141 REAL(wp), DIMENSION(jpij) :: zconv ! 1D rdg_conv (if EAP rheology) 140 142 ! 141 143 INTEGER, PARAMETER :: jp_itermax = 20 … … 174 176 175 177 ! just needed here 176 CALL tab_2d_1d( npti, nptidx(1:npti), zdelt (1:npti) , delta_i ) 178 CALL tab_2d_1d( npti, nptidx(1:npti), zdelt (1:npti) , delta_i ) 179 CALL tab_2d_1d( npti, nptidx(1:npti), zconv (1:npti) , rdg_conv ) 177 180 ! needed here and in the iteration loop 178 181 CALL tab_2d_1d( npti, nptidx(1:npti), zdivu (1:npti) , divu_i) ! zdivu is used as a work array here (no change in divu_i) … … 185 188 ! - ice area added in new ridges 186 189 closing_net(ji) = rn_csrdg * 0.5_wp * ( zdelt(ji) - ABS( zdivu(ji) ) ) - MIN( zdivu(ji), 0._wp ) 190 IF( ln_rhg_EVP ) closing_net(ji) = rn_csrdg * 0.5_wp * ( zdelt(ji) - ABS( zdivu(ji) ) ) - MIN( zdivu(ji), 0._wp ) 191 IF( ln_rhg_EAP ) closing_net(ji) = zconv(ji) 187 192 ! 188 193 IF( zdivu(ji) < 0._wp ) closing_net(ji) = MAX( closing_net(ji), -zdivu(ji) ) ! make sure the closing rate is large enough … … 768 773 ! !--------------------------------------------------! 769 774 strength(:,:) = rn_pstar * SUM( v_i(:,:,:), dim=3 ) * EXP( -rn_crhg * ( 1._wp - SUM( a_i(:,:,:), dim=3 ) ) ) 770 ismooth = 1 775 ! ismooth = 1 ! original code 776 ismooth = 0 ! try for EAP stability 771 777 ! !--------------------------------------------------! 772 778 ELSE ! Zero strength ! -
NEMO/branches/UKMO/NEMO_4.0.4_EAP_rheology/src/ICE/icedyn_rhg.F90
r14075 r14998 17 17 USE ice ! sea-ice: variables 18 18 USE icedyn_rhg_evp ! sea-ice: EVP rheology 19 USE icedyn_rhg_eap ! sea-ice: EAP rheology 19 20 USE icectl ! sea-ice: control prints 20 21 ! … … 33 34 ! ! associated indices: 34 35 INTEGER, PARAMETER :: np_rhgEVP = 1 ! EVP rheology 35 !!INTEGER, PARAMETER :: np_rhgEAP = 2 ! EAP rheology36 INTEGER, PARAMETER :: np_rhgEAP = 2 ! EAP rheology 36 37 37 38 ! ** namelist (namrhg) ** 38 LOGICAL :: ln_rhg_EVP ! EVP rheology 39 ! LOGICAL :: ln_rhg_EVP ! EVP rheology, used by icedyn_rdgrft 40 ! LOGICAL :: ln_rhg_EAP ! EAP rheology, used by icedyn_rdgrft 41 !LOGICAL, PUBLIC :: ln_rhg_EVP ! EVP rheology, used by icedyn_rdgrft 42 !LOGICAL, PUBLIC :: ln_rhg_EAP ! EAP rheology, used by icedyn_rdgrft 39 43 ! 40 44 !! * Substitutions … … 81 85 CALL ice_dyn_rhg_evp( kt, stress1_i, stress2_i, stress12_i, shear_i, divu_i, delta_i ) 82 86 ! 87 ! !----------------------------! 88 CASE( np_rhgEAP ) ! Elasto-Anisotropic-Plastic ! 89 ! !----------------------------! 90 CALL ice_dyn_rhg_eap( kt, stress1_i, stress2_i, stress12_i, shear_i, divu_i, delta_i, aniso_11, aniso_12, rdg_conv ) 83 91 END SELECT 84 92 ! 85 IF( lrst_ice ) THEN !* write EVP fields in the restart file 86 IF( ln_rhg_EVP ) CALL rhg_evp_rst( 'WRITE', kt ) 93 IF( lrst_ice ) THEN 94 IF( ln_rhg_EVP ) CALL rhg_evp_rst( 'WRITE', kt ) !* write EVP fields in the restart file 95 IF( ln_rhg_EAP ) CALL rhg_eap_rst( 'WRITE', kt ) !* write EAP fields in the restart file 87 96 ENDIF 88 97 ! … … 110 119 INTEGER :: ios, ioptio ! Local integer output status for namelist read 111 120 !! 112 NAMELIST/namdyn_rhg/ ln_rhg_EVP, ln_aEVP, rn_creepl, rn_ecc , nn_nevp, rn_relast, nn_rhg_chkcvg121 NAMELIST/namdyn_rhg/ ln_rhg_EVP, ln_aEVP, ln_rhg_EAP, rn_creepl, rn_ecc , nn_nevp, rn_relast, nn_rhg_chkcvg 113 122 !!------------------------------------------------------------------- 114 123 ! … … 137 146 ELSEIF( nn_rhg_chkcvg == 2 ) THEN ; WRITE(numout,*) ' check cvg at both main and rheology time steps' 138 147 ENDIF 148 WRITE(numout,*) ' rheology EAP (icedyn_rhg_eap) ln_rhg_EAP = ', ln_rhg_EAP 139 149 ENDIF 140 150 ! … … 142 152 ioptio = 0 143 153 IF( ln_rhg_EVP ) THEN ; ioptio = ioptio + 1 ; nice_rhg = np_rhgEVP ; ENDIF 144 !!IF( ln_rhg_EAP ) THEN ; ioptio = ioptio + 1 ; nice_rhg = np_rhgEAP ; ENDIF154 IF( ln_rhg_EAP ) THEN ; ioptio = ioptio + 1 ; nice_rhg = np_rhgEAP ; ENDIF 145 155 IF( ioptio /= 1 ) CALL ctl_stop( 'ice_dyn_rhg_init: choose one and only one ice rheology' ) 146 156 ! 147 157 IF( ln_rhg_EVP ) CALL rhg_evp_rst( 'READ' ) !* read or initialize all required files 158 IF( ln_rhg_EAP ) CALL rhg_eap_rst( 'READ' ) !* read or initialize all required files 148 159 ! 149 160 END SUBROUTINE ice_dyn_rhg_init -
NEMO/branches/UKMO/NEMO_4.0.4_EAP_rheology/src/ICE/icedyn_rhg_evp.F90
r14075 r14998 818 818 IF( iom_use('icediv') ) CALL iom_put( 'icediv' , pdivu_i * zmsk00 ) ! divergence 819 819 IF( iom_use('iceshe') ) CALL iom_put( 'iceshe' , pshear_i * zmsk00 ) ! shear 820 IF( iom_use('icedlt') ) CALL iom_put( 'icedlt' , pdelta_i * zmsk00 ) ! delta 820 821 IF( iom_use('icestr') ) CALL iom_put( 'icestr' , strength * zmsk00 ) ! strength 821 822 … … 856 857 ! Recommendation 1 : we use ice strength, not replacement pressure 857 858 ! Recommendation 2 : need to use deformations at PREVIOUS iterate for viscosities 858 !!$IF( iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN859 !!$!860 !!$ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) , zsig_I(jpi,jpj) , zsig_II(jpi,jpj) )861 !!$!862 !!$DO jj = 1, jpj863 !!$DO ji = 1, jpi864 !!$865 !!$! Ice stresses computed with **viscosities** (delta, p/delta) at **previous** iterates866 !!$! and **deformations** at current iterates867 !!$! 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-point874 !!$zsig_I(ji,jj) = zsig1 * 0.5_wp ! 1st stress invariant, aka average normal stress, aka negative pressure875 !!$zsig_II(ji,jj) = SQRT ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) ) ! 2nd '' '', aka maximum shear stress876 !!$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_strength880 !!$zsig2_p(ji,jj) = ( zsig_I(ji,jj) - zsig_II(ji,jj) ) * z1_strength881 !!$END DO882 !!$END DO883 !!$!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 !!$ENDIF859 IF( iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN 860 ! 861 ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) , zsig_I(jpi,jpj) , zsig_II(jpi,jpj) ) 862 ! 863 DO jj = 1, jpj 864 DO ji = 1, jpi 865 866 ! Ice stresses computed with **viscosities** (delta, p/delta) at **previous** iterates 867 ! and **deformations** at current iterates 868 ! following Lemieux & Dupont (2020) 869 zfac = zp_delt(ji,jj) 870 zsig1 = zfac * ( pdivu_i(ji,jj) - ( zdelta(ji,jj) + rn_creepl ) ) 871 zsig2 = zfac * z1_ecc2 * zten_i(ji,jj) 872 zsig12 = zfac * z1_ecc2 * pshear_i(ji,jj) 873 874 ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008), T-point 875 zsig_I(ji,jj) = zsig1 * 0.5_wp ! 1st stress invariant, aka average normal stress, aka negative pressure 876 zsig_II(ji,jj) = SQRT ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) ) ! 2nd '' '', aka maximum shear stress 877 878 ! Normalized principal stresses (used to display the ellipse) 879 z1_strength = 1._wp / MAX( 1._wp, strength(ji,jj) ) 880 zsig1_p(ji,jj) = ( zsig_I(ji,jj) + zsig_II(ji,jj) ) * z1_strength 881 zsig2_p(ji,jj) = ( zsig_I(ji,jj) - zsig_II(ji,jj) ) * z1_strength 882 END DO 883 END DO 884 ! 885 CALL iom_put( 'sig1_pnorm' , zsig1_p ) 886 CALL iom_put( 'sig2_pnorm' , zsig2_p ) 887 888 DEALLOCATE( zsig1_p , zsig2_p , zsig_I, zsig_II ) 889 890 ENDIF 890 891 891 892 ! --- SIMIP --- ! -
NEMO/branches/UKMO/NEMO_4.0.4_EAP_rheology/src/ICE/icewri.F90
r14075 r14998 254 254 CALL iom_rstput( 0, 0, kid, 'sivolu', vt_i ) ! Ice volume 255 255 CALL iom_rstput( 0, 0, kid, 'sidive', divu_i*1.0e8 ) ! Ice divergence 256 CALL iom_rstput( 0, 0, kid, 'sishea', shear_i*1.0e8 ) ! Ice shear 256 257 CALL iom_rstput( 0, 0, kid, 'si_amp', at_ip ) ! Melt pond fraction 257 258 CALL iom_rstput( 0, 0, kid, 'si_vmp', vt_ip ) ! Melt pond volume
Note: See TracChangeset
for help on using the changeset viewer.