Changeset 14072 for NEMO/trunk/src/OCE/TRA/traadv_fct.F90
- Timestamp:
- 2020-12-04T08:48:38+01:00 (3 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/OCE/TRA/traadv_fct.F90
r13982 r14072 10 10 !! tra_adv_fct : update the tracer trend with a 3D advective trends using a 2nd or 4th order FCT scheme 11 11 !! with sub-time-stepping in the vertical direction 12 !! nonosc : compute monotonic tracer fluxes by a non-oscillatory algorithm 12 !! nonosc : compute monotonic tracer fluxes by a non-oscillatory algorithm 13 13 !! interp_4th_cpt : 4th order compact scheme for the vertical component of the advection 14 14 !!---------------------------------------------------------------------- … … 24 24 ! 25 25 USE in_out_manager ! I/O manager 26 USE iom ! 26 USE iom ! 27 27 USE lib_mpp ! MPP library 28 USE lbclnk ! ocean lateral boundary condition (or mpp link) 29 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 28 USE lbclnk ! ocean lateral boundary condition (or mpp link) 29 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 30 30 31 31 IMPLICIT NONE … … 60 60 !!---------------------------------------------------------------------- 61 61 !! *** ROUTINE tra_adv_fct *** 62 !! 62 !! 63 63 !! ** Purpose : Compute the now trend due to total advection of tracers 64 64 !! and add it to the general trend of tracer equations … … 66 66 !! ** Method : - 2nd or 4th FCT scheme on the horizontal direction 67 67 !! (choice through the value of kn_fct) 68 !! - on the vertical the 4th order is a compact scheme 69 !! - corrected flux (monotonic correction) 68 !! - on the vertical the 4th order is a compact scheme 69 !! - corrected flux (monotonic correction) 70 70 !! 71 71 !! ** Action : - update pt(:,:,:,:,Krhs) with the now advective tracer trends … … 154 154 ! 155 155 ! !== upstream advection with initial mass fluxes & intermediate update ==! 156 ! !* upstream tracer flux in the i and j direction 156 ! !* upstream tracer flux in the i and j direction 157 157 DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 ) 158 158 ! upstream scheme … … 173 173 IF( ln_isfcav ) THEN ! top of the ice-shelf cavities and at the ocean surface 174 174 DO_2D( 1, 1, 1, 1 ) 175 zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb) ! linear free surface 175 zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb) ! linear free surface 176 176 END_2D 177 177 ELSE ! no cavities: only at the ocean surface … … 181 181 ENDIF 182 182 ENDIF 183 ! 183 ! 184 184 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) !* trend and after field with monotonic scheme 185 185 ! ! total intermediate advective trends … … 193 193 & / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk) 194 194 END_3D 195 195 196 196 IF ( ll_zAimp ) THEN 197 197 CALL tridia_solver( zwdia, zwsup, zwinf, zwi, zwi , 0 ) … … 215 215 END IF 216 216 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) 217 IF( l_ptr ) zptry(:,:,:) = zwy(:,:,:) 217 IF( l_ptr ) zptry(:,:,:) = zwy(:,:,:) 218 218 ! 219 219 ! !== anti-diffusive flux : high order minus low order ==! … … 268 268 zC4t_u = zC2t_u + r1_6 * ( ztu(ji-1,jj ,jk) - ztu(ji+1,jj ,jk) ) 269 269 zC4t_v = zC2t_v + r1_6 * ( ztv(ji ,jj-1,jk) - ztv(ji ,jj+1,jk) ) 270 ! ! C4 minus upstream advective fluxes 270 ! ! C4 minus upstream advective fluxes 271 271 zwx(ji,jj,jk) = 0.5_wp * pU(ji,jj,jk) * zC4t_u - zwx(ji,jj,jk) 272 272 zwy(ji,jj,jk) = 0.5_wp * pV(ji,jj,jk) * zC4t_v - zwy(ji,jj,jk) … … 275 275 ! 276 276 END SELECT 277 ! 277 ! 278 278 SELECT CASE( kn_fct_v ) !* vertical anti-diffusive fluxes (w-masked interior values) 279 279 ! … … 384 384 DEALLOCATE( ztrdx, ztrdy, ztrdz ) 385 385 ENDIF 386 IF( l_ptr ) THEN 386 IF( l_ptr ) THEN 387 387 DEALLOCATE( zptry ) 388 388 ENDIF … … 394 394 !!--------------------------------------------------------------------- 395 395 !! *** ROUTINE nonosc *** 396 !! 397 !! ** Purpose : compute monotonic tracer fluxes from the upstream 398 !! scheme and the before field by a nonoscillatory algorithm 396 !! 397 !! ** Purpose : compute monotonic tracer fluxes from the upstream 398 !! scheme and the before field by a nonoscillatory algorithm 399 399 !! 400 400 !! ** Method : ... ??? … … 492 492 !!---------------------------------------------------------------------- 493 493 !! *** ROUTINE interp_4th_cpt_org *** 494 !! 494 !! 495 495 !! ** Purpose : Compute the interpolation of tracer at w-point 496 496 !! … … 503 503 REAL(wp),DIMENSION(jpi,jpj,jpk) :: zwd, zwi, zws, zwrm, zwt 504 504 !!---------------------------------------------------------------------- 505 505 506 506 DO_3D( 1, 1, 1, 1, 3, jpkm1 ) !== build the three diagonal matrix ==! 507 507 zwd (ji,jj,jk) = 4._wp … … 514 514 zwi (ji,jj,jk) = 0._wp 515 515 zws (ji,jj,jk) = 0._wp 516 zwrm(ji,jj,jk) = 0.5 * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) ) 516 zwrm(ji,jj,jk) = 0.5 * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) ) 517 517 ENDIF 518 518 END_3D … … 538 538 END_2D 539 539 DO_3D( 1, 1, 1, 1, 3, jpkm1 ) 540 pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1) 540 pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1) 541 541 END_3D 542 542 … … 547 547 pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk) 548 548 END_3D 549 ! 549 ! 550 550 END SUBROUTINE interp_4th_cpt_org 551 551 552 552 553 553 SUBROUTINE interp_4th_cpt( pt_in, pt_out ) 554 554 !!---------------------------------------------------------------------- 555 555 !! *** ROUTINE interp_4th_cpt *** 556 !! 556 !! 557 557 !! ** Purpose : Compute the interpolation of tracer at w-point 558 558 !! … … 582 582 ! CASE( np_CEN2 ) ! 2nd order centered at top & bottom 583 583 ! END SELECT 584 !!gm 584 !!gm 585 585 ! 586 586 IF ( ln_isfcav ) THEN ! set level two values which may not be set in ISF case … … 600 600 zwi (ji,jj,ikb) = 0._wp 601 601 zws (ji,jj,ikb) = 0._wp 602 zwrm(ji,jj,ikb) = 0.5_wp * ( pt_in(ji,jj,ikb-1) + pt_in(ji,jj,ikb) ) 602 zwrm(ji,jj,ikb) = 0.5_wp * ( pt_in(ji,jj,ikb-1) + pt_in(ji,jj,ikb) ) 603 603 END_2D 604 604 ! … … 616 616 END_2D 617 617 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 3, jpkm1 ) 618 pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1) 618 pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1) 619 619 END_3D 620 620 … … 625 625 pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk) 626 626 END_3D 627 ! 627 ! 628 628 END SUBROUTINE interp_4th_cpt 629 629 … … 632 632 !!---------------------------------------------------------------------- 633 633 !! *** ROUTINE tridia_solver *** 634 !! 634 !! 635 635 !! ** Purpose : solve a symmetric 3diagonal system 636 636 !! 637 637 !! ** Method : solve M.t_out = RHS(t) where M is a tri diagonal matrix ( jpk*jpk ) 638 !! 638 !! 639 639 !! ( D_1 U_1 0 0 0 )( t_1 ) ( RHS_1 ) 640 640 !! ( L_2 D_2 U_2 0 0 )( t_2 ) ( RHS_2 ) … … 642 642 !! ( ... )( ... ) ( ... ) 643 643 !! ( 0 0 0 L_k D_k )( t_k ) ( RHS_k ) 644 !! 644 !! 645 645 !! M is decomposed in the product of an upper and lower triangular matrix. 646 !! The tri-diagonals matrix is given as input 3D arrays: pD, pU, pL 646 !! The tri-diagonals matrix is given as input 3D arrays: pD, pU, pL 647 647 !! (i.e. the Diagonal, the Upper diagonal, and the Lower diagonal). 648 648 !! The solution is pta. … … 672 672 END_2D 673 673 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, kstart+1, jpkm1 ) 674 pt_out(ji,jj,jk) = pRHS(ji,jj,jk) - pL(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1) 674 pt_out(ji,jj,jk) = pRHS(ji,jj,jk) - pL(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1) 675 675 END_3D 676 676
Note: See TracChangeset
for help on using the changeset viewer.