Changeset 791 for branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_cen2.F90
- Timestamp:
- 2008-01-12T21:33:34+01:00 (16 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_cen2.F90
r786 r791 4 4 !! Ocean active tracers: horizontal & vertical advective trend 5 5 !!============================================================================== 6 !! History : 8.2 ! 01-08 (G. Madec, E. Durand) trahad+trazad = traadv7 !! 8.5 ! 02-06 (G. Madec) F90: Free form and module8 !! NEMO 1.0 ! 05-11 (V. Garnier) Surface pressure gradient organization9 !! - ! 05-11 (V. Garnier) Surface pressure gradient organization10 !! - ! 06-04 (R. Benshila, G. Madec) Step reorganization11 !! 2.4 ! 08-01 (G. Madec) Merge TRA-TRC6 !! History : 8.2 ! 2001-08 (G. Madec, E. Durand) trahad+trazad = traadv 7 !! 8.5 ! 2002-06 (G. Madec) F90: Free form and module 8 !! NEMO 1.0 ! 2005-11 (V. Garnier) Surface pressure gradient organization 9 !! - ! 2005-11 (V. Garnier) Surface pressure gradient organization 10 !! - ! 2006-04 (R. Benshila, G. Madec) Step reorganization 11 !! 2.4 ! 2008-01 (G. Madec) merge TRC-TRA + switch from velocity to transport 12 12 !!---------------------------------------------------------------------- 13 13 … … 34 34 PUBLIC tra_adv_cen2 ! routine called by step.F90 35 35 36 REAL(wp), DIMENSION(jpi,jpj) :: btr2 ! inverse of T-point surface [1/(e1t*e2t)]37 38 36 !! * Substitutions 39 37 # include "domzgr_substitute.h90" … … 41 39 !!---------------------------------------------------------------------- 42 40 !! NEMO/OPA 2.4 , LOCEAN-IPSL (2008) 43 !! $Id :$41 !! $Id$ 44 42 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 45 43 !!---------------------------------------------------------------------- … … 62 60 !! Part 0 : compute the upstream / centered flag 63 61 !! (3D array, zind, defined at T-point (0<zind<1)) 64 !! Part I : horizontal advection 65 !! * centered flux: 66 !! zcenu = e2u*e3u un mi(ptn) 67 !! zcenv = e1v*e3v vn mj(ptn) 68 !! * upstream flux: 69 !! zupsu = e2u*e3u un (ptb(i) or ptb(i-1) ) [un>0 or <0] 70 !! zupsv = e1v*e3v vn (tb(j) or tb(j-1) ) [vn>0 or <0] 71 !! * mixed upstream / centered horizontal advection scheme 62 !! Part I : advective fluxes 63 !! * centered flux: zcenu = pun mi(ptn) (idem for v and w) 64 !! * upstream flux: zupsu = pun (ptb(i) or ptb(i-1) ) [un>0 or <0] 65 !! * mixed upstream / centered flux: 72 66 !! zcofi = max(zind(i+1), zind(i)) 73 !! zcofj = max(zind(j+1), zind(j))74 67 !! zwx = zcofi * zupsu + (1-zcofi) * zcenu 75 !! zwy = zcofj * zupsv + (1-zcofj) * zcenv68 !! Part III : advective trend 76 69 !! * horizontal advective trend (divergence of the fluxes) 77 !! zta = 1/(e1t*e2t*e3t) { di-1[zwx] + dj-1[zwy] } 78 !! * Add this trend now to the general trend of tracer (pta): 79 !! pta = pta + zta 80 !! * trend diagnostic (lk_trdtra=T): the trend is 81 !! saved for diagnostics. The trends saved is expressed as 82 !! Uh.gradh(T), i.e. 83 !! save trend = zta + ptn divn 84 !! In addition, the advective trend in the two horizontal direc- 85 !! tion is also re-computed as Uh gradh(T). Indeed hadt+ptn divn is 86 !! equal to (in s-coordinates, and similarly in z-coord.): 87 !! zta+ptn*divn=1/(e1t*e2t*e3t) { mi-1( e2u*e3u un di[ptn] ) 88 !! +mj-1( e1v*e3v vn mj[ptn] ) } 89 !! NB:in z-coordinate - full step (ln_zco=T) e3u=e3v=e3t, so 90 !! they vanish from the expression of the flux and divergence. 91 !! 92 !! Part II : vertical advection 93 !! the advective trend is computed as follows : 94 !! zta = 1/e3t dk+1[ zwx ] 95 !! where the vertical advective flux, zwx, is given by : 96 !! zwx = zcofk * zupst + (1-zcofk) * zcent 97 !! with 98 !! zupsv = upstream flux = wn * (tb(k) or tb(k-1) ) [wn>0 or <0] 99 !! zcenu = centered flux = wn * mk(ptn) 100 !! The surface boundary condition is : 101 !! rigid-lid (lk_dynspg_frd = T) : zero advective flux 102 !! free-surf (lk_dynspg_fsc = T) : wn(:,:,1) * ptn(:,:,1) 103 !! Add this trend now to the general trend of tracer (ta,sa): 104 !! (ta,sa) = (ta,sa) + ( zta , zsa ) 105 !! Trend diagnostic ('key_trdtra' defined): the trend is 106 !! saved for diagnostics. The trends saved is expressed as : 107 !! save trend = w.gradz(T) = zta - ptn divn. 108 !! 109 !! ** Action : - update (ta,sa) with the now advective tracer trends 70 !! added to the general trend of tracer (pta): 71 !! pta = pta - 1/(e1t*e2t*e3t) { di-1[zwx] + dj-1[zwy] + dk[zwz] } 72 !! * trend diagnostic (lk_trdtra=T): the advective fluxes are 73 !! where the 3 part of advective trend are recomputed as U.grad[T] 74 !! and then used in on-line diagnostics. 75 !! 76 !! NB: The surface boundary condition is no flux for both 77 !! rigid-lid (lk_dynspg_frd = T) and non-linear free-surface 78 !! but a flux (pwn(1) * ptn(1)) cross the surface in linear 79 !! free-surface case (lk_dynspg_fsc = T) 80 !! 81 !! ** Action : - update pta with the 2nd order centered advective trends 110 82 !! - trend diagnostics (lk_trdtra=T) 111 83 !!---------------------------------------------------------------------- … … 118 90 !! 119 91 INTEGER :: ji, jj, jk ! dummy loop indices 120 REAL(wp) :: zbtr , zta, zhw, ze3tr! temporary scalars121 REAL(wp) :: zcofi, z fui, zcenut, zupsut, zfp_ui, zfm_ui! " "122 REAL(wp) :: zcofj, z fvj, zcenvt, zupsvt, zfp_vj, zfm_vj! " "123 REAL(wp) :: zcofk, zcent , zupst , zfp_w , zfm_w! " "124 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwx, zwy, z ind ! 3D workspace92 REAL(wp) :: zbtr ! temporary scalars 93 REAL(wp) :: zcofi, zcenut, zupsut, zfp_ui, zfm_ui ! " " 94 REAL(wp) :: zcofj, zcenvt, zupsvt, zfp_vj, zfm_vj ! " " 95 REAL(wp) :: zcofk, zcent , zupst , zfp_w , zfm_w ! " " 96 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwx, zwy, zwz, zind ! 3D workspace 125 97 !!---------------------------------------------------------------------- 126 98 … … 129 101 IF(lwp) WRITE(numout,*) 'tra_adv_cen2 : 2nd order centered advection scheme' 130 102 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 131 !132 btr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) )133 103 ENDIF 134 104 … … 150 120 END DO 151 121 152 ! I. Horizontal advection 153 ! ----------------------- 154 155 ! ! =============== 156 DO jk = 1, jpkm1 ! Horizontal slab 157 ! ! =============== 158 ! 159 DO jj = 1, jpjm1 ! Horizontal advective fluxes 122 ! I. Advective fluxes 123 ! ------------------- 124 DO jk = 1, jpkm1 ! Horizontal advective fluxes 125 DO jj = 1, jpjm1 160 126 DO ji = 1, fs_jpim1 ! vector opt. 161 127 ! upstream indicator 162 128 zcofi = MAX( zind(ji+1,jj,jk), zind(ji,jj,jk) ) 163 129 zcofj = MAX( zind(ji,jj+1,jk), zind(ji,jj,jk) ) 164 ! volume fluxes * 1/2165 #if defined key_zco166 zfui = 0.5 * e2u(ji,jj) * pun(ji,jj,jk)167 zfvj = 0.5 * e1v(ji,jj) * pvn(ji,jj,jk)168 #else169 zfui = 0.5 * e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk)170 zfvj = 0.5 * e1v(ji,jj) * fse3v(ji,jj,jk) * pvn(ji,jj,jk)171 #endif172 130 ! upstream scheme 173 zfp_ui = zfui + ABS( zfui)174 zf p_vj = zfvj + ABS( zfvj)175 zf m_ui = zfui - ABS( zfui)176 zfm_vj = zfvj - ABS( zfvj)131 zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) ) 132 zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) ) 133 zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) ) 134 zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) ) 177 135 zupsut = zfp_ui * ptb(ji,jj,jk) + zfm_ui * ptb(ji+1,jj ,jk) 178 136 zupsvt = zfp_vj * ptb(ji,jj,jk) + zfm_vj * ptb(ji ,jj+1,jk) 179 137 ! centered scheme 180 zcenut = zfui* ( ptn(ji,jj,jk) + ptn(ji+1,jj ,jk) )181 zcenvt = zfvj* ( ptn(ji,jj,jk) + ptn(ji ,jj+1,jk) )138 zcenut = pun(ji,jj,jk) * ( ptn(ji,jj,jk) + ptn(ji+1,jj ,jk) ) 139 zcenvt = pvn(ji,jj,jk) * ( ptn(ji,jj,jk) + ptn(ji ,jj+1,jk) ) 182 140 ! mixed centered / upstream scheme 183 zwx(ji,jj,jk) = zcofi * zupsut + (1.-zcofi) * zcenut 184 zwy(ji,jj,jk) = zcofj * zupsvt + (1.-zcofj) * zcenvt 185 END DO 186 END DO 187 188 DO jj = 2, jpjm1 ! horizontal tracer flux divergence added to the general trend 189 DO ji = fs_2, fs_jpim1 ! vector opt. 190 #if defined key_zco 191 zbtr = btr2(ji,jj) 192 #else 193 zbtr = btr2(ji,jj) / fse3t(ji,jj,jk) 194 #endif 195 ! horizontal advective trends 196 zta = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk) & 197 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk) ) 198 ! add it to the general tracer trends 199 pta(ji,jj,jk) = pta(ji,jj,jk) + zta 200 END DO 201 END DO 202 ! ! =============== 203 END DO ! End of slab 204 ! ! =============== 205 206 ! Save the horizontal advective trends for diagnostic 207 ! ----------------------------------------------------- 208 IF( l_trdtra ) THEN 209 CALL trd_tra_adv( kt, ktra, jpt_trd_xad, cdtype, zwx, pun, ptn ) 210 CALL trd_tra_adv( kt, ktra, jpt_trd_yad, cdtype, zwy, pvn, ptn ) 211 ENDIF 212 213 IF(ln_ctl) CALL prt_ctl( tab3d_1=pta, clinfo1=' cen2 - had: ', mask1=tmask, clinfo3=cdtype ) 214 215 216 ! "Poleward" heat and salt transport 217 ! ---------------------------------- 218 IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN 219 IF( lk_zco ) THEN 220 DO jk = 1, jpkm1 221 DO jj = 2, jpjm1 222 DO ji = fs_2, fs_jpim1 ! vector opt. 223 zwy(ji,jj,jk) = zwy(ji,jj,jk) * fse3v(ji,jj,jk) 224 END DO 225 END DO 226 END DO 227 ENDIF 228 IF( ktra == jp_tem) pht_adv(:) = ptr_vj( zwy(:,:,:) ) 229 IF( ktra == jp_sal) pst_adv(:) = ptr_vj( zwy(:,:,:) ) 230 ENDIF 231 232 233 ! II. Vertical advection 234 ! ---------------------- 235 236 IF( lk_dynspg_rl .OR. lk_vvl ) THEN ! rigid lid or non-linear free surface 237 zwx(:,:, 1 ) = 0.e0 ! Surface value : zero flux 238 zwx(:,:,jpk) = 0.e0 ! Bottom value : flux set to zero 239 ELSE ! linear free surface 240 zwx(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1) ! Surface : : advection through z=0 241 zwx(:,:,jpk) = 0.e0 ! Bottom : flux set to zero 242 ENDIF 243 244 DO jk = 2, jpk ! Vertical advective fluxes (at w-point) 141 zwx(ji,jj,jk) = 0.5 * ( zcofi * zupsut + (1.-zcofi) * zcenut ) 142 zwy(ji,jj,jk) = 0.5 * ( zcofj * zupsvt + (1.-zcofj) * zcenvt ) 143 END DO 144 END DO 145 END DO 146 ! ! Vertical advective fluxes 147 zwz(:,:,jpk) = 0.e0 ! Bottom value : flux set to zero 148 ! ! Surface value : 149 IF( lk_dynspg_rl .OR. lk_vvl ) THEN ; zwz(:,:, 1 ) = 0.e0 ! rigid lid or non-linear fs 150 ELSE ; zwz(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1) ! linear free surface 151 ENDIF 152 ! 153 DO jk = 2, jpkm1 ! interior values 245 154 DO jj = 2, jpjm1 246 155 DO ji = fs_2, fs_jpim1 ! vector opt. 247 156 ! upstream indicator 248 157 zcofk = MAX( zind(ji,jj,jk-1), zind(ji,jj,jk) ) 249 ! velocity * 1/2250 zhw = 0.5 * pwn(ji,jj,jk)251 158 ! upstream scheme 252 zfp_w = zhw + ABS( zhw)253 zfm_w = zhw - ABS( zhw)159 zfp_w = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) ) 160 zfm_w = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) ) 254 161 zupst = zfp_w * ptb(ji,jj,jk) + zfm_w * ptb(ji,jj,jk-1) 255 162 ! centered scheme 256 zcent = zhw* ( ptn(ji,jj,jk) + ptn(ji,jj,jk-1) )163 zcent = pwn(ji,jj,jk) * ( ptn(ji,jj,jk) + ptn(ji,jj,jk-1) ) 257 164 ! mixed centered / upstream scheme 258 zwx(ji,jj,jk) = zcofk * zupst + (1.-zcofk) * zcent 259 END DO 260 END DO 261 END DO 262 263 DO jk = 1, jpkm1 ! Tracer flux divergence at t-point added to the general trend 264 DO jj = 2, jpjm1 165 zwz(ji,jj,jk) = 0.5 * ( zcofk * zupst + (1.-zcofk) * zcent ) 166 END DO 167 END DO 168 END DO 169 170 ! II. Divergence of advective fluxes 171 ! ---------------------------------- 172 DO jk = 1, jpkm1 173 DO jj = 2, jpjm1 265 174 DO ji = fs_2, fs_jpim1 ! vector opt. 266 ze3tr = 1. / fse3t(ji,jj,jk) 267 ! vertical advective trends 268 zta = - ze3tr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) 269 ! add it to the general tracer trends 270 pta(ji,jj,jk) = pta(ji,jj,jk) + zta 271 END DO 272 END DO 273 END DO 274 275 ! 3. Save the vertical advective trends for diagnostic 276 ! ---------------------------------------------------- 277 IF( l_trdtra ) CALL trd_tra_adv( kt, ktra, jpt_trd_zad, cdtype, zwx, pwn, ptn ) 278 279 IF(ln_ctl) CALL prt_ctl( tab3d_1=pta, clinfo1=' cen2 - zad : ', mask1=tmask, clinfo3=cdtype ) 175 zbtr = 1.e0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 176 ! advective trends added to the general tracer trends 177 pta(ji,jj,jk) = pta(ji,jj,jk) - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 178 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & 179 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) 180 END DO 181 END DO 182 END DO 183 184 185 ! III. Diagnostics and control print 186 ! ---------------------------------- 187 IF( l_trdtra ) THEN ! trend diagnostics 188 CALL trd_tra_adv( kt, ktra, jpt_trd_xad, cdtype, zwx, pun, ptn ) 189 CALL trd_tra_adv( kt, ktra, jpt_trd_yad, cdtype, zwy, pvn, ptn ) 190 CALL trd_tra_adv( kt, ktra, jpt_trd_zad, cdtype, zwz, pwn, ptn ) 191 ENDIF 192 ! ! "Poleward" heat and salt transport 193 IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN 194 IF( ktra == jp_tem) pht_adv(:) = ptr_vj( zwy(:,:,:) ) 195 IF( ktra == jp_sal) pst_adv(:) = ptr_vj( zwy(:,:,:) ) 196 ENDIF 197 ! ! control print 198 IF(ln_ctl) CALL prt_ctl( tab3d_1=pta, clinfo1=' cen2 - adv: ', mask1=tmask, clinfo3=cdtype ) 280 199 ! 281 200 END SUBROUTINE tra_adv_cen2
Note: See TracChangeset
for help on using the changeset viewer.
