- Timestamp:
- 2013-11-11T13:01:19+01:00 (11 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
r3896 r4178 33 33 USE diaar5, ONLY: lk_diaar5 34 34 USE iom 35 USE sbcrnf, ONLY: h_rnf, nk_rnf ! River runoff 35 USE sbcrnf, ONLY: h_rnf, nk_rnf, sbc_rnf_div ! River runoff 36 USE dynspg_ts, ONLY: un_b, vn_b, un_adv, vn_adv 37 USE dynspg_oce, ONLY: lk_dynspg_ts 36 38 #if defined key_agrif 37 39 USE agrif_opa_update … … 48 50 49 51 PUBLIC ssh_nxt ! called by step.F90 50 PUBLIC wzv ! called by step.F90 52 PUBLIC wzv_1 ! called by step.F90 53 PUBLIC wzv_2 ! called by step.F90 51 54 PUBLIC ssh_swp ! called by step.F90 52 55 … … 149 152 150 153 151 SUBROUTINE wzv ( kt )152 !!---------------------------------------------------------------------- 153 !! *** ROUTINE wzv ***154 SUBROUTINE wzv_1( kt ) 155 !!---------------------------------------------------------------------- 156 !! *** ROUTINE wzv_1 *** 154 157 !! 155 158 !! ** Purpose : compute the now vertical velocity … … 166 169 ! 167 170 INTEGER, INTENT(in) :: kt ! time step 168 REAL(wp), POINTER, DIMENSION(:,: ) :: z2d 169 REAL(wp), POINTER, DIMENSION(:,:,:) :: z3d, zhdiv 171 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhdiv 170 172 ! 171 173 INTEGER :: ji, jj, jk ! dummy loop indices … … 173 175 !!---------------------------------------------------------------------- 174 176 175 IF( nn_timing == 1 ) CALL timing_start('wzv') 176 ! 177 CALL wrk_alloc( jpi, jpj, z2d ) 178 CALL wrk_alloc( jpi, jpj, jpk, z3d, zhdiv ) 177 IF( nn_timing == 1 ) CALL timing_start('wzv_1') 178 ! 179 CALL wrk_alloc( jpi, jpj, jpk, zhdiv ) 179 180 ! 180 181 IF( kt == nit000 ) THEN 181 182 ! 182 183 IF(lwp) WRITE(numout,*) 183 IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity '184 IF(lwp) WRITE(numout,*) '~~~ '184 IF(lwp) WRITE(numout,*) 'wzv_1 : now vertical velocity ' 185 IF(lwp) WRITE(numout,*) '~~~~~ ' 185 186 ! 186 187 wn(:,:,jpk) = 0._wp ! bottom boundary condition: w=0 (set once for all) … … 224 225 #endif 225 226 227 ! 228 CALL wrk_dealloc( jpi, jpj, jpk, zhdiv ) 229 ! 230 IF( nn_timing == 1 ) CALL timing_stop('wzv_1') 231 232 233 END SUBROUTINE wzv_1 234 235 SUBROUTINE wzv_2( kt ) 236 !!---------------------------------------------------------------------- 237 !! *** ROUTINE wzv_2 *** 238 !! 239 !! ** Purpose : compute the now vertical velocity 240 !! 241 !! ** Method : - Using the incompressibility hypothesis, the vertical 242 !! velocity is computed by integrating the horizontal divergence 243 !! from the bottom to the surface minus the scale factor evolution. 244 !! The boundary conditions are w=0 at the bottom (no flux) and. 245 !! 246 !! ** action : wn : now vertical velocity 247 !! 248 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 249 !!---------------------------------------------------------------------- 250 ! 251 INTEGER, INTENT(in) :: kt ! time step 252 REAL(wp), POINTER, DIMENSION(:,: ) :: z2d 253 REAL(wp), POINTER, DIMENSION(:,:,:) :: z3d, zhdiv 254 ! 255 INTEGER :: ji, jj, jk ! dummy loop indices 256 REAL(wp) :: z1_2dt ! local scalars 257 !!---------------------------------------------------------------------- 258 259 IF( nn_timing == 1 ) CALL timing_start('wzv_2') 260 ! 261 CALL wrk_alloc( jpi, jpj, jpk, z3d, zhdiv ) 262 ! 263 IF( kt == nit000 ) THEN 264 ! 265 IF(lwp) WRITE(numout,*) 266 IF(lwp) WRITE(numout,*) 'wzv_2 : now vertical velocity ' 267 IF(lwp) WRITE(numout,*) '~~~~~ ' 268 ! 269 ENDIF 270 ! !------------------------------! 271 ! ! Now Vertical Velocity ! 272 ! !------------------------------! 273 z1_2dt = 1. / ( 2. * rdt ) ! set time step size (Euler/Leapfrog) 274 IF( neuler == 0 .AND. kt == nit000 ) z1_2dt = 1. / rdt 275 276 IF( lk_dynspg_ts ) THEN 277 ! At this stage: 278 ! 1) vertical scale factors have been updated. 279 ! 2) Time averaged barotropic velocities around step kt+1 (ua_b, va_b) as well as 280 ! "advective" barotropic velocities (un_adv, vn_adv) at step kt in agreement 281 ! with continuity equation are available. 282 ! 3) 3d velocities (ua, va) hold ua_b, va_b issued from time splitting as barotropic components. 283 ! Below => Update "now" velocities, divergence, then vertical velocity 284 ! NB: hdivn is NOT updated such that hdivb is not updated also. This means that horizontal 285 ! momentum diffusion is still performed on Instantaneous barotropic arrays. I experencied 286 ! some issues with UBS with the current method. See "divup" comments in development branch to 287 ! update the divergence fully if necessary (2013/dev_r3867_MERCATOR1_DYN). 288 ! 289 DO jk = 1, jpkm1 290 ! Correct velocities: 291 un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:) - un_b(:,:) )*umask(:,:,jk) 292 vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:) - vn_b(:,:) )*vmask(:,:,jk) 293 ! 294 ! Compute divergence: 295 DO jj = 2, jpjm1 296 DO ji = fs_2, fs_jpim1 ! vector opt. 297 z3d(ji,jj,jk) = & 298 ( e2u(ji,jj)*fse3u(ji,jj,jk) * un(ji,jj,jk) - e2u(ji-1,jj)*fse3u(ji-1,jj,jk) * un(ji-1,jj,jk) & 299 + e1v(ji,jj)*fse3v(ji,jj,jk) * vn(ji,jj,jk) - e1v(ji,jj-1)*fse3v(ji,jj-1,jk) * vn(ji,jj-1,jk) ) & 300 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 301 END DO 302 END DO 303 END DO 304 ! 305 IF( ln_rnf ) CALL sbc_rnf_div( z3d ) ! runoffs (update divergence) 306 ELSE 307 z3d(:,:,:) = hdivn(:,:,:) 308 ENDIF 309 ! 310 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde and layer cases 311 DO jk = 1, jpkm1 312 ! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t) 313 ! - ML - note: computation allready done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way) 314 DO jj = 2, jpjm1 315 DO ji = fs_2, fs_jpim1 ! vector opt. 316 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) ) 317 END DO 318 END DO 319 END DO 320 CALL lbc_lnk(zhdiv, 'T', 1.) ! - ML - Perhaps not necessary: not used for horizontal "connexions" 321 ! ! Is it problematic to have a wrong vertical velocity in boundary cells? 322 ! ! Same question holds for hdivn. Perhaps just for security 323 DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence 324 ! computation of w 325 wn(:,:,jk) = wn(:,:,jk+1) - ( fse3t_n(:,:,jk) * z3d(:,:,jk) + zhdiv(:,:,jk) & 326 & + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) ) * tmask(:,:,jk) 327 END DO 328 ! IF( ln_vvl_layer ) wn(:,:,:) = 0.e0 329 ELSE ! z_star and linear free surface cases 330 DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence 331 ! computation of w 332 wn(:,:,jk) = wn(:,:,jk+1) - ( fse3t_n(:,:,jk) * z3d(:,:,jk) & 333 & + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) ) * tmask(:,:,jk) 334 END DO 335 ENDIF 336 337 #if defined key_bdy 338 wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:) 339 #endif 340 226 341 ! !------------------------------! 227 342 ! ! outputs ! … … 229 344 CALL iom_put( "woce", wn ) ! vertical velocity 230 345 IF( lk_diaar5 ) THEN ! vertical mass transport & its square value 346 CALL wrk_alloc( jpi, jpj, z2d ) 231 347 ! Caution: in the VVL case, it only correponds to the baroclinic mass transport. 232 348 z2d(:,:) = rau0 * e12t(:,:) … … 236 352 CALL iom_put( "w_masstr" , z3d ) 237 353 CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) ) 238 ENDIF239 !240 CALL wrk_dealloc( jpi, jpj, z2d )354 CALL wrk_dealloc( jpi, jpj, z2d ) 355 ENDIF 356 ! 241 357 CALL wrk_dealloc( jpi, jpj, jpk, z3d, zhdiv ) 242 358 ! 243 IF( nn_timing == 1 ) CALL timing_stop('wzv') 244 245 246 END SUBROUTINE wzv 247 359 IF( nn_timing == 1 ) CALL timing_stop('wzv_2') 360 361 362 END SUBROUTINE wzv_2 248 363 249 364 SUBROUTINE ssh_swp( kt )
Note: See TracChangeset
for help on using the changeset viewer.
