Changeset 14072 for NEMO/trunk/src/OCE/TRA/traadv_qck.F90
- Timestamp:
- 2020-12-04T08:48:38+01:00 (3 years ago)
- File:
-
- 1 edited
-
NEMO/trunk/src/OCE/TRA/traadv_qck.F90 (modified) (15 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/OCE/TRA/traadv_qck.F90
r13982 r14072 19 19 USE trc_oce ! share passive tracers/Ocean variables 20 20 USE trd_oce ! trends: ocean variables 21 USE trdtra ! trends manager: tracers 21 USE trdtra ! trends manager: tracers 22 22 USE diaptr ! poleward transport diagnostics 23 23 USE iom … … 26 26 USE lib_mpp ! distribued memory computing 27 27 USE lbclnk ! ocean lateral boundary condition (or mpp link) 28 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 28 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 29 29 30 30 IMPLICIT NONE … … 112 112 ! 113 113 ! ! horizontal fluxes are computed with the QUICKEST + ULTIMATE scheme 114 CALL tra_adv_qck_i( kt, cdtype, p2dt, pU, Kbb, Kmm, pt, kjpt, Krhs ) 115 CALL tra_adv_qck_j( kt, cdtype, p2dt, pV, Kbb, Kmm, pt, kjpt, Krhs ) 114 CALL tra_adv_qck_i( kt, cdtype, p2dt, pU, Kbb, Kmm, pt, kjpt, Krhs ) 115 CALL tra_adv_qck_j( kt, cdtype, p2dt, pV, Kbb, Kmm, pt, kjpt, Krhs ) 116 116 117 117 ! ! vertical fluxes are computed with the 2nd order centered scheme … … 142 142 DO jn = 1, kjpt ! tracer loop 143 143 ! ! =========== 144 zfu(:,:,:) = 0._wp ; zfc(:,:,:) = 0._wp 145 zfd(:,:,:) = 0._wp ; zwx(:,:,:) = 0._wp 144 zfu(:,:,:) = 0._wp ; zfc(:,:,:) = 0._wp 145 zfd(:,:,:) = 0._wp ; zwx(:,:,:) = 0._wp 146 146 ! 147 147 !!gm why not using a SHIFT instruction... … … 151 151 END_3D 152 152 IF (nn_hls.EQ.1) CALL lbc_lnk_multi( 'traadv_qck', zfc(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp ) ! Lateral boundary conditions 153 153 154 154 ! 155 155 ! Horizontal advective fluxes 156 156 ! --------------------------- 157 157 DO_3D( 0, 0, nn_hls-1, 0, 1, jpkm1 ) 158 zdir = 0.5 + SIGN( 0.5_wp, pU(ji,jj,jk) ) ! if pU > 0 : zdir = 1 otherwise zdir = 0 159 zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji+1,jj,jk) ! FU in the x-direction for T 158 zdir = 0.5 + SIGN( 0.5_wp, pU(ji,jj,jk) ) ! if pU > 0 : zdir = 1 otherwise zdir = 0 159 zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji+1,jj,jk) ! FU in the x-direction for T 160 160 END_3D 161 161 ! 162 162 DO_3D( 0, 0, nn_hls-1, 0, 1, jpkm1 ) 163 zdir = 0.5 + SIGN( 0.5_wp, pU(ji,jj,jk) ) ! if pU > 0 : zdir = 1 otherwise zdir = 0 163 zdir = 0.5 + SIGN( 0.5_wp, pU(ji,jj,jk) ) ! if pU > 0 : zdir = 1 otherwise zdir = 0 164 164 zdx = ( zdir * e1t(ji,jj) + ( 1. - zdir ) * e1t(ji+1,jj) ) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm) 165 165 zwx(ji,jj,jk) = ABS( pU(ji,jj,jk) ) * p2dt / zdx ! (0<zc_cfl<1 : Courant number on x-direction) … … 167 167 zfd(ji,jj,jk) = zdir * pt(ji+1,jj,jk,jn,Kbb) + ( 1. - zdir ) * pt(ji ,jj,jk,jn,Kbb) ! FD in the x-direction for T 168 168 END_3D 169 !--- Lateral boundary conditions 169 !--- Lateral boundary conditions 170 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 ) 171 171 … … 227 227 DO jn = 1, kjpt ! tracer loop 228 228 ! ! =========== 229 zfu(:,:,:) = 0.0 ; zfc(:,:,:) = 0.0 230 zfd(:,:,:) = 0.0 ; zwy(:,:,:) = 0.0 231 ! 232 DO jk = 1, jpkm1 233 ! 229 zfu(:,:,:) = 0.0 ; zfc(:,:,:) = 0.0 230 zfd(:,:,:) = 0.0 ; zwy(:,:,:) = 0.0 231 ! 232 DO jk = 1, jpkm1 233 ! 234 234 !--- Computation of the ustream and downstream value of the tracer and the mask 235 235 DO_2D( nn_hls-1, nn_hls-1, 0, 0 ) … … 241 241 END DO 242 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 243 243 244 244 ! 245 245 ! Horizontal advective fluxes … … 247 247 ! 248 248 DO_3D( nn_hls-1, 0, 0, 0, 1, jpkm1 ) 249 zdir = 0.5 + SIGN( 0.5_wp, pV(ji,jj,jk) ) ! if pU > 0 : zdir = 1 otherwise zdir = 0 250 zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji,jj+1,jk) ! FU in the x-direction for T 249 zdir = 0.5 + SIGN( 0.5_wp, pV(ji,jj,jk) ) ! if pU > 0 : zdir = 1 otherwise zdir = 0 250 zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji,jj+1,jk) ! FU in the x-direction for T 251 251 END_3D 252 252 ! 253 253 DO_3D( nn_hls-1, 0, 0, 0, 1, jpkm1 ) 254 zdir = 0.5 + SIGN( 0.5_wp, pV(ji,jj,jk) ) ! if pU > 0 : zdir = 1 otherwise zdir = 0 254 zdir = 0.5 + SIGN( 0.5_wp, pV(ji,jj,jk) ) ! if pU > 0 : zdir = 1 otherwise zdir = 0 255 255 zdx = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * e3v(ji,jj,jk,Kmm) 256 256 zwy(ji,jj,jk) = ABS( pV(ji,jj,jk) ) * p2dt / zdx ! (0<zc_cfl<1 : Courant number on x-direction) … … 259 259 END_3D 260 260 261 !--- Lateral boundary conditions 261 !--- Lateral boundary conditions 262 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 ) 263 263 … … 328 328 IF( ln_isfcav ) THEN ! ice-shelf cavities (top of the ocean) 329 329 DO_2D( 0, 0, 0, 0 ) 330 zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kmm) ! linear free surface 330 zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kmm) ! linear free surface 331 331 END_2D 332 332 ELSE ! no ocean cavities (only ocean surface) … … 354 354 !! ** Purpose : Computation of advective flux with Quickest scheme 355 355 !! 356 !! ** Method : 356 !! ** Method : 357 357 !!---------------------------------------------------------------------- 358 358 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in ) :: pfu ! second upwind point … … 361 361 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(inout) :: puc ! input as Courant number ; output as flux 362 362 !! 363 INTEGER :: ji, jj, jk ! dummy loop indices 364 REAL(wp) :: zcoef1, zcoef2, zcoef3 ! local scalars 363 INTEGER :: ji, jj, jk ! dummy loop indices 364 REAL(wp) :: zcoef1, zcoef2, zcoef3 ! local scalars 365 365 REAL(wp) :: zc, zcurv, zfho ! - - 366 366 !---------------------------------------------------------------------- … … 372 372 zcoef2 = 0.5 * zc * ( pfd(ji,jj,jk) - pfc(ji,jj,jk) ) 373 373 zcoef3 = ( 1. - ( zc * zc ) ) * r1_6 * zcurv 374 zfho = zcoef1 - zcoef2 - zcoef3 ! phi_f QUICKEST 374 zfho = zcoef1 - zcoef2 - zcoef3 ! phi_f QUICKEST 375 375 ! 376 376 zcoef1 = pfd(ji,jj,jk) - pfu(ji,jj,jk) … … 378 378 zcoef3 = ABS( zcurv ) 379 379 IF( zcoef3 >= zcoef2 ) THEN 380 zfho = pfc(ji,jj,jk) 380 zfho = pfc(ji,jj,jk) 381 381 ELSE 382 382 zcoef3 = pfu(ji,jj,jk) + ( ( pfc(ji,jj,jk) - pfu(ji,jj,jk) ) / MAX( zc, 1.e-9 ) ) ! phi_REF 383 383 IF( zcoef1 >= 0. ) THEN 384 zfho = MAX( pfc(ji,jj,jk), zfho ) 385 zfho = MIN( zfho, MIN( zcoef3, pfd(ji,jj,jk) ) ) 384 zfho = MAX( pfc(ji,jj,jk), zfho ) 385 zfho = MIN( zfho, MIN( zcoef3, pfd(ji,jj,jk) ) ) 386 386 ELSE 387 zfho = MIN( pfc(ji,jj,jk), zfho ) 388 zfho = MAX( zfho, MAX( zcoef3, pfd(ji,jj,jk) ) ) 387 zfho = MIN( pfc(ji,jj,jk), zfho ) 388 zfho = MAX( zfho, MAX( zcoef3, pfd(ji,jj,jk) ) ) 389 389 ENDIF 390 390 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
