Changeset 5038 for branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
- Timestamp:
- 2015-01-20T15:26:13+01:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
r4333 r5038 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 … … 35 37 36 38 PUBLIC lim_thd_lac ! called by lim_thd 37 38 REAL(wp) :: epsi10 = 1.e-10_wp !39 REAL(wp) :: zzero = 0._wp !40 REAL(wp) :: zone = 1._wp !41 39 42 40 !!---------------------------------------------------------------------- … … 71 69 !! - Computation of variation of ice volume and mass 72 70 !! - Computation of frldb after lateral accretion and 73 !! update ht_s_ b, ht_i_band tbif_1d(:,:)71 !! update ht_s_1d, ht_i_1d and tbif_1d(:,:) 74 72 !!------------------------------------------------------------------------ 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 scalars73 INTEGER :: ji,jj,jk,jl ! dummy loop indices 74 INTEGER :: nbpac ! local integers 75 INTEGER :: ii, ij, iter ! - - 76 REAL(wp) :: ztmelts, zdv, zfrazb, zweight, zde ! local scalars 79 77 REAL(wp) :: zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf , zhicol_new ! - - 80 78 REAL(wp) :: ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, zsqcd , zhicrit ! - - 81 79 LOGICAL :: iterate_frazil ! iterate frazil ice collection thickness 82 80 CHARACTER (len = 15) :: fieldid 83 ! 84 INTEGER , POINTER, DIMENSION(:) :: zcatac ! indexes of categories where new ice grows 81 82 REAL(wp) :: zQm ! enthalpy exchanged with the ocean (J/m2, >0 towards ocean) 83 REAL(wp) :: zEi ! sea ice specific enthalpy (J/kg) 84 REAL(wp) :: zEw ! seawater specific enthalpy (J/kg) 85 REAL(wp) :: zfmdt ! mass flux x time step (kg/m2, >0 towards ocean) 86 87 REAL(wp) :: zv_newfra 88 89 INTEGER , POINTER, DIMENSION(:) :: jcat ! indexes of categories where new ice grows 85 90 REAL(wp), POINTER, DIMENSION(:) :: zswinew ! switch for new ice or not 86 91 … … 93 98 REAL(wp), POINTER, DIMENSION(:) :: zdv_res ! residual volume in case of excessive heat budget 94 99 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 100 REAL(wp), POINTER, DIMENSION(:) :: zat_i_1d ! total ice fraction 101 REAL(wp), POINTER, DIMENSION(:) :: zv_frazb ! accretion of frazil ice at the ice bottom 102 REAL(wp), POINTER, DIMENSION(:) :: zvrel_1d ! relative ice / frazil velocity (1D vector) 103 104 REAL(wp), POINTER, DIMENSION(:,:) :: zv_b ! old volume of ice in category jl 105 REAL(wp), POINTER, DIMENSION(:,:) :: za_b ! old area of ice in category jl 106 REAL(wp), POINTER, DIMENSION(:,:) :: za_i_1d ! 1-D version of a_i 107 REAL(wp), POINTER, DIMENSION(:,:) :: zv_i_1d ! 1-D version of v_i 108 REAL(wp), POINTER, DIMENSION(:,:) :: zoa_i_1d ! 1-D version of oa_i 109 REAL(wp), POINTER, DIMENSION(:,:) :: zsmv_i_1d ! 1-D version of smv_i 110 111 REAL(wp), POINTER, DIMENSION(:,:,:) :: ze_i_1d !: 1-D version of e_i 112 122 113 REAL(wp), POINTER, DIMENSION(:,:) :: zvrel ! relative ice / frazil velocity 123 114 !!-----------------------------------------------------------------------! 124 115 125 CALL wrk_alloc( jpij, zcatac) ! integer116 CALL wrk_alloc( jpij, jcat ) ! integer 126 117 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 118 CALL wrk_alloc( jpij, zdv_res, zda_res, zat_i_1d, zv_frazb, zvrel_1d ) 119 CALL wrk_alloc( jpij,jpl, zv_b, za_b, za_i_1d, zv_i_1d, zoa_i_1d, zsmv_i_1d ) 120 CALL wrk_alloc( jpij,nlay_i+1,jpl, ze_i_1d ) 121 CALL wrk_alloc( jpi,jpj, zvrel ) 147 122 148 123 !------------------------------------------------------------------------------| … … 154 129 DO ji = 1, jpi 155 130 !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 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 131 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , -v_i(ji,jj,jl) + epsi10 ) ) !0 if no ice and 1 if yes 132 e_i(ji,jj,jk,jl) = rswitch * e_i(ji,jj,jk,jl) & 133 & / ( area(ji,jj) * MAX( v_i(ji,jj,jl) , epsi10 ) ) * REAL( nlay_i, wp ) 134 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac 159 135 END DO 160 136 END DO … … 179 155 180 156 ! Default new ice thickness 181 hicol(:,:) = hiccrit (1)182 183 IF( fraz_swi == 1 ._wp) THEN157 hicol(:,:) = hiccrit 158 159 IF( fraz_swi == 1 ) THEN 184 160 185 161 !-------------------- … … 193 169 zgamafr = 0.03 194 170 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 171 DO jj = 2, jpj 172 DO ji = 2, jpi 173 IF ( qlead(ji,jj) < 0._wp ) THEN 199 174 !------------- 200 175 ! Wind stress … … 206 181 & + vtau_ice(ji ,jj ) * tmv(ji ,jj ) ) * 0.5_wp 207 182 ! Square root of wind stress 208 ztenagm = SQRT( SQRT( ztaux * ztaux + ztauy * ztauy) )183 ztenagm = SQRT( SQRT( ztaux**2 + ztauy**2 ) ) 209 184 210 185 !--------------------- 211 186 ! Frazil ice velocity 212 187 !--------------------- 213 zvfrx = zgamafr * zsqcd * ztaux / MAX(ztenagm,epsi10) 214 zvfry = zgamafr * zsqcd * ztauy / MAX(ztenagm,epsi10) 188 rswitch = MAX( 0._wp, SIGN( 1._wp , ztenagm - epsi10 ) ) 189 zvfrx = rswitch * zgamafr * zsqcd * ztaux / MAX( ztenagm, epsi10 ) 190 zvfry = rswitch * zgamafr * zsqcd * ztauy / MAX( ztenagm, epsi10 ) 215 191 216 192 !------------------- … … 218 194 !------------------- 219 195 ! C-grid ice velocity 220 zindb = MAX( 0._wp, SIGN( 1._wp , at_i(ji,jj) ) ) 221 zvgx = zindb * ( u_ice(ji-1,jj ) * tmu(ji-1,jj ) & 222 & + u_ice(ji,jj ) * tmu(ji ,jj ) ) * 0.5_wp 223 zvgy = zindb * ( v_ice(ji ,jj-1) * tmv(ji ,jj-1) & 224 & + v_ice(ji,jj ) * tmv(ji ,jj ) ) * 0.5_wp 196 rswitch = MAX( 0._wp, SIGN( 1._wp , at_i(ji,jj) ) ) 197 zvgx = rswitch * ( u_ice(ji-1,jj ) * tmu(ji-1,jj ) + u_ice(ji,jj) * tmu(ji,jj) ) * 0.5_wp 198 zvgy = rswitch * ( v_ice(ji ,jj-1) * tmv(ji ,jj-1) + v_ice(ji,jj) * tmv(ji,jj) ) * 0.5_wp 225 199 226 200 !----------------------------------- … … 264 238 END DO ! loop on ji ends 265 239 END DO ! loop on jj ends 240 ! 241 CALL lbc_lnk( zvrel(:,:), 'T', 1. ) 242 CALL lbc_lnk( hicol(:,:), 'T', 1. ) 266 243 267 244 ENDIF ! End of computation of frazil ice collection thickness … … 276 253 ! This occurs if open water energy budget is negative 277 254 nbpac = 0 255 npac(:) = 0 256 ! 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) ) 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) ) 299 CALL tab_2d_1d( nbpac, qlead_1d (1:nbpac) , qlead , jpi, jpj, npac(1:nbpac) ) 300 CALL tab_2d_1d( nbpac, t_bo_1d (1:nbpac) , t_bo , 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, hicol_1d (1:nbpac) , hicol , jpi, jpj, npac(1:nbpac) ) 304 CALL tab_2d_1d( nbpac, zvrel_1d (1:nbpac) , zvrel , jpi, jpj, npac(1:nbpac) ) 305 306 CALL tab_2d_1d( nbpac, hfx_thd_1d(1:nbpac) , hfx_thd, jpi, jpj, npac(1:nbpac) ) 307 CALL tab_2d_1d( nbpac, hfx_opw_1d(1:nbpac) , hfx_opw, jpi, jpj, npac(1:nbpac) ) 327 308 328 309 !------------------------------------------------------------------------------! … … 330 311 !------------------------------------------------------------------------------! 331 312 313 !----------------------------------------- 314 ! Keep old ice areas and volume in memory 315 !----------------------------------------- 316 zv_b(1:nbpac,:) = zv_i_1d(1:nbpac,:) 317 za_b(1:nbpac,:) = za_i_1d(1:nbpac,:) 332 318 !---------------------- 333 319 ! Thickness of new ice 334 320 !---------------------- 335 321 DO ji = 1, nbpac 336 zh_newice(ji) = hiccrit (1)337 END DO 338 IF( fraz_swi == 1 .0 ) zh_newice(:) = hicol_b(:)322 zh_newice(ji) = hiccrit 323 END DO 324 IF( fraz_swi == 1 ) zh_newice(1:nbpac) = hicol_1d(1:nbpac) 339 325 340 326 !---------------------- 341 327 ! Salinity of new ice 342 328 !---------------------- 343 344 329 SELECT CASE ( num_sal ) 345 330 CASE ( 1 ) ! Sice = constant 346 zs_newice( :) = bulk_sal331 zs_newice(1:nbpac) = bulk_sal 347 332 CASE ( 2 ) ! Sice = F(z,t) [Vancoppenolle et al (2005)] 348 333 DO ji = 1, nbpac … … 352 337 END DO 353 338 CASE ( 3 ) ! Sice = F(z) [multiyear ice] 354 zs_newice( :) = 2.3339 zs_newice(1:nbpac) = 2.3 355 340 END SELECT 356 357 341 358 342 !------------------------- … … 362 346 DO ji = 1, nbpac 363 347 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) ) &348 ze_newice(ji) = rhoic * ( cpic * ( ztmelts - t_bo_1d(ji) ) & 349 & + lfus * ( 1.0 - ( ztmelts - rtt ) / MIN( t_bo_1d(ji) - rtt, -epsi10 ) ) & 366 350 & - 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 351 END DO ! ji 352 370 353 !---------------- 371 354 ! Age of new ice … … 375 358 END DO ! ji 376 359 377 !--------------------------378 ! Open water energy budget379 !--------------------------380 DO ji = 1, nbpac381 zqbgow(ji) = qldif_1d(ji) - qcmif_1d(ji) !<0382 END DO ! ji383 384 360 !------------------- 385 361 ! Volume of new ice 386 362 !------------------- 387 363 DO ji = 1, nbpac 388 zv_newice(ji) = - zqbgow(ji) / ze_newice(ji) 364 365 zEi = - ze_newice(ji) / rhoic ! specific enthalpy of forming ice [J/kg] 366 367 zEw = rcp * ( t_bo_1d(ji) - rt0 ) ! specific enthalpy of seawater at t_bo_1d [J/kg] 368 ! clem: we suppose we are already at the freezing point (condition qlead<0 is satisfyied) 369 370 zdE = zEi - zEw ! specific enthalpy difference [J/kg] 371 372 zfmdt = - qlead_1d(ji) / zdE ! Fm.dt [kg/m2] (<0) 373 ! clem: we use qlead instead of zqld (limthd) because we suppose we are at the freezing point 374 zv_newice(ji) = - zfmdt / rhoic 375 376 zQm = zfmdt * zEw ! heat to the ocean >0 associated with mass flux 377 378 ! Contribution to heat flux to the ocean [W.m-2], >0 379 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * zEw * r1_rdtice 380 ! Total heat flux used in this process [W.m-2] 381 hfx_opw_1d(ji) = hfx_opw_1d(ji) - zfmdt * zdE * r1_rdtice 382 ! mass flux 383 wfx_opw_1d(ji) = wfx_opw_1d(ji) - zv_newice(ji) * rhoic * r1_rdtice 384 ! salt flux 385 sfx_opw_1d(ji) = sfx_opw_1d(ji) - zv_newice(ji) * rhoic * zs_newice(ji) * r1_rdtice 389 386 390 387 ! 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) 388 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp , - zat_i_1d(ji) ) ) 389 zfrazb = rswitch * ( TANH ( Cfrazb * ( zvrel_1d(ji) - vfrazb ) ) + 1.0 ) * 0.5 * maxfrazb 390 zv_frazb(ji) = zfrazb * zv_newice(ji) 393 391 zv_newice(ji) = ( 1.0 - zfrazb ) * zv_newice(ji) 394 392 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 393 413 394 !----------------- … … 415 396 !----------------- 416 397 DO ji = 1, nbpac 417 ii = MOD( npac(ji) - 1 , jpi ) + 1418 ij = ( npac(ji) - 1 ) / jpi + 1419 398 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 399 END DO 422 400 423 401 !------------------------------------------------------------------------------! … … 425 403 !------------------------------------------------------------------------------! 426 404 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 !------------------------------------------- 405 !------------------------ 406 ! 6.1) lateral ice growth 407 !------------------------ 436 408 ! If lateral ice growth gives an ice concentration gt 1, then 437 409 ! we keep the excessive volume in memory and attribute it later to bottom accretion 438 410 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) )411 IF ( za_newice(ji) > ( amax - zat_i_1d(ji) ) ) THEN 412 zda_res(ji) = za_newice(ji) - ( amax - zat_i_1d(ji) ) 441 413 zdv_res(ji) = zda_res (ji) * zh_newice(ji) 442 414 za_newice(ji) = za_newice(ji) - zda_res (ji) … … 446 418 zdv_res(ji) = 0._wp 447 419 ENDIF 448 END DO ! ji 449 450 !------------------------------------------------ 451 ! Laterally redistribute new ice volume and area 452 !------------------------------------------------ 453 zat_i_ac(:) = 0._wp 420 END DO 421 422 ! find which category to fill 423 zat_i_1d(:) = 0._wp 454 424 DO jl = 1, jpl 455 425 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 426 IF( zh_newice(ji) > hi_max(jl-1) .AND. zh_newice(ji) <= hi_max(jl) ) THEN 427 za_i_1d (ji,jl) = za_i_1d (ji,jl) + za_newice(ji) 428 zv_i_1d (ji,jl) = zv_i_1d (ji,jl) + zv_newice(ji) 429 jcat (ji) = jl 462 430 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 431 zat_i_1d(ji) = zat_i_1d(ji) + za_i_1d (ji,jl) 432 END DO 433 END DO 434 435 ! Heat content 436 DO ji = 1, nbpac 437 jl = jcat(ji) ! categroy in which new ice is put 438 zswinew (ji) = MAX( 0._wp , SIGN( 1._wp , - za_b(ji,jl) ) ) ! 0 if old ice 475 439 END DO 476 440 477 441 DO jk = 1, nlay_i 478 442 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) 443 jl = jcat(ji) 444 rswitch = MAX( 0._wp, SIGN( 1._wp , zv_i_1d(ji,jl) - epsi20 ) ) 445 ze_i_1d(ji,jk,jl) = zswinew(ji) * ze_newice(ji) + & 446 & ( 1.0 - zswinew(ji) ) * ( ze_newice(ji) * zv_newice(ji) + ze_i_1d(ji,jk,jl) * zv_b(ji,jl) ) & 447 & * rswitch / MAX( zv_i_1d(ji,jl), epsi20 ) 448 END DO 449 END DO 450 451 !------------------------------------------------ 452 ! 6.2) bottom ice growth + ice enthalpy remapping 453 !------------------------------------------------ 454 DO jl = 1, jpl 455 456 ! for remapping 457 h_i_old (1:nbpac,0:nlay_i+1) = 0._wp 458 qh_i_old(1:nbpac,0:nlay_i+1) = 0._wp 533 459 DO jk = 1, nlay_i 534 460 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)461 h_i_old (ji,jk) = zv_i_1d(ji,jl) / REAL( nlay_i ) 462 qh_i_old(ji,jk) = ze_i_1d(ji,jk,jl) * h_i_old(ji,jk) 537 463 END DO 538 464 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) 465 466 ! new volumes including lateral/bottom accretion + residual 542 467 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 468 rswitch = MAX( 0._wp, SIGN( 1._wp , zat_i_1d(ji) - epsi20 ) ) 469 zv_newfra = rswitch * ( zdv_res(ji) + zv_frazb(ji) ) * za_i_1d(ji,jl) / MAX( zat_i_1d(ji) , epsi20 ) 470 za_i_1d(ji,jl) = rswitch * za_i_1d(ji,jl) 471 zv_i_1d(ji,jl) = zv_i_1d(ji,jl) + zv_newfra 472 ! for remapping 473 h_i_old (ji,nlay_i+1) = zv_newfra 474 qh_i_old(ji,nlay_i+1) = ze_newice(ji) * zv_newfra 475 ENDDO 476 ! --- Ice enthalpy remapping --- ! 477 CALL lim_thd_ent( 1, nbpac, ze_i_1d(1:nbpac,:,jl) ) 478 ENDDO 574 479 575 480 !------------ … … 578 483 DO jl = 1, jpl 579 484 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 ) * zindb485 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_1d(ji,jl) + epsi20 ) ) ! 0 if no ice and 1 if yes 486 zoa_i_1d(ji,jl) = za_b(ji,jl) * zoa_i_1d(ji,jl) / MAX( za_i_1d(ji,jl) , epsi20 ) * rswitch 582 487 END DO 583 488 END DO … … 586 491 ! Update salinity 587 492 !----------------- 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 493 DO jl = 1, jpl 602 494 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 495 zdv = zv_i_1d(ji,jl) - zv_b(ji,jl) 496 zsmv_i_1d(ji,jl) = zsmv_i_1d(ji,jl) + zdv * zs_newice(ji) 497 END DO 608 498 END DO 609 499 610 500 !------------------------------------------------------------------------------! 611 ! 8) Change 2D vectors to 1D vectors501 ! 7) Change 2D vectors to 1D vectors 612 502 !------------------------------------------------------------------------------! 613 503 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 ) 504 CALL tab_1d_2d( nbpac, a_i (:,:,jl), npac(1:nbpac), za_i_1d (1:nbpac,jl), jpi, jpj ) 505 CALL tab_1d_2d( nbpac, v_i (:,:,jl), npac(1:nbpac), zv_i_1d (1:nbpac,jl), jpi, jpj ) 506 CALL tab_1d_2d( nbpac, oa_i(:,:,jl), npac(1:nbpac), zoa_i_1d(1:nbpac,jl), jpi, jpj ) 507 CALL tab_1d_2d( nbpac, smv_i (:,:,jl), npac(1:nbpac), zsmv_i_1d(1:nbpac,jl) , jpi, jpj ) 619 508 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 ) 509 CALL tab_1d_2d( nbpac, e_i(:,:,jk,jl), npac(1:nbpac), ze_i_1d(1:nbpac,jk,jl), jpi, jpj ) 510 END DO 511 END DO 512 CALL tab_1d_2d( nbpac, sfx_opw, npac(1:nbpac), sfx_opw_1d(1:nbpac), jpi, jpj ) 513 CALL tab_1d_2d( nbpac, wfx_opw, npac(1:nbpac), wfx_opw_1d(1:nbpac), jpi, jpj ) 514 515 CALL tab_1d_2d( nbpac, hfx_thd, npac(1:nbpac), hfx_thd_1d(1:nbpac), jpi, jpj ) 516 CALL tab_1d_2d( nbpac, hfx_opw, npac(1:nbpac), hfx_opw_1d(1:nbpac), jpi, jpj ) 625 517 ! 626 518 ENDIF ! nbpac > 0 627 519 628 520 !------------------------------------------------------------------------------! 629 ! 9) Change units for e_i521 ! 8) Change units for e_i 630 522 !------------------------------------------------------------------------------! 631 523 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 524 DO jk = 1, nlay_i 525 DO jj = 1, jpj 526 DO ji = 1, jpi 527 ! heat content in Joules 528 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 ) 529 END DO 530 END DO 634 531 END DO 635 532 END DO 636 533 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 534 ! 669 CALL wrk_dealloc( jpij, zcatac) ! integer535 CALL wrk_dealloc( jpij, jcat ) ! integer 670 536 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 ) 537 CALL wrk_dealloc( jpij, zdv_res, zda_res, zat_i_1d, zv_frazb, zvrel_1d ) 538 CALL wrk_dealloc( jpij,jpl, zv_b, za_b, za_i_1d, zv_i_1d, zoa_i_1d, zsmv_i_1d ) 539 CALL wrk_dealloc( jpij,nlay_i+1,jpl, ze_i_1d ) 540 CALL wrk_dealloc( jpi,jpj, zvrel ) 676 541 ! 677 542 END SUBROUTINE lim_thd_lac
Note: See TracChangeset
for help on using the changeset viewer.
