Changeset 15507 for NEMO/branches/UKMO/NEMO_4.0.4_ice_strength/src
- Timestamp:
- 2021-11-15T15:22:24+01:00 (3 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/UKMO/NEMO_4.0.4_ice_strength/src/ICE/icedyn_rdgrft.F90
r15346 r15507 43 43 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: opning ! rate of opening due to divergence/shear 44 44 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: closing_gross ! rate at which area removed, not counting area of new ridges 45 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: zaksum ! normalisation factor 45 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: zaksum ! normalisation factor 46 46 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: zasum ! sum of ice area plus open water 47 47 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: apartf ! participation function; fraction of ridging/closing associated w/ category n … … 49 49 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: hrmin ! minimum ridge thickness 50 50 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: hrmax ! maximum ridge thickness 51 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: hrexp ! e-folding ridge thickness52 51 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: hraft ! thickness of rafted ice 53 52 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: hi_hrdg ! thickness of ridging ice / mean ridge thickness … … 95 94 ALLOCATE( closing_net(jpij) , opning(jpij) , closing_gross(jpij) , & 96 95 & apartf(jpij,0:jpl) , hrmin(jpij,jpl) , hrmax(jpij,jpl) , zhi(jpij,jpl), & 97 & zaksum(jpij) , hraft(jpij,jpl) , hrexp(jpij,jpl) ,&98 & hi_hrdg(jpij,jpl) , araft(jpij,jpl) , aridge(jpij,jpl) , 96 & zaksum(jpij) , hraft(jpij,jpl) , & 97 & hi_hrdg(jpij,jpl) , araft(jpij,jpl) , aridge(jpij,jpl) , zasum(jpij), & 99 98 & ze_i_2d(jpij,nlay_i,jpl), ze_s_2d(jpij,nlay_s,jpl), STAT=ice_dyn_rdgrft_alloc ) 100 99 … … 146 145 INTEGER , DIMENSION(jpij) :: iptidx ! compute ridge/raft or not 147 146 REAL(wp), DIMENSION(jpij) :: zdivu, zdelt ! 1D divu_i & delta_i 148 REAL(wp), DIMENSION(jpij) :: zdivu_adv ! 1D divu_i as implied by transport scheme149 147 ! 150 148 INTEGER, PARAMETER :: jp_itermax = 20 … … 194 192 ! - ice area added in new ridges 195 193 closing_net(ji) = rn_csrdg * 0.5_wp * ( zdelt(ji) - ABS( zdivu(ji) ) ) - MIN( zdivu(ji), 0._wp ) 196 ! 197 ! zdivu_adv(ji) = (1._wp - zasum(ji)) / rdt_ice ! Need to do more testing/checking on this 198 199 ! WRITE(numout,*)'velocity div ',zdivu(ji), 'transport div ',zdivu_adv(ji) 200 194 ! 201 195 IF( zdivu(ji) < 0._wp ) closing_net(ji) = MAX( closing_net(ji), -zdivu(ji) ) ! make sure the closing rate is large enough 202 ! IF( zdivu_adv(ji) < 0._wp ) closing_net(ji) = MAX( closing_net(ji), -zdivu_adv(ji) )203 204 196 ! ! to give asum = 1.0 after ridging 205 197 ! Opening rate (non-negative) that will give asum = 1.0 after ridging. 206 198 opning(ji) = closing_net(ji) + zdivu(ji) 207 ! opning(ji) = closing_net(ji) + zdivu_adv(ji) 208 199 ! 209 200 END DO 210 201 ! … … 310 301 INTEGER :: ji, jl ! dummy loop indices 311 302 REAL(wp) :: z1_gstar, z1_astar, zhmean, zfac ! local scalar 312 ! REAL(wp), DIMENSION(jpij) :: zasum, z1_asum ! sum of a_i+ato_i and inverse313 303 REAL(wp), DIMENSION(jpij) :: z1_asum ! inverse of sum of a_i+ato_i 314 304 REAL(wp), DIMENSION(jpij,-1:jpl) :: zGsum ! zGsum(n) = sum of areas in categories 0 to n … … 324 314 ELSEWHERE ; zhi(1:npti,:) = 0._wp 325 315 END WHERE 326 ! Make sure ice thickness is not below minimum in any category (from iceitd.F90; no minimum was causing issues for R75 strength) 316 ! Make sure ice thickness is not below the minimum (from iceitd.F90; no minimum was causing issues for R75 317 ! strength) 327 318 WHERE( pa_i(1:npti,:) > epsi10 .AND. zhi(1:npti,:) < rn_himin ) 328 319 pa_i(1:npti,:) = pa_i(1:npti,:) * zhi(1:npti,:) / rn_himin 329 320 zhi(1:npti,:) = rn_himin 321 END WHERE 322 323 ! Set tiny values of pa_i to zero 324 ! If this isn't done, can end up with zhi=0 but apartf>0 325 WHERE( pa_i(1:npti,:) <= epsi10 ) ; pa_i(1:npti,:) = 0._wp 330 326 END WHERE 331 327 … … 344 340 ! Compute total area of ice plus open water. 345 341 ! This is in general not equal to one because of divergence during transport 346 ! zasum(1:npti) = pato_i(1:npti) + SUM( pa_i(1:npti,:), dim=2 ) 347 ! EF: to match CICE, this should be the following. What is pa_i(:,0)?? Very small. 348 zasum(1:npti) = pato_i(1:npti) + SUM( pa_i(1:npti,1:jpl), dim=2 ) 342 zasum(1:npti) = pato_i(1:npti) + SUM( pa_i(1:npti,1:jpl), dim=2 ) 349 343 350 344 WHERE( zasum(1:npti) > epsi10 ) ; z1_asum(1:npti) = 1._wp / zasum(1:npti) … … 359 353 zGsum(1:npti,0 ) = pato_i(1:npti) * z1_asum(1:npti) 360 354 DO jl = 1, jpl 361 zGsum(1:npti,jl) = ( pato_i(1:npti) + SUM( pa_i(1:npti,1:jl), dim=2 ) ) * z1_asum(1:npti) ! sum(1:jl) is ok(and not jpl)355 zGsum(1:npti,jl) = ( pato_i(1:npti) + SUM( pa_i(1:npti,1:jl), dim=2 ) ) * z1_asum(1:npti) ! sum(1:jl) is correct (and not jpl) 362 356 END DO 363 357 ! … … 457 451 hrmin (ji,jl) = MIN( 2._wp * zhi(ji,jl), 0.5_wp * ( zhmean + zhi(ji,jl) ) ) 458 452 hrmax (ji,jl) = 2._wp * zhmean - hrmin(ji,jl) 459 hrexp (ji,jl) = rn_murdg * SQRT( zhi(ji,jl) )460 453 hraft (ji,jl) = zhi(ji,jl) * zfac 461 hi_hrdg(ji,jl) = zhi(ji,jl) / MAX( zhmean, epsi20 )454 hi_hrdg(ji,jl) = zhi(ji,jl) / zhmean 462 455 ! 463 456 ! Normalization factor : zaksum, ensures mass conservation … … 467 460 hrmin (ji,jl) = 0._wp 468 461 hrmax (ji,jl) = 0._wp 469 hrexp (ji,jl) = 0._wp470 462 hraft (ji,jl) = 0._wp 471 463 hi_hrdg(ji,jl) = 1._wp … … 788 780 !! dissipated per unit area removed from the ice pack under compression, 789 781 !! and assumed proportional to the change in potential energy caused 790 !! by ridging. [**UPDATE THIS ---> Note that only Hibler's formulation is stable and that791 !! ice strength has to be smoothed**]782 !! by ridging. Note that ice strength using Hibler's formulation must be 783 !! smoothed. 792 784 !!---------------------------------------------------------------------- 793 785 INTEGER :: ji, jj, jl ! dummy loop indices … … 796 788 REAL(wp) :: zp, z1_3 ! local scalars 797 789 REAL(wp) :: Cp ! proportional constant for PE 798 !REAL(wp) :: hi ! ice thickness of category (m)799 790 REAL(wp) :: h2rdg ! mean value of h^2 for new ridge 800 791 REAL(wp) :: dh2rdg ! change in mean value of h^2 per unit area consumed by ridging 801 792 REAL(wp), DIMENSION(jpi,jpj) :: zworka ! temporary array used here 802 793 REAL(wp), DIMENSION(jpi,jpj) :: zstrp1, zstrp2 ! strength at previous time steps 803 REAL(wp), DIMENSION(jpij,jpl) :: krdg ! mean ridge thickness / thickness of ridging ice804 794 REAL(wp), DIMENSION(jpij) :: strength_1d ! 1-dimensional strength 805 795 REAL(wp), DIMENSION(jpij,jpl) :: strength_pe ! PE component of R75 strength … … 809 799 IF( ln_str_H79 ) THEN ! Ice strength => Hibler (1979) method ! 810 800 ! !--------------------------------------------------! 811 812 801 strength(:,:) = rn_pstar * SUM( v_i(:,:,:), dim=3 ) * EXP( -rn_crhg * ( 1._wp - SUM( a_i(:,:,:), dim=3 ) ) ) 802 ismooth = 1 813 803 ! !--------------------------------------------------! 814 804 ELSE IF( ln_str_R75 ) THEN ! Ice strength => Rothrock (1975) method ! 815 ! !--------------------------------------------------! 816 IF ( kt_ice == nit000 ) THEN 817 818 ! Replace this with input from start file 819 strength(:,:) = rn_pstar * SUM( v_i(:,:,:), dim=3 ) * EXP( -rn_crhg * ( 1._wp - SUM( a_i(:,:,:), dim=3 ) ) ) 820 ismooth = 1 821 822 ELSE 823 824 !-------------------------------- 825 ! 0) Identify grid cells with ice 826 !-------------------------------- 827 at_i(:,:) = SUM( a_i(:,:,:), dim=3 ) 828 vt_i(:,:) = SUM( v_i(:,:,:), dim=3 ) 829 ! 830 npti = 0 ; nptidx(:) = 0 831 DO jj = 1, jpj 832 DO ji = 1, jpi 833 IF ( at_i(ji,jj) > epsi10 ) THEN 834 npti = npti + 1 835 nptidx( npti ) = (jj - 1) * jpi + ji 836 ENDIF 837 END DO 838 END DO 839 840 CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i ) 841 CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i ) 842 CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti) , ato_i ) 843 844 845 !----------------------------------- 846 ! Get thickness distribution of ice 847 !----------------------------------- 848 849 CALL rdgrft_prep( a_i_2d, v_i_2d, ato_i_1d, closing_net ) 850 851 ! Don't need closing rate here (unless only calculating strength for ridging ice), could split this function up. 852 853 Cp = 0.5_wp * grav * (rhow-rhoi) * rhoi / rhow ! put this elsewhere, perhaps ice_prep? 854 855 856 IF ( npti > 0 ) THEN 857 858 DO jl = 1, jpl ! categories 859 DO ji = 1, npti ! grid cells 860 861 IF ( a_i_2d(ji,jl) > epsi10 .and. apartf(ji,jl) > 0._wp ) THEN 862 863 ! Move this to rdgrft_prep 864 krdg(ji,jl) = (hrmin(ji,jl) + hrmax(ji,jl))/2._wp * 1._wp/zhi(ji,jl) 865 866 867 ! Add a logical for distribution 868 ! Uniform distribution 869 h2rdg = (1._wp/3._wp) * (hrmax(ji,jl)**3._wp - hrmin(ji,jl)**3._wp) & 870 / (hrmax(ji,jl) - hrmin(ji,jl)) 871 872 873 ! Exponential distribution 874 ! h2rdg = hrmin(ji,jl) * hrmin(ji,jl) & 875 ! + 2._wp * hrmin(ji,jl) * hrexp(ji,jl) & 876 ! + 2._wp * hrexp(ji,jl) * hrexp(ji,jl) 877 878 dh2rdg = -zhi(ji,jl)*zhi(ji,jl) + h2rdg / krdg(ji,jl) 879 880 strength_pe(ji,jl) = apartf(ji,jl) * dh2rdg 881 882 ELSE 883 884 strength_pe(ji,jl) = 0._wp 885 886 END IF 887 888 END DO ! grid cells 889 END DO ! categories 890 891 892 strength_pe_sum(:) = SUM(strength_pe(:,:), dim=2) ! Sum over all categories 893 894 895 DO ji = 1, npti 896 IF ( zaksum(ji) > epsi10 ) THEN 897 strength_1d(ji) = rn_Cf * Cp * strength_pe_sum(ji) / zaksum(ji) 898 ELSE 899 strength_1d(ji) = 0._wp 900 END IF 901 END DO 902 903 ! Convert strength_1d to 2d 904 905 CALL tab_1d_2d( npti, nptidx(1:npti), strength_1d(1:npti), strength(:,:) ) 906 907 ! Enforce a maximum for strength (- Look into why numbers can be so large) 908 WHERE( strength(:,:) > 1.5e5) ; strength(:,:) = 1.5e5 909 END WHERE 910 911 ismooth = 0 ! Buffer size error with lbc_lnk... 912 913 END IF ! number of ice points 914 END IF ! timestep 805 ! !--------------------------------------------------! 806 ! 807 ! Initialise the strength field at each timestep 808 strength(:,:) = 0._wp 809 ismooth = 0 810 811 !-------------------------------- 812 ! 0) Identify grid cells with ice 813 !-------------------------------- 814 at_i(:,:) = SUM( a_i(:,:,:), dim=3 ) 815 vt_i(:,:) = SUM( v_i(:,:,:), dim=3 ) 816 ato_i(:,:) = 1 - at_i(:,:) 817 ! 818 npti = 0 ; nptidx(:) = 0 819 DO jj = 1, jpj 820 DO ji = 1, jpi 821 IF ( at_i(ji,jj) > epsi10 ) THEN 822 npti = npti + 1 823 nptidx( npti ) = (jj - 1) * jpi + ji 824 ENDIF 825 END DO 826 END DO 827 828 CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i ) 829 CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i ) 830 CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti) , ato_i ) 831 832 !----------------------------------- 833 ! Get thickness distribution of ice 834 !----------------------------------- 835 836 Cp = 0.5_wp * grav * (rhow-rhoi) * rhoi / rhow ! EF: This should go in an initialisation routine but unclear where, won't let me set as a parameter above 837 838 CALL rdgrft_prep( a_i_2d, v_i_2d, ato_i_1d, closing_net ) 839 ! EF: Note that closing rate is not needed here, could split this function up, and refactor as initial rdgrft_prep is 840 ! called here, before going into ridging iterations 841 842 ! Initialise 843 strength_pe_sum(:) = 0._wp 844 845 IF ( npti > 0 ) THEN 846 847 DO jl = 1, jpl ! categories 848 DO ji = 1, npti ! grid cells 849 IF ( apartf(ji,jl) > 0._wp ) THEN 850 851 ! Uniform distribution of ridged ice 852 h2rdg = (1._wp/3._wp) * (hrmax(ji,jl)**3._wp - hrmin(ji,jl)**3._wp) & 853 / (hrmax(ji,jl) - hrmin(ji,jl)) 854 855 dh2rdg = -zhi(ji,jl)*zhi(ji,jl) + h2rdg * hi_hrdg(ji,jl) 856 857 strength_pe_sum(ji) = strength_pe_sum(ji) + apartf(ji,jl) * dh2rdg 858 859 ELSE 860 861 strength_pe_sum(ji) = strength_pe_sum(ji) + 0._wp 862 863 END IF 864 865 END DO ! grid cells 866 END DO ! categories 867 868 869 DO ji = 1, npti 870 IF ( zaksum(ji) > epsi10 ) THEN ! Unclear why is this check should be needed? Occasional x10E-24... 871 strength_1d(ji) = rn_Cf * Cp * strength_pe_sum(ji) / zaksum(ji) 872 ELSE 873 strength_1d(ji) = 0._wp 874 END IF 875 END DO 876 877 878 ! Enforce a maximum for strength 879 WHERE( strength_1d(:) > 2.0e5) ; strength_1d(:) = 2.0e5 880 END WHERE 881 882 883 ! Convert strength_1d to 2d 884 CALL tab_1d_2d( npti, nptidx(1:npti), strength_1d(1:npti), strength(:,:) ) 885 886 END IF ! number of ice points 915 887 916 888 … … 920 892 strength(:,:) = 0._wp 921 893 ismooth = 0 922 ENDIF 894 ENDIF ! Ice strength formulation 923 895 ! 924 925 896 !--------------------------------------------------! 926 897 SELECT CASE( ismooth ) ! Smoothing ice strength ! … … 928 899 CASE( 0 ) !--- No smoothing 929 900 930 931 !CALL lbc_lnk( 'icedyn_rdgrft', strength, 'T', 1. ) 901 CALL lbc_lnk( 'icedyn_rdgrft', strength, 'T', 1. ) 932 902 933 903 CASE( 1 ) !--- Spatial smoothing … … 1113 1083 ENDIF 1114 1084 ! 1085 IF ( ( ln_str_H79 .AND. ln_str_R75 ) .OR. ( .NOT.ln_str_H79 .AND. .NOT.ln_str_R75 ) ) THEN 1086 CALL ctl_stop( 'ice_dyn_rdgrft_init: choose one and only one ice strength formulation (ln_str_H79 or ln_str_R75)' ) 1087 ENDIF 1088 ! 1115 1089 IF ( ( ln_partf_lin .AND. ln_partf_exp ) .OR. ( .NOT.ln_partf_lin .AND. .NOT.ln_partf_exp ) ) THEN 1116 1090 CALL ctl_stop( 'ice_dyn_rdgrft_init: choose one and only one participation function (ln_partf_lin or ln_partf_exp)' )
Note: See TracChangeset
for help on using the changeset viewer.
