Changeset 5954 for branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO
- Timestamp:
- 2015-11-30T15:04:08+01:00 (9 years ago)
- Location:
- branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC
- Files:
-
- 27 edited
- 1 copied
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90
r4990 r5954 509 509 ALLOCATE( dta_global(jpbdtau, 1, jpk) ) 510 510 IF ( icount>0 ) ALLOCATE( dta_global2(jpbdtas, nrimmax, jpk) ) 511 512 nbidta(:,:,:)=0._wp 513 nbjdta(:,:,:)=0._wp 514 nbrdta(:,:,:)=0._wp 511 515 ! 512 516 ENDIF -
branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90
r5123 r5954 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,limiters,and masks for Wetting and Drying 276 !! --------------------------------------------------------------------- 277 278 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wdmask !: u- and v- limiter 279 280 LOGICAL, PUBLIC, SAVE :: ln_wd !: key to turn on/off wetting/drying (T: on, F: off) 281 REAL(wp), PUBLIC, SAVE :: rn_wdmin1 !: minimum water depth on dried cells 282 REAL(wp), PUBLIC, SAVE :: rn_wdmin2 !: tolerrance of minimum water depth on dried cells 283 REAL(wp), PUBLIC, SAVE :: rn_wdld !: land elevation below which wetting/drying will be considered 284 INTEGER , PUBLIC, SAVE :: nn_wdit !: maximum number of iteration for W/D limiter 272 285 273 286 !!---------------------------------------------------------------------- … … 408 421 #endif 409 422 ! 423 IF(ln_wd) & 424 ALLOCATE( wdmask(jpi,jpj), STAT=ierr(12) ) 425 ! 410 426 dom_oce_alloc = MAXVAL(ierr) 411 427 ! -
branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90
r5363 r5954 86 86 CALL dom_zgr ! Vertical mesh and bathymetry 87 87 CALL dom_msk ! Masks 88 IF( ln_sco ) CALL dom_stiff ! Maximum stiffness ratio/hydrostatic consistency 88 IF( ln_sco.AND.(.NOT.ln_wd) ) & 89 & CALL dom_stiff ! Maximum stiffness ratio/hydrostatic consistency 89 90 ! 90 91 ht_0(:,:) = 0.0_wp ! Reference ocean depth at T-points … … 137 138 NAMELIST/namrun/ cn_ocerst_indir, cn_ocerst_outdir, nn_stocklist, ln_rst_list, & 138 139 & nn_no , cn_exp , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl, & 140 & ln_rstdate & 139 141 & nn_it000, nn_itend , nn_date0 , nn_leapy , nn_istate , nn_stock , & 140 142 & nn_write, ln_dimgnnn, ln_mskland , ln_cfmeta , ln_clobber, nn_chunksz, nn_euler … … 174 176 WRITE(numout,*) ' restart output directory cn_ocerst_outdir= ', cn_ocerst_outdir 175 177 WRITE(numout,*) ' restart logical ln_rstart = ', ln_rstart 178 WRITE(numout,*) ' datestamping of restarts ln_rstdate = ', ln_rstdate 179 WRITE(numout,*) ' restart input directory cn_ocerst_indir= ', cn_ocerst_indir 180 WRITE(numout,*) ' restart output directory cn_ocerst_outdir= ', cn_ocerst_outdir 176 181 WRITE(numout,*) ' start with forward time step nn_euler = ', nn_euler 177 182 WRITE(numout,*) ' control of time step nn_rstctl = ', nn_rstctl -
branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
r5506 r5954 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 ! ============================================================================= … … 717 729 DO jj = 1, jpjm1 718 730 DO ji = 1, fs_jpim1 ! vector opt. 719 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * r1_e12u(ji,jj) & 731 ! pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * r1_e12u(ji,jj) & 732 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * r1_e12u(ji,jj) & 720 733 & * ( e12t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) & 721 734 & + e12t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) ) … … 735 748 DO jj = 1, jpjm1 736 749 DO ji = 1, fs_jpim1 ! vector opt. 737 pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) * r1_e12v(ji,jj) & 738 & * ( e12t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) & 739 & + e12t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) ) 750 ! pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) * r1_e12v(ji,jj) & 751 pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) * r1_e12v(ji,jj) & 752 & * ( e12t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) & 753 & + e12t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) ) 740 754 END DO 741 755 END DO … … 753 767 DO jj = 1, jpjm1 754 768 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) & 769 ! pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) * r1_e12f(ji,jj) & 770 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) & 771 & * r1_e12f(ji,jj) & 756 772 & * ( e12u(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3u_0(ji,jj ,jk) ) & 757 773 & + e12u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) ) … … 817 833 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 818 834 !! * Local declarations 819 INTEGER :: j k835 INTEGER :: ji, jj, jk 820 836 INTEGER :: id1, id2, id3, id4, id5 ! local integers 821 837 !!---------------------------------------------------------------------- … … 905 921 fse3t_n(:,:,:) = e3t_0(:,:,:) 906 922 sshn(:,:) = 0.0_wp 923 924 IF(ln_wd) THEN 925 DO jj = 1, jpj 926 DO ji = 1, jpi 927 !IF(e3t_0(ji,jj,1) < 0._wp) THEN 928 !IF(mbathy(ji,jj) == 2 .AND. e3t_0(ji,jj,1) <= 0.5_wp * rn_wdmin1) THEN 929 IF( e3t_0(ji,jj,1) <= 0.5_wp * rn_wdmin1) THEN 930 fse3t_b(ji,jj,:) = 0.5_wp * rn_wdmin1 931 fse3t_n(ji,jj,:) = 0.5_wp * rn_wdmin1 932 fse3t_a(ji,jj,:) = 0.5_wp * rn_wdmin1 933 sshb(ji,jj) = rn_wdmin1 - bathy(ji,jj) 934 sshn(ji,jj) = rn_wdmin1 - bathy(ji,jj) 935 ssha(ji,jj) = rn_wdmin1 - bathy(ji,jj) 936 ENDIF 937 ENDDO 938 ENDDO 939 END IF 940 907 941 IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN 908 942 tilde_e3t_b(:,:,:) = 0.0_wp -
branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
r5506 r5954 17 17 !! 3.4 ! 2012-08 (J. Siddorn) added Siddorn and Furner stretching function 18 18 !! 3.4 ! 2012-12 (R. Bourdalle-Badie and G. Reffray) modify C1D case 19 !! 3.6 ! 2014-09 (H. Liu) Modifications for Wetting/Drying 19 20 !! 3.6 ! 2014-11 (P. Mathiot and C. Harris) add ice shelf capabilitye 20 21 !!---------------------------------------------------------------------- … … 316 317 zrefdep = 10._wp - 0.1_wp * MINVAL( e3w_1d ) ! ref. depth with tolerance (10% of minimum layer thickness) 317 318 nlb10 = MINLOC( gdepw_1d, mask = gdepw_1d > zrefdep, dim = 1 ) ! shallowest W level Below ~10m 319 nlb10 = MAX(nlb10, 2) ! prevent nla10 = 0 318 320 nla10 = nlb10 - 1 ! deepest W level Above ~10m 319 321 !!gm end bug … … 1796 1798 REAL(wp) :: zrmax, ztaper ! temporary scalars 1797 1799 REAL(wp) :: zrfact 1800 REAL(wp) :: zsmth 1798 1801 ! 1799 1802 REAL(wp), POINTER, DIMENSION(:,: ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2 … … 1851 1854 bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) ) 1852 1855 1853 DO jj = 1, jpj 1854 DO ji = 1, jpi 1855 IF( bathy(ji,jj) > 0._wp ) bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) ) 1856 END DO 1857 END DO 1856 IF (.NOT. ln_wd) THEN 1857 DO jj = 1, jpj 1858 DO ji = 1, jpi 1859 IF( bathy(ji,jj) > 0._wp ) bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) ) 1860 END DO 1861 END DO 1862 ENDIF 1858 1863 ! ! ============================= 1859 1864 ! ! Define the envelop bathymetry (hbatt) … … 1862 1867 zenv(:,:) = bathy(:,:) 1863 1868 ! 1864 ! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing 1865 DO jj = 1, jpj 1866 DO ji = 1, jpi 1867 IF( bathy(ji,jj) == 0._wp ) THEN 1868 iip1 = MIN( ji+1, jpi ) 1869 ijp1 = MIN( jj+1, jpj ) 1870 iim1 = MAX( ji-1, 1 ) 1871 ijm1 = MAX( jj-1, 1 ) 1872 IF( (bathy(iip1,jj) + bathy(iim1,jj) + bathy(ji,ijp1) + bathy(ji,ijm1) + & 1873 & bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0._wp ) THEN 1874 zenv(ji,jj) = rn_sbot_min 1869 IF (.NOT. ln_wd) THEN 1870 ! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing 1871 DO jj = 1, jpj 1872 DO ji = 1, jpi 1873 IF( bathy(ji,jj) == 0._wp ) THEN 1874 iip1 = MIN( ji+1, jpi ) 1875 ijp1 = MIN( jj+1, jpj ) 1876 iim1 = MAX( ji-1, 1 ) 1877 ijm1 = MAX( jj-1, 1 ) 1878 IF( (bathy(iip1,jj) + bathy(iim1,jj) + bathy(ji,ijp1) + bathy(ji,ijm1) + & 1879 & bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0._wp ) THEN 1880 zenv(ji,jj) = rn_sbot_min 1881 ENDIF 1875 1882 ENDIF 1876 END IF1877 1878 END DO1883 END DO 1884 END DO 1885 END IF 1879 1886 ! apply lateral boundary condition CAUTION: keep the value when the lbc field is zero 1880 1887 CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' ) … … 1916 1923 zri(:,:) = 0._wp 1917 1924 zrj(:,:) = 0._wp 1925 zsmth = 0._wp 1918 1926 DO jj = 1, nlcj 1919 1927 DO ji = 1, nlci 1920 1928 iip1 = MIN( ji+1, nlci ) ! force zri = 0 on last line (ji=ncli+1 to jpi) 1921 1929 ijp1 = MIN( jj+1, nlcj ) ! force zrj = 0 on last raw (jj=nclj+1 to jpj) 1922 IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(iip1,jj) > 0._wp)) THEN1930 IF( (zenv(ji,jj) > zsmth) .AND. (zenv(iip1,jj) > zsmth)) THEN 1923 1931 zri(ji,jj) = ( zenv(iip1,jj ) - zenv(ji,jj) ) / ( zenv(iip1,jj ) + zenv(ji,jj) ) 1924 1932 END IF 1925 IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(ji,ijp1) > 0._wp)) THEN1933 IF( (zenv(ji,jj) > zsmth) .AND. (zenv(ji,ijp1) > zsmth)) THEN 1926 1934 zrj(ji,jj) = ( zenv(ji ,ijp1) - zenv(ji,jj) ) / ( zenv(ji ,ijp1) + zenv(ji,jj) ) 1927 1935 END IF … … 1946 1954 END DO ! End loop ! 1947 1955 ! ! ================ ! 1948 DO jj = 1, jpj 1949 DO ji = 1, jpi 1950 zenv(ji,jj) = MAX( zenv(ji,jj), rn_sbot_min ) ! set all points to avoid undefined scale value warnings 1951 END DO 1952 END DO 1956 IF(ln_wd) THEN 1957 DO jj = 1, jpj 1958 DO ji = 1, jpi 1959 zenv(ji,jj) = MAX( zenv(ji,jj), -rn_wdld ) ! fill out land bathy data 1960 END DO 1961 END DO 1962 ELSE 1963 DO jj = 1, jpj 1964 DO ji = 1, jpi 1965 zenv(ji,jj) = MAX( zenv(ji,jj), rn_sbot_min ) ! set all points to avoid undefined scale value warnings 1966 END DO 1967 END DO 1968 END IF 1953 1969 ! 1954 1970 ! Envelope bathymetry saved in hbatt … … 1980 1996 IF(lwp) THEN 1981 1997 WRITE(numout,*) 1982 WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min 1998 IF (.NOT. ln_wd ) THEN 1999 WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min 2000 ELSE 2001 WRITE(numout,*) ' zgr_sco: minimum positive depth of the envelop topography set to : ', rn_sbot_min 2002 WRITE(numout,*) ' zgr_sco: minimum negative depth of the envelop topography set to : ', -rn_wdld 2003 ENDIF 1983 2004 ENDIF 1984 2005 hbatu(:,:) = rn_sbot_min … … 1993 2014 END DO 1994 2015 END DO 2016 2017 IF(ln_wd) THEN !avoid the zero depth on T- (U-,V-,F-) points 2018 DO jj = 1, jpj 2019 DO ji = 1, jpi 2020 IF(ABS(hbatt(ji,jj)) < rn_wdmin1) & 2021 & hbatt(ji,jj) = SIGN(1._wp, hbatt(ji,jj)) * rn_wdmin1 2022 IF(ABS(hbatu(ji,jj)) < rn_wdmin1) & 2023 & hbatu(ji,jj) = SIGN(1._wp, hbatu(ji,jj)) * rn_wdmin1 2024 IF(ABS(hbatv(ji,jj)) < rn_wdmin1) & 2025 & hbatv(ji,jj) = SIGN(1._wp, hbatv(ji,jj)) * rn_wdmin1 2026 IF(ABS(hbatf(ji,jj)) < rn_wdmin1) & 2027 & hbatf(ji,jj) = SIGN(1._wp, hbatf(ji,jj)) * rn_wdmin1 2028 END DO 2029 END DO 2030 END IF 1995 2031 ! 1996 2032 ! Apply lateral boundary condition … … 2000 2036 DO ji = 1, jpi 2001 2037 IF( hbatu(ji,jj) == 0._wp ) THEN 2038 !No worries about the following line when ln_wd == .true. 2002 2039 IF( zhbat(ji,jj) == 0._wp ) hbatu(ji,jj) = rn_sbot_min 2003 2040 IF( zhbat(ji,jj) /= 0._wp ) hbatu(ji,jj) = zhbat(ji,jj) … … 2025 2062 2026 2063 !!bug: key_helsinki a verifer 2027 hift(:,:) = MIN( hift(:,:), hbatt(:,:) ) 2028 hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) ) 2029 hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) ) 2030 hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) ) 2064 2065 IF (.NOT. ln_wd) THEN 2066 hift(:,:) = MIN( hift(:,:), hbatt(:,:) ) 2067 hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) ) 2068 hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) ) 2069 hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) ) 2070 END IF 2031 2071 2032 2072 IF( nprint == 1 .AND. lwp ) THEN … … 2070 2110 CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 2071 2111 2072 fsdepw(:,:,:) = gdepw_0 (:,:,:) 2073 fsde3w(:,:,:) = gdep3w_0(:,:,:) 2074 ! 2075 where (e3t_0 (:,:,:).eq.0.0) e3t_0(:,:,:) = 1.0 2076 where (e3u_0 (:,:,:).eq.0.0) e3u_0(:,:,:) = 1.0 2077 where (e3v_0 (:,:,:).eq.0.0) e3v_0(:,:,:) = 1.0 2078 where (e3f_0 (:,:,:).eq.0.0) e3f_0(:,:,:) = 1.0 2079 where (e3w_0 (:,:,:).eq.0.0) e3w_0(:,:,:) = 1.0 2080 where (e3uw_0 (:,:,:).eq.0.0) e3uw_0(:,:,:) = 1.0 2081 where (e3vw_0 (:,:,:).eq.0.0) e3vw_0(:,:,:) = 1.0 2112 !fsdepw(:,:,:) = gdepw_0 (:,:,:) 2113 !fsde3w(:,:,:) = gdep3w_0(:,:,:) 2114 ! 2115 IF (.NOT. ln_wd) THEN 2116 where (e3t_0 (:,:,:).eq.0.0) e3t_0(:,:,:) = 1.0 2117 where (e3u_0 (:,:,:).eq.0.0) e3u_0(:,:,:) = 1.0 2118 where (e3v_0 (:,:,:).eq.0.0) e3v_0(:,:,:) = 1.0 2119 where (e3f_0 (:,:,:).eq.0.0) e3f_0(:,:,:) = 1.0 2120 where (e3w_0 (:,:,:).eq.0.0) e3w_0(:,:,:) = 1.0 2121 where (e3uw_0 (:,:,:).eq.0.0) e3uw_0(:,:,:) = 1.0 2122 where (e3vw_0 (:,:,:).eq.0.0) e3vw_0(:,:,:) = 1.0 2123 END IF 2082 2124 2083 2125 #if defined key_agrif … … 2119 2161 DO ji = 1, jpi 2120 2162 DO jk = 1, jpkm1 2163 #if key_surge 2164 IF( scobot(ji,jj) >= fsdept(ji,jj,jk) ) mbathy(ji,jj) = MAX( 1, jk ) 2165 #else 2121 2166 IF( scobot(ji,jj) >= fsdept(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk ) 2122 END DO 2123 IF( scobot(ji,jj) == 0._wp ) mbathy(ji,jj) = 0 2167 #endif 2168 END DO 2169 IF(ln_wd) THEN 2170 IF( scobot(ji,jj) <= -(rn_wdld - rn_wdmin2)) THEN 2171 mbathy(ji,jj) = 0 2172 ELSEIF(scobot(ji,jj) <= rn_wdmin1) THEN 2173 #if key_surge 2174 mbathy(ji,jj) = 1 2175 #else 2176 mbathy(ji,jj) = 2 2177 #endif 2178 ENDIF 2179 ELSE 2180 IF( scobot(ji,jj) == 0._wp ) mbathy(ji,jj) = 0 2181 END IF 2124 2182 END DO 2125 2183 END DO … … 2179 2237 DO jj = 1, jpj 2180 2238 2181 IF( hbatt(ji,jj) > 0._wp) THEN2239 IF( bathy(ji,jj) > 0._wp) THEN 2182 2240 DO jk = 1, mbathy(ji,jj) 2183 2241 ! check coordinate is monotonically increasing 2184 IF (fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN 2185 WRITE(ctmp1,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk 2186 WRITE(numout,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk 2242 !RF test... 2243 !IF (fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN 2244 ! WRITE(ctmp1,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk 2245 ! WRITE(numout,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk 2246 IF (fse3w(ji,jj,jk) < 0._wp .OR. fse3t(ji,jj,jk) < 0._wp ) THEN 2247 WRITE(ctmp1,*) 'ERROR zgr_sco : e3w or e3t < 0 at point (i,j,k)= ', ji, jj, jk 2248 WRITE(numout,*) 'ERROR zgr_sco : e3w or e3t < 0 at point (i,j,k)= ', ji, jj, jk 2187 2249 WRITE(numout,*) 'e3w',fse3w(ji,jj,:) 2188 2250 WRITE(numout,*) 'e3t',fse3t(ji,jj,:) … … 2243 2305 INTEGER :: ji, jj, jk ! dummy loop argument 2244 2306 REAL(wp) :: zcoeft, zcoefw ! temporary scalars 2307 REAL(wp) :: ztmpu, ztmpv, ztmpf 2308 REAL(wp) :: ztmpu1, ztmpv1, ztmpf1 2245 2309 ! 2246 2310 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 … … 2296 2360 DO ji = 1, jpim1 2297 2361 DO jj = 1, jpjm1 2362 ! extended for Wetting/Drying case 2363 ztmpu = hbatt(ji,jj)+hbatt(ji+1,jj) 2364 ztmpv = hbatt(ji,jj)+hbatt(ji,jj+1) 2365 ztmpf = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) 2366 ztmpu1 = hbatt(ji,jj)*hbatt(ji+1,jj) 2367 ztmpv1 = hbatt(ji,jj)*hbatt(ji,jj+1) 2368 ztmpf1 = MIN(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) * & 2369 & MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) 2298 2370 DO jk = 1, jpk 2299 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) & 2300 & / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 2301 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) & 2302 & / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 2303 z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) & 2304 & + hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) & 2305 & / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 2306 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) & 2307 & / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 2308 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) & 2309 & / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 2371 IF((ztmpu1 < 0._wp.OR.ABS(ztmpu) < rn_wdmin1).AND.ln_wd) THEN 2372 z_esigtu3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) ) 2373 z_esigwu3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji+1,jj,jk) ) 2374 ELSE 2375 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) & 2376 & / ztmpu 2377 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) & 2378 & / ztmpu 2379 END IF 2380 2381 IF((ztmpv1 < 0._wp.OR.ABS(ztmpv) < rn_wdmin1).AND.ln_wd) THEN 2382 z_esigtv3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji,jj+1,jk) ) 2383 z_esigwv3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji,jj+1,jk) ) 2384 ELSE 2385 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) & 2386 & / ztmpv 2387 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) & 2388 & / ztmpv 2389 END IF 2390 2391 IF((ztmpf1 < 0._wp.OR.ABS(ztmpf) < rn_wdmin1).AND.ln_wd) THEN 2392 z_esigtf3(ji,jj,jk) = 0.25_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) & 2393 & + z_esigt3(ji,jj+1,jk) + z_esigt3(ji+1,jj+1,jk) ) 2394 ELSE 2395 z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) & 2396 & + hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) & 2397 & / ztmpf 2398 END IF 2399 2310 2400 ! 2311 2401 e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) … … 2347 2437 REAL(wp) :: zsmth ! smoothing around critical depth 2348 2438 REAL(wp) :: zzs, zzb ! Surface and bottom cell thickness in sigma space 2439 REAL(wp) :: ztmpu, ztmpv, ztmpf 2440 REAL(wp) :: ztmpu1, ztmpv1, ztmpf1 2349 2441 ! 2350 2442 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 … … 2425 2517 DO jj=1,jpj-1 2426 2518 2427 DO jk = 1, jpk 2428 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) / & 2429 ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 2430 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) / & 2431 ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 2432 z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) + & 2433 hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / & 2434 ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 2435 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) / & 2436 ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 2437 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) / & 2438 ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 2519 ! extend to suit for Wetting/Drying case 2520 ztmpu = hbatt(ji,jj)+hbatt(ji+1,jj) 2521 ztmpv = hbatt(ji,jj)+hbatt(ji,jj+1) 2522 ztmpf = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) 2523 ztmpu1 = hbatt(ji,jj)*hbatt(ji+1,jj) 2524 ztmpv1 = hbatt(ji,jj)*hbatt(ji,jj+1) 2525 ztmpf1 = MIN(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) * & 2526 & MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) 2527 DO jk = 1, jpk 2528 IF((ztmpu1 < 0._wp.OR.ABS(ztmpu) < rn_wdmin1).AND.ln_wd) THEN 2529 z_esigtu3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) ) 2530 z_esigwu3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji+1,jj,jk) ) 2531 ELSE 2532 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) & 2533 & / ztmpu 2534 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) & 2535 & / ztmpu 2536 END IF 2537 2538 IF((ztmpv1 < 0._wp.OR.ABS(ztmpv) < rn_wdmin1).AND.ln_wd) THEN 2539 z_esigtv3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji,jj+1,jk) ) 2540 z_esigwv3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji,jj+1,jk) ) 2541 ELSE 2542 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) & 2543 & / ztmpv 2544 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) & 2545 & / ztmpv 2546 END IF 2547 2548 IF((ztmpf1 < 0._wp.OR.ABS(ztmpf) < rn_wdmin1).AND.ln_wd) THEN 2549 z_esigtf3(ji,jj,jk) = 0.25_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) & 2550 & + z_esigt3(ji,jj+1,jk) + z_esigt3(ji+1,jj+1,jk) ) 2551 ELSE 2552 z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) & 2553 & + hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) & 2554 & / ztmpf 2555 END IF 2439 2556 2440 2557 e3t_0(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk) -
branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90
r5332 r5954 204 204 INTEGER :: ji, jj, jk 205 205 REAL(wp) :: zsal = 35.50 206 #if defined key_surge 207 REAL(wp) :: ztem = 10.0 208 #endif 206 209 !!---------------------------------------------------------------------- 207 210 ! … … 210 213 IF(lwp) WRITE(numout,*) '~~~~~~~~~~ and constant salinity (',zsal,' psu)' 211 214 ! 215 #if defined key_surge 216 tsn(:,:,:,jp_tem) = ztem * tmask(:,:,:) 217 tsb(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) 218 #else 212 219 DO jk = 1, jpk 213 220 tsn(:,:,jk,jp_tem) = ( ( ( 7.5 - 0. * ABS( gphit(:,:) )/30. ) * ( 1.-TANH((fsdept(:,:,jk)-80.)/30.) ) & … … 215 222 tsb(:,:,jk,jp_tem) = tsn(:,:,jk,jp_tem) 216 223 END DO 224 #endif 217 225 tsn(:,:,:,jp_sal) = zsal * tmask(:,:,:) 218 226 tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) -
branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
r5224 r5954 16 16 !! 3.4 ! 2011-11 (H. Liu) hpg_prj: Original code for s-coordinates 17 17 !! ! (A. Coward) suppression of hel, wdj and rot options 18 !! 3.6? ! 2014-09 (H. Liu) add Wetting/Drying pressure filter 18 19 !! 3.6 ! 2014-11 (P. Mathiot) hpg_isf: original code for ice shelf cavity 19 20 !!---------------------------------------------------------------------- … … 379 380 !! 380 381 INTEGER :: ji, jj, jk ! dummy loop indices 381 REAL(wp) :: zcoef0, zuap, zvap, znad ! temporary scalars 382 REAL(wp) :: zcoef0, zuap, zvap, znad, ztmp ! temporary scalars 383 LOGICAL :: ll_tmp1, ll_tmp2 ! 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 … … 398 402 ENDIF 399 403 404 IF(ln_wd) THEN 405 DO jj = 2, jpjm1 406 DO ji = 2, jpim1 407 IF ( tmask(ji+1,jj,1) == 0._wp) THEN 408 zcpx(ji,jj) = 1.0_wp 409 ELSE 410 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) 411 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) +& 412 & rn_wdmin1 + rn_wdmin2 413 414 IF(ll_tmp1) THEN 415 zcpx(ji,jj) = 1.0_wp 416 ELSE IF(ll_tmp2) THEN 417 ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 418 zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /& 419 & (sshn(ji+1,jj) - sshn(ji,jj))) 420 ELSE 421 zcpx(ji,jj) = 0._wp 422 END IF 423 ENDIF 424 425 IF ( tmask(ji,jj+1,1) == 0._wp) THEN 426 zcpy(ji,jj) = 1.0_wp 427 ELSE 428 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) 429 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) +& 430 & rn_wdmin1 + rn_wdmin2 431 432 IF(ll_tmp1) THEN 433 zcpy(ji,jj) = 1.0_wp 434 ELSE IF(ll_tmp2) THEN 435 ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 436 zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) /& 437 & (sshn(ji,jj+1) - sshn(ji,jj))) 438 ELSE 439 zcpy(ji,jj) = 0._wp 440 END IF 441 ENDIF 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 #if defined key_surge 448 ! Surface value 449 DO jj = 2, jpjm1 450 DO ji = fs_2, fs_jpim1 ! vector opt. 451 ! hydrostatic pressure gradient along s-surfaces 452 zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( fse3w(ji+1,jj ,1) * ( znad ) & 453 & - fse3w(ji ,jj ,1) * ( znad ) ) 454 zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( fse3w(ji ,jj+1,1) * ( znad ) & 455 & - fse3w(ji ,jj ,1) * ( znad ) ) 456 ! s-coordinate pressure gradient correction 457 zuap = -zcoef0 * ( 2._wp * znad ) & 458 & * ( fsde3w(ji+1,jj,1) - fsde3w(ji,jj,1) ) / e1u(ji,jj) 459 zvap = -zcoef0 * ( 2._wp * znad ) & 460 & * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj) 461 462 IF(ln_wd) THEN 463 zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 464 zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 465 zuap = zuap * zcpx(ji,jj) 466 zvap = zvap * zcpy(ji,jj) 467 ENDIF 468 469 ! add to the general momentum trend 470 ua(ji,jj,1) = ua(ji,jj,1) + ( zhpi(ji,jj,1) + zuap ) * umask(ji,jj,1) 471 va(ji,jj,1) = va(ji,jj,1) + ( zhpj(ji,jj,1) + zvap ) * vmask(ji,jj,1) 472 END DO 473 END DO 474 475 ! interior value (2=<jk=<jpkm1) 476 DO jk = 2, jpkm1 477 DO jj = 2, jpjm1 478 DO ji = fs_2, fs_jpim1 ! vector opt. 479 ! hydrostatic pressure gradient along s-surfaces 480 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj) & 481 & * ( fse3w(ji+1,jj,jk) * ( 2*znad ) & 482 & - fse3w(ji ,jj,jk) * ( 2*znad ) ) 483 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj) & 484 & * ( fse3w(ji,jj+1,jk) * ( 2*znad ) & 485 & - fse3w(ji,jj ,jk) * ( 2*znad ) ) 486 ! s-coordinate pressure gradient correction 487 zuap = -zcoef0 * ( 2._wp * znad ) & 488 & * ( fsde3w(ji+1,jj ,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj) 489 zvap = -zcoef0 * ( 2._wp * znad ) & 490 & * ( fsde3w(ji ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) 491 492 IF(ln_wd) THEN 493 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 494 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) 495 zuap = zuap * zcpx(ji,jj) 496 zvap = zvap * zcpy(ji,jj) 497 ENDIF 498 499 ! add to the general momentum trend 500 ua(ji,jj,jk) = ua(ji,jj,jk) + ( zhpi(ji,jj,jk) + zuap ) * umask(ji,jj,jk) 501 va(ji,jj,jk) = va(ji,jj,jk) + ( zhpj(ji,jj,jk) + zvap ) * vmask(ji,jj,jk) 502 END DO 503 END DO 504 END DO 505 ! 506 #else 400 507 ! Surface value 401 508 DO jj = 2, jpjm1 … … 411 518 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) + 2._wp * znad ) & 412 519 & * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj) 520 521 IF(ln_wd) THEN 522 zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 523 zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 524 zuap = zuap * zcpx(ji,jj) 525 zvap = zvap * zcpy(ji,jj) 526 ENDIF 527 413 528 ! add to the general momentum trend 414 529 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap … … 433 548 zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 434 549 & * ( fsde3w(ji ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) 550 551 IF(ln_wd) THEN 552 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 553 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) 554 zuap = zuap * zcpx(ji,jj) 555 zvap = zvap * zcpy(ji,jj) 556 ENDIF 557 435 558 ! add to the general momentum trend 436 559 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap … … 439 562 END DO 440 563 END DO 441 ! 564 #endif 565 ! 566 442 567 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj ) 568 IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 443 569 ! 444 570 END SUBROUTINE hpg_sco … … 719 845 REAL(wp) :: z1_10, cffu, cffx ! " " 720 846 REAL(wp) :: z1_12, cffv, cffy ! " " 847 LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables 721 848 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj 722 849 REAL(wp), POINTER, DIMENSION(:,:,:) :: dzx, dzy, dzz, dzu, dzv, dzw 723 850 REAL(wp), POINTER, DIMENSION(:,:,:) :: drhox, drhoy, drhoz, drhou, drhov, drhow 724 851 REAL(wp), POINTER, DIMENSION(:,:,:) :: rho_i, rho_j, rho_k 852 REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy !W/D pressure filter 725 853 !!---------------------------------------------------------------------- 726 854 ! … … 728 856 CALL wrk_alloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow ) 729 857 CALL wrk_alloc( jpi, jpj, jpk, rho_i, rho_j, rho_k, zhpi, zhpj ) 858 IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 730 859 ! 731 860 … … 734 863 IF(lwp) WRITE(numout,*) 'dyn:hpg_djc : hydrostatic pressure gradient trend' 735 864 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, density Jacobian with cubic polynomial scheme' 865 ENDIF 866 867 IF(ln_wd) THEN 868 DO jj = 2, jpjm1 869 DO ji = 2, jpim1 870 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) 871 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) +& 872 & rn_wdmin1 + rn_wdmin2 873 874 IF(ll_tmp1) THEN 875 zcpx(ji,jj) = 1.0_wp 876 ELSE IF(ll_tmp2) THEN 877 ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 878 zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /& 879 & (sshn(ji+1,jj) - sshn(ji,jj))) 880 ELSE 881 zcpx(ji,jj) = 0._wp 882 END IF 883 884 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) 885 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) +& 886 & rn_wdmin1 + rn_wdmin2 887 888 IF(ll_tmp1) THEN 889 zcpy(ji,jj) = 1.0_wp 890 ELSE IF(ll_tmp2) THEN 891 ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 892 zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) /& 893 & (sshn(ji,jj+1) - sshn(ji,jj))) 894 ELSE 895 zcpy(ji,jj) = 0._wp 896 END IF 897 END DO 898 END DO 899 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 736 900 ENDIF 737 901 … … 899 1063 zhpi(ji,jj,1) = ( rho_k(ji+1,jj ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) / e1u(ji,jj) 900 1064 zhpj(ji,jj,1) = ( rho_k(ji ,jj+1,1) - rho_k(ji,jj,1) - rho_j(ji,jj,1) ) / e2v(ji,jj) 1065 IF(ln_wd) THEN 1066 zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 1067 zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 1068 ENDIF 901 1069 ! add to the general momentum trend 902 1070 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) … … 918 1086 & + ( ( rho_k(ji,jj+1,jk) - rho_k(ji,jj,jk ) ) & 919 1087 & -( rho_j(ji,jj ,jk) - rho_j(ji,jj,jk-1) ) ) / e2v(ji,jj) 1088 IF(ln_wd) THEN 1089 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 1090 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) 1091 ENDIF 920 1092 ! add to the general momentum trend 921 1093 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) … … 928 1100 CALL wrk_dealloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow ) 929 1101 CALL wrk_dealloc( jpi, jpj, jpk, rho_i, rho_j, rho_k, zhpi, zhpj ) 1102 IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 930 1103 ! 931 1104 END SUBROUTINE hpg_djc … … 951 1124 !! The local variables for the correction term 952 1125 INTEGER :: jk1, jis, jid, jjs, jjd 1126 LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables 953 1127 REAL(wp) :: zuijk, zvijk, zpwes, zpwed, zpnss, zpnsd, zdeps 954 1128 REAL(wp) :: zrhdt1 … … 957 1131 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp 958 1132 REAL(wp), POINTER, DIMENSION(:,:) :: zsshu_n, zsshv_n 1133 REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy !W/D pressure filter 959 1134 !!---------------------------------------------------------------------- 960 1135 ! … … 962 1137 CALL wrk_alloc( jpi,jpj,jpk, zdept, zrhh ) 963 1138 CALL wrk_alloc( jpi,jpj, zsshu_n, zsshv_n ) 1139 IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 964 1140 ! 965 1141 IF( kt == nit000 ) THEN … … 974 1150 znad = 0.0_wp 975 1151 IF( lk_vvl ) znad = 1._wp 1152 1153 IF(ln_wd) THEN 1154 DO jj = 2, jpjm1 1155 DO ji = 2, jpim1 1156 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) 1157 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) +& 1158 & rn_wdmin1 + rn_wdmin2 1159 1160 IF(ll_tmp1) THEN 1161 zcpx(ji,jj) = 1.0_wp 1162 ELSE IF(ll_tmp2) THEN 1163 ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 1164 zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /& 1165 & (sshn(ji+1,jj) - sshn(ji,jj))) 1166 ELSE 1167 zcpx(ji,jj) = 0._wp 1168 END IF 1169 1170 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) 1171 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) +& 1172 & rn_wdmin1 + rn_wdmin2 1173 1174 IF(ll_tmp1.OR.ll_tmp2) THEN 1175 zcpy(ji,jj) = 1.0_wp 1176 ELSE IF(ll_tmp2) THEN 1177 ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 1178 zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) /& 1179 & (sshn(ji,jj+1) - sshn(ji,jj))) 1180 ELSE 1181 zcpy(ji,jj) = 0._wp 1182 END IF 1183 END DO 1184 END DO 1185 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 1186 ENDIF 976 1187 977 1188 ! Clean 3-D work arrays … … 1052 1263 END DO 1053 1264 END DO 1265 1266 CALL lbc_lnk (zsshu_n, 'U', 1.) 1267 CALL lbc_lnk (zsshv_n, 'V', 1.) 1054 1268 1055 1269 DO jj = 2, jpjm1 … … 1150 1364 ENDIF 1151 1365 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) 1366 IF(ln_wd) THEN 1367 zdpdx1 = zdpdx1 * zcpx(ji,jj) 1368 zdpdx2 = zdpdx2 * zcpx(ji,jj) 1369 ENDIF 1370 ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk) 1154 1371 ENDIF 1155 1372 … … 1207 1424 ENDIF 1208 1425 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) 1426 IF(ln_wd) THEN 1427 zdpdy1 = zdpdy1 * zcpy(ji,jj) 1428 zdpdy2 = zdpdy2 * zcpy(ji,jj) 1429 ENDIF 1430 1431 va(ji,jj,jk) = va(ji,jj,jk) + (zdpdy1 + zdpdy2) * vmask(ji,jj,jk) 1211 1432 ENDIF 1212 1433 … … 1219 1440 CALL wrk_dealloc( jpi,jpj,jpk, zdept, zrhh ) 1220 1441 CALL wrk_dealloc( jpi,jpj, zsshu_n, zsshv_n ) 1442 IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 1221 1443 ! 1222 1444 END SUBROUTINE hpg_prj -
branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90
r5328 r5954 137 137 DO jj = 2, jpjm1 138 138 DO ji = fs_2, fs_jpim1 ! vector opt. 139 ua(ji,jj,jk) = ua(ji,jj,jk) - ( zhke(ji+1,jj ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj) 140 va(ji,jj,jk) = va(ji,jj,jk) - ( zhke(ji ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj) 139 ua(ji,jj,jk) = ua(ji,jj,jk) - ( zhke(ji+1,jj ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj) * umask(ji,jj,jk) 140 va(ji,jj,jk) = va(ji,jj,jk) - ( zhke(ji ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj) * vmask(ji,jj,jk) 141 141 END DO 142 142 END DO -
branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilap.F90
r4990 r5954 140 140 141 141 ! Multiply by the eddy viscosity coef. (at u- and v-points) 142 zlu(:,:,jk) = zlu(:,:,jk) * ( fsahmu(:,:,jk) * (1-nkahm_smag) + nkahm_smag) 143 144 zlv(:,:,jk) = zlv(:,:,jk) * ( fsahmv(:,:,jk) * (1-nkahm_smag) + nkahm_smag) 142 zlu(:,:,jk) = zlu(:,:,jk) * ( fsahmu(:,:,jk) * (1-nkahm_smag) + nkahm_smag) * umask(:,:,jk) 143 144 zlv(:,:,jk) = zlv(:,:,jk) * ( fsahmv(:,:,jk) * (1-nkahm_smag) + nkahm_smag) * vmask(:,:,jk) 145 145 146 146 ! Contravariant "laplacian" … … 170 170 zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) 171 171 zut(ji,jj,jk) = ( zlu(ji,jj,jk) - zlu(ji-1,jj ,jk) & 172 & + zlv(ji,jj,jk) - zlv(ji ,jj-1,jk) ) / zbt 172 & + zlv(ji,jj,jk) - zlv(ji ,jj-1,jk) ) / zbt * tmask(ji,jj,jk) 173 173 END DO 174 174 END DO … … 197 197 & + ( zut(ji,jj+1,jk) - zut(ji ,jj,jk) ) / e2v(ji,jj) 198 198 ! add it to the general momentum trends 199 ua(ji,jj,jk) = ua(ji,jj,jk) + zua * ( fsahmu(ji,jj,jk)*nkahm_smag +(1 -nkahm_smag )) 200 va(ji,jj,jk) = va(ji,jj,jk) + zva * ( fsahmv(ji,jj,jk)*nkahm_smag +(1 -nkahm_smag )) 199 ua(ji,jj,jk) = ua(ji,jj,jk) + zua * ( fsahmu(ji,jj,jk)*nkahm_smag +(1 -nkahm_smag )) * umask(ji,jj,jk) 200 va(ji,jj,jk) = va(ji,jj,jk) + zva * ( fsahmv(ji,jj,jk)*nkahm_smag +(1 -nkahm_smag )) * vmask(ji,jj,jk) 201 201 END DO 202 202 END DO -
branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90
r5120 r5954 160 160 DO jj = 2, jpjm1 161 161 DO ji = fs_2, fs_jpim1 ! vector opt. 162 ua(ji,jj,jk) = ua(ji,jj,jk) + spgu(ji,jj) 163 va(ji,jj,jk) = va(ji,jj,jk) + spgv(ji,jj) 162 ua(ji,jj,jk) = ua(ji,jj,jk) + spgu(ji,jj) * umask(ji,jj,jk) 163 va(ji,jj,jk) = va(ji,jj,jk) + spgv(ji,jj) * vmask(ji,jj,jk) 164 164 END DO 165 165 END DO -
branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r5942 r5954 11 11 !! 3.5 ! 2013-07 (J. Chanut) Switch to Forward-backward time stepping 12 12 !! 3.6 ! 2013-11 (A. Coward) Update for z-tilde compatibility 13 !! 3.? ! 2014-09 (H. Liu) Add Wetting/Drying components 13 14 !!--------------------------------------------------------------------- 14 15 #if defined key_dynspg_ts || defined key_esopa … … 41 42 USE timing ! Timing 42 43 USE sbcapr ! surface boundary condition: atmospheric pressure 44 USE wadlmt ! wetting/drying flux limter 43 45 USE dynadv, ONLY: ln_dynadv_vec 44 46 #if defined key_agrif … … 140 142 INTEGER, INTENT(in) :: kt ! ocean time-step index 141 143 ! 142 LOGICAL :: ll_fw_start ! if true, forward integration 143 LOGICAL :: ll_init ! if true, special startup of 2d equations 144 INTEGER :: ji, jj, jk, jn ! dummy loop indices 145 INTEGER :: ikbu, ikbv, noffset ! local integers 146 REAL(wp) :: zraur, z1_2dt_b, z2dt_bf ! local scalars 147 REAL(wp) :: zx1, zy1, zx2, zy2 ! - - 148 REAL(wp) :: z1_12, z1_8, z1_4, z1_2 ! - - 149 REAL(wp) :: zu_spg, zv_spg ! - - 150 REAL(wp) :: zhura, zhvra ! - - 151 REAL(wp) :: za0, za1, za2, za3 ! - - 144 LOGICAL :: ll_fw_start ! if true, forward integration 145 LOGICAL :: ll_init ! if true, special startup of 2d equations 146 LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables used in W/D 147 INTEGER :: ji, jj, jk, jn ! dummy loop indices 148 INTEGER :: ikbu, ikbv, noffset ! local integers 149 REAL(wp) :: zraur, z1_2dt_b, z2dt_bf ! local scalars 150 REAL(wp) :: zx1, zy1, zx2, zy2 ! - - 151 REAL(wp) :: z1_12, z1_8, z1_4, z1_2 ! - - 152 REAL(wp) :: zu_spg, zv_spg ! - - 153 REAL(wp) :: zhura, zhvra ! - - 154 REAL(wp) :: za0, za1, za2, za3 ! - - 155 156 REAL(wp) :: ztmp ! temporary vaialbe 152 157 ! 153 158 REAL(wp), POINTER, DIMENSION(:,:) :: zun_e, zvn_e, zsshp2_e … … 157 162 REAL(wp), POINTER, DIMENSION(:,:) :: zsshu_a, zsshv_a 158 163 REAL(wp), POINTER, DIMENSION(:,:) :: zhf 164 REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy ! Wetting/Dying gravity filter coef. 159 165 !!---------------------------------------------------------------------- 160 166 ! … … 168 174 CALL wrk_alloc( jpi, jpj, zsshu_a, zsshv_a ) 169 175 CALL wrk_alloc( jpi, jpj, zhf ) 176 IF(ln_wd) CALL wrk_alloc( jpi, jpj, zcpx, zcpy ) 170 177 ! 171 178 ! !* Local constant initialization … … 176 183 zraur = 1._wp / rau0 177 184 ! 178 IF( kt == nit000 .AND. neuler == 0 ) THEN 185 IF( kt == nit000 .AND. neuler == 0 ) THEN ! reciprocal of baroclinic time step 179 186 z2dt_bf = rdt 180 187 ELSE … … 183 190 z1_2dt_b = 1.0_wp / z2dt_bf 184 191 ! 185 ll_init = ln_bt_av 192 ll_init = ln_bt_av ! if no time averaging, then no specific restart 186 193 ll_fw_start = .FALSE. 187 194 ! 188 189 IF (.NOT.ln_bt_fw) THEN ; noffset=-nn_baro ; ELSE ; noffset = 0 ; ENDIF 190 ! 191 IF( kt == nit000 ) THEN 195 ! time offset in steps for bdy data update 196 IF (.NOT.ln_bt_fw) THEN ; noffset=-nn_baro ; ELSE ; noffset = 0 ; ENDIF 197 ! 198 IF( kt == nit000 ) THEN !* initialisation 192 199 ! 193 200 IF(lwp) WRITE(numout,*) … … 378 385 ! !* Right-Hand-Side of the barotropic momentum equation 379 386 ! ! ---------------------------------------------------- 387 380 388 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 389 IF(ln_wd) THEN ! Calculating and applying W/D gravity filters 390 DO jj = 2, jpjm1 391 DO ji = 2, jpim1 392 IF ( tmask(ji+1,jj,1) == 0._wp ) THEN 393 zcpx = 1._wp 394 ELSE 395 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj)*tmask(ji,jj,1), -bathy(ji+1,jj)*tmask(ji+1,jj,1)) & 396 & .and. MAX(sshn(ji,jj) + bathy(ji,jj)*tmask(ji,jj,1), sshn(ji+1,jj) + bathy(ji+1,jj)*tmask(ji+1,jj,1)) & 397 & > rn_wdmin1 + rn_wdmin2 398 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj)*tmask(ji,jj,1), -bathy(ji+1,jj)*tmask(ji+1,jj,1)) + & 399 & rn_wdmin1 + rn_wdmin2 400 IF(ll_tmp1) THEN 401 zcpx(ji,jj) = 1.0_wp 402 ELSE IF(ll_tmp2) THEN 403 ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 404 zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj)*tmask(ji+1,jj,1) - sshn(ji,jj) - bathy(ji,jj)*tmask(ji,jj,1)) /& 405 & (sshn(ji+1,jj) - sshn(ji,jj))) 406 ELSE 407 zcpx(ji,jj) = 0._wp 408 END IF 409 ENDIF 410 411 IF ( tmask(ji,jj+1,1) == 0._wp ) THEN 412 zcpy = 1._wp 413 ELSE 414 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj)*tmask(ji,jj,1), -bathy(ji,jj+1)*tmask(ji,jj+1,1)) & 415 & .and. MAX(sshn(ji,jj) + bathy(ji,jj)*tmask(ji,jj,1), sshn(ji,jj+1) + bathy(ji,jj+1)*tmask(ji,jj+1,1)) & 416 & > rn_wdmin1 + rn_wdmin2 417 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj)*tmask(ji,jj,1), -bathy(ji,jj+1)*tmask(ji,jj+1,1)) + & 418 & rn_wdmin1 + rn_wdmin2 419 IF(ll_tmp1) THEN 420 zcpy(ji,jj) = 1.0_wp 421 ELSE IF(ll_tmp2) THEN 422 ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 423 zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1)*tmask(ji,jj+1,1) - sshn(ji,jj) - bathy(ji,jj)*tmask(ji,jj,1)) /& 424 & (sshn(ji,jj+1) - sshn(ji,jj))) 425 ELSE 426 zcpy(ji,jj) = 0._wp 427 END IF 428 ENDIF 429 END DO 430 END DO 431 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 432 ENDIF 433 434 IF(ln_wd) THEN ! Calculating and applying W/D gravity filters 435 DO jj = 2, jpjm1 436 DO ji = 2, jpim1 437 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) / & 438 & e1u(ji,jj) * zcpx(ji,jj) 439 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) / & 440 & e2v(ji,jj) * zcpy(ji,jj) 441 END DO 442 END DO 443 ELSE 444 DO jj = 2, jpjm1 445 DO ji = fs_2, fs_jpim1 ! vector opt. 446 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) / e1u(ji,jj) 447 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) / e2v(ji,jj) 448 END DO 449 END DO 450 END IF 387 451 ENDIF 388 452 … … 416 480 ! 417 481 ! 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(:,:) 482 IF(ln_wd) THEN 483 zu_frc(:,:) = zu_frc(:,:) + MAX(hur(:,:) * bfrua(:,:),-1._wp / rdtbt) * zwx(:,:) 484 zv_frc(:,:) = zv_frc(:,:) + MAX(hvr(:,:) * bfrva(:,:),-1._wp / rdtbt) * zwy(:,:) 485 ELSE 486 zu_frc(:,:) = zu_frc(:,:) + hur(:,:) * bfrua(:,:) * zwx(:,:) 487 zv_frc(:,:) = zv_frc(:,:) + hvr(:,:) * bfrva(:,:) * zwy(:,:) 488 END IF 420 489 ! 421 490 IF (ln_bt_fw) THEN ! Add wind forcing … … 488 557 ENDIF 489 558 ! 559 IF(ln_wd) THEN !preserve the positivity of water depth 560 !ssh[b,n,a] should have already been processed for this 561 sshbb_e(:,:jj) = MAX( sshbb_e(:,:), (rn_wdmin1 - bathy(:,:)) ) *tmask(:,:,1) 562 sshb_e(:,:jj) = MAX( sshb_e(:,:) , (rn_wdmin1 - bathy(:,:)) ) *tmask(:,:,1) 563 ENDIF 564 490 565 IF (ln_bt_fw) THEN ! FORWARD integration: start from NOW fields 491 566 sshn_e(:,:) = sshn (:,:) … … 523 598 ! Update only tidal forcing at open boundaries 524 599 #if defined key_tide 525 IF ( lk_bdy .AND. lk_tide ) CALL bdy_dta_tides( kt, kit=jn, time_offset=(noffset+1) ) 526 IF ( ln_tide_pot .AND. lk_tide ) CALL upd_tide( kt, kit=jn, time_offset=noffset ) 600 IF ( lk_bdy .AND. lk_tide ) CALL bdy_dta_tides( kt, kit=jn, time_offset=(noffset+1) ) 601 IF ( ln_tide_pot .AND. lk_tide ) CALL upd_tide( kt, kit=jn, time_offset=noffset ) 527 602 #endif 528 603 ! … … 561 636 zhup2_e (:,:) = hu_0(:,:) + zwx(:,:) ! Ocean depth at U- and V-points 562 637 zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:) 638 IF(ln_wd) THEN 639 zhup2_e(:,:) = MAX(zhup2_e (:,:), rn_wdmin1) 640 zhvp2_e(:,:) = MAX(zhvp2_e (:,:), rn_wdmin1) 641 END IF 563 642 ELSE 564 643 zhup2_e (:,:) = hu(:,:) … … 599 678 ENDIF 600 679 #endif 680 681 IF(ln_wd) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt) 601 682 ! 602 683 ! Sum over sub-time-steps to compute advective velocities … … 613 694 END DO 614 695 ssha_e(:,:) = ( sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) ) ) * tmask(:,:,1) 696 IF(ln_wd) ssha_e(:,:) = MAX(ssha_e(:,:), (rn_wdmin1 - bathy(:,:)) ) * tmask(:,:,1) 615 697 CALL lbc_lnk( ssha_e, 'T', 1._wp ) 616 698 … … 676 758 END DO 677 759 END DO 760 761 IF(ln_wd) THEN 762 zhust_e(:,:) = MAX(zhust_e (:,:), rn_wdmin1 ) 763 zhvst_e(:,:) = MAX(zhvst_e (:,:), rn_wdmin1 ) 764 END IF 765 !shall we call lbc_lnk for zhu(v)st_e() here? 766 678 767 ENDIF 679 768 ! … … 738 827 ! 739 828 ! Add bottom stresses: 740 zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * zun_e(:,:) * hur_e(:,:) 741 zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * zvn_e(:,:) * hvr_e(:,:) 742 ! 743 ! 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 829 IF(ln_wd) THEN 830 zu_trd(:,:) = zu_trd(:,:) + MAX(bfrua(:,:) * hur_e(:,:), -1._wp / rdtbt) * zun_e(:,:) 831 zv_trd(:,:) = zv_trd(:,:) + MAX(bfrva(:,:) * hvr_e(:,:), -1._wp / rdtbt) * zvn_e(:,:) 832 ELSE 833 zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * zun_e(:,:) * hur_e(:,:) 834 zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * zvn_e(:,:) * hvr_e(:,:) 835 END IF 836 837 IF(ln_wd) THEN ! Calculating and applying W/D gravity filters 838 DO jj = 2, jpjm1 839 DO ji = 2, jpim1 840 IF ( tmask(ji+1,jj,1) == 0._wp ) THEN 841 zcpx = 1._wp 842 ELSE 843 ll_tmp1 = MIN(zsshp2_e(ji,jj), zsshp2_e(ji+1,jj)) > MAX(-bathy(ji,jj)*tmask(ji,jj,1), -bathy(ji+1,jj)*tmask(ji+1,jj,1)) & 844 & .and. MAX(zsshp2_e(ji,jj) + bathy(ji,jj)*tmask(ji,jj,1), zsshp2_e(ji+1,jj) + bathy(ji+1,jj)*tmask(ji+1,jj,1)) & 845 & > rn_wdmin1 + rn_wdmin2 846 ll_tmp2 = MAX(zsshp2_e(ji,jj), zsshp2_e(ji+1,jj)) > MAX(-bathy(ji,jj)*tmask(ji,jj,1), -bathy(ji+1,jj)*tmask(ji+1,jj,1)) + & 847 & rn_wdmin1 + rn_wdmin2 848 IF(ll_tmp1) THEN 849 zcpx(ji,jj) = 1.0_wp 850 ELSE IF(ll_tmp2) THEN 851 ! no worries about zsshp2_e(ji+1,jj)-zsshp2_e(ji,jj) = 0, it won't happen ! here 852 zcpx(ji,jj) = ABS((zsshp2_e(ji+1,jj) + bathy(ji+1,jj)*tmask(ji+1,jj,1) - zsshp2_e(ji,jj) - bathy(ji,jj)*tmask(ji,jj,1)) /& 853 & (zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj))) 854 ELSE 855 zcpx(ji,jj) = 0._wp 856 END IF 857 ENDIF 858 859 IF ( tmask(ji,jj+1,1) == 0._wp ) THEN 860 zcpy = 1._wp 861 ELSE 862 ll_tmp1 = MIN(zsshp2_e(ji,jj), zsshp2_e(ji,jj+1)) > MAX(-bathy(ji,jj)*tmask(ji,jj,1), -bathy(ji,jj+1)*tmask(ji,jj+1,1))& 863 & .and. MAX(zsshp2_e(ji,jj) + bathy(ji,jj)*tmask(ji,jj,1), zsshp2_e(ji,jj+1) + bathy(ji,jj+1)*tmask(ji,jj+1,1)) & 864 & > rn_wdmin1 + rn_wdmin2 865 ll_tmp2 = MAX(zsshp2_e(ji,jj), zsshp2_e(ji,jj+1)) > MAX(-bathy(ji,jj)*tmask(ji,jj,1), -bathy(ji,jj+1)*tmask(ji,jj+1,1)) + & 866 & rn_wdmin1 + rn_wdmin2 867 IF(ll_tmp1) THEN 868 zcpy(ji,jj) = 1.0_wp 869 ELSE IF(ll_tmp2) THEN 870 ! no worries about zsshp2_e(ji,jj+1)-zsshp2_e(ji,jj) = 0, it won't happen ! here 871 zcpy(ji,jj) = ABS((zsshp2_e(ji,jj+1) + bathy(ji,jj+1)*tmask(ji,jj+1,1) - zsshp2_e(ji,jj) - bathy(ji,jj)*tmask(ji,jj,1)) /& 872 & (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj))) 873 ELSE 874 zcpy(ji,jj) = 0._wp 875 END IF 876 ENDIF 877 END DO 878 END DO 879 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 880 ENDIF 881 882 IF(ln_wd) THEN 883 DO jj = 2, jpjm1 884 DO ji = 2, jpim1 885 ! Add surface pressure gradient 886 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj) 887 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj) 888 zwx(ji,jj) = zu_spg * zcpx(ji,jj) 889 zwy(ji,jj) = zv_spg * zcpy(ji,jj) 890 END DO 891 END DO 892 ELSE 893 DO jj = 2, jpjm1 894 DO ji = fs_2, fs_jpim1 ! vector opt. 895 ! Add surface pressure gradient 896 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj) 897 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj) 898 zwx(ji,jj) = zu_spg 899 zwy(ji,jj) = zv_spg 900 END DO 901 END DO 902 END IF 753 903 ! 754 904 ! Set next velocities: … … 774 924 DO ji = fs_2, fs_jpim1 ! vector opt. 775 925 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) & 926 IF(ln_wd) THEN 927 zhura = MAX(hu_0(ji,jj) + zsshu_a(ji,jj), rn_wdmin1) 928 zhvra = MAX(hv_0(ji,jj) + zsshv_a(ji,jj), rn_wdmin1) 929 ELSE 930 zhura = hu_0(ji,jj) + zsshu_a(ji,jj) 931 zhvra = hv_0(ji,jj) + zsshv_a(ji,jj) 932 END IF 933 934 zhura = umask(ji,jj,1)/(zhura + 1._wp - umask(ji,jj,1)) 935 zhvra = vmask(ji,jj,1)/(zhvra + 1._wp - vmask(ji,jj,1)) 936 937 ua_e(ji,jj) = ( hu_e(ji,jj) * zun_e(ji,jj) & 938 & + rdtbt * ( zhust_e(ji,jj) * zwx(ji,jj) & 781 939 & + zhup2_e(ji,jj) * zu_trd(ji,jj) & 782 940 & + hu(ji,jj) * zu_frc(ji,jj) ) & … … 794 952 IF( lk_vvl ) THEN !* Update ocean depth (variable volume case only) 795 953 ! ! ---------------------------------------------- 796 hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 797 hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 954 IF(ln_wd) THEN 955 hu_e (:,:) = MAX(hu_0(:,:) + zsshu_a(:,:), rn_wdmin1 * umask(:,:,1) ) 956 hv_e (:,:) = MAX(hv_0(:,:) + zsshv_a(:,:), rn_wdmin1 * vmask(:,:,1) ) 957 ELSE 958 hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 959 hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 960 END IF 961 798 962 hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1._wp - umask(:,:,1) ) 799 963 hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1._wp - vmask(:,:,1) ) … … 926 1090 CALL wrk_dealloc( jpi, jpj, zsshu_a, zsshv_a ) 927 1091 CALL wrk_dealloc( jpi, jpj, zhf ) 1092 IF(ln_wd) CALL wrk_dealloc( jpi, jpj, zcpx, zcpy ) 928 1093 ! 929 1094 IF( nn_timing == 1 ) CALL timing_stop('dyn_spg_ts') -
branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90
r5029 r5954 285 285 zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1) 286 286 zx2 = zwx(ji ,jj) + zwx(ji ,jj+1) 287 pua(ji,jj,jk) = pua(ji,jj,jk) + zfact2 / e1u(ji,jj) * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 288 pva(ji,jj,jk) = pva(ji,jj,jk) - zfact2 / e2v(ji,jj) * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) 287 pua(ji,jj,jk) = pua(ji,jj,jk) + zfact2 / e1u(ji,jj) * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) * umask(ji,jj,jk) 288 pva(ji,jj,jk) = pva(ji,jj,jk) - zfact2 / e2v(ji,jj) * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) * vmask(ji,jj,jk) 289 289 END DO 290 290 END DO … … 404 404 zcva =-zfact2 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) 405 405 ! mixed vorticity trend added to the momentum trends 406 ua(ji,jj,jk) = ua(ji,jj,jk) + zcua + zua407 va(ji,jj,jk) = va(ji,jj,jk) + zcva + zva406 ua(ji,jj,jk) = ua(ji,jj,jk) + ( zcua + zua ) * umask(ji,jj,jk) 407 va(ji,jj,jk) = va(ji,jj,jk) + ( zcva + zva ) * vmask(ji,jj,jk) 408 408 END DO 409 409 END DO … … 522 522 zvau =-zfact1 / e2v(ji,jj) * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 523 523 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) 524 pua(ji,jj,jk) = pua(ji,jj,jk) + zuav * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 525 pva(ji,jj,jk) = pva(ji,jj,jk) + zvau * ( zwz(ji-1,jj ) + zwz(ji,jj) ) 524 pua(ji,jj,jk) = pua(ji,jj,jk) + zuav * ( zwz(ji ,jj-1) + zwz(ji,jj) ) * umask(ji,jj,jk) 525 pva(ji,jj,jk) = pva(ji,jj,jk) + zvau * ( zwz(ji-1,jj ) + zwz(ji,jj) ) * vmask(ji,jj,jk) 526 526 END DO 527 527 END DO … … 691 691 zva = - zfac12 / e2v(ji,jj) * ( ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji ,jj+1) & 692 692 & + ztnw(ji,jj ) * zwx(ji-1,jj ) + ztne(ji,jj ) * zwx(ji ,jj ) ) 693 pua(ji,jj,jk) = pua(ji,jj,jk) + zua 694 pva(ji,jj,jk) = pva(ji,jj,jk) + zva 693 pua(ji,jj,jk) = pua(ji,jj,jk) + zua * umask(ji,jj,jk) 694 pva(ji,jj,jk) = pva(ji,jj,jk) + zva * vmask(ji,jj,jk) 695 695 END DO 696 696 END DO -
branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90
r5120 r5954 124 124 zva = - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) ) 125 125 ! ! add the trends to the general momentum trends 126 ua(ji,jj,jk) = ua(ji,jj,jk) + zua 127 va(ji,jj,jk) = va(ji,jj,jk) + zva 126 ua(ji,jj,jk) = ua(ji,jj,jk) + zua * umask(ji,jj,jk) 127 va(ji,jj,jk) = va(ji,jj,jk) + zva * vmask(ji,jj,jk) 128 128 END DO 129 129 END DO -
branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
r5942 r5954 9 9 !! - ! 2010-09 (D.Storkey and E.O'Dea) bug fixes for BDY module 10 10 !! 3.3 ! 2011-10 (M. Leclair) split former ssh_wzv routine and remove all vvl related work 11 !! 3.? ! 2014-09 (H. Liu) add wetting and drying 11 12 !!---------------------------------------------------------------------- 12 13 … … 39 40 USE wrk_nemo ! Memory Allocation 40 41 USE timing ! Timing 42 USE wadlmt ! Wetting/Drying flux limting 41 43 42 44 IMPLICIT NONE … … 91 93 ENDIF 92 94 ! 93 CALL div_cur( kt ) ! Horizontal divergence & Relative vorticity94 !95 95 z2dt = 2._wp * rdt ! set time step size (Euler/Leapfrog) 96 96 IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt 97 97 ! 98 z1_rau0 = 0.5_wp * r1_rau0 99 ! 100 IF(ln_wd) CALL wad_lmt(sshb, z1_rau0 * (emp_b(:,:) + emp(:,:)), z2dt) 101 102 CALL div_cur( kt ) ! Horizontal divergence & Relative vorticity 103 ! 98 104 ! !------------------------------! 99 105 ! ! After Sea Surface Height ! … … 107 113 ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations. 108 114 ! 109 z1_rau0 = 0.5_wp * r1_rau0110 115 ssha(:,:) = ( sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * ssmask(:,:) 111 116 -
branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/ICB/icbrst.F90
r5341 r5954 24 24 USE icb_oce ! define iceberg arrays 25 25 USE icbutl ! iceberg utility routines 26 ! added for ln_rstdate, should include separate branch at later date... 27 USE phycst ! for rday 28 USE ioipsl, ONLY : ju2ymds ! for calendar 26 29 27 30 IMPLICIT NONE … … 235 238 TYPE(iceberg), POINTER :: this 236 239 TYPE(point) , POINTER :: pt 237 !!---------------------------------------------------------------------- 240 ! included for ln_rstdate.... 241 INTEGER :: iyear, imonth, iday 242 REAL (wp) :: zsec 243 CHARACTER(LEN=20) :: clkt ! ocean time-step deine as a character 244 !!---------------------------------------------------------------------- 245 246 IF ( ln_rstdate ) THEN 247 CALL ju2ymds( fjulday + rdttra(1) / rday, iyear, imonth, iday, zsec ) 248 WRITE(clkt, '(i4.4,2i2.2)') iyear, imonth, iday 249 ELSE 250 IF( kt > 999999999 ) THEN ; WRITE(clkt, * ) kt 251 ELSE ; WRITE(clkt, '(i8.8)') kt 252 ENDIF 253 ENDIF 238 254 239 255 ! Assume we write iceberg restarts to same directory as ocean restarts. 240 256 cl_path = TRIM(cn_ocerst_outdir) 241 257 IF( cl_path(LEN_TRIM(cl_path):) /= '/' ) cl_path = TRIM(cl_path) // '/' 258 259 ! changed for ln_rstdate.... 242 260 IF( lk_mpp ) THEN 243 WRITE(cl_filename,'(A,"_icebergs_", I8.8,"_restart_",I4.4,".nc")') TRIM(cexper), kt, narea-1261 WRITE(cl_filename,'(A,"_icebergs_",A,"_restart_",I4.4,".nc")') TRIM(cexper), TRIM(ADJUSTL(clkt)), narea-1 244 262 ELSE 245 WRITE(cl_filename,'(A,"_icebergs_", I8.8,"_restart.nc")') TRIM(cexper), kt263 WRITE(cl_filename,'(A,"_icebergs_",A,"_restart.nc")') TRIM(cexper), TRIM(ADJUSTL(clkt)) 246 264 ENDIF 247 265 IF (nn_verbose_level >= 0) WRITE(numout,'(2a)') 'icebergs, write_restart: creating ',TRIM(cl_path)//TRIM(cl_filename) -
branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/IOM/in_out_manager.F90
r5518 r5954 48 48 LOGICAL :: ln_clobber !: clobber (overwrite) an existing file 49 49 INTEGER :: nn_chunksz !: chunksize (bytes) for NetCDF file (works only with iom_nf90 routines) 50 ! included for rstdate and rstdir, should be moved to seperate branch! 51 LOGICAL :: ln_rstdate !: datestamping of restarts 52 CHARACTER(lc) :: cn_ocerst_indir !: restart input directory 53 CHARACTER(lc) :: cn_ocerst_outdir !: restart output directory 50 54 #if defined key_netcdf4 51 55 !!---------------------------------------------------------------------- -
branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90
r5407 r5954 24 24 USE trdmxl_oce ! ocean active mixed layer tracers trends variables 25 25 USE divcur ! hor. divergence and curl (div & cur routines) 26 ! for ln_rstdate.... 27 USE ioipsl, ONLY : ju2ymds ! for calendar 26 28 27 29 IMPLICIT NONE … … 58 60 CHARACTER(LEN=50) :: clname ! ocean output restart file name 59 61 CHARACTER(lc) :: clpath ! full path to ocean output restart file 62 ! Included for ln_rstdate 63 INTEGER :: iyear, imonth, iday 64 REAL (wp) :: zsec 65 60 66 !!---------------------------------------------------------------------- 61 67 ! … … 81 87 IF( kt == nitrst - 1 .OR. nstock == 1 .OR. ( kt == nitend .AND. .NOT. lrst_oce ) ) THEN 82 88 IF( nitrst <= nitend .AND. nitrst > 0 ) THEN 83 ! beware of the format used to write kt (default is i8.8, that should be large enough...) 84 IF( nitrst > 999999999 ) THEN ; WRITE(clkt, * ) nitrst 85 ELSE ; WRITE(clkt, '(i8.8)') nitrst 89 IF ( ln_rstdate ) THEN 90 CALL ju2ymds( fjulday + rdttra(1) / rday, iyear, imonth, iday, zsec ) 91 WRITE(clkt, '(i4.4,2i2.2)') iyear, imonth, iday 92 ELSE 93 ! beware of the format used to write kt (default is i8.8, that should be large enough...) 94 IF( nitrst > 999999999 ) THEN ; WRITE(clkt, * ) nitrst 95 ELSE ; WRITE(clkt, '(i8.8)') nitrst 96 ENDIF 86 97 ENDIF 87 98 ! create the file … … 89 100 clpath = TRIM(cn_ocerst_outdir) 90 101 IF( clpath(LEN_TRIM(clpath):) /= '/' ) clpath = TRIM(clpath) // '/' 102 ! 91 103 IF(lwp) THEN 92 104 WRITE(numout,*) … … 102 114 ENDIF 103 115 ENDIF 104 !116 ! 105 117 CALL iom_open( TRIM(clpath)//TRIM(clname), numrow, ldwrt = .TRUE., kiolib = jprstlib ) 106 118 lrst_oce = .TRUE. … … 172 184 LOGICAL :: llok 173 185 CHARACTER(lc) :: clpath ! full path to ocean output restart file 186 CHARACTER(lc) :: clpath ! full path to ocean output restart file 174 187 !!---------------------------------------------------------------------- 175 188 ! -
branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90
r5942 r5954 91 91 REAL(wp) :: rn_zqt ! z(q,t) : height of humidity and temperature measurements 92 92 REAL(wp) :: rn_zu ! z(u) : height of wind measurements 93 LOGICAL :: ln_charnock ! logical flag for charnock wind stress in surge model(true) or not(false) 93 94 94 95 !! * Substitutions … … 151 152 & sn_wndi, sn_wndj, sn_humi , sn_qsr , & 152 153 & sn_qlw , sn_tair, sn_prec , sn_snow, & 153 & sn_tdif, rn_zqt , rn_zu154 & sn_tdif, rn_zqt , rn_zu, ln_charnock 154 155 !!--------------------------------------------------------------------- 155 156 ! … … 247 248 INTEGER :: ji, jj ! dummy loop indices 248 249 REAL(wp) :: zcoef_qsatw, zztmp ! local variable 250 REAL(wp) :: z_z0, z_Cd1 ! local variable 251 REAL(wp) :: i ! local variable 252 REAL(wp) :: charn_const=0.0275 ! local variable 249 253 REAL(wp), DIMENSION(:,:), POINTER :: zwnd_i, zwnd_j ! wind speed components at T-point 250 254 REAL(wp), DIMENSION(:,:), POINTER :: zqsatw ! specific humidity at pst … … 300 304 ! I Radiative FLUXES ! 301 305 ! ----------------------------------------------------------------------------- ! 302 306 303 307 ! ocean albedo assumed to be constant + modify now Qsr to include the diurnal cycle ! Short Wave 304 308 zztmp = 1. - albo 309 #if defined key_surge 310 qsr(:,:)=0._wp 311 zqlw(:,:) = 0._wp 312 #else 305 313 IF( ln_dm2dc ) THEN ; qsr(:,:) = zztmp * sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) * tmask(:,:,1) 306 314 ELSE ; qsr(:,:) = zztmp * sf(jp_qsr)%fnow(:,:,1) * tmask(:,:,1) … … 308 316 309 317 zqlw(:,:) = ( sf(jp_qlw)%fnow(:,:,1) - Stef * zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:) ) * tmask(:,:,1) ! Long Wave 318 #endif 310 319 ! ----------------------------------------------------------------------------- ! 311 320 ! II Turbulent FLUXES ! 312 321 ! ----------------------------------------------------------------------------- ! 313 314 ! ... specific humidity at SST and IST 315 zqsatw(:,:) = zcoef_qsatw * EXP( -5107.4 / zst(:,:) ) 316 317 ! ... NCAR Bulk formulae, computation of Cd, Ch, Ce at T-point : 318 CALL turb_core_2z( rn_zqt, rn_zu, zst, sf(jp_tair)%fnow, zqsatw, sf(jp_humi)%fnow, wndm, & 319 & Cd, Ch, Ce, zt_zu, zq_zu ) 320 322 IF (ln_charnock) THEN 323 Cd(:,:)=0.0001_wp 324 DO jj = 1,jpj 325 DO ji = 1,jpi 326 z_Cd1=0._wp 327 i=1 328 !Iterate 329 DO WHILE((abs(Cd(ji,jj)-z_Cd1))>1E-6) 330 z_Cd1=Cd(ji,jj) 331 z_z0=charn_const*z_Cd1*wndm(ji,jj)**2/grav 332 Cd(ji,jj)=(0.41_wp/log(10._wp/z_z0))**2 333 i=i+1 334 ENDDO 335 ENDDO 336 ENDDO 337 ELSE 338 339 ! ... specific humidity at SST and IST 340 zqsatw(:,:) = zcoef_qsatw * EXP( -5107.4 / zst(:,:) ) 341 342 ! ... NCAR Bulk formulae, computation of Cd, Ch, Ce at T-point : 343 CALL turb_core_2z( rn_zqt, rn_zu, zst, sf(jp_tair)%fnow, zqsatw, sf(jp_humi)%fnow, wndm, & 344 & Cd, Ch, Ce, zt_zu, zq_zu ) 345 346 ENDIF 347 321 348 ! ... tau module, i and j component 322 349 DO jj = 1, jpj … … 352 379 ! Turbulent fluxes over ocean 353 380 ! ----------------------------- 381 #if ! defined key_surge 354 382 IF( ABS( rn_zu - rn_zqt) < 0.01_wp ) THEN 355 383 !! q_air and t_air are (or "are almost") given at 10m (wind reference height) … … 363 391 ENDIF 364 392 zqla (:,:) = Lv * zevap(:,:) ! Latent Heat 393 #endif 365 394 366 395 IF(ln_ctl) THEN … … 379 408 ! ----------------------------------------------------------------------------- ! 380 409 ! 410 #if defined key_surge 411 emp (:,:) = 0._wp 412 qns(:,:) = 0._wp 413 #else 381 414 emp (:,:) = ( zevap(:,:) & ! mass flux (evap. - precip.) 382 415 & - sf(jp_prec)%fnow(:,:,1) * rn_pfac ) * tmask(:,:,1) … … 389 422 & + sf(jp_snow)%fnow(:,:,1) * rn_pfac & ! add solid precip heat content at min(Tair,Tsnow) 390 423 & * ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) 424 #endif 391 425 ! 392 426 #if defined key_lim3 -
branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
r5942 r5954 250 250 nstop = nstop + 1 251 251 ENDIF 252 IF ( lk_surge .and. .not. ( ln_blk_core .or. ln_ana ) ) & 253 & CALL ctl_stop( ' surge model only compatible with analytical fluxes or core formulae' ) 252 254 IF(lwp) THEN 253 255 WRITE(numout,*) -
branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/SBC/tide.h90
r4292 r5954 28 28 Wave(18) = tide( 'L2' , 0.006694 , 2 , 2 , -1 , 2 , -1 , 0 , +180 , 2 , -2 , 0 , 0 , 0 , 215 ) 29 29 Wave(19) = tide( 'T2' , 0.006614 , 2 , 2 , 0 , -1 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ) 30 ! 31 ! !! name_tide , equitide , nutide , nt , ns , nh , np , np1 , shift , nksi , nnu0 , nnu1 , nnu2 , R , formula !! 32 Wave(20) = tide( 'MNS2' , 0.000000 , 2 , 2 , -5 , 4 , 1 , 0 , 0 , 4 , -4 , 0 , 0 , 0 , 6 ) 33 Wave(21) = tide( 'Lam2' , 0.001760 , 2 , 2 , -1 , 0 , 1 , 0 , +180 , 2 , -2 , 0 , 0 , 0 , 78 ) 34 Wave(22) = tide( 'MSN2' , 0.000000 , 2 , 2 , 1 , 0 , 1 , 0 , 0 , 2 , -2 , 0 , 2 , 0 , 6 ) 35 Wave(23) = tide( '2SM2' , 0.000000 , 2 , 2 , 2 , -2 , 0 , 0 , 0 , -2 , 2 , 0 , 0 , 0 , 16 ) 36 Wave(24) = tide( 'MO3' , 0.000000 , 3 , 3 , -4 , 1 , 0 , 0 , +90 , 2 , -2 , 0 , 0 , 0 , 13 ) 37 Wave(25) = tide( 'MK3' , 0.000000 , 3 , 3 , -2 , 3 , 0 , 0 , -90 , 2 , -2 , -1 , 0 , 0 , 10 ) 38 Wave(26) = tide( 'MN4' , 0.000000 , 4 , 4 , -5 , 4 , 1 , 0 , 0 , 4 , -4 , 0 , 0 , 0 , 1 ) 39 Wave(27) = tide( 'MS4' , 0.000000 , 4 , 4 , -2 , 2 , 0 , 0 , 0 , 2 , -2 , 0 , 0 , 0 , 2 ) 40 Wave(28) = tide( 'M6' , 0.000000 , 6 , 6 , -6 , 6 , 0 , 0 , 0 , 6 , -6 , 0 , 0 , 0 , 4 ) 41 Wave(29) = tide( '2MS6' , 0.000000 , 6 , 6 , -4 , 4 , 0 , 0 , 0 , 4 , -4 , 0 , 0 , 0 , 6 ) 42 Wave(30) = tide( '2MK6' , 0.000000 , 6 , 6 , -4 , 6 , 0 , 0 , 0 , 4 , -4 , 0 , -2 , 0 , 5 ) 43 Wave(31) = tide( '3M2S2' , 0.000000 , 2 , 2 , -6 , 6 , 0 , 0 , 0 , 6 , -6 , 0 , 0 , 0 , 12 ) -
branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/SBC/tide_mod.F90
r5215 r5954 16 16 PUBLIC tide_init_Wave ! called by tideini and diaharm modules 17 17 18 INTEGER, PUBLIC, PARAMETER :: jpmax_harmo = 19!: maximum number of harmonic18 INTEGER, PUBLIC, PARAMETER :: jpmax_harmo = 31 !: maximum number of harmonic 19 19 20 20 TYPE, PUBLIC :: tide 21 CHARACTER(LEN= 4) :: cname_tide21 CHARACTER(LEN=5) :: cname_tide 22 22 REAL(wp) :: equitide 23 23 INTEGER :: nutide -
branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/SBC/tideini.F90
r5215 r5954 48 48 INTEGER :: ji, jk 49 49 INTEGER, INTENT( in ) :: kt ! ocean time-step 50 CHARACTER(LEN= 4), DIMENSION(jpmax_harmo) :: clname50 CHARACTER(LEN=5), DIMENSION(jpmax_harmo) :: clname 51 51 INTEGER :: ios ! Local integer output status for namelist read 52 52 ! -
branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfini.F90
r5386 r5954 115 115 ioptio = ioptio+1 116 116 ENDIF 117 IF( ioptio == 0 .OR. ioptio > 1 .AND. .NOT. lk_esopa ) & 118 & CALL ctl_stop( ' one and only one vertical diffusion option has to be defined ' ) 117 IF( ioptio > 1 .AND. .NOT. lk_esopa ) & 118 & CALL ctl_stop( 'only one vertical diffusion option has to be defined ' ) 119 IF( ioptio == 0 .AND. .NOT. ( lk_esopa .OR. lk_surge ) ) & 120 & CALL ctl_stop( 'a vertical diffusion option has to be defined ' ) 119 121 IF( ( lk_zdfric .OR. lk_zdfgls .OR. lk_zdfkpp ) .AND. ln_isfcav ) & 120 122 & CALL ctl_stop( ' only zdfcst and zdftke were tested with ice shelves cavities ' ) … … 151 153 ENDIF 152 154 IF ( ioptio > 1 .AND. .NOT. lk_esopa ) CALL ctl_stop( ' chose between ln_zdfnpc and ln_zdfevd' ) 153 IF( ioptio == 0 .AND. .NOT.( lk_zdftke .OR. lk_zdfgls .OR. lk_zdfkpp ) ) &154 CALL ctl_stop( ' except for TKE, GLS or KPP physics , a convection scheme is', &155 IF( ioptio == 0 .AND. .NOT.( lk_zdftke .OR. lk_zdfgls .OR. lk_zdfkpp .OR. lk_surge ) ) & 156 CALL ctl_stop( ' except for TKE, GLS or KPP physics (or running in surge mode), a convection scheme is', & 155 157 & ' required: ln_zdfevd or ln_zdfnpc logicals' ) 156 158 -
branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
r5942 r5954 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.?! 2014-09 (H. Liu) Add Wetting and Drying 32 33 !!---------------------------------------------------------------------- 33 34 … … 228 229 NAMELIST/namcfg/ cp_cfg, cp_cfz, jp_cfg, jpidta, jpjdta, jpkdta, jpiglo, jpjglo, & 229 230 & jpizoom, jpjzoom, jperio, ln_use_jattr 231 NAMELIST/namwad/ ln_wd, rn_wdmin1, rn_wdmin2, rn_wdld, nn_wdit 230 232 !!---------------------------------------------------------------------- 231 233 ! … … 253 255 READ ( numnam_cfg, namcfg, IOSTAT = ios, ERR = 904 ) 254 256 904 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcfg in configuration namelist', .TRUE. ) 257 258 REWIND( numnam_ref ) ! Namelist namwad in reference namelist : Parameters for Wetting/Drying 259 READ ( numnam_ref, namwad, IOSTAT = ios, ERR = 905) 260 905 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namwad in reference namelist', lwp ) 261 262 REWIND( numnam_cfg ) ! Namelist namwad in configuration namelist : Parameters for Wetting/Drying 263 READ ( numnam_cfg, namwad, IOSTAT = ios, ERR = 906) 264 906 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namwad in configuration namelist', lwp ) 255 265 256 266 ! Force values for AGRIF zoom (cf. agrif_user.F90) … … 311 321 WRITE( numond, namctl ) 312 322 WRITE( numond, namcfg ) 323 WRITE( numond, namwad ) 313 324 ENDIF 314 325 … … 416 427 & CALL zdf_ddm_init ! double diffusive mixing 417 428 ! ! Lateral physics 429 #if ! defined key_surge 418 430 CALL ldf_tra_init ! Lateral ocean tracer physics 431 #endif 419 432 CALL ldf_dyn_init ! Lateral ocean momentum physics 420 433 IF( lk_ldfslp ) CALL ldf_slp_init ! slope of lateral mixing 421 434 435 #if ! defined key_surge 422 436 ! ! Active tracers 423 437 CALL tra_qsr_init ! penetrative solar radiation qsr … … 428 442 CALL tra_ldf_init ! lateral mixing 429 443 CALL tra_zdf_init ! vertical mixing and after tracer fields 444 #endif 430 445 431 446 ! ! Dynamics … … 435 450 CALL dyn_ldf_init ! lateral mixing 436 451 CALL dyn_hpg_init ! horizontal gradient of Hydrostatic pressure 452 #if ! defined key_surge 437 453 CALL dyn_zdf_init ! vertical diffusion 454 #endif 438 455 CALL dyn_spg_init ! surface pressure gradient 439 456 … … 520 537 WRITE(numout,*) ' use file attribute if exists as i/p j-start ln_use_jattr = ', ln_use_jattr 521 538 ENDIF 539 540 IF(lwp) THEN ! control print 541 WRITE(numout,*) 542 WRITE(numout,*) 'namwad : wetting and drying initialisation through namelist read' 543 WRITE(numout,*) '~~~~~~~ ' 544 WRITE(numout,*) ' Namelist namwad' 545 WRITE(numout,*) ' key to turn on/off wetting/drying ln_wd = ', ln_wd 546 WRITE(numout,*) ' minimum water depth on dried cells rn_wdmin1 = ', rn_wdmin1 547 WRITE(numout,*) ' tolerrance of min water depth on dried cells rn_wdmin2 = ', rn_wdmin2 548 WRITE(numout,*) ' land elevation below which wad considered rn_wdld = ', rn_wdld 549 WRITE(numout,*) ' maximum number of iteration for W/D limiter nn_wdit = ', nn_wdit 550 ENDIF 551 522 552 ! ! Parameter control 523 553 ! -
branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/par_oce.F90
r5118 r5954 101 101 #endif 102 102 103 #if defined key_vectopt_loop 104 LOGICAL, PUBLIC, PARAMETER :: lk_vopt_loop = .TRUE. !: vector optimization flag 105 #else 106 LOGICAL, PUBLIC, PARAMETER :: lk_vopt_loop = .FALSE. !: vector optimization flag 107 #endif 108 #if defined key_surge 109 LOGICAL, PUBLIC, PARAMETER :: lk_surge = .TRUE. ! flag for 2d baratropic modelling 110 #else 111 LOGICAL, PUBLIC, PARAMETER :: lk_surge = .FALSE. ! flag for 2d baratropic modelling 112 #endif 113 103 114 !!---------------------------------------------------------------------- 104 115 !! NEMO/OPA 3.3 , NEMO Consortium (2010) -
branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/step.F90
r5510 r5954 246 246 ! Active tracers (ua, va used as workspace) 247 247 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 248 #if ! defined key_surge 248 249 tsa(:,:,:,:) = 0.e0 ! set tracer trends to zero 249 250 … … 294 295 CALL tra_nxt( kstp ) ! tracer fields at next time step 295 296 ENDIF 296 297 #endif 297 298 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 298 299 ! Dynamics (tsa used as workspace) … … 309 310 310 311 CALL dyn_bfr( kstp ) ! bottom friction 312 #if ! defined key_surge 311 313 CALL dyn_zdf( kstp ) ! vertical diffusion 314 #endif 312 315 ELSE 313 316 ua(:,:,:) = 0.e0 ! set dynamics trends to zero … … 328 331 CALL dyn_hpg( kstp ) ! horizontal gradient of Hydrostatic pressure 329 332 CALL dyn_bfr( kstp ) ! bottom friction 333 #if ! defined key_surge 330 334 CALL dyn_zdf( kstp ) ! vertical diffusion 335 #endif 331 336 CALL dyn_spg( kstp, indic ) ! surface pressure gradient 332 337 ENDIF -
branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/step_oce.F90
r5501 r5954 116 116 USE trcstp ! passive tracer time-stepping (trc_stp routine) 117 117 #endif 118 USE wadlmt ! wetting/drying limiter (wad_lmt routine) 118 119 !!---------------------------------------------------------------------- 119 120 !! NEMO/OPA 3.3 , NEMO Consortium (2010)
Note: See TracChangeset
for help on using the changeset viewer.