Changeset 7646 for trunk/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90
- Timestamp:
- 2017-02-06T10:25:03+01:00 (7 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90 (modified) (17 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90
r6470 r7646 18 18 USE thd_ice ! LIM thermodynamics 19 19 USE ice ! LIM variables 20 USE dom_ice ! LIM domain21 20 USE limvar ! LIM 22 21 USE lbclnk ! lateral boundary condition - MPP exchanges 23 22 USE lib_mpp ! MPP library 24 23 USE wrk_nemo ! work arrays 25 USE prtctl ! Print control26 24 27 25 USE in_out_manager ! I/O manager 28 26 USE iom ! I/O manager 29 27 USE lib_fortran ! glob_sum 30 USE limdiahsb31 28 USE timing ! Timing 32 29 USE limcons ! conservation tests 30 USE limctl ! control prints 33 31 34 32 IMPLICIT NONE … … 70 68 !! *** ROUTINE lim_itd_me_alloc *** 71 69 !!---------------------------------------------------------------------! 72 ALLOCATE( &70 ALLOCATE( & 73 71 !* Variables shared among ridging subroutines 74 & asum (jpi,jpj) , athorn(jpi,jpj,0:jpl) , & 75 & aksum(jpi,jpj) , & 76 & hrmin(jpi,jpj,jpl) , hraft(jpi,jpj,jpl) , aridge(jpi,jpj,jpl) , & 77 & hrmax(jpi,jpj,jpl) , krdg (jpi,jpj,jpl) , araft (jpi,jpj,jpl) , STAT=lim_itd_me_alloc ) 72 & asum (jpi,jpj) , athorn(jpi,jpj,0:jpl) , aksum (jpi,jpj) , & 73 & hrmin(jpi,jpj,jpl) , hraft(jpi,jpj,jpl) , aridge(jpi,jpj,jpl) , & 74 & hrmax(jpi,jpj,jpl) , krdg (jpi,jpj,jpl) , araft (jpi,jpj,jpl) , STAT=lim_itd_me_alloc ) 78 75 ! 79 76 IF( lim_itd_me_alloc /= 0 ) CALL ctl_warn( 'lim_itd_me_alloc: failed to allocate arrays' ) … … 127 124 CALL wrk_alloc( jpi,jpj, closing_net, divu_adv, opning, closing_gross ) 128 125 129 IF(ln_ctl) THEN130 CALL prt_ctl(tab2d_1=ato_i , clinfo1=' lim_itd_me: ato_i : ', tab2d_2=at_i , clinfo2=' at_i : ')131 CALL prt_ctl(tab2d_1=divu_i, clinfo1=' lim_itd_me: divu_i : ', tab2d_2=delta_i, clinfo2=' delta_i : ')132 ENDIF133 134 IF( ln_limdyn ) THEN ! Start ridging and rafting !135 136 126 ! conservation test 137 IF( ln_limdia hsb) CALL lim_cons_hsm(0, 'limitd_me', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)127 IF( ln_limdiachk ) CALL lim_cons_hsm(0, 'limitd_me', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 138 128 139 129 !-----------------------------------------------------------------------------! … … 211 201 DO ji = 1, jpi 212 202 za = ( opning(ji,jj) - athorn(ji,jj,0) * closing_gross(ji,jj) ) * rdt_ice 213 IF ( za < 0._wp .AND. za > - ato_i(ji,jj) ) THEN! would lead to negative ato_i214 zfac = - ato_i(ji,jj) / za203 IF ( za < 0._wp .AND. za > - ato_i(ji,jj) ) THEN ! would lead to negative ato_i 204 zfac = - ato_i(ji,jj) / za 215 205 opning(ji,jj) = athorn(ji,jj,0) * closing_gross(ji,jj) - ato_i(ji,jj) * r1_rdtice 216 206 ELSEIF( za > 0._wp .AND. za > ( asum(ji,jj) - ato_i(ji,jj) ) ) THEN ! would lead to ato_i > asum 217 zfac = ( asum(ji,jj) - ato_i(ji,jj) ) / za207 zfac = ( asum(ji,jj) - ato_i(ji,jj) ) / za 218 208 opning(ji,jj) = athorn(ji,jj,0) * closing_gross(ji,jj) + ( asum(ji,jj) - ato_i(ji,jj) ) * r1_rdtice 219 209 ENDIF … … 259 249 closing_net(ji,jj) = 0._wp 260 250 opning (ji,jj) = 0._wp 251 ato_i (ji,jj) = MAX( 0._wp, 1._wp - SUM( a_i(ji,jj,:) ) ) 261 252 ELSE 262 253 iterate_ridging = 1 … … 292 283 ! control prints 293 284 !-----------------------------------------------------------------------------! 294 IF(ln_ctl) THEN295 CALL lim_var_glo2eqv296 297 CALL prt_ctl_info(' ')298 CALL prt_ctl_info(' - Cell values : ')299 CALL prt_ctl_info(' ~~~~~~~~~~~~~ ')300 CALL prt_ctl(tab2d_1=e1e2t , clinfo1=' lim_itd_me : cell area :')301 CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_itd_me : at_i :')302 CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_itd_me : vt_i :')303 CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_itd_me : vt_s :')304 DO jl = 1, jpl305 CALL prt_ctl_info(' ')306 CALL prt_ctl_info(' - Category : ', ivar1=jl)307 CALL prt_ctl_info(' ~~~~~~~~~~')308 CALL prt_ctl(tab2d_1=a_i (:,:,jl) , clinfo1= ' lim_itd_me : a_i : ')309 CALL prt_ctl(tab2d_1=ht_i (:,:,jl) , clinfo1= ' lim_itd_me : ht_i : ')310 CALL prt_ctl(tab2d_1=ht_s (:,:,jl) , clinfo1= ' lim_itd_me : ht_s : ')311 CALL prt_ctl(tab2d_1=v_i (:,:,jl) , clinfo1= ' lim_itd_me : v_i : ')312 CALL prt_ctl(tab2d_1=v_s (:,:,jl) , clinfo1= ' lim_itd_me : v_s : ')313 CALL prt_ctl(tab2d_1=e_s (:,:,1,jl) , clinfo1= ' lim_itd_me : e_s : ')314 CALL prt_ctl(tab2d_1=t_su (:,:,jl) , clinfo1= ' lim_itd_me : t_su : ')315 CALL prt_ctl(tab2d_1=t_s (:,:,1,jl) , clinfo1= ' lim_itd_me : t_snow : ')316 CALL prt_ctl(tab2d_1=sm_i (:,:,jl) , clinfo1= ' lim_itd_me : sm_i : ')317 CALL prt_ctl(tab2d_1=smv_i (:,:,jl) , clinfo1= ' lim_itd_me : smv_i : ')318 DO jk = 1, nlay_i319 CALL prt_ctl_info(' ')320 CALL prt_ctl_info(' - Layer : ', ivar1=jk)321 CALL prt_ctl_info(' ~~~~~~~')322 CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_itd_me : t_i : ')323 CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_itd_me : e_i : ')324 END DO325 END DO326 ENDIF327 328 285 ! conservation test 329 IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limitd_me', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 330 331 ENDIF ! ln_limdyn=.true. 332 ! 286 IF( ln_limdiachk ) CALL lim_cons_hsm(1, 'limitd_me', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 287 288 ! control prints 289 IF( ln_ctl ) CALL lim_prt3D( 'limitd_me' ) 290 333 291 CALL wrk_dealloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross ) 334 292 ! … … 368 326 rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) ) 369 327 ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch 370 328 END DO 371 329 END DO 372 330 END DO … … 438 396 ENDIF 439 397 440 IF( ln_rafting ) THEN ! Ridging and rafting ice participation functions 398 ! --- Ridging and rafting participation concentrations --- ! 399 IF( ln_rafting .AND. ln_ridging ) THEN 441 400 ! 442 401 DO jl = 1, jpl … … 445 404 zdummy = TANH ( rn_craft * ( ht_i(ji,jj,jl) - rn_hraft ) ) 446 405 aridge(ji,jj,jl) = ( 1._wp + zdummy ) * 0.5_wp * athorn(ji,jj,jl) 447 araft (ji,jj,jl) = ( 1._wp - zdummy ) * 0.5_wp * athorn(ji,jj,jl)406 araft (ji,jj,jl) = athorn(ji,jj,jl) - aridge(ji,jj,jl) 448 407 END DO 449 408 END DO 450 409 END DO 451 452 ELSE 410 ! 411 ELSEIF( ln_ridging .AND. .NOT. ln_rafting ) THEN 453 412 ! 454 413 DO jl = 1, jpl 455 414 aridge(:,:,jl) = athorn(:,:,jl) 415 END DO 416 ! 417 ELSEIF( ln_rafting .AND. .NOT. ln_ridging ) THEN 418 ! 419 DO jl = 1, jpl 420 araft(:,:,jl) = athorn(:,:,jl) 456 421 END DO 457 422 ! … … 657 622 & - sm_i(ji,jj,jl1) * vsw(ij) * rhoic * r1_rdtice ! and get sm_i from the ocean 658 623 ENDIF 659 624 660 625 !------------------------------------------ 661 626 ! 3.7 Put the snow somewhere in the ocean … … 795 760 INTEGER :: numts_rm ! number of time steps for the P smoothing 796 761 REAL(wp) :: zp, z1_3 ! local scalars 797 REAL(wp), POINTER, DIMENSION(:,:) :: zworka ! temporary array used here 762 REAL(wp), POINTER, DIMENSION(:,:) :: zworka ! temporary array used here 763 REAL(wp), POINTER, DIMENSION(:,:) :: zstrp1, zstrp2 ! strength at previous time steps 798 764 !!---------------------------------------------------------------------- 799 765 800 CALL wrk_alloc( jpi, jpj, zworka)766 CALL wrk_alloc( jpi,jpj, zworka, zstrp1, zstrp2 ) 801 767 802 768 !------------------------------------------------------------------------------! … … 844 810 END DO 845 811 846 strength(:,:) = rn_pe_rdg * Cp * strength(:,:) / aksum(:,:) 812 strength(:,:) = rn_pe_rdg * Cp * strength(:,:) / aksum(:,:) * tmask(:,:,1) 847 813 ! where Cp = (g/2)*(rhow-rhoi)*(rhoi/rhow) and rn_pe_rdg accounts for frictional dissipation 848 814 ksmooth = 1 849 815 850 851 852 816 !------------------------------------------------------------------------------! 817 ! 4) Hibler (1979)' method 818 !------------------------------------------------------------------------------! 853 819 ELSE ! kstrngth ne 1: Hibler (1979) form 854 820 ! 855 strength(:,:) = rn_pstar * vt_i(:,:) * EXP( - rn_crhg * ( 1._wp - at_i(:,:) ) ) 821 strength(:,:) = rn_pstar * vt_i(:,:) * EXP( - rn_crhg * ( 1._wp - at_i(:,:) ) ) * tmask(:,:,1) 856 822 ! 857 823 ksmooth = 1 … … 866 832 DO jj = 1, jpj 867 833 DO ji = 1, jpi 868 strength(ji,jj) = strength(ji,jj) * exp(-5.88*SQRT(MAX(bv _i(ji,jj),0.0)))834 strength(ji,jj) = strength(ji,jj) * exp(-5.88*SQRT(MAX(bvm_i(ji,jj),0.0))) 869 835 END DO 870 836 END DO … … 880 846 IF ( ksmooth == 1 ) THEN 881 847 882 CALL lbc_lnk( strength, 'T', 1. )883 884 848 DO jj = 2, jpjm1 885 849 DO ji = 2, jpim1 886 IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp ) THEN850 IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp ) THEN 887 851 zworka(ji,jj) = ( 4.0 * strength(ji,jj) & 888 852 & + strength(ji-1,jj) * tmask(ji-1,jj,1) + strength(ji+1,jj) * tmask(ji+1,jj,1) & … … 907 871 ! Temporal smoothing 908 872 !-------------------- 909 IF ( numit == nit000 + nn_fsbc - 1 ) THEN910 strp1(:,:) = 0.0911 strp2(:,:) = 0.0912 ENDIF913 914 873 IF ( ksmooth == 2 ) THEN 915 874 916 CALL lbc_lnk( strength, 'T', 1. ) 917 918 DO jj = 1, jpj - 1 919 DO ji = 1, jpi - 1 920 IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp) THEN 875 IF ( numit == nit000 + nn_fsbc - 1 ) THEN 876 zstrp1(:,:) = 0._wp 877 zstrp2(:,:) = 0._wp 878 ENDIF 879 880 DO jj = 2, jpjm1 881 DO ji = 2, jpim1 882 IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp ) THEN 921 883 numts_rm = 1 ! number of time steps for the running mean 922 IF ( strp1(ji,jj) > 0._wp ) numts_rm = numts_rm + 1923 IF ( strp2(ji,jj) > 0._wp ) numts_rm = numts_rm + 1924 zp = ( strength(ji,jj) + strp1(ji,jj) +strp2(ji,jj) ) / numts_rm925 strp2(ji,jj) =strp1(ji,jj)926 strp1(ji,jj) = strength(ji,jj)884 IF ( zstrp1(ji,jj) > 0._wp ) numts_rm = numts_rm + 1 885 IF ( zstrp2(ji,jj) > 0._wp ) numts_rm = numts_rm + 1 886 zp = ( strength(ji,jj) + zstrp1(ji,jj) + zstrp2(ji,jj) ) / numts_rm 887 zstrp2(ji,jj) = zstrp1(ji,jj) 888 zstrp1(ji,jj) = strength(ji,jj) 927 889 strength(ji,jj) = zp 928 929 890 ENDIF 930 891 END DO 931 892 END DO 932 893 894 CALL lbc_lnk( strength, 'T', 1. ) ! Boundary conditions 895 933 896 ENDIF ! ksmooth 934 897 935 CALL lbc_lnk( strength, 'T', 1. ) ! Boundary conditions 936 937 CALL wrk_dealloc( jpi, jpj, zworka ) 898 CALL wrk_dealloc( jpi,jpj, zworka, zstrp1, zstrp2 ) 938 899 ! 939 900 END SUBROUTINE lim_itd_me_icestrength … … 953 914 !!------------------------------------------------------------------- 954 915 INTEGER :: ios ! Local integer output status for namelist read 955 NAMELIST/namiceitdme/ rn_cs, rn_fsnowrdg, rn_fsnowrft, & 956 & rn_gstar, rn_astar, rn_hstar, ln_rafting, rn_hraft, rn_craft, rn_por_rdg, & 957 & nn_partfun 916 NAMELIST/namiceitdme/ rn_cs, nn_partfun, rn_gstar, rn_astar, & 917 & ln_ridging, rn_hstar, rn_por_rdg, rn_fsnowrdg, ln_rafting, rn_hraft, rn_craft, rn_fsnowrft 958 918 !!------------------------------------------------------------------- 959 919 ! … … 969 929 IF (lwp) THEN ! control print 970 930 WRITE(numout,*) 971 WRITE(numout,*)' 972 WRITE(numout,*)' 931 WRITE(numout,*)'lim_itd_me_init : ice parameters for mechanical ice redistribution ' 932 WRITE(numout,*)'~~~~~~~~~~~~~~~' 973 933 WRITE(numout,*)' Fraction of shear energy contributing to ridging rn_cs = ', rn_cs 974 WRITE(numout,*)' Fraction of snow volume conserved during ridging rn_fsnowrdg = ', rn_fsnowrdg 975 WRITE(numout,*)' Fraction of snow volume conserved during ridging rn_fsnowrft = ', rn_fsnowrft 934 WRITE(numout,*)' Switch for part. function (0) linear (1) exponential nn_partfun = ', nn_partfun 976 935 WRITE(numout,*)' Fraction of total ice coverage contributing to ridging rn_gstar = ', rn_gstar 977 936 WRITE(numout,*)' Equivalent to G* for an exponential part function rn_astar = ', rn_astar 937 WRITE(numout,*)' Ridging of ice sheets or not ln_ridging = ', ln_ridging 978 938 WRITE(numout,*)' Quantity playing a role in max ridged ice thickness rn_hstar = ', rn_hstar 939 WRITE(numout,*)' Initial porosity of ridges rn_por_rdg = ', rn_por_rdg 940 WRITE(numout,*)' Fraction of snow volume conserved during ridging rn_fsnowrdg = ', rn_fsnowrdg 979 941 WRITE(numout,*)' Rafting of ice sheets or not ln_rafting = ', ln_rafting 980 942 WRITE(numout,*)' Parmeter thickness (threshold between ridge-raft) rn_hraft = ', rn_hraft 981 943 WRITE(numout,*)' Rafting hyperbolic tangent coefficient rn_craft = ', rn_craft 982 WRITE(numout,*)' Initial porosity of ridges rn_por_rdg = ', rn_por_rdg 983 WRITE(numout,*)' Switch for part. function (0) linear (1) exponential nn_partfun = ', nn_partfun 944 WRITE(numout,*)' Fraction of snow volume conserved during ridging rn_fsnowrft = ', rn_fsnowrft 984 945 ENDIF 985 946 !
Note: See TracChangeset
for help on using the changeset viewer.
