Changeset 13497 for NEMO/trunk/src/OCE/TRA/traadv_fct.F90
- Timestamp:
- 2020-09-21T14:37:46+02:00 (4 years ago)
- File:
-
- 1 edited
-
NEMO/trunk/src/OCE/TRA/traadv_fct.F90 (modified) (24 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/OCE/TRA/traadv_fct.F90
r13295 r13497 160 160 zwy(ji,jj,jk) = 0.5 * ( zfp_vj * pt(ji,jj,jk,jn,Kbb) + zfm_vj * pt(ji ,jj+1,jk,jn,Kbb) ) 161 161 END_3D 162 ! !* upstream tracer flux in the k direction *!163 DO_3D( 1, 1, 1, 1, 2, jpkm1 ) 162 ! !* upstream tracer flux in the k direction *! 163 DO_3D( 1, 1, 1, 1, 2, jpkm1 ) ! Interior value ( multiplied by wmask) 164 164 zfp_wk = pW(ji,jj,jk) + ABS( pW(ji,jj,jk) ) 165 165 zfm_wk = pW(ji,jj,jk) - ABS( pW(ji,jj,jk) ) 166 166 zwz(ji,jj,jk) = 0.5 * ( zfp_wk * pt(ji,jj,jk,jn,Kbb) + zfm_wk * pt(ji,jj,jk-1,jn,Kbb) ) * wmask(ji,jj,jk) 167 167 END_3D 168 IF( ln_linssh ) THEN ! top ocean value (only in linear free surface as zwz has been w-masked)169 IF( ln_isfcav ) THEN ! top of the ice-shelf cavities and at the ocean surface168 IF( ln_linssh ) THEN ! top ocean value (only in linear free surface as zwz has been w-masked) 169 IF( ln_isfcav ) THEN ! top of the ice-shelf cavities and at the ocean surface 170 170 DO_2D( 1, 1, 1, 1 ) 171 171 zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb) ! linear free surface 172 172 END_2D 173 ELSE ! no cavities: only at the ocean surface173 ELSE ! no cavities: only at the ocean surface 174 174 DO_2D( 1, 1, 1, 1 ) 175 175 zwz(ji,jj,1) = pW(ji,jj,1) * pt(ji,jj,1,jn,Kbb) … … 178 178 ENDIF 179 179 ! 180 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 181 ! ! total intermediate advective trends180 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !* trend and after field with monotonic scheme 181 ! ! total intermediate advective trends 182 182 ztra = - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 183 183 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & 184 184 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) * r1_e1e2t(ji,jj) 185 ! ! update and guess with monotonic sheme185 ! ! update and guess with monotonic sheme 186 186 pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztra & 187 187 & / e3t(ji,jj,jk,Kmm ) * tmask(ji,jj,jk) … … 194 194 ! 195 195 ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp ; 196 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 196 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! Interior value ( multiplied by wmask) 197 197 zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) ) 198 198 zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) ) … … 227 227 zltv(:,:,jpk) = 0._wp 228 228 DO jk = 1, jpkm1 ! Laplacian 229 DO_2D( 1, 0, 1, 0 ) 229 DO_2D( 1, 0, 1, 0 ) ! 1st derivative (gradient) 230 230 ztu(ji,jj,jk) = ( pt(ji+1,jj ,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk) 231 231 ztv(ji,jj,jk) = ( pt(ji ,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk) 232 232 END_2D 233 DO_2D( 0, 0, 0, 0 ) 233 DO_2D( 0, 0, 0, 0 ) ! 2nd derivative * 1/ 6 234 234 zltu(ji,jj,jk) = ( ztu(ji,jj,jk) + ztu(ji-1,jj,jk) ) * r1_6 235 235 zltv(ji,jj,jk) = ( ztv(ji,jj,jk) + ztv(ji,jj-1,jk) ) * r1_6 … … 238 238 CALL lbc_lnk_multi( 'traadv_fct', zltu, 'T', 1.0_wp , zltv, 'T', 1.0_wp ) ! Lateral boundary cond. (unchanged sgn) 239 239 ! 240 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 240 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) ! Horizontal advective fluxes 241 241 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 242 242 zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji ,jj+1,jk,jn,Kmm) 243 ! ! C4 minus upstream advective fluxes243 ! ! C4 minus upstream advective fluxes 244 244 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) 245 245 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) … … 249 249 ztu(:,:,jpk) = 0._wp ! Bottom value : flux set to zero 250 250 ztv(:,:,jpk) = 0._wp 251 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 251 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) ! 1st derivative (gradient) 252 252 ztu(ji,jj,jk) = ( pt(ji+1,jj ,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk) 253 253 ztv(ji,jj,jk) = ( pt(ji ,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk) … … 255 255 CALL lbc_lnk_multi( 'traadv_fct', ztu, 'U', -1.0_wp , ztv, 'V', -1.0_wp ) ! Lateral boundary cond. (unchanged sgn) 256 256 ! 257 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 257 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! Horizontal advective fluxes 258 258 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 (x2) 259 259 zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji ,jj+1,jk,jn,Kmm) … … 288 288 ! 289 289 IF ( ll_zAimp ) THEN 290 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 291 ! ! total intermediate advective trends290 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !* trend and after field with monotonic scheme 291 ! ! total intermediate advective trends 292 292 ztra = - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 293 293 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & … … 298 298 CALL tridia_solver( zwdia, zwsup, zwinf, ztw, ztw , 0 ) 299 299 ! 300 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 300 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! Interior value ( multiplied by wmask) 301 301 zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) ) 302 302 zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) ) … … 324 324 ! 325 325 ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp 326 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 326 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! Interior value ( multiplied by wmask) 327 327 zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) ) 328 328 zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) ) … … 454 454 pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( 1._wp - zcv) * zbv ) 455 455 456 ! monotonic flux in the k direction, i.e. pcc457 ! -------------------------------------------456 ! monotonic flux in the k direction, i.e. pcc 457 ! ------------------------------------------- 458 458 za = MIN( 1., zbetdo(ji,jj,jk+1), zbetup(ji,jj,jk) ) 459 459 zb = MIN( 1., zbetup(ji,jj,jk+1), zbetdo(ji,jj,jk) ) … … 481 481 !!---------------------------------------------------------------------- 482 482 483 DO_3D( 1, 1, 1, 1, 3, jpkm1 ) 483 DO_3D( 1, 1, 1, 1, 3, jpkm1 ) !== build the three diagonal matrix ==! 484 484 zwd (ji,jj,jk) = 4._wp 485 485 zwi (ji,jj,jk) = 1._wp … … 495 495 END_3D 496 496 ! 497 jk = 2 497 jk = 2 ! Switch to second order centered at top 498 498 DO_2D( 1, 1, 1, 1 ) 499 499 zwd (ji,jj,jk) = 1._wp … … 504 504 ! 505 505 ! !== tridiagonal solve ==! 506 DO_2D( 1, 1, 1, 1 ) 506 DO_2D( 1, 1, 1, 1 ) ! first recurrence 507 507 zwt(ji,jj,2) = zwd(ji,jj,2) 508 508 END_2D … … 511 511 END_3D 512 512 ! 513 DO_2D( 1, 1, 1, 1 ) 513 DO_2D( 1, 1, 1, 1 ) ! second recurrence: Zk = Yk - Ik / Tk-1 Zk-1 514 514 pt_out(ji,jj,2) = zwrm(ji,jj,2) 515 515 END_2D … … 518 518 END_3D 519 519 520 DO_2D( 1, 1, 1, 1 ) 520 DO_2D( 1, 1, 1, 1 ) ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk 521 521 pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) 522 522 END_2D … … 546 546 ! !== build the three diagonal matrix & the RHS ==! 547 547 ! 548 DO_3D( 0, 0, 0, 0, 3, jpkm1 ) 548 DO_3D( 0, 0, 0, 0, 3, jpkm1 ) ! interior (from jk=3 to jpk-1) 549 549 zwd (ji,jj,jk) = 3._wp * wmask(ji,jj,jk) + 1._wp ! diagonal 550 550 zwi (ji,jj,jk) = wmask(ji,jj,jk) ! lower diagonal … … 565 565 END IF 566 566 ! 567 DO_2D( 0, 0, 0, 0 ) 567 DO_2D( 0, 0, 0, 0 ) ! 2nd order centered at top & bottom 568 568 ikt = mikt(ji,jj) + 1 ! w-point below the 1st wet point 569 569 ikb = MAX(mbkt(ji,jj), 2) ! - above the last wet point … … 582 582 ! !== tridiagonal solver ==! 583 583 ! 584 DO_2D( 0, 0, 0, 0 ) 584 DO_2D( 0, 0, 0, 0 ) !* 1st recurrence: Tk = Dk - Ik Sk-1 / Tk-1 585 585 zwt(ji,jj,2) = zwd(ji,jj,2) 586 586 END_2D … … 589 589 END_3D 590 590 ! 591 DO_2D( 0, 0, 0, 0 ) 591 DO_2D( 0, 0, 0, 0 ) !* 2nd recurrence: Zk = Yk - Ik / Tk-1 Zk-1 592 592 pt_out(ji,jj,2) = zwrm(ji,jj,2) 593 593 END_2D … … 596 596 END_3D 597 597 598 DO_2D( 0, 0, 0, 0 ) 598 DO_2D( 0, 0, 0, 0 ) !* 3d recurrence: Xk = (Zk - Sk Xk+1 ) / Tk 599 599 pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) 600 600 END_2D … … 638 638 kstart = 1 + klev 639 639 ! 640 DO_2D( 0, 0, 0, 0 ) 640 DO_2D( 0, 0, 0, 0 ) !* 1st recurrence: Tk = Dk - Ik Sk-1 / Tk-1 641 641 zwt(ji,jj,kstart) = pD(ji,jj,kstart) 642 642 END_2D … … 645 645 END_3D 646 646 ! 647 DO_2D( 0, 0, 0, 0 ) 647 DO_2D( 0, 0, 0, 0 ) !* 2nd recurrence: Zk = Yk - Ik / Tk-1 Zk-1 648 648 pt_out(ji,jj,kstart) = pRHS(ji,jj,kstart) 649 649 END_2D … … 652 652 END_3D 653 653 654 DO_2D( 0, 0, 0, 0 ) 654 DO_2D( 0, 0, 0, 0 ) !* 3d recurrence: Xk = (Zk - Sk Xk+1 ) / Tk 655 655 pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) 656 656 END_2D
Note: See TracChangeset
for help on using the changeset viewer.
