- Timestamp:
- 2016-07-19T10:38:35+02:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/NERC/dev_r5549_BDY_ZEROGRAD/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_ubs.F90
r5147 r6808 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 ! … … 35 35 36 36 !! * Substitutions 37 # include "domzgr_substitute.h90"38 37 # include "vectopt_loop_substitute.h90" 39 38 !!---------------------------------------------------------------------- 40 !! NEMO/OPA 3. 3 , NEMO Consortium (2010)39 !! NEMO/OPA 3.7 , NEMO Consortium (2015) 41 40 !! $Id$ 42 41 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 44 43 CONTAINS 45 44 46 SUBROUTINE tra_adv_ubs ( kt, kit000, cdtype, p2dt, pun, pvn, pwn,&47 & ptb, ptn, pta, kjpt)45 SUBROUTINE tra_adv_ubs( kt, kit000, cdtype, p2dt, pun, pvn, pwn, & 46 & ptb, ptn, pta, kjpt, kn_ubs_v ) 48 47 !!---------------------------------------------------------------------- 49 48 !! *** ROUTINE tra_adv_ubs *** … … 52 51 !! and add it to the general trend of passive tracer equations. 53 52 !! 54 !! ** Method : The upstream biased scheme (UBS) is based on a 3rd order53 !! ** Method : The 3rd order Upstream Biased Scheme (UBS) is based on an 55 54 !! upstream-biased parabolic interpolation (Shchepetkin and McWilliams 2005) 56 55 !! It is only used in the horizontal direction. … … 61 60 !! where zltu is the second derivative of the before temperature field: 62 61 !! zltu = 1/e3t di[ e2u e3u / e1u di[Tb] ] 63 !! This results in a dissipatively dominant (i.e. hyper-diffusive)62 !! This results in a dissipatively dominant (i.e. hyper-diffusive) 64 63 !! truncation error. The overall performance of the advection scheme 65 64 !! is similar to that reported in (Farrow and Stevens, 1995). 66 !! For stability reasons, the first term of the fluxes which corresponds65 !! For stability reasons, the first term of the fluxes which corresponds 67 66 !! to a second order centered scheme is evaluated using the now velocity 68 67 !! (centered in time) while the second term which is the diffusive part 69 68 !! of the scheme, is evaluated using the before velocity (forward in time). 70 69 !! 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. 70 !! On the vertical, the advection is evaluated using a FCT scheme, 71 !! as the UBS have been found to be too diffusive. 72 !! kn_ubs_v argument controles whether the FCT is based on 73 !! a 2nd order centrered scheme (kn_ubs_v=2) or on a 4th order compact 74 !! scheme (kn_ubs_v=4). 73 75 !! 74 !! ** Action : - update (pta) with the now advective tracer trends 76 !! ** Action : - update pta with the now advective tracer trends 77 !! - send trends to trdtra module for further diagnostcs (l_trdtra=T) 78 !! - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T) 75 79 !! 76 80 !! Reference : Shchepetkin, A. F., J. C. McWilliams, 2005, Ocean Modelling, 9, 347-404. … … 81 85 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 82 86 INTEGER , INTENT(in ) :: kjpt ! number of tracers 83 REAL(wp), DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile of tracer time-step 87 INTEGER , INTENT(in ) :: kn_ubs_v ! number of tracers 88 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 84 89 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun, pvn, pwn ! 3 ocean transport components 85 90 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before and now tracer fields … … 87 92 ! 88 93 INTEGER :: ji, jj, jk, jn ! dummy loop indices 89 REAL(wp) :: ztra, zbtr, zcoef , z2dtt! local scalars94 REAL(wp) :: ztra, zbtr, zcoef ! local scalars 90 95 REAL(wp) :: zfp_ui, zfm_ui, zcenut, ztak, zfp_wk, zfm_wk ! - - 91 96 REAL(wp) :: zfp_vj, zfm_vj, zcenvt, zeeu, zeev, z_hdivn ! - - … … 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 ztw (:,:, 1 ) = 0._wp ! surface & bottom value : set to zero for all tracers 114 zltu(:,:,jpk) = 0._wp ; zltv(:,:,jpk) = 0._wp 115 ztw (:,:,jpk) = 0._wp ; zti (:,:,jpk) = 0._wp 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) * e3u_n(ji,jj,jk) * umask(ji,jj,jk) 125 zeev = e1_e2v(ji,jj) * e3v_n(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 * e3t_n(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. 162 163 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) & 163 164 & - ( ztu(ji,jj,jk) - ztu(ji-1,jj ,jk) & 164 & + ztv(ji,jj,jk) - ztv(ji ,jj-1,jk) ) / ( e1e2t(ji,jj) * fse3t(ji,jj,jk))165 & + ztv(ji,jj,jk) - ztv(ji ,jj-1,jk) ) * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 165 166 END DO 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( ln_linssh ) 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 DO jj = 2, jpjm1 216 DO ji = fs_2, fs_jpim1 ! vector opt. 217 ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) / e3t_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) + p2dt * ( 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( ln_linssh ) 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( ln_linssh ) 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) ) * r1_e1e2t(ji,jj) / e3t_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 & * r1_e1e2t(ji,jj) / e3t_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') … … 281 293 !! in-space based differencing for fluid 282 294 !!---------------------------------------------------------------------- 283 REAL(wp), INTENT(in ) , DIMENSION(jpk) :: p2dt ! vertical profile oftracer time-step295 REAL(wp), INTENT(in ) :: p2dt ! tracer time-step 284 296 REAL(wp), DIMENSION (jpi,jpj,jpk) :: pbef ! before field 285 297 REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) :: paft ! after field … … 288 300 INTEGER :: ji, jj, jk ! dummy loop indices 289 301 INTEGER :: ikm1 ! local integer 290 REAL(wp) :: zpos, zneg, zbt, za, zb, zc, zbig, zrtrn , z2dtt! local scalars302 REAL(wp) :: zpos, zneg, zbt, za, zb, zc, zbig, zrtrn ! local scalars 291 303 REAL(wp), POINTER, DIMENSION(:,:,:) :: zbetup, zbetdo 292 304 !!---------------------------------------------------------------------- … … 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 z2dtt = p2dt(jk)343 351 DO jj = 2, jpjm1 344 352 DO ji = fs_2, fs_jpim1 ! vector opt. … … 347 355 zneg = MAX( 0., pcc(ji ,jj ,jk ) ) - MIN( 0., pcc(ji ,jj ,jk+1) ) 348 356 ! up & down beta terms 349 zbt = e1 t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) / z2dtt357 zbt = e1e2t(ji,jj) * e3t_n(ji,jj,jk) / p2dt 350 358 zbetup(ji,jj,jk) = ( zbetup(ji,jj,jk) - paft(ji,jj,jk) ) / (zpos+zrtrn) * zbt 351 359 zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zbetdo(ji,jj,jk) ) / (zneg+zrtrn) * zbt … … 353 361 END DO 354 362 END DO 363 ! 355 364 ! monotonic flux in the k direction, i.e. pcc 356 365 ! ------------------------------------------- … … 366 375 END DO 367 376 ! 368 CALL wrk_dealloc( jpi, jpj, jpk,zbetup, zbetdo )377 CALL wrk_dealloc( jpi,jpj,jpk, zbetup, zbetdo ) 369 378 ! 370 379 IF( nn_timing == 1 ) CALL timing_stop('nonosc_z')
Note: See TracChangeset
for help on using the changeset viewer.