Changeset 10799
- Timestamp:
- 2019-03-25T09:38:58+01:00 (6 years ago)
- Location:
- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE
- Files:
-
- 3 edited
-
DOM/domvvl.F90 (modified) (7 diffs)
-
DYN/dynnxt.F90 (modified) (17 diffs)
-
step.F90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM/domvvl.F90
r10795 r10799 563 563 564 564 565 SUBROUTINE dom_vvl_sf_swp( kt , kNm1, kNnn, kNm1_2lev, kNnn_2lev)565 SUBROUTINE dom_vvl_sf_swp( kt ) 566 566 !!---------------------------------------------------------------------- 567 567 !! *** ROUTINE dom_vvl_sf_swp *** … … 588 588 !!---------------------------------------------------------------------- 589 589 INTEGER, INTENT( in ) :: kt ! time step 590 INTEGER, INTENT( in ) :: kNm1, kNnn, kNm1_2lev, kNnn_2lev ! time level indices591 590 ! 592 591 INTEGER :: ji, jj, jk ! dummy loop indices … … 606 605 ! Time filter and swap of scale factors 607 606 ! ===================================== 608 ! - ML - e3(t/u/v)_b are al ready computed in dynnxt.607 ! - ML - e3(t/u/v)_b are allready computed in dynnxt. 609 608 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 610 609 IF( neuler == 0 .AND. kt == nit000 ) THEN … … 616 615 tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:) 617 616 ENDIF 617 gdept_b(:,:,:) = gdept_n(:,:,:) 618 gdepw_b(:,:,:) = gdepw_n(:,:,:) 619 620 e3t_n(:,:,:) = e3t_a(:,:,:) 621 e3u_n(:,:,:) = e3u_a(:,:,:) 622 e3v_n(:,:,:) = e3v_a(:,:,:) 618 623 619 624 ! Compute all missing vertical scale factor and depths … … 624 629 ! - JC - hu_b, hv_b, hur_b, hvr_b also 625 630 626 CALL dom_vvl_interpol( e3u (:,:,:,kNnn), e3f(:,:,:), 'F' )631 CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' ) 627 632 628 633 ! Vertical scale factor interpolations 629 CALL dom_vvl_interpol( e3t (:,:,:,kNnn), e3w(:,:,:,kNnn_2lev), 'W' )630 CALL dom_vvl_interpol( e3u (:,:,:,kNnn), e3uw(:,:,:,kNnn_2lev), 'UW' )631 CALL dom_vvl_interpol( e3v (:,:,:,kNnn), e3vw(:,:,:,kNnn_2lev), 'VW' )632 CALL dom_vvl_interpol( e3t (:,:,:,kNm1), e3w(:,:,:,kNm1_2lev), 'W' )633 CALL dom_vvl_interpol( e3u (:,:,:,kNm1), e3uw(:,:,:,kNm1_2lev), 'UW' )634 CALL dom_vvl_interpol( e3v (:,:,:,kNm1), e3vw(:,:,:,kNm1_2lev), 'VW' )634 CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n(:,:,:), 'W' ) 635 CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' ) 636 CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' ) 637 CALL dom_vvl_interpol( e3t_b(:,:,:), e3w_b(:,:,:), 'W' ) 638 CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' ) 639 CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' ) 635 640 636 641 ! t- and w- points depth (set the isf depth as it is in the initial step) 637 gdept (:,:,1,kNnn_2lev) = 0.5_wp * e3w(:,:,1,kNnn_2lev)638 gdepw (:,:,1,kNnn_2lev) = 0.0_wp639 gde3w (:,:,1) = gdept(:,:,1,kNnn_2lev) - sshn(:,:)642 gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1) 643 gdepw_n(:,:,1) = 0.0_wp 644 gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:) 640 645 DO jk = 2, jpk 641 646 DO jj = 1,jpj … … 644 649 ! 1 for jk = mikt 645 650 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 646 gdepw (ji,jj,jk,kNnn_2lev) = gdepw(ji,jj,jk-1,kNnn_2lev) + e3t(ji,jj,jk-1,kNnn)647 gdept (ji,jj,jk,kNnn_2lev) = zcoef * ( gdepw(ji,jj,jk,kNnn_2lev ) + 0.5 * e3w(ji,jj,jk,kNnn_2lev) ) &648 & + (1-zcoef) * ( gdept (ji,jj,jk-1,kNnn_2lev) + e3w(ji,jj,jk,kNnn_2lev) )649 gde3w (ji,jj,jk) = gdept(ji,jj,jk,kNnn_2lev) - sshn(ji,jj)651 gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1) 652 gdept_n(ji,jj,jk) = zcoef * ( gdepw_n(ji,jj,jk ) + 0.5 * e3w_n(ji,jj,jk) ) & 653 & + (1-zcoef) * ( gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk) ) 654 gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj) 650 655 END DO 651 656 END DO … … 654 659 ! Local depth and Inverse of the local depth of the water 655 660 ! ------------------------------------------------------- 656 ! 657 ht(:,:) = e3t(:,:,1,kNnn) * tmask(:,:,1) 661 hu_n(:,:) = hu_a(:,:) ; r1_hu_n(:,:) = r1_hu_a(:,:) 662 hv_n(:,:) = hv_a(:,:) ; r1_hv_n(:,:) = r1_hv_a(:,:) 663 ! 664 ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1) 658 665 DO jk = 2, jpkm1 659 ht (:,:) = ht(:,:) + e3t(:,:,jk,kNnn) * tmask(:,:,jk)666 ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 660 667 END DO 661 668 -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynnxt.F90
r10795 r10799 64 64 CONTAINS 65 65 66 SUBROUTINE dyn_nxt ( kt , kNm1, kNp1, puu, pvv, pe3t, pe3u, pe3v)66 SUBROUTINE dyn_nxt ( kt ) 67 67 !!---------------------------------------------------------------------- 68 68 !! *** ROUTINE dyn_nxt *** … … 92 92 !! un,vn now horizontal velocity of next time-step 93 93 !!---------------------------------------------------------------------- 94 INTEGER, INTENT( in ) :: kt ! ocean time-step index 95 INTEGER, INTENT( in ) :: kNm1, kNp1 ! before and after time level indices 96 REAL(wp), INTENT( inout), DIMENSION(jpi,jpj,jpk) :: puu, pvv ! now velocities to be time filtered 97 REAL(wp), INTENT( inout), DIMENSION(jpi,jpj,jpk) :: pe3t, pe3u, pe3v ! now scale factors to be time filtered 94 INTEGER, INTENT( in ) :: kt ! ocean time-step index 98 95 ! 99 96 INTEGER :: ji, jj, jk ! dummy loop indices 100 97 INTEGER :: ikt ! local integers 101 REAL(wp) :: zue3a, zue3n, zue3b, z coef ! local scalars102 REAL(wp) :: zve3a, zve3n, zve3b, z 1_2dt ! - -98 REAL(wp) :: zue3a, zue3n, zue3b, zuf, zcoef ! local scalars 99 REAL(wp) :: zve3a, zve3n, zve3b, zvf, z1_2dt ! - - 103 100 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zue, zve 104 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ze3 t_f, ze3u_f, ze3v_f, zua, zva101 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ze3u_f, ze3v_f, zua, zva 105 102 !!---------------------------------------------------------------------- 106 103 ! … … 118 115 ! Ensure below that barotropic velocities match time splitting estimate 119 116 ! Compute actual transport and replace it with ts estimate at "after" time step 120 zue(:,:) = e3u (:,:,1,kNp1) * uu(:,:,1,kNp1) * umask(:,:,1)121 zve(:,:) = e3v (:,:,1,kNp1) * vv(:,:,1,kNp1) * vmask(:,:,1)117 zue(:,:) = e3u_a(:,:,1) * ua(:,:,1) * umask(:,:,1) 118 zve(:,:) = e3v_a(:,:,1) * va(:,:,1) * vmask(:,:,1) 122 119 DO jk = 2, jpkm1 123 zue(:,:) = zue(:,:) + e3u (:,:,jk,kNp1) * uu(:,:,jk,kNp1) * umask(:,:,jk)124 zve(:,:) = zve(:,:) + e3v (:,:,jk,kNp1) * vv(:,:,jk,kNp1) * vmask(:,:,jk)120 zue(:,:) = zue(:,:) + e3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk) 121 zve(:,:) = zve(:,:) + e3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk) 125 122 END DO 126 123 DO jk = 1, jpkm1 127 u u(:,:,jk,kNp1) = ( uu(:,:,jk,kNp1) - zue(:,:) * r1_hu_a(:,:) + ua_b(:,:) ) * umask(:,:,jk)128 v v(:,:,jk,kNp1) = ( vv(:,:,jk,kNp1) - zve(:,:) * r1_hv_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) 129 126 END DO 130 127 ! … … 135 132 ! so that asselin contribution is removed at the same time 136 133 DO jk = 1, jpkm1 137 puu(:,:,jk) = ( puu(:,:,jk) - un_adv(:,:)*r1_hu_n(:,:) + un_b(:,:) )*umask(:,:,jk)138 pvv(:,:,jk) = ( pvv(:,:,jk) - vn_adv(:,:)*r1_hv_n(:,:) + vn_b(:,:) )*vmask(:,:,jk)134 un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:)*r1_hu_n(:,:) + un_b(:,:) )*umask(:,:,jk) 135 vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:)*r1_hv_n(:,:) + vn_b(:,:) )*vmask(:,:,jk) 139 136 END DO 140 137 ENDIF … … 147 144 # endif 148 145 ! 149 CALL lbc_lnk_multi( 'dynnxt', u u(:,:,:,kNp1), 'U', -1., vv(:,:,:,kNp1), 'V', -1. ) !* local domain boundaries146 CALL lbc_lnk_multi( 'dynnxt', ua, 'U', -1., va, 'V', -1. ) !* local domain boundaries 150 147 ! 151 148 ! !* BDY open boundaries … … 160 157 ! 161 158 ! ! Kinetic energy and Conversion 162 IF( ln_KE_trd ) CALL trd_dyn( u u(:,:,:,kNp1), vv(:,:,:,kNp1), jpdyn_ken, kt )159 IF( ln_KE_trd ) CALL trd_dyn( ua, va, jpdyn_ken, kt ) 163 160 ! 164 161 IF( ln_dyn_trd ) THEN ! 3D output: total momentum trends 165 zua(:,:,:) = ( u u(:,:,:,kNp1) - uu(:,:,:,kNm1) ) * z1_2dt166 zva(:,:,:) = ( v v(:,:,:,kNp1) - vv(:,:,:,kNm1) ) * z1_2dt162 zua(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * z1_2dt 163 zva(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * z1_2dt 167 164 CALL iom_put( "utrd_tot", zua ) ! total momentum trends, except the asselin time filter 168 165 CALL iom_put( "vtrd_tot", zva ) 169 166 ENDIF 170 167 ! 171 zua(:,:,:) = puu(:,:,:) ! save the now velocity before the asselin filter172 zva(:,:,:) = pvv(:,:,:) ! (caution: there will be a shift by 1 timestep in the168 zua(:,:,:) = un(:,:,:) ! save the now velocity before the asselin filter 169 zva(:,:,:) = vn(:,:,:) ! (caution: there will be a shift by 1 timestep in the 173 170 ! ! computation of the asselin filter trends) 174 171 ENDIF … … 176 173 ! Time filter and swap of dynamics arrays 177 174 ! ------------------------------------------ 175 IF( neuler == 0 .AND. kt == nit000 ) THEN !* Euler at first time-step: only swap 176 DO jk = 1, jpkm1 177 un(:,:,jk) = ua(:,:,jk) ! un <-- ua 178 vn(:,:,jk) = va(:,:,jk) 179 END DO 180 IF( .NOT.ln_linssh ) THEN ! e3._b <-- e3._n 181 !!gm BUG ???? I don't understand why it is not : e3._n <-- e3._a 182 DO jk = 1, jpkm1 183 ! e3t_b(:,:,jk) = e3t_n(:,:,jk) 184 ! e3u_b(:,:,jk) = e3u_n(:,:,jk) 185 ! e3v_b(:,:,jk) = e3v_n(:,:,jk) 186 ! 187 e3t_n(:,:,jk) = e3t_a(:,:,jk) 188 e3u_n(:,:,jk) = e3u_a(:,:,jk) 189 e3v_n(:,:,jk) = e3v_a(:,:,jk) 190 END DO 191 !!gm BUG end 192 ENDIF 193 ! 178 194 179 IF( .NOT.( neuler == 0 .AND. kt == nit000 ) ) THEN!* Leap-Frog : Asselin filter and swap195 ELSE !* Leap-Frog : Asselin filter and swap 180 196 ! ! =============! 181 197 IF( ln_linssh ) THEN ! Fixed volume ! … … 184 200 DO jj = 1, jpj 185 201 DO ji = 1, jpi 186 puu(ji,jj,jk) = puu(ji,jj,jk) + atfp * ( uu(ji,jj,jk,kNm1) - 2._wp * puu(ji,jj,jk) + uu(ji,jj,jk,kNp1) ) 187 pvv(ji,jj,jk) = pvv(ji,jj,jk) + atfp * ( vv(ji,jj,jk,kNm1) - 2._wp * pvv(ji,jj,jk) + vv(ji,jj,jk,kNp1) ) 202 zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) ) 203 zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) ) 204 ! 205 ub(ji,jj,jk) = zuf ! ub <-- filtered velocity 206 vb(ji,jj,jk) = zvf 207 un(ji,jj,jk) = ua(ji,jj,jk) ! un <-- ua 208 vn(ji,jj,jk) = va(ji,jj,jk) 188 209 END DO 189 210 END DO … … 192 213 ELSE ! Variable volume ! 193 214 ! ! ================! 194 ! Time-filtered scale factor at t-points 215 ! Before scale factor at t-points 216 ! (used as a now filtered scale factor until the swap) 195 217 ! ---------------------------------------------------- 196 ALLOCATE( ze3t_f(jpi,jpj,jpk) )197 218 DO jk = 1, jpkm1 198 ze3t_f(:,:,jk) = pe3t(:,:,jk) + atfp * ( e3t(:,:,jk,kNm1) - 2._wp * pe3t(:,:,jk) + e3t(:,:,jk,kNp1) )219 e3t_b(:,:,jk) = e3t_n(:,:,jk) + atfp * ( e3t_b(:,:,jk) - 2._wp * e3t_n(:,:,jk) + e3t_a(:,:,jk) ) 199 220 END DO 200 221 ! Add volume filter correction: compatibility with tracer advection scheme … … 202 223 zcoef = atfp * rdt * r1_rau0 203 224 204 ze3t_f(:,:,1) = ze3t_f(:,:,1) - zcoef * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)225 e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1) 205 226 206 227 IF ( ln_rnf ) THEN … … 210 231 DO ji = 1, jpi 211 232 IF( jk <= nk_rnf(ji,jj) ) THEN 212 ze3t_f(ji,jj,jk) = ze3t_f(ji,jj,jk) - zcoef * ( - rnf_b(ji,jj) + rnf(ji,jj) ) &213 & * ( pe3t(ji,jj,jk) / h_rnf(ji,jj) ) * tmask(ji,jj,jk)233 e3t_b(ji,jj,jk) = e3t_b(ji,jj,jk) - zcoef * ( - rnf_b(ji,jj) + rnf(ji,jj) ) & 234 & * ( e3t_n(ji,jj,jk) / h_rnf(ji,jj) ) * tmask(ji,jj,jk) 214 235 ENDIF 215 236 ENDDO … … 217 238 ENDDO 218 239 ELSE 219 ze3t_f(:,:,1) = ze3t_f(:,:,1) - zcoef * ( -rnf_b(:,:) + rnf(:,:))*tmask(:,:,1)240 e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef * ( -rnf_b(:,:) + rnf(:,:))*tmask(:,:,1) 220 241 ENDIF 221 242 END IF … … 226 247 DO ji = 1, jpi 227 248 IF( misfkt(ji,jj) <=jk .and. jk < misfkb(ji,jj) ) THEN 228 ze3t_f(ji,jj,jk) = ze3t_f(ji,jj,jk) - zcoef * ( fwfisf_b(ji,jj) - fwfisf(ji,jj) ) &229 & * ( pe3t(ji,jj,jk) * r1_hisf_tbl(ji,jj) ) * tmask(ji,jj,jk)249 e3t_b(ji,jj,jk) = e3t_b(ji,jj,jk) - zcoef * ( fwfisf_b(ji,jj) - fwfisf(ji,jj) ) & 250 & * ( e3t_n(ji,jj,jk) * r1_hisf_tbl(ji,jj) ) * tmask(ji,jj,jk) 230 251 ELSEIF ( jk==misfkb(ji,jj) ) THEN 231 ze3t_f(ji,jj,jk) = ze3t_f(ji,jj,jk) - zcoef * ( fwfisf_b(ji,jj) - fwfisf(ji,jj) ) &232 & * ( pe3t(ji,jj,jk) * r1_hisf_tbl(ji,jj) ) * ralpha(ji,jj) * tmask(ji,jj,jk)252 e3t_b(ji,jj,jk) = e3t_b(ji,jj,jk) - zcoef * ( fwfisf_b(ji,jj) - fwfisf(ji,jj) ) & 253 & * ( e3t_n(ji,jj,jk) * r1_hisf_tbl(ji,jj) ) * ralpha(ji,jj) * tmask(ji,jj,jk) 233 254 ENDIF 234 255 END DO … … 237 258 END IF 238 259 ! 239 pe3t(:,:,1:jpkm1) = ze3t_f(:,:,1:jpkm1) ! filtered scale factor at T-points240 !241 260 IF( ln_dynadv_vec ) THEN ! Asselin filter applied on velocity 242 261 ! Before filtered scale factor at (u/v)-points 243 CALL dom_vvl_interpol( pe3t(:,:,:), pe3u(:,:,:), 'U' )244 CALL dom_vvl_interpol( pe3t(:,:,:), pe3v(:,:,:), 'V' )262 CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' ) 263 CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' ) 245 264 DO jk = 1, jpkm1 246 265 DO jj = 1, jpj 247 266 DO ji = 1, jpi 248 puu(ji,jj,jk) = puu(ji,jj,jk) + atfp * ( uu(ji,jj,jk,kNm1) - 2._wp * puu(ji,jj,jk) + uu(ji,jj,jk,kNp1) ) 249 pvv(ji,jj,jk) = pvv(ji,jj,jk) + atfp * ( vv(ji,jj,jk,kNm1) - 2._wp * pvv(ji,jj,jk) + vv(ji,jj,jk,kNp1) ) 267 zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) ) 268 zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) ) 269 ! 270 ub(ji,jj,jk) = zuf ! ub <-- filtered velocity 271 vb(ji,jj,jk) = zvf 272 un(ji,jj,jk) = ua(ji,jj,jk) ! un <-- ua 273 vn(ji,jj,jk) = va(ji,jj,jk) 250 274 END DO 251 275 END DO … … 255 279 ! 256 280 ALLOCATE( ze3u_f(jpi,jpj,jpk) , ze3v_f(jpi,jpj,jpk) ) 257 ! Nowfiltered scale factor at (u/v)-points stored in ze3u_f, ze3v_f258 CALL dom_vvl_interpol( pe3t(:,:,:), ze3u_f, 'U' )259 CALL dom_vvl_interpol( pe3t(:,:,:), ze3v_f, 'V' )281 ! Before filtered scale factor at (u/v)-points stored in ze3u_f, ze3v_f 282 CALL dom_vvl_interpol( e3t_b(:,:,:), ze3u_f, 'U' ) 283 CALL dom_vvl_interpol( e3t_b(:,:,:), ze3v_f, 'V' ) 260 284 DO jk = 1, jpkm1 261 285 DO jj = 1, jpj 262 286 DO ji = 1, jpi 263 zue3a = e3u (ji,jj,jk,kNp1) * uu(ji,jj,jk,kNp1)264 zve3a = e3v (ji,jj,jk,kNp1) * vv(ji,jj,jk,kNp1)265 zue3n = pe3u(ji,jj,jk) * puu(ji,jj,jk)266 zve3n = pe3v(ji,jj,jk) * pvv(ji,jj,jk)267 zue3b = e3u (ji,jj,jk,kNm1) * uu(ji,jj,jk,kNm1)268 zve3b = e3v (ji,jj,jk,kNm1) * vv(ji,jj,jk,kNm1)287 zue3a = e3u_a(ji,jj,jk) * ua(ji,jj,jk) 288 zve3a = e3v_a(ji,jj,jk) * va(ji,jj,jk) 289 zue3n = e3u_n(ji,jj,jk) * un(ji,jj,jk) 290 zve3n = e3v_n(ji,jj,jk) * vn(ji,jj,jk) 291 zue3b = e3u_b(ji,jj,jk) * ub(ji,jj,jk) 292 zve3b = e3v_b(ji,jj,jk) * vb(ji,jj,jk) 269 293 ! 270 puu(ji,jj,jk) = ( zue3n + atfp * ( zue3b - 2._wp * zue3n + zue3a ) ) / ze3u_f(ji,jj,jk) 271 pvv(ji,jj,jk) = ( zve3n + atfp * ( zve3b - 2._wp * zve3n + zve3a ) ) / ze3v_f(ji,jj,jk) 294 zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n + zue3a ) ) / ze3u_f(ji,jj,jk) 295 zvf = ( zve3n + atfp * ( zve3b - 2._wp * zve3n + zve3a ) ) / ze3v_f(ji,jj,jk) 296 ! 297 ub(ji,jj,jk) = zuf ! ub <-- filtered velocity 298 vb(ji,jj,jk) = zvf 299 un(ji,jj,jk) = ua(ji,jj,jk) ! un <-- ua 300 vn(ji,jj,jk) = va(ji,jj,jk) 272 301 END DO 273 302 END DO 274 303 END DO 275 pe3u(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1)276 pe3v(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1)304 e3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1) ! e3u_b <-- filtered scale factor 305 e3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1) 277 306 ! 278 307 DEALLOCATE( ze3u_f , ze3v_f ) 279 308 ENDIF 280 309 ! 281 DEALLOCATE( ze3t_f )282 310 ENDIF 283 311 ! 284 312 IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN 285 ! Revert filtered "now" velocities to time split estimate313 ! Revert "before" velocities to time split estimate 286 314 ! Doing it here also means that asselin filter contribution is removed 287 zue(:,:) = pe3u(:,:,1) * puu(:,:,1) * umask(:,:,1)288 zve(:,:) = pe3v(:,:,1) * pvv(:,:,1) * vmask(:,:,1)315 zue(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1) 316 zve(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1) 289 317 DO jk = 2, jpkm1 290 zue(:,:) = zue(:,:) + pe3u(:,:,jk) * puu(:,:,jk) * umask(:,:,jk)291 zve(:,:) = zve(:,:) + pe3v(:,:,jk) * pvv(:,:,jk) * vmask(:,:,jk)318 zue(:,:) = zue(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk) 319 zve(:,:) = zve(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk) 292 320 END DO 293 321 DO jk = 1, jpkm1 294 puu(:,:,jk) = puu(:,:,jk) - (zue(:,:) * r1_hu_n(:,:) - un_b(:,:)) * umask(:,:,jk) 295 pvv(:,:,jk) = pvv(:,:,jk) - (zve(:,:) * r1_hv_n(:,:) - vn_b(:,:)) * vmask(:,:,jk) 296 END DO 297 ENDIF 298 ! 299 ikt = Nnn 300 ELSE 301 ikt = kNm1 302 puu(:,:,:) = uu(:,:,:,kNp1) 303 pvv(:,:,:) = vv(:,:,:,kNp1) 304 pe3t(:,:,:) = e3t(:,:,:,kNp1) 305 pe3u(:,:,:) = e3u(:,:,:,kNp1) 306 pe3v(:,:,:) = e3v(:,:,:,kNp1) 322 ub(:,:,jk) = ub(:,:,jk) - (zue(:,:) * r1_hu_n(:,:) - un_b(:,:)) * umask(:,:,jk) 323 vb(:,:,jk) = vb(:,:,jk) - (zve(:,:) * r1_hv_n(:,:) - vn_b(:,:)) * vmask(:,:,jk) 324 END DO 325 ENDIF 326 ! 307 327 ENDIF ! neuler =/0 308 328 ! … … 311 331 ! integration 312 332 ! 313 ! DS IMMERSE :: This is very confusing as it stands. But should314 ! hopefully be simpler when we do the time-level swapping for the315 ! external mode solution.316 333 ! 317 334 IF(.NOT.ln_linssh ) THEN 318 hu (:,:,ikt) = e3u(:,:,1,ikt) * umask(:,:,1)319 hv (:,:,ikt) = e3v(:,:,1,ikt) * vmask(:,:,1)335 hu_b(:,:) = e3u_b(:,:,1) * umask(:,:,1) 336 hv_b(:,:) = e3v_b(:,:,1) * vmask(:,:,1) 320 337 DO jk = 2, jpkm1 321 hu (:,:,ikt) = hu(:,:,ikt) + e3u(:,:,jk,ikt) * umask(:,:,jk)322 hv (:,:,ikt) = hv(:,:,ikt) + e3v(:,:,jk,ikt) * vmask(:,:,jk)338 hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk) 339 hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk) 323 340 END DO 324 r1_hu (:,:,ikt) = ssumask(:,:) / ( hu(:,:,ikt) + 1._wp - ssumask(:,:) )325 r1_hv (:,:,ikt) = ssvmask(:,:) / ( hv(:,:,ikt) + 1._wp - ssvmask(:,:) )326 ENDIF 327 ! 328 un_b(:,:) = e3u (:,:,1,kNp1) * uu(:,:,1,kNp1) * umask(:,:,1)329 ub_b(:,:) = e3u (:,:,1,ikt ) * uu(:,:,1,ikt) * umask(:,:,1)330 vn_b(:,:) = e3v (:,:,1,kNp1) * vv(:,:,1,kNp1) * vmask(:,:,1)331 vb_b(:,:) = e3v (:,:,1,ikt ) * vv(:,:,1,ikt) * vmask(:,:,1)341 r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) ) 342 r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) ) 343 ENDIF 344 ! 345 un_b(:,:) = e3u_a(:,:,1) * un(:,:,1) * umask(:,:,1) 346 ub_b(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1) 347 vn_b(:,:) = e3v_a(:,:,1) * vn(:,:,1) * vmask(:,:,1) 348 vb_b(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1) 332 349 DO jk = 2, jpkm1 333 un_b(:,:) = un_b(:,:) + e3u (:,:,jk,kNp1) * uu(:,:,jk,kNp1) * umask(:,:,jk)334 ub_b(:,:) = ub_b(:,:) + e3u (:,:,jk,ikt ) * uu(:,:,jk,ikt) * umask(:,:,jk)335 vn_b(:,:) = vn_b(:,:) + e3v (:,:,jk,kNp1) * vv(:,:,jk,kNp1) * vmask(:,:,jk)336 vb_b(:,:) = vb_b(:,:) + e3v (:,:,jk,ikt ) * vv(:,:,jk,ikt) * vmask(:,:,jk)350 un_b(:,:) = un_b(:,:) + e3u_a(:,:,jk) * un(:,:,jk) * umask(:,:,jk) 351 ub_b(:,:) = ub_b(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk) 352 vn_b(:,:) = vn_b(:,:) + e3v_a(:,:,jk) * vn(:,:,jk) * vmask(:,:,jk) 353 vb_b(:,:) = vb_b(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk) 337 354 END DO 338 un_b(:,:) = un_b(:,:) * r1_hu (:,:,kNp1)339 vn_b(:,:) = vn_b(:,:) * r1_hv (:,:,kNp1)340 ub_b(:,:) = ub_b(:,:) * r1_hu (:,:,ikt)341 vb_b(:,:) = vb_b(:,:) * r1_hv (:,:,ikt)355 un_b(:,:) = un_b(:,:) * r1_hu_a(:,:) 356 vn_b(:,:) = vn_b(:,:) * r1_hv_a(:,:) 357 ub_b(:,:) = ub_b(:,:) * r1_hu_b(:,:) 358 vb_b(:,:) = vb_b(:,:) * r1_hv_b(:,:) 342 359 ! 343 360 IF( .NOT.ln_dynspg_ts ) THEN ! output the barotropic currents … … 346 363 ENDIF 347 364 IF( l_trddyn ) THEN ! 3D output: asselin filter trends on momentum 348 zua(:,:,:) = ( puu(:,:,:) - zua(:,:,:) ) * z1_2dt349 zva(:,:,:) = ( pvv(:,:,:) - zva(:,:,:) ) * z1_2dt365 zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * z1_2dt 366 zva(:,:,:) = ( vb(:,:,:) - zva(:,:,:) ) * z1_2dt 350 367 CALL trd_dyn( zua, zva, jpdyn_atf, kt ) 351 368 ENDIF 352 369 ! 353 IF(ln_ctl) CALL prt_ctl( tab3d_1=u u(:,:,:,kNp1), clinfo1=' nxt - uu(:,:,:,kNp1): ', mask1=umask, &354 & tab3d_2=v v(:,:,:,kNp1), clinfo2=' vv(:,:,:,kNp1): ' , mask2=vmask )370 IF(ln_ctl) CALL prt_ctl( tab3d_1=un, clinfo1=' nxt - Un: ', mask1=umask, & 371 & tab3d_2=vn, clinfo2=' Vn: ' , mask2=vmask ) 355 372 ! 356 373 IF( ln_dynspg_ts ) DEALLOCATE( zue, zve ) -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/step.F90
r10789 r10799 277 277 !!jc2: dynnxt must be the latest call. e3t_b are indeed updated in that routine 278 278 CALL tra_nxt ( kstp ) ! finalize (bcs) tracer fields at next time step and swap 279 CALL dyn_nxt ( kstp , Nm1, Np1, uu(:,:,:,Nnn), vv(:,:,:,Nnn), e3t(:,:,:,Nnn), e3u(:,:,:,Nnn), e3v(:,:,:,Nnn)) !279 CALL dyn_nxt ( kstp ) ! 280 280 CALL ssh_swp ( kstp ) ! swap of sea surface height 281 ! 282 ! Swap time levels 283 IF( .NOT. (neuler == 0 .AND. kstp == nit000) ) THEN 284 Nrhs = Nm1 285 Nm1 = Nnn 286 Nnn = Np1 287 Np1 = Nrhs 288 ENDIF 289 ! 290 ! Update temporary pointers 291 CALL update_pointers() 292 ! 293 ! Note that 2-time-level indices don't need to be swapped because both "before" and "now" fields are derived in dom_vvl_sf_swp 294 IF(.NOT.ln_linssh) CALL dom_vvl_sf_swp( kstp, Nm1, Nnn, Nm1_2lev, Nnn_2lev ) ! interpolate vertical scale factors for Nnn time level 281 IF(.NOT.ln_linssh) CALL dom_vvl_sf_swp( kstp ) ! interpolate vertical scale factors for Nnn time level 295 282 ! 296 283 IF( ln_diahsb ) CALL dia_hsb ( kstp ) ! - ML - global conservation diagnostics
Note: See TracChangeset
for help on using the changeset viewer.
