- Timestamp:
- 2013-07-11T15:59:14+02:00 (11 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
r3764 r3970 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 !! 3.5 ! 2013-07 (J. Chanut) Compliant with time splitting changes 10 11 !!---------------------------------------------------------------------- 11 12 … … 28 29 USE obc_oce 29 30 USE bdy_oce 31 USE bdy_par 32 USE bdydyn2d ! bdy_ssh routine 30 33 USE diaar5, ONLY: lk_diaar5 31 34 USE iom 32 USE sbcrnf , ONLY: h_rnf, nk_rnf! River runoff35 USE sbcrnf ! River runoff 33 36 #if defined key_agrif 34 37 USE agrif_opa_update … … 40 43 USE wrk_nemo ! Memory Allocation 41 44 USE timing ! Timing 45 46 #if defined key_dynspg_ts 47 ! jchanut: Needed modules for dynamics update: 48 USE eosbn2 ! equation of state (eos_bn2 routine) 49 USE zpshde ! partial step: hor. derivative (zps_hde routine) 50 USE dynadv ! advection (dyn_adv routine) 51 USE dynvor ! vorticity term (dyn_vor routine) 52 USE dynhpg ! hydrostatic pressure grad. (dyn_hpg routine) 53 USE dynldf ! lateral momentum diffusion (dyn_ldf routine) 54 USE dynspg_oce ! surface pressure gradient (dyn_spg routine) 55 USE dynspg ! surface pressure gradient (dyn_spg routine) 56 USE dynnept ! simp. form of Neptune effect(dyn_nept_cor routine) 57 USE asminc ! assimilation increments (dyn_asm_inc routine) 58 USE bdydyn3d ! (bdy_dyn3d_dmp routine) 59 USE dynspg_ts ! for advective velocities issued from time splitting 60 USE zdfbfr ! Bottom friction 61 #if defined key_agrif 62 USE agrif_opa_sponge ! Momemtum and tracers sponges 63 #endif 64 #endif 42 65 43 66 IMPLICIT NONE … … 79 102 ! 80 103 INTEGER :: ji, jj, jk ! dummy loop indices 104 INTEGER :: indic 81 105 REAL(wp) :: zcoefu, zcoefv, zcoeff, z2dt, z1_2dt, z1_rau0 ! local scalars 82 106 REAL(wp), POINTER, DIMENSION(:,: ) :: z2d, zhdiv … … 146 170 hu(:,:) = hu_0(:,:) + sshu_n(:,:) ! now ocean depth (at u- and v-points) 147 171 hv(:,:) = hv_0(:,:) + sshv_n(:,:) 172 ! bg jchanut tschanges 173 #if defined key_dynspg_ts 174 ht(:,:) = ht_0(:,:) + sshn(:,:) ! now ocean depth (at t- and f-points) 175 hf(:,:) = hf_0(:,:) + sshf_n(:,:) 176 #endif 177 ! end jchanut tschanges 148 178 ! ! now masked inverse of the ocean depth (at u- and v-points) 149 179 hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1._wp - umask(:,:,1) ) … … 180 210 #endif 181 211 #if defined key_bdy 182 ssha(:,:) = ssha(:,:) * bdytmask(:,:) 183 CALL lbc_lnk( ssha, 'T', 1. ) ! absolutly compulsory !! (jmm) 184 #endif 212 ! bg jchanut tschanges 213 ! These lines are not necessary with time splitting since 214 ! boundary condition on sea level is set during ts loop 215 IF (lk_bdy) THEN 216 CALL lbc_lnk( ssha, 'T', 1. ) ! Not sure that's necessary 217 CALL bdy_ssh( ssha ) ! Duplicate sea level across open boundaries 218 ENDIF 219 #endif 220 ! end jchanut tschanges 185 221 #if defined key_asminc 186 222 ! ! Include the IAU weighted SSH increment … … 190 226 ENDIF 191 227 #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, jpjm1195 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) )202 END DO203 END DO204 CALL lbc_lnk( sshu_a, 'U', 1. ) ; CALL lbc_lnk( sshv_a, 'V', 1. ) ! Boundaries conditions205 ENDIF206 207 228 ! !------------------------------! 208 229 ! ! Now Vertical Velocity ! … … 219 240 END DO 220 241 242 ! bg jchanut tschanges 243 #if defined key_dynspg_ts 244 ! In case the time splitting case, update most of all momentum trends here: 245 ! Note that the computation of vertical velocity above, hence "after" sea level is necessary 246 ! to compute momentum advection for the rhs of barotropic loop: 247 CALL eos ( tsn, rhd, rhop ) ! now in situ density for hpg computation 248 249 IF( ln_zps ) CALL zps_hde( kt, jpts, tsn, gtsu, gtsv, & ! zps: now hor. derivative 250 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 251 CALL zdf_bfr( kt ) ! bottom friction (if quadratic) 252 253 ua(:,:,:) = 0.e0 ! set dynamics trends to zero 254 va(:,:,:) = 0.e0 255 IF( ln_asmiau .AND. & 256 & ln_dyninc ) CALL dyn_asm_inc( kt ) ! apply dynamics assimilation increment 257 IF( ln_neptsimp ) CALL dyn_nept_cor( kt ) ! subtract Neptune velocities (simplified) 258 IF( lk_bdy ) CALL bdy_dyn3d_dmp( kt ) ! bdy damping trends 259 CALL dyn_adv( kt ) ! advection (vector or flux form) 260 CALL dyn_vor( kt ) ! vorticity term including Coriolis 261 CALL dyn_ldf( kt ) ! lateral mixing 262 IF( ln_neptsimp ) CALL dyn_nept_cor( kt ) ! add Neptune velocities (simplified) 263 #if defined key_agrif 264 IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_dyn ! momentum sponge 265 #endif 266 CALL dyn_hpg( kt ) ! horizontal gradient of Hydrostatic pressure 267 CALL dyn_spg( kt, indic ) ! surface pressure gradient 268 269 ua_bak(:,:,:) = ua(:,:,:) ! save next velocities (not trends !) 270 va_bak(:,:,:) = va(:,:,:) 271 272 ! At this stage: 273 ! 1) ssha, sshu_a, sshv_a have been updated. 274 ! 2) Time averaged barotropic velocities around step kt+1 (ua_b, va_b) as well as 275 ! "advective" barotropic velocities (un_adv, vn_adv) at step kt in agreement 276 ! with continuity equation are available. 277 ! 3) 3d velocities (ua, va) hold ua_b, va_b issued from time splitting as barotropic components. 278 ! Below => Update "now" velocities, divergence, then vertical velocity 279 ! NB: hdivn is NOT updated such that hdivb is not updated also. This means that horizontal 280 ! momentum diffusion is still performed on Instantaneous barotropic arrays. I experencied 281 ! some issues with UBS with the current method. Uncomment "divup" below to update the divergence. 282 ! 283 CALL wrk_alloc( jpi,jpj,jpk, z3d ) 284 ! 285 DO jk = 1, jpkm1 286 ! Correct velocities: 287 un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:) - un_b(:,:) )*umask(:,:,jk) 288 vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:) - vn_b(:,:) )*vmask(:,:,jk) 289 ! 290 ! Compute divergence: 291 DO jj = 2, jpjm1 292 DO ji = fs_2, fs_jpim1 ! vector opt. 293 z3d(ji,jj,jk) = & 294 ( 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) & 295 + 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) ) & 296 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 297 END DO 298 END DO 299 END DO 300 ! 301 IF( ln_rnf ) CALL sbc_rnf_div( z3d ) ! runoffs (update divergence) 302 ! 303 ! 304 ! Set new vertical velocities: 305 DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence 306 ! - ML - need 3 lines here because replacement of fse3t by its expression yields too long lines otherwise 307 wn(:,:,jk) = wn(:,:,jk+1) - fse3t_n(:,:,jk) * z3d(:,:,jk) & 308 & - (fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) & 309 & * tmask(:,:,jk) * z1_2dt 310 #if defined key_bdy 311 ! JC: line below is purely cosmetic I guess: 312 wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:) 313 #endif 314 END DO 315 ! 316 CALL wrk_dealloc( jpi,jpj,jpk, z3d ) 317 318 !! bg divup 319 !! ! 320 !! DO jk = 1, jpkm1 321 !! ! Correct velocities: 322 !! un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:) - un_b(:,:) )*umask(:,:,jk) 323 !! vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:) - vn_b(:,:) )*vmask(:,:,jk) 324 !! ! 325 !! END DO 326 !! ! 327 !! ! 328 !! CALL div_cur( kt ) ! Horizontal divergence & Relative vorticity 329 !! ! 330 !! ! Set new vertical velocities: 331 !! DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence 332 !! ! - ML - need 3 lines here because replacement of fse3t by its expression yields too long lines otherwise 333 !! wn(:,:,jk) = wn(:,:,jk+1) - fse3t_n(:,:,jk) * hdivn(:,:,jk) & 334 !! & - (fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) & 335 !! & * tmask(:,:,jk) * z1_2dt 336 !!#if defined key_bdy 337 !! ! JC: line below is purely cosmetic I guess: 338 !! wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:) 339 !!#endif 340 !! END DO 341 !! end divup 342 ! 343 ! 344 ! End of time splitting case 345 #else 346 ! ! Sea Surface Height at u-,v- and f-points (vvl case only) 347 IF( lk_vvl ) THEN ! (required only in key_vvl case) 348 DO jj = 1, jpjm1 349 DO ji = 1, jpim1 ! NO Vector Opt. 350 sshu_a(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji ,jj) * e2u(ji ,jj) ) & 351 & * ( e1t(ji ,jj) * e2t(ji ,jj) * ssha(ji ,jj) & 352 & + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) ) 353 sshv_a(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj ) * e2v(ji,jj ) ) & 354 & * ( e1t(ji,jj ) * e2t(ji,jj ) * ssha(ji,jj ) & 355 & + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) ) 356 END DO 357 END DO 358 CALL lbc_lnk( sshu_a, 'U', 1. ) ; CALL lbc_lnk( sshv_a, 'V', 1. ) ! Boundaries conditions 359 ENDIF 360 #endif 221 361 ! !------------------------------! 222 362 ! ! outputs ! … … 283 423 ! !--------------------------! 284 424 ! 285 IF( neuler == 0 .AND. kt == nit000 ) THEN !** Euler time-stepping at first time-step : no filter 425 #if defined key_dynspg_ts 426 IF(( neuler == 0 .AND. kt == nit000 ).OR.(lk_dynspg_ts.AND.ln_bt_fw)) THEN !** Euler time-stepping: no filter 427 #else 428 IF ( neuler == 0 .AND. kt == nit000 ) THEN 429 #endif 430 sshb (:,:) = sshn (:,:) 286 431 sshn (:,:) = ssha (:,:) ! now <-- after (before already = now) 432 sshu_b(:,:) = sshu_n(:,:) 433 sshv_b(:,:) = sshv_n(:,:) 287 434 sshu_n(:,:) = sshu_a(:,:) 288 435 sshv_n(:,:) = sshv_a(:,:) … … 336 483 ! !--------------------------! 337 484 ! 338 IF( neuler == 0 .AND. kt == nit000 ) THEN !** Euler time-stepping at first time-step : no filter 485 #if defined key_dynspg_ts 486 IF(( neuler == 0 .AND. kt == nit000 ).OR.(lk_dynspg_ts.AND.ln_bt_fw)) THEN !** Euler time-stepping: no filter 487 #else 488 IF( neuler == 0 .AND. kt == nit000 ) THEN 489 #endif 490 sshb(:,:) = sshn(:,:) 339 491 sshn(:,:) = ssha(:,:) ! now <-- after (before already = now) 340 492 !
Note: See TracChangeset
for help on using the changeset viewer.
