Changeset 6152 for trunk/NEMOGCM/NEMO/OPA_SRC/DYN
- Timestamp:
- 2015-12-21T23:33:57+01:00 (8 years ago)
- Location:
- trunk/NEMOGCM/NEMO/OPA_SRC/DYN
- Files:
-
- 1 added
- 3 edited
-
dynhpg.F90 (modified) (19 diffs)
-
dynspg_ts.F90 (modified) (16 diffs)
-
sshwzv.F90 (modified) (2 diffs)
-
wet_dry.F90 (added)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
r6140 r6152 33 33 USE sbc_oce ! surface variable (only for the flag with ice shelf) 34 34 USE dom_oce ! ocean space and time domain 35 USE wet_dry ! wetting and drying 35 36 USE phycst ! physical constants 36 37 USE trd_oce ! trends: ocean variables … … 429 430 INTEGER, INTENT(in) :: kt ! ocean time-step index 430 431 !! 431 INTEGER :: ji, jj, jk ! dummy loop indices 432 REAL(wp) :: zcoef0, zuap, zvap, znad ! temporary scalars 432 INTEGER :: ji, jj, jk, jii, jjj ! dummy loop indices 433 REAL(wp) :: zcoef0, zuap, zvap, znad, ztmp ! temporary scalars 434 LOGICAL :: ll_tmp1, ll_tmp2, ll_tmp3 ! local logical variables 433 435 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj 436 REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy !W/D pressure filter 434 437 !!---------------------------------------------------------------------- 435 438 ! 436 439 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj ) 440 IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 437 441 ! 438 442 IF( kt == nit000 ) THEN … … 447 451 ENDIF 448 452 ! 453 IF(ln_wd) THEN 454 DO jj = 2, jpjm1 455 DO ji = 2, jpim1 456 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) 457 ll_tmp2 = MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) > rn_wdmin1 + rn_wdmin2 458 ll_tmp3 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) + & 459 & rn_wdmin1 + rn_wdmin2 460 461 IF(ll_tmp1.AND.ll_tmp2) THEN 462 zcpx(ji,jj) = 1.0_wp 463 wduflt(ji,jj) = 1.0_wp 464 ELSE IF(ll_tmp3) THEN 465 ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 466 zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) / & 467 & (sshn(ji+1,jj) - sshn(ji,jj))) 468 wduflt(ji,jj) = 1.0_wp 469 ELSE 470 zcpx(ji,jj) = 0._wp 471 wduflt(ji,jj) = 0.0_wp 472 END IF 473 474 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) 475 ll_tmp2 = MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) > rn_wdmin1 + rn_wdmin2 476 ll_tmp3 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) + & 477 & rn_wdmin1 + rn_wdmin2 478 479 IF(ll_tmp1.AND.ll_tmp2) THEN 480 zcpy(ji,jj) = 1.0_wp 481 wdvflt(ji,jj) = 1.0_wp 482 ELSE IF(ll_tmp3) THEN 483 ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 484 zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) / & 485 & (sshn(ji,jj+1) - sshn(ji,jj))) 486 wdvflt(ji,jj) = 1.0_wp 487 ELSE 488 zcpy(ji,jj) = 0._wp 489 wdvflt(ji,jj) = 0.0_wp 490 END IF 491 END DO 492 END DO 493 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 494 ENDIF 495 496 449 497 ! Surface value 450 498 DO jj = 2, jpjm1 … … 460 508 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) + 2._wp * znad ) & 461 509 & * ( gde3w_n(ji,jj+1,1) - gde3w_n(ji,jj,1) ) * r1_e2v(ji,jj) 510 511 512 IF(ln_wd) THEN 513 514 zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 515 zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 516 zuap = zuap * zcpx(ji,jj) 517 zvap = zvap * zcpy(ji,jj) 518 ENDIF 519 462 520 ! add to the general momentum trend 463 521 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap … … 482 540 zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 483 541 & * ( gde3w_n(ji ,jj+1,jk) - gde3w_n(ji,jj,jk) ) * r1_e2v(ji,jj) 542 543 IF(ln_wd) THEN 544 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 545 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) 546 zuap = zuap * zcpx(ji,jj) 547 zvap = zvap * zcpy(ji,jj) 548 ENDIF 549 484 550 ! add to the general momentum trend 485 551 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap … … 490 556 ! 491 557 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj ) 558 IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 492 559 ! 493 560 END SUBROUTINE hpg_sco … … 623 690 REAL(wp) :: z1_10, cffu, cffx ! " " 624 691 REAL(wp) :: z1_12, cffv, cffy ! " " 692 LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables 625 693 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj 626 694 REAL(wp), POINTER, DIMENSION(:,:,:) :: dzx, dzy, dzz, dzu, dzv, dzw 627 695 REAL(wp), POINTER, DIMENSION(:,:,:) :: drhox, drhoy, drhoz, drhou, drhov, drhow 628 696 REAL(wp), POINTER, DIMENSION(:,:,:) :: rho_i, rho_j, rho_k 697 REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy !W/D pressure filter 629 698 !!---------------------------------------------------------------------- 630 699 ! … … 632 701 CALL wrk_alloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow ) 633 702 CALL wrk_alloc( jpi, jpj, jpk, rho_i, rho_j, rho_k, zhpi, zhpj ) 634 ! 703 IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 704 ! 705 ! 706 IF(ln_wd) THEN 707 DO jj = 2, jpjm1 708 DO ji = 2, jpim1 709 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) & 710 & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) & 711 & > rn_wdmin1 + rn_wdmin2 712 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) +& 713 & rn_wdmin1 + rn_wdmin2 714 715 IF(ll_tmp1) THEN 716 zcpx(ji,jj) = 1.0_wp 717 ELSE IF(ll_tmp2) THEN 718 ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 719 zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /& 720 & (sshn(ji+1,jj) - sshn(ji,jj))) 721 ELSE 722 zcpx(ji,jj) = 0._wp 723 END IF 724 725 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) & 726 & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) & 727 & > rn_wdmin1 + rn_wdmin2 728 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) +& 729 & rn_wdmin1 + rn_wdmin2 730 731 IF(ll_tmp1) THEN 732 zcpy(ji,jj) = 1.0_wp 733 ELSE IF(ll_tmp2) THEN 734 ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 735 zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) /& 736 & (sshn(ji,jj+1) - sshn(ji,jj))) 737 ELSE 738 zcpy(ji,jj) = 0._wp 739 END IF 740 END DO 741 END DO 742 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 743 ENDIF 744 635 745 636 746 IF( kt == nit000 ) THEN … … 803 913 zhpi(ji,jj,1) = ( rho_k(ji+1,jj ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) * r1_e1u(ji,jj) 804 914 zhpj(ji,jj,1) = ( rho_k(ji ,jj+1,1) - rho_k(ji,jj,1) - rho_j(ji,jj,1) ) * r1_e2v(ji,jj) 915 IF(ln_wd) THEN 916 zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 917 zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 918 ENDIF 805 919 ! add to the general momentum trend 806 920 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) … … 822 936 & + ( ( rho_k(ji,jj+1,jk) - rho_k(ji,jj,jk ) ) & 823 937 & -( rho_j(ji,jj ,jk) - rho_j(ji,jj,jk-1) ) ) * r1_e2v(ji,jj) 938 IF(ln_wd) THEN 939 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 940 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) 941 ENDIF 824 942 ! add to the general momentum trend 825 943 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) … … 832 950 CALL wrk_dealloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow ) 833 951 CALL wrk_dealloc( jpi, jpj, jpk, rho_i, rho_j, rho_k, zhpi, zhpj ) 952 IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 834 953 ! 835 954 END SUBROUTINE hpg_djc … … 855 974 !! The local variables for the correction term 856 975 INTEGER :: jk1, jis, jid, jjs, jjd 976 LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables 857 977 REAL(wp) :: zuijk, zvijk, zpwes, zpwed, zpnss, zpnsd, zdeps 858 978 REAL(wp) :: zrhdt1 … … 861 981 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp 862 982 REAL(wp), POINTER, DIMENSION(:,:) :: zsshu_n, zsshv_n 983 REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy !W/D pressure filter 863 984 !!---------------------------------------------------------------------- 864 985 ! … … 866 987 CALL wrk_alloc( jpi,jpj,jpk, zdept, zrhh ) 867 988 CALL wrk_alloc( jpi,jpj, zsshu_n, zsshv_n ) 989 IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 868 990 ! 869 991 IF( kt == nit000 ) THEN … … 877 999 znad = 1._wp 878 1000 IF( ln_linssh ) znad = 0._wp 1001 1002 IF(ln_wd) THEN 1003 DO jj = 2, jpjm1 1004 DO ji = 2, jpim1 1005 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) & 1006 & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) & 1007 & > rn_wdmin1 + rn_wdmin2 1008 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) +& 1009 & rn_wdmin1 + rn_wdmin2 1010 1011 IF(ll_tmp1) THEN 1012 zcpx(ji,jj) = 1.0_wp 1013 ELSE IF(ll_tmp2) THEN 1014 ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 1015 zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /& 1016 & (sshn(ji+1,jj) - sshn(ji,jj))) 1017 ELSE 1018 zcpx(ji,jj) = 0._wp 1019 END IF 1020 1021 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) & 1022 & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) & 1023 & > rn_wdmin1 + rn_wdmin2 1024 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) +& 1025 & rn_wdmin1 + rn_wdmin2 1026 1027 IF(ll_tmp1.OR.ll_tmp2) THEN 1028 zcpy(ji,jj) = 1.0_wp 1029 ELSE IF(ll_tmp2) THEN 1030 ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 1031 zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) /& 1032 & (sshn(ji,jj+1) - sshn(ji,jj))) 1033 ELSE 1034 zcpy(ji,jj) = 0._wp 1035 END IF 1036 END DO 1037 END DO 1038 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 1039 ENDIF 879 1040 880 1041 ! Clean 3-D work arrays … … 960 1121 END DO 961 1122 END DO 1123 1124 CALL lbc_lnk (zsshu_n, 'U', 1.) 1125 CALL lbc_lnk (zsshv_n, 'V', 1.) 962 1126 963 1127 DO jj = 2, jpjm1 … … 1057 1221 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed) 1058 1222 ENDIF 1059 !!gm Since umask(ji,:,:) = tmask(ji,:,:) * tmask(ji+1,:,:) by definition 1060 !!gm in the line below only * umask(ji,jj,jk) is needed !! 1061 ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk) * tmask(ji,jj,jk) * tmask(ji+1,jj,jk) 1223 IF(ln_wd) THEN 1224 zdpdx1 = zdpdx1 * zcpx(ji,jj) 1225 zdpdx2 = zdpdx2 * zcpx(ji,jj) 1226 ENDIF 1227 ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk) 1062 1228 ENDIF 1063 1229 … … 1114 1280 zdpdy2 = zcoef0 * r1_e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd ) 1115 1281 ENDIF 1116 !!gm Since vmask(:,jj,:) = tmask(:,jj,:) * tmask(:,jj+1,:) by definition 1117 !!gm in the line below only * vmask(ji,jj,jk) is needed !! 1118 va(ji,jj,jk) = va(ji,jj,jk) + ( zdpdy1 + zdpdy2 ) * vmask(ji,jj,jk)*tmask(ji,jj,jk)*tmask(ji,jj+1,jk) 1282 IF(ln_wd) THEN 1283 zdpdy1 = zdpdy1 * zcpy(ji,jj) 1284 zdpdy2 = zdpdy2 * zcpy(ji,jj) 1285 ENDIF 1286 1287 va(ji,jj,jk) = va(ji,jj,jk) + (zdpdy1 + zdpdy2) * vmask(ji,jj,jk) 1119 1288 ENDIF 1120 1289 ! … … 1126 1295 CALL wrk_dealloc( jpi,jpj,jpk, zdept, zrhh ) 1127 1296 CALL wrk_dealloc( jpi,jpj, zsshu_n, zsshv_n ) 1297 IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 1128 1298 ! 1129 1299 END SUBROUTINE hpg_prj -
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 -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
r6140 r6152 39 39 USE wrk_nemo ! Memory Allocation 40 40 USE timing ! Timing 41 USE wet_dry ! Wetting/Drying flux limting 41 42 42 43 IMPLICIT NONE … … 104 105 ! 105 106 zcoef = 0.5_wp * r1_rau0 107 108 IF(ln_wd) CALL wad_lmt(sshb, zcoef * (emp_b(:,:) + emp(:,:)), z2dt) 109 106 110 ssha(:,:) = ( sshb(:,:) - z2dt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * ssmask(:,:) 107 111
Note: See TracChangeset
for help on using the changeset viewer.
