- Timestamp:
- 2015-12-04T17:05:58+01:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90
r5866 r6004 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 trends 20 !! 3.6 ! 2014-04 (G. Madec) add the diagnostic of the time filter trends 21 !! 3.7 ! 2015-11 (J. Chanut) Free surface simplification 21 22 !!------------------------------------------------------------------------- 22 23 23 24 !!------------------------------------------------------------------------- 24 !! dyn_nxt : obtain the next (after) horizontal velocity25 !! dyn_nxt : obtain the next (after) horizontal velocity 25 26 !!------------------------------------------------------------------------- 26 USE oce 27 USE dom_oce 28 USE sbc_oce 29 USE phycst 30 USE dyn spg_oce ! type of surface pressure gradient31 USE dyn adv ! dynamics: vector invariant versus flux form32 USE domvvl 33 USE bdy_oce 34 USE bdydta 35 USE bdydyn 36 USE bdyvol 37 USE trd_oce 38 USE trddyn 39 USE trdken 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 … … 66 67 !! *** ROUTINE dyn_nxt *** 67 68 !! 68 !! ** Purpose : Compute the after horizontal velocity. Apply the boundary69 !! condition on the after velocity, achieve dthe time stepping69 !! ** Purpose : Finalize after horizontal velocity. Apply the boundary 70 !! condition on the after velocity, achieve the time stepping 70 71 !! by applying the Asselin filter on now fields and swapping 71 72 !! the fields. 72 73 !! 73 !! ** Method : * After velocity is compute using a leap-frog scheme: 74 !! (ua,va) = (ub,vb) + 2 rdt (ua,va) 75 !! Note that with flux form advection and non linear free surface, 76 !! the leap-frog is applied on thickness weighted velocity. 77 !! Note also that in filtered free surface (lk_dynspg_flt=T), 78 !! the time stepping has already been done in dynspg module 74 !! ** Method : * Ensure after velocities transport matches time splitting 75 !! estimate (ln_dynspg_ts=T) 79 76 !! 80 77 !! * Apply lateral boundary conditions on after velocity … … 89 86 !! Note that with flux form advection and non linear free surface, 90 87 !! the time filter is applied on thickness weighted velocity. 88 !! As a result, dyn_nxt MUST be called after tra_nxt. 91 89 !! 92 90 !! ** Action : ub,vb filtered before horizontal velocity of next time-step … … 97 95 INTEGER :: ji, jj, jk ! dummy loop indices 98 96 INTEGER :: iku, ikv ! local integers 99 #if ! defined key_dynspg_flt 100 REAL(wp) :: z2dt ! temporary scalar 101 #endif 102 REAL(wp) :: zue3a, zue3n, zue3b, zuf, zec ! local scalars 97 REAL(wp) :: zue3a, zue3n, zue3b, zuf, zcoef ! local scalars 103 98 REAL(wp) :: zve3a, zve3n, zve3b, zvf, z1_2dt ! - - 104 99 REAL(wp), POINTER, DIMENSION(:,:) :: zue, zve … … 108 103 IF( nn_timing == 1 ) CALL timing_start('dyn_nxt') 109 104 ! 110 CALL wrk_alloc( jpi,jpj,jpk, ze3u_f, ze3v_f, zua, zva)111 IF( l k_dynspg_ts ) CALL wrk_alloc( jpi,jpj, zue, zve)105 IF( ln_dynspg_ts ) CALL wrk_alloc( jpi,jpj, zue, zve) 106 IF( l_trddyn ) CALL wrk_alloc( jpi,jpj,jpk, zua, zva) 112 107 ! 113 108 IF( kt == nit000 ) THEN … … 117 112 ENDIF 118 113 119 #if defined key_dynspg_flt 120 ! 121 ! Next velocity : Leap-frog time stepping already done in dynspg_flt.F routine 122 ! ------------- 123 124 ! Update after velocity on domain lateral boundaries (only local domain required) 125 ! -------------------------------------------------- 126 CALL lbc_lnk( ua, 'U', -1. ) ! local domain boundaries 127 CALL lbc_lnk( va, 'V', -1. ) 128 ! 129 #else 130 131 # if defined key_dynspg_exp 132 ! Next velocity : Leap-frog time stepping 133 ! ------------- 134 z2dt = 2. * rdt ! Euler or leap-frog time step 135 IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt 136 ! 137 IF( ln_dynadv_vec .OR. ln_linssh ) THEN !== applied on velocity ==! 114 IF ( ln_dynspg_ts ) THEN 115 ! Ensure below that barotropic velocities match time splitting estimate 116 ! Compute actual transport and replace it with ts estimate at "after" time step 117 zue(:,:) = e3u_a(:,:,1) * ua(:,:,1) * umask(:,:,1) 118 zve(:,:) = e3v_a(:,:,1) * va(:,:,1) * vmask(:,:,1) 119 DO jk = 2, jpkm1 120 zue(:,:) = zue(:,:) + e3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk) 121 zve(:,:) = zve(:,:) + e3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk) 122 END DO 138 123 DO jk = 1, jpkm1 139 ua(:,:,jk) = ( u b(:,:,jk) + z2dt * ua(:,:,jk) ) * umask(:,:,jk)140 va(:,:,jk) = ( v b(:,:,jk) + z2dt * va(:,:,jk) ) * 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) 141 126 END DO 142 ELSE !== applied on thickness weighted velocity ==! 143 DO jk = 1, jpkm1 144 ua(:,:,jk) = ( ub(:,:,jk) * e3u_b(:,:,jk) & 145 & + z2dt * ua(:,:,jk) * e3u_n(:,:,jk) ) / e3u_a(:,:,jk) * umask(:,:,jk) 146 va(:,:,jk) = ( vb(:,:,jk) * e3v_b(:,:,jk) & 147 & + z2dt * va(:,:,jk) * e3v_n(:,:,jk) ) / e3v_a(:,:,jk) * vmask(:,:,jk) 148 END DO 149 ENDIF 150 # endif 151 152 # if defined key_dynspg_ts 153 !!gm IF ( lk_dynspg_ts ) THEN .... 154 ! Ensure below that barotropic velocities match time splitting estimate 155 ! Compute actual transport and replace it with ts estimate at "after" time step 156 zue(:,:) = e3u_a(:,:,1) * ua(:,:,1) * umask(:,:,1) 157 zve(:,:) = e3v_a(:,:,1) * va(:,:,1) * vmask(:,:,1) 158 DO jk = 2, jpkm1 159 zue(:,:) = zue(:,:) + e3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk) 160 zve(:,:) = zve(:,:) + e3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk) 161 END DO 162 DO jk = 1, jpkm1 163 ua(:,:,jk) = ( ua(:,:,jk) - zue(:,:) * r1_hu_a(:,:) + ua_b(:,:) ) * umask(:,:,jk) 164 va(:,:,jk) = ( va(:,:,jk) - zve(:,:) * r1_hv_a(:,:) + va_b(:,:) ) * vmask(:,:,jk) 165 END DO 166 167 IF (lk_dynspg_ts.AND.(.NOT.ln_bt_fw)) THEN 168 ! Remove advective velocity from "now velocities" 169 ! prior to asselin filtering 170 ! In the forward case, this is done below after asselin filtering 171 ! so that asselin contribution is removed at the same time 172 DO jk = 1, jpkm1 173 un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:) + un_b(:,:) )*umask(:,:,jk) 174 vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:) + vn_b(:,:) )*vmask(:,:,jk) 175 END DO 176 ENDIF 177 !!gm ENDIF 178 # endif 127 ! 128 IF( .NOT.ln_bt_fw ) THEN 129 ! Remove advective velocity from "now velocities" 130 ! prior to asselin filtering 131 ! In the forward case, this is done below after asselin filtering 132 ! so that asselin contribution is removed at the same time 133 DO jk = 1, jpkm1 134 un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:) + un_b(:,:) )*umask(:,:,jk) 135 vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:) + vn_b(:,:) )*vmask(:,:,jk) 136 END DO 137 ENDIF 138 ENDIF 179 139 180 140 ! Update after velocity on domain lateral boundaries 181 141 ! -------------------------------------------------- 182 CALL lbc_lnk( ua, 'U', -1. ) !* local domain boundaries183 CALL lbc_lnk( va, 'V', -1. )184 !185 # if defined key_bdy186 ! !* BDY open boundaries187 IF( lk_bdy .AND. lk_dynspg_exp ) CALL bdy_dyn( kt )188 IF( lk_bdy .AND. lk_dynspg_ts ) CALL bdy_dyn( kt, dyn3d_only=.true. )189 190 !!$ Do we need a call to bdy_vol here??191 !192 # endif193 !194 142 # if defined key_agrif 195 143 CALL Agrif_dyn( kt ) !* AGRIF zoom boundaries 196 144 # endif 197 #endif 198 145 ! 146 CALL lbc_lnk( ua, 'U', -1. ) !* local domain boundaries 147 CALL lbc_lnk( va, 'V', -1. ) 148 ! 149 # if defined key_bdy 150 ! !* BDY open boundaries 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. ) 153 154 !!$ Do we need a call to bdy_vol here?? 155 ! 156 # endif 157 ! 199 158 IF( l_trddyn ) THEN ! prepare the atf trend computation + some diagnostics 200 159 z1_2dt = 1._wp / (2. * rdt) ! Euler or leap-frog time step … … 253 212 ! (used as a now filtered scale factor until the swap) 254 213 ! ---------------------------------------------------- 255 IF (lk_dynspg_ts.AND.ln_bt_fw) THEN 256 ! No asselin filtering on thicknesses if forward time splitting 257 e3t_b(:,:,:) = e3t_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) 258 216 ELSE 259 e3t_b(:,:,:) = e3t_n(:,:,:) + atfp * ( e3t_b(:,:,:) - 2._wp * e3t_n(:,:,:) + e3t_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 260 220 ! Add volume filter correction: compatibility with tracer advection scheme 261 221 ! => time filter + conservation correction (only at the first level) 262 IF ( nn_isf == 0) THEN ! if no ice shelf melting 263 e3t_b(:,:,1) = e3t_b(:,:,1) - atfp * rdt * r1_rau0 * ( emp_b(:,:) - emp(:,:) & 264 & -rnf_b(:,:) + rnf(:,:) ) * tmask(:,:,1) 222 IF( nn_isf == 0) THEN ! if no ice shelf melting 223 zcoef = atfp * rdt * r1_rau0 224 e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef * ( emp_b(:,:) - emp(:,:) & 225 & - rnf_b(:,:) + rnf(:,:) ) * tmask(:,:,1) 265 226 ELSE ! if ice shelf melting 266 DO jj = 1,jpj 267 DO ji = 1,jpi 227 zcoef = atfp * rdt * r1_rau0 228 DO jj = 1, jpj 229 DO ji = 1, jpi 268 230 jk = mikt(ji,jj) 269 e3t_b(ji,jj,jk) = e3t_b(ji,jj,jk) - atfp * rdt * r1_rau0 & 270 & * ( (emp_b(ji,jj) - emp(ji,jj) ) & 271 & - (rnf_b(ji,jj) - rnf(ji,jj) ) & 272 & + (fwfisf_b(ji,jj) - fwfisf(ji,jj)) ) * tmask(ji,jj,jk) 231 e3t_b(ji,jj,jk) = e3t_b(ji,jj,jk) - zcoef * ( 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) 273 234 END DO 274 235 END DO … … 276 237 ENDIF 277 238 ! 278 IF( ln_dynadv_vec ) THEN 279 ! Before scale factor at (u/v)-points 280 ! ----------------------------------- 239 IF( ln_dynadv_vec ) THEN ! Asselin filter applied on velocity 240 ! Before filtered scale factor at (u/v)-points 281 241 CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' ) 282 242 CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' ) 283 ! Leap-Frog - Asselin filter and swap: applied on velocity284 ! -----------------------------------285 243 DO jk = 1, jpkm1 286 244 DO jj = 1, jpj … … 297 255 END DO 298 256 ! 299 ELSE 300 ! Temporary filtered scale factor at (u/v)-points (will become before scale factor) 301 !------------------------------------------------ 257 ELSE ! Asselin filter applied on thickness weighted velocity 258 ! 259 CALL wrk_alloc( jpi,jpj,jpk, ze3u_f, ze3v_f ) 260 ! Before filtered scale factor at (u/v)-points stored in ze3u_f, ze3v_f 302 261 CALL dom_vvl_interpol( e3t_b(:,:,:), ze3u_f, 'U' ) 303 262 CALL dom_vvl_interpol( e3t_b(:,:,:), ze3v_f, 'V' ) 304 ! Leap-Frog - Asselin filter and swap: applied on thickness weighted velocity305 ! ----------------------------------- ===========================306 263 DO jk = 1, jpkm1 307 264 DO jj = 1, jpj 308 265 DO ji = 1, jpi 309 zue3a = ua(ji,jj,jk) * e3u_a(ji,jj,jk)310 zve3a = va(ji,jj,jk) * e3v_a(ji,jj,jk)311 zue3n = un(ji,jj,jk) * e3u_n(ji,jj,jk)312 zve3n = vn(ji,jj,jk) * e3v_n(ji,jj,jk)313 zue3b = ub(ji,jj,jk) * e3u_b(ji,jj,jk)314 zve3b = vb(ji,jj,jk) * e3v_b(ji,jj,jk)266 zue3a = e3u_a(ji,jj,jk) * ua(ji,jj,jk) 267 zve3a = e3v_a(ji,jj,jk) * va(ji,jj,jk) 268 zue3n = e3u_n(ji,jj,jk) * un(ji,jj,jk) 269 zve3n = e3v_n(ji,jj,jk) * vn(ji,jj,jk) 270 zue3b = e3u_b(ji,jj,jk) * ub(ji,jj,jk) 271 zve3b = e3v_b(ji,jj,jk) * vb(ji,jj,jk) 315 272 ! 316 273 zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n + zue3a ) ) / ze3u_f(ji,jj,jk) … … 324 281 END DO 325 282 END DO 326 e3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1) ! e3u_b <-- filtered scale factor283 e3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1) ! e3u_b <-- filtered scale factor 327 284 e3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1) 285 ! 286 CALL wrk_dealloc( jpi,jpj,jpk, ze3u_f, ze3v_f ) 328 287 ENDIF 329 288 ! 330 289 ENDIF 331 290 ! 332 IF (lk_dynspg_ts.AND.ln_bt_fw) THEN291 IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN 333 292 ! Revert "before" velocities to time split estimate 334 293 ! Doing it here also means that asselin filter contribution is removed … … 364 323 ENDIF 365 324 ! 366 un_b(:,:) = 0._wp ; vn_b(:,:) = 0._wp 367 ub_b(:,:) = 0._wp ; vb_b(:,:) = 0._wp 368 DO jk = 1, jpkm1 369 DO jj = 1, jpj 370 DO ji = 1, jpi 371 un_b(ji,jj) = un_b(ji,jj) + e3u_a(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk) 372 vn_b(ji,jj) = vn_b(ji,jj) + e3v_a(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk) 373 ! 374 ub_b(ji,jj) = ub_b(ji,jj) + e3u_b(ji,jj,jk) * ub(ji,jj,jk) * umask(ji,jj,jk) 375 vb_b(ji,jj) = vb_b(ji,jj) + e3v_b(ji,jj,jk) * vb(ji,jj,jk) * vmask(ji,jj,jk) 376 END DO 377 END DO 325 un_b(:,:) = e3u_a(:,:,jk) * un(:,:,1) * umask(:,:,1) 326 ub_b(:,:) = e3u_b(:,:,jk) * ub(:,:,1) * umask(:,:,1) 327 vn_b(:,:) = e3v_a(:,:,jk) * vn(:,:,1) * vmask(:,:,1) 328 vb_b(:,:) = e3v_b(:,:,jk) * vb(:,:,1) * vmask(:,:,1) 329 DO jk = 2, jpkm1 330 un_b(:,:) = un_b(:,:) + e3u_a(:,:,jk) * un(:,:,jk) * umask(ji,jj,jk) 331 ub_b(:,:) = ub_b(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(ji,jj,jk) 332 vn_b(:,:) = vn_b(:,:) + e3v_a(:,:,jk) * vn(:,:,jk) * vmask(ji,jj,jk) 333 vb_b(:,:) = vb_b(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(ji,jj,jk) 378 334 END DO 379 335 un_b(:,:) = un_b(:,:) * r1_hu_a(:,:) … … 382 338 vb_b(:,:) = vb_b(:,:) * r1_hv_b(:,:) 383 339 ! 340 IF( .NOT.ln_dynspg_ts ) THEN ! output the barotropic currents 341 CALL iom_put( "ubar", un_b(:,:) ) 342 CALL iom_put( "vbar", vn_b(:,:) ) 343 ENDIF 384 344 IF( l_trddyn ) THEN ! 3D output: asselin filter trends on momentum 385 345 zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * z1_2dt … … 391 351 & tab3d_2=vn, clinfo2=' Vn: ' , mask2=vmask ) 392 352 ! 393 CALL wrk_dealloc( jpi,jpj,jpk, ze3u_f, ze3v_f, zua, zva)394 IF( l k_dynspg_ts ) CALL wrk_dealloc( jpi,jpj, zue, zve)353 IF( ln_dynspg_ts ) CALL wrk_dealloc( jpi,jpj, zue, zve ) 354 IF( l_trddyn ) CALL wrk_dealloc( jpi,jpj,jpk, zua, zva ) 395 355 ! 396 356 IF( nn_timing == 1 ) CALL timing_stop('dyn_nxt')
Note: See TracChangeset
for help on using the changeset viewer.
