Changeset 11050 for NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/OCE/DYN/dynatf.F90
- Timestamp:
- 2019-05-24T10:27:58+02:00 (5 years ago)
- File:
-
- 1 moved
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/OCE/DYN/dynatf.F90
r11049 r11050 1 MODULE dyn nxt1 MODULE dynatf 2 2 !!========================================================================= 3 !! *** MODULE dyn nxt***4 !! Ocean dynamics: time stepping3 !! *** MODULE dynatf *** 4 !! Ocean dynamics: time filtering 5 5 !!========================================================================= 6 6 !! History : OPA ! 1987-02 (P. Andrich, D. L Hostis) Original code … … 20 20 !! 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 !! 4.1 ! 2019-05 (D. Storkey) Rename dynnxt -> dynatf. Now just does time filtering. 22 23 !!------------------------------------------------------------------------- 23 24 24 !!------------------------------------------------------------------------- 25 !! dyn_ nxt : obtain the next (after) horizontal velocity26 !!------------------------------------------------------------------------- 25 !!---------------------------------------------------------------------------------------------- 26 !! dyn_atf : apply Asselin time filtering to "now" velocities and vertical scale factors 27 !!---------------------------------------------------------------------------------------------- 27 28 USE oce ! ocean dynamics and tracers 28 29 USE dom_oce ! ocean space and time domain … … 50 51 #if defined key_agrif 51 52 USE agrif_oce_interp 52 #endif 53 #endif 53 54 54 55 IMPLICIT NONE 55 56 PRIVATE 56 57 57 PUBLIC dyn_ nxt! routine called by step.F9058 PUBLIC dyn_atf ! routine called by step.F90 58 59 59 60 !!---------------------------------------------------------------------- … … 64 65 CONTAINS 65 66 66 SUBROUTINE dyn_ nxt ( kt, Kbb, Kmm, Kaa)67 SUBROUTINE dyn_atf ( kt, Kbb, Kmm, Kaa, puu, pvv, pe3t, pe3u, pe3v ) 67 68 !!---------------------------------------------------------------------- 68 !! *** ROUTINE dyn_ nxt***69 !! *** ROUTINE dyn_atf *** 69 70 !! 71 !! !!!!!!!!!!!!!!!!!! REWRITE HEADER COMMENTS !!!!!!!!!!!!!!!!! 72 !! 70 73 !! ** Purpose : Finalize after horizontal velocity. Apply the boundary 71 74 !! condition on the after velocity, achieve the time stepping … … 83 86 !! * Apply the time filter applied and swap of the dynamics 84 87 !! arrays to start the next time step: 85 !! ( ub,vb) = (un,vn) + atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ]86 !! ( un,vn) = (ua,va).88 !! (puu(:,:,:,Kbb),pvv(:,:,:,Kbb)) = (puu(:,:,:,Kmm),pvv(:,:,:,Kmm)) + atfp [ (puu(:,:,:,Kbb),pvv(:,:,:,Kbb)) + (puu(:,:,:,Kaa),pvv(:,:,:,Kaa)) - 2 (puu(:,:,:,Kmm),pvv(:,:,:,Kmm)) ] 89 !! (puu(:,:,:,Kmm),pvv(:,:,:,Kmm)) = (puu(:,:,:,Kaa),pvv(:,:,:,Kaa)). 87 90 !! Note that with flux form advection and non linear free surface, 88 91 !! the time filter is applied on thickness weighted velocity. 89 !! As a result, dyn_nxt MUST be called after tra_nxt. 90 !! 91 !! ** Action : ub,vb filtered before horizontal velocity of next time-step 92 !! un,vn now horizontal velocity of next time-step 92 !! As a result, dyn_atf MUST be called after tra_nxt. 93 !! 94 !! ** Action : uu(Kmm),vv(Kmm) filtered now horizontal velocity 93 95 !!---------------------------------------------------------------------- 94 INTEGER, INTENT( in ) :: kt ! ocean time-step index 95 INTEGER, INTENT( in ) :: Kbb, Kmm, Kaa ! time level indices 96 INTEGER , INTENT(in ) :: kt ! ocean time-step index 97 INTEGER , INTENT(in ) :: Kbb, Kmm, Kaa ! before and after time level indices 98 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! velocities to be time filtered 99 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: pe3t, pe3u, pe3v ! scale factors to be time filtered 96 100 ! 97 101 INTEGER :: ji, jj, jk ! dummy loop indices 98 INTEGER :: ikt ! local integers 99 REAL(wp) :: zue3a, zue3n, zue3b, zuf, zcoef ! local scalars 100 REAL(wp) :: zve3a, zve3n, zve3b, zvf, z1_2dt ! - - 102 REAL(wp) :: zue3a, zue3n, zue3b, zcoef ! local scalars 103 REAL(wp) :: zve3a, zve3n, zve3b, z1_2dt ! - - 101 104 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zue, zve 102 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ze3 u_f, ze3v_f, zua, zva105 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ze3t_f, ze3u_f, ze3v_f, zua, zva 103 106 !!---------------------------------------------------------------------- 104 107 ! 105 IF( ln_timing ) CALL timing_start('dyn_ nxt')108 IF( ln_timing ) CALL timing_start('dyn_atf') 106 109 IF( ln_dynspg_ts ) ALLOCATE( zue(jpi,jpj) , zve(jpi,jpj) ) 107 110 IF( l_trddyn ) ALLOCATE( zua(jpi,jpj,jpk) , zva(jpi,jpj,jpk) ) … … 109 112 IF( kt == nit000 ) THEN 110 113 IF(lwp) WRITE(numout,*) 111 IF(lwp) WRITE(numout,*) 'dyn_ nxt : time stepping'114 IF(lwp) WRITE(numout,*) 'dyn_atf : Asselin time filtering' 112 115 IF(lwp) WRITE(numout,*) '~~~~~~~' 113 116 ENDIF … … 116 119 ! Ensure below that barotropic velocities match time splitting estimate 117 120 ! Compute actual transport and replace it with ts estimate at "after" time step 118 zue(:,:) = e3u_a(:,:,1) * ua(:,:,1) * umask(:,:,1)119 zve(:,:) = e3v_a(:,:,1) * va(:,:,1) * vmask(:,:,1)121 zue(:,:) = pe3u(:,:,1,Kaa) * puu(:,:,1,Kaa) * umask(:,:,1) 122 zve(:,:) = pe3v(:,:,1,Kaa) * pvv(:,:,1,Kaa) * vmask(:,:,1) 120 123 DO jk = 2, jpkm1 121 zue(:,:) = zue(:,:) + e3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)122 zve(:,:) = zve(:,:) + e3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)124 zue(:,:) = zue(:,:) + pe3u(:,:,jk,Kaa) * puu(:,:,jk,Kaa) * umask(:,:,jk) 125 zve(:,:) = zve(:,:) + pe3v(:,:,jk,Kaa) * pvv(:,:,jk,Kaa) * vmask(:,:,jk) 123 126 END DO 124 127 DO jk = 1, jpkm1 125 ua(:,:,jk) = ( ua(:,:,jk) - zue(:,:) * r1_hu_a(:,:) + ua_b(:,:) ) * umask(:,:,jk)126 va(:,:,jk) = ( va(:,:,jk) - zve(:,:) * r1_hv_a(:,:) + va_b(:,:) ) * vmask(:,:,jk)128 puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kaa) - zue(:,:) * r1_hu(:,:,Kaa) + uu_b(:,:,Kaa) ) * umask(:,:,jk) 129 pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kaa) - zve(:,:) * r1_hv(:,:,Kaa) + vv_b(:,:,Kaa) ) * vmask(:,:,jk) 127 130 END DO 128 131 ! … … 133 136 ! so that asselin contribution is removed at the same time 134 137 DO jk = 1, jpkm1 135 un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:)*r1_hu_n(:,:) + un_b(:,:) )*umask(:,:,jk)136 vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:)*r1_hv_n(:,:) + vn_b(:,:) )*vmask(:,:,jk)138 puu(:,:,jk,Kmm) = ( puu(:,:,jk,Kmm) - un_adv(:,:)*r1_hu(:,:,Kmm) + uu_b(:,:,Kmm) )*umask(:,:,jk) 139 pvv(:,:,jk,Kmm) = ( pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv(:,:,Kmm) + vv_b(:,:,Kmm) )*vmask(:,:,jk) 137 140 END DO 138 141 ENDIF … … 145 148 # endif 146 149 ! 147 CALL lbc_lnk_multi( 'dynnxt', ua, 'U', -1., va, 'V', -1. ) !* local domain boundaries150 CALL lbc_lnk_multi( 'dynnxt', puu(:,:,:,Kaa), 'U', -1., pvv(:,:,:,Kaa), 'V', -1. ) !* local domain boundaries 148 151 ! 149 152 ! !* BDY open boundaries 150 !! IMMERSE development : Following the general pattern for the new code we want to pass in the 151 !! velocities to bdy_dyn as arguments so here we use "uu" and "vv" even 152 !! though we haven't converted the velocity names in the rest of dynnxt.F90 153 !! because it will be completely rewritten. DS. 154 IF( ln_bdy .AND. ln_dynspg_exp ) CALL bdy_dyn( kt, Kbb, uu, vv, Kaa ) 155 IF( ln_bdy .AND. ln_dynspg_ts ) CALL bdy_dyn( kt, Kbb, uu, vv, Kaa, dyn3d_only=.true. ) 153 IF( ln_bdy .AND. ln_dynspg_exp ) CALL bdy_dyn( kt, Kbb, puu, pvv, Kaa ) 154 IF( ln_bdy .AND. ln_dynspg_ts ) CALL bdy_dyn( kt, Kbb, puu, pvv, Kaa, dyn3d_only=.true. ) 156 155 157 156 !!$ Do we need a call to bdy_vol here?? … … 162 161 ! 163 162 ! ! Kinetic energy and Conversion 164 IF( ln_KE_trd ) CALL trd_dyn( ua, va, jpdyn_ken, kt, Kmm )163 IF( ln_KE_trd ) CALL trd_dyn( puu(:,:,:,Kaa), pvv(:,:,:,Kaa), jpdyn_ken, kt, Kmm ) 165 164 ! 166 165 IF( ln_dyn_trd ) THEN ! 3D output: total momentum trends 167 zua(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * z1_2dt168 zva(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * z1_2dt166 zua(:,:,:) = ( puu(:,:,:,Kaa) - puu(:,:,:,Kbb) ) * z1_2dt 167 zva(:,:,:) = ( pvv(:,:,:,Kaa) - pvv(:,:,:,Kbb) ) * z1_2dt 169 168 CALL iom_put( "utrd_tot", zua ) ! total momentum trends, except the asselin time filter 170 169 CALL iom_put( "vtrd_tot", zva ) 171 170 ENDIF 172 171 ! 173 zua(:,:,:) = un(:,:,:) ! save the now velocity before the asselin filter174 zva(:,:,:) = vn(:,:,:) ! (caution: there will be a shift by 1 timestep in the172 zua(:,:,:) = puu(:,:,:,Kmm) ! save the now velocity before the asselin filter 173 zva(:,:,:) = pvv(:,:,:,Kmm) ! (caution: there will be a shift by 1 timestep in the 175 174 ! ! computation of the asselin filter trends) 176 175 ENDIF … … 178 177 ! Time filter and swap of dynamics arrays 179 178 ! ------------------------------------------ 180 IF( neuler == 0 .AND. kt == nit000 ) THEN !* Euler at first time-step: only swap181 DO jk = 1, jpkm1182 un(:,:,jk) = ua(:,:,jk) ! un <-- ua183 vn(:,:,jk) = va(:,:,jk)184 END DO185 IF( .NOT.ln_linssh ) THEN ! e3._b <-- e3._n186 !!gm BUG ???? I don't understand why it is not : e3._n <-- e3._a187 DO jk = 1, jpkm1188 ! e3t_b(:,:,jk) = e3t_n(:,:,jk)189 ! e3u_b(:,:,jk) = e3u_n(:,:,jk)190 ! e3v_b(:,:,jk) = e3v_n(:,:,jk)191 !192 e3t_n(:,:,jk) = e3t_a(:,:,jk)193 e3u_n(:,:,jk) = e3u_a(:,:,jk)194 e3v_n(:,:,jk) = e3v_a(:,:,jk)195 END DO196 !!gm BUG end197 ENDIF198 !199 179 200 ELSE !* Leap-Frog : Asselin filter and swap180 IF( .NOT.( neuler == 0 .AND. kt == nit000 ) ) THEN !* Leap-Frog : Asselin time filter 201 181 ! ! =============! 202 182 IF( ln_linssh ) THEN ! Fixed volume ! … … 205 185 DO jj = 1, jpj 206 186 DO ji = 1, jpi 207 zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) ) 208 zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) ) 209 ! 210 ub(ji,jj,jk) = zuf ! ub <-- filtered velocity 211 vb(ji,jj,jk) = zvf 212 un(ji,jj,jk) = ua(ji,jj,jk) ! un <-- ua 213 vn(ji,jj,jk) = va(ji,jj,jk) 187 puu(ji,jj,jk,Kmm) = puu(ji,jj,jk,Kmm) + atfp * ( puu(ji,jj,jk,Kbb) - 2._wp * puu(ji,jj,jk,Kmm) + puu(ji,jj,jk,Kaa) ) 188 pvv(ji,jj,jk,Kmm) = pvv(ji,jj,jk,Kmm) + atfp * ( pvv(ji,jj,jk,Kbb) - 2._wp * pvv(ji,jj,jk,Kmm) + pvv(ji,jj,jk,Kaa) ) 214 189 END DO 215 190 END DO … … 218 193 ELSE ! Variable volume ! 219 194 ! ! ================! 220 ! Before scale factor at t-points 221 ! (used as a now filtered scale factor until the swap) 195 ! Time-filtered scale factor at t-points 222 196 ! ---------------------------------------------------- 197 ALLOCATE( ze3t_f(jpi,jpj,jpk) ) 223 198 DO jk = 1, jpkm1 224 e3t_b(:,:,jk) = e3t_n(:,:,jk) + atfp * ( e3t_b(:,:,jk) - 2._wp * e3t_n(:,:,jk) + e3t_a(:,:,jk) )199 ze3t_f(:,:,jk) = pe3t(:,:,jk,Kmm) + atfp * ( pe3t(:,:,jk,Kbb) - 2._wp * pe3t(:,:,jk,Kmm) + pe3t(:,:,jk,Kaa) ) 225 200 END DO 226 201 ! Add volume filter correction: compatibility with tracer advection scheme … … 228 203 zcoef = atfp * rdt * r1_rau0 229 204 230 e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)205 ze3t_f(:,:,1) = ze3t_f(:,:,1) - zcoef * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1) 231 206 232 207 IF ( ln_rnf ) THEN 233 208 IF( ln_rnf_depth ) THEN 234 DO jk = 1, jpkm1 ! Deal with Rivers separ etely, as can be through depth too209 DO jk = 1, jpkm1 ! Deal with Rivers separately, as can be through depth too 235 210 DO jj = 1, jpj 236 211 DO ji = 1, jpi 237 212 IF( jk <= nk_rnf(ji,jj) ) THEN 238 e3t_b(ji,jj,jk) = e3t_b(ji,jj,jk) - zcoef * ( - rnf_b(ji,jj) + rnf(ji,jj) ) &239 & * ( e3t_n(ji,jj,jk) / h_rnf(ji,jj) ) * tmask(ji,jj,jk)213 ze3t_f(ji,jj,jk) = ze3t_f(ji,jj,jk) - zcoef * ( - rnf_b(ji,jj) + rnf(ji,jj) ) & 214 & * ( pe3t(ji,jj,jk,Kmm) / h_rnf(ji,jj) ) * tmask(ji,jj,jk) 240 215 ENDIF 241 216 ENDDO … … 243 218 ENDDO 244 219 ELSE 245 e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef * ( -rnf_b(:,:) + rnf(:,:))*tmask(:,:,1)220 ze3t_f(:,:,1) = ze3t_f(:,:,1) - zcoef * ( -rnf_b(:,:) + rnf(:,:))*tmask(:,:,1) 246 221 ENDIF 247 222 END IF … … 252 227 DO ji = 1, jpi 253 228 IF( misfkt(ji,jj) <=jk .and. jk < misfkb(ji,jj) ) THEN 254 e3t_b(ji,jj,jk) = e3t_b(ji,jj,jk) - zcoef * ( fwfisf_b(ji,jj) - fwfisf(ji,jj) ) &255 & * ( e3t_n(ji,jj,jk) * r1_hisf_tbl(ji,jj) ) * tmask(ji,jj,jk)229 ze3t_f(ji,jj,jk) = ze3t_f(ji,jj,jk) - zcoef * ( fwfisf_b(ji,jj) - fwfisf(ji,jj) ) & 230 & * ( pe3t(ji,jj,jk,Kmm) * r1_hisf_tbl(ji,jj) ) * tmask(ji,jj,jk) 256 231 ELSEIF ( jk==misfkb(ji,jj) ) THEN 257 e3t_b(ji,jj,jk) = e3t_b(ji,jj,jk) - zcoef * ( fwfisf_b(ji,jj) - fwfisf(ji,jj) ) &258 & * ( e3t_n(ji,jj,jk) * r1_hisf_tbl(ji,jj) ) * ralpha(ji,jj) * tmask(ji,jj,jk)232 ze3t_f(ji,jj,jk) = ze3t_f(ji,jj,jk) - zcoef * ( fwfisf_b(ji,jj) - fwfisf(ji,jj) ) & 233 & * ( pe3t(ji,jj,jk,Kmm) * r1_hisf_tbl(ji,jj) ) * ralpha(ji,jj) * tmask(ji,jj,jk) 259 234 ENDIF 260 235 END DO … … 263 238 END IF 264 239 ! 240 pe3t(:,:,1:jpkm1,Kmm) = ze3t_f(:,:,1:jpkm1) ! filtered scale factor at T-points 241 ! 265 242 IF( ln_dynadv_vec ) THEN ! Asselin filter applied on velocity 266 243 ! Before filtered scale factor at (u/v)-points 267 CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )268 CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )244 CALL dom_vvl_interpol( pe3t(:,:,:,Kmm), pe3u(:,:,:,Kmm), 'U' ) 245 CALL dom_vvl_interpol( pe3t(:,:,:,Kmm), pe3v(:,:,:,Kmm), 'V' ) 269 246 DO jk = 1, jpkm1 270 247 DO jj = 1, jpj 271 248 DO ji = 1, jpi 272 zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) ) 273 zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) ) 274 ! 275 ub(ji,jj,jk) = zuf ! ub <-- filtered velocity 276 vb(ji,jj,jk) = zvf 277 un(ji,jj,jk) = ua(ji,jj,jk) ! un <-- ua 278 vn(ji,jj,jk) = va(ji,jj,jk) 249 puu(ji,jj,jk,Kmm) = puu(ji,jj,jk,Kmm) + atfp * ( puu(ji,jj,jk,Kbb) - 2._wp * puu(ji,jj,jk,Kmm) + puu(ji,jj,jk,Kaa) ) 250 pvv(ji,jj,jk,Kmm) = pvv(ji,jj,jk,Kmm) + atfp * ( pvv(ji,jj,jk,Kbb) - 2._wp * pvv(ji,jj,jk,Kmm) + pvv(ji,jj,jk,Kaa) ) 279 251 END DO 280 252 END DO … … 284 256 ! 285 257 ALLOCATE( ze3u_f(jpi,jpj,jpk) , ze3v_f(jpi,jpj,jpk) ) 286 ! Beforefiltered scale factor at (u/v)-points stored in ze3u_f, ze3v_f287 CALL dom_vvl_interpol( e3t_b(:,:,:), ze3u_f, 'U' )288 CALL dom_vvl_interpol( e3t_b(:,:,:), ze3v_f, 'V' )258 ! Now filtered scale factor at (u/v)-points stored in ze3u_f, ze3v_f 259 CALL dom_vvl_interpol( pe3t(:,:,:,Kmm), ze3u_f, 'U' ) 260 CALL dom_vvl_interpol( pe3t(:,:,:,Kmm), ze3v_f, 'V' ) 289 261 DO jk = 1, jpkm1 290 262 DO jj = 1, jpj 291 263 DO ji = 1, jpi 292 zue3a = e3u_a(ji,jj,jk) * ua(ji,jj,jk)293 zve3a = e3v_a(ji,jj,jk) * va(ji,jj,jk)294 zue3n = e3u_n(ji,jj,jk) * un(ji,jj,jk)295 zve3n = e3v_n(ji,jj,jk) * vn(ji,jj,jk)296 zue3b = e3u_b(ji,jj,jk) * ub(ji,jj,jk)297 zve3b = e3v_b(ji,jj,jk) * vb(ji,jj,jk)264 zue3a = pe3u(ji,jj,jk,Kaa) * puu(ji,jj,jk,Kaa) 265 zve3a = pe3v(ji,jj,jk,Kaa) * pvv(ji,jj,jk,Kaa) 266 zue3n = pe3u(ji,jj,jk,Kmm) * puu(ji,jj,jk,Kmm) 267 zve3n = pe3v(ji,jj,jk,Kmm) * pvv(ji,jj,jk,Kmm) 268 zue3b = pe3u(ji,jj,jk,Kbb) * puu(ji,jj,jk,Kbb) 269 zve3b = pe3v(ji,jj,jk,Kbb) * pvv(ji,jj,jk,Kbb) 298 270 ! 299 zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n + zue3a ) ) / ze3u_f(ji,jj,jk) 300 zvf = ( zve3n + atfp * ( zve3b - 2._wp * zve3n + zve3a ) ) / ze3v_f(ji,jj,jk) 301 ! 302 ub(ji,jj,jk) = zuf ! ub <-- filtered velocity 303 vb(ji,jj,jk) = zvf 304 un(ji,jj,jk) = ua(ji,jj,jk) ! un <-- ua 305 vn(ji,jj,jk) = va(ji,jj,jk) 271 puu(ji,jj,jk,Kmm) = ( zue3n + atfp * ( zue3b - 2._wp * zue3n + zue3a ) ) / ze3u_f(ji,jj,jk) 272 pvv(ji,jj,jk,Kmm) = ( zve3n + atfp * ( zve3b - 2._wp * zve3n + zve3a ) ) / ze3v_f(ji,jj,jk) 306 273 END DO 307 274 END DO 308 275 END DO 309 e3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1) ! e3u_b <-- filtered scale factor310 e3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1)276 pe3u(:,:,1:jpkm1,Kmm) = ze3u_f(:,:,1:jpkm1) 277 pe3v(:,:,1:jpkm1,Kmm) = ze3v_f(:,:,1:jpkm1) 311 278 ! 312 279 DEALLOCATE( ze3u_f , ze3v_f ) 313 280 ENDIF 314 281 ! 282 DEALLOCATE( ze3t_f ) 315 283 ENDIF 316 284 ! 317 285 IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN 318 ! Revert "before" velocities to time split estimate286 ! Revert filtered "now" velocities to time split estimate 319 287 ! Doing it here also means that asselin filter contribution is removed 320 zue(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1)321 zve(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)288 zue(:,:) = pe3u(:,:,1,Kmm) * puu(:,:,1,Kmm) * umask(:,:,1) 289 zve(:,:) = pe3v(:,:,1,Kmm) * pvv(:,:,1,Kmm) * vmask(:,:,1) 322 290 DO jk = 2, jpkm1 323 zue(:,:) = zue(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)324 zve(:,:) = zve(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)291 zue(:,:) = zue(:,:) + pe3u(:,:,jk,Kmm) * puu(:,:,jk,Kmm) * umask(:,:,jk) 292 zve(:,:) = zve(:,:) + pe3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm) * vmask(:,:,jk) 325 293 END DO 326 294 DO jk = 1, jpkm1 327 ub(:,:,jk) = ub(:,:,jk) - (zue(:,:) * r1_hu_n(:,:) - un_b(:,:)) * umask(:,:,jk)328 vb(:,:,jk) = vb(:,:,jk) - (zve(:,:) * r1_hv_n(:,:) - vn_b(:,:)) * vmask(:,:,jk)295 puu(:,:,jk,Kmm) = puu(:,:,jk,Kmm) - (zue(:,:) * r1_hu(:,:,Kmm) - uu_b(:,:,Kmm)) * umask(:,:,jk) 296 pvv(:,:,jk,Kmm) = pvv(:,:,jk,Kmm) - (zve(:,:) * r1_hv(:,:,Kmm) - vv_b(:,:,Kmm)) * vmask(:,:,jk) 329 297 END DO 330 298 ENDIF 331 299 ! 332 ENDIF ! neuler =/0300 ENDIF ! neuler /= 0 333 301 ! 334 302 ! Set "now" and "before" barotropic velocities for next time step: … … 336 304 ! integration 337 305 ! 338 !339 306 IF(.NOT.ln_linssh ) THEN 340 hu _b(:,:) = e3u_b(:,:,1) * umask(:,:,1)341 hv _b(:,:) = e3v_b(:,:,1) * vmask(:,:,1)307 hu(:,:,Kmm) = pe3u(:,:,1,Kmm ) * umask(:,:,1) 308 hv(:,:,Kmm) = pe3v(:,:,1,Kmm ) * vmask(:,:,1) 342 309 DO jk = 2, jpkm1 343 hu _b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk)344 hv _b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk)310 hu(:,:,Kmm) = hu(:,:,Kmm) + pe3u(:,:,jk,Kmm ) * umask(:,:,jk) 311 hv(:,:,Kmm) = hv(:,:,Kmm) + pe3v(:,:,jk,Kmm ) * vmask(:,:,jk) 345 312 END DO 346 r1_hu _b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) )347 r1_hv _b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) )348 ENDIF 349 ! 350 u n_b(:,:) = e3u_a(:,:,1) * un(:,:,1) * umask(:,:,1)351 u b_b(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1)352 v n_b(:,:) = e3v_a(:,:,1) * vn(:,:,1) * vmask(:,:,1)353 v b_b(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)313 r1_hu(:,:,Kmm) = ssumask(:,:) / ( hu(:,:,Kmm) + 1._wp - ssumask(:,:) ) 314 r1_hv(:,:,Kmm) = ssvmask(:,:) / ( hv(:,:,Kmm) + 1._wp - ssvmask(:,:) ) 315 ENDIF 316 ! 317 uu_b(:,:,Kaa) = pe3u(:,:,1,Kaa) * puu(:,:,1,Kaa) * umask(:,:,1) 318 uu_b(:,:,Kmm) = pe3u(:,:,1,Kmm) * puu(:,:,1,Kmm) * umask(:,:,1) 319 vv_b(:,:,Kaa) = pe3v(:,:,1,Kaa) * pvv(:,:,1,Kaa) * vmask(:,:,1) 320 vv_b(:,:,Kmm) = pe3v(:,:,1,Kmm) * pvv(:,:,1,Kmm) * vmask(:,:,1) 354 321 DO jk = 2, jpkm1 355 u n_b(:,:) = un_b(:,:) + e3u_a(:,:,jk) * un(:,:,jk) * umask(:,:,jk)356 u b_b(:,:) = ub_b(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)357 v n_b(:,:) = vn_b(:,:) + e3v_a(:,:,jk) * vn(:,:,jk) * vmask(:,:,jk)358 v b_b(:,:) = vb_b(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)322 uu_b(:,:,Kaa) = uu_b(:,:,Kaa) + pe3u(:,:,jk,Kaa) * puu(:,:,jk,Kaa) * umask(:,:,jk) 323 uu_b(:,:,Kmm) = uu_b(:,:,Kmm) + pe3u(:,:,jk,Kmm) * puu(:,:,jk,Kmm) * umask(:,:,jk) 324 vv_b(:,:,Kaa) = vv_b(:,:,Kaa) + pe3v(:,:,jk,Kaa) * pvv(:,:,jk,Kaa) * vmask(:,:,jk) 325 vv_b(:,:,Kmm) = vv_b(:,:,Kmm) + pe3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm) * vmask(:,:,jk) 359 326 END DO 360 u n_b(:,:) = un_b(:,:) * r1_hu_a(:,:)361 v n_b(:,:) = vn_b(:,:) * r1_hv_a(:,:)362 u b_b(:,:) = ub_b(:,:) * r1_hu_b(:,:)363 v b_b(:,:) = vb_b(:,:) * r1_hv_b(:,:)327 uu_b(:,:,Kaa) = uu_b(:,:,Kaa) * r1_hu(:,:,Kaa) 328 vv_b(:,:,Kaa) = vv_b(:,:,Kaa) * r1_hv(:,:,Kaa) 329 uu_b(:,:,Kmm) = uu_b(:,:,Kmm) * r1_hu(:,:,Kmm) 330 vv_b(:,:,Kmm) = vv_b(:,:,Kmm) * r1_hv(:,:,Kmm) 364 331 ! 365 332 IF( .NOT.ln_dynspg_ts ) THEN ! output the barotropic currents 366 CALL iom_put( "ubar", u n_b(:,:) )367 CALL iom_put( "vbar", v n_b(:,:) )333 CALL iom_put( "ubar", uu_b(:,:,Kmm) ) 334 CALL iom_put( "vbar", vv_b(:,:,Kmm) ) 368 335 ENDIF 369 336 IF( l_trddyn ) THEN ! 3D output: asselin filter trends on momentum 370 zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * z1_2dt371 zva(:,:,:) = ( vb(:,:,:) - zva(:,:,:) ) * z1_2dt337 zua(:,:,:) = ( puu(:,:,:,Kmm) - zua(:,:,:) ) * z1_2dt 338 zva(:,:,:) = ( pvv(:,:,:,Kmm) - zva(:,:,:) ) * z1_2dt 372 339 CALL trd_dyn( zua, zva, jpdyn_atf, kt, Kmm ) 373 340 ENDIF 374 341 ! 375 IF(ln_ctl) CALL prt_ctl( tab3d_1= un, clinfo1=' nxt - Un: ', mask1=umask, &376 & tab3d_2= vn, clinfo2=' Vn: ' , mask2=vmask )342 IF(ln_ctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Kaa), clinfo1=' nxt - puu(:,:,:,Kaa): ', mask1=umask, & 343 & tab3d_2=pvv(:,:,:,Kaa), clinfo2=' pvv(:,:,:,Kaa): ' , mask2=vmask ) 377 344 ! 378 345 IF( ln_dynspg_ts ) DEALLOCATE( zue, zve ) 379 346 IF( l_trddyn ) DEALLOCATE( zua, zva ) 380 IF( ln_timing ) CALL timing_stop('dyn_ nxt')381 ! 382 END SUBROUTINE dyn_ nxt347 IF( ln_timing ) CALL timing_stop('dyn_atf') 348 ! 349 END SUBROUTINE dyn_atf 383 350 384 351 !!========================================================================= 385 END MODULE dyn nxt352 END MODULE dynatf
Note: See TracChangeset
for help on using the changeset viewer.