- Timestamp:
- 2016-07-19T10:38:35+02:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/NERC/dev_r5549_BDY_ZEROGRAD/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90
r5147 r6808 20 20 USE trd_oce ! trends: ocean variables 21 21 USE trdtra ! trends manager: tracers 22 USE dynspg_oce ! surface pressure gradient variables23 22 USE diaptr ! poleward transport diagnostics 24 23 ! … … 39 38 40 39 !! * Substitutions 41 # include "domzgr_substitute.h90"42 40 # include "vectopt_loop_substitute.h90" 43 41 !!---------------------------------------------------------------------- … … 79 77 !! prevent the appearance of spurious numerical oscillations 80 78 !! 81 !! ** Action : - update (pta) with the now advective tracer trends 82 !! - save the trends 79 !! ** Action : - update pta with the now advective tracer trends 80 !! - send trends to trdtra module for further diagnostcs (l_trdtra=T) 81 !! - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T) 83 82 !! 84 83 !! ** Reference : Leonard (1979, 1991) … … 88 87 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 89 88 INTEGER , INTENT(in ) :: kjpt ! number of tracers 90 REAL(wp) , DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile oftracer time-step89 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 91 90 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun, pvn, pwn ! 3 ocean velocity components 92 91 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before and now tracer fields … … 102 101 IF(lwp) WRITE(numout,*) 103 102 ENDIF 103 ! 104 104 l_trd = .FALSE. 105 105 IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 106 106 ! 107 ! I. Thehorizontal fluxes are computed with the QUICKEST + ULTIMATE scheme107 ! ! horizontal fluxes are computed with the QUICKEST + ULTIMATE scheme 108 108 CALL tra_adv_qck_i( kt, cdtype, p2dt, pun, ptb, ptn, pta, kjpt ) 109 109 CALL tra_adv_qck_j( kt, cdtype, p2dt, pvn, ptb, ptn, pta, kjpt ) 110 110 111 ! II. Thevertical fluxes are computed with the 2nd order centered scheme111 ! ! vertical fluxes are computed with the 2nd order centered scheme 112 112 CALL tra_adv_cen2_k( kt, cdtype, pwn, ptn, pta, kjpt ) 113 113 ! … … 125 125 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 126 126 INTEGER , INTENT(in ) :: kjpt ! number of tracers 127 REAL(wp) , DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile oftracer time-step127 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 128 128 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun ! i-velocity components 129 129 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before and now tracer fields 130 130 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 131 131 !! 132 INTEGER :: ji, jj, jk, jn ! dummy loop indices133 REAL(wp) :: ztra, zbtr, zdir, zdx, zdt, zmsk ! local scalars134 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwx, zfu, zfc, zfd132 INTEGER :: ji, jj, jk, jn ! dummy loop indices 133 REAL(wp) :: ztra, zbtr, zdir, zdx, zmsk ! local scalars 134 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwx, zfu, zfc, zfd 135 135 !---------------------------------------------------------------------- 136 136 ! … … 139 139 DO jn = 1, kjpt ! tracer loop 140 140 ! ! =========== 141 zfu(:,:,:) = 0.0 ; zfc(:,:,:) = 0.0 142 zfd(:,:,:) = 0.0 ; zwx(:,:,:) = 0.0 143 ! 144 DO jk = 1, jpkm1 145 ! 146 !--- Computation of the ustream and downstream value of the tracer and the mask 141 zfu(:,:,:) = 0._wp ; zfc(:,:,:) = 0._wp 142 zfd(:,:,:) = 0._wp ; zwx(:,:,:) = 0._wp 143 ! 144 !!gm why not using a SHIFT instruction... 145 DO jk = 1, jpkm1 !--- Computation of the ustream and downstream value of the tracer and the mask 147 146 DO jj = 2, jpjm1 148 147 DO ji = fs_2, fs_jpim1 ! vector opt. 149 ! Upstream in the x-direction for the tracer 150 zfc(ji,jj,jk) = ptb(ji-1,jj,jk,jn) 151 ! Downstream in the x-direction for the tracer 152 zfd(ji,jj,jk) = ptb(ji+1,jj,jk,jn) 148 zfc(ji,jj,jk) = ptb(ji-1,jj,jk,jn) ! Upstream in the x-direction for the tracer 149 zfd(ji,jj,jk) = ptb(ji+1,jj,jk,jn) ! Downstream in the x-direction for the tracer 153 150 END DO 154 151 END DO … … 159 156 ! Horizontal advective fluxes 160 157 ! --------------------------- 161 !162 158 DO jk = 1, jpkm1 163 159 DO jj = 2, jpjm1 … … 170 166 ! 171 167 DO jk = 1, jpkm1 172 zdt = p2dt(jk)173 168 DO jj = 2, jpjm1 174 169 DO ji = fs_2, fs_jpim1 ! vector opt. 175 170 zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0 176 zdx = ( zdir * e1t(ji,jj) + ( 1. - zdir ) * e1t(ji+1,jj) ) * e2u(ji,jj) * fse3u(ji,jj,jk)177 zwx(ji,jj,jk) = ABS( pun(ji,jj,jk) ) * zdt / zdx ! (0<zc_cfl<1 : Courant number on x-direction)171 zdx = ( zdir * e1t(ji,jj) + ( 1. - zdir ) * e1t(ji+1,jj) ) * e2u(ji,jj) * e3u_n(ji,jj,jk) 172 zwx(ji,jj,jk) = ABS( pun(ji,jj,jk) ) * p2dt / zdx ! (0<zc_cfl<1 : Courant number on x-direction) 178 173 zfc(ji,jj,jk) = zdir * ptb(ji ,jj,jk,jn) + ( 1. - zdir ) * ptb(ji+1,jj,jk,jn) ! FC in the x-direction for T 179 174 zfd(ji,jj,jk) = zdir * ptb(ji+1,jj,jk,jn) + ( 1. - zdir ) * ptb(ji ,jj,jk,jn) ! FD in the x-direction for T … … 220 215 DO jj = 2, jpjm1 221 216 DO ji = fs_2, fs_jpim1 ! vector opt. 222 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk))217 zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 223 218 ! horizontal advective trends 224 219 ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj,jk) ) … … 228 223 END DO 229 224 END DO 230 ! ! trend diagnostics (contribution of upstream fluxes)225 ! ! trend diagnostics 231 226 IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) ) 232 227 ! … … 246 241 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 247 242 INTEGER , INTENT(in ) :: kjpt ! number of tracers 248 REAL(wp) , DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile oftracer time-step243 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 249 244 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pvn ! j-velocity components 250 245 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before and now tracer fields … … 252 247 !! 253 248 INTEGER :: ji, jj, jk, jn ! dummy loop indices 254 REAL(wp) :: ztra, zbtr, zdir, zdx, z dt, zmsk ! local scalars249 REAL(wp) :: ztra, zbtr, zdir, zdx, zmsk ! local scalars 255 250 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwy, zfu, zfc, zfd 256 251 !---------------------------------------------------------------------- … … 293 288 ! 294 289 DO jk = 1, jpkm1 295 zdt = p2dt(jk)296 290 DO jj = 2, jpjm1 297 291 DO ji = fs_2, fs_jpim1 ! vector opt. 298 292 zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0 299 zdx = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * fse3v(ji,jj,jk)300 zwy(ji,jj,jk) = ABS( pvn(ji,jj,jk) ) * zdt / zdx ! (0<zc_cfl<1 : Courant number on x-direction)293 zdx = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * e3v_n(ji,jj,jk) 294 zwy(ji,jj,jk) = ABS( pvn(ji,jj,jk) ) * p2dt / zdx ! (0<zc_cfl<1 : Courant number on x-direction) 301 295 zfc(ji,jj,jk) = zdir * ptb(ji,jj ,jk,jn) + ( 1. - zdir ) * ptb(ji,jj+1,jk,jn) ! FC in the x-direction for T 302 296 zfd(ji,jj,jk) = zdir * ptb(ji,jj+1,jk,jn) + ( 1. - zdir ) * ptb(ji,jj ,jk,jn) ! FD in the x-direction for T … … 344 338 DO jj = 2, jpjm1 345 339 DO ji = fs_2, fs_jpim1 ! vector opt. 346 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk))340 zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 347 341 ! horizontal advective trends 348 342 ztra = - zbtr * ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) ) … … 352 346 END DO 353 347 END DO 354 ! ! trend diagnostics (contribution of upstream fluxes)348 ! ! trend diagnostics 355 349 IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) ) 356 350 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) … … 380 374 ! 381 375 INTEGER :: ji, jj, jk, jn ! dummy loop indices 382 REAL(wp) :: zbtr , ztra ! local scalars383 376 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwz 384 377 !!---------------------------------------------------------------------- 385 378 ! 386 CALL wrk_alloc( jpi, jpj, jpk, zwz ) 379 CALL wrk_alloc( jpi,jpj,jpk, zwz ) 380 ! 381 zwz(:,:, 1 ) = 0._wp ! surface & bottom values set to zero for all tracers 382 zwz(:,:,jpk) = 0._wp 383 ! 387 384 ! ! =========== 388 385 DO jn = 1, kjpt ! tracer loop 389 386 ! ! =========== 390 ! 1. Bottom value : flux set to zero 391 zwz(:,:,jpk) = 0.e0 ! Bottom value : flux set to zero 392 ! 393 ! ! Surface value 394 IF( lk_vvl ) THEN ; zwz(:,:, 1 ) = 0.e0 ! Variable volume : flux set to zero 395 ELSE ; zwz(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1,jn) ! Constant volume : advective flux through the surface 387 ! 388 DO jk = 2, jpkm1 !* Interior point (w-masked 2nd order centered flux) 389 DO jj = 2, jpjm1 390 DO ji = fs_2, fs_jpim1 ! vector opt. 391 zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk-1,jn) + ptn(ji,jj,jk,jn) ) * wmask(ji,jj,jk) 392 END DO 393 END DO 394 END DO 395 IF( ln_linssh ) THEN !* top value (only in linear free surf. as zwz is multiplied by wmask) 396 IF( ln_isfcav ) THEN ! ice-shelf cavities (top of the ocean) 397 DO jj = 1, jpj 398 DO ji = 1, jpi 399 zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptn(ji,jj,mikt(ji,jj),jn) ! linear free surface 400 END DO 401 END DO 402 ELSE ! no ocean cavities (only ocean surface) 403 zwz(:,:,1) = pwn(:,:,1) * ptn(:,:,1,jn) 404 ENDIF 396 405 ENDIF 397 406 ! 398 DO jk = 2, jpkm1 ! Interior point: second order centered tracer flux at w-point407 DO jk = 1, jpkm1 !== Tracer flux divergence added to the general trend ==! 399 408 DO jj = 2, jpjm1 400 409 DO ji = fs_2, fs_jpim1 ! vector opt. 401 zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk-1,jn) + ptn(ji,jj,jk,jn) ) 402 END DO 403 END DO 404 END DO 405 ! 406 DO jk = 1, jpkm1 !== Tracer flux divergence added to the general trend ==! 407 DO jj = 2, jpjm1 408 DO ji = fs_2, fs_jpim1 ! vector opt. 409 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 410 ! k- vertical advective trends 411 ztra = - zbtr * ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) ) 412 ! added to the general tracer trends 413 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 414 END DO 415 END DO 416 END DO 417 ! ! Save the vertical advective trends for diagnostic 410 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) ) & 411 & * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 412 END DO 413 END DO 414 END DO 415 ! ! Send trends for diagnostic 418 416 IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pwn, ptn(:,:,:,jn) ) 419 417 ! 420 418 END DO 421 419 ! 422 CALL wrk_dealloc( jpi, jpj, jpk,zwz )420 CALL wrk_dealloc( jpi,jpj,jpk, zwz ) 423 421 ! 424 422 END SUBROUTINE tra_adv_cen2_k
Note: See TracChangeset
for help on using the changeset viewer.
