Changeset 13892 for NEMO/branches/2020/dev_r12985_TOP-04_IMMERSE_BGC_interface/src/OCE/TRA/traadv_fct.F90
- Timestamp:
- 2020-11-26T17:47:20+01:00 (3 years ago)
- Location:
- NEMO/branches/2020/dev_r12985_TOP-04_IMMERSE_BGC_interface
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/OCE/TRA/traadv_fct.F90 (modified) (28 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r12985_TOP-04_IMMERSE_BGC_interface
- Property svn:externals
-
old new 3 3 ^/utils/build/mk@HEAD mk 4 4 ^/utils/tools@HEAD tools 5 ^/vendors/AGRIF/dev @HEADext/AGRIF5 ^/vendors/AGRIF/dev_r12970_AGRIF_CMEMS ext/AGRIF 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 8 9 9 # SETTE 10 ^/utils/CI/sette@1 2931sette10 ^/utils/CI/sette@13559 sette
-
- Property svn:externals
-
NEMO/branches/2020/dev_r12985_TOP-04_IMMERSE_BGC_interface/src/OCE/TRA/traadv_fct.F90
r12489 r13892 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) … … 96 97 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 97 98 ENDIF 99 !! -- init to 0 100 zwi(:,:,:) = 0._wp 101 zwx(:,:,:) = 0._wp 102 zwy(:,:,:) = 0._wp 103 zwz(:,:,:) = 0._wp 104 ztu(:,:,:) = 0._wp 105 ztv(:,:,:) = 0._wp 106 zltu(:,:,:) = 0._wp 107 zltv(:,:,:) = 0._wp 108 ztw(:,:,:) = 0._wp 98 109 ! 99 110 l_trd = .FALSE. ! set local switches … … 128 139 IF( ll_zAimp ) THEN 129 140 ALLOCATE(zwdia(jpi,jpj,jpk), zwinf(jpi,jpj,jpk),zwsup(jpi,jpj,jpk)) 130 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) 141 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 142 zwdia(ji,jj,jk) = 1._wp + p2dt * ( MAX( wi(ji,jj,jk) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) ) & 143 & / e3t(ji,jj,jk,Krhs) 132 144 zwinf(ji,jj,jk) = p2dt * MIN( wi(ji,jj,jk ) , 0._wp ) / e3t(ji,jj,jk,Krhs) 133 145 zwsup(ji,jj,jk) = -p2dt * MAX( wi(ji,jj,jk+1) , 0._wp ) / e3t(ji,jj,jk,Krhs) … … 139 151 ! !== upstream advection with initial mass fluxes & intermediate update ==! 140 152 ! !* upstream tracer flux in the i and j direction 141 DO_3D _10_10(1, jpkm1 )153 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 142 154 ! upstream scheme 143 155 zfp_ui = pU(ji,jj,jk) + ABS( pU(ji,jj,jk) ) … … 148 160 zwy(ji,jj,jk) = 0.5 * ( zfp_vj * pt(ji,jj,jk,jn,Kbb) + zfm_vj * pt(ji ,jj+1,jk,jn,Kbb) ) 149 161 END_3D 150 ! !* upstream tracer flux in the k direction *!151 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) 152 164 zfp_wk = pW(ji,jj,jk) + ABS( pW(ji,jj,jk) ) 153 165 zfm_wk = pW(ji,jj,jk) - ABS( pW(ji,jj,jk) ) 154 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) 155 167 END_3D 156 IF( ln_linssh ) THEN ! top ocean value (only in linear free surface as zwz has been w-masked)157 IF( ln_isfcav ) THEN ! top of the ice-shelf cavities and at the ocean surface158 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 ) 159 171 zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb) ! linear free surface 160 172 END_2D 161 ELSE ! no cavities: only at the ocean surface 162 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 163 177 ENDIF 164 178 ENDIF 165 179 ! 166 DO_3D _00_00( 1, jpkm1 )167 ! ! total intermediate advective trends180 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !* trend and after field with monotonic scheme 181 ! ! total intermediate advective trends 168 182 ztra = - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 169 183 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & 170 184 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) * r1_e1e2t(ji,jj) 171 ! ! 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) 185 ! ! update and guess with monotonic sheme 186 pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztra & 187 & / e3t(ji,jj,jk,Kmm ) * tmask(ji,jj,jk) 188 zwi(ji,jj,jk) = ( e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) + p2dt * ztra ) & 189 & / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk) 174 190 END_3D 175 191 … … 178 194 ! 179 195 ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp ; 180 DO_3D _00_00( 2, jpkm1)196 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! Interior value ( multiplied by wmask) 181 197 zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) ) 182 198 zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) ) … … 184 200 zwz(ji,jj,jk) = zwz(ji,jj,jk) + ztw(ji,jj,jk) ! update vertical fluxes 185 201 END_3D 186 DO_3D _00_00(1, jpkm1 )202 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 187 203 pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ztw(ji,jj,jk) - ztw(ji ,jj ,jk+1) ) & 188 204 & * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) … … 202 218 ! 203 219 CASE( 2 ) !- 2nd order centered 204 DO_3D _10_10(1, jpkm1 )220 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 205 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) 206 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) … … 211 227 zltv(:,:,jpk) = 0._wp 212 228 DO jk = 1, jpkm1 ! Laplacian 213 DO_2D _10_10229 DO_2D( 1, 0, 1, 0 ) ! 1st derivative (gradient) 214 230 ztu(ji,jj,jk) = ( pt(ji+1,jj ,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk) 215 231 ztv(ji,jj,jk) = ( pt(ji ,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk) 216 232 END_2D 217 DO_2D _00_00233 DO_2D( 0, 0, 0, 0 ) ! 2nd derivative * 1/ 6 218 234 zltu(ji,jj,jk) = ( ztu(ji,jj,jk) + ztu(ji-1,jj,jk) ) * r1_6 219 235 zltv(ji,jj,jk) = ( ztv(ji,jj,jk) + ztv(ji,jj-1,jk) ) * r1_6 220 236 END_2D 221 237 END DO 222 CALL lbc_lnk_multi( 'traadv_fct', zltu, 'T', 1. , zltv, 'T', 1.) ! Lateral boundary cond. (unchanged sgn)223 ! 224 DO_3D _10_10( 1, jpkm1 )238 CALL lbc_lnk_multi( 'traadv_fct', zltu, 'T', 1.0_wp , zltv, 'T', 1.0_wp ) ! Lateral boundary cond. (unchanged sgn) 239 ! 240 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) ! Horizontal advective fluxes 225 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 226 242 zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji ,jj+1,jk,jn,Kmm) 227 ! ! C4 minus upstream advective fluxes243 ! ! C4 minus upstream advective fluxes 228 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) 229 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) … … 233 249 ztu(:,:,jpk) = 0._wp ! Bottom value : flux set to zero 234 250 ztv(:,:,jpk) = 0._wp 235 DO_3D _10_10( 1, jpkm1)251 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) ! 1st derivative (gradient) 236 252 ztu(ji,jj,jk) = ( pt(ji+1,jj ,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk) 237 253 ztv(ji,jj,jk) = ( pt(ji ,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk) 238 254 END_3D 239 CALL lbc_lnk_multi( 'traadv_fct', ztu, 'U', -1. , ztv, 'V', -1.) ! Lateral boundary cond. (unchanged sgn)240 ! 241 DO_3D _00_00( 1, jpkm1 )255 CALL lbc_lnk_multi( 'traadv_fct', ztu, 'U', -1.0_wp , ztv, 'V', -1.0_wp ) ! Lateral boundary cond. (unchanged sgn) 256 ! 257 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! Horizontal advective fluxes 242 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) 243 259 zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji ,jj+1,jk,jn,Kmm) … … 255 271 ! 256 272 CASE( 2 ) !- 2nd order centered 257 DO_3D _00_00(2, jpkm1 )273 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 258 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) ) & 259 275 & - zwz(ji,jj,jk) ) * wmask(ji,jj,jk) … … 262 278 CASE( 4 ) !- 4th order COMPACT 263 279 CALL interp_4th_cpt( pt(:,:,:,jn,Kmm) , ztw ) ! zwt = COMPACT interpolation of T at w-point 264 DO_3D _00_00(2, jpkm1 )280 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 265 281 zwz(ji,jj,jk) = ( pW(ji,jj,jk) * ztw(ji,jj,jk) - zwz(ji,jj,jk) ) * wmask(ji,jj,jk) 266 282 END_3D … … 272 288 ! 273 289 IF ( ll_zAimp ) THEN 274 DO_3D _00_00( 1, jpkm1 )275 ! ! total intermediate advective trends290 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !* trend and after field with monotonic scheme 291 ! ! total intermediate advective trends 276 292 ztra = - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 277 293 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & … … 282 298 CALL tridia_solver( zwdia, zwsup, zwinf, ztw, ztw , 0 ) 283 299 ! 284 DO_3D _00_00( 2, jpkm1)300 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! Interior value ( multiplied by wmask) 285 301 zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) ) 286 302 zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) ) … … 289 305 END IF 290 306 ! 291 CALL lbc_lnk_multi( 'traadv_fct', zwi, 'T', 1. , zwx, 'U', -1. , zwy, 'V', -1., zwz, 'W', 1.)307 CALL lbc_lnk_multi( 'traadv_fct', zwi, 'T', 1.0_wp, zwx, 'U', -1.0_wp , zwy, 'V', -1.0_wp, zwz, 'W', 1.0_wp ) 292 308 ! 293 309 ! !== monotonicity algorithm ==! … … 297 313 ! !== final trend with corrected fluxes ==! 298 314 ! 299 DO_3D _00_00(1, jpkm1 )315 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 300 316 ztra = - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 301 317 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & … … 308 324 ! 309 325 ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp 310 DO_3D _00_00( 2, jpkm1)326 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! Interior value ( multiplied by wmask) 311 327 zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) ) 312 328 zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) ) … … 314 330 zwz(ji,jj,jk) = zwz(ji,jj,jk) + ztw(ji,jj,jk) ! Update vertical fluxes for trend diagnostic 315 331 END_3D 316 DO_3D _00_00(1, jpkm1 )332 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 317 333 pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ztw(ji,jj,jk) - ztw(ji ,jj ,jk+1) ) & 318 334 & * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) … … 374 390 INTEGER :: ji, jj, jk ! dummy loop indices 375 391 INTEGER :: ikm1 ! local integer 376 REAL( wp) :: zpos, zneg, zbt, za, zb, zc, zbig, zrtrn ! local scalars377 REAL( wp) :: zau, zbu, zcu, zav, zbv, zcv, zup, zdo ! - -378 REAL( wp), DIMENSION(jpi,jpj,jpk) :: zbetup, zbetdo, zbup, zbdo379 !!---------------------------------------------------------------------- 380 ! 381 zbig = 1.e+40_ wp382 zrtrn = 1.e-15_ wp383 zbetup(:,:,:) = 0._ wp ; zbetdo(:,:,:) = 0._wp392 REAL(dp) :: zpos, zneg, zbt, za, zb, zc, zbig, zrtrn ! local scalars 393 REAL(dp) :: zau, zbu, zcu, zav, zbv, zcv, zup, zdo ! - - 394 REAL(dp), DIMENSION(jpi,jpj,jpk) :: zbetup, zbetdo, zbup, zbdo 395 !!---------------------------------------------------------------------- 396 ! 397 zbig = 1.e+40_dp 398 zrtrn = 1.e-15_dp 399 zbetup(:,:,:) = 0._dp ; zbetdo(:,:,:) = 0._dp 384 400 385 401 ! Search local extrema … … 393 409 DO jk = 1, jpkm1 394 410 ikm1 = MAX(jk-1,1) 395 DO_2D _00_00411 DO_2D( 0, 0, 0, 0 ) 396 412 397 413 ! search maximum in neighbourhood … … 423 439 END_2D 424 440 END DO 425 CALL lbc_lnk_multi( 'traadv_fct', zbetup, 'T', 1. , zbetdo, 'T', 1.) ! lateral boundary cond. (unchanged sign)441 CALL lbc_lnk_multi( 'traadv_fct', zbetup, 'T', 1.0_wp , zbetdo, 'T', 1.0_wp ) ! lateral boundary cond. (unchanged sign) 426 442 427 443 ! 3. monotonic flux in the i & j direction (paa & pbb) 428 444 ! ---------------------------------------- 429 DO_3D _00_00(1, jpkm1 )445 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 430 446 zau = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) ) 431 447 zbu = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) ) 432 zcu = ( 0.5 + SIGN( 0.5 , paa(ji,jj,jk) ) )448 zcu = ( 0.5 + SIGN( 0.5_wp , paa(ji,jj,jk) ) ) 433 449 paa(ji,jj,jk) = paa(ji,jj,jk) * ( zcu * zau + ( 1._wp - zcu) * zbu ) 434 450 435 451 zav = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) ) 436 452 zbv = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) ) 437 zcv = ( 0.5 + SIGN( 0.5 , pbb(ji,jj,jk) ) )453 zcv = ( 0.5 + SIGN( 0.5_wp , pbb(ji,jj,jk) ) ) 438 454 pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( 1._wp - zcv) * zbv ) 439 455 440 ! monotonic flux in the k direction, i.e. pcc441 ! -------------------------------------------456 ! monotonic flux in the k direction, i.e. pcc 457 ! ------------------------------------------- 442 458 za = MIN( 1., zbetdo(ji,jj,jk+1), zbetup(ji,jj,jk) ) 443 459 zb = MIN( 1., zbetup(ji,jj,jk+1), zbetdo(ji,jj,jk) ) 444 zc = ( 0.5 + SIGN( 0.5 , pcc(ji,jj,jk+1) ) )460 zc = ( 0.5 + SIGN( 0.5_wp , pcc(ji,jj,jk+1) ) ) 445 461 pcc(ji,jj,jk+1) = pcc(ji,jj,jk+1) * ( zc * za + ( 1._wp - zc) * zb ) 446 462 END_3D 447 CALL lbc_lnk_multi( 'traadv_fct', paa, 'U', -1. , pbb, 'V', -1.) ! lateral boundary condition (changed sign)463 CALL lbc_lnk_multi( 'traadv_fct', paa, 'U', -1.0_wp , pbb, 'V', -1.0_wp ) ! lateral boundary condition (changed sign) 448 464 ! 449 465 END SUBROUTINE nonosc … … 465 481 !!---------------------------------------------------------------------- 466 482 467 DO_3D _11_11( 3, jpkm1 )483 DO_3D( 1, 1, 1, 1, 3, jpkm1 ) !== build the three diagonal matrix ==! 468 484 zwd (ji,jj,jk) = 4._wp 469 485 zwi (ji,jj,jk) = 1._wp … … 479 495 END_3D 480 496 ! 481 jk = 2 482 DO_2D _11_11497 jk = 2 ! Switch to second order centered at top 498 DO_2D( 1, 1, 1, 1 ) 483 499 zwd (ji,jj,jk) = 1._wp 484 500 zwi (ji,jj,jk) = 0._wp … … 488 504 ! 489 505 ! !== tridiagonal solve ==! 490 DO_2D _11_11506 DO_2D( 1, 1, 1, 1 ) ! first recurrence 491 507 zwt(ji,jj,2) = zwd(ji,jj,2) 492 508 END_2D 493 DO_3D _11_11(3, jpkm1 )509 DO_3D( 1, 1, 1, 1, 3, jpkm1 ) 494 510 zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1) 495 511 END_3D 496 512 ! 497 DO_2D _11_11513 DO_2D( 1, 1, 1, 1 ) ! second recurrence: Zk = Yk - Ik / Tk-1 Zk-1 498 514 pt_out(ji,jj,2) = zwrm(ji,jj,2) 499 515 END_2D 500 DO_3D _11_11(3, jpkm1 )516 DO_3D( 1, 1, 1, 1, 3, jpkm1 ) 501 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) 502 518 END_3D 503 519 504 DO_2D _11_11520 DO_2D( 1, 1, 1, 1 ) ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk 505 521 pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) 506 522 END_2D 507 DO_3DS _11_11(jpk-2, 2, -1 )523 DO_3DS( 1, 1, 1, 1, jpk-2, 2, -1 ) 508 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) 509 525 END_3D … … 530 546 ! !== build the three diagonal matrix & the RHS ==! 531 547 ! 532 DO_3D _00_00( 3, jpkm1)548 DO_3D( 0, 0, 0, 0, 3, jpkm1 ) ! interior (from jk=3 to jpk-1) 533 549 zwd (ji,jj,jk) = 3._wp * wmask(ji,jj,jk) + 1._wp ! diagonal 534 550 zwi (ji,jj,jk) = wmask(ji,jj,jk) ! lower diagonal … … 549 565 END IF 550 566 ! 551 DO_2D _00_00567 DO_2D( 0, 0, 0, 0 ) ! 2nd order centered at top & bottom 552 568 ikt = mikt(ji,jj) + 1 ! w-point below the 1st wet point 553 569 ikb = MAX(mbkt(ji,jj), 2) ! - above the last wet point … … 566 582 ! !== tridiagonal solver ==! 567 583 ! 568 DO_2D _00_00584 DO_2D( 0, 0, 0, 0 ) !* 1st recurrence: Tk = Dk - Ik Sk-1 / Tk-1 569 585 zwt(ji,jj,2) = zwd(ji,jj,2) 570 586 END_2D 571 DO_3D _00_00(3, jpkm1 )587 DO_3D( 0, 0, 0, 0, 3, jpkm1 ) 572 588 zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1) 573 589 END_3D 574 590 ! 575 DO_2D _00_00591 DO_2D( 0, 0, 0, 0 ) !* 2nd recurrence: Zk = Yk - Ik / Tk-1 Zk-1 576 592 pt_out(ji,jj,2) = zwrm(ji,jj,2) 577 593 END_2D 578 DO_3D _00_00(3, jpkm1 )594 DO_3D( 0, 0, 0, 0, 3, jpkm1 ) 579 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) 580 596 END_3D 581 597 582 DO_2D _00_00598 DO_2D( 0, 0, 0, 0 ) !* 3d recurrence: Xk = (Zk - Sk Xk+1 ) / Tk 583 599 pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) 584 600 END_2D 585 DO_3DS _00_00(jpk-2, 2, -1 )601 DO_3DS( 0, 0, 0, 0, jpk-2, 2, -1 ) 586 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) 587 603 END_3D … … 622 638 kstart = 1 + klev 623 639 ! 624 DO_2D _00_00640 DO_2D( 0, 0, 0, 0 ) !* 1st recurrence: Tk = Dk - Ik Sk-1 / Tk-1 625 641 zwt(ji,jj,kstart) = pD(ji,jj,kstart) 626 642 END_2D 627 DO_3D _00_00(kstart+1, jpkm1 )643 DO_3D( 0, 0, 0, 0, kstart+1, jpkm1 ) 628 644 zwt(ji,jj,jk) = pD(ji,jj,jk) - pL(ji,jj,jk) * pU(ji,jj,jk-1) /zwt(ji,jj,jk-1) 629 645 END_3D 630 646 ! 631 DO_2D _00_00647 DO_2D( 0, 0, 0, 0 ) !* 2nd recurrence: Zk = Yk - Ik / Tk-1 Zk-1 632 648 pt_out(ji,jj,kstart) = pRHS(ji,jj,kstart) 633 649 END_2D 634 DO_3D _00_00(kstart+1, jpkm1 )650 DO_3D( 0, 0, 0, 0, kstart+1, jpkm1 ) 635 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) 636 652 END_3D 637 653 638 DO_2D _00_00654 DO_2D( 0, 0, 0, 0 ) !* 3d recurrence: Xk = (Zk - Sk Xk+1 ) / Tk 639 655 pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) 640 656 END_2D 641 DO_3DS _00_00(jpk-2, kstart, -1 )657 DO_3DS( 0, 0, 0, 0, jpk-2, kstart, -1 ) 642 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) 643 659 END_3D
Note: See TracChangeset
for help on using the changeset viewer.
