Changeset 4736 for branches/2014/dev_r4743_NOC2_ZTS
- Timestamp:
- 2014-08-05T12:24:00+02:00 (10 years ago)
- Location:
- branches/2014/dev_r4743_NOC2_ZTS/NEMOGCM
- Files:
-
- 5 edited
-
CONFIG/SHARED/namelist_ref (modified) (2 diffs)
-
NEMO/OPA_SRC/DYN/dynadv.F90 (modified) (6 diffs)
-
NEMO/OPA_SRC/DYN/dynzad.F90 (modified) (4 diffs)
-
NEMO/OPA_SRC/TRA/traadv.F90 (modified) (7 diffs)
-
NEMO/OPA_SRC/TRA/traadv_tvd.F90 (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4743_NOC2_ZTS/NEMOGCM/CONFIG/SHARED/namelist_ref
r4699 r4736 684 684 ln_traadv_qck = .false. ! QUICKEST scheme 685 685 ln_traadv_msc_ups= .false. ! use upstream scheme within muscl 686 ln_traadv_tvd_zts= .false. ! TVD scheme with sub-timestepping of vertical tracer advection 686 687 / 687 688 !----------------------------------------------------------------------- … … 758 759 ln_dynadv_cen2= .false. ! flux form - 2nd order centered scheme 759 760 ln_dynadv_ubs = .false. ! flux form - 3rd order UBS scheme 761 ln_dynzad_zts = .false. ! Use (T) sub timestepping for vertical momentum advection 760 762 / 761 763 !----------------------------------------------------------------------- -
branches/2014/dev_r4743_NOC2_ZTS/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv.F90
r4624 r4736 30 30 LOGICAL, PUBLIC :: ln_dynadv_cen2 !: flux form - 2nd order centered scheme flag 31 31 LOGICAL, PUBLIC :: ln_dynadv_ubs !: flux form - 3rd order UBS scheme flag 32 LOGICAL, PUBLIC :: ln_dynzad_zts !: vertical advection with sub-timestepping (requires vector form) 32 33 33 34 INTEGER :: nadv ! choice of the formulation and scheme for the advection … … 64 65 CALL dyn_keg ( kt ) ! vector form : horizontal gradient of kinetic energy 65 66 CALL dyn_zad ( kt ) ! vector form : vertical advection 66 CASE ( 1 ) 67 CASE ( 1 ) 68 CALL dyn_keg ( kt ) ! vector form : horizontal gradient of kinetic energy 69 CALL dyn_zad_zts ( kt ) ! vector form : vertical advection with sub-timestepping 70 CASE ( 2 ) 67 71 CALL dyn_adv_cen2( kt ) ! 2nd order centered scheme 68 CASE ( 2)72 CASE ( 3 ) 69 73 CALL dyn_adv_ubs ( kt ) ! 3rd order UBS scheme 70 74 ! … … 91 95 INTEGER :: ios ! Local integer output status for namelist read 92 96 !! 93 NAMELIST/namdyn_adv/ ln_dynadv_vec, ln_dynadv_cen2 , ln_dynadv_ubs 97 NAMELIST/namdyn_adv/ ln_dynadv_vec, ln_dynadv_cen2 , ln_dynadv_ubs, ln_dynzad_zts 94 98 !!---------------------------------------------------------------------- 95 99 … … 111 115 WRITE(numout,*) ' 2nd order centred advection scheme ln_dynadv_cen2 = ', ln_dynadv_cen2 112 116 WRITE(numout,*) ' 3rd order UBS advection scheme ln_dynadv_ubs = ', ln_dynadv_ubs 117 WRITE(numout,*) ' Sub timestepping of vertical advection ln_dynzad_zts = ', ln_dynzad_zts 113 118 ENDIF 114 119 … … 120 125 121 126 IF( ioptio /= 1 ) CALL ctl_stop( 'Choose ONE advection scheme in namelist namdyn_adv' ) 127 IF( ln_dynzad_zts .AND. .NOT. ln_dynadv_vec ) & 128 CALL ctl_stop( 'Sub timestepping of vertical advection requires vector form; set ln_dynadv_vec = .TRUE.' ) 122 129 123 130 ! ! Set nadv 124 131 IF( ln_dynadv_vec ) nadv = 0 125 IF( ln_dynadv_cen2 ) nadv = 1 126 IF( ln_dynadv_ubs ) nadv = 2 132 IF( ln_dynzad_zts ) nadv = 1 133 IF( ln_dynadv_cen2 ) nadv = 2 134 IF( ln_dynadv_ubs ) nadv = 3 127 135 IF( lk_esopa ) nadv = -1 128 136 … … 130 138 WRITE(numout,*) 131 139 IF( nadv == 0 ) WRITE(numout,*) ' vector form : keg + zad + vor is used' 132 IF( nadv == 1 ) WRITE(numout,*) ' flux form : 2nd order scheme is used' 133 IF( nadv == 2 ) WRITE(numout,*) ' flux form : UBS scheme is used' 140 IF( nadv == 1 ) WRITE(numout,*) ' vector form : keg + zad_zts + vor is used' 141 IF( nadv == 2 ) WRITE(numout,*) ' flux form : 2nd order scheme is used' 142 IF( nadv == 3 ) WRITE(numout,*) ' flux form : UBS scheme is used' 134 143 IF( nadv == -1 ) WRITE(numout,*) ' esopa test: use all advection formulation' 135 144 ENDIF -
branches/2014/dev_r4743_NOC2_ZTS/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90
r3294 r4736 27 27 PRIVATE 28 28 29 PUBLIC dyn_zad ! routine called by step.F90 29 PUBLIC dyn_zad ! routine called by dynadv.F90 30 PUBLIC dyn_zad_zts ! routine called by dynadv.F90 30 31 31 32 !! * Substitutions … … 83 84 DO jj = 2, jpj ! vertical fluxes 84 85 DO ji = fs_2, jpi ! vector opt. 85 zww(ji,jj) = 0.25 * e1t(ji,jj) * e2t(ji,jj) * wn(ji,jj,jk)86 zww(ji,jj) = 0.25_wp * e1t(ji,jj) * e2t(ji,jj) * wn(ji,jj,jk) 86 87 END DO 87 88 END DO … … 95 96 DO jj = 2, jpjm1 ! Surface and bottom values set to zero 96 97 DO ji = fs_2, fs_jpim1 ! vector opt. 97 zwuw(ji,jj, 1 ) = 0. e098 zwvw(ji,jj, 1 ) = 0. e099 zwuw(ji,jj,jpk) = 0. e0100 zwvw(ji,jj,jpk) = 0. e098 zwuw(ji,jj, 1 ) = 0._wp 99 zwvw(ji,jj, 1 ) = 0._wp 100 zwuw(ji,jj,jpk) = 0._wp 101 zwvw(ji,jj,jpk) = 0._wp 101 102 END DO 102 103 END DO … … 132 133 END SUBROUTINE dyn_zad 133 134 135 SUBROUTINE dyn_zad_zts ( kt ) 136 !!---------------------------------------------------------------------- 137 !! *** ROUTINE dynzad_zts *** 138 !! 139 !! ** Purpose : Compute the now vertical momentum advection trend and 140 !! add it to the general trend of momentum equation. This version 141 !! uses sub-timesteps for improved numerical stability with small 142 !! vertical grid sizes. This is especially relevant when using 143 !! embedded ice with thin surface boxes. 144 !! 145 !! ** Method : The now vertical advection of momentum is given by: 146 !! w dz(u) = ua + 1/(e1u*e2u*e3u) mk+1[ mi(e1t*e2t*wn) dk(un) ] 147 !! w dz(v) = va + 1/(e1v*e2v*e3v) mk+1[ mj(e1t*e2t*wn) dk(vn) ] 148 !! Add this trend to the general trend (ua,va): 149 !! (ua,va) = (ua,va) + w dz(u,v) 150 !! 151 !! ** Action : - Update (ua,va) with the vert. momentum adv. trends 152 !! - Save the trends in (ztrdu,ztrdv) ('key_trddyn') 153 !!---------------------------------------------------------------------- 154 INTEGER, INTENT(in) :: kt ! ocean time-step inedx 155 ! 156 INTEGER :: ji, jj, jk, jl ! dummy loop indices 157 INTEGER :: jnzts = 5 ! number of sub-timesteps for vertical advection 158 INTEGER :: jtb, jtn, jta ! sub timestep pointers for leap-frog/euler forward steps 159 REAL(wp) :: zua, zva ! temporary scalars 160 REAL(wp) :: zr_rdt ! temporary scalar 161 REAL(wp) :: z2dtzts ! length of Euler forward sub-timestep for vertical advection 162 REAL(wp) :: zts ! length of sub-timestep for vertical advection 163 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwuw , zwvw, zww 164 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdu, ztrdv 165 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zus , zvs 166 !!---------------------------------------------------------------------- 167 ! 168 IF( nn_timing == 1 ) CALL timing_start('dyn_zad_zts') 169 ! 170 CALL wrk_alloc( jpi,jpj,jpk, zwuw , zwvw, zww ) 171 CALL wrk_alloc( jpi,jpj,jpk,3, zus, zvs ) 172 ! 173 IF( kt == nit000 ) THEN 174 IF(lwp)WRITE(numout,*) 175 IF(lwp)WRITE(numout,*) 'dyn_zad_zts : arakawa advection scheme with sub-timesteps' 176 ENDIF 177 178 IF( l_trddyn ) THEN ! Save ua and va trends 179 CALL wrk_alloc( jpi, jpj, jpk, ztrdu, ztrdv ) 180 ztrdu(:,:,:) = ua(:,:,:) 181 ztrdv(:,:,:) = va(:,:,:) 182 ENDIF 183 184 IF( neuler == 0 .AND. kt == nit000 ) THEN 185 z2dtzts = rdt / REAL( jnzts, wp ) ! = rdt (restart with Euler time stepping) 186 ELSE 187 z2dtzts = 2._wp * rdt / REAL( jnzts, wp ) ! = 2 rdt (leapfrog) 188 ENDIF 189 190 DO jk = 2, jpkm1 ! Calculate and store vertical fluxes 191 DO jj = 2, jpj 192 DO ji = fs_2, jpi ! vector opt. 193 zww(ji,jj,jk) = 0.25_wp * e1t(ji,jj) * e2t(ji,jj) * wn(ji,jj,jk) 194 END DO 195 END DO 196 END DO 197 198 DO jj = 2, jpjm1 ! Surface and bottom advective fluxes set to zero 199 DO ji = fs_2, fs_jpim1 ! vector opt. 200 zwuw(ji,jj, 1 ) = 0._wp 201 zwvw(ji,jj, 1 ) = 0._wp 202 zwuw(ji,jj,jpk) = 0._wp 203 zwvw(ji,jj,jpk) = 0._wp 204 END DO 205 END DO 206 207 ! Start with before values and use sub timestepping to reach after values 208 209 zus(:,:,:,1) = ub(:,:,:) 210 zvs(:,:,:,1) = vb(:,:,:) 211 212 DO jl = 1, jnzts ! Start of sub timestepping loop 213 214 IF( jl == 1 ) THEN ! Euler forward to kick things off 215 jtb = 1 ; jtn = 1 ; jta = 2 216 zts = z2dtzts 217 ELSEIF( jl == 2 ) THEN ! First leapfrog step 218 jtb = 1 ; jtn = 2 ; jta = 3 219 zts = 2._wp * z2dtzts 220 ELSE ! Shuffle pointers for subsequent leapfrog steps 221 jtb = MOD(jtb,3) + 1 222 jtn = MOD(jtn,3) + 1 223 jta = MOD(jta,3) + 1 224 ENDIF 225 226 DO jk = 2, jpkm1 ! Vertical momentum advection at level w and u- and v- vertical 227 DO jj = 2, jpjm1 ! vertical momentum advection at w-point 228 DO ji = fs_2, fs_jpim1 ! vector opt. 229 zwuw(ji,jj,jk) = ( zww(ji+1,jj ,jk) + zww(ji,jj,jk) ) * ( zus(ji,jj,jk-1,jtn)-zus(ji,jj,jk,jtn) ) 230 zwvw(ji,jj,jk) = ( zww(ji ,jj+1,jk) + zww(ji,jj,jk) ) * ( zvs(ji,jj,jk-1,jtn)-zvs(ji,jj,jk,jtn) ) 231 END DO 232 END DO 233 END DO 234 DO jk = 1, jpkm1 ! Vertical momentum advection at u- and v-points 235 DO jj = 2, jpjm1 236 DO ji = fs_2, fs_jpim1 ! vector opt. 237 ! ! vertical momentum advective trends 238 zua = - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) ) 239 zva = - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) ) 240 zus(ji,jj,jk,jta) = zus(ji,jj,jk,jtb) + zua * zts 241 zvs(ji,jj,jk,jta) = zvs(ji,jj,jk,jtb) + zva * zts 242 END DO 243 END DO 244 END DO 245 246 END DO ! End of sub timestepping loop 247 248 zr_rdt = 1._wp / ( REAL( jnzts, wp ) * z2dtzts ) 249 DO jk = 1, jpkm1 ! Recover trends over the outer timestep 250 DO jj = 2, jpjm1 251 DO ji = fs_2, fs_jpim1 ! vector opt. 252 ! ! vertical momentum advective trends 253 ! ! add the trends to the general momentum trends 254 ua(ji,jj,jk) = ua(ji,jj,jk) + ( zus(ji,jj,jk,jta) - ub(ji,jj,jk)) * zr_rdt 255 va(ji,jj,jk) = va(ji,jj,jk) + ( zvs(ji,jj,jk,jta) - vb(ji,jj,jk)) * zr_rdt 256 END DO 257 END DO 258 END DO 259 260 IF( l_trddyn ) THEN ! save the vertical advection trends for diagnostic 261 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 262 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 263 CALL trd_mod(ztrdu, ztrdv, jpdyn_trd_zad, 'DYN', kt) 264 CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv ) 265 ENDIF 266 ! ! Control print 267 IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' zad - Ua: ', mask1=umask, & 268 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 269 ! 270 CALL wrk_dealloc( jpi,jpj,jpk, zwuw , zwvw, zww ) 271 CALL wrk_dealloc( jpi,jpj,jpk,3, zus, zvs ) 272 ! 273 IF( nn_timing == 1 ) CALL timing_stop('dyn_zad_zts') 274 ! 275 END SUBROUTINE dyn_zad_zts 276 134 277 !!====================================================================== 135 278 END MODULE dynzad -
branches/2014/dev_r4743_NOC2_ZTS/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90
r4624 r4736 43 43 LOGICAL :: ln_traadv_cen2 ! 2nd order centered scheme flag 44 44 LOGICAL :: ln_traadv_tvd ! TVD scheme flag 45 LOGICAL :: ln_traadv_tvd_zts ! TVD scheme flag with vertical sub time-stepping 45 46 LOGICAL :: ln_traadv_muscl ! MUSCL scheme flag 46 47 LOGICAL :: ln_traadv_muscl2 ! MUSCL2 scheme flag … … 120 121 CASE ( 5 ) ; CALL tra_adv_ubs ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts ) ! UBS 121 122 CASE ( 6 ) ; CALL tra_adv_qck ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts ) ! QUICKEST 123 CASE ( 7 ) ; CALL tra_adv_tvd_zts( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts ) ! TVD ZTS 122 124 ! 123 125 CASE (-1 ) !== esopa: test all possibility with control print ==! … … 166 168 & ln_traadv_muscl, ln_traadv_muscl2, & 167 169 & ln_traadv_ubs , ln_traadv_qck, & 168 & ln_traadv_msc_ups 170 & ln_traadv_msc_ups, ln_traadv_tvd_zts 169 171 !!---------------------------------------------------------------------- 170 172 … … 190 192 WRITE(numout,*) ' QUICKEST advection scheme ln_traadv_qck = ', ln_traadv_qck 191 193 WRITE(numout,*) ' upstream scheme within muscl ln_traadv_msc_ups = ', ln_traadv_msc_ups 194 WRITE(numout,*) ' TVD advection scheme with zts ln_traadv_tvd_zts = ', ln_traadv_tvd_zts 192 195 ENDIF 193 196 … … 199 202 IF( ln_traadv_ubs ) ioptio = ioptio + 1 200 203 IF( ln_traadv_qck ) ioptio = ioptio + 1 204 IF( ln_traadv_tvd_zts) ioptio = ioptio + 1 201 205 IF( lk_esopa ) ioptio = 1 202 206 … … 210 214 IF( ln_traadv_ubs ) nadv = 5 211 215 IF( ln_traadv_qck ) nadv = 6 216 IF( ln_traadv_tvd_zts) nadv = 7 212 217 IF( lk_esopa ) nadv = -1 213 218 … … 220 225 IF( nadv == 5 ) WRITE(numout,*) ' UBS scheme is used' 221 226 IF( nadv == 6 ) WRITE(numout,*) ' QUICKEST scheme is used' 227 IF( nadv == 7 ) WRITE(numout,*) ' TVD ZTS scheme is used' 222 228 IF( nadv == -1 ) WRITE(numout,*) ' esopa test: use all advection scheme' 223 229 ENDIF -
branches/2014/dev_r4743_NOC2_ZTS/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90
r4499 r4736 37 37 PRIVATE 38 38 39 PUBLIC tra_adv_tvd ! routine called by step.F90 39 PUBLIC tra_adv_tvd ! routine called by traadv.F90 40 PUBLIC tra_adv_tvd_zts ! routine called by traadv.F90 40 41 41 42 LOGICAL :: l_trd ! flag to compute trends … … 247 248 END SUBROUTINE tra_adv_tvd 248 249 250 SUBROUTINE tra_adv_tvd_zts ( kt, kit000, cdtype, p2dt, pun, pvn, pwn, & 251 & ptb, ptn, pta, kjpt ) 252 !!---------------------------------------------------------------------- 253 !! *** ROUTINE tra_adv_tvd_zts *** 254 !! 255 !! ** Purpose : Compute the now trend due to total advection of 256 !! tracers and add it to the general trend of tracer equations 257 !! 258 !! ** Method : TVD ZTS scheme, i.e. 2nd order centered scheme with 259 !! corrected flux (monotonic correction). This version use sub- 260 !! timestepping for the vertical advection which increases stability 261 !! when vertical metrics are small. 262 !! note: - this advection scheme needs a leap-frog time scheme 263 !! 264 !! ** Action : - update (pta) with the now advective tracer trends 265 !! - save the trends 266 !!---------------------------------------------------------------------- 267 USE oce , ONLY: zwx => ua , zwy => va ! (ua,va) used as workspace 268 ! 269 INTEGER , INTENT(in ) :: kt ! ocean time-step index 270 INTEGER , INTENT(in ) :: kit000 ! first time step index 271 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 272 INTEGER , INTENT(in ) :: kjpt ! number of tracers 273 REAL(wp), DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile of tracer time-step 274 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun, pvn, pwn ! 3 ocean velocity components 275 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before and now tracer fields 276 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 277 ! 278 REAL(wp), DIMENSION( jpk ) :: zts ! length of sub-timestep for vertical advection 279 REAL(wp), DIMENSION( jpk ) :: zr_p2dt ! reciprocal of tracer timestep 280 INTEGER :: ji, jj, jk, jl, jn ! dummy loop indices 281 INTEGER :: jnzts = 5 ! number of sub-timesteps for vertical advection 282 INTEGER :: jtb, jtn, jta ! sub timestep pointers for leap-frog/euler forward steps 283 INTEGER :: jtaken ! toggle for collecting appropriate fluxes from sub timesteps 284 REAL(wp) :: z_rzts ! Fractional length of Euler forward sub-timestep for vertical advection 285 REAL(wp) :: z2dtt, zbtr, ztra ! local scalar 286 REAL(wp) :: zfp_ui, zfp_vj, zfp_wk ! - - 287 REAL(wp) :: zfm_ui, zfm_vj, zfm_wk ! - - 288 REAL(wp), POINTER, DIMENSION(:,: ) :: zwx_sav , zwy_sav 289 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwi, zwz, zhdiv, zwz_sav, zwzts 290 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdx, ztrdy, ztrdz 291 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: ztrs 292 !!---------------------------------------------------------------------- 293 ! 294 IF( nn_timing == 1 ) CALL timing_start('tra_adv_tvd_zts') 295 ! 296 CALL wrk_alloc( jpi, jpj, zwx_sav, zwy_sav ) 297 CALL wrk_alloc( jpi, jpj, jpk, zwi, zwz , zhdiv, zwz_sav, zwzts ) 298 CALL wrk_alloc( jpi, jpj, jpk, 3, ztrs ) 299 ! 300 IF( kt == kit000 ) THEN 301 IF(lwp) WRITE(numout,*) 302 IF(lwp) WRITE(numout,*) 'tra_adv_tvd_zts : TVD ZTS advection scheme on ', cdtype 303 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 304 ENDIF 305 ! 306 l_trd = .FALSE. 307 IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 308 ! 309 IF( l_trd ) THEN 310 CALL wrk_alloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz ) 311 ztrdx(:,:,:) = 0._wp ; ztrdy(:,:,:) = 0._wp ; ztrdz(:,:,:) = 0._wp 312 ENDIF 313 ! 314 zwi(:,:,:) = 0._wp 315 z_rzts = 1._wp / REAL( jnzts, wp ) 316 zr_p2dt(:) = 1._wp / p2dt(:) 317 ! 318 ! ! =========== 319 DO jn = 1, kjpt ! tracer loop 320 ! ! =========== 321 ! 1. Bottom value : flux set to zero 322 ! ---------------------------------- 323 zwx(:,:,jpk) = 0._wp ; zwz(:,:,jpk) = 0._wp 324 zwy(:,:,jpk) = 0._wp ; zwi(:,:,jpk) = 0._wp 325 326 ! 2. upstream advection with initial mass fluxes & intermediate update 327 ! -------------------------------------------------------------------- 328 ! upstream tracer flux in the i and j direction 329 DO jk = 1, jpkm1 330 DO jj = 1, jpjm1 331 DO ji = 1, fs_jpim1 ! vector opt. 332 ! upstream scheme 333 zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) ) 334 zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) ) 335 zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) ) 336 zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) ) 337 zwx(ji,jj,jk) = 0.5_wp * ( zfp_ui * ptb(ji,jj,jk,jn) + zfm_ui * ptb(ji+1,jj ,jk,jn) ) 338 zwy(ji,jj,jk) = 0.5_wp * ( zfp_vj * ptb(ji,jj,jk,jn) + zfm_vj * ptb(ji ,jj+1,jk,jn) ) 339 END DO 340 END DO 341 END DO 342 343 ! upstream tracer flux in the k direction 344 ! Surface value 345 IF( lk_vvl ) THEN ; zwz(:,:, 1 ) = 0._wp ! volume variable 346 ELSE ; zwz(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1,jn) ! linear free surface 347 ENDIF 348 ! Interior value 349 DO jk = 2, jpkm1 350 DO jj = 1, jpj 351 DO ji = 1, jpi 352 zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) ) 353 zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) ) 354 zwz(ji,jj,jk) = 0.5_wp * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) ) 355 END DO 356 END DO 357 END DO 358 359 ! total advective trend 360 DO jk = 1, jpkm1 361 z2dtt = p2dt(jk) 362 DO jj = 2, jpjm1 363 DO ji = fs_2, fs_jpim1 ! vector opt. 364 zbtr = 1._wp / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 365 ! total intermediate advective trends 366 ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 367 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & 368 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) 369 ! update and guess with monotonic sheme 370 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 371 zwi(ji,jj,jk) = ( ptb(ji,jj,jk,jn) + z2dtt * ztra ) * tmask(ji,jj,jk) 372 END DO 373 END DO 374 END DO 375 ! ! Lateral boundary conditions on zwi (unchanged sign) 376 CALL lbc_lnk( zwi, 'T', 1. ) 377 378 ! ! trend diagnostics (contribution of upstream fluxes) 379 IF( l_trd ) THEN 380 ! store intermediate advective trends 381 ztrdx(:,:,:) = zwx(:,:,:) ; ztrdy(:,:,:) = zwy(:,:,:) ; ztrdz(:,:,:) = zwz(:,:,:) 382 END IF 383 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) 384 IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN 385 IF( jn == jp_tem ) htr_adv(:) = ptr_vj( zwy(:,:,:) ) 386 IF( jn == jp_sal ) str_adv(:) = ptr_vj( zwy(:,:,:) ) 387 ENDIF 388 389 ! 3. antidiffusive flux : high order minus low order 390 ! -------------------------------------------------- 391 ! antidiffusive flux on i and j 392 393 394 DO jk = 1, jpkm1 395 396 DO jj = 1, jpjm1 397 DO ji = 1, fs_jpim1 ! vector opt. 398 zwx_sav(ji,jj) = zwx(ji,jj,jk) 399 zwy_sav(ji,jj) = zwy(ji,jj,jk) 400 401 zwx(ji,jj,jk) = 0.5_wp * pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj,jk,jn) ) 402 zwy(ji,jj,jk) = 0.5_wp * pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj+1,jk,jn) ) 403 END DO 404 END DO 405 406 DO jj = 2, jpjm1 ! partial horizontal divergence 407 DO ji = fs_2, fs_jpim1 408 zhdiv(ji,jj,jk) = ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk) & 409 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk) ) 410 END DO 411 END DO 412 413 DO jj = 1, jpjm1 414 DO ji = 1, fs_jpim1 ! vector opt. 415 zwx(ji,jj,jk) = zwx(ji,jj,jk) - zwx_sav(ji,jj) 416 zwy(ji,jj,jk) = zwy(ji,jj,jk) - zwy_sav(ji,jj) 417 END DO 418 END DO 419 END DO 420 421 ! antidiffusive flux on k 422 zwz(:,:,1) = 0._wp ! Surface value 423 zwz_sav(:,:,:) = zwz(:,:,:) 424 ! 425 ztrs(:,:,:,1) = ptb(:,:,:,jn) 426 zwzts(:,:,:) = 0._wp 427 428 DO jl = 1, jnzts ! Start of sub timestepping loop 429 430 IF( jl == 1 ) THEN ! Euler forward to kick things off 431 jtb = 1 ; jtn = 1 ; jta = 2 432 zts(:) = p2dt(:) * z_rzts 433 jtaken = MOD( jnzts + 1 , 2) ! Toggle to collect every second flux 434 ! starting at jl =1 if jnzts is odd; 435 ! starting at jl =2 otherwise 436 ELSEIF( jl == 2 ) THEN ! First leapfrog step 437 jtb = 1 ; jtn = 2 ; jta = 3 438 zts(:) = 2._wp * p2dt(:) * z_rzts 439 ELSE ! Shuffle pointers for subsequent leapfrog steps 440 jtb = MOD(jtb,3) + 1 441 jtn = MOD(jtn,3) + 1 442 jta = MOD(jta,3) + 1 443 ENDIF 444 DO jk = 2, jpkm1 ! Interior value 445 DO jj = 2, jpjm1 446 DO ji = fs_2, fs_jpim1 447 zwz(ji,jj,jk) = 0.5_wp * pwn(ji,jj,jk) * ( ztrs(ji,jj,jk,jtn) + ztrs(ji,jj,jk-1,jtn) ) 448 IF( jtaken == 0 ) zwzts(ji,jj,jk) = zwzts(ji,jj,jk) + zwz(ji,jj,jk)*zts(jk) ! Accumulate time-weighted vertcal flux 449 END DO 450 END DO 451 END DO 452 453 jtaken = MOD( jtaken + 1 , 2 ) 454 455 DO jk = 2, jpkm1 ! Interior value 456 DO jj = 2, jpjm1 457 DO ji = fs_2, fs_jpim1 458 zbtr = 1._wp / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 459 ! total advective trends 460 ztra = - zbtr * ( zhdiv(ji,jj,jk) + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) 461 ztrs(ji,jj,jk,jta) = ztrs(ji,jj,jk,jtb) + zts(jk) * ztra 462 END DO 463 END DO 464 END DO 465 466 END DO 467 468 DO jk = 2, jpkm1 ! Anti-diffusive vertical flux using average flux from the sub-timestepping 469 DO jj = 2, jpjm1 470 DO ji = fs_2, fs_jpim1 471 zwz(ji,jj,jk) = zwzts(ji,jj,jk) * zr_p2dt(jk) - zwz_sav(ji,jj,jk) 472 END DO 473 END DO 474 END DO 475 CALL lbc_lnk( zwx, 'U', -1. ) ; CALL lbc_lnk( zwy, 'V', -1. ) ! Lateral bondary conditions 476 CALL lbc_lnk( zwz, 'W', 1. ) 477 478 ! 4. monotonicity algorithm 479 ! ------------------------- 480 CALL nonosc( ptb(:,:,:,jn), zwx, zwy, zwz, zwi, p2dt ) 481 482 483 ! 5. final trend with corrected fluxes 484 ! ------------------------------------ 485 DO jk = 1, jpkm1 486 DO jj = 2, jpjm1 487 DO ji = fs_2, fs_jpim1 ! vector opt. 488 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 489 ! total advective trends 490 ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 491 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & 492 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) 493 ! add them to the general tracer trends 494 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 495 END DO 496 END DO 497 END DO 498 499 ! ! trend diagnostics (contribution of upstream fluxes) 500 IF( l_trd ) THEN 501 ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:) ! <<< Add to previously computed 502 ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:) ! <<< Add to previously computed 503 ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:) ! <<< Add to previously computed 504 505 CALL trd_tra( kt, cdtype, jn, jptra_trd_xad, ztrdx, pun, ptn(:,:,:,jn) ) 506 CALL trd_tra( kt, cdtype, jn, jptra_trd_yad, ztrdy, pvn, ptn(:,:,:,jn) ) 507 CALL trd_tra( kt, cdtype, jn, jptra_trd_zad, ztrdz, pwn, ptn(:,:,:,jn) ) 508 END IF 509 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) 510 IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN 511 IF( jn == jp_tem ) htr_adv(:) = ptr_vj( zwy(:,:,:) ) + htr_adv(:) 512 IF( jn == jp_sal ) str_adv(:) = ptr_vj( zwy(:,:,:) ) + str_adv(:) 513 ENDIF 514 ! 515 END DO 516 ! 517 CALL wrk_dealloc( jpi, jpj, jpk, zwi, zwz, zhdiv, zwz_sav, zwzts ) 518 CALL wrk_dealloc( jpi, jpj, jpk, 3, ztrs ) 519 CALL wrk_dealloc( jpi, jpj, zwx_sav, zwy_sav ) 520 IF( l_trd ) CALL wrk_dealloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz ) 521 ! 522 IF( nn_timing == 1 ) CALL timing_stop('tra_adv_tvd_zts') 523 ! 524 END SUBROUTINE tra_adv_tvd_zts 249 525 250 526 SUBROUTINE nonosc( pbef, paa, pbb, pcc, paft, p2dt )
Note: See TracChangeset
for help on using the changeset viewer.
