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