- Timestamp:
- 2010-07-29T12:57:35+02:00 (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traadv_qck.F90
r1559 r2024 2 2 !!============================================================================== 3 3 !! *** MODULE traadv_qck *** 4 !! Ocean activetracers: horizontal & vertical advective trend4 !! Ocean tracers: horizontal & vertical advective trend 5 5 !!============================================================================== 6 6 !! History : 3.0 ! 2008-07 (G. Reffray) Original code 7 !! 3.3 ! 2010-05 (C.Ethe, G. Madec) merge TRC-TRA + switch from velocity to transport 7 8 !!---------------------------------------------------------------------- 8 9 … … 16 17 USE oce ! ocean dynamics and active tracers 17 18 USE dom_oce ! ocean space and time domain 18 USE trdmod ! ocean active tracers trends19 USE trd mod_oce ! ocean variables trends19 USE trdmod_oce ! ocean space and time domain 20 USE trdtra ! ocean tracers trends 20 21 USE trabbl ! advective term in the BBL 21 22 USE lib_mpp ! distribued memory computing … … 24 25 USE in_out_manager ! I/O manager 25 26 USE diaptr ! poleward transport diagnostics 26 USE prtctl ! Print control27 27 28 28 IMPLICIT NONE … … 31 31 PUBLIC tra_adv_qck ! routine called by step.F90 32 32 33 REAL(wp) , DIMENSION(jpi,jpj) :: btr234 REAL(wp) :: r1_633 REAL(wp) :: r1_6 = 1./ 6. 34 LOGICAL :: l_trd ! flag to compute trends 35 35 36 36 !! * Substitutions … … 45 45 CONTAINS 46 46 47 SUBROUTINE tra_adv_qck( kt, pun, pvn, pwn ) 47 SUBROUTINE tra_adv_qck ( kt , cdtype, pun , pvn, pwn, & 48 & ptrab, ptran , ptraa, kjpt ) 48 49 !!---------------------------------------------------------------------- 49 50 !! *** ROUTINE tra_adv_qck *** … … 69 70 !! dt = 2*rdtra and the scalar values are tb and sb 70 71 !! 71 !! On the vertical, the simple centered scheme used tn and sn72 !! On the vertical, the simple centered scheme used ptran 72 73 !! 73 74 !! The fluxes are bounded by the ULTIMATE limiter to … … 75 76 !! prevent the appearance of spurious numerical oscillations 76 77 !! 77 !! ** Action : - update ( ta,sa) with the now advective tracer trends78 !! - save the trends ('key_trdtra')78 !! ** Action : - update (ptraa) with the now advective tracer trends 79 !! - save the trends 79 80 !! 80 81 !! ** Reference : Leonard (1979, 1991) 81 82 !!---------------------------------------------------------------------- 82 USE oce, ONLY : ztrdt => ua ! use ua as workspace 83 USE oce, ONLY : ztrds => va ! use va as workspace 84 !! 85 INTEGER , INTENT(in) :: kt ! ocean time-step index 86 REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pun ! effective ocean velocity, u_component 87 REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pvn ! effective ocean velocity, v_component 88 REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pwn ! effective ocean velocity, w_component 89 !! 90 INTEGER :: ji, jj, jk ! dummy loop indices 91 REAL(wp) :: z_hdivn_x, z_hdivn_y, z_hdivn ! temporary scalars 92 REAL(wp) :: zbtr, z2 ! " " 83 !!* Arguments 84 INTEGER , INTENT(in ) :: kt ! ocean time-step index 85 CHARACTER(len=3), INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 86 INTEGER , INTENT(in ) :: kjpt ! number of tracers 87 REAL(wp) , INTENT(in ), DIMENSION(jpi,jpj,jpk) :: pun, pvn, pwn ! 3 ocean velocity components 88 REAL(wp) , INTENT(in ), DIMENSION(jpi,jpj,jpk,kjpt) :: ptrab, ptran ! before and now tracer fields 89 REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: ptraa ! tracer trend 90 !!* Local declarations 91 REAL(wp) :: z2 ! temporary scalar 93 92 !!---------------------------------------------------------------------- 94 93 … … 98 97 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 99 98 IF(lwp) WRITE(numout,*) 100 btr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) ) 101 r1_6 = 1. / 6. 99 ! 100 l_trd = .FALSE. 101 IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 102 102 ENDIF 103 103 … … 109 109 !--------------------------------------------------------------------------- 110 110 111 CALL tra_adv_qck_i( pun, tb, tn, ta, ztrdt, z2) 112 CALL tra_adv_qck_i( pun, sb, sn, sa, ztrds, z2) 113 114 IF( l_trdtra ) CALL trd_mod(ztrdt, ztrds, jptra_trd_xad, 'TRA', kt) 115 116 CALL tra_adv_qck_j( kt, pvn, tb, tn, ta, ztrdt, pht_adv, z2) 117 CALL tra_adv_qck_j( kt, pvn, sb, sn, sa, ztrds, pst_adv, z2) 118 119 IF( l_trdtra ) THEN 120 CALL trd_mod(ztrdt, ztrds, jptra_trd_yad, 'TRA', kt) 121 ! 122 ztrdt(:,:,:) = ta(:,:,:) ! Save the horizontal up-to-date ta/sa trends 123 ztrds(:,:,:) = sa(:,:,:) 124 END IF 125 126 IF(ln_ctl) CALL prt_ctl( tab3d_1=ta, clinfo1=' qck had - Ta: ', mask1=tmask, & 127 & tab3d_2=sa, clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' ) 111 CALL tra_adv_qck_i( kt , cdtype, pun , z2, & 112 & ptrab, ptran , ptraa, kjpt ) 113 CALL tra_adv_qck_j( kt , cdtype, pvn , z2, & 114 & ptrab, ptran , ptraa, kjpt ) 128 115 129 116 ! II. The vertical fluxes are computed with the 2nd order centered scheme 130 117 !------------------------------------------------------------------------- 131 118 ! 132 CALL tra_adv_cen2_k( pwn, tn, ta ) 133 CALL tra_adv_cen2_k( pwn, sn, sa ) 134 ! 135 !Save the vertical advective trends for diagnostic 136 ! ---------------------------------------------------- 137 IF( l_trdtra ) THEN 138 ! Recompute the vertical advection zta & zsa trends computed 139 ! at the step 2. above in making the difference between the new 140 ! trends and the previous one: ta()/sa - ztrdt()/ztrds() and substract 141 ! the term tn()/sn()*hdivn() to recover the W gradz(T/S) trends 142 143 DO jk = 1, jpkm1 119 CALL tra_adv_cen2_k( kt , cdtype, pwn, & 120 & ptran, ptraa , kjpt ) 121 ! 122 END SUBROUTINE tra_adv_qck 123 124 SUBROUTINE tra_adv_qck_i( kt , cdtype, pun , pz2, & 125 & ptrab, ptran , ptraa, kjpt ) 126 !!---------------------------------------------------------------------- 127 !! 128 !!---------------------------------------------------------------------- 129 !!* Module used 130 USE oce , zwx => ua ! use ua as workspace 131 !!* Arguments 132 INTEGER , INTENT(in ) :: kt ! ocean time-step index 133 CHARACTER(len=3), INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 134 INTEGER , INTENT(in ) :: kjpt ! number of tracers 135 REAL(wp) , INTENT(in ), DIMENSION(jpi,jpj,jpk) :: pun ! zonal velocity component 136 REAL(wp) , INTENT(in ), DIMENSION(jpi,jpj,jpk,kjpt) :: ptrab, ptran ! before tracer fields 137 REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: ptraa ! tracer trend 138 REAL(wp) , INTENT(in ) :: pz2 139 !!* Local declarations 140 INTEGER :: ji, jj, jk, jn ! dummy loop indices 141 REAL(wp) :: ztra, zbtr ! temporary scalars 142 REAL(wp) :: zdir, zdx, zdt, zmsk ! temporary scalars 143 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zfu, zfc, zfd 144 !---------------------------------------------------------------------- 145 146 147 DO jn = 1, kjpt ! tracer loop 148 ! ! =========== 149 zfu(:,:,:) = 0.0 ; zfc(:,:,:) = 0.0 150 zfd(:,:,:) = 0.0 ; zwx(:,:,:) = 0.0 151 ! 152 DO jk = 1, jpkm1 153 ! 154 !--- Computation of the ustream and downstream value of the tracer and the mask 144 155 DO jj = 2, jpjm1 145 156 DO ji = fs_2, fs_jpim1 ! vector opt. 146 #if defined key_zco 147 zbtr = btr2(ji,jj) 148 z_hdivn_x = e2u(ji,jj)*pun(ji,jj,jk) - e2u(ji-1,jj)*pun(ji-1,jj,jk) 149 z_hdivn_y = e1v(ji,jj)*pvn(ji,jj,jk) - e1v(ji,jj-1)*pvn(ji,jj-1,jk) 150 #else 151 zbtr = btr2(ji,jj) / fse3t(ji,jj,jk) 152 z_hdivn_x = e2u(ji,jj)*fse3u(ji,jj,jk)*pun(ji,jj,jk) - e2u(ji-1,jj)*fse3u(ji-1,jj,jk)*pun(ji-1,jj,jk) 153 z_hdivn_y = e1v(ji,jj)*fse3v(ji,jj,jk)*pvn(ji,jj,jk) - e1v(ji,jj-1)*fse3v(ji,jj-1,jk)*pvn(ji,jj-1,jk) 154 #endif 155 z_hdivn = (z_hdivn_x + z_hdivn_y) * zbtr 156 ztrdt(ji,jj,jk) = ta(ji,jj,jk) - ztrdt(ji,jj,jk) - tn(ji,jj,jk) * z_hdivn 157 ztrds(ji,jj,jk) = sa(ji,jj,jk) - ztrds(ji,jj,jk) - sn(ji,jj,jk) * z_hdivn 158 END DO 159 END DO 160 END DO 161 CALL trd_mod(ztrdt, ztrds, jptra_trd_zad, 'TRA', kt) 162 ENDIF 163 164 IF(ln_ctl) CALL prt_ctl( tab3d_1=ta, clinfo1=' qck zad - Ta: ', mask1=tmask, & 165 & tab3d_2=sa, clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' ) 166 ! 167 END SUBROUTINE tra_adv_qck 168 169 170 SUBROUTINE tra_adv_qck_i ( pun, tra, tran, traa, ztrdtra, z2 ) 171 !!---------------------------------------------------------------------- 172 !! 173 !!---------------------------------------------------------------------- 174 REAL, INTENT(in) :: z2 175 REAL(wp), INTENT(in) , DIMENSION(jpi,jpj,jpk) :: pun, tra, tran ! horizontal effective velocity 176 REAL(wp), INTENT(out) , DIMENSION(jpi,jpj,jpk) :: ztrdtra 177 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: traa 178 ! 179 INTEGER :: ji, jj, jk 180 REAL(wp) :: za, zbtr, dir, dx, dt ! temporary scalars 181 REAL(wp) :: z_hdivn_x 182 REAL(wp), DIMENSION(jpi,jpj) :: zmask, zupst, zdwst, zc_cfl 183 REAL(wp), DIMENSION(jpi,jpj) :: zfu, zfc, zfd, zfho, zmskl, zsc_e 184 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zflux 185 !---------------------------------------------------------------------- 186 187 zfu (:,jpj) = 0.e0 ; zfc (:,jpj) = 0.e0 188 zfd (:,jpj) = 0.e0 ; zc_cfl(:,jpj) = 0.e0 189 zsc_e (:,jpj) = 0.e0 ; zmskl (:,jpj) = 0.e0 190 zfho (:,jpj) = 0.e0 191 ! =============== 192 DO jk = 1, jpkm1 ! Horizontal slab 193 ! ! =============== 194 !--- Computation of the ustream and downstream value of the tracer and the mask 195 DO jj = 2, jpjm1 196 DO ji = 2, fs_jpim1 ! vector opt. 197 ! Upstream in the x-direction for the tracer 198 zupst(ji,jj)=tra(ji-1,jj,jk) 199 ! Downstream in the x-direction for the tracer 200 zdwst(ji,jj)=tra(ji+1,jj,jk) 201 ! Mask at the T-points in the x-direction (mask=0 or mask=1) 202 zmask(ji,jj)=tmask(ji-1,jj,jk)+tmask(ji,jj,jk)+tmask(ji+1,jj,jk)-2 203 END DO 204 END DO 205 ! 206 !--- Lateral boundary conditions 207 CALL lbc_lnk( zupst(:,:), 'T', 1. ) 208 CALL lbc_lnk( zdwst(:,:), 'T', 1. ) 209 CALL lbc_lnk( zmask(:,:), 'T', 1. ) 157 ! Upstream in the x-direction for the tracer 158 zfc(ji,jj,jk) = ptrab(ji-1,jj,jk,jn) 159 ! Downstream in the x-direction for the tracer 160 zfd(ji,jj,jk) = ptrab(ji+1,jj,jk,jn) 161 END DO 162 END DO 163 END DO 164 ! 165 !--- Lateral boundary conditions 166 CALL lbc_lnk( zfc(:,:,:), 'T', 1. ) ; CALL lbc_lnk( zfd(:,:,:), 'T', 1. ) 167 210 168 ! 211 169 ! Horizontal advective fluxes 212 170 ! --------------------------- 213 171 ! 214 dt = z2 * rdttra(jk) 215 !--- tracer flux at u-points 216 DO jj = 1, jpjm1 217 DO ji = 1, jpi 218 #if defined key_zco 219 zsc_e(ji,jj) = e2u(ji,jj) 220 #else 221 zsc_e(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) 222 #endif 223 dir = 0.5 + sign(0.5,pun(ji,jj,jk)) ! if pun>0 : dir = 1 otherwise dir = 0 224 dx = dir * e1t(ji,jj) + (1-dir)* e1t(ji+1,jj) 225 zc_cfl (ji,jj) = ABS(pun(ji,jj,jk))*dt/dx ! (0<zc_cfl<1 : Courant number on x-direction) 226 227 zfu(ji,jj) = dir*zupst(ji ,jj )+(1-dir)*zdwst(ji+1,jj ) ! FU in the x-direction for T 228 zfc(ji,jj) = dir*tra (ji ,jj,jk)+(1-dir)*tra (ji+1,jj,jk) ! FC in the x-direction for T 229 zfd(ji,jj) = dir*tra (ji+1,jj,jk)+(1-dir)*tra (ji ,jj,jk) ! FD in the x-direction for T 230 zmskl(ji,jj) = dir*zmask(ji ,jj) +(1-dir)*zmask(ji+1,jj) 231 END DO 232 END DO 233 ! 172 DO jk = 1, jpkm1 173 DO jj = 2, jpjm1 174 DO ji = fs_2, fs_jpim1 ! vector opt. 175 zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0 176 zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji+1,jj,jk) ! FU in the x-direction for T 177 END DO 178 END DO 179 END DO 180 ! 181 DO jk = 1, jpkm1 182 zdt = pz2 * rdttra(jk) 183 DO jj = 2, jpjm1 184 DO ji = fs_2, fs_jpim1 ! vector opt. 185 zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0 186 zdx = ( zdir * e1t(ji,jj) + ( 1. - zdir ) * e1t(ji+1,jj) ) * e2u(ji,jj) * fse3u(ji,jj,jk) 187 zwx(ji,jj,jk) = ABS( pun(ji,jj,jk) ) * zdt / zdx ! (0<zc_cfl<1 : Courant number on x-direction) 188 zfc(ji,jj,jk) = zdir * ptrab(ji ,jj,jk,jn) + ( 1. - zdir ) * ptrab(ji+1,jj,jk,jn) ! FC in the x-direction for T 189 zfd(ji,jj,jk) = zdir * ptrab(ji+1,jj,jk,jn) + ( 1. - zdir ) * ptrab(ji ,jj,jk,jn) ! FD in the x-direction for T 190 END DO 191 END DO 192 END DO ! 193 194 !--- Lateral boundary conditions 195 CALL lbc_lnk( zfu(:,:,:), 'T', 1. ) ; CALL lbc_lnk( zfd(:,:,:), 'T', 1. ) 196 CALL lbc_lnk( zfc(:,:,:), 'T', 1. ) ; CALL lbc_lnk( zwx(:,:,:), 'T', 1. ) 197 234 198 !--- QUICKEST scheme 199 CALL quickest( zfu, zfd, zfc, zwx ) 200 ! 201 ! Mask at the T-points in the x-direction (mask=0 or mask=1) 202 DO jk = 1, jpkm1 203 DO jj = 2, jpjm1 204 DO ji = fs_2, fs_jpim1 ! vector opt. 205 zfu(ji,jj,jk) = tmask(ji-1,jj,jk) + tmask(ji,jj,jk) + tmask(ji+1,jj,jk) - 2. 206 ENDDO 207 END DO 208 END DO 209 !--- Lateral boundary conditions 210 CALL lbc_lnk( zfu(:,:,:), 'T', 1. ) 211 ! 235 212 ! Tracer flux on the x-direction 236 CALL quickest(zfu,zfd,zfc,zfho,zc_cfl) 237 !--- If the second ustream point is a land point 238 !--- the flux is computed by the 1st order UPWIND scheme 239 zfho(:,:) = zmskl(:,:)*zfho(:,:) + (1.-zmskl(:,:))*zfc(:,:) 240 !--- Computation of fluxes 241 zflux(:,:,jk) = zsc_e(:,:)*pun(:,:,jk)*zfho(:,:) 242 ! 243 !--- Tracer flux divergence at t-point added to the general trend 244 DO jj = 2, jpjm1 245 DO ji = fs_2, fs_jpim1 ! vector opt. 246 !--- horizontal advective trends 247 #if defined key_zco 248 zbtr = btr2(ji,jj) 249 #else 250 zbtr = btr2(ji,jj) / fse3t(ji,jj,jk) 251 #endif 252 za = - zbtr * ( zflux(ji,jj,jk) - zflux(ji-1,jj,jk) ) 253 !--- add it to the general tracer trends 254 traa(ji,jj,jk) = traa(ji,jj,jk) + za 255 END DO 256 END DO 257 ! ! =============== 258 END DO ! End of slab 259 ! ! =============== 260 ! 261 ! Save the horizontal advective trends for diagnostic 262 ! ----------------------------------------------------- 263 IF( l_trdtra ) THEN 264 ! T/S ZONAL advection trends 265 ztrdtra(:,:,:) = 0.e0 266 ! 267 DO jk = 1, jpkm1 213 DO jk = 1, jpkm1 214 ! 215 DO jj = 2, jpjm1 216 DO ji = fs_2, fs_jpim1 ! vector opt. 217 zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0 218 !--- If the second ustream point is a land point 219 !--- the flux is computed by the 1st order UPWIND scheme 220 zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji+1,jj,jk) 221 zwx(ji,jj,jk) = zmsk * zwx(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk) 222 zwx(ji,jj,jk) = zwx(ji,jj,jk) * pun(ji,jj,jk) 223 END DO 224 END DO 225 ! 226 ! Computation of the trend 227 DO jj = 2, jpjm1 228 DO ji = fs_2, fs_jpim1 ! vector opt. 229 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 230 ! horizontal advective trends 231 ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj,jk) ) 232 !--- add it to the general tracer trends 233 ptraa(ji,jj,jk,jn) = ptraa(ji,jj,jk,jn) + ztra 234 END DO 235 END DO 236 ! 237 END DO 238 ! ! trend diagnostics (contribution of upstream fluxes) 239 IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_trd_xad, zwx, pun, ptran(:,:,:,jn) ) 240 ! 241 END DO 242 ! 243 END SUBROUTINE tra_adv_qck_i 244 245 SUBROUTINE tra_adv_qck_j( kt , cdtype, pvn , pz2, & 246 & ptrab, ptran , ptraa, kjpt ) 247 !!---------------------------------------------------------------------- 248 !! 249 !!---------------------------------------------------------------------- 250 !!* Module used 251 USE oce , zwy => ua ! use ua as workspace 252 !!* Arguments 253 INTEGER , INTENT(in ) :: kt ! ocean time-step index 254 CHARACTER(len=3), INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 255 INTEGER , INTENT(in ) :: kjpt ! number of tracers 256 REAL(wp) , INTENT(in ), DIMENSION(jpi,jpj,jpk) :: pvn ! meridional velocity component 257 REAL(wp) , INTENT(in ), DIMENSION(jpi,jpj,jpk,kjpt) :: ptrab, ptran ! before tracer fields 258 REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: ptraa ! tracer trend 259 REAL(wp) , INTENT(in ) :: pz2 260 !!* Local declarations 261 INTEGER :: ji, jj, jk, jn ! dummy loop indices 262 REAL(wp) :: ztra, zbtr ! temporary scalars 263 REAL(wp) :: zdir, zdx, zdt, zmsk ! temporary scalars 264 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zfu, zfc, zfd 265 !---------------------------------------------------------------------- 266 267 DO jn = 1, kjpt ! tracer loop 268 ! ! =========== 269 zfu(:,:,:) = 0.0 ; zfc(:,:,:) = 0.0 270 zfd(:,:,:) = 0.0 ; zwy(:,:,:) = 0.0 271 ! 272 DO jk = 1, jpkm1 273 ! 274 !--- Computation of the ustream and downstream value of the tracer and the mask 268 275 DO jj = 2, jpjm1 269 276 DO ji = fs_2, fs_jpim1 ! vector opt. 270 !-- Compute zonal divergence by splitting hdivn (see divcur.F90) 271 ! N.B. This computation is not valid along OBCs (if any) 272 #if defined key_zco 273 zbtr = btr2(ji,jj) 274 z_hdivn_x = ( e2u(ji ,jj) * pun(ji ,jj,jk) & 275 & - e2u(ji-1,jj) * pun(ji-1,jj,jk) ) * zbtr 276 #else 277 zbtr = btr2(ji,jj) / fse3t(ji,jj,jk) 278 z_hdivn_x = ( e2u(ji ,jj) * fse3u(ji ,jj,jk) * pun(ji ,jj,jk) & 279 & - e2u(ji-1,jj) * fse3u(ji-1,jj,jk) * pun(ji-1,jj,jk) ) * zbtr 280 #endif 281 ztrdtra(ji,jj,jk) = - zbtr * ( zflux(ji,jj,jk) - zflux(ji-1,jj,jk) ) + tran(ji,jj,jk) * z_hdivn_x 282 END DO 283 END DO 284 END DO 285 END IF 286 287 END SUBROUTINE tra_adv_qck_i 288 289 290 SUBROUTINE tra_adv_qck_j ( kt, pvn, tra, tran, traa, ztrdtra, trd_adv, z2 ) 291 !!---------------------------------------------------------------------- 292 !! 293 !!---------------------------------------------------------------------- 294 INTEGER, INTENT(in) :: kt ! ocean time-step index 295 REAL, INTENT(in) :: z2 296 REAL(wp), INTENT(in) , DIMENSION(jpi,jpj,jpk) :: pvn, tra, tran ! horizontal effective velocity 297 REAL(wp), INTENT(out) , DIMENSION(jpj) :: trd_adv 298 REAL(wp), INTENT(out) , DIMENSION(jpi,jpj,jpk) :: ztrdtra 299 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: traa 300 !! 301 INTEGER :: ji, jj, jk 302 REAL(wp) :: za, zbtr, dir, dx, dt ! temporary scalars 303 REAL(wp) :: z_hdivn_y 304 REAL(wp), DIMENSION(jpi,jpj) :: zmask, zupst, zdwst, zc_cfl 305 REAL(wp), DIMENSION(jpi,jpj) :: zfu, zfc, zfd, zfho, zmskl, zsc_e 306 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zflux 307 !---------------------------------------------------------------------- 308 ! II. Part 2 : y-direction 309 !---------------------------------------------------------------------- 310 311 zfu (:,jpj) = 0.e0 ; zfc (:,jpj) = 0.e0 312 zfd (:,jpj) = 0.e0 ; zc_cfl(:,jpj) = 0.e0 313 zsc_e (:,jpj) = 0.e0 ; zmskl (:,jpj) = 0.e0 314 zfho (:,jpj) = 0.e0 315 316 ! =============== 317 DO jk = 1, jpkm1 ! Horizontal slab 318 ! ! =============== 319 !--- Computation of the ustream and downstream value of the tracer and the mask 320 DO jj = 2, jpjm1 321 DO ji = 2, fs_jpim1 ! vector opt. 322 ! Upstream in the x-direction for the tracer 323 zupst(ji,jj)=tra(ji,jj-1,jk) 324 ! Downstream in the x-direction for the tracer 325 zdwst(ji,jj)=tra(ji,jj+1,jk) 326 ! Mask at the T-points in the x-direction (mask=0 or mask=1) 327 zmask(ji,jj)=tmask(ji,jj-1,jk)+tmask(ji,jj,jk)+tmask(ji,jj+1,jk)-2 328 END DO 329 END DO 330 ! 331 !--- Lateral boundary conditions 332 CALL lbc_lnk( zupst(:,:), 'T', 1. ) 333 CALL lbc_lnk( zdwst(:,:), 'T', 1. ) 334 CALL lbc_lnk( zmask(:,:), 'T', 1. ) 277 ! Upstream in the x-direction for the tracer 278 zfc(ji,jj,jk) = ptrab(ji,jj-1,jk,jn) 279 ! Downstream in the x-direction for the tracer 280 zfd(ji,jj,jk) = ptrab(ji,jj+1,jk,jn) 281 END DO 282 END DO 283 END DO 284 ! 285 !--- Lateral boundary conditions 286 CALL lbc_lnk( zfc(:,:,:), 'T', 1. ) ; CALL lbc_lnk( zfd(:,:,:), 'T', 1. ) 287 335 288 ! 336 289 ! Horizontal advective fluxes 337 290 ! --------------------------- 338 291 ! 339 dt = z2 * rdttra(jk) 340 !--- tracer flux at v-points 341 DO jj = 1, jpjm1 292 DO jk = 1, jpkm1 293 DO jj = 2, jpjm1 294 DO ji = fs_2, fs_jpim1 ! vector opt. 295 zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0 296 zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji,jj+1,jk) ! FU in the x-direction for T 297 END DO 298 END DO 299 END DO 300 ! 301 DO jk = 1, jpkm1 302 zdt = pz2 * rdttra(jk) 303 DO jj = 2, jpjm1 304 DO ji = fs_2, fs_jpim1 ! vector opt. 305 zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0 306 zdx = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * fse3v(ji,jj,jk) 307 zwy(ji,jj,jk) = ABS( pvn(ji,jj,jk) ) * zdt / zdx ! (0<zc_cfl<1 : Courant number on x-direction) 308 zfc(ji,jj,jk) = zdir * ptrab(ji,jj ,jk,jn) + ( 1. - zdir ) * ptrab(ji,jj+1,jk,jn) ! FC in the x-direction for T 309 zfd(ji,jj,jk) = zdir * ptrab(ji,jj+1,jk,jn) + ( 1. - zdir ) * ptrab(ji,jj ,jk,jn) ! FD in the x-direction for T 310 END DO 311 END DO 312 END DO ! 313 314 !--- Lateral boundary conditions 315 CALL lbc_lnk( zfu(:,:,:), 'T', 1. ) ; CALL lbc_lnk( zfd(:,:,:), 'T', 1. ) 316 CALL lbc_lnk( zfc(:,:,:), 'T', 1. ) ; CALL lbc_lnk( zwy(:,:,:), 'T', 1. ) 317 318 !--- QUICKEST scheme 319 CALL quickest( zfu, zfd, zfc, zwy ) 320 ! 321 ! Mask at the T-points in the x-direction (mask=0 or mask=1) 322 DO jk = 1, jpkm1 323 DO jj = 2, jpjm1 324 DO ji = fs_2, fs_jpim1 ! vector opt. 325 zfu(ji,jj,jk) = tmask(ji,jj-1,jk) + tmask(ji,jj,jk) + tmask(ji,jj+1,jk) - 2. 326 ENDDO 327 END DO 328 END DO 329 !--- Lateral boundary conditions 330 CALL lbc_lnk( zfu(:,:,:), 'T', 1. ) 331 ! 332 ! Tracer flux on the x-direction 333 DO jk = 1, jpkm1 334 ! 335 DO jj = 2, jpjm1 336 DO ji = fs_2, fs_jpim1 ! vector opt. 337 zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0 338 !--- If the second ustream point is a land point 339 !--- the flux is computed by the 1st order UPWIND scheme 340 zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji,jj+1,jk) 341 zwy(ji,jj,jk) = zmsk * zwy(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk) 342 zwy(ji,jj,jk) = zwy(ji,jj,jk) * pvn(ji,jj,jk) 343 END DO 344 END DO 345 ! 346 ! Computation of the trend 347 DO jj = 2, jpjm1 348 DO ji = fs_2, fs_jpim1 ! vector opt. 349 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 350 ! horizontal advective trends 351 ztra = - zbtr * ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) ) 352 !--- add it to the general tracer trends 353 ptraa(ji,jj,jk,jn) = ptraa(ji,jj,jk,jn) + ztra 354 END DO 355 END DO 356 ! 357 END DO 358 ! ! trend diagnostics (contribution of upstream fluxes) 359 IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_trd_yad, zwy, pvn, ptran(:,:,:,jn) ) 360 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) 361 IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN 362 IF( jn == jp_tem ) pht_adv(:) = ptr_vj( zwy(:,:,:) ) 363 IF( jn == jp_sal ) pst_adv(:) = ptr_vj( zwy(:,:,:) ) 364 ENDIF 365 ! 366 END DO 367 368 END SUBROUTINE tra_adv_qck_j 369 370 SUBROUTINE tra_adv_cen2_k( kt , cdtype, pwn, & 371 & ptran, ptraa , kjpt ) 372 !!---------------------------------------------------------------------- 373 !! 374 !!---------------------------------------------------------------------- 375 !!* Module used 376 USE oce , zwz => ua ! use ua as workspace 377 !!* Arguments 378 INTEGER , INTENT(in ) :: kt ! ocean time-step index 379 CHARACTER(len=3), INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 380 INTEGER , INTENT(in ) :: kjpt ! number of tracers 381 REAL(wp) , INTENT(in ), DIMENSION(jpi,jpj,jpk) :: pwn ! vertical velocity component 382 REAL(wp) , INTENT(in ), DIMENSION(jpi,jpj,jpk,kjpt) :: ptran ! now tracer field 383 REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: ptraa ! tracer trend 384 !!* Local declarations 385 INTEGER :: ji, jj, jk, jn ! dummy loop indices 386 REAL(wp) :: zbtr , ztra ! temporary scalars 387 !!---------------------------------------------------------------------- 388 389 ! 390 DO jn = 1, kjpt ! tracer loop 391 ! ! =========== 392 ! 1. Bottom value : flux set to zero 393 zwz(:,:,jpk) = 0.e0 ! Bottom value : flux set to zero 394 ! 395 ! ! Surface value 396 IF( lk_vvl ) THEN ; zwz(:,:, 1 ) = 0.e0 ! Variable volume : flux set to zero 397 ELSE ; zwz(:,:, 1 ) = pwn(:,:,1) * ptran(:,:,1,jn) ! Constant volume : advective flux through the surface 398 ENDIF 399 ! 400 DO jk = 2, jpkm1 ! Interior point: second order centered tracer flux at w-point 401 DO jj = 2, jpjm1 402 DO ji = fs_2, fs_jpim1 ! vector opt. 403 zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptran(ji,jj,jk-1,jn) + ptran(ji,jj,jk,jn) ) 404 END DO 405 END DO 406 END DO 407 ! 408 DO jk = 1, jpkm1 !== Tracer flux divergence added to the general trend ==! 409 DO jj = 2, jpjm1 410 DO ji = fs_2, fs_jpim1 ! vector opt. 411 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 412 ! k- vertical advective trends 413 ztra = - zbtr * ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) ) 414 ! added to the general tracer trends 415 ptraa(ji,jj,jk,jn) = ptraa(ji,jj,jk,jn) + ztra 416 END DO 417 END DO 418 END DO 419 ! ! Save the vertical advective trends for diagnostic 420 IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_trd_zad, zwz, pwn, ptran(:,:,:,jn) ) 421 ! 422 END DO 423 ! 424 END SUBROUTINE tra_adv_cen2_k 425 426 427 SUBROUTINE quickest( pfu, pfd, pfc, puc ) 428 !!---------------------------------------------------------------------- 429 !! 430 !! ** Purpose : Computation of advective flux with Quickest scheme 431 !! 432 !! ** Method : 433 !!---------------------------------------------------------------------- 434 REAL(wp), INTENT(in) , DIMENSION(jpi,jpj,jpk) :: pfu ! second upwind point 435 REAL(wp), INTENT(in) , DIMENSION(jpi,jpj,jpk) :: pfd ! first douwning point 436 REAL(wp), INTENT(in) , DIMENSION(jpi,jpj,jpk) :: pfc ! the central point (or the first upwind point) 437 REAL(wp), INTENT(inout) , DIMENSION(jpi,jpj,jpk) :: puc ! input as Courant number ; output as flux 438 !! 439 INTEGER :: ji, jj, jk ! dummy loop indices 440 REAL(wp) :: zcoef1, zcoef2, zcoef3 ! temporary scalars 441 REAL(wp) :: zc, zcurv, zfho ! 442 !---------------------------------------------------------------------- 443 444 DO jk = 1, jpkm1 445 DO jj = 1, jpj 342 446 DO ji = 1, jpi 343 #if defined key_zco 344 zsc_e(ji,jj) = e1v(ji,jj) 345 #else 346 zsc_e(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) 347 #endif 348 dir = 0.5 + sign(0.5,pvn(ji,jj,jk)) ! if pvn>0 : dir = 1 otherwise dir = 0 349 dx = dir * e2t(ji,jj) + (1-dir)* e2t(ji,jj+1) 350 zc_cfl(ji,jj) = ABS(pvn(ji,jj,jk))*dt/dx ! (0<zc_cfl<1 : Courant number on y-direction) 351 352 zfu(ji,jj) = dir*zupst(ji,jj )+(1-dir)*zdwst(ji,jj+1 ) ! FU in the y-direction for T 353 zfc(ji,jj) = dir*tra (ji,jj ,jk)+(1-dir)*tra (ji,jj+1,jk) ! FC in the y-direction for T 354 zfd(ji,jj) = dir*tra (ji,jj+1,jk)+(1-dir)*tra (ji,jj ,jk) ! FD in the y-direction for T 355 zmskl(ji,jj) = dir*zmask(ji,jj )+(1-dir)*zmask(ji,jj+1) 356 END DO 357 END DO 358 ! 359 !--- QUICKEST scheme 360 ! Tracer flux on the y-direction 361 CALL quickest(zfu,zfd,zfc,zfho,zc_cfl) 362 !--- If the second ustream point is a land point 363 !--- the flux is computed by the 1st order UPWIND scheme 364 zfho(:,:) = zmskl(:,:)*zfho(:,:) + (1.-zmskl(:,:))*zfc(:,:) 365 !--- Computation of fluxes 366 zflux(:,:,jk) = zsc_e(:,:)*pvn(:,:,jk)*zfho(:,:) 367 ! 368 !--- Tracer flux divergence at t-point added to the general trend 369 DO jj = 2, jpjm1 370 DO ji = fs_2, fs_jpim1 ! vector opt. 371 !--- horizontal advective trends 372 #if defined key_zco 373 zbtr = btr2(ji,jj) 374 #else 375 zbtr = btr2(ji,jj) / fse3t(ji,jj,jk) 376 #endif 377 za = - zbtr * ( zflux(ji,jj,jk) - zflux(ji,jj-1,jk) ) 378 !--- add it to the general tracer trends 379 traa(ji,jj,jk) = traa(ji,jj,jk) + za 380 END DO 381 END DO 382 ! ! =============== 383 END DO ! End of slab 384 ! ! =============== 385 ! 386 ! Save the horizontal advective trends for diagnostic 387 ! ----------------------------------------------------- 388 IF( l_trdtra ) THEN 389 ! T/S MERIDIONAL advection trends 390 DO jk = 1, jpkm1 391 DO jj = 2, jpjm1 392 DO ji = fs_2, fs_jpim1 ! vector opt. 393 !-- Compute merid. divergence by splitting hdivn (see divcur.F90) 394 ! N.B. This computation is not valid along OBCs (if any) 395 #if defined key_zco 396 zbtr = btr2(ji,jj) 397 z_hdivn_y = ( e1v(ji,jj ) * pvn(ji,jj ,jk) & 398 & - e1v(ji,jj-1) * pvn(ji,jj-1,jk) ) * zbtr 399 #else 400 zbtr = btr2(ji,jj) / fse3t(ji,jj,jk) 401 z_hdivn_y = ( e1v(ji, jj) * fse3v(ji,jj ,jk) * pvn(ji,jj ,jk) & 402 & - e1v(ji,jj-1) * fse3v(ji,jj-1,jk) * pvn(ji,jj-1,jk) ) * zbtr 403 #endif 404 ztrdtra(ji,jj,jk) = - zbtr * ( zflux(ji,jj,jk) - zflux(ji,jj-1,jk) ) + tran(ji,jj,jk) * z_hdivn_y 405 END DO 406 END DO 407 END DO 408 END IF 409 410 ! "zonal" mean advective heat and salt transport 411 ! ---------------------------------------------- 412 413 IF( ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN 414 IF( lk_zco ) THEN 415 DO jk = 1, jpkm1 416 DO jj = 2, jpjm1 417 DO ji = fs_2, fs_jpim1 ! vector opt. 418 zflux(ji,jj,jk) = zflux(ji,jj,jk) * fse3v(ji,jj,jk) 419 END DO 420 END DO 421 END DO 422 ENDIF 423 trd_adv(:) = ptr_vj( zflux(:,:,:) ) 424 ENDIF 425 426 END SUBROUTINE tra_adv_qck_j 427 428 429 SUBROUTINE tra_adv_cen2_k ( pwn, ptn, pta ) 430 !!---------------------------------------------------------------------- 431 !! 432 !!---------------------------------------------------------------------- 433 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj,jpk) :: pwn ! vertical effective velocity 434 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj,jpk) :: ptn ! now tracer 435 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pta ! tracer general trend 436 !! 437 INTEGER :: ji, jj, jk ! dummy loop indices 438 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zflux ! 3D workspace 439 !!---------------------------------------------------------------------- 440 ! 441 ! !== Vertical advective fluxes ==! 442 zflux(:,:,jpk) = 0.e0 ! Bottom value : flux set to zero 443 ! 444 ! ! Surface value 445 IF( lk_vvl ) THEN ; zflux(:,:, 1 ) = 0.e0 ! Variable volume : flux set to zero 446 ELSE ; zflux(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1) ! Constant volume : advective flux through the surface 447 ENDIF 448 ! 449 DO jk = 2, jpkm1 ! Interior point: second order centered tracer flux at w-point 450 DO jj = 2, jpjm1 451 DO ji = fs_2, fs_jpim1 ! vector opt. 452 zflux(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk-1) + ptn(ji,jj,jk) ) 453 END DO 454 END DO 455 END DO 456 ! 457 DO jk = 1, jpkm1 !== Tracer flux divergence added to the general trend ==! 458 DO jj = 2, jpjm1 459 DO ji = fs_2, fs_jpim1 ! vector opt. 460 pta(ji,jj,jk) = pta(ji,jj,jk) - ( zflux(ji,jj,jk) - zflux(ji,jj,jk+1) ) & 461 & / fse3t(ji,jj,jk) 462 END DO 463 END DO 464 END DO 465 ! 466 END SUBROUTINE tra_adv_cen2_k 467 468 469 SUBROUTINE quickest( fu, fd, fc, fho, fc_cfl ) 470 !!---------------------------------------------------------------------- 471 !! 472 !!---------------------------------------------------------------------- 473 REAL(wp), INTENT(in) , DIMENSION(jpi,jpj) :: fu, fd, fc, fc_cfl 474 REAL(wp), INTENT(out) , DIMENSION(jpi,jpj) :: fho 475 REAL(wp) , DIMENSION(jpi,jpj) :: zcurv, zcoef1, zcoef2, zcoef3 ! temporary scalars 476 ! 477 zcurv (:,:) = fd(:,:) + fu(:,:) - 2.*fc(:,:) 478 zcoef1(:,:) = 0.5*( fc(:,:) + fd(:,:) ) 479 zcoef2(:,:) = 0.5*fc_cfl(:,:)*( fd(:,:) - fc(:,:) ) 480 zcoef3(:,:) = ( ( 1. - ( fc_cfl(:,:)*fc_cfl(:,:) ) )*r1_6 )*zcurv(:,:) 481 fho (:,:) = zcoef1(:,:) - zcoef2(:,:) - zcoef3(:,:) ! phi_f QUICKEST 482 ! 483 zcoef1(:,:) = fd(:,:) - fu(:,:) ! DEL 484 zcoef2(:,:) = ABS( zcoef1(:,:) ) ! ABS(DEL) 485 zcoef3(:,:) = ABS( zcurv(:,:) ) ! ABS(CURV) 486 ! 487 WHERE ( zcoef3(:,:) >= zcoef2(:,:) ) 488 fho(:,:) = fc(:,:) 489 ELSEWHERE 490 zcoef3(:,:) = fu(:,:) + ( ( fc(:,:) - fu(:,:) )/MAX(fc_cfl(:,:),1.e-9) ) ! phi_REF 491 WHERE ( zcoef1(:,:) >= 0.e0 ) 492 fho(:,:) = MAX(fc(:,:),fho(:,:)) 493 fho(:,:) = MIN(fho(:,:),MIN(zcoef3(:,:),fd(:,:))) 494 ELSEWHERE 495 fho(:,:) = MIN(fc(:,:),fho(:,:)) 496 fho(:,:) = MAX(fho(:,:),MAX(zcoef3(:,:),fd(:,:))) 497 ENDWHERE 498 ENDWHERE 447 zc = puc(ji,jj,jk) ! Courant number 448 zcurv = pfd(ji,jj,jk) + pfu(ji,jj,jk) - 2. * pfc(ji,jj,jk) 449 zcoef1 = 0.5 * ( pfc(ji,jj,jk) + pfd(ji,jj,jk) ) 450 zcoef2 = 0.5 * zc * ( pfd(ji,jj,jk) - pfc(ji,jj,jk) ) 451 zcoef3 = ( 1. - ( zc * zc ) ) * r1_6 * zcurv 452 zfho = zcoef1 - zcoef2 - zcoef3 ! phi_f QUICKEST 453 ! 454 zcoef1 = pfd(ji,jj,jk) - pfu(ji,jj,jk) 455 zcoef2 = ABS( zcoef1 ) 456 zcoef3 = ABS( zcurv ) 457 IF( zcoef3 >= zcoef2 ) THEN 458 zfho = pfc(ji,jj,jk) 459 ELSE 460 zcoef3 = pfu(ji,jj,jk) + ( ( pfc(ji,jj,jk) - pfu(ji,jj,jk) ) / MAX( zc, 1.e-9 ) ) ! phi_REF 461 IF( zcoef1 >= 0. ) THEN 462 zfho = MAX( pfc(ji,jj,jk), zfho ) 463 zfho = MIN( zfho, MIN( zcoef3, pfd(ji,jj,jk) ) ) 464 ELSE 465 zfho = MIN( pfc(ji,jj,jk), zfho ) 466 zfho = MAX( zfho, MAX( zcoef3, pfd(ji,jj,jk) ) ) 467 ENDIF 468 ENDIF 469 puc(ji,jj,jk) = zfho 470 ENDDO 471 ENDDO 472 ENDDO 499 473 ! 500 474 END SUBROUTINE quickest
Note: See TracChangeset
for help on using the changeset viewer.
