- Timestamp:
- 2015-10-30T16:34:24+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r5656 r5842 41 41 USE timing ! Timing 42 42 USE sbcapr ! surface boundary condition: atmospheric pressure 43 USE wadlmt ! wetting/drying flux limter 43 44 USE dynadv, ONLY: ln_dynadv_vec 44 45 #if defined key_agrif … … 142 143 LOGICAL :: ll_fw_start ! if true, forward integration 143 144 LOGICAL :: ll_init ! if true, special startup of 2d equations 145 LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables used in W/D 144 146 INTEGER :: ji, jj, jk, jn ! dummy loop indices 145 147 INTEGER :: ikbu, ikbv, noffset ! local integers … … 157 159 REAL(wp), POINTER, DIMENSION(:,:) :: zsshu_a, zsshv_a 158 160 REAL(wp), POINTER, DIMENSION(:,:) :: zhf 161 REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy ! Wetting/Dying gravity filter coef. 162 REAL(wp), POINTER, DIMENSION(:,:) :: wduflt1, wdvflt1 ! Wetting/Dying velocity filter coef. 159 163 !!---------------------------------------------------------------------- 164 165 ! tempoary debugging integer 166 INTEGER :: jidbg, jjdbg 167 jidbg = 163; jjdbg = 179 160 168 ! 161 169 IF( nn_timing == 1 ) CALL timing_start('dyn_spg_ts') … … 168 176 CALL wrk_alloc( jpi, jpj, zsshu_a, zsshv_a ) 169 177 CALL wrk_alloc( jpi, jpj, zhf ) 178 IF(ln_wd) CALL wrk_alloc( jpi, jpj, zcpx, zcpy, wduflt1, wdvflt1 ) 170 179 ! 171 180 ! !* Local constant initialization … … 379 388 ! ! ---------------------------------------------------- 380 389 IF( lk_vvl ) THEN ! Variable volume : remove surface pressure gradient 381 DO jj = 2, jpjm1 382 DO ji = fs_2, fs_jpim1 ! vector opt. 383 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) / e1u(ji,jj) 384 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) / e2v(ji,jj) 385 END DO 386 END DO 387 ENDIF 388 390 IF(ln_wd) THEN ! Calculating and applying W/D gravity filters 391 DO jj = 2, jpjm1 392 DO ji = 2, jpim1 393 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) & 394 & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) & 395 & > rn_wdmin1 + rn_wdmin2 396 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) + & 397 & rn_wdmin1 + rn_wdmin2 398 IF(ll_tmp1) THEN 399 zcpx(ji,jj) = 1.0_wp 400 wduflt1(ji,jj) = 1.0_wp 401 ELSE IF(ll_tmp2) THEN 402 ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 403 zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /& 404 & (sshn(ji+1,jj) - sshn(ji,jj))) 405 wduflt1(ji,jj) = 1.0_wp 406 ELSE 407 zcpx(ji,jj) = 0._wp 408 wduflt1(ji,jj) = 0.0_wp 409 END IF 410 411 !for w/d debugging, delete it when finished 412 ! zcpx(ji,jj) = 0._wp 413 ! end debugging part 414 415 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) & 416 & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) & 417 & > rn_wdmin1 + rn_wdmin2 418 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) + & 419 & rn_wdmin1 + rn_wdmin2 420 IF(ll_tmp1) THEN 421 zcpy(ji,jj) = 1.0_wp 422 wdvflt1(ji,jj) = 1.0_wp 423 ELSE IF(ll_tmp2) THEN 424 ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 425 zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) /& 426 & (sshn(ji,jj+1) - sshn(ji,jj))) 427 wdvflt1(ji,jj) = 1.0_wp 428 ELSE 429 zcpy(ji,jj) = 0._wp 430 wdvflt1(ji,jj) = 0.0_wp 431 END IF 432 433 !for w/d debugging, delete it when finished 434 ! zcpy(ji,jj) = 0._wp 435 ! end debugging part 436 END DO 437 END DO 438 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 439 ENDIF 440 441 IF(ln_wd) THEN ! Calculating and applying W/D gravity filters 442 DO jj = 2, jpjm1 443 DO ji = 2, jpim1 444 zu_trd(ji,jj) = (zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) / & 445 & e1u(ji,jj)) * zcpx(ji,jj) 446 zv_trd(ji,jj) = (zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) / & 447 & e2v(ji,jj)) * zcpy(ji,jj) 448 IF(ji+nimpp-1 == jidbg .and. jj+njmpp-1 == jjdbg) THEN 449 write(333,*) 'in spt_ts: zutrdx, zutrdy, hpgx, hpgy, zcpx, zcpy' 450 write(333,*) zu_trd(ji,jj), zv_trd(ji,jj), & 451 & -grav*(sshn(ji+1,jj)-sshn(ji,jj))/e1u(ji,jj)*zcpx(ji,jj), & 452 & -grav*(sshn(ji,jj+1)-sshn(ji,jj))/e2v(ji,jj)*zcpy(ji,jj), & 453 & zcpx(ji,jj), zcpy(ji,jj) 454 write(334,*) 'in spt_ts: sshn_ij, sshn_ip1j, sshn_ijp1, bathy_ij, bathy_ip1j, bathy_ijp1' 455 write(334,*) sshn(ji,jj), sshn(ji+1,jj), sshn(ji,jj+1), bathy(ji,jj), bathy(ji+1,jj), bathy(ji,jj+1) 456 END IF 457 END DO 458 END DO 459 ELSE 460 DO jj = 2, jpjm1 461 DO ji = fs_2, fs_jpim1 ! vector opt. 462 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) / e1u(ji,jj) 463 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) / e2v(ji,jj) 464 END DO 465 END DO 466 END IF 467 468 ENDIF 469 470 IF(ln_wd) THEN 471 DO jj = 2, jpjm1 ! Remove coriolis term (and possibly spg) from barotropic trend 472 DO ji = fs_2, fs_jpim1 473 zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * umask(ji,jj,1) * wduflt1(ji,jj) 474 zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * vmask(ji,jj,1) * wdvflt1(ji,jj) 475 END DO 476 END DO 477 ELSE 389 478 DO jj = 2, jpjm1 ! Remove coriolis term (and possibly spg) from barotropic trend 390 479 DO ji = fs_2, fs_jpim1 … … 392 481 zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * vmask(ji,jj,1) 393 482 END DO 394 END DO 483 END DO 484 END IF 395 485 ! 396 486 ! ! Add bottom stress contribution from baroclinic velocities: … … 416 506 ! 417 507 ! Note that the "unclipped" bottom friction parameter is used even with explicit drag 418 zu_frc(:,:) = zu_frc(:,:) + hur(:,:) * bfrua(:,:) * zwx(:,:) 419 zv_frc(:,:) = zv_frc(:,:) + hvr(:,:) * bfrva(:,:) * zwy(:,:) 420 ! 508 IF(ln_wd) THEN 509 zu_frc(:,:) = zu_frc(:,:) + MAX(hur(:,:) * bfrua(:,:),-1._wp / rdtbt) * zwx(:,:) 510 zv_frc(:,:) = zv_frc(:,:) + MAX(hvr(:,:) * bfrva(:,:),-1._wp / rdtbt) * zwy(:,:) 511 ELSE 512 zu_frc(:,:) = zu_frc(:,:) + hur(:,:) * bfrua(:,:) * zwx(:,:) 513 zv_frc(:,:) = zv_frc(:,:) + hvr(:,:) * bfrva(:,:) * zwy(:,:) 514 END IF 515 ! 421 516 IF (ln_bt_fw) THEN ! Add wind forcing 422 517 zu_frc(:,:) = zu_frc(:,:) + zraur * utau(:,:) * hur(:,:) … … 487 582 vb_e (:,:) = 0._wp 488 583 ENDIF 489 ! 490 IF (ln_bt_fw) THEN ! FORWARD integration: start from NOW fields 491 sshn_e(:,:) = sshn (:,:) 492 zun_e (:,:) = un_b (:,:) 584 585 IF(ln_wd) THEN !preserve the positivity of water depth 586 !ssh[b,n,a] should have already been processed for this 587 sshbb_e(:,:) = MAX(sshbb_e(:,:), rn_wdmin1 - bathy(:,:)) 588 sshb_e(:,:) = MAX(sshb_e(:,:) , rn_wdmin1 - bathy(:,:)) 589 ENDIF 590 ! 591 IF (ln_bt_fw) THEN ! FORWARD integration: start from NOW fields 592 sshn_e(:,:) = sshn (:,:) 593 zun_e (:,:) = un_b (:,:) 493 594 zvn_e (:,:) = vn_b (:,:) 494 595 ! … … 561 662 zhup2_e (:,:) = hu_0(:,:) + zwx(:,:) ! Ocean depth at U- and V-points 562 663 zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:) 664 IF(ln_wd) THEN 665 zhup2_e(:,:) = MAX(zhup2_e (:,:), rn_wdmin1) 666 zhvp2_e(:,:) = MAX(zhvp2_e (:,:), rn_wdmin1) 667 END IF 563 668 ELSE 564 669 zhup2_e (:,:) = hu(:,:) … … 599 704 ENDIF 600 705 #endif 706 IF(ln_wd) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt) 601 707 ! 602 708 ! Sum over sub-time-steps to compute advective velocities … … 613 719 END DO 614 720 ssha_e(:,:) = ( sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) ) ) * tmask(:,:,1) 721 IF(ln_wd) ssha_e(:,:) = MAX(ssha_e(:,:), rn_wdmin1 - bathy(:,:)) 615 722 CALL lbc_lnk( ssha_e, 'T', 1._wp ) 616 723 … … 660 767 & + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:) 661 768 769 IF(ln_wd) THEN ! Calculating and applying W/D gravity filters 770 DO jj = 2, jpjm1 771 DO ji = 2, jpim1 772 ll_tmp1 = MIN(zsshp2_e(ji,jj), zsshp2_e(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj))& 773 & .and. MAX(zsshp2_e(ji,jj) + bathy(ji,jj), zsshp2_e(ji+1,jj) + bathy(ji+1,jj)) & 774 & > rn_wdmin1 + rn_wdmin2 775 ll_tmp2 = MAX(zsshp2_e(ji,jj), zsshp2_e(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) + & 776 & rn_wdmin1 + rn_wdmin2 777 IF(ll_tmp1) THEN 778 zcpx(ji,jj) = 1.0_wp 779 wduflt1(ji,jj) = 1.0_wp 780 ELSE IF(ll_tmp2) THEN 781 ! no worries about zsshp2_e(ji+1,jj)-zsshp2_e(ji,jj) = 0, it won't happen ! here 782 zcpx(ji,jj) = ABS((zsshp2_e(ji+1,jj) + bathy(ji+1,jj) - zsshp2_e(ji,jj) - bathy(ji,jj)) /& 783 & (zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj))) 784 wduflt1(ji,jj) = 1.0_wp 785 ELSE 786 zcpx(ji,jj) = 0._wp 787 wduflt1(ji,jj) = 0.0_wp 788 END IF 789 790 !for w/d debugging, delete it when finished 791 ! zcpx(ji,jj) = 0._wp 792 ! end debugging part 793 794 ll_tmp1 = MIN(zsshp2_e(ji,jj), zsshp2_e(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1))& 795 & .and. MAX(zsshp2_e(ji,jj) + bathy(ji,jj), zsshp2_e(ji,jj+1) + bathy(ji,jj+1)) & 796 & > rn_wdmin1 + rn_wdmin2 797 ll_tmp2 = MAX(zsshp2_e(ji,jj), zsshp2_e(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) + & 798 & rn_wdmin1 + rn_wdmin2 799 IF(ll_tmp1) THEN 800 zcpy(ji,jj) = 1.0_wp 801 wdvflt1(ji,jj) = 1.0_wp 802 ELSE IF(ll_tmp2) THEN 803 ! no worries about zsshp2_e(ji,jj+1)-zsshp2_e(ji,jj) = 0, it won't happen ! here 804 zcpy(ji,jj) = ABS((zsshp2_e(ji,jj+1) + bathy(ji,jj+1) - zsshp2_e(ji,jj) - bathy(ji,jj)) /& 805 & (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj))) 806 wdvflt1(ji,jj) = 1.0_wp 807 ELSE 808 zcpy(ji,jj) = 0._wp 809 wdvflt1(ji,jj) = 0.0_wp 810 END IF 811 812 !for w/d debugging, delete it when finished 813 ! zcpy(ji,jj) = 0._wp 814 ! end debugging part 815 END DO 816 END DO 817 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 818 ENDIF 662 819 ! 663 820 ! Compute associated depths at U and V points: … … 676 833 END DO 677 834 END DO 835 836 IF(ln_wd) THEN 837 zhust_e(:,:) = MAX(zhust_e (:,:), rn_wdmin1 ) 838 zhvst_e(:,:) = MAX(zhvst_e (:,:), rn_wdmin1 ) 839 END IF 840 !shall we call lbc_lnk for zhu(v)st_e() here? 841 678 842 ENDIF 679 843 ! … … 742 906 ! 743 907 ! Surface pressure trend: 744 DO jj = 2, jpjm1 745 DO ji = fs_2, fs_jpim1 ! vector opt. 746 ! Add surface pressure gradient 747 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj) 748 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj) 749 zwx(ji,jj) = zu_spg 750 zwy(ji,jj) = zv_spg 751 END DO 752 END DO 908 909 IF(ln_wd) THEN 910 DO jj = 2, jpjm1 911 DO ji = 2, jpim1 912 ! Add surface pressure gradient 913 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj) 914 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj) 915 ! for W/D debugging 916 !zwx(ji,jj) = zu_spg * zcpx(ji,jj) * 0._wp 917 !zwy(ji,jj) = zv_spg * zcpy(ji,jj) * 0._wp 918 IF(ji+nimpp-1 == jidbg .and. jj+njmpp-1 == jjdbg) THEN 919 write(444,*) 'in spt_ts_intg: zu_spg, zv_spg, zcpx, zcpy' 920 write(444,*) zu_spg, zv_spg, zcpx(ji,jj), zcpy(ji,jj) 921 write(445,*) 'in spt_ts_intg: sshn_ij, sshn_ip1j, sshn_ijp1, bathy_ij, bathy_ip1j, bathy_ijp1' 922 write(445,*) zsshp2_e(ji,jj), zsshp2_e(ji+1,jj), zsshp2_e(ji,jj+1), & 923 & bathy(ji,jj), bathy(ji+1,jj), bathy(ji,jj+1) 924 END IF 925 zwx(ji,jj) = zu_spg * zcpx(ji,jj) 926 zwy(ji,jj) = zv_spg * zcpy(ji,jj) 927 END DO 928 END DO 929 ELSE 930 DO jj = 2, jpjm1 931 DO ji = fs_2, fs_jpim1 ! vector opt. 932 ! Add surface pressure gradient 933 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj) 934 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj) 935 zwx(ji,jj) = zu_spg 936 zwy(ji,jj) = zv_spg 937 END DO 938 END DO 939 END IF 940 753 941 ! 754 942 ! Set next velocities: … … 774 962 DO ji = fs_2, fs_jpim1 ! vector opt. 775 963 776 zhura = umask(ji,jj,1)/(hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - umask(ji,jj,1)) 777 zhvra = vmask(ji,jj,1)/(hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - vmask(ji,jj,1)) 778 779 ua_e(ji,jj) = ( hu_e(ji,jj) * zun_e(ji,jj) & 780 & + rdtbt * ( zhust_e(ji,jj) * zwx(ji,jj) & 964 IF(ln_wd) THEN 965 zhura = MAX(hu_0(ji,jj) + zsshu_a(ji,jj), rn_wdmin1) 966 zhvra = MAX(hv_0(ji,jj) + zsshv_a(ji,jj), rn_wdmin1) 967 ELSE 968 zhura = hu_0(ji,jj) + zsshu_a(ji,jj) 969 zhvra = hv_0(ji,jj) + zsshv_a(ji,jj) 970 END IF 971 972 zhura = umask(ji,jj,1)/(zhura + 1._wp - umask(ji,jj,1)) 973 zhvra = vmask(ji,jj,1)/(zhvra + 1._wp - vmask(ji,jj,1)) 974 975 ua_e(ji,jj) = ( hu_e(ji,jj) * zun_e(ji,jj) & 976 & + rdtbt * ( zhust_e(ji,jj) * zwx(ji,jj) & 781 977 & + zhup2_e(ji,jj) * zu_trd(ji,jj) & 782 978 & + hu(ji,jj) * zu_frc(ji,jj) ) & … … 793 989 ! 794 990 IF( lk_vvl ) THEN !* Update ocean depth (variable volume case only) 795 ! ! ---------------------------------------------- 796 hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 797 hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 991 ! ! ---------------------------------------------- 992 IF(ln_wd) THEN 993 hu_e (:,:) = MAX(hu_0(:,:) + zsshu_a(:,:), rn_wdmin1) 994 hv_e (:,:) = MAX(hv_0(:,:) + zsshv_a(:,:), rn_wdmin1) 995 ELSE 996 hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 997 hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 998 END IF 999 798 1000 hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1._wp - umask(:,:,1) ) 799 1001 hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1._wp - vmask(:,:,1) ) … … 925 1127 CALL wrk_dealloc( jpi, jpj, zsshu_a, zsshv_a ) 926 1128 CALL wrk_dealloc( jpi, jpj, zhf ) 1129 IF(ln_wd) CALL wrk_dealloc( jpi, jpj, zcpx, zcpy, wduflt1, wdvflt1 ) 927 1130 ! 928 1131 IF( nn_timing == 1 ) CALL timing_stop('dyn_spg_ts')
Note: See TracChangeset
for help on using the changeset viewer.
