Changeset 2528 for trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_ubs.F90
- Timestamp:
- 2010-12-27T18:33:53+01:00 (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_ubs.F90
- Property svn:eol-style deleted
r1528 r2528 4 4 !! Ocean active tracers: horizontal & vertical advective trend 5 5 !!============================================================================== 6 !! History : 9.0 ! 06-08 (L. Debreu, R. Benshila) Original code 6 !! History : 1.0 ! 2006-08 (L. Debreu, R. Benshila) Original code 7 !! 3.3 ! 2010-05 (C. Ethe, G. Madec) merge TRC-TRA + switch from velocity to transport 7 8 !!---------------------------------------------------------------------- 8 9 … … 13 14 USE oce ! ocean dynamics and active tracers 14 15 USE dom_oce ! ocean space and time domain 15 USE trdmod 16 USE trd mod_oce16 USE trdmod_oce ! ocean space and time domain 17 USE trdtra 17 18 USE lib_mpp 18 19 USE lbclnk ! ocean lateral boundary condition (or mpp link) … … 20 21 USE diaptr ! poleward transport diagnostics 21 22 USE dynspg_oce ! choice/control of key cpp for surface pressure gradient 22 USE prtctl23 USE trc_oce ! share passive tracers/Ocean variables 23 24 24 25 IMPLICIT NONE … … 27 28 PUBLIC tra_adv_ubs ! routine called by traadv module 28 29 29 REAL(wp), DIMENSION(jpi,jpj) :: e1e2tr ! = 1/(e1t * e2t)30 LOGICAL :: l_trd ! flag to compute trends or not 30 31 31 32 !! * Substitutions … … 33 34 # include "vectopt_loop_substitute.h90" 34 35 !!---------------------------------------------------------------------- 35 !! OPA 9.0 , LOCEAN-IPSL (2006)36 !! NEMO/OPA 3.3 , NEMO Consortium (2010) 36 37 !! $Id$ 37 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 38 !!---------------------------------------------------------------------- 39 38 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 39 !!---------------------------------------------------------------------- 40 40 CONTAINS 41 41 42 SUBROUTINE tra_adv_ubs( kt, pun, pvn, pwn ) 42 SUBROUTINE tra_adv_ubs ( kt, cdtype, p2dt, pun, pvn, pwn, & 43 & ptb, ptn, pta, kjpt ) 43 44 !!---------------------------------------------------------------------- 44 45 !! *** ROUTINE tra_adv_ubs *** … … 67 68 !! the UBS have been found to be too diffusive. 68 69 !! 69 !! ** Action : - update ( ta,sa) with the now advective tracer trends70 !! ** Action : - update (pta) with the now advective tracer trends 70 71 !! 71 72 !! Reference : Shchepetkin, A. F., J. C. McWilliams, 2005, Ocean Modelling, 9, 347-404. 72 73 !! Farrow, D.E., Stevens, D.P., 1995, J. Phys. Ocean. 25, 1731Ð1741. 73 74 !!---------------------------------------------------------------------- 74 USE oce, ONLY : zwx => ua ! use ua as workspace 75 USE oce, ONLY : zwy => va ! use va as workspace 76 !! 77 INTEGER , INTENT(in) :: kt ! ocean time-step index 78 REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pun ! effective ocean velocity, u_component 79 REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pvn ! effective ocean velocity, v_component 80 REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pwn ! effective ocean velocity, w_component 81 !! 82 INTEGER :: ji, jj, jk ! dummy loop indices 83 REAL(wp) :: zta, zsa, zbtr, zcoef ! temporary scalars 84 REAL(wp) :: zfui, zfp_ui, zfm_ui, zcenut, zcenus ! " " 85 REAL(wp) :: zfvj, zfp_vj, zfm_vj, zcenvt, zcenvs ! " " 86 REAL(wp) :: z_hdivn_x, z_hdivn_y, z_hdivn ! " " 87 REAL(wp), DIMENSION(jpi,jpj) :: zeeu, zeev ! temporary 2D workspace 88 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwz , zww ! temporary 3D workspace 89 REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztu , ztv , zltu , zltv, ztrdt ! " " 90 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zsu , zsv , zlsu , zlsv, ztrds ! " " 91 !!---------------------------------------------------------------------- 92 93 zltu(:,:,:) = 0.e0 94 zltv(:,:,:) = 0.e0 95 zlsu(:,:,:) = 0.e0 96 zlsv(:,:,:) = 0.e0 97 98 IF( kt == nit000 ) THEN 75 USE oce , zwx => ua ! use ua as workspace 76 USE oce , zwy => va ! use va as workspace 77 !! 78 INTEGER , INTENT(in ) :: kt ! ocean time-step index 79 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 80 INTEGER , INTENT(in ) :: kjpt ! number of tracers 81 REAL(wp), DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile of tracer time-step 82 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun, pvn, pwn ! 3 ocean velocity components 83 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before and now tracer fields 84 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 85 !! 86 INTEGER :: ji, jj, jk, jn ! dummy loop indices 87 REAL(wp) :: ztra, zbtr, zcoef ! local scalars 88 REAL(wp) :: zfp_ui, zfm_ui, zcenut ! - - 89 REAL(wp) :: zfp_vj, zfm_vj, zcenvt ! - - 90 REAL(wp) :: z2dtt ! - - 91 REAL(wp) :: ztak, zfp_wk, zfm_wk ! - - 92 REAL(wp) :: zeeu, zeev, z_hdivn ! - - 93 REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztu, ztv, zltu , zltv ! 3D workspace 94 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zti, ztw ! - - 95 !!---------------------------------------------------------------------- 96 97 IF( kt == nit000 ) THEN 99 98 IF(lwp) WRITE(numout,*) 100 IF(lwp) WRITE(numout,*) 'tra_adv_ubs : horizontal UBS advection scheme '99 IF(lwp) WRITE(numout,*) 'tra_adv_ubs : horizontal UBS advection scheme on ', cdtype 101 100 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 102 101 ! 103 e1e2tr(:,:) = 1. / ( e1t(:,:) * e2t(:,:) ) 102 l_trd = .FALSE. 103 IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 104 104 ENDIF 105 106 ! Save ta and sa trends 107 ztrdt(:,:,:) = ta(:,:,:) 108 ztrds(:,:,:) = sa(:,:,:) 109 110 zcoef = 1./6. 111 ! ! =============== 112 DO jk = 1, jpkm1 ! Horizontal slab 113 ! ! =============== 114 115 ! Initialization of metric arrays (for z- or s-coordinates) 116 DO jj = 1, jpjm1 117 DO ji = 1, fs_jpim1 ! vector opt. 118 #if defined key_zco 119 ! z-coordinates, no vertical scale factors 120 zeeu(ji,jj) = e2u(ji,jj) / e1u(ji,jj) * umask(ji,jj,jk) 121 zeev(ji,jj) = e1v(ji,jj) / e2v(ji,jj) * vmask(ji,jj,jk) 122 #else 123 ! s-coordinates, vertical scale factor are used 124 zeeu(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) * umask(ji,jj,jk) 125 zeev(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) * vmask(ji,jj,jk) 126 #endif 127 END DO 128 END DO 129 130 ! Laplacian 131 ! First derivative (gradient) 132 DO jj = 1, jpjm1 133 DO ji = 1, fs_jpim1 ! vector opt. 134 ztu(ji,jj,jk) = zeeu(ji,jj) * ( tb(ji+1,jj ,jk) - tb(ji,jj,jk) ) 135 zsu(ji,jj,jk) = zeeu(ji,jj) * ( sb(ji+1,jj ,jk) - sb(ji,jj,jk) ) 136 ztv(ji,jj,jk) = zeev(ji,jj) * ( tb(ji ,jj+1,jk) - tb(ji,jj,jk) ) 137 zsv(ji,jj,jk) = zeev(ji,jj) * ( sb(ji ,jj+1,jk) - sb(ji,jj,jk) ) 138 END DO 139 END DO 140 ! Second derivative (divergence) 141 DO jj = 2, jpjm1 142 DO ji = fs_2, fs_jpim1 ! vector opt. 143 #if ! defined key_zco 144 zcoef = 1. / ( 6. * fse3t(ji,jj,jk) ) 145 #endif 146 zltu(ji,jj,jk) = ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) ) * zcoef 147 zlsu(ji,jj,jk) = ( zsu(ji,jj,jk) - zsu(ji-1,jj,jk) ) * zcoef 148 zltv(ji,jj,jk) = ( ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) * zcoef 149 zlsv(ji,jj,jk) = ( zsv(ji,jj,jk) - zsv(ji,jj-1,jk) ) * zcoef 150 END DO 151 END DO 152 ! ! ================= 153 END DO ! End of slab 154 ! ! ================= 155 156 ! Lateral boundary conditions on the laplacian (zlt,zls) (unchanged sgn) 157 CALL lbc_lnk( zltu, 'T', 1. ) ; CALL lbc_lnk( zlsu, 'T', 1. ) 158 CALL lbc_lnk( zltv, 'T', 1. ) ; CALL lbc_lnk( zlsv, 'T', 1. ) 159 160 ! ! =============== 161 DO jk = 1, jpkm1 ! Horizontal slab 162 ! ! =============== 163 ! Horizontal advective fluxes 164 DO jj = 1, jpjm1 165 DO ji = 1, fs_jpim1 ! vector opt. 166 ! volume fluxes * 1/2 167 #if defined key_zco 168 zfui = 0.5 * e2u(ji,jj) * pun(ji,jj,jk) 169 zfvj = 0.5 * e1v(ji,jj) * pvn(ji,jj,jk) 170 #else 171 zfui = 0.5 * e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk) 172 zfvj = 0.5 * e1v(ji,jj) * fse3v(ji,jj,jk) * pvn(ji,jj,jk) 173 #endif 174 ! upstream scheme 175 zfp_ui = zfui + ABS( zfui ) 176 zfp_vj = zfvj + ABS( zfvj ) 177 zfm_ui = zfui - ABS( zfui ) 178 zfm_vj = zfvj - ABS( zfvj ) 179 ! centered scheme 180 zcenut = zfui * ( tn(ji,jj,jk) + tn(ji+1,jj ,jk) ) 181 zcenvt = zfvj * ( tn(ji,jj,jk) + tn(ji ,jj+1,jk) ) 182 zcenus = zfui * ( sn(ji,jj,jk) + sn(ji+1,jj ,jk) ) 183 zcenvs = zfvj * ( sn(ji,jj,jk) + sn(ji ,jj+1,jk) ) 184 ! mixed centered / upstream scheme 185 zwx(ji,jj,jk) = zcenut - zfp_ui * zltu(ji,jj,jk) -zfm_ui * zltu(ji+1,jj,jk) 186 zwy(ji,jj,jk) = zcenvt - zfp_vj * zltv(ji,jj,jk) -zfm_vj * zltv(ji,jj+1,jk) 187 zww(ji,jj,jk) = zcenus - zfp_ui * zlsu(ji,jj,jk) -zfm_ui * zlsu(ji+1,jj,jk) 188 zwz(ji,jj,jk) = zcenvs - zfp_vj * zlsv(ji,jj,jk) -zfm_vj * zlsv(ji,jj+1,jk) 189 END DO 190 END DO 191 192 ! Tracer flux divergence at t-point added to the general trend 193 DO jj = 2, jpjm1 194 DO ji = fs_2, fs_jpim1 ! vector opt. 195 ! horizontal advective trends 196 #if defined key_zco 197 zbtr = e1e2tr(ji,jj) 198 #else 199 zbtr = e1e2tr(ji,jj) / fse3t(ji,jj,jk) 200 #endif 201 zta = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk) & 202 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk) ) 203 zsa = - zbtr * ( zww(ji,jj,jk) - zww(ji-1,jj ,jk) & 204 & + zwz(ji,jj,jk) - zwz(ji ,jj-1,jk) ) 205 ! add it to the general tracer trends 206 ta(ji,jj,jk) = ta(ji,jj,jk) + zta 207 sa(ji,jj,jk) = sa(ji,jj,jk) + zsa 208 END DO 209 END DO 210 ! ! =============== 211 END DO ! End of slab 212 ! ! =============== 213 214 ! Horizontal trend used in tra_adv_ztvd subroutine 215 zltu(:,:,:) = ta(:,:,:) - ztrdt(:,:,:) 216 zlsu(:,:,:) = sa(:,:,:) - ztrds(:,:,:) 217 218 ! 3. Save the horizontal advective trends for diagnostic 219 ! ------------------------------------------------------ 220 IF( l_trdtra ) THEN 221 ! Recompute the hoizontal advection zta & zsa trends computed 222 ! at the step 2. above in making the difference between the new 223 ! trends and the previous one ta()/sa - ztrdt()/ztrds() and add 224 ! the term tn()/sn()*hdivn() to recover the Uh gradh(T/S) trends 225 ztrdt(:,:,:) = 0.e0 ; ztrds(:,:,:) = 0.e0 226 ! 227 ! T/S ZONAL advection trends 105 ! 106 ! ! =========== 107 DO jn = 1, kjpt ! tracer loop 108 ! ! =========== 109 ! 1. Bottom value : flux set to zero 110 ! ---------------------------------- 111 zltu(:,:,jpk) = 0.e0 ; zltv(:,:,jpk) = 0.e0 112 ! 113 DO jk = 1, jpkm1 ! Horizontal slab 114 ! 115 ! Laplacian 116 DO jj = 1, jpjm1 ! First derivative (gradient) 117 DO ji = 1, fs_jpim1 ! vector opt. 118 zeeu = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) * umask(ji,jj,jk) 119 zeev = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) * vmask(ji,jj,jk) 120 ztu(ji,jj,jk) = zeeu * ( ptb(ji+1,jj ,jk,jn) - ptb(ji,jj,jk,jn) ) 121 ztv(ji,jj,jk) = zeev * ( ptb(ji ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) 122 END DO 123 END DO 124 DO jj = 2, jpjm1 ! Second derivative (divergence) 125 DO ji = fs_2, fs_jpim1 ! vector opt. 126 zcoef = 1. / ( 6. * fse3t(ji,jj,jk) ) 127 zltu(ji,jj,jk) = ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) ) * zcoef 128 zltv(ji,jj,jk) = ( ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) * zcoef 129 END DO 130 END DO 131 ! 132 END DO ! End of slab 133 CALL lbc_lnk( zltu, 'T', 1. ) ; CALL lbc_lnk( zltv, 'T', 1. ) ! Lateral boundary cond. (unchanged sgn) 134 135 ! 136 ! Horizontal advective fluxes 137 DO jk = 1, jpkm1 ! Horizontal slab 138 DO jj = 1, jpjm1 139 DO ji = 1, fs_jpim1 ! vector opt. 140 ! upstream transport 141 zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) ) 142 zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) ) 143 zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) ) 144 zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) ) 145 ! centered scheme 146 zcenut = 0.5 * pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj ,jk,jn) ) 147 zcenvt = 0.5 * pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji ,jj+1,jk,jn) ) 148 ! UBS scheme 149 zwx(ji,jj,jk) = zcenut - zfp_ui * zltu(ji,jj,jk) - zfm_ui * zltu(ji+1,jj,jk) 150 zwy(ji,jj,jk) = zcenvt - zfp_vj * zltv(ji,jj,jk) - zfm_vj * zltv(ji,jj+1,jk) 151 END DO 152 END DO 153 END DO ! End of slab 154 155 zltu(:,:,:) = pta(:,:,:,jn) ! store pta trends 156 157 ! Horizontal advective trends 228 158 DO jk = 1, jpkm1 159 ! Tracer flux divergence at t-point added to the general trend 229 160 DO jj = 2, jpjm1 230 161 DO ji = fs_2, fs_jpim1 ! vector opt. 231 !-- Compute zonal divergence by splitting hdivn (see divcur.F90) 232 #if defined key_zco 233 zbtr = e1e2tr(ji,jj) 234 z_hdivn_x = ( e2u(ji ,jj) * pun(ji ,jj,jk) & 235 & - e2u(ji-1,jj) * pun(ji-1,jj,jk) ) * zbtr 236 #else 237 zbtr = e1e2tr(ji,jj) / fse3t(ji,jj,jk) 238 z_hdivn_x = ( e2u(ji ,jj) * fse3u(ji ,jj,jk) * pun(ji ,jj,jk) & 239 & - e2u(ji-1,jj) * fse3u(ji-1,jj,jk) * pun(ji-1,jj,jk) ) * zbtr 240 #endif 241 ztrdt(ji,jj,jk) = - ( zwx(ji,jj,jk) - zwx(ji-1,jj,jk) ) * zbtr + tn(ji,jj,jk) * z_hdivn_x 242 ztrds(ji,jj,jk) = - ( zww(ji,jj,jk) - zww(ji-1,jj,jk) ) * zbtr + sn(ji,jj,jk) * z_hdivn_x 243 END DO 244 END DO 245 END DO 246 CALL trd_mod(ztrdt, ztrds, jptra_trd_xad, 'TRA', kt) ! save the trends 247 ! 248 ! T/S MERIDIONAL advection trends 162 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 163 ! horizontal advective 164 ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk) & 165 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk) ) 166 ! add it to the general tracer trends 167 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 168 END DO 169 END DO 170 ! 171 END DO ! End of slab 172 173 ! Horizontal trend used in tra_adv_ztvd subroutine 174 zltu(:,:,:) = pta(:,:,:,jn) - zltu(:,:,:) 175 176 ! 3. Save the horizontal advective trends for diagnostic 177 ! ------------------------------------------------------ 178 ! ! trend diagnostics (contribution of upstream fluxes) 179 IF( l_trd ) THEN 180 CALL trd_tra( kt, cdtype, jn, jptra_trd_xad, zwx, pun, ptn(:,:,:,jn) ) 181 CALL trd_tra( kt, cdtype, jn, jptra_trd_yad, zwy, pvn, ptn(:,:,:,jn) ) 182 END IF 183 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) 184 IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN 185 IF( jn == jp_tem ) htr_adv(:) = ptr_vj( zwy(:,:,:) ) 186 IF( jn == jp_sal ) str_adv(:) = ptr_vj( zwy(:,:,:) ) 187 ENDIF 188 189 ! TVD scheme for the vertical direction 190 ! ---------------------- 191 IF( l_trd ) zltv(:,:,:) = pta(:,:,:,jn) ! store pta if trend diag. 192 193 ! Bottom value : flux set to zero 194 ztw(:,:,jpk) = 0.e0 ; zti(:,:,jpk) = 0.e0 195 196 ! Surface value 197 IF( lk_vvl ) THEN ; ztw(:,:,1) = 0.e0 ! variable volume : flux set to zero 198 ELSE ; ztw(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn) ! free constant surface 199 ENDIF 200 ! upstream advection with initial mass fluxes & intermediate update 201 ! ------------------------------------------------------------------- 202 ! Interior value 203 DO jk = 2, jpkm1 204 DO jj = 1, jpj 205 DO ji = 1, jpi 206 zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) ) 207 zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) ) 208 ztw(ji,jj,jk) = 0.5 * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) ) 209 END DO 210 END DO 211 END DO 212 ! update and guess with monotonic sheme 249 213 DO jk = 1, jpkm1 214 z2dtt = p2dt(jk) 250 215 DO jj = 2, jpjm1 251 216 DO ji = fs_2, fs_jpim1 ! vector opt. 252 !-- Compute merid. divergence by splitting hdivn (see divcur.F90) 253 #if defined key_zco 254 zbtr = e1e2tr(ji,jj) 255 z_hdivn_y = ( e1v(ji, jj) * pvn(ji,jj ,jk) & 256 & - e1v(ji,jj-1) * pvn(ji,jj-1,jk) ) * zbtr 257 #else 258 zbtr = e1e2tr(ji,jj) / fse3t(ji,jj,jk) 259 z_hdivn_y = ( e1v(ji, jj) * fse3v(ji,jj ,jk) * pvn(ji,jj ,jk) & 260 & - e1v(ji,jj-1) * fse3v(ji,jj-1,jk) * pvn(ji,jj-1,jk) ) * zbtr 261 #endif 262 ztrdt(ji,jj,jk) = - ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) ) * zbtr + tn(ji,jj,jk) * z_hdivn_y 263 ztrds(ji,jj,jk) = - ( zwz(ji,jj,jk) - zwz(ji,jj-1,jk) ) * zbtr + sn(ji,jj,jk) * z_hdivn_y 264 END DO 265 END DO 266 END DO 267 CALL trd_mod(ztrdt, ztrds, jptra_trd_yad, 'TRA', kt) ! save the trends 217 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 218 ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * zbtr 219 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztak 220 zti(ji,jj,jk) = ( ptb(ji,jj,jk,jn) + z2dtt * ( ztak + zltu(ji,jj,jk) ) ) * tmask(ji,jj,jk) 221 END DO 222 END DO 223 END DO 268 224 ! 269 ENDIF 270 271 IF(ln_ctl) CALL prt_ctl( tab3d_1=ta, clinfo1=' ubs had - Ta: ', mask1=tmask, & 272 & tab3d_2=sa, clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' ) 273 274 ! "zonal" mean advective heat and salt transport 275 IF( ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN 276 IF( lk_zco ) THEN 277 DO jk = 1, jpkm1 225 CALL lbc_lnk( zti, 'T', 1. ) ! Lateral boundary conditions on zti, zsi (unchanged sign) 226 227 ! antidiffusive flux : high order minus low order 228 ztw(:,:,1) = 0.e0 ! Surface value 229 DO jk = 2, jpkm1 ! Interior value 230 DO jj = 1, jpj 231 DO ji = 1, jpi 232 ztw(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) ) - ztw(ji,jj,jk) 233 END DO 234 END DO 235 END DO 236 ! 237 CALL nonosc_z( ptb(:,:,:,jn), ztw, zti, p2dt ) ! monotonicity algorithm 238 239 ! final trend with corrected fluxes 240 DO jk = 1, jpkm1 241 DO jj = 2, jpjm1 242 DO ji = fs_2, fs_jpim1 ! vector opt. 243 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 244 ! k- vertical advective trends 245 ztra = - zbtr * ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) 246 ! added to the general tracer trends 247 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 248 END DO 249 END DO 250 END DO 251 252 ! Save the final vertical advective trends 253 IF( l_trd ) THEN ! vertical advective trend diagnostics 254 DO jk = 1, jpkm1 ! (compute -w.dk[ptn]= -dk[w.ptn] + ptn.dk[w]) 278 255 DO jj = 2, jpjm1 279 256 DO ji = fs_2, fs_jpim1 ! vector opt. 280 zwy(ji,jj,jk) = zwy(ji,jj,jk) * fse3v(ji,jj,jk) 281 zwz(ji,jj,jk) = zwz(ji,jj,jk) * fse3v(ji,jj,jk) 257 zbtr = 1.e0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 258 z_hdivn = ( pwn(ji,jj,jk) - pwn(ji,jj,jk+1) ) * zbtr 259 zltv(ji,jj,jk) = pta(ji,jj,jk,jn) - zltv(ji,jj,jk) + ptn(ji,jj,jk,jn) * z_hdivn 282 260 END DO 283 261 END DO 284 262 END DO 263 CALL trd_tra( kt, cdtype, jn, jptra_trd_zad, zltv ) 285 264 ENDIF 286 pht_adv(:) = ptr_vj( zwy(:,:,:) )287 pst_adv(:) = ptr_vj( zwz(:,:,:) )288 ENDIF289 290 ! II. Vertical advection291 ! ----------------------292 IF( l_trdtra ) THEN ! Save ta and sa trends293 ztrdt(:,:,:) = ta(:,:,:)294 ztrds(:,:,:) = sa(:,:,:)295 ENDIF296 297 ! TVD scheme the vertical direction298 CALL tra_adv_ztvd(kt, pwn, zltu, zlsu)299 300 IF( l_trdtra ) THEN ! Save the final vertical advective trends301 DO jk = 1, jpkm1302 DO jj = 2, jpjm1303 DO ji = fs_2, fs_jpim1 ! vector opt.304 #if defined key_zco305 zbtr = e1e2tr(ji,jj)306 z_hdivn_x = e2u(ji,jj)*pun(ji,jj,jk) - e2u(ji-1,jj)*pun(ji-1,jj,jk)307 z_hdivn_y = e1v(ji,jj)*pvn(ji,jj,jk) - e1v(ji,jj-1)*pvn(ji,jj-1,jk)308 #else309 zbtr = e1e2tr(ji,jj) / fse3t(ji,jj,jk)310 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)311 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)312 #endif313 z_hdivn = (z_hdivn_x + z_hdivn_y) * zbtr314 zbtr = e1e2tr(ji,jj) / fse3t(ji,jj,jk)315 ztrdt(ji,jj,jk) = ta(ji,jj,jk) - ztrdt(ji,jj,jk) - tn(ji,jj,jk) * z_hdivn316 ztrds(ji,jj,jk) = sa(ji,jj,jk) - ztrds(ji,jj,jk) - sn(ji,jj,jk) * z_hdivn317 END DO318 END DO319 END DO320 CALL trd_mod(ztrdt, ztrds, jptra_trd_zad, 'TRA', kt) ! <<< ADD TO PREVIOUSLY COMPUTED321 265 ! 322 ENDIF 323 324 IF(ln_ctl) CALL prt_ctl( tab3d_1=ta, clinfo1=' ubs zad - Ta: ', mask1=tmask, & 325 & tab3d_2=sa, clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra') 266 ENDDO 326 267 ! 327 268 END SUBROUTINE tra_adv_ubs 328 269 329 270 330 SUBROUTINE tra_adv_ztvd( kt, pwn, zttrd, zstrd ) 331 !!---------------------------------------------------------------------- 332 !! *** ROUTINE tra_adv_ztvd *** 333 !! 334 !! ** Purpose : Compute the now trend due to total advection of 335 !! tracers and add it to the general trend of tracer equations 336 !! 337 !! ** Method : TVD scheme, i.e. 2nd order centered scheme with 338 !! corrected flux (monotonic correction) 339 !! note: - this advection scheme needs a leap-frog time scheme 340 !! 341 !! ** Action : - update (ta,sa) with the now advective tracer trends 342 !! - save the trends in (ztrdt,ztrds) ('key_trdtra') 343 !!---------------------------------------------------------------------- 344 INTEGER , INTENT(in) :: kt ! ocean time-step 345 REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pwn ! verical effective velocity 346 REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: zttrd, zstrd ! lateral advective trends on T & S 347 !! 348 INTEGER :: ji, jj, jk ! dummy loop indices 349 REAL(wp) :: z2dtt, zbtr, zew, z2 ! temporary scalar 350 REAL(wp) :: ztak, zfp_wk ! " " 351 REAL(wp) :: zsak, zfm_wk ! " " 352 REAL(wp), DIMENSION (jpi,jpj,jpk) :: zti, ztw ! temporary 3D workspace 353 REAL(wp), DIMENSION (jpi,jpj,jpk) :: zsi, zsw ! " " 354 !!---------------------------------------------------------------------- 355 356 IF( kt == nit000 .AND. lwp ) THEN 357 WRITE(numout,*) 358 WRITE(numout,*) 'tra_adv_ztvd : vertical TVD advection scheme' 359 WRITE(numout,*) '~~~~~~~~~~~~' 360 ENDIF 361 362 IF( neuler == 0 .AND. kt == nit000 ) THEN ; z2 = 1. 363 ELSE ; z2 = 2. 364 ENDIF 365 366 ! Bottom value : flux set to zero 367 ! -------------- 368 ztw(:,:,jpk) = 0.e0 ; zsw(:,:,jpk) = 0.e0 369 zti (:,:,:) = 0.e0 ; zsi (:,:,:) = 0.e0 370 371 372 ! upstream advection with initial mass fluxes & intermediate update 373 ! ------------------------------------------------------------------- 374 ! Surface value 375 IF( lk_vvl ) THEN 376 ! variable volume : flux set to zero 377 ztw(:,:,1) = 0.e0 378 zsw(:,:,1) = 0.e0 379 ELSE 380 ! free surface-constant volume 381 DO jj = 1, jpj 382 DO ji = 1, jpi 383 zew = e1t(ji,jj) * e2t(ji,jj) * pwn(ji,jj,1) 384 ztw(ji,jj,1) = zew * tb(ji,jj,1) 385 zsw(ji,jj,1) = zew * sb(ji,jj,1) 386 END DO 387 END DO 388 ENDIF 389 390 ! Interior value 391 DO jk = 2, jpkm1 392 DO jj = 1, jpj 393 DO ji = 1, jpi 394 zew = 0.5 * e1t(ji,jj) * e2t(ji,jj) * pwn(ji,jj,jk) 395 zfp_wk = zew + ABS( zew ) 396 zfm_wk = zew - ABS( zew ) 397 ztw(ji,jj,jk) = zfp_wk * tb(ji,jj,jk) + zfm_wk * tb(ji,jj,jk-1) 398 zsw(ji,jj,jk) = zfp_wk * sb(ji,jj,jk) + zfm_wk * sb(ji,jj,jk-1) 399 END DO 400 END DO 401 END DO 402 403 ! update and guess with monotonic sheme 404 DO jk = 1, jpkm1 405 z2dtt = z2 * rdttra(jk) 406 DO jj = 2, jpjm1 407 DO ji = fs_2, fs_jpim1 ! vector opt. 408 zbtr = 1./ ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 409 ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * zbtr 410 zsak = - ( zsw(ji,jj,jk) - zsw(ji,jj,jk+1) ) * zbtr 411 ta(ji,jj,jk) = ta(ji,jj,jk) + ztak 412 sa(ji,jj,jk) = sa(ji,jj,jk) + zsak 413 zti (ji,jj,jk) = ( tb(ji,jj,jk) + z2dtt * ( ztak + zttrd(ji,jj,jk) ) ) * tmask(ji,jj,jk) 414 zsi (ji,jj,jk) = ( sb(ji,jj,jk) + z2dtt * ( zsak + zstrd(ji,jj,jk) ) ) * tmask(ji,jj,jk) 415 END DO 416 END DO 417 END DO 418 419 ! Lateral boundary conditions on zti, zsi (unchanged sign) 420 CALL lbc_lnk( zti, 'T', 1. ) 421 CALL lbc_lnk( zsi, 'T', 1. ) 422 423 424 ! antidiffusive flux : high order minus low order 425 ! ------------------------------------------------- 426 ! Surface value 427 ztw(:,:,1) = 0.e0 ; zsw(:,:,1) = 0.e0 428 429 ! Interior value 430 DO jk = 2, jpkm1 431 DO jj = 1, jpj 432 DO ji = 1, jpi 433 zew = 0.5 * e1t(ji,jj) * e2t(ji,jj) * pwn(ji,jj,jk) 434 ztw(ji,jj,jk) = zew * ( tn(ji,jj,jk) + tn(ji,jj,jk-1) ) - ztw(ji,jj,jk) 435 zsw(ji,jj,jk) = zew * ( sn(ji,jj,jk) + sn(ji,jj,jk-1) ) - zsw(ji,jj,jk) 436 END DO 437 END DO 438 END DO 439 440 ! monotonicity algorithm 441 ! ------------------------ 442 CALL nonosc_z( tb, ztw, zti, z2 ) 443 CALL nonosc_z( sb, zsw, zsi, z2 ) 444 445 446 ! final trend with corrected fluxes 447 ! ----------------------------------- 448 DO jk = 1, jpkm1 449 DO jj = 2, jpjm1 450 DO ji = fs_2, fs_jpim1 ! vector opt. 451 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 452 ! k- vertical advective trends 453 ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * zbtr 454 zsak = - ( zsw(ji,jj,jk) - zsw(ji,jj,jk+1) ) * zbtr 455 ! add them to the general tracer trends 456 ta(ji,jj,jk) = ta(ji,jj,jk) + ztak 457 sa(ji,jj,jk) = sa(ji,jj,jk) + zsak 458 END DO 459 END DO 460 END DO 461 ! 462 END SUBROUTINE tra_adv_ztvd 463 464 465 SUBROUTINE nonosc_z( pbef, pcc, paft, prdt ) 271 SUBROUTINE nonosc_z( pbef, pcc, paft, p2dt ) 466 272 !!--------------------------------------------------------------------- 467 273 !! *** ROUTINE nonosc_z *** … … 476 282 !! in-space based differencing for fluid 477 283 !!---------------------------------------------------------------------- 478 REAL(wp), INTENT(in ) :: prdt ! ???479 REAL(wp), INTENT(inout),DIMENSION (jpi,jpj,jpk) :: pbef ! before field284 REAL(wp), INTENT(in ), DIMENSION(jpk) :: p2dt ! vertical profile of tracer time-step 285 REAL(wp), DIMENSION (jpi,jpj,jpk) :: pbef ! before field 480 286 REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) :: paft ! after field 481 287 REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) :: pcc ! monotonic flux in the k direction … … 525 331 pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) 526 332 paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) 527 333 528 334 529 335 ! 2. Positive and negative part of fluxes and beta terms … … 531 337 532 338 DO jk = 1, jpkm1 533 z2dtt = p rdt * rdttra(jk)339 z2dtt = p2dt(jk) 534 340 DO jj = 2, jpjm1 535 341 DO ji = fs_2, fs_jpim1 ! vector opt. … … 544 350 END DO 545 351 END DO 546 547 352 ! monotonic flux in the k direction, i.e. pcc 548 353 ! ------------------------------------------- … … 550 355 DO jj = 2, jpjm1 551 356 DO ji = fs_2, fs_jpim1 ! vector opt. 552 553 357 za = MIN( 1., zbetdo(ji,jj,jk), zbetup(ji,jj,jk-1) ) 554 358 zb = MIN( 1., zbetup(ji,jj,jk), zbetdo(ji,jj,jk-1) )
Note: See TracChangeset
for help on using the changeset viewer.
