- Timestamp:
- 2016-07-19T10:38:35+02:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/NERC/dev_r5549_BDY_ZEROGRAD/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90
r5467 r6808 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 … … 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) … … 68 67 !! *** ROUTINE dyn_nxt *** 69 68 !! 70 !! ** Purpose : Compute the after horizontal velocity. Apply the boundary71 !! 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 72 71 !! by applying the Asselin filter on now fields and swapping 73 72 !! the fields. 74 73 !! 75 !! ** Method : * After velocity is compute using a leap-frog scheme: 76 !! (ua,va) = (ub,vb) + 2 rdt (ua,va) 77 !! Note that with flux form advection and variable volume layer 78 !! (lk_vvl=T), the leap-frog is applied on thickness weighted 79 !! velocity. 80 !! Note also that in filtered free surface (lk_dynspg_flt=T), 81 !! 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) 82 76 !! 83 77 !! * Apply lateral boundary conditions on after velocity … … 90 84 !! (ub,vb) = (un,vn) + atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ] 91 85 !! (un,vn) = (ua,va). 92 !! Note that with flux form advection and variable volume layer93 !! (lk_vvl=T), the time filter is applied on thickness weighted94 !! 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. 95 89 !! 96 90 !! ** Action : ub,vb filtered before horizontal velocity of next time-step … … 100 94 ! 101 95 INTEGER :: ji, jj, jk ! dummy loop indices 102 INTEGER :: iku, ikv ! local integers 103 #if ! defined key_dynspg_flt 104 REAL(wp) :: z2dt ! temporary scalar 105 #endif 106 REAL(wp) :: zue3a, zue3n, zue3b, zuf, zec ! local scalars 96 INTEGER :: ikt ! local integers 97 REAL(wp) :: zue3a, zue3n, zue3b, zuf, zcoef ! local scalars 107 98 REAL(wp) :: zve3a, zve3n, zve3b, zvf, z1_2dt ! - - 108 99 REAL(wp), POINTER, DIMENSION(:,:) :: zue, zve … … 112 103 IF( nn_timing == 1 ) CALL timing_start('dyn_nxt') 113 104 ! 114 CALL wrk_alloc( jpi,jpj,jpk, ze3u_f, ze3v_f, zua, zva)115 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) 116 107 ! 117 108 IF( kt == nit000 ) THEN … … 121 112 ENDIF 122 113 123 #if defined key_dynspg_flt 124 ! 125 ! Next velocity : Leap-frog time stepping already done in dynspg_flt.F routine 126 ! ------------- 127 128 ! Update after velocity on domain lateral boundaries (only local domain required) 129 ! -------------------------------------------------- 130 CALL lbc_lnk( ua, 'U', -1. ) ! local domain boundaries 131 CALL lbc_lnk( va, 'V', -1. ) 132 ! 133 #else 134 135 # if defined key_dynspg_exp 136 ! Next velocity : Leap-frog time stepping 137 ! ------------- 138 z2dt = 2. * rdt ! Euler or leap-frog time step 139 IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt 140 ! 141 IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) 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 142 123 DO jk = 1, jpkm1 143 ua(:,:,jk) = ( u b(:,:,jk) + z2dt * ua(:,:,jk) ) * umask(:,:,jk)144 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) 145 126 END DO 146 ELSE ! applied on thickness weighted velocity 147 DO jk = 1, jpkm1 148 ua(:,:,jk) = ( ub(:,:,jk) * fse3u_b(:,:,jk) & 149 & + z2dt * ua(:,:,jk) * fse3u_n(:,:,jk) ) & 150 & / fse3u_a(:,:,jk) * umask(:,:,jk) 151 va(:,:,jk) = ( vb(:,:,jk) * fse3v_b(:,:,jk) & 152 & + z2dt * va(:,:,jk) * fse3v_n(:,:,jk) ) & 153 & / fse3v_a(:,:,jk) * vmask(:,:,jk) 154 END DO 155 ENDIF 156 # endif 157 158 # if defined key_dynspg_ts 159 !!gm IF ( lk_dynspg_ts ) THEN .... 160 ! Ensure below that barotropic velocities match time splitting estimate 161 ! Compute actual transport and replace it with ts estimate at "after" time step 162 zue(:,:) = fse3u_a(:,:,1) * ua(:,:,1) * umask(:,:,1) 163 zve(:,:) = fse3v_a(:,:,1) * va(:,:,1) * vmask(:,:,1) 164 DO jk = 2, jpkm1 165 zue(:,:) = zue(:,:) + fse3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk) 166 zve(:,:) = zve(:,:) + fse3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk) 167 END DO 168 DO jk = 1, jpkm1 169 ua(:,:,jk) = ( ua(:,:,jk) - zue(:,:) * hur_a(:,:) + ua_b(:,:) ) * umask(:,:,jk) 170 va(:,:,jk) = ( va(:,:,jk) - zve(:,:) * hvr_a(:,:) + va_b(:,:) ) * vmask(:,:,jk) 171 END DO 172 173 IF (lk_dynspg_ts.AND.(.NOT.ln_bt_fw)) THEN 174 ! Remove advective velocity from "now velocities" 175 ! prior to asselin filtering 176 ! In the forward case, this is done below after asselin filtering 177 ! so that asselin contribution is removed at the same time 178 DO jk = 1, jpkm1 179 un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:) + un_b(:,:) )*umask(:,:,jk) 180 vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:) + vn_b(:,:) )*vmask(:,:,jk) 181 END DO 182 ENDIF 183 !!gm ENDIF 184 # 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 185 139 186 140 ! Update after velocity on domain lateral boundaries 187 141 ! -------------------------------------------------- 188 CALL lbc_lnk( ua, 'U', -1. ) !* local domain boundaries189 CALL lbc_lnk( va, 'V', -1. )190 !191 # if defined key_bdy192 ! !* BDY open boundaries193 IF( lk_bdy .AND. lk_dynspg_exp ) CALL bdy_dyn( kt )194 IF( lk_bdy .AND. lk_dynspg_ts ) CALL bdy_dyn( kt, dyn3d_only=.true. )195 196 !!$ Do we need a call to bdy_vol here??197 !198 # endif199 !200 142 # if defined key_agrif 201 143 CALL Agrif_dyn( kt ) !* AGRIF zoom boundaries 202 144 # endif 203 #endif 204 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 ! 205 158 IF( l_trddyn ) THEN ! prepare the atf trend computation + some diagnostics 206 159 z1_2dt = 1._wp / (2. * rdt) ! Euler or leap-frog time step … … 229 182 vn(:,:,jk) = va(:,:,jk) 230 183 END DO 231 IF (lk_vvl) THEN184 IF(.NOT.ln_linssh ) THEN 232 185 DO jk = 1, jpkm1 233 fse3t_b(:,:,jk) = fse3t_n(:,:,jk)234 fse3u_b(:,:,jk) = fse3u_n(:,:,jk)235 fse3v_b(:,:,jk) = fse3v_n(:,:,jk)236 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 237 190 ENDIF 238 191 ELSE !* Leap-Frog : Asselin filter and swap 239 192 ! ! =============! 240 IF( .NOT. lk_vvl ) THEN! Fixed volume !193 IF( ln_linssh ) THEN ! Fixed volume ! 241 194 ! ! =============! 242 195 DO jk = 1, jpkm1 … … 259 212 ! (used as a now filtered scale factor until the swap) 260 213 ! ---------------------------------------------------- 261 IF (lk_dynspg_ts.AND.ln_bt_fw) THEN 262 ! No asselin filtering on thicknesses if forward time splitting 263 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) 264 216 ELSE 265 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 266 220 ! Add volume filter correction: compatibility with tracer advection scheme 267 221 ! => time filter + conservation correction (only at the first level) 268 fse3t_b(:,:,1) = fse3t_b(:,:,1) - atfp * rdt * r1_rau0 * ( emp_b(:,:) - emp(:,:) & 269 & -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 ELSE ! if ice shelf melting 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) 233 END DO 234 END DO 235 END IF 270 236 ENDIF 271 237 ! 272 IF( ln_dynadv_vec ) THEN 273 ! Before scale factor at (u/v)-points 274 ! ----------------------------------- 275 CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3u_b(:,:,:), 'U' ) 276 CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3v_b(:,:,:), 'V' ) 277 ! Leap-Frog - Asselin filter and swap: applied on velocity 278 ! ----------------------------------- 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' ) 279 242 DO jk = 1, jpkm1 280 243 DO jj = 1, jpj … … 291 254 END DO 292 255 ! 293 ELSE 294 ! Temporary filtered scale factor at (u/v)-points (will become before scale factor) 295 !------------------------------------------------ 296 CALL dom_vvl_interpol( fse3t_b(:,:,:), ze3u_f, 'U' ) 297 CALL dom_vvl_interpol( fse3t_b(:,:,:), ze3v_f, 'V' ) 298 ! Leap-Frog - Asselin filter and swap: applied on thickness weighted velocity 299 ! ----------------------------------- =========================== 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' ) 300 262 DO jk = 1, jpkm1 301 263 DO jj = 1, jpj 302 264 DO ji = 1, jpi 303 zue3a = ua(ji,jj,jk) * fse3u_a(ji,jj,jk)304 zve3a = va(ji,jj,jk) * fse3v_a(ji,jj,jk)305 zue3n = un(ji,jj,jk) * fse3u_n(ji,jj,jk)306 zve3n = vn(ji,jj,jk) * fse3v_n(ji,jj,jk)307 zue3b = ub(ji,jj,jk) * fse3u_b(ji,jj,jk)308 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) 309 271 ! 310 272 zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n + zue3a ) ) / ze3u_f(ji,jj,jk) … … 318 280 END DO 319 281 END DO 320 fse3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1) ! e3u_b <-- filtered scale factor 321 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 ) 322 286 ENDIF 323 287 ! 324 288 ENDIF 325 289 ! 326 IF (lk_dynspg_ts.AND.ln_bt_fw) THEN290 IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN 327 291 ! Revert "before" velocities to time split estimate 328 292 ! Doing it here also means that asselin filter contribution is removed 329 zue(:,:) = fse3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1)330 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) 331 295 DO jk = 2, jpkm1 332 zue(:,:) = zue(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)333 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) 334 298 END DO 335 299 DO jk = 1, jpkm1 336 ub(:,:,jk) = ub(:,:,jk) - (zue(:,:) * hur(:,:) - un_b(:,:)) * umask(:,:,jk)337 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) 338 302 END DO 339 303 ENDIF … … 346 310 ! 347 311 ! 348 IF (lk_vvl) THEN349 hu_b(:,:) = 0.350 hv_b(:,:) = 0.351 DO jk = 1, jpkm1352 hu_b(:,:) = hu_b(:,:) + fse3u_b(:,:,jk) * umask(:,:,jk)353 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) 354 318 END DO 355 hur_b(:,:) = umask_i(:,:) / ( hu_b(:,:) + 1._wp - umask_i(:,:) ) 356 hvr_b(:,:) = vmask_i(:,:) / ( hv_b(:,:) + 1._wp - vmask_i(:,:) ) 357 ENDIF 358 ! 359 un_b(:,:) = 0._wp ; vn_b(:,:) = 0._wp 360 ub_b(:,:) = 0._wp ; vb_b(:,:) = 0._wp 361 ! 362 DO jk = 1, jpkm1 363 DO jj = 1, jpj 364 DO ji = 1, jpi 365 un_b(ji,jj) = un_b(ji,jj) + fse3u_a(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk) 366 vn_b(ji,jj) = vn_b(ji,jj) + fse3v_a(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk) 367 ! 368 ub_b(ji,jj) = ub_b(ji,jj) + fse3u_b(ji,jj,jk) * ub(ji,jj,jk) * umask(ji,jj,jk) 369 vb_b(ji,jj) = vb_b(ji,jj) + fse3v_b(ji,jj,jk) * vb(ji,jj,jk) * vmask(ji,jj,jk) 370 END DO 371 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) 372 332 END DO 373 !374 !375 u n_b(:,:) = un_b(:,:) * hur_a(:,:)376 v n_b(:,:) = vn_b(:,:) * hvr_a(:,:)377 ub_b(:,:) = ub_b(:,:) * hur_b(:,:)378 vb_b(:,:) = vb_b(:,:) * hvr_b(:,:)379 !380 !381 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 382 342 IF( l_trddyn ) THEN ! 3D output: asselin filter trends on momentum 383 343 zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * z1_2dt … … 389 349 & tab3d_2=vn, clinfo2=' Vn: ' , mask2=vmask ) 390 350 ! 391 CALL wrk_dealloc( jpi,jpj,jpk, ze3u_f, ze3v_f, zua, zva)392 IF( l k_dynspg_ts ) CALL wrk_dealloc( jpi,jpj, zue, zve)351 IF( ln_dynspg_ts ) CALL wrk_dealloc( jpi,jpj, zue, zve ) 352 IF( l_trddyn ) CALL wrk_dealloc( jpi,jpj,jpk, zua, zva ) 393 353 ! 394 354 IF( nn_timing == 1 ) CALL timing_stop('dyn_nxt')
Note: See TracChangeset
for help on using the changeset viewer.
