Changeset 11560
- Timestamp:
- 2019-09-18T13:31:36+02:00 (5 years ago)
- Location:
- NEMO/trunk/src/ICE
- Files:
-
- 2 edited
-
icedyn_rdgrft.F90 (modified) (1 diff)
-
icevar.F90 (modified) (5 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/ICE/icedyn_rdgrft.F90
r11536 r11560 278 278 ! controls 279 279 IF( ln_ctl ) CALL ice_prt3D ('icedyn_rdgrft') ! prints 280 IF( ln_icectl ) CALL ice_prt (kt, iiceprt, jiceprt,-1, ' - ice dyn rdgrft - ') ! prints 280 281 IF( ln_icediachk ) CALL ice_cons_hsm(1, 'icedyn_rdgrft', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 281 282 IF( ln_icediachk ) CALL ice_cons2D (1, 'icedyn_rdgrft', diag_v, diag_s, diag_t, diag_fv, diag_fs, diag_ft) ! conservation -
NEMO/trunk/src/ICE/icevar.F90
r11536 r11560 868 868 !! ** Purpose : converting 1-cat ice to jpl ice categories 869 869 !! 870 !! ice thickness distribution follows a gaussian law 871 !! around the concentration of the most likely ice thickness 872 !! (similar as iceistate.F90) 873 !! 874 !! ** Method: Iterative procedure 875 !! 876 !! 1) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian 877 !! 878 !! 2) Check whether the distribution conserves area and volume, positivity and 879 !! category boundaries 870 !! 871 !! ** Method: ice thickness distribution follows a gamma function from Abraham et al. (2015) 872 !! it has the property of conserving total concentration and volume 880 873 !! 881 !! 3) If not (input ice is too thin), the last category is empty and882 !! the number of categories is reduced (jpl-1)883 !!884 !! 4) Iterate until ok (SUM(itest(:) = 4)885 874 !! 886 875 !! ** Arguments : phti: 1-cat ice thickness … … 890 879 !! ** Output : jpl-cat 891 880 !! 892 !! (Example of application: BDY forcings when input are cell averaged) 881 !! Abraham, C., Steiner, N., Monahan, A. and Michel, C., 2015. 882 !! Effects of subgrid‐scale snow thickness variability on radiative transfer in sea ice. 883 !! Journal of Geophysical Research: Oceans, 120(8), pp.5597-5614 893 884 !!------------------------------------------------------------------- 894 885 REAL(wp), DIMENSION(:), INTENT(in) :: phti, phts, pati ! input ice/snow variables … … 897 888 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ! output ice/snow temp & sal & ponds 898 889 ! 899 INTEGER , DIMENSION(4) :: itest 900 REAL(wp), ALLOCATABLE, DIMENSION(:) :: zfra 890 REAL(wp), ALLOCATABLE, DIMENSION(:) :: zfra, z1_hti 901 891 INTEGER :: ji, jk, jl 902 INTEGER :: idim, i_fill, jl0 903 REAL(wp) :: zarg, zV, zconv, zdh, zdv 904 !!------------------------------------------------------------------- 905 ! 906 ! == thickness and concentration == ! 907 ! distribution over the jpl ice categories 908 ! a gaussian distribution for ice concentration is used 909 ! then we check whether the distribution fullfills 910 ! volume and area conservation, positivity and ice categories bounds 892 INTEGER :: idim 893 REAL(wp) :: zv, zdh 894 !!------------------------------------------------------------------- 895 ! 911 896 idim = SIZE( phti , 1 ) 912 897 ! … … 915 900 pa_i(1:idim,1:jpl) = 0._wp 916 901 ! 902 ALLOCATE( z1_hti(idim) ) 903 WHERE( phti(:) /= 0._wp ) ; z1_hti(:) = 1._wp / phti(:) 904 ELSEWHERE ; z1_hti(:) = 0._wp 905 END WHERE 906 ! 907 ! == thickness and concentration == ! 908 ! for categories 1:jpl-1, integrate the gamma function from hi_max(jl-1) to hi_max(jl) 909 DO jl = 1, jpl-1 910 DO ji = 1, idim 911 ! 912 IF( phti(ji) > 0._wp ) THEN 913 ! concentration : integrate ((4A/H^2)xexp(-2x/H))dx from x=hi_max(jl-1) to hi_max(jl) 914 pa_i(ji,jl) = pati(ji) * z1_hti(ji) * ( ( phti(ji) + 2.*hi_max(jl-1) ) * EXP( -2.*hi_max(jl-1)*z1_hti(ji) ) & 915 & - ( phti(ji) + 2.*hi_max(jl ) ) * EXP( -2.*hi_max(jl )*z1_hti(ji) ) ) 916 ! 917 ! volume : integrate ((4A/H^2)x^2exp(-2x/H))dx from x=hi_max(jl-1) to hi_max(jl) 918 zv = pati(ji) * z1_hti(ji) * ( ( phti(ji)*phti(ji) + 2.*phti(ji)*hi_max(jl-1) + 2.*hi_max(jl-1)*hi_max(jl-1) ) & 919 & * EXP( -2.*hi_max(jl-1)*z1_hti(ji) ) & 920 & - ( phti(ji)*phti(ji) + 2.*phti(ji)*hi_max(jl) + 2.*hi_max(jl)*hi_max(jl) ) & 921 & * EXP(-2.*hi_max(jl)*z1_hti(ji)) ) 922 ! thickness 923 IF( pa_i(ji,jl) > epsi06 ) THEN 924 ph_i(ji,jl) = zv / pa_i(ji,jl) 925 ELSE 926 ph_i(ji,jl) = 0. 927 pa_i(ji,jl) = 0. 928 ENDIF 929 ENDIF 930 ! 931 ENDDO 932 ENDDO 933 ! 934 ! for the last category (jpl), integrate the gamma function from hi_max(jpl-1) to infinity 917 935 DO ji = 1, idim 918 936 ! 919 937 IF( phti(ji) > 0._wp ) THEN 920 ! 921 ! find which category (jl0) the input ice thickness falls into 922 jl0 = jpl 923 DO jl = 1, jpl 924 IF ( ( phti(ji) >= hi_max(jl-1) ) .AND. ( phti(ji) < hi_max(jl) ) ) THEN 925 jl0 = jl 926 CYCLE 927 ENDIF 928 END DO 929 ! 930 itest(:) = 0 931 i_fill = jpl + 1 !------------------------------------ 932 DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) ) ! iterative loop on i_fill categories 933 ! !------------------------------------ 934 i_fill = i_fill - 1 935 ! 936 ph_i(ji,1:jpl) = 0._wp 937 pa_i(ji,1:jpl) = 0._wp 938 itest(:) = 0 939 ! 940 IF ( i_fill == 1 ) THEN !-- case very thin ice: fill only category 1 941 ph_i(ji,1) = phti(ji) 942 pa_i(ji,1) = pati (ji) 943 ELSE !-- case ice is thicker: fill categories >1 944 ! thickness 945 DO jl = 1, i_fill - 1 946 ph_i(ji,jl) = hi_mean(jl) 947 END DO 948 ! 949 ! concentration 950 pa_i(ji,jl0) = pati(ji) / SQRT(REAL(jpl)) 951 DO jl = 1, i_fill - 1 952 IF ( jl /= jl0 ) THEN 953 zarg = ( ph_i(ji,jl) - phti(ji) ) / ( phti(ji) * 0.5_wp ) 954 pa_i(ji,jl) = pa_i (ji,jl0) * EXP(-zarg**2) 955 ENDIF 956 END DO 957 ! 958 ! last category 959 pa_i(ji,i_fill) = pati(ji) - SUM( pa_i(ji,1:i_fill-1) ) 960 zV = SUM( pa_i(ji,1:i_fill-1) * ph_i(ji,1:i_fill-1) ) 961 ph_i(ji,i_fill) = ( phti(ji) * pati(ji) - zV ) / MAX( pa_i(ji,i_fill), epsi10 ) 962 ! 963 ! correction if concentration of upper cat is greater than lower cat 964 ! (it should be a gaussian around jl0 but sometimes it is not) 965 IF ( jl0 /= jpl ) THEN 966 DO jl = jpl, jl0+1, -1 967 IF ( pa_i(ji,jl) > pa_i(ji,jl-1) ) THEN 968 zdv = ph_i(ji,jl) * pa_i(ji,jl) 969 ph_i(ji,jl ) = 0._wp 970 pa_i (ji,jl ) = 0._wp 971 pa_i (ji,1:jl-1) = pa_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * phti(ji), epsi10 ) 972 END IF 973 END DO 974 ENDIF 975 ! 976 ENDIF 977 ! 978 ! Compatibility tests 979 zconv = ABS( pati(ji) - SUM( pa_i(ji,1:jpl) ) ) 980 IF ( zconv < epsi06 ) itest(1) = 1 ! Test 1: area conservation 981 ! 982 zconv = ABS( phti(ji)*pati(ji) - SUM( pa_i(ji,1:jpl)*ph_i(ji,1:jpl) ) ) 983 IF ( zconv < epsi06 ) itest(2) = 1 ! Test 2: volume conservation 984 ! 985 IF ( ph_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1 ! Test 3: thickness of the last category is in-bounds ? 986 ! 987 itest(4) = 1 988 DO jl = 1, i_fill 989 IF ( pa_i(ji,jl) < 0._wp ) itest(4) = 0 ! Test 4: positivity of ice concentrations 990 END DO 991 ! !---------------------------- 992 END DO ! end iteration on categories 993 ! !---------------------------- 938 ! concentration : integrate ((4A/H^2)xexp(-2x/H))dx from x=hi_max(jpl-1) to infinity 939 pa_i(ji,jpl) = pati(ji) * z1_hti(ji) * ( phti(ji) + 2.*hi_max(jpl-1) ) * EXP( -2.*hi_max(jpl-1)*z1_hti(ji) ) 940 941 ! volume : integrate ((4A/H^2)x^2exp(-2x/H))dx from x=hi_max(jpl-1) to infinity 942 zv = pati(ji) * z1_hti(ji) * ( phti(ji)*phti(ji) + 2.*phti(ji)*hi_max(jpl-1) + 2.*hi_max(jpl-1)*hi_max(jpl-1) ) & 943 & * EXP( -2.*hi_max(jpl-1)*z1_hti(ji) ) 944 ! thickness 945 IF( pa_i(ji,jpl) > epsi06 ) THEN 946 ph_i(ji,jpl) = zv / pa_i(ji,jpl) 947 else 948 ph_i(ji,jpl) = 0. 949 pa_i(ji,jpl) = 0. 950 ENDIF 994 951 ENDIF 995 END DO 996 952 ! 953 ENDDO 954 ! 997 955 ! Add Snow in each category where pa_i is not 0 998 956 DO jl = 1, jpl 999 957 DO ji = 1, idim 1000 958 IF( pa_i(ji,jl) > 0._wp ) THEN 1001 ph_s(ji,jl) = ph_i(ji,jl) * ( phts(ji) / phti(ji))959 ph_s(ji,jl) = ph_i(ji,jl) * phts(ji) * z1_hti(ji) 1002 960 ! In case snow load is in excess that would lead to transformation from snow to ice 1003 961 ! Then, transfer the snow excess into the ice (different from icethd_dh) … … 1009 967 END DO 1010 968 END DO 969 ! 970 DEALLOCATE( z1_hti ) 1011 971 ! 1012 972 ! == temperature and salinity == !
Note: See TracChangeset
for help on using the changeset viewer.
