Changeset 791 for branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_muscl.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_muscl.F90
r786 r791 4 4 !! Ocean active tracers: horizontal & vertical advective trend 5 5 !!====================================================================== 6 !! History : !06-00 (A.Estublier) for passive tracers7 !! !01-08 (E.Durand, G.Madec) adapted for T & S8 !! NEMO 1.0 ! 02-06 (G. Madec) F90: Free form and module9 !! 2.4 ! 08-01 (G. Madec) Merge TRA-TRC6 !! History : OPA ! 2006-00 (A.Estublier) for passive tracers 7 !! 8.2 ! 2001-08 (E.Durand, G.Madec) adapted for T & S 8 !! NEMO 1.0 ! 2002-06 (G. Madec) F90: Free form and module 9 !! 2.4 ! 2008-01 (G. Madec) merge TRC-TRA + switch from velocity to transport 10 10 !!---------------------------------------------------------------------- 11 11 … … 28 28 PRIVATE 29 29 30 PUBLIC tra_adv_muscl ! routine called by step.F9030 PUBLIC tra_adv_muscl ! routine called by step.F90 31 31 32 32 !! * Substitutions … … 35 35 !!---------------------------------------------------------------------- 36 36 !! NEMO/OPA 2.4 , LOCEAN-IPSL (2008) 37 !! $Id :$37 !! $Id$ 38 38 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 39 39 !!---------------------------------------------------------------------- … … 66 66 !! 67 67 INTEGER :: ji, jj, jk ! dummy loop indices 68 REAL(wp) :: zu, zv, zw, zeu, zev 69 REAL(wp) :: zew, zbtr, z2, zstep 70 REAL(wp) :: z0u, z0v, z0w 71 REAL(wp) :: zzwx, zzwy, zalpha 72 REAL(wp) :: zta 68 REAL(wp) :: zu, z0u, zzwx 69 REAL(wp) :: zv, z0v, zzwy 70 REAL(wp) :: zw, z0w 71 REAL(wp) :: zbtr, z2, zdt, zalpha 73 72 REAL(wp), DIMENSION (jpi,jpj,jpk) :: zwx, zwy, zslpx, zslpy ! 3D workspace 74 73 !!---------------------------------------------------------------------- … … 80 79 ENDIF 81 80 82 IF( neuler == 0 .AND. kt == nit000 ) THEN ; z2 = 1. 83 ELSE ; z2 = 2. 81 IF( neuler == 0 .AND. kt == nit000 ) THEN ; z2 = 1. ! euler time-stepping 82 ELSE ; z2 = 2. ! leap-frog time-stepping 84 83 ENDIF 85 84 86 85 ! I. Horizontal advective fluxes 87 86 ! ------------------------------ 88 ! first guess of the slopes89 ! interiorvalues90 DO jk = 1, jpkm1 87 ! !-- first guess of the slopes 88 zwx(:,:,jpk) = 0.e0 ; zwy(:,:,jpk) = 0.e0 ! bottom values 89 DO jk = 1, jpkm1 ! interior values 91 90 DO jj = 1, jpjm1 92 91 DO ji = 1, fs_jpim1 ! vector opt. … … 96 95 END DO 97 96 END DO 98 ! bottom values 99 zwx(:,:,jpk) = 0.e0 ; zwy(:,:,jpk) = 0.e0 100 101 ! lateral boundary conditions on zwx, zwy (changed sign) 102 CALL lbc_lnk( zwx, 'U', -1. ) ; CALL lbc_lnk( zwy, 'V', -1. ) 103 104 ! Slopes 105 ! interior values 106 DO jk = 1, jpkm1 107 DO jj = 2, jpj 108 DO ji = fs_2, jpi ! vector opt. 97 !!gm optimisation (lbc to be suppressed) 98 CALL lbc_lnk( zwx, 'U', -1. ) ! lateral boundary conditions on zwx, zwy (changed sign) 99 CALL lbc_lnk( zwy, 'V', -1. ) 100 !!gm 101 102 ! !-- Slopes of tracer 103 zslpx(:,:,jpk) = 0.e0 ; zslpy(:,:,jpk) = 0.e0 ! bottom values 104 DO jk = 1, jpkm1 ! interior values 105 DO jj = 2, jpjm1 106 DO ji = fs_2, jpim1 ! vector opt. 109 107 zslpx(ji,jj,jk) = ( zwx(ji,jj,jk) + zwx(ji-1,jj ,jk) ) & 110 & * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji-1,jj ,jk) ) )108 & * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji-1,jj ,jk) ) ) 111 109 zslpy(ji,jj,jk) = ( zwy(ji,jj,jk) + zwy(ji ,jj-1,jk) ) & 112 & * ( 0.25 + SIGN( 0.25, zwy(ji,jj,jk) * zwy(ji ,jj-1,jk) ) ) 113 END DO 114 END DO 115 END DO 116 ! bottom values 117 zslpx(:,:,jpk) = 0.e0 ; zslpy(:,:,jpk) = 0.e0 118 119 ! Slopes limitation 120 DO jk = 1, jpkm1 121 DO jj = 2, jpj 122 DO ji = fs_2, jpi ! vector opt. 123 zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) & 124 & * MIN( ABS( zslpx(ji ,jj,jk) ), & 125 & 2.*ABS( zwx (ji-1,jj,jk) ), & 126 & 2.*ABS( zwx (ji ,jj,jk) ) ) 127 zslpy(ji,jj,jk) = SIGN( 1., zslpy(ji,jj,jk) ) & 128 & * MIN( ABS( zslpy(ji,jj ,jk) ), & 129 & 2.*ABS( zwy (ji,jj-1,jk) ), & 130 & 2.*ABS( zwy (ji,jj ,jk) ) ) 110 & * ( 0.25 + SIGN( 0.25, zwy(ji,jj,jk) * zwy(ji ,jj-1,jk) ) ) 111 END DO 112 END DO 113 END DO 114 !!gm :merge the 2 loop (above & below) 115 DO jk = 1, jpkm1 ! Slopes limitation 116 DO jj = 2, jpjm1 117 DO ji = fs_2, jpim1 ! vector opt. 118 zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN( ABS( zslpx(ji ,jj,jk) ), & 119 & 2.*ABS( zwx (ji-1,jj,jk) ), & 120 & 2.*ABS( zwx (ji ,jj,jk) ) ) 121 zslpy(ji,jj,jk) = SIGN( 1., zslpy(ji,jj,jk) ) * MIN( ABS( zslpy(ji,jj ,jk) ), & 122 & 2.*ABS( zwy (ji,jj-1,jk) ), & 123 & 2.*ABS( zwy (ji,jj ,jk) ) ) 131 124 END DO 132 125 END DO 133 126 END DO 134 135 ! Advection terms 136 ! interior values 137 DO jk = 1, jpkm1 138 zstep = z2 * rdttra(jk) 139 DO jj = 2, jpjm1 140 DO ji = fs_2, fs_jpim1 ! vector opt. 141 ! volume fluxes 142 #if defined key_zco 143 zeu = e2u(ji,jj) * pun(ji,jj,jk) 144 zev = e1v(ji,jj) * pvn(ji,jj,jk) 145 #else 146 zeu = e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk) 147 zev = e1v(ji,jj) * fse3v(ji,jj,jk) * pvn(ji,jj,jk) 148 #endif 149 ! MUSCL fluxes 127 !!gm optimisation: call lbclnk just 1 time here on zslpx & zslpy 128 !! CALL lbc_lnk( zslpx, 'U', -1. ) ! lateral boundary conditions on zwx, zwy (changed sign) 129 !! CALL lbc_lnk( zslpy, 'V', -1. ) 130 !! 131 !! and loop below from 1 to jpjm1 and to jpim1 132 !!gm 133 134 ! !-- MUSCL horizontal advective fluxes 135 DO jk = 1, jpkm1 ! interior values 136 zdt = z2 * rdttra(jk) 137 DO jj = 2, jpjm1 138 DO ji = fs_2, fs_jpim1 ! vector opt. 150 139 z0u = SIGN( 0.5, pun(ji,jj,jk) ) 151 140 zalpha = 0.5 - z0u 152 zu = z0u - 0.5 * pun(ji,jj,jk) * z step / e1u(ji,jj)141 zu = z0u - 0.5 * pun(ji,jj,jk) * zdt / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) ) 153 142 zzwx = ptb(ji+1,jj,jk) + zu * zslpx(ji+1,jj,jk) 154 143 zzwy = ptb(ji ,jj,jk) + zu * zslpx(ji ,jj,jk) 155 zwx(ji,jj,jk) = zeu* ( zalpha * zzwx + (1.-zalpha) * zzwy )144 zwx(ji,jj,jk) = pun(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 156 145 ! 157 146 z0v = SIGN( 0.5, pvn(ji,jj,jk) ) 158 147 zalpha = 0.5 - z0v 159 zv = z0v - 0.5 * pvn(ji,jj,jk) * z step / e2v(ji,jj)148 zv = z0v - 0.5 * pvn(ji,jj,jk) * zdt / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) ) 160 149 zzwx = ptb(ji,jj+1,jk) + zv * zslpy(ji,jj+1,jk) 161 150 zzwy = ptb(ji,jj ,jk) + zv * zslpy(ji,jj ,jk) 162 zwy(ji,jj,jk) = zev * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 163 END DO 164 END DO 165 END DO 166 167 !!gm bug? there is too many lbc: this have to be checked 168 ! lateral boundary conditions on zwx, zwy (changed sign) 151 zwy(ji,jj,jk) = pvn(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 152 END DO 153 END DO 154 END DO 155 !!gm optimisation (lbc to be suppressed) 156 ! ! lateral boundary conditions on zwx, zwy (changed sign) 169 157 CALL lbc_lnk( zwx, 'U', -1. ) ; CALL lbc_lnk( zwy, 'V', -1. ) 170 158 !!gm 159 160 ! !-- MUSCL advective trend 171 161 ! Tracer flux divergence at t-point added to the general trend 172 162 DO jk = 1, jpkm1 173 163 DO jj = 2, jpjm1 174 164 DO ji = fs_2, fs_jpim1 ! vector opt. 175 #if defined key_zco176 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) )177 #else178 165 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 179 #endif 180 ! horizontal advective trends 181 zta = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 182 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) ) 183 ! add it to the general tracer trends 184 pta(ji,jj,jk) = pta(ji,jj,jk) + zta 166 ! horizontal advective trends added to the general tracer trends 167 pta(ji,jj,jk) = pta(ji,jj,jk) - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 168 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) ) 185 169 END DO 186 170 END DO 187 171 END DO 188 172 189 IF(ln_ctl) CALL prt_ctl( tab3d_1=pta, clinfo1=' muscl - had: ', mask1=tmask, clinfo3=cdtype ) 190 191 ! Save the horizontal advective trends for diagnostics 192 IF( l_trdtra ) THEN 173 IF( l_trdtra ) THEN !-- trend diagnostics 193 174 CALL trd_tra_adv( kt, ktra, jpt_trd_xad, cdtype, zwx, pun, ptb ) 194 175 CALL trd_tra_adv( kt, ktra, jpt_trd_yad, cdtype, zwy, pvn, ptb ) 195 176 ENDIF 196 177 ! !-- "Poleward" heat or salt transports 197 178 IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN 198 IF( lk_zco ) THEN199 DO jk = 1, jpkm1200 DO jj = 2, jpjm1201 DO ji = fs_2, fs_jpim1 ! vector opt.202 zwy(ji,jj,jk) = zwy(ji,jj,jk) * fse3v(ji,jj,jk)203 END DO204 END DO205 END DO206 ENDIF207 179 IF( ktra == jp_tem) pht_adv(:) = ptr_vj( zwy(:,:,:) ) 208 180 IF( ktra == jp_sal) pst_adv(:) = ptr_vj( zwy(:,:,:) ) 209 181 ENDIF 182 ! !-- control print 183 IF(ln_ctl) CALL prt_ctl( tab3d_1=pta, clinfo1=' muscl - had: ', mask1=tmask, clinfo3=cdtype ) 210 184 211 185 … … 213 187 ! ----------------------------- 214 188 215 ! First guess of the slope216 ! interior values217 DO jk = 2, jpkm1 189 ! !-- first guess of the slopes 190 zwx (:,:, 1 ) = 0.e0 ; zwx (:,:,jpk) = 0.e0 ! surface & bottom boundary conditions 191 DO jk = 2, jpkm1 ! interior values 218 192 zwx(:,:,jk) = tmask(:,:,jk) * ( ptb(:,:,jk-1) - ptb(:,:,jk) ) 219 193 END DO 220 ! surface & bottom boundary conditions 221 zwx (:,:, 1 ) = 0.e0 ; zwx (:,:,jpk) = 0.e0 222 223 ! Slopes 224 DO jk = 2, jpkm1 194 195 ! !-- Slopes of tracer 196 zslpx(:,:,1) = 0.e0 ! surface values 197 DO jk = 2, jpkm1 ! interior value 225 198 DO jj = 1, jpj 226 199 DO ji = 1, jpi … … 230 203 END DO 231 204 END DO 232 233 ! Slopes limitation234 DO jk = 2, jpkm1 ! interior values205 !!gm : merge the above and below loops 206 ! !-- Slopes limitation 207 DO jk = 2, jpkm1 ! interior values 235 208 DO jj = 1, jpj 236 209 DO ji = 1, jpi 237 zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) & 238 & * MIN( ABS( zslpx(ji,jj,jk ) ), & 239 & 2.*ABS( zwx (ji,jj,jk+1) ), & 240 & 2.*ABS( zwx (ji,jj,jk ) ) ) 241 END DO 242 END DO 243 END DO 244 zslpx(:,:,1) = 0.e0 ! surface values 245 246 ! vertical advective flux 247 DO jk = 1, jpkm1 ! interior values 248 zstep = z2 * rdttra(jk) 249 DO jj = 2, jpjm1 250 DO ji = fs_2, fs_jpim1 ! vector opt. 251 zew = pwn(ji,jj,jk+1) 210 zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN( ABS( zslpx(ji,jj,jk ) ), & 211 & 2.*ABS( zwx (ji,jj,jk+1) ), & 212 & 2.*ABS( zwx (ji,jj,jk ) ) ) 213 END DO 214 END DO 215 END DO 216 217 ! !-- vertical advective flux 218 ! ! surface values (bottom already set to zero) 219 IF( lk_dynspg_rl .OR. lk_vvl) THEN ; zwx(:,:, 1 ) = 0.e0 ! rigid lid or variable volume 220 ELSE ; zwx(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1) ! linear free surface 221 ENDIF 222 DO jk = 1, jpkm1 ! interior values 223 zdt = z2 * rdttra(jk) 224 DO jj = 2, jpjm1 225 DO ji = fs_2, fs_jpim1 ! vector opt. 252 226 z0w = SIGN( 0.5, pwn(ji,jj,jk+1) ) 253 227 zalpha = 0.5 + z0w 254 zw = z0w - 0.5 * pwn(ji,jj,jk+1) * z step / fse3w(ji,jj,jk+1)228 zw = z0w - 0.5 * pwn(ji,jj,jk+1) * zdt / (e1t(ji,jj) * e2t(ji,jj) * fse3w(ji,jj,jk+1) ) 255 229 zzwx = ptb(ji,jj,jk+1) + zw * zslpx(ji,jj,jk+1) 256 230 zzwy = ptb(ji,jj,jk ) + zw * zslpx(ji,jj,jk ) 257 zwx(ji,jj,jk+1) = zew * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 258 END DO 259 END DO 260 END DO 261 ! ! surface values 262 IF( lk_dynspg_rl .OR. lk_vvl) THEN ! rigid lid or variable volume: flux set to zero 263 zwx(:,:, 1 ) = 0.e0 ! surface 264 zwx(:,:,jpk) = 0.e0 ! bottom 265 ELSE ! free surface 266 zwx(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1) ! Surface 267 zwx(:,:,jpk) = 0.e0 ! bottom 268 269 ENDIF 270 271 ! Compute & add the vertical advective trend 231 zwx(ji,jj,jk+1) = pwn(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 232 END DO 233 END DO 234 END DO 235 236 ! !-- vertical advective trend 272 237 DO jk = 1, jpkm1 273 238 DO jj = 2, jpjm1 274 239 DO ji = fs_2, fs_jpim1 ! vector opt. 275 zbtr = 1. / fse3t(ji,jj,jk) 276 ! horizontal advective trends 277 zta = - zbtr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) 278 ! add it to the general tracer trends 279 pta(ji,jj,jk) = pta(ji,jj,jk) + zta 280 END DO 281 END DO 282 END DO 283 284 ! Save the vertical advective trends for diagnostic 285 ! ------------------------------------------------- 240 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 241 ! horizontal advective trends added to the general tracer trends 242 pta(ji,jj,jk) = pta(ji,jj,jk) - zbtr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) 243 END DO 244 END DO 245 END DO 246 247 ! !-- trend diagnostic 286 248 IF( l_trdtra ) CALL trd_tra_adv( kt, ktra, jpt_trd_zad, cdtype, zwx, pwn, ptb ) 287 249 ! !-- control print 288 250 IF(ln_ctl) CALL prt_ctl( tab3d_1=pta, clinfo1=' muscl - zad: ', mask1=tmask, clinfo3=cdtype ) 289 251 !
Note: See TracChangeset
for help on using the changeset viewer.
