- Timestamp:
- 2013-06-26T09:54:16+02:00 (11 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
r3294 r3938 29 29 USE lib_mpp ! MPP library 30 30 USE wrk_nemo ! work arrays 31 USE lib_fortran ! to use key_nosignedzero 31 32 32 33 IMPLICIT NONE … … 159 160 !Energy of melting q(S,T) [J.m-3] 160 161 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / & 161 MAX( area(ji,jj) * v_i(ji,jj,jl) , epsi10 ) * & 162 nlay_i 162 MAX( area(ji,jj) * v_i(ji,jj,jl) , epsi10 ) * REAL( nlay_i ) 163 163 zindb = 1.0-MAX(0.0,SIGN(1.0,-v_i(ji,jj,jl))) !0 if no ice and 1 if yes 164 164 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl)*unit_fac*zindb … … 185 185 186 186 ! Default new ice thickness 187 DO jj = 1, jpj 188 DO ji = 1, jpi 189 hicol(ji,jj) = hiccrit(1) 190 END DO 191 END DO 187 hicol(:,:) = hiccrit(1) 192 188 193 189 IF (fraz_swi.eq.1.0) THEN … … 335 331 CALL tab_2d_1d( nbpac, zvrel_ac (1:nbpac) , zvrel , & 336 332 jpi, jpj, npac(1:nbpac) ) 333 CALL tab_2d_1d( nbpac, rdmicif_1d (1:nbpac) , rdmicif, jpi, jpj, npac(1:nbpac) ) !martin 337 334 338 335 !------------------------------------------------------------------------------! … … 410 407 zdh_frazb(ji) = zfrazb*zv_newice(ji) 411 408 zv_newice(ji) = ( 1.0 - zfrazb ) * zv_newice(ji) 409 ! 410 !clem@mv rdmicif_1d(ji) = rdmicif_1d(ji) + zv_newice(ji) * rhoic !martin 412 411 END DO 413 412 … … 415 414 ! Salt flux due to new ice growth 416 415 !--------------------------------- 417 IF ( ( num_sal .EQ. 4 ) ) THEN 418 DO ji = 1, nbpac 419 zji = MOD( npac(ji) - 1, jpi ) + 1 420 zjj = ( npac(ji) - 1 ) / jpi + 1 421 fseqv_1d(ji) = fseqv_1d(ji) + & 422 ( sss_m(zji,zjj) - bulk_sal ) * rhoic * & 423 zv_newice(ji) / rdt_ice 424 END DO 425 ELSE 426 DO ji = 1, nbpac 427 zji = MOD( npac(ji) - 1, jpi ) + 1 428 zjj = ( npac(ji) - 1 ) / jpi + 1 429 fseqv_1d(ji) = fseqv_1d(ji) + & 430 ( sss_m(zji,zjj) - zs_newice(ji) ) * rhoic * & 431 zv_newice(ji) / rdt_ice 432 END DO ! ji 433 ENDIF 416 !clem@mv IF ( ( num_sal .EQ. 4 ) ) THEN 417 !clem@mv DO ji = 1, nbpac 418 !clem@mv ! IOVINO 419 !clem@mv fseqv_1d(ji) = fseqv_1d(ji) - bulk_sal * rhoic * & 420 !clem@mv zv_newice(ji) / rdt_ice 421 !clem@mv END DO 422 !clem@mv ELSE 423 !clem@mv DO ji = 1, nbpac 424 !clem@mv ! IOVINO 425 !clem@mv fseqv_1d(ji) = fseqv_1d(ji) - zs_newice(ji) * rhoic * & 426 !clem@mv zv_newice(ji) / rdt_ice 427 !clem@mv END DO ! ji 428 !clem@mv ENDIF 434 429 435 430 !------------------------------------ … … 460 455 zjj = ( npac(ji) - 1 ) / jpi + 1 461 456 diag_lat_gr(zji,zjj) = zv_newice(ji) / rdt_ice 462 END DO !ji457 END DO !ji 463 458 464 459 !------------------------------------------------------------------------------! … … 480 475 DO ji = 1, nbpac 481 476 ! vectorize 482 IF ( za_newice(ji) .GT. ( 1.0- zat_i_ac(ji) ) ) THEN483 zda_res(ji) = za_newice(ji) - ( 1.0- zat_i_ac(ji) )477 IF ( za_newice(ji) .GT. ( amax - zat_i_ac(ji) ) ) THEN 478 zda_res(ji) = za_newice(ji) - ( amax - zat_i_ac(ji) ) 484 479 zdv_res(ji) = zda_res(ji) * zh_newice(ji) 485 480 za_newice(ji) = za_newice(ji) - zda_res(ji) 486 481 zv_newice(ji) = zv_newice(ji) - zdv_res(ji) 487 482 ELSE 488 483 zda_res(ji) = 0.0 489 484 zdv_res(ji) = 0.0 … … 512 507 DO ji = 1, nbpac 513 508 jl = zcatac(ji) ! categroy in which new ice is put 514 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , -za_old(ji,jl) ) ) ! zindb=1 if ice =0 otherwise509 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , -za_old(ji,jl) + epsi10 ) ) ! zindb=1 if ice =0 otherwise 515 510 zhice_old(ji,jl) = zv_old(ji,jl) / MAX( za_old(ji,jl) , epsi10 ) * zindb ! old ice thickness 516 511 zdhex (ji) = MAX( 0._wp , zh_newice(ji) - zhice_old(ji,jl) ) ! difference in thickness 517 zswinew (ji) = MAX( 0._wp , SIGN( 1._wp , - za_old(ji,jl) + epsi1 1) ) ! ice totally new in jl category512 zswinew (ji) = MAX( 0._wp , SIGN( 1._wp , - za_old(ji,jl) + epsi10 ) ) ! ice totally new in jl category 518 513 END DO 519 514 … … 522 517 jl = zcatac(ji) 523 518 zqold = ze_i_ac(ji,jk,jl) ! [ J.m-3 ] 524 zalphai = MIN( zhice_old(ji,jl) * jk / nlay_i, zh_newice(ji) ) &525 & - MIN( zhice_old(ji,jl) * ( jk - 1 ) / nlay_i, zh_newice(ji) )519 zalphai = MIN( zhice_old(ji,jl) * REAL( jk ) / REAL( nlay_i ), zh_newice(ji) ) & 520 & - MIN( zhice_old(ji,jl) * REAL( jk - 1 ) / REAL( nlay_i ), zh_newice(ji) ) 526 521 ze_i_ac(ji,jk,jl) = zswinew(ji) * ze_newice(ji) & 527 + ( 1.0 - zswinew(ji) ) * ( za_old(ji,jl) * zqold * zhice_old(ji,jl) / nlay_i&522 + ( 1.0 - zswinew(ji) ) * ( za_old(ji,jl) * zqold * zhice_old(ji,jl) / REAL( nlay_i ) & 528 523 + za_newice(ji) * ze_newice(ji) * zalphai & 529 + za_newice(ji) * ze_newice(ji) * zdhex(ji) / nlay_i ) / ( ( zv_i_ac(ji,jl) ) / nlay_i)524 + za_newice(ji) * ze_newice(ji) * zdhex(ji) / REAL( nlay_i ) ) / ( ( zv_i_ac(ji,jl) ) / REAL( nlay_i ) ) 530 525 END DO 531 526 END DO … … 563 558 DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) 564 559 DO ji = 1, nbpac 565 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl ) ) ) ! zindb=1 if ice =0 otherwise560 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl ) + epsi10 ) ) ! zindb=1 if ice =0 otherwise 566 561 zhice_old(ji,jl) = zv_i_ac(ji,jl) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb 567 562 zdhicbot (ji,jl) = zdv_res(ji) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb & … … 575 570 DO jk = 1, nlay_i 576 571 DO ji = 1, nbpac 577 zthick0(ji,jk,jl) = zhice_old(ji,jl) / nlay_i572 zthick0(ji,jk,jl) = zhice_old(ji,jl) / REAL( nlay_i ) 578 573 zqm0 (ji,jk,jl) = ze_i_ac(ji,jk,jl) * zthick0(ji,jk,jl) 579 574 END DO … … 594 589 DO layer = 1, nlay_i + 1 595 590 DO ji = 1, nbpac 596 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl) ) )591 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl) + epsi10 ) ) 597 592 ! Redistributing energy on the new grid 598 zweight = MAX ( MIN( zhice_old(ji,jl) * layer , zdummy(ji,jl) * jk) &599 & - MAX( zhice_old(ji,jl) * ( layer - 1 ) , zdummy(ji,jl) *( jk - 1 ) ) , 0._wp ) &600 & /( MAX( nlay_i* zthick0(ji,layer,jl),epsi10) ) * zindb593 zweight = MAX ( MIN( zhice_old(ji,jl) * REAL( layer ), zdummy(ji,jl) * REAL( jk ) ) & 594 & - MAX( zhice_old(ji,jl) * REAL( layer - 1 ) , zdummy(ji,jl) * REAL( jk - 1 ) ) , 0._wp ) & 595 & /( MAX(REAL(nlay_i) * zthick0(ji,layer,jl),epsi10) ) * zindb 601 596 ze_i_ac(ji,jk,jl) = ze_i_ac(ji,jk,jl) + zweight * zqm0(ji,layer,jl) 602 597 END DO ! ji … … 608 603 DO jk = 1, nlay_i 609 604 DO ji = 1, nbpac 610 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) ) )605 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) ) 611 606 ze_i_ac(ji,jk,jl) = ze_i_ac(ji,jk,jl) & 612 & / MAX( zv_i_ac(ji,jl) , epsi10) * za_i_ac(ji,jl) * nlay_i* zindb607 & / MAX( zv_i_ac(ji,jl) , epsi10) * za_i_ac(ji,jl) * REAL( nlay_i ) * zindb 613 608 END DO 614 609 END DO … … 620 615 DO jl = 1, jpl 621 616 DO ji = 1, nbpac 622 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl) ) ) ! 0 if no ice and 1 if yes617 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl) + epsi10 ) ) ! 0 if no ice and 1 if yes 623 618 zoa_i_ac(ji,jl) = za_old(ji,jl) * zoa_i_ac(ji,jl) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb 624 619 END DO … … 631 626 DO jl = 1, jpl 632 627 DO ji = 1, nbpac 633 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) ) ) ! 0 if no ice and 1 if yes628 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) ) ! 0 if no ice and 1 if yes 634 629 zdv = zv_i_ac(ji,jl) - zv_old(ji,jl) 635 630 zsmv_i_ac(ji,jl) = ( zsmv_i_ac(ji,jl) + zdv * zs_newice(ji) ) * zindb … … 637 632 END DO 638 633 ENDIF 634 635 !-------------------------------- 636 ! Update mass/salt fluxes (clem) 637 !-------------------------------- 638 DO jl = 1, jpl 639 DO ji = 1, nbpac 640 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) ) ! 0 if no ice and 1 if yes 641 zdv = zv_i_ac(ji,jl) - zv_old(ji,jl) 642 rdmicif_1d(ji) = rdmicif_1d(ji) + zdv * rhoic !* zindb 643 fseqv_1d(ji) = fseqv_1d(ji) - zdv * rhoic * zs_newice(ji) / rdt_ice * zindb 644 END DO 645 END DO 639 646 640 647 !------------------------------------------------------------------------------! … … 652 659 END DO 653 660 CALL tab_1d_2d( nbpac, fseqv , npac(1:nbpac), fseqv_1d (1:nbpac) , jpi, jpj ) 661 CALL tab_1d_2d( nbpac, rdmicif , npac(1:nbpac), rdmicif_1d (1:nbpac) , jpi, jpj ) 654 662 ! 655 663 ENDIF ! nbpac > 0 … … 660 668 DO jl = 1, jpl 661 669 DO jk = 1, nlay_i ! heat content in 10^9 Joules 662 e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * area(:,:) * v_i(:,:,jl) / nlay_i/ unit_fac670 e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * area(:,:) * v_i(:,:,jl) / REAL( nlay_i ) / unit_fac 663 671 END DO 664 672 END DO
Note: See TracChangeset
for help on using the changeset viewer.
