Changeset 791 for branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_muscl2.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_muscl2.F90
r786 r791 4 4 !! Ocean active tracers: horizontal & vertical advective trend 5 5 !!============================================================================== 6 !! History : 1.0 ! 02-06 (G. Madec) from traadv_muscl7 !! 2.4 ! 08-01 (G. Madec) Merge TRA-TRC6 !! History : 1.0 ! 2002-06 (G. Madec) from traadv_muscl 7 !! 2.4 ! 2008-01 (G. Madec) merge TRC-TRA + switch from velocity to transport 8 8 !!---------------------------------------------------------------------- 9 9 … … 26 26 PRIVATE 27 27 28 !! * Accessibility29 28 PUBLIC tra_adv_muscl2 ! routine called by step.F90 30 29 … … 34 33 !!---------------------------------------------------------------------- 35 34 !! NEMO/OPA 2.4 , LOCEAN-IPSL (2008) 36 !! $Id :$35 !! $Id$ 37 36 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 38 37 !!---------------------------------------------------------------------- … … 66 65 !! 67 66 INTEGER :: ji, jj, jk ! dummy loop indices 68 REAL(wp) :: zu, zv, zw , zeu, zev69 REAL(wp) :: z ew, zbtr, zstep, z267 REAL(wp) :: zu, zv, zw 68 REAL(wp) :: zbtr, zstep, z2 70 69 REAL(wp) :: z0u, z0v, z0w 71 REAL(wp) :: zzwx, zzwy, zalpha , zta70 REAL(wp) :: zzwx, zzwy, zalpha 72 71 REAL(wp), DIMENSION (jpi,jpj,jpk) :: zwx, zwy, zslpx, zslpy ! 3D workspace 73 72 !!---------------------------------------------------------------------- … … 98 97 zwx(:,:,jpk) = 0.e0 ; zwy(:,:,jpk) = 0.e0 ! bottom values 99 98 99 !!gm this lbc_lnk can be omitted 100 100 ! lateral boundary conditions on zwx, zwy (changed sign) 101 101 CALL lbc_lnk( zwx, 'U', -1. ) ; CALL lbc_lnk( zwy, 'V', -1. ) … … 117 117 DO jj = 2, jpj 118 118 DO ji = fs_2, jpi ! vector opt. 119 zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) & 120 & * MIN( ABS( zslpx(ji ,jj,jk) ), & 121 & 2.*ABS( zwx (ji-1,jj,jk) ), & 122 & 2.*ABS( zwx (ji ,jj,jk) ) ) 123 zslpy(ji,jj,jk) = SIGN( 1., zslpy(ji,jj,jk) ) & 124 & * MIN( ABS( zslpy(ji,jj ,jk) ), & 125 & 2.*ABS( zwy (ji,jj-1,jk) ), & 126 & 2.*ABS( zwy (ji,jj ,jk) ) ) 119 zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN( ABS( zslpx(ji ,jj,jk) ), & 120 & 2.*ABS( zwx (ji-1,jj,jk) ), & 121 & 2.*ABS( zwx (ji ,jj,jk) ) ) 122 zslpy(ji,jj,jk) = SIGN( 1., zslpy(ji,jj,jk) ) * MIN( ABS( zslpy(ji,jj ,jk) ), & 123 & 2.*ABS( zwy (ji,jj-1,jk) ), & 124 & 2.*ABS( zwy (ji,jj ,jk) ) ) 127 125 END DO 128 126 END DO … … 134 132 DO jj = 2, jpjm1 135 133 DO ji = fs_2, fs_jpim1 ! vector opt. 136 ! volume fluxes 137 #if defined key_zco 138 zeu = e2u(ji,jj) * pun(ji,jj,jk) 139 zev = e1v(ji,jj) * pvn(ji,jj,jk) 140 #else 141 zeu = e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk) 142 zev = e1v(ji,jj) * fse3v(ji,jj,jk) * pvn(ji,jj,jk) 143 #endif 144 ! MUSCL fluxes 145 z0u = SIGN( 0.5, pun(ji,jj,jk) ) 134 z0u = SIGN( 0.5, pun(ji,jj,jk) ) 146 135 zalpha = 0.5 - z0u 147 zu = z0u - 0.5 * pun(ji,jj,jk) * zstep / e1u(ji,jj)148 zzwx = ptb(ji+1,jj,jk) + zu *zslpx(ji+1,jj,jk)149 zzwy = ptb(ji ,jj,jk) + zu *zslpx(ji ,jj,jk)150 zwx(ji,jj,jk) = zeu* ( zalpha * zzwx + (1.-zalpha) * zzwy )136 zu = z0u - 0.5 * pun(ji,jj,jk) * zstep / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) ) 137 zzwx = ptb(ji+1,jj,jk) + zu * zslpx(ji+1,jj,jk) 138 zzwy = ptb(ji ,jj,jk) + zu * zslpx(ji ,jj,jk) 139 zwx(ji,jj,jk) = pun(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 151 140 ! 152 141 z0v = SIGN( 0.5, pvn(ji,jj,jk) ) 153 142 zalpha = 0.5 - z0v 154 zv = z0v - 0.5 * pvn(ji,jj,jk) * zstep / e2v(ji,jj)155 zzwx = ptb(ji,jj+1,jk) + zv *zslpy(ji,jj+1,jk)156 zzwy = ptb(ji,jj ,jk) + zv *zslpy(ji,jj ,jk)157 zwy(ji,jj,jk) = zev* ( zalpha * zzwx + (1.-zalpha) * zzwy )143 zv = z0v - 0.5 * pvn(ji,jj,jk) * zstep / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) ) 144 zzwx = ptb(ji,jj+1,jk) + zv * zslpy(ji,jj+1,jk) 145 zzwy = ptb(ji,jj ,jk) + zv * zslpy(ji,jj ,jk) 146 zwy(ji,jj,jk) = pvn(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 158 147 END DO 159 148 END DO … … 166 155 DO jj = 2, jpjm1 167 156 DO ji = fs_2, fs_jpim1 ! vector opt. 168 #if defined key_zco169 157 IF( umask(ji,jj,jk) == 0. ) THEN 170 158 IF( pun(ji+1,jj,jk) > 0. .AND. ji /= jpi ) THEN 171 zwx(ji+1,jj,jk) = e2u(ji+1,jj) *pun(ji+1,jj,jk) * ( ptn(ji+1,jj,jk) + ptn(ji+2,jj,jk) ) * 0.5159 zwx(ji+1,jj,jk) = pun(ji+1,jj,jk) * ( ptn(ji+1,jj,jk) + ptn(ji+2,jj,jk) ) * 0.5 172 160 ENDIF 173 161 IF( pun(ji-1,jj,jk) < 0. ) THEN 174 zwx(ji-1,jj,jk) = e2u(ji-1,jj) *pun(ji-1,jj,jk) * ( ptn(ji-1,jj,jk) + ptn(ji ,jj,jk) ) * 0.5162 zwx(ji-1,jj,jk) = pun(ji-1,jj,jk) * ( ptn(ji-1,jj,jk) + ptn(ji ,jj,jk) ) * 0.5 175 163 ENDIF 176 164 ENDIF 177 165 IF( vmask(ji,jj,jk) == 0. ) THEN 178 166 IF( pvn(ji,jj+1,jk) > 0. .AND. jj /= jpj ) THEN 179 zwy(ji,jj+1,jk) = e1v(ji,jj+1) *pvn(ji,jj+1,jk) * ( ptn(ji,jj+1,jk) + ptn(ji,jj+2,jk) ) * 0.5167 zwy(ji,jj+1,jk) = pvn(ji,jj+1,jk) * ( ptn(ji,jj+1,jk) + ptn(ji,jj+2,jk) ) * 0.5 180 168 ENDIF 181 169 IF( pvn(ji,jj-1,jk) < 0. ) THEN 182 zwy(ji,jj-1,jk) = e1v(ji,jj-1) *pvn(ji,jj-1,jk) * ( ptn(ji,jj-1,jk) + ptn(ji ,jj,jk) ) * 0.5170 zwy(ji,jj-1,jk) = pvn(ji,jj-1,jk) * ( ptn(ji,jj-1,jk) + ptn(ji ,jj,jk) ) * 0.5 183 171 ENDIF 184 172 ENDIF 185 #else186 IF( umask(ji,jj,jk) == 0. ) THEN187 IF( pun(ji+1,jj,jk) > 0. .AND. ji /= jpi ) THEN188 zwx(ji+1,jj,jk) = e2u(ji+1,jj)* fse3u(ji+1,jj,jk) &189 & * pun(ji+1,jj,jk) * ( ptn(ji+1,jj,jk) + ptn(ji+2,jj,jk) ) * 0.5190 ENDIF191 IF( pun(ji-1,jj,jk) < 0. ) THEN192 zwx(ji-1,jj,jk) = e2u(ji-1,jj)* fse3u(ji-1,jj,jk) &193 & * pun(ji-1,jj,jk) * ( ptn(ji-1,jj,jk) + ptn(ji ,jj,jk) ) * 0.5194 ENDIF195 ENDIF196 IF( vmask(ji,jj,jk) == 0. ) THEN197 IF( pvn(ji,jj+1,jk) > 0. .AND. jj /= jpj ) THEN198 zwy(ji,jj+1,jk) = e1v(ji,jj+1) * fse3v(ji,jj+1,jk) &199 & * pvn(ji,jj+1,jk) * ( ptn(ji,jj+1,jk) + ptn(ji,jj+2,jk) ) * 0.5200 ENDIF201 IF( pvn(ji,jj-1,jk) < 0. ) THEN202 zwy(ji,jj-1,jk) = e1v(ji,jj-1)* fse3v(ji,jj-1,jk) &203 & * pvn(ji,jj-1,jk) * ( ptn(ji,jj-1,jk) + ptn(ji ,jj,jk) ) * 0.5204 ENDIF205 ENDIF206 #endif207 173 END DO 208 174 END DO … … 217 183 DO jj = 2, jpjm1 218 184 DO ji = fs_2, fs_jpim1 ! vector opt. 219 #if defined key_zco220 zbtr = 1. / ( e1t(ji,jj)*e2t(ji,jj) )221 #else222 185 zbtr = 1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) ) 223 #endif 224 ! horizontal advective trends 225 zta = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 226 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) ) 227 ! add it to the general tracer trends 228 pta(ji,jj,jk) = pta(ji,jj,jk) + zta 186 ! horizontal advective trends added to the general tracer trends 187 pta(ji,jj,jk) = pta(ji,jj,jk) - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 188 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) ) 229 189 END DO 230 190 END DO … … 270 230 DO ji = 1, jpi 271 231 zslpx(ji,jj,jk) = ( zwx(ji,jj,jk) + zwx(ji,jj,jk+1) ) & 272 & * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji,jj,jk+1) ) )232 & * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji,jj,jk+1) ) ) 273 233 END DO 274 234 END DO … … 279 239 DO jj = 1, jpj 280 240 DO ji = 1, jpi 281 zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) & 282 & * MIN( ABS( zslpx(ji,jj,jk ) ), & 283 & 2.*ABS( zwx (ji,jj,jk+1) ), & 284 & 2.*ABS( zwx (ji,jj,jk ) ) ) 241 zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN( ABS( zslpx(ji,jj,jk ) ), & 242 & 2.*ABS( zwx (ji,jj,jk+1) ), & 243 & 2.*ABS( zwx (ji,jj,jk ) ) ) 285 244 END DO 286 245 END DO … … 293 252 DO jj = 2, jpjm1 294 253 DO ji = fs_2, fs_jpim1 ! vector opt. 295 zew = pwn(ji,jj,jk+1)296 254 z0w = SIGN( 0.5, pwn(ji,jj,jk+1) ) 297 255 zalpha = 0.5 + z0w 298 zw = z0w - 0.5 * pwn(ji,jj,jk+1)*zstep / fse3w(ji,jj,jk+1)299 zzwx = ptb(ji,jj,jk+1) + zw *zslpx(ji,jj,jk+1)300 zzwy = ptb(ji,jj,jk ) + zw *zslpx(ji,jj,jk )301 zwx(ji,jj,jk+1) = zew* ( zalpha * zzwx + (1.-zalpha)*zzwy )256 zw = z0w - 0.5 * pwn(ji,jj,jk+1)*zstep / ( e1t(ji,jj) *e2t(ji,jj) * fse3w(ji,jj,jk+1) ) 257 zzwx = ptb(ji,jj,jk+1) + zw * zslpx(ji,jj,jk+1) 258 zzwy = ptb(ji,jj,jk ) + zw * zslpx(ji,jj,jk ) 259 zwx(ji,jj,jk+1) = pwn(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha)*zzwy ) 302 260 END DO 303 261 END DO … … 315 273 END DO 316 274 317 ! surface values 318 IF( lk_dynspg_rl .OR. lk_vvl ) THEN ! rigid lid or variable volume: flux set to zero 319 zwx(:,:, 1 ) = 0.e0 ! surface 320 zwx(:,:,jpk) = 0.e0 ! bottom 321 ELSE ! free surface 322 zwx(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1) ! surface 323 zwx(:,:,jpk) = 0.e0 ! bottom 324 ENDIF 275 zwx(:,:,jpk) = 0.e0 ! bottom value 276 ! ! surface values 277 IF( lk_dynspg_rl .OR. lk_vvl) THEN ; zwx(:,:, 1 ) = 0.e0 ! rigid lid or variable volume 278 ELSE ; zwx(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1) ! linear free surface 279 ENDIF 325 280 326 281 … … 330 285 DO ji = fs_2, fs_jpim1 ! vector opt. 331 286 zbtr = 1. / fse3t(ji,jj,jk) 332 ! horizontal advective trends 333 zta = - zbtr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) 334 ! add it to the general tracer trends 335 pta(ji,jj,jk) = pta(ji,jj,jk) + zta 287 ! horizontal advective trends added to the general tracer trends 288 pta(ji,jj,jk) = pta(ji,jj,jk) - zbtr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) 336 289 END DO 337 290 END DO
Note: See TracChangeset
for help on using the changeset viewer.
