Changeset 9604 for NEMO/trunk/src/ICE/icethd_do.F90
- Timestamp:
- 2018-05-18T09:53:22+02:00 (6 years ago)
- File:
-
- 1 edited
-
NEMO/trunk/src/ICE/icethd_do.F90 (modified) (23 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/ICE/icethd_do.F90
r9598 r9604 4 4 !! sea-ice: sea ice growth in the leads (open water) 5 5 !!====================================================================== 6 !! History : LIM ! 2005-12 (M. Vancoppenolle) Original code 7 !! - ! 2006-01 (M. Vancoppenolle) add ITD 8 !! 3.0 ! 2007-07 (M. Vancoppenolle) Mass and energy conservation tested 9 !! 4.0 ! 2011-02 (G. Madec) dynamical allocation 6 !! History : ! 2005-12 (M. Vancoppenolle) Original code 7 !! 4.0 ! 2018 (many people) SI3 [aka Sea Ice cube] 10 8 !!---------------------------------------------------------------------- 11 9 #if defined key_si3 … … 58 56 !! 59 57 !! ** Purpose : Computation of the evolution of the ice thickness and 60 !! concentration as a function of the heat balance in the leads.58 !! concentration as a function of the heat balance in the leads 61 59 !! 62 !! ** Method : Ice is formed in the open water when ocean lose heat 63 !! (heat budget of open water Bl is negative) . 64 !! Computation of the increase of 1-A (ice concentration) fol- 65 !! lowing the law : 66 !! (dA/dt)acc = F[ (1-A)/(1-a) ] * [ Bl / (Li*h0) ] 67 !! where - h0 is the thickness of ice created in the lead 68 !! - a is a minimum fraction for leads 69 !! - F is a monotonic non-increasing function defined as: 60 !! ** Method : Ice is formed in the open water when ocean looses heat 61 !! (heat budget of open water is negative) following 62 !! 63 !! (dA/dt)acc = F[ (1-A)/(1-a) ] * [ Bl / (Li*h0) ] 64 !! where - h0 is the thickness of ice created in the lead 65 !! - a is a minimum fraction for leads 66 !! - F is a monotonic non-increasing function defined as: 70 67 !! F(X)=( 1 - X**exld )**(1.0/exld) 71 !! - exld is the exponent closure rate (=2 default val.)68 !! - exld is the exponent closure rate (=2 default val.) 72 69 !! 73 70 !! ** Action : - Adjustment of snow and ice thicknesses and heat … … 79 76 !!------------------------------------------------------------------------ 80 77 INTEGER :: ji, jj, jk, jl ! dummy loop indices 81 INTEGER :: iter ! - -82 REAL(wp) :: ztmelts, zfrazb, zweight, zde ! local scalars78 INTEGER :: iter ! - - 79 REAL(wp) :: ztmelts, zfrazb, zweight, zde ! local scalars 83 80 REAL(wp) :: zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf ! - - 84 81 REAL(wp) :: ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, zsqcd , zhicrit ! - - … … 105 102 REAL(wp), DIMENSION(jpij) :: zvrel_1d ! relative ice / frazil velocity (1D vector) 106 103 ! 107 REAL(wp), DIMENSION(jpij,jpl) :: zv_b 108 REAL(wp), DIMENSION(jpij,jpl) :: za_b 104 REAL(wp), DIMENSION(jpij,jpl) :: zv_b ! old volume of ice in category jl 105 REAL(wp), DIMENSION(jpij,jpl) :: za_b ! old area of ice in category jl 109 106 ! 110 107 REAL(wp), DIMENSION(jpij,nlay_i,jpl) :: ze_i_2d !: 1-D version of e_i 111 108 ! 112 REAL(wp), DIMENSION(jpi,jpj) :: zvrel 113 ! 114 REAL(wp) :: zcai = 1.4e-3_wp 109 REAL(wp), DIMENSION(jpi,jpj) :: zvrel ! relative ice / frazil velocity 110 ! 111 REAL(wp) :: zcai = 1.4e-3_wp ! ice-air drag (clem: should be dependent on coupling/forcing used) 115 112 !!-----------------------------------------------------------------------! 116 113 … … 121 118 122 119 !------------------------------------------------------------------------------! 123 ! 3) Collection thickness of ice formed in leads and polynyas120 ! 1) Collection thickness of ice formed in leads and polynyas 124 121 !------------------------------------------------------------------------------! 125 122 ! ht_i_new is the thickness of new ice formed in open water 126 ! ht_i_new can be either prescribed ( frazswi = 0) or computed (frazswi = 1)123 ! ht_i_new can be either prescribed (ln_frazil=F) or computed (ln_frazil=T) 127 124 ! Frazil ice forms in open water, is transported by wind 128 125 ! accumulates at the edge of the consolidated ice edge … … 130 127 ! collection thickness. 131 128 132 ! Note : the following algorithm currently breaks vectorization133 !134 135 129 zvrel(:,:) = 0._wp 136 130 … … 142 136 IF( ln_frazil ) THEN 143 137 ! 144 !--------------------145 ! Physical constants146 !--------------------147 138 ht_i_new(:,:) = 0._wp 148 139 ! 149 zhicrit = 0.04 ! frazil ice thickness 140 ! Physical constants 141 zhicrit = 0.04 ! frazil ice thickness 150 142 ztwogp = 2. * rau0 / ( grav * 0.3 * ( rau0 - rhoic ) ) ! reduced grav 151 zsqcd = 1.0 / SQRT( 1.3 * zcai ) ! 1/SQRT(airdensity*drag)143 zsqcd = 1.0 / SQRT( 1.3 * zcai ) ! 1/SQRT(airdensity*drag) 152 144 zgamafr = 0.03 153 145 ! … … 155 147 DO ji = 2, jpim1 156 148 IF ( qlead(ji,jj) < 0._wp .AND. tau_icebfr(ji,jj) == 0._wp ) THEN ! activated if cooling and no landfast 157 !------------- 158 ! Wind stress 159 !------------- 160 ! C-grid wind stress components 149 ! -- Wind stress -- ! 161 150 ztaux = ( utau_ice(ji-1,jj ) * umask(ji-1,jj ,1) & 162 151 & + utau_ice(ji ,jj ) * umask(ji ,jj ,1) ) * 0.5_wp … … 166 155 ztenagm = SQRT( SQRT( ztaux * ztaux + ztauy * ztauy ) ) 167 156 168 !--------------------- 169 ! Frazil ice velocity 170 !--------------------- 157 ! -- Frazil ice velocity -- ! 171 158 rswitch = MAX( 0._wp, SIGN( 1._wp , ztenagm - epsi10 ) ) 172 159 zvfrx = rswitch * zgamafr * zsqcd * ztaux / MAX( ztenagm, epsi10 ) 173 160 zvfry = rswitch * zgamafr * zsqcd * ztauy / MAX( ztenagm, epsi10 ) 174 161 175 !------------------- 176 ! Pack ice velocity 177 !------------------- 178 ! C-grid ice velocity 162 ! -- Pack ice velocity -- ! 179 163 zvgx = ( u_ice(ji-1,jj ) * umask(ji-1,jj ,1) + u_ice(ji,jj) * umask(ji,jj,1) ) * 0.5_wp 180 164 zvgy = ( v_ice(ji ,jj-1) * vmask(ji ,jj-1,1) + v_ice(ji,jj) * vmask(ji,jj,1) ) * 0.5_wp 181 165 182 !----------------------------------- 183 ! Relative frazil/pack ice velocity 184 !----------------------------------- 185 ! absolute relative velocity 166 ! -- Relative frazil/pack ice velocity -- ! 186 167 rswitch = MAX( 0._wp, SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) 187 168 zvrel2 = MAX( ( zvfrx - zvgx ) * ( zvfrx - zvgx ) & … … 189 170 zvrel(ji,jj) = SQRT( zvrel2 ) 190 171 191 !--------------------- 192 ! Iterative procedure 193 !--------------------- 172 ! -- new ice thickness (iterative loop) -- ! 194 173 ht_i_new(ji,jj) = zhicrit + ( zhicrit + 0.1 ) & 195 174 & / ( ( zhicrit + 0.1 ) * ( zhicrit + 0.1 ) - zhicrit * zhicrit ) * ztwogp * zvrel2 … … 205 184 END DO 206 185 ! 207 ENDIF ! end of selection of pixels where ice forms186 ENDIF 208 187 ! 209 188 END DO … … 212 191 CALL lbc_lnk_multi( zvrel, 'T', 1., ht_i_new, 'T', 1. ) 213 192 214 ENDIF ! End of computation of frazil ice collection thickness193 ENDIF 215 194 216 195 !------------------------------------------------------------------------------! 217 ! 4) Identify grid points where new ice forms196 ! 2) Compute thickness, salinity, enthalpy, age, area and volume of new ice 218 197 !------------------------------------------------------------------------------! 219 220 !-------------------------------------221 ! Select points for new ice formation222 !-------------------------------------223 198 ! This occurs if open water energy budget is negative (cooling) and there is no landfast ice 199 200 ! Identify grid points where new ice forms 224 201 npti = 0 ; nptidx(:) = 0 225 202 DO jj = 1, jpj … … 232 209 END DO 233 210 234 !------------------------------235 211 ! Move from 2-D to 1-D vectors 236 !------------------------------237 ! If ocean gains heat do nothing. Otherwise compute new ice formation238 239 212 IF ( npti > 0 ) THEN 240 213 … … 248 221 END DO 249 222 END DO 250 CALL tab_2d_1d( npti, nptidx(1:npti), qlead_1d (1:npti) , qlead ) 251 CALL tab_2d_1d( npti, nptidx(1:npti), t_bo_1d (1:npti) , t_bo ) 252 CALL tab_2d_1d( npti, nptidx(1:npti), sfx_opw_1d(1:npti) , sfx_opw ) 253 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_opw_1d(1:npti) , wfx_opw ) 254 CALL tab_2d_1d( npti, nptidx(1:npti), zh_newice (1:npti) , ht_i_new ) 255 CALL tab_2d_1d( npti, nptidx(1:npti), zvrel_1d (1:npti) , zvrel ) 256 257 CALL tab_2d_1d( npti, nptidx(1:npti), hfx_thd_1d(1:npti) , hfx_thd ) 258 CALL tab_2d_1d( npti, nptidx(1:npti), hfx_opw_1d(1:npti) , hfx_opw ) 259 CALL tab_2d_1d( npti, nptidx(1:npti), rn_amax_1d(1:npti) , rn_amax_2d ) 260 CALL tab_2d_1d( npti, nptidx(1:npti), sss_1d (1:npti) , sss_m ) 261 262 !------------------------------------------------------------------------------| 263 ! 2) Convert units for ice internal energy 264 !------------------------------------------------------------------------------| 223 CALL tab_2d_1d( npti, nptidx(1:npti), qlead_1d (1:npti) , qlead ) 224 CALL tab_2d_1d( npti, nptidx(1:npti), t_bo_1d (1:npti) , t_bo ) 225 CALL tab_2d_1d( npti, nptidx(1:npti), sfx_opw_1d(1:npti) , sfx_opw ) 226 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_opw_1d(1:npti) , wfx_opw ) 227 CALL tab_2d_1d( npti, nptidx(1:npti), zh_newice (1:npti) , ht_i_new ) 228 CALL tab_2d_1d( npti, nptidx(1:npti), zvrel_1d (1:npti) , zvrel ) 229 230 CALL tab_2d_1d( npti, nptidx(1:npti), hfx_thd_1d(1:npti) , hfx_thd ) 231 CALL tab_2d_1d( npti, nptidx(1:npti), hfx_opw_1d(1:npti) , hfx_opw ) 232 CALL tab_2d_1d( npti, nptidx(1:npti), rn_amax_1d(1:npti) , rn_amax_2d ) 233 CALL tab_2d_1d( npti, nptidx(1:npti), sss_1d (1:npti) , sss_m ) 234 235 ! Convert units for ice internal energy 265 236 DO jl = 1, jpl 266 237 DO jk = 1, nlay_i … … 272 243 END DO 273 244 END DO 274 !------------------------------------------------------------------------------! 275 ! 5) Compute thickness, salinity, enthalpy, age, area and volume of new ice 276 !------------------------------------------------------------------------------! 277 278 !----------------------------------------- 245 279 246 ! Keep old ice areas and volume in memory 280 !-----------------------------------------281 247 zv_b(1:npti,:) = v_i_2d(1:npti,:) 282 248 za_b(1:npti,:) = a_i_2d(1:npti,:) 283 249 284 !---------------------- 285 ! Salinity of new ice 286 !---------------------- 250 ! --- Salinity of new ice --- ! 287 251 SELECT CASE ( nn_icesal ) 288 252 CASE ( 1 ) ! Sice = constant … … 296 260 END SELECT 297 261 298 !------------------------- 299 ! Heat content of new ice 300 !------------------------- 262 ! --- Heat content of new ice --- ! 301 263 ! We assume that new ice is formed at the seawater freezing point 302 264 DO ji = 1, npti … … 307 269 END DO 308 270 309 !---------------- 310 ! Age of new ice 311 !---------------- 271 ! --- Age of new ice --- ! 312 272 zo_newice(1:npti) = 0._wp 313 273 314 !------------------- 315 ! Volume of new ice 316 !------------------- 274 ! --- Volume of new ice --- ! 317 275 DO ji = 1, npti 318 276 … … 351 309 END IF 352 310 353 !----------------- 354 ! Area of new ice 355 !----------------- 311 ! --- Area of new ice --- ! 356 312 DO ji = 1, npti 357 313 za_newice(ji) = zv_newice(ji) / zh_newice(ji) … … 359 315 360 316 !------------------------------------------------------------------------------! 361 ! 6) Redistribute new ice area and volume into ice categories !317 ! 3) Redistribute new ice area and volume into ice categories ! 362 318 !------------------------------------------------------------------------------! 363 319 364 !------------------------ 365 ! 6.1) lateral ice growth 366 !------------------------ 320 ! --- lateral ice growth --- ! 367 321 ! If lateral ice growth gives an ice concentration gt 1, then 368 322 ! we keep the excessive volume in memory and attribute it later to bottom accretion … … 407 361 END DO 408 362 409 !------------------------------------------------ 410 ! 6.2) bottom ice growth + ice enthalpy remapping 411 !------------------------------------------------ 363 ! --- bottom ice growth + ice enthalpy remapping --- ! 412 364 DO jl = 1, jpl 413 365 … … 436 388 END DO 437 389 438 !----------------- 439 ! Update salinity 440 !----------------- 390 ! --- Update salinity --- ! 441 391 DO jl = 1, jpl 442 392 DO ji = 1, npti … … 445 395 END DO 446 396 447 !------------------------------------------------------------------------------! 448 ! 8) Change units for e_i 449 !------------------------------------------------------------------------------! 397 ! Change units for e_i 450 398 DO jl = 1, jpl 451 399 DO jk = 1, nlay_i … … 453 401 END DO 454 402 END DO 455 !------------------------------------------------------------------------------! 456 ! 7) Change 2D vectors to 1D vectors 457 !------------------------------------------------------------------------------! 403 404 ! Move 2D vectors to 1D vectors 458 405 CALL tab_2d_3d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i (:,:,:) ) 459 406 CALL tab_2d_3d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:) )
Note: See TracChangeset
for help on using the changeset viewer.
