Changeset 6152 for trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
- Timestamp:
- 2015-12-21T23:33:57+01:00 (8 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 (modified) (16 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r6140 r6152 32 32 USE phycst ! physical constants 33 33 USE dynvor ! vorticity term 34 USE wet_dry ! wetting/drying flux limter 34 35 USE bdy_par ! for lk_bdy 35 36 USE bdytides ! open boundary condition data … … 136 137 LOGICAL :: ll_fw_start ! if true, forward integration 137 138 LOGICAL :: ll_init ! if true, special startup of 2d equations 139 LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables used in W/D 138 140 INTEGER :: ji, jj, jk, jn ! dummy loop indices 139 141 INTEGER :: ikbu, ikbv, noffset ! local integers … … 153 155 REAL(wp), POINTER, DIMENSION(:,:) :: zsshu_a, zsshv_a 154 156 REAL(wp), POINTER, DIMENSION(:,:) :: zhf 157 REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy ! Wetting/Dying gravity filter coef. 158 REAL(wp), POINTER, DIMENSION(:,:) :: wduflt1, wdvflt1 ! Wetting/Dying velocity filter coef. 155 159 !!---------------------------------------------------------------------- 156 160 ! … … 162 166 CALL wrk_alloc( jpi,jpj, zwx, zwy, zssh_frc, zu_frc, zv_frc) 163 167 CALL wrk_alloc( jpi,jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e) 164 CALL wrk_alloc( jpi,jpj, zsshu_a, zsshv_a 168 CALL wrk_alloc( jpi,jpj, zsshu_a, zsshv_a ) 165 169 CALL wrk_alloc( jpi,jpj, zhf ) 170 IF( ln_wd ) CALL wrk_alloc( jpi, jpj, zcpx, zcpy, wduflt1, wdvflt1 ) 166 171 ! 167 172 zmdi=1.e+20 ! missing data indicator for masking … … 368 373 ! ! ---------------------------------------------------- 369 374 IF( .NOT.ln_linssh ) THEN ! Variable volume : remove surface pressure gradient 370 DO jj = 2, jpjm1 371 DO ji = fs_2, fs_jpim1 ! vector opt. 372 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) * r1_e1u(ji,jj) 373 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) * r1_e2v(ji,jj) 374 END DO 375 END DO 375 IF( ln_wd ) THEN ! Calculating and applying W/D gravity filters 376 wduflt1(:,:) = 1.0_wp 377 wdvflt1(:,:) = 1.0_wp 378 DO jj = 2, jpjm1 379 DO ji = 2, jpim1 380 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) & 381 & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) & 382 & > rn_wdmin1 + rn_wdmin2 383 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) & 384 & + rn_wdmin1 + rn_wdmin2 385 IF(ll_tmp1) THEN 386 zcpx(ji,jj) = 1.0_wp 387 ELSEIF(ll_tmp2) THEN 388 ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen here 389 zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) & 390 & /(sshn(ji+1,jj) - sshn(ji,jj))) 391 ELSE 392 zcpx(ji,jj) = 0._wp 393 wduflt1(ji,jj) = 0.0_wp 394 END IF 395 396 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) & 397 & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) & 398 & > rn_wdmin1 + rn_wdmin2 399 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) & 400 & + rn_wdmin1 + rn_wdmin2 401 IF(ll_tmp1) THEN 402 zcpy(ji,jj) = 1.0_wp 403 ELSEIF(ll_tmp2) THEN 404 ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen here 405 zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) & 406 & /(sshn(ji,jj+1) - sshn(ji,jj))) 407 ELSE 408 zcpy(ji,jj) = 0._wp 409 wdvflt1(ji,jj) = 0.0_wp 410 ENDIF 411 412 END DO 413 END DO 414 415 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 416 417 DO jj = 2, jpjm1 418 DO ji = 2, jpim1 419 zu_trd(ji,jj) = ( zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) & 420 & * r1_e1u(ji,jj) ) * zcpx(ji,jj) * wduflt1(ji,jj) 421 zv_trd(ji,jj) = ( zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) & 422 & * r1_e2v(ji,jj) ) * zcpy(ji,jj) * wdvflt1(ji,jj) 423 END DO 424 END DO 425 426 ELSE 427 428 DO jj = 2, jpjm1 429 DO ji = fs_2, fs_jpim1 ! vector opt. 430 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) * r1_e1u(ji,jj) 431 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) * r1_e2v(ji,jj) 432 END DO 433 END DO 434 ENDIF 435 376 436 ENDIF 377 437 … … 405 465 ! 406 466 ! Note that the "unclipped" bottom friction parameter is used even with explicit drag 407 zu_frc(:,:) = zu_frc(:,:) + r1_hu_n(:,:) * bfrua(:,:) * zwx(:,:) 408 zv_frc(:,:) = zv_frc(:,:) + r1_hv_n(:,:) * bfrva(:,:) * zwy(:,:) 409 ! 467 IF( ln_wd ) THEN 468 zu_frc(:,:) = zu_frc(:,:) + MAX(r1_hu_n(:,:) * bfrua(:,:),-1._wp / rdtbt) * zwx(:,:) 469 zv_frc(:,:) = zv_frc(:,:) + MAX(r1_hv_n(:,:) * bfrva(:,:),-1._wp / rdtbt) * zwy(:,:) 470 ELSE 471 zu_frc(:,:) = zu_frc(:,:) + r1_hu_n(:,:) * bfrua(:,:) * zwx(:,:) 472 zv_frc(:,:) = zv_frc(:,:) + r1_hv_n(:,:) * bfrva(:,:) * zwy(:,:) 473 END IF 474 ! 410 475 ! ! Add top stress contribution from baroclinic velocities: 411 476 IF (ln_bt_fw) THEN … … 500 565 ub_e (:,:) = 0._wp 501 566 vb_e (:,:) = 0._wp 567 ENDIF 568 569 IF( ln_wd ) THEN !preserve the positivity of water depth 570 !ssh[b,n,a] should have already been processed for this 571 sshbb_e(:,:) = MAX(sshbb_e(:,:), rn_wdmin1 - bathy(:,:)) 572 sshb_e(:,:) = MAX(sshb_e(:,:) , rn_wdmin1 - bathy(:,:)) 502 573 ENDIF 503 574 ! … … 575 646 zhup2_e (:,:) = hu_0(:,:) + zwx(:,:) ! Ocean depth at U- and V-points 576 647 zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:) 648 IF( ln_wd ) THEN 649 zhup2_e(:,:) = MAX(zhup2_e (:,:), rn_wdmin1) 650 zhvp2_e(:,:) = MAX(zhvp2_e (:,:), rn_wdmin1) 651 END IF 577 652 ELSE 578 653 zhup2_e (:,:) = hu_n(:,:) … … 612 687 ENDIF 613 688 #endif 689 IF( ln_wd ) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt) 614 690 ! 615 691 ! Sum over sub-time-steps to compute advective velocities … … 626 702 END DO 627 703 ssha_e(:,:) = ( sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) ) ) * ssmask(:,:) 704 IF( ln_wd ) ssha_e(:,:) = MAX(ssha_e(:,:), rn_wdmin1 - bathy(:,:)) 628 705 CALL lbc_lnk( ssha_e, 'T', 1._wp ) 629 706 … … 672 749 zsshp2_e(:,:) = za0 * ssha_e(:,:) + za1 * sshn_e (:,:) & 673 750 & + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:) 751 IF( ln_wd ) THEN ! Calculating and applying W/D gravity filters 752 wduflt1(:,:) = 1._wp 753 wdvflt1(:,:) = 1._wp 754 DO jj = 2, jpjm1 755 DO ji = 2, jpim1 756 ll_tmp1 = MIN( zsshp2_e(ji,jj), zsshp2_e(ji+1,jj) ) > MAX( -bathy(ji,jj), -bathy(ji+1,jj) ) & 757 & .AND. MAX( zsshp2_e(ji,jj) + bathy(ji,jj), zsshp2_e(ji+1,jj) + bathy(ji+1,jj) ) & 758 & > rn_wdmin1 + rn_wdmin2 759 ll_tmp2 = MAX( zsshp2_e(ji,jj), zsshp2_e(ji+1,jj) ) > MAX( -bathy(ji,jj), -bathy(ji+1,jj) ) & 760 & + rn_wdmin1 + rn_wdmin2 761 IF(ll_tmp1) THEN 762 zcpx(ji,jj) = 1._wp 763 ELSE IF(ll_tmp2) THEN 764 ! no worries about zsshp2_e(ji+1,jj)-zsshp2_e(ji,jj) = 0, it won't happen here 765 zcpx(ji,jj) = ABS( (zsshp2_e(ji+1,jj) + bathy(ji+1,jj) - zsshp2_e(ji,jj) - bathy(ji,jj)) & 766 & / (zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj)) ) 767 ELSE 768 zcpx(ji,jj) = 0._wp 769 wduflt1(ji,jj) = 0._wp 770 END IF 771 772 ll_tmp1 = MIN( zsshp2_e(ji,jj), zsshp2_e(ji,jj+1) ) > MAX( -bathy(ji,jj), -bathy(ji,jj+1) ) & 773 & .AND. MAX( zsshp2_e(ji,jj) + bathy(ji,jj), zsshp2_e(ji,jj+1) + bathy(ji,jj+1) ) & 774 & > rn_wdmin1 + rn_wdmin2 775 ll_tmp2 = MAX( zsshp2_e(ji,jj), zsshp2_e(ji,jj+1) ) > MAX( -bathy(ji,jj), -bathy(ji,jj+1) ) & 776 & + rn_wdmin1 + rn_wdmin2 777 IF(ll_tmp1) THEN 778 zcpy(ji,jj) = 1._wp 779 ELSE IF(ll_tmp2) THEN 780 ! no worries about zsshp2_e(ji,jj+1)-zsshp2_e(ji,jj) = 0, it won't happen here 781 zcpy(ji,jj) = ABS( (zsshp2_e(ji,jj+1) + bathy(ji,jj+1) - zsshp2_e(ji,jj) - bathy(ji,jj)) & 782 & / (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj)) ) 783 ELSE 784 zcpy(ji,jj) = 0._wp 785 wdvflt1(ji,jj) = 0._wp 786 END IF 787 END DO 788 END DO 789 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 790 ENDIF 674 791 ! 675 792 ! Compute associated depths at U and V points: … … 688 805 END DO 689 806 END DO 807 808 IF( ln_wd ) THEN 809 zhust_e(:,:) = MAX(zhust_e (:,:), rn_wdmin1 ) 810 zhvst_e(:,:) = MAX(zhvst_e (:,:), rn_wdmin1 ) 811 END IF 812 690 813 ENDIF 691 814 ! … … 758 881 ! 759 882 ! Surface pressure trend: 760 DO jj = 2, jpjm1 761 DO ji = fs_2, fs_jpim1 ! vector opt. 762 ! Add surface pressure gradient 763 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj) 764 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj) 765 zwx(ji,jj) = zu_spg 766 zwy(ji,jj) = zv_spg 767 END DO 768 END DO 883 884 IF( ln_wd ) THEN 885 DO jj = 2, jpjm1 886 DO ji = 2, jpim1 887 ! Add surface pressure gradient 888 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj) 889 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj) 890 zwx(ji,jj) = zu_spg * zcpx(ji,jj) 891 zwy(ji,jj) = zv_spg * zcpy(ji,jj) 892 END DO 893 END DO 894 ELSE 895 DO jj = 2, jpjm1 896 DO ji = fs_2, fs_jpim1 ! vector opt. 897 ! Add surface pressure gradient 898 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj) 899 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj) 900 zwx(ji,jj) = zu_spg 901 zwy(ji,jj) = zv_spg 902 END DO 903 END DO 904 END IF 905 769 906 ! 770 907 ! Set next velocities: … … 789 926 DO jj = 2, jpjm1 790 927 DO ji = fs_2, fs_jpim1 ! vector opt. 791 zhura = ssumask(ji,jj)/(hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - ssumask(ji,jj)) 792 zhvra = ssvmask(ji,jj)/(hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - ssvmask(ji,jj)) 928 929 IF( ln_wd ) THEN 930 zhura = MAX(hu_0(ji,jj) + zsshu_a(ji,jj), rn_wdmin1) 931 zhvra = MAX(hv_0(ji,jj) + zsshv_a(ji,jj), rn_wdmin1) 932 ELSE 933 zhura = hu_0(ji,jj) + zsshu_a(ji,jj) 934 zhvra = hv_0(ji,jj) + zsshv_a(ji,jj) 935 END IF 936 zhura = ssumask(ji,jj)/(zhura + 1._wp - ssumask(ji,jj)) 937 zhvra = ssvmask(ji,jj)/(zhvra + 1._wp - ssvmask(ji,jj)) 793 938 794 939 ua_e(ji,jj) = ( hu_e(ji,jj) * un_e(ji,jj) & … … 808 953 ! 809 954 IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only) 810 hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 811 hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 955 IF( ln_wd ) THEN 956 hu_e (:,:) = MAX(hu_0(:,:) + zsshu_a(:,:), rn_wdmin1) 957 hv_e (:,:) = MAX(hv_0(:,:) + zsshv_a(:,:), rn_wdmin1) 958 ELSE 959 hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 960 hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 961 END IF 812 962 hur_e(:,:) = ssumask(:,:) / ( hu_e(:,:) + 1._wp - ssumask(:,:) ) 813 963 hvr_e(:,:) = ssvmask(:,:) / ( hv_e(:,:) + 1._wp - ssvmask(:,:) ) … … 936 1086 CALL wrk_dealloc( jpi,jpj, zsshu_a, zsshv_a ) 937 1087 CALL wrk_dealloc( jpi,jpj, zhf ) 1088 IF( ln_wd ) CALL wrk_dealloc( jpi, jpj, zcpx, zcpy, wduflt1, wdvflt1 ) 938 1089 ! 939 1090 IF ( ln_diatmb ) THEN
Note: See TracChangeset
for help on using the changeset viewer.
