Changeset 5842 for branches/2015
- Timestamp:
- 2015-10-30T16:34:24+01:00 (8 years ago)
- Location:
- branches/2015/dev_r5803_NOC_WAD/NEMOGCM
- Files:
-
- 1 added
- 9 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5803_NOC_WAD/NEMOGCM/CONFIG/AMM12/cpp_AMM12.fcm
r4245 r5842 1 bld::tool::fppkeys key_bdy key_tide key_dynspg_ts key_ldfslp key_zdfgls key_vvl key_diainstant key_mpp_mpi key_iomput1 bld::tool::fppkeys key_bdy key_tide key_dynspg_ts key_ldfslp key_zdfgls key_vvl key_diainstant key_mpp_mpi -
branches/2015/dev_r5803_NOC_WAD/NEMOGCM/CONFIG/cfg.txt
r5656 r5842 5 5 C1D_PAPA OPA_SRC 6 6 GYRE_BFM OPA_SRC TOP_SRC 7 AMM12 OPA_SRC8 7 ORCA2_LIM_PISCES OPA_SRC LIM_SRC_2 NST_SRC TOP_SRC 9 8 GYRE OPA_SRC … … 11 10 ORCA2_LIM OPA_SRC LIM_SRC_2 NST_SRC 12 11 ORCA2_OFF_PISCES OPA_SRC OFF_SRC TOP_SRC 12 AMM12 OPA_SRC 13 IRS18WD OPA_SRC -
branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90
r5123 r5842 10 10 !! 3.5 ! 2012 (S. Mocavero, I. Epicoco) Add arrays associated 11 11 !! to the optimization of BDY communications 12 !! 3.6.?! 2014 (H. Liu) Add arrays associated with Wetting and Drying 12 13 !!---------------------------------------------------------------------- 13 14 … … 270 271 INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: nicoa, njcoa !: ??? 271 272 #endif 273 274 !!---------------------------------------------------------------------- 275 !! critical depths,filters, limiters,and masks for Wetting and Drying 276 !! --------------------------------------------------------------------- 277 278 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wduflt, wdvflt !: u- and v- filter 279 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wdmask !: u- and v- limiter 280 281 LOGICAL, PUBLIC, SAVE :: ln_wd !: key to turn on/off wetting/drying (T: on, F: off) 282 REAL(wp), PUBLIC, SAVE :: rn_wdmin1 !: minimum water depth on dried cells 283 REAL(wp), PUBLIC, SAVE :: rn_wdmin2 !: tolerrance of minimum water depth on dried cells 284 REAL(wp), PUBLIC, SAVE :: rn_wdld !: land elevation below which wetting/drying will be considered 285 INTEGER , PUBLIC, SAVE :: nn_wdit !: maximum number of iteration for W/D limiter 286 287 !LOGICAL, PUBLIC, SAVE :: ln_wd = .FALSE. !: turn on wetting/drying (T: on, F: off) 288 !REAL(wp), PUBLIC, SAVE :: rn_wdmin1 = 0.10_wp !: minimum water depth on dried cells 289 !REAL(wp), PUBLIC, SAVE :: rn_wdmin2 = 0.01_wp !: tolerrance of minimum water depth on dried cells 290 !REAL(wp), PUBLIC, SAVE :: rn_wdld = 20.0_wp !: land elevation below which wetting/drying will be considered 291 !INTEGER , PUBLIC, SAVE :: nn_wdit = 10 !: maximum number of iteration for W/D limiter 272 292 273 293 !!---------------------------------------------------------------------- … … 407 427 ALLOCATE( npcoa(4,jpk), nicoa(2*(jpi+jpj),4,jpk), njcoa(2*(jpi+jpj),4,jpk), STAT=ierr(12) ) 408 428 #endif 429 430 IF(ln_wd) & 431 ALLOCATE( wduflt(jpi,jpj), wdvflt(jpi,jpj), wdmask(jpi,jpj), STAT=ierr(12) ) 409 432 ! 410 433 dom_oce_alloc = MAXVAL(ierr) -
branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
r5506 r5842 147 147 fse3t_a(:,:,jpk) = e3t_0(:,:,jpk) 148 148 149 !IF(ln_wd) THEN 150 ! DO jj = 1, jpj 151 ! DO ji = 1, jpi 152 ! IF(mbathy(ji,jj) == 2 .AND. e3t_0(ji,jj,1) <= 0.5_wp * rn_wdmin1) THEN 153 ! fse3t_a(ji,jj,1:2) = 0.5_wp * rn_wdmin1 154 ! fse3t_n(ji,jj,1:2) = 0.5_wp * rn_wdmin1 155 ! fse3t_b(ji,jj,1:2) = 0.5_wp * rn_wdmin1 156 ! END IF 157 ! ENDDO 158 ! ENDDO 159 !END IF 160 149 161 ! Reconstruction of all vertical scale factors at now and before time steps 150 162 ! ============================================================================= … … 701 713 ! ! = 'U', 'V', 'W, 'F', 'UW' or 'VW' 702 714 !! * Local declarations 715 REAL(wp) :: zwad ! = 1.0 when ln_wd = .true. 716 ! = 0.0 when ln_wd = .false. 717 ! 703 718 INTEGER :: ji, jj, jk ! dummy loop indices 704 719 LOGICAL :: l_is_orca ! local logical … … 706 721 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_interpol') 707 722 ! 723 IF(ln_wd) THEN 724 zwad = 1.0_wp 725 ELSE 726 zwad = 0.0_wp 727 END IF 728 708 729 l_is_orca = .FALSE. 709 730 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) l_is_orca = .TRUE. ! ORCA R2 configuration - will need to correct some locations … … 717 738 DO jj = 1, jpjm1 718 739 DO ji = 1, fs_jpim1 ! vector opt. 719 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * r1_e12u(ji,jj) & 740 !pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * r1_e12u(ji,jj) & 741 pe3_out(ji,jj,jk) = 0.5_wp * (umask(ji,jj,jk) * (1.0_wp - zwad) + zwad) * r1_e12u(ji,jj) & 720 742 & * ( e12t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) & 721 743 & + e12t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) ) … … 735 757 DO jj = 1, jpjm1 736 758 DO ji = 1, fs_jpim1 ! vector opt. 737 pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) * r1_e12v(ji,jj) & 759 !pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) * r1_e12v(ji,jj) & 760 pe3_out(ji,jj,jk) = 0.5_wp * (vmask(ji,jj,jk) * (1.0_wp - zwad) + zwad) * r1_e12v(ji,jj) & 738 761 & * ( e12t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) & 739 762 & + e12t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) ) … … 753 776 DO jj = 1, jpjm1 754 777 DO ji = 1, fs_jpim1 ! vector opt. 755 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) * r1_e12f(ji,jj) & 778 !pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) * r1_e12f(ji,jj) & 779 pe3_out(ji,jj,jk) = 0.5_wp * (umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zwad) + zwad) & 780 & * r1_e12f(ji,jj) & 756 781 & * ( e12u(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3u_0(ji,jj ,jk) ) & 757 782 & + e12u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) ) … … 771 796 ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing 772 797 DO jk = 2, jpk 773 pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * tmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) ) & 774 & + 0.5_wp * tmask(:,:,jk) * ( pe3_in(:,:,jk ) - e3t_0(:,:,jk ) ) 798 !pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * tmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) ) & 799 ! & + 0.5_wp * tmask(:,:,jk) * ( pe3_in(:,:,jk ) - e3t_0(:,:,jk ) ) 800 pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * (tmask(:,:,jk) * (1.0_wp - zwad) + zwad) ) & 801 & * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) ) & 802 & + 0.5_wp * (tmask(:,:,jk) * (1.0_wp - zwad) + zwad) & 803 & * ( pe3_in(:,:,jk ) - e3t_0(:,:,jk ) ) 775 804 END DO 776 805 ! ! -------------------------------------- ! … … 781 810 ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing 782 811 DO jk = 2, jpk 783 pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * umask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) ) & 784 & + 0.5_wp * umask(:,:,jk) * ( pe3_in(:,:,jk ) - e3u_0(:,:,jk ) ) 812 !pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * umask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) ) & 813 ! & + 0.5_wp * umask(:,:,jk) * ( pe3_in(:,:,jk ) - e3u_0(:,:,jk ) ) 814 pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * (umask(:,:,jk) * (1.0_wp - zwad) + zwad) ) & 815 & * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) ) & 816 & + 0.5_wp * (umask(:,:,jk) * (1.0_wp - zwad) + zwad) & 817 & * ( pe3_in(:,:,jk ) - e3u_0(:,:,jk ) ) 785 818 END DO 786 819 ! ! -------------------------------------- ! … … 791 824 ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing 792 825 DO jk = 2, jpk 793 pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * vmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) ) & 794 & + 0.5_wp * vmask(:,:,jk) * ( pe3_in(:,:,jk ) - e3v_0(:,:,jk ) ) 826 !pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * vmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) ) & 827 ! & + 0.5_wp * vmask(:,:,jk) * ( pe3_in(:,:,jk ) - e3v_0(:,:,jk ) ) 828 pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * (vmask(:,:,jk) * (1.0_wp - zwad) + zwad) ) & 829 & * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) ) & 830 & + 0.5_wp * (vmask(:,:,jk) * (1.0_wp - zwad) + zwad) & 831 & * ( pe3_in(:,:,jk ) - e3v_0(:,:,jk ) ) 795 832 END DO 796 833 END SELECT … … 817 854 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 818 855 !! * Local declarations 819 INTEGER :: j k856 INTEGER :: ji, jj, jk 820 857 INTEGER :: id1, id2, id3, id4, id5 ! local integers 821 858 !!---------------------------------------------------------------------- … … 905 942 fse3t_n(:,:,:) = e3t_0(:,:,:) 906 943 sshn(:,:) = 0.0_wp 944 945 IF(ln_wd) THEN 946 DO jj = 1, jpj 947 DO ji = 1, jpi 948 !IF(e3t_0(ji,jj,1) < 0._wp) THEN 949 !IF(mbathy(ji,jj) == 2 .AND. e3t_0(ji,jj,1) <= 0.5_wp * rn_wdmin1) THEN 950 IF( e3t_0(ji,jj,1) <= 0.5_wp * rn_wdmin1) THEN 951 fse3t_b(ji,jj,:) = 0.5_wp * rn_wdmin1 952 fse3t_n(ji,jj,:) = 0.5_wp * rn_wdmin1 953 fse3t_a(ji,jj,:) = 0.5_wp * rn_wdmin1 954 sshb(ji,jj) = rn_wdmin1 - bathy(ji,jj) 955 sshn(ji,jj) = rn_wdmin1 - bathy(ji,jj) 956 ssha(ji,jj) = rn_wdmin1 - bathy(ji,jj) 957 ENDIF 958 ENDDO 959 ENDDO 960 END IF 961 907 962 IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN 908 963 tilde_e3t_b(:,:,:) = 0.0_wp -
branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
r5656 r5842 18 18 !! 3.4 ! 2012-12 (R. Bourdalle-Badie and G. Reffray) modify C1D case 19 19 !! 3.6 ! 2014-11 (P. Mathiot and C. Harris) add ice shelf capabilitye 20 !! 3.? ! 2015-11 (H. Liu) Modifications for Wetting/Drying 20 21 !!---------------------------------------------------------------------- 21 22 … … 1810 1811 REAL(wp) :: zrmax, ztaper ! temporary scalars 1811 1812 REAL(wp) :: zrfact 1813 REAL(wp) :: zsmth 1812 1814 ! 1813 1815 REAL(wp), POINTER, DIMENSION(:,: ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2 … … 1865 1867 bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) ) 1866 1868 1869 IF(.NOT.ln_wd) THEN 1867 1870 DO jj = 1, jpj 1868 1871 DO ji = 1, jpi … … 1870 1873 END DO 1871 1874 END DO 1875 END IF 1872 1876 ! ! ============================= 1873 1877 ! ! Define the envelop bathymetry (hbatt) … … 1876 1880 zenv(:,:) = bathy(:,:) 1877 1881 ! 1882 IF(.NOT.ln_wd) THEN 1878 1883 ! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing 1879 1884 DO jj = 1, jpj … … 1891 1896 END DO 1892 1897 END DO 1898 END IF 1899 1893 1900 ! apply lateral boundary condition CAUTION: keep the value when the lbc field is zero 1894 1901 CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' ) … … 1930 1937 zri(:,:) = 0._wp 1931 1938 zrj(:,:) = 0._wp 1939 1940 !IF(ln_wd) THEN !extend the smoothed region to cover the W/D zones 1941 ! zsmth = -rn_wdld 1942 !ELSE 1943 zsmth = 0._wp ! The original form (only smooth ocean points) 1944 !ENDIF 1945 1932 1946 DO jj = 1, nlcj 1933 1947 DO ji = 1, nlci 1934 1948 iip1 = MIN( ji+1, nlci ) ! force zri = 0 on last line (ji=ncli+1 to jpi) 1935 1949 ijp1 = MIN( jj+1, nlcj ) ! force zrj = 0 on last raw (jj=nclj+1 to jpj) 1936 IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(iip1,jj) > 0._wp)) THEN1950 IF( (zenv(ji,jj) > zsmth) .AND. (zenv(iip1,jj) > zsmth)) THEN 1937 1951 zri(ji,jj) = ( zenv(iip1,jj ) - zenv(ji,jj) ) / ( zenv(iip1,jj ) + zenv(ji,jj) ) 1938 1952 END IF 1939 IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(ji,ijp1) > 0._wp)) THEN1953 IF( (zenv(ji,jj) > zsmth) .AND. (zenv(ji,ijp1) > zsmth)) THEN 1940 1954 zrj(ji,jj) = ( zenv(ji ,ijp1) - zenv(ji,jj) ) / ( zenv(ji ,ijp1) + zenv(ji,jj) ) 1941 1955 END IF … … 1960 1974 END DO ! End loop ! 1961 1975 ! ! ================ ! 1962 DO jj = 1, jpj 1963 DO ji = 1, jpi 1964 zenv(ji,jj) = MAX( zenv(ji,jj), rn_sbot_min ) ! set all points to avoid undefined scale value warnings 1965 END DO 1966 END DO 1976 IF(ln_wd) THEN 1977 DO jj = 1, jpj 1978 DO ji = 1, jpi 1979 !zenv(ji,jj) = MAX( zenv(ji,jj), -rn_wdld ) ! filt out land bathy data 1980 END DO 1981 END DO 1982 ELSE 1983 DO jj = 1, jpj 1984 DO ji = 1, jpi 1985 zenv(ji,jj) = MAX( zenv(ji,jj), rn_sbot_min ) ! set all points to avoid undefined scale value warnings 1986 END DO 1987 END DO 1988 ENDIF 1967 1989 ! 1968 1990 ! Envelope bathymetry saved in hbatt … … 1994 2016 IF(lwp) THEN 1995 2017 WRITE(numout,*) 1996 WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min 2018 IF(.NOT.ln_wd) THEN 2019 WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min 2020 ELSE 2021 WRITE(numout,*) ' zgr_sco: minimum positive depth of the envelop topography set to : ', rn_sbot_min 2022 WRITE(numout,*) ' zgr_sco: minimum negative depth of the envelop topography set to : ', -rn_wdld 2023 ENDIF 1997 2024 ENDIF 1998 2025 hbatu(:,:) = rn_sbot_min … … 2007 2034 END DO 2008 2035 END DO 2036 2037 IF(ln_wd) THEN !avoid the zero depth on T- (U-,V-,F-) points 2038 DO jj = 1, jpj 2039 DO ji = 1, jpi 2040 IF(ABS(hbatt(ji,jj)) < rn_wdmin1) & 2041 & hbatt(ji,jj) = SIGN(1._wp, hbatt(ji,jj)) * rn_wdmin1 2042 IF(ABS(hbatu(ji,jj)) < rn_wdmin1) & 2043 & hbatu(ji,jj) = SIGN(1._wp, hbatu(ji,jj)) * rn_wdmin1 2044 IF(ABS(hbatv(ji,jj)) < rn_wdmin1) & 2045 & hbatv(ji,jj) = SIGN(1._wp, hbatv(ji,jj)) * rn_wdmin1 2046 IF(ABS(hbatf(ji,jj)) < rn_wdmin1) & 2047 & hbatf(ji,jj) = SIGN(1._wp, hbatf(ji,jj)) * rn_wdmin1 2048 END DO 2049 END DO 2050 END IF 2009 2051 ! 2010 2052 ! Apply lateral boundary condition … … 2014 2056 DO ji = 1, jpi 2015 2057 IF( hbatu(ji,jj) == 0._wp ) THEN 2058 !No worries about the following line when ln_wd == .true. 2016 2059 IF( zhbat(ji,jj) == 0._wp ) hbatu(ji,jj) = rn_sbot_min 2017 2060 IF( zhbat(ji,jj) /= 0._wp ) hbatu(ji,jj) = zhbat(ji,jj) … … 2039 2082 2040 2083 !!bug: key_helsinki a verifer 2041 hift(:,:) = MIN( hift(:,:), hbatt(:,:) ) 2042 hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) ) 2043 hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) ) 2044 hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) ) 2084 IF(.NOT.ln_wd) THEN 2085 hift(:,:) = MIN( hift(:,:), hbatt(:,:) ) 2086 hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) ) 2087 hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) ) 2088 hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) ) 2089 END IF 2045 2090 2046 2091 IF( nprint == 1 .AND. lwp ) THEN … … 2087 2132 fsde3w(:,:,:) = gdep3w_0(:,:,:) 2088 2133 ! 2089 where (e3t_0 (:,:,:).eq.0.0) e3t_0(:,:,:) = 1.0 2090 where (e3u_0 (:,:,:).eq.0.0) e3u_0(:,:,:) = 1.0 2091 where (e3v_0 (:,:,:).eq.0.0) e3v_0(:,:,:) = 1.0 2092 where (e3f_0 (:,:,:).eq.0.0) e3f_0(:,:,:) = 1.0 2093 where (e3w_0 (:,:,:).eq.0.0) e3w_0(:,:,:) = 1.0 2094 where (e3uw_0 (:,:,:).eq.0.0) e3uw_0(:,:,:) = 1.0 2095 where (e3vw_0 (:,:,:).eq.0.0) e3vw_0(:,:,:) = 1.0 2134 IF(.NOT.ln_wd) THEN 2135 where (e3t_0 (:,:,:).eq.0.0) e3t_0(:,:,:) = 1.0 2136 where (e3u_0 (:,:,:).eq.0.0) e3u_0(:,:,:) = 1.0 2137 where (e3v_0 (:,:,:).eq.0.0) e3v_0(:,:,:) = 1.0 2138 where (e3f_0 (:,:,:).eq.0.0) e3f_0(:,:,:) = 1.0 2139 where (e3w_0 (:,:,:).eq.0.0) e3w_0(:,:,:) = 1.0 2140 where (e3uw_0 (:,:,:).eq.0.0) e3uw_0(:,:,:) = 1.0 2141 where (e3vw_0 (:,:,:).eq.0.0) e3vw_0(:,:,:) = 1.0 2142 END IF 2096 2143 2097 2144 #if defined key_agrif … … 2135 2182 IF( scobot(ji,jj) >= fsdept(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk ) 2136 2183 END DO 2137 IF( scobot(ji,jj) == 0._wp ) mbathy(ji,jj) = 0 2138 END DO 2139 END DO 2184 2185 IF(ln_wd) THEN 2186 IF( scobot(ji,jj) <= -(rn_wdld - rn_wdmin2)) THEN 2187 mbathy(ji,jj) = 0 2188 ELSEIF(scobot(ji,jj) <= rn_wdmin1) THEN 2189 mbathy(ji,jj) = 2 2190 ENDIF 2191 ELSE 2192 IF( scobot(ji,jj) == 0._wp ) mbathy(ji,jj) = 0 2193 ENDIF 2194 END DO 2195 END DO 2196 2140 2197 IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ), & 2141 2198 & ' MAX ', MAXVAL( mbathy(:,:) ) … … 2257 2314 INTEGER :: ji, jj, jk ! dummy loop argument 2258 2315 REAL(wp) :: zcoeft, zcoefw ! temporary scalars 2316 REAL(wp) :: ztmpu, ztmpv, ztmpf 2317 REAL(wp) :: ztmpu1, ztmpv1, ztmpf1 2259 2318 ! 2260 2319 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 … … 2310 2369 DO ji = 1, jpim1 2311 2370 DO jj = 1, jpjm1 2371 ! extended for Wetting/Drying case 2372 ztmpu = hbatt(ji,jj)+hbatt(ji+1,jj) 2373 ztmpv = hbatt(ji,jj)+hbatt(ji,jj+1) 2374 ztmpf = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) 2375 ztmpu1 = hbatt(ji,jj)*hbatt(ji+1,jj) 2376 ztmpv1 = hbatt(ji,jj)*hbatt(ji,jj+1) 2377 ztmpf1 = MIN(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) * & 2378 & MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) 2312 2379 DO jk = 1, jpk 2313 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) & 2314 & / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 2315 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) & 2316 & / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 2317 z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) & 2318 & + hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) & 2319 & / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 2320 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) & 2321 & / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 2322 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) & 2323 & / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 2380 IF((ztmpu1 < 0._wp.OR.ABS(ztmpu) < rn_wdmin1).AND.ln_wd) THEN 2381 z_esigtu3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) ) 2382 z_esigwu3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji+1,jj,jk) ) 2383 ELSE 2384 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) & 2385 & / ztmpu 2386 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) & 2387 & / ztmpu 2388 END IF 2389 2390 IF((ztmpv1 < 0._wp.OR.ABS(ztmpv) < rn_wdmin1).AND.ln_wd) THEN 2391 z_esigtv3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji,jj+1,jk) ) 2392 z_esigwv3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji,jj+1,jk) ) 2393 ELSE 2394 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) & 2395 & / ztmpv 2396 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) & 2397 & / ztmpv 2398 END IF 2399 2400 IF((ztmpf1 < 0._wp.OR.ABS(ztmpf) < rn_wdmin1).AND.ln_wd) THEN 2401 z_esigtf3(ji,jj,jk) = 0.25_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) & 2402 & + z_esigt3(ji,jj+1,jk) + z_esigt3(ji+1,jj+1,jk) ) 2403 ELSE 2404 z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) & 2405 & + hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) & 2406 & / ztmpf 2407 END IF 2408 2409 ! 2410 !z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) & 2411 ! & / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 2412 !z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) & 2413 ! & / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 2414 !z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) & 2415 ! & + hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) & 2416 ! & / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 2417 !z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) & 2418 ! & / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 2419 !z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) & 2420 ! & / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 2324 2421 ! 2325 2422 e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) … … 2361 2458 REAL(wp) :: zsmth ! smoothing around critical depth 2362 2459 REAL(wp) :: zzs, zzb ! Surface and bottom cell thickness in sigma space 2460 REAL(wp) :: ztmpu, ztmpv, ztmpf 2461 REAL(wp) :: ztmpu1, ztmpv1, ztmpf1 2363 2462 ! 2364 2463 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 … … 2439 2538 DO jj=1,jpj-1 2440 2539 2441 DO jk = 1, jpk 2442 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) / & 2443 ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 2444 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) / & 2445 ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 2446 z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) + & 2447 hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / & 2448 ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 2449 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) / & 2450 ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 2451 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) / & 2452 ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 2540 ! extend to suit for Wetting/Drying case 2541 ztmpu = hbatt(ji,jj)+hbatt(ji+1,jj) 2542 ztmpv = hbatt(ji,jj)+hbatt(ji,jj+1) 2543 ztmpf = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) 2544 ztmpu1 = hbatt(ji,jj)*hbatt(ji+1,jj) 2545 ztmpv1 = hbatt(ji,jj)*hbatt(ji,jj+1) 2546 ztmpf1 = MIN(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) * & 2547 & MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) 2548 DO jk = 1, jpk 2549 IF((ztmpu1 < 0._wp.OR.ABS(ztmpu) < rn_wdmin1).AND.ln_wd) THEN 2550 z_esigtu3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) ) 2551 z_esigwu3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji+1,jj,jk) ) 2552 ELSE 2553 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) & 2554 & / ztmpu 2555 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) & 2556 & / ztmpu 2557 END IF 2558 2559 IF((ztmpv1 < 0._wp.OR.ABS(ztmpv) < rn_wdmin1).AND.ln_wd) THEN 2560 z_esigtv3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji,jj+1,jk) ) 2561 z_esigwv3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji,jj+1,jk) ) 2562 ELSE 2563 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) & 2564 & / ztmpv 2565 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) & 2566 & / ztmpv 2567 END IF 2568 2569 IF((ztmpf1 < 0._wp.OR.ABS(ztmpf) < rn_wdmin1).AND.ln_wd) THEN 2570 z_esigtf3(ji,jj,jk) = 0.25_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) & 2571 & + z_esigt3(ji,jj+1,jk) + z_esigt3(ji+1,jj+1,jk) ) 2572 ELSE 2573 z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) & 2574 & + hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) & 2575 & / ztmpf 2576 END IF 2577 2578 !z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) / & 2579 ! ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 2580 !z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) / & 2581 ! ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 2582 !z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) + & 2583 ! hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / & 2584 ! ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 2585 !z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) / & 2586 ! ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 2587 !z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) / & 2588 ! ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 2453 2589 2454 2590 e3t_0(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk) -
branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
r5224 r5842 17 17 !! ! (A. Coward) suppression of hel, wdj and rot options 18 18 !! 3.6 ! 2014-11 (P. Mathiot) hpg_isf: original code for ice shelf cavity 19 !! 3.6? ! 2015-11 (H. Liu) add Wetting/Drying pressure filter 19 20 !!---------------------------------------------------------------------- 20 21 … … 378 379 INTEGER, INTENT(in) :: kt ! ocean time-step index 379 380 !! 380 INTEGER :: ji, jj, jk ! dummy loop indices 381 REAL(wp) :: zcoef0, zuap, zvap, znad ! temporary scalars 381 INTEGER :: ji, jj, jk, jii, jjj ! dummy loop indices 382 REAL(wp) :: zcoef0, zuap, zvap, znad, ztmp ! temporary scalars 383 LOGICAL :: ll_tmp1, ll_tmp2, ll_tmp3 ! local logical variables 382 384 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj 385 REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy !W/D pressure filter 383 386 !!---------------------------------------------------------------------- 384 387 ! 385 388 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj ) 389 IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 386 390 ! 387 391 IF( kt == nit000 ) THEN … … 397 401 ELSE ; znad = 0._wp ! Fixed volume 398 402 ENDIF 403 404 IF(ln_wd) THEN 405 DO jj = 2, jpjm1 406 DO ji = 2, jpim1 407 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) 408 ll_tmp2 = MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) > rn_wdmin1 + rn_wdmin2 409 ll_tmp3 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) + & 410 & rn_wdmin1 + rn_wdmin2 411 412 IF(ll_tmp1.AND.ll_tmp2) THEN 413 zcpx(ji,jj) = 1.0_wp 414 wduflt(ji,jj) = 1.0_wp 415 ELSE IF(ll_tmp3) THEN 416 ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 417 zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) / & 418 & (sshn(ji+1,jj) - sshn(ji,jj))) 419 wduflt(ji,jj) = 1.0_wp 420 ELSE 421 zcpx(ji,jj) = 0._wp 422 wduflt(ji,jj) = 0.0_wp 423 END IF 424 425 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) 426 ll_tmp2 = MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) > rn_wdmin1 + rn_wdmin2 427 ll_tmp3 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) + & 428 & rn_wdmin1 + rn_wdmin2 429 430 IF(ll_tmp1.AND.ll_tmp2) THEN 431 zcpy(ji,jj) = 1.0_wp 432 wdvflt(ji,jj) = 1.0_wp 433 ELSE IF(ll_tmp3) THEN 434 ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 435 zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) / & 436 & (sshn(ji,jj+1) - sshn(ji,jj))) 437 wdvflt(ji,jj) = 1.0_wp 438 ELSE 439 zcpy(ji,jj) = 0._wp 440 wdvflt(ji,jj) = 0.0_wp 441 END IF 442 END DO 443 END DO 444 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 445 ENDIF 446 447 448 !jii=jidbg-nimpp+1;jjj=jjdbg-njmpp+1 399 449 400 450 ! Surface value … … 411 461 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) + 2._wp * znad ) & 412 462 & * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj) 463 464 465 IF(ln_wd) THEN 466 467 zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 468 zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 469 zuap = zuap * zcpx(ji,jj) 470 zvap = zvap * zcpy(ji,jj) 471 ENDIF 472 413 473 ! add to the general momentum trend 414 474 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap … … 433 493 zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 434 494 & * ( fsde3w(ji ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) 495 496 IF(ln_wd) THEN 497 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 498 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) 499 zuap = zuap * zcpx(ji,jj) 500 zvap = zvap * zcpy(ji,jj) 501 ENDIF 502 435 503 ! add to the general momentum trend 436 504 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap … … 441 509 ! 442 510 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj ) 511 IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 443 512 ! 444 513 END SUBROUTINE hpg_sco … … 719 788 REAL(wp) :: z1_10, cffu, cffx ! " " 720 789 REAL(wp) :: z1_12, cffv, cffy ! " " 790 LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables 721 791 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj 722 792 REAL(wp), POINTER, DIMENSION(:,:,:) :: dzx, dzy, dzz, dzu, dzv, dzw 723 793 REAL(wp), POINTER, DIMENSION(:,:,:) :: drhox, drhoy, drhoz, drhou, drhov, drhow 724 794 REAL(wp), POINTER, DIMENSION(:,:,:) :: rho_i, rho_j, rho_k 795 REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy !W/D pressure filter 725 796 !!---------------------------------------------------------------------- 726 797 ! … … 728 799 CALL wrk_alloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow ) 729 800 CALL wrk_alloc( jpi, jpj, jpk, rho_i, rho_j, rho_k, zhpi, zhpj ) 730 ! 801 IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 802 ! 803 !!---------------------------------------------------------------------- 804 ! 805 ! 806 IF(ln_wd) THEN 807 DO jj = 2, jpjm1 808 DO ji = 2, jpim1 809 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) & 810 & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) & 811 & > rn_wdmin1 + rn_wdmin2 812 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) +& 813 & rn_wdmin1 + rn_wdmin2 814 815 IF(ll_tmp1) THEN 816 zcpx(ji,jj) = 1.0_wp 817 ELSE IF(ll_tmp2) THEN 818 ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 819 zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /& 820 & (sshn(ji+1,jj) - sshn(ji,jj))) 821 ELSE 822 zcpx(ji,jj) = 0._wp 823 END IF 824 825 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) & 826 & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) & 827 & > rn_wdmin1 + rn_wdmin2 828 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) +& 829 & rn_wdmin1 + rn_wdmin2 830 831 IF(ll_tmp1) THEN 832 zcpy(ji,jj) = 1.0_wp 833 ELSE IF(ll_tmp2) THEN 834 ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 835 zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) /& 836 & (sshn(ji,jj+1) - sshn(ji,jj))) 837 ELSE 838 zcpy(ji,jj) = 0._wp 839 END IF 840 END DO 841 END DO 842 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 843 ENDIF 844 731 845 732 846 IF( kt == nit000 ) THEN … … 899 1013 zhpi(ji,jj,1) = ( rho_k(ji+1,jj ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) / e1u(ji,jj) 900 1014 zhpj(ji,jj,1) = ( rho_k(ji ,jj+1,1) - rho_k(ji,jj,1) - rho_j(ji,jj,1) ) / e2v(ji,jj) 1015 IF(ln_wd) THEN 1016 zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 1017 zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 1018 ENDIF 901 1019 ! add to the general momentum trend 902 1020 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) … … 918 1036 & + ( ( rho_k(ji,jj+1,jk) - rho_k(ji,jj,jk ) ) & 919 1037 & -( rho_j(ji,jj ,jk) - rho_j(ji,jj,jk-1) ) ) / e2v(ji,jj) 1038 IF(ln_wd) THEN 1039 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 1040 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) 1041 ENDIF 920 1042 ! add to the general momentum trend 921 1043 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) … … 928 1050 CALL wrk_dealloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow ) 929 1051 CALL wrk_dealloc( jpi, jpj, jpk, rho_i, rho_j, rho_k, zhpi, zhpj ) 1052 IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 930 1053 ! 931 1054 END SUBROUTINE hpg_djc … … 951 1074 !! The local variables for the correction term 952 1075 INTEGER :: jk1, jis, jid, jjs, jjd 1076 LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables 953 1077 REAL(wp) :: zuijk, zvijk, zpwes, zpwed, zpnss, zpnsd, zdeps 954 1078 REAL(wp) :: zrhdt1 … … 957 1081 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp 958 1082 REAL(wp), POINTER, DIMENSION(:,:) :: zsshu_n, zsshv_n 1083 REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy !W/D pressure filter 959 1084 !!---------------------------------------------------------------------- 960 1085 ! … … 962 1087 CALL wrk_alloc( jpi,jpj,jpk, zdept, zrhh ) 963 1088 CALL wrk_alloc( jpi,jpj, zsshu_n, zsshv_n ) 1089 IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 964 1090 ! 965 1091 IF( kt == nit000 ) THEN … … 974 1100 znad = 0.0_wp 975 1101 IF( lk_vvl ) znad = 1._wp 1102 1103 IF(ln_wd) THEN 1104 DO jj = 2, jpjm1 1105 DO ji = 2, jpim1 1106 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) & 1107 & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) & 1108 & > rn_wdmin1 + rn_wdmin2 1109 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) +& 1110 & rn_wdmin1 + rn_wdmin2 1111 1112 IF(ll_tmp1) THEN 1113 zcpx(ji,jj) = 1.0_wp 1114 ELSE IF(ll_tmp2) THEN 1115 ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 1116 zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /& 1117 & (sshn(ji+1,jj) - sshn(ji,jj))) 1118 ELSE 1119 zcpx(ji,jj) = 0._wp 1120 END IF 1121 1122 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) & 1123 & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) & 1124 & > rn_wdmin1 + rn_wdmin2 1125 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) +& 1126 & rn_wdmin1 + rn_wdmin2 1127 1128 IF(ll_tmp1.OR.ll_tmp2) THEN 1129 zcpy(ji,jj) = 1.0_wp 1130 ELSE IF(ll_tmp2) THEN 1131 ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 1132 zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) /& 1133 & (sshn(ji,jj+1) - sshn(ji,jj))) 1134 ELSE 1135 zcpy(ji,jj) = 0._wp 1136 END IF 1137 END DO 1138 END DO 1139 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 1140 ENDIF 976 1141 977 1142 ! Clean 3-D work arrays … … 1052 1217 END DO 1053 1218 END DO 1219 1220 CALL lbc_lnk (zsshu_n, 'U', 1.) 1221 CALL lbc_lnk (zsshv_n, 'V', 1.) 1054 1222 1055 1223 DO jj = 2, jpjm1 … … 1150 1318 ENDIF 1151 1319 1152 ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * & 1153 & umask(ji,jj,jk) * tmask(ji,jj,jk) * tmask(ji+1,jj,jk) 1320 IF(ln_wd) THEN 1321 zdpdx1 = zdpdx1 * zcpx(ji,jj) 1322 zdpdx2 = zdpdx2 * zcpx(ji,jj) 1323 ENDIF 1324 ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk) 1154 1325 ENDIF 1155 1326 … … 1207 1378 ENDIF 1208 1379 1209 va(ji,jj,jk) = va(ji,jj,jk) + (zdpdy1 + zdpdy2)*& 1210 & vmask(ji,jj,jk)*tmask(ji,jj,jk)*tmask(ji,jj+1,jk) 1380 IF(ln_wd) THEN 1381 zdpdy1 = zdpdy1 * zcpy(ji,jj) 1382 zdpdy2 = zdpdy2 * zcpy(ji,jj) 1383 ENDIF 1384 1385 va(ji,jj,jk) = va(ji,jj,jk) + (zdpdy1 + zdpdy2) * vmask(ji,jj,jk) 1211 1386 ENDIF 1212 1387 … … 1219 1394 CALL wrk_dealloc( jpi,jpj,jpk, zdept, zrhh ) 1220 1395 CALL wrk_dealloc( jpi,jpj, zsshu_n, zsshv_n ) 1396 IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 1221 1397 ! 1222 1398 END SUBROUTINE hpg_prj -
branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r5656 r5842 41 41 USE timing ! Timing 42 42 USE sbcapr ! surface boundary condition: atmospheric pressure 43 USE wadlmt ! wetting/drying flux limter 43 44 USE dynadv, ONLY: ln_dynadv_vec 44 45 #if defined key_agrif … … 142 143 LOGICAL :: ll_fw_start ! if true, forward integration 143 144 LOGICAL :: ll_init ! if true, special startup of 2d equations 145 LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables used in W/D 144 146 INTEGER :: ji, jj, jk, jn ! dummy loop indices 145 147 INTEGER :: ikbu, ikbv, noffset ! local integers … … 157 159 REAL(wp), POINTER, DIMENSION(:,:) :: zsshu_a, zsshv_a 158 160 REAL(wp), POINTER, DIMENSION(:,:) :: zhf 161 REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy ! Wetting/Dying gravity filter coef. 162 REAL(wp), POINTER, DIMENSION(:,:) :: wduflt1, wdvflt1 ! Wetting/Dying velocity filter coef. 159 163 !!---------------------------------------------------------------------- 164 165 ! tempoary debugging integer 166 INTEGER :: jidbg, jjdbg 167 jidbg = 163; jjdbg = 179 160 168 ! 161 169 IF( nn_timing == 1 ) CALL timing_start('dyn_spg_ts') … … 168 176 CALL wrk_alloc( jpi, jpj, zsshu_a, zsshv_a ) 169 177 CALL wrk_alloc( jpi, jpj, zhf ) 178 IF(ln_wd) CALL wrk_alloc( jpi, jpj, zcpx, zcpy, wduflt1, wdvflt1 ) 170 179 ! 171 180 ! !* Local constant initialization … … 379 388 ! ! ---------------------------------------------------- 380 389 IF( lk_vvl ) THEN ! Variable volume : remove surface pressure gradient 381 DO jj = 2, jpjm1 382 DO ji = fs_2, fs_jpim1 ! vector opt. 383 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) / e1u(ji,jj) 384 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) / e2v(ji,jj) 385 END DO 386 END DO 387 ENDIF 388 390 IF(ln_wd) THEN ! Calculating and applying W/D gravity filters 391 DO jj = 2, jpjm1 392 DO ji = 2, jpim1 393 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) & 394 & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) & 395 & > rn_wdmin1 + rn_wdmin2 396 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) + & 397 & rn_wdmin1 + rn_wdmin2 398 IF(ll_tmp1) THEN 399 zcpx(ji,jj) = 1.0_wp 400 wduflt1(ji,jj) = 1.0_wp 401 ELSE IF(ll_tmp2) THEN 402 ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 403 zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /& 404 & (sshn(ji+1,jj) - sshn(ji,jj))) 405 wduflt1(ji,jj) = 1.0_wp 406 ELSE 407 zcpx(ji,jj) = 0._wp 408 wduflt1(ji,jj) = 0.0_wp 409 END IF 410 411 !for w/d debugging, delete it when finished 412 ! zcpx(ji,jj) = 0._wp 413 ! end debugging part 414 415 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) & 416 & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) & 417 & > rn_wdmin1 + rn_wdmin2 418 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) + & 419 & rn_wdmin1 + rn_wdmin2 420 IF(ll_tmp1) THEN 421 zcpy(ji,jj) = 1.0_wp 422 wdvflt1(ji,jj) = 1.0_wp 423 ELSE IF(ll_tmp2) THEN 424 ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 425 zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) /& 426 & (sshn(ji,jj+1) - sshn(ji,jj))) 427 wdvflt1(ji,jj) = 1.0_wp 428 ELSE 429 zcpy(ji,jj) = 0._wp 430 wdvflt1(ji,jj) = 0.0_wp 431 END IF 432 433 !for w/d debugging, delete it when finished 434 ! zcpy(ji,jj) = 0._wp 435 ! end debugging part 436 END DO 437 END DO 438 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 439 ENDIF 440 441 IF(ln_wd) THEN ! Calculating and applying W/D gravity filters 442 DO jj = 2, jpjm1 443 DO ji = 2, jpim1 444 zu_trd(ji,jj) = (zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) / & 445 & e1u(ji,jj)) * zcpx(ji,jj) 446 zv_trd(ji,jj) = (zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) / & 447 & e2v(ji,jj)) * zcpy(ji,jj) 448 IF(ji+nimpp-1 == jidbg .and. jj+njmpp-1 == jjdbg) THEN 449 write(333,*) 'in spt_ts: zutrdx, zutrdy, hpgx, hpgy, zcpx, zcpy' 450 write(333,*) zu_trd(ji,jj), zv_trd(ji,jj), & 451 & -grav*(sshn(ji+1,jj)-sshn(ji,jj))/e1u(ji,jj)*zcpx(ji,jj), & 452 & -grav*(sshn(ji,jj+1)-sshn(ji,jj))/e2v(ji,jj)*zcpy(ji,jj), & 453 & zcpx(ji,jj), zcpy(ji,jj) 454 write(334,*) 'in spt_ts: sshn_ij, sshn_ip1j, sshn_ijp1, bathy_ij, bathy_ip1j, bathy_ijp1' 455 write(334,*) sshn(ji,jj), sshn(ji+1,jj), sshn(ji,jj+1), bathy(ji,jj), bathy(ji+1,jj), bathy(ji,jj+1) 456 END IF 457 END DO 458 END DO 459 ELSE 460 DO jj = 2, jpjm1 461 DO ji = fs_2, fs_jpim1 ! vector opt. 462 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) / e1u(ji,jj) 463 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) / e2v(ji,jj) 464 END DO 465 END DO 466 END IF 467 468 ENDIF 469 470 IF(ln_wd) THEN 471 DO jj = 2, jpjm1 ! Remove coriolis term (and possibly spg) from barotropic trend 472 DO ji = fs_2, fs_jpim1 473 zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * umask(ji,jj,1) * wduflt1(ji,jj) 474 zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * vmask(ji,jj,1) * wdvflt1(ji,jj) 475 END DO 476 END DO 477 ELSE 389 478 DO jj = 2, jpjm1 ! Remove coriolis term (and possibly spg) from barotropic trend 390 479 DO ji = fs_2, fs_jpim1 … … 392 481 zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * vmask(ji,jj,1) 393 482 END DO 394 END DO 483 END DO 484 END IF 395 485 ! 396 486 ! ! Add bottom stress contribution from baroclinic velocities: … … 416 506 ! 417 507 ! Note that the "unclipped" bottom friction parameter is used even with explicit drag 418 zu_frc(:,:) = zu_frc(:,:) + hur(:,:) * bfrua(:,:) * zwx(:,:) 419 zv_frc(:,:) = zv_frc(:,:) + hvr(:,:) * bfrva(:,:) * zwy(:,:) 420 ! 508 IF(ln_wd) THEN 509 zu_frc(:,:) = zu_frc(:,:) + MAX(hur(:,:) * bfrua(:,:),-1._wp / rdtbt) * zwx(:,:) 510 zv_frc(:,:) = zv_frc(:,:) + MAX(hvr(:,:) * bfrva(:,:),-1._wp / rdtbt) * zwy(:,:) 511 ELSE 512 zu_frc(:,:) = zu_frc(:,:) + hur(:,:) * bfrua(:,:) * zwx(:,:) 513 zv_frc(:,:) = zv_frc(:,:) + hvr(:,:) * bfrva(:,:) * zwy(:,:) 514 END IF 515 ! 421 516 IF (ln_bt_fw) THEN ! Add wind forcing 422 517 zu_frc(:,:) = zu_frc(:,:) + zraur * utau(:,:) * hur(:,:) … … 487 582 vb_e (:,:) = 0._wp 488 583 ENDIF 489 ! 490 IF (ln_bt_fw) THEN ! FORWARD integration: start from NOW fields 491 sshn_e(:,:) = sshn (:,:) 492 zun_e (:,:) = un_b (:,:) 584 585 IF(ln_wd) THEN !preserve the positivity of water depth 586 !ssh[b,n,a] should have already been processed for this 587 sshbb_e(:,:) = MAX(sshbb_e(:,:), rn_wdmin1 - bathy(:,:)) 588 sshb_e(:,:) = MAX(sshb_e(:,:) , rn_wdmin1 - bathy(:,:)) 589 ENDIF 590 ! 591 IF (ln_bt_fw) THEN ! FORWARD integration: start from NOW fields 592 sshn_e(:,:) = sshn (:,:) 593 zun_e (:,:) = un_b (:,:) 493 594 zvn_e (:,:) = vn_b (:,:) 494 595 ! … … 561 662 zhup2_e (:,:) = hu_0(:,:) + zwx(:,:) ! Ocean depth at U- and V-points 562 663 zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:) 664 IF(ln_wd) THEN 665 zhup2_e(:,:) = MAX(zhup2_e (:,:), rn_wdmin1) 666 zhvp2_e(:,:) = MAX(zhvp2_e (:,:), rn_wdmin1) 667 END IF 563 668 ELSE 564 669 zhup2_e (:,:) = hu(:,:) … … 599 704 ENDIF 600 705 #endif 706 IF(ln_wd) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt) 601 707 ! 602 708 ! Sum over sub-time-steps to compute advective velocities … … 613 719 END DO 614 720 ssha_e(:,:) = ( sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) ) ) * tmask(:,:,1) 721 IF(ln_wd) ssha_e(:,:) = MAX(ssha_e(:,:), rn_wdmin1 - bathy(:,:)) 615 722 CALL lbc_lnk( ssha_e, 'T', 1._wp ) 616 723 … … 660 767 & + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:) 661 768 769 IF(ln_wd) THEN ! Calculating and applying W/D gravity filters 770 DO jj = 2, jpjm1 771 DO ji = 2, jpim1 772 ll_tmp1 = MIN(zsshp2_e(ji,jj), zsshp2_e(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj))& 773 & .and. MAX(zsshp2_e(ji,jj) + bathy(ji,jj), zsshp2_e(ji+1,jj) + bathy(ji+1,jj)) & 774 & > rn_wdmin1 + rn_wdmin2 775 ll_tmp2 = MAX(zsshp2_e(ji,jj), zsshp2_e(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) + & 776 & rn_wdmin1 + rn_wdmin2 777 IF(ll_tmp1) THEN 778 zcpx(ji,jj) = 1.0_wp 779 wduflt1(ji,jj) = 1.0_wp 780 ELSE IF(ll_tmp2) THEN 781 ! no worries about zsshp2_e(ji+1,jj)-zsshp2_e(ji,jj) = 0, it won't happen ! here 782 zcpx(ji,jj) = ABS((zsshp2_e(ji+1,jj) + bathy(ji+1,jj) - zsshp2_e(ji,jj) - bathy(ji,jj)) /& 783 & (zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj))) 784 wduflt1(ji,jj) = 1.0_wp 785 ELSE 786 zcpx(ji,jj) = 0._wp 787 wduflt1(ji,jj) = 0.0_wp 788 END IF 789 790 !for w/d debugging, delete it when finished 791 ! zcpx(ji,jj) = 0._wp 792 ! end debugging part 793 794 ll_tmp1 = MIN(zsshp2_e(ji,jj), zsshp2_e(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1))& 795 & .and. MAX(zsshp2_e(ji,jj) + bathy(ji,jj), zsshp2_e(ji,jj+1) + bathy(ji,jj+1)) & 796 & > rn_wdmin1 + rn_wdmin2 797 ll_tmp2 = MAX(zsshp2_e(ji,jj), zsshp2_e(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) + & 798 & rn_wdmin1 + rn_wdmin2 799 IF(ll_tmp1) THEN 800 zcpy(ji,jj) = 1.0_wp 801 wdvflt1(ji,jj) = 1.0_wp 802 ELSE IF(ll_tmp2) THEN 803 ! no worries about zsshp2_e(ji,jj+1)-zsshp2_e(ji,jj) = 0, it won't happen ! here 804 zcpy(ji,jj) = ABS((zsshp2_e(ji,jj+1) + bathy(ji,jj+1) - zsshp2_e(ji,jj) - bathy(ji,jj)) /& 805 & (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj))) 806 wdvflt1(ji,jj) = 1.0_wp 807 ELSE 808 zcpy(ji,jj) = 0._wp 809 wdvflt1(ji,jj) = 0.0_wp 810 END IF 811 812 !for w/d debugging, delete it when finished 813 ! zcpy(ji,jj) = 0._wp 814 ! end debugging part 815 END DO 816 END DO 817 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 818 ENDIF 662 819 ! 663 820 ! Compute associated depths at U and V points: … … 676 833 END DO 677 834 END DO 835 836 IF(ln_wd) THEN 837 zhust_e(:,:) = MAX(zhust_e (:,:), rn_wdmin1 ) 838 zhvst_e(:,:) = MAX(zhvst_e (:,:), rn_wdmin1 ) 839 END IF 840 !shall we call lbc_lnk for zhu(v)st_e() here? 841 678 842 ENDIF 679 843 ! … … 742 906 ! 743 907 ! Surface pressure trend: 744 DO jj = 2, jpjm1 745 DO ji = fs_2, fs_jpim1 ! vector opt. 746 ! Add surface pressure gradient 747 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj) 748 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj) 749 zwx(ji,jj) = zu_spg 750 zwy(ji,jj) = zv_spg 751 END DO 752 END DO 908 909 IF(ln_wd) THEN 910 DO jj = 2, jpjm1 911 DO ji = 2, jpim1 912 ! Add surface pressure gradient 913 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj) 914 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj) 915 ! for W/D debugging 916 !zwx(ji,jj) = zu_spg * zcpx(ji,jj) * 0._wp 917 !zwy(ji,jj) = zv_spg * zcpy(ji,jj) * 0._wp 918 IF(ji+nimpp-1 == jidbg .and. jj+njmpp-1 == jjdbg) THEN 919 write(444,*) 'in spt_ts_intg: zu_spg, zv_spg, zcpx, zcpy' 920 write(444,*) zu_spg, zv_spg, zcpx(ji,jj), zcpy(ji,jj) 921 write(445,*) 'in spt_ts_intg: sshn_ij, sshn_ip1j, sshn_ijp1, bathy_ij, bathy_ip1j, bathy_ijp1' 922 write(445,*) zsshp2_e(ji,jj), zsshp2_e(ji+1,jj), zsshp2_e(ji,jj+1), & 923 & bathy(ji,jj), bathy(ji+1,jj), bathy(ji,jj+1) 924 END IF 925 zwx(ji,jj) = zu_spg * zcpx(ji,jj) 926 zwy(ji,jj) = zv_spg * zcpy(ji,jj) 927 END DO 928 END DO 929 ELSE 930 DO jj = 2, jpjm1 931 DO ji = fs_2, fs_jpim1 ! vector opt. 932 ! Add surface pressure gradient 933 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj) 934 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj) 935 zwx(ji,jj) = zu_spg 936 zwy(ji,jj) = zv_spg 937 END DO 938 END DO 939 END IF 940 753 941 ! 754 942 ! Set next velocities: … … 774 962 DO ji = fs_2, fs_jpim1 ! vector opt. 775 963 776 zhura = umask(ji,jj,1)/(hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - umask(ji,jj,1)) 777 zhvra = vmask(ji,jj,1)/(hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - vmask(ji,jj,1)) 778 779 ua_e(ji,jj) = ( hu_e(ji,jj) * zun_e(ji,jj) & 780 & + rdtbt * ( zhust_e(ji,jj) * zwx(ji,jj) & 964 IF(ln_wd) THEN 965 zhura = MAX(hu_0(ji,jj) + zsshu_a(ji,jj), rn_wdmin1) 966 zhvra = MAX(hv_0(ji,jj) + zsshv_a(ji,jj), rn_wdmin1) 967 ELSE 968 zhura = hu_0(ji,jj) + zsshu_a(ji,jj) 969 zhvra = hv_0(ji,jj) + zsshv_a(ji,jj) 970 END IF 971 972 zhura = umask(ji,jj,1)/(zhura + 1._wp - umask(ji,jj,1)) 973 zhvra = vmask(ji,jj,1)/(zhvra + 1._wp - vmask(ji,jj,1)) 974 975 ua_e(ji,jj) = ( hu_e(ji,jj) * zun_e(ji,jj) & 976 & + rdtbt * ( zhust_e(ji,jj) * zwx(ji,jj) & 781 977 & + zhup2_e(ji,jj) * zu_trd(ji,jj) & 782 978 & + hu(ji,jj) * zu_frc(ji,jj) ) & … … 793 989 ! 794 990 IF( lk_vvl ) THEN !* Update ocean depth (variable volume case only) 795 ! ! ---------------------------------------------- 796 hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 797 hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 991 ! ! ---------------------------------------------- 992 IF(ln_wd) THEN 993 hu_e (:,:) = MAX(hu_0(:,:) + zsshu_a(:,:), rn_wdmin1) 994 hv_e (:,:) = MAX(hv_0(:,:) + zsshv_a(:,:), rn_wdmin1) 995 ELSE 996 hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 997 hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 998 END IF 999 798 1000 hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1._wp - umask(:,:,1) ) 799 1001 hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1._wp - vmask(:,:,1) ) … … 925 1127 CALL wrk_dealloc( jpi, jpj, zsshu_a, zsshv_a ) 926 1128 CALL wrk_dealloc( jpi, jpj, zhf ) 1129 IF(ln_wd) CALL wrk_dealloc( jpi, jpj, zcpx, zcpy, wduflt1, wdvflt1 ) 927 1130 ! 928 1131 IF( nn_timing == 1 ) CALL timing_stop('dyn_spg_ts') -
branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
r5656 r5842 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.? ! 2015-11 (H. Liu) add wetting and drying 11 12 !!---------------------------------------------------------------------- 12 13 … … 38 39 USE wrk_nemo ! Memory Allocation 39 40 USE timing ! Timing 41 USE wadlmt ! Wetting/Drying flux limting 40 42 41 43 IMPLICIT NONE … … 90 92 ENDIF 91 93 ! 92 CALL div_cur( kt ) ! Horizontal divergence & Relative vorticity93 !94 94 z2dt = 2._wp * rdt ! set time step size (Euler/Leapfrog) 95 95 IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt 96 ! 97 z1_rau0 = 0.5_wp * r1_rau0 98 ! 99 IF(ln_wd) CALL wad_lmt(sshb, z1_rau0 * (emp_b(:,:) + emp(:,:)), z2dt) 100 101 CALL div_cur( kt ) ! Horizontal divergence & Relative vorticity 102 ! 96 103 97 104 ! !------------------------------! … … 106 113 ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations. 107 114 ! 108 z1_rau0 = 0.5_wp * r1_rau0109 115 ssha(:,:) = ( sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * ssmask(:,:) 110 116 -
branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
r5656 r5842 30 30 !! 3.4 ! 2011-11 (C. Harris) decomposition changes for running with CICE 31 31 !! ! 2012-05 (C. Calone, J. Simeon, G. Madec, C. Ethe) Add grid coarsening 32 !! 3.6.?! 2015-11 (H. Liu) Add Wetting and Drying 32 33 !!---------------------------------------------------------------------- 33 34 … … 232 233 NAMELIST/namcfg/ cp_cfg, cp_cfz, jp_cfg, jpidta, jpjdta, jpkdta, jpiglo, jpjglo, & 233 234 & jpizoom, jpjzoom, jperio, ln_use_jattr 235 NAMELIST/namwad/ ln_wd, rn_wdmin1, rn_wdmin2, rn_wdld, nn_wdit 234 236 !!---------------------------------------------------------------------- 235 237 ! … … 257 259 READ ( numnam_cfg, namcfg, IOSTAT = ios, ERR = 904 ) 258 260 904 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcfg in configuration namelist', .TRUE. ) 261 262 REWIND( numnam_ref ) ! Namelist namwad in reference namelist : Parameters for Wetting/Drying 263 READ ( numnam_ref, namwad, IOSTAT = ios, ERR = 905) 264 905 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namwad in reference namelist', .TRUE. ) 265 266 REWIND( numnam_cfg ) ! Namelist namwad in configuration namelist : Parameters for Wetting/Drying 267 READ ( numnam_cfg, namwad, IOSTAT = ios, ERR = 906) 268 906 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namwad in configuration namelist', .TRUE. ) 269 270 259 271 260 272 ! Force values for AGRIF zoom (cf. agrif_user.F90) … … 315 327 WRITE( numond, namctl ) 316 328 WRITE( numond, namcfg ) 329 WRITE( numond, namwad ) 317 330 ENDIF 318 331
Note: See TracChangeset
for help on using the changeset viewer.