- Timestamp:
- 2020-11-27T17:26:33+01:00 (4 years ago)
- Location:
- NEMO/branches/2020/tickets_icb_1900
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/OCE/TRA/traadv_fct.F90 (modified) (27 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/tickets_icb_1900
- Property svn:externals
-
NEMO/branches/2020/tickets_icb_1900/src/OCE/TRA/traadv_fct.F90
r13237 r13899 139 139 IF( ll_zAimp ) THEN 140 140 ALLOCATE(zwdia(jpi,jpj,jpk), zwinf(jpi,jpj,jpk),zwsup(jpi,jpj,jpk)) 141 DO_3D _00_00(1, jpkm1 )141 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 142 142 zwdia(ji,jj,jk) = 1._wp + p2dt * ( MAX( wi(ji,jj,jk) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) ) & 143 143 & / e3t(ji,jj,jk,Krhs) … … 151 151 ! !== upstream advection with initial mass fluxes & intermediate update ==! 152 152 ! !* upstream tracer flux in the i and j direction 153 DO_3D _10_10(1, jpkm1 )153 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 154 154 ! upstream scheme 155 155 zfp_ui = pU(ji,jj,jk) + ABS( pU(ji,jj,jk) ) … … 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 _11_11( 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 surface170 DO_2D _11_11168 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 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 surface 174 zwz(:,:,1) = pW(:,:,1) * pt(:,:,1,jn,Kbb) 173 ELSE ! no cavities: only at the ocean surface 174 DO_2D( 1, 1, 1, 1 ) 175 zwz(ji,jj,1) = pW(ji,jj,1) * pt(ji,jj,1,jn,Kbb) 176 END_2D 175 177 ENDIF 176 178 ENDIF 177 179 ! 178 DO_3D _00_00( 1, jpkm1 )179 ! ! total intermediate advective trends180 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !* trend and after field with monotonic scheme 181 ! ! total intermediate advective trends 180 182 ztra = - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 181 183 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & 182 184 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) * r1_e1e2t(ji,jj) 183 ! ! update and guess with monotonic sheme185 ! ! update and guess with monotonic sheme 184 186 pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztra & 185 187 & / e3t(ji,jj,jk,Kmm ) * tmask(ji,jj,jk) … … 192 194 ! 193 195 ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp ; 194 DO_3D _00_00( 2, jpkm1)196 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! Interior value ( multiplied by wmask) 195 197 zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) ) 196 198 zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) ) … … 198 200 zwz(ji,jj,jk) = zwz(ji,jj,jk) + ztw(ji,jj,jk) ! update vertical fluxes 199 201 END_3D 200 DO_3D _00_00(1, jpkm1 )202 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 201 203 pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ztw(ji,jj,jk) - ztw(ji ,jj ,jk+1) ) & 202 204 & * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) … … 216 218 ! 217 219 CASE( 2 ) !- 2nd order centered 218 DO_3D _10_10(1, jpkm1 )220 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 219 221 zwx(ji,jj,jk) = 0.5_wp * pU(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj,jk,jn,Kmm) ) - zwx(ji,jj,jk) 220 222 zwy(ji,jj,jk) = 0.5_wp * pV(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj+1,jk,jn,Kmm) ) - zwy(ji,jj,jk) … … 225 227 zltv(:,:,jpk) = 0._wp 226 228 DO jk = 1, jpkm1 ! Laplacian 227 DO_2D _10_10229 DO_2D( 1, 0, 1, 0 ) ! 1st derivative (gradient) 228 230 ztu(ji,jj,jk) = ( pt(ji+1,jj ,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk) 229 231 ztv(ji,jj,jk) = ( pt(ji ,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk) 230 232 END_2D 231 DO_2D _00_00233 DO_2D( 0, 0, 0, 0 ) ! 2nd derivative * 1/ 6 232 234 zltu(ji,jj,jk) = ( ztu(ji,jj,jk) + ztu(ji-1,jj,jk) ) * r1_6 233 235 zltv(ji,jj,jk) = ( ztv(ji,jj,jk) + ztv(ji,jj-1,jk) ) * r1_6 … … 236 238 CALL lbc_lnk_multi( 'traadv_fct', zltu, 'T', 1.0_wp , zltv, 'T', 1.0_wp ) ! Lateral boundary cond. (unchanged sgn) 237 239 ! 238 DO_3D _10_10( 1, jpkm1 )240 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) ! Horizontal advective fluxes 239 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 240 242 zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji ,jj+1,jk,jn,Kmm) 241 ! ! C4 minus upstream advective fluxes243 ! ! C4 minus upstream advective fluxes 242 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) 243 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) … … 247 249 ztu(:,:,jpk) = 0._wp ! Bottom value : flux set to zero 248 250 ztv(:,:,jpk) = 0._wp 249 DO_3D _10_10( 1, jpkm1)251 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) ! 1st derivative (gradient) 250 252 ztu(ji,jj,jk) = ( pt(ji+1,jj ,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk) 251 253 ztv(ji,jj,jk) = ( pt(ji ,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk) … … 253 255 CALL lbc_lnk_multi( 'traadv_fct', ztu, 'U', -1.0_wp , ztv, 'V', -1.0_wp ) ! Lateral boundary cond. (unchanged sgn) 254 256 ! 255 DO_3D _00_00( 1, jpkm1 )257 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! Horizontal advective fluxes 256 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) 257 259 zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji ,jj+1,jk,jn,Kmm) … … 269 271 ! 270 272 CASE( 2 ) !- 2nd order centered 271 DO_3D _00_00(2, jpkm1 )273 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 272 274 zwz(ji,jj,jk) = ( pW(ji,jj,jk) * 0.5_wp * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj,jk-1,jn,Kmm) ) & 273 275 & - zwz(ji,jj,jk) ) * wmask(ji,jj,jk) … … 276 278 CASE( 4 ) !- 4th order COMPACT 277 279 CALL interp_4th_cpt( pt(:,:,:,jn,Kmm) , ztw ) ! zwt = COMPACT interpolation of T at w-point 278 DO_3D _00_00(2, jpkm1 )280 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 279 281 zwz(ji,jj,jk) = ( pW(ji,jj,jk) * ztw(ji,jj,jk) - zwz(ji,jj,jk) ) * wmask(ji,jj,jk) 280 282 END_3D … … 286 288 ! 287 289 IF ( ll_zAimp ) THEN 288 DO_3D _00_00( 1, jpkm1 )289 ! ! total intermediate advective trends290 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !* trend and after field with monotonic scheme 291 ! ! total intermediate advective trends 290 292 ztra = - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 291 293 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & … … 296 298 CALL tridia_solver( zwdia, zwsup, zwinf, ztw, ztw , 0 ) 297 299 ! 298 DO_3D _00_00( 2, jpkm1)300 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! Interior value ( multiplied by wmask) 299 301 zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) ) 300 302 zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) ) … … 311 313 ! !== final trend with corrected fluxes ==! 312 314 ! 313 DO_3D _00_00(1, jpkm1 )315 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 314 316 ztra = - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 315 317 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & … … 322 324 ! 323 325 ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp 324 DO_3D _00_00( 2, jpkm1)326 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! Interior value ( multiplied by wmask) 325 327 zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) ) 326 328 zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) ) … … 328 330 zwz(ji,jj,jk) = zwz(ji,jj,jk) + ztw(ji,jj,jk) ! Update vertical fluxes for trend diagnostic 329 331 END_3D 330 DO_3D _00_00(1, jpkm1 )332 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 331 333 pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ztw(ji,jj,jk) - ztw(ji ,jj ,jk+1) ) & 332 334 & * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) … … 407 409 DO jk = 1, jpkm1 408 410 ikm1 = MAX(jk-1,1) 409 DO_2D _00_00411 DO_2D( 0, 0, 0, 0 ) 410 412 411 413 ! search maximum in neighbourhood … … 441 443 ! 3. monotonic flux in the i & j direction (paa & pbb) 442 444 ! ---------------------------------------- 443 DO_3D _00_00(1, jpkm1 )445 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 444 446 zau = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) ) 445 447 zbu = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) ) … … 452 454 pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( 1._wp - zcv) * zbv ) 453 455 454 ! monotonic flux in the k direction, i.e. pcc455 ! -------------------------------------------456 ! monotonic flux in the k direction, i.e. pcc 457 ! ------------------------------------------- 456 458 za = MIN( 1., zbetdo(ji,jj,jk+1), zbetup(ji,jj,jk) ) 457 459 zb = MIN( 1., zbetup(ji,jj,jk+1), zbetdo(ji,jj,jk) ) … … 479 481 !!---------------------------------------------------------------------- 480 482 481 DO_3D _11_11( 3, jpkm1 )483 DO_3D( 1, 1, 1, 1, 3, jpkm1 ) !== build the three diagonal matrix ==! 482 484 zwd (ji,jj,jk) = 4._wp 483 485 zwi (ji,jj,jk) = 1._wp … … 493 495 END_3D 494 496 ! 495 jk = 2 496 DO_2D _11_11497 jk = 2 ! Switch to second order centered at top 498 DO_2D( 1, 1, 1, 1 ) 497 499 zwd (ji,jj,jk) = 1._wp 498 500 zwi (ji,jj,jk) = 0._wp … … 502 504 ! 503 505 ! !== tridiagonal solve ==! 504 DO_2D _11_11506 DO_2D( 1, 1, 1, 1 ) ! first recurrence 505 507 zwt(ji,jj,2) = zwd(ji,jj,2) 506 508 END_2D 507 DO_3D _11_11(3, jpkm1 )509 DO_3D( 1, 1, 1, 1, 3, jpkm1 ) 508 510 zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1) 509 511 END_3D 510 512 ! 511 DO_2D _11_11513 DO_2D( 1, 1, 1, 1 ) ! second recurrence: Zk = Yk - Ik / Tk-1 Zk-1 512 514 pt_out(ji,jj,2) = zwrm(ji,jj,2) 513 515 END_2D 514 DO_3D _11_11(3, jpkm1 )516 DO_3D( 1, 1, 1, 1, 3, jpkm1 ) 515 517 pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1) 516 518 END_3D 517 519 518 DO_2D _11_11520 DO_2D( 1, 1, 1, 1 ) ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk 519 521 pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) 520 522 END_2D 521 DO_3DS _11_11(jpk-2, 2, -1 )523 DO_3DS( 1, 1, 1, 1, jpk-2, 2, -1 ) 522 524 pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk) 523 525 END_3D … … 544 546 ! !== build the three diagonal matrix & the RHS ==! 545 547 ! 546 DO_3D _00_00( 3, jpkm1)548 DO_3D( 0, 0, 0, 0, 3, jpkm1 ) ! interior (from jk=3 to jpk-1) 547 549 zwd (ji,jj,jk) = 3._wp * wmask(ji,jj,jk) + 1._wp ! diagonal 548 550 zwi (ji,jj,jk) = wmask(ji,jj,jk) ! lower diagonal … … 563 565 END IF 564 566 ! 565 DO_2D _00_00567 DO_2D( 0, 0, 0, 0 ) ! 2nd order centered at top & bottom 566 568 ikt = mikt(ji,jj) + 1 ! w-point below the 1st wet point 567 569 ikb = MAX(mbkt(ji,jj), 2) ! - above the last wet point … … 580 582 ! !== tridiagonal solver ==! 581 583 ! 582 DO_2D _00_00584 DO_2D( 0, 0, 0, 0 ) !* 1st recurrence: Tk = Dk - Ik Sk-1 / Tk-1 583 585 zwt(ji,jj,2) = zwd(ji,jj,2) 584 586 END_2D 585 DO_3D _00_00(3, jpkm1 )587 DO_3D( 0, 0, 0, 0, 3, jpkm1 ) 586 588 zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1) 587 589 END_3D 588 590 ! 589 DO_2D _00_00591 DO_2D( 0, 0, 0, 0 ) !* 2nd recurrence: Zk = Yk - Ik / Tk-1 Zk-1 590 592 pt_out(ji,jj,2) = zwrm(ji,jj,2) 591 593 END_2D 592 DO_3D _00_00(3, jpkm1 )594 DO_3D( 0, 0, 0, 0, 3, jpkm1 ) 593 595 pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1) 594 596 END_3D 595 597 596 DO_2D _00_00598 DO_2D( 0, 0, 0, 0 ) !* 3d recurrence: Xk = (Zk - Sk Xk+1 ) / Tk 597 599 pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) 598 600 END_2D 599 DO_3DS _00_00(jpk-2, 2, -1 )601 DO_3DS( 0, 0, 0, 0, jpk-2, 2, -1 ) 600 602 pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk) 601 603 END_3D … … 636 638 kstart = 1 + klev 637 639 ! 638 DO_2D _00_00640 DO_2D( 0, 0, 0, 0 ) !* 1st recurrence: Tk = Dk - Ik Sk-1 / Tk-1 639 641 zwt(ji,jj,kstart) = pD(ji,jj,kstart) 640 642 END_2D 641 DO_3D _00_00(kstart+1, jpkm1 )643 DO_3D( 0, 0, 0, 0, kstart+1, jpkm1 ) 642 644 zwt(ji,jj,jk) = pD(ji,jj,jk) - pL(ji,jj,jk) * pU(ji,jj,jk-1) /zwt(ji,jj,jk-1) 643 645 END_3D 644 646 ! 645 DO_2D _00_00647 DO_2D( 0, 0, 0, 0 ) !* 2nd recurrence: Zk = Yk - Ik / Tk-1 Zk-1 646 648 pt_out(ji,jj,kstart) = pRHS(ji,jj,kstart) 647 649 END_2D 648 DO_3D _00_00(kstart+1, jpkm1 )650 DO_3D( 0, 0, 0, 0, kstart+1, jpkm1 ) 649 651 pt_out(ji,jj,jk) = pRHS(ji,jj,jk) - pL(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1) 650 652 END_3D 651 653 652 DO_2D _00_00654 DO_2D( 0, 0, 0, 0 ) !* 3d recurrence: Xk = (Zk - Sk Xk+1 ) / Tk 653 655 pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) 654 656 END_2D 655 DO_3DS _00_00(jpk-2, kstart, -1 )657 DO_3DS( 0, 0, 0, 0, jpk-2, kstart, -1 ) 656 658 pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - pU(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk) 657 659 END_3D
Note: See TracChangeset
for help on using the changeset viewer.
