Skip to content
Snippets Groups Projects
Commit 4a17a1f8 authored by Martin Vancoppenolle's avatar Martin Vancoppenolle Committed by Sebastien Masson
Browse files

ADD all files related to commit

parent 2c2c8a48
No related branches found
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
cfgs/ORCA2_SAS_ICE/EXPREF/file_def_nemo-ice.xml
+ 8
− 1
View file @ 4a17a1f8
......@@ -48,11 +48,18 @@
<!-- rheology -->
<field field_ref="icediv" name="sidive" />
<field field_ref="iceshe" name="sishea" />
<field field_ref="icedlt" name="sidelt" />
<field field_ref="icestr" name="sistre" />
<field field_ref="normstr" name="normstr" />
<field field_ref="sheastr" name="sheastr" />
<field field_ref="sig1_pnorm" name="sig1_pnorm"/>
<field field_ref="sig2_pnorm" name="sig2_pnorm"/>
<field field_ref="velo_res" name="velo_res" />
<field field_ref="velo_ero" name="velo_ero" />
<field field_ref="uice_eri" name="uice_eri" />
<field field_ref="vice_eri" name="vice_eri" />
<field field_ref="uice_cvg" name="uice_cvg" />
<!-- heat fluxes -->
<field field_ref="qt_oce_ai" name="qt_oce_ai" />
......@@ -64,7 +71,7 @@
<field field_ref="qns_ice" name="qns_ice" />
<field field_ref="qemp_ice" name="qemp_ice" />
<field field_ref="albedo" name="albedo" />
<field field_ref="hfxcndtop" name="hfxcndtop" />
<field field_ref="hfxcndbot" name="hfxcndbot" />
<field field_ref="hfxsensib" name="hfxsensib" />
......
cfgs/SHARED/field_def_nemo-ice.xml
+ 7
− 1
View file @ 4a17a1f8
......@@ -200,6 +200,12 @@
<!-- rheology convergence tests -->
<field id="uice_cvg" long_name="sea ice velocity convergence" standard_name="sea_ice_velocity_convergence" unit="m/s" />
<!-- vp rheology convergence tests -->
<field id="velo_res" long_name="sea ice velocity residual" standard_name="sea_ice_velocity_residual" unit="N/m2" />
<field id="velo_ero" long_name="sea ice velocity error last outer iteration" standard_name="sea_ice_velocity_outer_error" unit="m/s" />
<field id="uice_eri" long_name="uice velocity error last inner iteration" standard_name="sea_ice_u_velocity_inner_error" unit="m/s" />
<field id="vice_eri" long_name="vice velocity error last inner iteration" standard_name="sea_ice_v_velocity_inner_error" unit="m/s" />
<!-- ================= -->
<!-- Add-ons for SIMIP -->
<!-- ================= -->
......@@ -429,7 +435,7 @@
<field field_ref="sheastr" name="sheastr" />
<field field_ref="sig1_pnorm" name="sig1_pnorm"/>
<field field_ref="sig2_pnorm" name="sig2_pnorm"/>
<field field_ref="icedlt" name="sidelta" />
<field field_ref="icedlt" name="sidelt" />
<!-- heat fluxes -->
<field field_ref="qt_oce_ai" name="qt_oce_ai" />
......
src/ICE/icedyn_rhg_evp.F90
+ 31
− 27
View file @ 4a17a1f8
......@@ -145,7 +145,7 @@ CONTAINS
REAL(wp), DIMENSION(jpi,jpj) :: v_oceU, u_oceV, v_iceU, u_iceV ! ocean/ice u/v component on V/U points
!
REAL(wp), DIMENSION(jpi,jpj) :: zds ! shear
REAL(wp), DIMENSION(jpi,jpj) :: zten_i ! tension
REAL(wp), DIMENSION(jpi,jpj) :: zten_i, zshear ! tension, shear
REAL(wp), DIMENSION(jpi,jpj) :: zs1, zs2, zs12 ! stress tensor components
REAL(wp), DIMENSION(jpi,jpj) :: zsshdyn ! array used for the calculation of ice surface slope:
! ! ocean surface (ssh_m) if ice is not embedded
......@@ -749,8 +749,11 @@ CONTAINS
& + zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1) &
& ) * 0.25_wp * r1_e1e2t(ji,jj)
! shear at T points
pshear_i(ji,jj) = SQRT( zdt2 + zds2 ) * zmsk(ji,jj)
! maximum shear rate at T points (includes tension, output only)
pshear_i(ji,jj) = SQRT( zdt2 + zds2 ) * zmsk(ji,jj) !
! shear at T-points
zshear(ji,jj) = SQRT( zds2 ) * zmsk(ji,jj)
! divergence at T points
pdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) &
......@@ -758,21 +761,25 @@ CONTAINS
& ) * r1_e1e2t(ji,jj) * zmsk(ji,jj)
! delta at T points
zfac = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) * zmsk(ji,jj) ! delta
rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zfac ) ) ! 0 if delta=0
pdelta_i(ji,jj) = zfac + rn_creepl * rswitch ! delta+creepl
zdelta(ji,jj) = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) * zmsk(ji,jj) ! delta
! delta* at T points (pdelta_i)
rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta(ji,jj) ) ) ! 0 if delta=0
pdelta_i(ji,jj) = zdelta(ji,jj) + rn_creepl * rswitch
! it seems that deformation used for advection and mech redistribution is delta*
! MV in principle adding creep limit is a regularization for viscosity not for delta
! delta_star should not (in my view) be used in a replacement for delta
END_2D
CALL lbc_lnk( 'icedyn_rhg_evp', pshear_i, 'T', 1._wp, pdivu_i, 'T', 1._wp, pdelta_i, 'T', 1._wp, zten_i, 'T', 1._wp, &
& zs1 , 'T', 1._wp, zs2 , 'T', 1._wp, zs12 , 'F', 1._wp )
& zdelta , 'T', 1._wp, zs1 , 'T', 1._wp, zs2 , 'T', 1._wp, zs12 , 'F', 1._wp )
! --- Store the stress tensor for the next time step --- !
pstress1_i (:,:) = zs1 (:,:)
pstress2_i (:,:) = zs2 (:,:)
pstress12_i(:,:) = zs12(:,:)
!
!------------------------------------------------------------------------------!
! 5) diagnostics
!------------------------------------------------------------------------------!
! --- ice-ocean, ice-atm. & ice-oceanbottom(landfast) stresses --- !
......@@ -795,7 +802,7 @@ CONTAINS
IF( iom_use('icediv') ) CALL iom_put( 'icediv' , pdivu_i * zmsk00 ) ! divergence
IF( iom_use('iceshe') ) CALL iom_put( 'iceshe' , pshear_i * zmsk00 ) ! shear
IF( iom_use('icestr') ) CALL iom_put( 'icestr' , strength * zmsk00 ) ! strength
IF( iom_use('icedlt') ) CALL iom_put( 'icedlt' , pdelta_i * zmsk00 ) ! delta
IF( iom_use('icedlt') ) CALL iom_put( 'icedlt' , zdelta * zmsk00 ) ! delta
! --- Stress tensor invariants (SIMIP diags) --- !
IF( iom_use('normstr') .OR. iom_use('sheastr') ) THEN
......@@ -805,21 +812,19 @@ CONTAINS
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
! Ice stresses
! sigma1, sigma2, sigma12 are some useful recombination of the stresses (Hunke and Dukowicz MWR 2002, Bouillon et al., OM2013)
! These are NOT stress tensor components, neither stress invariants, neither stress principal components
! I know, this can be confusing...
zfac = strength(ji,jj) / ( pdelta_i(ji,jj) + rn_creepl )
zsig1 = zfac * ( pdivu_i(ji,jj) - pdelta_i(ji,jj) )
! sigma1, sigma2, sigma12 are some recombination of the stresses (HD MWR002, Bouillon et al., OM2013)
! not to be confused with stress tensor components, stress invariants, or stress principal components
zfac = strength(ji,jj) / pdelta_i(ji,jj) ! viscosity
zsig1 = zfac * ( pdivu_i(ji,jj) - zdelta(ji,jj) )
zsig2 = zfac * z1_ecc2 * zten_i(ji,jj)
zsig12 = zfac * z1_ecc2 * pshear_i(ji,jj)
zsig12 = zfac * z1_ecc2 * zshear(ji,jj) * 0.5_wp
! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008)
zsig_I (ji,jj) = zsig1 * 0.5_wp ! 1st stress invariant, aka average normal stress, aka negative pressure
zsig_II(ji,jj) = SQRT ( zsig2 * zsig2 * 0.25_wp + zsig12 * zsig12 ) ! 2nd '' '' , aka maximum shear stress
zsig_I (ji,jj) = 0.5_wp * zsig1
zsig_II(ji,jj) = 0.5_wp * SQRT ( zsig2 * zsig2 + 4. * zsig12 * zsig12 )
END_2D
!
! Stress tensor invariants (normal and shear stress N/m) - SIMIP diags - definitions following Coon (1974) and Feltham (2008)
IF( iom_use('normstr') ) CALL iom_put( 'normstr', zsig_I (:,:) * zmsk00(:,:) ) ! Normal stress
IF( iom_use('sheastr') ) CALL iom_put( 'sheastr', zsig_II(:,:) * zmsk00(:,:) ) ! Maximum shear stress
......@@ -830,24 +835,23 @@ CONTAINS
! --- Normalized stress tensor principal components --- !
! This are used to plot the normalized yield curve, see Lemieux & Dupont, 2020
! Recommendation 1 : we use ice strength, not replacement pressure
! Recommendation 2 : need to use deformations at PREVIOUS iterate for viscosities
! Recommendation 2 : for EVP, no need to use viscosities at last iteration (stress is properly iterated)
IF( iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN
!
ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) , zsig_I(jpi,jpj) , zsig_II(jpi,jpj) )
!
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
! Ice stresses computed with **viscosities** (delta, p/delta) at **previous** iterates
! and **deformations** at current iterates
! For EVP solvers, ice stresses at current iterates can be used
! following Lemieux & Dupont (2020)
zfac = zp_delt(ji,jj)
zsig1 = zfac * ( pdivu_i(ji,jj) - ( zdelta(ji,jj) + rn_creepl ) )
zfac = strength(ji,jj) / pdelta_i(ji,jj)
zsig1 = zfac * ( pdivu_i(ji,jj) - zdelta(ji,jj) )
zsig2 = zfac * z1_ecc2 * zten_i(ji,jj)
zsig12 = zfac * z1_ecc2 * pshear_i(ji,jj)
zsig12 = zfac * z1_ecc2 * zshear(ji,jj) * 0.5_wp
! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008), T-point
zsig_I(ji,jj) = zsig1 * 0.5_wp ! 1st stress invariant, aka average normal stress, aka negative pressure
zsig_II(ji,jj) = SQRT ( zsig2 * zsig2 * 0.25_wp + zsig12 * zsig12 ) ! 2nd '' '' , aka maximum shear stress
zsig_I(ji,jj) = 0.5_wp * zsig1 ! normal stress
zsig_II(ji,jj) = 0.5_wp * SQRT ( zsig2 * zsig2 + 4. * zsig12 * zsig12 ) ! max shear stress
! Normalized principal stresses (used to display the ellipse)
z1_strength = 1._wp / MAX( 1._wp, strength(ji,jj) )
......
src/ICE/icedyn_rhg_vp.F90
+ 325
− 289
View file @ 4a17a1f8
This diff is collapsed.
tests/ICE_ADV2D/EXPREF/README
+ 7
− 2
View file @ 4a17a1f8
!! comment from MV
!! not clear why one has both EXP00 and EXPREF directories
!! end MV
WARNING: For now, the test case ICE_ADV2D only works if the logical "ll_neg" is set to FALSE in the routine icedyn_adv_umx.F90
(it is still unclear why)
-------
......@@ -28,8 +33,8 @@ How to run
----------
a) Compile and run the model once to get a mesh_mask.nc file with the following command:
../../../makenemo -r ICE_ADV2D -n ICE_ADV2D -m X64_ADA -j 4
mpirun ./nemo -np 1
./makenemo -n ICE_ADV2D -a ICE_ADV2D -m X64_JEANZAY -j8
./nemo.exe (on jeanzaypp)
b) Create the initial condition file for sea-ice (initice.nc) by running this python script:
python ./make_INITICE.py
......
tests/ICE_ADV2D/EXPREF/make_INITICE.py
+ 7
− 7
View file @ 4a17a1f8
......@@ -15,18 +15,18 @@ fcoord='mesh_mask.nc'
# output file
fflx='initice.nc'
print ' creating init ice file ' +fflx
print (' creating init ice file ' +fflx)
# Reading coordinates file
nccoord=netcdf(fcoord,'r')
nav_lon=nccoord.variables['nav_lon']
nav_lat=nccoord.variables['nav_lat']
nav_lon=nccoord.variables['glamt']
nav_lat=nccoord.variables['gphit']
time_counter=1
LON1= nav_lon.shape[1]
LAT1= nav_lon.shape[0]
print 'nav_lon.shape[1]' ,nav_lon.shape[1]
print 'LON1 ', LON1
print 'LAT1 ', LAT1
LAT1= nav_lon.shape[1]
print ('nav_lon.shape[1]' ,nav_lon.shape[1])
print ('LON1 ', LON1)
print ('LAT1 ', LAT1)
# Creating INITICE netcdf file
nc=netcdf(fflx,'w')
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment