Changeset 503 for trunk/NEMO/OPA_SRC/TRA/traadv_muscl.F90
- Timestamp:
- 2006-09-27T10:52:29+02:00 (18 years ago)
- File:
-
- 1 edited
-
trunk/NEMO/OPA_SRC/TRA/traadv_muscl.F90 (modified) (11 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/OPA_SRC/TRA/traadv_muscl.F90
r457 r503 1 1 MODULE traadv_muscl 2 !!====================================================================== ========2 !!====================================================================== 3 3 !! *** MODULE traadv_muscl *** 4 4 !! Ocean active tracers: horizontal & vertical advective trend 5 !!============================================================================== 5 !!====================================================================== 6 !! History : ! 06-00 (A.Estublier) for passive tracers 7 !! ! 01-08 (E.Durand, G.Madec) adapted for T & S 8 !! 8.5 ! 02-06 (G. Madec) F90: Free form and module 9 !!---------------------------------------------------------------------- 6 10 7 11 !!---------------------------------------------------------------------- … … 9 13 !! and vertical advection trends using MUSCL scheme 10 14 !!---------------------------------------------------------------------- 11 !! * Modules used12 15 USE oce ! ocean dynamics and active tracers 13 16 USE dom_oce ! ocean space and time domain … … 25 28 PRIVATE 26 29 27 !! * Accessibility 28 PUBLIC tra_adv_muscl ! routine called by step.F90 30 PUBLIC tra_adv_muscl ! routine called by step.F90 29 31 30 32 !! * Substitutions … … 32 34 # include "vectopt_loop_substitute.h90" 33 35 !!---------------------------------------------------------------------- 34 !! OPA 9.0 , LOCEAN-IPSL (200 5)36 !! OPA 9.0 , LOCEAN-IPSL (2006) 35 37 !! $Header$ 36 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt38 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 37 39 !!---------------------------------------------------------------------- 38 40 … … 50 52 !! 51 53 !! ** Action : - update (ta,sa) with the now advective tracer trends 52 !! - save trends in (ztdta,ztdsa) ('key_trdtra') 53 !! 54 !! References : 55 !! Estubier, A., and M. Levy, Notes Techn. Pole de Modelisation 56 !! IPSL, Sept. 2000 (http://www.lodyc.jussieu.fr/opa) 57 !! 58 !! History : 59 !! ! 06-00 (A.Estublier) for passive tracers 60 !! ! 01-08 (E.Durand G.Madec) adapted for T & S 61 !! 8.5 ! 02-06 (G. Madec) F90: Free form and module 62 !! 9.0 ! 04-08 (C. Talandier) New trends organization 54 !! - save trends in (ztrdt,ztrds) ('key_trdtra') 55 !! 56 !! References : Estubier, A., and M. Levy, Notes Techn. Pole de Modelisation 57 !! IPSL, Sept. 2000 (http://www.lodyc.jussieu.fr/opa) 63 58 !!---------------------------------------------------------------------- 64 !! * Arguments 65 INTEGER , INTENT( in ) :: kt ! ocean time-step index 66 REAL(wp), INTENT( in ), DIMENSION(jpi,jpj,jpk) :: & 67 pun, pvn, pwn ! now ocean velocity fields 68 69 !! * Local declarations 70 INTEGER :: ji, jj, jk ! dummy loop indices 59 USE oce, ONLY : ztrdt => ua ! use ua as workspace 60 USE oce, ONLY : ztrds => va ! use va as workspace 61 !! 62 INTEGER , INTENT(in) :: kt ! ocean time-step index 63 REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pun ! ocean velocity u-component 64 REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pvn ! ocean velocity v-component 65 REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pwn ! ocean velocity w-component 66 !! 67 INTEGER :: ji, jj, jk ! dummy loop indices 71 68 REAL(wp) :: & 72 69 zu, zv, zw, zeu, zev, & … … 75 72 zzt1, zzt2, zalpha, & 76 73 zzs1, zzs2, z2, & 77 zta, zsa 78 REAL(wp), DIMENSION (jpi,jpj,jpk) :: & 79 zt1, zt2, ztp1, ztp2, & 80 zs1, zs2, zsp1, zsp2, & 81 ztdta, ztdsa 74 zta, zsa, & 75 z_hdivn_x, z_hdivn_y, z_hdivn 76 REAL(wp), DIMENSION (jpi,jpj,jpk) :: zt1, zt2, ztp1, ztp2 ! 3D workspace 77 REAL(wp), DIMENSION (jpi,jpj,jpk) :: zs1, zs2, zsp1, zsp2 ! " " 82 78 !!---------------------------------------------------------------------- 83 79 … … 88 84 ENDIF 89 85 90 IF( neuler == 0 .AND. kt == nit000 ) THEN 91 z2=1. 92 ELSE 93 z2=2. 94 ENDIF 95 96 ! Save ta and sa trends 97 IF( l_trdtra ) THEN 98 ztdta(:,:,:) = ta(:,:,:) 99 ztdsa(:,:,:) = sa(:,:,:) 100 l_adv = 'mus' 101 ENDIF 102 86 IF( neuler == 0 .AND. kt == nit000 ) THEN ; z2 = 1. 87 ELSE ; z2 = 2. 88 ENDIF 103 89 104 90 ! I. Horizontal advective fluxes 105 91 ! ------------------------------ 106 107 92 ! first guess of the slopes 108 93 ! interior values … … 193 178 zt1(ji,jj,jk) = zeu * ( zalpha * zzt1 + (1.-zalpha) * zzt2 ) 194 179 zs1(ji,jj,jk) = zeu * ( zalpha * zzs1 + (1.-zalpha) * zzs2 ) 195 180 ! 196 181 z0v = SIGN( 0.5, pvn(ji,jj,jk) ) 197 182 zalpha = 0.5 - z0v … … 211 196 CALL lbc_lnk( zt2, 'V', -1. ) ; CALL lbc_lnk( zs2, 'V', -1. ) 212 197 213 ! Save MUSCL fluxes to compute i- & j- horizontal 214 ! advection trends in the MLD 215 IF( l_trdtra ) THEN 216 ! save i- terms 217 tladi(:,:,:) = zt1(:,:,:) 218 sladi(:,:,:) = zs1(:,:,:) 219 ! save j- terms 220 tladj(:,:,:) = zt2(:,:,:) 221 sladj(:,:,:) = zs2(:,:,:) 222 ENDIF 223 224 ! Compute & add the horizontal advective trend 225 198 ! Tracer flux divergence at t-point added to the general trend 226 199 DO jk = 1, jpkm1 227 200 DO jj = 2, jpjm1 … … 244 217 END DO 245 218 246 ! Save the horizontal advective trends for diagnostic 247 248 IF( l_trdtra ) THEN 249 ! Recompute the horizontal advection zta & zsa trends computed 250 ! at the step 2. above in making the difference between the new 251 ! trends and the previous one ta()/sa - ztdta()/ztdsa() and add 252 ! the term tn()/sn()*hdivn() to recover the Uh gradh(T/S) trends 253 ztdta(:,:,:) = ta(:,:,:) - ztdta(:,:,:) + tn(:,:,:) * hdivn(:,:,:) 254 ztdsa(:,:,:) = sa(:,:,:) - ztdsa(:,:,:) + sn(:,:,:) * hdivn(:,:,:) 255 256 CALL trd_mod(ztdta, ztdsa, jpttdlad, 'TRA', kt) 257 258 ! Save the new ta and sa trends 259 ztdta(:,:,:) = ta(:,:,:) 260 ztdsa(:,:,:) = sa(:,:,:) 261 262 ENDIF 263 264 IF(ln_ctl) THEN 265 CALL prt_ctl(tab3d_1=ta, clinfo1=' muscl had - Ta: ', mask1=tmask , & 266 & tab3d_2=sa, clinfo2=' Sa: ', mask2=tmask, clinfo3='tra' ) 219 IF(ln_ctl) CALL prt_ctl( tab3d_1=ta, clinfo1=' muscl had - Ta: ', mask1=tmask , & 220 & tab3d_2=sa, clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' ) 221 222 ! Save the horizontal advective trends for diagnostics 223 IF( l_trdtra ) THEN 224 ztrdt(:,:,:) = 0.e0 ; ztrds(:,:,:) = 0.e0 225 ! 226 ! T/S ZONAL advection trends 227 DO jk = 1, jpkm1 228 DO jj = 2, jpjm1 229 DO ji = fs_2, fs_jpim1 ! vector opt. 230 !-- Compute zonal divergence by splitting hdivn (see divcur.F90) 231 ! N.B. This computation is not valid along OBCs (if any) 232 #if defined key_zco 233 zbtr = 1. / ( e1t(ji,jj) * e2t(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 = 1. / ( e1t(ji,jj) * e2t(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) = - zbtr * ( zt1(ji,jj,jk) - zt1(ji-1,jj,jk) ) + tn(ji,jj,jk) * z_hdivn_x 242 ztrds(ji,jj,jk) = - zbtr * ( zs1(ji,jj,jk) - zs1(ji-1,jj,jk) ) + 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) 247 248 ! T/S MERIDIONAL advection trends 249 DO jk = 1, jpkm1 250 DO jj = 2, jpjm1 251 DO ji = fs_2, fs_jpim1 ! vector opt. 252 !-- Compute merid. divergence by splitting hdivn (see divcur.F90) 253 ! N.B. This computation is not valid along OBCs (if any) 254 #if defined key_zco 255 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) ) 256 z_hdivn_y = ( e1v(ji,jj ) * pvn(ji,jj ,jk) & 257 & - e1v(ji,jj-1) * pvn(ji,jj-1,jk) ) * zbtr 258 #else 259 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 260 z_hdivn_y = ( e1v(ji, jj) * fse3v(ji,jj ,jk) * pvn(ji,jj ,jk) & 261 & - e1v(ji,jj-1) * fse3v(ji,jj-1,jk) * pvn(ji,jj-1,jk) ) * zbtr 262 #endif 263 ztrdt(ji,jj,jk) = - zbtr * ( zt2(ji,jj,jk) - zt2(ji,jj-1,jk) ) + tn(ji,jj,jk) * z_hdivn_y 264 ztrds(ji,jj,jk) = - zbtr * ( zs2(ji,jj,jk) - zs2(ji,jj-1,jk) ) + sn(ji,jj,jk) * z_hdivn_y 265 END DO 266 END DO 267 END DO 268 CALL trd_mod(ztrdt, ztrds, jptra_trd_yad, 'TRA', kt) 269 270 ! Save the up-to-date ta and sa trends 271 ztrdt(:,:,:) = ta(:,:,:) 272 ztrds(:,:,:) = sa(:,:,:) 273 ! 267 274 ENDIF 268 275 … … 378 385 379 386 ! Save the vertical advective trends for diagnostic 380 387 ! ------------------------------------------------- 381 388 IF( l_trdtra ) THEN 382 389 ! Recompute the vertical advection zta & zsa trends computed 383 390 ! at the step 2. above in making the difference between the new 384 ! trends and the previous one: ta()/sa - zt dta()/ztdsa() and substract391 ! trends and the previous one: ta()/sa - ztrdt()/ztrds() and substract 385 392 ! the term tn()/sn()*hdivn() to recover the W gradz(T/S) trends 386 ztdta(:,:,:) = ta(:,:,:) - ztdta(:,:,:) - tn(:,:,:) * hdivn(:,:,:) 387 ztdsa(:,:,:) = sa(:,:,:) - ztdsa(:,:,:) - sn(:,:,:) * hdivn(:,:,:) 388 389 CALL trd_mod(ztdta, ztdsa, jpttdzad, 'TRA', kt) 390 ENDIF 391 392 IF(ln_ctl) THEN 393 CALL prt_ctl(tab3d_1=ta, clinfo1=' muscl zad - Ta: ', mask1=tmask , & 394 & tab3d_2=sa, clinfo2=' Sa: ', mask2=tmask, clinfo3='tra') 395 ENDIF 396 393 394 DO jk = 1, jpkm1 395 DO jj = 2, jpjm1 396 DO ji = fs_2, fs_jpim1 ! vector opt. 397 #if defined key_zco 398 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) ) 399 z_hdivn_x = e2u(ji,jj)*pun(ji,jj,jk) - e2u(ji-1,jj)*pun(ji-1,jj,jk) 400 z_hdivn_y = e1v(ji,jj)*pvn(ji,jj,jk) - e1v(ji,jj-1)*pvn(ji,jj-1,jk) 401 #else 402 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 403 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) 404 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) 405 #endif 406 z_hdivn = (z_hdivn_x + z_hdivn_y) * zbtr 407 ztrdt(ji,jj,jk) = ta(ji,jj,jk) - ztrdt(ji,jj,jk) - tn(ji,jj,jk) * z_hdivn 408 ztrds(ji,jj,jk) = sa(ji,jj,jk) - ztrds(ji,jj,jk) - sn(ji,jj,jk) * z_hdivn 409 END DO 410 END DO 411 END DO 412 CALL trd_mod(ztrdt, ztrds, jptra_trd_zad, 'TRA', kt) 413 ! 414 ENDIF 415 416 IF(ln_ctl) CALL prt_ctl( tab3d_1=ta, clinfo1=' muscl zad - Ta: ', mask1=tmask , & 417 & tab3d_2=sa, clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' ) 418 ! 397 419 END SUBROUTINE tra_adv_muscl 398 420
Note: See TracChangeset
for help on using the changeset viewer.
