Changeset 14015 for NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/TRA/traadv_fct.F90
- Timestamp:
- 2020-12-02T16:26:43+01:00 (4 years ago)
- Location:
- NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source
- Property svn:externals
-
old new 8 8 9 9 # SETTE 10 ^/utils/CI/sette @13559sette10 ^/utils/CI/sette_wave@13990 sette
-
- Property svn:externals
-
NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/TRA/traadv_fct.F90
r13497 r14015 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.F90 37 PUBLIC nonosc ! called by traadv_fct_lf.F90 - key_agrif 36 38 37 39 LOGICAL :: l_trd ! flag to compute trends … … 79 81 INTEGER , INTENT(in ) :: kn_fct_v ! order of the FCT scheme (=2 or 4) 80 82 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 83 ! TEMP: [tiling] This can be A2D(nn_hls) if using XIOS (subdomain support) 84 ! NOTE: [tiling-comms-merge] These were changed to INTENT(inout) but they are not modified, so it is reverted 81 85 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pU, pV, pW ! 3 ocean volume flux components 82 86 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! tracers and RHS of tracer equation 83 87 ! 84 INTEGER :: ji, jj, jk, jn ! dummy loop indices 88 INTEGER :: ji, jj, jk, jn ! dummy loop indices 85 89 REAL(wp) :: ztra ! local scalar 86 90 REAL(wp) :: zfp_ui, zfp_vj, zfp_wk, zC2t_u, zC4t_u ! - - 87 91 REAL(wp) :: zfm_ui, zfm_vj, zfm_wk, zC2t_v, zC4t_v ! - - 88 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw92 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw 89 93 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdx, ztrdy, ztrdz, zptry 90 94 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zwinf, zwdia, zwsup … … 92 96 !!---------------------------------------------------------------------- 93 97 ! 94 IF( kt == kit000 ) THEN 95 IF(lwp) WRITE(numout,*) 96 IF(lwp) WRITE(numout,*) 'tra_adv_fct : FCT advection scheme on ', cdtype 97 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 98 IF( ntile == 0 .OR. ntile == 1 ) THEN ! Do only on the first tile 99 IF( kt == kit000 ) THEN 100 IF(lwp) WRITE(numout,*) 101 IF(lwp) WRITE(numout,*) 'tra_adv_fct : FCT advection scheme on ', cdtype 102 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 103 ENDIF 104 ! NOTE: [tiling-comms-merge] Bug fix- move array zeroing out of this IF block 105 ! 106 l_trd = .FALSE. ! set local switches 107 l_hst = .FALSE. 108 l_ptr = .FALSE. 109 ll_zAimp = .FALSE. 110 IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype =='TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 111 IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) ) l_ptr = .TRUE. 112 IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. & 113 & iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) ) l_hst = .TRUE. 114 ! 98 115 ENDIF 116 99 117 !! -- init to 0 100 118 zwi(:,:,:) = 0._wp … … 108 126 ztw(:,:,:) = 0._wp 109 127 ! 110 l_trd = .FALSE. ! set local switches111 l_hst = .FALSE.112 l_ptr = .FALSE.113 ll_zAimp = .FALSE.114 IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype =='TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.115 IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) ) l_ptr = .TRUE.116 IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &117 & iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) ) l_hst = .TRUE.118 !119 128 IF( l_trd .OR. l_hst ) THEN 120 ALLOCATE( ztrdx( jpi,jpj,jpk), ztrdy(jpi,jpj,jpk), ztrdz(jpi,jpj,jpk) )129 ALLOCATE( ztrdx(A2D(nn_hls),jpk), ztrdy(A2D(nn_hls),jpk), ztrdz(A2D(nn_hls),jpk) ) 121 130 ztrdx(:,:,:) = 0._wp ; ztrdy(:,:,:) = 0._wp ; ztrdz(:,:,:) = 0._wp 122 131 ENDIF 123 132 ! 124 IF( l_ptr ) THEN 125 ALLOCATE( zptry( jpi,jpj,jpk) )133 IF( l_ptr ) THEN 134 ALLOCATE( zptry(A2D(nn_hls),jpk) ) 126 135 zptry(:,:,:) = 0._wp 127 136 ENDIF 128 ! ! surface & bottom value : flux set to zero one for all129 zwz(:,:, 1 ) = 0._wp130 zwx(:,:,jpk) = 0._wp ; zwy(:,:,jpk) = 0._wp ; zwz(:,:,jpk) = 0._wp131 !132 zwi(:,:,:) = 0._wp133 137 ! 134 138 ! If adaptive vertical advection, check if it is needed on this PE at this time 135 139 IF( ln_zad_Aimp ) THEN 136 IF( MAXVAL( ABS( wi( :,:,:) ) ) > 0._wp ) ll_zAimp = .TRUE.140 IF( MAXVAL( ABS( wi(A2D(nn_hls),:) ) ) > 0._wp ) ll_zAimp = .TRUE. 137 141 END IF 138 142 ! If active adaptive vertical advection, build tridiagonal matrix 139 143 IF( ll_zAimp ) THEN 140 ALLOCATE(zwdia( jpi,jpj,jpk), zwinf(jpi,jpj,jpk),zwsup(jpi,jpj,jpk))141 DO_3D( 0, 0, 0, 0, 1, jpkm1 )144 ALLOCATE(zwdia(A2D(nn_hls),jpk), zwinf(A2D(nn_hls),jpk), zwsup(A2D(nn_hls),jpk)) 145 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 142 146 zwdia(ji,jj,jk) = 1._wp + p2dt * ( MAX( wi(ji,jj,jk) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) ) & 143 147 & / e3t(ji,jj,jk,Krhs) … … 151 155 ! !== upstream advection with initial mass fluxes & intermediate update ==! 152 156 ! !* upstream tracer flux in the i and j direction 153 DO_3D( 1, 0, 1, 0, 1, jpkm1 )157 DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 ) 154 158 ! upstream scheme 155 159 zfp_ui = pU(ji,jj,jk) + ABS( pU(ji,jj,jk) ) … … 178 182 ENDIF 179 183 ! 180 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !* trend and after field with monotonic scheme184 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) !* trend and after field with monotonic scheme 181 185 ! ! total intermediate advective trends 182 186 ztra = - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & … … 194 198 ! 195 199 ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp ; 196 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! Interior value ( multiplied by wmask)200 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! Interior value ( multiplied by wmask) 197 201 zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) ) 198 202 zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) ) … … 206 210 ! 207 211 END IF 208 ! 212 ! 209 213 IF( l_trd .OR. l_hst ) THEN ! trend diagnostics (contribution of upstream fluxes) 210 214 ztrdx(:,:,:) = zwx(:,:,:) ; ztrdy(:,:,:) = zwy(:,:,:) ; ztrdz(:,:,:) = zwz(:,:,:) … … 218 222 ! 219 223 CASE( 2 ) !- 2nd order centered 220 DO_3D( 1, 0, 1, 0, 1, jpkm1 )224 DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 ) 221 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) 222 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) … … 238 242 CALL lbc_lnk_multi( 'traadv_fct', zltu, 'T', 1.0_wp , zltv, 'T', 1.0_wp ) ! Lateral boundary cond. (unchanged sgn) 239 243 ! 240 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) ! Horizontal advective fluxes244 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 241 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 242 246 zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji ,jj+1,jk,jn,Kmm) 243 ! ! C4 minus upstream advective fluxes 247 ! ! C4 minus upstream advective fluxes 244 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) 245 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) 246 250 END_3D 251 IF (nn_hls.EQ.2) CALL lbc_lnk_multi( 'traadv_fct', zwx, 'U', -1.0_wp, zwy, 'V', -1.0_wp ) ! Lateral boundary cond. (unchanged sgn) 247 252 ! 248 253 CASE( 41 ) !- 4th order centered ==>> !!gm coding attempt need to be tested 249 254 ztu(:,:,jpk) = 0._wp ! Bottom value : flux set to zero 250 255 ztv(:,:,jpk) = 0._wp 251 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) ! 1st derivative (gradient)256 DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 ) ! 1st derivative (gradient) 252 257 ztu(ji,jj,jk) = ( pt(ji+1,jj ,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk) 253 258 ztv(ji,jj,jk) = ( pt(ji ,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk) 254 259 END_3D 255 CALL lbc_lnk_multi( 'traadv_fct', ztu, 'U', -1.0_wp , ztv, 'V', -1.0_wp ) ! Lateral boundary cond. (unchanged sgn) 260 IF (nn_hls.EQ.1) CALL lbc_lnk_multi( 'traadv_fct', ztu, 'U', -1.0_wp , ztv, 'V', -1.0_wp ) ! Lateral boundary cond. (unchanged sgn) 261 ! 262 IF (nn_hls.EQ.1) CALL lbc_lnk_multi( 'traadv_fct', ztu, 'U', -1.0_wp , ztv, 'V', -1.0_wp ) ! Lateral boundary cond. (unchanged sgn) 256 263 ! 257 264 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! Horizontal advective fluxes … … 265 272 zwy(ji,jj,jk) = 0.5_wp * pV(ji,jj,jk) * zC4t_v - zwy(ji,jj,jk) 266 273 END_3D 274 IF (nn_hls.EQ.2) CALL lbc_lnk_multi( 'traadv_fct', zwx, 'U', -1.0_wp , zwy, 'V', -1.0_wp ) ! Lateral boundary cond. (unchanged sgn) 267 275 ! 268 276 END SELECT … … 271 279 ! 272 280 CASE( 2 ) !- 2nd order centered 273 DO_3D( 0, 0, 0, 0, 2, jpkm1 )281 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 274 282 zwz(ji,jj,jk) = ( pW(ji,jj,jk) * 0.5_wp * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj,jk-1,jn,Kmm) ) & 275 283 & - zwz(ji,jj,jk) ) * wmask(ji,jj,jk) … … 278 286 CASE( 4 ) !- 4th order COMPACT 279 287 CALL interp_4th_cpt( pt(:,:,:,jn,Kmm) , ztw ) ! zwt = COMPACT interpolation of T at w-point 280 DO_3D( 0, 0, 0, 0, 2, jpkm1 )288 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 281 289 zwz(ji,jj,jk) = ( pW(ji,jj,jk) * ztw(ji,jj,jk) - zwz(ji,jj,jk) ) * wmask(ji,jj,jk) 282 290 END_3D … … 286 294 zwz(:,:,1) = 0._wp ! only ocean surface as interior zwz values have been w-masked 287 295 ENDIF 288 ! 296 ! 297 IF (nn_hls.EQ.1) THEN 298 CALL lbc_lnk_multi( 'traadv_fct', zwi, 'T', 1.0_wp, zwx, 'U', -1.0_wp , zwy, 'V', -1.0_wp, zwz, 'T', 1.0_wp ) 299 ELSE 300 CALL lbc_lnk( 'traadv_fct', zwi, 'T', 1.0_wp) 301 END IF 302 ! 303 IF (nn_hls.EQ.1) THEN 304 CALL lbc_lnk_multi( 'traadv_fct', zwi, 'T', 1.0_wp, 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 ! 289 309 IF ( ll_zAimp ) THEN 290 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !* trend and after field with monotonic scheme310 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) !* trend and after field with monotonic scheme 291 311 ! ! total intermediate advective trends 292 312 ztra = - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 293 313 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & 294 314 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) * r1_e1e2t(ji,jj) 295 ztw(ji,jj,jk) 315 ztw(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk) 296 316 END_3D 297 317 ! 298 318 CALL tridia_solver( zwdia, zwsup, zwinf, ztw, ztw , 0 ) 299 319 ! 300 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! Interior value ( multiplied by wmask)320 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) ! Interior value ( multiplied by wmask) 301 321 zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) ) 302 322 zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) ) 303 zwz(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) 304 324 END_3D 305 325 END IF 306 !307 CALL lbc_lnk_multi( 'traadv_fct', zwi, 'T', 1.0_wp, zwx, 'U', -1.0_wp , zwy, 'V', -1.0_wp, zwz, 'W', 1.0_wp )308 326 ! 309 327 ! !== monotonicity algorithm ==! … … 334 352 & * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 335 353 END_3D 336 END IF 337 ! 354 END IF 355 ! NOTE: [tiling-comms-merge] I tested this 356 ! NOT TESTED - NEED l_trd OR l_hst TRUE 338 357 IF( l_trd .OR. l_hst ) THEN ! trend diagnostics // heat/salt transport 339 ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:) ! <<< add anti-diffusive fluxes 358 ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:) ! <<< add anti-diffusive fluxes 340 359 ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:) ! to upstream fluxes 341 360 ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:) ! … … 350 369 ! 351 370 ENDIF 371 ! NOTE: [tiling-comms-merge] I tested this 372 ! NOT TESTED - NEED l_ptr TRUE 352 373 IF( l_ptr ) THEN ! "Poleward" transports 353 374 zptry(:,:,:) = zptry(:,:,:) + zwy(:,:,:) ! <<< add anti-diffusive fluxes … … 360 381 DEALLOCATE( zwdia, zwinf, zwsup ) 361 382 ENDIF 362 IF( l_trd .OR. l_hst ) THEN 383 IF( l_trd .OR. l_hst ) THEN 363 384 DEALLOCATE( ztrdx, ztrdy, ztrdz ) 364 385 ENDIF … … 383 404 !! in-space based differencing for fluid 384 405 !!---------------------------------------------------------------------- 385 INTEGER , INTENT(in ) :: Kmm ! time level index 386 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 387 REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(in ) :: pbef, paft ! before & after field 388 REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(inout) :: paa, pbb, pcc ! monotonic fluxes in the 3 directions 406 INTEGER , INTENT(in ) :: Kmm ! time level index 407 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 408 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pbef ! before field 409 REAL(wp), DIMENSION(A2D(nn_hls) ,jpk), INTENT(in ) :: paft ! after field 410 REAL(wp), DIMENSION(A2D(nn_hls) ,jpk), INTENT(inout) :: paa, pbb, pcc ! monotonic fluxes in the 3 directions 389 411 ! 390 412 INTEGER :: ji, jj, jk ! dummy loop indices … … 392 414 REAL(dp) :: zpos, zneg, zbt, za, zb, zc, zbig, zrtrn ! local scalars 393 415 REAL(dp) :: zau, zbu, zcu, zav, zbv, zcv, zup, zdo ! - - 394 REAL(dp), DIMENSION( jpi,jpj,jpk) :: zbetup, zbetdo, zbup, zbdo416 REAL(dp), DIMENSION(A2D(nn_hls),jpk) :: zbetup, zbetdo, zbup, zbdo 395 417 !!---------------------------------------------------------------------- 396 418 ! … … 402 424 ! -------------------- 403 425 ! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land 404 zbup = MAX( pbef * tmask - zbig * ( 1._wp - tmask ), & 405 & paft * tmask - zbig * ( 1._wp - tmask ) ) 406 zbdo = MIN( pbef * tmask + zbig * ( 1._wp - tmask ), & 407 & paft * tmask + zbig * ( 1._wp - tmask ) ) 426 DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk ) 427 zbup(ji,jj,jk) = MAX( pbef(ji,jj,jk) * tmask(ji,jj,jk) - zbig * ( 1._wp - tmask(ji,jj,jk) ), & 428 & paft(ji,jj,jk) * tmask(ji,jj,jk) - zbig * ( 1._wp - tmask(ji,jj,jk) ) ) 429 zbdo(ji,jj,jk) = MIN( pbef(ji,jj,jk) * tmask(ji,jj,jk) + zbig * ( 1._wp - tmask(ji,jj,jk) ), & 430 & paft(ji,jj,jk) * tmask(ji,jj,jk) + zbig * ( 1._wp - tmask(ji,jj,jk) ) ) 431 END_3D 408 432 409 433 DO jk = 1, jpkm1 410 434 ikm1 = MAX(jk-1,1) 411 DO_2D( 0, 0, 0, 0)435 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 412 436 413 437 ! search maximum in neighbourhood … … 439 463 END_2D 440 464 END DO 441 CALL lbc_lnk_multi( 'traadv_fct', zbetup, 'T', 1.0_wp , zbetdo, 'T', 1.0_wp ) ! lateral boundary cond. (unchanged sign)465 IF (nn_hls.EQ.1) CALL lbc_lnk_multi( 'traadv_fct', zbetup, 'T', 1.0_wp , zbetdo, 'T', 1.0_wp ) ! lateral boundary cond. (unchanged sign) 442 466 443 467 ! 3. monotonic flux in the i & j direction (paa & pbb) 444 468 ! ---------------------------------------- 445 DO_3D( 0, 0, 0, 0, 1, jpkm1 )469 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 446 470 zau = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) ) 447 471 zbu = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) ) … … 461 485 pcc(ji,jj,jk+1) = pcc(ji,jj,jk+1) * ( zc * za + ( 1._wp - zc) * zb ) 462 486 END_3D 463 CALL lbc_lnk_multi( 'traadv_fct', paa, 'U', -1.0_wp , pbb, 'V', -1.0_wp ) ! lateral boundary condition (changed sign)464 487 ! 465 488 END SUBROUTINE nonosc … … 537 560 !!---------------------------------------------------------------------- 538 561 REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pt_in ! field at t-point 539 REAL(wp),DIMENSION( jpi,jpj,jpk), INTENT( out) :: pt_out ! field interpolated at w-point562 REAL(wp),DIMENSION(A2D(nn_hls) ,jpk), INTENT( out) :: pt_out ! field interpolated at w-point 540 563 ! 541 564 INTEGER :: ji, jj, jk ! dummy loop integers 542 565 INTEGER :: ikt, ikb ! local integers 543 REAL(wp),DIMENSION( jpi,jpj,jpk) :: zwd, zwi, zws, zwrm, zwt566 REAL(wp),DIMENSION(A2D(nn_hls),jpk) :: zwd, zwi, zws, zwrm, zwt 544 567 !!---------------------------------------------------------------------- 545 568 ! 546 569 ! !== build the three diagonal matrix & the RHS ==! 547 570 ! 548 DO_3D( 0, 0, 0, 0, 3, jpkm1 ) ! interior (from jk=3 to jpk-1)571 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 3, jpkm1 ) ! interior (from jk=3 to jpk-1) 549 572 zwd (ji,jj,jk) = 3._wp * wmask(ji,jj,jk) + 1._wp ! diagonal 550 573 zwi (ji,jj,jk) = wmask(ji,jj,jk) ! lower diagonal … … 565 588 END IF 566 589 ! 567 DO_2D( 0, 0, 0, 0) ! 2nd order centered at top & bottom590 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! 2nd order centered at top & bottom 568 591 ikt = mikt(ji,jj) + 1 ! w-point below the 1st wet point 569 592 ikb = MAX(mbkt(ji,jj), 2) ! - above the last wet point … … 582 605 ! !== tridiagonal solver ==! 583 606 ! 584 DO_2D( 0, 0, 0, 0) !* 1st recurrence: Tk = Dk - Ik Sk-1 / Tk-1607 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) !* 1st recurrence: Tk = Dk - Ik Sk-1 / Tk-1 585 608 zwt(ji,jj,2) = zwd(ji,jj,2) 586 609 END_2D 587 DO_3D( 0, 0, 0, 0, 3, jpkm1 )610 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 3, jpkm1 ) 588 611 zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1) 589 612 END_3D 590 613 ! 591 DO_2D( 0, 0, 0, 0) !* 2nd recurrence: Zk = Yk - Ik / Tk-1 Zk-1614 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) !* 2nd recurrence: Zk = Yk - Ik / Tk-1 Zk-1 592 615 pt_out(ji,jj,2) = zwrm(ji,jj,2) 593 616 END_2D 594 DO_3D( 0, 0, 0, 0, 3, jpkm1 )617 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 3, jpkm1 ) 595 618 pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1) 596 619 END_3D 597 620 598 DO_2D( 0, 0, 0, 0) !* 3d recurrence: Xk = (Zk - Sk Xk+1 ) / Tk621 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) !* 3d recurrence: Xk = (Zk - Sk Xk+1 ) / Tk 599 622 pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) 600 623 END_2D 601 DO_3DS( 0, 0, 0, 0, jpk-2, 2, -1 )624 DO_3DS( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpk-2, 2, -1 ) 602 625 pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk) 603 626 END_3D … … 626 649 !! The 3d array zwt is used as a work space array. 627 650 !!---------------------------------------------------------------------- 628 REAL(wp),DIMENSION( :,:,:), INTENT(in ) :: pD, pU, PL ! 3-diagonal matrix629 REAL(wp),DIMENSION( :,:,:), INTENT(in ) :: pRHS ! Right-Hand-Side630 REAL(wp),DIMENSION( :,:,:), INTENT( out) :: pt_out !!gm field at level=F(klev)631 INTEGER , INTENT(in ) :: klev ! =1 pt_out at w-level632 ! ! =0 pt at t-level651 REAL(wp),DIMENSION(A2D(nn_hls),jpk), INTENT(in ) :: pD, pU, PL ! 3-diagonal matrix 652 REAL(wp),DIMENSION(A2D(nn_hls),jpk), INTENT(in ) :: pRHS ! Right-Hand-Side 653 REAL(wp),DIMENSION(A2D(nn_hls),jpk), INTENT( out) :: pt_out !!gm field at level=F(klev) 654 INTEGER , INTENT(in ) :: klev ! =1 pt_out at w-level 655 ! ! =0 pt at t-level 633 656 INTEGER :: ji, jj, jk ! dummy loop integers 634 657 INTEGER :: kstart ! local indices 635 REAL(wp),DIMENSION( jpi,jpj,jpk) :: zwt ! 3D work array658 REAL(wp),DIMENSION(A2D(nn_hls),jpk) :: zwt ! 3D work array 636 659 !!---------------------------------------------------------------------- 637 660 ! 638 661 kstart = 1 + klev 639 662 ! 640 DO_2D( 0, 0, 0, 0) !* 1st recurrence: Tk = Dk - Ik Sk-1 / Tk-1663 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) !* 1st recurrence: Tk = Dk - Ik Sk-1 / Tk-1 641 664 zwt(ji,jj,kstart) = pD(ji,jj,kstart) 642 665 END_2D 643 DO_3D( 0, 0, 0, 0, kstart+1, jpkm1 )666 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, kstart+1, jpkm1 ) 644 667 zwt(ji,jj,jk) = pD(ji,jj,jk) - pL(ji,jj,jk) * pU(ji,jj,jk-1) /zwt(ji,jj,jk-1) 645 668 END_3D 646 669 ! 647 DO_2D( 0, 0, 0, 0) !* 2nd recurrence: Zk = Yk - Ik / Tk-1 Zk-1670 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) !* 2nd recurrence: Zk = Yk - Ik / Tk-1 Zk-1 648 671 pt_out(ji,jj,kstart) = pRHS(ji,jj,kstart) 649 672 END_2D 650 DO_3D( 0, 0, 0, 0, kstart+1, jpkm1 )673 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, kstart+1, jpkm1 ) 651 674 pt_out(ji,jj,jk) = pRHS(ji,jj,jk) - pL(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1) 652 675 END_3D 653 676 654 DO_2D( 0, 0, 0, 0) !* 3d recurrence: Xk = (Zk - Sk Xk+1 ) / Tk677 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) !* 3d recurrence: Xk = (Zk - Sk Xk+1 ) / Tk 655 678 pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) 656 679 END_2D 657 DO_3DS( 0, 0, 0, 0, jpk-2, kstart, -1 )680 DO_3DS( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpk-2, kstart, -1 ) 658 681 pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - pU(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk) 659 682 END_3D
Note: See TracChangeset
for help on using the changeset viewer.