- Timestamp:
- 2017-07-13T11:29:29+02:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
r8325 r8327 71 71 !!------------------------------------------------------------------------ 72 72 INTEGER :: ji,jj,jk,jl ! dummy loop indices 73 INTEGER :: n bpac! local integers74 INTEGER :: i i, ij, iter ! - -73 INTEGER :: nidx ! local integers 74 INTEGER :: iter ! - - 75 75 REAL(wp) :: ztmelts, zdv, zfrazb, zweight, zde ! local scalars 76 76 REAL(wp) :: zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf ! - - 77 77 REAL(wp) :: ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, zsqcd , zhicrit ! - - 78 CHARACTER (len = 15) :: fieldid79 78 80 79 REAL(wp) :: zQm ! enthalpy exchanged with the ocean (J/m2, >0 towards ocean) … … 122 121 CALL lim_var_agg(1) 123 122 CALL lim_var_glo2eqv 124 !------------------------------------------------------------------------------|125 ! 2) Convert units for ice internal energy126 !------------------------------------------------------------------------------|127 DO jl = 1, jpl128 DO jk = 1, nlay_i129 DO jj = 1, jpj130 DO ji = 1, jpi131 !Energy of melting q(S,T) [J.m-3]132 rswitch = MAX( 0._wp , SIGN( 1._wp , v_i(ji,jj,jl) - epsi20 ) ) !0 if no ice133 e_i(ji,jj,jk,jl) = rswitch * e_i(ji,jj,jk,jl) / MAX( v_i(ji,jj,jl), epsi20 ) * REAL( nlay_i, wp )134 END DO135 END DO136 END DO137 END DO138 123 139 124 !------------------------------------------------------------------------------! … … 240 225 !------------------------------------- 241 226 ! This occurs if open water energy budget is negative (cooling) and there is no landfast ice 242 nbpac = 0 243 npac(:) = 0 244 ! 245 DO jj = 1, jpj 246 DO ji = 1, jpi 227 nidx = 0 ; idxice(:) = 0 228 DO jj = 2, jpjm1 229 DO ji = 2, jpim1 247 230 IF ( qlead(ji,jj) < 0._wp .AND. tau_icebfr(ji,jj) == 0._wp ) THEN 248 n bpac = nbpac+ 1249 npac( nbpac) = (jj - 1) * jpi + ji231 nidx = nidx + 1 232 idxice( nidx ) = (jj - 1) * jpi + ji 250 233 ENDIF 251 234 END DO … … 264 247 ENDIF 265 248 266 IF( ln_limctl ) WRITE(numout,*) 'lim_thd_lac : n bpac = ', nbpac249 IF( ln_limctl ) WRITE(numout,*) 'lim_thd_lac : nidx = ', nidx 267 250 268 251 !------------------------------ … … 271 254 ! If ocean gains heat do nothing. Otherwise compute new ice formation 272 255 273 IF ( n bpac> 0 ) THEN274 275 CALL tab_2d_1d( n bpac, zat_i_1d (1:nbpac) , at_i , jpi, jpj, npac(1:nbpac) )276 DO jl = 1, jpl 277 CALL tab_2d_1d( n bpac, za_i_1d (1:nbpac,jl), a_i (:,:,jl), jpi, jpj, npac(1:nbpac) )278 CALL tab_2d_1d( n bpac, zv_i_1d (1:nbpac,jl), v_i (:,:,jl), jpi, jpj, npac(1:nbpac) )279 CALL tab_2d_1d( n bpac, zsmv_i_1d(1:nbpac,jl), smv_i(:,:,jl), jpi, jpj, npac(1:nbpac) )256 IF ( nidx > 0 ) THEN 257 258 CALL tab_2d_1d( nidx, zat_i_1d (1:nidx) , at_i , jpi, jpj, idxice(1:nidx) ) 259 DO jl = 1, jpl 260 CALL tab_2d_1d( nidx, za_i_1d (1:nidx,jl), a_i (:,:,jl), jpi, jpj, idxice(1:nidx) ) 261 CALL tab_2d_1d( nidx, zv_i_1d (1:nidx,jl), v_i (:,:,jl), jpi, jpj, idxice(1:nidx) ) 262 CALL tab_2d_1d( nidx, zsmv_i_1d(1:nidx,jl), smv_i(:,:,jl), jpi, jpj, idxice(1:nidx) ) 280 263 DO jk = 1, nlay_i 281 CALL tab_2d_1d( nbpac, ze_i_1d(1:nbpac,jk,jl), e_i(:,:,jk,jl) , jpi, jpj, npac(1:nbpac) ) 282 END DO 283 END DO 284 285 CALL tab_2d_1d( nbpac, qlead_1d (1:nbpac) , qlead , jpi, jpj, npac(1:nbpac) ) 286 CALL tab_2d_1d( nbpac, t_bo_1d (1:nbpac) , t_bo , jpi, jpj, npac(1:nbpac) ) 287 CALL tab_2d_1d( nbpac, sfx_opw_1d(1:nbpac) , sfx_opw , jpi, jpj, npac(1:nbpac) ) 288 CALL tab_2d_1d( nbpac, wfx_opw_1d(1:nbpac) , wfx_opw , jpi, jpj, npac(1:nbpac) ) 289 CALL tab_2d_1d( nbpac, hicol_1d (1:nbpac) , hicol , jpi, jpj, npac(1:nbpac) ) 290 CALL tab_2d_1d( nbpac, zvrel_1d (1:nbpac) , zvrel , jpi, jpj, npac(1:nbpac) ) 291 292 CALL tab_2d_1d( nbpac, hfx_thd_1d(1:nbpac) , hfx_thd , jpi, jpj, npac(1:nbpac) ) 293 CALL tab_2d_1d( nbpac, hfx_opw_1d(1:nbpac) , hfx_opw , jpi, jpj, npac(1:nbpac) ) 294 CALL tab_2d_1d( nbpac, rn_amax_1d(1:nbpac) , rn_amax_2d, jpi, jpj, npac(1:nbpac) ) 295 264 CALL tab_2d_1d( nidx, ze_i_1d(1:nidx,jk,jl), e_i(:,:,jk,jl) , jpi, jpj, idxice(1:nidx) ) 265 END DO 266 END DO 267 268 CALL tab_2d_1d( nidx, qlead_1d (1:nidx) , qlead , jpi, jpj, idxice(1:nidx) ) 269 CALL tab_2d_1d( nidx, t_bo_1d (1:nidx) , t_bo , jpi, jpj, idxice(1:nidx) ) 270 CALL tab_2d_1d( nidx, sfx_opw_1d(1:nidx) , sfx_opw , jpi, jpj, idxice(1:nidx) ) 271 CALL tab_2d_1d( nidx, wfx_opw_1d(1:nidx) , wfx_opw , jpi, jpj, idxice(1:nidx) ) 272 CALL tab_2d_1d( nidx, hicol_1d (1:nidx) , hicol , jpi, jpj, idxice(1:nidx) ) 273 CALL tab_2d_1d( nidx, zvrel_1d (1:nidx) , zvrel , jpi, jpj, idxice(1:nidx) ) 274 275 CALL tab_2d_1d( nidx, hfx_thd_1d(1:nidx) , hfx_thd , jpi, jpj, idxice(1:nidx) ) 276 CALL tab_2d_1d( nidx, hfx_opw_1d(1:nidx) , hfx_opw , jpi, jpj, idxice(1:nidx) ) 277 CALL tab_2d_1d( nidx, rn_amax_1d(1:nidx) , rn_amax_2d, jpi, jpj, idxice(1:nidx) ) 278 CALL tab_2d_1d( nidx, sss_1d (1:nidx) , sss_m , jpi, jpj, idxice(1:nidx) ) 279 280 !------------------------------------------------------------------------------| 281 ! 2) Convert units for ice internal energy 282 !------------------------------------------------------------------------------| 283 DO jl = 1, jpl 284 DO jk = 1, nlay_i 285 DO ji = 1, nidx 286 IF( zv_i_1d(ji,jl) > 0._wp ) ze_i_1d(ji,jk,jl) = ze_i_1d(ji,jk,jl) / zv_i_1d(ji,jl) * REAL( nlay_i ) 287 END DO 288 END DO 289 END DO 296 290 !------------------------------------------------------------------------------! 297 291 ! 5) Compute thickness, salinity, enthalpy, age, area and volume of new ice … … 301 295 ! Keep old ice areas and volume in memory 302 296 !----------------------------------------- 303 zv_b(1:n bpac,:) = zv_i_1d(1:nbpac,:)304 za_b(1:n bpac,:) = za_i_1d(1:nbpac,:)297 zv_b(1:nidx,:) = zv_i_1d(1:nidx,:) 298 za_b(1:nidx,:) = za_i_1d(1:nidx,:) 305 299 306 300 !---------------------- 307 301 ! Thickness of new ice 308 302 !---------------------- 309 zh_newice(1:n bpac) = hicol_1d(1:nbpac)303 zh_newice(1:nidx) = hicol_1d(1:nidx) 310 304 311 305 !---------------------- … … 314 308 SELECT CASE ( nn_icesal ) 315 309 CASE ( 1 ) ! Sice = constant 316 zs_newice(1:n bpac) = rn_icesal310 zs_newice(1:nidx) = rn_icesal 317 311 CASE ( 2 ) ! Sice = F(z,t) [Vancoppenolle et al (2005)] 318 DO ji = 1, nbpac 319 ii = MOD( npac(ji) - 1 , jpi ) + 1 320 ij = ( npac(ji) - 1 ) / jpi + 1 321 zs_newice(ji) = MIN( 4.606 + 0.91 / zh_newice(ji) , rn_simax , 0.5 * sss_m(ii,ij) ) 312 DO ji = 1, nidx 313 zs_newice(ji) = MIN( 4.606 + 0.91 / zh_newice(ji) , rn_simax , 0.5 * sss_1d(ji) ) 322 314 END DO 323 315 CASE ( 3 ) ! Sice = F(z) [multiyear ice] 324 zs_newice(1:n bpac) = 2.3316 zs_newice(1:nidx) = 2.3 325 317 END SELECT 326 318 … … 329 321 !------------------------- 330 322 ! We assume that new ice is formed at the seawater freezing point 331 DO ji = 1, n bpac323 DO ji = 1, nidx 332 324 ztmelts = - tmut * zs_newice(ji) + rt0 ! Melting point (K) 333 325 ze_newice(ji) = rhoic * ( cpic * ( ztmelts - t_bo_1d(ji) ) & … … 339 331 ! Age of new ice 340 332 !---------------- 341 DO ji = 1, n bpac333 DO ji = 1, nidx 342 334 zo_newice(ji) = 0._wp 343 335 END DO … … 346 338 ! Volume of new ice 347 339 !------------------- 348 DO ji = 1, n bpac340 DO ji = 1, nidx 349 341 350 342 zEi = - ze_newice(ji) * r1_rhoic ! specific enthalpy of forming ice [J/kg] … … 374 366 IF( ln_frazil ) THEN 375 367 ! A fraction zfrazb of frazil ice is accreted at the ice bottom 376 DO ji = 1, n bpac368 DO ji = 1, nidx 377 369 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp , - zat_i_1d(ji) ) ) 378 370 zfrazb = rswitch * ( TANH( rn_Cfrazb * ( zvrel_1d(ji) - rn_vfrazb ) ) + 1.0 ) * 0.5 * rn_maxfrazb … … 385 377 ! Area of new ice 386 378 !----------------- 387 DO ji = 1, n bpac379 DO ji = 1, nidx 388 380 za_newice(ji) = zv_newice(ji) / zh_newice(ji) 389 381 END DO … … 398 390 ! If lateral ice growth gives an ice concentration gt 1, then 399 391 ! we keep the excessive volume in memory and attribute it later to bottom accretion 400 DO ji = 1, n bpac392 DO ji = 1, nidx 401 393 IF ( za_newice(ji) > ( rn_amax_1d(ji) - zat_i_1d(ji) ) ) THEN 402 394 zda_res(ji) = za_newice(ji) - ( rn_amax_1d(ji) - zat_i_1d(ji) ) … … 413 405 zat_i_1d(:) = 0._wp 414 406 DO jl = 1, jpl 415 DO ji = 1, n bpac407 DO ji = 1, nidx 416 408 IF( zh_newice(ji) > hi_max(jl-1) .AND. zh_newice(ji) <= hi_max(jl) ) THEN 417 409 za_i_1d (ji,jl) = za_i_1d (ji,jl) + za_newice(ji) … … 424 416 425 417 ! Heat content 426 DO ji = 1, n bpac418 DO ji = 1, nidx 427 419 jl = jcat(ji) ! categroy in which new ice is put 428 420 zswinew (ji) = MAX( 0._wp , SIGN( 1._wp , - za_b(ji,jl) ) ) ! 0 if old ice … … 430 422 431 423 DO jk = 1, nlay_i 432 DO ji = 1, n bpac424 DO ji = 1, nidx 433 425 jl = jcat(ji) 434 426 rswitch = MAX( 0._wp, SIGN( 1._wp , zv_i_1d(ji,jl) - epsi20 ) ) … … 445 437 446 438 ! for remapping 447 h_i_old (1:n bpac,0:nlay_i+1) = 0._wp448 eh_i_old(1:n bpac,0:nlay_i+1) = 0._wp439 h_i_old (1:nidx,0:nlay_i+1) = 0._wp 440 eh_i_old(1:nidx,0:nlay_i+1) = 0._wp 449 441 DO jk = 1, nlay_i 450 DO ji = 1, n bpac442 DO ji = 1, nidx 451 443 h_i_old (ji,jk) = zv_i_1d(ji,jl) * r1_nlay_i 452 444 eh_i_old(ji,jk) = ze_i_1d(ji,jk,jl) * h_i_old(ji,jk) … … 455 447 456 448 ! new volumes including lateral/bottom accretion + residual 457 DO ji = 1, n bpac449 DO ji = 1, nidx 458 450 rswitch = MAX( 0._wp, SIGN( 1._wp , zat_i_1d(ji) - epsi20 ) ) 459 451 zv_newfra = rswitch * ( zdv_res(ji) + zv_frazb(ji) ) * za_i_1d(ji,jl) / MAX( zat_i_1d(ji) , epsi20 ) … … 465 457 ENDDO 466 458 ! --- Ice enthalpy remapping --- ! 467 CALL lim_thd_ent( 1, n bpac, ze_i_1d(1:nbpac,:,jl) )459 CALL lim_thd_ent( 1, nidx, ze_i_1d(1:nidx,:,jl) ) 468 460 ENDDO 469 461 … … 472 464 !----------------- 473 465 DO jl = 1, jpl 474 DO ji = 1, n bpac466 DO ji = 1, nidx 475 467 zdv = zv_i_1d(ji,jl) - zv_b(ji,jl) 476 468 zsmv_i_1d(ji,jl) = zsmv_i_1d(ji,jl) + zdv * zs_newice(ji) … … 479 471 480 472 !------------------------------------------------------------------------------! 473 ! 8) Change units for e_i 474 !------------------------------------------------------------------------------! 475 DO jl = 1, jpl 476 DO jk = 1, nlay_i 477 DO ji = 1, nidx 478 ze_i_1d(ji,jk,jl) = ze_i_1d(ji,jk,jl) * zv_i_1d(ji,jl) * r1_nlay_i 479 END DO 480 END DO 481 END DO 482 !------------------------------------------------------------------------------! 481 483 ! 7) Change 2D vectors to 1D vectors 482 484 !------------------------------------------------------------------------------! 483 485 DO jl = 1, jpl 484 CALL tab_1d_2d( n bpac, a_i (:,:,jl), npac(1:nbpac), za_i_1d (1:nbpac,jl), jpi, jpj )485 CALL tab_1d_2d( n bpac, v_i (:,:,jl), npac(1:nbpac), zv_i_1d (1:nbpac,jl), jpi, jpj )486 CALL tab_1d_2d( n bpac, smv_i (:,:,jl), npac(1:nbpac), zsmv_i_1d(1:nbpac,jl) , jpi, jpj )486 CALL tab_1d_2d( nidx, a_i (:,:,jl), idxice(1:nidx), za_i_1d (1:nidx,jl), jpi, jpj ) 487 CALL tab_1d_2d( nidx, v_i (:,:,jl), idxice(1:nidx), zv_i_1d (1:nidx,jl), jpi, jpj ) 488 CALL tab_1d_2d( nidx, smv_i (:,:,jl), idxice(1:nidx), zsmv_i_1d(1:nidx,jl) , jpi, jpj ) 487 489 DO jk = 1, nlay_i 488 CALL tab_1d_2d( nbpac, e_i(:,:,jk,jl), npac(1:nbpac), ze_i_1d(1:nbpac,jk,jl), jpi, jpj ) 489 END DO 490 END DO 491 CALL tab_1d_2d( nbpac, sfx_opw, npac(1:nbpac), sfx_opw_1d(1:nbpac), jpi, jpj ) 492 CALL tab_1d_2d( nbpac, wfx_opw, npac(1:nbpac), wfx_opw_1d(1:nbpac), jpi, jpj ) 493 494 CALL tab_1d_2d( nbpac, hfx_thd, npac(1:nbpac), hfx_thd_1d(1:nbpac), jpi, jpj ) 495 CALL tab_1d_2d( nbpac, hfx_opw, npac(1:nbpac), hfx_opw_1d(1:nbpac), jpi, jpj ) 490 CALL tab_1d_2d( nidx, e_i(:,:,jk,jl), idxice(1:nidx), ze_i_1d(1:nidx,jk,jl), jpi, jpj ) 491 END DO 492 END DO 493 CALL tab_1d_2d( nidx, sfx_opw, idxice(1:nidx), sfx_opw_1d(1:nidx), jpi, jpj ) 494 CALL tab_1d_2d( nidx, wfx_opw, idxice(1:nidx), wfx_opw_1d(1:nidx), jpi, jpj ) 495 CALL tab_1d_2d( nidx, hfx_thd, idxice(1:nidx), hfx_thd_1d(1:nidx), jpi, jpj ) 496 CALL tab_1d_2d( nidx, hfx_opw, idxice(1:nidx), hfx_opw_1d(1:nidx), jpi, jpj ) 496 497 ! 497 ENDIF ! nbpac > 0 498 499 !------------------------------------------------------------------------------! 500 ! 8) Change units for e_i 501 !------------------------------------------------------------------------------! 502 DO jl = 1, jpl 503 DO jk = 1, nlay_i 504 DO jj = 1, jpj 505 DO ji = 1, jpi 506 ! heat content in J/m2 507 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i 508 END DO 509 END DO 510 END DO 511 END DO 512 498 ENDIF ! nidx > 0 513 499 ! 514 500 CALL wrk_dealloc( jpij, jcat ) ! integer
Note: See TracChangeset
for help on using the changeset viewer.
