- Timestamp:
- 2020-03-23T22:16:19+01:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/traadv_fct.F90
r12377 r12590 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 … … 46 46 !! * Substitutions 47 47 # include "do_loop_substitute.h90" 48 # include "domzgr_substitute.h90" 48 49 !!---------------------------------------------------------------------- 49 50 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 57 58 !!---------------------------------------------------------------------- 58 59 !! *** ROUTINE tra_adv_fct *** 59 !! 60 !! 60 61 !! ** Purpose : Compute the now trend due to total advection of tracers 61 62 !! and add it to the general trend of tracer equations … … 63 64 !! ** Method : - 2nd or 4th FCT scheme on the horizontal direction 64 65 !! (choice through the value of kn_fct) 65 !! - on the vertical the 4th order is a compact scheme 66 !! - corrected flux (monotonic correction) 66 !! - on the vertical the 4th order is a compact scheme 67 !! - corrected flux (monotonic correction) 67 68 !! 68 69 !! ** Action : - update pt(:,:,:,:,Krhs) with the now advective tracer trends … … 81 82 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! tracers and RHS of tracer equation 82 83 ! 83 INTEGER :: ji, jj, jk, jn ! dummy loop indices 84 INTEGER :: ji, jj, jk, jn ! dummy loop indices 84 85 REAL(wp) :: ztra ! local scalar 85 86 REAL(wp) :: zfp_ui, zfp_vj, zfp_wk, zC2t_u, zC4t_u ! - - … … 102 103 ll_zAimp = .FALSE. 103 104 IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype =='TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 104 IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) ) l_ptr = .TRUE. 105 IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) ) l_ptr = .TRUE. 105 106 IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. & 106 107 & iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) ) l_hst = .TRUE. … … 111 112 ENDIF 112 113 ! 113 IF( l_ptr ) THEN 114 IF( l_ptr ) THEN 114 115 ALLOCATE( zptry(jpi,jpj,jpk) ) 115 116 zptry(:,:,:) = 0._wp 116 117 ENDIF 117 118 ! ! surface & bottom value : flux set to zero one for all 118 zwz(:,:, 1 ) = 0._wp 119 zwz(:,:, 1 ) = 0._wp 119 120 zwx(:,:,jpk) = 0._wp ; zwy(:,:,jpk) = 0._wp ; zwz(:,:,jpk) = 0._wp 120 121 ! 121 zwi(:,:,:) = 0._wp 122 zwi(:,:,:) = 0._wp 122 123 ! 123 124 ! If adaptive vertical advection, check if it is needed on this PE at this time … … 129 130 ALLOCATE(zwdia(jpi,jpj,jpk), zwinf(jpi,jpj,jpk),zwsup(jpi,jpj,jpk)) 130 131 DO_3D_00_00( 1, jpkm1 ) 131 zwdia(ji,jj,jk) = 1._wp + p2dt * ( MAX( wi(ji,jj,jk ) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) ) / e3t(ji,jj,jk,Krhs) 132 zwdia(ji,jj,jk) = 1._wp + p2dt * ( MAX( wi(ji,jj,jk) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) ) & 133 & / e3t(ji,jj,jk,Krhs) 132 134 zwinf(ji,jj,jk) = p2dt * MIN( wi(ji,jj,jk ) , 0._wp ) / e3t(ji,jj,jk,Krhs) 133 135 zwsup(ji,jj,jk) = -p2dt * MAX( wi(ji,jj,jk+1) , 0._wp ) / e3t(ji,jj,jk,Krhs) … … 138 140 ! 139 141 ! !== upstream advection with initial mass fluxes & intermediate update ==! 140 ! !* upstream tracer flux in the i and j direction 142 ! !* upstream tracer flux in the i and j direction 141 143 DO_3D_10_10( 1, jpkm1 ) 142 144 ! upstream scheme … … 157 159 IF( ln_isfcav ) THEN ! top of the ice-shelf cavities and at the ocean surface 158 160 DO_2D_11_11 159 zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb) ! linear free surface 161 zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb) ! linear free surface 160 162 END_2D 161 163 ELSE ! no cavities: only at the ocean surface … … 163 165 ENDIF 164 166 ENDIF 165 ! 167 ! 166 168 DO_3D_00_00( 1, jpkm1 ) 167 169 ! ! total intermediate advective trends … … 170 172 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) * r1_e1e2t(ji,jj) 171 173 ! ! update and guess with monotonic sheme 172 pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztra / e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk) 173 zwi(ji,jj,jk) = ( e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) + p2dt * ztra ) / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk) 174 pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztra & 175 & / e3t(ji,jj,jk,Kmm ) * tmask(ji,jj,jk) 176 zwi(ji,jj,jk) = ( e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) + p2dt * ztra ) & 177 & / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk) 174 178 END_3D 175 179 176 180 IF ( ll_zAimp ) THEN 177 181 CALL tridia_solver( zwdia, zwsup, zwinf, zwi, zwi , 0 ) … … 186 190 DO_3D_00_00( 1, jpkm1 ) 187 191 pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ztw(ji,jj,jk) - ztw(ji ,jj ,jk+1) ) & 188 & 192 & * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 189 193 END_3D 190 194 ! 191 195 END IF 192 ! 196 ! 193 197 IF( l_trd .OR. l_hst ) THEN ! trend diagnostics (contribution of upstream fluxes) 194 198 ztrdx(:,:,:) = zwx(:,:,:) ; ztrdy(:,:,:) = zwy(:,:,:) ; ztrdz(:,:,:) = zwz(:,:,:) 195 199 END IF 196 200 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) 197 IF( l_ptr ) zptry(:,:,:) = zwy(:,:,:) 201 IF( l_ptr ) zptry(:,:,:) = zwy(:,:,:) 198 202 ! 199 203 ! !== anti-diffusive flux : high order minus low order ==! … … 225 229 zC2t_u = pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj ,jk,jn,Kmm) ! 2 x C2 interpolation of T at u- & v-points 226 230 zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji ,jj+1,jk,jn,Kmm) 227 ! ! C4 minus upstream advective fluxes 231 ! ! C4 minus upstream advective fluxes 228 232 zwx(ji,jj,jk) = 0.5_wp * pU(ji,jj,jk) * ( zC2t_u + zltu(ji,jj,jk) - zltu(ji+1,jj,jk) ) - zwx(ji,jj,jk) 229 233 zwy(ji,jj,jk) = 0.5_wp * pV(ji,jj,jk) * ( zC2t_v + zltv(ji,jj,jk) - zltv(ji,jj+1,jk) ) - zwy(ji,jj,jk) … … 245 249 zC4t_u = zC2t_u + r1_6 * ( ztu(ji-1,jj ,jk) - ztu(ji+1,jj ,jk) ) 246 250 zC4t_v = zC2t_v + r1_6 * ( ztv(ji ,jj-1,jk) - ztv(ji ,jj+1,jk) ) 247 ! ! C4 minus upstream advective fluxes 251 ! ! C4 minus upstream advective fluxes 248 252 zwx(ji,jj,jk) = 0.5_wp * pU(ji,jj,jk) * zC4t_u - zwx(ji,jj,jk) 249 253 zwy(ji,jj,jk) = 0.5_wp * pV(ji,jj,jk) * zC4t_v - zwy(ji,jj,jk) … … 251 255 ! 252 256 END SELECT 253 ! 257 ! 254 258 SELECT CASE( kn_fct_v ) !* vertical anti-diffusive fluxes (w-masked interior values) 255 259 ! … … 270 274 zwz(:,:,1) = 0._wp ! only ocean surface as interior zwz values have been w-masked 271 275 ENDIF 272 ! 276 ! 273 277 IF ( ll_zAimp ) THEN 274 278 DO_3D_00_00( 1, jpkm1 ) … … 277 281 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & 278 282 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) * r1_e1e2t(ji,jj) 279 ztw(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t(ji,jj,jk,Krhs) *tmask(ji,jj,jk)283 ztw(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t(ji,jj,jk,Krhs)*tmask(ji,jj,jk) 280 284 END_3D 281 285 ! … … 316 320 DO_3D_00_00( 1, jpkm1 ) 317 321 pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ztw(ji,jj,jk) - ztw(ji ,jj ,jk+1) ) & 318 & 319 END_3D 320 END IF 322 & * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 323 END_3D 324 END IF 321 325 ! 322 326 IF( l_trd .OR. l_hst ) THEN ! trend diagnostics // heat/salt transport 323 ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:) ! <<< add anti-diffusive fluxes 327 ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:) ! <<< add anti-diffusive fluxes 324 328 ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:) ! to upstream fluxes 325 329 ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:) ! … … 344 348 DEALLOCATE( zwdia, zwinf, zwsup ) 345 349 ENDIF 346 IF( l_trd .OR. l_hst ) THEN 350 IF( l_trd .OR. l_hst ) THEN 347 351 DEALLOCATE( ztrdx, ztrdy, ztrdz ) 348 352 ENDIF 349 IF( l_ptr ) THEN 353 IF( l_ptr ) THEN 350 354 DEALLOCATE( zptry ) 351 355 ENDIF … … 357 361 !!--------------------------------------------------------------------- 358 362 !! *** ROUTINE nonosc *** 359 !! 360 !! ** Purpose : compute monotonic tracer fluxes from the upstream 361 !! scheme and the before field by a nonoscillatory algorithm 363 !! 364 !! ** Purpose : compute monotonic tracer fluxes from the upstream 365 !! scheme and the before field by a nonoscillatory algorithm 362 366 !! 363 367 !! ** Method : ... ??? … … 367 371 !! in-space based differencing for fluid 368 372 !!---------------------------------------------------------------------- 369 INTEGER , INTENT(in ) :: Kmm ! time level index 373 INTEGER , INTENT(in ) :: Kmm ! time level index 370 374 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 371 375 REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(in ) :: pbef, paft ! before & after field … … 453 457 !!---------------------------------------------------------------------- 454 458 !! *** ROUTINE interp_4th_cpt_org *** 455 !! 459 !! 456 460 !! ** Purpose : Compute the interpolation of tracer at w-point 457 461 !! … … 464 468 REAL(wp),DIMENSION(jpi,jpj,jpk) :: zwd, zwi, zws, zwrm, zwt 465 469 !!---------------------------------------------------------------------- 466 470 467 471 DO_3D_11_11( 3, jpkm1 ) 468 472 zwd (ji,jj,jk) = 4._wp … … 475 479 zwi (ji,jj,jk) = 0._wp 476 480 zws (ji,jj,jk) = 0._wp 477 zwrm(ji,jj,jk) = 0.5 * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) ) 481 zwrm(ji,jj,jk) = 0.5 * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) ) 478 482 ENDIF 479 483 END_3D … … 499 503 END_2D 500 504 DO_3D_11_11( 3, jpkm1 ) 501 pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1) 505 pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1) 502 506 END_3D 503 507 … … 508 512 pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk) 509 513 END_3D 510 ! 514 ! 511 515 END SUBROUTINE interp_4th_cpt_org 512 516 513 517 514 518 SUBROUTINE interp_4th_cpt( pt_in, pt_out ) 515 519 !!---------------------------------------------------------------------- 516 520 !! *** ROUTINE interp_4th_cpt *** 517 !! 521 !! 518 522 !! ** Purpose : Compute the interpolation of tracer at w-point 519 523 !! … … 543 547 ! CASE( np_CEN2 ) ! 2nd order centered at top & bottom 544 548 ! END SELECT 545 !!gm 549 !!gm 546 550 ! 547 551 IF ( ln_isfcav ) THEN ! set level two values which may not be set in ISF case … … 561 565 zwi (ji,jj,ikb) = 0._wp 562 566 zws (ji,jj,ikb) = 0._wp 563 zwrm(ji,jj,ikb) = 0.5_wp * ( pt_in(ji,jj,ikb-1) + pt_in(ji,jj,ikb) ) 567 zwrm(ji,jj,ikb) = 0.5_wp * ( pt_in(ji,jj,ikb-1) + pt_in(ji,jj,ikb) ) 564 568 END_2D 565 569 ! … … 577 581 END_2D 578 582 DO_3D_00_00( 3, jpkm1 ) 579 pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1) 583 pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1) 580 584 END_3D 581 585 … … 586 590 pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk) 587 591 END_3D 588 ! 592 ! 589 593 END SUBROUTINE interp_4th_cpt 590 594 … … 593 597 !!---------------------------------------------------------------------- 594 598 !! *** ROUTINE tridia_solver *** 595 !! 599 !! 596 600 !! ** Purpose : solve a symmetric 3diagonal system 597 601 !! 598 602 !! ** Method : solve M.t_out = RHS(t) where M is a tri diagonal matrix ( jpk*jpk ) 599 !! 603 !! 600 604 !! ( D_1 U_1 0 0 0 )( t_1 ) ( RHS_1 ) 601 605 !! ( L_2 D_2 U_2 0 0 )( t_2 ) ( RHS_2 ) … … 603 607 !! ( ... )( ... ) ( ... ) 604 608 !! ( 0 0 0 L_k D_k )( t_k ) ( RHS_k ) 605 !! 609 !! 606 610 !! M is decomposed in the product of an upper and lower triangular matrix. 607 !! The tri-diagonals matrix is given as input 3D arrays: pD, pU, pL 611 !! The tri-diagonals matrix is given as input 3D arrays: pD, pU, pL 608 612 !! (i.e. the Diagonal, the Upper diagonal, and the Lower diagonal). 609 613 !! The solution is pta. … … 613 617 REAL(wp),DIMENSION(:,:,:), INTENT(in ) :: pRHS ! Right-Hand-Side 614 618 REAL(wp),DIMENSION(:,:,:), INTENT( out) :: pt_out !!gm field at level=F(klev) 615 INTEGER , INTENT(in ) :: klev ! =1 pt_out at w-level 619 INTEGER , INTENT(in ) :: klev ! =1 pt_out at w-level 616 620 ! ! =0 pt at t-level 617 621 INTEGER :: ji, jj, jk ! dummy loop integers … … 633 637 END_2D 634 638 DO_3D_00_00( kstart+1, jpkm1 ) 635 pt_out(ji,jj,jk) = pRHS(ji,jj,jk) - pL(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1) 639 pt_out(ji,jj,jk) = pRHS(ji,jj,jk) - pL(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1) 636 640 END_3D 637 641
Note: See TracChangeset
for help on using the changeset viewer.
