Changeset 15647
- Timestamp:
- 2022-01-17T11:04:18+01:00 (2 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/UKMO/NEMO_4.0.4_ice_strength_EAP/src/ICE/icedyn_rdgrft.F90
r15570 r15647 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 REAL(wp) :: rn_pstar ! determines ice strength, Hibler JPO79 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 INTEGER :: ismooth ! smoothing the resistance to deformation (strength) 68 LOGICAL :: ln_redist_ufm ! uniform redistribution of ridged ice (Hibler, 1980) 69 LOGICAL :: ln_redist_exp ! exponential redistribution of ridged ice 70 REAL(wp) :: rn_murdg ! gives e-folding scale of ridged ice (m^.5) for ln_redist_exp 61 71 REAL(wp) :: rn_csrdg ! fraction of shearing energy contributing to ridging 62 72 LOGICAL :: ln_partf_lin ! participation function linear (Thorndike et al. (1975)) … … 86 96 !! *** ROUTINE ice_dyn_rdgrft_alloc *** 87 97 !!------------------------------------------------------------------- 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) , & 98 ALLOCATE( closing_net(jpij) , opning(jpij) , closing_gross(jpij) , & 99 & apartf(jpij,0:jpl) , hrmin(jpij,jpl) , hrmax(jpij,jpl) , zhi(jpij,jpl), & 100 & zaksum(jpij) , hraft(jpij,jpl) , hrexp(jpij,jpl) , & 101 & hi_hrdg(jpij,jpl) , araft(jpij,jpl) , aridge(jpij,jpl) , zasum(jpij), & 91 102 & ze_i_2d(jpij,nlay_i,jpl), ze_s_2d(jpij,nlay_s,jpl), STAT=ice_dyn_rdgrft_alloc ) 92 103 … … 188 199 IF( ln_rhg_EVP ) closing_net(ji) = rn_csrdg * 0.5_wp * ( zdelt(ji) - ABS( zdivu(ji) ) ) - MIN( zdivu(ji), 0._wp ) 189 200 IF( ln_rhg_EAP ) closing_net(ji) = zconv(ji) 190 ! 201 ! 191 202 IF( zdivu(ji) < 0._wp ) closing_net(ji) = MAX( closing_net(ji), -zdivu(ji) ) ! make sure the closing rate is large enough 192 203 ! ! to give asum = 1.0 after ridging 193 204 ! Opening rate (non-negative) that will give asum = 1.0 after ridging. 194 205 opning(ji) = closing_net(ji) + zdivu(ji) 206 ! 195 207 END DO 196 208 ! … … 280 292 281 293 282 SUBROUTINE rdgrft_prep( pa_i, pv_i, pato_i, pclosing_net ) 294 295 SUBROUTINE rdgrft_prep( pa_i_in, pv_i, pato_i, pclosing_net ) 283 296 !!------------------------------------------------------------------- 284 297 !! *** ROUTINE rdgrft_prep *** … … 290 303 !!------------------------------------------------------------------- 291 304 REAL(wp), DIMENSION(:) , INTENT(in) :: pato_i, pclosing_net 292 REAL(wp), DIMENSION(:,:), INTENT(in) :: pa_i, pv_i 305 REAL(wp), DIMENSION(:,:), INTENT(in) :: pa_i_in, pv_i 306 REAL(wp), DIMENSION(jpij,jpl) :: pa_i ! area adjusted for minimum zhi 293 307 !! 294 308 INTEGER :: ji, jl ! dummy loop indices 295 309 REAL(wp) :: z1_gstar, z1_astar, zhmean, zfac ! local scalar 296 REAL(wp), DIMENSION(jpij) :: zasum, z1_asum, zaksum ! sum of a_i+ato_i and reverse 297 REAL(wp), DIMENSION(jpij,jpl) :: zhi ! ice thickness 310 REAL(wp), DIMENSION(jpij) :: z1_asum ! inverse of sum of a_i+ato_i 298 311 REAL(wp), DIMENSION(jpij,-1:jpl) :: zGsum ! zGsum(n) = sum of areas in categories 0 to n 299 312 !-------------------------------------------------------------------- … … 301 314 z1_gstar = 1._wp / rn_gstar 302 315 z1_astar = 1._wp / rn_astar 316 317 pa_i(:,:) = pa_i_in(:,:) 303 318 304 319 ! ! Ice thickness needed for rafting … … 306 321 ELSEWHERE ; zhi(1:npti,:) = 0._wp 307 322 END WHERE 323 ! Make sure ice thickness is not below the minimum (from iceitd.F90; no minimum was causing issues for R75 324 ! strength) 325 WHERE( pa_i(1:npti,:) > epsi10 .AND. zhi(1:npti,:) < rn_himin ) 326 pa_i(1:npti,:) = pa_i(1:npti,:) * zhi(1:npti,:) / rn_himin 327 zhi(1:npti,:) = rn_himin 328 END WHERE 329 330 ! Set tiny values of pa_i to zero 331 ! If this isn't done, can end up with zhi=0 but apartf>0 332 WHERE( pa_i(1:npti,:) <= epsi10 ) ; pa_i(1:npti,:) = 0._wp 333 END WHERE 334 308 335 309 336 ! 1) Participation function (apartf): a(h) = b(h).g(h) … … 320 347 ! Compute total area of ice plus open water. 321 348 ! This is in general not equal to one because of divergence during transport 322 zasum(1:npti) = pato_i(1:npti) + SUM( pa_i(1:npti, :), dim=2 )323 ! 349 zasum(1:npti) = pato_i(1:npti) + SUM( pa_i(1:npti,1:jpl), dim=2 ) 350 324 351 WHERE( zasum(1:npti) > epsi10 ) ; z1_asum(1:npti) = 1._wp / zasum(1:npti) 325 352 ELSEWHERE ; z1_asum(1:npti) = 0._wp 326 353 END WHERE 327 ! 354 328 355 ! Compute cumulative thickness distribution function 329 356 ! Compute the cumulative thickness distribution function zGsum, … … 333 360 zGsum(1:npti,0 ) = pato_i(1:npti) * z1_asum(1:npti) 334 361 DO jl = 1, jpl 335 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)362 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) 336 363 END DO 337 364 ! … … 398 425 ENDIF 399 426 427 428 400 429 ! 2) Transfer function 401 430 !----------------------------------------------------------------- … … 416 445 ! (relative to the Hibler formulation) when very thick ice is ridging. 417 446 ! 418 ! zaksum = net area removed/ total area removed419 ! where total area removed= area of ice that ridges420 ! net area removed = total area removed- area of new ridges447 ! zaksum = net area removed/ total area participating 448 ! where total area participating = area of ice that ridges 449 ! net area removed = total area participating - area of new ridges 421 450 !----------------------------------------------------------------- 422 451 zfac = 1._wp / hi_hrft … … 428 457 zhmean = MAX( SQRT( rn_hstar * zhi(ji,jl) ), zhi(ji,jl) * hrdg_hi_min ) 429 458 hrmin (ji,jl) = MIN( 2._wp * zhi(ji,jl), 0.5_wp * ( zhmean + zhi(ji,jl) ) ) 430 hrmax (ji,jl) = 2._wp * zhmean - hrmin(ji,jl)431 459 hraft (ji,jl) = zhi(ji,jl) * zfac 432 hi_hrdg(ji,jl) = zhi(ji,jl) / MAX( zhmean, epsi20 ) 460 ! 461 IF ( ln_redist_ufm ) THEN 462 hrmax (ji,jl) = 2._wp * zhmean - hrmin(ji,jl) 463 hi_hrdg(ji,jl) = zhi(ji,jl) / zhmean 464 ELSE IF ( ln_redist_exp ) THEN 465 hrexp (ji,jl) = rn_murdg * SQRT( zhi(ji,jl) ) 466 hi_hrdg(ji,jl) = zhi(ji,jl) / ( hrmin(ji,jl) + hrexp(ji,jl) ) 467 END IF 433 468 ! 434 469 ! Normalization factor : zaksum, ensures mass conservation … … 439 474 hrmax (ji,jl) = 0._wp 440 475 hraft (ji,jl) = 0._wp 476 hrexp (ji,jl) = 0._wp 441 477 hi_hrdg(ji,jl) = 1._wp 478 zaksum (ji) = 0._wp 442 479 ENDIF 443 480 END DO … … 493 530 INTEGER :: ji, jj, jl, jl1, jl2, jk ! dummy loop indices 494 531 REAL(wp) :: hL, hR, farea ! left and right limits of integration and new area going to jl2 532 REAL(wp) :: expL, expR ! exponentials involving hL, hR 495 533 REAL(wp) :: vsw ! vol of water trapped into ridges 496 534 REAL(wp) :: afrdg, afrft ! fraction of category area ridged/rafted … … 661 699 itest_rdg(1:npti) = 0 662 700 itest_rft(1:npti) = 0 701 ! 663 702 DO jl2 = 1, jpl 664 703 ! 665 704 DO ji = 1, npti 666 705 667 IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN 668 669 ! Compute the fraction of ridged ice area and volume going to thickness category jl2 706 IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN 707 ! 708 ! Compute the fraction of ridged ice area and volume going to thickness category jl2 709 ! 710 IF ( ln_redist_ufm ) THEN ! Hibler 1980 uniform formulation 711 ! 670 712 IF( hrmin(ji,jl1) <= hi_max(jl2) .AND. hrmax(ji,jl1) > hi_max(jl2-1) ) THEN 671 713 hL = MAX( hrmin(ji,jl1), hi_max(jl2-1) ) … … 679 721 fvol(ji) = 0._wp 680 722 ENDIF 723 ! 724 ELSE IF ( ln_redist_exp ) THEN ! exponential formulation 725 ! 726 IF (jl2 < jpl) THEN 727 ! 728 IF (hrmin(ji,jl1) <= hi_max(jl2)) THEN 729 hL = MAX( hrmin(ji,jl1), hi_max(jl2-1) ) 730 hR = hi_max(jl2) 731 expL = EXP( -( hL - hrmin(ji,jl1) ) / hrexp(ji,jl1) ) 732 expR = EXP( -( hR - hrmin(ji,jl1) ) / hrexp(ji,jl1) ) 733 farea = expL - expR 734 fvol(ji) = ( ( hL + hrexp(ji,jl1) ) * expL & 735 - ( hR + hrexp(ji,jl1) ) * expR ) / ( hrmin(ji,jl1) + hrexp(ji,jl1) ) 736 ELSE 737 farea = 0._wp 738 fvol(ji) = 0._wp 739 END IF 740 ! 741 ELSE ! jl2 = jpl 742 ! 743 hL = MAX( hrmin(ji,jl1), hi_max(jl2-1) ) 744 expL = EXP(-( hL - hrmin(ji,jl1) ) / hrexp(ji,jl1) ) 745 farea = expL 746 fvol(ji) = ( hL + hrexp(ji,jl1) ) * expL / ( hrmin(ji,jl1) + hrexp(ji,jl1) ) 747 ! 748 END IF ! jl2 < jpl 749 ! 750 END IF ! ridge redistribution 681 751 682 752 ! Compute the fraction of rafted ice area and volume going to thickness category jl2 … … 757 827 !! dissipated per unit area removed from the ice pack under compression, 758 828 !! and assumed proportional to the change in potential energy caused 759 !! by ridging. Note that only Hibler's formulation is stable and that760 !! ice strength has to be smoothed829 !! by ridging. Note that ice strength using Hibler's formulation must be 830 !! smoothed. 761 831 !!---------------------------------------------------------------------- 762 832 INTEGER :: ji, jj, jl ! dummy loop indices 763 INTEGER :: ismooth ! smoothing the resistance to deformation764 833 INTEGER :: itframe ! number of time steps for the P smoothing 765 834 REAL(wp) :: zp, z1_3 ! local scalars 835 REAL(wp) :: Cp ! proportional constant for PE 836 REAL(wp) :: h2rdg ! mean value of h^2 for new ridge 837 REAL(wp) :: dh2rdg ! change in mean value of h^2 per unit area consumed by ridging 766 838 REAL(wp), DIMENSION(jpi,jpj) :: zworka ! temporary array used here 767 839 REAL(wp), DIMENSION(jpi,jpj) :: zstrp1, zstrp2 ! strength at previous time steps 840 REAL(wp), DIMENSION(jpij) :: strength_1d ! 1-dimensional strength 841 REAL(wp), DIMENSION(jpij,jpl) :: strength_pe ! PE component of R75 strength 842 REAL(wp), DIMENSION(jpij) :: strength_pe_sum ! PE component of R75 strength 768 843 !!---------------------------------------------------------------------- 769 844 ! !--------------------------------------------------! 770 845 IF( ln_str_H79 ) THEN ! Ice strength => Hibler (1979) method ! 771 846 ! !--------------------------------------------------! 772 strength(:,:) = rn_pstar * SUM( v_i(:,:,:), dim=3 ) * EXP( -rn_crhg * ( 1._wp - SUM( a_i(:,:,:), dim=3 ) ) ) 773 ismooth = 1 847 strength(:,:) = rn_pstar * SUM( v_i(:,:,:), dim=3 ) * EXP( -rn_crhg * ( 1._wp - SUM( a_i(:,:,:), dim=3 ) ) ) 848 ! Recommend ismooth = 1 849 ! !--------------------------------------------------! 850 ELSE IF( ln_str_R75 ) THEN ! Ice strength => Rothrock (1975) method ! 851 ! !--------------------------------------------------! 852 ! 853 ! Initialise the strength field at each timestep 854 strength(:,:) = 0._wp 855 ! Recommend ismooth = 0 856 857 !-------------------------------- 858 ! 0) Identify grid cells with ice 859 !-------------------------------- 860 at_i(:,:) = SUM( a_i(:,:,:), dim=3 ) 861 vt_i(:,:) = SUM( v_i(:,:,:), dim=3 ) 862 ato_i(:,:) = 1 - at_i(:,:) 863 ! 864 npti = 0 ; nptidx(:) = 0 865 DO jj = 1, jpj 866 DO ji = 1, jpi 867 IF ( at_i(ji,jj) > epsi10 ) THEN 868 npti = npti + 1 869 nptidx( npti ) = (jj - 1) * jpi + ji 870 ENDIF 871 END DO 872 END DO 873 874 CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i ) 875 CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i ) 876 CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti) , ato_i ) 877 878 !----------------------------------- 879 ! Get thickness distribution of ice 880 !----------------------------------- 881 882 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 883 884 CALL rdgrft_prep( a_i_2d, v_i_2d, ato_i_1d, closing_net ) 885 ! EF: Note that closing rate is not needed here, could split this function up, and refactor as initial rdgrft_prep is 886 ! called here, before going into ridging iterations 887 888 ! Initialise 889 strength_pe_sum(:) = 0._wp 890 891 IF ( npti > 0 ) THEN 892 893 DO jl = 1, jpl ! categories 894 DO ji = 1, npti ! grid cells 895 IF ( apartf(ji,jl) > 0._wp ) THEN 896 ! 897 IF ( ln_redist_ufm ) THEN ! Uniform distribution of ridged ice 898 h2rdg = (1._wp/3._wp) * (hrmax(ji,jl)**3._wp - hrmin(ji,jl)**3._wp) & 899 / (hrmax(ji,jl) - hrmin(ji,jl)) 900 ! 901 ELSE IF ( ln_redist_exp ) THEN ! Exponential distribution of ridged ice 902 h2rdg = hrmin(ji,jl) * hrmin(ji,jl) & 903 + 2._wp * hrmin(ji,jl) * hrexp(ji,jl) & 904 + 2._wp * hrexp(ji,jl) * hrexp(ji,jl) 905 END IF 906 907 dh2rdg = -zhi(ji,jl)*zhi(ji,jl) + h2rdg * hi_hrdg(ji,jl) 908 909 strength_pe_sum(ji) = strength_pe_sum(ji) + apartf(ji,jl) * dh2rdg 910 911 ELSE 912 913 strength_pe_sum(ji) = strength_pe_sum(ji) + 0._wp 914 915 END IF 916 917 END DO ! grid cells 918 END DO ! categories 919 920 921 DO ji = 1, npti 922 IF ( zaksum(ji) > epsi10 ) THEN ! Unclear why is this check should be needed? Occasional x10E-24... 923 strength_1d(ji) = rn_Cf * Cp * strength_pe_sum(ji) / zaksum(ji) 924 ELSE 925 strength_1d(ji) = 0._wp 926 END IF 927 END DO 928 929 930 ! Enforce a maximum for strength 931 WHERE( strength_1d(:) > 2.0e5) ; strength_1d(:) = 2.0e5 932 END WHERE 933 934 935 ! Convert strength_1d to 2d 936 CALL tab_1d_2d( npti, nptidx(1:npti), strength_1d(1:npti), strength(:,:) ) 937 938 END IF ! number of ice points 939 940 774 941 ! !--------------------------------------------------! 775 ELSE ! Zero strength 942 ELSE ! Zero strength (this just crashes) ! 776 943 ! !--------------------------------------------------! 777 944 strength(:,:) = 0._wp 778 ismooth = 0 779 ENDIF 780 ! !--------------------------------------------------! 945 ! 946 ENDIF ! Ice strength formulation 947 ! 948 !--------------------------------------------------! 781 949 SELECT CASE( ismooth ) ! Smoothing ice strength ! 782 950 ! !--------------------------------------------------! 951 CASE( 0 ) !--- No smoothing 952 953 CALL lbc_lnk( 'icedyn_rdgrft', strength, 'T', 1. ) 954 783 955 CASE( 1 ) !--- Spatial smoothing 784 956 DO jj = 2, jpjm1 … … 918 1090 INTEGER :: ios ! Local integer output status for namelist read 919 1091 !! 920 NAMELIST/namdyn_rdgrft/ ln_str_H79, rn_pstar, rn_crhg, & 921 & rn_csrdg , & 1092 NAMELIST/namdyn_rdgrft/ ismooth, ln_str_H79, rn_pstar, & 1093 & ln_str_R75, rn_Cf, & 1094 & ln_redist_ufm, ln_redist_exp, rn_murdg, & 1095 & rn_crhg, rn_csrdg, & 922 1096 & ln_partf_lin, rn_gstar, & 923 1097 & ln_partf_exp, rn_astar, & … … 939 1113 WRITE(numout,*) '~~~~~~~~~~~~~~~~~~' 940 1114 WRITE(numout,*) ' Namelist namdyn_rdgrft:' 941 WRITE(numout,*) ' ice strength parameterization Hibler (1979) ln_str_H79 = ', ln_str_H79 1115 WRITE(numout,*) ' ice strength parameterization Hibler (1979) ln_str_H79 = ', ln_str_H79 942 1116 WRITE(numout,*) ' 1st bulk-rheology parameter rn_pstar = ', rn_pstar 943 1117 WRITE(numout,*) ' 2nd bulk-rhelogy parameter rn_crhg = ', rn_crhg 1118 WRITE(numout,*) ' ice strength parameterization Rothrock (1975) ln_str_R75 = ', ln_str_R75 1119 WRITE(numout,*) ' ratio of ridging work to PE change in ridging rn_Cf = ', rn_Cf 1120 WRITE(numout,*) ' smoothing the resistance to deformation (strength) ismooth = ', ismooth 1121 WRITE(numout,*) ' uniform redistribution of ridged ice (Hibler, 1980) ln_redist_ufm = ', ln_redist_ufm 1122 WRITE(numout,*) ' exponential redistribution of ridged ice ln_redist_exp = ', ln_redist_exp 1123 WRITE(numout,*) ' e-folding scale of ridged ice for ln_redist_exp rn_murdg = ', rn_murdg 944 1124 WRITE(numout,*) ' Fraction of shear energy contributing to ridging rn_csrdg = ', rn_csrdg 945 1125 WRITE(numout,*) ' linear ridging participation function ln_partf_lin = ', ln_partf_lin … … 959 1139 ENDIF 960 1140 ! 1141 IF ( ( ln_str_H79 .AND. ln_str_R75 ) .OR. ( .NOT.ln_str_H79 .AND. .NOT.ln_str_R75 ) ) THEN 1142 CALL ctl_stop( 'ice_dyn_rdgrft_init: choose one and only one ice strength formulation (ln_str_H79 or ln_str_R75)' ) 1143 ENDIF 1144 ! 961 1145 IF ( ( ln_partf_lin .AND. ln_partf_exp ) .OR. ( .NOT.ln_partf_lin .AND. .NOT.ln_partf_exp ) ) THEN 962 1146 CALL ctl_stop( 'ice_dyn_rdgrft_init: choose one and only one participation function (ln_partf_lin or ln_partf_exp)' ) 963 1147 ENDIF 964 1148 ! 1149 IF ( ( ln_redist_ufm .AND. ln_redist_exp ) .OR. ( .NOT.ln_redist_ufm .AND. .NOT.ln_redist_exp ) ) THEN 1150 CALL ctl_stop( 'ice_dyn_rdgrft_init: choose one and only one ice redistribution function (ln_redist_ufm or ln_redist_exp)' ) 1151 ENDIF 1152 ! 1153 !! Add check for ismooth too, =0, 1 or 2 1154 965 1155 IF( .NOT. ln_icethd ) THEN 966 1156 rn_porordg = 0._wp
Note: See TracChangeset
for help on using the changeset viewer.
