Changeset 1073 for trunk/NEMO/TOP_SRC/PISCES/p4zsink.F90
- Timestamp:
- 2008-06-05T14:15:34+02:00 (16 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/TOP_SRC/PISCES/p4zsink.F90
r935 r1073 13 13 USE oce_trc ! 14 14 USE trp_trc 15 USE sms 16 USE p4zsink2 ! 15 USE sms_pisces 17 16 USE prtctl_trc 18 17 … … 615 614 #endif 616 615 616 SUBROUTINE p4z_sink2( pwsink, psinkflx, jp_tra ) 617 !!--------------------------------------------------------------------- 618 !! *** ROUTINE p4z_sink2 *** 619 !! 620 !! ** Purpose : Compute the sedimentation terms for the various sinking 621 !! particles. The scheme used to compute the trends is based 622 !! on MUSCL. 623 !! 624 !! ** Method : - this ROUTINE compute not exactly the advection but the 625 !! transport term, i.e. div(u*tra). 626 !!--------------------------------------------------------------------- 627 INTEGER , INTENT(in ) :: jp_tra ! tracer index index 628 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj,jpk) :: pwsink ! sinking speed 629 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: psinkflx ! sinking fluxe 630 !! 631 INTEGER :: ji, jj, jk, jn 632 REAL(wp) :: zigma,zew,zstep,zign, zflx 633 REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztraz, zakz 634 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwsink2 635 !!--------------------------------------------------------------------- 636 637 zstep = rfact2 / 2. 638 639 ztraz(:,:,:) = 0.e0 640 zakz (:,:,:) = 0.e0 641 642 DO jk = 1, jpkm1 643 # if defined key_off_degrad 644 zwsink2(:,:,jk+1) = -pwsink(:,:,jk) / rjjss * tmask(:,:,jk+1) * facvol(:,:,jk) 645 # else 646 zwsink2(:,:,jk+1) = -pwsink(:,:,jk) / rjjss * tmask(:,:,jk+1) 647 # endif 648 END DO 649 zwsink2(:,:,1) = 0.e0 650 651 652 ! Vertical advective flux 653 DO jn = 1, 2 654 ! first guess of the slopes interior values 655 DO jk = 2, jpkm1 656 ztraz(:,:,jk) = ( trn(:,:,jk-1,jp_tra) - trn(:,:,jk,jp_tra) ) * tmask(:,:,jk) 657 END DO 658 ztraz(:,:,1 ) = 0.0 659 ztraz(:,:,jpk) = 0.0 660 661 ! slopes 662 DO jk = 2, jpkm1 663 DO jj = 1,jpj 664 DO ji = 1, jpi 665 zign = 0.25 + SIGN( 0.25, ztraz(ji,jj,jk) * ztraz(ji,jj,jk+1) ) 666 zakz(ji,jj,jk) = ( ztraz(ji,jj,jk) + ztraz(ji,jj,jk+1) ) * zign 667 END DO 668 END DO 669 END DO 670 671 ! Slopes limitation 672 DO jk = 2, jpkm1 673 DO jj = 1, jpj 674 DO ji = 1, jpi 675 zakz(ji,jj,jk) = SIGN( 1., zakz(ji,jj,jk) ) * & 676 & MIN( ABS( zakz(ji,jj,jk) ), 2. * ABS(ztraz(ji,jj,jk+1)), 2. * ABS(ztraz(ji,jj,jk) ) ) 677 END DO 678 END DO 679 END DO 680 681 ! vertical advective flux 682 DO jk = 1, jpkm1 683 DO jj = 1, jpj 684 DO ji = 1, jpi 685 zigma = zwsink2(ji,jj,jk+1) * zstep / fse3w(ji,jj,jk+1) 686 zew = zwsink2(ji,jj,jk+1) 687 psinkflx(ji,jj,jk+1) = -zew * ( trn(ji,jj,jk,jp_tra) - 0.5 * ( 1 + zigma ) * zakz(ji,jj,jk) ) * zstep 688 END DO 689 END DO 690 END DO 691 ! 692 ! Boundary conditions 693 psinkflx(:,:,1 ) = 0.e0 694 psinkflx(:,:,jpk) = 0.e0 695 696 DO jk=1,jpkm1 697 DO jj = 1,jpj 698 DO ji = 1, jpi 699 zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / fse3t(ji,jj,jk) 700 trn(ji,jj,jk,jp_tra) = trn(ji,jj,jk,jp_tra) + zflx 701 END DO 702 END DO 703 END DO 704 705 ENDDO 706 707 DO jk=1,jpkm1 708 DO jj = 1,jpj 709 DO ji = 1, jpi 710 zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / fse3t(ji,jj,jk) 711 trb(ji,jj,jk,jp_tra) = trb(ji,jj,jk,jp_tra) + 2. * zflx 712 END DO 713 END DO 714 END DO 715 716 trn(:,:,:,jp_tra) = trb(:,:,:,jp_tra) 717 psinkflx(:,:,:) = 2. * psinkflx(:,:,:) 718 719 ! 720 END SUBROUTINE p4z_sink2 721 617 722 #else 618 723 !!======================================================================
Note: See TracChangeset
for help on using the changeset viewer.