- Timestamp:
- 2013-11-20T17:28:04+01:00 (10 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
r3764 r4292 8 8 !! - ! 2010-05 (K. Mogensen, A. Weaver, M. Martin, D. Lea) Assimilation interface 9 9 !! - ! 2010-09 (D.Storkey and E.O'Dea) bug fixes for BDY module 10 !!---------------------------------------------------------------------- 11 12 !!---------------------------------------------------------------------- 13 !! ssh_wzv : after ssh & now vertical velocity 14 !! ssh_nxt : filter ans swap the ssh arrays 10 !! 3.3 ! 2011-10 (M. Leclair) split former ssh_wzv routine and remove all vvl related work 11 !!---------------------------------------------------------------------- 12 13 !!---------------------------------------------------------------------- 14 !! ssh_nxt : after ssh 15 !! ssh_swp : filter ans swap the ssh arrays 16 !! wzv : compute now vertical velocity 15 17 !!---------------------------------------------------------------------- 16 18 USE oce ! ocean dynamics and tracers variables … … 20 22 USE divcur ! hor. divergence and curl (div & cur routines) 21 23 USE iom ! I/O library 24 USE restart ! only for lrst_oce 22 25 USE in_out_manager ! I/O manager 23 26 USE prtctl ! Print control … … 28 31 USE obc_oce 29 32 USE bdy_oce 33 USE bdy_par 34 USE bdydyn2d ! bdy_ssh routine 30 35 USE diaar5, ONLY: lk_diaar5 31 36 USE iom 32 USE sbcrnf, ONLY: h_rnf, nk_rnf ! River runoff 37 USE sbcrnf, ONLY: h_rnf, nk_rnf, sbc_rnf_div ! River runoff 38 USE dynspg_ts, ONLY: ln_bt_fw 39 USE dynspg_oce, ONLY: lk_dynspg_ts 33 40 #if defined key_agrif 34 41 USE agrif_opa_update … … 44 51 PRIVATE 45 52 46 PUBLIC ssh_wzv ! called by step.F9047 53 PUBLIC ssh_nxt ! called by step.F90 54 PUBLIC wzv ! called by step.F90 55 PUBLIC ssh_swp ! called by step.F90 48 56 49 57 !! * Substitutions … … 57 65 CONTAINS 58 66 59 SUBROUTINE ssh_ wzv( kt )60 !!---------------------------------------------------------------------- 61 !! *** ROUTINE ssh_ wzv***67 SUBROUTINE ssh_nxt( kt ) 68 !!---------------------------------------------------------------------- 69 !! *** ROUTINE ssh_nxt *** 62 70 !! 63 !! ** Purpose : compute the after ssh (ssha), the now vertical velocity 64 !! and update the now vertical coordinate (lk_vvl=T). 65 !! 66 !! ** Method : - Using the incompressibility hypothesis, the vertical 67 !! velocity is computed by integrating the horizontal divergence 68 !! from the bottom to the surface minus the scale factor evolution. 69 !! The boundary conditions are w=0 at the bottom (no flux) and. 71 !! ** Purpose : compute the after ssh (ssha) 72 !! 73 !! ** Method : - Using the incompressibility hypothesis, the ssh increment 74 !! is computed by integrating the horizontal divergence and multiply by 75 !! by the time step. 70 76 !! 71 77 !! ** action : ssha : after sea surface height 72 !! wn : now vertical velocity73 !! sshu_a, sshv_a, sshf_a : after sea surface height (lk_vvl=T)74 !! hu, hv, hur, hvr : ocean depth and its inverse at u-,v-points75 78 !! 76 79 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 77 80 !!---------------------------------------------------------------------- 78 INTEGER, INTENT(in) :: kt ! time step79 !80 INTEGER :: ji, jj, jk ! dummy loop indices81 REAL(wp) :: zcoefu, zcoefv, zcoeff, z2dt, z1_2dt, z1_rau0 ! local scalars82 REAL(wp), POINTER, DIMENSION(:,: ) :: z2d, zhdiv83 REAL(wp) , POINTER, DIMENSION(:,:,:) :: z3d84 !!---------------------------------------------------------------------- 85 ! 86 IF( nn_timing == 1 ) CALL timing_start('ssh_ wzv')87 ! 88 CALL wrk_alloc( jpi, jpj, z 2d, zhdiv )81 ! 82 REAL(wp), POINTER, DIMENSION(:,: ) :: zhdiv 83 INTEGER, INTENT(in) :: kt ! time step 84 ! 85 INTEGER :: jk ! dummy loop indice 86 REAL(wp) :: z2dt, z1_rau0 ! local scalars 87 !!---------------------------------------------------------------------- 88 ! 89 IF( nn_timing == 1 ) CALL timing_start('ssh_nxt') 90 ! 91 CALL wrk_alloc( jpi, jpj, zhdiv ) 89 92 ! 90 93 IF( kt == nit000 ) THEN 91 94 ! 92 95 IF(lwp) WRITE(numout,*) 93 IF(lwp) WRITE(numout,*) 'ssh_ wzv : after sea surface height and now vertical velocity'96 IF(lwp) WRITE(numout,*) 'ssh_nxt : after sea surface height' 94 97 IF(lwp) WRITE(numout,*) '~~~~~~~ ' 95 98 ! 96 wn(:,:,jpk) = 0._wp ! bottom boundary condition: w=0 (set once for all)97 !98 IF( lk_vvl ) THEN ! before and now Sea SSH at u-, v-, f-points (vvl case only)99 DO jj = 1, jpjm1100 DO ji = 1, jpim1 ! caution: use of Vector Opt. not possible101 zcoefu = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )102 zcoefv = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )103 zcoeff = 0.25 * umask(ji,jj,1) * umask(ji,jj+1,1)104 sshu_b(ji,jj) = zcoefu * ( e1t(ji ,jj) * e2t(ji ,jj) * sshb(ji ,jj) &105 & + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )106 sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj ) * e2t(ji,jj ) * sshb(ji,jj ) &107 & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )108 sshu_n(ji,jj) = zcoefu * ( e1t(ji ,jj) * e2t(ji ,jj) * sshn(ji ,jj) &109 & + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) )110 sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj ) * e2t(ji,jj ) * sshn(ji,jj ) &111 & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) )112 END DO113 END DO114 CALL lbc_lnk( sshu_b, 'U', 1. ) ; CALL lbc_lnk( sshu_n, 'U', 1. )115 CALL lbc_lnk( sshv_b, 'V', 1. ) ; CALL lbc_lnk( sshv_n, 'V', 1. )116 DO jj = 1, jpjm1117 DO ji = 1, jpim1 ! NO Vector Opt.118 sshf_n(ji,jj) = 0.5 * umask(ji,jj,1) * umask(ji,jj+1,1) &119 & / ( e1f(ji,jj ) * e2f(ji,jj ) ) &120 & * ( e1u(ji,jj ) * e2u(ji,jj ) * sshu_n(ji,jj ) &121 & + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )122 END DO123 END DO124 CALL lbc_lnk( sshf_n, 'F', 1. )125 ENDIF126 !127 ENDIF128 129 ! !------------------------------------------!130 IF( lk_vvl ) THEN ! Regridding: Update Now Vertical coord. ! (only in vvl case)131 ! !------------------------------------------!132 DO jk = 1, jpkm1133 fsdept(:,:,jk) = fsdept_n(:,:,jk) ! now local depths stored in fsdep. arrays134 fsdepw(:,:,jk) = fsdepw_n(:,:,jk)135 fsde3w(:,:,jk) = fsde3w_n(:,:,jk)136 !137 fse3t (:,:,jk) = fse3t_n (:,:,jk) ! vertical scale factors stored in fse3. arrays138 fse3u (:,:,jk) = fse3u_n (:,:,jk)139 fse3v (:,:,jk) = fse3v_n (:,:,jk)140 fse3f (:,:,jk) = fse3f_n (:,:,jk)141 fse3w (:,:,jk) = fse3w_n (:,:,jk)142 fse3uw(:,:,jk) = fse3uw_n(:,:,jk)143 fse3vw(:,:,jk) = fse3vw_n(:,:,jk)144 END DO145 !146 hu(:,:) = hu_0(:,:) + sshu_n(:,:) ! now ocean depth (at u- and v-points)147 hv(:,:) = hv_0(:,:) + sshv_n(:,:)148 ! ! now masked inverse of the ocean depth (at u- and v-points)149 hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1._wp - umask(:,:,1) )150 hvr(:,:) = vmask(:,:,1) / ( hv(:,:) + 1._wp - vmask(:,:,1) )151 !152 99 ENDIF 153 100 ! … … 162 109 zhdiv(:,:) = 0._wp 163 110 DO jk = 1, jpkm1 ! Horizontal divergence of barotropic transports 164 zhdiv(:,:) = zhdiv(:,:) + fse3t (:,:,jk) * hdivn(:,:,jk)111 zhdiv(:,:) = zhdiv(:,:) + fse3t_n(:,:,jk) * hdivn(:,:,jk) 165 112 END DO 166 113 ! ! Sea surface elevation time stepping 167 114 ! In forward Euler time stepping case, the same formulation as in the leap-frog case can be used 168 115 ! because emp_b field is initialized with the vlaues of emp field. Hence, 0.5 * ( emp + emp_b ) = emp 169 z1_rau0 = 0.5 /rau0116 z1_rau0 = 0.5_wp * r1_rau0 170 117 ssha(:,:) = ( sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * tmask(:,:,1) 171 118 … … 180 127 #endif 181 128 #if defined key_bdy 182 ssha(:,:) = ssha(:,:) * bdytmask(:,:) 183 CALL lbc_lnk( ssha, 'T', 1. ) ! absolutly compulsory !! (jmm) 184 #endif 129 ! bg jchanut tschanges 130 ! These lines are not necessary with time splitting since 131 ! boundary condition on sea level is set during ts loop 132 IF (lk_bdy) THEN 133 CALL lbc_lnk( ssha, 'T', 1. ) ! Not sure that's necessary 134 CALL bdy_ssh( ssha ) ! Duplicate sea level across open boundaries 135 ENDIF 136 #endif 137 ! end jchanut tschanges 185 138 #if defined key_asminc 186 139 ! ! Include the IAU weighted SSH increment … … 190 143 ENDIF 191 144 #endif 192 ! ! Sea Surface Height at u-,v- and f-points (vvl case only) 193 IF( lk_vvl ) THEN ! (required only in key_vvl case) 194 DO jj = 1, jpjm1 195 DO ji = 1, jpim1 ! NO Vector Opt. 196 sshu_a(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji ,jj) * e2u(ji ,jj) ) & 197 & * ( e1t(ji ,jj) * e2t(ji ,jj) * ssha(ji ,jj) & 198 & + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) ) 199 sshv_a(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj ) * e2v(ji,jj ) ) & 200 & * ( e1t(ji,jj ) * e2t(ji,jj ) * ssha(ji,jj ) & 201 & + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) ) 145 146 ! !------------------------------! 147 ! ! outputs ! 148 ! !------------------------------! 149 CALL iom_put( "ssh" , sshn ) ! sea surface height 150 CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) ) ! square of sea surface height 151 ! 152 IF(ln_ctl) CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha - : ', mask1=tmask, ovlap=1 ) 153 ! 154 CALL wrk_dealloc( jpi, jpj, zhdiv ) 155 ! 156 IF( nn_timing == 1 ) CALL timing_stop('ssh_nxt') 157 ! 158 END SUBROUTINE ssh_nxt 159 160 161 SUBROUTINE wzv( kt ) 162 !!---------------------------------------------------------------------- 163 !! *** ROUTINE wzv *** 164 !! 165 !! ** Purpose : compute the now vertical velocity 166 !! 167 !! ** Method : - Using the incompressibility hypothesis, the vertical 168 !! velocity is computed by integrating the horizontal divergence 169 !! from the bottom to the surface minus the scale factor evolution. 170 !! The boundary conditions are w=0 at the bottom (no flux) and. 171 !! 172 !! ** action : wn : now vertical velocity 173 !! 174 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 175 !!---------------------------------------------------------------------- 176 ! 177 INTEGER, INTENT(in) :: kt ! time step 178 REAL(wp), POINTER, DIMENSION(:,: ) :: z2d 179 REAL(wp), POINTER, DIMENSION(:,:,:) :: z3d, zhdiv 180 ! 181 INTEGER :: ji, jj, jk ! dummy loop indices 182 REAL(wp) :: z1_2dt ! local scalars 183 !!---------------------------------------------------------------------- 184 185 IF( nn_timing == 1 ) CALL timing_start('wzv') 186 ! 187 IF( kt == nit000 ) THEN 188 ! 189 IF(lwp) WRITE(numout,*) 190 IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity ' 191 IF(lwp) WRITE(numout,*) '~~~~~ ' 192 ! 193 wn(:,:,jpk) = 0._wp ! bottom boundary condition: w=0 (set once for all) 194 ! 195 ENDIF 196 ! !------------------------------! 197 ! ! Now Vertical Velocity ! 198 ! !------------------------------! 199 z1_2dt = 1. / ( 2. * rdt ) ! set time step size (Euler/Leapfrog) 200 IF( neuler == 0 .AND. kt == nit000 ) z1_2dt = 1. / rdt 201 ! 202 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde and layer cases 203 CALL wrk_alloc( jpi, jpj, jpk, zhdiv ) 204 ! 205 DO jk = 1, jpkm1 206 ! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t) 207 ! - ML - note: computation allready done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way) 208 DO jj = 2, jpjm1 209 DO ji = fs_2, fs_jpim1 ! vector opt. 210 zhdiv(ji,jj,jk) = r1_e12t(ji,jj) * ( un_td(ji,jj,jk) - un_td(ji-1,jj,jk) + vn_td(ji,jj,jk) - vn_td(ji,jj-1,jk) ) 211 END DO 202 212 END DO 203 213 END DO 204 CALL lbc_lnk( sshu_a, 'U', 1. ) ; CALL lbc_lnk( sshv_a, 'V', 1. ) ! Boundaries conditions 205 ENDIF 206 207 ! !------------------------------! 208 ! ! Now Vertical Velocity ! 209 ! !------------------------------! 210 z1_2dt = 1.e0 / z2dt 211 DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence 212 ! - ML - need 3 lines here because replacement of fse3t by its expression yields too long lines otherwise 213 wn(:,:,jk) = wn(:,:,jk+1) - fse3t_n(:,:,jk) * hdivn(:,:,jk) & 214 & - ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) & 215 & * tmask(:,:,jk) * z1_2dt 214 CALL lbc_lnk(zhdiv, 'T', 1.) ! - ML - Perhaps not necessary: not used for horizontal "connexions" 215 ! ! Is it problematic to have a wrong vertical velocity in boundary cells? 216 ! ! Same question holds for hdivn. Perhaps just for security 217 DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence 218 ! computation of w 219 wn(:,:,jk) = wn(:,:,jk+1) - ( fse3t_n(:,:,jk) * hdivn(:,:,jk) + zhdiv(:,:,jk) & 220 & + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) ) * tmask(:,:,jk) 221 END DO 222 ! IF( ln_vvl_layer ) wn(:,:,:) = 0.e0 223 CALL wrk_dealloc( jpi, jpj, jpk, zhdiv ) 224 ELSE ! z_star and linear free surface cases 225 DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence 226 ! computation of w 227 wn(:,:,jk) = wn(:,:,jk+1) - ( fse3t_n(:,:,jk) * hdivn(:,:,jk) & 228 & + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) ) * tmask(:,:,jk) 229 END DO 230 ENDIF 231 216 232 #if defined key_bdy 217 233 wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:) 218 234 #endif 219 END DO 220 235 ! 221 236 ! !------------------------------! 222 237 ! ! outputs ! 223 238 ! !------------------------------! 224 CALL iom_put( "woce", wn ) ! vertical velocity 225 CALL iom_put( "ssh" , sshn ) ! sea surface height 226 CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) ) ! square of sea surface height 239 CALL iom_put( "woce", wn ) ! vertical velocity 227 240 IF( lk_diaar5 ) THEN ! vertical mass transport & its square value 241 CALL wrk_alloc( jpi, jpj, z2d ) 242 CALL wrk_alloc( jpi, jpj, jpk, z3d ) 228 243 ! Caution: in the VVL case, it only correponds to the baroclinic mass transport. 229 CALL wrk_alloc( jpi,jpj,jpk, z3d ) 230 z2d(:,:) = rau0 * e1t(:,:) * e2t(:,:) 244 z2d(:,:) = rau0 * e12t(:,:) 231 245 DO jk = 1, jpk 232 246 z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:) … … 234 248 CALL iom_put( "w_masstr" , z3d ) 235 249 CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) ) 236 CALL wrk_dealloc( jpi,jpj,jpk, z3d ) 237 ENDIF 238 ! 239 IF(ln_ctl) CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha - : ', mask1=tmask, ovlap=1 ) 240 ! 241 CALL wrk_dealloc( jpi, jpj, z2d, zhdiv ) 242 ! 243 IF( nn_timing == 1 ) CALL timing_stop('ssh_wzv') 244 ! 245 END SUBROUTINE ssh_wzv 246 247 248 SUBROUTINE ssh_nxt( kt ) 250 CALL wrk_dealloc( jpi, jpj, z2d ) 251 CALL wrk_dealloc( jpi, jpj, jpk, z3d ) 252 ENDIF 253 ! 254 IF( nn_timing == 1 ) CALL timing_stop('wzv') 255 256 257 END SUBROUTINE wzv 258 259 SUBROUTINE ssh_swp( kt ) 249 260 !!---------------------------------------------------------------------- 250 261 !! *** ROUTINE ssh_nxt *** … … 252 263 !! ** Purpose : achieve the sea surface height time stepping by 253 264 !! applying Asselin time filter and swapping the arrays 254 !! ssha already computed in ssh_ wzv265 !! ssha already computed in ssh_nxt 255 266 !! 256 267 !! ** Method : - apply Asselin time fiter to now ssh (excluding the forcing … … 266 277 !!---------------------------------------------------------------------- 267 278 INTEGER, INTENT(in) :: kt ! ocean time-step index 268 !! 269 INTEGER :: ji, jj ! dummy loop indices 270 REAL(wp) :: zec ! temporary scalar 271 !!---------------------------------------------------------------------- 272 ! 273 IF( nn_timing == 1 ) CALL timing_start('ssh_nxt') 279 !!---------------------------------------------------------------------- 280 ! 281 IF( nn_timing == 1 ) CALL timing_start('ssh_swp') 274 282 ! 275 283 IF( kt == nit000 ) THEN 276 284 IF(lwp) WRITE(numout,*) 277 IF(lwp) WRITE(numout,*) 'ssh_ nxt : next sea surface height (Asselin time filter + swap)'285 IF(lwp) WRITE(numout,*) 'ssh_swp : Asselin time filter and swap of sea surface height' 278 286 IF(lwp) WRITE(numout,*) '~~~~~~~ ' 279 287 ENDIF 280 288 281 ! !--------------------------! 282 IF( lk_vvl ) THEN ! Variable volume levels ! (ssh at t-, u-, v, f-points) 283 ! !--------------------------! 284 ! 285 IF( neuler == 0 .AND. kt == nit000 ) THEN !** Euler time-stepping at first time-step : no filter 286 sshn (:,:) = ssha (:,:) ! now <-- after (before already = now) 287 sshu_n(:,:) = sshu_a(:,:) 288 sshv_n(:,:) = sshv_a(:,:) 289 DO jj = 1, jpjm1 ! ssh now at f-point 290 DO ji = 1, jpim1 ! NO Vector Opt. 291 sshf_n(ji,jj) = 0.5 * umask(ji,jj,1) * umask(ji,jj+1,1) & 292 & / ( e1f(ji,jj ) * e2f(ji,jj ) ) & 293 & * ( e1u(ji,jj ) * e2u(ji,jj ) * sshu_n(ji,jj ) & 294 & + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) ) 295 END DO 296 END DO 297 CALL lbc_lnk( sshf_n, 'F', 1. ) ! Boundaries conditions 298 ! 299 ELSE !** Leap-Frog time-stepping: Asselin filter + swap 300 zec = atfp * rdt / rau0 301 DO jj = 1, jpj 302 DO ji = 1, jpi ! before <-- now filtered 303 sshb (ji,jj) = sshn (ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) ) & 304 & - zec * ( emp_b(ji,jj) - emp(ji,jj) ) * tmask(ji,jj,1) 305 sshn (ji,jj) = ssha (ji,jj) ! now <-- after 306 sshu_n(ji,jj) = sshu_a(ji,jj) 307 sshv_n(ji,jj) = sshv_a(ji,jj) 308 END DO 309 END DO 310 DO jj = 1, jpjm1 ! ssh now at f-point 311 DO ji = 1, jpim1 ! NO Vector Opt. 312 sshf_n(ji,jj) = 0.5 * umask(ji,jj,1) * umask(ji,jj+1,1) & 313 & / ( e1f(ji,jj ) * e2f(ji,jj ) ) & 314 & * ( e1u(ji,jj ) * e2u(ji,jj ) * sshu_n(ji,jj ) & 315 & + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) ) 316 END DO 317 END DO 318 CALL lbc_lnk( sshf_n, 'F', 1. ) ! Boundaries conditions 319 ! 320 DO jj = 1, jpjm1 ! ssh before at u- & v-points 321 DO ji = 1, jpim1 ! NO Vector Opt. 322 sshu_b(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji ,jj) * e2u(ji ,jj) ) & 323 & * ( e1t(ji ,jj) * e2t(ji ,jj) * sshb(ji ,jj) & 324 & + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) ) 325 sshv_b(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj ) * e2v(ji,jj ) ) & 326 & * ( e1t(ji,jj ) * e2t(ji,jj ) * sshb(ji,jj ) & 327 & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) ) 328 END DO 329 END DO 330 CALL lbc_lnk( sshu_b, 'U', 1. ) 331 CALL lbc_lnk( sshv_b, 'V', 1. ) ! Boundaries conditions 332 ! 333 ENDIF 334 ! !--------------------------! 335 ELSE ! fixed levels ! (ssh at t-point only) 336 ! !--------------------------! 337 ! 338 IF( neuler == 0 .AND. kt == nit000 ) THEN !** Euler time-stepping at first time-step : no filter 339 sshn(:,:) = ssha(:,:) ! now <-- after (before already = now) 340 ! 341 ELSE ! Leap-Frog time-stepping: Asselin filter + swap 342 DO jj = 1, jpj 343 DO ji = 1, jpi ! before <-- now filtered 344 sshb(ji,jj) = sshn(ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) ) 345 sshn(ji,jj) = ssha(ji,jj) ! now <-- after 346 END DO 347 END DO 348 ENDIF 349 ! 289 # if defined key_dynspg_ts 290 IF( ( neuler == 0 .AND. kt == nit000 ) .OR. ln_bt_fw ) THEN !** Euler time-stepping: no filter 291 # else 292 IF ( neuler == 0 .AND. kt == nit000 ) THEN !** Euler time-stepping at first time-step : no filter 293 #endif 294 sshb(:,:) = sshn(:,:) ! before <-- now 295 sshn(:,:) = ssha(:,:) ! now <-- after (before already = now) 296 ELSE !** Leap-Frog time-stepping: Asselin filter + swap 297 sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) ) ! before <-- now filtered 298 IF( lk_vvl ) sshb(:,:) = sshb(:,:) - atfp * rdt / rau0 * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1) 299 sshn(:,:) = ssha(:,:) ! now <-- after 350 300 ENDIF 351 301 ! … … 357 307 IF(ln_ctl) CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb - : ', mask1=tmask, ovlap=1 ) 358 308 ! 359 IF( nn_timing == 1 ) CALL timing_stop('ssh_ nxt')360 ! 361 END SUBROUTINE ssh_ nxt309 IF( nn_timing == 1 ) CALL timing_stop('ssh_swp') 310 ! 311 END SUBROUTINE ssh_swp 362 312 363 313 !!======================================================================
Note: See TracChangeset
for help on using the changeset viewer.
