Changeset 13982 for NEMO/trunk/src/OCE/TRA/traadv_qck.F90
- Timestamp:
- 2020-12-02T11:57:05+01:00 (3 years ago)
- File:
-
- 1 edited
-
NEMO/trunk/src/OCE/TRA/traadv_qck.F90 (modified) (17 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/OCE/TRA/traadv_qck.F90
r13497 r13982 91 91 INTEGER , INTENT(in ) :: kjpt ! number of tracers 92 92 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 93 ! TEMP: [tiling] This can be A2D(nn_hls) if using XIOS (subdomain support) 94 ! NOTE: [tiling-comms-merge] These were changed to INTENT(inout) but they are not modified, so it is reverted 93 95 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pU, pV, pW ! 3 ocean volume transport components 94 96 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! tracers and RHS of tracer equation 95 97 !!---------------------------------------------------------------------- 96 98 ! 97 IF( kt == kit000 ) THEN 98 IF(lwp) WRITE(numout,*) 99 IF(lwp) WRITE(numout,*) 'tra_adv_qck : 3rd order quickest advection scheme on ', cdtype 100 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 101 IF(lwp) WRITE(numout,*) 99 IF( ntile == 0 .OR. ntile == 1 ) THEN ! Do only on the first tile 100 IF( kt == kit000 ) THEN 101 IF(lwp) WRITE(numout,*) 102 IF(lwp) WRITE(numout,*) 'tra_adv_qck : 3rd order quickest advection scheme on ', cdtype 103 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 104 IF(lwp) WRITE(numout,*) 105 ENDIF 106 ! 107 l_trd = .FALSE. 108 l_ptr = .FALSE. 109 IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 110 IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) ) l_ptr = .TRUE. 102 111 ENDIF 103 !104 l_trd = .FALSE.105 l_ptr = .FALSE.106 IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.107 IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) ) l_ptr = .TRUE.108 !109 112 ! 110 113 ! ! horizontal fluxes are computed with the QUICKEST + ULTIMATE scheme … … 127 130 INTEGER , INTENT(in ) :: kjpt ! number of tracers 128 131 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 132 ! TEMP: [tiling] This can be A2D(nn_hls) if using XIOS (subdomain support) 129 133 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pU ! i-velocity components 130 134 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! active tracers and RHS of tracer equation … … 132 136 INTEGER :: ji, jj, jk, jn ! dummy loop indices 133 137 REAL(wp) :: ztra, zbtr, zdir, zdx, zmsk ! local scalars 134 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zwx, zfu, zfc, zfd138 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zwx, zfu, zfc, zfd 135 139 !---------------------------------------------------------------------- 136 140 ! … … 142 146 ! 143 147 !!gm why not using a SHIFT instruction... 144 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !--- Computation of the ustream and downstream value of the tracer and the mask148 DO_3D( 0, 0, nn_hls-1, nn_hls-1, 1, jpkm1 ) !--- Computation of the ustream and downstream value of the tracer and the mask 145 149 zfc(ji,jj,jk) = pt(ji-1,jj,jk,jn,Kbb) ! Upstream in the x-direction for the tracer 146 150 zfd(ji,jj,jk) = pt(ji+1,jj,jk,jn,Kbb) ! Downstream in the x-direction for the tracer 147 151 END_3D 148 CALL lbc_lnk_multi( 'traadv_qck', zfc(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp ) ! Lateral boundary conditions152 IF (nn_hls.EQ.1) CALL lbc_lnk_multi( 'traadv_qck', zfc(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp ) ! Lateral boundary conditions 149 153 150 154 ! 151 155 ! Horizontal advective fluxes 152 156 ! --------------------------- 153 DO_3D( 0, 0, 0, 0, 1, jpkm1 )157 DO_3D( 0, 0, nn_hls-1, 0, 1, jpkm1 ) 154 158 zdir = 0.5 + SIGN( 0.5_wp, pU(ji,jj,jk) ) ! if pU > 0 : zdir = 1 otherwise zdir = 0 155 159 zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji+1,jj,jk) ! FU in the x-direction for T 156 160 END_3D 157 161 ! 158 DO_3D( 0, 0, 0, 0, 1, jpkm1 )162 DO_3D( 0, 0, nn_hls-1, 0, 1, jpkm1 ) 159 163 zdir = 0.5 + SIGN( 0.5_wp, pU(ji,jj,jk) ) ! if pU > 0 : zdir = 1 otherwise zdir = 0 160 164 zdx = ( zdir * e1t(ji,jj) + ( 1. - zdir ) * e1t(ji+1,jj) ) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm) … … 164 168 END_3D 165 169 !--- Lateral boundary conditions 166 CALL lbc_lnk_multi( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp, zfc(:,:,:), 'T', 1.0_wp, zwx(:,:,:), 'T', 1.0_wp )170 IF (nn_hls.EQ.1) CALL lbc_lnk_multi( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp, zfc(:,:,:), 'T', 1.0_wp, zwx(:,:,:), 'T', 1.0_wp ) 167 171 168 172 !--- QUICKEST scheme … … 170 174 ! 171 175 ! Mask at the T-points in the x-direction (mask=0 or mask=1) 172 DO_3D( 0, 0, 0, 0, 1, jpkm1 )176 DO_3D( 0, 0, nn_hls-1, nn_hls-1, 1, jpkm1 ) 173 177 zfu(ji,jj,jk) = tmask(ji-1,jj,jk) + tmask(ji,jj,jk) + tmask(ji+1,jj,jk) - 2. 174 178 END_3D 175 CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp ) ! Lateral boundary conditions179 IF (nn_hls.EQ.1) CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp ) ! Lateral boundary conditions 176 180 177 181 ! 178 182 ! Tracer flux on the x-direction 179 DO jk = 1, jpkm1 180 ! 181 DO_2D( 0, 0, 0, 0 ) 182 zdir = 0.5 + SIGN( 0.5_wp, pU(ji,jj,jk) ) ! if pU > 0 : zdir = 1 otherwise zdir = 0 183 !--- If the second ustream point is a land point 184 !--- the flux is computed by the 1st order UPWIND scheme 185 zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji+1,jj,jk) 186 zwx(ji,jj,jk) = zmsk * zwx(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk) 187 zwx(ji,jj,jk) = zwx(ji,jj,jk) * pU(ji,jj,jk) 188 END_2D 189 END DO 190 ! 191 CALL lbc_lnk( 'traadv_qck', zwx(:,:,:), 'T', 1.0_wp ) ! Lateral boundary conditions 183 DO_3D( 0, 0, 1, 0, 1, jpkm1 ) 184 zdir = 0.5 + SIGN( 0.5_wp, pU(ji,jj,jk) ) ! if pU > 0 : zdir = 1 otherwise zdir = 0 185 !--- If the second ustream point is a land point 186 !--- the flux is computed by the 1st order UPWIND scheme 187 zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji+1,jj,jk) 188 zwx(ji,jj,jk) = zmsk * zwx(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk) 189 zwx(ji,jj,jk) = zwx(ji,jj,jk) * pU(ji,jj,jk) 190 END_3D 192 191 ! 193 192 ! Computation of the trend … … 216 215 INTEGER , INTENT(in ) :: kjpt ! number of tracers 217 216 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 217 ! TEMP: [tiling] This can be A2D(nn_hls) if using XIOS (subdomain support) 218 218 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pV ! j-velocity components 219 219 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! active tracers and RHS of tracer equation … … 221 221 INTEGER :: ji, jj, jk, jn ! dummy loop indices 222 222 REAL(wp) :: ztra, zbtr, zdir, zdx, zmsk ! local scalars 223 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zwy, zfu, zfc, zfd ! 3D workspace223 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zwy, zfu, zfc, zfd ! 3D workspace 224 224 !---------------------------------------------------------------------- 225 225 ! … … 233 233 ! 234 234 !--- Computation of the ustream and downstream value of the tracer and the mask 235 DO_2D( 0, 0, 0, 0 )235 DO_2D( nn_hls-1, nn_hls-1, 0, 0 ) 236 236 ! Upstream in the x-direction for the tracer 237 237 zfc(ji,jj,jk) = pt(ji,jj-1,jk,jn,Kbb) … … 240 240 END_2D 241 241 END DO 242 CALL lbc_lnk_multi( 'traadv_qck', zfc(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp ) ! Lateral boundary conditions 243 242 IF (nn_hls.EQ.1) CALL lbc_lnk_multi( 'traadv_qck', zfc(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp ) ! Lateral boundary conditions 244 243 245 244 ! … … 247 246 ! --------------------------- 248 247 ! 249 DO_3D( 0, 0, 0, 0, 1, jpkm1 )248 DO_3D( nn_hls-1, 0, 0, 0, 1, jpkm1 ) 250 249 zdir = 0.5 + SIGN( 0.5_wp, pV(ji,jj,jk) ) ! if pU > 0 : zdir = 1 otherwise zdir = 0 251 250 zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji,jj+1,jk) ! FU in the x-direction for T 252 251 END_3D 253 252 ! 254 DO_3D( 0, 0, 0, 0, 1, jpkm1 )253 DO_3D( nn_hls-1, 0, 0, 0, 1, jpkm1 ) 255 254 zdir = 0.5 + SIGN( 0.5_wp, pV(ji,jj,jk) ) ! if pU > 0 : zdir = 1 otherwise zdir = 0 256 255 zdx = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * e3v(ji,jj,jk,Kmm) … … 261 260 262 261 !--- Lateral boundary conditions 263 CALL lbc_lnk_multi( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp, zfc(:,:,:), 'T', 1.0_wp, zwy(:,:,:), 'T', 1.0_wp )262 IF (nn_hls.EQ.1) CALL lbc_lnk_multi( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp, zfc(:,:,:), 'T', 1.0_wp, zwy(:,:,:), 'T', 1.0_wp ) 264 263 265 264 !--- QUICKEST scheme … … 267 266 ! 268 267 ! Mask at the T-points in the x-direction (mask=0 or mask=1) 269 DO_3D( 0, 0, 0, 0, 1, jpkm1 )268 DO_3D( nn_hls-1, nn_hls-1, 0, 0, 1, jpkm1 ) 270 269 zfu(ji,jj,jk) = tmask(ji,jj-1,jk) + tmask(ji,jj,jk) + tmask(ji,jj+1,jk) - 2. 271 270 END_3D 272 CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp ) !--- Lateral boundary conditions271 IF (nn_hls.EQ.1) CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp ) !--- Lateral boundary conditions 273 272 ! 274 273 ! Tracer flux on the x-direction 275 DO jk = 1, jpkm1 276 ! 277 DO_2D( 0, 0, 0, 0 ) 278 zdir = 0.5 + SIGN( 0.5_wp, pV(ji,jj,jk) ) ! if pU > 0 : zdir = 1 otherwise zdir = 0 279 !--- If the second ustream point is a land point 280 !--- the flux is computed by the 1st order UPWIND scheme 281 zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji,jj+1,jk) 282 zwy(ji,jj,jk) = zmsk * zwy(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk) 283 zwy(ji,jj,jk) = zwy(ji,jj,jk) * pV(ji,jj,jk) 284 END_2D 285 END DO 286 ! 287 CALL lbc_lnk( 'traadv_qck', zwy(:,:,:), 'T', 1.0_wp ) ! Lateral boundary conditions 274 DO_3D( 1, 0, 0, 0, 1, jpkm1 ) 275 zdir = 0.5 + SIGN( 0.5_wp, pV(ji,jj,jk) ) ! if pU > 0 : zdir = 1 otherwise zdir = 0 276 !--- If the second ustream point is a land point 277 !--- the flux is computed by the 1st order UPWIND scheme 278 zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji,jj+1,jk) 279 zwy(ji,jj,jk) = zmsk * zwy(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk) 280 zwy(ji,jj,jk) = zwy(ji,jj,jk) * pV(ji,jj,jk) 281 END_3D 288 282 ! 289 283 ! Computation of the trend … … 313 307 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 314 308 INTEGER , INTENT(in ) :: kjpt ! number of tracers 315 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pW ! vertical velocity 309 ! TEMP: [tiling] This can be A2D(nn_hls) if using XIOS (subdomain support) 310 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pW ! vertical velocity 316 311 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! active tracers and RHS of tracer equation 317 312 ! 318 313 INTEGER :: ji, jj, jk, jn ! dummy loop indices 319 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zwz ! 3D workspace314 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zwz ! 3D workspace 320 315 !!---------------------------------------------------------------------- 321 316 ! … … 332 327 IF( ln_linssh ) THEN !* top value (only in linear free surf. as zwz is multiplied by wmask) 333 328 IF( ln_isfcav ) THEN ! ice-shelf cavities (top of the ocean) 334 DO_2D( 1, 1, 1, 1)329 DO_2D( 0, 0, 0, 0 ) 335 330 zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kmm) ! linear free surface 336 331 END_2D 337 332 ELSE ! no ocean cavities (only ocean surface) 338 zwz(:,:,1) = pW(:,:,1) * pt(:,:,1,jn,Kmm) 333 DO_2D( 0, 0, 0, 0 ) 334 zwz(ji,jj,1) = pW(ji,jj,1) * pt(ji,jj,1,jn,Kmm) 335 END_2D 339 336 ENDIF 340 337 ENDIF … … 359 356 !! ** Method : 360 357 !!---------------------------------------------------------------------- 361 REAL(wp), DIMENSION( jpi,jpj,jpk), INTENT(in ) :: pfu ! second upwind point362 REAL(wp), DIMENSION( jpi,jpj,jpk), INTENT(in ) :: pfd ! first douwning point363 REAL(wp), DIMENSION( jpi,jpj,jpk), INTENT(in ) :: pfc ! the central point (or the first upwind point)364 REAL(wp), DIMENSION( jpi,jpj,jpk), INTENT(inout) :: puc ! input as Courant number ; output as flux358 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in ) :: pfu ! second upwind point 359 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in ) :: pfd ! first douwning point 360 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in ) :: pfc ! the central point (or the first upwind point) 361 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(inout) :: puc ! input as Courant number ; output as flux 365 362 !! 366 363 INTEGER :: ji, jj, jk ! dummy loop indices … … 369 366 !---------------------------------------------------------------------- 370 367 ! 371 DO_3D( 1, 1, 1, 1, 1, jpkm1 )368 DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) 372 369 zc = puc(ji,jj,jk) ! Courant number 373 370 zcurv = pfd(ji,jj,jk) + pfu(ji,jj,jk) - 2. * pfc(ji,jj,jk)
Note: See TracChangeset
for help on using the changeset viewer.
