- Timestamp:
- 2020-10-06T18:17:44+02:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/TRA/traadv_fct.F90
r13570 r13571 161 161 zwy(ji,jj,jk) = 0.5 * ( zfp_vj * pt(ji,jj,jk,jn,Kbb) + zfm_vj * pt(ji ,jj+1,jk,jn,Kbb) ) 162 162 END_3D 163 ! !* upstream tracer flux in the k direction *!164 DO_3D( 1, 1, 1, 1, 2, jpkm1 ) 163 ! !* upstream tracer flux in the k direction *! 164 DO_3D( 1, 1, 1, 1, 2, jpkm1 ) ! Interior value ( multiplied by wmask) 165 165 zfp_wk = pW(ji,jj,jk) + ABS( pW(ji,jj,jk) ) 166 166 zfm_wk = pW(ji,jj,jk) - ABS( pW(ji,jj,jk) ) 167 167 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) 168 168 END_3D 169 IF( ln_linssh ) THEN ! top ocean value (only in linear free surface as zwz has been w-masked)170 IF( ln_isfcav ) THEN ! top of the ice-shelf cavities and at the ocean surface169 IF( ln_linssh ) THEN ! top ocean value (only in linear free surface as zwz has been w-masked) 170 IF( ln_isfcav ) THEN ! top of the ice-shelf cavities and at the ocean surface 171 171 DO_2D( 1, 1, 1, 1 ) 172 172 zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb) ! linear free surface 173 173 END_2D 174 ELSE ! no cavities: only at the ocean surface174 ELSE ! no cavities: only at the ocean surface 175 175 DO_2D( 1, 1, 1, 1 ) 176 176 zwz(ji,jj,1) = pW(ji,jj,1) * pt(ji,jj,1,jn,Kbb) … … 179 179 ENDIF 180 180 ! 181 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 182 ! ! total intermediate advective trends181 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !* trend and after field with monotonic scheme 182 ! ! total intermediate advective trends 183 183 ztra = - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 184 184 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & 185 185 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) * r1_e1e2t(ji,jj) 186 ! ! update and guess with monotonic sheme186 ! ! update and guess with monotonic sheme 187 187 pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztra & 188 188 & / e3t(ji,jj,jk,Kmm ) * tmask(ji,jj,jk) … … 195 195 ! 196 196 ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp ; 197 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 197 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! Interior value ( multiplied by wmask) 198 198 zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) ) 199 199 zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) ) … … 228 228 zltv(:,:,jpk) = 0._wp 229 229 DO jk = 1, jpkm1 ! Laplacian 230 DO_2D( 1, 0, 1, 0 ) 230 DO_2D( 1, 0, 1, 0 ) ! 1st derivative (gradient) 231 231 ztu(ji,jj,jk) = ( pt(ji+1,jj ,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk) 232 232 ztv(ji,jj,jk) = ( pt(ji ,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk) 233 233 END_2D 234 DO_2D( 0, 0, 0, 0 ) 234 DO_2D( 0, 0, 0, 0 ) ! 2nd derivative * 1/ 6 235 235 zltu(ji,jj,jk) = ( ztu(ji,jj,jk) + ztu(ji-1,jj,jk) ) * r1_6 236 236 zltv(ji,jj,jk) = ( ztv(ji,jj,jk) + ztv(ji,jj-1,jk) ) * r1_6 … … 243 243 #endif 244 244 ! 245 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 245 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) ! Horizontal advective fluxes 246 246 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 247 247 zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji ,jj+1,jk,jn,Kmm) 248 ! ! C4 minus upstream advective fluxes248 ! ! C4 minus upstream advective fluxes 249 249 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) 250 250 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) … … 254 254 ztu(:,:,jpk) = 0._wp ! Bottom value : flux set to zero 255 255 ztv(:,:,jpk) = 0._wp 256 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 256 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) ! 1st derivative (gradient) 257 257 ztu(ji,jj,jk) = ( pt(ji+1,jj ,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk) 258 258 ztv(ji,jj,jk) = ( pt(ji ,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk) … … 264 264 #endif 265 265 ! 266 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 266 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! Horizontal advective fluxes 267 267 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) 268 268 zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji ,jj+1,jk,jn,Kmm) … … 297 297 ! 298 298 IF ( ll_zAimp ) THEN 299 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 300 ! ! total intermediate advective trends299 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !* trend and after field with monotonic scheme 300 ! ! total intermediate advective trends 301 301 ztra = - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 302 302 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & … … 307 307 CALL tridia_solver( zwdia, zwsup, zwinf, ztw, ztw , 0 ) 308 308 ! 309 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 309 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! Interior value ( multiplied by wmask) 310 310 zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) ) 311 311 zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) ) … … 337 337 ! 338 338 ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp 339 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 339 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! Interior value ( multiplied by wmask) 340 340 zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) ) 341 341 zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) ) … … 471 471 pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( 1._wp - zcv) * zbv ) 472 472 473 ! monotonic flux in the k direction, i.e. pcc474 ! -------------------------------------------473 ! monotonic flux in the k direction, i.e. pcc 474 ! ------------------------------------------- 475 475 za = MIN( 1., zbetdo(ji,jj,jk+1), zbetup(ji,jj,jk) ) 476 476 zb = MIN( 1., zbetup(ji,jj,jk+1), zbetdo(ji,jj,jk) ) … … 502 502 !!---------------------------------------------------------------------- 503 503 504 DO_3D( 1, 1, 1, 1, 3, jpkm1 ) 504 DO_3D( 1, 1, 1, 1, 3, jpkm1 ) !== build the three diagonal matrix ==! 505 505 zwd (ji,jj,jk) = 4._wp 506 506 zwi (ji,jj,jk) = 1._wp … … 516 516 END_3D 517 517 ! 518 jk = 2 518 jk = 2 ! Switch to second order centered at top 519 519 DO_2D( 1, 1, 1, 1 ) 520 520 zwd (ji,jj,jk) = 1._wp … … 525 525 ! 526 526 ! !== tridiagonal solve ==! 527 DO_2D( 1, 1, 1, 1 ) 527 DO_2D( 1, 1, 1, 1 ) ! first recurrence 528 528 zwt(ji,jj,2) = zwd(ji,jj,2) 529 529 END_2D … … 532 532 END_3D 533 533 ! 534 DO_2D( 1, 1, 1, 1 ) 534 DO_2D( 1, 1, 1, 1 ) ! second recurrence: Zk = Yk - Ik / Tk-1 Zk-1 535 535 pt_out(ji,jj,2) = zwrm(ji,jj,2) 536 536 END_2D … … 539 539 END_3D 540 540 541 DO_2D( 1, 1, 1, 1 ) 541 DO_2D( 1, 1, 1, 1 ) ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk 542 542 pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) 543 543 END_2D … … 567 567 ! !== build the three diagonal matrix & the RHS ==! 568 568 ! 569 DO_3D( 0, 0, 0, 0, 3, jpkm1 ) 569 DO_3D( 0, 0, 0, 0, 3, jpkm1 ) ! interior (from jk=3 to jpk-1) 570 570 zwd (ji,jj,jk) = 3._wp * wmask(ji,jj,jk) + 1._wp ! diagonal 571 571 zwi (ji,jj,jk) = wmask(ji,jj,jk) ! lower diagonal … … 586 586 END IF 587 587 ! 588 DO_2D( 0, 0, 0, 0 ) 588 DO_2D( 0, 0, 0, 0 ) ! 2nd order centered at top & bottom 589 589 ikt = mikt(ji,jj) + 1 ! w-point below the 1st wet point 590 590 ikb = MAX(mbkt(ji,jj), 2) ! - above the last wet point … … 603 603 ! !== tridiagonal solver ==! 604 604 ! 605 DO_2D( 0, 0, 0, 0 ) 605 DO_2D( 0, 0, 0, 0 ) !* 1st recurrence: Tk = Dk - Ik Sk-1 / Tk-1 606 606 zwt(ji,jj,2) = zwd(ji,jj,2) 607 607 END_2D … … 610 610 END_3D 611 611 ! 612 DO_2D( 0, 0, 0, 0 ) 612 DO_2D( 0, 0, 0, 0 ) !* 2nd recurrence: Zk = Yk - Ik / Tk-1 Zk-1 613 613 pt_out(ji,jj,2) = zwrm(ji,jj,2) 614 614 END_2D … … 617 617 END_3D 618 618 619 DO_2D( 0, 0, 0, 0 ) 619 DO_2D( 0, 0, 0, 0 ) !* 3d recurrence: Xk = (Zk - Sk Xk+1 ) / Tk 620 620 pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) 621 621 END_2D … … 659 659 kstart = 1 + klev 660 660 ! 661 DO_2D( 0, 0, 0, 0 ) 661 DO_2D( 0, 0, 0, 0 ) !* 1st recurrence: Tk = Dk - Ik Sk-1 / Tk-1 662 662 zwt(ji,jj,kstart) = pD(ji,jj,kstart) 663 663 END_2D … … 666 666 END_3D 667 667 ! 668 DO_2D( 0, 0, 0, 0 ) 668 DO_2D( 0, 0, 0, 0 ) !* 2nd recurrence: Zk = Yk - Ik / Tk-1 Zk-1 669 669 pt_out(ji,jj,kstart) = pRHS(ji,jj,kstart) 670 670 END_2D … … 673 673 END_3D 674 674 675 DO_2D( 0, 0, 0, 0 ) 675 DO_2D( 0, 0, 0, 0 ) !* 3d recurrence: Xk = (Zk - Sk Xk+1 ) / Tk 676 676 pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) 677 677 END_2D
Note: See TracChangeset
for help on using the changeset viewer.
