- Timestamp:
- 2016-11-06T17:31:33+01:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_fct.F90
r6900 r7200 38 38 REAL(wp) :: r1_6 = 1._wp / 6._wp ! =1/6 39 39 40 ! ! tridiag solver associated indices: 41 INTEGER, PARAMETER :: np_NH = 0 ! Neumann homogeneous boundary condition 42 INTEGER, PARAMETER :: np_CEN2 = 1 ! 2nd order centered boundary condition 43 40 44 !! * Substitutions 41 45 # include "vectopt_loop_substitute.h90" … … 706 710 707 711 708 SUBROUTINE interp_4th_cpt ( pt_in, pt_out )709 !!---------------------------------------------------------------------- 710 !! *** ROUTINE interp_4th_cpt ***712 SUBROUTINE interp_4th_cpt_org( pt_in, pt_out ) 713 !!---------------------------------------------------------------------- 714 !! *** ROUTINE interp_4th_cpt_org *** 711 715 !! 712 716 !! ** Purpose : Compute the interpolation of tracer at w-point … … 739 743 END DO 740 744 ! 741 jk =2! Switch to second order centered at top742 DO jj =1,jpj743 DO ji =1,jpi745 jk = 2 ! Switch to second order centered at top 746 DO jj = 1, jpj 747 DO ji = 1, jpi 744 748 zwd (ji,jj,jk) = 1._wp 745 749 zwi (ji,jj,jk) = 0._wp … … 789 793 END DO 790 794 ! 795 END SUBROUTINE interp_4th_cpt_org 796 797 798 SUBROUTINE interp_4th_cpt( pt_in, pt_out ) 799 !!---------------------------------------------------------------------- 800 !! *** ROUTINE interp_4th_cpt *** 801 !! 802 !! ** Purpose : Compute the interpolation of tracer at w-point 803 !! 804 !! ** Method : 4th order compact interpolation 805 !!---------------------------------------------------------------------- 806 REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pt_in ! field at t-point 807 REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT( out) :: pt_out ! field interpolated at w-point 808 ! 809 INTEGER :: ji, jj, jk ! dummy loop integers 810 INTEGER :: ikt, ikb ! local integers 811 REAL(wp),DIMENSION(jpi,jpj,jpk) :: zwd, zwi, zws, zwrm, zwt 812 !!---------------------------------------------------------------------- 813 ! 814 ! !== build the three diagonal matrix & the RHS ==! 815 ! 816 DO jk = 3, jpkm1 ! interior (from jk=3 to jpk-1) 817 DO jj = 2, jpjm1 818 DO ji = fs_2, fs_jpim1 819 zwd (ji,jj,jk) = 3._wp * wmask(ji,jj,jk) + 1._wp ! diagonal 820 zwi (ji,jj,jk) = wmask(ji,jj,jk) ! lower diagonal 821 zws (ji,jj,jk) = wmask(ji,jj,jk) ! upper diagonal 822 zwrm(ji,jj,jk) = 3._wp * wmask(ji,jj,jk) & ! RHS 823 & * ( pt_in(ji,jj,jk) + pt_in(ji,jj,jk-1) ) 824 END DO 825 END DO 826 END DO 827 ! 828 !!gm 829 ! SELECT CASE( kbc ) !* boundary condition 830 ! CASE( np_NH ) ! Neumann homogeneous at top & bottom 831 ! CASE( np_CEN2 ) ! 2nd order centered at top & bottom 832 ! END SELECT 833 !!gm 834 ! 835 DO jj = 2, jpjm1 ! 2nd order centered at top & bottom 836 DO ji = fs_2, fs_jpim1 837 ikt = mikt(ji,jj) + 1 ! w-point below the 1st wet point 838 ikb = mbkt(ji,jj) ! - above the last wet point 839 ! 840 zwd (ji,jj,ikt) = 1._wp ! top 841 zwi (ji,jj,ikt) = 0._wp 842 zws (ji,jj,ikt) = 0._wp 843 zwrm(ji,jj,ikt) = 0.5_wp * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) ) 844 ! 845 zwd (ji,jj,ikb) = 1._wp ! bottom 846 zwi (ji,jj,ikb) = 0._wp 847 zws (ji,jj,ikb) = 0._wp 848 zwrm(ji,jj,ikb) = 0.5_wp * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) ) 849 END DO 850 END DO 851 ! 852 ! !== tridiagonal solver ==! 853 ! 854 DO jj = 2, jpjm1 !* 1st recurrence: Tk = Dk - Ik Sk-1 / Tk-1 855 DO ji = fs_2, fs_jpim1 856 zwt(ji,jj,2) = zwd(ji,jj,2) 857 END DO 858 END DO 859 DO jk = 3, jpkm1 860 DO jj = 2, jpjm1 861 DO ji = fs_2, fs_jpim1 862 zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1) 863 END DO 864 END DO 865 END DO 866 ! 867 DO jj = 2, jpjm1 !* 2nd recurrence: Zk = Yk - Ik / Tk-1 Zk-1 868 DO ji = fs_2, fs_jpim1 869 pt_out(ji,jj,2) = zwrm(ji,jj,2) 870 END DO 871 END DO 872 DO jk = 3, jpkm1 873 DO jj = 2, jpjm1 874 DO ji = fs_2, fs_jpim1 875 pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1) 876 END DO 877 END DO 878 END DO 879 880 DO jj = 2, jpjm1 !* 3d recurrence: Xk = (Zk - Sk Xk+1 ) / Tk 881 DO ji = fs_2, fs_jpim1 882 pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) 883 END DO 884 END DO 885 DO jk = jpk-2, 2, -1 886 DO jj = 2, jpjm1 887 DO ji = fs_2, fs_jpim1 888 pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk) 889 END DO 890 END DO 891 END DO 892 ! 791 893 END SUBROUTINE interp_4th_cpt 792 894 895 896 SUBROUTINE tridia_solver( pD, pU, pL, pRHS, pt_out , klev ) 897 !!---------------------------------------------------------------------- 898 !! *** ROUTINE tridia_solver *** 899 !! 900 !! ** Purpose : solve a symmetric 3diagonal system 901 !! 902 !! ** Method : solve M.t_out = RHS(t) where M is a tri diagonal matrix ( jpk*jpk ) 903 !! 904 !! ( D_1 U_1 0 0 0 )( t_1 ) ( RHS_1 ) 905 !! ( L_2 D_2 U_2 0 0 )( t_2 ) ( RHS_2 ) 906 !! ( 0 L_3 D_3 U_3 0 )( t_3 ) = ( RHS_3 ) 907 !! ( ... )( ... ) ( ... ) 908 !! ( 0 0 0 L_k D_k )( t_k ) ( RHS_k ) 909 !! 910 !! M is decomposed in the product of an upper and lower triangular matrix. 911 !! The tri-diagonals matrix is given as input 3D arrays: pD, pU, pL 912 !! (i.e. the Diagonal, the Upper diagonal, and the Lower diagonal). 913 !! The solution is pta. 914 !! The 3d array zwt is used as a work space array. 915 !!---------------------------------------------------------------------- 916 REAL(wp),DIMENSION(:,:,:), INTENT(in ) :: pD, pU, PL ! 3-diagonal matrix 917 REAL(wp),DIMENSION(:,:,:), INTENT(in ) :: pRHS ! Right-Hand-Side 918 REAL(wp),DIMENSION(:,:,:), INTENT( out) :: pt_out !!gm field at level=F(klev) 919 INTEGER , INTENT(in ) :: klev ! =1 pt_out at w-level 920 ! ! =0 pt at t-level 921 INTEGER :: ji, jj, jk ! dummy loop integers 922 INTEGER :: kstart ! local indices 923 REAL(wp),DIMENSION(jpi,jpj,jpk) :: zwt ! 3D work array 924 !!---------------------------------------------------------------------- 925 ! 926 kstart = 1 + klev 927 ! 928 DO jj = 2, jpjm1 !* 1st recurrence: Tk = Dk - Ik Sk-1 / Tk-1 929 DO ji = fs_2, fs_jpim1 930 zwt(ji,jj,kstart) = pD(ji,jj,kstart) 931 END DO 932 END DO 933 DO jk = kstart+1, jpkm1 934 DO jj = 2, jpjm1 935 DO ji = fs_2, fs_jpim1 936 zwt(ji,jj,jk) = pD(ji,jj,jk) - pL(ji,jj,jk) * pU(ji,jj,jk-1) /zwt(ji,jj,jk-1) 937 END DO 938 END DO 939 END DO 940 ! 941 DO jj = 2, jpjm1 !* 2nd recurrence: Zk = Yk - Ik / Tk-1 Zk-1 942 DO ji = fs_2, fs_jpim1 943 pt_out(ji,jj,kstart) = pRHS(ji,jj,kstart) 944 END DO 945 END DO 946 DO jk = kstart+1, jpkm1 947 DO jj = 2, jpjm1 948 DO ji = fs_2, fs_jpim1 949 pt_out(ji,jj,jk) = pRHS(ji,jj,jk) - pL(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1) 950 END DO 951 END DO 952 END DO 953 954 DO jj = 2, jpjm1 !* 3d recurrence: Xk = (Zk - Sk Xk+1 ) / Tk 955 DO ji = fs_2, fs_jpim1 956 pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) 957 END DO 958 END DO 959 DO jk = jpk-2, kstart, -1 960 DO jj = 2, jpjm1 961 DO ji = fs_2, fs_jpim1 962 pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - pU(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk) 963 END DO 964 END DO 965 END DO 966 ! 967 END SUBROUTINE tridia_solver 968 793 969 !!====================================================================== 794 970 END MODULE traadv_fct
Note: See TracChangeset
for help on using the changeset viewer.
