- Timestamp:
- 2013-11-07T11:01:27+01:00 (10 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
r3625 r4161 46 46 47 47 !!---------------------------------------------------------------------- 48 !! NEMO/LIM3 3.4, UCL - NEMO Consortium (2011)48 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 49 49 !! $Id$ 50 50 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 78 78 !! update ht_s_b, ht_i_b and tbif_1d(:,:) 79 79 !!------------------------------------------------------------------------ 80 INTEGER 81 INTEGER 82 INTEGER :: zji, zjj, iter ! - -83 REAL(wp) :: ztmelts, zdv, zqold, zfrazb, zweight, zalphai, zindb, zde! local scalars80 INTEGER :: ji,jj,jk,jl,jm ! dummy loop indices 81 INTEGER :: layer, nbpac ! local integers 82 INTEGER :: ii, ij, iter ! - - 83 REAL(wp) :: ztmelts, zdv, zqold, zfrazb, zweight, zalphai, zindb, zinda, zde ! local scalars 84 84 REAL(wp) :: zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf , zhicol_new ! - - 85 85 REAL(wp) :: ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, zsqcd , zhicrit ! - - 86 REAL(wp) :: zcoef ! - -87 86 LOGICAL :: iterate_frazil ! iterate frazil ice collection thickness 88 87 CHARACTER (len = 15) :: fieldid … … 160 159 DO ji = 1, jpi 161 160 !Energy of melting q(S,T) [J.m-3] 162 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / MAX( area(ji,jj) * v_i(ji,jj,jl) , epsi10 ) * nlay_i161 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 ) 163 162 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , -v_i(ji,jj,jl) ) ) !0 if no ice and 1 if yes 164 163 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac * zindb … … 342 341 CASE ( 2 ) ! Sice = F(z,t) [Vancoppenolle et al (2005)] 343 342 DO ji = 1, nbpac 344 zji = MOD( npac(ji) - 1 , jpi ) + 1345 zjj = ( npac(ji) - 1 ) / jpi + 1346 zs_newice(ji) = MIN( 4.606 + 0.91 / zh_newice(ji) , s_i_max , 0.5 * sss_m( zji,zjj) )343 ii = MOD( npac(ji) - 1 , jpi ) + 1 344 ij = ( npac(ji) - 1 ) / jpi + 1 345 zs_newice(ji) = MIN( 4.606 + 0.91 / zh_newice(ji) , s_i_max , 0.5 * sss_m(ii,ij) ) 347 346 END DO 348 347 CASE ( 3 ) ! Sice = F(z) [multiyear ice] … … 389 388 END DO 390 389 391 !---------------------------------392 ! Salt flux due to new ice growth393 !---------------------------------394 ! note that for constant salinity zs_newice() = bulk_sal (see top of the subroutine)395 DO ji = 1, nbpac396 sfx_thd_1d(ji) = sfx_thd_1d(ji) - zs_newice(ji) * rhoic * zv_newice(ji) * r1_rdtice397 rdm_ice_1d(ji) = rdm_ice_1d(ji) + rhoic * zv_newice(ji)398 END DO ! ji399 400 390 !------------------------------------ 401 391 ! Diags for energy conservation test 402 392 !------------------------------------ 403 393 DO ji = 1, nbpac 404 zji = MOD( npac(ji) - 1 , jpi ) + 1405 zjj = ( npac(ji) - 1 ) / jpi + 1394 ii = MOD( npac(ji) - 1 , jpi ) + 1 395 ij = ( npac(ji) - 1 ) / jpi + 1 406 396 ! 407 zde = ze_newice(ji) / unit_fac * area( zji,zjj) * zv_newice(ji)397 zde = ze_newice(ji) / unit_fac * area(ii,ij) * zv_newice(ji) 408 398 ! 409 vt_i_init( zji,zjj) = vt_i_init(zji,zjj) + zv_newice(ji) ! volume410 et_i_init( zji,zjj) = et_i_init(zji,zjj) + zde ! Energy399 vt_i_init(ii,ij) = vt_i_init(ii,ij) + zv_newice(ji) ! volume 400 et_i_init(ii,ij) = et_i_init(ii,ij) + zde ! Energy 411 401 412 402 END DO … … 419 409 !----------------- 420 410 DO ji = 1, nbpac 421 zji = MOD( npac(ji) - 1 , jpi ) + 1422 zjj = ( npac(ji) - 1 ) / jpi + 1411 ii = MOD( npac(ji) - 1 , jpi ) + 1 412 ij = ( npac(ji) - 1 ) / jpi + 1 423 413 za_newice(ji) = zv_newice(ji) / zh_newice(ji) 424 diag_lat_gr( zji,zjj) = zv_newice(ji) * r1_rdtice414 diag_lat_gr(ii,ij) = diag_lat_gr(ii,ij) + zv_newice(ji) * r1_rdtice ! clem 425 415 END DO !ji 426 416 … … 441 431 ! we keep the excessive volume in memory and attribute it later to bottom accretion 442 432 DO ji = 1, nbpac 443 IF ( za_newice(ji) > ( 1._wp- zat_i_ac(ji) ) ) THEN444 zda_res(ji) = za_newice(ji) - ( 1.0- zat_i_ac(ji) )433 IF ( za_newice(ji) > ( amax - zat_i_ac(ji) ) ) THEN 434 zda_res(ji) = za_newice(ji) - ( amax - zat_i_ac(ji) ) 445 435 zdv_res(ji) = zda_res (ji) * zh_newice(ji) 446 436 za_newice(ji) = za_newice(ji) - zda_res (ji) … … 473 463 DO ji = 1, nbpac 474 464 jl = zcatac(ji) ! categroy in which new ice is put 475 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , -za_old(ji,jl) ) ) ! zindb=1 if ice =0 otherwise465 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , -za_old(ji,jl) + epsi10 ) ) ! zindb=1 if ice =0 otherwise 476 466 zhice_old(ji,jl) = zv_old(ji,jl) / MAX( za_old(ji,jl) , epsi10 ) * zindb ! old ice thickness 477 467 zdhex (ji) = MAX( 0._wp , zh_newice(ji) - zhice_old(ji,jl) ) ! difference in thickness 478 zswinew (ji) = MAX( 0._wp , SIGN( 1._wp , - za_old(ji,jl) + epsi1 1) ) ! ice totally new in jl category468 zswinew (ji) = MAX( 0._wp , SIGN( 1._wp , - za_old(ji,jl) + epsi10 ) ) ! ice totally new in jl category 479 469 END DO 480 470 … … 482 472 DO ji = 1, nbpac 483 473 jl = zcatac(ji) 484 zqold = ze_i_ac(ji,jk,jl) 485 zalphai = MIN( zhice_old(ji,jl) * jk / nlay_i, zh_newice(ji) ) &486 & - MIN( zhice_old(ji,jl) * ( jk - 1 ) / nlay_i, zh_newice(ji) )474 zqold = ze_i_ac(ji,jk,jl) ! [ J.m-3 ] 475 zalphai = MIN( zhice_old(ji,jl) * REAL( jk ) / REAL( nlay_i ), zh_newice(ji) ) & 476 & - MIN( zhice_old(ji,jl) * REAL( jk - 1 ) / REAL( nlay_i ), zh_newice(ji) ) 487 477 ze_i_ac(ji,jk,jl) = zswinew(ji) * ze_newice(ji) & 488 + ( 1.0 - zswinew(ji) ) * ( za_old(ji,jl) * zqold * zhice_old(ji,jl) / nlay_i&478 + ( 1.0 - zswinew(ji) ) * ( za_old(ji,jl) * zqold * zhice_old(ji,jl) / REAL( nlay_i ) & 489 479 + za_newice(ji) * ze_newice(ji) * zalphai & 490 + za_newice(ji) * ze_newice(ji) * zdhex(ji) / nlay_i ) / ( zv_i_ac(ji,jl) / nlay_i)480 + za_newice(ji) * ze_newice(ji) * zdhex(ji) / REAL( nlay_i ) ) / ( ( zv_i_ac(ji,jl) ) / REAL( nlay_i ) ) 491 481 END DO 492 482 END DO … … 513 503 DO ji = 1, nbpac 514 504 zindb = MAX( 0._wp, SIGN( 1._wp , zdv_res(ji) ) ) 515 zv_i_ac(ji,jl) = zv_i_ac(ji,jl) + zindb * zdv_res(ji) * za_i_ac(ji,jl) / MAX( zat_i_lev(ji) , epsi06 ) 505 zinda = MAX( 0._wp, SIGN( 1._wp , zat_i_lev(ji) - epsi06 ) ) ! clem 506 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) , epsi06 ) 516 507 END DO 517 508 END DO … … 524 515 DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) 525 516 DO ji = 1, nbpac 526 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl ) ) ) ! zindb=1 if ice =0 otherwise517 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl ) + epsi10 ) ) ! zindb=1 if ice =0 otherwise 527 518 zhice_old(ji,jl) = zv_i_ac(ji,jl) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb 528 519 zdhicbot (ji,jl) = zdv_res(ji) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb & … … 536 527 DO jk = 1, nlay_i 537 528 DO ji = 1, nbpac 538 zthick0(ji,jk,jl) = zhice_old(ji,jl) / nlay_i529 zthick0(ji,jk,jl) = zhice_old(ji,jl) / REAL( nlay_i ) 539 530 zqm0 (ji,jk,jl) = ze_i_ac(ji,jk,jl) * zthick0(ji,jk,jl) 540 531 END DO … … 555 546 DO layer = 1, nlay_i + 1 556 547 DO ji = 1, nbpac 557 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl) ) )548 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl) + epsi10 ) ) 558 549 ! Redistributing energy on the new grid 559 zweight = MAX ( MIN( zhice_old(ji,jl) * layer , zdummy(ji,jl) * jk) &560 & - MAX( zhice_old(ji,jl) * ( layer - 1 ) , zdummy(ji,jl) *( jk - 1 ) ) , 0._wp ) &561 & /( MAX( nlay_i* zthick0(ji,layer,jl),epsi10) ) * zindb550 zweight = MAX ( MIN( zhice_old(ji,jl) * REAL( layer ), zdummy(ji,jl) * REAL( jk ) ) & 551 & - MAX( zhice_old(ji,jl) * REAL( layer - 1 ) , zdummy(ji,jl) * REAL( jk - 1 ) ) , 0._wp ) & 552 & /( MAX(REAL(nlay_i) * zthick0(ji,layer,jl),epsi10) ) * zindb 562 553 ze_i_ac(ji,jk,jl) = ze_i_ac(ji,jk,jl) + zweight * zqm0(ji,layer,jl) 563 554 END DO ! ji … … 569 560 DO jk = 1, nlay_i 570 561 DO ji = 1, nbpac 571 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) ) )562 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) ) 572 563 ze_i_ac(ji,jk,jl) = ze_i_ac(ji,jk,jl) & 573 & / MAX( zv_i_ac(ji,jl) , epsi10) * za_i_ac(ji,jl) * nlay_i* zindb564 & / MAX( zv_i_ac(ji,jl) , epsi10) * za_i_ac(ji,jl) * REAL( nlay_i ) * zindb 574 565 END DO 575 566 END DO … … 581 572 DO jl = 1, jpl 582 573 DO ji = 1, nbpac 583 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl) ) ) ! 0 if no ice and 1 if yes574 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl) + epsi10 ) ) ! 0 if no ice and 1 if yes 584 575 zoa_i_ac(ji,jl) = za_old(ji,jl) * zoa_i_ac(ji,jl) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb 585 576 END DO … … 589 580 ! Update salinity 590 581 !----------------- 591 IF( num_sal == 2 ) THEN ! Sice = F(z,t)582 !clem IF( num_sal == 2 ) THEN 592 583 DO jl = 1, jpl 593 584 DO ji = 1, nbpac 594 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) ) )! 0 if no ice and 1 if yes585 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) ) ! 0 if no ice and 1 if yes 595 586 zdv = zv_i_ac(ji,jl) - zv_old(ji,jl) 596 zsmv_i_ac(ji,jl) = ( zsmv_i_ac(ji,jl) + zdv * zs_newice(ji) ) * zindb587 zsmv_i_ac(ji,jl) = zsmv_i_ac(ji,jl) + zdv * zs_newice(ji) * zindb ! clem modif 597 588 END DO 598 589 END DO 599 ENDIF 590 !clem ENDIF 591 592 !-------------------------------- 593 ! Update mass/salt fluxes (clem) 594 !-------------------------------- 595 DO jl = 1, jpl 596 DO ji = 1, nbpac 597 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) ) ! 0 if no ice and 1 if yes 598 zdv = zv_i_ac(ji,jl) - zv_old(ji,jl) 599 rdm_ice_1d(ji) = rdm_ice_1d(ji) + zdv * rhoic * zindb 600 sfx_thd_1d(ji) = sfx_thd_1d(ji) - zdv * rhoic * zs_newice(ji) * r1_rdtice * zindb 601 END DO 602 END DO 600 603 601 604 !------------------------------------------------------------------------------! … … 606 609 CALL tab_1d_2d( nbpac, v_i (:,:,jl), npac(1:nbpac), zv_i_ac (1:nbpac,jl), jpi, jpj ) 607 610 CALL tab_1d_2d( nbpac, oa_i(:,:,jl), npac(1:nbpac), zoa_i_ac(1:nbpac,jl), jpi, jpj ) 608 IF ( num_sal == 2 ) &611 !clem IF ( num_sal == 2 ) & 609 612 CALL tab_1d_2d( nbpac, smv_i (:,:,jl), npac(1:nbpac), zsmv_i_ac(1:nbpac,jl) , jpi, jpj ) 610 613 DO jk = 1, nlay_i … … 622 625 DO jl = 1, jpl 623 626 DO jk = 1, nlay_i ! heat content in 10^9 Joules 624 e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * area(:,:) * v_i(:,:,jl) / nlay_i/ unit_fac627 e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * area(:,:) * v_i(:,:,jl) / REAL( nlay_i ) / unit_fac 625 628 END DO 626 629 END DO
Note: See TracChangeset
for help on using the changeset viewer.