- Timestamp:
- 2015-02-06T17:02:20+01:00 (9 years ago)
- Location:
- branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN
- Files:
-
- 1 added
- 3 edited
-
dynhpg.F90 (modified) (20 diffs)
-
dynspg_ts.F90 (modified) (17 diffs)
-
sshwzv.F90 (modified) (4 diffs)
-
wadlmt.F90 (added)
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
r4794 r5066 16 16 !! 3.4 ! 2011-11 (H. Liu) hpg_prj: Original code for s-coordinates 17 17 !! ! (A. Coward) suppression of hel, wdj and rot options 18 !! 3.6? ! 2014-09 (H. Liu) add Wetting/Drying pressure filter 18 19 !!---------------------------------------------------------------------- 19 20 … … 369 370 !! 370 371 INTEGER :: ji, jj, jk ! dummy loop indices 371 REAL(wp) :: zcoef0, zuap, zvap, znad ! temporary scalars 372 REAL(wp) :: zcoef0, zuap, zvap, znad, ztmp ! temporary scalars 373 LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables 372 374 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj 375 REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy !W/D pressure filter 373 376 !!---------------------------------------------------------------------- 374 377 ! 375 378 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj ) 379 IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 376 380 ! 377 381 IF( kt == nit000 ) THEN … … 386 390 IF ( lk_vvl ) THEN ; znad = 1._wp ! Variable volume 387 391 ELSE ; znad = 0._wp ! Fixed volume 392 ENDIF 393 394 IF(ln_wd) THEN 395 DO jj = 2, jpjm1 396 DO ji = 2, jpim1 397 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) 398 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) +& 399 & rn_wdmin1 + rn_wdmin2 400 401 IF(ll_tmp1) THEN 402 zcpx(ji,jj) = 1.0_wp 403 ELSE IF(ll_tmp2) THEN 404 ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 405 zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /& 406 & (sshn(ji+1,jj) - sshn(ji,jj))) 407 ELSE 408 zcpx(ji,jj) = 0._wp 409 END IF 410 411 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) 412 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) +& 413 & rn_wdmin1 + rn_wdmin2 414 415 IF(ll_tmp1) THEN 416 zcpy(ji,jj) = 1.0_wp 417 ELSE IF(ll_tmp2) THEN 418 ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 419 zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) /& 420 & (sshn(ji,jj+1) - sshn(ji,jj))) 421 ELSE 422 zcpy(ji,jj) = 0._wp 423 END IF 424 END DO 425 END DO 426 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 388 427 ENDIF 389 428 … … 402 441 zvap = -zcoef0 * ( 2._wp * znad ) & 403 442 & * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj) 443 444 IF(ln_wd) THEN 445 zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 446 zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 447 zuap = zuap * zcpx(ji,jj) 448 zvap = zvap * zcpy(ji,jj) 449 ENDIF 450 404 451 ! add to the general momentum trend 405 452 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap … … 424 471 zvap = -zcoef0 * ( 2._wp * znad ) & 425 472 & * ( fsde3w(ji ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) 473 474 IF(ln_wd) THEN 475 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 476 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) 477 zuap = zuap * zcpx(ji,jj) 478 zvap = zvap * zcpy(ji,jj) 479 ENDIF 480 426 481 ! add to the general momentum trend 427 482 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap … … 445 500 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) + 2._wp * znad ) & 446 501 & * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj) 502 503 IF(ln_wd) THEN 504 zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 505 zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 506 zuap = zuap * zcpx(ji,jj) 507 zvap = zvap * zcpy(ji,jj) 508 ENDIF 509 447 510 ! add to the general momentum trend 448 511 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap … … 467 530 zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 468 531 & * ( fsde3w(ji ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) 532 533 IF(ln_wd) THEN 534 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 535 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) 536 zuap = zuap * zcpx(ji,jj) 537 zvap = zvap * zcpy(ji,jj) 538 ENDIF 539 469 540 ! add to the general momentum trend 470 541 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap … … 476 547 ! 477 548 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj ) 549 IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 478 550 ! 479 551 END SUBROUTINE hpg_sco … … 493 565 REAL(wp) :: z1_10, cffu, cffx ! " " 494 566 REAL(wp) :: z1_12, cffv, cffy ! " " 567 LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables 495 568 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj 496 569 REAL(wp), POINTER, DIMENSION(:,:,:) :: dzx, dzy, dzz, dzu, dzv, dzw 497 570 REAL(wp), POINTER, DIMENSION(:,:,:) :: drhox, drhoy, drhoz, drhou, drhov, drhow 498 571 REAL(wp), POINTER, DIMENSION(:,:,:) :: rho_i, rho_j, rho_k 572 REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy !W/D pressure filter 499 573 !!---------------------------------------------------------------------- 500 574 ! … … 502 576 CALL wrk_alloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow ) 503 577 CALL wrk_alloc( jpi, jpj, jpk, rho_i, rho_j, rho_k, zhpi, zhpj ) 578 IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 504 579 ! 505 580 … … 508 583 IF(lwp) WRITE(numout,*) 'dyn:hpg_djc : hydrostatic pressure gradient trend' 509 584 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, density Jacobian with cubic polynomial scheme' 585 ENDIF 586 587 IF(ln_wd) THEN 588 DO jj = 2, jpjm1 589 DO ji = 2, jpim1 590 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) 591 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) +& 592 & rn_wdmin1 + rn_wdmin2 593 594 IF(ll_tmp1) THEN 595 zcpx(ji,jj) = 1.0_wp 596 ELSE IF(ll_tmp2) THEN 597 ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 598 zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /& 599 & (sshn(ji+1,jj) - sshn(ji,jj))) 600 ELSE 601 zcpx(ji,jj) = 0._wp 602 END IF 603 604 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) 605 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) +& 606 & rn_wdmin1 + rn_wdmin2 607 608 IF(ll_tmp1) THEN 609 zcpy(ji,jj) = 1.0_wp 610 ELSE IF(ll_tmp2) THEN 611 ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 612 zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) /& 613 & (sshn(ji,jj+1) - sshn(ji,jj))) 614 ELSE 615 zcpy(ji,jj) = 0._wp 616 END IF 617 END DO 618 END DO 619 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 510 620 ENDIF 511 621 … … 673 783 zhpi(ji,jj,1) = ( rho_k(ji+1,jj ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) / e1u(ji,jj) 674 784 zhpj(ji,jj,1) = ( rho_k(ji ,jj+1,1) - rho_k(ji,jj,1) - rho_j(ji,jj,1) ) / e2v(ji,jj) 785 IF(ln_wd) THEN 786 zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 787 zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 788 ENDIF 675 789 ! add to the general momentum trend 676 790 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) … … 692 806 & + ( ( rho_k(ji,jj+1,jk) - rho_k(ji,jj,jk ) ) & 693 807 & -( rho_j(ji,jj ,jk) - rho_j(ji,jj,jk-1) ) ) / e2v(ji,jj) 808 IF(ln_wd) THEN 809 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 810 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) 811 ENDIF 694 812 ! add to the general momentum trend 695 813 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) … … 702 820 CALL wrk_dealloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow ) 703 821 CALL wrk_dealloc( jpi, jpj, jpk, rho_i, rho_j, rho_k, zhpi, zhpj ) 822 IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 704 823 ! 705 824 END SUBROUTINE hpg_djc … … 727 846 !! The local variables for the correction term 728 847 INTEGER :: jk1, jis, jid, jjs, jjd 848 LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables 729 849 REAL(wp) :: zuijk, zvijk, zpwes, zpwed, zpnss, zpnsd, zdeps 730 850 REAL(wp) :: zrhdt1 … … 732 852 REAL(wp), POINTER, DIMENSION(:,:,:) :: zdept, zrhh 733 853 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp 854 REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy !W/D pressure filter 734 855 !!---------------------------------------------------------------------- 735 856 ! 736 857 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp ) 737 858 CALL wrk_alloc( jpi,jpj,jpk, zdept, zrhh ) 859 IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 738 860 ! 739 861 IF( kt == nit000 ) THEN … … 748 870 znad = 0.0_wp 749 871 IF( lk_vvl ) znad = 1._wp 872 873 IF(ln_wd) THEN 874 DO jj = 2, jpjm1 875 DO ji = 2, jpim1 876 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) 877 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) +& 878 & rn_wdmin1 + rn_wdmin2 879 880 IF(ll_tmp1) THEN 881 zcpx(ji,jj) = 1.0_wp 882 ELSE IF(ll_tmp2) THEN 883 ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 884 zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /& 885 & (sshn(ji+1,jj) - sshn(ji,jj))) 886 ELSE 887 zcpx(ji,jj) = 0._wp 888 END IF 889 890 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) 891 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) +& 892 & rn_wdmin1 + rn_wdmin2 893 894 IF(ll_tmp1.OR.ll_tmp2) THEN 895 zcpy(ji,jj) = 1.0_wp 896 ELSE IF(ll_tmp2) THEN 897 ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 898 zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) /& 899 & (sshn(ji,jj+1) - sshn(ji,jj))) 900 ELSE 901 zcpy(ji,jj) = 0._wp 902 END IF 903 END DO 904 END DO 905 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 906 ENDIF 750 907 751 908 ! Clean 3-D work arrays … … 907 1064 ENDIF 908 1065 909 ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * & 910 & umask(ji,jj,jk) * tmask(ji,jj,jk) * tmask(ji+1,jj,jk) 1066 IF(ln_wd) THEN 1067 zdpdx1 = zdpdx1 * zcpx(ji,jj) 1068 zdpdx2 = zdpdx2 * zcpx(ji,jj) 1069 ENDIF 1070 ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk) 911 1071 ENDIF 912 1072 … … 964 1124 ENDIF 965 1125 966 va(ji,jj,jk) = va(ji,jj,jk) + (zdpdy1 + zdpdy2)*& 967 & vmask(ji,jj,jk)*tmask(ji,jj,jk)*tmask(ji,jj+1,jk) 1126 IF(ln_wd) THEN 1127 zdpdy1 = zdpdy1 * zcpy(ji,jj) 1128 zdpdy2 = zdpdy2 * zcpy(ji,jj) 1129 ENDIF 1130 1131 va(ji,jj,jk) = va(ji,jj,jk) + (zdpdy1 + zdpdy2) * vmask(ji,jj,jk) 968 1132 ENDIF 969 1133 … … 975 1139 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp ) 976 1140 CALL wrk_dealloc( jpi,jpj,jpk, zdept, zrhh ) 1141 IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 977 1142 ! 978 1143 END SUBROUTINE hpg_prj -
branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r4624 r5066 11 11 !! 3.5 ! 2013-07 (J. Chanut) Switch to Forward-backward time stepping 12 12 !! 3.6 ! 2013-11 (A. Coward) Update for z-tilde compatibility 13 !! 3.? ! 2014-09 (H. Liu) Add Wetting/Drying components 13 14 !!--------------------------------------------------------------------- 14 15 #if defined key_dynspg_ts || defined key_esopa … … 40 41 USE timing ! Timing 41 42 USE sbcapr ! surface boundary condition: atmospheric pressure 43 USE wadlmt ! wetting/drying flux limter 42 44 USE dynadv, ONLY: ln_dynadv_vec 43 45 #if defined key_agrif … … 137 139 INTEGER, INTENT(in) :: kt ! ocean time-step index 138 140 ! 139 LOGICAL :: ll_fw_start ! if true, forward integration 140 LOGICAL :: ll_init ! if true, special startup of 2d equations 141 INTEGER :: ji, jj, jk, jn ! dummy loop indices 142 INTEGER :: ikbu, ikbv, noffset ! local integers 143 REAL(wp) :: zraur, z1_2dt_b, z2dt_bf ! local scalars 144 REAL(wp) :: zx1, zy1, zx2, zy2 ! - - 145 REAL(wp) :: z1_12, z1_8, z1_4, z1_2 ! - - 146 REAL(wp) :: zu_spg, zv_spg ! - - 147 REAL(wp) :: zhura, zhvra ! - - 148 REAL(wp) :: za0, za1, za2, za3 ! - - 141 LOGICAL :: ll_fw_start ! if true, forward integration 142 LOGICAL :: ll_init ! if true, special startup of 2d equations 143 LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables used in W/D 144 INTEGER :: ji, jj, jk, jn ! dummy loop indices 145 INTEGER :: ikbu, ikbv, noffset ! local integers 146 REAL(wp) :: zraur, z1_2dt_b, z2dt_bf ! local scalars 147 REAL(wp) :: zx1, zy1, zx2, zy2 ! - - 148 REAL(wp) :: z1_12, z1_8, z1_4, z1_2 ! - - 149 REAL(wp) :: zu_spg, zv_spg ! - - 150 REAL(wp) :: zhura, zhvra ! - - 151 REAL(wp) :: za0, za1, za2, za3 ! - - 152 153 REAL(wp) :: ztmp ! temporary vaialbe 149 154 ! 150 155 REAL(wp), POINTER, DIMENSION(:,:) :: zun_e, zvn_e, zsshp2_e … … 154 159 REAL(wp), POINTER, DIMENSION(:,:) :: zsshu_a, zsshv_a 155 160 REAL(wp), POINTER, DIMENSION(:,:) :: zhf 161 REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy ! Wetting/Dying gravity filter coef. 156 162 !!---------------------------------------------------------------------- 157 163 ! … … 165 171 CALL wrk_alloc( jpi, jpj, zsshu_a, zsshv_a ) 166 172 CALL wrk_alloc( jpi, jpj, zhf ) 173 IF(ln_wd) CALL wrk_alloc( jpi, jpj, zcpx, zcpy ) 167 174 ! 168 175 ! !* Local constant initialization … … 173 180 zraur = 1._wp / rau0 174 181 ! 175 IF( kt == nit000 .AND. neuler == 0 ) THEN 182 IF( kt == nit000 .AND. neuler == 0 ) THEN ! reciprocal of baroclinic time step 176 183 z2dt_bf = rdt 177 184 ELSE … … 180 187 z1_2dt_b = 1.0_wp / z2dt_bf 181 188 ! 182 ll_init = ln_bt_av 189 ll_init = ln_bt_av ! if no time averaging, then no specific restart 183 190 ll_fw_start = .FALSE. 184 191 ! 185 192 ! time offset in steps for bdy data update 186 193 IF (.NOT.ln_bt_fw) THEN ; noffset=-2*nn_baro ; ELSE ; noffset = 0 ; ENDIF 187 194 ! 188 IF( kt == nit000 ) THEN 195 IF( kt == nit000 ) THEN !* initialisation 189 196 ! 190 197 IF(lwp) WRITE(numout,*) … … 365 372 ! ! ---------------------------------------------------- 366 373 IF( lk_vvl ) THEN ! Variable volume : remove surface pressure gradient 367 DO jj = 2, jpjm1 368 DO ji = fs_2, fs_jpim1 ! vector opt. 369 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) / e1u(ji,jj) 370 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) / e2v(ji,jj) 371 END DO 372 END DO 374 IF(ln_wd) THEN ! Calculating and applying W/D gravity filters 375 DO jj = 2, jpjm1 376 DO ji = 2, jpim1 377 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) 378 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) + & 379 & rn_wdmin1 + rn_wdmin2 380 IF(ll_tmp1) THEN 381 zcpx(ji,jj) = 1.0_wp 382 ELSE IF(ll_tmp2) THEN 383 ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 384 zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /& 385 & (sshn(ji+1,jj) - sshn(ji,jj))) 386 ELSE 387 zcpx(ji,jj) = 0._wp 388 END IF 389 390 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) 391 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) + & 392 & rn_wdmin1 + rn_wdmin2 393 IF(ll_tmp1) THEN 394 zcpy(ji,jj) = 1.0_wp 395 ELSE IF(ll_tmp2) THEN 396 ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 397 zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) + bathy(ji,jj)) /& 398 & (sshn(ji,jj+1) - sshn(ji,jj))) 399 ELSE 400 zcpy(ji,jj) = 0._wp 401 END IF 402 END DO 403 END DO 404 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 405 ENDIF 406 407 IF(ln_wd) THEN ! Calculating and applying W/D gravity filters 408 DO jj = 2, jpjm1 409 DO ji = 2, jpim1 410 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) / & 411 & e1u(ji,jj) * zcpx(ji,jj) 412 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) / & 413 & e2v(ji,jj) * zcpy(ji,jj) 414 END DO 415 END DO 416 ELSE 417 DO jj = 2, jpjm1 418 DO ji = fs_2, fs_jpim1 ! vector opt. 419 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) / e1u(ji,jj) 420 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) / e2v(ji,jj) 421 END DO 422 END DO 423 END IF 373 424 ENDIF 374 425 … … 464 515 ! ! ==================== ! 465 516 ! Initialize barotropic variables: 517 518 IF(ln_wd) THEN !preserve the positivity of water depth 519 !ssh[b,n,a] should have already been processed for this 520 sshbb_e(:,:) = MAX(sshbb_e(:,:), rn_wdmin1 - bathy(:,:)) 521 sshb_e(:,:) = MAX(sshb_e(:,:) , rn_wdmin1 - bathy(:,:)) 522 ENDIF 523 466 524 IF (ln_bt_fw) THEN ! FORWARD integration: start from NOW fields 467 525 sshn_e(:,:) = sshn (:,:) … … 537 595 zhup2_e (:,:) = hu_0(:,:) + zwx(:,:) ! Ocean depth at U- and V-points 538 596 zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:) 597 IF(ln_wd) THEN 598 zhup2_e(:,:) = MAX(zhup2_e (:,:), rn_wdmin1) 599 zhvp2_e(:,:) = MAX(zhvp2_e (:,:), rn_wdmin1) 600 END IF 539 601 ELSE 540 602 zhup2_e (:,:) = hu(:,:) … … 575 637 ENDIF 576 638 #endif 639 IF(ln_wd) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt) 577 640 ! 578 641 ! Sum over sub-time-steps to compute advective velocities … … 589 652 END DO 590 653 ssha_e(:,:) = ( sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) ) ) * tmask(:,:,1) 654 IF(ln_wd) ssha_e(:,:) = MAX(ssha_e(:,:), rn_wdmin1 - bathy(:,:)) 591 655 CALL lbc_lnk( ssha_e, 'T', 1._wp ) 592 656 … … 652 716 END DO 653 717 END DO 718 719 IF(ln_wd) THEN 720 zhust_e(:,:) = MAX(zhust_e (:,:), rn_wdmin1 ) 721 zhvst_e(:,:) = MAX(zhvst_e (:,:), rn_wdmin1 ) 722 END IF 723 !shall we call lbc_lnk for zhu(v)st_e() here? 724 654 725 ENDIF 655 726 ! … … 717 788 zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * zvn_e(:,:) * hvr_e(:,:) 718 789 ! 719 ! Surface pressure trend: 720 DO jj = 2, jpjm1 721 DO ji = fs_2, fs_jpim1 ! vector opt. 722 ! Add surface pressure gradient 723 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj) 724 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj) 725 zwx(ji,jj) = zu_spg 726 zwy(ji,jj) = zv_spg 727 END DO 728 END DO 790 IF(ln_wd) THEN ! Calculating and applying W/D gravity filters 791 DO jj = 2, jpjm1 792 DO ji = 2, jpim1 793 ll_tmp1 = MIN(zsshp2_e(ji,jj), zsshp2_e(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) 794 ll_tmp2 = MAX(zsshp2_e(ji,jj), zsshp2_e(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) + & 795 & rn_wdmin1 + rn_wdmin2 796 IF(ll_tmp1) THEN 797 zcpx(ji,jj) = 1.0_wp 798 ELSE IF(ll_tmp2) THEN 799 ! no worries about zsshp2_e(ji+1,jj)-zsshp2_e(ji,jj) = 0, it won't happen ! here 800 zcpx(ji,jj) = ABS((zsshp2_e(ji+1,jj) + bathy(ji+1,jj) - zsshp2_e(ji,jj) - bathy(ji,jj)) /& 801 & (zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj))) 802 ELSE 803 zcpx(ji,jj) = 0._wp 804 END IF 805 806 ll_tmp1 = MIN(zsshp2_e(ji,jj), zsshp2_e(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) 807 ll_tmp2 = MAX(zsshp2_e(ji,jj), zsshp2_e(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) + & 808 & rn_wdmin1 + rn_wdmin2 809 IF(ll_tmp1) THEN 810 zcpy(ji,jj) = 1.0_wp 811 ELSE IF(ll_tmp2) THEN 812 ! no worries about zsshp2_e(ji,jj+1)-zsshp2_e(ji,jj) = 0, it won't happen ! here 813 zcpy(ji,jj) = ABS((zsshp2_e(ji,jj+1) + bathy(ji,jj+1) - zsshp2_e(ji,jj) - bathy(ji,jj)) /& 814 & (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj))) 815 ELSE 816 zcpy(ji,jj) = 0._wp 817 END IF 818 END DO 819 END DO 820 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 821 ENDIF 822 823 IF(ln_wd) THEN 824 DO jj = 2, jpjm1 825 DO ji = 2, jpim1 826 ! Add surface pressure gradient 827 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj) 828 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj) 829 zwx(ji,jj) = zu_spg * zcpx(ji,jj) 830 zwy(ji,jj) = zv_spg * zcpy(ji,jj) 831 END DO 832 END DO 833 ELSE 834 DO jj = 2, jpjm1 835 DO ji = fs_2, fs_jpim1 ! vector opt. 836 ! Add surface pressure gradient 837 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj) 838 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj) 839 zwx(ji,jj) = zu_spg 840 zwy(ji,jj) = zv_spg 841 END DO 842 END DO 843 END IF 729 844 ! 730 845 ! Set next velocities: … … 750 865 DO ji = fs_2, fs_jpim1 ! vector opt. 751 866 752 zhura = umask(ji,jj,1)/(hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - umask(ji,jj,1)) 753 zhvra = vmask(ji,jj,1)/(hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - vmask(ji,jj,1)) 754 755 ua_e(ji,jj) = ( hu_e(ji,jj) * zun_e(ji,jj) & 756 & + rdtbt * ( zhust_e(ji,jj) * zwx(ji,jj) & 867 IF(ln_wd) THEN 868 zhura = MAX(hu_0(ji,jj) + zsshu_a(ji,jj), rn_wdmin1) 869 zhvra = MAX(hv_0(ji,jj) + zsshv_a(ji,jj), rn_wdmin1) 870 ELSE 871 zhura = hu_0(ji,jj) + zsshu_a(ji,jj) 872 zhvra = hv_0(ji,jj) + zsshv_a(ji,jj) 873 END IF 874 875 zhura = umask(ji,jj,1)/(zhura + 1._wp - umask(ji,jj,1)) 876 zhvra = vmask(ji,jj,1)/(zhvra + 1._wp - vmask(ji,jj,1)) 877 878 ua_e(ji,jj) = ( hu_e(ji,jj) * zun_e(ji,jj) & 879 & + rdtbt * ( zhust_e(ji,jj) * zwx(ji,jj) & 757 880 & + zhup2_e(ji,jj) * zu_trd(ji,jj) & 758 881 & + hu(ji,jj) * zu_frc(ji,jj) ) & … … 770 893 IF( lk_vvl ) THEN !* Update ocean depth (variable volume case only) 771 894 ! ! ---------------------------------------------- 772 hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 773 hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 895 IF(ln_wd) THEN 896 hu_e (:,:) = MAX(hu_0(:,:) + zsshu_a(:,:), rn_wdmin1) 897 hv_e (:,:) = MAX(hv_0(:,:) + zsshv_a(:,:), rn_wdmin1) 898 ELSE 899 hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 900 hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 901 END IF 902 774 903 hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1._wp - umask(:,:,1) ) 775 904 hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1._wp - vmask(:,:,1) ) … … 903 1032 CALL wrk_dealloc( jpi, jpj, zsshu_a, zsshv_a ) 904 1033 CALL wrk_dealloc( jpi, jpj, zhf ) 1034 IF(ln_wd) CALL wrk_alloc( jpi, jpj, zcpx, zcpy ) 905 1035 ! 906 1036 IF( nn_timing == 1 ) CALL timing_stop('dyn_spg_ts') -
branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
r4486 r5066 9 9 !! - ! 2010-09 (D.Storkey and E.O'Dea) bug fixes for BDY module 10 10 !! 3.3 ! 2011-10 (M. Leclair) split former ssh_wzv routine and remove all vvl related work 11 !! 3.? ! 2014-09 (H. Liu) add wetting and drying 11 12 !!---------------------------------------------------------------------- 12 13 … … 42 43 USE wrk_nemo ! Memory Allocation 43 44 USE timing ! Timing 45 USE wadlmt ! Wetting/Drying flux limting 44 46 45 47 IMPLICIT NONE … … 94 96 ENDIF 95 97 ! 96 CALL div_cur( kt ) ! Horizontal divergence & Relative vorticity97 !98 98 z2dt = 2._wp * rdt ! set time step size (Euler/Leapfrog) 99 99 IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt 100 100 ! 101 z1_rau0 = 0.5_wp * r1_rau0 102 ! 103 IF(ln_wd) CALL wad_lmt(sshb, z1_rau0 * (emp_b(:,:) + emp(:,:)), z2dt) 104 105 CALL div_cur( kt ) ! Horizontal divergence & Relative vorticity 106 ! 101 107 ! !------------------------------! 102 108 ! ! After Sea Surface Height ! … … 110 116 ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations. 111 117 ! 112 z1_rau0 = 0.5_wp * r1_rau0113 118 ssha(:,:) = ( sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * tmask(:,:,1) 114 119
Note: See TracChangeset
for help on using the changeset viewer.
