diff --git a/ICE_RHEO/EXPREF/README b/ICE_RHEO/EXPREF/README index 7c535ba..100a105 100644 --- a/ICE_RHEO/EXPREF/README +++ b/ICE_RHEO/EXPREF/README @@ -23,8 +23,6 @@ A range of experiments can be configured: - The ratio of the windstress component can be changed to cause different modes of failure (Rwind in usrdef_sbc.F90, default is -0.8) -It takes about 10 minutes to run a simulation on 16 processors - ---------- How to run ---------- diff --git a/ICE_RHEO/MY_SRC/icedyn_rhg_eap.F90 b/ICE_RHEO/MY_SRC/icedyn_rhg_eap.F90 index 115b6ff..da650a5 100644 --- a/ICE_RHEO/MY_SRC/icedyn_rhg_eap.F90 +++ b/ICE_RHEO/MY_SRC/icedyn_rhg_eap.F90 @@ -57,16 +57,15 @@ MODULE icedyn_rhg_eap INTEGER, PARAMETER :: na_yield = 21 REAL(wp), DIMENSION(nx_yield, ny_yield, na_yield) :: s11r, s12r, s22r, s11s, s12s, s22s - - !! * Substitutions -# include "do_loop_substitute.h90" -# include "domzgr_substitute.h90" + REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fimask ! mask at F points for the ice !! for convergence tests INTEGER :: ncvgid ! netcdf file id INTEGER :: nvarid ! netcdf variable id - REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: aimsk00 - REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: eap_res , aimsk15 + + !! * Substitutions +# include "do_loop_substitute.h90" +# include "domzgr_substitute.h90" !!---------------------------------------------------------------------- !! NEMO/ICE 4.0 , NEMO Consortium (2018) !! $Id: icedyn_rhg_eap.F90 11536 2019-09-11 13:54:18Z smasson $ @@ -154,7 +153,6 @@ SUBROUTINE ice_dyn_rhg_eap( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear REAL(wp) :: zalphar, zalphas ! for mechanical redistribution REAL(wp) :: zmresult11, zmresult12, z1dtevpkth, zp5kth, z1_dtevp_A ! for structure tensor evolution REAL(wp) :: zinvw ! for test case - ! REAL(wp), DIMENSION(jpi,jpj) :: zstress12tmp ! anisotropic stress tensor component for regridding REAL(wp), DIMENSION(jpi,jpj) :: zyield11, zyield22, zyield12 ! yield surface tensor for history @@ -181,9 +179,9 @@ SUBROUTINE ice_dyn_rhg_eap( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear REAL(wp), DIMENSION(jpi,jpj) :: ztaux_bi, ztauy_bi ! ice-OceanBottom stress at U-V points (landfast) REAL(wp), DIMENSION(jpi,jpj) :: ztaux_base, ztauy_base ! ice-bottom stress at U-V points (landfast) ! + REAL(wp), DIMENSION(jpi,jpj) :: zmsk, zmsk00, zmsk15 REAL(wp), DIMENSION(jpi,jpj) :: zmsk01x, zmsk01y ! dummy arrays REAL(wp), DIMENSION(jpi,jpj) :: zmsk00x, zmsk00y ! mask for ice presence - REAL(wp), DIMENSION(jpi,jpj) :: zfmask ! mask at F points for the ice REAL(wp), PARAMETER :: zepsi = 1.0e-20_wp ! tolerance parameter REAL(wp), PARAMETER :: zmmin = 1._wp ! ice mass (kg/m2) below which ice velocity becomes very small @@ -204,63 +202,45 @@ SUBROUTINE ice_dyn_rhg_eap( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_dyn_rhg_eap: EAP sea-ice rheology' ! - IF( kt == nit000 ) THEN - ! - ! for diagnostics - ALLOCATE( aimsk00(jpi,jpj) ) - ! for convergence tests - IF( nn_rhg_chkcvg > 0 ) ALLOCATE( eap_res(jpi,jpj), aimsk15(jpi,jpj) ) - ENDIF - ! + ! for diagnostics and convergence tests DO_2D( 1, 1, 1, 1 ) - aimsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice , 0 if no ice + zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice , 0 if no ice + zmsk (ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) ! 1 if ice , 0 if no ice END_2D IF( nn_rhg_chkcvg > 0 ) THEN DO_2D( 1, 1, 1, 1 ) - aimsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less + zmsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less END_2D ENDIF ! -!!gm for Clem: OPTIMIZATION: I think zfmask can be computed one for all at the initialization.... !------------------------------------------------------------------------------! ! 0) mask at F points for the ice !------------------------------------------------------------------------------! - ! ocean/land mask - DO_2D( 1, 0, 1, 0 ) - zfmask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1) - END_2D - CALL lbc_lnk( 'icedyn_rhg_eap', zfmask, 'F', 1._wp ) - - ! Lateral boundary conditions on velocity (modify zfmask) - DO_2D( 0, 0, 0, 0 ) - IF( zfmask(ji,jj) == 0._wp ) THEN - zfmask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(ji,jj,1), umask(ji,jj+1,1), & - & vmask(ji,jj,1), vmask(ji+1,jj,1) ) ) - ENDIF - END_2D - DO jj = 2, jpjm1 - IF( zfmask(1,jj) == 0._wp ) THEN - zfmask(1 ,jj) = rn_ishlat * MIN( 1._wp , MAX( vmask(2,jj,1), umask(1,jj+1,1), umask(1,jj,1) ) ) - ENDIF - IF( zfmask(jpi,jj) == 0._wp ) THEN - zfmask(jpi,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(jpi,jj+1,1), vmask(jpim1,jj,1), umask(jpi,jj-1,1) ) ) - ENDIF - END DO - DO ji = 2, jpim1 - IF( zfmask(ji,1) == 0._wp ) THEN - zfmask(ji, 1 ) = rn_ishlat * MIN( 1._wp , MAX( vmask(ji+1,1,1), umask(ji,2,1), vmask(ji,1,1) ) ) - ENDIF - IF( zfmask(ji,jpj) == 0._wp ) THEN - zfmask(ji,jpj) = rn_ishlat * MIN( 1._wp , MAX( vmask(ji+1,jpj,1), vmask(ji-1,jpj,1), umask(ji,jpjm1,1) ) ) + IF( kt == nit000 ) THEN + ! ocean/land mask + ALLOCATE( fimask(jpi,jpj) ) + IF( rn_ishlat == 0._wp ) THEN + DO_2D( 0, 0, 0, 0 ) + fimask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1) + END_2D + ELSE + DO_2D( 0, 0, 0, 0 ) + fimask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1) + ! Lateral boundary conditions on velocity (modify fimask) + IF( fimask(ji,jj) == 0._wp ) THEN + fimask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(ji,jj,1), umask(ji,jj+1,1), & + & vmask(ji,jj,1), vmask(ji+1,jj,1) ) ) + ENDIF + END_2D ENDIF - END DO - CALL lbc_lnk( 'icedyn_rhg_eap', zfmask, 'F', 1.0_wp ) + CALL lbc_lnk( 'icedyn_rhg_eap', fimask, 'F', 1.0_wp ) + ENDIF !------------------------------------------------------------------------------! ! 1) define some variables and initialize arrays !------------------------------------------------------------------------------! zrhoco = rho0 * rn_cio -!extra code for test case boundary conditions + !extra code for test case boundary conditions zinvw=1._wp/(zrhoco*0.5_wp) ! ecc2: square of yield ellipse eccenticrity @@ -404,7 +384,7 @@ SUBROUTINE ice_dyn_rhg_eap( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear ! shear at F points zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) & & + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & - & ) * r1_e1e2f(ji,jj) * zfmask(ji,jj) + & ) * r1_e1e2f(ji,jj) * fimask(ji,jj) END_2D @@ -428,7 +408,7 @@ SUBROUTINE ice_dyn_rhg_eap( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear zdt2 = zdt * zdt ! delta at T points - zdelta(ji,jj) = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 ) + zdelta(ji,jj) = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 ) * zmsk(ji,jj) ! zmsk is for reducing cpu END_2D CALL lbc_lnk( 'icedyn_rhg_eap', zdelta, 'T', 1.0_wp ) @@ -462,8 +442,8 @@ SUBROUTINE ice_dyn_rhg_eap( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear ! structure tensor update CALL calc_ffrac(zstressptmp, zstressmtmp, zstress12tmp(ji,jj), paniso_11(ji,jj), paniso_12(ji,jj), zmresult11, zmresult12) - paniso_11(ji,jj) = (paniso_11(ji,jj) + 0.5*2.e-5*zdtevp + zmresult11*zdtevp) / (1. + 2.e-5*zdtevp) ! implicit - paniso_12(ji,jj) = (paniso_12(ji,jj) + zmresult12*zdtevp) / (1. + 2.e-5*zdtevp) ! implicit + paniso_11(ji,jj) = zmsk00(ji,jj) * (paniso_11(ji,jj) + 0.5*2.e-5*zdtevp + zmresult11*zdtevp) / (1. + 2.e-5*zdtevp) ! implicit + paniso_12(ji,jj) = zmsk00(ji,jj) * (paniso_12(ji,jj) + zmresult12*zdtevp) / (1. + 2.e-5*zdtevp) ! implicit IF (jter == nn_nevp) THEN zyield11(ji,jj) = 0.5_wp * (zstressptmp + zstressmtmp) @@ -488,8 +468,8 @@ SUBROUTINE ice_dyn_rhg_eap( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear ENDIF ! stress at T points (zkt/=0 if landfast) - zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zstressptmp ) * z1_alph1 - zs2(ji,jj) = ( zs2(ji,jj) * zalph1 + zstressmtmp ) * z1_alph1 + zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zstressptmp ) * z1_alph1 * zmsk(ji,jj) ! zmsk is for reducing cpu + zs2(ji,jj) = ( zs2(ji,jj) * zalph1 + zstressmtmp ) * z1_alph1 * zmsk(ji,jj) ! zmsk is for reducing cpu END_2D CALL lbc_lnk( 'icedyn_rhg_eap', zstress12tmp, 'T', 1.0_wp , paniso_11, 'T', 1.0_wp , paniso_12, 'T', 1.0_wp) @@ -593,11 +573,10 @@ SUBROUTINE ice_dyn_rhg_eap( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear & ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin & ) * zmsk00y(ji,jj) ENDIF -!extra code for test case boundary conditions + !extra code for test case boundary conditions IF (mjg(jj)<25 .or. mjg(jj)>975 .or. mig(ji)<25 .or. mig(ji)>975) THEN v_ice(ji,jj) = zinvw*(ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)) END IF - END_2D CALL lbc_lnk( 'icedyn_rhg_eap', v_ice, 'V', -1.0_wp ) ! @@ -650,11 +629,10 @@ SUBROUTINE ice_dyn_rhg_eap( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin & ) * zmsk00x(ji,jj) ENDIF -!extra code for test case boundary conditions + !extra code for test case boundary conditions IF (mjg(jj)<25 .or. mjg(jj)>975 .or. mig(ji)<25 .or. mig(ji)>975) THEN u_ice(ji,jj) = zinvw*(ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)) END IF - END_2D CALL lbc_lnk( 'icedyn_rhg_eap', u_ice, 'U', -1.0_wp ) ! @@ -709,7 +687,7 @@ SUBROUTINE ice_dyn_rhg_eap( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin & ) * zmsk00x(ji,jj) ENDIF -!extra code for test case boundary conditions + !extra code for test case boundary conditions IF (mjg(jj)<25 .or. mjg(jj)>975 .or. mig(ji)<25 .or. mig(ji)>975) THEN u_ice(ji,jj) = zinvw*(ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)) END IF @@ -765,11 +743,11 @@ SUBROUTINE ice_dyn_rhg_eap( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear & ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin & ) * zmsk00y(ji,jj) ENDIF -!extra code for test case boundary conditions + !extra code for test case boundary conditions IF (mjg(jj)<25 .or. mjg(jj)>975 .or. mig(ji)<25 .or. mig(ji)>975) THEN v_ice(ji,jj) = zinvw*(ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)) END IF - END_2D + END_2D CALL lbc_lnk( 'icedyn_rhg_eap', v_ice, 'V', -1.0_wp ) ! #if defined key_agrif @@ -781,7 +759,7 @@ SUBROUTINE ice_dyn_rhg_eap( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear ENDIF ! convergence test - IF( nn_rhg_chkcvg == 2 ) CALL rhg_cvg_eap( kt, jter, nn_nevp, u_ice, v_ice, zu_ice, zv_ice ) + IF( nn_rhg_chkcvg == 2 ) CALL rhg_cvg_eap( kt, jter, nn_nevp, u_ice, v_ice, zu_ice, zv_ice, zmsk15 ) ! ! ! ==================== ! END DO ! end loop over jter ! @@ -798,7 +776,7 @@ SUBROUTINE ice_dyn_rhg_eap( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear ! shear at F points zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) & & + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & - & ) * r1_e1e2f(ji,jj) * zfmask(ji,jj) + & ) * r1_e1e2f(ji,jj) * fimask(ji,jj) END_2D @@ -818,15 +796,15 @@ SUBROUTINE ice_dyn_rhg_eap( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear & ) * 0.25_wp * r1_e1e2t(ji,jj) ! shear at T points - pshear_i(ji,jj) = SQRT( zdt2 + zds2 ) + pshear_i(ji,jj) = SQRT( zdt2 + 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) & & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) & - & ) * r1_e1e2t(ji,jj) + & ) * 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 ) ! delta + 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 @@ -851,19 +829,19 @@ SUBROUTINE ice_dyn_rhg_eap( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear CALL lbc_lnk( 'icedyn_rhg_eap', ztaux_oi, 'U', -1.0_wp, ztauy_oi, 'V', -1.0_wp, ztaux_ai, 'U', -1.0_wp, & & ztauy_ai, 'V', -1.0_wp, ztaux_bi, 'U', -1.0_wp, ztauy_bi, 'V', -1.0_wp ) ! - CALL iom_put( 'utau_oi' , ztaux_oi * aimsk00 ) - CALL iom_put( 'vtau_oi' , ztauy_oi * aimsk00 ) - CALL iom_put( 'utau_ai' , ztaux_ai * aimsk00 ) - CALL iom_put( 'vtau_ai' , ztauy_ai * aimsk00 ) - CALL iom_put( 'utau_bi' , ztaux_bi * aimsk00 ) - CALL iom_put( 'vtau_bi' , ztauy_bi * aimsk00 ) + CALL iom_put( 'utau_oi' , ztaux_oi * zmsk00 ) + CALL iom_put( 'vtau_oi' , ztauy_oi * zmsk00 ) + CALL iom_put( 'utau_ai' , ztaux_ai * zmsk00 ) + CALL iom_put( 'vtau_ai' , ztauy_ai * zmsk00 ) + CALL iom_put( 'utau_bi' , ztaux_bi * zmsk00 ) + CALL iom_put( 'vtau_bi' , ztauy_bi * zmsk00 ) ENDIF ! --- divergence, shear and strength --- ! - IF( iom_use('icediv') ) CALL iom_put( 'icediv' , pdivu_i * aimsk00 ) ! divergence - IF( iom_use('iceshe') ) CALL iom_put( 'iceshe' , pshear_i * aimsk00 ) ! shear - IF( iom_use('icedlt') ) CALL iom_put( 'icedlt' , pdelta_i * aimsk00 ) ! delta - IF( iom_use('icestr') ) CALL iom_put( 'icestr' , strength * aimsk00 ) ! strength + 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('icedlt') ) CALL iom_put( 'icedlt' , pdelta_i * zmsk00 ) ! delta + IF( iom_use('icestr') ) CALL iom_put( 'icestr' , strength * zmsk00 ) ! strength ! --- Stress tensor invariants (SIMIP diags) --- ! IF( iom_use('normstr') .OR. iom_use('sheastr') ) THEN @@ -882,14 +860,14 @@ SUBROUTINE ice_dyn_rhg_eap( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear zsig12 = zfac * z1_ecc2 * pshear_i(ji,jj) ! 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 ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) ) ! 2nd '' '', aka maximum shear stress + 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 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 (:,:) * aimsk00(:,:) ) ! Normal stress - IF( iom_use('sheastr') ) CALL iom_put( 'sheastr', zsig_II(:,:) * aimsk00(:,:) ) ! Maximum shear stress + 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 DEALLOCATE ( zsig_I, zsig_II ) @@ -914,8 +892,8 @@ SUBROUTINE ice_dyn_rhg_eap( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear zsig12 = zfac * z1_ecc2 * pshear_i(ji,jj) ! 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 ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) ) ! 2nd '' '', aka maximum shear stress + 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 ! Normalized principal stresses (used to display the ellipse) z1_strength = 1._wp / MAX( 1._wp, strength(ji,jj) ) @@ -923,8 +901,8 @@ SUBROUTINE ice_dyn_rhg_eap( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear zsig2_p(ji,jj) = ( zsig_I(ji,jj) - zsig_II(ji,jj) ) * z1_strength END_2D ! - CALL iom_put( 'sig1_pnorm' , zsig1_p ) - CALL iom_put( 'sig2_pnorm' , zsig2_p ) + IF( iom_use('sig1_pnorm') ) CALL iom_put( 'sig1_pnorm' , zsig1_p(:,:) * zmsk00(:,:) ) + IF( iom_use('sig2_pnorm') ) CALL iom_put( 'sig2_pnorm' , zsig2_p(:,:) * zmsk00(:,:) ) DEALLOCATE( zsig1_p , zsig2_p , zsig_I, zsig_II ) @@ -935,15 +913,15 @@ SUBROUTINE ice_dyn_rhg_eap( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear CALL lbc_lnk( 'icedyn_rhg_eap', zyield11, 'T', 1.0_wp, zyield22, 'T', 1.0_wp, zyield12, 'T', 1.0_wp ) - CALL iom_put( 'yield11', zyield11 * aimsk00 ) - CALL iom_put( 'yield22', zyield22 * aimsk00 ) - CALL iom_put( 'yield12', zyield12 * aimsk00 ) + CALL iom_put( 'yield11', zyield11 * zmsk00 ) + CALL iom_put( 'yield22', zyield22 * zmsk00 ) + CALL iom_put( 'yield12', zyield12 * zmsk00 ) ENDIF ! --- anisotropy tensor --- ! IF( iom_use('aniso') ) THEN CALL lbc_lnk( 'icedyn_rhg_eap', paniso_11, 'T', 1.0_wp ) - CALL iom_put( 'aniso' , paniso_11 * aimsk00 ) + CALL iom_put( 'aniso' , paniso_11 * zmsk00 ) ENDIF ! --- SIMIP --- ! @@ -954,12 +932,12 @@ SUBROUTINE ice_dyn_rhg_eap( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear & zCorU, 'U', -1.0_wp, zCorV, 'V', -1.0_wp, & & zfU, 'U', -1.0_wp, zfV, 'V', -1.0_wp ) - CALL iom_put( 'dssh_dx' , zspgU * aimsk00 ) ! Sea-surface tilt term in force balance (x) - CALL iom_put( 'dssh_dy' , zspgV * aimsk00 ) ! Sea-surface tilt term in force balance (y) - CALL iom_put( 'corstrx' , zCorU * aimsk00 ) ! Coriolis force term in force balance (x) - CALL iom_put( 'corstry' , zCorV * aimsk00 ) ! Coriolis force term in force balance (y) - CALL iom_put( 'intstrx' , zfU * aimsk00 ) ! Internal force term in force balance (x) - CALL iom_put( 'intstry' , zfV * aimsk00 ) ! Internal force term in force balance (y) + CALL iom_put( 'dssh_dx' , zspgU * zmsk00 ) ! Sea-surface tilt term in force balance (x) + CALL iom_put( 'dssh_dy' , zspgV * zmsk00 ) ! Sea-surface tilt term in force balance (y) + CALL iom_put( 'corstrx' , zCorU * zmsk00 ) ! Coriolis force term in force balance (x) + CALL iom_put( 'corstry' , zCorV * zmsk00 ) ! Coriolis force term in force balance (y) + CALL iom_put( 'intstrx' , zfU * zmsk00 ) ! Internal force term in force balance (x) + CALL iom_put( 'intstry' , zfV * zmsk00 ) ! Internal force term in force balance (y) ENDIF IF( iom_use('xmtrpice') .OR. iom_use('ymtrpice') .OR. & @@ -970,8 +948,8 @@ SUBROUTINE ice_dyn_rhg_eap( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear ! DO_2D( 0, 0, 0, 0 ) ! 2D ice mass, snow mass, area transport arrays (X, Y) - zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * aimsk00(ji,jj) - zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * aimsk00(ji,jj) + zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * zmsk00(ji,jj) + zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * zmsk00(ji,jj) zdiag_xmtrp_ice(ji,jj) = rhoi * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component zdiag_ymtrp_ice(ji,jj) = rhoi * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) ! '' Y- '' @@ -1005,10 +983,10 @@ SUBROUTINE ice_dyn_rhg_eap( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear IF( iom_use('uice_cvg') ) THEN IF( ln_aEVP ) THEN ! output: beta * ( u(t=nn_nevp) - u(t=nn_nevp-1) ) CALL iom_put( 'uice_cvg', MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * zbeta(:,:) * umask(:,:,1) , & - & ABS( v_ice(:,:) - zv_ice(:,:) ) * zbeta(:,:) * vmask(:,:,1) ) * aimsk15(:,:) ) + & ABS( v_ice(:,:) - zv_ice(:,:) ) * zbeta(:,:) * vmask(:,:,1) ) * zmsk15(:,:) ) ELSE ! output: nn_nevp * ( u(t=nn_nevp) - u(t=nn_nevp-1) ) CALL iom_put( 'uice_cvg', REAL( nn_nevp ) * MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * umask(:,:,1) , & - & ABS( v_ice(:,:) - zv_ice(:,:) ) * vmask(:,:,1) ) * aimsk15(:,:) ) + & ABS( v_ice(:,:) - zv_ice(:,:) ) * vmask(:,:,1) ) * zmsk15(:,:) ) ENDIF ENDIF ENDIF @@ -1016,7 +994,7 @@ SUBROUTINE ice_dyn_rhg_eap( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear END SUBROUTINE ice_dyn_rhg_eap - SUBROUTINE rhg_cvg_eap( kt, kiter, kitermax, pu, pv, pub, pvb ) + SUBROUTINE rhg_cvg_eap( kt, kiter, kitermax, pu, pv, pub, pvb, pmsk15 ) !!---------------------------------------------------------------------- !! *** ROUTINE rhg_cvg_eap *** !! @@ -1031,6 +1009,7 @@ SUBROUTINE rhg_cvg_eap( kt, kiter, kitermax, pu, pv, pub, pvb ) !!---------------------------------------------------------------------- INTEGER , INTENT(in) :: kt, kiter, kitermax ! ocean time-step index REAL(wp), DIMENSION(:,:), INTENT(in) :: pu, pv, pub, pvb ! now and before velocities + REAL(wp), DIMENSION(:,:), INTENT(in) :: pmsk15 !! INTEGER :: it, idtime, istatus INTEGER :: ji, jj ! dummy loop indices @@ -1059,22 +1038,20 @@ SUBROUTINE rhg_cvg_eap( kt, kiter, kitermax, pu, pv, pub, pvb ) ENDIF ! time - it = ( kt - 1 ) * kitermax + kiter + it = ( kt - nit000 ) * kitermax + kiter ! convergence IF( kiter == 1 ) THEN ! remove the first iteration for calculations of convergence (always very large) zresm = 0._wp ELSE - DO_2D( 1, 1, 1, 1 ) - eap_res(ji,jj) = MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), & - & ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * aimsk15(ji,jj) + zresm = 0._wp + DO_2D( 0, 0, 0, 0 ) ! cut of the boundary of the box (forced velocities) - IF (mjg(jj)<=30 .or. mjg(jj)>970 .or. mig(ji)<=30 .or. mig(ji)>970) THEN - eap_res(ji,jj) = 0._wp - END IF - END_2D - - zresm = MAXVAL( eap_res ) + IF (mjg0(jj)>30 .AND. mjg0(jj)<=970 .AND. mig0(ji)>30 .AND. mig0(ji)<=970) THEN + zresm = MAX( zresm, MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), & + & ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * pmsk15(ji,jj) ) + ENDIF + END_2D CALL mpp_max( 'icedyn_rhg_evp', zresm ) ! max over the global domain ENDIF @@ -1082,7 +1059,7 @@ SUBROUTINE rhg_cvg_eap( kt, kiter, kitermax, pu, pv, pub, pvb ) ! write variables istatus = NF90_PUT_VAR( ncvgid, nvarid, (/zresm/), (/it/), (/1/) ) ! close file - IF( kt == nitend ) istatus = NF90_CLOSE(ncvgid) + IF( kt == nitend - nn_fsbc + 1 .AND. kiter == kitermax ) istatus = NF90_CLOSE(ncvgid) ENDIF END SUBROUTINE rhg_cvg_eap diff --git a/ICE_RHEO/MY_SRC/icedyn_rhg_evp.F90 b/ICE_RHEO/MY_SRC/icedyn_rhg_evp.F90 index d632307..82e1432 100644 --- a/ICE_RHEO/MY_SRC/icedyn_rhg_evp.F90 +++ b/ICE_RHEO/MY_SRC/icedyn_rhg_evp.F90 @@ -5,9 +5,9 @@ MODULE icedyn_rhg_evp !!====================================================================== !! History : - ! 2007-03 (M.A. Morales Maqueda, S. Bouillon) Original code !! 3.0 ! 2008-03 (M. Vancoppenolle) adaptation to new model - !! - ! 2008-11 (M. Vancoppenolle, S. Bouillon, Y. Aksenov) add surface tilt in ice rheolohy + !! - ! 2008-11 (M. Vancoppenolle, S. Bouillon, Y. Aksenov) add surface tilt in ice rheolohy !! 3.3 ! 2009-05 (G.Garric) addition of the evp case - !! 3.4 ! 2011-01 (A. Porter) dynamical allocation + !! 3.4 ! 2011-01 (A. Porter) dynamical allocation !! 3.5 ! 2012-08 (R. Benshila) AGRIF !! 3.6 ! 2016-06 (C. Rousset) Rewriting + landfast ice + mEVP (Bouillon 2013) !! 3.7 ! 2017 (C. Rousset) add aEVP (Kimmritz 2016-2017) @@ -27,8 +27,8 @@ MODULE icedyn_rhg_evp USE ice ! sea-ice: ice variables USE icevar ! ice_var_sshdyn USE icedyn_rdgrft ! sea-ice: ice strength - USE bdy_oce , ONLY : ln_bdy - USE bdyice + USE bdy_oce , ONLY : ln_bdy + USE bdyice #if defined key_agrif USE agrif_ice_interp #endif @@ -47,17 +47,17 @@ MODULE icedyn_rhg_evp PUBLIC ice_dyn_rhg_evp ! called by icedyn_rhg.F90 PUBLIC rhg_evp_rst ! called by icedyn_rhg.F90 - !! * Substitutions -# include "do_loop_substitute.h90" -# include "domzgr_substitute.h90" - !! for convergence tests INTEGER :: ncvgid ! netcdf file id INTEGER :: nvarid ! netcdf variable id - REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zmsk00, zmsk15 + REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fimask ! mask at F points for the ice + + !! * Substitutions +# include "do_loop_substitute.h90" +# include "domzgr_substitute.h90" !!---------------------------------------------------------------------- !! NEMO/ICE 4.0 , NEMO Consortium (2018) - !! $Id: icedyn_rhg_evp.F90 13612 2020-10-14 17:18:19Z clem $ + !! $Id: icedyn_rhg_evp.F90 15550 2021-11-28 20:02:31Z clem $ !! Software governed by the CeCILL license (see ./LICENSE) !!---------------------------------------------------------------------- CONTAINS @@ -68,9 +68,9 @@ SUBROUTINE ice_dyn_rhg_evp( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear !! EVP-C-grid !! !! ** purpose : determines sea ice drift from wind stress, ice-ocean - !! stress and sea-surface slope. Ice-ice interaction is described by - !! a non-linear elasto-viscous-plastic (EVP) law including shear - !! strength and a bulk rheology (Hunke and Dukowicz, 2002). + !! stress and sea-surface slope. Ice-ice interaction is described by + !! a non-linear elasto-viscous-plastic (EVP) law including shear + !! strength and a bulk rheology (Hunke and Dukowicz, 2002). !! !! The points in the C-grid look like this, dear reader !! @@ -78,24 +78,24 @@ SUBROUTINE ice_dyn_rhg_evp( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear !! | !! | !! (ji-1,jj) | (ji,jj) - !! --------- + !! --------- !! | | !! | (ji,jj) |------(ji,jj) !! | | - !! --------- + !! --------- !! (ji-1,jj-1) (ji,jj-1) !! !! ** Inputs : - wind forcing (stress), oceanic currents !! ice total volume (vt_i) per unit area !! snow total volume (vt_s) per unit area !! - !! ** Action : - compute u_ice, v_ice : the components of the + !! ** Action : - compute u_ice, v_ice : the components of the !! sea-ice velocity vector !! - compute delta_i, shear_i, divu_i, which are inputs !! of the ice thickness distribution !! !! ** Steps : 0) compute mask at F point - !! 1) Compute ice snow mass, ice strength + !! 1) Compute ice snow mass, ice strength !! 2) Compute wind, oceanic stresses, mass terms and !! coriolis terms of the momentum equation !! 3) Solve the momentum equation (iterative procedure) @@ -133,10 +133,8 @@ SUBROUTINE ice_dyn_rhg_evp( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear REAL(wp) :: zkt ! isotropic tensile strength for landfast ice REAL(wp) :: zvCr ! critical ice volume above which ice is landfast ! - REAL(wp) :: zintb, zintn ! dummy argument REAL(wp) :: zfac_x, zfac_y - REAL(wp) :: zshear, zdum1, zdum2 - REAL(wp) :: zinvw ! for test case + REAL(wp) :: zinvw ! for test case ! REAL(wp), DIMENSION(jpi,jpj) :: zdelta, zp_delt ! delta and P/delta at T points @@ -153,7 +151,7 @@ SUBROUTINE ice_dyn_rhg_evp( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear 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 - ! ! ice bottom surface if ice is embedded + ! ! ice bottom surface if ice is embedded REAL(wp), DIMENSION(jpi,jpj) :: zfU , zfV ! internal stresses REAL(wp), DIMENSION(jpi,jpj) :: zspgU, zspgV ! surface pressure gradient at U/V points REAL(wp), DIMENSION(jpi,jpj) :: zCorU, zCorV ! Coriolis stress array @@ -162,9 +160,9 @@ SUBROUTINE ice_dyn_rhg_evp( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear REAL(wp), DIMENSION(jpi,jpj) :: ztaux_bi, ztauy_bi ! ice-OceanBottom stress at U-V points (landfast) REAL(wp), DIMENSION(jpi,jpj) :: ztaux_base, ztauy_base ! ice-bottom stress at U-V points (landfast) ! + REAL(wp), DIMENSION(jpi,jpj) :: zmsk, zmsk00, zmsk15 REAL(wp), DIMENSION(jpi,jpj) :: zmsk01x, zmsk01y ! dummy arrays REAL(wp), DIMENSION(jpi,jpj) :: zmsk00x, zmsk00y ! mask for ice presence - REAL(wp), DIMENSION(jpi,jpj) :: zfmask ! mask at F points for the ice REAL(wp), PARAMETER :: zepsi = 1.0e-20_wp ! tolerance parameter REAL(wp), PARAMETER :: zmmin = 1._wp ! ice mass (kg/m2) below which ice velocity becomes very small @@ -173,65 +171,61 @@ SUBROUTINE ice_dyn_rhg_evp( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear REAL(wp), DIMENSION(jpi,jpj) :: zu_ice, zv_ice !! --- diags REAL(wp) :: zsig1, zsig2, zsig12, zfac, z1_strength - REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zsig_I, zsig_II, zsig1_p, zsig2_p + REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zsig_I, zsig_II, zsig1_p, zsig2_p !! --- SIMIP diags REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s) REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_ymtrp_ice ! Y-component of ice mass transport (kg/s) REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_xmtrp_snw ! X-component of snow mass transport (kg/s) REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_ymtrp_snw ! Y-component of snow mass transport (kg/s) REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_xatrp ! X-component of area transport (m2/s) - REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_yatrp ! Y-component of area transport (m2/s) + REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_yatrp ! Y-component of area transport (m2/s) + !! -- advect fields at the rheology time step for the calculation of strength + !! it seems that convergence is worse when ll_advups=true. So it is not really a good idea + LOGICAL :: ll_advups = .FALSE. + REAL(wp) :: zdt_ups + REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: za_i_ups, zv_i_ups ! tracers advected upstream !!------------------------------------------------------------------- IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_dyn_rhg_evp: EVP sea-ice rheology' ! ! for diagnostics and convergence tests - ALLOCATE( zmsk00(jpi,jpj), zmsk15(jpi,jpj) ) - DO_2D( 1, 1, 1, 1 ) + DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice , 0 if no ice - zmsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less + zmsk (ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) ! 1 if ice , 0 if no ice END_2D + IF( nn_rhg_chkcvg > 0 ) THEN + DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) + zmsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less + END_2D + ENDIF ! - !!gm for Clem: OPTIMIZATION: I think zfmask can be computed one for all at the initialization.... !------------------------------------------------------------------------------! ! 0) mask at F points for the ice !------------------------------------------------------------------------------! - ! ocean/land mask - DO_2D( 1, 0, 1, 0 ) - zfmask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1) - END_2D - CALL lbc_lnk( 'icedyn_rhg_evp', zfmask, 'F', 1._wp ) - - ! Lateral boundary conditions on velocity (modify zfmask) - DO_2D( 0, 0, 0, 0 ) - IF( zfmask(ji,jj) == 0._wp ) THEN - zfmask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(ji,jj,1), umask(ji,jj+1,1), & - & vmask(ji,jj,1), vmask(ji+1,jj,1) ) ) - ENDIF - END_2D - DO jj = 2, jpjm1 - IF( zfmask(1,jj) == 0._wp ) THEN - zfmask(1 ,jj) = rn_ishlat * MIN( 1._wp , MAX( vmask(2,jj,1), umask(1,jj+1,1), umask(1,jj,1) ) ) - ENDIF - IF( zfmask(jpi,jj) == 0._wp ) THEN - zfmask(jpi,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(jpi,jj+1,1), vmask(jpim1,jj,1), umask(jpi,jj-1,1) ) ) - ENDIF - END DO - DO ji = 2, jpim1 - IF( zfmask(ji,1) == 0._wp ) THEN - zfmask(ji, 1 ) = rn_ishlat * MIN( 1._wp , MAX( vmask(ji+1,1,1), umask(ji,2,1), vmask(ji,1,1) ) ) - ENDIF - IF( zfmask(ji,jpj) == 0._wp ) THEN - zfmask(ji,jpj) = rn_ishlat * MIN( 1._wp , MAX( vmask(ji+1,jpj,1), vmask(ji-1,jpj,1), umask(ji,jpjm1,1) ) ) + IF( kt == nit000 ) THEN + ! ocean/land mask + ALLOCATE( fimask(jpi,jpj) ) + IF( rn_ishlat == 0._wp ) THEN + DO_2D( 0, 0, 0, 0 ) + fimask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1) + END_2D + ELSE + DO_2D( 0, 0, 0, 0 ) + fimask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1) + ! Lateral boundary conditions on velocity (modify fimask) + IF( fimask(ji,jj) == 0._wp ) THEN + fimask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(ji,jj,1), umask(ji,jj+1,1), & + & vmask(ji,jj,1), vmask(ji+1,jj,1) ) ) + ENDIF + END_2D ENDIF - END DO - CALL lbc_lnk( 'icedyn_rhg_evp', zfmask, 'F', 1._wp ) - + CALL lbc_lnk( 'icedyn_rhg_evp', fimask, 'F', 1._wp ) + ENDIF !------------------------------------------------------------------------------! ! 1) define some variables and initialize arrays !------------------------------------------------------------------------------! - zrhoco = rho0 * rn_cio -!extra code for test case boundary conditions + zrhoco = rho0 * rn_cio + !extra code for test case boundary conditions zinvw=1._wp/(zrhoco*0.5_wp) ! ecc2: square of yield ellipse eccenticrity @@ -247,13 +241,13 @@ SUBROUTINE ice_dyn_rhg_evp( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear z1_alph1 = 1._wp / ( zalph1 + 1._wp ) z1_alph2 = 1._wp / ( zalph2 + 1._wp ) ELSE - zdtevp = rdt_ice + zdtevp = rDt_ice ! zalpha parameters set later on adaptatively ENDIF z1_dtevp = 1._wp / zdtevp - - ! Initialise stress tensor - zs1 (:,:) = pstress1_i (:,:) + + ! Initialise stress tensor + zs1 (:,:) = pstress1_i (:,:) zs2 (:,:) = pstress2_i (:,:) zs12(:,:) = pstress12_i(:,:) @@ -273,7 +267,13 @@ SUBROUTINE ice_dyn_rhg_evp( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear ! non-embedded sea ice: use ocean surface for slope calculation zsshdyn(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b) - DO_2D( 0, 0, 0, 0 ) + DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) + zm1 = ( rhos * vt_s(ji,jj) + rhoi * vt_i(ji,jj) ) ! Ice/snow mass at U-V points + zmf (ji,jj) = zm1 * ff_t(ji,jj) ! Coriolis at T points (m*f) + zdt_m(ji,jj) = zdtevp / MAX( zm1, zmmin ) ! dt/m at T points (for alpha and beta coefficients) + END_2D + + DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! ice fraction at U-V points zaU(ji,jj) = 0.5_wp * ( at_i(ji,jj) * e1e2t(ji,jj) + at_i(ji+1,jj) * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1) @@ -287,19 +287,14 @@ SUBROUTINE ice_dyn_rhg_evp( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear zmassV = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1) ! Ocean currents at U-V points - v_oceU(ji,jj) = 0.25_wp * ( v_oce(ji,jj) + v_oce(ji,jj-1) + v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * umask(ji,jj,1) - u_oceV(ji,jj) = 0.25_wp * ( u_oce(ji,jj) + u_oce(ji-1,jj) + u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * vmask(ji,jj,1) + ! (brackets added to fix the order of floating point operations for halo 1 - halo 2 compatibility) + v_oceU(ji,jj) = 0.25_wp * ( (v_oce(ji,jj) + v_oce(ji,jj-1)) + (v_oce(ji+1,jj) + v_oce(ji+1,jj-1)) ) * umask(ji,jj,1) + u_oceV(ji,jj) = 0.25_wp * ( (u_oce(ji,jj) + u_oce(ji-1,jj)) + (u_oce(ji,jj+1) + u_oce(ji-1,jj+1)) ) * vmask(ji,jj,1) - ! Coriolis at T points (m*f) - zmf(ji,jj) = zm1 * ff_t(ji,jj) - - ! dt/m at T points (for alpha and beta coefficients) - zdt_m(ji,jj) = zdtevp / MAX( zm1, zmmin ) - ! m/dt zmU_t(ji,jj) = zmassU * z1_dtevp zmV_t(ji,jj) = zmassV * z1_dtevp - + ! Drag ice-atm. ztaux_ai(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj) ztauy_ai(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj) @@ -319,29 +314,28 @@ SUBROUTINE ice_dyn_rhg_evp( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear ELSE ; zmsk01y(ji,jj) = 1._wp ; ENDIF END_2D - CALL lbc_lnk( 'icedyn_rhg_evp', zmf, 'T', 1.0_wp, zdt_m, 'T', 1.0_wp ) ! ! !== Landfast ice parameterization ==! ! IF( ln_landfast_L16 ) THEN !-- Lemieux 2016 - DO_2D( 0, 0, 0, 0 ) + DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! ice thickness at U-V points zvU = 0.5_wp * ( vt_i(ji,jj) * e1e2t(ji,jj) + vt_i(ji+1,jj) * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1) zvV = 0.5_wp * ( vt_i(ji,jj) * e1e2t(ji,jj) + vt_i(ji,jj+1) * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1) ! ice-bottom stress at U points - zvCr = zaU(ji,jj) * rn_lf_depfra * hu(ji,jj,Kmm) + zvCr = zaU(ji,jj) * rn_lf_depfra * hu(ji,jj,Kmm) * ( 1._wp - icb_mask(ji,jj) ) ! if grounded icebergs are read: ocean depth = 0 ztaux_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) ) ! ice-bottom stress at V points - zvCr = zaV(ji,jj) * rn_lf_depfra * hv(ji,jj,Kmm) + zvCr = zaV(ji,jj) * rn_lf_depfra * hv(ji,jj,Kmm) * ( 1._wp - icb_mask(ji,jj) ) ! if grounded icebergs are read: ocean depth = 0 ztauy_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) ) ! ice_bottom stress at T points - zvCr = at_i(ji,jj) * rn_lf_depfra * ht(ji,jj) + zvCr = at_i(ji,jj) * rn_lf_depfra * ht(ji,jj) * ( 1._wp - icb_mask(ji,jj) ) ! if grounded icebergs are read: ocean depth = 0 tau_icebfr(ji,jj) = - rn_lf_bfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) ) END_2D CALL lbc_lnk( 'icedyn_rhg_evp', tau_icebfr(:,:), 'T', 1.0_wp ) ! ELSE !-- no landfast - DO_2D( 0, 0, 0, 0 ) + DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ztaux_base(ji,jj) = 0._wp ztauy_base(ji,jj) = 0._wp END_2D @@ -353,7 +347,7 @@ SUBROUTINE ice_dyn_rhg_evp( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear ! ! ! ==================== ! DO jter = 1 , nn_nevp ! loop over jter ! - ! ! ==================== ! + ! ! ==================== ! l_full_nf_update = jter == nn_nevp ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1 ! ! convergence test @@ -365,12 +359,12 @@ SUBROUTINE ice_dyn_rhg_evp( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear ENDIF ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- ! - DO_2D( 1, 0, 1, 0 ) + DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) ! shear at F points zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) & & + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & - & ) * r1_e1e2f(ji,jj) * zfmask(ji,jj) + & ) * r1_e1e2f(ji,jj) * fimask(ji,jj) END_2D @@ -380,42 +374,42 @@ SUBROUTINE ice_dyn_rhg_evp( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear zds2 = ( zds(ji,jj ) * zds(ji,jj ) * e1e2f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e1e2f(ji-1,jj ) & & + 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) - + ! divergence at T points zdiv = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) & & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) & & ) * r1_e1e2t(ji,jj) zdiv2 = zdiv * zdiv - + ! tension at T points zdt = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) & & - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) & & ) * r1_e1e2t(ji,jj) zdt2 = zdt * zdt - + ! delta at T points - zdelta(ji,jj) = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 ) + zdelta(ji,jj) = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 ) * zmsk(ji,jj) ! zmsk is for reducing cpu - END_2D - CALL lbc_lnk( 'icedyn_rhg_evp', zdelta, 'T', 1.0_wp ) + ! P/delta at T points + zp_delt(ji,jj) = strength(ji,jj) / ( zdelta(ji,jj) + rn_creepl ) * zmsk(ji,jj) ! zmsk is for reducing cpu - ! P/delta at T points - DO_2D( 1, 1, 1, 1 ) - zp_delt(ji,jj) = strength(ji,jj) / ( zdelta(ji,jj) + rn_creepl ) END_2D + CALL lbc_lnk( 'icedyn_rhg_evp', zdelta, 'T', 1.0_wp, zp_delt, 'T', 1.0_wp ) - DO_2D( 0, 1, 0, 1 ) ! loop ends at jpi,jpj so that no lbc_lnk are needed for zs1 and zs2 + ! + DO_2D( nn_hls-1, nn_hls, nn_hls-1, nn_hls ) ! loop ends at jpi,jpj so that no lbc_lnk are needed for zs1 and zs2 ! divergence at T points (duplication to avoid communications) - zdiv = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) & - & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) & + ! (brackets added to fix the order of floating point operations for halo 1 - halo 2 compatibility) + zdiv = ( (e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)) & + & + (e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)) & & ) * r1_e1e2t(ji,jj) - + ! tension at T points (duplication to avoid communications) zdt = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) & & - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) & & ) * r1_e1e2t(ji,jj) - + ! alpha for aEVP ! gamma = 0.5*P/(delta+creepl) * (c*pi)**2/Area * dt/m ! alpha = beta = sqrt(4*gamma) @@ -430,21 +424,23 @@ SUBROUTINE ice_dyn_rhg_evp( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear ! zalph1 = zalph1 - 1._wp ! zalph2 = zalph1 ENDIF - + ! stress at T points (zkt/=0 if landfast) - zs1(ji,jj) = ( zs1(ji,jj)*zalph1 + zp_delt(ji,jj) * ( zdiv*(1._wp + zkt) - zdelta(ji,jj)*(1._wp - zkt) ) ) * z1_alph1 - zs2(ji,jj) = ( zs2(ji,jj)*zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2 - + zs1(ji,jj) = ( zs1(ji,jj)*zalph1 + zp_delt(ji,jj) * ( zdiv*(1._wp + zkt) - zdelta(ji,jj)*(1._wp - zkt) ) ) & + & * z1_alph1 * zmsk(ji,jj) ! zmsk is for reducing cpu + zs2(ji,jj) = ( zs2(ji,jj)*zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) & + & * z1_alph2 * zmsk(ji,jj) ! zmsk is for reducing cpu + END_2D ! Save beta at T-points for further computations IF( ln_aEVP ) THEN - DO_2D( 1, 1, 1, 1 ) + DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) zbeta(ji,jj) = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) ) END_2D ENDIF - - DO_2D( 1, 0, 1, 0 ) + + DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) ! alpha for aEVP IF( ln_aEVP ) THEN @@ -454,36 +450,39 @@ SUBROUTINE ice_dyn_rhg_evp( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear ! z1_alph2 = 1._wp / zalph2 ! zalph2 = zalph2 - 1._wp ENDIF - + ! P/delta at F points - zp_delf = 0.25_wp * ( zp_delt(ji,jj) + zp_delt(ji+1,jj) + zp_delt(ji,jj+1) + zp_delt(ji+1,jj+1) ) - + ! (brackets added to fix the order of floating point operations for halo 1 - halo 2 compatibility) + zp_delf = 0.25_wp * ( (zp_delt(ji,jj) + zp_delt(ji+1,jj)) + (zp_delt(ji,jj+1) + zp_delt(ji+1,jj+1)) ) + ! stress at F points (zkt/=0 if landfast) - zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 * (1._wp + zkt) ) * 0.5_wp ) * z1_alph2 + zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 * (1._wp + zkt) ) * 0.5_wp ) & + & * z1_alph2 END_2D ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- ! - DO_2D( 0, 0, 0, 0 ) + ! (brackets added to fix the order of floating point operations for halo 1 - halo 2 compatibility) + DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! !--- U points - zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj) & + zfU(ji,jj) = 0.5_wp * ( (( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj) & & + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj) & - & ) * r1_e2u(ji,jj) & + & ) * r1_e2u(ji,jj)) & & + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) & & ) * 2._wp * r1_e1u(ji,jj) & & ) * r1_e1e2u(ji,jj) ! ! !--- V points - zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj) & + zfV(ji,jj) = 0.5_wp * ( (( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj) & & - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj) & - & ) * r1_e1v(ji,jj) & + & ) * r1_e1v(ji,jj)) & & + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) & & ) * 2._wp * r1_e2v(ji,jj) & & ) * r1_e1e2v(ji,jj) ! ! !--- ice currents at U-V point - v_iceU(ji,jj) = 0.25_wp * ( v_ice(ji,jj) + v_ice(ji,jj-1) + v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * umask(ji,jj,1) - u_iceV(ji,jj) = 0.25_wp * ( u_ice(ji,jj) + u_ice(ji-1,jj) + u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * vmask(ji,jj,1) + v_iceU(ji,jj) = 0.25_wp * ( (v_ice(ji,jj) + v_ice(ji,jj-1)) + (v_ice(ji+1,jj) + v_ice(ji+1,jj-1)) ) * umask(ji,jj,1) + u_iceV(ji,jj) = 0.25_wp * ( (u_ice(ji,jj) + u_ice(ji-1,jj)) + (u_ice(ji,jj+1) + u_ice(ji-1,jj+1)) ) * vmask(ji,jj,1) ! END_2D ! @@ -492,7 +491,7 @@ SUBROUTINE ice_dyn_rhg_evp( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear ! Bouillon et al. 2009 (eq 34-35) => stable IF( MOD(jter,2) == 0 ) THEN ! even iterations ! - DO_2D( 0, 0, 0, 0 ) + DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! !--- tau_io/(v_oce - v_ice) zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) ) & & + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) ) @@ -522,7 +521,7 @@ SUBROUTINE ice_dyn_rhg_evp( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * ( zbetav * v_ice(ji,jj) + v_ice_b(ji,jj) ) & ! previous velocity & + zRHS + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) & ) / MAX( zepsi, zmV_t(ji,jj) * ( zbetav + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast - & + ( 1._wp - rswitch ) * ( v_ice_b(ji,jj) & + & + ( 1._wp - rswitch ) * ( v_ice_b(ji,jj) & & + v_ice (ji,jj) * MAX( 0._wp, zbetav - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 & ) / ( zbetav + 1._wp ) & & ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin @@ -535,18 +534,12 @@ SUBROUTINE ice_dyn_rhg_evp( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear & ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin & ) * zmsk00y(ji,jj) ENDIF -!extra code for test case boundary conditions + !extra code for test case boundary conditions IF (mjg(jj)<25 .or. mjg(jj)>975 .or. mig(ji)<25 .or. mig(ji)>975) THEN v_ice(ji,jj) = zinvw*(ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)) END IF - END_2D - CALL lbc_lnk( 'icedyn_rhg_evp', v_ice, 'V', -1.0_wp ) - ! -#if defined key_agrif -!! CALL agrif_interp_ice( 'V', jter, nn_nevp ) - CALL agrif_interp_ice( 'V' ) -#endif - IF( ln_bdy ) CALL bdy_ice_dyn( 'V' ) + END_2D + IF( nn_hls == 1 ) CALL lbc_lnk( 'icedyn_rhg_evp', v_ice, 'V', -1.0_wp ) ! DO_2D( 0, 0, 0, 0 ) ! !--- tau_io/(u_oce - u_ice) @@ -581,32 +574,28 @@ SUBROUTINE ice_dyn_rhg_evp( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear & + ( 1._wp - rswitch ) * ( u_ice_b(ji,jj) & & + u_ice (ji,jj) * MAX( 0._wp, zbetau - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 & ) / ( zbetau + 1._wp ) & - & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin + & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin & ) * zmsk00x(ji,jj) ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) & ! previous velocity & + zRHS + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) & ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 - & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin + & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin & ) * zmsk00x(ji,jj) ENDIF -!extra code for test case boundary conditions + !extra code for test case boundary conditions IF (mjg(jj)<25 .or. mjg(jj)>975 .or. mig(ji)<25 .or. mig(ji)>975) THEN u_ice(ji,jj) = zinvw*(ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)) END IF - END_2D - CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1.0_wp ) - ! -#if defined key_agrif -!! CALL agrif_interp_ice( 'U', jter, nn_nevp ) - CALL agrif_interp_ice( 'U' ) -#endif - IF( ln_bdy ) CALL bdy_ice_dyn( 'U' ) + END_2D + IF( nn_hls == 1 ) THEN ; CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1.0_wp ) + ELSE ; CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1.0_wp, v_ice, 'V', -1.0_wp ) + ENDIF ! ELSE ! odd iterations ! - DO_2D( 0, 0, 0, 0 ) + DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! !--- tau_io/(u_oce - u_ice) zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) ) & & + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) ) @@ -639,7 +628,7 @@ SUBROUTINE ice_dyn_rhg_evp( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear & + ( 1._wp - rswitch ) * ( u_ice_b(ji,jj) & & + u_ice (ji,jj) * MAX( 0._wp, zbetau - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 & ) / ( zbetau + 1._wp ) & - & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin + & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin & ) * zmsk00x(ji,jj) ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) & ! previous velocity @@ -649,18 +638,12 @@ SUBROUTINE ice_dyn_rhg_evp( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear & ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin & ) * zmsk00x(ji,jj) ENDIF -!extra code for test case boundary conditions + !extra code for test case boundary conditions IF (mjg(jj)<25 .or. mjg(jj)>975 .or. mig(ji)<25 .or. mig(ji)>975) THEN u_ice(ji,jj) = zinvw*(ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)) END IF - END_2D - CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1.0_wp ) - ! -#if defined key_agrif -!! CALL agrif_interp_ice( 'U', jter, nn_nevp ) - CALL agrif_interp_ice( 'U' ) -#endif - IF( ln_bdy ) CALL bdy_ice_dyn( 'U' ) + END_2D + IF( nn_hls == 1 ) CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1.0_wp ) ! DO_2D( 0, 0, 0, 0 ) ! !--- tau_io/(v_oce - v_ice) @@ -694,7 +677,7 @@ SUBROUTINE ice_dyn_rhg_evp( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear & ) / MAX( zepsi, zmV_t(ji,jj) * ( zbetav + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast & + ( 1._wp - rswitch ) * ( v_ice_b(ji,jj) & & + v_ice (ji,jj) * MAX( 0._wp, zbetav - zdtevp * rn_lf_relax ) & ! static friction => slow decrease to v=0 - & ) / ( zbetav + 1._wp ) & + & ) / ( zbetav + 1._wp ) & & ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin & ) * zmsk00y(ji,jj) ELSE !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) @@ -705,43 +688,74 @@ SUBROUTINE ice_dyn_rhg_evp( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear & ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) ) & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin & ) * zmsk00y(ji,jj) ENDIF -!extra code for test case boundary conditions + !extra code for test case boundary conditions IF (mjg(jj)<25 .or. mjg(jj)>975 .or. mig(ji)<25 .or. mig(ji)>975) THEN v_ice(ji,jj) = zinvw*(ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)) END IF - END_2D - CALL lbc_lnk( 'icedyn_rhg_evp', v_ice, 'V', -1.0_wp ) + END_2D + IF( nn_hls == 1 ) THEN ; CALL lbc_lnk( 'icedyn_rhg_evp', v_ice, 'V', -1.0_wp ) + ELSE ; CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1.0_wp, v_ice, 'V', -1.0_wp ) + ENDIF ! + ENDIF + ! #if defined key_agrif -!! CALL agrif_interp_ice( 'V', jter, nn_nevp ) - CALL agrif_interp_ice( 'V' ) +!! CALL agrif_interp_ice( 'U', jter, nn_nevp ) +!! CALL agrif_interp_ice( 'V', jter, nn_nevp ) + CALL agrif_interp_ice( 'U' ) + CALL agrif_interp_ice( 'V' ) #endif - IF( ln_bdy ) CALL bdy_ice_dyn( 'V' ) - ! - ENDIF - + IF( ln_bdy ) CALL bdy_ice_dyn( 'U' ) + IF( ln_bdy ) CALL bdy_ice_dyn( 'V' ) + ! ! convergence test - IF( nn_rhg_chkcvg == 2 ) CALL rhg_cvg( kt, jter, nn_nevp, u_ice, v_ice, zu_ice, zv_ice ) + IF( nn_rhg_chkcvg == 2 ) CALL rhg_cvg( kt, jter, nn_nevp, u_ice, v_ice, zu_ice, zv_ice, zmsk15 ) + ! ! + ! --- change strength according to advected a_i and v_i (upstream for now) --- ! + IF( ll_advups .AND. ln_str_H79 ) THEN + ! + IF( jter == 1 ) THEN ! init + ALLOCATE( za_i_ups(jpi,jpj,jpl), zv_i_ups(jpi,jpj,jpl) ) + zdt_ups = rDt_ice / REAL( nn_nevp ) + za_i_ups(:,:,:) = a_i(:,:,:) + zv_i_ups(:,:,:) = v_i(:,:,:) + ELSE + CALL lbc_lnk( 'icedyn_rhg_evp', za_i_ups, 'T', 1.0_wp, zv_i_ups, 'T', 1.0_wp ) + ENDIF + ! + CALL rhg_upstream( jter, zdt_ups, u_ice, v_ice, za_i_ups ) ! upstream advection: a_i + CALL rhg_upstream( jter, zdt_ups, u_ice, v_ice, zv_i_ups ) ! upstream advection: v_i + ! + DO_2D( 0, 0, 0, 0 ) ! strength + strength(ji,jj) = rn_pstar * SUM( zv_i_ups(ji,jj,:) ) * EXP( -rn_crhg * ( 1._wp - SUM( za_i_ups(ji,jj,:) ) ) ) + END_2D + ! + IF( jter == nn_nevp ) THEN + DEALLOCATE( za_i_ups, zv_i_ups ) + ENDIF + ENDIF ! ! ==================== ! END DO ! end loop over jter ! ! ! ==================== ! IF( ln_aEVP ) CALL iom_put( 'beta_evp' , zbeta ) ! + IF( ll_advups .AND. ln_str_H79 ) CALL lbc_lnk( 'icedyn_rhg_evp', strength, 'T', 1.0_wp ) + ! !------------------------------------------------------------------------------! - ! 4) Recompute delta, shear and div (inputs for mechanical redistribution) + ! 4) Recompute delta, shear and div (inputs for mechanical redistribution) !------------------------------------------------------------------------------! - DO_2D( 1, 0, 1, 0 ) + DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) ! shear at F points zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) & & + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & - & ) * r1_e1e2f(ji,jj) * zfmask(ji,jj) + & ) * r1_e1e2f(ji,jj) * fimask(ji,jj) END_2D - + DO_2D( 0, 0, 0, 0 ) ! no vector loop - + ! tension**2 at T points zdt = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) & & - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) & @@ -749,29 +763,29 @@ SUBROUTINE ice_dyn_rhg_evp( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear zdt2 = zdt * zdt zten_i(ji,jj) = zdt - + ! shear**2 at T points (doc eq. A16) zds2 = ( zds(ji,jj ) * zds(ji,jj ) * e1e2f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e1e2f(ji-1,jj ) & & + 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 ) + pshear_i(ji,jj) = SQRT( zdt2 + 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) & & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) & - & ) * r1_e1e2t(ji,jj) - + & ) * 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 ) ! delta + 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 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 ) - + ! --- Store the stress tensor for the next time step --- ! pstress1_i (:,:) = zs1 (:,:) pstress2_i (:,:) = zs2 (:,:) @@ -785,8 +799,8 @@ SUBROUTINE ice_dyn_rhg_evp( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear IF( iom_use('utau_oi') .OR. iom_use('vtau_oi') .OR. iom_use('utau_ai') .OR. iom_use('vtau_ai') .OR. & & iom_use('utau_bi') .OR. iom_use('vtau_bi') ) THEN ! - CALL lbc_lnk( 'icedyn_rhg_evp', ztaux_oi, 'U', -1.0_wp, ztauy_oi, 'V', -1.0_wp, & - & ztaux_ai, 'U', -1.0_wp, ztauy_ai, 'V', -1.0_wp, & + CALL lbc_lnk( 'icedyn_rhg_evp', ztaux_oi, 'U', -1.0_wp, ztauy_oi, 'V', -1.0_wp, & + & ztaux_ai, 'U', -1.0_wp, ztauy_ai, 'V', -1.0_wp, & & ztaux_bi, 'U', -1.0_wp, ztauy_bi, 'V', -1.0_wp ) ! CALL iom_put( 'utau_oi' , ztaux_oi * zmsk00 ) @@ -796,40 +810,41 @@ SUBROUTINE ice_dyn_rhg_evp( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear CALL iom_put( 'utau_bi' , ztaux_bi * zmsk00 ) CALL iom_put( 'vtau_bi' , ztauy_bi * zmsk00 ) ENDIF - + ! --- divergence, shear and strength --- ! 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 ! --- Stress tensor invariants (SIMIP diags) --- ! IF( iom_use('normstr') .OR. iom_use('sheastr') ) THEN ! ALLOCATE( zsig_I(jpi,jpj) , zsig_II(jpi,jpj) ) - ! - DO_2D( 1, 1, 1, 1 ) - + ! + 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 ) + zfac = strength(ji,jj) / ( pdelta_i(ji,jj) + rn_creepl ) zsig1 = zfac * ( pdivu_i(ji,jj) - pdelta_i(ji,jj) ) zsig2 = zfac * z1_ecc2 * zten_i(ji,jj) zsig12 = zfac * z1_ecc2 * pshear_i(ji,jj) - + ! 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 ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) ) ! 2nd '' '', aka maximum shear stress - - END_2D + 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 + + 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 - + DEALLOCATE ( zsig_I, zsig_II ) - + ENDIF ! --- Normalized stress tensor principal components --- ! @@ -838,33 +853,33 @@ SUBROUTINE ice_dyn_rhg_evp( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear ! Recommendation 2 : need to use deformations at PREVIOUS iterate for viscosities 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( 1, 1, 1, 1 ) - - ! Ice stresses computed with **viscosities** (delta, p/delta) at **previous** iterates + 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 ! following Lemieux & Dupont (2020) zfac = zp_delt(ji,jj) zsig1 = zfac * ( pdivu_i(ji,jj) - ( zdelta(ji,jj) + rn_creepl ) ) zsig2 = zfac * z1_ecc2 * zten_i(ji,jj) zsig12 = zfac * z1_ecc2 * pshear_i(ji,jj) - + ! 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 ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) ) ! 2nd '' '', aka maximum shear stress - + 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 + ! Normalized principal stresses (used to display the ellipse) z1_strength = 1._wp / MAX( 1._wp, strength(ji,jj) ) zsig1_p(ji,jj) = ( zsig_I(ji,jj) + zsig_II(ji,jj) ) * z1_strength zsig2_p(ji,jj) = ( zsig_I(ji,jj) - zsig_II(ji,jj) ) * z1_strength - END_2D + END_2D ! - CALL iom_put( 'sig1_pnorm' , zsig1_p ) - CALL iom_put( 'sig2_pnorm' , zsig2_p ) + IF( iom_use('sig1_pnorm') ) CALL iom_put( 'sig1_pnorm' , zsig1_p(:,:) * zmsk00(:,:) ) + IF( iom_use('sig2_pnorm') ) CALL iom_put( 'sig2_pnorm' , zsig2_p(:,:) * zmsk00(:,:) ) DEALLOCATE( zsig1_p , zsig2_p , zsig_I, zsig_II ) - + ENDIF ! --- SIMIP --- ! @@ -909,7 +924,7 @@ SUBROUTINE ice_dyn_rhg_evp( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear & zdiag_xatrp , 'U', -1.0_wp, zdiag_yatrp , 'V', -1.0_wp ) CALL iom_put( 'xmtrpice' , zdiag_xmtrp_ice ) ! X-component of sea-ice mass transport (kg/s) - CALL iom_put( 'ymtrpice' , zdiag_ymtrp_ice ) ! Y-component of sea-ice mass transport + CALL iom_put( 'ymtrpice' , zdiag_ymtrp_ice ) ! Y-component of sea-ice mass transport CALL iom_put( 'xmtrpsnw' , zdiag_xmtrp_snw ) ! X-component of snow mass transport (kg/s) CALL iom_put( 'ymtrpsnw' , zdiag_ymtrp_snw ) ! Y-component of snow mass transport CALL iom_put( 'xatrp' , zdiag_xatrp ) ! X-component of ice area transport @@ -931,17 +946,15 @@ SUBROUTINE ice_dyn_rhg_evp( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear & ABS( v_ice(:,:) - zv_ice(:,:) ) * vmask(:,:,1) ) * zmsk15(:,:) ) ENDIF ENDIF - ENDIF - ! - DEALLOCATE( zmsk00, zmsk15 ) + ENDIF ! END SUBROUTINE ice_dyn_rhg_evp - SUBROUTINE rhg_cvg( kt, kiter, kitermax, pu, pv, pub, pvb ) + SUBROUTINE rhg_cvg( kt, kiter, kitermax, pu, pv, pub, pvb, pmsk15 ) !!---------------------------------------------------------------------- !! *** ROUTINE rhg_cvg *** - !! + !! !! ** Purpose : check convergence of oce rheology !! !! ** Method : create a file ice_cvg.nc containing the convergence of ice velocity @@ -949,18 +962,22 @@ SUBROUTINE rhg_cvg( kt, kiter, kitermax, pu, pv, pub, pvb ) !! uice_cvg = MAX( u(t+1) - u(t) , v(t+1) - v(t) ) !! This routine is called every sub-iteration, so it is cpu expensive !! - !! ** Note : for the first sub-iteration, uice_cvg is set to 0 (too large otherwise) + !! ** Note : for the first sub-iteration, uice_cvg is set to 0 (too large otherwise) !!---------------------------------------------------------------------- INTEGER , INTENT(in) :: kt, kiter, kitermax ! ocean time-step index REAL(wp), DIMENSION(:,:), INTENT(in) :: pu, pv, pub, pvb ! now and before velocities + REAL(wp), DIMENSION(:,:), INTENT(in) :: pmsk15 !! INTEGER :: it, idtime, istatus INTEGER :: ji, jj ! dummy loop indices - REAL(wp) :: zresm ! local real + REAL(wp) :: zresm ! local real CHARACTER(len=20) :: clname - REAL(wp), DIMENSION(jpi,jpj) :: zres ! check convergence + LOGICAL :: ll_maxcvg + REAL(wp), DIMENSION(jpi,jpj,2) :: zres + REAL(wp), DIMENSION(2) :: ztmp !!---------------------------------------------------------------------- - + ll_maxcvg = .FALSE. + ! ! create file IF( kt == nit000 .AND. kiter == 1 ) THEN ! @@ -975,47 +992,51 @@ SUBROUTINE rhg_cvg( kt, kiter, kitermax, pu, pv, pub, pvb ) IF( .NOT. Agrif_Root() ) clname = TRIM(Agrif_CFixed())//"_"//TRIM(clname) istatus = NF90_CREATE( TRIM(clname), NF90_CLOBBER, ncvgid ) istatus = NF90_DEF_DIM( ncvgid, 'time' , NF90_UNLIMITED, idtime ) - istatus = NF90_DEF_VAR( ncvgid, 'uice_cvg', NF90_DOUBLE , (/ idtime /), nvarid ) + istatus = NF90_DEF_VAR( ncvgid, 'uice_cvg', NF90_DOUBLE , (/ idtime /), nvarid ) istatus = NF90_ENDDEF(ncvgid) ENDIF ! ENDIF ! time - it = ( kt - 1 ) * kitermax + kiter - + it = ( kt - nit000 ) * kitermax + kiter + ! convergence IF( kiter == 1 ) THEN ! remove the first iteration for calculations of convergence (always very large) zresm = 0._wp ELSE - DO_2D( 1, 1, 1, 1 ) - zres(ji,jj) = MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), & - & ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * zmsk15(ji,jj) - END_2D - - ! cut of the boundary of the box (forced velocities) - IF (mjg(jj)<=30 .or. mjg(jj)>970 .or. mig(ji)<=30 .or. mig(ji)>970) THEN - zres(ji,jj) = 0._wp - END IF - - zresm = MAXVAL( zres ) - CALL mpp_max( 'icedyn_rhg_evp', zresm ) ! max over the global domain + zresm = 0._wp + IF( ll_maxcvg ) THEN ! error max over the domain + DO_2D( 0, 0, 0, 0 ) + zresm = MAX( zresm, MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), & + & ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * pmsk15(ji,jj) ) + END_2D + CALL mpp_max( 'icedyn_rhg_evp', zresm ) + ELSE ! error averaged over the domain + DO_2D( 0, 0, 0, 0 ) + zres(ji,jj,1) = MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), & + & ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * pmsk15(ji,jj) + zres(ji,jj,2) = pmsk15(ji,jj) + END_2D + ztmp(:) = glob_sum_vec( 'icedyn_rhg_evp', zres ) + IF( ztmp(2) /= 0._wp ) zresm = ztmp(1) / ztmp(2) + ENDIF ENDIF IF( lwm ) THEN ! write variables istatus = NF90_PUT_VAR( ncvgid, nvarid, (/zresm/), (/it/), (/1/) ) ! close file - IF( kt == nitend - nn_fsbc + 1 ) istatus = NF90_CLOSE(ncvgid) + IF( kt == nitend - nn_fsbc + 1 .AND. kiter == kitermax ) istatus = NF90_CLOSE(ncvgid) ENDIF - + END SUBROUTINE rhg_cvg SUBROUTINE rhg_evp_rst( cdrw, kt ) !!--------------------------------------------------------------------- !! *** ROUTINE rhg_evp_rst *** - !! + !! !! ** Purpose : Read or write RHG file in restart file !! !! ** Method : use of IOM library @@ -1059,15 +1080,86 @@ SUBROUTINE rhg_evp_rst( cdrw, kt ) IF(lwp) WRITE(numout,*) '---- rhg-rst ----' iter = kt + nn_fsbc - 1 ! ice restarts are written at kt == nitrst - nn_fsbc + 1 ! - CALL iom_rstput( iter, nitrst, numriw, 'stress1_i' , stress1_i ) - CALL iom_rstput( iter, nitrst, numriw, 'stress2_i' , stress2_i ) + CALL iom_rstput( iter, nitrst, numriw, 'stress1_i' , stress1_i ) + CALL iom_rstput( iter, nitrst, numriw, 'stress2_i' , stress2_i ) CALL iom_rstput( iter, nitrst, numriw, 'stress12_i', stress12_i ) ! ENDIF ! END SUBROUTINE rhg_evp_rst - + SUBROUTINE rhg_upstream( jter, pdt, pu, pv, pt ) + !!--------------------------------------------------------------------- + !! *** ROUTINE rhg_upstream *** + !! + !! ** Purpose : compute the upstream fluxes and upstream guess of tracer + !!---------------------------------------------------------------------- + INTEGER , INTENT(in ) :: jter + REAL(wp) , INTENT(in ) :: pdt ! tracer time-step + REAL(wp), DIMENSION(:,: ) , INTENT(in ) :: pu, pv ! 2 ice velocity components + REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pt ! tracer fields + ! + INTEGER :: ji, jj, jl ! dummy loop indices + REAL(wp) :: ztra ! local scalar + LOGICAL :: ll_upsxy = .TRUE. + REAL(wp), DIMENSION(jpi,jpj) :: zfu_ups, zfv_ups, zpt ! upstream fluxes and tracer guess + !!---------------------------------------------------------------------- + DO jl = 1, jpl + IF( .NOT. ll_upsxy ) THEN !** no alternate directions **! + ! + DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) + zfu_ups(ji,jj) = MAX(pu(ji,jj)*e2u(ji,jj), 0._wp) * pt(ji,jj,jl) + MIN(pu(ji,jj)*e2u(ji,jj), 0._wp) * pt(ji+1,jj,jl) + zfv_ups(ji,jj) = MAX(pv(ji,jj)*e1v(ji,jj), 0._wp) * pt(ji,jj,jl) + MIN(pv(ji,jj)*e1v(ji,jj), 0._wp) * pt(ji,jj+1,jl) + END_2D + ! + ELSE !** alternate directions **! + ! + IF( MOD(jter,2) == 1 ) THEN !== odd ice time step: adv_x then adv_y ==! + ! + DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls ) !-- flux in x-direction + zfu_ups(ji,jj) = MAX( pu(ji,jj)*e2u(ji,jj), 0._wp ) * pt(ji ,jj,jl) + & + & MIN( pu(ji,jj)*e2u(ji,jj), 0._wp ) * pt(ji+1,jj,jl) + END_2D + ! + DO_2D( nn_hls-1, nn_hls-1, nn_hls, nn_hls ) !-- first guess of tracer from u-flux + ztra = - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) ) + zpt(ji,jj) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) + END_2D + ! + DO_2D( nn_hls-1, nn_hls-1, nn_hls, nn_hls-1 ) !-- flux in y-direction + zfv_ups(ji,jj) = MAX( pv(ji,jj)*e1v(ji,jj), 0._wp ) * zpt(ji,jj ) + & + & MIN( pv(ji,jj)*e1v(ji,jj), 0._wp ) * zpt(ji,jj+1) + END_2D + ! + ELSE !== even ice time step: adv_y then adv_x ==! + ! + DO_2D( nn_hls, nn_hls, nn_hls, nn_hls-1 ) !-- flux in y-direction + zfv_ups(ji,jj) = MAX( pv(ji,jj)*e1v(ji,jj), 0._wp ) * pt(ji,jj ,jl) + & + & MIN( pv(ji,jj)*e1v(ji,jj), 0._wp ) * pt(ji,jj+1,jl) + END_2D + ! + DO_2D( nn_hls, nn_hls, nn_hls-1, nn_hls-1 ) !-- first guess of tracer from v-flux + ztra = - ( zfv_ups(ji,jj) - zfv_ups(ji,jj-1) ) + zpt(ji,jj) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) + END_2D + ! + DO_2D( nn_hls, nn_hls-1, nn_hls-1, nn_hls-1 ) !-- flux in x-direction + zfu_ups(ji,jj) = MAX( pu(ji,jj)*e2u(ji,jj), 0._wp ) * zpt(ji ,jj) + & + & MIN( pu(ji,jj)*e2u(ji,jj), 0._wp ) * zpt(ji+1,jj) + END_2D + ! + ENDIF + ! + ENDIF + ! + DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) + ztra = - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) + zfv_ups(ji,jj) - zfv_ups(ji,jj-1) ) + pt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) + END_2D + END DO + ! + END SUBROUTINE rhg_upstream + #else !!---------------------------------------------------------------------- !! Default option Empty module NO SI3 sea-ice model diff --git a/ICE_RHEO/MY_SRC/usrdef_zgr.F90 b/ICE_RHEO/MY_SRC/usrdef_zgr.F90 index 096e524..7c295aa 100644 --- a/ICE_RHEO/MY_SRC/usrdef_zgr.F90 +++ b/ICE_RHEO/MY_SRC/usrdef_zgr.F90 @@ -14,6 +14,7 @@ MODULE usrdef_zgr !!--------------------------------------------------------------------- USE oce ! ocean variables USE usrdef_nam ! User defined : namelist variables + USE dom_oce , ONLY : mig, mjg ! USE in_out_manager ! I/O manager USE lbclnk ! ocean lateral boundary conditions (or mpp link) @@ -24,6 +25,8 @@ MODULE usrdef_zgr PUBLIC usr_def_zgr ! called by domzgr.F90 + !! * Substitutions +# include "do_loop_substitute.h90" !!---------------------------------------------------------------------- !! NEMO/OCE 4.0 , NEMO Consortium (2018) !! $Id: usrdef_zgr.F90 10074 2018-08-28 16:15:49Z nicolasmartin $ @@ -53,6 +56,7 @@ SUBROUTINE usr_def_zgr( ld_zco , ld_zps , ld_sco , ld_isfcav, & ! type o INTEGER , DIMENSION(:,:) , INTENT(out) :: k_top, k_bot ! first & last ocean level ! INTEGER :: jk, k_dz ! dummy indices + INTEGER :: ji, jj ! dummy loop indices !!---------------------------------------------------------------------- ! IF(lwp) WRITE(numout,*) @@ -85,6 +89,14 @@ SUBROUTINE usr_def_zgr( ld_zco , ld_zps , ld_sco , ld_isfcav, & ! type o ! ! no ocean cavities : top ocean level is ONE, except over land k_top(:,:) = 1 + + ! add island + DO_2D( 1, 1, 1, 1 ) + IF (mjg(jj)>=900000/rn_dy .and. mjg(jj)<1100000/rn_dy .and. mig(ji)>=900000/rn_dx .and. mig(ji)<1100000/rn_dx) THEN + k_top(ji,jj) = 0 + END IF + END_2D + ! ! !== z-coordinate ==! (step-like topography) ! !* bottom ocean compute from the depth of grid-points