Changeset 11822 for NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/ICE/icedyn.F90
- Timestamp:
- 2019-10-29T11:41:36+01:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/ICE/icedyn.F90
r11480 r11822 75 75 INTEGER, INTENT(in) :: Kmm ! ocean time level index 76 76 !! 77 INTEGER :: ji, jj , jl! dummy loop indices77 INTEGER :: ji, jj ! dummy loop indices 78 78 REAL(wp) :: zcoefu, zcoefv 79 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zhi_max, zhs_max, zhip_max 80 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdivu_i 79 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdivu_i 81 80 !!-------------------------------------------------------------------- 82 81 ! … … 90 89 ENDIF 91 90 ! 92 IF( ln_landfast_home ) THEN !-- Landfast ice parameterization 93 tau_icebfr(:,:) = 0._wp 94 DO jl = 1, jpl 95 WHERE( h_i_b(:,:,jl) > ht(:,:) * rn_depfra ) tau_icebfr(:,:) = tau_icebfr(:,:) + a_i(:,:,jl) * rn_icebfr 96 END DO 97 ENDIF 98 ! 99 ! !-- Record max of the surrounding 9-pts ice thick. (for CALL Hbig) 100 DO jl = 1, jpl 101 DO jj = 2, jpjm1 102 DO ji = fs_2, fs_jpim1 103 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), & 104 & h_ip_b(ji-1,jj ,jl), h_ip_b(ji ,jj-1,jl), & 105 & h_ip_b(ji+1,jj+1,jl), h_ip_b(ji-1,jj-1,jl), & 106 & h_ip_b(ji+1,jj-1,jl), h_ip_b(ji-1,jj+1,jl) ) 107 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), & 108 & h_i_b (ji-1,jj ,jl), h_i_b (ji ,jj-1,jl), & 109 & h_i_b (ji+1,jj+1,jl), h_i_b (ji-1,jj-1,jl), & 110 & h_i_b (ji+1,jj-1,jl), h_i_b (ji-1,jj+1,jl) ) 111 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), & 112 & h_s_b (ji-1,jj ,jl), h_s_b (ji ,jj-1,jl), & 113 & h_s_b (ji+1,jj+1,jl), h_s_b (ji-1,jj-1,jl), & 114 & h_s_b (ji+1,jj-1,jl), h_s_b (ji-1,jj+1,jl) ) 115 END DO 116 END DO 117 END DO 118 CALL lbc_lnk_multi( 'icedyn', zhi_max(:,:,:), 'T', 1., zhs_max(:,:,:), 'T', 1., zhip_max(:,:,:), 'T', 1. ) 119 ! 120 ! 121 SELECT CASE( nice_dyn ) !-- Set which dynamics is running 91 ! retrieve thickness from volume for landfast param. and UMx advection scheme 92 WHERE( a_i(:,:,:) >= epsi20 ) 93 h_i(:,:,:) = v_i(:,:,:) / a_i_b(:,:,:) 94 h_s(:,:,:) = v_s(:,:,:) / a_i_b(:,:,:) 95 ELSEWHERE 96 h_i(:,:,:) = 0._wp 97 h_s(:,:,:) = 0._wp 98 END WHERE 99 ! 100 WHERE( a_ip(:,:,:) >= epsi20 ) 101 h_ip(:,:,:) = v_ip(:,:,:) / a_ip(:,:,:) 102 ELSEWHERE 103 h_ip(:,:,:) = 0._wp 104 END WHERE 105 ! 106 ! 107 SELECT CASE( nice_dyn ) !-- Set which dynamics is running 122 108 123 109 CASE ( np_dynALL ) !== all dynamical processes ==! 124 CALL ice_dyn_rhg ( kt, Kmm ) ! -- rheology 125 CALL ice_dyn_adv ( kt ) ; CALL Hbig( zhi_max, zhs_max, zhip_max ) ! -- advection of ice + correction on ice thickness 126 CALL ice_dyn_rdgrft( kt ) ! -- ridging/rafting 127 CALL ice_cor ( kt , 1 ) ! -- Corrections 128 110 ! 111 CALL ice_dyn_rhg ( kt, Kmm ) ! -- rheology 112 CALL ice_dyn_adv ( kt ) ! -- advection of ice 113 CALL ice_dyn_rdgrft( kt ) ! -- ridging/rafting 114 CALL ice_cor ( kt , 1 ) ! -- Corrections 115 ! 129 116 CASE ( np_dynRHGADV ) !== no ridge/raft & no corrections ==! 130 CALL ice_dyn_rhg ( kt, Kmm ) ! -- rheology 131 CALL ice_dyn_adv ( kt ) ; CALL Hbig( zhi_max, zhs_max, zhip_max ) ! -- advection of ice + correction on ice thickness 132 CALL Hpiling ! -- simple pile-up (replaces ridging/rafting) 133 117 ! 118 CALL ice_dyn_rhg ( kt, Kmm ) ! -- rheology 119 CALL ice_dyn_adv ( kt ) ! -- advection of ice 120 CALL Hpiling ! -- simple pile-up (replaces ridging/rafting) 121 CALL ice_var_zapsmall ! -- zap small areas 122 ! 134 123 CASE ( np_dynADV1D ) !== pure advection ==! (1D) 135 ALLOCATE( zdivu_i(jpi,jpj) )124 ! 136 125 ! --- monotonicity test from Schar & Smolarkiewicz 1996 --- ! 137 126 ! CFL = 0.5 at a distance from the bound of 1/6 of the basin length … … 146 135 END DO 147 136 ! --- 148 CALL ice_dyn_adv ( kt ) ! -- advection of ice 149 150 ! diagnostics: divergence at T points 151 DO jj = 2, jpjm1 152 DO ji = 2, jpim1 153 zdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) & 154 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj) 155 END DO 156 END DO 157 CALL lbc_lnk( 'icedyn', zdivu_i, 'T', 1. ) 158 IF( iom_use('icediv') ) CALL iom_put( "icediv" , zdivu_i(:,:) ) 159 160 DEALLOCATE( zdivu_i ) 161 137 CALL ice_dyn_adv ( kt ) ! -- advection of ice 138 ! 162 139 CASE ( np_dynADV2D ) !== pure advection ==! (2D w prescribed velocities) 163 ALLOCATE( zdivu_i(jpi,jpj) )140 ! 164 141 u_ice(:,:) = rn_uice * umask(:,:,1) 165 142 v_ice(:,:) = rn_vice * vmask(:,:,1) … … 167 144 !CALL RANDOM_NUMBER(v_ice(:,:)) ; v_ice(:,:) = v_ice(:,:) * 0.1 + rn_vice * 0.9 * vmask(:,:,1) 168 145 ! --- 169 CALL ice_dyn_adv ( kt ) ! -- advection of ice 170 171 ! diagnostics: divergence at T points 172 DO jj = 2, jpjm1 173 DO ji = 2, jpim1 174 zdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) & 175 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj) 146 CALL ice_dyn_adv ( kt ) ! -- advection of ice 147 148 END SELECT 149 ! 150 ! 151 ! diagnostics: divergence at T points 152 IF( iom_use('icediv') ) THEN 153 ! 154 SELECT CASE( nice_dyn ) 155 156 CASE ( np_dynADV1D , np_dynADV2D ) 157 158 ALLOCATE( zdivu_i(jpi,jpj) ) 159 DO jj = 2, jpjm1 160 DO ji = 2, jpim1 161 zdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) & 162 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj) 163 END DO 176 164 END DO 177 END DO 178 CALL lbc_lnk( 'icedyn', zdivu_i, 'T', 1. ) 179 IF( iom_use('icediv') ) CALL iom_put( "icediv" , zdivu_i(:,:) ) 180 181 DEALLOCATE( zdivu_i ) 182 183 END SELECT 165 CALL lbc_lnk( 'icedyn', zdivu_i, 'T', 1. ) 166 ! output 167 CALL iom_put( 'icediv' , zdivu_i ) 168 169 DEALLOCATE( zdivu_i ) 170 171 END SELECT 172 ! 173 ENDIF 184 174 ! 185 175 ! controls … … 189 179 190 180 191 SUBROUTINE Hbig( phi_max, phs_max, phip_max )192 !!-------------------------------------------------------------------193 !! *** ROUTINE Hbig ***194 !!195 !! ** Purpose : Thickness correction in case advection scheme creates196 !! abnormally tick ice or snow197 !!198 !! ** Method : 1- check whether ice thickness is larger than the surrounding 9-points199 !! (before advection) and reduce it by adapting ice concentration200 !! 2- check whether snow thickness is larger than the surrounding 9-points201 !! (before advection) and reduce it by sending the excess in the ocean202 !! 3- check whether snow load deplets the snow-ice interface below sea level$203 !! and reduce it by sending the excess in the ocean204 !! 4- correct pond fraction to avoid a_ip > a_i205 !!206 !! ** input : Max thickness of the surrounding 9-points207 !!-------------------------------------------------------------------208 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: phi_max, phs_max, phip_max ! max ice thick from surrounding 9-pts209 !210 INTEGER :: ji, jj, jl ! dummy loop indices211 REAL(wp) :: zhip, zhi, zhs, zvs_excess, zfra212 !!-------------------------------------------------------------------213 ! controls214 IF( ln_icediachk ) CALL ice_cons_hsm(0, 'Hbig', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation215 !216 CALL ice_var_zapsmall !-- zap small areas217 !218 DO jl = 1, jpl219 DO jj = 1, jpj220 DO ji = 1, jpi221 IF ( v_i(ji,jj,jl) > 0._wp ) THEN222 !223 ! ! -- check h_ip -- !224 ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip225 IF( ln_pnd_H12 .AND. v_ip(ji,jj,jl) > 0._wp ) THEN226 zhip = v_ip(ji,jj,jl) / MAX( epsi20, a_ip(ji,jj,jl) )227 IF( zhip > phip_max(ji,jj,jl) .AND. a_ip(ji,jj,jl) < 0.15 ) THEN228 a_ip(ji,jj,jl) = v_ip(ji,jj,jl) / phip_max(ji,jj,jl)229 ENDIF230 ENDIF231 !232 ! ! -- check h_i -- !233 ! if h_i is larger than the surrounding 9 pts => reduce h_i and increase a_i234 zhi = v_i(ji,jj,jl) / a_i(ji,jj,jl)235 IF( zhi > phi_max(ji,jj,jl) .AND. a_i(ji,jj,jl) < 0.15 ) THEN236 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)237 ENDIF238 !239 ! ! -- check h_s -- !240 ! if h_s is larger than the surrounding 9 pts => put the snow excess in the ocean241 zhs = v_s(ji,jj,jl) / a_i(ji,jj,jl)242 IF( v_s(ji,jj,jl) > 0._wp .AND. zhs > phs_max(ji,jj,jl) .AND. a_i(ji,jj,jl) < 0.15 ) THEN243 zfra = phs_max(ji,jj,jl) / MAX( zhs, epsi20 )244 !245 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_rdtice246 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 <0247 !248 e_s(ji,jj,1:nlay_s,jl) = e_s(ji,jj,1:nlay_s,jl) * zfra249 v_s(ji,jj,jl) = a_i(ji,jj,jl) * phs_max(ji,jj,jl)250 ENDIF251 !252 ! ! -- check snow load -- !253 ! if snow load makes snow-ice interface to deplet below the ocean surface => put the snow excess in the ocean254 ! this correction is crucial because of the call to routine icecor afterwards which imposes a mini of ice thick. (rn_himin)255 ! this imposed mini can artificially make the snow very thick (if concentration decreases drastically)256 zvs_excess = MAX( 0._wp, v_s(ji,jj,jl) - v_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos )257 IF( zvs_excess > 0._wp ) THEN258 zfra = ( v_s(ji,jj,jl) - zvs_excess ) / MAX( v_s(ji,jj,jl), epsi20 )259 wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * r1_rdtice260 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 <0261 !262 e_s(ji,jj,1:nlay_s,jl) = e_s(ji,jj,1:nlay_s,jl) * zfra263 v_s(ji,jj,jl) = v_s(ji,jj,jl) - zvs_excess264 ENDIF265 266 ENDIF267 END DO268 END DO269 END DO270 ! !-- correct pond fraction to avoid a_ip > a_i271 WHERE( a_ip(:,:,:) > a_i(:,:,:) ) a_ip(:,:,:) = a_i(:,:,:)272 !273 ! controls274 IF( ln_icediachk ) CALL ice_cons_hsm(1, 'Hbig', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation275 !276 END SUBROUTINE Hbig277 278 279 181 SUBROUTINE Hpiling 280 182 !!------------------------------------------------------------------- … … 291 193 ! controls 292 194 IF( ln_icediachk ) CALL ice_cons_hsm(0, 'Hpiling', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 293 !294 CALL ice_var_zapsmall !-- zap small areas295 195 ! 296 196 at_i(:,:) = SUM( a_i(:,:,:), dim=3 ) … … 322 222 NAMELIST/namdyn/ ln_dynALL, ln_dynRHGADV, ln_dynADV1D, ln_dynADV2D, rn_uice, rn_vice, & 323 223 & rn_ishlat , & 324 & ln_landfast_L16, ln_landfast_home,rn_depfra, rn_icebfr, rn_lfrelax, rn_tensile224 & ln_landfast_L16, rn_depfra, rn_icebfr, rn_lfrelax, rn_tensile 325 225 !!------------------------------------------------------------------- 326 226 ! 327 227 REWIND( numnam_ice_ref ) ! Namelist namdyn in reference namelist : Ice dynamics 328 228 READ ( numnam_ice_ref, namdyn, IOSTAT = ios, ERR = 901) 329 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn in reference namelist' , lwp)229 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn in reference namelist' ) 330 230 REWIND( numnam_ice_cfg ) ! Namelist namdyn in configuration namelist : Ice dynamics 331 231 READ ( numnam_ice_cfg, namdyn, IOSTAT = ios, ERR = 902 ) 332 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namdyn in configuration namelist' , lwp)232 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namdyn in configuration namelist' ) 333 233 IF(lwm) WRITE( numoni, namdyn ) 334 234 ! … … 338 238 WRITE(numout,*) '~~~~~~~~~~~~' 339 239 WRITE(numout,*) ' Namelist namdyn:' 340 WRITE(numout,*) ' Full ice dynamics (rhg + adv + ridge/raft + corr) ln_dynALL = ', ln_dynALL 341 WRITE(numout,*) ' No ridge/raft & No cor (rhg + adv) ln_dynRHGADV = ', ln_dynRHGADV 342 WRITE(numout,*) ' Advection 1D only (Schar & Smolarkiewicz 1996) ln_dynADV1D = ', ln_dynADV1D 343 WRITE(numout,*) ' Advection 2D only (rn_uvice + adv) ln_dynADV2D = ', ln_dynADV2D 344 WRITE(numout,*) ' with prescribed velocity given by (u,v)_ice = (rn_uice,rn_vice) = (', rn_uice,',', rn_vice,')' 345 WRITE(numout,*) ' lateral boundary condition for sea ice dynamics rn_ishlat = ', rn_ishlat 346 WRITE(numout,*) ' Landfast: param from Lemieux 2016 ln_landfast_L16 = ', ln_landfast_L16 347 WRITE(numout,*) ' Landfast: param from home made ln_landfast_home= ', ln_landfast_home 348 WRITE(numout,*) ' fraction of ocean depth that ice must reach rn_depfra = ', rn_depfra 349 WRITE(numout,*) ' maximum bottom stress per unit area of contact rn_icebfr = ', rn_icebfr 350 WRITE(numout,*) ' relax time scale (s-1) to reach static friction rn_lfrelax = ', rn_lfrelax 351 WRITE(numout,*) ' isotropic tensile strength rn_tensile = ', rn_tensile 240 WRITE(numout,*) ' Full ice dynamics (rhg + adv + ridge/raft + corr) ln_dynALL = ', ln_dynALL 241 WRITE(numout,*) ' No ridge/raft & No cor (rhg + adv) ln_dynRHGADV = ', ln_dynRHGADV 242 WRITE(numout,*) ' Advection 1D only (Schar & Smolarkiewicz 1996) ln_dynADV1D = ', ln_dynADV1D 243 WRITE(numout,*) ' Advection 2D only (rn_uvice + adv) ln_dynADV2D = ', ln_dynADV2D 244 WRITE(numout,*) ' with prescribed velocity given by (u,v)_ice = (rn_uice,rn_vice) = (', rn_uice,',',rn_vice,')' 245 WRITE(numout,*) ' lateral boundary condition for sea ice dynamics rn_ishlat = ', rn_ishlat 246 WRITE(numout,*) ' Landfast: param from Lemieux 2016 ln_landfast_L16 = ', ln_landfast_L16 247 WRITE(numout,*) ' fraction of ocean depth that ice must reach rn_depfra = ', rn_depfra 248 WRITE(numout,*) ' maximum bottom stress per unit area of contact rn_icebfr = ', rn_icebfr 249 WRITE(numout,*) ' relax time scale (s-1) to reach static friction rn_lfrelax = ', rn_lfrelax 250 WRITE(numout,*) ' isotropic tensile strength rn_tensile = ', rn_tensile 352 251 WRITE(numout,*) 353 252 ENDIF … … 372 271 ENDIF 373 272 ! !--- Landfast ice 374 IF( .NOT.ln_landfast_L16 .AND. .NOT.ln_landfast_home ) tau_icebfr(:,:) = 0._wp 375 ! 376 IF ( ln_landfast_L16 .AND. ln_landfast_home ) THEN 377 CALL ctl_stop( 'ice_dyn_init: choose one and only one landfast parameterization (ln_landfast_L16 or ln_landfast_home)' ) 378 ENDIF 273 IF( .NOT.ln_landfast_L16 ) tau_icebfr(:,:) = 0._wp 379 274 ! 380 275 CALL ice_dyn_rdgrft_init ! set ice ridging/rafting parameters
Note: See TracChangeset
for help on using the changeset viewer.
