- Timestamp:
- 2018-07-13T09:28:50+02:00 (6 years ago)
- Location:
- NEMO/branches/2018/dev_r9838_ENHANCE04_RK3
- Files:
-
- 1 edited
- 1 copied
-
. (copied) (copied from NEMO/trunk)
-
src/OCE/DYN/dynnxt.F90 (modified) (13 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DYN/dynnxt.F90
r9598 r9939 64 64 CONTAINS 65 65 66 SUBROUTINE dyn_nxt 66 SUBROUTINE dyn_nxt( kt ) 67 67 !!---------------------------------------------------------------------- 68 68 !! *** ROUTINE dyn_nxt *** … … 83 83 !! * Apply the time filter applied and swap of the dynamics 84 84 !! arrays to start the next time step: 85 !! (ub,vb) = (un,vn) + atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ]85 !! (ub,vb) = (un,vn) + rn_atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ] 86 86 !! (un,vn) = (ua,va). 87 87 !! Note that with flux form advection and non linear free surface, … … 92 92 !! un,vn now horizontal velocity of next time-step 93 93 !!---------------------------------------------------------------------- 94 INTEGER, INTENT( in ) :: kt 94 INTEGER, INTENT( in ) :: kt ! ocean time-step index 95 95 ! 96 96 INTEGER :: ji, jj, jk ! dummy loop indices 97 97 INTEGER :: ikt ! local integers 98 REAL(wp) :: zue3a, zue3n, zue3b, zuf, zcoef 99 REAL(wp) :: zve3a, zve3n, zve3b, zvf , z1_2dt! - -98 REAL(wp) :: zue3a, zue3n, zue3b, zuf, zcoef ! local scalars 99 REAL(wp) :: zve3a, zve3n, zve3b, zvf ! - - 100 100 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zue, zve 101 101 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ze3u_f, ze3v_f, zua, zva … … 132 132 ! so that asselin contribution is removed at the same time 133 133 DO jk = 1, jpkm1 134 un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:)*r1_hu_n(:,:) + un_b(:,:) ) *umask(:,:,jk)135 vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:)*r1_hv_n(:,:) + vn_b(:,:) ) *vmask(:,:,jk)134 un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:)*r1_hu_n(:,:) + un_b(:,:) ) * umask(:,:,jk) 135 vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:)*r1_hv_n(:,:) + vn_b(:,:) ) * vmask(:,:,jk) 136 136 END DO 137 137 ENDIF … … 152 152 !!$ Do we need a call to bdy_vol here?? 153 153 ! 154 IF( l_trddyn ) THEN ! prepare the atf trend computation + some diagnostics155 z1_2dt = 1._wp / (2. * rdt) ! Euler or leap-frog time step156 IF( neuler == 0 .AND. kt == nit000 ) z1_2dt = 1._wp / rdt157 !158 ! ! Kinetic energy and Conversion159 IF( ln_KE_trd ) CALL trd_dyn( ua, va, jpdyn_ken, kt )160 !161 IF( ln_dyn_trd ) THEN ! 3D output: total momentum trends162 zua(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * z1_2dt163 zva(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * z1_2dt164 CALL iom_put( "utrd_tot", zua ) ! total momentum trends, except the asselin timefilter154 IF( l_trddyn ) THEN !* prepare the atf trend computation + some diagnostics 155 IF( ln_KE_trd ) CALL trd_dyn( ua, va, jpdyn_ken, kt ) ! Kinetic energy and Conversion 156 IF( ln_dyn_trd ) THEN ! 3D output: total momentum trends 157 IF( ln_dynadv_vec ) THEN 158 zua(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * r1_Dt 159 zva(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * r1_Dt 160 ELSE 161 zua(:,:,:) = ( e3u_a(:,:,:)*ua(:,:,:) - e3u_b(:,:,:)*ub(:,:,:) ) / e3u_n(:,:,:) * r1_Dt 162 zva(:,:,:) = ( e3v_a(:,:,:)*va(:,:,:) - e3v_b(:,:,:)*vb(:,:,:) ) / e3v_n(:,:,:) * r1_Dt 163 ENDIF 164 CALL iom_put( "utrd_tot", zua ) ! total momentum trends, except the asselin filter 165 165 CALL iom_put( "vtrd_tot", zva ) 166 166 ENDIF 167 ! 168 zua(:,:,:) = un(:,:,:) ! save the now velocity before the asselin filter 169 zva(:,:,:) = vn(:,:,:) ! (caution: there will be a shift by 1 timestep in the 170 ! ! computation of the asselin filter trends) 167 zua(:,:,:) = un(:,:,:) ! save the now velocity before the asselin filter 168 zva(:,:,:) = vn(:,:,:) ! (caution: the Asselin filter trends computation will be shifted by 1 timestep) 171 169 ENDIF 172 170 173 171 ! Time filter and swap of dynamics arrays 174 ! --------------------------------------- ---175 IF( neuler == 0 .AND. kt == nit000 ) THEN !* Euler at first time-step: only swap172 ! --------------------------------------- 173 IF( l_1st_euler ) THEN !== Euler at 1st time-step ==! (swap only) 176 174 DO jk = 1, jpkm1 177 175 un(:,:,jk) = ua(:,:,jk) ! un <-- ua 178 176 vn(:,:,jk) = va(:,:,jk) 179 177 END DO 180 IF( .NOT.ln_linssh ) THEN ! e3._b <-- e3._n 181 !!gm BUG ???? I don't understand why it is not : e3._n <-- e3._a 178 IF( .NOT.ln_linssh ) THEN ! e3._n <-- e3._a 182 179 DO jk = 1, jpkm1 183 ! e3t_b(:,:,jk) = e3t_n(:,:,jk)184 ! e3u_b(:,:,jk) = e3u_n(:,:,jk)185 ! e3v_b(:,:,jk) = e3v_n(:,:,jk)186 !187 180 e3t_n(:,:,jk) = e3t_a(:,:,jk) 188 181 e3u_n(:,:,jk) = e3u_a(:,:,jk) 189 182 e3v_n(:,:,jk) = e3v_a(:,:,jk) 190 183 END DO 191 !!gm BUG end 192 ENDIF 193 ! 194 195 ELSE !* Leap-Frog : Asselin filter and swap 184 ENDIF 185 ! 186 ELSE !== Leap-Frog ==! (Asselin filter and swap) 187 ! 196 188 ! ! =============! 197 189 IF( ln_linssh ) THEN ! Fixed volume ! … … 200 192 DO jj = 1, jpj 201 193 DO ji = 1, jpi 202 zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )203 zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )194 zuf = un(ji,jj,jk) + rn_atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) ) 195 zvf = vn(ji,jj,jk) + rn_atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) ) 204 196 ! 205 197 ub(ji,jj,jk) = zuf ! ub <-- filtered velocity … … 213 205 ELSE ! Variable volume ! 214 206 ! ! ================! 215 ! Before scale factor at t-points 216 ! (used as a now filtered scale factor until the swap) 217 ! ---------------------------------------------------- 207 ! Before scale factor at t-points (used as a now filtered scale factor until the swap) 218 208 DO jk = 1, jpkm1 219 e3t_b(:,:,jk) = e3t_n(:,:,jk) + atfp * ( e3t_b(:,:,jk) - 2._wp * e3t_n(:,:,jk) + e3t_a(:,:,jk) )209 e3t_b(:,:,jk) = e3t_n(:,:,jk) + rn_atfp * ( e3t_b(:,:,jk) - 2._wp * e3t_n(:,:,jk) + e3t_a(:,:,jk) ) 220 210 END DO 221 211 ! Add volume filter correction: compatibility with tracer advection scheme 222 ! => time filter + conservation correction (only at the first level)223 zcoef = atfp * rdt * r1_rau0212 ! => time filter + conservation correction (only at the first level) 213 zcoef = rn_atfp * rn_Dt * r1_rho0 224 214 225 215 e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1) … … 232 222 IF( jk <= nk_rnf(ji,jj) ) THEN 233 223 e3t_b(ji,jj,jk) = e3t_b(ji,jj,jk) - zcoef * ( - rnf_b(ji,jj) + rnf(ji,jj) ) & 234 &* ( e3t_n(ji,jj,jk) / h_rnf(ji,jj) ) * tmask(ji,jj,jk)224 & * ( e3t_n(ji,jj,jk) / h_rnf(ji,jj) ) * tmask(ji,jj,jk) 235 225 ENDIF 236 END DO237 END DO238 END DO226 END DO 227 END DO 228 END DO 239 229 ELSE 240 230 e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef * ( -rnf_b(:,:) + rnf(:,:))*tmask(:,:,1) 241 231 ENDIF 242 END 232 ENDIF 243 233 244 234 IF ( ln_isf ) THEN ! if ice shelf melting … … 253 243 END DO 254 244 END DO 255 END 245 ENDIF 256 246 ! 257 247 IF( ln_dynadv_vec ) THEN ! Asselin filter applied on velocity … … 262 252 DO jj = 1, jpj 263 253 DO ji = 1, jpi 264 zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )265 zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )254 zuf = un(ji,jj,jk) + rn_atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) ) 255 zvf = vn(ji,jj,jk) + rn_atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) ) 266 256 ! 267 257 ub(ji,jj,jk) = zuf ! ub <-- filtered velocity … … 289 279 zve3b = e3v_b(ji,jj,jk) * vb(ji,jj,jk) 290 280 ! 291 zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n + zue3a ) ) / ze3u_f(ji,jj,jk)292 zvf = ( zve3n + atfp * ( zve3b - 2._wp * zve3n + zve3a ) ) / ze3v_f(ji,jj,jk)281 zuf = ( zue3n + rn_atfp * ( zue3b - 2._wp * zue3n + zue3a ) ) / ze3u_f(ji,jj,jk) 282 zvf = ( zve3n + rn_atfp * ( zve3b - 2._wp * zve3n + zve3a ) ) / ze3v_f(ji,jj,jk) 293 283 ! 294 284 ub(ji,jj,jk) = zuf ! ub <-- filtered velocity … … 322 312 ENDIF 323 313 ! 324 ENDIF ! neuler =/0314 ENDIF ! end Leap-Frog time stepping 325 315 ! 326 316 ! Set "now" and "before" barotropic velocities for next time step: 327 ! JC: Would be more clever to swap variables than to make a full vertical 328 ! integration 317 ! JC: Would be more clever to swap variables than to make a full vertical integration 329 318 ! 330 319 ! … … 360 349 ENDIF 361 350 IF( l_trddyn ) THEN ! 3D output: asselin filter trends on momentum 362 zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * z1_2dt363 zva(:,:,:) = ( vb(:,:,:) - zva(:,:,:) ) * z1_2dt351 zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * r1_Dt 352 zva(:,:,:) = ( vb(:,:,:) - zva(:,:,:) ) * r1_Dt 364 353 CALL trd_dyn( zua, zva, jpdyn_atf, kt ) 365 354 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
