Changeset 592 for trunk/NEMO/OPA_SRC/DYN/dynspg_ts.F90
- Timestamp:
- 2007-02-09T10:15:25+01:00 (17 years ago)
- File:
-
- 1 edited
-
trunk/NEMO/OPA_SRC/DYN/dynspg_ts.F90 (modified) (12 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r575 r592 34 34 USE iom 35 35 USE restart ! only for lrst_oce 36 USE domvvl ! variable volume 36 37 37 38 IMPLICIT NONE … … 102 103 zsshb_e, zub_e, zvb_e, & ! " " 103 104 zun_e, zvn_e ! " " 105 !! Variable volume 106 REAL(wp), DIMENSION(jpi,jpj) :: & 107 zspgu_1, zspgv_1, zsshun_e, zsshvn_e, hu_e, hv_e ! 2D workspace 108 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zfse3un_e, zfse3vn_e ! 3D workspace 104 109 !!---------------------------------------------------------------------- 105 110 … … 137 142 ENDIF 138 143 ! 144 ENDIF 145 146 ! Substract the surface pressure gradient from the momentum 147 ! --------------------------------------------------------- 148 IF( lk_vvl ) THEN ! Variable volume 149 150 ! 0. Local initialization 151 ! ----------------------- 152 zspgu_1(:,:) = 0.e0 153 zspgv_1(:,:) = 0.e0 154 155 ! 1. Surface pressure gradient (now) 156 ! ---------------------------- 157 DO jj = 2, jpjm1 158 DO ji = fs_2, fs_jpim1 ! vector opt. 159 zspgu_1(ji,jj) = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn(ji+1,jj ) & 160 & - ( rhd(ji ,jj ,1) + 1 ) * sshn(ji ,jj ) ) / e1u(ji,jj) 161 zspgv_1(ji,jj) = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn(ji ,jj+1) & 162 & - ( rhd(ji ,jj ,1) + 1 ) * sshn(ji ,jj ) ) / e2v(ji,jj) 163 END DO 164 END DO 165 166 ! 2. Substract the surface pressure trend from the general trend 167 ! ------------------------------------------------------ 168 DO jk = 1, jpkm1 169 DO jj = 2, jpjm1 170 DO ji = fs_2, fs_jpim1 ! vector opt. 171 ua(ji,jj,jk) = ua(ji,jj,jk) - zspgu_1(ji,jj) 172 va(ji,jj,jk) = va(ji,jj,jk) - zspgv_1(ji,jj) 173 END DO 174 END DO 175 END DO 176 139 177 ENDIF 140 178 … … 269 307 zua_b (:,:) = un_b (:,:) 270 308 zva_b (:,:) = vn_b (:,:) 309 IF( lk_vvl ) THEN 310 zsshun_e(:,:) = sshu (:,:) ! (barotropic) sea surface height (now) at u-point 311 zsshvn_e(:,:) = sshv (:,:) ! (barotropic) sea surface height (now) at v-point 312 hu_e(:,:) = hu (:,:) ! (barotropic) ocean depth at u-point 313 hv_e(:,:) = hv (:,:) ! (barotropic) ocean depth at v-point 314 zfse3un_e(:,:,:) = fse3u(:,:,:) ! (barotropic) scale factors at u-point 315 zfse3un_e(:,:,:) = fse3v(:,:,:) ! (barotropic) scale factors at v-point 316 ENDIF 271 317 272 318 ! set ssh corrections to 0 … … 330 376 DO ji = fs_2, fs_jpim1 ! vector opt. 331 377 ! surface pressure gradient 332 zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj) 333 zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj) 378 IF( lk_vvl) THEN 379 zspgu = -grav * ( ( rhd(ji+1,jj,1) + 1 ) * sshn_e(ji+1,jj) & 380 & - ( rhd(ji,jj,1) + 1 ) * sshn_e(ji,jj) ) * hu_e(ji,jj) / e1u(ji,jj) 381 zspgv = -grav * ( ( rhd(ji,jj+1,1) + 1 ) * sshn_e(ji,jj+1) & 382 & - ( rhd(ji,jj,1) + 1 ) * sshn_e(ji,jj) ) * hv_e(ji,jj) / e2v(ji,jj) 383 ELSE 384 zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj) 385 zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj) 386 ENDIF 334 387 ! energy conserving formulation for planetary vorticity term 335 388 zy1 = ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj) … … 349 402 DO ji = fs_2, fs_jpim1 ! vector opt. 350 403 ! surface pressure gradient 351 zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj) 352 zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj) 404 IF( lk_vvl) THEN 405 zspgu = -grav * ( ( rhd(ji+1,jj,1) + 1 ) * sshn_e(ji+1,jj) & 406 & - ( rhd(ji,jj,1) + 1 ) * sshn_e(ji,jj) ) * hu_e(ji,jj) / e1u(ji,jj) 407 zspgv = -grav * ( ( rhd(ji,jj+1,1) + 1 ) * sshn_e(ji,jj+1) & 408 & - ( rhd(ji,jj,1) + 1 ) * sshn_e(ji,jj) ) * hv_e(ji,jj) / e2v(ji,jj) 409 ELSE 410 zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj) 411 zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj) 412 ENDIF 353 413 ! enstrophy conserving formulation for planetary vorticity term 354 414 zy1 = zfact1 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & … … 369 429 DO ji = fs_2, fs_jpim1 ! vector opt. 370 430 ! surface pressure gradient 371 zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj) 372 zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj) 431 IF( lk_vvl) THEN 432 zspgu = -grav * ( ( rhd(ji+1,jj,1) + 1 ) * sshn_e(ji+1,jj) & 433 & - ( rhd(ji,jj,1) + 1 ) * sshn_e(ji,jj) ) * hu_e(ji,jj) / e1u(ji,jj) 434 zspgv = -grav * ( ( rhd(ji,jj+1,1) + 1 ) * sshn_e(ji,jj+1) & 435 & - ( rhd(ji,jj,1) + 1 ) * sshn_e(ji,jj) ) * hv_e(ji,jj) / e2v(ji,jj) 436 ELSE 437 zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj) 438 zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj) 439 ENDIF 373 440 ! energy/enstrophy conserving formulation for planetary vorticity term 374 441 zcubt = + zfac25 / e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) + ftnw(ji+1,jj) * zwy(ji+1,jj ) & … … 403 470 ! Time filter and swap of dynamics arrays 404 471 ! --------------------------------------- 405 IF( neuler == 0 .AND. kt == nit000) THEN ! Euler (forward) time stepping472 IF( neuler == 0 .AND. jit == 1 ) THEN ! Euler (forward) time stepping 406 473 zsshb_e(:,:) = sshn_e(:,:) 407 474 zub_e (:,:) = zun_e (:,:) … … 419 486 ENDIF 420 487 488 IF( lk_vvl ) THEN ! Variable volume 489 490 ! Sea surface elevation 491 ! --------------------- 492 DO jj = 1, jpjm1 493 DO ji = 1,jpim1 494 495 ! Sea Surface Height at u-point before 496 zsshun_e(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) & 497 & * ( e1t(ji ,jj) * e2t(ji ,jj) * sshn_e(ji ,jj) & 498 & + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn_e(ji+1,jj) ) 499 500 ! Sea Surface Height at v-point before 501 zsshvn_e(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) & 502 & * ( e1t(ji,jj ) * e2t(ji,jj ) * sshn_e(ji,jj ) & 503 & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn_e(ji,jj+1) ) 504 505 END DO 506 END DO 507 508 ! Boundaries conditions 509 CALL lbc_lnk( zsshun_e, 'U', 1. ) ; CALL lbc_lnk( zsshvn_e, 'V', 1. ) 510 511 ! Scale factors at before and after time step 512 ! ------------------------------------------- 513 DO jk = 1, jpkm1 514 zfse3un_e(:,:,jk) = fsve3u(:,:,jk) * ( 1 + zsshun_e(:,:) * muu(:,:,jk) ) 515 zfse3vn_e(:,:,jk) = fsve3v(:,:,jk) * ( 1 + zsshvn_e(:,:) * muv(:,:,jk) ) 516 END DO 517 518 ! Ocean depth at U- and V-points 519 hu_e(:,:) = 0.e0 520 hv_e(:,:) = 0.e0 521 522 DO jk = 1, jpk 523 hu_e(:,:) = hu_e(:,:) + zfse3un_e(:,:,jk) * umask(:,:,jk) 524 hv_e(:,:) = hv_e(:,:) + zfse3vn_e(:,:,jk) * vmask(:,:,jk) 525 END DO 526 527 ENDIF ! End variable volume case 528 421 529 ! ! ==================== ! 422 530 END DO ! end loop ! … … 444 552 ! Horizontal divergence of time averaged barotropic transports 445 553 !------------------------------------------------------------- 446 DO jj = 2, jpjm1 447 DO ji = fs_2, fs_jpim1 ! vector opt. 448 zhdiv(ji,jj) = ( e2u(ji,jj) * un_b(ji,jj) - e2u(ji-1,jj ) * un_b(ji-1,jj ) & 449 & +e1v(ji,jj) * vn_b(ji,jj) - e1v(ji ,jj-1) * vn_b(ji ,jj-1) ) & 450 & / ( e1t(ji,jj) * e2t(ji,jj) ) 451 END DO 452 END DO 453 454 #if defined key_obc 554 IF( .NOT. lk_vvl ) THEN 555 DO jj = 2, jpjm1 556 DO ji = fs_2, fs_jpim1 ! vector opt. 557 zhdiv(ji,jj) = ( e2u(ji,jj) * un_b(ji,jj) - e2u(ji-1,jj ) * un_b(ji-1,jj ) & 558 & +e1v(ji,jj) * vn_b(ji,jj) - e1v(ji ,jj-1) * vn_b(ji ,jj-1) ) & 559 & / ( e1t(ji,jj) * e2t(ji,jj) ) 560 END DO 561 END DO 562 ENDIF 563 564 #if defined key_obc && ! defined key_vvl 455 565 ! open boundaries (div must be zero behind the open boundary) 456 566 ! mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column … … 463 573 ! sea surface height 464 574 !------------------- 465 sshb(:,:) = sshn(:,:) 466 sshn(:,:) = ( sshb_b(:,:) - z2dt_b * ( zraur * emp(:,:) + zhdiv(:,:) ) ) * tmask(:,:,1) 575 IF( lk_vvl ) THEN 576 sshbb(:,:) = sshb(:,:) 577 sshb (:,:) = sshn(:,:) 578 sshn (:,:) = ssha(:,:) 579 ELSE 580 sshb (:,:) = sshn(:,:) 581 sshn(:,:) = ( sshb_b(:,:) - z2dt_b * ( zraur * emp(:,:) + zhdiv(:,:) ) ) * tmask(:,:,1) 582 ENDIF 467 583 468 584 ! ... Boundary conditions on sshn … … 477 593 ! ------------------------------------------ 478 594 sshb_b(:,:) = sshn_b (:,:) 595 IF ( kt == nit000 ) zssha_b(:,:) = sshn(:,:) 479 596 sshn_b(:,:) = zssha_b(:,:) 480 597 un_b (:,:) = zua_b (:,:)
Note: See TracChangeset
for help on using the changeset viewer.
