- Timestamp:
- 2019-07-18T17:12:52+02:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/UKMO/NEMO_4.0_GC_couple_pkg/src/ICE/icedyn.F90
r10888 r11303 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 CALL ice_var_zapsmall ! -- zap small areas 121 ! 133 122 CASE ( np_dynADV1D ) !== pure advection ==! (1D) 134 ALLOCATE( zdivu_i(jpi,jpj) )123 ! 135 124 ! --- monotonicity test from Schar & Smolarkiewicz 1996 --- ! 136 125 ! CFL = 0.5 at a distance from the bound of 1/6 of the basin length … … 145 134 END DO 146 135 ! --- 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 136 CALL ice_dyn_adv ( kt ) ! -- advection of ice 137 ! 161 138 CASE ( np_dynADV2D ) !== pure advection ==! (2D w prescribed velocities) 162 ALLOCATE( zdivu_i(jpi,jpj) )139 ! 163 140 u_ice(:,:) = rn_uice * umask(:,:,1) 164 141 v_ice(:,:) = rn_vice * vmask(:,:,1) … … 166 143 !CALL RANDOM_NUMBER(v_ice(:,:)) ; v_ice(:,:) = v_ice(:,:) * 0.1 + rn_vice * 0.9 * vmask(:,:,1) 167 144 ! --- 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) 145 CALL ice_dyn_adv ( kt ) ! -- advection of ice 146 147 END SELECT 148 ! 149 ! 150 ! diagnostics: divergence at T points 151 IF( iom_use('icediv') ) THEN 152 ! 153 SELECT CASE( nice_dyn ) 154 155 CASE ( np_dynADV1D , np_dynADV2D ) 156 157 ALLOCATE( zdivu_i(jpi,jpj) ) 158 DO jj = 2, jpjm1 159 DO ji = 2, jpim1 160 zdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) & 161 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj) 162 END DO 175 163 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 SELECT164 CALL lbc_lnk( 'icedyn', zdivu_i, 'T', 1. ) 165 CALL iom_put( "icediv" , zdivu_i(:,:) ) 166 DEALLOCATE( zdivu_i ) 167 168 END SELECT 169 ! 170 ENDIF 183 171 ! 184 172 ! controls … … 188 176 189 177 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 277 278 178 SUBROUTINE Hpiling 279 179 !!------------------------------------------------------------------- … … 290 190 ! controls 291 191 IF( ln_icediachk ) CALL ice_cons_hsm(0, 'Hpiling', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 292 !293 CALL ice_var_zapsmall !-- zap small areas294 192 ! 295 193 at_i(:,:) = SUM( a_i(:,:,:), dim=3 ) … … 337 235 WRITE(numout,*) '~~~~~~~~~~~~' 338 236 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 237 WRITE(numout,*) ' Full ice dynamics (rhg + adv + ridge/raft + corr) ln_dynALL = ', ln_dynALL 238 WRITE(numout,*) ' No ridge/raft & No cor (rhg + adv) ln_dynRHGADV = ', ln_dynRHGADV 239 WRITE(numout,*) ' Advection 1D only (Schar & Smolarkiewicz 1996) ln_dynADV1D = ', ln_dynADV1D 240 WRITE(numout,*) ' Advection 2D only (rn_uvice + adv) ln_dynADV2D = ', ln_dynADV2D 241 WRITE(numout,*) ' with prescribed velocity given by (u,v)_ice = (rn_uice,rn_vice) = (', rn_uice,',',rn_vice,')' 242 WRITE(numout,*) ' lateral boundary condition for sea ice dynamics rn_ishlat = ', rn_ishlat 243 WRITE(numout,*) ' Landfast: param from Lemieux 2016 ln_landfast_L16 = ', ln_landfast_L16 244 WRITE(numout,*) ' Landfast: param from home made ln_landfast_home= ', ln_landfast_home 245 WRITE(numout,*) ' fraction of ocean depth that ice must reach rn_depfra = ', rn_depfra 246 WRITE(numout,*) ' maximum bottom stress per unit area of contact rn_icebfr = ', rn_icebfr 247 WRITE(numout,*) ' relax time scale (s-1) to reach static friction rn_lfrelax = ', rn_lfrelax 248 WRITE(numout,*) ' isotropic tensile strength rn_tensile = ', rn_tensile 351 249 WRITE(numout,*) 352 250 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
