Changeset 15344
- Timestamp:
- 2021-10-07T17:06:57+02:00 (20 months ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/UKMO/NEMO_4.0.4_ice_strength/src/ICE/icedyn_rdgrft.F90
r14075 r15344 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 46 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: zasum ! sum of ice area plus open water 45 47 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: apartf ! participation function; fraction of ridging/closing associated w/ category n 48 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zhi ! ice thickness in category n 46 49 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: hrmin ! minimum ridge thickness 47 50 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: hrmax ! maximum ridge thickness 51 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: hrexp ! e-folding ridge thickness 48 52 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: hraft ! thickness of rafted ice 49 53 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: hi_hrdg ! thickness of ridging ice / mean ridge thickness … … 54 58 ! 55 59 REAL(wp), PARAMETER :: hrdg_hi_min = 1.1_wp ! min ridge thickness multiplier: min(hrdg/hi) 56 REAL(wp), PARAMETER :: hi_hrft = 0.5_wp ! rafting multipli yer: (hi/hraft)60 REAL(wp), PARAMETER :: hi_hrft = 0.5_wp ! rafting multiplier: (hi/hraft) 57 61 ! 58 62 ! ** namelist (namdyn_rdgrft) ** 59 63 LOGICAL :: ln_str_H79 ! ice strength parameterization (Hibler79) 60 64 REAL(wp) :: rn_pstar ! determines ice strength, Hibler JPO79 65 LOGICAL :: ln_str_R75 ! ice strength parameterization (Rothrock75) 66 REAL(wp) :: rn_Cf ! ratio of ridging work to PE change in ridging, R75 67 REAL(wp) :: rn_murdg ! gives e-folding scale of ridged ice (m^.5), R75 61 68 REAL(wp) :: rn_csrdg ! fraction of shearing energy contributing to ridging 62 69 LOGICAL :: ln_partf_lin ! participation function linear (Thorndike et al. (1975)) … … 86 93 !! *** ROUTINE ice_dyn_rdgrft_alloc *** 87 94 !!------------------------------------------------------------------- 88 ALLOCATE( closing_net(jpij) , opning(jpij) , closing_gross(jpij) , & 89 & apartf(jpij,0:jpl) , hrmin (jpij,jpl) , hraft(jpij,jpl) , aridge(jpij,jpl), & 90 & hrmax (jpij,jpl) , hi_hrdg(jpij,jpl) , araft(jpij,jpl) , & 95 ALLOCATE( closing_net(jpij) , opning(jpij) , closing_gross(jpij) , & 96 & 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) , zasum(jpij), & 91 99 & ze_i_2d(jpij,nlay_i,jpl), ze_s_2d(jpij,nlay_s,jpl), STAT=ice_dyn_rdgrft_alloc ) 92 100 … … 138 146 INTEGER , DIMENSION(jpij) :: iptidx ! compute ridge/raft or not 139 147 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 scheme 140 149 ! 141 150 INTEGER, PARAMETER :: jp_itermax = 20 … … 186 195 closing_net(ji) = rn_csrdg * 0.5_wp * ( zdelt(ji) - ABS( zdivu(ji) ) ) - MIN( zdivu(ji), 0._wp ) 187 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 188 201 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 189 204 ! ! to give asum = 1.0 after ridging 190 205 ! Opening rate (non-negative) that will give asum = 1.0 after ridging. 191 206 opning(ji) = closing_net(ji) + zdivu(ji) 207 ! opning(ji) = closing_net(ji) + zdivu_adv(ji) 208 192 209 END DO 193 210 ! … … 277 294 278 295 279 SUBROUTINE rdgrft_prep( pa_i, pv_i, pato_i, pclosing_net ) 296 297 SUBROUTINE rdgrft_prep( pa_i_in, pv_i, pato_i, pclosing_net ) 280 298 !!------------------------------------------------------------------- 281 299 !! *** ROUTINE rdgrft_prep *** … … 287 305 !!------------------------------------------------------------------- 288 306 REAL(wp), DIMENSION(:) , INTENT(in) :: pato_i, pclosing_net 289 REAL(wp), DIMENSION(:,:), INTENT(in) :: pa_i, pv_i 307 REAL(wp), DIMENSION(:,:), INTENT(in) :: pa_i_in, pv_i 308 REAL(wp), DIMENSION(jpij,jpl) :: pa_i ! area adjusted for minimum zhi 290 309 !! 291 310 INTEGER :: ji, jl ! dummy loop indices 292 311 REAL(wp) :: z1_gstar, z1_astar, zhmean, zfac ! local scalar 293 REAL(wp), DIMENSION(jpij) :: zasum, z1_asum, zaksum ! sum of a_i+ato_i and reverse294 REAL(wp), DIMENSION(jpij ,jpl) :: zhi ! ice thickness312 ! REAL(wp), DIMENSION(jpij) :: zasum, z1_asum ! sum of a_i+ato_i and inverse 313 REAL(wp), DIMENSION(jpij) :: z1_asum ! inverse of sum of a_i+ato_i 295 314 REAL(wp), DIMENSION(jpij,-1:jpl) :: zGsum ! zGsum(n) = sum of areas in categories 0 to n 296 315 !-------------------------------------------------------------------- … … 298 317 z1_gstar = 1._wp / rn_gstar 299 318 z1_astar = 1._wp / rn_astar 319 320 pa_i(:,:) = pa_i_in(:,:) 300 321 301 322 ! ! Ice thickness needed for rafting … … 303 324 ELSEWHERE ; zhi(1:npti,:) = 0._wp 304 325 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) 327 WHERE( pa_i(1:npti,:) > epsi10 .AND. zhi(1:npti,:) < rn_himin ) 328 pa_i(1:npti,:) = pa_i(1:npti,:) * zhi(1:npti,:) / rn_himin 329 zhi(1:npti,:) = rn_himin 330 END WHERE 331 305 332 306 333 ! 1) Participation function (apartf): a(h) = b(h).g(h) … … 317 344 ! Compute total area of ice plus open water. 318 345 ! This is in general not equal to one because of divergence during transport 319 zasum(1:npti) = pato_i(1:npti) + SUM( pa_i(1:npti,:), dim=2 ) 320 ! 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 ) 349 321 350 WHERE( zasum(1:npti) > epsi10 ) ; z1_asum(1:npti) = 1._wp / zasum(1:npti) 322 351 ELSEWHERE ; z1_asum(1:npti) = 0._wp 323 352 END WHERE 324 ! 353 325 354 ! Compute cumulative thickness distribution function 326 355 ! Compute the cumulative thickness distribution function zGsum, … … 330 359 zGsum(1:npti,0 ) = pato_i(1:npti) * z1_asum(1:npti) 331 360 DO jl = 1, jpl 332 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) 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) 333 362 END DO 334 363 ! … … 395 424 ENDIF 396 425 426 427 397 428 ! 2) Transfer function 398 429 !----------------------------------------------------------------- … … 413 444 ! (relative to the Hibler formulation) when very thick ice is ridging. 414 445 ! 415 ! zaksum = net area removed/ total area removed416 ! where total area removed= area of ice that ridges417 ! net area removed = total area removed- area of new ridges446 ! zaksum = net area removed/ total area participating 447 ! where total area participating = area of ice that ridges 448 ! net area removed = total area participating - area of new ridges 418 449 !----------------------------------------------------------------- 419 450 zfac = 1._wp / hi_hrft … … 422 453 DO jl = 1, jpl 423 454 DO ji = 1, npti 424 IF ( apartf(ji,jl) > 0._wp ) THEN 455 IF ( apartf(ji,jl) > 0._wp ) THEN 456 ! IF ( pa_i(ji,jl) > epsi10 .and. apartf(ji,jl) > 0._wp ) THEN 425 457 zhmean = MAX( SQRT( rn_hstar * zhi(ji,jl) ), zhi(ji,jl) * hrdg_hi_min ) 426 458 hrmin (ji,jl) = MIN( 2._wp * zhi(ji,jl), 0.5_wp * ( zhmean + zhi(ji,jl) ) ) 427 459 hrmax (ji,jl) = 2._wp * zhmean - hrmin(ji,jl) 460 hrexp (ji,jl) = rn_murdg * SQRT( zhi(ji,jl) ) 428 461 hraft (ji,jl) = zhi(ji,jl) * zfac 429 462 hi_hrdg(ji,jl) = zhi(ji,jl) / MAX( zhmean, epsi20 ) … … 435 468 hrmin (ji,jl) = 0._wp 436 469 hrmax (ji,jl) = 0._wp 470 hrexp (ji,jl) = 0._wp 437 471 hraft (ji,jl) = 0._wp 438 472 hi_hrdg(ji,jl) = 1._wp 473 zaksum (ji) = 0._wp 439 474 ENDIF 440 475 END DO … … 754 789 !! dissipated per unit area removed from the ice pack under compression, 755 790 !! and assumed proportional to the change in potential energy caused 756 !! by ridging. Note that only Hibler's formulation is stable and that757 !! ice strength has to be smoothed 791 !! by ridging. [**UPDATE THIS ---> Note that only Hibler's formulation is stable and that 792 !! ice strength has to be smoothed**] 758 793 !!---------------------------------------------------------------------- 759 794 INTEGER :: ji, jj, jl ! dummy loop indices … … 761 796 INTEGER :: itframe ! number of time steps for the P smoothing 762 797 REAL(wp) :: zp, z1_3 ! local scalars 798 REAL(wp) :: Cp ! proportional constant for PE 799 !REAL(wp) :: hi ! ice thickness of category (m) 800 REAL(wp) :: h2rdg ! mean value of h^2 for new ridge 801 REAL(wp) :: dh2rdg ! change in mean value of h^2 per unit area consumed by ridging 763 802 REAL(wp), DIMENSION(jpi,jpj) :: zworka ! temporary array used here 764 803 REAL(wp), DIMENSION(jpi,jpj) :: zstrp1, zstrp2 ! strength at previous time steps 804 REAL(wp), DIMENSION(jpij,jpl) :: krdg ! mean ridge thickness / thickness of ridging ice 805 REAL(wp), DIMENSION(jpij) :: strength_1d ! 1-dimensional strength 806 REAL(wp), DIMENSION(jpij,jpl) :: strength_pe ! PE component of R75 strength 807 REAL(wp), DIMENSION(jpij) :: strength_pe_sum ! PE component of R75 strength 765 808 !!---------------------------------------------------------------------- 766 809 ! !--------------------------------------------------! … … 769 812 strength(:,:) = rn_pstar * SUM( v_i(:,:,:), dim=3 ) * EXP( -rn_crhg * ( 1._wp - SUM( a_i(:,:,:), dim=3 ) ) ) 770 813 ismooth = 1 814 ! !--------------------------------------------------! 815 ELSE IF( ln_str_R75 ) THEN ! Ice strength => Rothrock (1975) method ! 816 ! !--------------------------------------------------! 817 IF ( kt_ice == nit000 ) THEN 818 819 ! Replace this with input from start file 820 strength(:,:) = rn_pstar * SUM( v_i(:,:,:), dim=3 ) * EXP( -rn_crhg * ( 1._wp - SUM( a_i(:,:,:), dim=3 ) ) ) 821 ismooth = 1 822 823 ELSE 824 825 !-------------------------------- 826 ! 0) Identify grid cells with ice 827 !-------------------------------- 828 at_i(:,:) = SUM( a_i(:,:,:), dim=3 ) 829 vt_i(:,:) = SUM( v_i(:,:,:), dim=3 ) 830 831 ! ato_i(:,:) = 1._wp - at_i(:,:) 832 833 ! DO jj = 1, jpj 834 ! ato_i_1d (ji) = MAX( 0._wp, 1._wp - SUM( at_i(ji,:) ) ) 835 ! END DO 836 ! 837 npti = 0 ; nptidx(:) = 0 838 DO jj = 1, jpj 839 DO ji = 1, jpi 840 IF ( at_i(ji,jj) > epsi10 ) THEN 841 npti = npti + 1 842 nptidx( npti ) = (jj - 1) * jpi + ji 843 ENDIF 844 END DO 845 END DO 846 847 CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i ) 848 CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i ) 849 CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti) , ato_i ) 850 851 852 !----------------------------------- 853 ! Get thickness distribution of ice 854 !----------------------------------- 855 856 CALL rdgrft_prep( a_i_2d, v_i_2d, ato_i_1d, closing_net ) 857 858 ! Don't need closing rate here (unless only calculating strength for ridging ice), could split this function up. 859 860 Cp = 0.5_wp * grav * (rhow-rhoi) * rhoi / rhow ! put this elsewhere, perhaps ice_prep? 861 862 863 IF ( npti > 0 ) THEN 864 865 DO jl = 1, jpl ! categories 866 DO ji = 1, npti ! grid cells 867 868 IF ( a_i_2d(ji,jl) > epsi10 .and. apartf(ji,jl) > 0._wp ) THEN 869 870 ! Move this to rdgrft_prep 871 krdg(ji,jl) = (hrmin(ji,jl) + hrmax(ji,jl))/2._wp * 1._wp/zhi(ji,jl) 872 873 874 ! Add a logical for distribution 875 ! Uniform distribution 876 h2rdg = (1._wp/3._wp) * (hrmax(ji,jl)**3._wp - hrmin(ji,jl)**3._wp) & 877 / (hrmax(ji,jl) - hrmin(ji,jl)) 878 879 880 ! Exponential distribution 881 ! h2rdg = hrmin(ji,jl) * hrmin(ji,jl) & 882 ! + 2._wp * hrmin(ji,jl) * hrexp(ji,jl) & 883 ! + 2._wp * hrexp(ji,jl) * hrexp(ji,jl) 884 885 dh2rdg = -zhi(ji,jl)*zhi(ji,jl) + h2rdg / krdg(ji,jl) 886 887 strength_pe(ji,jl) = apartf(ji,jl) * dh2rdg 888 889 ! strength_pe(ji,jl) = strength_pe(ji,jl) + apartf(ji,jl) * dh2rdg !CICE format 890 891 ELSE 892 893 strength_pe(ji,jl) = 0._wp 894 ! strength_pe(ji,jl) = strength_pe(ji,jl) + 0._wp !CICE format 895 896 END IF 897 898 END DO ! grid cells 899 END DO ! categories 900 901 902 strength_pe_sum(:) = SUM(strength_pe(:,:), dim=2) ! Sum over all categories 903 ! strength_pe_1d(:) = ??? ! CICE format, already summed so how to convert to 1D?? 904 905 DO ji = 1, npti 906 IF ( zaksum(ji) > epsi10 ) THEN 907 strength_1d(ji) = rn_Cf * Cp * strength_pe_sum(ji) / zaksum(ji) 908 ELSE 909 strength_1d(ji) = 0._wp 910 END IF 911 END DO 912 913 ! Convert strength_1d to 2d 914 915 CALL tab_1d_2d( npti, nptidx(1:npti), strength_1d(1:npti), strength(:,:) ) 916 917 ! Enforce a maximum for strength (- Look into why numbers can be so large) 918 WHERE( strength(:,:) > 1.5e5) ; strength(:,:) = 1.5e5 919 END WHERE 920 921 ismooth = 0 ! Buffer size error with lbc_lnk... 922 923 END IF ! number of ice points 924 END IF ! timestep 925 926 771 927 ! !--------------------------------------------------! 772 ELSE ! Zero strength 928 ELSE ! Zero strength (this just crashes) ! 773 929 ! !--------------------------------------------------! 774 930 strength(:,:) = 0._wp 775 931 ismooth = 0 776 932 ENDIF 777 ! !--------------------------------------------------! 933 ! 934 935 !--------------------------------------------------! 778 936 SELECT CASE( ismooth ) ! Smoothing ice strength ! 779 937 ! !--------------------------------------------------! 938 CASE( 0 ) !--- No smoothing 939 940 941 !CALL lbc_lnk( 'icedyn_rdgrft', strength, 'T', 1. ) 942 780 943 CASE( 1 ) !--- Spatial smoothing 781 944 DO jj = 2, jpjm1 … … 915 1078 INTEGER :: ios ! Local integer output status for namelist read 916 1079 !! 917 NAMELIST/namdyn_rdgrft/ ln_str_H79, rn_pstar, rn_crhg, & 918 & rn_csrdg , & 1080 NAMELIST/namdyn_rdgrft/ ln_str_H79, rn_pstar, & 1081 & ln_str_R75, rn_Cf, rn_murdg, & 1082 & rn_crhg, rn_csrdg, & 919 1083 & ln_partf_lin, rn_gstar, & 920 1084 & ln_partf_exp, rn_astar, & … … 936 1100 WRITE(numout,*) '~~~~~~~~~~~~~~~~~~' 937 1101 WRITE(numout,*) ' Namelist namdyn_rdgrft:' 938 WRITE(numout,*) ' ice strength parameterization Hibler (1979) ln_str_H79 = ', ln_str_H79 1102 WRITE(numout,*) ' ice strength parameterization Hibler (1979) ln_str_H79 = ', ln_str_H79 939 1103 WRITE(numout,*) ' 1st bulk-rheology parameter rn_pstar = ', rn_pstar 940 1104 WRITE(numout,*) ' 2nd bulk-rhelogy parameter rn_crhg = ', rn_crhg 1105 WRITE(numout,*) ' ice strength parameterization Rothrock (1975) ln_str_R75 = ', ln_str_R75 1106 WRITE(numout,*) ' ratio of ridging work to PE change in ridging rn_Cf = ', rn_Cf 1107 WRITE(numout,*) ' gives e-folding scale of ridged ice (m^.5) rn_murdg = ', rn_murdg 941 1108 WRITE(numout,*) ' Fraction of shear energy contributing to ridging rn_csrdg = ', rn_csrdg 942 1109 WRITE(numout,*) ' linear ridging participation function ln_partf_lin = ', ln_partf_lin
Note: See TracChangeset
for help on using the changeset viewer.
