- Timestamp:
- 2015-12-16T10:25:22+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90
r5930 r6060 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 … … 96 95 INTEGER :: ji, jj, jk ! dummy loop indices 97 96 INTEGER :: iku, ikv ! local integers 98 REAL(wp) :: zue3a, zue3n, zue3b, zuf 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 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) 226 226 ELSE ! if ice shelf melting 227 DO jj = 1,jpj 228 DO ji = 1,jpi 227 zcoef = atfp * rdt * r1_rau0 228 DO jj = 1, jpj 229 DO ji = 1, jpi 229 230 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) 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) 234 234 END DO 235 235 END DO … … 237 237 ENDIF 238 238 ! 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 ! ----------------------------------- 239 IF( ln_dynadv_vec ) THEN ! Asselin filter applied on velocity 240 ! Before filtered scale factor at (u/v)-points 241 CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' ) 242 CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' ) 246 243 DO jk = 1, jpkm1 247 244 DO jj = 1, jpj … … 258 255 END DO 259 256 ! 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 ! ----------------------------------- =========================== 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 261 CALL dom_vvl_interpol( e3t_b(:,:,:), ze3u_f, 'U' ) 262 CALL dom_vvl_interpol( e3t_b(:,:,:), ze3v_f, 'V' ) 267 263 DO jk = 1, jpkm1 268 264 DO jj = 1, jpj 269 265 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)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) 276 272 ! 277 273 zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n + zue3a ) ) / ze3u_f(ji,jj,jk) … … 285 281 END DO 286 282 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) 283 e3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1) ! e3u_b <-- filtered scale factor 284 e3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1) 285 ! 286 CALL wrk_dealloc( jpi,jpj,jpk, ze3u_f, ze3v_f ) 289 287 ENDIF 290 288 ! 291 289 ENDIF 292 290 ! 293 IF (ln_dynspg_ts.AND.ln_bt_fw) THEN291 IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN 294 292 ! Revert "before" velocities to time split estimate 295 293 ! 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)294 zue(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1) 295 zve(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1) 298 296 DO jk = 2, jpkm1 299 zue(:,:) = zue(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)300 zve(:,:) = zve(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)297 zue(:,:) = zue(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk) 298 zve(:,:) = zve(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk) 301 299 END DO 302 300 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)301 ub(:,:,jk) = ub(:,:,jk) - (zue(:,:) * r1_hu_n(:,:) - un_b(:,:)) * umask(:,:,jk) 302 vb(:,:,jk) = vb(:,:,jk) - (zve(:,:) * r1_hv_n(:,:) - vn_b(:,:)) * vmask(:,:,jk) 305 303 END DO 306 304 ENDIF … … 313 311 ! 314 312 ! 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)313 IF(.NOT.ln_linssh ) THEN 314 hu_b(:,:) = e3u_b(:,:,1) * umask(:,:,1) 315 hv_b(:,:) = e3v_b(:,:,1) * vmask(:,:,1) 316 DO jk = 2, jpkm1 317 hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk) 318 hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk) 321 319 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 320 !!gm don't understand the use of umask_i .... 321 r1_hu_b(:,:) = umask_i(:,:) / ( hu_b(:,:) + 1._wp - umask_i(:,:) ) 322 r1_hv_b(:,:) = vmask_i(:,:) / ( hv_b(:,:) + 1._wp - vmask_i(:,:) ) 323 ENDIF 324 ! 325 un_b(:,:) = e3u_a(:,:,1) * un(:,:,1) * umask(:,:,1) 326 ub_b(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1) 327 vn_b(:,:) = e3v_a(:,:,1) * vn(:,:,1) * vmask(:,:,1) 328 vb_b(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1) 329 DO jk = 2, jpkm1 330 un_b(:,:) = un_b(:,:) + e3u_a(:,:,jk) * un(:,:,jk) * umask(:,:,jk) 331 ub_b(:,:) = ub_b(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk) 332 vn_b(:,:) = vn_b(:,:) + e3v_a(:,:,jk) * vn(:,:,jk) * vmask(:,:,jk) 333 vb_b(:,:) = vb_b(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk) 339 334 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 335 un_b(:,:) = un_b(:,:) * r1_hu_a(:,:) 336 vn_b(:,:) = vn_b(:,:) * r1_hv_a(:,:) 337 ub_b(:,:) = ub_b(:,:) * r1_hu_b(:,:) 338 vb_b(:,:) = vb_b(:,:) * r1_hv_b(:,:) 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 349 344 IF( l_trddyn ) THEN ! 3D output: asselin filter trends on momentum 350 345 zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * z1_2dt … … 355 350 IF(ln_ctl) CALL prt_ctl( tab3d_1=un, clinfo1=' nxt - Un: ', mask1=umask, & 356 351 & 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 ) 352 ! 353 IF( ln_dynspg_ts ) CALL wrk_dealloc( jpi,jpj, zue, zve ) 354 IF( l_trddyn ) CALL wrk_dealloc( jpi,jpj,jpk, zua, zva ) 361 355 ! 362 356 IF( nn_timing == 1 ) CALL timing_stop('dyn_nxt')
Note: See TracChangeset
for help on using the changeset viewer.
