Changeset 5901 for branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_ubs.F90
- Timestamp:
- 2015-11-20T09:39:06+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_ubs.F90
r5620 r5901 16 16 USE trc_oce ! share passive tracers/Ocean variables 17 17 USE trd_oce ! trends: ocean variables 18 USE traadv_fct ! acces to routine interp_4th_cpt 18 19 USE trdtra ! trends manager: tracers 19 20 USE dynspg_oce ! choice/control of key cpp for surface pressure gradient … … 38 39 # include "vectopt_loop_substitute.h90" 39 40 !!---------------------------------------------------------------------- 40 !! NEMO/OPA 3. 3 , NEMO Consortium (2010)41 !! NEMO/OPA 3.7 , NEMO Consortium (2015) 41 42 !! $Id$ 42 43 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 44 45 CONTAINS 45 46 46 SUBROUTINE tra_adv_ubs ( kt, kit000, cdtype, p2dt, pun, pvn, pwn,&47 & ptb, ptn, pta, kjpt)47 SUBROUTINE tra_adv_ubs( kt, kit000, cdtype, p2dt, pun, pvn, pwn, & 48 & ptb, ptn, pta, kjpt, kn_ubs_v ) 48 49 !!---------------------------------------------------------------------- 49 50 !! *** ROUTINE tra_adv_ubs *** … … 52 53 !! and add it to the general trend of passive tracer equations. 53 54 !! 54 !! ** Method : The upstream biased scheme (UBS) is based on a 3rd order55 !! ** Method : The 3rd order Upstream Biased Scheme (UBS) is based on an 55 56 !! upstream-biased parabolic interpolation (Shchepetkin and McWilliams 2005) 56 57 !! It is only used in the horizontal direction. … … 61 62 !! where zltu is the second derivative of the before temperature field: 62 63 !! zltu = 1/e3t di[ e2u e3u / e1u di[Tb] ] 63 !! This results in a dissipatively dominant (i.e. hyper-diffusive)64 !! This results in a dissipatively dominant (i.e. hyper-diffusive) 64 65 !! truncation error. The overall performance of the advection scheme 65 66 !! is similar to that reported in (Farrow and Stevens, 1995). 66 !! For stability reasons, the first term of the fluxes which corresponds67 !! For stability reasons, the first term of the fluxes which corresponds 67 68 !! to a second order centered scheme is evaluated using the now velocity 68 69 !! (centered in time) while the second term which is the diffusive part 69 70 !! of the scheme, is evaluated using the before velocity (forward in time). 70 71 !! Note that UBS is not positive. Do not use it on passive tracers. 71 !! On the vertical, the advection is evaluated using a TVD scheme, 72 !! as the UBS have been found to be too diffusive. 72 !! On the vertical, the advection is evaluated using a FCT scheme, 73 !! as the UBS have been found to be too diffusive. 74 !!gm !! kn_ubs_v argument (not coded for the moment) 75 !! controles whether the FCT is based on a 2nd order centrered scheme (kn_ubs_v=2) 76 !! or on a 4th order compact scheme (kn_ubs_v=4). 73 77 !! 74 78 !! ** Action : - update (pta) with the now advective tracer trends … … 81 85 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 82 86 INTEGER , INTENT(in ) :: kjpt ! number of tracers 87 INTEGER , INTENT(in ) :: kn_ubs_v ! number of tracers 83 88 REAL(wp), DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile of tracer time-step 84 89 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun, pvn, pwn ! 3 ocean transport components … … 95 100 IF( nn_timing == 1 ) CALL timing_start('tra_adv_ubs') 96 101 ! 97 CALL wrk_alloc( jpi, jpj, jpk,ztu, ztv, zltu, zltv, zti, ztw )102 CALL wrk_alloc( jpi,jpj,jpk, ztu, ztv, zltu, zltv, zti, ztw ) 98 103 ! 99 104 IF( kt == kit000 ) THEN … … 106 111 IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 107 112 ! 113 zltu(:,:,jpk) = 0._wp ; zltv(:,:,jpk) = 0._wp ! Bottom value : set to zero one for all 114 ztw (:,:,jpk) = 0._wp ; zti (:,:,jpk) = 0._wp 115 IF( lk_vvl ) ztw(:,:, 1 ) = 0._wp ! surface value: set to zero only in vvl case 116 ! 108 117 ! ! =========== 109 118 DO jn = 1, kjpt ! tracer loop 110 119 ! ! =========== 111 ! 1. Bottom value : flux set to zero112 ! ----------------------------------113 zltu(:,:,jpk) = 0.e0 ; zltv(:,:,jpk) = 0.e0114 120 ! 115 DO jk = 1, jpkm1 ! Horizontal slab 116 ! 117 ! Laplacian 118 DO jj = 1, jpjm1 ! First derivative (gradient) 121 DO jk = 1, jpkm1 !== horizontal laplacian of before tracer ==! 122 DO jj = 1, jpjm1 ! First derivative (masked gradient) 119 123 DO ji = 1, fs_jpim1 ! vector opt. 120 zeeu = e2 u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) * umask(ji,jj,jk)121 zeev = e1 v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) * vmask(ji,jj,jk)124 zeeu = e2_e1u(ji,jj) * fse3u(ji,jj,jk) * umask(ji,jj,jk) 125 zeev = e1_e2v(ji,jj) * fse3v(ji,jj,jk) * vmask(ji,jj,jk) 122 126 ztu(ji,jj,jk) = zeeu * ( ptb(ji+1,jj ,jk,jn) - ptb(ji,jj,jk,jn) ) 123 127 ztv(ji,jj,jk) = zeev * ( ptb(ji ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) 124 128 END DO 125 129 END DO 126 DO jj = 2, jpjm1 ! Second derivative (divergence)130 DO jj = 2, jpjm1 ! Second derivative (divergence) 127 131 DO ji = fs_2, fs_jpim1 ! vector opt. 128 zcoef = 1. / ( 6.* fse3t(ji,jj,jk) )132 zcoef = 1._wp / ( 6._wp * fse3t(ji,jj,jk) ) 129 133 zltu(ji,jj,jk) = ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) ) * zcoef 130 134 zltv(ji,jj,jk) = ( ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) * zcoef … … 132 136 END DO 133 137 ! 134 END DO ! End of slab138 END DO 135 139 CALL lbc_lnk( zltu, 'T', 1. ) ; CALL lbc_lnk( zltv, 'T', 1. ) ! Lateral boundary cond. (unchanged sgn) 136 137 140 ! 138 ! Horizontal advective fluxes 139 DO jk = 1, jpkm1 ! Horizontal slab 141 DO jk = 1, jpkm1 !== Horizontal advective fluxes ==! (UBS) 140 142 DO jj = 1, jpjm1 141 143 DO ji = 1, fs_jpim1 ! vector opt. 142 ! upstream transport (x2) 143 zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) ) 144 zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) ) ! upstream transport (x2) 144 145 zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) ) 145 146 zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) ) 146 147 zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) ) 147 ! 2nd order centered advective fluxes (x2)148 ! ! 2nd order centered advective fluxes (x2) 148 149 zcenut = pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj ,jk,jn) ) 149 150 zcenvt = pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji ,jj+1,jk,jn) ) 150 ! UBS advective fluxes151 ! ! UBS advective fluxes 151 152 ztu(ji,jj,jk) = 0.5 * ( zcenut - zfp_ui * zltu(ji,jj,jk) - zfm_ui * zltu(ji+1,jj,jk) ) 152 153 ztv(ji,jj,jk) = 0.5 * ( zcenvt - zfp_vj * zltv(ji,jj,jk) - zfm_vj * zltv(ji,jj+1,jk) ) 153 154 END DO 154 155 END DO 155 END DO ! End of slab156 157 zltu(:,:,:) = pta(:,:,:,jn) ! store pta trends158 159 DO jk = 1, jpkm1 ! Horizontal advective trends156 END DO 157 ! 158 zltu(:,:,:) = pta(:,:,:,jn) ! store the initial trends before its update 159 ! 160 DO jk = 1, jpkm1 !== add the horizontal advective trend ==! 160 161 DO jj = 2, jpjm1 161 162 DO ji = fs_2, fs_jpim1 ! vector opt. … … 166 167 END DO 167 168 ! 168 END DO ! End of slab 169 170 ! Horizontal trend used in tra_adv_ztvd subroutine 171 zltu(:,:,:) = pta(:,:,:,jn) - zltu(:,:,:) 172 169 END DO 170 ! 171 zltu(:,:,:) = pta(:,:,:,jn) - zltu(:,:,:) ! Horizontal advective trend used in vertical 2nd order FCT case 172 ! ! and/or in trend diagnostic (l_trd=T) 173 173 ! 174 174 IF( l_trd ) THEN ! trend diagnostics … … 181 181 IF( jn == jp_sal ) str_adv(:) = ptr_sj( ztv(:,:,:) ) 182 182 ENDIF 183 184 ! TVD scheme for the vertical direction 185 ! ---------------------- 186 IF( l_trd ) zltv(:,:,:) = pta(:,:,:,jn) ! store pta if trend diag. 187 188 ! Bottom value : flux set to zero 189 ztw(:,:,jpk) = 0.e0 ; zti(:,:,jpk) = 0.e0 190 191 ! Surface value 192 IF( lk_vvl ) THEN ; ztw(:,:,1) = 0.e0 ! variable volume : flux set to zero 193 ELSE ; ztw(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn) ! free constant surface 194 ENDIF 195 ! upstream advection with initial mass fluxes & intermediate update 196 ! ------------------------------------------------------------------- 197 ! Interior value 198 DO jk = 2, jpkm1 199 DO jj = 1, jpj 200 DO ji = 1, jpi 201 zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) ) 202 zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) ) 203 ztw(ji,jj,jk) = 0.5 * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) ) 204 END DO 205 END DO 206 END DO 207 ! update and guess with monotonic sheme 208 DO jk = 1, jpkm1 209 z2dtt = p2dt(jk) 210 DO jj = 2, jpjm1 211 DO ji = fs_2, fs_jpim1 ! vector opt. 212 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 213 ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * zbtr 214 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztak 215 zti(ji,jj,jk) = ( ptb(ji,jj,jk,jn) + z2dtt * ( ztak + zltu(ji,jj,jk) ) ) * tmask(ji,jj,jk) 216 END DO 217 END DO 218 END DO 219 ! 220 CALL lbc_lnk( zti, 'T', 1. ) ! Lateral boundary conditions on zti, zsi (unchanged sign) 221 222 ! antidiffusive flux : high order minus low order 223 ztw(:,:,1) = 0.e0 ! Surface value 224 DO jk = 2, jpkm1 ! Interior value 225 DO jj = 1, jpj 226 DO ji = 1, jpi 227 ztw(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) ) - ztw(ji,jj,jk) 228 END DO 229 END DO 230 END DO 231 ! 232 CALL nonosc_z( ptb(:,:,:,jn), ztw, zti, p2dt ) ! monotonicity algorithm 233 234 ! final trend with corrected fluxes 235 DO jk = 1, jpkm1 183 ! 184 ! !== vertical advective trend ==! 185 ! 186 SELECT CASE( kn_ubs_v ) ! select the vertical advection scheme 187 ! 188 CASE( 2 ) ! 2nd order FCT 189 ! 190 IF( l_trd ) zltv(:,:,:) = pta(:,:,:,jn) ! store pta if trend diag. 191 ! 192 ! !* upstream advection with initial mass fluxes & intermediate update ==! 193 DO jk = 2, jpkm1 ! Interior value (w-masked) 194 DO jj = 1, jpj 195 DO ji = 1, jpi 196 zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) ) 197 zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) ) 198 ztw(ji,jj,jk) = 0.5_wp * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) ) * wmask(ji,jj,jk) 199 END DO 200 END DO 201 END DO 202 IF(.NOT.lk_vvl ) THEN ! top ocean value (only in linear free surface as ztw has been w-masked) 203 IF( ln_isfcav ) THEN ! top of the ice-shelf cavities and at the ocean surface 204 DO jj = 1, jpj 205 DO ji = 1, jpi 206 ztw(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn) ! linear free surface 207 END DO 208 END DO 209 ELSE ! no cavities: only at the ocean surface 210 ztw(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn) 211 ENDIF 212 ENDIF 213 ! 214 DO jk = 1, jpkm1 !* trend and after field with monotonic scheme 215 z2dtt = p2dt(jk) 216 DO jj = 2, jpjm1 217 DO ji = fs_2, fs_jpim1 ! vector opt. 218 ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk) ) 219 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztak 220 zti(ji,jj,jk) = ( ptb(ji,jj,jk,jn) + z2dtt * ( ztak + zltu(ji,jj,jk) ) ) * tmask(ji,jj,jk) 221 END DO 222 END DO 223 END DO 224 CALL lbc_lnk( zti, 'T', 1. ) ! Lateral boundary conditions on zti, zsi (unchanged sign) 225 ! 226 ! !* anti-diffusive flux : high order minus low order 227 DO jk = 2, jpkm1 ! Interior value (w-masked) 228 DO jj = 1, jpj 229 DO ji = 1, jpi 230 ztw(ji,jj,jk) = ( 0.5_wp * pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) ) & 231 & - ztw(ji,jj,jk) ) * wmask(ji,jj,jk) 232 END DO 233 END DO 234 END DO 235 ! ! top ocean value: high order == upstream ==>> zwz=0 236 IF(.NOT.lk_vvl ) ztw(:,:, 1 ) = 0._wp ! only ocean surface as interior zwz values have been w-masked 237 ! 238 CALL nonosc_z( ptb(:,:,:,jn), ztw, zti, p2dt ) ! monotonicity algorithm 239 ! 240 CASE( 4 ) ! 4th order COMPACT 241 CALL interp_4th_cpt( ptn(:,:,:,jn) , ztw ) ! 4th order compact interpolation of T at w-point 242 DO jk = 2, jpkm1 243 DO jj = 2, jpjm1 244 DO ji = fs_2, fs_jpim1 245 ztw(ji,jj,jk) = pwn(ji,jj,jk) * ztw(ji,jj,jk) * wmask(ji,jj,jk) 246 END DO 247 END DO 248 END DO 249 IF(.NOT.lk_vvl ) ztw(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1,jn) !!gm ISF & 4th COMPACT doesn't work 250 ! 251 END SELECT 252 ! 253 DO jk = 1, jpkm1 ! final trend with corrected fluxes 236 254 DO jj = 2, jpjm1 237 255 DO ji = fs_2, fs_jpim1 ! vector opt. 238 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 239 ! k- vertical advective trends 240 ztra = - zbtr * ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) 241 ! added to the general tracer trends 242 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 243 END DO 244 END DO 245 END DO 246 247 ! Save the final vertical advective trends 248 IF( l_trd ) THEN ! vertical advective trend diagnostics 256 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk) ) 257 END DO 258 END DO 259 END DO 260 ! 261 IF( l_trd ) THEN ! vertical advective trend diagnostics 249 262 DO jk = 1, jpkm1 ! (compute -w.dk[ptn]= -dk[w.ptn] + ptn.dk[w]) 250 263 DO jj = 2, jpjm1 251 264 DO ji = fs_2, fs_jpim1 ! vector opt. 252 z btr = 1.e0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )253 z_hdivn = ( pwn(ji,jj,jk) - pwn(ji,jj,jk+1) ) * zbtr254 zltv(ji,jj,jk) = pta(ji,jj,jk,jn) - zltv(ji,jj,jk) + ptn(ji,jj,jk,jn) * z_hdivn265 zltv(ji,jj,jk) = pta(ji,jj,jk,jn) - zltv(ji,jj,jk) & 266 & + ptn(ji,jj,jk,jn) * ( pwn(ji,jj,jk) - pwn(ji,jj,jk+1) ) & 267 & / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk) ) 255 268 END DO 256 269 END DO … … 261 274 END DO 262 275 ! 263 CALL wrk_dealloc( jpi, jpj, jpk,ztu, ztv, zltu, zltv, zti, ztw )276 CALL wrk_dealloc( jpi,jpj,jpk, ztu, ztv, zltu, zltv, zti, ztw ) 264 277 ! 265 278 IF( nn_timing == 1 ) CALL timing_stop('tra_adv_ubs') … … 294 307 IF( nn_timing == 1 ) CALL timing_start('nonosc_z') 295 308 ! 296 CALL wrk_alloc( jpi, jpj, jpk,zbetup, zbetdo )309 CALL wrk_alloc( jpi,jpj,jpk, zbetup, zbetdo ) 297 310 ! 298 311 zbig = 1.e+40_wp 299 312 zrtrn = 1.e-15_wp 300 313 zbetup(:,:,:) = 0._wp ; zbetdo(:,:,:) = 0._wp 301 314 ! 302 315 ! Search local extrema 303 316 ! -------------------- 304 ! large negative value (-zbig) inside land317 ! ! large negative value (-zbig) inside land 305 318 pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) ) 306 319 paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) ) 307 ! search maximum in neighbourhood308 DO jk = 1, jpkm1 320 ! 321 DO jk = 1, jpkm1 ! search maximum in neighbourhood 309 322 ikm1 = MAX(jk-1,1) 310 323 DO jj = 2, jpjm1 … … 316 329 END DO 317 330 END DO 318 ! large positive value (+zbig) inside land331 ! ! large positive value (+zbig) inside land 319 332 pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) ) 320 333 paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) ) 321 ! search minimum in neighbourhood322 DO jk = 1, jpkm1 334 ! 335 DO jk = 1, jpkm1 ! search minimum in neighbourhood 323 336 ikm1 = MAX(jk-1,1) 324 337 DO jj = 2, jpjm1 … … 330 343 END DO 331 344 END DO 332 333 ! restore masked values to zero 345 ! ! restore masked values to zero 334 346 pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) 335 347 paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) 336 337 338 ! 2. Positive and negative part of fluxes and beta terms 339 ! ------------------------------------------------------ 340 348 ! 349 ! Positive and negative part of fluxes and beta terms 350 ! --------------------------------------------------- 341 351 DO jk = 1, jpkm1 342 352 z2dtt = p2dt(jk) … … 347 357 zneg = MAX( 0., pcc(ji ,jj ,jk ) ) - MIN( 0., pcc(ji ,jj ,jk+1) ) 348 358 ! up & down beta terms 349 zbt = e1 t(ji,jj) *e2t(ji,jj) * fse3t(ji,jj,jk) / z2dtt359 zbt = e1e2t(ji,jj) * fse3t(ji,jj,jk) / z2dtt 350 360 zbetup(ji,jj,jk) = ( zbetup(ji,jj,jk) - paft(ji,jj,jk) ) / (zpos+zrtrn) * zbt 351 361 zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zbetdo(ji,jj,jk) ) / (zneg+zrtrn) * zbt … … 353 363 END DO 354 364 END DO 365 ! 355 366 ! monotonic flux in the k direction, i.e. pcc 356 367 ! ------------------------------------------- … … 366 377 END DO 367 378 ! 368 CALL wrk_dealloc( jpi, jpj, jpk,zbetup, zbetdo )379 CALL wrk_dealloc( jpi,jpj,jpk, zbetup, zbetdo ) 369 380 ! 370 381 IF( nn_timing == 1 ) CALL timing_stop('nonosc_z')
Note: See TracChangeset
for help on using the changeset viewer.
