Changeset 13536 for NEMO/branches/2020
- Timestamp:
- 2020-09-29T11:22:04+02:00 (4 years ago)
- Location:
- NEMO/branches/2020/SI3-03_VP_rheology
- Files:
-
- 6 edited
-
cfgs/SHARED/field_def_nemo-ice.xml (modified) (2 diffs)
-
cfgs/SHARED/namelist_ice_ref (modified) (2 diffs)
-
src/ICE/ice.F90 (modified) (1 diff)
-
src/ICE/icedyn_rhg.F90 (modified) (5 diffs)
-
src/ICE/icedyn_rhg_evp.F90 (modified) (10 diffs)
-
src/ICE/icedyn_rhg_vp.F90 (modified) (29 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/SI3-03_VP_rheology/cfgs/SHARED/field_def_nemo-ice.xml
r13064 r13536 75 75 <field id="utau_bi" long_name="X-component of ocean bottom stress on sea ice -landfast" standard_name="ocean_bottom_upward_x_stress" unit="N/m2" /> 76 76 <field id="vtau_bi" long_name="Y-component of ocean bottom stress on sea ice -landfast" standard_name="ocean_bottom_upward_y_stress" unit="N/m2" /> 77 <field id="isig1" long_name="1st principal stress component for EVP rhg" unit="" /> 78 <field id="isig2" long_name="2nd principal stress component for EVP rhg" unit="" /> 79 <field id="isig3" long_name="convergence measure for EVP rheology (must be around 1)" unit="" /> 77 <field id="sig1_pnorm" long_name="P-normalized 1st principal stress component" unit="" /> 78 <field id="sig2_pnorm" long_name="P-normalized 2nd principal stress component" unit="" /> 80 79 <field id="normstr" long_name="Average normal stress in sea ice" standard_name="average_normal_stress" unit="N/m" /> 81 80 <field id="sheastr" long_name="Maximum shear stress in sea ice" standard_name="maximum_shear_stress" unit="N/m" /> … … 398 397 <field field_ref="normstr" name="normstr" /> 399 398 <field field_ref="sheastr" name="sheastr" /> 400 <field field_ref="isig1" name="isig1" /> 401 <field field_ref="isig2" name="isig2" /> 402 <field field_ref="isig3" name="isig3" /> 399 <field field_ref="sig1_pnorm" name="sig1_pnorm" /> 400 <field field_ref="sig2_pnorm" name="sig2_pnorm" /> 403 401 404 402 <!-- heat fluxes --> -
NEMO/branches/2020/SI3-03_VP_rheology/cfgs/SHARED/namelist_ice_ref
r12920 r13536 91 91 &namdyn_rhg ! Ice rheology 92 92 !------------------------------------------------------------------------------ 93 ln_rhg_EVP = . true.! EVP rheology93 ln_rhg_EVP = .false. ! EVP rheology 94 94 ln_aEVP = .true. ! adaptive rheology (Kimmritz et al. 2016 & 2017) 95 95 rn_creepl = 2.0e-9 ! creep limit [1/s] … … 99 99 ! advised value: 1/3 (rn_nevp=120) or 1/9 (rn_nevp=300) 100 100 ln_rhg_chkcvg = .false. ! check convergence of rheology (outputs: file ice_cvg.nc & variable uice_cvg) 101 ln_rhg_VP = .true. ! VP rheology 102 nn_nout_vp = 10 ! number of outer iterations 103 nn_ninn_vp = 1500 ! number of inner iterations 104 ln_zebra_vp = .false. ! activate zebra solver 105 rn_relaxu_vp = 0.95 ! relaxation factor for u-velocity 106 rn_relaxv_vp = 0.95 ! relaxation factor for v-velocity 107 rn_uerr_max_vp = 0.80 ! maximum error on velocity 108 rn_uerr_min_vp = 1.0e-04 ! velocity to decide convergence 109 nn_cvgchk_vp = 5 ! iteration step for convergence check 101 110 / 102 111 !------------------------------------------------------------------------------ -
NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/ice.F90
r13522 r13536 160 160 INTEGER , PUBLIC :: nn_nout_vp !: Number of outer iterations 161 161 INTEGER , PUBLIC :: nn_ninn_vp !: Number of inner iterations (linear system solver) 162 LOGICAL , PUBLIC :: ln_smooth_vp !: zeta viscosity smoothing (if yes, tanh function of Lemieux and Tremblay JGR 2009)163 162 LOGICAL , PUBLIC :: ln_zebra_vp !: activate zebra (solve the linear system problem every odd j-band, then one every even one) 164 163 REAL(wp), PUBLIC :: rn_relaxu_vp !: U-relaxation factor (MV: can probably be merged with V-factor once ok) -
NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg.F90
r13522 r13536 56 56 !! 57 57 !! ** Action : comupte - ice velocity (u_ice, v_ice) 58 !! - 3 components of the stress tensor (stress1_i, stress2_i, stress12_i) 58 !! - 3 components of the stress tensor (stress1_i, stress2_i, stress12_i) - EVP only 59 59 !! - shear, divergence and delta (shear_i, divu_i, delta_i) 60 60 !!-------------------------------------------------------------------- … … 88 88 CASE( np_rhgVP ) ! Viscous-Plastic rheology ! 89 89 ! !---------------------------------------------! 90 CALL ice_dyn_rhg_vp ( kt, s tress1_i, stress2_i, stress12_i, shear_i, divu_i, delta_i )90 CALL ice_dyn_rhg_vp ( kt, shear_i, divu_i, delta_i ) 91 91 ! 92 92 END SELECT … … 94 94 IF( lrst_ice ) THEN !* write EVP fields in the restart file 95 95 IF( ln_rhg_EVP ) CALL rhg_evp_rst( 'WRITE', kt ) 96 IF( ln_rhg_VP ) CALL rhg_lsr_rst( 'WRITE', kt )96 ! MV note: no restart needed for VP as there is no time equation for stress tensor 97 97 ENDIF 98 98 ! … … 121 121 !! 122 122 NAMELIST/namdyn_rhg/ ln_rhg_EVP, ln_aEVP, rn_creepl, rn_ecc , nn_nevp, rn_relast, ln_rhg_chkcvg, & !- evp 123 & ln_rhg_VP, nn_nout_vp, nn_ninn_vp, ln_ smooth_vp, ln_zebra_vp, rn_relaxu_vp, rn_relaxv_vp, rn_uerr_max_vp, rn_uerr_min_vp, nn_cvgchk_vp !-vp123 & ln_rhg_VP, nn_nout_vp, nn_ninn_vp, ln_zebra_vp, rn_relaxu_vp, rn_relaxv_vp, rn_uerr_max_vp, rn_uerr_min_vp, nn_cvgchk_vp !-vp 124 124 !!------------------------------------------------------------------- 125 125 … … 149 149 WRITE(numout,*) ' number of outer iterations nn_nout_vp = ', nn_nout_vp 150 150 WRITE(numout,*) ' number of inner iterations nn_ninn_vp = ', nn_ninn_vp 151 WRITE(numout,*) ' tanh viscosity smooothing ln_smooth_vp = ', ln_smooth_vp152 151 WRITE(numout,*) ' activate zebra solver ln_zebra_vp = ', ln_zebra_vp 153 152 WRITE(numout,*) ' relaxation factor for u rn_relaxu_vp = ', rn_relaxu_vp -
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 -
NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg_vp.F90
r13522 r13536 1 ! TODOLIST 2 ! 3 ! Define all symbols 4 ! - Do viscosity smoothing with sum (differentiable rheology) 5 ! - Re-calculate internal force diagnostic (end of the routine) 6 ! - Do proper masking of output fileds with ice mass (zmsk00 see EVP routine) 7 ! 8 ! quality control - compare code to mitGCM 9 ! quality control - comparing EVP and VP codes 10 ! 1 ! TODOLOIST 11 2 ! 12 3 ! Commit 13 4 ! Compile 14 5 ! 6 ! - Re-calculate internal force diagnostic (currently incorrect, see end of the routine) 7 ! - QC:compare code to mitGCM 8 ! - QC: comparing EVP and VP codes 9 ! 15 10 ! Clarify 16 11 ! --- Boundary conditions --> how to enforce them - is the fmask strategy sufficient ? 17 ! --- Swap of independent term in the U-V solvers ? 18 ! --- Is stress tensor for restart needed ? 19 ! --- Is stress tensor calculated at the end of the time step 12 ! --- Swap of independent term in the U-V solvers ? -> JF Lemieux says maybe 13 ! --- Is stress tensor for restart needed ? No!!! 20 14 ! 21 15 ! Test 22 ! --- Can we add the 15% mask in the convergence criteria ?23 16 ! --- Try ADI for u-v solver 24 17 ! 25 18 ! Add 26 19 ! - Charge ellipse as an output (good one from Lemieux and Dupont) 27 ! - Bottom drag 20 ! - Bottom drag (landfast ice) 28 21 ! - Tensile strength 29 22 ! - agrif … … 52 45 !!---------------------------------------------------------------------- 53 46 !! ice_dyn_rhg_vp : computes ice velocities from VP rheolog with LSR solvery 54 !! rhg_vp_rst : read/write EVP fields in ice restart55 47 !!---------------------------------------------------------------------- 56 48 USE phycst ! Physical constants … … 79 71 80 72 PUBLIC ice_dyn_rhg_vp ! called by icedyn_rhg.F90 81 PUBLIC rhg_vp_rst ! called by icedyn_rhg.F9082 73 83 74 !! for convergence tests … … 99 90 CONTAINS 100 91 101 SUBROUTINE ice_dyn_rhg_vp( kt, ps tress1_i, pstress2_i, pstress12_i, pshear_i, pdivu_i, pdelta_i )92 SUBROUTINE ice_dyn_rhg_vp( kt, pshear_i, pdivu_i, pdelta_i ) 102 93 !!------------------------------------------------------------------- 103 94 !! … … 164 155 !! 165 156 INTEGER , INTENT(in ) :: kt ! time step 166 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pstress1_i, pstress2_i, pstress12_i !167 157 REAL(wp), DIMENSION(:,:), INTENT( out) :: pshear_i , pdivu_i , pdelta_i ! 168 158 !! … … 180 170 REAL(wp) :: zkt ! isotropic tensile strength for landfast ice 181 171 REAL(wp) :: zm2, zm3, zmassU, zmassV ! ice/snow mass and volume 182 REAL(wp) :: zdeltat, zd eltat_star, zds2, zdt, zdt2, zdiv, zdiv2! temporary scalars172 REAL(wp) :: zdeltat, zds2, zdt, zdt2, zdiv, zdiv2 ! temporary scalars 183 173 REAL(wp) :: zp_deltastar_f 184 174 REAL(wp) :: zfac, zfac1, zfac2, zfac3 … … 192 182 REAL(wp), DIMENSION(jpi,jpj) :: zmassU_t, zmassV_t ! Mass per unit area divided by time step 193 183 ! 194 REAL(wp), DIMENSION(jpi,jpj) :: zp_deltastar_t ! P/delta at T points 184 REAL(wp), DIMENSION(jpi,jpj) :: zdeltastar_t ! Delta* at T-points 185 REAL(wp), DIMENSION(jpi,jpj) :: zp_deltastar_t ! P/delta* at T points 195 186 REAL(wp), DIMENSION(jpi,jpj) :: zzt, zet ! Viscosity pre-factors at T points 196 187 REAL(wp), DIMENSION(jpi,jpj) :: zef ! Viscosity pre-factor at F point … … 205 196 ! 206 197 REAL(wp), DIMENSION(jpi,jpj) :: zds ! shear 207 REAL(wp), DIMENSION(jpi,jpj) :: zs1, zs2, zs12 ! stress tensor components208 198 REAL(wp), DIMENSION(jpi,jpj) :: zsshdyn ! array used for the calculation of ice surface slope: 209 199 ! ! ocean surface (ssh_m) if ice is not embedded … … 225 215 !!! REAL(wp), DIMENSION(jpi,jpj) :: ztaux_base, ztauy_base ! ice-bottom stress at U-V points (landfast) 226 216 ! 227 REAL(wp), DIMENSION(jpi,jpj) :: zmsk01x, zmsk01y ! dummy arrays228 REAL(wp), DIMENSION(jpi,jpj) :: zmsk00x, zmsk00y ! mask for ice presence 217 REAL(wp), DIMENSION(jpi,jpj) :: zmsk01x, zmsk01y ! mask for lots of ice (1), little ice (0) 218 REAL(wp), DIMENSION(jpi,jpj) :: zmsk00x, zmsk00y ! mask for ice presence (1), no ice (0) 229 219 ! 230 220 REAL(wp), PARAMETER :: zepsi = 1.0e-20_wp ! tolerance parameter … … 232 222 REAL(wp), PARAMETER :: zamin = 0.001_wp ! ice concentration below which ice velocity becomes very small 233 223 !! --- diags 234 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zsig1, zsig2, zsig3 235 !! --- SIMIP diags 236 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s) 237 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_ymtrp_ice ! Y-component of ice mass transport (kg/s) 238 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_xmtrp_snw ! X-component of snow mass transport (kg/s) 239 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_ymtrp_snw ! Y-component of snow mass transport (kg/s) 240 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_xatrp ! X-component of area transport (m2/s) 241 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_yatrp ! Y-component of area transport (m2/s) 224 REAL(wp) :: zsig1, zsig2, zsig12 225 REAL(wp), DIMENSION(jpi,jpj) :: zs1, zs2, zs12, zs12f, zten_i ! stress tensor components 226 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zsig_I, zsig_II, zsig1_p, zsig2_p 227 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s, SIMIP) 228 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_ymtrp_ice ! Y-component of ice mass transport (kg/s, SIMIP) 229 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_xmtrp_snw ! X-component of snow mass transport (kg/s, SIMIP) 230 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_ymtrp_snw ! Y-component of snow mass transport (kg/s, SIMIP) 231 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_xatrp ! X-component of area transport (m2/s, SIMIP) 232 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_yatrp ! Y-component of area transport (m2/s, SIMIP) 242 233 243 234 !!---------------------------------------------------------------------------------------------------------------------- … … 269 260 ecc2 = rn_ecc * rn_ecc 270 261 z1_ecc2 = 1._wp / ecc2 271 272 ! Initialise stress tensor from previous time step 273 zs1 (:,:) = pstress1_i (:,:) 274 zs2 (:,:) = pstress2_i (:,:) 275 zs12(:,:) = pstress12_i(:,:) 276 262 277 263 ! Initialise convergence checks 278 264 IF( ln_rhg_chkcvg ) THEN … … 389 375 zspgV(ji,jj) = - zmassV * grav * ( zsshdyn(ji,jj+1) - zsshdyn(ji,jj) ) * r1_e2v(ji,jj) 390 376 391 ! masks377 ! Mask for ice presence (1) or absence (0) 392 378 zmsk00x(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) ) ! 0 if no ice 393 379 zmsk00y(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) ) ! 0 if no ice 394 380 395 ! switches381 ! Mask for lots of ice (1) or little ice (0) 396 382 IF( zmassU <= zmmin .AND. zaU(ji,jj) <= zamin ) THEN ; zmsk01x(ji,jj) = 0._wp 397 383 ELSE ; zmsk01x(ji,jj) = 1._wp ; ENDIF 398 384 IF( zmassV <= zmmin .AND. zaV(ji,jj) <= zamin ) THEN ; zmsk01y(ji,jj) = 0._wp 399 ELSE ; zmsk01y(ji,jj) = 1._wp ; ENDIF 385 ELSE ; zmsk01y(ji,jj) = 1._wp ; ENDIF 400 386 401 387 END DO … … 473 459 474 460 ! delta* at T points (following Lemieux and Dupont, GMD 2020) 475 zdeltat_star = MAX( zdeltat, rn_delta_min ) 476 477 IF ( ln_smooth_vp ) zdelta_star = zdeltat + rn_delta_min 478 461 zdeltastar_t(ji,jj) = zdeltat + rn_creepl 462 479 463 ! P/delta at T-points 480 zp_deltastar_t(ji,jj) = strength(ji,jj) / zdelta t_star464 zp_deltastar_t(ji,jj) = strength(ji,jj) / zdeltastar_t(ji,jj) 481 465 482 466 ! Temporary zzt and zet factors at T-points … … 486 470 END DO 487 471 END DO 488 472 489 473 CALL lbc_lnk_multi( 'icedyn_rhg_vp', zp_deltastar_t , 'T', 1. , zzt , 'T', 1., zet, 'T', 1. ) 490 474 … … 502 486 503 487 CALL lbc_lnk( 'icedyn_rhg_vp', zef, 'F', 1. ) 504 488 505 489 !--------------------------------------------------- 506 490 ! -- Ocean-ice drag and Coriolis RHS contributions … … 799 783 !--------------- 800 784 DO ji = 2, jpi - 1 785 801 786 u_ice(ji,jj) = zu_b(ji,jj) + rn_relaxu_vp * ( u_ice(ji,jj) - zu_b(ji,jj) ) 787 788 ! velocity masking for little-ice and no-ice cases 789 ! put 0.01 of ocean current, appropriate choice to avoid too much spreading of the ice edge 790 u_ice(ji,jj) = zmsk00x(ji,jj) & 791 & * ( zmsk01x(ji,jj) * u_ice(ji,jj) & 792 + ( 1._wp - zmsk01x(ji,jj) ) * u_oce(ji,jj) * 0.01_wp ) 793 802 794 END DO 803 795 … … 871 863 DO jj = 2, jpj - 1 872 864 v_ice(ji,jj) = zv_b(ji,jj) + rn_relaxv_vp * ( v_ice(ji,jj) - zv_b(ji,jj) ) 865 866 ! mask velocity for no-ice and little-ice cases 867 v_ice(ji,jj) = zmsk00y(ji,jj) & 868 & * ( zmsk01y(ji,jj) * v_ice(ji,jj) & 869 + ( 1._wp - zmsk01y(ji,jj) ) * v_oce(ji,jj) * 0.01_wp ) 870 873 871 END DO 874 872 … … 897 895 DO jj = 2, jpj - 1 898 896 DO ji = 2, jpi - 1 899 zuerr(ji,jj) = ABS ( ( u_ice(ji,jj) - zu_b(ji,jj) ) ) * umask(ji,jj,1) ! * zmask15 ---> MV TEST mask at 15% concentration897 zuerr(ji,jj) = ABS ( ( u_ice(ji,jj) - zu_b(ji,jj) ) ) * umask(ji,jj,1) 900 898 END DO 901 899 END DO … … 975 973 !------------------------------------------------------------------------------! 976 974 ! 977 ! OPT: subroutinize ?975 ! MV OPT: subroutinize ? 978 976 979 977 DO jj = 1, jpj - 1 … … 997 995 zdt2 = zdt * zdt 998 996 997 zten_i(ji,jj) = zdt 998 999 999 ! shear**2 at T points (doc eq. A16) 1000 1000 zds2 = ( zds(ji,jj ) * zds(ji,jj ) * e1e2f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e1e2f(ji-1,jj ) & … … 1013 1013 zdelta = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) 1014 1014 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0 1015 pdelta_i(ji,jj) = zdelta + rn_creepl * rswitch 1015 1016 !pdelta_i(ji,jj) = zdelta + rn_creepl * rswitch 1017 pdelta_i(ji,jj) = zdelta + rn_creepl 1016 1018 1017 1019 END DO … … 1019 1021 1020 1022 CALL lbc_lnk_multi( 'icedyn_rhg_vp', pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1. ) 1021 CALL lbc_lnk_multi( 'icedyn_rhg_vp', zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. ) 1022 1023 ! --- Store the stress tensor for the next time step --- ! 1024 ! MV OPT: Are they computed at the end of the tucking fime step ??? 1025 ! MV Is stress at previous time step needed for VP, normally no, because they equation is not tucking fime dependent!!! 1026 ! 1027 pstress1_i (:,:) = zs1 (:,:) 1028 pstress2_i (:,:) = zs2 (:,:) 1029 pstress12_i(:,:) = zs12(:,:) 1030 ! 1031 1023 1032 1024 !------------------------------------------------------------------------------! 1033 1025 ! … … 1035 1027 ! 1036 1028 !------------------------------------------------------------------------------! 1037 ! MV OPT: subroutinize 1038 ! 1039 ! --- Ice-ocean, ice-atm. & ice-oceanbottom(landfast) stresses --- ! 1029 ! 1030 ! MV OPT: subroutinize ? 1031 ! 1032 !---------------------------------- 1033 ! --- Recompute stresses if needed 1034 !---------------------------------- 1035 ! 1036 ! ---- Sea ice stresses at T-points 1037 IF ( iom_use('normstr') .OR. iom_use('sheastr') .OR. iom_use('intstrx') .OR. iom_use('intstry') ) THEN 1038 1039 DO jj = 2, jpj - 1 1040 DO ji = 2, jpi - 1 1041 zp_deltastar_t(ji,jj) = strength(ji,jj) / pdelta_i(ji,jj) 1042 zfac = zp_deltastar_t(ji,jj) 1043 zs1(ji,jj) = zfac * ( pdivu_i(ji,jj) - pdelta_i(ji,jj) ) 1044 zs2(ji,jj) = zfac * z1_ecc2 * zten_i(ji,jj) 1045 zs12(ji,jj) = zfac * z1_ecc2 * pshear_i(ji,jj) 1046 END DO 1047 END DO 1048 1049 CALL lbc_lnk_multi( 'icedyn_rhg_vp', zs1, 'T', 1., zs2, 'T', 1., zs12, 'T', 1. ) 1050 1051 ENDIF 1052 1053 ! ---- s12 at F-points 1054 IF ( iom_use('intstrx') .OR. iom_use('intstry') ) THEN 1055 1056 DO jj = 1, jpj - 1 1057 DO ji = 1, jpi - 1 1058 1059 ! P/delta* at F points 1060 zp_deltastar_f = 0.25_wp * ( zp_deltastar_t(ji,jj) + zp_deltastar_t(ji+1,jj) + zp_deltastar_t(ji,jj+1) + zp_deltastar_t(ji+1,jj+1) ) 1061 1062 ! s12 at F-points 1063 zs12f(ji,jj) = zp_deltastar_f * z1_ecc2 * zds(ji,jj) 1064 1065 END DO 1066 END DO 1067 1068 CALL lbc_lnk( 'icedyn_rhg_vp', zs12f, 'F', 1. ) 1069 1070 ENDIF 1071 ! 1072 !----------------------- 1073 ! --- Store diagnostics 1074 !----------------------- 1075 ! 1076 ! --- Ice-ocean, ice-atm. & ice-ocean bottom (landfast) stresses --- ! 1040 1077 IF( iom_use('utau_oi') .OR. iom_use('vtau_oi') .OR. iom_use('utau_ai') .OR. iom_use('vtau_ai') .OR. & 1041 1078 & iom_use('utau_bi') .OR. iom_use('vtau_bi') ) THEN 1042 1043 1079 ! 1044 1080 CALL lbc_lnk_multi( 'icedyn_rhg_vp', ztaux_oi, 'U', -1., ztauy_oi, 'V', -1., ztaux_ai, 'U', -1., ztauy_ai, 'V', -1., & 1045 & 1081 & ztaux_bi, 'U', -1., ztauy_bi, 'V', -1. ) 1046 1082 ! 1047 1083 CALL iom_put( 'utau_oi' , ztaux_oi * zmsk00 ) … … 1059 1095 IF( iom_use('icestr') ) CALL iom_put( 'icestr' , strength * zmsk00 ) ! strength 1060 1096 1061 ! --- Stress tensor --- ! 1062 IF( iom_use('isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') .OR. iom_use('normstr') .OR. iom_use('sheastr') ) THEN 1063 ! 1064 ALLOCATE( zsig1(jpi,jpj) , zsig2(jpi,jpj) , zsig3(jpi,jpj) ) 1097 ! --- Stress tensor invariants (SIMIP diags) --- ! 1098 IF( iom_use('normstr') .OR. iom_use('sheastr') ) THEN 1099 ! 1100 ! Stress tensor invariants (normal and shear stress N/m) - SIMIP diags. 1101 ! Definitions following Coon (1974) and Feltham (2008) 1102 ! 1103 ! sigma1, sigma2, sigma12 are useful (Hunke and Dukowicz MWR 2002, Bouillon et al., OM2013) 1104 ! however these are NOT stress tensor components, neither stress invariants, nor stress principal components 1105 ! 1106 ALLOCATE( zsig_I(jpi,jpj) , zsig_II(jpi,jpj) ) 1065 1107 ! 1066 1108 DO jj = 2, jpj - 1 1067 1109 DO ji = 2, jpi - 1 1068 1069 zdum1 = ( zmsk00(ji-1,jj) * pstress12_i(ji-1,jj) + zmsk00(ji ,jj-1) * pstress12_i(ji ,jj-1) + & ! stress12_i at T-point 1070 & zmsk00(ji ,jj) * pstress12_i(ji ,jj) + zmsk00(ji-1,jj-1) * pstress12_i(ji-1,jj-1) ) & 1071 & / MAX( 1._wp, zmsk00(ji-1,jj) + zmsk00(ji,jj-1) + zmsk00(ji,jj) + zmsk00(ji-1,jj-1) ) 1072 1073 zshear = SQRT( pstress2_i(ji,jj) * pstress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress 1074 1075 zdum2 = zmsk00(ji,jj) / MAX( 1._wp, strength(ji,jj) ) 1076 1077 zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) ) ! compressive stress, see Bouillon et al. 2015 1078 zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear ) ! shear stress 1079 zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) 1080 1110 ! Stress invariants 1111 zsig_I(ji,jj) = zs1(ji,jj) * 0.5_wp ! 1st invariant, aka average normal stress aka negative pressure 1112 zsig_II(ji,jj) = SQRT ( zs2(ji,jj) * zs2(ji,jj) * 0.25_wp + zs12(ji,jj) ) ! 2nd invariant, aka maximum shear stress 1081 1113 END DO 1082 1114 END DO 1083 1115 1084 CALL lbc_lnk_multi( 'icedyn_rhg_vp', zsig1, 'T', 1., zsig2, 'T', 1., zsig3, 'T', 1. ) 1085 ! 1086 CALL iom_put( 'isig1' , zsig1 ) 1087 CALL iom_put( 'isig2' , zsig2 ) 1088 CALL iom_put( 'isig3' , zsig3 ) 1089 ! 1090 ! Stress tensor invariants (normal and shear stress N/m) 1091 IF( iom_use('normstr') ) CALL iom_put( 'normstr' , ( zs1(:,:) + zs2(:,:) ) * zmsk00(:,:) ) ! Normal stress 1092 IF( iom_use('sheastr') ) CALL iom_put( 'sheastr' , SQRT( ( zs1(:,:) - zs2(:,:) )**2 + 4*zs12(:,:)**2 ) * zmsk00(:,:) ) ! Shear stress 1093 1094 DEALLOCATE( zsig1 , zsig2 , zsig3 ) 1116 CALL lbc_lnk_multi( 'icedyn_rhg_vp', zsig_I, 'T', 1., zsig_II, 'T', 1.) 1117 1118 IF( iom_use('normstr') ) CALL iom_put( 'normstr' , zsig_I(:,:) * zmsk00(:,:) ) ! Normal stress 1119 IF( iom_use('sheastr') ) CALL iom_put( 'sheastr' , zsig_II(:,:) * zmsk00(:,:) ) ! Maximum shear stress 1120 1121 DEALLOCATE ( zsig_I, zsig_II ) 1122 1123 ENDIF 1124 1125 ! --- Normalized stress tensor principal components --- ! 1126 ! These are used to plot the normalized yield curve (Lemieux & Dupont, GMD 2020) 1127 ! To plot the yield curve and evaluate physical convergence, they have two recommendations 1128 ! Recommendation 1 : Use ice strength, not replacement pressure 1129 ! Recommendation 2 : Need to use deformations at PREVIOUS iterate for viscosities (see p. 1765) 1130 ! R2 means we need to recompute stresses 1131 IF( iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN 1132 ! 1133 ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) ) 1134 ! 1135 DO jj = 2, jpj - 1 1136 DO ji = 2, jpi - 1 1137 ! Ice stresses computed with **viscosities** (delta, p/delta) at **previous** iterates 1138 ! and **deformations** at current iterates 1139 ! following Lemieux & Dupont (2020) 1140 zfac = zp_deltastar_t(ji,jj) 1141 zsig1 = zfac * ( pdivu_i(ji,jj) - zdeltastar_t(ji,jj) ) 1142 zsig2 = zfac * z1_ecc2 * zten_i(ji,jj) 1143 zsig12 = zfac * z1_ecc2 * pshear_i(ji,jj) 1144 1145 ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008), T-point 1146 zsig_I(ji,jj) = zsig1 * 0.5_wp ! 1st invariant 1147 zsig_II(ji,jj) = SQRT ( zsig2 * zsig2 * 0.25_wp + zsig12 ) ! 2nd invariant 1148 1149 ! Normalized principal stresses (used to display the ellipse) 1150 z1_strength = 1._wp / strength(ji,jj) 1151 zsig1_p(ji,jj) = ( zsig_I(ji,jj) + zsig_II(ji,jj) ) * z1_strength 1152 zsig2_p(ji,jj) = ( zsig_II(ji,jj) - zsig_II(ji,jj) ) * z1_strength 1153 END DO 1154 END DO 1155 ! 1156 CALL lbc_lnk_multi( 'icedyn_rhg_vp', zsig1_p, 'T', 1., zsig2_p, 'T', 1.) 1157 ! 1158 ! 1159 CALL iom_put( 'sig1_pnorm' , zsig1_p ) 1160 CALL iom_put( 'sig2_pnorm' , zsig2_p ) 1161 1162 DEALLOCATE( zsig1_p , zsig2_p ) 1095 1163 1096 1164 ENDIF … … 1098 1166 ! --- SIMIP, terms of tendency for momentum equation --- ! 1099 1167 IF( iom_use('dssh_dx') .OR. iom_use('dssh_dy') .OR. & 1100 & iom_use('corstrx') .OR. iom_use('corstry') .OR. iom_use('intstrx') .OR. iom_use('intstry') ) THEN 1101 ! 1102 !!!!!!!!! ATTENTION LA FORCE INTERNE DOIT ETTRE RECALCIULEEE ICCI !!!!!!!!!!!!!!!! 1168 & iom_use('corstrx') .OR. iom_use('corstry') ) THEN 1169 ! 1103 1170 CALL lbc_lnk_multi( 'icedyn_rhg_vp', zspgU, 'U', -1., zspgV, 'V', -1., & 1104 & zCorU, 'U', -1., zCorV, 'V', -1., zfU, 'U', -1., zfV, 'V', -1. )1171 & zCorU, 'U', -1., zCorV, 'V', -1. ) 1105 1172 1106 1173 CALL iom_put( 'dssh_dx' , zspgU * zmsk00 ) ! Sea-surface tilt term in force balance (x) … … 1108 1175 CALL iom_put( 'corstrx' , zCorU * zmsk00 ) ! Coriolis force term in force balance (x) 1109 1176 CALL iom_put( 'corstry' , zCorV * zmsk00 ) ! Coriolis force term in force balance (y) 1177 ENDIF 1178 1179 IF ( iom_use('intstrx') .OR. iom_use('intstry') ) THEN 1180 1181 1182 ! Recalculate internal forces (divergence of stress tensor) 1183 DO jj = 2, jpj - 1 1184 DO ji = 2, jpi - 1 1185 zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj) & 1186 & + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj) & 1187 & ) * r1_e2u(ji,jj) & 1188 & + ( zs12f(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12f(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) & 1189 & ) * 2._wp * r1_e1u(ji,jj) & 1190 & ) * r1_e1e2u(ji,jj) 1191 zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj) & 1192 & - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj) & 1193 & ) * r1_e1v(ji,jj) & 1194 & + ( zs12f(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12f(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) & 1195 & ) * 2._wp * r1_e2v(ji,jj) & 1196 & ) * r1_e1e2v(ji,jj) 1197 END DO 1198 END DO 1199 1200 CALL lbc_lnk_multi( 'icedyn_rhg_vp', zfU, 'U', -1., zfV, 'V', -1. ) 1201 1110 1202 CALL iom_put( 'intstrx' , zfU * zmsk00 ) ! Internal force term in force balance (x) 1111 1203 CALL iom_put( 'intstry' , zfV * zmsk00 ) ! Internal force term in force balance (y) 1204 1112 1205 ENDIF 1113 1206 1207 ! --- Ice & snow mass and ice area transports 1114 1208 IF( iom_use('xmtrpice') .OR. iom_use('ymtrpice') .OR. & 1115 1209 & iom_use('xmtrpsnw') .OR. iom_use('ymtrpsnw') .OR. iom_use('xatrp') .OR. iom_use('yatrp') ) THEN … … 1137 1231 1138 1232 CALL lbc_lnk_multi( 'icedyn_rhg_vp', zdiag_xmtrp_ice, 'U', -1., zdiag_ymtrp_ice, 'V', -1., & 1139 & 1140 & 1233 & zdiag_xmtrp_snw, 'U', -1., zdiag_ymtrp_snw, 'V', -1., & 1234 & zdiag_xatrp , 'U', -1., zdiag_yatrp , 'V', -1. ) 1141 1235 1142 1236 CALL iom_put( 'xmtrpice' , zdiag_xmtrp_ice ) ! X-component of sea-ice mass transport (kg/s) … … 1172 1266 & prhsu, pAU, pBU, pCU, pDU, pEU, prhsv, pAV, pBV, pCV, pDV, pEV ) 1173 1267 1174 ! CALL rhg_cvg_vp( kt, jter, nn_nvp, u_ice, v_ice, zmt, zuerr_max, zverr_max, zglob_area, &1175 ! & zrhsu, zAU, zBU, zCU, zDU, zEU, zrhsv, zAV, zBV, zCV, zDV, zEV )1176 1268 !!---------------------------------------------------------------------- 1177 1269 !! *** ROUTINE rhg_cvg_vp *** … … 1302 1394 1303 1395 1304 1305 SUBROUTINE rhg_vp_rst( cdrw, kt )1306 !!---------------------------------------------------------------------1307 !! *** ROUTINE rhg_vp_rst ***1308 !!1309 !! ** Purpose : Read or write RHG file in restart file1310 !!1311 !! ** Method : use of IOM library1312 !!----------------------------------------------------------------------1313 CHARACTER(len=*) , INTENT(in) :: cdrw ! "READ"/"WRITE" flag1314 INTEGER, OPTIONAL, INTENT(in) :: kt ! ice time-step1315 !1316 INTEGER :: iter ! local integer1317 INTEGER :: id1, id2, id3 ! local integers1318 !!----------------------------------------------------------------------1319 !1320 IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialize1321 ! ! ---------------1322 IF( ln_rstart ) THEN !* Read the restart file1323 !1324 id1 = iom_varid( numrir, 'stress1_i' , ldstop = .FALSE. )1325 id2 = iom_varid( numrir, 'stress2_i' , ldstop = .FALSE. )1326 id3 = iom_varid( numrir, 'stress12_i', ldstop = .FALSE. )1327 !1328 IF( MIN( id1, id2, id3 ) > 0 ) THEN ! fields exist1329 CALL iom_get( numrir, jpdom_autoglo, 'stress1_i' , stress1_i )1330 CALL iom_get( numrir, jpdom_autoglo, 'stress2_i' , stress2_i )1331 CALL iom_get( numrir, jpdom_autoglo, 'stress12_i', stress12_i )1332 ELSE ! start rheology from rest1333 IF(lwp) WRITE(numout,*)1334 IF(lwp) WRITE(numout,*) ' ==>>> previous run without rheology, set stresses to 0'1335 stress1_i (:,:) = 0._wp1336 stress2_i (:,:) = 0._wp1337 stress12_i(:,:) = 0._wp1338 ENDIF1339 ELSE !* Start from rest1340 IF(lwp) WRITE(numout,*)1341 IF(lwp) WRITE(numout,*) ' ==>>> start from rest: set stresses to 0'1342 stress1_i (:,:) = 0._wp1343 stress2_i (:,:) = 0._wp1344 stress12_i(:,:) = 0._wp1345 ENDIF1346 !1347 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file1348 ! ! -------------------1349 IF(lwp) WRITE(numout,*) '---- rhg-rst ----'1350 iter = kt + nn_fsbc - 1 ! ice restarts are written at kt == nitrst - nn_fsbc + 11351 !1352 CALL iom_rstput( iter, nitrst, numriw, 'stress1_i' , stress1_i )1353 CALL iom_rstput( iter, nitrst, numriw, 'stress2_i' , stress2_i )1354 CALL iom_rstput( iter, nitrst, numriw, 'stress12_i', stress12_i )1355 !1356 ENDIF1357 !1358 END SUBROUTINE rhg_vp_rst1359 1360 1396 1361 1397 #else
Note: See TracChangeset
for help on using the changeset viewer.
