- Timestamp:
- 2021-11-28T18:59:49+01:00 (3 years ago)
- Location:
- NEMO/branches/2021/ticket2632_r14588_theta_sbcblk
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/OCE/TRA/traadv_fct.F90 (modified) (17 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2021/ticket2632_r14588_theta_sbcblk
- Property svn:externals
-
old new 9 9 10 10 # SETTE 11 ^/utils/CI/sette@14244 sette 11 ^/utils/CI/sette@HEAD sette 12
-
- Property svn:externals
-
NEMO/branches/2021/ticket2632_r14588_theta_sbcblk/src/OCE/TRA/traadv_fct.F90
r14433 r15548 34 34 PUBLIC tra_adv_fct ! called by traadv.F90 35 35 PUBLIC interp_4th_cpt ! called by traadv_cen.F90 36 PUBLIC tridia_solver ! called by traadv_fct_lf.F9037 PUBLIC nonosc ! called by traadv_fct_lf.F90 - key_agrif38 36 39 37 LOGICAL :: l_trd ! flag to compute trends … … 81 79 INTEGER , INTENT(in ) :: kn_fct_v ! order of the FCT scheme (=2 or 4) 82 80 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 83 ! TEMP: [tiling] This can be A2D(nn_hls) if using XIOS (subdomain support)81 ! TEMP: [tiling] This can be A2D(nn_hls) after all lbc_lnks removed in the nn_hls = 2 case 84 82 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pU, pV, pW ! 3 ocean volume flux components 85 83 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! tracers and RHS of tracer equation … … 95 93 !!---------------------------------------------------------------------- 96 94 ! 97 IF( ntile == 0 .OR. ntile == 1 ) THEN ! Do only on the first tile 95 #if defined key_loop_fusion 96 CALL tra_adv_fct_lf ( kt, nit000, cdtype, p2dt, pU, pV, pW, Kbb, Kmm, pt, kjpt, Krhs, kn_fct_h, kn_fct_v ) 97 #else 98 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 98 99 IF( kt == kit000 ) THEN 99 100 IF(lwp) WRITE(numout,*) … … 136 137 ! If adaptive vertical advection, check if it is needed on this PE at this time 137 138 IF( ln_zad_Aimp ) THEN 138 IF( MAXVAL( ABS( wi(A2D( nn_hls),:) ) ) > 0._wp ) ll_zAimp = .TRUE.139 IF( MAXVAL( ABS( wi(A2D(1),:) ) ) > 0._wp ) ll_zAimp = .TRUE. 139 140 END IF 140 141 ! If active adaptive vertical advection, build tridiagonal matrix … … 162 163 zwy(ji,jj,jk) = 0.5 * ( zfp_vj * pt(ji,jj,jk,jn,Kbb) + zfm_vj * pt(ji ,jj+1,jk,jn,Kbb) ) 163 164 END_3D 165 ! !* upstream tracer flux in the k direction *! 166 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! Interior value ( multiplied by wmask) 167 zfp_wk = pW(ji,jj,jk) + ABS( pW(ji,jj,jk) ) 168 zfm_wk = pW(ji,jj,jk) - ABS( pW(ji,jj,jk) ) 169 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) 170 END_3D 171 IF( ln_linssh ) THEN ! top ocean value (only in linear free surface as zwz has been w-masked) 172 IF( ln_isfcav ) THEN ! top of the ice-shelf cavities and at the ocean surface 173 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 174 zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb) ! linear free surface 175 END_2D 176 ELSE ! no cavities: only at the ocean surface 177 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 178 zwz(ji,jj,1) = pW(ji,jj,1) * pt(ji,jj,1,jn,Kbb) 179 END_2D 180 ENDIF 181 ENDIF 182 ! 183 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) !* trend and after field with monotonic scheme 184 ! ! total intermediate advective trends 185 ztra = - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 186 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & 187 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) * r1_e1e2t(ji,jj) 188 ! ! update and guess with monotonic sheme 189 pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztra & 190 & / e3t(ji,jj,jk,Kmm ) * tmask(ji,jj,jk) 191 zwi(ji,jj,jk) = ( e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) + p2dt * ztra ) & 192 & / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk) 193 END_3D 194 195 IF ( ll_zAimp ) THEN 196 CALL tridia_solver( zwdia, zwsup, zwinf, zwi, zwi , 0 ) 197 ! 198 ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp ; 199 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! Interior value ( multiplied by wmask) 200 zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) ) 201 zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) ) 202 ztw(ji,jj,jk) = 0.5 * e1e2t(ji,jj) * ( zfp_wk * zwi(ji,jj,jk) + zfm_wk * zwi(ji,jj,jk-1) ) * wmask(ji,jj,jk) 203 zwz(ji,jj,jk) = zwz(ji,jj,jk) + ztw(ji,jj,jk) ! update vertical fluxes 204 END_3D 205 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 206 pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ztw(ji,jj,jk) - ztw(ji ,jj ,jk+1) ) & 207 & * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 208 END_3D 209 ! 210 END IF 211 ! 212 IF( l_trd .OR. l_hst ) THEN ! trend diagnostics (contribution of upstream fluxes) 213 ztrdx(:,:,:) = zwx(:,:,:) ; ztrdy(:,:,:) = zwy(:,:,:) ; ztrdz(:,:,:) = zwz(:,:,:) 214 END IF 215 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) 216 IF( l_ptr ) zptry(:,:,:) = zwy(:,:,:) 217 ! 218 ! !== anti-diffusive flux : high order minus low order ==! 219 ! 220 SELECT CASE( kn_fct_h ) !* horizontal anti-diffusive fluxes 221 ! 222 CASE( 2 ) !- 2nd order centered 223 DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 ) 224 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) 225 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) 226 END_3D 227 ! 228 CASE( 4 ) !- 4th order centered 229 zltu(:,:,jpk) = 0._wp ! Bottom value : flux set to zero 230 zltv(:,:,jpk) = 0._wp 231 DO jk = 1, jpkm1 ! Laplacian 232 DO_2D( 1, 0, 1, 0 ) ! 1st derivative (gradient) 233 ztu(ji,jj,jk) = ( pt(ji+1,jj ,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk) 234 ztv(ji,jj,jk) = ( pt(ji ,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk) 235 END_2D 236 DO_2D( 0, 0, 0, 0 ) ! 2nd derivative * 1/ 6 237 zltu(ji,jj,jk) = ( ztu(ji,jj,jk) + ztu(ji-1,jj,jk) ) * r1_6 238 zltv(ji,jj,jk) = ( ztv(ji,jj,jk) + ztv(ji,jj-1,jk) ) * r1_6 239 END_2D 240 END DO 241 ! NOTE [ comm_cleanup ] : need to change sign to ensure halo 1 - halo 2 compatibility 242 CALL lbc_lnk( 'traadv_fct', zltu, 'T', -1.0_wp , zltv, 'T', -1.0_wp, ld4only= .TRUE. ) ! Lateral boundary cond. (unchanged sgn) 243 ! 244 DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 ) 245 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 246 zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji ,jj+1,jk,jn,Kmm) 247 ! ! C4 minus upstream advective fluxes 248 ! round brackets added to fix the order of floating point operations 249 ! needed to ensure halo 1 - halo 2 compatibility 250 zwx(ji,jj,jk) = 0.5_wp * pU(ji,jj,jk) * ( zC2t_u + ( zltu(ji,jj,jk) - zltu(ji+1,jj,jk) & 251 & ) & ! bracket for halo 1 - halo 2 compatibility 252 & ) - zwx(ji,jj,jk) 253 zwy(ji,jj,jk) = 0.5_wp * pV(ji,jj,jk) * ( zC2t_v + ( zltv(ji,jj,jk) - zltv(ji,jj+1,jk) & 254 & ) & ! bracket for halo 1 - halo 2 compatibility 255 & ) - zwy(ji,jj,jk) 256 END_3D 257 ! 258 CASE( 41 ) !- 4th order centered ==>> !!gm coding attempt need to be tested 259 ztu(:,:,jpk) = 0._wp ! Bottom value : flux set to zero 260 ztv(:,:,jpk) = 0._wp 261 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) ! 1st derivative (gradient) 262 ztu(ji,jj,jk) = ( pt(ji+1,jj ,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk) 263 ztv(ji,jj,jk) = ( pt(ji ,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk) 264 END_3D 265 IF (nn_hls==1) CALL lbc_lnk( 'traadv_fct', ztu, 'U', -1.0_wp , ztv, 'V', -1.0_wp, ld4only= .TRUE. ) ! Lateral boundary cond. (unchanged sgn) 266 ! 267 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! Horizontal advective fluxes 268 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) 269 zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji ,jj+1,jk,jn,Kmm) 270 ! ! C4 interpolation of T at u- & v-points (x2) 271 zC4t_u = zC2t_u + r1_6 * ( ztu(ji-1,jj ,jk) - ztu(ji+1,jj ,jk) ) 272 zC4t_v = zC2t_v + r1_6 * ( ztv(ji ,jj-1,jk) - ztv(ji ,jj+1,jk) ) 273 ! ! C4 minus upstream advective fluxes 274 zwx(ji,jj,jk) = 0.5_wp * pU(ji,jj,jk) * zC4t_u - zwx(ji,jj,jk) 275 zwy(ji,jj,jk) = 0.5_wp * pV(ji,jj,jk) * zC4t_v - zwy(ji,jj,jk) 276 END_3D 277 IF (nn_hls==2) CALL lbc_lnk( 'traadv_fct', zwx, 'U', -1.0_wp , zwy, 'V', -1.0_wp ) ! Lateral boundary cond. (unchanged sgn) 278 ! 279 END SELECT 280 ! 281 SELECT CASE( kn_fct_v ) !* vertical anti-diffusive fluxes (w-masked interior values) 282 ! 283 CASE( 2 ) !- 2nd order centered 284 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 285 zwz(ji,jj,jk) = ( pW(ji,jj,jk) * 0.5_wp * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj,jk-1,jn,Kmm) ) & 286 & - zwz(ji,jj,jk) ) * wmask(ji,jj,jk) 287 END_3D 288 ! 289 CASE( 4 ) !- 4th order COMPACT 290 CALL interp_4th_cpt( pt(:,:,:,jn,Kmm) , ztw ) ! zwt = COMPACT interpolation of T at w-point 291 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 292 zwz(ji,jj,jk) = ( pW(ji,jj,jk) * ztw(ji,jj,jk) - zwz(ji,jj,jk) ) * wmask(ji,jj,jk) 293 END_3D 294 ! 295 END SELECT 296 IF( ln_linssh ) THEN ! top ocean value: high order = upstream ==>> zwz=0 297 zwz(:,:,1) = 0._wp ! only ocean surface as interior zwz values have been w-masked 298 ENDIF 299 ! 300 IF (nn_hls==1) THEN 301 CALL lbc_lnk( 'traadv_fct', zwi, 'T', 1.0_wp, zwx, 'U', -1.0_wp , zwy, 'V', -1.0_wp, zwz, 'T', 1.0_wp ) 302 ELSE 303 CALL lbc_lnk( 'traadv_fct', zwi, 'T', 1.0_wp) 304 END IF 305 ! 306 IF ( ll_zAimp ) THEN 307 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) !* trend and after field with monotonic scheme 308 ! ! total intermediate advective trends 309 ztra = - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 310 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & 311 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) * r1_e1e2t(ji,jj) 312 ztw(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk) 313 END_3D 314 ! 315 CALL tridia_solver( zwdia, zwsup, zwinf, ztw, ztw , 0 ) 316 ! 317 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! Interior value ( multiplied by wmask) 318 zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) ) 319 zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) ) 320 zwz(ji,jj,jk) = zwz(ji,jj,jk) + 0.5 * e1e2t(ji,jj) * ( zfp_wk * ztw(ji,jj,jk) + zfm_wk * ztw(ji,jj,jk-1) ) * wmask(ji,jj,jk) 321 END_3D 322 END IF 323 ! 324 ! !== monotonicity algorithm ==! 325 ! 326 CALL nonosc( Kmm, pt(:,:,:,jn,Kbb), zwx, zwy, zwz, zwi, p2dt ) 327 ! 328 ! !== final trend with corrected fluxes ==! 329 ! 330 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 331 ztra = - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 332 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & 333 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) * r1_e1e2t(ji,jj) 334 pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztra / e3t(ji,jj,jk,Kmm) 335 zwi(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk) 336 END_3D 337 ! 338 IF ( ll_zAimp ) THEN 339 ! 340 ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp 341 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! Interior value ( multiplied by wmask) 342 zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) ) 343 zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) ) 344 ztw(ji,jj,jk) = - 0.5 * e1e2t(ji,jj) * ( zfp_wk * zwi(ji,jj,jk) + zfm_wk * zwi(ji,jj,jk-1) ) * wmask(ji,jj,jk) 345 zwz(ji,jj,jk) = zwz(ji,jj,jk) + ztw(ji,jj,jk) ! Update vertical fluxes for trend diagnostic 346 END_3D 347 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 348 pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ztw(ji,jj,jk) - ztw(ji ,jj ,jk+1) ) & 349 & * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 350 END_3D 351 END IF 352 IF( l_trd .OR. l_hst ) THEN ! trend diagnostics // heat/salt transport 353 ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:) ! <<< add anti-diffusive fluxes 354 ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:) ! to upstream fluxes 355 ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:) ! 356 ! 357 IF( l_trd ) THEN ! trend diagnostics 358 CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_xad, ztrdx, pU, pt(:,:,:,jn,Kmm) ) 359 CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, ztrdy, pV, pt(:,:,:,jn,Kmm) ) 360 CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, ztrdz, pW, pt(:,:,:,jn,Kmm) ) 361 ENDIF 362 ! ! heat/salt transport 363 IF( l_hst ) CALL dia_ar5_hst( jn, 'adv', ztrdx(:,:,:), ztrdy(:,:,:) ) 364 ! 365 ENDIF 366 IF( l_ptr ) THEN ! "Poleward" transports 367 zptry(:,:,:) = zptry(:,:,:) + zwy(:,:,:) ! <<< add anti-diffusive fluxes 368 CALL dia_ptr_hst( jn, 'adv', zptry(:,:,:) ) 369 ENDIF 370 ! 371 END DO ! end of tracer loop 372 ! 373 IF ( ll_zAimp ) THEN 374 DEALLOCATE( zwdia, zwinf, zwsup ) 375 ENDIF 376 IF( l_trd .OR. l_hst ) THEN 377 DEALLOCATE( ztrdx, ztrdy, ztrdz ) 378 ENDIF 379 IF( l_ptr ) THEN 380 DEALLOCATE( zptry ) 381 ENDIF 382 ! 383 #endif 384 END SUBROUTINE tra_adv_fct 385 386 387 SUBROUTINE nonosc( Kmm, pbef, paa, pbb, pcc, paft, p2dt ) 388 !!--------------------------------------------------------------------- 389 !! *** ROUTINE nonosc *** 390 !! 391 !! ** Purpose : compute monotonic tracer fluxes from the upstream 392 !! scheme and the before field by a nonoscillatory algorithm 393 !! 394 !! ** Method : ... ??? 395 !! warning : pbef and paft must be masked, but the boundaries 396 !! conditions on the fluxes are not necessary zalezak (1979) 397 !! drange (1995) multi-dimensional forward-in-time and upstream- 398 !! in-space based differencing for fluid 399 !!---------------------------------------------------------------------- 400 INTEGER , INTENT(in ) :: Kmm ! time level index 401 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 402 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pbef ! before field 403 REAL(wp), DIMENSION(A2D(nn_hls) ,jpk), INTENT(in ) :: paft ! after field 404 REAL(wp), DIMENSION(A2D(nn_hls) ,jpk), INTENT(inout) :: paa, pbb, pcc ! monotonic fluxes in the 3 directions 405 ! 406 INTEGER :: ji, jj, jk ! dummy loop indices 407 INTEGER :: ikm1 ! local integer 408 REAL(dp) :: zpos, zneg, zbt, za, zb, zc, zbig, zrtrn ! local scalars 409 REAL(dp) :: zau, zbu, zcu, zav, zbv, zcv, zup, zdo ! - - 410 REAL(dp), DIMENSION(A2D(nn_hls),jpk) :: zbetup, zbetdo, zbup, zbdo 411 !!---------------------------------------------------------------------- 412 ! 413 zbig = 1.e+40_dp 414 zrtrn = 1.e-15_dp 415 zbetup(:,:,:) = 0._dp ; zbetdo(:,:,:) = 0._dp 416 417 ! Search local extrema 418 ! -------------------- 419 ! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land 420 DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk ) 421 zbup(ji,jj,jk) = MAX( pbef(ji,jj,jk) * tmask(ji,jj,jk) - zbig * ( 1._wp - tmask(ji,jj,jk) ), & 422 & paft(ji,jj,jk) * tmask(ji,jj,jk) - zbig * ( 1._wp - tmask(ji,jj,jk) ) ) 423 zbdo(ji,jj,jk) = MIN( pbef(ji,jj,jk) * tmask(ji,jj,jk) + zbig * ( 1._wp - tmask(ji,jj,jk) ), & 424 & paft(ji,jj,jk) * tmask(ji,jj,jk) + zbig * ( 1._wp - tmask(ji,jj,jk) ) ) 425 END_3D 426 427 DO jk = 1, jpkm1 428 ikm1 = MAX(jk-1,1) 429 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 430 431 ! search maximum in neighbourhood 432 zup = MAX( zbup(ji ,jj ,jk ), & 433 & zbup(ji-1,jj ,jk ), zbup(ji+1,jj ,jk ), & 434 & zbup(ji ,jj-1,jk ), zbup(ji ,jj+1,jk ), & 435 & zbup(ji ,jj ,ikm1), zbup(ji ,jj ,jk+1) ) 436 437 ! search minimum in neighbourhood 438 zdo = MIN( zbdo(ji ,jj ,jk ), & 439 & zbdo(ji-1,jj ,jk ), zbdo(ji+1,jj ,jk ), & 440 & zbdo(ji ,jj-1,jk ), zbdo(ji ,jj+1,jk ), & 441 & zbdo(ji ,jj ,ikm1), zbdo(ji ,jj ,jk+1) ) 442 443 ! positive part of the flux 444 zpos = MAX( 0., paa(ji-1,jj ,jk ) ) - MIN( 0., paa(ji ,jj ,jk ) ) & 445 & + MAX( 0., pbb(ji ,jj-1,jk ) ) - MIN( 0., pbb(ji ,jj ,jk ) ) & 446 & + MAX( 0., pcc(ji ,jj ,jk+1) ) - MIN( 0., pcc(ji ,jj ,jk ) ) 447 448 ! negative part of the flux 449 zneg = MAX( 0., paa(ji ,jj ,jk ) ) - MIN( 0., paa(ji-1,jj ,jk ) ) & 450 & + MAX( 0., pbb(ji ,jj ,jk ) ) - MIN( 0., pbb(ji ,jj-1,jk ) ) & 451 & + MAX( 0., pcc(ji ,jj ,jk ) ) - MIN( 0., pcc(ji ,jj ,jk+1) ) 452 453 ! up & down beta terms 454 zbt = e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) / p2dt 455 zbetup(ji,jj,jk) = ( zup - paft(ji,jj,jk) ) / ( zpos + zrtrn ) * zbt 456 zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zdo ) / ( zneg + zrtrn ) * zbt 457 END_2D 458 END DO 459 IF (nn_hls==1) CALL lbc_lnk( 'traadv_fct', zbetup, 'T', 1.0_wp , zbetdo, 'T', 1.0_wp, ld4only= .TRUE. ) ! lateral boundary cond. (unchanged sign) 460 461 ! 3. monotonic flux in the i & j direction (paa & pbb) 462 ! ---------------------------------------- 463 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 464 zau = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) ) 465 zbu = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) ) 466 zcu = ( 0.5 + SIGN( 0.5_wp , paa(ji,jj,jk) ) ) 467 paa(ji,jj,jk) = paa(ji,jj,jk) * ( zcu * zau + ( 1._wp - zcu) * zbu ) 468 469 zav = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) ) 470 zbv = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) ) 471 zcv = ( 0.5 + SIGN( 0.5_wp , pbb(ji,jj,jk) ) ) 472 pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( 1._wp - zcv) * zbv ) 473 474 ! monotonic flux in the k direction, i.e. pcc 475 ! ------------------------------------------- 476 za = MIN( 1., zbetdo(ji,jj,jk+1), zbetup(ji,jj,jk) ) 477 zb = MIN( 1., zbetup(ji,jj,jk+1), zbetdo(ji,jj,jk) ) 478 zc = ( 0.5 + SIGN( 0.5_wp , pcc(ji,jj,jk+1) ) ) 479 pcc(ji,jj,jk+1) = pcc(ji,jj,jk+1) * ( zc * za + ( 1._wp - zc) * zb ) 480 END_3D 481 ! 482 END SUBROUTINE nonosc 483 484 485 SUBROUTINE interp_4th_cpt_org( pt_in, pt_out ) 486 !!---------------------------------------------------------------------- 487 !! *** ROUTINE interp_4th_cpt_org *** 488 !! 489 !! ** Purpose : Compute the interpolation of tracer at w-point 490 !! 491 !! ** Method : 4th order compact interpolation 492 !!---------------------------------------------------------------------- 493 REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pt_in ! now tracer fields 494 REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT( out) :: pt_out ! now tracer field interpolated at w-pts 495 ! 496 INTEGER :: ji, jj, jk ! dummy loop integers 497 REAL(wp),DIMENSION(jpi,jpj,jpk) :: zwd, zwi, zws, zwrm, zwt 498 !!---------------------------------------------------------------------- 499 500 DO_3D( 1, 1, 1, 1, 3, jpkm1 ) !== build the three diagonal matrix ==! 501 zwd (ji,jj,jk) = 4._wp 502 zwi (ji,jj,jk) = 1._wp 503 zws (ji,jj,jk) = 1._wp 504 zwrm(ji,jj,jk) = 3._wp * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) ) 505 ! 506 IF( tmask(ji,jj,jk+1) == 0._wp) THEN ! Switch to second order centered at bottom 507 zwd (ji,jj,jk) = 1._wp 508 zwi (ji,jj,jk) = 0._wp 509 zws (ji,jj,jk) = 0._wp 510 zwrm(ji,jj,jk) = 0.5 * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) ) 511 ENDIF 512 END_3D 513 ! 514 jk = 2 ! Switch to second order centered at top 515 DO_2D( 1, 1, 1, 1 ) 516 zwd (ji,jj,jk) = 1._wp 517 zwi (ji,jj,jk) = 0._wp 518 zws (ji,jj,jk) = 0._wp 519 zwrm(ji,jj,jk) = 0.5 * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) ) 520 END_2D 521 ! 522 ! !== tridiagonal solve ==! 523 DO_2D( 1, 1, 1, 1 ) ! first recurrence 524 zwt(ji,jj,2) = zwd(ji,jj,2) 525 END_2D 526 DO_3D( 1, 1, 1, 1, 3, jpkm1 ) 527 zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1) 528 END_3D 529 ! 530 DO_2D( 1, 1, 1, 1 ) ! second recurrence: Zk = Yk - Ik / Tk-1 Zk-1 531 pt_out(ji,jj,2) = zwrm(ji,jj,2) 532 END_2D 533 DO_3D( 1, 1, 1, 1, 3, jpkm1 ) 534 pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1) 535 END_3D 536 537 DO_2D( 1, 1, 1, 1 ) ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk 538 pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) 539 END_2D 540 DO_3DS( 1, 1, 1, 1, jpk-2, 2, -1 ) 541 pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk) 542 END_3D 543 ! 544 END SUBROUTINE interp_4th_cpt_org 545 546 547 SUBROUTINE interp_4th_cpt( pt_in, pt_out ) 548 !!---------------------------------------------------------------------- 549 !! *** ROUTINE interp_4th_cpt *** 550 !! 551 !! ** Purpose : Compute the interpolation of tracer at w-point 552 !! 553 !! ** Method : 4th order compact interpolation 554 !!---------------------------------------------------------------------- 555 REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pt_in ! field at t-point 556 REAL(wp),DIMENSION(A2D(nn_hls) ,jpk), INTENT( out) :: pt_out ! field interpolated at w-point 557 ! 558 INTEGER :: ji, jj, jk ! dummy loop integers 559 INTEGER :: ikt, ikb ! local integers 560 REAL(wp),DIMENSION(A2D(nn_hls),jpk) :: zwd, zwi, zws, zwrm, zwt 561 !!---------------------------------------------------------------------- 562 ! 563 ! !== build the three diagonal matrix & the RHS ==! 564 ! 565 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 3, jpkm1 ) ! interior (from jk=3 to jpk-1) 566 zwd (ji,jj,jk) = 3._wp * wmask(ji,jj,jk) + 1._wp ! diagonal 567 zwi (ji,jj,jk) = wmask(ji,jj,jk) ! lower diagonal 568 zws (ji,jj,jk) = wmask(ji,jj,jk) ! upper diagonal 569 zwrm(ji,jj,jk) = 3._wp * wmask(ji,jj,jk) & ! RHS 570 & * ( pt_in(ji,jj,jk) + pt_in(ji,jj,jk-1) ) 571 END_3D 572 ! 573 !!gm 574 ! SELECT CASE( kbc ) !* boundary condition 575 ! CASE( np_NH ) ! Neumann homogeneous at top & bottom 576 ! CASE( np_CEN2 ) ! 2nd order centered at top & bottom 577 ! END SELECT 578 !!gm 579 ! 580 IF ( ln_isfcav ) THEN ! set level two values which may not be set in ISF case 581 zwd(:,:,2) = 1._wp ; zwi(:,:,2) = 0._wp ; zws(:,:,2) = 0._wp ; zwrm(:,:,2) = 0._wp 582 END IF 583 ! 584 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! 2nd order centered at top & bottom 585 ikt = mikt(ji,jj) + 1 ! w-point below the 1st wet point 586 ikb = MAX(mbkt(ji,jj), 2) ! - above the last wet point 587 ! 588 zwd (ji,jj,ikt) = 1._wp ! top 589 zwi (ji,jj,ikt) = 0._wp 590 zws (ji,jj,ikt) = 0._wp 591 zwrm(ji,jj,ikt) = 0.5_wp * ( pt_in(ji,jj,ikt-1) + pt_in(ji,jj,ikt) ) 592 ! 593 zwd (ji,jj,ikb) = 1._wp ! bottom 594 zwi (ji,jj,ikb) = 0._wp 595 zws (ji,jj,ikb) = 0._wp 596 zwrm(ji,jj,ikb) = 0.5_wp * ( pt_in(ji,jj,ikb-1) + pt_in(ji,jj,ikb) ) 597 END_2D 598 ! 599 ! !== tridiagonal solver ==! 600 ! 601 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) !* 1st recurrence: Tk = Dk - Ik Sk-1 / Tk-1 602 zwt(ji,jj,2) = zwd(ji,jj,2) 603 END_2D 604 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 3, jpkm1 ) 605 zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1) 606 END_3D 607 ! 608 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) !* 2nd recurrence: Zk = Yk - Ik / Tk-1 Zk-1 609 pt_out(ji,jj,2) = zwrm(ji,jj,2) 610 END_2D 611 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 3, jpkm1 ) 612 pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1) 613 END_3D 614 615 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) !* 3d recurrence: Xk = (Zk - Sk Xk+1 ) / Tk 616 pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) 617 END_2D 618 DO_3DS( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpk-2, 2, -1 ) 619 pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk) 620 END_3D 621 ! 622 END SUBROUTINE interp_4th_cpt 623 624 625 SUBROUTINE tridia_solver( pD, pU, pL, pRHS, pt_out , klev ) 626 !!---------------------------------------------------------------------- 627 !! *** ROUTINE tridia_solver *** 628 !! 629 !! ** Purpose : solve a symmetric 3diagonal system 630 !! 631 !! ** Method : solve M.t_out = RHS(t) where M is a tri diagonal matrix ( jpk*jpk ) 632 !! 633 !! ( D_1 U_1 0 0 0 )( t_1 ) ( RHS_1 ) 634 !! ( L_2 D_2 U_2 0 0 )( t_2 ) ( RHS_2 ) 635 !! ( 0 L_3 D_3 U_3 0 )( t_3 ) = ( RHS_3 ) 636 !! ( ... )( ... ) ( ... ) 637 !! ( 0 0 0 L_k D_k )( t_k ) ( RHS_k ) 638 !! 639 !! M is decomposed in the product of an upper and lower triangular matrix. 640 !! The tri-diagonals matrix is given as input 3D arrays: pD, pU, pL 641 !! (i.e. the Diagonal, the Upper diagonal, and the Lower diagonal). 642 !! The solution is pta. 643 !! The 3d array zwt is used as a work space array. 644 !!---------------------------------------------------------------------- 645 REAL(wp),DIMENSION(A2D(nn_hls),jpk), INTENT(in ) :: pD, pU, PL ! 3-diagonal matrix 646 REAL(wp),DIMENSION(A2D(nn_hls),jpk), INTENT(in ) :: pRHS ! Right-Hand-Side 647 REAL(wp),DIMENSION(A2D(nn_hls),jpk), INTENT( out) :: pt_out !!gm field at level=F(klev) 648 INTEGER , INTENT(in ) :: klev ! =1 pt_out at w-level 649 ! ! =0 pt at t-level 650 INTEGER :: ji, jj, jk ! dummy loop integers 651 INTEGER :: kstart ! local indices 652 REAL(wp),DIMENSION(A2D(nn_hls),jpk) :: zwt ! 3D work array 653 !!---------------------------------------------------------------------- 654 ! 655 kstart = 1 + klev 656 ! 657 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) !* 1st recurrence: Tk = Dk - Ik Sk-1 / Tk-1 658 zwt(ji,jj,kstart) = pD(ji,jj,kstart) 659 END_2D 660 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, kstart+1, jpkm1 ) 661 zwt(ji,jj,jk) = pD(ji,jj,jk) - pL(ji,jj,jk) * pU(ji,jj,jk-1) /zwt(ji,jj,jk-1) 662 END_3D 663 ! 664 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) !* 2nd recurrence: Zk = Yk - Ik / Tk-1 Zk-1 665 pt_out(ji,jj,kstart) = pRHS(ji,jj,kstart) 666 END_2D 667 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, kstart+1, jpkm1 ) 668 pt_out(ji,jj,jk) = pRHS(ji,jj,jk) - pL(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1) 669 END_3D 670 671 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) !* 3d recurrence: Xk = (Zk - Sk Xk+1 ) / Tk 672 pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) 673 END_2D 674 DO_3DS( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpk-2, kstart, -1 ) 675 pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - pU(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk) 676 END_3D 677 ! 678 END SUBROUTINE tridia_solver 679 680 #if defined key_loop_fusion 681 #define tracer_flux_i(out,zfp,zfm,ji,jj,jk) \ 682 zfp = pU(ji,jj,jk) + ABS( pU(ji,jj,jk) ) ; \ 683 zfm = pU(ji,jj,jk) - ABS( pU(ji,jj,jk) ) ; \ 684 out = 0.5 * ( zfp * pt(ji,jj,jk,jn,Kbb) + zfm * pt(ji+1,jj,jk,jn,Kbb) ) 685 686 #define tracer_flux_j(out,zfp,zfm,ji,jj,jk) \ 687 zfp = pV(ji,jj,jk) + ABS( pV(ji,jj,jk) ) ; \ 688 zfm = pV(ji,jj,jk) - ABS( pV(ji,jj,jk) ) ; \ 689 out = 0.5 * ( zfp * pt(ji,jj,jk,jn,Kbb) + zfm * pt(ji,jj+1,jk,jn,Kbb) ) 690 691 SUBROUTINE tra_adv_fct_lf( kt, kit000, cdtype, p2dt, pU, pV, pW, & 692 & Kbb, Kmm, pt, kjpt, Krhs, kn_fct_h, kn_fct_v ) 693 !!---------------------------------------------------------------------- 694 !! *** ROUTINE tra_adv_fct *** 695 !! 696 !! ** Purpose : Compute the now trend due to total advection of tracers 697 !! and add it to the general trend of tracer equations 698 !! 699 !! ** Method : - 2nd or 4th FCT scheme on the horizontal direction 700 !! (choice through the value of kn_fct) 701 !! - on the vertical the 4th order is a compact scheme 702 !! - corrected flux (monotonic correction) 703 !! 704 !! ** Action : - update pt(:,:,:,:,Krhs) with the now advective tracer trends 705 !! - send trends to trdtra module for further diagnostics (l_trdtra=T) 706 !! - poleward advective heat and salt transport (ln_diaptr=T) 707 !!---------------------------------------------------------------------- 708 INTEGER , INTENT(in ) :: kt ! ocean time-step index 709 INTEGER , INTENT(in ) :: Kbb, Kmm, Krhs ! ocean time level indices 710 INTEGER , INTENT(in ) :: kit000 ! first time step index 711 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 712 INTEGER , INTENT(in ) :: kjpt ! number of tracers 713 INTEGER , INTENT(in ) :: kn_fct_h ! order of the FCT scheme (=2 or 4) 714 INTEGER , INTENT(in ) :: kn_fct_v ! order of the FCT scheme (=2 or 4) 715 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 716 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pU, pV, pW ! 3 ocean volume flux components 717 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! tracers and RHS of tracer equation 718 ! 719 INTEGER :: ji, jj, jk, jn ! dummy loop indices 720 REAL(wp) :: ztra ! local scalar 721 REAL(wp) :: zwx_im1, zfp_ui, zfp_ui_m1, zfp_vj, zfp_vj_m1, zfp_wk, zC2t_u, zC4t_u ! - - 722 REAL(wp) :: zwy_jm1, zfm_ui, zfm_ui_m1, zfm_vj, zfm_vj_m1, zfm_wk, zC2t_v, zC4t_v ! - - 723 REAL(wp) :: ztu, ztv, ztu_im1, ztu_ip1, ztv_jm1, ztv_jp1 724 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwi, zwx_3d, zwy_3d, zwz, ztw, zltu_3d, zltv_3d 725 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdx, ztrdy, ztrdz, zptry 726 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zwinf, zwdia, zwsup 727 LOGICAL :: ll_zAimp ! flag to apply adaptive implicit vertical advection 728 !!---------------------------------------------------------------------- 729 ! 730 IF( kt == kit000 ) THEN 731 IF(lwp) WRITE(numout,*) 732 IF(lwp) WRITE(numout,*) 'tra_adv_fct_lf : FCT advection scheme on ', cdtype 733 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 734 ENDIF 735 !! -- init to 0 736 zwx_3d(:,:,:) = 0._wp 737 zwy_3d(:,:,:) = 0._wp 738 zwz(:,:,:) = 0._wp 739 zwi(:,:,:) = 0._wp 740 ! 741 l_trd = .FALSE. ! set local switches 742 l_hst = .FALSE. 743 l_ptr = .FALSE. 744 ll_zAimp = .FALSE. 745 IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype =='TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 746 IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) ) l_ptr = .TRUE. 747 IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. & 748 & iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) ) l_hst = .TRUE. 749 ! 750 IF( l_trd .OR. l_hst ) THEN 751 ALLOCATE( ztrdx(jpi,jpj,jpk), ztrdy(jpi,jpj,jpk), ztrdz(jpi,jpj,jpk) ) 752 ztrdx(:,:,:) = 0._wp ; ztrdy(:,:,:) = 0._wp ; ztrdz(:,:,:) = 0._wp 753 ENDIF 754 ! 755 IF( l_ptr ) THEN 756 ALLOCATE( zptry(jpi,jpj,jpk) ) 757 zptry(:,:,:) = 0._wp 758 ENDIF 759 ! 760 ! If adaptive vertical advection, check if it is needed on this PE at this time 761 IF( ln_zad_Aimp ) THEN 762 IF( MAXVAL( ABS( wi(:,:,:) ) ) > 0._wp ) ll_zAimp = .TRUE. 763 END IF 764 ! If active adaptive vertical advection, build tridiagonal matrix 765 IF( ll_zAimp ) THEN 766 ALLOCATE(zwdia(jpi,jpj,jpk), zwinf(jpi,jpj,jpk),zwsup(jpi,jpj,jpk)) 767 DO_3D( 1, 1, 1, 1, 1, jpkm1 ) 768 zwdia(ji,jj,jk) = 1._wp + p2dt * ( MAX( wi(ji,jj,jk) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) ) & 769 & / e3t(ji,jj,jk,Krhs) 770 zwinf(ji,jj,jk) = p2dt * MIN( wi(ji,jj,jk ) , 0._wp ) / e3t(ji,jj,jk,Krhs) 771 zwsup(ji,jj,jk) = -p2dt * MAX( wi(ji,jj,jk+1) , 0._wp ) / e3t(ji,jj,jk,Krhs) 772 END_3D 773 END IF 774 ! 775 DO jn = 1, kjpt !== loop over the tracers ==! 776 ! 777 ! !== upstream advection with initial mass fluxes & intermediate update ==! 164 778 ! !* upstream tracer flux in the k direction *! 165 779 DO_3D( 1, 1, 1, 1, 2, jpkm1 ) ! Interior value ( multiplied by wmask) … … 180 794 ENDIF 181 795 ! 182 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) !* trend and after field with monotonic scheme 183 ! ! total intermediate advective trends 184 ztra = - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 185 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & 186 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) * r1_e1e2t(ji,jj) 187 ! ! update and guess with monotonic sheme 188 pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztra & 189 & / e3t(ji,jj,jk,Kmm ) * tmask(ji,jj,jk) 190 zwi(ji,jj,jk) = ( e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) + p2dt * ztra ) & 191 & / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk) 192 END_3D 796 ! !* upstream tracer flux in the i and j direction 797 DO jk = 1, jpkm1 798 DO jj = 1, jpj-1 799 tracer_flux_i(zwx_3d(1,jj,jk),zfp_ui,zfm_ui,1,jj,jk) 800 tracer_flux_j(zwy_3d(1,jj,jk),zfp_vj,zfm_vj,1,jj,jk) 801 END DO 802 DO ji = 1, jpi-1 803 tracer_flux_i(zwx_3d(ji,1,jk),zfp_ui,zfm_ui,ji,1,jk) 804 tracer_flux_j(zwy_3d(ji,1,jk),zfp_vj,zfm_vj,ji,1,jk) 805 END DO 806 DO_2D( 1, 1, 1, 1 ) 807 tracer_flux_i(zwx_3d(ji,jj,jk),zfp_ui,zfm_ui,ji,jj,jk) 808 tracer_flux_i(zwx_im1,zfp_ui_m1,zfm_ui_m1,ji-1,jj,jk) 809 tracer_flux_j(zwy_3d(ji,jj,jk),zfp_vj,zfm_vj,ji,jj,jk) 810 tracer_flux_j(zwy_jm1,zfp_vj_m1,zfm_vj_m1,ji,jj-1,jk) 811 ztra = - ( zwx_3d(ji,jj,jk) - zwx_im1 + zwy_3d(ji,jj,jk) - zwy_jm1 + zwz(ji,jj,jk) - zwz(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) 812 ! ! update and guess with monotonic sheme 813 pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztra & 814 & / e3t(ji,jj,jk,Kmm ) * tmask(ji,jj,jk) 815 zwi(ji,jj,jk) = ( e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) + p2dt * ztra ) & 816 & / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk) 817 END_2D 818 END DO 193 819 194 820 IF ( ll_zAimp ) THEN … … 196 822 ! 197 823 ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp ; 198 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! Interior value ( multiplied by wmask)824 DO_3D( 1, 1, 1, 1, 2, jpkm1 ) ! Interior value ( multiplied by wmask) 199 825 zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) ) 200 826 zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) ) … … 210 836 ! 211 837 IF( l_trd .OR. l_hst ) THEN ! trend diagnostics (contribution of upstream fluxes) 212 ztrdx(:,:,:) = zwx (:,:,:) ; ztrdy(:,:,:) = zwy(:,:,:) ; ztrdz(:,:,:) = zwz(:,:,:)838 ztrdx(:,:,:) = zwx_3d(:,:,:) ; ztrdy(:,:,:) = zwy_3d(:,:,:) ; ztrdz(:,:,:) = zwz(:,:,:) 213 839 END IF 214 840 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) 215 IF( l_ptr ) zptry(:,:,:) = zwy (:,:,:)841 IF( l_ptr ) zptry(:,:,:) = zwy_3d(:,:,:) 216 842 ! 217 843 ! !== anti-diffusive flux : high order minus low order ==! … … 220 846 ! 221 847 CASE( 2 ) !- 2nd order centered 222 DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 )223 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)224 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)848 DO_3D( 2, 1, 2, 1, 1, jpkm1 ) 849 zwx_3d(ji,jj,jk) = 0.5_wp * pU(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj,jk,jn,Kmm) ) - zwx_3d(ji,jj,jk) 850 zwy_3d(ji,jj,jk) = 0.5_wp * pV(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj+1,jk,jn,Kmm) ) - zwy_3d(ji,jj,jk) 225 851 END_3D 226 852 ! 227 853 CASE( 4 ) !- 4th order centered 228 zltu(:,:,jpk) = 0._wp ! Bottom value : flux set to zero 229 zltv(:,:,jpk) = 0._wp 230 DO jk = 1, jpkm1 ! Laplacian 231 DO_2D( 1, 0, 1, 0 ) ! 1st derivative (gradient) 232 ztu(ji,jj,jk) = ( pt(ji+1,jj ,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk) 233 ztv(ji,jj,jk) = ( pt(ji ,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk) 234 END_2D 235 DO_2D( 0, 0, 0, 0 ) ! 2nd derivative * 1/ 6 236 zltu(ji,jj,jk) = ( ztu(ji,jj,jk) + ztu(ji-1,jj,jk) ) * r1_6 237 zltv(ji,jj,jk) = ( ztv(ji,jj,jk) + ztv(ji,jj-1,jk) ) * r1_6 854 zltu_3d(:,:,jpk) = 0._wp ! Bottom value : flux set to zero 855 zltv_3d(:,:,jpk) = 0._wp 856 ! ! Laplacian 857 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! 2nd derivative * 1/ 6 858 ! ! 1st derivative (gradient) 859 ztu = ( pt(ji+1,jj,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk) 860 ztu_im1 = ( pt(ji,jj,jk,jn,Kmm) - pt(ji-1,jj,jk,jn,Kmm) ) * umask(ji-1,jj,jk) 861 ztv = ( pt(ji,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk) 862 ztv_jm1 = ( pt(ji,jj,jk,jn,Kmm) - pt(ji,jj-1,jk,jn,Kmm) ) * vmask(ji,jj-1,jk) 863 ! ! 2nd derivative * 1/ 6 864 zltu_3d(ji,jj,jk) = ( ztu + ztu_im1 ) * r1_6 865 zltv_3d(ji,jj,jk) = ( ztv + ztv_jm1 ) * r1_6 238 866 END_2D 239 867 END DO 240 CALL lbc_lnk( 'traadv_fct', zltu, 'T', 1.0_wp , zltv, 'T', 1.0_wp ) ! Lateral boundary cond. (unchanged sgn) 241 ! 242 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 868 ! NOTE [ comm_cleanup ] : need to change sign to ensure halo 1 - halo 2 compatibility 869 CALL lbc_lnk( 'traadv_fct', zltu_3d, 'T', -1.0_wp , zltv_3d, 'T', -1.0_wp ) ! Lateral boundary cond. (unchanged sgn) 870 ! 871 DO_3D( 2, 1, 2, 1, 1, jpkm1 ) 243 872 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 244 873 zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji ,jj+1,jk,jn,Kmm) 245 874 ! ! C4 minus upstream advective fluxes 246 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) 247 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) 248 END_3D 249 IF (nn_hls.EQ.2) CALL lbc_lnk( 'traadv_fct', zwx, 'U', -1.0_wp, zwy, 'V', -1.0_wp ) ! Lateral boundary cond. (unchanged sgn) 875 ! round brackets added to fix the order of floating point operations 876 ! needed to ensure halo 1 - halo 2 compatibility 877 zwx_3d(ji,jj,jk) = 0.5_wp * pU(ji,jj,jk) * ( zC2t_u + ( zltu_3d(ji,jj,jk) - zltu_3d(ji+1,jj,jk) & 878 & ) & ! bracket for halo 1 - halo 2 compatibility 879 & ) - zwx_3d(ji,jj,jk) 880 zwy_3d(ji,jj,jk) = 0.5_wp * pV(ji,jj,jk) * ( zC2t_v + ( zltv_3d(ji,jj,jk) - zltv_3d(ji,jj+1,jk) & 881 & ) & ! bracket for halo 1 - halo 2 compatibility 882 & ) - zwy_3d(ji,jj,jk) 883 END_3D 250 884 ! 251 885 CASE( 41 ) !- 4th order centered ==>> !!gm coding attempt need to be tested 252 ztu(:,:,jpk) = 0._wp ! Bottom value : flux set to zero253 ztv(:,:,jpk) = 0._wp254 DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 ) ! 1st derivative (gradient)255 ztu(ji,jj,jk) = ( pt(ji+1,jj ,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk)256 ztv(ji,jj,jk) = ( pt(ji ,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk)257 END_3D258 IF (nn_hls.EQ.1) CALL lbc_lnk( 'traadv_fct', ztu, 'U', -1.0_wp , ztv, 'V', -1.0_wp ) ! Lateral boundary cond. (unchanged sgn)259 !260 886 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! Horizontal advective fluxes 887 ztu_im1 = ( pt(ji ,jj ,jk,jn,Kmm) - pt(ji-1,jj,jk,jn,Kmm) ) * umask(ji-1,jj,jk) 888 ztu_ip1 = ( pt(ji+2,jj ,jk,jn,Kmm) - pt(ji+1,jj,jk,jn,Kmm) ) * umask(ji+1,jj,jk) 889 890 ztv_jm1 = ( pt(ji,jj ,jk,jn,Kmm) - pt(ji,jj-1,jk,jn,Kmm) ) * vmask(ji,jj-1,jk) 891 ztv_jp1 = ( pt(ji,jj+2,jk,jn,Kmm) - pt(ji,jj+1,jk,jn,Kmm) ) * vmask(ji,jj+1,jk) 261 892 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) 262 893 zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji ,jj+1,jk,jn,Kmm) 263 894 ! ! C4 interpolation of T at u- & v-points (x2) 264 zC4t_u = zC2t_u + r1_6 * ( ztu (ji-1,jj ,jk) - ztu(ji+1,jj ,jk))265 zC4t_v = zC2t_v + r1_6 * ( ztv (ji ,jj-1,jk) - ztv(ji ,jj+1,jk))895 zC4t_u = zC2t_u + r1_6 * ( ztu_im1 - ztu_ip1 ) 896 zC4t_v = zC2t_v + r1_6 * ( ztv_jm1 - ztv_jp1 ) 266 897 ! ! C4 minus upstream advective fluxes 267 zwx (ji,jj,jk) = 0.5_wp * pU(ji,jj,jk) * zC4t_u - zwx(ji,jj,jk)268 zwy (ji,jj,jk) = 0.5_wp * pV(ji,jj,jk) * zC4t_v - zwy(ji,jj,jk)269 END_3D 270 IF (nn_hls.EQ.2) CALL lbc_lnk( 'traadv_fct', zwx, 'U', -1.0_wp , zwy, 'V', -1.0_wp ) ! Lateral boundary cond. (unchanged sgn)898 zwx_3d(ji,jj,jk) = 0.5_wp * pU(ji,jj,jk) * zC4t_u - zwx_3d(ji,jj,jk) 899 zwy_3d(ji,jj,jk) = 0.5_wp * pV(ji,jj,jk) * zC4t_v - zwy_3d(ji,jj,jk) 900 END_3D 901 CALL lbc_lnk( 'traadv_fct', zwx_3d, 'U', -1.0_wp , zwy_3d, 'V', -1.0_wp ) ! Lateral boundary cond. (unchanged sgn) 271 902 ! 272 903 END SELECT … … 275 906 ! 276 907 CASE( 2 ) !- 2nd order centered 277 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )908 DO_3D( 1, 1, 1, 1, 2, jpkm1 ) 278 909 zwz(ji,jj,jk) = ( pW(ji,jj,jk) * 0.5_wp * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj,jk-1,jn,Kmm) ) & 279 910 & - zwz(ji,jj,jk) ) * wmask(ji,jj,jk) … … 282 913 CASE( 4 ) !- 4th order COMPACT 283 914 CALL interp_4th_cpt( pt(:,:,:,jn,Kmm) , ztw ) ! zwt = COMPACT interpolation of T at w-point 284 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )915 DO_3D( 1, 1, 1, 1, 2, jpkm1 ) 285 916 zwz(ji,jj,jk) = ( pW(ji,jj,jk) * ztw(ji,jj,jk) - zwz(ji,jj,jk) ) * wmask(ji,jj,jk) 286 917 END_3D … … 291 922 ENDIF 292 923 ! 293 IF (nn_hls.EQ.1) THEN 294 CALL lbc_lnk( 'traadv_fct', zwi, 'T', 1.0_wp, zwx, 'U', -1.0_wp , zwy, 'V', -1.0_wp, zwz, 'T', 1.0_wp ) 295 ELSE 296 CALL lbc_lnk( 'traadv_fct', zwi, 'T', 1.0_wp) 297 END IF 924 CALL lbc_lnk( 'traadv_fct', zwi, 'T', 1.0_wp) 298 925 ! 299 926 IF ( ll_zAimp ) THEN 300 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) !* trend and after field with monotonic scheme927 DO_3D( 1, 1, 1, 1, 1, jpkm1 ) !* trend and after field with monotonic scheme 301 928 ! ! total intermediate advective trends 302 ztra = - ( zwx (ji,jj,jk) - zwx(ji-1,jj ,jk ) &303 & + zwy (ji,jj,jk) - zwy(ji ,jj-1,jk ) &929 ztra = - ( zwx_3d(ji,jj,jk) - zwx_3d(ji-1,jj ,jk ) & 930 & + zwy_3d(ji,jj,jk) - zwy_3d(ji ,jj-1,jk ) & 304 931 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) * r1_e1e2t(ji,jj) 305 932 ztw(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk) … … 308 935 CALL tridia_solver( zwdia, zwsup, zwinf, ztw, ztw , 0 ) 309 936 ! 310 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! Interior value ( multiplied by wmask)937 DO_3D( 1, 1, 1, 1, 2, jpkm1 ) ! Interior value ( multiplied by wmask) 311 938 zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) ) 312 939 zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) ) … … 317 944 ! !== monotonicity algorithm ==! 318 945 ! 319 CALL nonosc( Kmm, pt(:,:,:,jn,Kbb), zwx , zwy, zwz, zwi, p2dt )946 CALL nonosc( Kmm, pt(:,:,:,jn,Kbb), zwx_3d, zwy_3d, zwz, zwi, p2dt ) 320 947 ! 321 948 ! !== final trend with corrected fluxes ==! 322 949 ! 323 950 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 324 ztra = - ( zwx (ji,jj,jk) - zwx(ji-1,jj ,jk ) &325 & + zwy (ji,jj,jk) - zwy(ji ,jj-1,jk ) &951 ztra = - ( zwx_3d(ji,jj,jk) - zwx_3d(ji-1,jj ,jk ) & 952 & + zwy_3d(ji,jj,jk) - zwy_3d(ji ,jj-1,jk ) & 326 953 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) * r1_e1e2t(ji,jj) 327 954 pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztra / e3t(ji,jj,jk,Kmm) … … 343 970 END_3D 344 971 END IF 972 ! NOT TESTED - NEED l_trd OR l_hst TRUE 345 973 IF( l_trd .OR. l_hst ) THEN ! trend diagnostics // heat/salt transport 346 ztrdx(:,:,:) = ztrdx(:,:,:) + zwx (:,:,:) ! <<< add anti-diffusive fluxes347 ztrdy(:,:,:) = ztrdy(:,:,:) + zwy (:,:,:) ! to upstream fluxes974 ztrdx(:,:,:) = ztrdx(:,:,:) + zwx_3d(:,:,:) ! <<< add anti-diffusive fluxes 975 ztrdy(:,:,:) = ztrdy(:,:,:) + zwy_3d(:,:,:) ! to upstream fluxes 348 976 ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:) ! 349 977 ! … … 357 985 ! 358 986 ENDIF 987 ! NOT TESTED - NEED l_ptr TRUE 359 988 IF( l_ptr ) THEN ! "Poleward" transports 360 zptry(:,:,:) = zptry(:,:,:) + zwy (:,:,:) ! <<< add anti-diffusive fluxes989 zptry(:,:,:) = zptry(:,:,:) + zwy_3d(:,:,:) ! <<< add anti-diffusive fluxes 361 990 CALL dia_ptr_hst( jn, 'adv', zptry(:,:,:) ) 362 991 ENDIF … … 374 1003 ENDIF 375 1004 ! 376 END SUBROUTINE tra_adv_fct 377 378 379 SUBROUTINE nonosc( Kmm, pbef, paa, pbb, pcc, paft, p2dt ) 380 !!--------------------------------------------------------------------- 381 !! *** ROUTINE nonosc *** 382 !! 383 !! ** Purpose : compute monotonic tracer fluxes from the upstream 384 !! scheme and the before field by a nonoscillatory algorithm 385 !! 386 !! ** Method : ... ??? 387 !! warning : pbef and paft must be masked, but the boundaries 388 !! conditions on the fluxes are not necessary zalezak (1979) 389 !! drange (1995) multi-dimensional forward-in-time and upstream- 390 !! in-space based differencing for fluid 391 !!---------------------------------------------------------------------- 392 INTEGER , INTENT(in ) :: Kmm ! time level index 393 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 394 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pbef ! before field 395 REAL(wp), DIMENSION(A2D(nn_hls) ,jpk), INTENT(in ) :: paft ! after field 396 REAL(wp), DIMENSION(A2D(nn_hls) ,jpk), INTENT(inout) :: paa, pbb, pcc ! monotonic fluxes in the 3 directions 397 ! 398 INTEGER :: ji, jj, jk ! dummy loop indices 399 INTEGER :: ikm1 ! local integer 400 REAL(dp) :: zpos, zneg, zbt, za, zb, zc, zbig, zrtrn ! local scalars 401 REAL(dp) :: zau, zbu, zcu, zav, zbv, zcv, zup, zdo ! - - 402 REAL(dp), DIMENSION(A2D(nn_hls),jpk) :: zbetup, zbetdo, zbup, zbdo 403 !!---------------------------------------------------------------------- 404 ! 405 zbig = 1.e+40_dp 406 zrtrn = 1.e-15_dp 407 zbetup(:,:,:) = 0._dp ; zbetdo(:,:,:) = 0._dp 408 409 ! Search local extrema 410 ! -------------------- 411 ! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land 412 DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk ) 413 zbup(ji,jj,jk) = MAX( pbef(ji,jj,jk) * tmask(ji,jj,jk) - zbig * ( 1._wp - tmask(ji,jj,jk) ), & 414 & paft(ji,jj,jk) * tmask(ji,jj,jk) - zbig * ( 1._wp - tmask(ji,jj,jk) ) ) 415 zbdo(ji,jj,jk) = MIN( pbef(ji,jj,jk) * tmask(ji,jj,jk) + zbig * ( 1._wp - tmask(ji,jj,jk) ), & 416 & paft(ji,jj,jk) * tmask(ji,jj,jk) + zbig * ( 1._wp - tmask(ji,jj,jk) ) ) 417 END_3D 418 419 DO jk = 1, jpkm1 420 ikm1 = MAX(jk-1,1) 421 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 422 423 ! search maximum in neighbourhood 424 zup = MAX( zbup(ji ,jj ,jk ), & 425 & zbup(ji-1,jj ,jk ), zbup(ji+1,jj ,jk ), & 426 & zbup(ji ,jj-1,jk ), zbup(ji ,jj+1,jk ), & 427 & zbup(ji ,jj ,ikm1), zbup(ji ,jj ,jk+1) ) 428 429 ! search minimum in neighbourhood 430 zdo = MIN( zbdo(ji ,jj ,jk ), & 431 & zbdo(ji-1,jj ,jk ), zbdo(ji+1,jj ,jk ), & 432 & zbdo(ji ,jj-1,jk ), zbdo(ji ,jj+1,jk ), & 433 & zbdo(ji ,jj ,ikm1), zbdo(ji ,jj ,jk+1) ) 434 435 ! positive part of the flux 436 zpos = MAX( 0., paa(ji-1,jj ,jk ) ) - MIN( 0., paa(ji ,jj ,jk ) ) & 437 & + MAX( 0., pbb(ji ,jj-1,jk ) ) - MIN( 0., pbb(ji ,jj ,jk ) ) & 438 & + MAX( 0., pcc(ji ,jj ,jk+1) ) - MIN( 0., pcc(ji ,jj ,jk ) ) 439 440 ! negative part of the flux 441 zneg = MAX( 0., paa(ji ,jj ,jk ) ) - MIN( 0., paa(ji-1,jj ,jk ) ) & 442 & + MAX( 0., pbb(ji ,jj ,jk ) ) - MIN( 0., pbb(ji ,jj-1,jk ) ) & 443 & + MAX( 0., pcc(ji ,jj ,jk ) ) - MIN( 0., pcc(ji ,jj ,jk+1) ) 444 445 ! up & down beta terms 446 zbt = e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) / p2dt 447 zbetup(ji,jj,jk) = ( zup - paft(ji,jj,jk) ) / ( zpos + zrtrn ) * zbt 448 zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zdo ) / ( zneg + zrtrn ) * zbt 449 END_2D 450 END DO 451 IF (nn_hls.EQ.1) CALL lbc_lnk( 'traadv_fct', zbetup, 'T', 1.0_wp , zbetdo, 'T', 1.0_wp ) ! lateral boundary cond. (unchanged sign) 452 453 ! 3. monotonic flux in the i & j direction (paa & pbb) 454 ! ---------------------------------------- 455 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 456 zau = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) ) 457 zbu = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) ) 458 zcu = ( 0.5 + SIGN( 0.5_wp , paa(ji,jj,jk) ) ) 459 paa(ji,jj,jk) = paa(ji,jj,jk) * ( zcu * zau + ( 1._wp - zcu) * zbu ) 460 461 zav = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) ) 462 zbv = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) ) 463 zcv = ( 0.5 + SIGN( 0.5_wp , pbb(ji,jj,jk) ) ) 464 pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( 1._wp - zcv) * zbv ) 465 466 ! monotonic flux in the k direction, i.e. pcc 467 ! ------------------------------------------- 468 za = MIN( 1., zbetdo(ji,jj,jk+1), zbetup(ji,jj,jk) ) 469 zb = MIN( 1., zbetup(ji,jj,jk+1), zbetdo(ji,jj,jk) ) 470 zc = ( 0.5 + SIGN( 0.5_wp , pcc(ji,jj,jk+1) ) ) 471 pcc(ji,jj,jk+1) = pcc(ji,jj,jk+1) * ( zc * za + ( 1._wp - zc) * zb ) 472 END_3D 473 ! 474 END SUBROUTINE nonosc 475 476 477 SUBROUTINE interp_4th_cpt_org( pt_in, pt_out ) 478 !!---------------------------------------------------------------------- 479 !! *** ROUTINE interp_4th_cpt_org *** 480 !! 481 !! ** Purpose : Compute the interpolation of tracer at w-point 482 !! 483 !! ** Method : 4th order compact interpolation 484 !!---------------------------------------------------------------------- 485 REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pt_in ! now tracer fields 486 REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT( out) :: pt_out ! now tracer field interpolated at w-pts 487 ! 488 INTEGER :: ji, jj, jk ! dummy loop integers 489 REAL(wp),DIMENSION(jpi,jpj,jpk) :: zwd, zwi, zws, zwrm, zwt 490 !!---------------------------------------------------------------------- 491 492 DO_3D( 1, 1, 1, 1, 3, jpkm1 ) !== build the three diagonal matrix ==! 493 zwd (ji,jj,jk) = 4._wp 494 zwi (ji,jj,jk) = 1._wp 495 zws (ji,jj,jk) = 1._wp 496 zwrm(ji,jj,jk) = 3._wp * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) ) 497 ! 498 IF( tmask(ji,jj,jk+1) == 0._wp) THEN ! Switch to second order centered at bottom 499 zwd (ji,jj,jk) = 1._wp 500 zwi (ji,jj,jk) = 0._wp 501 zws (ji,jj,jk) = 0._wp 502 zwrm(ji,jj,jk) = 0.5 * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) ) 503 ENDIF 504 END_3D 505 ! 506 jk = 2 ! Switch to second order centered at top 507 DO_2D( 1, 1, 1, 1 ) 508 zwd (ji,jj,jk) = 1._wp 509 zwi (ji,jj,jk) = 0._wp 510 zws (ji,jj,jk) = 0._wp 511 zwrm(ji,jj,jk) = 0.5 * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) ) 512 END_2D 513 ! 514 ! !== tridiagonal solve ==! 515 DO_2D( 1, 1, 1, 1 ) ! first recurrence 516 zwt(ji,jj,2) = zwd(ji,jj,2) 517 END_2D 518 DO_3D( 1, 1, 1, 1, 3, jpkm1 ) 519 zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1) 520 END_3D 521 ! 522 DO_2D( 1, 1, 1, 1 ) ! second recurrence: Zk = Yk - Ik / Tk-1 Zk-1 523 pt_out(ji,jj,2) = zwrm(ji,jj,2) 524 END_2D 525 DO_3D( 1, 1, 1, 1, 3, jpkm1 ) 526 pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1) 527 END_3D 528 529 DO_2D( 1, 1, 1, 1 ) ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk 530 pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) 531 END_2D 532 DO_3DS( 1, 1, 1, 1, jpk-2, 2, -1 ) 533 pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk) 534 END_3D 535 ! 536 END SUBROUTINE interp_4th_cpt_org 537 538 539 SUBROUTINE interp_4th_cpt( pt_in, pt_out ) 540 !!---------------------------------------------------------------------- 541 !! *** ROUTINE interp_4th_cpt *** 542 !! 543 !! ** Purpose : Compute the interpolation of tracer at w-point 544 !! 545 !! ** Method : 4th order compact interpolation 546 !!---------------------------------------------------------------------- 547 REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pt_in ! field at t-point 548 REAL(wp),DIMENSION(A2D(nn_hls) ,jpk), INTENT( out) :: pt_out ! field interpolated at w-point 549 ! 550 INTEGER :: ji, jj, jk ! dummy loop integers 551 INTEGER :: ikt, ikb ! local integers 552 REAL(wp),DIMENSION(A2D(nn_hls),jpk) :: zwd, zwi, zws, zwrm, zwt 553 !!---------------------------------------------------------------------- 554 ! 555 ! !== build the three diagonal matrix & the RHS ==! 556 ! 557 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 3, jpkm1 ) ! interior (from jk=3 to jpk-1) 558 zwd (ji,jj,jk) = 3._wp * wmask(ji,jj,jk) + 1._wp ! diagonal 559 zwi (ji,jj,jk) = wmask(ji,jj,jk) ! lower diagonal 560 zws (ji,jj,jk) = wmask(ji,jj,jk) ! upper diagonal 561 zwrm(ji,jj,jk) = 3._wp * wmask(ji,jj,jk) & ! RHS 562 & * ( pt_in(ji,jj,jk) + pt_in(ji,jj,jk-1) ) 563 END_3D 564 ! 565 !!gm 566 ! SELECT CASE( kbc ) !* boundary condition 567 ! CASE( np_NH ) ! Neumann homogeneous at top & bottom 568 ! CASE( np_CEN2 ) ! 2nd order centered at top & bottom 569 ! END SELECT 570 !!gm 571 ! 572 IF ( ln_isfcav ) THEN ! set level two values which may not be set in ISF case 573 zwd(:,:,2) = 1._wp ; zwi(:,:,2) = 0._wp ; zws(:,:,2) = 0._wp ; zwrm(:,:,2) = 0._wp 574 END IF 575 ! 576 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! 2nd order centered at top & bottom 577 ikt = mikt(ji,jj) + 1 ! w-point below the 1st wet point 578 ikb = MAX(mbkt(ji,jj), 2) ! - above the last wet point 579 ! 580 zwd (ji,jj,ikt) = 1._wp ! top 581 zwi (ji,jj,ikt) = 0._wp 582 zws (ji,jj,ikt) = 0._wp 583 zwrm(ji,jj,ikt) = 0.5_wp * ( pt_in(ji,jj,ikt-1) + pt_in(ji,jj,ikt) ) 584 ! 585 zwd (ji,jj,ikb) = 1._wp ! bottom 586 zwi (ji,jj,ikb) = 0._wp 587 zws (ji,jj,ikb) = 0._wp 588 zwrm(ji,jj,ikb) = 0.5_wp * ( pt_in(ji,jj,ikb-1) + pt_in(ji,jj,ikb) ) 589 END_2D 590 ! 591 ! !== tridiagonal solver ==! 592 ! 593 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) !* 1st recurrence: Tk = Dk - Ik Sk-1 / Tk-1 594 zwt(ji,jj,2) = zwd(ji,jj,2) 595 END_2D 596 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 3, jpkm1 ) 597 zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1) 598 END_3D 599 ! 600 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) !* 2nd recurrence: Zk = Yk - Ik / Tk-1 Zk-1 601 pt_out(ji,jj,2) = zwrm(ji,jj,2) 602 END_2D 603 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 3, jpkm1 ) 604 pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1) 605 END_3D 606 607 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) !* 3d recurrence: Xk = (Zk - Sk Xk+1 ) / Tk 608 pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) 609 END_2D 610 DO_3DS( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpk-2, 2, -1 ) 611 pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk) 612 END_3D 613 ! 614 END SUBROUTINE interp_4th_cpt 615 616 617 SUBROUTINE tridia_solver( pD, pU, pL, pRHS, pt_out , klev ) 618 !!---------------------------------------------------------------------- 619 !! *** ROUTINE tridia_solver *** 620 !! 621 !! ** Purpose : solve a symmetric 3diagonal system 622 !! 623 !! ** Method : solve M.t_out = RHS(t) where M is a tri diagonal matrix ( jpk*jpk ) 624 !! 625 !! ( D_1 U_1 0 0 0 )( t_1 ) ( RHS_1 ) 626 !! ( L_2 D_2 U_2 0 0 )( t_2 ) ( RHS_2 ) 627 !! ( 0 L_3 D_3 U_3 0 )( t_3 ) = ( RHS_3 ) 628 !! ( ... )( ... ) ( ... ) 629 !! ( 0 0 0 L_k D_k )( t_k ) ( RHS_k ) 630 !! 631 !! M is decomposed in the product of an upper and lower triangular matrix. 632 !! The tri-diagonals matrix is given as input 3D arrays: pD, pU, pL 633 !! (i.e. the Diagonal, the Upper diagonal, and the Lower diagonal). 634 !! The solution is pta. 635 !! The 3d array zwt is used as a work space array. 636 !!---------------------------------------------------------------------- 637 REAL(wp),DIMENSION(A2D(nn_hls),jpk), INTENT(in ) :: pD, pU, PL ! 3-diagonal matrix 638 REAL(wp),DIMENSION(A2D(nn_hls),jpk), INTENT(in ) :: pRHS ! Right-Hand-Side 639 REAL(wp),DIMENSION(A2D(nn_hls),jpk), INTENT( out) :: pt_out !!gm field at level=F(klev) 640 INTEGER , INTENT(in ) :: klev ! =1 pt_out at w-level 641 ! ! =0 pt at t-level 642 INTEGER :: ji, jj, jk ! dummy loop integers 643 INTEGER :: kstart ! local indices 644 REAL(wp),DIMENSION(A2D(nn_hls),jpk) :: zwt ! 3D work array 645 !!---------------------------------------------------------------------- 646 ! 647 kstart = 1 + klev 648 ! 649 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) !* 1st recurrence: Tk = Dk - Ik Sk-1 / Tk-1 650 zwt(ji,jj,kstart) = pD(ji,jj,kstart) 651 END_2D 652 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, kstart+1, jpkm1 ) 653 zwt(ji,jj,jk) = pD(ji,jj,jk) - pL(ji,jj,jk) * pU(ji,jj,jk-1) /zwt(ji,jj,jk-1) 654 END_3D 655 ! 656 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) !* 2nd recurrence: Zk = Yk - Ik / Tk-1 Zk-1 657 pt_out(ji,jj,kstart) = pRHS(ji,jj,kstart) 658 END_2D 659 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, kstart+1, jpkm1 ) 660 pt_out(ji,jj,jk) = pRHS(ji,jj,jk) - pL(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1) 661 END_3D 662 663 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) !* 3d recurrence: Xk = (Zk - Sk Xk+1 ) / Tk 664 pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) 665 END_2D 666 DO_3DS( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpk-2, kstart, -1 ) 667 pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - pU(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk) 668 END_3D 669 ! 670 END SUBROUTINE tridia_solver 671 1005 END SUBROUTINE tra_adv_fct_lf 1006 #endif 672 1007 !!====================================================================== 673 1008 END MODULE traadv_fct
Note: See TracChangeset
for help on using the changeset viewer.
