Changeset 7412 for branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/DYN
- Timestamp:
- 2016-12-01T11:30:29+01:00 (8 years ago)
- Location:
- branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/DYN
- Files:
-
- 7 edited
-
dynhpg.F90 (modified) (17 diffs)
-
dynkeg.F90 (modified) (4 diffs)
-
dynnxt.F90 (modified) (3 diffs)
-
dynspg.F90 (modified) (2 diffs)
-
dynspg_ts.F90 (modified) (17 diffs)
-
sshwzv.F90 (modified) (5 diffs)
-
wet_dry.F90 (modified) (18 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
r6152 r7412 432 432 INTEGER :: ji, jj, jk, jii, jjj ! dummy loop indices 433 433 REAL(wp) :: zcoef0, zuap, zvap, znad, ztmp ! temporary scalars 434 LOGICAL :: ll_tmp1, ll_tmp2 , ll_tmp3! local logical variables434 LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables 435 435 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj 436 436 REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy !W/D pressure filter … … 438 438 ! 439 439 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj ) 440 IF( ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy )440 IF( ln_wd ) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 441 441 ! 442 442 IF( kt == nit000 ) THEN … … 451 451 ENDIF 452 452 ! 453 IF( ln_wd) THEN453 IF( ln_wd ) THEN 454 454 DO jj = 2, jpjm1 455 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 456 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > & 457 & MAX( -bathy(ji,jj) , -bathy(ji+1,jj) ) .AND. & 458 & MAX( sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj) ) & 459 & > rn_wdmin1 + rn_wdmin2 460 ll_tmp2 = MAX( sshn(ji,jj) , sshn(ji+1,jj) ) > & 461 & MAX( -bathy(ji,jj) , -bathy(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 462 463 IF(ll_tmp1) THEN 462 464 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 465 ELSE IF(ll_tmp2) THEN 466 ! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here 467 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) & 468 & / (sshn(ji+1,jj) - sshn(ji ,jj)) ) 469 469 ELSE 470 470 zcpx(ji,jj) = 0._wp 471 wduflt(ji,jj) = 0.0_wp472 471 END IF 473 472 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 473 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > & 474 & MAX( -bathy(ji,jj) , -bathy(ji,jj+1) ) .AND. & 475 & MAX( sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1) ) & 476 & > rn_wdmin1 + rn_wdmin2 477 ll_tmp2 = MAX( sshn(ji,jj) , sshn(ji,jj+1) ) > & 478 & MAX( -bathy(ji,jj) , -bathy(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 479 480 IF(ll_tmp1) THEN 480 481 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 482 ELSE IF(ll_tmp2) 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 )) ) 487 486 ELSE 488 487 zcpy(ji,jj) = 0._wp 489 wdvflt(ji,jj) = 0.0_wp490 488 END IF 491 489 END DO 492 490 END DO 493 491 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 494 ENDIF 495 492 END IF 496 493 497 494 ! Surface value … … 510 507 511 508 512 IF( ln_wd) THEN509 IF( ln_wd ) THEN 513 510 514 511 zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) … … 541 538 & * ( gde3w_n(ji ,jj+1,jk) - gde3w_n(ji,jj,jk) ) * r1_e2v(ji,jj) 542 539 543 IF( ln_wd) THEN540 IF( ln_wd ) THEN 544 541 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 545 542 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) … … 556 553 ! 557 554 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj ) 558 IF( ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy )555 IF( ln_wd ) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 559 556 ! 560 557 END SUBROUTINE hpg_sco … … 701 698 CALL wrk_alloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow ) 702 699 CALL wrk_alloc( jpi, jpj, jpk, rho_i, rho_j, rho_k, zhpi, zhpj ) 703 IF( ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy )704 ! 705 ! 706 IF( ln_wd) THEN700 IF( ln_wd ) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 701 ! 702 ! 703 IF( ln_wd ) THEN 707 704 DO jj = 2, jpjm1 708 705 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 706 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > & 707 & MAX( -bathy(ji,jj) , -bathy(ji+1,jj) ) .AND. & 708 & MAX( sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj) ) & 709 & > rn_wdmin1 + rn_wdmin2 710 ll_tmp2 = MAX( sshn(ji,jj) , sshn(ji+1,jj) ) > & 711 & MAX( -bathy(ji,jj) , -bathy(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 714 712 715 713 IF(ll_tmp1) THEN 716 714 zcpx(ji,jj) = 1.0_wp 717 715 ELSE IF(ll_tmp2) THEN 718 ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here719 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)))716 ! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here 717 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) & 718 & / (sshn(ji+1,jj) - sshn(ji ,jj)) ) 721 719 ELSE 722 720 zcpx(ji,jj) = 0._wp 723 721 END IF 724 722 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 723 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > & 724 & MAX( -bathy(ji,jj) , -bathy(ji,jj+1) ) .AND. & 725 & MAX( sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1) ) & 726 & > rn_wdmin1 + rn_wdmin2 727 ll_tmp2 = MAX( sshn(ji,jj) , sshn(ji,jj+1) ) > & 728 & MAX( -bathy(ji,jj) , -bathy(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 730 729 731 730 IF(ll_tmp1) THEN 732 731 zcpy(ji,jj) = 1.0_wp 733 732 ELSE IF(ll_tmp2) THEN 734 ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here735 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)))733 ! no worries about sshn(ji,jj+1) - sshn(ji,jj ) = 0, it won't happen ! here 734 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) & 735 & / (sshn(ji,jj+1) - sshn(ji,jj )) ) 737 736 ELSE 738 737 zcpy(ji,jj) = 0._wp … … 741 740 END DO 742 741 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 743 ENDIF 744 742 END IF 745 743 746 744 IF( kt == nit000 ) THEN … … 913 911 zhpi(ji,jj,1) = ( rho_k(ji+1,jj ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) * r1_e1u(ji,jj) 914 912 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) THEN913 IF( ln_wd ) THEN 916 914 zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 917 915 zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) … … 936 934 & + ( ( rho_k(ji,jj+1,jk) - rho_k(ji,jj,jk ) ) & 937 935 & -( rho_j(ji,jj ,jk) - rho_j(ji,jj,jk-1) ) ) * r1_e2v(ji,jj) 938 IF( ln_wd) THEN936 IF( ln_wd ) THEN 939 937 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 940 938 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) … … 950 948 CALL wrk_dealloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow ) 951 949 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 )950 IF( ln_wd ) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 953 951 ! 954 952 END SUBROUTINE hpg_djc … … 987 985 CALL wrk_alloc( jpi,jpj,jpk, zdept, zrhh ) 988 986 CALL wrk_alloc( jpi,jpj, zsshu_n, zsshv_n ) 989 IF( ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy )987 IF( ln_wd ) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 990 988 ! 991 989 IF( kt == nit000 ) THEN … … 1000 998 IF( ln_linssh ) znad = 0._wp 1001 999 1002 IF( ln_wd) THEN1000 IF( ln_wd ) THEN 1003 1001 DO jj = 2, jpjm1 1004 1002 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 1003 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > & 1004 & MAX( -bathy(ji,jj) , -bathy(ji+1,jj) ) .AND. & 1005 & MAX( sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj) ) & 1006 & > rn_wdmin1 + rn_wdmin2 1007 ll_tmp2 = MAX( sshn(ji,jj) , sshn(ji+1,jj) ) > & 1008 & MAX( -bathy(ji,jj) , -bathy(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 1010 1009 1011 1010 IF(ll_tmp1) THEN 1012 1011 zcpx(ji,jj) = 1.0_wp 1013 1012 ELSE IF(ll_tmp2) THEN 1014 ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here1015 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)))1013 ! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here 1014 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) & 1015 & / (sshn(ji+1,jj) - sshn(ji ,jj)) ) 1017 1016 ELSE 1018 1017 zcpx(ji,jj) = 0._wp 1019 1018 END IF 1020 1019 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 1020 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > & 1021 & MAX( -bathy(ji,jj) , -bathy(ji,jj+1) ) .AND. & 1022 & 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) ) > & 1025 & MAX( -bathy(ji,jj) , -bathy(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 1026 1027 IF(ll_tmp1) THEN 1028 1028 zcpy(ji,jj) = 1.0_wp 1029 1029 ELSE IF(ll_tmp2) THEN 1030 ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here1031 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)))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 1033 ELSE 1034 1034 zcpy(ji,jj) = 0._wp … … 1037 1037 END DO 1038 1038 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 1039 END IF1039 END IF 1040 1040 1041 1041 ! Clean 3-D work arrays … … 1221 1221 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed) 1222 1222 ENDIF 1223 IF( ln_wd) THEN1223 IF( ln_wd ) THEN 1224 1224 zdpdx1 = zdpdx1 * zcpx(ji,jj) 1225 1225 zdpdx2 = zdpdx2 * zcpx(ji,jj) … … 1280 1280 zdpdy2 = zcoef0 * r1_e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd ) 1281 1281 ENDIF 1282 IF( ln_wd) THEN1282 IF( ln_wd ) THEN 1283 1283 zdpdy1 = zdpdy1 * zcpy(ji,jj) 1284 1284 zdpdy2 = zdpdy2 * zcpy(ji,jj) … … 1295 1295 CALL wrk_dealloc( jpi,jpj,jpk, zdept, zrhh ) 1296 1296 CALL wrk_dealloc( jpi,jpj, zsshu_n, zsshv_n ) 1297 IF( ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy )1297 IF( ln_wd ) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 1298 1298 ! 1299 1299 END SUBROUTINE hpg_prj -
branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90
r5328 r7412 24 24 USE wrk_nemo ! Memory Allocation 25 25 USE timing ! Timing 26 USE bdy_oce ! ocean open boundary conditions 26 27 27 28 IMPLICIT NONE … … 78 79 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhke 79 80 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdu, ztrdv 81 INTEGER :: jb ! dummy loop indices 82 INTEGER :: ii, ij, igrd, ib_bdy ! local integers 83 INTEGER :: fu, fv 80 84 !!---------------------------------------------------------------------- 81 85 ! … … 98 102 zhke(:,:,jpk) = 0._wp 99 103 104 IF (ln_bdy) THEN 105 ! Maria Luneva & Fred Wobus: July-2016 106 ! compensate for lack of turbulent kinetic energy on liquid bdy points 107 DO ib_bdy = 1, nb_bdy 108 IF( cn_dyn3d(ib_bdy) /= 'none' ) THEN 109 igrd = 2 ! Copying normal velocity into points outside bdy 110 DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 111 DO jk = 1, jpkm1 112 ii = idx_bdy(ib_bdy)%nbi(jb,igrd) 113 ij = idx_bdy(ib_bdy)%nbj(jb,igrd) 114 fu = NINT( idx_bdy(ib_bdy)%flagu(jb,igrd) ) 115 un(ii-fu,ij,jk) = un(ii,ij,jk) * umask(ii,ij,jk) 116 END DO 117 END DO 118 ! 119 igrd = 3 ! Copying normal velocity into points outside bdy 120 DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 121 DO jk = 1, jpkm1 122 ii = idx_bdy(ib_bdy)%nbi(jb,igrd) 123 ij = idx_bdy(ib_bdy)%nbj(jb,igrd) 124 fv = NINT( idx_bdy(ib_bdy)%flagv(jb,igrd) ) 125 vn(ii,ij-fv,jk) = vn(ii,ij,jk) * vmask(ii,ij,jk) 126 END DO 127 END DO 128 ENDIF 129 ENDDO 130 ENDIF 131 100 132 SELECT CASE ( kscheme ) !== Horizontal kinetic energy at T-point ==! 101 133 ! … … 133 165 ! 134 166 END SELECT 167 168 IF (ln_bdy) THEN 169 ! restore velocity masks at points outside boundary 170 un(:,:,:) = un(:,:,:) * umask(:,:,:) 171 vn(:,:,:) = vn(:,:,:) * vmask(:,:,:) 172 ENDIF 173 174 135 175 ! 136 176 DO jk = 1, jpkm1 !== grad( KE ) added to the general momentum trends ==! -
branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90
r6140 r7412 32 32 USE dynspg_ts ! surface pressure gradient: split-explicit scheme 33 33 USE domvvl ! variable volume 34 USE bdy_oce ! ocean open boundary conditions34 USE bdy_oce , ONLY: ln_bdy 35 35 USE bdydta ! ocean open boundary conditions 36 36 USE bdydyn ! ocean open boundary conditions … … 77 77 !! * Apply lateral boundary conditions on after velocity 78 78 !! at the local domain boundaries through lbc_lnk call, 79 !! at the one-way open boundaries (l k_bdy=T),79 !! at the one-way open boundaries (ln_bdy=T), 80 80 !! at the AGRIF zoom boundaries (lk_agrif=T) 81 81 !! … … 147 147 CALL lbc_lnk( va, 'V', -1. ) 148 148 ! 149 # if defined key_bdy150 149 ! !* BDY open boundaries 151 IF( l k_bdy .AND. ln_dynspg_exp ) CALL bdy_dyn( kt )152 IF( l k_bdy .AND. ln_dynspg_ts ) CALL bdy_dyn( kt, dyn3d_only=.true. )150 IF( ln_bdy .AND. ln_dynspg_exp ) CALL bdy_dyn( kt ) 151 IF( ln_bdy .AND. ln_dynspg_ts ) CALL bdy_dyn( kt, dyn3d_only=.true. ) 153 152 154 153 !!$ Do we need a call to bdy_vol here?? 155 !156 # endif157 154 ! 158 155 IF( l_trddyn ) THEN ! prepare the atf trend computation + some diagnostics -
branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90
r6981 r7412 88 88 ! 89 89 IF( ln_apr_dyn & ! atmos. pressure 90 .OR. ( .NOT.ln_dynspg_ts .AND. (ln_tide_pot .AND. l k_tide) ) & ! tide potential (no time slitting)90 .OR. ( .NOT.ln_dynspg_ts .AND. (ln_tide_pot .AND. ln_tide) ) & ! tide potential (no time slitting) 91 91 .OR. nn_ice_embd == 2 ) THEN ! embedded sea-ice 92 92 ! … … 111 111 ! 112 112 ! !== tide potential forcing term ==! 113 IF( .NOT.ln_dynspg_ts .AND. ( ln_tide_pot .AND. l k_tide ) ) THEN ! N.B. added directly at sub-time-step in ts-case113 IF( .NOT.ln_dynspg_ts .AND. ( ln_tide_pot .AND. ln_tide ) ) THEN ! N.B. added directly at sub-time-step in ts-case 114 114 ! 115 115 CALL upd_tide( kt ) ! update tide potential -
branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r6152 r7412 33 33 USE dynvor ! vorticity term 34 34 USE wet_dry ! wetting/drying flux limter 35 USE bdy_ par ! for lk_bdy35 USE bdy_oce , ONLY: ln_bdy 36 36 USE bdytides ! open boundary condition data 37 37 USE bdydyn2d ! open boundary conditions on barotropic variables … … 156 156 REAL(wp), POINTER, DIMENSION(:,:) :: zhf 157 157 REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy ! Wetting/Dying gravity filter coef. 158 REAL(wp), POINTER, DIMENSION(:,:) :: wduflt1, wdvflt1 ! Wetting/Dying velocity filter coef.159 158 !!---------------------------------------------------------------------- 160 159 ! … … 168 167 CALL wrk_alloc( jpi,jpj, zsshu_a, zsshv_a ) 169 168 CALL wrk_alloc( jpi,jpj, zhf ) 170 IF( ln_wd ) CALL wrk_alloc( jpi, jpj, zcpx, zcpy , wduflt1, wdvflt1)169 IF( ln_wd ) CALL wrk_alloc( jpi, jpj, zcpx, zcpy ) 171 170 ! 172 171 zmdi=1.e+20 ! missing data indicator for masking … … 374 373 IF( .NOT.ln_linssh ) THEN ! Variable volume : remove surface pressure gradient 375 374 IF( ln_wd ) THEN ! Calculating and applying W/D gravity filters 376 wduflt1(:,:) = 1.0_wp377 wdvflt1(:,:) = 1.0_wp378 DO jj = 2, jpjm1379 DO ji = 2, jpim1380 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_wdmin2383 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) &384 & + rn_wdmin1 + rn_wdmin2375 DO jj = 2, jpjm1 376 DO ji = 2, jpim1 377 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > & 378 & MAX( -bathy(ji,jj) , -bathy(ji+1,jj) ) .AND. & 379 & MAX( sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj) ) & 380 & > rn_wdmin1 + rn_wdmin2 381 ll_tmp2 = MAX( sshn(ji,jj) , sshn(ji+1,jj) ) > & 382 & MAX( -bathy(ji,jj) , -bathy(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 383 385 384 IF(ll_tmp1) THEN 386 zcpx(ji,jj) 387 ELSE IF(ll_tmp2) THEN388 ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happenhere389 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)))385 zcpx(ji,jj) = 1.0_wp 386 ELSE IF(ll_tmp2) THEN 387 ! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here 388 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) & 389 & / (sshn(ji+1,jj) - sshn(ji ,jj)) ) 391 390 ELSE 392 zcpx(ji,jj) = 0._wp 393 wduflt1(ji,jj) = 0.0_wp 391 zcpx(ji,jj) = 0._wp 394 392 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 393 394 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > & 395 & MAX( -bathy(ji,jj) , -bathy(ji,jj+1) ) .AND. & 396 & MAX( sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1) ) & 397 & > rn_wdmin1 + rn_wdmin2 398 ll_tmp2 = MAX( sshn(ji,jj) , sshn(ji,jj+1) ) > & 399 & MAX( -bathy(ji,jj) , -bathy(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 400 401 401 IF(ll_tmp1) THEN 402 zcpy(ji,jj)= 1.0_wp403 ELSE IF(ll_tmp2) THEN404 ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happenhere405 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)))402 zcpy(ji,jj) = 1.0_wp 403 ELSE IF(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 407 ELSE 408 zcpy(ji,jj) = 0._wp 409 wdvflt1(ji,jj) = 0.0_wp 410 ENDIF 411 412 END DO 408 zcpy(ji,jj) = 0._wp 409 END IF 410 END DO 413 411 END DO 414 415 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 416 412 417 413 DO jj = 2, jpjm1 418 414 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)415 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) & 416 & * r1_e1u(ji,jj) * zcpx(ji,jj) 417 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) & 418 & * r1_e2v(ji,jj) * zcpy(ji,jj) 423 419 END DO 424 420 END DO … … 567 563 ENDIF 568 564 569 IF( ln_wd ) THEN !preserve the positivity of water depth570 !ssh[b,n,a] should have already been processed for this571 sshbb_e(:,:) = MAX(sshbb_e(:,:), rn_wdmin1 - bathy(:,:))572 sshb_e(:,:) = MAX(sshb_e(:,:) , rn_wdmin1 - bathy(:,:))573 ENDIF574 565 ! 575 566 IF (ln_bt_fw) THEN ! FORWARD integration: start from NOW fields … … 607 598 ! ! ------------------ 608 599 ! Update only tidal forcing at open boundaries 609 #if defined key_tide 610 IF( lk_bdy .AND. lk_tide ) CALL bdy_dta_tides( kt, kit=jn, time_offset= noffset+1 ) 611 IF( ln_tide_pot .AND. lk_tide ) CALL upd_tide ( kt, kit=jn, time_offset= noffset ) 612 #endif 600 IF( ln_bdy .AND. ln_tide ) CALL bdy_dta_tides( kt, kit=jn, time_offset= noffset+1 ) 601 IF( ln_tide_pot .AND. ln_tide ) CALL upd_tide ( kt, kit=jn, time_offset= noffset ) 613 602 ! 614 603 ! Set extrapolation coefficients for predictor step: … … 646 635 zhup2_e (:,:) = hu_0(:,:) + zwx(:,:) ! Ocean depth at U- and V-points 647 636 zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:) 648 IF( ln_wd ) THEN649 zhup2_e(:,:) = MAX(zhup2_e (:,:), rn_wdmin1)650 zhvp2_e(:,:) = MAX(zhvp2_e (:,:), rn_wdmin1)651 END IF652 637 ELSE 653 638 zhup2_e (:,:) = hu_n(:,:) … … 701 686 END DO 702 687 END DO 688 703 689 ssha_e(:,:) = ( sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) ) ) * ssmask(:,:) 704 IF( ln_wd ) ssha_e(:,:) = MAX(ssha_e(:,:), rn_wdmin1 - bathy(:,:))690 705 691 CALL lbc_lnk( ssha_e, 'T', 1._wp ) 706 692 707 #if defined key_bdy708 693 ! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T) 709 IF( l k_bdy ) CALL bdy_ssh( ssha_e )710 #endif 694 IF( ln_bdy ) CALL bdy_ssh( ssha_e ) 695 711 696 #if defined key_agrif 712 697 IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jn ) … … 749 734 zsshp2_e(:,:) = za0 * ssha_e(:,:) + za1 * sshn_e (:,:) & 750 735 & + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:) 736 751 737 IF( ln_wd ) THEN ! Calculating and applying W/D gravity filters 752 wduflt1(:,:) = 1._wp753 wdvflt1(:,:) = 1._wp754 738 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 739 DO ji = 2, jpim1 740 ll_tmp1 = MIN( zsshp2_e(ji,jj) , zsshp2_e(ji+1,jj) ) > & 741 & MAX( -bathy(ji,jj) , -bathy(ji+1,jj) ) .AND. & 742 & MAX( zsshp2_e(ji,jj) + bathy(ji,jj), zsshp2_e(ji+1,jj) + bathy(ji+1,jj) ) & 743 & > rn_wdmin1 + rn_wdmin2 744 ll_tmp2 = MAX( zsshp2_e(ji,jj) , zsshp2_e(ji+1,jj) ) > & 745 & MAX( -bathy(ji,jj) , -bathy(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 746 747 IF(ll_tmp1) THEN 748 zcpx(ji,jj) = 1.0_wp 749 ELSE IF(ll_tmp2) THEN 750 ! no worries about zsshp2_e(ji+1,jj) - zsshp2_e(ji ,jj) = 0, it won't happen ! here 751 zcpx(ji,jj) = ABS( (zsshp2_e(ji+1,jj) + bathy(ji+1,jj) - zsshp2_e(ji,jj) - bathy(ji,jj)) & 752 & / (zsshp2_e(ji+1,jj) - zsshp2_e(ji ,jj)) ) 753 ELSE 754 zcpx(ji,jj) = 0._wp 755 END IF 756 757 ll_tmp1 = MIN( zsshp2_e(ji,jj) , zsshp2_e(ji,jj+1) ) > & 758 & MAX( -bathy(ji,jj) , -bathy(ji,jj+1) ) .AND. & 759 & MAX( zsshp2_e(ji,jj) + bathy(ji,jj), zsshp2_e(ji,jj+1) + bathy(ji,jj+1) ) & 760 & > rn_wdmin1 + rn_wdmin2 761 ll_tmp2 = MAX( zsshp2_e(ji,jj) , zsshp2_e(ji,jj+1) ) > & 762 & MAX( -bathy(ji,jj) , -bathy(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 763 764 IF(ll_tmp1) THEN 765 zcpy(ji,jj) = 1.0_wp 766 ELSE IF(ll_tmp2) THEN 767 ! no worries about zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj ) = 0, it won't happen ! here 768 zcpy(ji,jj) = ABS( (zsshp2_e(ji,jj+1) + bathy(ji,jj+1) - zsshp2_e(ji,jj) - bathy(ji,jj)) & 769 & / (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj )) ) 770 ELSE 771 zcpy(ji,jj) = 0._wp 772 END IF 787 773 END DO 788 END DO 789 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 790 ENDIF 774 END DO 775 END IF 791 776 ! 792 777 ! Compute associated depths at U and V points: … … 806 791 END DO 807 792 808 IF( ln_wd ) THEN809 zhust_e(:,:) = MAX(zhust_e (:,:), rn_wdmin1 )810 zhvst_e(:,:) = MAX(zhvst_e (:,:), rn_wdmin1 )811 END IF812 813 793 ENDIF 814 794 ! … … 861 841 ! 862 842 ! Add tidal astronomical forcing if defined 863 IF ( l k_tide.AND.ln_tide_pot ) THEN843 IF ( ln_tide.AND.ln_tide_pot ) THEN 864 844 DO jj = 2, jpjm1 865 845 DO ji = fs_2, fs_jpim1 ! vector opt. … … 888 868 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj) 889 869 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) 870 zwx(ji,jj) = zu_spg * zcpx(ji,jj) * wdmask(ji,jj) * wdmask(ji+1, jj) 871 zwy(ji,jj) = zv_spg * zcpy(ji,jj) * wdmask(ji,jj) * wdmask(ji, jj+1) 892 872 END DO 893 873 END DO … … 927 907 DO ji = fs_2, fs_jpim1 ! vector opt. 928 908 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 909 zhura = hu_0(ji,jj) + zsshu_a(ji,jj) 910 zhvra = hv_0(ji,jj) + zsshv_a(ji,jj) 936 911 zhura = ssumask(ji,jj)/(zhura + 1._wp - ssumask(ji,jj)) 937 912 zhvra = ssvmask(ji,jj)/(zhvra + 1._wp - ssvmask(ji,jj)) … … 953 928 ! 954 929 IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only) 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 930 hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 931 hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 962 932 hur_e(:,:) = ssumask(:,:) / ( hu_e(:,:) + 1._wp - ssumask(:,:) ) 963 933 hvr_e(:,:) = ssvmask(:,:) / ( hv_e(:,:) + 1._wp - ssvmask(:,:) ) … … 967 937 CALL lbc_lnk_multi( ua_e, 'U', -1._wp, va_e , 'V', -1._wp ) 968 938 ! 969 #if defined key_bdy970 939 ! ! open boundaries 971 IF( l k_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e )972 #endif 940 IF( ln_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e ) 941 973 942 #if defined key_agrif 974 943 IF( .NOT.Agrif_Root() ) CALL agrif_dyn_ts( jn ) ! Agrif … … 1024 993 ! 1025 994 ! Update barotropic trend: 1026 IF( ln_dynadv_vec .OR. ln_linssh ) THEN 1027 DO jk=1,jpkm1 1028 ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b 1029 va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * z1_2dt_b 1030 END DO 995 IF(ln_wd) THEN 996 IF( ln_dynadv_vec .OR. ln_linssh ) THEN 997 DO jk=1,jpkm1 998 ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b * wdmask(:,:) 999 va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * z1_2dt_b * wdmask(:,:) 1000 END DO 1001 ELSE 1002 ! At this stage, ssha has been corrected: compute new depths at velocity points 1003 DO jj = 1, jpjm1 1004 DO ji = 1, jpim1 ! NO Vector Opt. 1005 zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e1e2u(ji,jj) & 1006 & * ( e1e2t(ji ,jj) * ssha(ji ,jj) & 1007 & + e1e2t(ji+1,jj) * ssha(ji+1,jj) ) 1008 zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e1e2v(ji,jj) & 1009 & * ( e1e2t(ji,jj ) * ssha(ji,jj ) & 1010 & + e1e2t(ji,jj+1) * ssha(ji,jj+1) ) 1011 END DO 1012 END DO 1013 CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions 1014 ! 1015 DO jk=1,jpkm1 1016 ua(:,:,jk) = ua(:,:,jk) + r1_hu_n(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * z1_2dt_b * wdmask(:,:) 1017 va(:,:,jk) = va(:,:,jk) + r1_hv_n(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * z1_2dt_b * wdmask(:,:) 1018 END DO 1019 ! Save barotropic velocities not transport: 1020 ua_b(:,:) = ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) ) 1021 va_b(:,:) = va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) ) 1022 ENDIF 1031 1023 ELSE 1032 ! At this stage, ssha has been corrected: compute new depths at velocity points 1033 DO jj = 1, jpjm1 1034 DO ji = 1, jpim1 ! NO Vector Opt. 1035 zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e1e2u(ji,jj) & 1036 & * ( e1e2t(ji ,jj) * ssha(ji ,jj) & 1037 & + e1e2t(ji+1,jj) * ssha(ji+1,jj) ) 1038 zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e1e2v(ji,jj) & 1039 & * ( e1e2t(ji,jj ) * ssha(ji,jj ) & 1040 & + e1e2t(ji,jj+1) * ssha(ji,jj+1) ) 1041 END DO 1042 END DO 1043 CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions 1044 ! 1045 DO jk=1,jpkm1 1046 ua(:,:,jk) = ua(:,:,jk) + r1_hu_n(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * z1_2dt_b 1047 va(:,:,jk) = va(:,:,jk) + r1_hv_n(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * z1_2dt_b 1048 END DO 1049 ! Save barotropic velocities not transport: 1050 ua_b(:,:) = ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) ) 1051 va_b(:,:) = va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) ) 1052 ENDIF 1024 IF( ln_dynadv_vec .OR. ln_linssh ) THEN 1025 DO jk=1,jpkm1 1026 ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b 1027 va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * z1_2dt_b 1028 END DO 1029 ELSE 1030 ! At this stage, ssha has been corrected: compute new depths at velocity points 1031 DO jj = 1, jpjm1 1032 DO ji = 1, jpim1 ! NO Vector Opt. 1033 zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e1e2u(ji,jj) & 1034 & * ( e1e2t(ji ,jj) * ssha(ji ,jj) & 1035 & + e1e2t(ji+1,jj) * ssha(ji+1,jj) ) 1036 zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e1e2v(ji,jj) & 1037 & * ( e1e2t(ji,jj ) * ssha(ji,jj ) & 1038 & + e1e2t(ji,jj+1) * ssha(ji,jj+1) ) 1039 END DO 1040 END DO 1041 CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions 1042 ! 1043 DO jk=1,jpkm1 1044 ua(:,:,jk) = ua(:,:,jk) + r1_hu_n(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * z1_2dt_b 1045 va(:,:,jk) = va(:,:,jk) + r1_hv_n(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * z1_2dt_b 1046 END DO 1047 ! Save barotropic velocities not transport: 1048 ua_b(:,:) = ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) ) 1049 va_b(:,:) = va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) ) 1050 ENDIF 1051 1052 END IF 1053 1053 ! 1054 1054 DO jk = 1, jpkm1 … … 1086 1086 CALL wrk_dealloc( jpi,jpj, zsshu_a, zsshv_a ) 1087 1087 CALL wrk_dealloc( jpi,jpj, zhf ) 1088 IF( ln_wd ) CALL wrk_dealloc( jpi, jpj, zcpx, zcpy , wduflt1, wdvflt1)1088 IF( ln_wd ) CALL wrk_dealloc( jpi, jpj, zcpx, zcpy ) 1089 1089 ! 1090 1090 IF ( ln_diatmb ) THEN -
branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
r6152 r7412 22 22 USE divhor ! horizontal divergence 23 23 USE phycst ! physical constants 24 USE bdy_oce ! 25 USE bdy_par ! 24 USE bdy_oce , ONLY: ln_bdy, bdytmask 26 25 USE bdydyn2d ! bdy_ssh routine 27 26 #if defined key_agrif … … 88 87 ENDIF 89 88 ! 90 CALL div_hor( kt ) ! Horizontal divergence 91 ! 92 z2dt = 2._wp * rdt ! set time step size (Euler/Leapfrog) 89 z2dt = 2._wp * rdt ! set time step size (Euler/Leapfrog) 93 90 IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt 91 zcoef = 0.5_wp * r1_rau0 94 92 95 93 ! !------------------------------! 96 94 ! ! After Sea Surface Height ! 97 95 ! !------------------------------! 96 IF(ln_wd) THEN 97 CALL wad_lmt(sshb, zcoef * (emp_b(:,:) + emp(:,:)), z2dt) 98 ENDIF 99 100 CALL div_hor( kt ) ! Horizontal divergence 101 ! 98 102 zhdiv(:,:) = 0._wp 99 103 DO jk = 1, jpkm1 ! Horizontal divergence of barotropic transports … … 104 108 ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations. 105 109 ! 106 zcoef = 0.5_wp * r1_rau0107 108 IF(ln_wd) CALL wad_lmt(sshb, zcoef * (emp_b(:,:) + emp(:,:)), z2dt)109 110 110 ssha(:,:) = ( sshb(:,:) - z2dt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * ssmask(:,:) 111 111 … … 116 116 CALL agrif_ssh( kt ) 117 117 # endif 118 # if defined key_bdy 119 IF( lk_bdy ) THEN 118 IF( ln_bdy ) THEN 120 119 CALL lbc_lnk( ssha, 'T', 1. ) ! Not sure that's necessary 121 120 CALL bdy_ssh( ssha ) ! Duplicate sea level across open boundaries 122 121 ENDIF 123 # endif124 122 ENDIF 125 123 … … 211 209 ENDIF 212 210 213 #if defined key_bdy 214 IF( lk_bdy ) THEN 211 IF( ln_bdy ) THEN 215 212 DO jk = 1, jpkm1 216 213 wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:) 217 214 END DO 218 215 ENDIF 219 #endif220 216 ! 221 217 IF( nn_timing == 1 ) CALL timing_stop('wzv') -
branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/DYN/wet_dry.F90
r6152 r7412 33 33 !! --------------------------------------------------------------------- 34 34 35 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: wduflt, wdvflt !: u- and v- filter36 35 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: wdmask !: u- and v- limiter 37 36 … … 46 45 PUBLIC wad_lmt ! routine called by sshwzv.F90 47 46 PUBLIC wad_lmt_bt ! routine called by dynspg_ts.F90 47 PUBLIC wad_istate ! routine called by istate.F90 and domvvl.F90 48 48 49 49 !! * Substitutions … … 87 87 88 88 IF(ln_wd) THEN 89 ALLOCATE( wd uflt(jpi,jpj), wdvflt(jpi,jpj), wdmask(jpi,jpj), STAT=ierr )89 ALLOCATE( wdmask(jpi,jpj), STAT=ierr ) 90 90 IF( ierr /= 0 ) CALL ctl_stop('STOP', 'wad_init : Array allocation error') 91 91 ENDIF … … 145 145 ! Horizontal Flux in u and v direction 146 146 DO jk = 1, jpkm1 147 DO jj = 1, jpj m1148 DO ji = 1, jpi m1147 DO jj = 1, jpj 148 DO ji = 1, jpi 149 149 zflxu(ji,jj) = zflxu(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk) 150 150 zflxv(ji,jj) = zflxv(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk) … … 156 156 zflxv(:,:) = zflxv(:,:) * e1v(:,:) 157 157 158 DO jj = 2, jpjm1 159 DO ji = 2, jpim1 158 wdmask(:,:) = 1 159 DO jj = 2, jpj 160 DO ji = 2, jpi 160 161 161 162 IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE ! we don't care about land cells … … 168 169 169 170 zdep2 = bathy(ji,jj) + sshb1(ji,jj) - rn_wdmin1 170 IF(zdep2 <0._wp) THEN !add more safty, but not necessary171 IF(zdep2 .le. 0._wp) THEN !add more safty, but not necessary 171 172 !zdep2 = 0._wp 172 173 sshb1(ji,jj) = rn_wdmin1 - bathy(ji,jj) 174 wdmask(ji,jj) = 0._wp 173 175 END IF 174 176 ENDDO … … 183 185 zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:) 184 186 185 DO jj = 2, jpj m1186 DO ji = 2, jpi m1187 DO jj = 2, jpj 188 DO ji = 2, jpi 187 189 188 wdmask(ji,jj) = 0189 190 IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE 190 191 IF(bathy(ji,jj) > zdepwd) CYCLE … … 202 203 IF(zdep1 > zdep2) THEN 203 204 zflag = 1 204 wdmask(ji, jj) = 1205 wdmask(ji, jj) = 0 205 206 zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt ) 206 207 zcoef = max(zcoef, 0._wp) … … 209 210 IF(zflxu1(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = zcoef 210 211 IF(zflxv1(ji, jj) > 0._wp) zwdlmtv(ji ,jj) = zcoef 211 IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji -1,jj) = zcoef212 IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = zcoef 212 213 END IF 213 214 END DO ! ji loop … … 231 232 CALL lbc_lnk( un, 'U', -1. ) 232 233 CALL lbc_lnk( vn, 'V', -1. ) 234 ! 235 un_b(:,:) = un_b(:,:) * zwdlmtu(:, :) 236 vn_b(:,:) = vn_b(:,:) * zwdlmtv(:, :) 237 CALL lbc_lnk( un_b, 'U', -1. ) 238 CALL lbc_lnk( vn_b, 'V', -1. ) 233 239 234 240 IF(zflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt!!!' … … 291 297 zflxp(:,:) = 0._wp 292 298 zflxn(:,:) = 0._wp 293 !zflxu(:,:) = 0._wp294 !zflxv(:,:) = 0._wp295 299 296 300 zwdlmtu(:,:) = 1._wp … … 299 303 ! Horizontal Flux in u and v direction 300 304 301 !zflxu(:,:) = zflxu(:,:) * e2u(:,:) 302 !zflxv(:,:) = zflxv(:,:) * e1v(:,:) 303 304 DO jj = 2, jpjm1 305 DO ji = 2, jpim1 305 DO jj = 2, jpj 306 DO ji = 2, jpi 306 307 307 308 IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE ! we don't care about land cells … … 314 315 315 316 zdep2 = bathy(ji,jj) + sshn_e(ji,jj) - rn_wdmin1 316 IF(zdep2 < 0._wp) THEN !add more safty, but not necessary317 !zdep2 = 0._wp318 sshn_e(ji,jj) = rn_wdmin1 - bathy(ji,jj)319 END IF320 317 ENDDO 321 318 END DO … … 329 326 zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:) 330 327 331 DO jj = 2, jpj m1332 DO ji = 2, jpi m1328 DO jj = 2, jpj 329 DO ji = 2, jpi 333 330 334 wdmask(ji,jj) = 0335 331 IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE 336 332 IF(bathy(ji,jj) > zdepwd) CYCLE … … 349 345 IF(zdep1 > zdep2) THEN 350 346 zflag = 1 351 !wdmask(ji, jj) = 1352 347 zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt ) 353 348 zcoef = max(zcoef, 0._wp) … … 356 351 IF(zflxu1(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = zcoef 357 352 IF(zflxv1(ji, jj) > 0._wp) zwdlmtv(ji ,jj) = zcoef 358 IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji -1,jj) = zcoef353 IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = zcoef 359 354 END IF 360 355 END DO ! ji loop … … 379 374 IF(zflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt_bt!!!' 380 375 381 !IF( ln_rnf ) CALL sbc_rnf_div( hdivn ) ! runoffs (update hdivn field)382 !IF( nn_cla == 1 ) CALL cla_div ( kt ) ! Cross Land Advection (update hdivn field)383 376 ! 384 377 ! … … 390 383 IF( nn_timing == 1 ) CALL timing_stop('wad_lmt') 391 384 END SUBROUTINE wad_lmt_bt 385 386 SUBROUTINE wad_istate 387 !!---------------------------------------------------------------------- 388 !! *** ROUTINE wad_istate *** 389 !! 390 !! ** Purpose : Initialization of the dynamics and tracers for WAD test 391 !! configurations (channels or bowls with initial ssh gradients) 392 !! 393 !! ** Method : - set temperature field 394 !! - set salinity field 395 !! - set ssh slope (needs to be repeated in domvvl_rst_init to 396 !! set vertical metrics ) 397 !!---------------------------------------------------------------------- 398 ! 399 INTEGER :: ji, jj ! dummy loop indices 400 REAL(wp) :: zi, zj 401 !!---------------------------------------------------------------------- 402 ! 403 ! Uniform T & S in all test cases 404 tsn(:,:,:,jp_tem) = 10._wp 405 tsb(:,:,:,jp_tem) = 10._wp 406 tsn(:,:,:,jp_sal) = 35._wp 407 tsb(:,:,:,jp_sal) = 35._wp 408 SELECT CASE ( jp_cfg ) 409 ! ! ==================== 410 CASE ( 1 ) ! WAD 1 configuration 411 ! ! ==================== 412 ! 413 IF(lwp) WRITE(numout,*) 414 IF(lwp) WRITE(numout,*) 'istate_wad : Closed box with EW linear bottom slope' 415 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 416 ! 417 do ji = 1,jpi 418 sshn(ji,:) = ( -5.5_wp + 5.5_wp*FLOAT(mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1) 419 end do 420 ! ! ==================== 421 CASE ( 2 ) ! WAD 2 configuration 422 ! ! ==================== 423 ! 424 IF(lwp) WRITE(numout,*) 425 IF(lwp) WRITE(numout,*) 'istate_wad : Parobolic EW channel, mid-range initial ssh slope' 426 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 427 ! 428 do ji = 1,jpi 429 sshn(ji,:) = ( -5.5_wp + 3.9_wp*FLOAT(jpidta - mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1) 430 end do 431 ! ! ==================== 432 CASE ( 3 ) ! WAD 3 configuration 433 ! ! ==================== 434 ! 435 IF(lwp) WRITE(numout,*) 436 IF(lwp) WRITE(numout,*) 'istate_wad : Parobolic EW channel, extreme initial ssh slope' 437 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 438 ! 439 do ji = 1,jpi 440 sshn(ji,:) = ( -7.5_wp + 6.9_wp*FLOAT(jpidta - mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1) 441 end do 442 443 ! 444 ! ! ==================== 445 CASE ( 4 ) ! WAD 4 configuration 446 ! ! ==================== 447 ! 448 IF(lwp) WRITE(numout,*) 449 IF(lwp) WRITE(numout,*) 'istate_wad : Parobolic bowl, mid-range initial ssh slope' 450 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 451 ! 452 DO ji = 1, jpi 453 zi = MAX(1.0-FLOAT((mig(ji)-25)**2)/400.0, 0.0 ) 454 DO jj = 1, jpj 455 zj = MAX(1.0-FLOAT((mjg(jj)-17)**2)/144.0, 0.0 ) 456 sshn(ji,jj) = -8.5_wp + 8.5_wp*zi*zj 457 END DO 458 END DO 459 460 ! 461 ! ! =========================== 462 CASE ( 5 ) ! WAD 5 configuration 463 ! ! ==================== 464 ! 465 IF(lwp) WRITE(numout,*) 466 IF(lwp) WRITE(numout,*) 'istate_wad : Double slope with shelf' 467 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 468 ! 469 ! Needed rn_wdmin2 increased to 0.01 for this case? 470 do ji = 1,jpi 471 sshn(ji,:) = ( -5.5_wp + 9.0_wp*FLOAT(mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1) 472 end do 473 474 ! 475 ! ! =========================== 476 CASE ( 6 ) ! WAD 6 configuration 477 ! ! ==================== 478 ! 479 IF(lwp) WRITE(numout,*) 480 IF(lwp) WRITE(numout,*) 'istate_wad : Parobolic EW channel with gaussian ridge' 481 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 482 ! 483 do ji = 1,jpi 484 !6a 485 sshn(ji,:) = ( -5.5_wp + 9.0_wp*FLOAT(jpidta - mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1) 486 !Some variations in initial slope that have been tested 487 !6b 488 !sshn(ji,:) = ( -5.5_wp + 6.5_wp*FLOAT(jpidta - mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1) 489 !6c 490 !sshn(ji,:) = ( -5.5_wp + 7.5_wp*FLOAT(jpidta - mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1) 491 !6d 492 !sshn(ji,:) = ( -4.5_wp + 8.0_wp*FLOAT(jpidta - mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1) 493 end do 494 495 ! 496 ! ! =========================== 497 CASE DEFAULT ! NONE existing configuration 498 ! ! =========================== 499 WRITE(ctmp1,*) 'WAD test with a ', jp_cfg,' option is not coded' 500 ! 501 CALL ctl_stop( ctmp1 ) 502 ! 503 END SELECT 504 ! 505 ! Apply minimum wetdepth criterion 506 ! 507 do jj = 1,jpj 508 do ji = 1,jpi 509 IF( bathy(ji,jj) + sshn(ji,jj) < rn_wdmin1 ) THEN 510 sshn(ji,jj) = tmask(ji,jj,1)*( rn_wdmin1 - bathy(ji,jj) ) 511 ENDIF 512 end do 513 end do 514 sshb = sshn 515 ssha = sshn 516 ! 517 END SUBROUTINE wad_istate 518 519 !!===================================================================== 392 520 END MODULE wet_dry
Note: See TracChangeset
for help on using the changeset viewer.
