Changeset 14072 for NEMO/trunk/src/ICE/icevar.F90
- Timestamp:
- 2020-12-04T08:48:38+01:00 (4 years ago)
- File:
-
- 1 edited
-
NEMO/trunk/src/ICE/icevar.F90 (modified) (47 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/ICE/icevar.F90
r14005 r14072 34 34 !! - st_i(jpi,jpj) 35 35 !! - et_s(jpi,jpj) total snow heat content 36 !! - et_i(jpi,jpj) total ice thermal content 36 !! - et_i(jpi,jpj) total ice thermal content 37 37 !! - sm_i(jpi,jpj) mean ice salinity 38 38 !! - tm_i(jpi,jpj) mean ice temperature … … 55 55 !!---------------------------------------------------------------------- 56 56 USE dom_oce ! ocean space and time domain 57 USE phycst ! physical constants (ocean directory) 57 USE phycst ! physical constants (ocean directory) 58 58 USE sbc_oce , ONLY : sss_m, ln_ice_embd, nn_fsbc 59 59 USE ice ! sea-ice: variables … … 67 67 PRIVATE 68 68 69 PUBLIC ice_var_agg 70 PUBLIC ice_var_glo2eqv 71 PUBLIC ice_var_eqv2glo 72 PUBLIC ice_var_salprof 73 PUBLIC ice_var_salprof1d 69 PUBLIC ice_var_agg 70 PUBLIC ice_var_glo2eqv 71 PUBLIC ice_var_eqv2glo 72 PUBLIC ice_var_salprof 73 PUBLIC ice_var_salprof1d 74 74 PUBLIC ice_var_zapsmall 75 75 PUBLIC ice_var_zapneg 76 76 PUBLIC ice_var_roundoff 77 PUBLIC ice_var_bv 78 PUBLIC ice_var_enthalpy 77 PUBLIC ice_var_bv 78 PUBLIC ice_var_enthalpy 79 79 PUBLIC ice_var_sshdyn 80 80 PUBLIC ice_var_itd … … 108 108 !! *** ROUTINE ice_var_agg *** 109 109 !! 110 !! ** Purpose : aggregates ice-thickness-category variables to 110 !! ** Purpose : aggregates ice-thickness-category variables to 111 111 !! all-ice variables, i.e. it turns VGLO into VAGG 112 112 !!------------------------------------------------------------------- … … 130 130 vt_il(:,:) = SUM( v_il(:,:,:), dim=3 ) 131 131 ! 132 ato_i(:,:) = 1._wp - at_i(:,:) ! open water fraction 132 ato_i(:,:) = 1._wp - at_i(:,:) ! open water fraction 133 133 ! 134 134 !!GS: tm_su always needed by ABL over sea-ice … … 155 155 hm_i(:,:) = vt_i(:,:) * z1_at_i(:,:) 156 156 hm_s(:,:) = vt_s(:,:) * z1_at_i(:,:) 157 ! 157 ! 158 158 ! ! mean temperature (K), salinity and age 159 159 tm_si(:,:) = SUM( t_si(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:) … … 182 182 WHERE( at_ip(:,:) > epsi20 ) ; hm_ip(:,:) = vt_ip(:,:) / at_ip(:,:) ; hm_il(:,:) = vt_il(:,:) / at_ip(:,:) 183 183 ELSEWHERE ; hm_ip(:,:) = 0._wp ; hm_il(:,:) = 0._wp 184 END WHERE 184 END WHERE 185 185 ! 186 186 DEALLOCATE( z1_vt_i , z1_vt_s ) … … 197 197 !! *** ROUTINE ice_var_glo2eqv *** 198 198 !! 199 !! ** Purpose : computes equivalent variables as function of 199 !! ** Purpose : computes equivalent variables as function of 200 200 !! global variables, i.e. it turns VGLO into VEQV 201 201 !!------------------------------------------------------------------- … … 210 210 !!------------------------------------------------------------------- 211 211 212 !!gm Question 2: It is possible to define existence of sea-ice in a common way between 212 !!gm Question 2: It is possible to define existence of sea-ice in a common way between 213 213 !! ice area and ice volume ? 214 214 !! the idea is to be able to define one for all at the begining of this routine … … 234 234 235 235 zhmax = hi_max(jpl) 236 z1_zhmax = 1._wp / hi_max(jpl) 236 z1_zhmax = 1._wp / hi_max(jpl) 237 237 WHERE( h_i(:,:,jpl) > zhmax ) ! bound h_i by hi_max (i.e. 99 m) with associated update of ice area 238 238 h_i (:,:,jpl) = zhmax 239 a_i (:,:,jpl) = v_i(:,:,jpl) * z1_zhmax 239 a_i (:,:,jpl) = v_i(:,:,jpl) * z1_zhmax 240 240 z1_a_i(:,:,jpl) = zhmax * z1_v_i(:,:,jpl) 241 241 END WHERE 242 242 ! !--- snow thickness 243 243 h_s(:,:,:) = v_s (:,:,:) * z1_a_i(:,:,:) 244 ! !--- ice age 244 ! !--- ice age 245 245 o_i(:,:,:) = oa_i(:,:,:) * z1_a_i(:,:,:) 246 ! !--- pond and lid thickness 246 ! !--- pond and lid thickness 247 247 h_ip(:,:,:) = v_ip(:,:,:) * z1_a_ip(:,:,:) 248 248 h_il(:,:,:) = v_il(:,:,:) * z1_a_ip(:,:,:) … … 258 258 a_ip_eff = MIN( a_ip_eff, 1._wp - za_s_fra ) ! make sure (a_ip_eff + a_s_fra) <= 1 259 259 ! 260 ! !--- salinity (with a minimum value imposed everywhere) 260 ! !--- salinity (with a minimum value imposed everywhere) 261 261 IF( nn_icesal == 2 ) THEN 262 262 WHERE( v_i(:,:,:) > epsi20 ) ; s_i(:,:,:) = MAX( rn_simin , MIN( rn_simax, sv_i(:,:,:) * z1_v_i(:,:,:) ) ) … … 272 272 DO jl = 1, jpl 273 273 DO_3D( 1, 1, 1, 1, 1, nlay_i ) 274 IF ( v_i(ji,jj,jl) > epsi20 ) THEN !--- icy area 274 IF ( v_i(ji,jj,jl) > epsi20 ) THEN !--- icy area 275 275 ! 276 276 ze_i = e_i (ji,jj,jk,jl) * z1_v_i(ji,jj,jl) * zlay_i ! Energy of melting e(S,T) [J.m-3] … … 300 300 END DO 301 301 ! 302 ! integrated values 302 ! integrated values 303 303 vt_i (:,:) = SUM( v_i , dim=3 ) 304 304 vt_s (:,:) = SUM( v_s , dim=3 ) … … 312 312 !! *** ROUTINE ice_var_eqv2glo *** 313 313 !! 314 !! ** Purpose : computes global variables as function of 314 !! ** Purpose : computes global variables as function of 315 315 !! equivalent variables, i.e. it turns VEQV into VGLO 316 316 !!------------------------------------------------------------------- … … 329 329 !! *** ROUTINE ice_var_salprof *** 330 330 !! 331 !! ** Purpose : computes salinity profile in function of bulk salinity 332 !! 333 !! ** Method : If bulk salinity greater than zsi1, 331 !! ** Purpose : computes salinity profile in function of bulk salinity 332 !! 333 !! ** Method : If bulk salinity greater than zsi1, 334 334 !! the profile is assumed to be constant (S_inf) 335 335 !! If bulk salinity lower than zsi0, … … 348 348 !!------------------------------------------------------------------- 349 349 350 !!gm Question: Remove the option 3 ? How many years since it last use ? 350 !!gm Question: Remove the option 3 ? How many years since it last use ? 351 351 352 352 SELECT CASE ( nn_icesal ) … … 369 369 END DO 370 370 END DO 371 ! ! Slope of the linear profile 371 ! ! Slope of the linear profile 372 372 WHERE( h_i(:,:,:) > epsi20 ) ; z_slope_s(:,:,:) = 2._wp * s_i(:,:,:) / h_i(:,:,:) 373 373 ELSEWHERE ; z_slope_s(:,:,:) = 0._wp … … 379 379 zalpha(ji,jj,jl) = MAX( 0._wp , MIN( ( zsi1 - s_i(ji,jj,jl) ) * z1_dS , 1._wp ) ) 380 380 ! ! force a constant profile when SSS too low (Baltic Sea) 381 IF( 2._wp * s_i(ji,jj,jl) >= sss_m(ji,jj) ) zalpha(ji,jj,jl) = 0._wp 381 IF( 2._wp * s_i(ji,jj,jl) >= sss_m(ji,jj) ) zalpha(ji,jj,jl) = 0._wp 382 382 END_2D 383 383 END DO … … 448 448 ALLOCATE( z_slope_s(jpij), zalpha(jpij) ) 449 449 ! 450 ! ! Slope of the linear profile 450 ! ! Slope of the linear profile 451 451 WHERE( h_i_1d(1:npti) > epsi20 ) ; z_slope_s(1:npti) = 2._wp * s_i_1d(1:npti) / h_i_1d(1:npti) 452 452 ELSEWHERE ; z_slope_s(1:npti) = 0._wp 453 453 END WHERE 454 454 455 455 z1_dS = 1._wp / ( zsi1 - zsi0 ) 456 456 DO ji = 1, npti … … 557 557 END_2D 558 558 ! 559 END DO 559 END DO 560 560 561 561 ! to be sure that at_i is the sum of a_i(jl) … … 648 648 END_2D 649 649 ! 650 END DO 650 END DO 651 651 ! 652 652 WHERE( pato_i(:,:) < 0._wp ) pato_i(:,:) = 0._wp … … 693 693 ! 694 694 END SUBROUTINE ice_var_roundoff 695 695 696 696 697 697 SUBROUTINE ice_var_bv … … 713 713 DO jl = 1, jpl 714 714 DO jk = 1, nlay_i 715 WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 ) 715 WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 ) 716 716 bv_i(:,:,jl) = bv_i(:,:,jl) - rTmlt * sz_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 ) 717 717 END WHERE … … 727 727 SUBROUTINE ice_var_enthalpy 728 728 !!------------------------------------------------------------------- 729 !! *** ROUTINE ice_var_enthalpy *** 730 !! 729 !! *** ROUTINE ice_var_enthalpy *** 730 !! 731 731 !! ** Purpose : Computes sea ice energy of melting q_i (J.m-3) from temperature 732 732 !! … … 734 734 !!------------------------------------------------------------------- 735 735 INTEGER :: ji, jk ! dummy loop indices 736 REAL(wp) :: ztmelts ! local scalar 736 REAL(wp) :: ztmelts ! local scalar 737 737 !!------------------------------------------------------------------- 738 738 ! … … 741 741 ztmelts = - rTmlt * sz_i_1d(ji,jk) 742 742 t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts + rt0 ) ! Force t_i_1d to be lower than melting point => likely conservation issue 743 ! (sometimes zdf scheme produces abnormally high temperatures) 743 ! (sometimes zdf scheme produces abnormally high temperatures) 744 744 e_i_1d(ji,jk) = rhoi * ( rcpi * ( ztmelts - ( t_i_1d(ji,jk) - rt0 ) ) & 745 745 & + rLfus * ( 1._wp - ztmelts / ( t_i_1d(ji,jk) - rt0 ) ) & … … 755 755 END SUBROUTINE ice_var_enthalpy 756 756 757 757 758 758 FUNCTION ice_var_sshdyn(pssh, psnwice_mass, psnwice_mass_b) 759 759 !!--------------------------------------------------------------------- 760 760 !! *** ROUTINE ice_var_sshdyn *** 761 !! 761 !! 762 762 !! ** Purpose : compute the equivalent ssh in lead when sea ice is embedded 763 763 !! … … 765 765 !! 766 766 !! ** Reference : Jean-Michel Campin, John Marshall, David Ferreira, 767 !! Sea ice-ocean coupling using a rescaled vertical coordinate z*, 767 !! Sea ice-ocean coupling using a rescaled vertical coordinate z*, 768 768 !! Ocean Modelling, Volume 24, Issues 1-2, 2008 769 769 !!---------------------------------------------------------------------- … … 783 783 ! compute ice load used to define the equivalent ssh in lead 784 784 IF( ln_ice_embd ) THEN 785 ! 785 ! 786 786 ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1} 787 787 ! = (1/nn_fsbc)^2 * {SUM[n] , n=0,nn_fsbc-1} … … 802 802 END FUNCTION ice_var_sshdyn 803 803 804 804 805 805 !!------------------------------------------------------------------- 806 806 !! *** INTERFACE ice_var_itd *** … … 831 831 ph_ip(:) = phtip(:) 832 832 ph_il(:) = phtil(:) 833 833 834 834 END SUBROUTINE ice_var_itd_1c1c 835 835 … … 846 846 REAL(wp), ALLOCATABLE, DIMENSION(:) :: z1_ai, z1_vi, z1_vs 847 847 ! 848 INTEGER :: idim 848 INTEGER :: idim 849 849 !!------------------------------------------------------------------- 850 850 ! … … 888 888 ! 889 889 END SUBROUTINE ice_var_itd_Nc1c 890 890 891 891 SUBROUTINE ice_var_itd_1cMc( phti, phts, pati , ph_i, ph_s, pa_i, & 892 892 & ptmi, ptms, ptmsu, psmi, patip, phtip, phtil, pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il ) … … 898 898 !! ** Method: ice thickness distribution follows a gamma function from Abraham et al. (2015) 899 899 !! it has the property of conserving total concentration and volume 900 !! 900 !! 901 901 !! 902 902 !! ** Arguments : phti: 1-cat ice thickness … … 904 904 !! pati: 1-cat ice concentration 905 905 !! 906 !! ** Output : jpl-cat 906 !! ** Output : jpl-cat 907 907 !! 908 908 !! Abraham, C., Steiner, N., Monahan, A. and Michel, C., 2015. 909 909 !! Effects of subgrid‐scale snow thickness variability on radiative transfer in sea ice. 910 !! Journal of Geophysical Research: Oceans, 120(8), pp.5597-5614 910 !! Journal of Geophysical Research: Oceans, 120(8), pp.5597-5614 911 911 !!------------------------------------------------------------------- 912 912 REAL(wp), DIMENSION(:), INTENT(in) :: phti, phts, pati ! input ice/snow variables … … 987 987 ! In case snow load is in excess that would lead to transformation from snow to ice 988 988 ! Then, transfer the snow excess into the ice (different from icethd_dh) 989 zdh = MAX( 0._wp, ( rhos * ph_s(ji,jl) + ( rhoi - rho0 ) * ph_i(ji,jl) ) * r1_rho0 ) 989 zdh = MAX( 0._wp, ( rhos * ph_s(ji,jl) + ( rhoi - rho0 ) * ph_i(ji,jl) ) * r1_rho0 ) 990 990 ! recompute h_i, h_s avoiding out of bounds values 991 991 ph_i(ji,jl) = MIN( hi_max(jl), ph_i(ji,jl) + zdh ) … … 1047 1047 !! 1048 1048 !! ** Method: Iterative procedure 1049 !! 1049 !! 1050 1050 !! 1) Fill ice cat that correspond to input thicknesses 1051 1051 !! Find the lowest(jlmin) and highest(jlmax) cat that are filled 1052 1052 !! 1053 1053 !! 2) Expand the filling to the cat jlmin-1 and jlmax+1 1054 !! by removing 25% ice area from jlmin and jlmax (resp.) 1055 !! 1056 !! 3) Expand the filling to the empty cat between jlmin and jlmax 1054 !! by removing 25% ice area from jlmin and jlmax (resp.) 1055 !! 1056 !! 3) Expand the filling to the empty cat between jlmin and jlmax 1057 1057 !! by a) removing 25% ice area from the lower cat (ascendant loop jlmin=>jlmax) 1058 1058 !! b) removing 25% ice area from the higher cat (descendant loop jlmax=>jlmin) … … 1062 1062 !! pati: N-cat ice concentration 1063 1063 !! 1064 !! ** Output : jpl-cat 1065 !! 1066 !! (Example of application: BDY forcings when inputs have N-cat /= jpl) 1064 !! ** Output : jpl-cat 1065 !! 1066 !! (Example of application: BDY forcings when inputs have N-cat /= jpl) 1067 1067 !!------------------------------------------------------------------- 1068 1068 REAL(wp), DIMENSION(:,:), INTENT(in) :: phti, phts, pati ! input ice/snow variables … … 1077 1077 REAL(wp), PARAMETER :: ztrans = 0.25_wp 1078 1078 INTEGER :: ji, jl, jl1, jl2 1079 INTEGER :: idim, icat 1079 INTEGER :: idim, icat 1080 1080 !!------------------------------------------------------------------- 1081 1081 ! … … 1116 1116 ELSE ! input cat /= output cat ! 1117 1117 ! ! ----------------------- ! 1118 1118 1119 1119 ALLOCATE( jlfil(idim,jpl), jlfil2(idim,jpl) ) ! allocate arrays 1120 1120 ALLOCATE( jlmin(idim), jlmax(idim) ) … … 1126 1126 ! 1127 1127 ! --- fill the categories --- ! 1128 ! find where cat-input = cat-output and fill cat-output fields 1128 ! find where cat-input = cat-output and fill cat-output fields 1129 1129 jlmax(:) = 0 1130 1130 jlmin(:) = 999 … … 1147 1147 END DO 1148 1148 ! 1149 ! --- fill the gaps between categories --- ! 1149 ! --- fill the gaps between categories --- ! 1150 1150 ! transfer from categories filled at the previous step to the empty ones in between 1151 1151 DO ji = 1, idim … … 1168 1168 END DO 1169 1169 ! 1170 jlfil2(:,:) = jlfil(:,:) 1170 jlfil2(:,:) = jlfil(:,:) 1171 1171 ! fill categories from low to high 1172 1172 DO jl = 2, jpl-1 … … 1189 1189 ! fill low 1190 1190 pa_i(ji,jl) = pa_i(ji,jl) + ztrans * pa_i(ji,jl+1) 1191 ph_i(ji,jl) = hi_mean(jl) 1191 ph_i(ji,jl) = hi_mean(jl) 1192 1192 jlfil2(ji,jl) = jl 1193 1193 ! remove high … … 1279 1279 !! we argue that snow does not cover the whole ice because 1280 1280 !! of wind blowing... 1281 !! 1281 !! 1282 1282 !! ** Arguments : ph_s: snow thickness 1283 !! 1283 !! 1284 1284 !! ** Output : pa_s_fra: fraction of ice covered by snow 1285 1285 !! … … 1326 1326 ENDIF 1327 1327 END SUBROUTINE ice_var_snwfra_1d 1328 1328 1329 1329 !!-------------------------------------------------------------------------- 1330 1330 !! INTERFACE ice_var_snwblow … … 1336 1336 !! If snow fall was uniform, a fraction (1-at_i) would fall into leads 1337 1337 !! but because of the winds, more snow falls on leads than on sea ice 1338 !! and a greater fraction (1-at_i)^beta of the total mass of snow 1338 !! and a greater fraction (1-at_i)^beta of the total mass of snow 1339 1339 !! (beta < 1) falls in leads. 1340 !! In reality, beta depends on wind speed, 1341 !! and should decrease with increasing wind speed but here, it is 1340 !! In reality, beta depends on wind speed, 1341 !! and should decrease with increasing wind speed but here, it is 1342 1342 !! considered as a constant. an average value is 0.66 1343 1343 !!--------------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.
