- Timestamp:
- 2021-12-03T20:32:50+01:00 (3 years ago)
- Location:
- NEMO/branches/2021/dev_r14318_RK3_stage1
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/ICE/icedyn_rdgrft.F90 (modified) (16 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2021/dev_r14318_RK3_stage1
- Property svn:externals
-
old new 9 9 10 10 # SETTE 11 ^/utils/CI/sette@14244 sette 11 ^/utils/CI/sette@HEAD sette 12
-
- Property svn:externals
-
NEMO/branches/2021/dev_r14318_RK3_stage1/src/ICE/icedyn_rdgrft.F90
r14072 r15574 38 38 PUBLIC ice_strength ! called by icedyn_rhg_evp 39 39 40 INTEGER :: nice_str ! choice of the type of strength 41 ! ! associated indices: 42 INTEGER, PARAMETER :: np_strh79 = 1 ! Hibler 79 43 INTEGER, PARAMETER :: np_strr75 = 2 ! Rothrock 75 44 INTEGER, PARAMETER :: np_strcst = 3 ! Constant value 45 40 46 ! Variables shared among ridging subroutines 41 47 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: closing_net ! net rate at which area is removed (1/s) … … 57 63 ! 58 64 ! ** namelist (namdyn_rdgrft) ** 59 LOGICAL :: ln_str_H79 ! ice strength parameterization (Hibler79) 60 REAL(wp) :: rn_pstar ! determines ice strength, Hibler JPO79 65 LOGICAL :: ln_str_R75 ! ice strength parameterization: Rothrock 75 66 REAL(wp) :: rn_pe_rdg ! coef accounting for frictional dissipation 67 LOGICAL :: ln_str_CST ! ice strength parameterization: Constant 68 REAL(wp) :: rn_str ! constant value of ice strength 69 LOGICAL :: ln_str_smooth ! ice strength spatial smoothing 70 LOGICAL :: ln_distf_lin ! redistribution of ridged ice: linear (Hibler 1980) 71 LOGICAL :: ln_distf_exp ! redistribution of ridged ice: exponential 72 REAL(wp) :: rn_murdg ! gives e-folding scale of ridged ice (m^.5) 61 73 REAL(wp) :: rn_csrdg ! fraction of shearing energy contributing to ridging 62 74 LOGICAL :: ln_partf_lin ! participation function linear (Thorndike et al. (1975)) … … 162 174 npti = 0 ; nptidx(:) = 0 163 175 ipti = 0 ; iptidx(:) = 0 164 DO_2D( 1, 1, 1, 1)176 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 165 177 IF ( at_i(ji,jj) > epsi10 ) THEN 166 178 npti = npti + 1 … … 272 284 273 285 ! controls 274 IF( sn_cfctl%l_prtctl ) CALL ice_prt3D ('icedyn_rdgrft')! prints286 IF( sn_cfctl%l_prtctl ) CALL ice_prt3D('icedyn_rdgrft') ! prints 275 287 IF( ln_icectl ) CALL ice_prt (kt, iiceprt, jiceprt,-1, ' - ice dyn rdgrft - ') ! prints 276 288 IF( ln_icediachk ) CALL ice_cons_hsm(1, 'icedyn_rdgrft', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation … … 520 532 ! 521 533 INTEGER , DIMENSION(jpij) :: itest_rdg, itest_rft ! test for conservation 534 LOGICAL , DIMENSION(jpij) :: ll_shift ! logical for doing calculation or not 522 535 !!------------------------------------------------------------------- 523 536 ! … … 540 553 DO ji = 1, npti 541 554 542 IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN ! only if ice is ridging 555 ! set logical to true when ridging 556 IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN ; ll_shift(ji) = .TRUE. 557 ELSE ; ll_shift(ji) = .FALSE. 558 ENDIF 559 560 IF( ll_shift(ji) ) THEN ! only if ice is ridging 543 561 544 562 IF( a_i_2d(ji,jl1) > epsi10 ) THEN ; z1_ai(ji) = 1._wp / a_i_2d(ji,jl1) … … 630 648 DO jk = 1, nlay_s 631 649 DO ji = 1, npti 632 IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp) THEN650 IF( ll_shift(ji) ) THEN 633 651 ! Compute ridging /rafting fractions 634 652 afrdg = aridge(ji,jl1) * closing_gross(ji) * rDt_ice * z1_ai(ji) … … 651 669 DO jk = 1, nlay_i 652 670 DO ji = 1, npti 653 IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp) THEN671 IF( ll_shift(ji) ) THEN 654 672 ! Compute ridging /rafting fractions 655 673 afrdg = aridge(ji,jl1) * closing_gross(ji) * rDt_ice * z1_ai(ji) … … 674 692 DO ji = 1, npti 675 693 676 IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp) THEN694 IF( ll_shift(ji) ) THEN 677 695 678 696 ! Compute the fraction of ridged ice area and volume going to thickness category jl2 … … 731 749 DO jk = 1, nlay_s 732 750 DO ji = 1, npti 733 IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp) &751 IF( ll_shift(ji) ) & 734 752 & ze_s_2d(ji,jk,jl2) = ze_s_2d(ji,jk,jl2) + ( esrdg(ji,jk) * rn_fsnwrdg * fvol(ji) + & 735 753 & esrft(ji,jk) * rn_fsnwrft * zswitch(ji) ) … … 740 758 DO jk = 1, nlay_i 741 759 DO ji = 1, npti 742 IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp) &760 IF( ll_shift(ji) ) & 743 761 & ze_i_2d(ji,jk,jl2) = ze_i_2d(ji,jk,jl2) + eirdg(ji,jk) * fvol(ji) + eirft(ji,jk) * zswitch(ji) 744 762 END DO … … 770 788 !!---------------------------------------------------------------------- 771 789 INTEGER :: ji, jj, jl ! dummy loop indices 772 INTEGER :: ismooth ! smoothing the resistance to deformation 773 INTEGER :: itframe ! number of time steps for the P smoothing 774 REAL(wp) :: zp, z1_3 ! local scalars 790 REAL(wp) :: z1_3 ! local scalars 775 791 REAL(wp), DIMENSION(jpi,jpj) :: zworka ! temporary array used here 776 REAL(wp), DIMENSION(jpi,jpj) :: zstrp1, zstrp2 ! strength at previous time steps 792 !!clem 793 LOGICAL :: ln_str_R75 794 REAL(wp) :: zhi, zcp, rn_pe_rdg 795 REAL(wp), DIMENSION(jpij) :: zstrength, zaksum ! strength in 1D 777 796 !!---------------------------------------------------------------------- 778 ! !--------------------------------------------------! 779 IF( ln_str_H79 ) THEN ! Ice strength => Hibler (1979) method ! 780 ! !--------------------------------------------------! 797 ! 798 SELECT CASE( nice_str ) !--- Set which ice strength is chosen 799 800 CASE ( np_strr75 ) !== Rothrock(1975)'s method ==! 801 802 ! these 2 param should be defined once for all at the 1st time step 803 zcp = 0.5_wp * grav * (rho0-rhoi) * rhoi * r1_rho0 ! proport const for PE 804 ! 805 strength(:,:) = 0._wp 806 ! 807 ! Identify grid cells with ice 808 at_i(:,:) = SUM( a_i, dim=3 ) 809 npti = 0 ; nptidx(:) = 0 810 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 811 IF ( at_i(ji,jj) > epsi10 ) THEN 812 npti = npti + 1 813 nptidx( npti ) = (jj - 1) * jpi + ji 814 ENDIF 815 END_2D 816 817 IF( npti > 0 ) THEN 818 CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i ) 819 CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i ) 820 CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti) , ato_i ) 821 CALL tab_2d_1d( npti, nptidx(1:npti), zstrength(1:npti) , strength ) 822 823 CALL rdgrft_prep( a_i_2d, v_i_2d, ato_i_1d, closing_net ) 824 ! 825 zaksum(1:npti) = apartf(1:npti,0) !clem: aksum should be defined in the header => local to module 826 DO jl = 1, jpl 827 DO ji = 1, npti 828 IF ( apartf(ji,jl) > 0._wp ) THEN 829 zaksum(ji) = zaksum(ji) + aridge(ji,jl) * ( 1._wp - hi_hrdg(ji,jl) ) & 830 & + araft (ji,jl) * ( 1._wp - hi_hrft ) 831 ENDIF 832 END DO 833 END DO 834 ! 835 z1_3 = 1._wp / 3._wp 836 DO jl = 1, jpl 837 DO ji = 1, npti 838 ! 839 IF( apartf(ji,jl) > 0._wp ) THEN 840 ! 841 IF( a_i_2d(ji,jl) > epsi10 ) THEN ; zhi = v_i_2d(ji,jl) / a_i_2d(ji,jl) 842 ELSE ; zhi = 0._wp 843 ENDIF 844 zstrength(ji) = zstrength(ji) - apartf(ji,jl) * zhi * zhi ! PE loss from deforming ice 845 zstrength(ji) = zstrength(ji) + 2._wp * araft (ji,jl) * zhi * zhi ! PE gain from rafting ice 846 zstrength(ji) = zstrength(ji) + aridge(ji,jl) * hi_hrdg(ji,jl) * z1_3 * & ! PE gain from ridging ice 847 & ( hrmax(ji,jl) * hrmax(ji,jl) + & ! (a**3-b**3)/(a-b) = a*a+ab+b*b 848 & hrmin(ji,jl) * hrmin(ji,jl) + & 849 & hrmax(ji,jl) * hrmin(ji,jl) ) 850 ENDIF 851 ! 852 END DO 853 END DO 854 ! 855 zstrength(1:npti) = rn_pe_rdg * zcp * zstrength(1:npti) / zaksum(1:npti) 856 ! 857 CALL tab_1d_2d( npti, nptidx(1:npti), zstrength(1:npti), strength ) 858 ! 859 ENDIF 860 ! 861 CASE ( np_strh79 ) !== Hibler(1979)'s method ==! 781 862 strength(:,:) = rn_pstar * SUM( v_i(:,:,:), dim=3 ) * EXP( -rn_crhg * ( 1._wp - SUM( a_i(:,:,:), dim=3 ) ) ) 782 ismooth = 1 ! original code 783 ! ismooth = 0 ! try for EAP stability 784 ! !--------------------------------------------------! 785 ELSE ! Zero strength ! 786 ! !--------------------------------------------------! 787 strength(:,:) = 0._wp 788 ismooth = 0 789 ENDIF 790 ! !--------------------------------------------------! 791 SELECT CASE( ismooth ) ! Smoothing ice strength ! 792 ! !--------------------------------------------------! 793 CASE( 1 ) !--- Spatial smoothing 863 ! 864 CASE ( np_strcst ) !== Constant strength ==! 865 strength(:,:) = rn_str 866 ! 867 END SELECT 868 ! 869 IF( ln_str_smooth ) THEN !--- Spatial smoothing 794 870 DO_2D( 0, 0, 0, 0 ) 795 871 IF ( SUM( a_i(ji,jj,:) ) > 0._wp ) THEN 796 zworka(ji,jj) = ( 4. 0* strength(ji,jj) &797 & + strength(ji-1,jj) * tmask(ji-1,jj,1) + strength(ji+1,jj) * tmask(ji+1,jj,1) &798 & + strength(ji,jj-1) * tmask(ji,jj-1,1) + strength(ji,jj+1) * tmask(ji,jj+1,1) &799 & ) / ( 4. 0+ tmask(ji-1,jj,1) + tmask(ji+1,jj,1) + tmask(ji,jj-1,1) + tmask(ji,jj+1,1) )872 zworka(ji,jj) = ( 4._wp * strength(ji,jj) & 873 & + strength(ji-1,jj) * tmask(ji-1,jj,1) + strength(ji+1,jj) * tmask(ji+1,jj,1) & 874 & + strength(ji,jj-1) * tmask(ji,jj-1,1) + strength(ji,jj+1) * tmask(ji,jj+1,1) & 875 & ) / ( 4._wp + tmask(ji-1,jj,1) + tmask(ji+1,jj,1) + tmask(ji,jj-1,1) + tmask(ji,jj+1,1) ) 800 876 ELSE 801 877 zworka(ji,jj) = 0._wp … … 808 884 CALL lbc_lnk( 'icedyn_rdgrft', strength, 'T', 1.0_wp ) 809 885 ! 810 CASE( 2 ) !--- Temporal smoothing 811 IF ( kt_ice == nit000 ) THEN 812 zstrp1(:,:) = 0._wp 813 zstrp2(:,:) = 0._wp 814 ENDIF 815 ! 816 DO_2D( 0, 0, 0, 0 ) 817 IF ( SUM( a_i(ji,jj,:) ) > 0._wp ) THEN 818 itframe = 1 ! number of time steps for the running mean 819 IF ( zstrp1(ji,jj) > 0._wp ) itframe = itframe + 1 820 IF ( zstrp2(ji,jj) > 0._wp ) itframe = itframe + 1 821 zp = ( strength(ji,jj) + zstrp1(ji,jj) + zstrp2(ji,jj) ) / itframe 822 zstrp2 (ji,jj) = zstrp1 (ji,jj) 823 zstrp1 (ji,jj) = strength(ji,jj) 824 strength(ji,jj) = zp 825 ENDIF 826 END_2D 827 CALL lbc_lnk( 'icedyn_rdgrft', strength, 'T', 1.0_wp ) 828 ! 829 END SELECT 886 ENDIF 830 887 ! 831 888 END SUBROUTINE ice_strength … … 920 977 !! ** input : Namelist namdyn_rdgrft 921 978 !!------------------------------------------------------------------- 922 INTEGER :: ios ! Local integer output status for namelist read 923 !! 924 NAMELIST/namdyn_rdgrft/ ln_str_H79, rn_pstar, rn_crhg, & 925 & rn_csrdg , & 926 & ln_partf_lin, rn_gstar, & 927 & ln_partf_exp, rn_astar, & 979 INTEGER :: ios, ioptio ! Local integer output status for namelist read 980 !! 981 NAMELIST/namdyn_rdgrft/ ln_str_H79, rn_pstar, rn_crhg, ln_str_R75, rn_pe_rdg, ln_str_CST, rn_str, ln_str_smooth, & 982 & ln_distf_lin, ln_distf_exp, rn_murdg, rn_csrdg, & 983 & ln_partf_lin, rn_gstar, ln_partf_exp, rn_astar, & 928 984 & ln_ridging, rn_hstar, rn_porordg, rn_fsnwrdg, rn_fpndrdg, & 929 985 & ln_rafting, rn_hraft, rn_craft , rn_fsnwrft, rn_fpndrft … … 944 1000 WRITE(numout,*) ' 1st bulk-rheology parameter rn_pstar = ', rn_pstar 945 1001 WRITE(numout,*) ' 2nd bulk-rhelogy parameter rn_crhg = ', rn_crhg 1002 WRITE(numout,*) ' ice strength parameterization Rothrock (1975) ln_str_R75 = ', ln_str_R75 1003 WRITE(numout,*) ' coef accounting for frictional dissipation rn_pe_rdg = ', rn_pe_rdg 1004 WRITE(numout,*) ' ice strength parameterization Constant ln_str_CST = ', ln_str_CST 1005 WRITE(numout,*) ' ice strength value rn_str = ', rn_str 1006 WRITE(numout,*) ' spatial smoothing of the strength ln_str_smooth= ', ln_str_smooth 1007 WRITE(numout,*) ' redistribution of ridged ice: linear (Hibler 1980) ln_distf_lin = ', ln_distf_lin 1008 WRITE(numout,*) ' redistribution of ridged ice: exponential ln_distf_exp = ', ln_distf_exp 1009 WRITE(numout,*) ' e-folding scale of ridged ice rn_murdg = ', rn_murdg 946 1010 WRITE(numout,*) ' Fraction of shear energy contributing to ridging rn_csrdg = ', rn_csrdg 947 1011 WRITE(numout,*) ' linear ridging participation function ln_partf_lin = ', ln_partf_lin … … 961 1025 ENDIF 962 1026 ! 1027 ioptio = 0 1028 IF( ln_str_H79 ) THEN ; ioptio = ioptio + 1 ; nice_str = np_strh79 ; ENDIF 1029 IF( ln_str_R75 ) THEN ; ioptio = ioptio + 1 ; nice_str = np_strr75 ; ENDIF 1030 IF( ln_str_CST ) THEN ; ioptio = ioptio + 1 ; nice_str = np_strcst ; ENDIF 1031 IF( ioptio /= 1 ) CALL ctl_stop( 'ice_dyn_rdgrft_init: one and only one ice strength option has to be defined ' ) 1032 ! 1033 IF ( ( ln_distf_lin .AND. ln_distf_exp ) .OR. ( .NOT.ln_distf_lin .AND. .NOT.ln_distf_exp ) ) THEN 1034 CALL ctl_stop( 'ice_dyn_rdgrft_init: choose one and only one redistribution function (ln_distf_lin or ln_distf_exp)' ) 1035 ENDIF 1036 !!clem 1037 IF( ln_distf_exp ) CALL ctl_stop( 'ice_dyn_rdgrft_init: exponential redistribution function not coded yet (ln_distf_exp)' ) 1038 ! 963 1039 IF ( ( ln_partf_lin .AND. ln_partf_exp ) .OR. ( .NOT.ln_partf_lin .AND. .NOT.ln_partf_exp ) ) THEN 964 1040 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.
