Changeset 6140 for trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90
- Timestamp:
- 2015-12-21T12:35:23+01:00 (8 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90 (modified) (14 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90
r5930 r6140 18 18 !! 3.3 ! 2011-03 (P. Oddo) Bug fix for time-splitting+(BDY-OBC) and not VVL 19 19 !! 3.5 ! 2013-07 (J. Chanut) Compliant with time splitting changes 20 !! 3. 7! 2014-04 (G. Madec) add the diagnostic of the time filter trends20 !! 3.6 ! 2014-04 (G. Madec) add the diagnostic of the time filter trends 21 21 !! 3.7 ! 2015-11 (J. Chanut) Free surface simplification 22 22 !!------------------------------------------------------------------------- 23 23 24 24 !!------------------------------------------------------------------------- 25 !! dyn_nxt : obtain the next (after) horizontal velocity25 !! dyn_nxt : obtain the next (after) horizontal velocity 26 26 !!------------------------------------------------------------------------- 27 USE oce ! ocean dynamics and tracers 28 USE dom_oce ! ocean space and time domain 29 USE sbc_oce ! Surface boundary condition: ocean fields 30 USE phycst ! physical constants 31 USE dynadv ! dynamics: vector invariant versus flux form 32 USE domvvl ! variable volume 33 USE bdy_oce ! ocean open boundary conditions 34 USE bdydta ! ocean open boundary conditions 35 USE bdydyn ! ocean open boundary conditions 36 USE bdyvol ! ocean open boundary condition (bdy_vol routines) 37 USE trd_oce ! trends: ocean variables 38 USE trddyn ! trend manager: dynamics 39 USE trdken ! trend manager: kinetic energy 27 USE oce ! ocean dynamics and tracers 28 USE dom_oce ! ocean space and time domain 29 USE sbc_oce ! Surface boundary condition: ocean fields 30 USE phycst ! physical constants 31 USE dynadv ! dynamics: vector invariant versus flux form 32 USE dynspg_ts ! surface pressure gradient: split-explicit scheme 33 USE domvvl ! variable volume 34 USE bdy_oce ! ocean open boundary conditions 35 USE bdydta ! ocean open boundary conditions 36 USE bdydyn ! ocean open boundary conditions 37 USE bdyvol ! ocean open boundary condition (bdy_vol routines) 38 USE trd_oce ! trends: ocean variables 39 USE trddyn ! trend manager: dynamics 40 USE trdken ! trend manager: kinetic energy 40 41 ! 41 USE in_out_manager 42 USE iom 43 USE lbclnk 44 USE lib_mpp 45 USE wrk_nemo 46 USE prtctl 47 USE timing 42 USE in_out_manager ! I/O manager 43 USE iom ! I/O manager library 44 USE lbclnk ! lateral boundary condition (or mpp link) 45 USE lib_mpp ! MPP library 46 USE wrk_nemo ! Memory Allocation 47 USE prtctl ! Print control 48 USE timing ! Timing 48 49 #if defined key_agrif 49 50 USE agrif_opa_interp … … 55 56 PUBLIC dyn_nxt ! routine called by step.F90 56 57 57 !! * Substitutions58 # include "domzgr_substitute.h90"59 58 !!---------------------------------------------------------------------- 60 59 !! NEMO/OPA 3.3 , NEMO Consortium (2010) … … 85 84 !! (ub,vb) = (un,vn) + atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ] 86 85 !! (un,vn) = (ua,va). 87 !! Note that with flux form advection and variable volume layer88 !! (lk_vvl=T), the time filter is applied on thickness weighted89 !! velocity.86 !! Note that with flux form advection and non linear free surface, 87 !! the time filter is applied on thickness weighted velocity. 88 !! As a result, dyn_nxt MUST be called after tra_nxt. 90 89 !! 91 90 !! ** Action : ub,vb filtered before horizontal velocity of next time-step … … 95 94 ! 96 95 INTEGER :: ji, jj, jk ! dummy loop indices 97 INTEGER :: ik u, ikv! local integers98 REAL(wp) :: zue3a, zue3n, zue3b, zuf 96 INTEGER :: ikt ! local integers 97 REAL(wp) :: zue3a, zue3n, zue3b, zuf, zcoef ! local scalars 99 98 REAL(wp) :: zve3a, zve3n, zve3b, zvf, z1_2dt ! - - 100 99 REAL(wp), POINTER, DIMENSION(:,:) :: zue, zve … … 104 103 IF( nn_timing == 1 ) CALL timing_start('dyn_nxt') 105 104 ! 106 IF( ln_dynspg_ts ) CALL wrk_alloc( jpi,jpj, zue, zve ) 107 IF( lk_vvl.AND.(.NOT.ln_dynadv_vec ) ) CALL wrk_alloc( jpi,jpj,jpk, ze3u_f, ze3v_f ) 108 IF( l_trddyn ) CALL wrk_alloc( jpi,jpj,jpk, zua, zva ) 105 IF( ln_dynspg_ts ) CALL wrk_alloc( jpi,jpj, zue, zve) 106 IF( l_trddyn ) CALL wrk_alloc( jpi,jpj,jpk, zua, zva) 109 107 ! 110 108 IF( kt == nit000 ) THEN … … 117 115 ! Ensure below that barotropic velocities match time splitting estimate 118 116 ! Compute actual transport and replace it with ts estimate at "after" time step 119 zue(:,:) = fse3u_a(:,:,1) * ua(:,:,1) * umask(:,:,1)120 zve(:,:) = fse3v_a(:,:,1) * va(:,:,1) * vmask(:,:,1)117 zue(:,:) = e3u_a(:,:,1) * ua(:,:,1) * umask(:,:,1) 118 zve(:,:) = e3v_a(:,:,1) * va(:,:,1) * vmask(:,:,1) 121 119 DO jk = 2, jpkm1 122 zue(:,:) = zue(:,:) + fse3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)123 zve(:,:) = zve(:,:) + fse3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)120 zue(:,:) = zue(:,:) + e3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk) 121 zve(:,:) = zve(:,:) + e3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk) 124 122 END DO 125 123 DO jk = 1, jpkm1 126 ua(:,:,jk) = ( ua(:,:,jk) - zue(:,:) * hur_a(:,:) + ua_b(:,:) ) * umask(:,:,jk)127 va(:,:,jk) = ( va(:,:,jk) - zve(:,:) * hvr_a(:,:) + va_b(:,:) ) * vmask(:,:,jk)124 ua(:,:,jk) = ( ua(:,:,jk) - zue(:,:) * r1_hu_a(:,:) + ua_b(:,:) ) * umask(:,:,jk) 125 va(:,:,jk) = ( va(:,:,jk) - zve(:,:) * r1_hv_a(:,:) + va_b(:,:) ) * vmask(:,:,jk) 128 126 END DO 129 130 IF 127 ! 128 IF( .NOT.ln_bt_fw ) THEN 131 129 ! Remove advective velocity from "now velocities" 132 130 ! prior to asselin filtering … … 151 149 # if defined key_bdy 152 150 ! !* BDY open boundaries 153 IF( lk_bdy .AND. ln_dynspg_exp ) CALL bdy_dyn( kt )154 IF( lk_bdy .AND. ln_dynspg_ts ) CALL bdy_dyn( kt, dyn3d_only=.true. )151 IF( lk_bdy .AND. ln_dynspg_exp ) CALL bdy_dyn( kt ) 152 IF( lk_bdy .AND. ln_dynspg_ts ) CALL bdy_dyn( kt, dyn3d_only=.true. ) 155 153 156 154 !!$ Do we need a call to bdy_vol here?? … … 184 182 vn(:,:,jk) = va(:,:,jk) 185 183 END DO 186 IF (lk_vvl) THEN184 IF(.NOT.ln_linssh ) THEN 187 185 DO jk = 1, jpkm1 188 fse3t_b(:,:,jk) = fse3t_n(:,:,jk)189 fse3u_b(:,:,jk) = fse3u_n(:,:,jk)190 fse3v_b(:,:,jk) = fse3v_n(:,:,jk)191 END DO186 e3t_b(:,:,jk) = e3t_n(:,:,jk) 187 e3u_b(:,:,jk) = e3u_n(:,:,jk) 188 e3v_b(:,:,jk) = e3v_n(:,:,jk) 189 END DO 192 190 ENDIF 193 191 ELSE !* Leap-Frog : Asselin filter and swap 194 192 ! ! =============! 195 IF( .NOT. lk_vvl ) THEN! Fixed volume !193 IF( ln_linssh ) THEN ! Fixed volume ! 196 194 ! ! =============! 197 195 DO jk = 1, jpkm1 … … 214 212 ! (used as a now filtered scale factor until the swap) 215 213 ! ---------------------------------------------------- 216 IF (ln_dynspg_ts.AND.ln_bt_fw) THEN 217 ! No asselin filtering on thicknesses if forward time splitting 218 fse3t_b(:,:,:) = fse3t_n(:,:,:) 214 IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN ! No asselin filtering on thicknesses if forward time splitting 215 e3t_b(:,:,1:jpkm1) = e3t_n(:,:,1:jpkm1) 219 216 ELSE 220 fse3t_b(:,:,:) = fse3t_n(:,:,:) + atfp * ( fse3t_b(:,:,:) - 2._wp * fse3t_n(:,:,:) + fse3t_a(:,:,:) ) 217 DO jk = 1, jpkm1 218 e3t_b(:,:,jk) = e3t_n(:,:,jk) + atfp * ( e3t_b(:,:,jk) - 2._wp * e3t_n(:,:,jk) + e3t_a(:,:,jk) ) 219 END DO 221 220 ! Add volume filter correction: compatibility with tracer advection scheme 222 221 ! => time filter + conservation correction (only at the first level) 223 IF ( nn_isf == 0) THEN ! if no ice shelf melting 224 fse3t_b(:,:,1) = fse3t_b(:,:,1) - atfp * rdt * r1_rau0 * ( emp_b(:,:) - emp(:,:) & 225 & -rnf_b(:,:) + rnf(:,:) ) * tmask(:,:,1) 222 zcoef = atfp * rdt * r1_rau0 223 IF ( .NOT. ln_isf ) THEN ! if no ice shelf melting 224 e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef * ( emp_b(:,:) - emp(:,:) & 225 & - rnf_b(:,:) + rnf(:,:) ) * tmask(:,:,1) 226 226 ELSE ! if ice shelf melting 227 DO jj = 1,jpj 228 DO ji = 1,jpi 229 jk = mikt(ji,jj) 230 fse3t_b(ji,jj,jk) = fse3t_b(ji,jj,jk) - atfp * rdt * r1_rau0 & 231 & * ( (emp_b(ji,jj) - emp(ji,jj) ) & 232 & - (rnf_b(ji,jj) - rnf(ji,jj) ) & 233 & + (fwfisf_b(ji,jj) - fwfisf(ji,jj)) ) * tmask(ji,jj,jk) 227 DO jj = 1, jpj 228 DO ji = 1, jpi 229 ikt = mikt(ji,jj) 230 e3t_b(ji,jj,ikt) = e3t_b(ji,jj,ikt) - zcoef * ( emp_b (ji,jj) - emp (ji,jj) & 231 & - rnf_b (ji,jj) + rnf (ji,jj) & 232 & + fwfisf_b(ji,jj) - fwfisf(ji,jj) ) * tmask(ji,jj,ikt) 234 233 END DO 235 234 END DO … … 237 236 ENDIF 238 237 ! 239 IF( ln_dynadv_vec ) THEN 240 ! Before scale factor at (u/v)-points 241 ! ----------------------------------- 242 CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3u_b(:,:,:), 'U' ) 243 CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3v_b(:,:,:), 'V' ) 244 ! Leap-Frog - Asselin filter and swap: applied on velocity 245 ! ----------------------------------- 238 IF( ln_dynadv_vec ) THEN ! Asselin filter applied on velocity 239 ! Before filtered scale factor at (u/v)-points 240 CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' ) 241 CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' ) 246 242 DO jk = 1, jpkm1 247 243 DO jj = 1, jpj … … 258 254 END DO 259 255 ! 260 ELSE 261 ! Temporary filtered scale factor at (u/v)-points (will become before scale factor) 262 !------------------------------------------------ 263 CALL dom_vvl_interpol( fse3t_b(:,:,:), ze3u_f, 'U' ) 264 CALL dom_vvl_interpol( fse3t_b(:,:,:), ze3v_f, 'V' ) 265 ! Leap-Frog - Asselin filter and swap: applied on thickness weighted velocity 266 ! ----------------------------------- =========================== 256 ELSE ! Asselin filter applied on thickness weighted velocity 257 ! 258 CALL wrk_alloc( jpi,jpj,jpk, ze3u_f, ze3v_f ) 259 ! Before filtered scale factor at (u/v)-points stored in ze3u_f, ze3v_f 260 CALL dom_vvl_interpol( e3t_b(:,:,:), ze3u_f, 'U' ) 261 CALL dom_vvl_interpol( e3t_b(:,:,:), ze3v_f, 'V' ) 267 262 DO jk = 1, jpkm1 268 263 DO jj = 1, jpj 269 264 DO ji = 1, jpi 270 zue3a = ua(ji,jj,jk) * fse3u_a(ji,jj,jk)271 zve3a = va(ji,jj,jk) * fse3v_a(ji,jj,jk)272 zue3n = un(ji,jj,jk) * fse3u_n(ji,jj,jk)273 zve3n = vn(ji,jj,jk) * fse3v_n(ji,jj,jk)274 zue3b = ub(ji,jj,jk) * fse3u_b(ji,jj,jk)275 zve3b = vb(ji,jj,jk) * fse3v_b(ji,jj,jk)265 zue3a = e3u_a(ji,jj,jk) * ua(ji,jj,jk) 266 zve3a = e3v_a(ji,jj,jk) * va(ji,jj,jk) 267 zue3n = e3u_n(ji,jj,jk) * un(ji,jj,jk) 268 zve3n = e3v_n(ji,jj,jk) * vn(ji,jj,jk) 269 zue3b = e3u_b(ji,jj,jk) * ub(ji,jj,jk) 270 zve3b = e3v_b(ji,jj,jk) * vb(ji,jj,jk) 276 271 ! 277 272 zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n + zue3a ) ) / ze3u_f(ji,jj,jk) … … 285 280 END DO 286 281 END DO 287 fse3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1) ! e3u_b <-- filtered scale factor 288 fse3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1) 282 e3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1) ! e3u_b <-- filtered scale factor 283 e3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1) 284 ! 285 CALL wrk_dealloc( jpi,jpj,jpk, ze3u_f, ze3v_f ) 289 286 ENDIF 290 287 ! 291 288 ENDIF 292 289 ! 293 IF (ln_dynspg_ts.AND.ln_bt_fw) THEN290 IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN 294 291 ! Revert "before" velocities to time split estimate 295 292 ! Doing it here also means that asselin filter contribution is removed 296 zue(:,:) = fse3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1)297 zve(:,:) = fse3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)293 zue(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1) 294 zve(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1) 298 295 DO jk = 2, jpkm1 299 zue(:,:) = zue(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)300 zve(:,:) = zve(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)296 zue(:,:) = zue(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk) 297 zve(:,:) = zve(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk) 301 298 END DO 302 299 DO jk = 1, jpkm1 303 ub(:,:,jk) = ub(:,:,jk) - (zue(:,:) * hur(:,:) - un_b(:,:)) * umask(:,:,jk)304 vb(:,:,jk) = vb(:,:,jk) - (zve(:,:) * hvr(:,:) - vn_b(:,:)) * vmask(:,:,jk)300 ub(:,:,jk) = ub(:,:,jk) - (zue(:,:) * r1_hu_n(:,:) - un_b(:,:)) * umask(:,:,jk) 301 vb(:,:,jk) = vb(:,:,jk) - (zve(:,:) * r1_hv_n(:,:) - vn_b(:,:)) * vmask(:,:,jk) 305 302 END DO 306 303 ENDIF … … 313 310 ! 314 311 ! 315 IF (lk_vvl) THEN316 hu_b(:,:) = 0.317 hv_b(:,:) = 0.318 DO jk = 1, jpkm1319 hu_b(:,:) = hu_b(:,:) + fse3u_b(:,:,jk) * umask(:,:,jk)320 hv_b(:,:) = hv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk)312 IF(.NOT.ln_linssh ) THEN 313 hu_b(:,:) = e3u_b(:,:,1) * umask(:,:,1) 314 hv_b(:,:) = e3v_b(:,:,1) * vmask(:,:,1) 315 DO jk = 2, jpkm1 316 hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk) 317 hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk) 321 318 END DO 322 hur_b(:,:) = umask_i(:,:) / ( hu_b(:,:) + 1._wp - umask_i(:,:) ) 323 hvr_b(:,:) = vmask_i(:,:) / ( hv_b(:,:) + 1._wp - vmask_i(:,:) ) 324 ENDIF 325 ! 326 un_b(:,:) = 0._wp ; vn_b(:,:) = 0._wp 327 ub_b(:,:) = 0._wp ; vb_b(:,:) = 0._wp 328 ! 329 DO jk = 1, jpkm1 330 DO jj = 1, jpj 331 DO ji = 1, jpi 332 un_b(ji,jj) = un_b(ji,jj) + fse3u_a(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk) 333 vn_b(ji,jj) = vn_b(ji,jj) + fse3v_a(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk) 334 ! 335 ub_b(ji,jj) = ub_b(ji,jj) + fse3u_b(ji,jj,jk) * ub(ji,jj,jk) * umask(ji,jj,jk) 336 vb_b(ji,jj) = vb_b(ji,jj) + fse3v_b(ji,jj,jk) * vb(ji,jj,jk) * vmask(ji,jj,jk) 337 END DO 338 END DO 319 r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) ) 320 r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) ) 321 ENDIF 322 ! 323 un_b(:,:) = e3u_a(:,:,1) * un(:,:,1) * umask(:,:,1) 324 ub_b(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1) 325 vn_b(:,:) = e3v_a(:,:,1) * vn(:,:,1) * vmask(:,:,1) 326 vb_b(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1) 327 DO jk = 2, jpkm1 328 un_b(:,:) = un_b(:,:) + e3u_a(:,:,jk) * un(:,:,jk) * umask(:,:,jk) 329 ub_b(:,:) = ub_b(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk) 330 vn_b(:,:) = vn_b(:,:) + e3v_a(:,:,jk) * vn(:,:,jk) * vmask(:,:,jk) 331 vb_b(:,:) = vb_b(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk) 339 332 END DO 340 !341 !342 u n_b(:,:) = un_b(:,:) * hur_a(:,:)343 v n_b(:,:) = vn_b(:,:) * hvr_a(:,:)344 ub_b(:,:) = ub_b(:,:) * hur_b(:,:)345 vb_b(:,:) = vb_b(:,:) * hvr_b(:,:)346 !347 !348 333 un_b(:,:) = un_b(:,:) * r1_hu_a(:,:) 334 vn_b(:,:) = vn_b(:,:) * r1_hv_a(:,:) 335 ub_b(:,:) = ub_b(:,:) * r1_hu_b(:,:) 336 vb_b(:,:) = vb_b(:,:) * r1_hv_b(:,:) 337 ! 338 IF( .NOT.ln_dynspg_ts ) THEN ! output the barotropic currents 339 CALL iom_put( "ubar", un_b(:,:) ) 340 CALL iom_put( "vbar", vn_b(:,:) ) 341 ENDIF 349 342 IF( l_trddyn ) THEN ! 3D output: asselin filter trends on momentum 350 343 zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * z1_2dt … … 355 348 IF(ln_ctl) CALL prt_ctl( tab3d_1=un, clinfo1=' nxt - Un: ', mask1=umask, & 356 349 & tab3d_2=vn, clinfo2=' Vn: ' , mask2=vmask ) 357 ! 358 IF( ln_dynspg_ts ) CALL wrk_dealloc( jpi,jpj, zue, zve ) 359 IF( lk_vvl.AND.(.NOT.ln_dynadv_vec ) ) CALL wrk_dealloc( jpi,jpj,jpk, ze3u_f, ze3v_f ) 360 IF( l_trddyn ) CALL wrk_dealloc( jpi,jpj,jpk, zua, zva ) 350 ! 351 IF( ln_dynspg_ts ) CALL wrk_dealloc( jpi,jpj, zue, zve ) 352 IF( l_trddyn ) CALL wrk_dealloc( jpi,jpj,jpk, zua, zva ) 361 353 ! 362 354 IF( nn_timing == 1 ) CALL timing_stop('dyn_nxt')
Note: See TracChangeset
for help on using the changeset viewer.
