Changeset 503 for trunk/NEMO/OPA_SRC/TRA/traadv_tvd.F90
- Timestamp:
- 2006-09-27T10:52:29+02:00 (18 years ago)
- File:
-
- 1 edited
-
trunk/NEMO/OPA_SRC/TRA/traadv_tvd.F90 (modified) (15 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/OPA_SRC/TRA/traadv_tvd.F90
r457 r503 4 4 !! Ocean active tracers: horizontal & vertical advective trend 5 5 !!============================================================================== 6 !! History : ! 95-12 (L. Mortier) Original code 7 !! ! 00-01 (H. Loukos) adapted to ORCA 8 !! ! 00-10 (MA Foujols E.Kestenare) include file not routine 9 !! ! 00-12 (E. Kestenare M. Levy) fix bug in trtrd indexes 10 !! ! 01-07 (E. Durand G. Madec) adaptation to ORCA config 11 !! 8.5 ! 02-06 (G. Madec) F90: Free form and module 12 !! 9.0 ! 04-01 (A. de Miranda, G. Madec, J.M. Molines ): advective bbl 13 !! 9.0 ! 08-04 (S. Cravatte) add the i-, j- & k- trends computation 14 !! " " ! 05-11 (V. Garnier) Surface pressure gradient organization 15 !!---------------------------------------------------------------------- 16 6 17 7 18 !!---------------------------------------------------------------------- … … 11 22 !! algorithm 12 23 !!---------------------------------------------------------------------- 13 !! * Modules used14 24 USE oce ! ocean dynamics and active tracers 15 25 USE dom_oce ! ocean space and time domain … … 28 38 PRIVATE 29 39 30 !! * Accessibility 31 PUBLIC tra_adv_tvd ! routine called by step.F90 40 PUBLIC tra_adv_tvd ! routine called by step.F90 32 41 33 42 !! * Substitutions … … 35 44 # include "vectopt_loop_substitute.h90" 36 45 !!---------------------------------------------------------------------- 37 !! OPA 9.0 , LOCEAN-IPSL (200 5)38 !! $Header$ 39 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt46 !! OPA 9.0 , LOCEAN-IPSL (2006) 47 !! $Header$ 48 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 40 49 !!---------------------------------------------------------------------- 41 50 … … 54 63 !! 55 64 !! ** Action : - update (ta,sa) with the now advective tracer trends 56 !! - save the trends in (ttrdh,strdh) ('key_trdtra') 65 !! - save the trends in (ztrdt,ztrds) ('key_trdtra') 66 !!---------------------------------------------------------------------- 67 USE oce , ztrdt => ua ! use ua as workspace 68 USE oce , ztrds => va ! use va as workspace 57 69 !! 58 !! History : 59 !! ! 95-12 (L. Mortier) Original code 60 !! ! 00-01 (H. Loukos) adapted to ORCA 61 !! ! 00-10 (MA Foujols E.Kestenare) include file not routine 62 !! ! 00-12 (E. Kestenare M. Levy) fix bug in trtrd indexes 63 !! ! 01-07 (E. Durand G. Madec) adaptation to ORCA config 64 !! 8.5 ! 02-06 (G. Madec) F90: Free form and module 65 !! 9.0 ! 04-01 (A. de Miranda, G. Madec, J.M. Molines ): advective bbl 66 !! 9.0 ! 08-04 (S. Cravatte) add the i-, j- & k- trends computation 67 !! " ! 05-11 (V. Garnier) Surface pressure gradient organization 68 !!---------------------------------------------------------------------- 69 !! * Modules used 70 USE trdmod_oce , ztay => tladj, & ! use tladj latter 71 & zsay => sladj, & ! use sladj latter 72 & ztaz => tladi, & ! use ua as workspace 73 & zsaz => sladi ! use ua as workspace 74 75 !! * Arguments 76 INTEGER , INTENT( in ) :: kt ! ocean time-step index 77 REAL(wp), INTENT( in ), DIMENSION(jpi,jpj,jpk) :: & 78 pun, pvn, pwn ! now ocean velocity fields 79 80 !! * Local declarations 70 INTEGER , INTENT(in) :: kt ! ocean time-step index 71 REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pun ! ocean velocity u-component 72 REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pvn ! ocean velocity v-component 73 REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pwn ! ocean velocity w-component 74 !! 81 75 INTEGER :: ji, jj, jk ! dummy loop indices 82 76 REAL(wp) :: & ! temporary scalar 83 77 ztai, ztaj, ztak, & ! " " 84 zsai, zsaj, zsak ! " " 85 REAL(wp), DIMENSION (jpi,jpj,jpk) :: & 86 zti, ztu, ztv, ztw, & ! temporary workspace 87 zsi, zsu, zsv, zsw, & ! " " 88 ztdta, ztdsa ! " " 78 zsai, zsaj, zsak, & ! " " 79 z_hdivn_x, z_hdivn_y, z_hdivn 89 80 REAL(wp) :: & 90 z2dtt, zbtr, zeu, zev, zew, z2, & ! temporary scalar 81 z2dtt, zbtr, zeu, zev, & ! temporary scalar 82 zew, z2, zbtr1, & ! temporary scalar 91 83 zfp_ui, zfp_vj, zfp_wk, & ! " " 92 84 zfm_ui, zfm_vj, zfm_wk ! " " 85 REAL(wp), DIMENSION (jpi,jpj,jpk) :: zti, ztu, ztv, ztw ! temporary workspace 86 REAL(wp), DIMENSION (jpi,jpj,jpk) :: zsi, zsu, zsv, zsw ! " " 93 87 !!---------------------------------------------------------------------- 94 88 … … 101 95 ENDIF 102 96 103 IF( neuler == 0 .AND. kt == nit000 ) THEN 104 z2=1. 105 ELSE 106 z2=2. 97 IF( neuler == 0 .AND. kt == nit000 ) THEN ; z2 = 1. 98 ELSE ; z2 = 2. 107 99 ENDIF 108 109 ! Save ta and sa trends110 IF( l_trdtra ) THEN111 ztdta(:,:,:) = ta(:,:,:)112 ztdsa(:,:,:) = sa(:,:,:)113 l_adv = 'tvd'114 ENDIF115 116 100 117 101 ! 1. Bottom value : flux set to zero … … 177 161 DO ji = fs_2, fs_jpim1 ! vector opt. 178 162 zbtr = 1./ ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 179 180 163 ! i- j- horizontal & k- vertical advective trends 181 164 ztai = - ( ztu(ji,jj,jk) - ztu(ji-1,jj ,jk ) ) * zbtr … … 185 168 zsaj = - ( zsv(ji,jj,jk) - zsv(ji ,jj-1,jk ) ) * zbtr 186 169 zsak = - ( zsw(ji,jj,jk) - zsw(ji ,jj ,jk+1) ) * zbtr 187 188 170 ! total intermediate advective trends 189 171 zti(ji,jj,jk) = ztai + ztaj + ztak … … 193 175 END DO 194 176 195 ! Save the intermediate vertical & j- horizontal advection trends 196 IF( l_trdtra ) THEN 177 178 ! Save the intermediate i / j / k advective trends for diagnostics 179 ! ------------------------------------------------------------------- 180 ! Warning : We should use zun instead of un in the computations below, but we 181 ! also use hdivn which is computed with un, vn (check ???). So we use un, vn 182 ! for consistency. Results are therefore approximate with key_trabbl_adv. 183 184 IF( l_trdtra ) THEN 185 ztrdt(:,:,:) = 0.e0 ; ztrds(:,:,:) = 0.e0 186 ! 187 ! T/S ZONAL advection trends 197 188 DO jk = 1, jpkm1 198 189 DO jj = 2, jpjm1 199 190 DO ji = fs_2, fs_jpim1 ! vector opt. 200 zbtr = 1./ ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 201 ztay(ji,jj,jk) = - ( ztv(ji,jj,jk) - ztv(ji ,jj-1,jk ) ) * zbtr 202 zsay(ji,jj,jk) = - ( zsv(ji,jj,jk) - zsv(ji ,jj-1,jk ) ) * zbtr 203 ztaz(ji,jj,jk) = - ( ztw(ji,jj,jk) - ztw(ji ,jj ,jk+1) ) * zbtr 204 zsaz(ji,jj,jk) = - ( zsw(ji,jj,jk) - zsw(ji ,jj ,jk+1) ) * zbtr 191 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 192 ztrdt(ji,jj,jk) = - ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) ) * zbtr 193 ztrds(ji,jj,jk) = - ( zsu(ji,jj,jk) - zsu(ji-1,jj,jk) ) * zbtr 205 194 END DO 206 195 END DO 207 196 END DO 197 CALL trd_mod(ztrdt, ztrds, jptra_trd_xad, 'TRA', kt) ! save the trends 198 ! 199 ! T/S MERIDIONAL advection trends 200 DO jk = 1, jpkm1 201 DO jj = 2, jpjm1 202 DO ji = fs_2, fs_jpim1 ! vector opt. 203 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 204 ztrdt(ji,jj,jk) = - ( ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) * zbtr 205 ztrds(ji,jj,jk) = - ( zsv(ji,jj,jk) - zsv(ji,jj-1,jk) ) * zbtr 206 END DO 207 END DO 208 END DO 209 CALL trd_mod(ztrdt, ztrds, jptra_trd_yad, 'TRA', kt) ! save the trends 210 ! 211 ! T/S VERTICAL advection trends 212 DO jk = 1, jpkm1 213 DO jj = 2, jpjm1 214 DO ji = fs_2, fs_jpim1 ! vector opt. 215 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 216 ztrdt(ji,jj,jk) = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * zbtr 217 ztrds(ji,jj,jk) = - ( zsw(ji,jj,jk) - zsw(ji,jj,jk+1) ) * zbtr 218 END DO 219 END DO 220 END DO 221 CALL trd_mod(ztrdt, ztrds, jptra_trd_zad, 'TRA', kt) ! save the trends 222 ! 208 223 ENDIF 209 224 … … 244 259 ! antidiffusive flux on k 245 260 ! Surface value 246 ztw(:,:,1) = 0. 247 zsw(:,:,1) = 0. 261 ztw(:,:,1) = 0.e0 262 zsw(:,:,1) = 0.e0 248 263 249 264 ! Interior value … … 290 305 END DO 291 306 292 ! save the advective trends for diagnostic 293 ! tracers trends 294 IF( l_trdtra ) THEN 295 ! Compute the final vertical & j- horizontal advection trends 307 308 ! Save the advective trends for diagnostics 309 ! -------------------------------------------- 310 311 IF( l_trdtra ) THEN 312 ztrdt(:,:,:) = 0.e0 ; ztrds(:,:,:) = 0.e0 313 ! 314 ! T/S ZONAL advection trends 296 315 DO jk = 1, jpkm1 297 316 DO jj = 2, jpjm1 298 317 DO ji = fs_2, fs_jpim1 ! vector opt. 299 zbtr = 1./ ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 300 ztay(ji,jj,jk) = - ( ztv(ji,jj,jk) - ztv(ji ,jj-1,jk ) ) * zbtr & 301 & + ztay(ji,jj,jk) 302 zsay(ji,jj,jk) = - ( zsv(ji,jj,jk) - zsv(ji ,jj-1,jk ) ) * zbtr & 303 & + zsay(ji,jj,jk) 304 ztaz(ji,jj,jk) = - ( ztw(ji,jj,jk) - ztw(ji ,jj ,jk+1) ) * zbtr & 305 & + ztaz(ji,jj,jk) 306 zsaz(ji,jj,jk) = - ( zsw(ji,jj,jk) - zsw(ji ,jj ,jk+1) ) * zbtr & 307 & + zsaz(ji,jj,jk) 318 !-- Compute zonal divergence by splitting hdivn (see divcur.F90) 319 ! N.B. This computation is not valid along OBCs (if any) 320 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 321 z_hdivn_x = ( e2u(ji ,jj) * fse3u(ji ,jj,jk) * pun(ji ,jj,jk) & 322 & - e2u(ji-1,jj) * fse3u(ji-1,jj,jk) * pun(ji-1,jj,jk) ) * zbtr 323 !-- Compute T/S zonal advection trends 324 ztrdt(ji,jj,jk) = - ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) ) * zbtr + tn(ji,jj,jk) * z_hdivn_x 325 ztrds(ji,jj,jk) = - ( zsu(ji,jj,jk) - zsu(ji-1,jj,jk) ) * zbtr + sn(ji,jj,jk) * z_hdivn_x 308 326 END DO 309 327 END DO 310 328 END DO 311 312 ! horizontal advection: 313 ! make the difference between the new trends ta()/sa() and the 314 ! previous one ztdta()/ztdsa() to have the total advection trends 315 ! to which we substract the vertical trends ztaz()/zsaz() 316 ztdta(:,:,:) = ta(:,:,:) - ztdta(:,:,:) - ztaz(:,:,:) 317 ztdsa(:,:,:) = sa(:,:,:) - ztdsa(:,:,:) - zsaz(:,:,:) 318 319 ! Add the term tn()/sn()*hdivn() to recover the Uh gradh(T/S) trends 320 ztdta(:,:,:) = ztdta(:,:,:) + tn(:,:,:) * hdivn(:,:,:) 321 ztdsa(:,:,:) = ztdsa(:,:,:) + sn(:,:,:) * hdivn(:,:,:) 322 323 CALL trd_mod(ztdta, ztdsa, jpttdlad, 'TRA', kt) 324 325 ! vertical advection: 326 ! Substract the term tn()/sn()*hdivn() to recover the W gradz(T/S) trends 327 ztaz(:,:,:) = ztaz(:,:,:) - tn(:,:,:) * hdivn(:,:,:) 328 zsaz(:,:,:) = zsaz(:,:,:) - sn(:,:,:) * hdivn(:,:,:) 329 330 CALL trd_mod(ztaz, zsaz, jpttdzad, 'TRA', kt) 331 329 CALL trd_mod(ztrdt, ztrds, jptra_trd_xad, 'TRA', kt, cnbpas='bis') ! <<< ADD TO PREVIOUSLY COMPUTED 330 ! 331 ! T/S MERIDIONAL advection trends 332 DO jk = 1, jpkm1 333 DO jj = 2, jpjm1 334 DO ji = fs_2, fs_jpim1 ! vector opt. 335 !-- Compute merid. divergence by splitting hdivn (see divcur.F90) 336 ! N.B. This computation is not valid along OBCs (if any) 337 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 338 z_hdivn_y = ( e1v(ji, jj) * fse3v(ji,jj ,jk) * pvn(ji,jj ,jk) & 339 & - e1v(ji,jj-1) * fse3v(ji,jj-1,jk) * pvn(ji,jj-1,jk) ) * zbtr 340 !-- Compute T/S meridional advection trends 341 ztrdt(ji,jj,jk) = - ( ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) * zbtr + tn(ji,jj,jk) * z_hdivn_y 342 ztrds(ji,jj,jk) = - ( zsv(ji,jj,jk) - zsv(ji,jj-1,jk) ) * zbtr + sn(ji,jj,jk) * z_hdivn_y 343 END DO 344 END DO 345 END DO 346 CALL trd_mod(ztrdt, ztrds, jptra_trd_yad, 'TRA', kt, cnbpas='bis') ! <<< ADD TO PREVIOUSLY COMPUTED 347 ! 348 ! T/S VERTICAL advection trends 349 DO jk = 1, jpkm1 350 DO jj = 2, jpjm1 351 DO ji = fs_2, fs_jpim1 ! vector opt. 352 zbtr1 = 1. / ( e1t(ji,jj) * e2t(ji,jj) ) 353 #if defined key_zco 354 zbtr = zbtr1 355 z_hdivn_x = e2u(ji,jj)*pun(ji,jj,jk) - e2u(ji-1,jj)*pun(ji-1,jj,jk) 356 z_hdivn_y = e1v(ji,jj)*pvn(ji,jj,jk) - e1v(ji,jj-1)*pvn(ji,jj-1,jk) 357 #else 358 zbtr = zbtr1 / fse3t(ji,jj,jk) 359 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) 360 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) 361 #endif 362 z_hdivn = (z_hdivn_x + z_hdivn_y) * zbtr 363 zbtr = zbtr1 / fse3t(ji,jj,jk) 364 ztrdt(ji,jj,jk) = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * zbtr - tn(ji,jj,jk) * z_hdivn 365 ztrds(ji,jj,jk) = - ( zsw(ji,jj,jk) - zsw(ji,jj,jk+1) ) * zbtr - sn(ji,jj,jk) * z_hdivn 366 END DO 367 END DO 368 END DO 369 CALL trd_mod(ztrdt, ztrds, jptra_trd_zad, 'TRA', kt, cnbpas='bis') ! <<< ADD TO PREVIOUSLY COMPUTED 370 ! 332 371 ENDIF 333 372 334 IF(ln_ctl) THEN 335 CALL prt_ctl(tab3d_1=ta, clinfo1=' tvd adv - Ta: ', mask1=tmask, & 336 & tab3d_2=sa, clinfo2=' Sa: ', mask2=tmask, clinfo3='tra') 337 ENDIF 373 IF(ln_ctl) CALL prt_ctl( tab3d_1=ta, clinfo1=' tvd adv - Ta: ', mask1=tmask, & 374 & tab3d_2=sa, clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' ) 338 375 339 376 ! "zonal" mean advective heat and salt transport … … 342 379 pst_adv(:) = ptr_vj( zsv(:,:,:) ) 343 380 ENDIF 344 381 ! 345 382 END SUBROUTINE tra_adv_tvd 346 383 … … 358 395 !! drange (1995) multi-dimensional forward-in-time and upstream- 359 396 !! in-space based differencing for fluid 360 !!361 !! History :362 !! ! 97-04 (L. Mortier) Original code363 !! ! 00-02 (H. Loukos) rewritting for opa8364 !! ! 00-10 (M.A Foujols, E. Kestenare) lateral b.c.365 !! ! 01-03 (E. Kestenare) add key_passivetrc366 !! ! 01-07 (E. Durand G. Madec) adapted for T & S367 !! 8.5 ! 02-06 (G. Madec) F90: Free form and module368 397 !!---------------------------------------------------------------------- 369 !! * Arguments 370 REAL(wp), INTENT( in ) :: & 371 prdt ! ??? 398 REAL(wp), INTENT( in ) :: prdt ! ??? 372 399 REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT( inout ) :: & 373 400 pbef, & ! before field … … 376 403 pbb, & ! monotonic flux in the j direction 377 404 pcc ! monotonic flux in the k direction 378 379 !! * Local declarations 405 !! 380 406 INTEGER :: ji, jj, jk ! dummy loop indices 381 407 INTEGER :: ikm1 … … 496 522 CALL lbc_lnk( pbb, 'V', -1. ) ! changed sign 497 523 CALL lbc_lnk( pcc, 'W', 1. ) ! NO changed sign 498 524 ! 499 525 END SUBROUTINE nonosc 500 526
Note: See TracChangeset
for help on using the changeset viewer.
