Changeset 10910
- Timestamp:
- 2019-04-29T18:10:38+02:00 (5 years ago)
- Location:
- NEMO/releases/release-4.0
- Files:
-
- 11 edited
-
cfgs/ORCA2_ICE_PISCES/EXPREF/file_def_nemo-ice.xml (modified) (1 diff)
-
cfgs/SHARED/field_def_nemo-ice.xml (modified) (1 diff)
-
cfgs/SHARED/namelist_ice_ref (modified) (1 diff)
-
cfgs/SPITZ12/EXPREF/file_def_nemo-ice.xml (modified) (1 diff)
-
cfgs/SPITZ12/EXPREF/namelist_cfg (modified) (10 diffs)
-
cfgs/SPITZ12/EXPREF/namelist_ice_cfg (modified) (4 diffs)
-
src/ICE/icedyn.F90 (modified) (6 diffs)
-
src/ICE/icedyn_adv.F90 (modified) (2 diffs)
-
src/ICE/icedyn_adv_umx.F90 (modified) (35 diffs)
-
src/ICE/icedyn_rhg.F90 (modified) (2 diffs)
-
src/ICE/icewri.F90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/releases/release-4.0/cfgs/ORCA2_ICE_PISCES/EXPREF/file_def_nemo-ice.xml
r9943 r10910 79 79 <field field_ref="vfxsnw" name="vfxsnw" /> 80 80 81 <!-- diag error for negative ice volume after advection -->82 <field field_ref="iceneg_pres" name="sineg_pres" />83 <field field_ref="iceneg_volu" name="sineg_volu" />84 <field field_ref="iceneg_hfx" name="sineg_hfx" />85 86 81 <!-- categories --> 87 82 <field field_ref="icemask_cat" name="simskcat"/> -
NEMO/releases/release-4.0/cfgs/SHARED/field_def_nemo-ice.xml
r10578 r10910 160 160 <field id="hfxdhc" long_name="Heat content variation in snow and ice (neg = ice cooling)" unit="W/m2" /> 161 161 162 <!-- diagnostics of the negative values resulting from the advection scheme -->163 <field id="iceneg_pres" long_name="Fraction of time steps with negative sea ice volume" unit="" />164 <field id="iceneg_volu" long_name="Negative sea ice volume per area arising from advection" unit="m" />165 <field id="iceneg_hfx" long_name="Negative sea ice heat content (eq. heat flux) arising from advection" unit="W/m2" />166 167 162 <!-- sbcssm variables --> 168 163 <field id="sst_m" unit="degC" /> -
NEMO/releases/release-4.0/cfgs/SHARED/namelist_ice_ref
r10612 r10910 154 154 &namthd_do ! Ice growth in open water 155 155 !------------------------------------------------------------------------------ 156 rn_hinew = 0.1 ! thickness for new ice formation in open water (m), must be larger than rn_h newice156 rn_hinew = 0.1 ! thickness for new ice formation in open water (m), must be larger than rn_himin 157 157 ln_frazil = .false. ! Frazil ice parameterization (ice collection as a function of wind) 158 158 rn_maxfraz = 1.0 ! maximum fraction of frazil ice collecting at the ice base -
NEMO/releases/release-4.0/cfgs/SPITZ12/EXPREF/file_def_nemo-ice.xml
r9943 r10910 79 79 <field field_ref="vfxsnw" name="vfxsnw" /> 80 80 81 <!-- diag error for negative ice volume after advection -->82 <field field_ref="iceneg_pres" name="sineg_pres" />83 <field field_ref="iceneg_volu" name="sineg_volu" />84 <field field_ref="iceneg_hfx" name="sineg_hfx" />85 86 81 <!-- categories --> 87 82 <field field_ref="icemask_cat" name="simskcat"/> -
NEMO/releases/release-4.0/cfgs/SPITZ12/EXPREF/namelist_cfg
r10075 r10910 80 80 ! Sea-ice : 81 81 nn_ice = 2 ! SI3 82 ln_ice_embd = .false. ! =T embedded sea-ice (pressure + mass and salt exchanges) 83 ! ! =F levitating ice (no pressure, mass and salt exchanges) 82 84 ! Misc. options of sbc : 83 85 ln_traqsr = .true. ! Light penetration in the ocean (T => fill namtra_qsr ) … … 90 92 ! 91 93 ln_Cd_L12 = .false. ! air-ice drags = F(ice concentration) (Lupkes et al. 2012) 92 ln_Cd_L15 = . false.! air-ice drags = F(ice concentration) (Lupkes et al. 2015)94 ln_Cd_L15 = .true. ! air-ice drags = F(ice concentration) (Lupkes et al. 2015) 93 95 ! 94 96 cn_dir = './' ! root directory for the bulk data location … … 96 98 ! ! file name ! frequency (hours) ! variable ! time interp.! clim ! 'yearly'/ ! weights filename ! rotation ! land/sea mask ! 97 99 ! ! ! (if <0 months) ! name ! (logical) ! (T/F) ! 'monthly' ! ! pairing ! filename ! 98 !! ERAI99 ! sn_wndi = 'u10_era_spitz' , 6 , 'u10' , .true. , .false. , 'yearly' , 'weights_bicub', 'Uwnd' , ''100 ! sn_wndj = 'v10_era_spitz' , 6 , 'v10' , .true. , .false. , 'yearly' , 'weights_bicub', 'Vwnd' , ''101 ! sn_qsr = 'ssrd_era_spitz' , 6 , 'ssrd' , .true. , .false. , 'yearly' , 'weights_bilin', '' , ''102 ! sn_qlw = 'strd_era_spitz' , 6 , 'strd' , .true. , .false. , 'yearly' , 'weights_bilin', '' , ''103 ! sn_tair = 't2_era_spitz' , 6 , 't2' , .true. , .false. , 'yearly' , 'weights_bilin', '' , ''104 ! sn_humi = 'humi_era_spitz' , 6 , 'humi' , .true. , .false. , 'yearly' , 'weights_bilin', '' , ''105 ! sn_prec = 'precip_era_spitz' , 6 , 'precip' , .true. , .false. , 'yearly' , 'weights_bilin', '' , ''106 ! sn_snow = 'snow_era_spitz' , 6 , 'snow' , .true. , .false. , 'yearly' , 'weights_bilin', '' , ''107 ! sn_tdif = 'taudif' , 6 , 'taudif' , .true. , .false. , 'yearly' , '' , '' , ''108 ! rn_zqt = 10. ! Air temperature & humidity reference height (m)109 !! MAR110 100 sn_wndi = 'MARv3.6-9km-Svalbard-2hourly_spitz' , 2 , 'u10' , .true. , .false. , 'yearly' , 'weights_bicub', 'Uwnd' , '' 111 101 sn_wndj = 'MARv3.6-9km-Svalbard-2hourly_spitz' , 2 , 'v10' , .true. , .false. , 'yearly' , 'weights_bicub', 'Vwnd' , '' … … 129 119 ! ! type of penetration (default: NO selection) 130 120 ln_qsr_rgb = .true. ! RGB light penetration (Red-Green-Blue) 131 /132 !-----------------------------------------------------------------------133 &namsbc_ssr ! surface boundary condition : sea surface restoring (ln_ssr =T)134 !-----------------------------------------------------------------------135 /136 !-----------------------------------------------------------------------137 &namsbc_rnf ! runoffs (ln_rnf =T)138 !-----------------------------------------------------------------------139 121 / 140 122 !!====================================================================== … … 236 218 !----------------------------------------------------------------------- 237 219 ln_loglayer = .true. ! logarithmic drag: Cd = vkarmn/log(z/z0) |U| 220 ln_drgimp = .true. ! implicit top/bottom friction flag 238 221 / 239 222 !----------------------------------------------------------------------- 240 223 &namdrg_bot ! BOTTOM friction (ln_OFF =F) 241 224 !----------------------------------------------------------------------- 242 rn_Cd0 = 2.5e-3 ! !1.e-3 !drag coefficient [-]243 rn_Cdmax = 0. 01 !!0.1 ! drag value maximum [-] (logarithmic drag)244 rn_ke0 = 0. ! !2.5e-3 !background kinetic energy [m2/s2] (non-linear cases)225 rn_Cd0 = 2.5e-3 ! drag coefficient [-] 226 rn_Cdmax = 0.1 ! drag value maximum [-] (logarithmic drag) 227 rn_ke0 = 0. ! background kinetic energy [m2/s2] (non-linear cases) 245 228 rn_z0 = 3.e-3 ! roughness [m] (ln_loglayer=T) 246 ln_boost = .false. ! =T regional boost of Cd0 ; =F constant247 rn_boost = 50. ! local boost factor [-]248 /249 !-----------------------------------------------------------------------250 &nambbc ! bottom temperature boundary condition (default: OFF)251 !-----------------------------------------------------------------------252 ln_trabbc = .false. ! Apply a geothermal heating at the ocean bottom253 229 / 254 230 !----------------------------------------------------------------------- … … 258 234 nn_bbl_ldf = 1 ! diffusive bbl (=1) or not (=0) 259 235 nn_bbl_adv = 0 ! advective bbl (=1/2) or not (=0) 260 rn_ahtbbl = 1000. ! lateral mixing coefficient in the bbl [m2/s]261 rn_gambbl = 10. ! advective bbl coefficient [s]262 236 / 263 237 !!====================================================================== … … 293 267 ! ! Coefficients: 294 268 nn_aht_ijk_t = 31 ! space/time variation of eddy coefficient: 295 !! = 31 F(i,j,k,t)=F(local velocity and grid-spacing)269 ! ! = 31 F(i,j,k,t)=F(local velocity and grid-spacing) 296 270 / 297 271 !!====================================================================== … … 320 294 &namdyn_vor ! Vorticity / Coriolis scheme (default: NO selection) 321 295 !----------------------------------------------------------------------- 322 ln_dynvor_e nT = .true. ! energy conserving scheme (T-point)296 ln_dynvor_eeT = .true. ! energy conserving scheme (een using e3t) 323 297 / 324 298 !----------------------------------------------------------------------- … … 352 326 &namzdf ! vertical physics manager (default: NO selection) 353 327 !----------------------------------------------------------------------- 328 ! ! adaptive-implicit vertical advection 329 ln_zad_Aimp = .true. ! Courant number dependent scheme (Shchepetkin 2015) 354 330 ! ! type of vertical closure (required) 355 331 ln_zdftke = .true. ! Turbulent Kinetic Energy closure (T => fill namzdf_tke) 356 332 ln_zdfgls = .false. ! Generic Length Scale closure (T => fill namzdf_gls) 333 ! ! convection 334 ln_zdfevd = .true. ! enhanced vertical diffusion 335 ! 357 336 ln_zdfddm = .true. ! double diffusive mixing 337 ! 358 338 ! ! Coefficients 359 339 rn_avm0 = 1.2e-4 ! vertical eddy viscosity [m2/s] (background Kz if ln_zdfcst=F) … … 384 364 / 385 365 !----------------------------------------------------------------------- 386 &nam_diaharm ! Harmonic analysis of tidal constituents ('key_diaharm')387 !-----------------------------------------------------------------------388 /389 !-----------------------------------------------------------------------390 366 &namnc4 ! netcdf4 chunking and compression settings ("key_netcdf4") 391 367 !----------------------------------------------------------------------- -
NEMO/releases/release-4.0/cfgs/SPITZ12/EXPREF/namelist_ice_cfg
r10535 r10910 35 35 &namdyn ! Ice dynamics 36 36 !------------------------------------------------------------------------------ 37 ln_landfast_L16 = .true. ! landfast: parameterization from Lemieux 2016 37 38 / 38 39 !------------------------------------------------------------------------------ … … 43 44 &namdyn_rhg ! Ice rheology 44 45 !------------------------------------------------------------------------------ 46 ln_rhg_EVP = .true. ! EVP rheology 47 ln_aEVP = .true. ! adaptive rheology (Kimmritz et al. 2016 & 2017) 45 48 / 46 49 !------------------------------------------------------------------------------ … … 67 70 &namthd_do ! Ice growth in open water 68 71 !------------------------------------------------------------------------------ 69 rn_hinew = 0.02 ! thickness for new ice formation in open water (m), must be larger than rn_h newice72 rn_hinew = 0.02 ! thickness for new ice formation in open water (m), must be larger than rn_himin 70 73 ln_frazil = .true. ! Frazil ice parameterization (ice collection as a function of wind) 71 74 / … … 77 80 &namthd_pnd ! Melt ponds 78 81 !------------------------------------------------------------------------------ 79 ln_pnd_H12 = .false. ! activate evolutive melt ponds (from Holland et al 2012) 80 ln_pnd_CST = .false. ! activate constant melt ponds 81 ln_pnd_alb = .false. ! melt ponds affect albedo or not 82 ln_pnd_H12 = .true. ! activate evolutive melt ponds (from Holland et al 2012) 83 ln_pnd_alb = .true. ! melt ponds affect albedo or not 82 84 / 83 85 -
NEMO/releases/release-4.0/src/ICE/icedyn.F90
r10564 r10910 74 74 INTEGER, INTENT(in) :: kt ! ice time step 75 75 !! 76 INTEGER :: ji, jj , jl! dummy loop indices76 INTEGER :: ji, jj ! dummy loop indices 77 77 REAL(wp) :: zcoefu, zcoefv 78 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zhi_max, zhs_max, zhip_max 79 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdivu_i 78 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdivu_i 80 79 !!-------------------------------------------------------------------- 81 80 ! … … 89 88 ENDIF 90 89 ! 91 IF( ln_landfast_home ) THEN !-- Landfast ice parameterization 92 tau_icebfr(:,:) = 0._wp 93 DO jl = 1, jpl 94 WHERE( h_i_b(:,:,jl) > ht_n(:,:) * rn_depfra ) tau_icebfr(:,:) = tau_icebfr(:,:) + a_i(:,:,jl) * rn_icebfr 95 END DO 96 ENDIF 97 ! 98 ! !-- Record max of the surrounding 9-pts ice thick. (for CALL Hbig) 99 DO jl = 1, jpl 100 DO jj = 2, jpjm1 101 DO ji = fs_2, fs_jpim1 102 zhip_max(ji,jj,jl) = MAX( epsi20, h_ip_b(ji,jj,jl), h_ip_b(ji+1,jj ,jl), h_ip_b(ji ,jj+1,jl), & 103 & h_ip_b(ji-1,jj ,jl), h_ip_b(ji ,jj-1,jl), & 104 & h_ip_b(ji+1,jj+1,jl), h_ip_b(ji-1,jj-1,jl), & 105 & h_ip_b(ji+1,jj-1,jl), h_ip_b(ji-1,jj+1,jl) ) 106 zhi_max (ji,jj,jl) = MAX( epsi20, h_i_b (ji,jj,jl), h_i_b (ji+1,jj ,jl), h_i_b (ji ,jj+1,jl), & 107 & h_i_b (ji-1,jj ,jl), h_i_b (ji ,jj-1,jl), & 108 & h_i_b (ji+1,jj+1,jl), h_i_b (ji-1,jj-1,jl), & 109 & h_i_b (ji+1,jj-1,jl), h_i_b (ji-1,jj+1,jl) ) 110 zhs_max (ji,jj,jl) = MAX( epsi20, h_s_b (ji,jj,jl), h_s_b (ji+1,jj ,jl), h_s_b (ji ,jj+1,jl), & 111 & h_s_b (ji-1,jj ,jl), h_s_b (ji ,jj-1,jl), & 112 & h_s_b (ji+1,jj+1,jl), h_s_b (ji-1,jj-1,jl), & 113 & h_s_b (ji+1,jj-1,jl), h_s_b (ji-1,jj+1,jl) ) 114 END DO 115 END DO 116 END DO 117 CALL lbc_lnk_multi( 'icedyn', zhi_max(:,:,:), 'T', 1., zhs_max(:,:,:), 'T', 1., zhip_max(:,:,:), 'T', 1. ) 118 ! 119 ! 120 SELECT CASE( nice_dyn ) !-- Set which dynamics is running 90 ! retrieve thickness from volume for landfast param. and UMx advection scheme 91 WHERE( a_i(:,:,:) >= epsi20 ) 92 h_i(:,:,:) = v_i(:,:,:) / a_i_b(:,:,:) 93 h_s(:,:,:) = v_s(:,:,:) / a_i_b(:,:,:) 94 ELSEWHERE 95 h_i(:,:,:) = 0._wp 96 h_s(:,:,:) = 0._wp 97 END WHERE 98 ! 99 WHERE( a_ip(:,:,:) >= epsi20 ) 100 h_ip(:,:,:) = v_ip(:,:,:) / a_ip(:,:,:) 101 ELSEWHERE 102 h_ip(:,:,:) = 0._wp 103 END WHERE 104 ! 105 ! 106 SELECT CASE( nice_dyn ) !-- Set which dynamics is running 121 107 122 108 CASE ( np_dynALL ) !== all dynamical processes ==! 123 CALL ice_dyn_rhg ( kt ) ! -- rheology 124 CALL ice_dyn_adv ( kt ) ; CALL Hbig( zhi_max, zhs_max, zhip_max ) ! -- advection of ice + correction on ice thickness 125 CALL ice_dyn_rdgrft( kt ) ! -- ridging/rafting 126 CALL ice_cor ( kt , 1 ) ! -- Corrections 127 109 ! 110 CALL ice_dyn_rhg ( kt ) ! -- rheology 111 CALL ice_dyn_adv ( kt ) ! -- advection of ice 112 CALL ice_dyn_rdgrft( kt ) ! -- ridging/rafting 113 CALL ice_cor ( kt , 1 ) ! -- Corrections 114 ! 128 115 CASE ( np_dynRHGADV ) !== no ridge/raft & no corrections ==! 129 CALL ice_dyn_rhg ( kt ) ! -- rheology 130 CALL ice_dyn_adv ( kt ) ; CALL Hbig( zhi_max, zhs_max, zhip_max ) ! -- advection of ice + correction on ice thickness 131 CALL Hpiling ! -- simple pile-up (replaces ridging/rafting) 132 116 ! 117 CALL ice_dyn_rhg ( kt ) ! -- rheology 118 CALL ice_dyn_adv ( kt ) ! -- advection of ice 119 CALL Hpiling ! -- simple pile-up (replaces ridging/rafting) 120 ! 133 121 CASE ( np_dynADV1D ) !== pure advection ==! (1D) 134 ALLOCATE( zdivu_i(jpi,jpj) )122 ! 135 123 ! --- monotonicity test from Schar & Smolarkiewicz 1996 --- ! 136 124 ! CFL = 0.5 at a distance from the bound of 1/6 of the basin length … … 145 133 END DO 146 134 ! --- 147 CALL ice_dyn_adv ( kt ) ! -- advection of ice 148 149 ! diagnostics: divergence at T points 150 DO jj = 2, jpjm1 151 DO ji = 2, jpim1 152 zdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) & 153 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj) 154 END DO 155 END DO 156 CALL lbc_lnk( 'icedyn', zdivu_i, 'T', 1. ) 157 IF( iom_use('icediv') ) CALL iom_put( "icediv" , zdivu_i(:,:) ) 158 159 DEALLOCATE( zdivu_i ) 160 135 CALL ice_dyn_adv ( kt ) ! -- advection of ice 136 ! 161 137 CASE ( np_dynADV2D ) !== pure advection ==! (2D w prescribed velocities) 162 ALLOCATE( zdivu_i(jpi,jpj) )138 ! 163 139 u_ice(:,:) = rn_uice * umask(:,:,1) 164 140 v_ice(:,:) = rn_vice * vmask(:,:,1) … … 166 142 !CALL RANDOM_NUMBER(v_ice(:,:)) ; v_ice(:,:) = v_ice(:,:) * 0.1 + rn_vice * 0.9 * vmask(:,:,1) 167 143 ! --- 168 CALL ice_dyn_adv ( kt ) ! -- advection of ice 169 170 ! diagnostics: divergence at T points 171 DO jj = 2, jpjm1 172 DO ji = 2, jpim1 173 zdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) & 174 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj) 144 CALL ice_dyn_adv ( kt ) ! -- advection of ice 145 146 END SELECT 147 ! 148 ! 149 ! diagnostics: divergence at T points 150 IF( iom_use('icediv') ) THEN 151 ! 152 SELECT CASE( nice_dyn ) 153 154 CASE ( np_dynADV1D , np_dynADV2D ) 155 156 ALLOCATE( zdivu_i(jpi,jpj) ) 157 DO jj = 2, jpjm1 158 DO ji = 2, jpim1 159 zdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) & 160 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj) 161 END DO 175 162 END DO 176 END DO177 CALL lbc_lnk( 'icedyn', zdivu_i, 'T', 1.)178 IF( iom_use('icediv') ) CALL iom_put( "icediv" , zdivu_i(:,:))179 180 DEALLOCATE( zdivu_i )181 182 END SELECT163 CALL lbc_lnk( 'icedyn', zdivu_i, 'T', 1. ) 164 CALL iom_put( "icediv" , zdivu_i(:,:) ) 165 DEALLOCATE( zdivu_i ) 166 167 END SELECT 168 ! 169 ENDIF 183 170 ! 184 171 ! controls … … 186 173 ! 187 174 END SUBROUTINE ice_dyn 188 189 190 SUBROUTINE Hbig( phi_max, phs_max, phip_max )191 !!-------------------------------------------------------------------192 !! *** ROUTINE Hbig ***193 !!194 !! ** Purpose : Thickness correction in case advection scheme creates195 !! abnormally tick ice or snow196 !!197 !! ** Method : 1- check whether ice thickness is larger than the surrounding 9-points198 !! (before advection) and reduce it by adapting ice concentration199 !! 2- check whether snow thickness is larger than the surrounding 9-points200 !! (before advection) and reduce it by sending the excess in the ocean201 !! 3- check whether snow load deplets the snow-ice interface below sea level$202 !! and reduce it by sending the excess in the ocean203 !! 4- correct pond fraction to avoid a_ip > a_i204 !!205 !! ** input : Max thickness of the surrounding 9-points206 !!-------------------------------------------------------------------207 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: phi_max, phs_max, phip_max ! max ice thick from surrounding 9-pts208 !209 INTEGER :: ji, jj, jl ! dummy loop indices210 REAL(wp) :: zhip, zhi, zhs, zvs_excess, zfra211 !!-------------------------------------------------------------------212 ! controls213 IF( ln_icediachk ) CALL ice_cons_hsm(0, 'Hbig', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation214 !215 CALL ice_var_zapsmall !-- zap small areas216 !217 DO jl = 1, jpl218 DO jj = 1, jpj219 DO ji = 1, jpi220 IF ( v_i(ji,jj,jl) > 0._wp ) THEN221 !222 ! ! -- check h_ip -- !223 ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip224 IF( ln_pnd_H12 .AND. v_ip(ji,jj,jl) > 0._wp ) THEN225 zhip = v_ip(ji,jj,jl) / MAX( epsi20, a_ip(ji,jj,jl) )226 IF( zhip > phip_max(ji,jj,jl) .AND. a_ip(ji,jj,jl) < 0.15 ) THEN227 a_ip(ji,jj,jl) = v_ip(ji,jj,jl) / phip_max(ji,jj,jl)228 ENDIF229 ENDIF230 !231 ! ! -- check h_i -- !232 ! if h_i is larger than the surrounding 9 pts => reduce h_i and increase a_i233 zhi = v_i(ji,jj,jl) / a_i(ji,jj,jl)234 IF( zhi > phi_max(ji,jj,jl) .AND. a_i(ji,jj,jl) < 0.15 ) THEN235 a_i(ji,jj,jl) = v_i(ji,jj,jl) / MIN( phi_max(ji,jj,jl), hi_max(jpl) ) !-- bound h_i to hi_max (99 m)236 ENDIF237 !238 ! ! -- check h_s -- !239 ! if h_s is larger than the surrounding 9 pts => put the snow excess in the ocean240 zhs = v_s(ji,jj,jl) / a_i(ji,jj,jl)241 IF( v_s(ji,jj,jl) > 0._wp .AND. zhs > phs_max(ji,jj,jl) .AND. a_i(ji,jj,jl) < 0.15 ) THEN242 zfra = phs_max(ji,jj,jl) / MAX( zhs, epsi20 )243 !244 wfx_res(ji,jj) = wfx_res(ji,jj) + ( v_s(ji,jj,jl) - a_i(ji,jj,jl) * phs_max(ji,jj,jl) ) * rhos * r1_rdtice245 hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( e_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * r1_rdtice ! W.m-2 <0246 !247 e_s(ji,jj,1:nlay_s,jl) = e_s(ji,jj,1:nlay_s,jl) * zfra248 v_s(ji,jj,jl) = a_i(ji,jj,jl) * phs_max(ji,jj,jl)249 ENDIF250 !251 ! ! -- check snow load -- !252 ! if snow load makes snow-ice interface to deplet below the ocean surface => put the snow excess in the ocean253 ! this correction is crucial because of the call to routine icecor afterwards which imposes a mini of ice thick. (rn_himin)254 ! this imposed mini can artificially make the snow very thick (if concentration decreases drastically)255 zvs_excess = MAX( 0._wp, v_s(ji,jj,jl) - v_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos )256 IF( zvs_excess > 0._wp ) THEN257 zfra = ( v_s(ji,jj,jl) - zvs_excess ) / MAX( v_s(ji,jj,jl), epsi20 )258 wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * r1_rdtice259 hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( e_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * r1_rdtice ! W.m-2 <0260 !261 e_s(ji,jj,1:nlay_s,jl) = e_s(ji,jj,1:nlay_s,jl) * zfra262 v_s(ji,jj,jl) = v_s(ji,jj,jl) - zvs_excess263 ENDIF264 265 ENDIF266 END DO267 END DO268 END DO269 ! !-- correct pond fraction to avoid a_ip > a_i270 WHERE( a_ip(:,:,:) > a_i(:,:,:) ) a_ip(:,:,:) = a_i(:,:,:)271 !272 ! controls273 IF( ln_icediachk ) CALL ice_cons_hsm(1, 'Hbig', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation274 !275 END SUBROUTINE Hbig276 175 277 176 … … 337 236 WRITE(numout,*) '~~~~~~~~~~~~' 338 237 WRITE(numout,*) ' Namelist namdyn:' 339 WRITE(numout,*) ' Full ice dynamics (rhg + adv + ridge/raft + corr) 340 WRITE(numout,*) ' No ridge/raft & No cor (rhg + adv) 341 WRITE(numout,*) ' Advection 1D only (Schar & Smolarkiewicz 1996) 342 WRITE(numout,*) ' Advection 2D only (rn_uvice + adv) 343 WRITE(numout,*) ' with prescribed velocity given by (u,v)_ice = (rn_uice,rn_vice) = (', rn_uice,',',rn_vice,')'344 WRITE(numout,*) ' lateral boundary condition for sea ice dynamics 345 WRITE(numout,*) ' Landfast: param from Lemieux 2016 346 WRITE(numout,*) ' Landfast: param from home made 347 WRITE(numout,*) ' fraction of ocean depth that ice must reach 348 WRITE(numout,*) ' maximum bottom stress per unit area of contact 349 WRITE(numout,*) ' relax time scale (s-1) to reach static friction 350 WRITE(numout,*) ' isotropic tensile strength 238 WRITE(numout,*) ' Full ice dynamics (rhg + adv + ridge/raft + corr) ln_dynALL = ', ln_dynALL 239 WRITE(numout,*) ' No ridge/raft & No cor (rhg + adv) ln_dynRHGADV = ', ln_dynRHGADV 240 WRITE(numout,*) ' Advection 1D only (Schar & Smolarkiewicz 1996) ln_dynADV1D = ', ln_dynADV1D 241 WRITE(numout,*) ' Advection 2D only (rn_uvice + adv) ln_dynADV2D = ', ln_dynADV2D 242 WRITE(numout,*) ' with prescribed velocity given by (u,v)_ice = (rn_uice,rn_vice) = (', rn_uice,',',rn_vice,')' 243 WRITE(numout,*) ' lateral boundary condition for sea ice dynamics rn_ishlat = ', rn_ishlat 244 WRITE(numout,*) ' Landfast: param from Lemieux 2016 ln_landfast_L16 = ', ln_landfast_L16 245 WRITE(numout,*) ' Landfast: param from home made ln_landfast_home= ', ln_landfast_home 246 WRITE(numout,*) ' fraction of ocean depth that ice must reach rn_depfra = ', rn_depfra 247 WRITE(numout,*) ' maximum bottom stress per unit area of contact rn_icebfr = ', rn_icebfr 248 WRITE(numout,*) ' relax time scale (s-1) to reach static friction rn_lfrelax = ', rn_lfrelax 249 WRITE(numout,*) ' isotropic tensile strength rn_tensile = ', rn_tensile 351 250 WRITE(numout,*) 352 251 ENDIF -
NEMO/releases/release-4.0/src/ICE/icedyn_adv.F90
r10413 r10910 64 64 !!---------------------------------------------------------------------- 65 65 INTEGER, INTENT(in) :: kt ! number of iteration 66 !67 INTEGER :: jl ! dummy loop indice68 REAL(wp), DIMENSION(jpi,jpj) :: zmask ! fraction of time step with v_i < 069 66 !!--------------------------------------------------------------------- 70 67 ! 71 IF( ln_timing ) CALL timing_start('icedyn_adv') 68 ! controls 69 IF( ln_timing ) CALL timing_start('icedyn_adv') ! timing 70 IF( ln_icediachk ) CALL ice_cons_hsm(0, 'icedyn_adv', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 72 71 ! 73 72 IF( kt == nit000 .AND. lwp ) THEN … … 76 75 WRITE(numout,*) '~~~~~~~~~~~' 77 76 ENDIF 78 79 ! conservation test 80 IF( ln_icediachk ) CALL ice_cons_hsm(0, 'icedyn_adv', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 81 82 !---------- 83 ! Advection 84 !---------- 77 ! 78 !---------------! 79 !== Advection ==! 80 !---------------! 85 81 SELECT CASE( nice_adv ) 86 82 ! !-----------------------! 87 83 CASE( np_advUMx ) ! ULTIMATE-MACHO scheme ! 88 84 ! !-----------------------! 89 CALL ice_dyn_adv_umx( nn_UMx, kt, u_ice, v_ice, ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i ) 90 ! !-----------------------! 85 CALL ice_dyn_adv_umx( nn_UMx, kt, u_ice, v_ice, h_i, h_s, h_ip, & 86 & ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i ) 87 ! !-----------------------! 91 88 CASE( np_advPRA ) ! PRATHER scheme ! 92 89 ! !-----------------------! 93 CALL ice_dyn_adv_pra( kt, u_ice, v_ice, ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i ) 90 CALL ice_dyn_adv_pra( kt, u_ice, v_ice, & 91 & ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i ) 94 92 END SELECT 95 96 !----------------------------97 ! Debug the advection schemes98 !----------------------------99 ! clem: At least one advection scheme above is not strictly positive => UMx100 ! In Prather, I am not sure if the fields are bounded by 0 or not (it seems yes)101 ! In UMx , advected fields are not perfectly bounded and negative values can appear.102 ! These values are usually very small but in some occasions they can also be non-negligible103 ! Therefore one needs to bound the advected fields by 0 (though this is not a clean fix)104 !105 ! record the negative values resulting from UMx106 zmask(:,:) = 0._wp ! keep the init to 0 here107 DO jl = 1, jpl108 WHERE( v_i(:,:,jl) < 0._wp ) zmask(:,:) = 1._wp109 END DO110 IF( iom_use('iceneg_pres') ) CALL iom_put("iceneg_pres", zmask ) ! fraction of time step with v_i < 0111 IF( iom_use('iceneg_volu') ) CALL iom_put("iceneg_volu", SUM(MIN( v_i, 0. ), dim=3 ) ) ! negative ice volume (only)112 IF( iom_use('iceneg_hfx' ) ) CALL iom_put("iceneg_hfx" , ( SUM(SUM( MIN( e_i(:,:,1:nlay_i,:), 0. ) & ! negative ice heat content (only)113 & , dim=4 ), dim=3 ) )* r1_rdtice ) ! -- eq. heat flux [W/m2]114 !115 ! ==> conservation is ensured by calling this routine below,116 ! however the global ice volume is then changed by advection (but errors are small)117 CALL ice_var_zapneg( ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i )118 93 119 94 !------------ -
NEMO/releases/release-4.0/src/ICE/icedyn_adv_umx.F90
r10785 r10910 11 11 !! 'key_si3' SI3 sea-ice model 12 12 !!---------------------------------------------------------------------- 13 !! ice_dyn_adv_umx : update the tracer trend with the 3D advection trends using a TVD scheme13 !! ice_dyn_adv_umx : update the tracer fields 14 14 !! ultimate_x(_y) : compute a tracer value at velocity points using ULTIMATE scheme at various orders 15 !! macho : ???16 !! nonosc_ice : compute monotonic tracer fluxes bya non-oscillatory algorithm15 !! macho : compute the fluxes 16 !! nonosc_ice : limit the fluxes using a non-oscillatory algorithm 17 17 !!---------------------------------------------------------------------- 18 18 USE phycst ! physical constant … … 39 39 INTEGER :: kn_limiter = 1 40 40 41 ! if T interpolated at u/v points is negative, then interpolate T at u/v points using the upstream scheme 42 ! clem: if set to true, the 2D test case "diagonal advection" does not work (I do not understand why) 43 ! but in realistic cases, it avoids having very negative ice temperature (-50) at low ice concentration 41 ! if there is an outward velocity in a grid cell where there is no ice initially (typically at the ice edge), 42 ! interpolated T at u/v points can be non-zero while it should 43 ! (because of the high order of the advection scheme). Thus set it to 0 in this case 44 LOGICAL :: ll_icedge = .TRUE. 45 46 ! if T interpolated at u/v points is negative or v_i < 1.e-6 47 ! then interpolate T at u/v points using the upstream scheme 44 48 LOGICAL :: ll_neg = .TRUE. 45 49 … … 51 55 52 56 ! prelimiter: use it to avoid overshoot in H 53 ! clem: if set to true, the 2D test case "diagnoal advection" does not work (I do not understand why)54 57 LOGICAL :: ll_prelimiter_zalesak = .FALSE. ! from: Zalesak(1979) eq. 14 => better for 1D. Not well defined in 2D 55 58 59 ! advection for S and T: dVS/dt = -div( uA * uHS / u ) => ll_advS = F 60 ! or dVS/dt = -div( uV * uS / u ) => ll_advS = T 61 LOGICAL :: ll_advS = .FALSE. 56 62 57 63 !! * Substitutions … … 64 70 CONTAINS 65 71 66 SUBROUTINE ice_dyn_adv_umx( kn_umx, kt, pu_ice, pv_ice, &72 SUBROUTINE ice_dyn_adv_umx( kn_umx, kt, pu_ice, pv_ice, ph_i, ph_s, ph_ip, & 67 73 & pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 68 74 !!---------------------------------------------------------------------- … … 79 85 REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pu_ice ! ice i-velocity 80 86 REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pv_ice ! ice j-velocity 87 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: ph_i ! ice thickness 88 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: ph_s ! snw thickness 89 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: ph_ip ! ice pond thickness 81 90 REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pato_i ! open water area 82 91 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i ! ice volume … … 94 103 REAL(wp) :: zamsk ! 1 if advection of concentration, 0 if advection of other tracers 95 104 REAL(wp) :: zdt 96 REAL(wp), DIMENSION(1) :: zcflprv, zcflnow ! send zcflnow and receive zcflprv105 REAL(wp), DIMENSION(1) :: zcflprv, zcflnow ! for global communication 97 106 REAL(wp), DIMENSION(jpi,jpj) :: zudy, zvdx, zcu_box, zcv_box 98 107 REAL(wp), DIMENSION(jpi,jpj) :: zati1, zati2 99 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zua_ho, zva_ho 100 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_ai, z1_aip 108 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zua_ho, zva_ho, zua_ups, zva_ups 109 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zuv_ho, zvv_ho, zuv_ups, zvv_ups 110 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_ai, z1_aip, z1_vi, z1_vs 101 111 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zhvar 112 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zhi_max, zhs_max, zhip_max 102 113 !!---------------------------------------------------------------------- 103 114 ! 104 115 IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_dyn_adv_umx: Ultimate-Macho advection scheme' 105 116 ! 106 ! --- If ice drift field is too fast, use an appropriate time step for advection (CFL test for stability) --- ! 107 ! When needed, the advection split is applied at the next time-step in order to avoid blocking global comm. 108 ! ...this should not affect too much the stability... Was ok on the tests we did... 117 ! --- Record max of the surrounding 9-pts ice thick. (for call Hbig) --- ! 118 DO jl = 1, jpl 119 DO jj = 2, jpjm1 120 DO ji = fs_2, fs_jpim1 121 zhip_max(ji,jj,jl) = MAX( epsi20, ph_ip(ji,jj,jl), ph_ip(ji+1,jj ,jl), ph_ip(ji ,jj+1,jl), & 122 & ph_ip(ji-1,jj ,jl), ph_ip(ji ,jj-1,jl), & 123 & ph_ip(ji+1,jj+1,jl), ph_ip(ji-1,jj-1,jl), & 124 & ph_ip(ji+1,jj-1,jl), ph_ip(ji-1,jj+1,jl) ) 125 zhi_max (ji,jj,jl) = MAX( epsi20, ph_i (ji,jj,jl), ph_i (ji+1,jj ,jl), ph_i (ji ,jj+1,jl), & 126 & ph_i (ji-1,jj ,jl), ph_i (ji ,jj-1,jl), & 127 & ph_i (ji+1,jj+1,jl), ph_i (ji-1,jj-1,jl), & 128 & ph_i (ji+1,jj-1,jl), ph_i (ji-1,jj+1,jl) ) 129 zhs_max (ji,jj,jl) = MAX( epsi20, ph_s (ji,jj,jl), ph_s (ji+1,jj ,jl), ph_s (ji ,jj+1,jl), & 130 & ph_s (ji-1,jj ,jl), ph_s (ji ,jj-1,jl), & 131 & ph_s (ji+1,jj+1,jl), ph_s (ji-1,jj-1,jl), & 132 & ph_s (ji+1,jj-1,jl), ph_s (ji-1,jj+1,jl) ) 133 END DO 134 END DO 135 END DO 136 CALL lbc_lnk_multi( 'icedyn_adv_umx', zhi_max(:,:,:), 'T', 1., zhs_max(:,:,:), 'T', 1., zhip_max(:,:,:), 'T', 1. ) 137 ! 138 ! 139 ! --- If ice drift is too fast, use subtime steps for advection (CFL test for stability) --- ! 140 ! Note: the advection split is applied at the next time-step in order to avoid blocking global comm. 141 ! this should not affect too much the stability 109 142 zcflnow(1) = MAXVAL( ABS( pu_ice(:,:) ) * rdt_ice * r1_e1u(:,:) ) 110 143 zcflnow(1) = MAX( zcflnow(1), MAXVAL( ABS( pv_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) ) … … 116 149 ELSE ; icycle = 1 117 150 ENDIF 118 119 151 zdt = rdt_ice / REAL(icycle) 120 152 … … 154 186 END WHERE 155 187 ! 156 ! set u* a=u for advection of A only188 ! set u*A=u for advection of A only 157 189 DO jl = 1, jpl 158 190 zua_ho(:,:,jl) = zudy(:,:) 159 191 zva_ho(:,:,jl) = zvdx(:,:) 160 192 END DO 161 193 194 !== Ice area ==! 162 195 zamsk = 1._wp 163 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, pa_i, pa_i, zua_ho, zva_ho ) !-- Ice area 196 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho , zva_ho , zcu_box, zcv_box, & 197 & pa_i, pa_i, zua_ups, zva_ups, zua_ho , zva_ho ) 164 198 zamsk = 0._wp 165 ! 166 zhvar(:,:,:) = pv_i(:,:,:) * z1_ai(:,:,:) 167 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pv_i ) !-- Ice volume 168 ! 169 zhvar(:,:,:) = pv_s(:,:,:) * z1_ai(:,:,:) 170 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pv_s ) !-- Snw volume 171 ! 172 zhvar(:,:,:) = psv_i(:,:,:) * z1_ai(:,:,:) 173 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, psv_i ) !-- Salt content 174 ! 175 DO jk = 1, nlay_i 176 zhvar(:,:,:) = pe_i(:,:,jk,:) * z1_ai(:,:,:) 177 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pe_i(:,:,jk,:) ) !-- Ice heat content 178 END DO 179 ! 180 DO jk = 1, nlay_s 181 zhvar(:,:,:) = pe_s(:,:,jk,:) * z1_ai(:,:,:) 182 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pe_s(:,:,jk,:) ) !-- Snw heat content 183 END DO 184 ! 185 IF( iom_use('iceage') .OR. iom_use('iceage_cat') ) THEN !-- Ice Age 199 200 IF( .NOT. ll_advS ) THEN !-- advection form: uA * uHS / u 201 !== Ice volume ==! 202 zhvar(:,:,:) = pv_i(:,:,:) * z1_ai(:,:,:) 203 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, & 204 & zhvar, pv_i, zua_ups, zva_ups ) 205 ! 206 !== Snw volume ==! 207 zhvar(:,:,:) = pv_s(:,:,:) * z1_ai(:,:,:) 208 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, & 209 & zhvar, pv_s, zua_ups, zva_ups ) 210 ! 211 !== Salt content ==! 212 zhvar(:,:,:) = psv_i(:,:,:) * z1_ai(:,:,:) 213 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, & 214 & zhvar, psv_i, zua_ups, zva_ups ) 215 ! 216 !== Ice heat content ==! 217 DO jk = 1, nlay_i 218 zhvar(:,:,:) = pe_i(:,:,jk,:) * z1_ai(:,:,:) 219 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho, zva_ho, zcu_box, zcv_box, & 220 & zhvar, pe_i(:,:,jk,:), zua_ups, zva_ups ) 221 END DO 222 ! 223 !== Snw heat content ==! 224 DO jk = 1, nlay_s 225 zhvar(:,:,:) = pe_s(:,:,jk,:) * z1_ai(:,:,:) 226 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho, zva_ho, zcu_box, zcv_box, & 227 & zhvar, pe_s(:,:,jk,:), zua_ups, zva_ups ) 228 END DO 229 ! 230 ELSE !-- advection form: uV * uS / u 231 ! inverse of Vi 232 WHERE( pv_i(:,:,:) >= epsi20 ) ; z1_vi(:,:,:) = 1._wp / pv_i(:,:,:) 233 ELSEWHERE ; z1_vi(:,:,:) = 0. 234 END WHERE 235 ! inverse of Vs 236 WHERE( pv_s(:,:,:) >= epsi20 ) ; z1_vs(:,:,:) = 1._wp / pv_s(:,:,:) 237 ELSEWHERE ; z1_vs(:,:,:) = 0. 238 END WHERE 239 ! 240 ! It is important to first calculate the ice fields and then the snow fields (because we use the same arrays) 241 ! 242 !== Ice volume ==! 243 zuv_ups = zua_ups 244 zvv_ups = zva_ups 245 zhvar(:,:,:) = pv_i(:,:,:) * z1_ai(:,:,:) 246 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, & 247 & zhvar, pv_i, zuv_ups, zvv_ups, zuv_ho , zvv_ho ) 248 ! 249 !== Salt content ==! 250 zhvar(:,:,:) = psv_i(:,:,:) * z1_vi(:,:,:) 251 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zuv_ho , zvv_ho , zcu_box, zcv_box, & 252 & zhvar, psv_i, zuv_ups, zvv_ups ) 253 ! 254 !== Ice heat content ==! 255 DO jk = 1, nlay_i 256 zhvar(:,:,:) = pe_i(:,:,jk,:) * z1_vi(:,:,:) 257 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zuv_ho, zvv_ho, zcu_box, zcv_box, & 258 & zhvar, pe_i(:,:,jk,:), zuv_ups, zvv_ups ) 259 END DO 260 ! 261 !== Snow volume ==! 262 zuv_ups = zua_ups 263 zvv_ups = zva_ups 264 zhvar(:,:,:) = pv_s(:,:,:) * z1_ai(:,:,:) 265 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zua_ho , zva_ho , zcu_box, zcv_box, & 266 & zhvar, pv_s, zuv_ups, zvv_ups, zuv_ho , zvv_ho ) 267 ! 268 !== Snw heat content ==! 269 DO jk = 1, nlay_s 270 zhvar(:,:,:) = pe_s(:,:,jk,:) * z1_vs(:,:,:) 271 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx, zuv_ho, zvv_ho, zcu_box, zcv_box, & 272 & zhvar, pe_s(:,:,jk,:), zuv_ups, zvv_ups ) 273 END DO 274 ! 275 ENDIF 276 ! 277 !== Ice age ==! 278 IF( iom_use('iceage') .OR. iom_use('iceage_cat') ) THEN 186 279 ! clem: in theory we should use the formulation below to advect the ice age, but the code is unable to deal with 187 280 ! fields that do not depend on volume (here oa_i depends on concentration). It creates abnormal ages that … … 189 282 !!zhvar(:,:,:) = poa_i(:,:,:) * z1_ai(:,:,:) 190 283 !!CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, poa_i ) 191 ! set u* a=u for advection of ice age284 ! set u*A=u for advection of ice age 192 285 DO jl = 1, jpl 193 286 zua_ho(:,:,jl) = zudy(:,:) … … 195 288 END DO 196 289 zamsk = 1._wp 197 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, poa_i, poa_i ) 290 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho, zva_ho, zcu_box, zcv_box, & 291 & poa_i, poa_i ) 198 292 zamsk = 0._wp 199 293 ENDIF 200 294 ! 201 IF ( ln_pnd_H12 ) THEN !-- melt ponds 202 ! set u*a=u for advection of Ap only 295 !== melt ponds ==! 296 IF ( ln_pnd_H12 ) THEN 297 ! set u*A=u for advection of Ap only 203 298 DO jl = 1, jpl 204 299 zua_ho(:,:,jl) = zudy(:,:) 205 300 zva_ho(:,:,jl) = zvdx(:,:) 206 301 END DO 207 ! 302 ! fraction 208 303 zamsk = 1._wp 209 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, pa_ip, pa_ip, zua_ho, zva_ho ) ! fraction 304 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, & 305 & pa_ip, pa_ip, zua_ups, zva_ups, zua_ho , zva_ho ) 210 306 zamsk = 0._wp 211 ! 307 ! volume 212 308 zhvar(:,:,:) = pv_ip(:,:,:) * z1_aip(:,:,:) 213 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy, zvdx, zua_ho, zva_ho, zcu_box, zcv_box, zhvar, pv_ip ) ! volume 309 CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, & 310 & zhvar, pv_ip, zua_ups, zva_ups ) 214 311 ENDIF 215 312 ! 313 !== Open water area ==! 216 314 zati2(:,:) = SUM( pa_i(:,:,:), dim=3 ) 217 315 DO jj = 2, jpjm1 218 316 DO ji = fs_2, fs_jpim1 219 pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) & !-- Open water area317 pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) & 220 318 & - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt 221 319 END DO … … 223 321 CALL lbc_lnk( 'icedyn_adv_umx', pato_i(:,:), 'T', 1. ) 224 322 ! 323 ! 324 ! --- Ensure non-negative fields and in-bound thicknesses --- ! 325 ! 326 ! Remove negative values (conservation is ensured) 327 ! (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20) 328 CALL ice_var_zapneg( pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 329 ! 330 ! Make sure ice thickness is not too big 331 ! (because ice thickness can be too large where ice concentration is very small) 332 CALL Hbig( zhi_max, zhs_max, zhip_max, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 333 225 334 END DO 226 335 ! … … 228 337 229 338 230 SUBROUTINE adv_umx( pamsk, kn_umx, jt, kt, pdt, pu, pv, puc, pvc, pubox, pvbox, pt, ptc, pua_ho, pva_ho ) 339 SUBROUTINE adv_umx( pamsk, kn_umx, jt, kt, pdt, pu, pv, puc, pvc, pubox, pvbox, & 340 & pt, ptc, pua_ups, pva_ups, pua_ho, pva_ho ) 231 341 !!---------------------------------------------------------------------- 232 342 !! *** ROUTINE adv_umx *** … … 235 345 !! tracers and add it to the general trend of tracer equations 236 346 !! 237 !! ** Method : - calculate upstream fluxes and upstream solution for tracer H347 !! ** Method : - calculate upstream fluxes and upstream solution for tracers V/A(=H) etc 238 348 !! - calculate tracer H at u and v points (Ultimate) 239 !! - calculate the high order fluxes using alterning directions (Macho ?)349 !! - calculate the high order fluxes using alterning directions (Macho) 240 350 !! - apply a limiter on the fluxes (nonosc_ice) 241 !! - convert this tracer flux to a tracer content flux (uH -> uV) 242 !! - calculate the high order solution for tracer content V 351 !! - convert this tracer flux to a "volume" flux (uH -> uV) 352 !! - apply a limiter a second time on the volumes fluxes (nonosc_ice) 353 !! - calculate the high order solution for V 243 354 !! 244 !! ** Action : solve 2 equations => a) da/dt = -div(ua) 245 !! b) dV/dt = -div(uV) using dH/dt = -u.grad(H) 246 !! in eq. b), - fluxes uH are evaluated (with UMx) and limited (with nonosc_ice). This step is necessary to get a good H. 247 !! - then we convert this flux to a "volume" flux this way => uH*ua/u 248 !! where ua is the flux from eq. a) 249 !! - at last we estimate dV/dt = -div(uH*ua/u) 355 !! ** Action : solve 3 equations => a) dA/dt = -div(uA) 356 !! b) dV/dt = -div(uV) using dH/dt = -u.grad(H) 357 !! c) dVS/dt = -div(uVS) using either dHS/dt = -u.grad(HS) or dS/dt = -u.grad(S) 250 358 !! 251 !! ** Note : - this method can lead to small negative V (since we only limit H) => corrected in icedyn_adv.F90 conserving mass etc. 252 !! - negative tracers at u-v points can also occur from the Ultimate scheme (usually at the ice edge) and the solution for now 253 !! is to apply an upstream scheme when it occurs. A better solution would be to degrade the order of 254 !! the scheme automatically by applying a mask of the ice cover inside Ultimate (not done). 359 !! in eq. b), - fluxes uH are evaluated (with UMx) and limited with nonosc_ice. This step is necessary to get a good H. 360 !! - then we convert this flux to a "volume" flux this way => uH * uA / u 361 !! where uA is the flux from eq. a) 362 !! this "volume" flux is also limited with nonosc_ice (otherwise overshoots can occur) 363 !! - at last we estimate dV/dt = -div(uH * uA / u) 364 !! 365 !! in eq. c), one can solve the equation for S (ln_advS=T), then dVS/dt = -div(uV * uS / u) 366 !! or for HS (ln_advS=F), then dVS/dt = -div(uA * uHS / u) 367 !! 368 !! ** Note : - this method can lead to tiny negative V (-1.e-20) => set it to 0 while conserving mass etc. 369 !! - At the ice edge, Ultimate scheme can lead to: 370 !! 1) negative interpolated tracers at u-v points 371 !! 2) non-zero interpolated tracers at u-v points eventhough there is no ice and velocity is outward 372 !! Solution for 1): apply an upstream scheme when it occurs. A better solution would be to degrade the order of 373 !! the scheme automatically by applying a mask of the ice cover inside Ultimate (not done). 374 !! Solution for 2): we set it to 0 in this case 255 375 !! - Eventhough 1D tests give very good results (typically the one from Schar & Smolarkiewiecz), the 2D is less good. 256 376 !! Large values of H can appear for very small ice concentration, and when it does it messes the things up since we 257 !! work on H (and not V). It probably comes from the prelimiter of zalesak which is coded for 1D and not 2D.377 !! work on H (and not V). It is partly related to the multi-category approach 258 378 !! Therefore, after advection we limit the thickness to the largest value of the 9-points around (only if ice 259 !! concentration is small). 260 !! To-do: expand the prelimiter from zalesak to make it work in 2D 261 !!---------------------------------------------------------------------- 262 REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0) 263 INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2) 264 INTEGER , INTENT(in ) :: jt ! number of sub-iteration 265 INTEGER , INTENT(in ) :: kt ! number of iteration 266 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 267 REAL(wp), DIMENSION(:,: ) , INTENT(in ) :: pu , pv ! 2 ice velocity components => u*e2 268 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: puc , pvc ! 2 ice velocity components => u*e2 or u*a*e2u 269 REAL(wp), DIMENSION(:,: ) , INTENT(in ) :: pubox, pvbox ! upstream velocity 270 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pt ! tracer field 271 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: ptc ! tracer content field 272 REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out), OPTIONAL :: pua_ho, pva_ho ! high order u*a fluxes 379 !! concentration is small). Since we do not limit S and T, large values can occur at the edge but it does not really matter 380 !! since sv_i and e_i are still good. 381 !!---------------------------------------------------------------------- 382 REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0) 383 INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2) 384 INTEGER , INTENT(in ) :: jt ! number of sub-iteration 385 INTEGER , INTENT(in ) :: kt ! number of iteration 386 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 387 REAL(wp), DIMENSION(:,: ) , INTENT(in ) :: pu , pv ! 2 ice velocity components => u*e2 388 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: puc , pvc ! 2 ice velocity components => u*e2 or u*a*e2u 389 REAL(wp), DIMENSION(:,: ) , INTENT(in ) :: pubox, pvbox ! upstream velocity 390 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pt ! tracer field 391 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: ptc ! tracer content field 392 REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(inout), OPTIONAL :: pua_ups, pva_ups ! upstream u*a fluxes 393 REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out), OPTIONAL :: pua_ho, pva_ho ! high order u*a fluxes 273 394 ! 274 395 INTEGER :: ji, jj, jl ! dummy loop indices … … 303 424 DO jj = 1, jpjm1 304 425 DO ji = 1, fs_jpim1 305 IF( ABS( pu c(ji,jj,jl) ) > 0._wp .AND. ABS( pu(ji,jj) ) > 0._wp) THEN306 zfu_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) * puc (ji,jj,jl) / pu(ji,jj)307 zfu_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) * pu c(ji,jj,jl) / pu(ji,jj)426 IF( ABS( pu(ji,jj) ) > epsi10 ) THEN 427 zfu_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) * puc (ji,jj,jl) / pu(ji,jj) 428 zfu_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) * pua_ups(ji,jj,jl) / pu(ji,jj) 308 429 ELSE 309 430 zfu_ho (ji,jj,jl) = 0._wp … … 311 432 ENDIF 312 433 ! 313 IF( ABS( pv c(ji,jj,jl) ) > 0._wp .AND. ABS( pv(ji,jj) ) > 0._wp) THEN314 zfv_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) * pvc (ji,jj,jl) / pv(ji,jj)315 zfv_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) * pv c(ji,jj,jl) / pv(ji,jj)434 IF( ABS( pv(ji,jj) ) > epsi10 ) THEN 435 zfv_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) * pvc (ji,jj,jl) / pv(ji,jj) 436 zfv_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) * pva_ups(ji,jj,jl) / pv(ji,jj) 316 437 ELSE 317 438 zfv_ho (ji,jj,jl) = 0._wp … … 321 442 END DO 322 443 END DO 444 445 ! the new "volume" fluxes must also be "flux corrected" 446 ! thus we calculate the upstream solution and apply a limiter again 447 DO jl = 1, jpl 448 DO jj = 2, jpjm1 449 DO ji = fs_2, fs_jpim1 450 ztra = - ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj,jl) + zfv_ups(ji,jj,jl) - zfv_ups(ji,jj-1,jl) ) 451 ! 452 zt_ups(ji,jj,jl) = ( ptc(ji,jj,jl) + ztra * r1_e1e2t(ji,jj) * pdt ) * tmask(ji,jj,1) 453 END DO 454 END DO 455 END DO 456 CALL lbc_lnk( 'icedyn_adv_umx', zt_ups, 'T', 1. ) 457 ! 458 IF ( kn_limiter == 1 ) THEN 459 CALL nonosc_ice( 1._wp, pdt, pu, pv, ptc, zt_ups, zfu_ups, zfv_ups, zfu_ho, zfv_ho ) 460 ELSEIF( kn_limiter == 2 .OR. kn_limiter == 3 ) THEN 461 CALL limiter_x( pdt, pu, ptc, zfu_ups, zfu_ho ) 462 CALL limiter_y( pdt, pv, ptc, zfv_ups, zfv_ho ) 463 ENDIF 464 ! 323 465 ENDIF 324 ! --ho 325 ! in case of advection of A: output u*a 326 ! ------------------------------------- 466 ! --ho --ups 467 ! in case of advection of A: output u*a and u*a 468 ! ----------------------------------------------- 327 469 IF( PRESENT( pua_ho ) ) THEN 328 470 DO jl = 1, jpl 329 471 DO jj = 1, jpjm1 330 472 DO ji = 1, fs_jpim1 331 pua_ho (ji,jj,jl) = zfu_ho(ji,jj,jl)332 p va_ho(ji,jj,jl) = zfv_ho(ji,jj,jl)333 473 pua_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) ; pva_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) 474 pua_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) ; pva_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) 475 END DO 334 476 END DO 335 477 END DO … … 609 751 ! 610 752 ! !-- ultimate interpolation of pt at u-point --! 611 CALL ultimate_x( kn_umx, pdt, pt, pu, zt_u, pfu_ho )753 CALL ultimate_x( pamsk, kn_umx, pdt, pt, pu, zt_u, pfu_ho ) 612 754 ! !-- limiter in x --! 613 755 IF( kn_limiter == 2 .OR. kn_limiter == 3 ) CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) … … 619 761 & + pt (ji,jj,jl) * ( pu (ji,jj ) - pu (ji-1,jj ) ) * r1_e1e2t(ji,jj) & 620 762 & * pamsk & 621 & ) * pdt ) * tmask(ji,jj,1) 763 & ) * pdt ) * tmask(ji,jj,1) 764 !!clem test 765 !!zpt(ji,jj,jl) = MAX( 0._wp, zpt(ji,jj,jl) ) 766 !!clem test 622 767 END DO 623 768 END DO … … 627 772 ! !-- ultimate interpolation of pt at v-point --! 628 773 IF( ll_hoxy ) THEN 629 CALL ultimate_y( kn_umx, pdt, zpt, pv, zt_v, pfv_ho )774 CALL ultimate_y( pamsk, kn_umx, pdt, zpt, pv, zt_v, pfv_ho ) 630 775 ELSE 631 CALL ultimate_y( kn_umx, pdt, pt , pv, zt_v, pfv_ho )776 CALL ultimate_y( pamsk, kn_umx, pdt, pt , pv, zt_v, pfv_ho ) 632 777 ENDIF 633 778 ! !-- limiter in y --! … … 638 783 ! 639 784 ! !-- ultimate interpolation of pt at v-point --! 640 CALL ultimate_y( kn_umx, pdt, pt, pv, zt_v, pfv_ho )785 CALL ultimate_y( pamsk, kn_umx, pdt, pt, pv, zt_v, pfv_ho ) 641 786 ! !-- limiter in y --! 642 787 IF( kn_limiter == 2 .OR. kn_limiter == 3 ) CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) … … 649 794 & * pamsk & 650 795 & ) * pdt ) * tmask(ji,jj,1) 796 !!clem test 797 !!zpt(ji,jj,jl) = MAX( 0._wp, zpt(ji,jj,jl) ) 798 !!clem test 651 799 END DO 652 800 END DO … … 656 804 ! !-- ultimate interpolation of pt at u-point --! 657 805 IF( ll_hoxy ) THEN 658 CALL ultimate_x( kn_umx, pdt, zpt, pu, zt_u, pfu_ho )806 CALL ultimate_x( pamsk, kn_umx, pdt, zpt, pu, zt_u, pfu_ho ) 659 807 ELSE 660 CALL ultimate_x( kn_umx, pdt, pt , pu, zt_u, pfu_ho )808 CALL ultimate_x( pamsk, kn_umx, pdt, pt , pu, zt_u, pfu_ho ) 661 809 ENDIF 662 810 ! !-- limiter in x --! … … 670 818 671 819 672 SUBROUTINE ultimate_x( kn_umx, pdt, pt, pu, pt_u, pfu_ho )820 SUBROUTINE ultimate_x( pamsk, kn_umx, pdt, pt, pu, pt_u, pfu_ho ) 673 821 !!--------------------------------------------------------------------- 674 822 !! *** ROUTINE ultimate_x *** … … 680 828 !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 681 829 !!---------------------------------------------------------------------- 830 REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0) 682 831 INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2) 683 832 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step … … 688 837 ! 689 838 INTEGER :: ji, jj, jl ! dummy loop indices 690 REAL(wp) :: zcu, zdx2, zdx4 ! - -839 REAL(wp) :: zcu, zdx2, zdx4, zvi_cen2 ! - - 691 840 REAL(wp), DIMENSION(jpi,jpj,jpl) :: ztu1, ztu2, ztu3, ztu4 692 841 !!---------------------------------------------------------------------- … … 799 948 END SELECT 800 949 ! 950 ! if there is an outward velocity in a grid cell where there is no ice initially (typically at the ice edge), 951 ! interpolated T at u/v points can be non-zero while it should 952 ! (because of the high order of the advection scheme). Thus set it to 0 in this case 953 IF( ll_icedge ) THEN 954 DO jl = 1, jpl 955 DO jj = 1, jpjm1 956 DO ji = 1, fs_jpim1 957 IF( pt(ji,jj,jl) <= 0._wp .AND. pu(ji,jj) >= 0._wp ) THEN 958 pt_u(ji,jj,jl) = 0._wp 959 ENDIF 960 END DO 961 END DO 962 END DO 963 ENDIF 964 ! 801 965 ! if pt at u-point is negative then use the upstream value 802 966 ! this should not be necessary if a proper sea-ice mask is set in Ultimate … … 806 970 DO jj = 1, jpjm1 807 971 DO ji = 1, fs_jpim1 808 IF( pt_u(ji,jj,jl) < 0._wp ) THEN 972 zvi_cen2 = 0.5_wp * ( v_i(ji+1,jj,jl) + v_i(ji,jj,jl) ) 973 IF( pt_u(ji,jj,jl) < 0._wp .OR. ( zvi_cen2 < epsi06 .AND. pamsk == 0._wp ) ) THEN 809 974 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) & 810 975 & - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) … … 826 991 827 992 828 SUBROUTINE ultimate_y( kn_umx, pdt, pt, pv, pt_v, pfv_ho )993 SUBROUTINE ultimate_y( pamsk, kn_umx, pdt, pt, pv, pt_v, pfv_ho ) 829 994 !!--------------------------------------------------------------------- 830 995 !! *** ROUTINE ultimate_y *** … … 836 1001 !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 837 1002 !!---------------------------------------------------------------------- 1003 REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0) 838 1004 INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2) 839 1005 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step … … 844 1010 ! 845 1011 INTEGER :: ji, jj, jl ! dummy loop indices 846 REAL(wp) :: zcv, zdy2, zdy4 ! - -1012 REAL(wp) :: zcv, zdy2, zdy4, zvi_cen2 ! - - 847 1013 REAL(wp), DIMENSION(jpi,jpj,jpl) :: ztv1, ztv2, ztv3, ztv4 848 1014 !!---------------------------------------------------------------------- … … 952 1118 END SELECT 953 1119 ! 1120 ! if there is an outward velocity in a grid cell where there is no ice initially (typically at the ice edge), 1121 ! interpolated T at u/v points can be non-zero while it should 1122 ! (because of the high order of the advection scheme). Thus set it to 0 in this case 1123 IF( ll_icedge ) THEN 1124 DO jl = 1, jpl 1125 DO jj = 1, jpjm1 1126 DO ji = 1, fs_jpim1 1127 IF( pt(ji,jj,jl) <= 0._wp .AND. pv(ji,jj) >= 0._wp ) THEN 1128 pt_v(ji,jj,jl) = 0._wp 1129 ENDIF 1130 END DO 1131 END DO 1132 END DO 1133 ENDIF 1134 ! 954 1135 ! if pt at v-point is negative then use the upstream value 955 1136 ! this should not be necessary if a proper sea-ice mask is set in Ultimate … … 959 1140 DO jj = 1, jpjm1 960 1141 DO ji = 1, fs_jpim1 961 IF( pt_v(ji,jj,jl) < 0._wp ) THEN 1142 zvi_cen2 = 0.5_wp * ( v_i(ji,jj+1,jl) + v_i(ji,jj,jl) ) 1143 IF( pt_v(ji,jj,jl) < 0._wp .OR. ( zvi_cen2 < epsi06 .AND. pamsk == 0._wp ) ) THEN 962 1144 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt(ji,jj+1,jl) + pt(ji,jj,jl) ) & 963 1145 & - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) ) … … 1102 1284 ! 1103 1285 ! ! up & down beta terms 1104 IF( zpos > 0._wp ) THEN ; zbetup(ji,jj,jl) = MAX( 0._wp, zup - pt_ups(ji,jj,jl) ) / zpos * e1e2t(ji,jj) * z1_dt 1105 ELSE ; zbetup(ji,jj,jl) = 0._wp ! zbig 1286 ! clem: zbetup and zbetdo must be 0 for zpos>1.e-10 & zneg>1.e-10 (do not put 0 instead of 1.e-10 !!!) 1287 IF( zpos > epsi10 ) THEN ; zbetup(ji,jj,jl) = MAX( 0._wp, zup - pt_ups(ji,jj,jl) ) / zpos * e1e2t(ji,jj) * z1_dt 1288 ELSE ; zbetup(ji,jj,jl) = 0._wp ! zbig 1106 1289 ENDIF 1107 1290 ! 1108 IF( zneg > 0._wp) THEN ; zbetdo(ji,jj,jl) = MAX( 0._wp, pt_ups(ji,jj,jl) - zdo ) / zneg * e1e2t(ji,jj) * z1_dt1109 ELSE ; zbetdo(ji,jj,jl) = 0._wp ! zbig1291 IF( zneg > epsi10 ) THEN ; zbetdo(ji,jj,jl) = MAX( 0._wp, pt_ups(ji,jj,jl) - zdo ) / zneg * e1e2t(ji,jj) * z1_dt 1292 ELSE ; zbetdo(ji,jj,jl) = 0._wp ! zbig 1110 1293 ENDIF 1111 1294 ! … … 1148 1331 END DO 1149 1332 END DO 1150 1151 ! clem test1152 !! DO jj = 2, jpjm11153 !! DO ji = 2, fs_jpim1 ! vector opt.1154 !! zzt = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) &1155 !! & - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) &1156 !! & + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &1157 !! & + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) &1158 !! & ) * tmask(ji,jj,1)1159 !! IF( zzt < -epsi20 ) THEN1160 !! WRITE(numout,*) 'T<0 nonosc_ice',zzt1161 !! ENDIF1162 !! END DO1163 !! END DO1164 1333 1165 1334 END DO … … 1358 1527 END SUBROUTINE limiter_y 1359 1528 1529 1530 SUBROUTINE Hbig( phi_max, phs_max, phip_max, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 1531 !!------------------------------------------------------------------- 1532 !! *** ROUTINE Hbig *** 1533 !! 1534 !! ** Purpose : Thickness correction in case advection scheme creates 1535 !! abnormally tick ice or snow 1536 !! 1537 !! ** Method : 1- check whether ice thickness is larger than the surrounding 9-points 1538 !! (before advection) and reduce it by adapting ice concentration 1539 !! 2- check whether snow thickness is larger than the surrounding 9-points 1540 !! (before advection) and reduce it by sending the excess in the ocean 1541 !! 3- check whether snow load deplets the snow-ice interface below sea level$ 1542 !! and reduce it by sending the excess in the ocean 1543 !! 4- correct pond fraction to avoid a_ip > a_i 1544 !! 1545 !! ** input : Max thickness of the surrounding 9-points 1546 !!------------------------------------------------------------------- 1547 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: phi_max, phs_max, phip_max ! max ice thick from surrounding 9-pts 1548 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip 1549 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s 1550 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_i 1551 ! 1552 INTEGER :: ji, jj, jk, jl ! dummy loop indices 1553 REAL(wp) :: zhip, zhi, zhs, zvs_excess, zfra 1554 REAL(wp), DIMENSION(jpi,jpj) :: zswitch 1555 !!------------------------------------------------------------------- 1556 ! 1557 ! 1558 DO jl = 1, jpl 1559 1560 DO jj = 1, jpj 1561 DO ji = 1, jpi 1562 IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 1563 ! 1564 ! ! -- check h_ip -- ! 1565 ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip 1566 IF( ln_pnd_H12 .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN 1567 zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) ) 1568 IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN 1569 pa_ip(ji,jj,jl) = pv_ip(ji,jj,jl) / phip_max(ji,jj,jl) 1570 ENDIF 1571 ENDIF 1572 ! 1573 ! ! -- check h_i -- ! 1574 ! if h_i is larger than the surrounding 9 pts => reduce h_i and increase a_i 1575 zhi = pv_i(ji,jj,jl) / pa_i(ji,jj,jl) 1576 IF( zhi > phi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 1577 pa_i(ji,jj,jl) = pv_i(ji,jj,jl) / MIN( phi_max(ji,jj,jl), hi_max(jpl) ) !-- bound h_i to hi_max (99 m) 1578 ENDIF 1579 ! 1580 ! ! -- check h_s -- ! 1581 ! if h_s is larger than the surrounding 9 pts => put the snow excess in the ocean 1582 zhs = pv_s(ji,jj,jl) / pa_i(ji,jj,jl) 1583 IF( pv_s(ji,jj,jl) > 0._wp .AND. zhs > phs_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 1584 zfra = phs_max(ji,jj,jl) / MAX( zhs, epsi20 ) 1585 ! 1586 wfx_res(ji,jj) = wfx_res(ji,jj) + ( pv_s(ji,jj,jl) - pa_i(ji,jj,jl) * phs_max(ji,jj,jl) ) * rhos * r1_rdtice 1587 hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * r1_rdtice ! W.m-2 <0 1588 ! 1589 pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra 1590 pv_s(ji,jj,jl) = pa_i(ji,jj,jl) * phs_max(ji,jj,jl) 1591 ENDIF 1592 ! 1593 ! ! -- check snow load -- ! 1594 ! if snow load makes snow-ice interface to deplet below the ocean surface => put the snow excess in the ocean 1595 ! this correction is crucial because of the call to routine icecor afterwards which imposes a mini of ice thick. (rn_himin) 1596 ! this imposed mini can artificially make the snow very thick (if concentration decreases drastically) 1597 zvs_excess = MAX( 0._wp, pv_s(ji,jj,jl) - pv_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos ) 1598 IF( zvs_excess > 0._wp ) THEN 1599 zfra = ( pv_s(ji,jj,jl) - zvs_excess ) / MAX( pv_s(ji,jj,jl), epsi20 ) 1600 wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * r1_rdtice 1601 hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * r1_rdtice ! W.m-2 <0 1602 ! 1603 pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra 1604 pv_s(ji,jj,jl) = pv_s(ji,jj,jl) - zvs_excess 1605 ENDIF 1606 1607 ENDIF 1608 END DO 1609 END DO 1610 END DO 1611 ! !-- correct pond fraction to avoid a_ip > a_i 1612 WHERE( pa_ip(:,:,:) > pa_i(:,:,:) ) pa_ip(:,:,:) = pa_i(:,:,:) 1613 ! 1614 ! 1615 END SUBROUTINE Hbig 1616 1360 1617 #else 1361 1618 !!---------------------------------------------------------------------- -
NEMO/releases/release-4.0/src/ICE/icedyn_rhg.F90
r10413 r10910 58 58 !!-------------------------------------------------------------------- 59 59 INTEGER, INTENT(in) :: kt ! ice time step 60 ! 61 INTEGER :: jl ! dummy loop indices 60 62 !!-------------------------------------------------------------------- 61 63 ! controls … … 68 70 WRITE(numout,*)'~~~~~~~~~~~' 69 71 ENDIF 70 71 ! -------- 72 ! Rheology 73 ! -------- 72 ! 73 IF( ln_landfast_home ) THEN !-- Landfast ice parameterization 74 tau_icebfr(:,:) = 0._wp 75 DO jl = 1, jpl 76 WHERE( h_i(:,:,jl) > ht_n(:,:) * rn_depfra ) tau_icebfr(:,:) = tau_icebfr(:,:) + a_i(:,:,jl) * rn_icebfr 77 END DO 78 ENDIF 79 ! 80 !--------------! 81 !== Rheology ==! 82 !--------------! 74 83 SELECT CASE( nice_rhg ) 75 84 ! !------------------------! -
NEMO/releases/release-4.0/src/ICE/icewri.F90
r10785 r10910 145 145 146 146 ! record presence of fast ice 147 WHERE( z2d(:,:) < 5.e-04_wp .AND. zmsk 00(:,:) == 1._wp ) ; zfast(:,:) = 1._wp147 WHERE( z2d(:,:) < 5.e-04_wp .AND. zmsk15(:,:) == 1._wp ) ; zfast(:,:) = 1._wp 148 148 ELSEWHERE ; zfast(:,:) = 0._wp 149 149 END WHERE
Note: See TracChangeset
for help on using the changeset viewer.
