Changeset 7646 for trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_fct.F90
- Timestamp:
- 2017-02-06T10:25:03+01:00 (7 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_fct.F90 (modified) (13 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_fct.F90
r6771 r7646 20 20 USE trdtra ! tracers trends 21 21 USE diaptr ! poleward transport diagnostics 22 USE diaar5 ! AR5 diagnostics 23 USE phycst, ONLY: rau0_rcp 22 24 ! 23 25 USE in_out_manager ! I/O manager 26 USE iom 24 27 USE lib_mpp ! MPP library 25 28 USE lbclnk ! ocean lateral boundary condition (or mpp link) … … 36 39 37 40 LOGICAL :: l_trd ! flag to compute trends 41 LOGICAL :: l_ptr ! flag to compute poleward transport 42 LOGICAL :: l_hst ! flag to compute heat/salt transport 38 43 REAL(wp) :: r1_6 = 1._wp / 6._wp ! =1/6 44 45 ! ! tridiag solver associated indices: 46 INTEGER, PARAMETER :: np_NH = 0 ! Neumann homogeneous boundary condition 47 INTEGER, PARAMETER :: np_CEN2 = 1 ! 2nd order centered boundary condition 39 48 40 49 !! * Substitutions … … 80 89 REAL(wp) :: zfm_ui, zfm_vj, zfm_wk, zC2t_v, zC4t_v ! - - 81 90 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw 82 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdx, ztrdy, ztrdz 91 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdx, ztrdy, ztrdz, zptry 92 REAL(wp), POINTER, DIMENSION(:,:) :: z2d 83 93 !!---------------------------------------------------------------------- 84 94 ! … … 94 104 ! 95 105 l_trd = .FALSE. 96 IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 97 ! 98 IF( l_trd ) THEN 106 l_hst = .FALSE. 107 l_ptr = .FALSE. 108 IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 109 IF( cdtype == 'TRA' .AND. ln_diaptr ) l_ptr = .TRUE. 110 IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. & 111 & iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) ) l_hst = .TRUE. 112 ! 113 IF( l_trd .OR. l_hst ) THEN 99 114 CALL wrk_alloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz ) 100 115 ztrdx(:,:,:) = 0._wp ; ztrdy(:,:,:) = 0._wp ; ztrdz(:,:,:) = 0._wp 101 116 ENDIF 102 117 ! 118 IF( l_ptr ) THEN 119 CALL wrk_alloc( jpi, jpj, jpk, zptry ) 120 zptry(:,:,:) = 0._wp 121 ENDIF 103 122 ! ! surface & bottom value : flux set to zero one for all 104 123 zwz(:,:, 1 ) = 0._wp … … 161 180 CALL lbc_lnk( zwi, 'T', 1. ) ! Lateral boundary conditions on zwi (unchanged sign) 162 181 ! 163 IF( l_trd ) THEN ! trend diagnostics (contribution of upstream fluxes)182 IF( l_trd .OR. l_hst ) THEN ! trend diagnostics (contribution of upstream fluxes) 164 183 ztrdx(:,:,:) = zwx(:,:,:) ; ztrdy(:,:,:) = zwy(:,:,:) ; ztrdz(:,:,:) = zwz(:,:,:) 165 184 END IF 166 185 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) 167 IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 168 IF( jn == jp_tem ) htr_adv(:) = ptr_sj( zwy(:,:,:) ) 169 IF( jn == jp_sal ) str_adv(:) = ptr_sj( zwy(:,:,:) ) 170 ENDIF 186 IF( l_ptr ) zptry(:,:,:) = zwy(:,:,:) 171 187 ! 172 188 ! !== anti-diffusive flux : high order minus low order ==! … … 292 308 END DO 293 309 ! 294 IF( l_trd ) THEN ! trend diagnostics (contribution of upstream fluxes)310 IF( l_trd .OR. l_hst ) THEN ! trend diagnostics (contribution of upstream fluxes) 295 311 ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:) ! <<< Add to previously computed 296 312 ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:) ! <<< Add to previously computed 297 313 ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:) ! <<< Add to previously computed 298 ! 314 ENDIF 315 ! 316 IF( l_trd ) THEN 299 317 CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pun, ptn(:,:,:,jn) ) 300 318 CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pvn, ptn(:,:,:,jn) ) 301 319 CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pwn, ptn(:,:,:,jn) ) 302 320 ! 303 CALL wrk_dealloc( jpi,jpj,jpk, ztrdx, ztrdy, ztrdz )304 321 END IF 305 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) 306 IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 307 IF( jn == jp_tem ) htr_adv(:) = htr_adv(:) + ptr_sj( zwy(:,:,:) ) 308 IF( jn == jp_sal ) str_adv(:) = str_adv(:) + ptr_sj( zwy(:,:,:) ) 322 ! ! heat/salt transport 323 IF( l_hst ) CALL dia_ar5_hst( jn, 'adv', ztrdx(:,:,:), ztrdy(:,:,:) ) 324 325 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) 326 IF( l_ptr ) THEN 327 zptry(:,:,:) = zptry(:,:,:) + zwy(:,:,:) ! <<< Add to previously computed 328 CALL dia_ptr_hst( jn, 'adv', zptry(:,:,:) ) 309 329 ENDIF 310 330 ! 311 331 END DO ! end of tracer loop 312 332 ! 313 CALL wrk_dealloc( jpi,jpj,jpk, zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw ) 333 CALL wrk_dealloc( jpi,jpj,jpk, zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw ) 334 IF( l_trd .OR. l_hst ) CALL wrk_dealloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz ) 335 IF( l_ptr ) CALL wrk_dealloc( jpi, jpj, jpk, zptry ) 314 336 ! 315 337 IF( nn_timing == 1 ) CALL timing_stop('tra_adv_fct') … … 357 379 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwi, zwx, zwy, zwz, zhdiv, zwzts, zwz_sav 358 380 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdx, ztrdy, ztrdz 381 REAL(wp), POINTER, DIMENSION(:,:,:) :: zptry 359 382 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: ztrs 360 383 !!---------------------------------------------------------------------- … … 373 396 ! 374 397 l_trd = .FALSE. 375 IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 376 ! 377 IF( l_trd ) THEN 398 l_hst = .FALSE. 399 l_ptr = .FALSE. 400 IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 401 IF( cdtype == 'TRA' .AND. ln_diaptr ) l_ptr = .TRUE. 402 IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. & 403 & iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) ) l_hst = .TRUE. 404 ! 405 IF( l_trd .OR. l_hst ) THEN 378 406 CALL wrk_alloc( jpi,jpj,jpk, ztrdx, ztrdy, ztrdz ) 379 407 ztrdx(:,:,:) = 0._wp ; ztrdy(:,:,:) = 0._wp ; ztrdz(:,:,:) = 0._wp 380 408 ENDIF 381 409 ! 410 IF( l_ptr ) THEN 411 CALL wrk_alloc( jpi, jpj,jpk, zptry ) 412 zptry(:,:,:) = 0._wp 413 ENDIF 382 414 zwi(:,:,:) = 0._wp 383 415 z_rzts = 1._wp / REAL( kn_fct_zts, wp ) … … 445 477 CALL lbc_lnk( zwi, 'T', 1. ) ! Lateral boundary conditions on zwi (unchanged sign) 446 478 ! 447 IF( l_trd ) THEN ! trend diagnostics (contribution of upstream fluxes)479 IF( l_trd .OR. l_hst ) THEN ! trend diagnostics (contribution of upstream fluxes) 448 480 ztrdx(:,:,:) = zwx(:,:,:) ; ztrdy(:,:,:) = zwy(:,:,:) ; ztrdz(:,:,:) = zwz(:,:,:) 449 481 END IF 450 482 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) 451 IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 452 IF( jn == jp_tem ) htr_adv(:) = ptr_sj( zwy(:,:,:) ) 453 IF( jn == jp_sal ) str_adv(:) = ptr_sj( zwy(:,:,:) ) 454 ENDIF 483 IF( l_ptr ) zptry(:,:,:) = zwy(:,:,:) 455 484 456 485 ! 3. anti-diffusive flux : high order minus low order … … 568 597 END DO 569 598 570 ! ! trend diagnostics (contribution of upstream fluxes)571 IF( l_trd ) THEN599 ! 600 IF( l_trd .OR. l_hst ) THEN ! trend diagnostics (contribution of upstream fluxes) 572 601 ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:) ! <<< Add to previously computed 573 602 ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:) ! <<< Add to previously computed 574 603 ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:) ! <<< Add to previously computed 575 ! 576 CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pun, ptn(:,:,:,jn) ) 577 CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pvn, ptn(:,:,:,jn) ) 578 CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pwn, ptn(:,:,:,jn) ) 579 ! 580 CALL wrk_dealloc( jpi,jpj,jpk, ztrdx, ztrdy, ztrdz ) 604 ENDIF 605 ! 606 IF( l_trd ) THEN 607 CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pun, ptn(:,:,:,jn) ) 608 CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pvn, ptn(:,:,:,jn) ) 609 CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pwn, ptn(:,:,:,jn) ) 610 ! 581 611 END IF 582 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) 583 IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 584 IF( jn == jp_tem ) htr_adv(:) = ptr_sj( zwy(:,:,:) ) + htr_adv(:) 585 IF( jn == jp_sal ) str_adv(:) = ptr_sj( zwy(:,:,:) ) + str_adv(:) 612 ! ! heat/salt transport 613 IF( l_hst ) CALL dia_ar5_hst( jn, 'adv', ztrdx(:,:,:), ztrdy(:,:,:) ) 614 615 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) 616 IF( l_ptr ) THEN 617 zptry(:,:,:) = zptry(:,:,:) + zwy(:,:,:) ! <<< Add to previously computed 618 CALL dia_ptr_hst( jn, 'adv', zptry(:,:,:) ) 586 619 ENDIF 587 620 ! 588 621 END DO 589 622 ! 590 CALL wrk_alloc( jpi,jpj, zwx_sav, zwy_sav ) 591 CALL wrk_alloc( jpi,jpj, jpk, zwx, zwy, zwz, zwi, zhdiv, zwzts, zwz_sav ) 592 CALL wrk_alloc( jpi,jpj,jpk,kjpt+1, ztrs ) 623 CALL wrk_alloc( jpi,jpj, zwx_sav, zwy_sav ) 624 CALL wrk_alloc( jpi,jpj, jpk, zwx, zwy, zwz, zwi, zhdiv, zwzts, zwz_sav ) 625 CALL wrk_alloc( jpi,jpj,jpk,kjpt+1, ztrs ) 626 IF( l_trd .OR. l_hst ) CALL wrk_dealloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz ) 627 IF( l_ptr ) CALL wrk_dealloc( jpi, jpj, jpk, zptry ) 593 628 ! 594 629 IF( nn_timing == 1 ) CALL timing_stop('tra_adv_fct_zts') … … 706 741 707 742 708 SUBROUTINE interp_4th_cpt ( pt_in, pt_out )709 !!---------------------------------------------------------------------- 710 !! *** ROUTINE interp_4th_cpt ***743 SUBROUTINE interp_4th_cpt_org( pt_in, pt_out ) 744 !!---------------------------------------------------------------------- 745 !! *** ROUTINE interp_4th_cpt_org *** 711 746 !! 712 747 !! ** Purpose : Compute the interpolation of tracer at w-point … … 739 774 END DO 740 775 ! 741 jk =2! Switch to second order centered at top742 DO jj =1,jpj743 DO ji =1,jpi776 jk = 2 ! Switch to second order centered at top 777 DO jj = 1, jpj 778 DO ji = 1, jpi 744 779 zwd (ji,jj,jk) = 1._wp 745 780 zwi (ji,jj,jk) = 0._wp … … 789 824 END DO 790 825 ! 826 END SUBROUTINE interp_4th_cpt_org 827 828 829 SUBROUTINE interp_4th_cpt( pt_in, pt_out ) 830 !!---------------------------------------------------------------------- 831 !! *** ROUTINE interp_4th_cpt *** 832 !! 833 !! ** Purpose : Compute the interpolation of tracer at w-point 834 !! 835 !! ** Method : 4th order compact interpolation 836 !!---------------------------------------------------------------------- 837 REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pt_in ! field at t-point 838 REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT( out) :: pt_out ! field interpolated at w-point 839 ! 840 INTEGER :: ji, jj, jk ! dummy loop integers 841 INTEGER :: ikt, ikb ! local integers 842 REAL(wp),DIMENSION(jpi,jpj,jpk) :: zwd, zwi, zws, zwrm, zwt 843 !!---------------------------------------------------------------------- 844 ! 845 ! !== build the three diagonal matrix & the RHS ==! 846 ! 847 DO jk = 3, jpkm1 ! interior (from jk=3 to jpk-1) 848 DO jj = 2, jpjm1 849 DO ji = fs_2, fs_jpim1 850 zwd (ji,jj,jk) = 3._wp * wmask(ji,jj,jk) + 1._wp ! diagonal 851 zwi (ji,jj,jk) = wmask(ji,jj,jk) ! lower diagonal 852 zws (ji,jj,jk) = wmask(ji,jj,jk) ! upper diagonal 853 zwrm(ji,jj,jk) = 3._wp * wmask(ji,jj,jk) & ! RHS 854 & * ( pt_in(ji,jj,jk) + pt_in(ji,jj,jk-1) ) 855 END DO 856 END DO 857 END DO 858 ! 859 !!gm 860 ! SELECT CASE( kbc ) !* boundary condition 861 ! CASE( np_NH ) ! Neumann homogeneous at top & bottom 862 ! CASE( np_CEN2 ) ! 2nd order centered at top & bottom 863 ! END SELECT 864 !!gm 865 ! 866 DO jj = 2, jpjm1 ! 2nd order centered at top & bottom 867 DO ji = fs_2, fs_jpim1 868 ikt = mikt(ji,jj) + 1 ! w-point below the 1st wet point 869 ikb = mbkt(ji,jj) ! - above the last wet point 870 ! 871 zwd (ji,jj,ikt) = 1._wp ! top 872 zwi (ji,jj,ikt) = 0._wp 873 zws (ji,jj,ikt) = 0._wp 874 zwrm(ji,jj,ikt) = 0.5_wp * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) ) 875 ! 876 zwd (ji,jj,ikb) = 1._wp ! bottom 877 zwi (ji,jj,ikb) = 0._wp 878 zws (ji,jj,ikb) = 0._wp 879 zwrm(ji,jj,ikb) = 0.5_wp * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) ) 880 END DO 881 END DO 882 ! 883 ! !== tridiagonal solver ==! 884 ! 885 DO jj = 2, jpjm1 !* 1st recurrence: Tk = Dk - Ik Sk-1 / Tk-1 886 DO ji = fs_2, fs_jpim1 887 zwt(ji,jj,2) = zwd(ji,jj,2) 888 END DO 889 END DO 890 DO jk = 3, jpkm1 891 DO jj = 2, jpjm1 892 DO ji = fs_2, fs_jpim1 893 zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1) 894 END DO 895 END DO 896 END DO 897 ! 898 DO jj = 2, jpjm1 !* 2nd recurrence: Zk = Yk - Ik / Tk-1 Zk-1 899 DO ji = fs_2, fs_jpim1 900 pt_out(ji,jj,2) = zwrm(ji,jj,2) 901 END DO 902 END DO 903 DO jk = 3, jpkm1 904 DO jj = 2, jpjm1 905 DO ji = fs_2, fs_jpim1 906 pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1) 907 END DO 908 END DO 909 END DO 910 911 DO jj = 2, jpjm1 !* 3d recurrence: Xk = (Zk - Sk Xk+1 ) / Tk 912 DO ji = fs_2, fs_jpim1 913 pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) 914 END DO 915 END DO 916 DO jk = jpk-2, 2, -1 917 DO jj = 2, jpjm1 918 DO ji = fs_2, fs_jpim1 919 pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk) 920 END DO 921 END DO 922 END DO 923 ! 791 924 END SUBROUTINE interp_4th_cpt 792 925 926 927 SUBROUTINE tridia_solver( pD, pU, pL, pRHS, pt_out , klev ) 928 !!---------------------------------------------------------------------- 929 !! *** ROUTINE tridia_solver *** 930 !! 931 !! ** Purpose : solve a symmetric 3diagonal system 932 !! 933 !! ** Method : solve M.t_out = RHS(t) where M is a tri diagonal matrix ( jpk*jpk ) 934 !! 935 !! ( D_1 U_1 0 0 0 )( t_1 ) ( RHS_1 ) 936 !! ( L_2 D_2 U_2 0 0 )( t_2 ) ( RHS_2 ) 937 !! ( 0 L_3 D_3 U_3 0 )( t_3 ) = ( RHS_3 ) 938 !! ( ... )( ... ) ( ... ) 939 !! ( 0 0 0 L_k D_k )( t_k ) ( RHS_k ) 940 !! 941 !! M is decomposed in the product of an upper and lower triangular matrix. 942 !! The tri-diagonals matrix is given as input 3D arrays: pD, pU, pL 943 !! (i.e. the Diagonal, the Upper diagonal, and the Lower diagonal). 944 !! The solution is pta. 945 !! The 3d array zwt is used as a work space array. 946 !!---------------------------------------------------------------------- 947 REAL(wp),DIMENSION(:,:,:), INTENT(in ) :: pD, pU, PL ! 3-diagonal matrix 948 REAL(wp),DIMENSION(:,:,:), INTENT(in ) :: pRHS ! Right-Hand-Side 949 REAL(wp),DIMENSION(:,:,:), INTENT( out) :: pt_out !!gm field at level=F(klev) 950 INTEGER , INTENT(in ) :: klev ! =1 pt_out at w-level 951 ! ! =0 pt at t-level 952 INTEGER :: ji, jj, jk ! dummy loop integers 953 INTEGER :: kstart ! local indices 954 REAL(wp),DIMENSION(jpi,jpj,jpk) :: zwt ! 3D work array 955 !!---------------------------------------------------------------------- 956 ! 957 kstart = 1 + klev 958 ! 959 DO jj = 2, jpjm1 !* 1st recurrence: Tk = Dk - Ik Sk-1 / Tk-1 960 DO ji = fs_2, fs_jpim1 961 zwt(ji,jj,kstart) = pD(ji,jj,kstart) 962 END DO 963 END DO 964 DO jk = kstart+1, jpkm1 965 DO jj = 2, jpjm1 966 DO ji = fs_2, fs_jpim1 967 zwt(ji,jj,jk) = pD(ji,jj,jk) - pL(ji,jj,jk) * pU(ji,jj,jk-1) /zwt(ji,jj,jk-1) 968 END DO 969 END DO 970 END DO 971 ! 972 DO jj = 2, jpjm1 !* 2nd recurrence: Zk = Yk - Ik / Tk-1 Zk-1 973 DO ji = fs_2, fs_jpim1 974 pt_out(ji,jj,kstart) = pRHS(ji,jj,kstart) 975 END DO 976 END DO 977 DO jk = kstart+1, jpkm1 978 DO jj = 2, jpjm1 979 DO ji = fs_2, fs_jpim1 980 pt_out(ji,jj,jk) = pRHS(ji,jj,jk) - pL(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1) 981 END DO 982 END DO 983 END DO 984 985 DO jj = 2, jpjm1 !* 3d recurrence: Xk = (Zk - Sk Xk+1 ) / Tk 986 DO ji = fs_2, fs_jpim1 987 pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) 988 END DO 989 END DO 990 DO jk = jpk-2, kstart, -1 991 DO jj = 2, jpjm1 992 DO ji = fs_2, fs_jpim1 993 pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - pU(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk) 994 END DO 995 END DO 996 END DO 997 ! 998 END SUBROUTINE tridia_solver 999 793 1000 !!====================================================================== 794 1001 END MODULE traadv_fct
Note: See TracChangeset
for help on using the changeset viewer.
