Changeset 4688 for trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
- Timestamp:
- 2014-06-25T01:39:59+02:00 (10 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90 (modified) (20 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
r4333 r4688 30 30 USE wrk_nemo ! work arrays 31 31 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 32 USE limthd_ent 32 33 33 34 IMPLICIT NONE … … 37 38 38 39 REAL(wp) :: epsi10 = 1.e-10_wp ! 39 REAL(wp) :: zzero = 0._wp ! 40 REAL(wp) :: zone = 1._wp ! 40 REAL(wp) :: epsi20 = 1.e-20_wp ! 41 41 42 42 !!---------------------------------------------------------------------- … … 76 76 INTEGER :: layer, nbpac ! local integers 77 77 INTEGER :: ii, ij, iter ! - - 78 REAL(wp) :: ztmelts, zdv, z qold, zfrazb, zweight, zalphai, zindb, zinda, zde ! local scalars78 REAL(wp) :: ztmelts, zdv, zfrazb, zweight, zindb, zinda, zde ! local scalars 79 79 REAL(wp) :: zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf , zhicol_new ! - - 80 80 REAL(wp) :: ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, zsqcd , zhicrit ! - - 81 81 LOGICAL :: iterate_frazil ! iterate frazil ice collection thickness 82 82 CHARACTER (len = 15) :: fieldid 83 ! 84 INTEGER , POINTER, DIMENSION(:) :: zcatac ! indexes of categories where new ice grows 83 84 REAL(wp) :: zQm ! enthalpy exchanged with the ocean (J/m2, >0 towards ocean) 85 REAL(wp) :: zEi ! sea ice specific enthalpy (J/kg) 86 REAL(wp) :: zEw ! seawater specific enthalpy (J/kg) 87 REAL(wp) :: zfmdt ! mass flux x time step (kg/m2, >0 towards ocean) 88 89 REAL(wp) :: zv_newfra 90 91 INTEGER , POINTER, DIMENSION(:) :: jcat ! indexes of categories where new ice grows 85 92 REAL(wp), POINTER, DIMENSION(:) :: zswinew ! switch for new ice or not 86 93 … … 93 100 REAL(wp), POINTER, DIMENSION(:) :: zdv_res ! residual volume in case of excessive heat budget 94 101 REAL(wp), POINTER, DIMENSION(:) :: zda_res ! residual area in case of excessive heat budget 95 REAL(wp), POINTER, DIMENSION(:) :: zat_i_ ac! total ice fraction102 REAL(wp), POINTER, DIMENSION(:) :: zat_i_1d ! total ice fraction 96 103 REAL(wp), POINTER, DIMENSION(:) :: zat_i_lev ! total ice fraction for level ice only (type 1) 97 REAL(wp), POINTER, DIMENSION(:) :: zdh_frazb ! accretion of frazil ice at the ice bottom 98 REAL(wp), POINTER, DIMENSION(:) :: zvrel_ac ! relative ice / frazil velocity (1D vector) 99 100 REAL(wp), POINTER, DIMENSION(:,:) :: zhice_old ! previous ice thickness 101 REAL(wp), POINTER, DIMENSION(:,:) :: zdummy ! dummy thickness of new ice 102 REAL(wp), POINTER, DIMENSION(:,:) :: zdhicbot ! thickness of new ice which is accreted vertically 104 REAL(wp), POINTER, DIMENSION(:) :: zv_frazb ! accretion of frazil ice at the ice bottom 105 REAL(wp), POINTER, DIMENSION(:) :: zvrel_1d ! relative ice / frazil velocity (1D vector) 106 103 107 REAL(wp), POINTER, DIMENSION(:,:) :: zv_old ! old volume of ice in category jl 104 108 REAL(wp), POINTER, DIMENSION(:,:) :: za_old ! old area of ice in category jl 105 REAL(wp), POINTER, DIMENSION(:,:) :: za_i_ac ! 1-D version of a_i 106 REAL(wp), POINTER, DIMENSION(:,:) :: zv_i_ac ! 1-D version of v_i 107 REAL(wp), POINTER, DIMENSION(:,:) :: zoa_i_ac ! 1-D version of oa_i 108 REAL(wp), POINTER, DIMENSION(:,:) :: zsmv_i_ac ! 1-D version of smv_i 109 110 REAL(wp), POINTER, DIMENSION(:,:,:) :: ze_i_ac !: 1-D version of e_i 111 112 REAL(wp), POINTER, DIMENSION(:) :: zqbgow ! heat budget of the open water (negative) 113 REAL(wp), POINTER, DIMENSION(:) :: zdhex ! excessively thick accreted sea ice (hlead-hice) 114 115 REAL(wp), POINTER, DIMENSION(:,:,:) :: zqm0 ! old layer-system heat content 116 REAL(wp), POINTER, DIMENSION(:,:,:) :: zthick0 ! old ice thickness 117 118 REAL(wp), POINTER, DIMENSION(:,:) :: vt_i_init, vt_i_final ! ice volume summed over categories 119 REAL(wp), POINTER, DIMENSION(:,:) :: vt_s_init, vt_s_final ! snow volume summed over categories 120 REAL(wp), POINTER, DIMENSION(:,:) :: et_i_init, et_i_final ! ice energy summed over categories 121 REAL(wp), POINTER, DIMENSION(:,:) :: et_s_init ! snow energy summed over categories 109 REAL(wp), POINTER, DIMENSION(:,:) :: za_i_1d ! 1-D version of a_i 110 REAL(wp), POINTER, DIMENSION(:,:) :: zv_i_1d ! 1-D version of v_i 111 REAL(wp), POINTER, DIMENSION(:,:) :: zoa_i_1d ! 1-D version of oa_i 112 REAL(wp), POINTER, DIMENSION(:,:) :: zsmv_i_1d ! 1-D version of smv_i 113 114 REAL(wp), POINTER, DIMENSION(:,:,:) :: ze_i_1d !: 1-D version of e_i 115 122 116 REAL(wp), POINTER, DIMENSION(:,:) :: zvrel ! relative ice / frazil velocity 123 117 !!-----------------------------------------------------------------------! 124 118 125 CALL wrk_alloc( jpij, zcatac) ! integer119 CALL wrk_alloc( jpij, jcat ) ! integer 126 120 CALL wrk_alloc( jpij, zswinew, zv_newice, za_newice, zh_newice, ze_newice, zs_newice, zo_newice ) 127 CALL wrk_alloc( jpij, zdv_res, zda_res, zat_i_ac, zat_i_lev, zdh_frazb, zvrel_ac, zqbgow, zdhex ) 128 CALL wrk_alloc( jpij,jpl, zhice_old, zdummy, zdhicbot, zv_old, za_old, za_i_ac, zv_i_ac, zoa_i_ac, zsmv_i_ac ) 129 CALL wrk_alloc( jpij,jkmax,jpl, ze_i_ac ) 130 CALL wrk_alloc( jpij,jkmax+1,jpl, zqm0, zthick0 ) 131 CALL wrk_alloc( jpi,jpj, vt_i_init, vt_i_final, vt_s_init, vt_s_final, et_i_init, et_i_final, et_s_init, zvrel ) 132 133 et_i_init(:,:) = 0._wp 134 et_s_init(:,:) = 0._wp 135 vt_i_init(:,:) = 0._wp 136 vt_s_init(:,:) = 0._wp 137 138 !------------------------------------------------------------------------------! 139 ! 1) Conservation check and changes in each ice category 140 !------------------------------------------------------------------------------! 141 IF( con_i ) THEN 142 CALL lim_column_sum ( jpl, v_i , vt_i_init) 143 CALL lim_column_sum ( jpl, v_s , vt_s_init) 144 CALL lim_column_sum_energy ( jpl, nlay_i , e_i , et_i_init) 145 CALL lim_column_sum ( jpl, e_s(:,:,1,:) , et_s_init) 146 ENDIF 121 CALL wrk_alloc( jpij, zdv_res, zda_res, zat_i_1d, zat_i_lev, zv_frazb, zvrel_1d ) 122 CALL wrk_alloc( jpij,jpl, zv_old, za_old, za_i_1d, zv_i_1d, zoa_i_1d, zsmv_i_1d ) 123 CALL wrk_alloc( jpij,jkmax,jpl, ze_i_1d ) 124 CALL wrk_alloc( jpi,jpj, zvrel ) 147 125 148 126 !------------------------------------------------------------------------------| … … 154 132 DO ji = 1, jpi 155 133 !Energy of melting q(S,T) [J.m-3] 156 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / MAX( area(ji,jj) * v_i(ji,jj,jl) , epsi10 ) * REAL( nlay_i )157 134 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , -v_i(ji,jj,jl) + epsi10 ) ) !0 if no ice and 1 if yes 158 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac * zindb 135 e_i(ji,jj,jk,jl) = zindb * e_i(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_i(ji,jj,jl) , epsi10 ) ) * REAL( nlay_i ) 136 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac 159 137 END DO 160 138 END DO … … 179 157 180 158 ! Default new ice thickness 181 hicol(:,:) = hiccrit (1)182 183 IF( fraz_swi == 1 ._wp) THEN159 hicol(:,:) = hiccrit 160 161 IF( fraz_swi == 1 ) THEN 184 162 185 163 !-------------------- … … 196 174 DO ji = 1, jpi 197 175 198 IF ( tms(ji,jj) * ( qcmif(ji,jj) - qldif(ji,jj) ) > 0.e0) THEN176 IF ( qlead(ji,jj) < 0._wp ) THEN 199 177 !------------- 200 178 ! Wind stress … … 206 184 & + vtau_ice(ji ,jj ) * tmv(ji ,jj ) ) * 0.5_wp 207 185 ! Square root of wind stress 208 ztenagm = SQRT( SQRT( ztaux * ztaux + ztauy * ztauy) )186 ztenagm = SQRT( SQRT( ztaux**2 + ztauy**2 ) ) 209 187 210 188 !--------------------- 211 189 ! Frazil ice velocity 212 190 !--------------------- 213 zvfrx = zgamafr * zsqcd * ztaux / MAX(ztenagm,epsi10) 214 zvfry = zgamafr * zsqcd * ztauy / MAX(ztenagm,epsi10) 191 zindb = MAX( 0._wp, SIGN( 1._wp , ztenagm - epsi10 ) ) 192 zvfrx = zindb * zgamafr * zsqcd * ztaux / MAX( ztenagm, epsi10 ) 193 zvfry = zindb * zgamafr * zsqcd * ztauy / MAX( ztenagm, epsi10 ) 215 194 216 195 !------------------- … … 278 257 DO jj = 1, jpj 279 258 DO ji = 1, jpi 280 IF ( tms(ji,jj) * ( qcmif(ji,jj) - qldif(ji,jj) ) >0._wp ) THEN259 IF ( qlead(ji,jj) < 0._wp ) THEN 281 260 nbpac = nbpac + 1 282 261 npac( nbpac ) = (jj - 1) * jpi + ji … … 290 269 DO ji = mi0(jiindx), mi1(jiindx) 291 270 DO jj = mj0(jjindx), mj1(jjindx) 292 IF ( tms(ji,jj) * ( qcmif(ji,jj) - qldif(ji,jj) ) >0._wp ) THEN271 IF ( qlead(ji,jj) < 0._wp ) THEN 293 272 jiindex_1d = (jj - 1) * jpi + ji 294 273 ENDIF … … 307 286 IF ( nbpac > 0 ) THEN 308 287 309 CALL tab_2d_1d( nbpac, zat_i_ ac(1:nbpac) , at_i , jpi, jpj, npac(1:nbpac) )288 CALL tab_2d_1d( nbpac, zat_i_1d (1:nbpac) , at_i , jpi, jpj, npac(1:nbpac) ) 310 289 DO jl = 1, jpl 311 CALL tab_2d_1d( nbpac, za_i_ ac(1:nbpac,jl), a_i (:,:,jl), jpi, jpj, npac(1:nbpac) )312 CALL tab_2d_1d( nbpac, zv_i_ ac(1:nbpac,jl), v_i (:,:,jl), jpi, jpj, npac(1:nbpac) )313 CALL tab_2d_1d( nbpac, zoa_i_ ac(1:nbpac,jl), oa_i (:,:,jl), jpi, jpj, npac(1:nbpac) )314 CALL tab_2d_1d( nbpac, zsmv_i_ ac(1:nbpac,jl), smv_i(:,:,jl), jpi, jpj, npac(1:nbpac) )290 CALL tab_2d_1d( nbpac, za_i_1d (1:nbpac,jl), a_i (:,:,jl), jpi, jpj, npac(1:nbpac) ) 291 CALL tab_2d_1d( nbpac, zv_i_1d (1:nbpac,jl), v_i (:,:,jl), jpi, jpj, npac(1:nbpac) ) 292 CALL tab_2d_1d( nbpac, zoa_i_1d (1:nbpac,jl), oa_i (:,:,jl), jpi, jpj, npac(1:nbpac) ) 293 CALL tab_2d_1d( nbpac, zsmv_i_1d(1:nbpac,jl), smv_i(:,:,jl), jpi, jpj, npac(1:nbpac) ) 315 294 DO jk = 1, nlay_i 316 CALL tab_2d_1d( nbpac, ze_i_ ac(1:nbpac,jk,jl), e_i(:,:,jk,jl) , jpi, jpj, npac(1:nbpac) )295 CALL tab_2d_1d( nbpac, ze_i_1d(1:nbpac,jk,jl), e_i(:,:,jk,jl) , jpi, jpj, npac(1:nbpac) ) 317 296 END DO ! jk 318 297 END DO ! jl 319 298 320 CALL tab_2d_1d( nbpac, qldif_1d (1:nbpac) , qldif , jpi, jpj, npac(1:nbpac) ) 321 CALL tab_2d_1d( nbpac, qcmif_1d (1:nbpac) , qcmif , jpi, jpj, npac(1:nbpac) ) 299 CALL tab_2d_1d( nbpac, qlead_1d (1:nbpac) , qlead , jpi, jpj, npac(1:nbpac) ) 322 300 CALL tab_2d_1d( nbpac, t_bo_b (1:nbpac) , t_bo , jpi, jpj, npac(1:nbpac) ) 323 CALL tab_2d_1d( nbpac, sfx_thd_1d(1:nbpac) , sfx_thd, jpi, jpj, npac(1:nbpac) ) 324 CALL tab_2d_1d( nbpac, rdm_ice_1d(1:nbpac) , rdm_ice, jpi, jpj, npac(1:nbpac) ) 301 CALL tab_2d_1d( nbpac, sfx_opw_1d(1:nbpac) , sfx_opw, jpi, jpj, npac(1:nbpac) ) 302 CALL tab_2d_1d( nbpac, wfx_opw_1d(1:nbpac) , wfx_opw, jpi, jpj, npac(1:nbpac) ) 303 CALL tab_2d_1d( nbpac, wfx_opw_1d(1:nbpac) , wfx_opw, jpi, jpj, npac(1:nbpac) ) 325 304 CALL tab_2d_1d( nbpac, hicol_b (1:nbpac) , hicol , jpi, jpj, npac(1:nbpac) ) 326 CALL tab_2d_1d( nbpac, zvrel_ac (1:nbpac) , zvrel , jpi, jpj, npac(1:nbpac) ) 305 CALL tab_2d_1d( nbpac, zvrel_1d (1:nbpac) , zvrel , jpi, jpj, npac(1:nbpac) ) 306 307 CALL tab_2d_1d( nbpac, hfx_thd_1d(1:nbpac) , hfx_thd, jpi, jpj, npac(1:nbpac) ) 308 CALL tab_2d_1d( nbpac, hfx_opw_1d(1:nbpac) , hfx_opw, jpi, jpj, npac(1:nbpac) ) 327 309 328 310 !------------------------------------------------------------------------------! … … 330 312 !------------------------------------------------------------------------------! 331 313 314 !----------------------------------------- 315 ! Keep old ice areas and volume in memory 316 !----------------------------------------- 317 zv_old(:,:) = zv_i_1d(:,:) 318 za_old(:,:) = za_i_1d(:,:) 319 332 320 !---------------------- 333 321 ! Thickness of new ice 334 322 !---------------------- 335 323 DO ji = 1, nbpac 336 zh_newice(ji) = hiccrit (1)337 END DO 338 IF( fraz_swi == 1 .0 )zh_newice(:) = hicol_b(:)324 zh_newice(ji) = hiccrit 325 END DO 326 IF( fraz_swi == 1 ) zh_newice(:) = hicol_b(:) 339 327 340 328 !---------------------- 341 329 ! Salinity of new ice 342 330 !---------------------- 343 344 331 SELECT CASE ( num_sal ) 345 332 CASE ( 1 ) ! Sice = constant … … 355 342 END SELECT 356 343 357 358 344 !------------------------- 359 345 ! Heat content of new ice … … 363 349 ztmelts = - tmut * zs_newice(ji) + rtt ! Melting point (K) 364 350 ze_newice(ji) = rhoic * ( cpic * ( ztmelts - t_bo_b(ji) ) & 365 & + lfus * ( 1.0 - ( ztmelts - rtt ) / ( t_bo_b(ji) - rtt) ) &351 & + lfus * ( 1.0 - ( ztmelts - rtt ) / MIN( t_bo_b(ji) - rtt, -epsi10 ) ) & 366 352 & - rcp * ( ztmelts - rtt ) ) 367 ze_newice(ji) = MAX( ze_newice(ji) , 0._wp ) &368 & + MAX( 0.0 , SIGN( 1.0 , - ze_newice(ji) ) ) * rhoic * lfus369 353 END DO ! ji 354 370 355 !---------------- 371 356 ! Age of new ice … … 375 360 END DO ! ji 376 361 377 !--------------------------378 ! Open water energy budget379 !--------------------------380 DO ji = 1, nbpac381 zqbgow(ji) = qldif_1d(ji) - qcmif_1d(ji) !<0382 END DO ! ji383 384 362 !------------------- 385 363 ! Volume of new ice 386 364 !------------------- 387 365 DO ji = 1, nbpac 388 zv_newice(ji) = - zqbgow(ji) / ze_newice(ji) 366 367 zEi = - ze_newice(ji) / rhoic ! specific enthalpy of forming ice [J/kg] 368 369 zEw = rcp * ( t_bo_b(ji) - rt0 ) ! specific enthalpy of seawater at t_bo_b [J/kg] 370 ! clem: we suppose we are already at the freezing point (condition qlead<0 is satisfyied) 371 372 zdE = zEi - zEw ! specific enthalpy difference [J/kg] 373 374 zfmdt = - qlead_1d(ji) / zdE ! Fm.dt [kg/m2] (<0) 375 ! clem: we use qlead instead of zqld (limthd) because we suppose we are at the freezing point 376 zv_newice(ji) = - zfmdt / rhoic 377 378 zQm = zfmdt * zEw ! heat to the ocean >0 associated with mass flux 379 380 ! Contribution to heat flux to the ocean [W.m-2], >0 381 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * zEw * r1_rdtice 382 ! Total heat flux used in this process [W.m-2] 383 hfx_opw_1d(ji) = hfx_opw_1d(ji) - zfmdt * zdE * r1_rdtice 384 ! mass flux 385 wfx_opw_1d(ji) = wfx_opw_1d(ji) - zv_newice(ji) * rhoic * r1_rdtice 386 ! salt flux 387 sfx_opw_1d(ji) = sfx_opw_1d(ji) - zv_newice(ji) * rhoic * zs_newice(ji) * r1_rdtice 389 388 390 389 ! A fraction zfrazb of frazil ice is accreted at the ice bottom 391 zfrazb = ( TANH ( Cfrazb * ( zvrel_ac(ji) - vfrazb ) ) + 1.0 ) * 0.5 * maxfrazb 392 zdh_frazb(ji) = zfrazb * zv_newice(ji) 390 zinda = 1._wp - MAX( 0._wp, SIGN( 1._wp , - zat_i_1d(ji) ) ) 391 zfrazb = zinda * ( TANH ( Cfrazb * ( zvrel_1d(ji) - vfrazb ) ) + 1.0 ) * 0.5 * maxfrazb 392 zv_frazb(ji) = zfrazb * zv_newice(ji) 393 393 zv_newice(ji) = ( 1.0 - zfrazb ) * zv_newice(ji) 394 394 END DO 395 396 !------------------------------------397 ! Diags for energy conservation test398 !------------------------------------399 DO ji = 1, nbpac400 ii = MOD( npac(ji) - 1 , jpi ) + 1401 ij = ( npac(ji) - 1 ) / jpi + 1402 !403 zde = ze_newice(ji) / unit_fac * area(ii,ij) * zv_newice(ji)404 !405 vt_i_init(ii,ij) = vt_i_init(ii,ij) + zv_newice(ji) ! volume406 et_i_init(ii,ij) = et_i_init(ii,ij) + zde ! Energy407 408 END DO409 410 ! keep new ice volume in memory411 CALL tab_1d_2d( nbpac, v_newice , npac(1:nbpac), zv_newice(1:nbpac) , jpi, jpj )412 395 413 396 !----------------- … … 415 398 !----------------- 416 399 DO ji = 1, nbpac 417 ii = MOD( npac(ji) - 1 , jpi ) + 1418 ij = ( npac(ji) - 1 ) / jpi + 1419 400 za_newice(ji) = zv_newice(ji) / zh_newice(ji) 420 diag_lat_gr(ii,ij) = diag_lat_gr(ii,ij) + zv_newice(ji) * r1_rdtice ! clem 421 END DO !ji 401 END DO 422 402 423 403 !------------------------------------------------------------------------------! … … 425 405 !------------------------------------------------------------------------------! 426 406 427 !----------------------------------------- 428 ! Keep old ice areas and volume in memory 429 !----------------------------------------- 430 zv_old(:,:) = zv_i_ac(:,:) 431 za_old(:,:) = za_i_ac(:,:) 432 433 !------------------------------------------- 434 ! Compute excessive new ice area and volume 435 !------------------------------------------- 407 !------------------------ 408 ! 6.1) lateral ice growth 409 !------------------------ 436 410 ! If lateral ice growth gives an ice concentration gt 1, then 437 411 ! we keep the excessive volume in memory and attribute it later to bottom accretion 438 412 DO ji = 1, nbpac 439 IF ( za_newice(ji) > ( amax - zat_i_ ac(ji) ) ) THEN440 zda_res(ji) = za_newice(ji) - ( amax - zat_i_ ac(ji) )413 IF ( za_newice(ji) > ( amax - zat_i_1d(ji) ) ) THEN 414 zda_res(ji) = za_newice(ji) - ( amax - zat_i_1d(ji) ) 441 415 zdv_res(ji) = zda_res (ji) * zh_newice(ji) 442 416 za_newice(ji) = za_newice(ji) - zda_res (ji) … … 446 420 zdv_res(ji) = 0._wp 447 421 ENDIF 448 END DO ! ji 449 450 !------------------------------------------------ 451 ! Laterally redistribute new ice volume and area 452 !------------------------------------------------ 453 zat_i_ac(:) = 0._wp 422 END DO 423 424 ! find which category to fill 425 zat_i_1d(:) = 0._wp 454 426 DO jl = 1, jpl 455 427 DO ji = 1, nbpac 456 IF( hi_max (jl-1) < zh_newice(ji) .AND. & 457 & zh_newice(ji) <= hi_max (jl) ) THEN 458 za_i_ac (ji,jl) = za_i_ac (ji,jl) + za_newice(ji) 459 zv_i_ac (ji,jl) = zv_i_ac (ji,jl) + zv_newice(ji) 460 zat_i_ac(ji) = zat_i_ac(ji) + za_i_ac (ji,jl) 461 zcatac (ji) = jl 428 IF( zh_newice(ji) > hi_max(jl-1) .AND. zh_newice(ji) <= hi_max(jl) ) THEN 429 za_i_1d (ji,jl) = za_i_1d (ji,jl) + za_newice(ji) 430 zv_i_1d (ji,jl) = zv_i_1d (ji,jl) + zv_newice(ji) 431 jcat (ji) = jl 462 432 ENDIF 463 END DO 464 END DO 465 466 !---------------------------------- 467 ! Heat content - lateral accretion 468 !---------------------------------- 469 DO ji = 1, nbpac 470 jl = zcatac(ji) ! categroy in which new ice is put 471 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , -za_old(ji,jl) + epsi10 ) ) ! zindb=1 if ice =0 otherwise 472 zhice_old(ji,jl) = zv_old(ji,jl) / MAX( za_old(ji,jl) , epsi10 ) * zindb ! old ice thickness 473 zdhex (ji) = MAX( 0._wp , zh_newice(ji) - zhice_old(ji,jl) ) ! difference in thickness 474 zswinew (ji) = MAX( 0._wp , SIGN( 1._wp , - za_old(ji,jl) + epsi10 ) ) ! ice totally new in jl category 433 zat_i_1d(ji) = zat_i_1d(ji) + za_i_1d (ji,jl) 434 END DO 435 END DO 436 437 ! Heat content 438 DO ji = 1, nbpac 439 jl = jcat(ji) ! categroy in which new ice is put 440 zswinew (ji) = MAX( 0._wp , SIGN( 1._wp , - za_old(ji,jl) ) ) ! 0 if old ice 475 441 END DO 476 442 477 443 DO jk = 1, nlay_i 478 444 DO ji = 1, nbpac 479 jl = zcatac(ji) 480 zqold = ze_i_ac(ji,jk,jl) ! [ J.m-3 ] 481 zalphai = MIN( zhice_old(ji,jl) * REAL( jk ) / REAL( nlay_i ), zh_newice(ji) ) & 482 & - MIN( zhice_old(ji,jl) * REAL( jk - 1 ) / REAL( nlay_i ), zh_newice(ji) ) 483 ze_i_ac(ji,jk,jl) = zswinew(ji) * ze_newice(ji) & 484 + ( 1.0 - zswinew(ji) ) * ( za_old(ji,jl) * zqold * zhice_old(ji,jl) / REAL( nlay_i ) & 485 + za_newice(ji) * ze_newice(ji) * zalphai & 486 + za_newice(ji) * ze_newice(ji) * zdhex(ji) / REAL( nlay_i ) ) / ( ( zv_i_ac(ji,jl) ) / REAL( nlay_i ) ) 487 END DO 488 END DO 489 490 !----------------------------------------------- 491 ! Add excessive volume of new ice at the bottom 492 !----------------------------------------------- 493 ! If the ice concentration exceeds 1, the remaining volume of new ice 494 ! is equally redistributed among all ice categories in which there is 495 ! ice 496 497 ! Fraction of level ice 498 jm = 1 499 zat_i_lev(:) = 0._wp 500 501 DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) 502 DO ji = 1, nbpac 503 zat_i_lev(ji) = zat_i_lev(ji) + za_i_ac(ji,jl) 504 END DO 505 END DO 506 507 IF( ln_nicep .AND. jiindex_1d > 0 ) WRITE(numout,*) ' zv_i_ac : ', zv_i_ac(jiindex_1d, 1:jpl) 508 DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) 509 DO ji = 1, nbpac 510 zindb = MAX( 0._wp, SIGN( 1._wp , zdv_res(ji) ) ) 511 zinda = MAX( 0._wp, SIGN( 1._wp , zat_i_lev(ji) - epsi10 ) ) ! clem 512 zv_i_ac(ji,jl) = zv_i_ac(ji,jl) + zindb * zinda * zdv_res(ji) * za_i_ac(ji,jl) / MAX( zat_i_lev(ji) , epsi10 ) 513 END DO 514 END DO 515 IF( ln_nicep .AND. jiindex_1d > 0 ) WRITE(numout,*) ' zv_i_ac : ', zv_i_ac(jiindex_1d, 1:jpl) 516 517 !--------------------------------- 518 ! Heat content - bottom accretion 519 !--------------------------------- 520 jm = 1 521 DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) 522 DO ji = 1, nbpac 523 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl ) + epsi10 ) ) ! zindb=1 if ice =0 otherwise 524 zhice_old(ji,jl) = zv_i_ac(ji,jl) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb 525 zdhicbot (ji,jl) = zdv_res(ji) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb & 526 & + zindb * zdh_frazb(ji) ! frazil ice may coalesce 527 zdummy(ji,jl) = zv_i_ac(ji,jl) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb ! thickness of residual ice 528 END DO 529 END DO 530 531 ! old layers thicknesses and enthalpies 532 DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) 445 jl = jcat(ji) 446 zinda = MAX( 0._wp, SIGN( 1._wp , zv_i_1d(ji,jl) - epsi20 ) ) 447 ze_i_1d(ji,jk,jl) = zswinew(ji) * ze_newice(ji) + & 448 & ( 1.0 - zswinew(ji) ) * ( ze_newice(ji) * zv_newice(ji) + ze_i_1d(ji,jk,jl) * zv_old(ji,jl) ) & 449 & * zinda / MAX( zv_i_1d(ji,jl), epsi20 ) 450 END DO 451 END DO 452 453 !------------------------------------------------ 454 ! 6.2) bottom ice growth + ice enthalpy remapping 455 !------------------------------------------------ 456 DO jl = 1, jpl 457 458 ! for remapping 459 h_i_old (1:nbpac,0:nlay_i+1) = 0._wp 460 qh_i_old(1:nbpac,0:nlay_i+1) = 0._wp 533 461 DO jk = 1, nlay_i 534 462 DO ji = 1, nbpac 535 zthick0(ji,jk,jl) = zhice_old(ji,jl) / REAL( nlay_i )536 zqm0 (ji,jk,jl) = ze_i_ac(ji,jk,jl) * zthick0(ji,jk,jl)463 h_i_old (ji,jk) = zv_i_1d(ji,jl) / REAL( nlay_i ) 464 qh_i_old(ji,jk) = ze_i_1d(ji,jk,jl) * h_i_old(ji,jk) 537 465 END DO 538 466 END DO 539 END DO 540 !!gm ??? why the previous do loop if ocerwriten by the following one ? 541 DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) 467 468 ! new volumes including lateral/bottom accretion + residual 542 469 DO ji = 1, nbpac 543 zthick0(ji,nlay_i+1,jl) = zdhicbot(ji,jl) 544 zqm0 (ji,nlay_i+1,jl) = ze_newice(ji) * zdhicbot(ji,jl) 545 END DO ! ji 546 END DO ! jl 547 548 ! Redistributing energy on the new grid 549 ze_i_ac(:,:,:) = 0._wp 550 DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) 551 DO jk = 1, nlay_i 552 DO layer = 1, nlay_i + 1 553 DO ji = 1, nbpac 554 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl) + epsi10 ) ) 555 ! Redistributing energy on the new grid 556 zweight = MAX ( MIN( zhice_old(ji,jl) * REAL( layer ), zdummy(ji,jl) * REAL( jk ) ) & 557 & - MAX( zhice_old(ji,jl) * REAL( layer - 1 ) , zdummy(ji,jl) * REAL( jk - 1 ) ) , 0._wp ) & 558 & /( MAX(REAL(nlay_i) * zthick0(ji,layer,jl),epsi10) ) * zindb 559 ze_i_ac(ji,jk,jl) = ze_i_ac(ji,jk,jl) + zweight * zqm0(ji,layer,jl) 560 END DO ! ji 561 END DO ! layer 562 END DO ! jk 563 END DO ! jl 564 565 DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) 566 DO jk = 1, nlay_i 567 DO ji = 1, nbpac 568 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) ) 569 ze_i_ac(ji,jk,jl) = ze_i_ac(ji,jk,jl) & 570 & / MAX( zv_i_ac(ji,jl) , epsi10) * za_i_ac(ji,jl) * REAL( nlay_i ) * zindb 571 END DO 572 END DO 573 END DO 470 zinda = MAX( 0._wp, SIGN( 1._wp , zat_i_1d(ji) - epsi20 ) ) 471 zv_newfra = zinda * ( zdv_res(ji) + zv_frazb(ji) ) * za_i_1d(ji,jl) / MAX( zat_i_1d(ji) , epsi20 ) 472 za_i_1d(ji,jl) = zinda * za_i_1d(ji,jl) 473 zv_i_1d(ji,jl) = zv_i_1d(ji,jl) + zv_newfra 474 475 ! for remapping 476 h_i_old (ji,nlay_i+1) = zv_newfra 477 qh_i_old(ji,nlay_i+1) = ze_newice(ji) * zv_newfra 478 ENDDO 479 480 ! --- Ice enthalpy remapping --- ! 481 IF( zv_newfra > 0._wp ) THEN 482 CALL lim_thd_ent( 1, nbpac, ze_i_1d(1:nbpac,:,jl) ) 483 ENDIF 484 485 ENDDO 574 486 575 487 !------------ … … 578 490 DO jl = 1, jpl 579 491 DO ji = 1, nbpac 580 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ ac(ji,jl) + epsi10 ) ) ! 0 if no ice and 1 if yes581 zoa_i_ ac(ji,jl) = za_old(ji,jl) * zoa_i_ac(ji,jl) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb492 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_1d(ji,jl) + epsi20 ) ) ! 0 if no ice and 1 if yes 493 zoa_i_1d(ji,jl) = za_old(ji,jl) * zoa_i_1d(ji,jl) / MAX( za_i_1d(ji,jl) , epsi20 ) * zindb 582 494 END DO 583 495 END DO … … 586 498 ! Update salinity 587 499 !----------------- 588 !clem IF( num_sal == 2 ) THEN589 DO jl = 1, jpl590 DO ji = 1, nbpac591 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) ) ! 0 if no ice and 1 if yes592 zdv = zv_i_ac(ji,jl) - zv_old(ji,jl)593 zsmv_i_ac(ji,jl) = zsmv_i_ac(ji,jl) + zdv * zs_newice(ji) * zindb ! clem modif594 END DO595 END DO596 !clem ENDIF597 598 !--------------------------------599 ! Update mass/salt fluxes (clem)600 !--------------------------------601 500 DO jl = 1, jpl 602 501 DO ji = 1, nbpac 603 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) ) ! 0 if no ice and 1 if yes 604 zdv = zv_i_ac(ji,jl) - zv_old(ji,jl) 605 rdm_ice_1d(ji) = rdm_ice_1d(ji) + zdv * rhoic * zindb 606 sfx_thd_1d(ji) = sfx_thd_1d(ji) - zdv * rhoic * zs_newice(ji) * r1_rdtice * zindb 607 END DO 502 zdv = zv_i_1d(ji,jl) - zv_old(ji,jl) 503 zsmv_i_1d(ji,jl) = zsmv_i_1d(ji,jl) + zdv * zs_newice(ji) 504 END DO 608 505 END DO 609 506 610 507 !------------------------------------------------------------------------------! 611 ! 8) Change 2D vectors to 1D vectors508 ! 7) Change 2D vectors to 1D vectors 612 509 !------------------------------------------------------------------------------! 613 510 DO jl = 1, jpl 614 CALL tab_1d_2d( nbpac, a_i (:,:,jl), npac(1:nbpac), za_i_ac (1:nbpac,jl), jpi, jpj ) 615 CALL tab_1d_2d( nbpac, v_i (:,:,jl), npac(1:nbpac), zv_i_ac (1:nbpac,jl), jpi, jpj ) 616 CALL tab_1d_2d( nbpac, oa_i(:,:,jl), npac(1:nbpac), zoa_i_ac(1:nbpac,jl), jpi, jpj ) 617 !clem IF ( num_sal == 2 ) & 618 CALL tab_1d_2d( nbpac, smv_i (:,:,jl), npac(1:nbpac), zsmv_i_ac(1:nbpac,jl) , jpi, jpj ) 511 CALL tab_1d_2d( nbpac, a_i (:,:,jl), npac(1:nbpac), za_i_1d (1:nbpac,jl), jpi, jpj ) 512 CALL tab_1d_2d( nbpac, v_i (:,:,jl), npac(1:nbpac), zv_i_1d (1:nbpac,jl), jpi, jpj ) 513 CALL tab_1d_2d( nbpac, oa_i(:,:,jl), npac(1:nbpac), zoa_i_1d(1:nbpac,jl), jpi, jpj ) 514 CALL tab_1d_2d( nbpac, smv_i (:,:,jl), npac(1:nbpac), zsmv_i_1d(1:nbpac,jl) , jpi, jpj ) 619 515 DO jk = 1, nlay_i 620 CALL tab_1d_2d( nbpac, e_i(:,:,jk,jl), npac(1:nbpac), ze_i_ac(1:nbpac,jk,jl), jpi, jpj ) 621 END DO 622 END DO 623 CALL tab_1d_2d( nbpac, sfx_thd, npac(1:nbpac), sfx_thd_1d(1:nbpac), jpi, jpj ) 624 CALL tab_1d_2d( nbpac, rdm_ice, npac(1:nbpac), rdm_ice_1d(1:nbpac), jpi, jpj ) 516 CALL tab_1d_2d( nbpac, e_i(:,:,jk,jl), npac(1:nbpac), ze_i_1d(1:nbpac,jk,jl), jpi, jpj ) 517 END DO 518 END DO 519 CALL tab_1d_2d( nbpac, sfx_opw, npac(1:nbpac), sfx_opw_1d(1:nbpac), jpi, jpj ) 520 CALL tab_1d_2d( nbpac, wfx_opw, npac(1:nbpac), wfx_opw_1d(1:nbpac), jpi, jpj ) 521 CALL tab_1d_2d( nbpac, wfx_opw, npac(1:nbpac), wfx_opw_1d(1:nbpac), jpi, jpj ) 522 523 CALL tab_1d_2d( nbpac, hfx_thd, npac(1:nbpac), hfx_thd_1d(1:nbpac), jpi, jpj ) 524 CALL tab_1d_2d( nbpac, hfx_opw, npac(1:nbpac), hfx_opw_1d(1:nbpac), jpi, jpj ) 625 525 ! 626 526 ENDIF ! nbpac > 0 627 527 628 528 !------------------------------------------------------------------------------! 629 ! 9) Change units for e_i529 ! 8) Change units for e_i 630 530 !------------------------------------------------------------------------------! 631 531 DO jl = 1, jpl 632 DO jk = 1, nlay_i ! heat content in 10^9 Joules 633 e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * area(:,:) * v_i(:,:,jl) / REAL( nlay_i ) / unit_fac 532 DO jk = 1, nlay_i 533 DO jj = 1, jpj 534 DO ji = 1, jpi 535 ! heat content in Joules 536 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * area(ji,jj) * v_i(ji,jj,jl) / ( REAL( nlay_i ) * unit_fac ) 537 END DO 538 END DO 634 539 END DO 635 540 END DO 636 541 637 !------------------------------------------------------------------------------|638 ! 10) Conservation check and changes in each ice category639 !------------------------------------------------------------------------------|640 IF( con_i ) THEN641 CALL lim_column_sum (jpl, v_i, vt_i_final)642 fieldid = 'v_i, limthd_lac'643 CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid)644 !645 CALL lim_column_sum_energy(jpl, nlay_i, e_i, et_i_final)646 fieldid = 'e_i, limthd_lac'647 CALL lim_cons_check (et_i_final, et_i_final, 1.0e-3, fieldid)648 !649 CALL lim_column_sum (jpl, v_s, vt_s_final)650 fieldid = 'v_s, limthd_lac'651 CALL lim_cons_check (vt_s_init, vt_s_final, 1.0e-6, fieldid)652 !653 ! CALL lim_column_sum (jpl, e_s(:,:,1,:) , et_s_init)654 ! fieldid = 'e_s, limthd_lac'655 ! CALL lim_cons_check (et_s_init, et_s_final, 1.0e-3, fieldid)656 IF( ln_nicep ) THEN657 DO ji = mi0(jiindx), mi1(jiindx)658 DO jj = mj0(jjindx), mj1(jjindx)659 WRITE(numout,*) ' vt_i_init : ', vt_i_init (ji,jj)660 WRITE(numout,*) ' vt_i_final: ', vt_i_final(ji,jj)661 WRITE(numout,*) ' et_i_init : ', et_i_init (ji,jj)662 WRITE(numout,*) ' et_i_final: ', et_i_final(ji,jj)663 END DO664 END DO665 ENDIF666 !667 ENDIF668 542 ! 669 CALL wrk_dealloc( jpij, zcatac) ! integer543 CALL wrk_dealloc( jpij, jcat ) ! integer 670 544 CALL wrk_dealloc( jpij, zswinew, zv_newice, za_newice, zh_newice, ze_newice, zs_newice, zo_newice ) 671 CALL wrk_dealloc( jpij, zdv_res, zda_res, zat_i_ac, zat_i_lev, zdh_frazb, zvrel_ac, zqbgow, zdhex ) 672 CALL wrk_dealloc( jpij,jpl, zhice_old, zdummy, zdhicbot, zv_old, za_old, za_i_ac, zv_i_ac, zoa_i_ac, zsmv_i_ac ) 673 CALL wrk_dealloc( jpij,jkmax,jpl, ze_i_ac ) 674 CALL wrk_dealloc( jpij,jkmax+1,jpl, zqm0, zthick0 ) 675 CALL wrk_dealloc( jpi,jpj, vt_i_init, vt_i_final, vt_s_init, vt_s_final, et_i_init, et_i_final, et_s_init, zvrel ) 545 CALL wrk_dealloc( jpij, zdv_res, zda_res, zat_i_1d, zat_i_lev, zv_frazb, zvrel_1d ) 546 CALL wrk_dealloc( jpij,jpl, zv_old, za_old, za_i_1d, zv_i_1d, zoa_i_1d, zsmv_i_1d ) 547 CALL wrk_dealloc( jpij,jkmax,jpl, ze_i_1d ) 548 CALL wrk_dealloc( jpi,jpj, zvrel ) 676 549 ! 677 550 END SUBROUTINE lim_thd_lac
Note: See TracChangeset
for help on using the changeset viewer.
