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