Changeset 4921 for branches/2014/dev_r4650_UKMO13_CICE_changes_take2/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
- Timestamp:
- 2014-11-28T14:59:01+01:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4650_UKMO13_CICE_changes_take2/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
r4333 r4921 29 29 USE lib_mpp ! MPP library 30 30 USE wrk_nemo ! work arrays 31 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 31 32 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 33 USE limthd_ent 32 34 33 35 IMPLICIT NONE … … 37 39 38 40 REAL(wp) :: epsi10 = 1.e-10_wp ! 39 REAL(wp) :: zzero = 0._wp ! 40 REAL(wp) :: zone = 1._wp ! 41 REAL(wp) :: epsi20 = 1.e-20_wp ! 41 42 42 43 !!---------------------------------------------------------------------- … … 71 72 !! - Computation of variation of ice volume and mass 72 73 !! - Computation of frldb after lateral accretion and 73 !! update ht_s_ b, ht_i_band tbif_1d(:,:)74 !! update ht_s_1d, ht_i_1d and tbif_1d(:,:) 74 75 !!------------------------------------------------------------------------ 75 INTEGER :: ji,jj,jk,jl ,jm! dummy loop indices76 INTEGER :: layer, nbpac! local integers77 INTEGER :: ii, ij, iter ! - -78 REAL(wp) :: ztmelts, zdv, z qold, zfrazb, zweight, zalphai, zindb, zinda, zde ! local scalars76 INTEGER :: ji,jj,jk,jl ! dummy loop indices 77 INTEGER :: nbpac ! local integers 78 INTEGER :: ii, ij, iter ! - - 79 REAL(wp) :: ztmelts, zdv, zfrazb, zweight, zindb, zinda, zde ! local scalars 79 80 REAL(wp) :: zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf , zhicol_new ! - - 80 81 REAL(wp) :: ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, zsqcd , zhicrit ! - - 81 82 LOGICAL :: iterate_frazil ! iterate frazil ice collection thickness 82 83 CHARACTER (len = 15) :: fieldid 83 ! 84 INTEGER , POINTER, DIMENSION(:) :: zcatac ! indexes of categories where new ice grows 84 85 REAL(wp) :: zQm ! enthalpy exchanged with the ocean (J/m2, >0 towards ocean) 86 REAL(wp) :: zEi ! sea ice specific enthalpy (J/kg) 87 REAL(wp) :: zEw ! seawater specific enthalpy (J/kg) 88 REAL(wp) :: zfmdt ! mass flux x time step (kg/m2, >0 towards ocean) 89 90 REAL(wp) :: zv_newfra 91 92 INTEGER , POINTER, DIMENSION(:) :: jcat ! indexes of categories where new ice grows 85 93 REAL(wp), POINTER, DIMENSION(:) :: zswinew ! switch for new ice or not 86 94 … … 93 101 REAL(wp), POINTER, DIMENSION(:) :: zdv_res ! residual volume in case of excessive heat budget 94 102 REAL(wp), POINTER, DIMENSION(:) :: zda_res ! residual area in case of excessive heat budget 95 REAL(wp), POINTER, DIMENSION(:) :: zat_i_ac ! total ice fraction 96 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 103 REAL(wp), POINTER, DIMENSION(:,:) :: zv_old ! old volume of ice in category jl 104 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 103 REAL(wp), POINTER, DIMENSION(:) :: zat_i_1d ! total ice fraction 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 107 REAL(wp), POINTER, DIMENSION(:,:) :: zv_b ! old volume of ice in category jl 108 REAL(wp), POINTER, DIMENSION(:,:) :: za_b ! old area of ice in category jl 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, zv_frazb, zvrel_1d ) 122 CALL wrk_alloc( jpij,jpl, zv_b, za_b, za_i_1d, zv_i_1d, zoa_i_1d, zsmv_i_1d ) 123 CALL wrk_alloc( jpij,nlay_i+1,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, wp ) 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 !-------------------- … … 193 171 zgamafr = 0.03 194 172 195 DO jj = 1, jpj 196 DO ji = 1, jpi 197 198 IF ( tms(ji,jj) * ( qcmif(ji,jj) - qldif(ji,jj) ) > 0.e0 ) THEN 173 DO jj = 2, jpj 174 DO ji = 2, jpi 175 IF ( qlead(ji,jj) < 0._wp ) THEN 199 176 !------------- 200 177 ! Wind stress … … 206 183 & + vtau_ice(ji ,jj ) * tmv(ji ,jj ) ) * 0.5_wp 207 184 ! Square root of wind stress 208 ztenagm = SQRT( SQRT( ztaux * ztaux + ztauy * ztauy) )185 ztenagm = SQRT( SQRT( ztaux**2 + ztauy**2 ) ) 209 186 210 187 !--------------------- 211 188 ! Frazil ice velocity 212 189 !--------------------- 213 zvfrx = zgamafr * zsqcd * ztaux / MAX(ztenagm,epsi10) 214 zvfry = zgamafr * zsqcd * ztauy / MAX(ztenagm,epsi10) 190 zindb = MAX( 0._wp, SIGN( 1._wp , ztenagm - epsi10 ) ) 191 zvfrx = zindb * zgamafr * zsqcd * ztaux / MAX( ztenagm, epsi10 ) 192 zvfry = zindb * zgamafr * zsqcd * ztauy / MAX( ztenagm, epsi10 ) 215 193 216 194 !------------------- … … 264 242 END DO ! loop on ji ends 265 243 END DO ! loop on jj ends 244 ! 245 CALL lbc_lnk( zvrel(:,:), 'T', 1. ) 246 CALL lbc_lnk( hicol(:,:), 'T', 1. ) 266 247 267 248 ENDIF ! End of computation of frazil ice collection thickness … … 276 257 ! This occurs if open water energy budget is negative 277 258 nbpac = 0 259 npac(:) = 0 260 ! 278 261 DO jj = 1, jpj 279 262 DO ji = 1, jpi 280 IF ( tms(ji,jj) * ( qcmif(ji,jj) - qldif(ji,jj) ) >0._wp ) THEN263 IF ( qlead(ji,jj) < 0._wp ) THEN 281 264 nbpac = nbpac + 1 282 265 npac( nbpac ) = (jj - 1) * jpi + ji … … 290 273 DO ji = mi0(jiindx), mi1(jiindx) 291 274 DO jj = mj0(jjindx), mj1(jjindx) 292 IF ( tms(ji,jj) * ( qcmif(ji,jj) - qldif(ji,jj) ) >0._wp ) THEN275 IF ( qlead(ji,jj) < 0._wp ) THEN 293 276 jiindex_1d = (jj - 1) * jpi + ji 294 277 ENDIF … … 307 290 IF ( nbpac > 0 ) THEN 308 291 309 CALL tab_2d_1d( nbpac, zat_i_ ac(1:nbpac) , at_i , jpi, jpj, npac(1:nbpac) )292 CALL tab_2d_1d( nbpac, zat_i_1d (1:nbpac) , at_i , jpi, jpj, npac(1:nbpac) ) 310 293 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) )294 CALL tab_2d_1d( nbpac, za_i_1d (1:nbpac,jl), a_i (:,:,jl), jpi, jpj, npac(1:nbpac) ) 295 CALL tab_2d_1d( nbpac, zv_i_1d (1:nbpac,jl), v_i (:,:,jl), jpi, jpj, npac(1:nbpac) ) 296 CALL tab_2d_1d( nbpac, zoa_i_1d (1:nbpac,jl), oa_i (:,:,jl), jpi, jpj, npac(1:nbpac) ) 297 CALL tab_2d_1d( nbpac, zsmv_i_1d(1:nbpac,jl), smv_i(:,:,jl), jpi, jpj, npac(1:nbpac) ) 315 298 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) )299 CALL tab_2d_1d( nbpac, ze_i_1d(1:nbpac,jk,jl), e_i(:,:,jk,jl) , jpi, jpj, npac(1:nbpac) ) 317 300 END DO ! jk 318 301 END DO ! jl 319 302 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) ) 322 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) ) 325 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) ) 303 CALL tab_2d_1d( nbpac, qlead_1d (1:nbpac) , qlead , jpi, jpj, npac(1:nbpac) ) 304 CALL tab_2d_1d( nbpac, t_bo_1d (1:nbpac) , t_bo , jpi, jpj, npac(1:nbpac) ) 305 CALL tab_2d_1d( nbpac, sfx_opw_1d(1:nbpac) , sfx_opw, jpi, jpj, npac(1:nbpac) ) 306 CALL tab_2d_1d( nbpac, wfx_opw_1d(1:nbpac) , wfx_opw, jpi, jpj, npac(1:nbpac) ) 307 CALL tab_2d_1d( nbpac, hicol_1d (1:nbpac) , hicol , jpi, jpj, npac(1:nbpac) ) 308 CALL tab_2d_1d( nbpac, zvrel_1d (1:nbpac) , zvrel , jpi, jpj, npac(1:nbpac) ) 309 310 CALL tab_2d_1d( nbpac, hfx_thd_1d(1:nbpac) , hfx_thd, jpi, jpj, npac(1:nbpac) ) 311 CALL tab_2d_1d( nbpac, hfx_opw_1d(1:nbpac) , hfx_opw, jpi, jpj, npac(1:nbpac) ) 327 312 328 313 !------------------------------------------------------------------------------! … … 330 315 !------------------------------------------------------------------------------! 331 316 317 !----------------------------------------- 318 ! Keep old ice areas and volume in memory 319 !----------------------------------------- 320 zv_b(1:nbpac,:) = zv_i_1d(1:nbpac,:) 321 za_b(1:nbpac,:) = za_i_1d(1:nbpac,:) 332 322 !---------------------- 333 323 ! Thickness of new ice 334 324 !---------------------- 335 325 DO ji = 1, nbpac 336 zh_newice(ji) = hiccrit (1)337 END DO 338 IF( fraz_swi == 1 .0 ) zh_newice(:) = hicol_b(:)326 zh_newice(ji) = hiccrit 327 END DO 328 IF( fraz_swi == 1 ) zh_newice(1:nbpac) = hicol_1d(1:nbpac) 339 329 340 330 !---------------------- 341 331 ! Salinity of new ice 342 332 !---------------------- 343 344 333 SELECT CASE ( num_sal ) 345 334 CASE ( 1 ) ! Sice = constant 346 zs_newice( :) = bulk_sal335 zs_newice(1:nbpac) = bulk_sal 347 336 CASE ( 2 ) ! Sice = F(z,t) [Vancoppenolle et al (2005)] 348 337 DO ji = 1, nbpac … … 352 341 END DO 353 342 CASE ( 3 ) ! Sice = F(z) [multiyear ice] 354 zs_newice( :) = 2.3343 zs_newice(1:nbpac) = 2.3 355 344 END SELECT 356 357 345 358 346 !------------------------- … … 362 350 DO ji = 1, nbpac 363 351 ztmelts = - tmut * zs_newice(ji) + rtt ! Melting point (K) 364 ze_newice(ji) = rhoic * ( cpic * ( ztmelts - t_bo_ b(ji) ) &365 & + lfus * ( 1.0 - ( ztmelts - rtt ) / ( t_bo_b(ji) - rtt) ) &352 ze_newice(ji) = rhoic * ( cpic * ( ztmelts - t_bo_1d(ji) ) & 353 & + lfus * ( 1.0 - ( ztmelts - rtt ) / MIN( t_bo_1d(ji) - rtt, -epsi10 ) ) & 366 354 & - 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 355 END DO ! ji 356 370 357 !---------------- 371 358 ! Age of new ice … … 375 362 END DO ! ji 376 363 377 !--------------------------378 ! Open water energy budget379 !--------------------------380 DO ji = 1, nbpac381 zqbgow(ji) = qldif_1d(ji) - qcmif_1d(ji) !<0382 END DO ! ji383 384 364 !------------------- 385 365 ! Volume of new ice 386 366 !------------------- 387 367 DO ji = 1, nbpac 388 zv_newice(ji) = - zqbgow(ji) / ze_newice(ji) 368 369 zEi = - ze_newice(ji) / rhoic ! specific enthalpy of forming ice [J/kg] 370 371 zEw = rcp * ( t_bo_1d(ji) - rt0 ) ! specific enthalpy of seawater at t_bo_1d [J/kg] 372 ! clem: we suppose we are already at the freezing point (condition qlead<0 is satisfyied) 373 374 zdE = zEi - zEw ! specific enthalpy difference [J/kg] 375 376 zfmdt = - qlead_1d(ji) / zdE ! Fm.dt [kg/m2] (<0) 377 ! clem: we use qlead instead of zqld (limthd) because we suppose we are at the freezing point 378 zv_newice(ji) = - zfmdt / rhoic 379 380 zQm = zfmdt * zEw ! heat to the ocean >0 associated with mass flux 381 382 ! Contribution to heat flux to the ocean [W.m-2], >0 383 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * zEw * r1_rdtice 384 ! Total heat flux used in this process [W.m-2] 385 hfx_opw_1d(ji) = hfx_opw_1d(ji) - zfmdt * zdE * r1_rdtice 386 ! mass flux 387 wfx_opw_1d(ji) = wfx_opw_1d(ji) - zv_newice(ji) * rhoic * r1_rdtice 388 ! salt flux 389 sfx_opw_1d(ji) = sfx_opw_1d(ji) - zv_newice(ji) * rhoic * zs_newice(ji) * r1_rdtice 389 390 390 391 ! 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) 392 zinda = 1._wp - MAX( 0._wp, SIGN( 1._wp , - zat_i_1d(ji) ) ) 393 zfrazb = zinda * ( TANH ( Cfrazb * ( zvrel_1d(ji) - vfrazb ) ) + 1.0 ) * 0.5 * maxfrazb 394 zv_frazb(ji) = zfrazb * zv_newice(ji) 393 395 zv_newice(ji) = ( 1.0 - zfrazb ) * zv_newice(ji) 394 396 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 397 413 398 !----------------- … … 415 400 !----------------- 416 401 DO ji = 1, nbpac 417 ii = MOD( npac(ji) - 1 , jpi ) + 1418 ij = ( npac(ji) - 1 ) / jpi + 1419 402 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 403 END DO 422 404 423 405 !------------------------------------------------------------------------------! … … 425 407 !------------------------------------------------------------------------------! 426 408 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 !------------------------------------------- 409 !------------------------ 410 ! 6.1) lateral ice growth 411 !------------------------ 436 412 ! If lateral ice growth gives an ice concentration gt 1, then 437 413 ! we keep the excessive volume in memory and attribute it later to bottom accretion 438 414 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) )415 IF ( za_newice(ji) > ( amax - zat_i_1d(ji) ) ) THEN 416 zda_res(ji) = za_newice(ji) - ( amax - zat_i_1d(ji) ) 441 417 zdv_res(ji) = zda_res (ji) * zh_newice(ji) 442 418 za_newice(ji) = za_newice(ji) - zda_res (ji) … … 446 422 zdv_res(ji) = 0._wp 447 423 ENDIF 448 END DO ! ji 449 450 !------------------------------------------------ 451 ! Laterally redistribute new ice volume and area 452 !------------------------------------------------ 453 zat_i_ac(:) = 0._wp 424 END DO 425 426 ! find which category to fill 427 zat_i_1d(:) = 0._wp 454 428 DO jl = 1, jpl 455 429 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 430 IF( zh_newice(ji) > hi_max(jl-1) .AND. zh_newice(ji) <= hi_max(jl) ) THEN 431 za_i_1d (ji,jl) = za_i_1d (ji,jl) + za_newice(ji) 432 zv_i_1d (ji,jl) = zv_i_1d (ji,jl) + zv_newice(ji) 433 jcat (ji) = jl 462 434 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 435 zat_i_1d(ji) = zat_i_1d(ji) + za_i_1d (ji,jl) 436 END DO 437 END DO 438 439 ! Heat content 440 DO ji = 1, nbpac 441 jl = jcat(ji) ! categroy in which new ice is put 442 zswinew (ji) = MAX( 0._wp , SIGN( 1._wp , - za_b(ji,jl) ) ) ! 0 if old ice 475 443 END DO 476 444 477 445 DO jk = 1, nlay_i 478 446 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) 447 jl = jcat(ji) 448 zinda = MAX( 0._wp, SIGN( 1._wp , zv_i_1d(ji,jl) - epsi20 ) ) 449 ze_i_1d(ji,jk,jl) = zswinew(ji) * ze_newice(ji) + & 450 & ( 1.0 - zswinew(ji) ) * ( ze_newice(ji) * zv_newice(ji) + ze_i_1d(ji,jk,jl) * zv_b(ji,jl) ) & 451 & * zinda / MAX( zv_i_1d(ji,jl), epsi20 ) 452 END DO 453 END DO 454 455 !------------------------------------------------ 456 ! 6.2) bottom ice growth + ice enthalpy remapping 457 !------------------------------------------------ 458 DO jl = 1, jpl 459 460 ! for remapping 461 h_i_old (1:nbpac,0:nlay_i+1) = 0._wp 462 qh_i_old(1:nbpac,0:nlay_i+1) = 0._wp 533 463 DO jk = 1, nlay_i 534 464 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)465 h_i_old (ji,jk) = zv_i_1d(ji,jl) / REAL( nlay_i ) 466 qh_i_old(ji,jk) = ze_i_1d(ji,jk,jl) * h_i_old(ji,jk) 537 467 END DO 538 468 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) 469 470 ! new volumes including lateral/bottom accretion + residual 542 471 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 472 zinda = MAX( 0._wp, SIGN( 1._wp , zat_i_1d(ji) - epsi20 ) ) 473 zv_newfra = zinda * ( zdv_res(ji) + zv_frazb(ji) ) * za_i_1d(ji,jl) / MAX( zat_i_1d(ji) , epsi20 ) 474 za_i_1d(ji,jl) = zinda * za_i_1d(ji,jl) 475 zv_i_1d(ji,jl) = zv_i_1d(ji,jl) + zv_newfra 476 ! for remapping 477 h_i_old (ji,nlay_i+1) = zv_newfra 478 qh_i_old(ji,nlay_i+1) = ze_newice(ji) * zv_newfra 479 ENDDO 480 481 ! --- Ice enthalpy remapping --- ! 482 CALL lim_thd_ent( 1, nbpac, ze_i_1d(1:nbpac,:,jl) ) 483 ENDDO 574 484 575 485 !------------ … … 578 488 DO jl = 1, jpl 579 489 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 ) * zindb490 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_1d(ji,jl) + epsi20 ) ) ! 0 if no ice and 1 if yes 491 zoa_i_1d(ji,jl) = za_b(ji,jl) * zoa_i_1d(ji,jl) / MAX( za_i_1d(ji,jl) , epsi20 ) * zindb 582 492 END DO 583 493 END DO … … 586 496 ! Update salinity 587 497 !----------------- 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 498 DO jl = 1, jpl 602 499 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 500 zdv = zv_i_1d(ji,jl) - zv_b(ji,jl) 501 zsmv_i_1d(ji,jl) = zsmv_i_1d(ji,jl) + zdv * zs_newice(ji) 502 END DO 608 503 END DO 609 504 610 505 !------------------------------------------------------------------------------! 611 ! 8) Change 2D vectors to 1D vectors506 ! 7) Change 2D vectors to 1D vectors 612 507 !------------------------------------------------------------------------------! 613 508 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 ) 509 CALL tab_1d_2d( nbpac, a_i (:,:,jl), npac(1:nbpac), za_i_1d (1:nbpac,jl), jpi, jpj ) 510 CALL tab_1d_2d( nbpac, v_i (:,:,jl), npac(1:nbpac), zv_i_1d (1:nbpac,jl), jpi, jpj ) 511 CALL tab_1d_2d( nbpac, oa_i(:,:,jl), npac(1:nbpac), zoa_i_1d(1:nbpac,jl), jpi, jpj ) 512 CALL tab_1d_2d( nbpac, smv_i (:,:,jl), npac(1:nbpac), zsmv_i_1d(1:nbpac,jl) , jpi, jpj ) 619 513 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 ) 514 CALL tab_1d_2d( nbpac, e_i(:,:,jk,jl), npac(1:nbpac), ze_i_1d(1:nbpac,jk,jl), jpi, jpj ) 515 END DO 516 END DO 517 CALL tab_1d_2d( nbpac, sfx_opw, npac(1:nbpac), sfx_opw_1d(1:nbpac), jpi, jpj ) 518 CALL tab_1d_2d( nbpac, wfx_opw, npac(1:nbpac), wfx_opw_1d(1:nbpac), jpi, jpj ) 519 520 CALL tab_1d_2d( nbpac, hfx_thd, npac(1:nbpac), hfx_thd_1d(1:nbpac), jpi, jpj ) 521 CALL tab_1d_2d( nbpac, hfx_opw, npac(1:nbpac), hfx_opw_1d(1:nbpac), jpi, jpj ) 625 522 ! 626 523 ENDIF ! nbpac > 0 627 524 628 525 !------------------------------------------------------------------------------! 629 ! 9) Change units for e_i526 ! 8) Change units for e_i 630 527 !------------------------------------------------------------------------------! 631 528 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 529 DO jk = 1, nlay_i 530 DO jj = 1, jpj 531 DO ji = 1, jpi 532 ! heat content in Joules 533 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * area(ji,jj) * v_i(ji,jj,jl) / ( REAL( nlay_i ,wp ) * unit_fac ) 534 END DO 535 END DO 634 536 END DO 635 537 END DO 636 538 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 539 ! 669 CALL wrk_dealloc( jpij, zcatac) ! integer540 CALL wrk_dealloc( jpij, jcat ) ! integer 670 541 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 ) 542 CALL wrk_dealloc( jpij, zdv_res, zda_res, zat_i_1d, zv_frazb, zvrel_1d ) 543 CALL wrk_dealloc( jpij,jpl, zv_b, za_b, za_i_1d, zv_i_1d, zoa_i_1d, zsmv_i_1d ) 544 CALL wrk_dealloc( jpij,nlay_i+1,jpl, ze_i_1d ) 545 CALL wrk_dealloc( jpi,jpj, zvrel ) 676 546 ! 677 547 END SUBROUTINE lim_thd_lac
Note: See TracChangeset
for help on using the changeset viewer.
