New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 6152 for trunk/NEMOGCM – NEMO

Changeset 6152 for trunk/NEMOGCM


Ignore:
Timestamp:
2015-12-21T23:33:57+01:00 (8 years ago)
Author:
acc
Message:

Add wetting and drying option from dev_r5803_NOC_WAD branch. Logically isolated code changes in domvvl.F90, domzgr.F90, dynhpg.F90, dynspg_ts.F90, sshwzv.F90 and nemogcm.F90. New module wet_dry.F90 in DYN. Fully SETTE tested with code deactivated (ln_wad=.false.). No test case yet available to justify activating option (still under development)

Location:
trunk/NEMOGCM
Files:
1 added
9 edited

Legend:

Unmodified
Added
Removed
  • trunk/NEMOGCM/ARCH/arch-X64_MOBILIS.fcm

    r6140 r6152  
    3838%NCDF_HOME           /home/acc/shared 
    3939%HDF5_HOME           /home/acc/shared 
    40 %XIOS_HOME           /home/acc/XIOS_1.0_r777 
     40%XIOS_HOME           /home/acc/XIOS_1.0_r803 
    4141%OASIS_HOME           
    4242 
  • trunk/NEMOGCM/CONFIG/SHARED/namelist_ref

    r6140 r6152  
    12511251/ 
    12521252!----------------------------------------------------------------------- 
     1253&namwad  !   Wetting and drying 
     1254!----------------------------------------------------------------------- 
     1255   ln_wd             = .false.  ! T/F activation of wetting and drying 
     1256   rn_wdmin1         =  0.1     ! Minimum wet depth on dried cells 
     1257   rn_wdmin2         =  0.01    ! Tolerance of min wet depth on dried cells 
     1258   rn_wdld           =  20.0    ! Land elevation below which wetting/drying is allowed 
     1259   nn_wdit           =  10      ! Max iterations for W/D limiter 
     1260/ 
     1261!----------------------------------------------------------------------- 
    12531262&nam_dia25h  !  25h Mean Output 
    12541263!----------------------------------------------------------------------- 
  • trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90

    r6140 r6152  
    2323   USE dom_oce         ! ocean space and time domain 
    2424   USE sbc_oce         ! ocean surface boundary condition 
     25   USE wet_dry         ! wetting and drying 
    2526   USE restart         ! ocean restart 
    2627   ! 
     
    687688      !                                                             !   =  'U', 'V', 'W, 'F', 'UW' or 'VW' 
    688689      ! 
    689       INTEGER ::   ji, jj, jk                                          ! dummy loop indices 
     690      INTEGER ::   ji, jj, jk                                       ! dummy loop indices 
     691      REAL(wp) ::  zlnwd                                            ! =1./0. when ln_wd = T/F 
    690692      !!---------------------------------------------------------------------- 
    691693      ! 
    692694      IF( nn_timing == 1 )   CALL timing_start('dom_vvl_interpol') 
     695      ! 
     696      IF(ln_wd) THEN 
     697        zlnwd = 1.0_wp 
     698      ELSE 
     699        zlnwd = 0.0_wp 
     700      END IF 
    693701      ! 
    694702      SELECT CASE ( pout )    !==  type of interpolation  ==! 
     
    698706            DO jj = 1, jpjm1 
    699707               DO ji = 1, fs_jpim1   ! vector opt. 
    700                   pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * r1_e1e2u(ji,jj)                                   & 
     708                  pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj)   & 
    701709                     &                       * (   e1e2t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     & 
    702710                     &                           + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) ) 
     
    711719            DO jj = 1, jpjm1 
    712720               DO ji = 1, fs_jpim1   ! vector opt. 
    713                   pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) * r1_e1e2v(ji,jj)                                   & 
     721                  pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk)  * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj)   & 
    714722                     &                       * (   e1e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     & 
    715723                     &                           + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) ) 
     
    724732            DO jj = 1, jpjm1 
    725733               DO ji = 1, fs_jpim1   ! vector opt. 
    726                   pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) * r1_e1e2f(ji,jj)               & 
     734                  pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) & 
     735                     &                       *    r1_e1e2f(ji,jj)                                                  & 
    727736                     &                       * (   e1e2u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     & 
    728737                     &                           + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) ) 
     
    739748!!gm BUG? use here wmask in case of ISF ?  to be checked 
    740749         DO jk = 2, jpk 
    741             pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * tmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) )   & 
    742                &                            +            0.5_wp * tmask(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3t_0(:,:,jk  ) ) 
     750            pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )   & 
     751               &                            * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) )                               & 
     752               &                            +            0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )     & 
     753               &                            * ( pe3_in(:,:,jk  ) - e3t_0(:,:,jk  ) ) 
    743754         END DO 
    744755         ! 
     
    749760!!gm BUG? use here wumask in case of ISF ?  to be checked 
    750761         DO jk = 2, jpk 
    751             pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * umask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) )   & 
    752                &                             +            0.5_wp * umask(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3u_0(:,:,jk  ) ) 
     762            pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )  & 
     763               &                             * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) )                              & 
     764               &                             +            0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )    & 
     765               &                             * ( pe3_in(:,:,jk  ) - e3u_0(:,:,jk  ) ) 
    753766         END DO 
    754767         ! 
     
    759772!!gm BUG? use here wvmask in case of ISF ?  to be checked 
    760773         DO jk = 2, jpk 
    761             pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * vmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) )   & 
    762                &                             +            0.5_wp * vmask(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3v_0(:,:,jk  ) ) 
     774            pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )  & 
     775               &                             * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) )                              & 
     776               &                             +            0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )    & 
     777               &                             * ( pe3_in(:,:,jk  ) - e3v_0(:,:,jk  ) ) 
    763778         END DO 
    764779      END SELECT 
     
    784799      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag 
    785800      ! 
    786       INTEGER ::   jk 
     801      INTEGER ::   ji, jj, jk 
    787802      INTEGER ::   id1, id2, id3, id4, id5     ! local integers 
    788803      !!---------------------------------------------------------------------- 
     
    872887            e3t_n(:,:,:) = e3t_0(:,:,:) 
    873888            sshn(:,:) = 0.0_wp 
     889 
     890            IF( ln_wd ) THEN 
     891              DO jj = 1, jpj 
     892                DO ji = 1, jpi 
     893                  IF( e3t_0(ji,jj,1) <= 0.5_wp * rn_wdmin1 ) THEN 
     894                     e3t_b(ji,jj,:) = 0.5_wp * rn_wdmin1  
     895                     e3t_n(ji,jj,:) = 0.5_wp * rn_wdmin1  
     896                     e3t_a(ji,jj,:) = 0.5_wp * rn_wdmin1  
     897                     sshb(ji,jj) = rn_wdmin1 - bathy(ji,jj) 
     898                     sshn(ji,jj) = rn_wdmin1 - bathy(ji,jj) 
     899                     ssha(ji,jj) = rn_wdmin1 - bathy(ji,jj) 
     900                  ENDIF 
     901                ENDDO 
     902              ENDDO 
     903            END IF 
     904 
    874905            IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN 
    875906               tilde_e3t_b(:,:,:) = 0.0_wp 
  • trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90

    r6140 r6152  
    1818   !!            3.4  ! 2012-12  (R. Bourdalle-Badie and G. Reffray)  modify C1D case   
    1919   !!            3.6  ! 2014-11  (P. Mathiot and C. Harris) add ice shelf capabilitye   
     20   !!            3.?  ! 2015-11  (H. Liu) Modifications for Wetting/Drying 
    2021   !!---------------------------------------------------------------------- 
    2122 
     
    3637   USE oce               ! ocean variables 
    3738   USE dom_oce           ! ocean domain 
     39   USE wet_dry           ! wetting and drying 
    3840   USE closea            ! closed seas 
    3941   USE c1d               ! 1D vertical configuration 
     
    19401942      bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) ) 
    19411943 
    1942       DO jj = 1, jpj 
    1943          DO ji = 1, jpi 
    1944            IF( bathy(ji,jj) > 0._wp )   bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) ) 
    1945          END DO 
    1946       END DO 
     1944      IF( .NOT.ln_wd ) THEN                       
     1945         DO jj = 1, jpj 
     1946            DO ji = 1, jpi 
     1947              IF( bathy(ji,jj) > 0._wp )   bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) ) 
     1948            END DO 
     1949         END DO 
     1950      END IF 
    19471951      !                                        ! ============================= 
    19481952      !                                        ! Define the envelop bathymetry   (hbatt) 
     
    19511955      zenv(:,:) = bathy(:,:) 
    19521956      ! 
     1957      IF( .NOT.ln_wd ) THEN     
    19531958      ! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing 
    1954       DO jj = 1, jpj 
    1955          DO ji = 1, jpi 
    1956             IF( bathy(ji,jj) == 0._wp ) THEN 
    1957                iip1 = MIN( ji+1, jpi ) 
    1958                ijp1 = MIN( jj+1, jpj ) 
    1959                iim1 = MAX( ji-1, 1 ) 
    1960                ijm1 = MAX( jj-1, 1 ) 
     1959         DO jj = 1, jpj 
     1960            DO ji = 1, jpi 
     1961               IF( bathy(ji,jj) == 0._wp ) THEN 
     1962                  iip1 = MIN( ji+1, jpi ) 
     1963                  ijp1 = MIN( jj+1, jpj ) 
     1964                  iim1 = MAX( ji-1, 1 ) 
     1965                  ijm1 = MAX( jj-1, 1 ) 
    19611966!!gm BUG fix see ticket #1617 
    1962                IF( ( + bathy(iim1,ijm1) + bathy(ji,ijp1) + bathy(iip1,ijp1)            & 
    1963                   &  + bathy(iim1,jj  )                  + bathy(iip1,jj  )            & 
    1964                   &  + bathy(iim1,ijm1) + bathy(ji,ijm1) + bathy(iip1,ijp1)  ) > 0._wp )   zenv(ji,jj) = rn_sbot_min 
     1967                  IF( ( + bathy(iim1,ijm1) + bathy(ji,ijp1) + bathy(iip1,ijp1)              & 
     1968                     &  + bathy(iim1,jj  )                  + bathy(iip1,jj  )              & 
     1969                     &  + bathy(iim1,ijm1) + bathy(ji,ijm1) + bathy(iip1,ijp1)  ) > 0._wp ) & 
     1970                     &    zenv(ji,jj) = rn_sbot_min 
    19651971!!gm 
    19661972!!gm               IF( ( bathy(iip1,jj  ) + bathy(iim1,jj  ) + bathy(ji,ijp1  ) + bathy(ji,ijm1) +         & 
     
    19691975!!gm             ENDIF 
    19701976!!gm end 
    1971            ENDIF 
    1972          END DO 
    1973       END DO 
     1977              ENDIF 
     1978            END DO 
     1979         END DO 
     1980      END IF 
     1981 
    19741982      ! apply lateral boundary condition   CAUTION: keep the value when the lbc field is zero 
    19751983      CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' ) 
     
    20642072      IF(lwp) THEN 
    20652073         WRITE(numout,*) 
    2066          WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min 
     2074         IF( .NOT.ln_wd ) THEN 
     2075           WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min 
     2076         ELSE 
     2077           WRITE(numout,*) ' zgr_sco: minimum positive depth of the envelop topography set to : ', rn_sbot_min 
     2078           WRITE(numout,*) ' zgr_sco: minimum negative depth of the envelop topography set to : ', -rn_wdld 
     2079         ENDIF 
    20672080      ENDIF 
    20682081      hbatu(:,:) = rn_sbot_min 
     
    20772090        END DO 
    20782091      END DO 
     2092 
     2093      IF( ln_wd ) THEN               !avoid the zero depth on T- (U-,V-,F-) points 
     2094        DO jj = 1, jpj 
     2095          DO ji = 1, jpi 
     2096            IF(ABS(hbatt(ji,jj)) < rn_wdmin1) & 
     2097              & hbatt(ji,jj) = SIGN(1._wp, hbatt(ji,jj)) * rn_wdmin1 
     2098            IF(ABS(hbatu(ji,jj)) < rn_wdmin1) & 
     2099              & hbatu(ji,jj) = SIGN(1._wp, hbatu(ji,jj)) * rn_wdmin1 
     2100            IF(ABS(hbatv(ji,jj)) < rn_wdmin1) & 
     2101              & hbatv(ji,jj) = SIGN(1._wp, hbatv(ji,jj)) * rn_wdmin1 
     2102            IF(ABS(hbatf(ji,jj)) < rn_wdmin1) & 
     2103              & hbatf(ji,jj) = SIGN(1._wp, hbatf(ji,jj)) * rn_wdmin1 
     2104          END DO 
     2105        END DO 
     2106      END IF 
    20792107      !  
    20802108      ! Apply lateral boundary condition 
     
    20842112         DO ji = 1, jpi 
    20852113            IF( hbatu(ji,jj) == 0._wp ) THEN 
     2114               !No worries about the following line when ln_wd == .true. 
    20862115               IF( zhbat(ji,jj) == 0._wp )   hbatu(ji,jj) = rn_sbot_min 
    20872116               IF( zhbat(ji,jj) /= 0._wp )   hbatu(ji,jj) = zhbat(ji,jj) 
     
    21092138 
    21102139!!bug:  key_helsinki a verifer 
    2111       hift(:,:) = MIN( hift(:,:), hbatt(:,:) ) 
    2112       hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) ) 
    2113       hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) ) 
    2114       hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) ) 
     2140      IF( .NOT.ln_wd ) THEN 
     2141        hift(:,:) = MIN( hift(:,:), hbatt(:,:) ) 
     2142        hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) ) 
     2143        hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) ) 
     2144        hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) ) 
     2145      END IF 
    21152146 
    21162147      IF( nprint == 1 .AND. lwp )   THEN 
     
    21542185      CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 
    21552186      ! 
    2156       WHERE( e3t_0 (:,:,:) == 0._wp )   e3t_0 (:,:,:) = 1._wp 
    2157       WHERE( e3u_0 (:,:,:) == 0._wp )   e3u_0 (:,:,:) = 1._wp 
    2158       WHERE( e3v_0 (:,:,:) == 0._wp )   e3v_0 (:,:,:) = 1._wp 
    2159       WHERE( e3f_0 (:,:,:) == 0._wp )   e3f_0 (:,:,:) = 1._wp 
    2160       WHERE( e3w_0 (:,:,:) == 0._wp )   e3w_0 (:,:,:) = 1._wp 
    2161       WHERE( e3uw_0(:,:,:) == 0._wp )   e3uw_0(:,:,:) = 1._wp 
    2162       WHERE( e3vw_0(:,:,:) == 0._wp )   e3vw_0(:,:,:) = 1._wp 
     2187      IF( .NOT.ln_wd ) THEN 
     2188        WHERE( e3t_0 (:,:,:) == 0._wp )   e3t_0 (:,:,:) = 1._wp 
     2189        WHERE( e3u_0 (:,:,:) == 0._wp )   e3u_0 (:,:,:) = 1._wp 
     2190        WHERE( e3v_0 (:,:,:) == 0._wp )   e3v_0 (:,:,:) = 1._wp 
     2191        WHERE( e3f_0 (:,:,:) == 0._wp )   e3f_0 (:,:,:) = 1._wp 
     2192        WHERE( e3w_0 (:,:,:) == 0._wp )   e3w_0 (:,:,:) = 1._wp 
     2193        WHERE( e3uw_0(:,:,:) == 0._wp )   e3uw_0(:,:,:) = 1._wp 
     2194        WHERE( e3vw_0(:,:,:) == 0._wp )   e3vw_0(:,:,:) = 1._wp 
     2195      END IF 
    21632196 
    21642197#if defined key_agrif 
     
    21932226               IF( scobot(ji,jj) >= gdept_n(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk ) 
    21942227            END DO 
    2195             IF( scobot(ji,jj) == 0._wp                )   mbathy(ji,jj) = 0 
     2228            IF( ln_wd ) THEN 
     2229              IF( scobot(ji,jj) <= -(rn_wdld - rn_wdmin2)) THEN 
     2230                mbathy(ji,jj) = 0 
     2231              ELSEIF(scobot(ji,jj) <= rn_wdmin1) THEN 
     2232                mbathy(ji,jj) = 2 
     2233              ENDIF 
     2234            ELSE 
     2235              IF( scobot(ji,jj) == 0._wp )   mbathy(ji,jj) = 0 
     2236            ENDIF 
    21962237         END DO 
    21972238      END DO 
     
    23112352      INTEGER  ::   ji, jj, jk           ! dummy loop argument 
    23122353      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars 
     2354      REAL(wp) ::   ztmpu,  ztmpv,  ztmpf 
     2355      REAL(wp) ::   ztmpu1, ztmpv1, ztmpf1 
    23132356      ! 
    23142357      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 
     
    23652408      DO ji = 1, jpim1 
    23662409         DO jj = 1, jpjm1 
     2410            ! extended for Wetting/Drying case 
     2411            ztmpu  = hbatt(ji,jj)+hbatt(ji+1,jj) 
     2412            ztmpv  = hbatt(ji,jj)+hbatt(ji,jj+1) 
     2413            ztmpf  = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) 
     2414            ztmpu1 = hbatt(ji,jj)*hbatt(ji+1,jj) 
     2415            ztmpv1 = hbatt(ji,jj)*hbatt(ji,jj+1) 
     2416            ztmpf1 = MIN(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) * & 
     2417                   & MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) 
    23672418            DO jk = 1, jpk 
    2368                z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) )   & 
    2369                   &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
    2370                z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) )   & 
    2371                   &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
    2372                z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk)     & 
    2373                   &                + hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) )   & 
    2374                   &              / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 
    2375                z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) )   & 
    2376                   &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
    2377                z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) )   & 
    2378                   &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
     2419               IF( ln_wd .AND. (ztmpu1 < 0._wp.OR.ABS(ztmpu) < rn_wdmin1) ) THEN 
     2420                 z_esigtu3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) ) 
     2421                 z_esigwu3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji+1,jj,jk) ) 
     2422               ELSE 
     2423                 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) & 
     2424                        &              / ztmpu 
     2425                 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) & 
     2426                        &              / ztmpu 
     2427               END IF 
     2428 
     2429               IF( ln_wd .AND. (ztmpv1 < 0._wp.OR.ABS(ztmpv) < rn_wdmin1) ) THEN 
     2430                 z_esigtv3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji,jj+1,jk) ) 
     2431                 z_esigwv3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji,jj+1,jk) ) 
     2432               ELSE 
     2433                 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) & 
     2434                        &              / ztmpv 
     2435                 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) & 
     2436                        &              / ztmpv 
     2437               END IF 
     2438 
     2439               IF( ln_wd .AND. (ztmpf1 < 0._wp.OR.ABS(ztmpf) < rn_wdmin1) ) THEN 
     2440                 z_esigtf3(ji,jj,jk) = 0.25_wp * ( z_esigt3(ji,jj  ,jk) + z_esigt3(ji+1,jj  ,jk)               & 
     2441                        &                        + z_esigt3(ji,jj+1,jk) + z_esigt3(ji+1,jj+1,jk) ) 
     2442               ELSE 
     2443                 z_esigtf3(ji,jj,jk) = ( hbatt(ji  ,jj  )*z_esigt3(ji  ,jj  ,jk)                               & 
     2444                        &            +   hbatt(ji+1,jj  )*z_esigt3(ji+1,jj  ,jk)                               & 
     2445                        &            +   hbatt(ji  ,jj+1)*z_esigt3(ji  ,jj+1,jk)                               & 
     2446                        &            +   hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / ztmpf 
     2447               END IF 
     2448 
    23792449               ! 
    23802450               e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 
     
    24152485      REAL(wp) ::   zsmth               ! smoothing around critical depth 
    24162486      REAL(wp) ::   zzs, zzb           ! Surface and bottom cell thickness in sigma space 
     2487      REAL(wp) ::   ztmpu, ztmpv, ztmpf 
     2488      REAL(wp) ::   ztmpu1, ztmpv1, ztmpf1 
    24172489      ! 
    24182490      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 
     
    24932565        DO jj=1,jpj-1 
    24942566 
    2495           DO jk = 1, jpk 
    2496                 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) / & 
    2497                                     ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
    2498                 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) / & 
    2499                                     ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
    2500                 z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) +  & 
    2501                                       hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / & 
    2502                                     ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 
    2503                 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) / & 
    2504                                     ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
    2505                 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) / & 
    2506                                     ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
     2567           ! extend to suit for Wetting/Drying case 
     2568           ztmpu  = hbatt(ji,jj)+hbatt(ji+1,jj) 
     2569           ztmpv  = hbatt(ji,jj)+hbatt(ji,jj+1) 
     2570           ztmpf  = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) 
     2571           ztmpu1 = hbatt(ji,jj)*hbatt(ji+1,jj) 
     2572           ztmpv1 = hbatt(ji,jj)*hbatt(ji,jj+1) 
     2573           ztmpf1 = MIN(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) * & 
     2574                  & MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) 
     2575           DO jk = 1, jpk 
     2576              IF( ln_wd .AND. (ztmpu1 < 0._wp.OR.ABS(ztmpu) < rn_wdmin1) ) THEN 
     2577                z_esigtu3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) ) 
     2578                z_esigwu3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji+1,jj,jk) ) 
     2579              ELSE 
     2580                z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) & 
     2581                       &              / ztmpu 
     2582                z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) & 
     2583                       &              / ztmpu 
     2584              END IF 
     2585 
     2586              IF( ln_wd .AND. (ztmpv1 < 0._wp.OR.ABS(ztmpv) < rn_wdmin1) ) THEN 
     2587                z_esigtv3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji,jj+1,jk) ) 
     2588                z_esigwv3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji,jj+1,jk) ) 
     2589              ELSE 
     2590                z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) & 
     2591                       &              / ztmpv 
     2592                z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) & 
     2593                       &              / ztmpv 
     2594              END IF 
     2595 
     2596              IF( ln_wd .AND. (ztmpf1 < 0._wp.OR.ABS(ztmpf) < rn_wdmin1) ) THEN 
     2597                z_esigtf3(ji,jj,jk) = 0.25_wp * ( z_esigt3(ji,jj,jk)   + z_esigt3(ji+1,jj,jk)                 & 
     2598                       &                        + z_esigt3(ji,jj+1,jk) + z_esigt3(ji+1,jj+1,jk) ) 
     2599              ELSE 
     2600                z_esigtf3(ji,jj,jk) = ( hbatt(ji  ,jj  )*z_esigt3(ji  ,jj  ,jk)                               & 
     2601                       &              + hbatt(ji+1,jj  )*z_esigt3(ji+1,jj  ,jk)                               & 
     2602                       &              + hbatt(ji  ,jj+1)*z_esigt3(ji  ,jj+1,jk)                               & 
     2603                       &              + hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / ztmpf 
     2604              END IF 
     2605 
     2606             ! Code prior to wetting and drying option (for reference) 
     2607             !z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) )   & 
     2608             !                     /( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
     2609             ! 
     2610             !z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) )   & 
     2611             !                     /( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
     2612             ! 
     2613             !z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) )   & 
     2614             !                     /( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
     2615             ! 
     2616             !z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) )   & 
     2617             !                     /( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
     2618             ! 
     2619             !z_esigtf3(ji,jj,jk) = ( hbatt(ji  ,jj  )*z_esigt3(ji  ,jj  ,jk)                                 & 
     2620             !                    &  +hbatt(ji+1,jj  )*z_esigt3(ji+1,jj  ,jk)                                 & 
     2621             !                       +hbatt(ji  ,jj+1)*z_esigt3(ji  ,jj+1,jk)                                 & 
     2622             !                    &  +hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) )                               & 
     2623             !                     /( hbatt(ji  ,jj  )+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 
    25072624 
    25082625             e3t_0(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk) 
  • trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90

    r6140 r6152  
    3333   USE sbc_oce         ! surface variable (only for the flag with ice shelf) 
    3434   USE dom_oce         ! ocean space and time domain 
     35   USE wet_dry         ! wetting and drying 
    3536   USE phycst          ! physical constants 
    3637   USE trd_oce         ! trends: ocean variables 
     
    429430      INTEGER, INTENT(in) ::   kt    ! ocean time-step index 
    430431      !! 
    431       INTEGER  ::   ji, jj, jk                 ! dummy loop indices 
    432       REAL(wp) ::   zcoef0, zuap, zvap, znad   ! temporary scalars 
     432      INTEGER  ::   ji, jj, jk, jii, jjj                 ! dummy loop indices 
     433      REAL(wp) ::   zcoef0, zuap, zvap, znad, ztmp       ! temporary scalars 
     434      LOGICAL  ::   ll_tmp1, ll_tmp2, ll_tmp3            ! local logical variables 
    433435      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zhpi, zhpj 
     436      REAL(wp), POINTER, DIMENSION(:,:)   ::  zcpx, zcpy !W/D pressure filter 
    434437      !!---------------------------------------------------------------------- 
    435438      ! 
    436439      CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj ) 
     440      IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 
    437441      ! 
    438442      IF( kt == nit000 ) THEN 
     
    447451      ENDIF 
    448452      ! 
     453      IF(ln_wd) THEN 
     454        DO jj = 2, jpjm1 
     455           DO ji = 2, jpim1  
     456             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj))  
     457             ll_tmp2 = MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) > rn_wdmin1 + rn_wdmin2 
     458             ll_tmp3 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) + & 
     459                                                       & rn_wdmin1 + rn_wdmin2 
     460 
     461             IF(ll_tmp1.AND.ll_tmp2) THEN 
     462               zcpx(ji,jj) = 1.0_wp 
     463               wduflt(ji,jj) = 1.0_wp 
     464             ELSE IF(ll_tmp3) THEN 
     465               ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 
     466               zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) / & 
     467                           &     (sshn(ji+1,jj) - sshn(ji,jj))) 
     468               wduflt(ji,jj) = 1.0_wp 
     469             ELSE 
     470               zcpx(ji,jj) = 0._wp 
     471               wduflt(ji,jj) = 0.0_wp 
     472             END IF 
     473       
     474             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1))  
     475             ll_tmp2 = MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) > rn_wdmin1 + rn_wdmin2 
     476             ll_tmp3 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) + & 
     477                                                       & rn_wdmin1 + rn_wdmin2 
     478 
     479             IF(ll_tmp1.AND.ll_tmp2) THEN 
     480               zcpy(ji,jj) = 1.0_wp 
     481               wdvflt(ji,jj) = 1.0_wp 
     482             ELSE IF(ll_tmp3) THEN 
     483               ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 
     484               zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) / & 
     485                           &     (sshn(ji,jj+1) - sshn(ji,jj))) 
     486               wdvflt(ji,jj) = 1.0_wp 
     487             ELSE 
     488               zcpy(ji,jj) = 0._wp 
     489               wdvflt(ji,jj) = 0.0_wp 
     490             END IF 
     491           END DO 
     492        END DO 
     493        CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp ) 
     494      ENDIF 
     495 
     496 
    449497      ! Surface value 
    450498      DO jj = 2, jpjm1 
     
    460508            zvap = -zcoef0 * ( rhd    (ji,jj+1,1) + rhd    (ji,jj,1) + 2._wp * znad )   & 
    461509               &           * ( gde3w_n(ji,jj+1,1) - gde3w_n(ji,jj,1) ) * r1_e2v(ji,jj) 
     510 
     511 
     512            IF(ln_wd) THEN 
     513 
     514              zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 
     515              zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj)  
     516              zuap = zuap * zcpx(ji,jj) 
     517              zvap = zvap * zcpy(ji,jj) 
     518            ENDIF 
     519 
    462520            ! add to the general momentum trend 
    463521            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap 
     
    482540               zvap = -zcoef0 * ( rhd    (ji  ,jj+1,jk) + rhd    (ji,jj,jk) + 2._wp * znad )   & 
    483541                  &           * ( gde3w_n(ji  ,jj+1,jk) - gde3w_n(ji,jj,jk) ) * r1_e2v(ji,jj) 
     542 
     543               IF(ln_wd) THEN 
     544                 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 
     545                 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj)  
     546                 zuap = zuap * zcpx(ji,jj) 
     547                 zvap = zvap * zcpy(ji,jj) 
     548               ENDIF 
     549 
    484550               ! add to the general momentum trend 
    485551               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap 
     
    490556      ! 
    491557      CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj ) 
     558      IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 
    492559      ! 
    493560   END SUBROUTINE hpg_sco 
     
    623690      REAL(wp) ::   z1_10, cffu, cffx   !    "         " 
    624691      REAL(wp) ::   z1_12, cffv, cffy   !    "         " 
     692      LOGICAL  ::   ll_tmp1, ll_tmp2    ! local logical variables 
    625693      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zhpi, zhpj 
    626694      REAL(wp), POINTER, DIMENSION(:,:,:) ::  dzx, dzy, dzz, dzu, dzv, dzw 
    627695      REAL(wp), POINTER, DIMENSION(:,:,:) ::  drhox, drhoy, drhoz, drhou, drhov, drhow 
    628696      REAL(wp), POINTER, DIMENSION(:,:,:) ::  rho_i, rho_j, rho_k 
     697      REAL(wp), POINTER, DIMENSION(:,:)   ::  zcpx, zcpy    !W/D pressure filter 
    629698      !!---------------------------------------------------------------------- 
    630699      ! 
     
    632701      CALL wrk_alloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow ) 
    633702      CALL wrk_alloc( jpi, jpj, jpk, rho_i, rho_j, rho_k,  zhpi,  zhpj        ) 
    634       ! 
     703      IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 
     704      ! 
     705      ! 
     706      IF(ln_wd) THEN 
     707        DO jj = 2, jpjm1 
     708           DO ji = 2, jpim1  
     709             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) & 
     710                     & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) & 
     711                     &  > rn_wdmin1 + rn_wdmin2 
     712             ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) +& 
     713                                                       & rn_wdmin1 + rn_wdmin2 
     714 
     715             IF(ll_tmp1) THEN 
     716               zcpx(ji,jj) = 1.0_wp 
     717             ELSE IF(ll_tmp2) THEN 
     718               ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 
     719               zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /& 
     720                           &     (sshn(ji+1,jj) - sshn(ji,jj))) 
     721             ELSE 
     722               zcpx(ji,jj) = 0._wp 
     723             END IF 
     724       
     725             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) & 
     726                     & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) & 
     727                     &  > rn_wdmin1 + rn_wdmin2 
     728             ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) +& 
     729                                                       & rn_wdmin1 + rn_wdmin2 
     730 
     731             IF(ll_tmp1) THEN 
     732               zcpy(ji,jj) = 1.0_wp 
     733             ELSE IF(ll_tmp2) THEN 
     734               ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 
     735               zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) /& 
     736                           &     (sshn(ji,jj+1) - sshn(ji,jj))) 
     737             ELSE 
     738               zcpy(ji,jj) = 0._wp 
     739             END IF 
     740           END DO 
     741        END DO 
     742        CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp ) 
     743      ENDIF 
     744 
    635745 
    636746      IF( kt == nit000 ) THEN 
     
    803913            zhpi(ji,jj,1) = ( rho_k(ji+1,jj  ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) * r1_e1u(ji,jj) 
    804914            zhpj(ji,jj,1) = ( rho_k(ji  ,jj+1,1) - rho_k(ji,jj,1) - rho_j(ji,jj,1) ) * r1_e2v(ji,jj) 
     915            IF(ln_wd) THEN 
     916              zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 
     917              zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj)  
     918            ENDIF 
    805919            ! add to the general momentum trend 
    806920            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) 
     
    822936                  &           + (  ( rho_k(ji,jj+1,jk) - rho_k(ji,jj,jk  ) )    & 
    823937                  &               -( rho_j(ji,jj  ,jk) - rho_j(ji,jj,jk-1) )  ) * r1_e2v(ji,jj) 
     938               IF(ln_wd) THEN 
     939                 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 
     940                 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj)  
     941               ENDIF 
    824942               ! add to the general momentum trend 
    825943               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) 
     
    832950      CALL wrk_dealloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow ) 
    833951      CALL wrk_dealloc( jpi, jpj, jpk, rho_i, rho_j, rho_k,  zhpi,  zhpj        ) 
     952      IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 
    834953      ! 
    835954   END SUBROUTINE hpg_djc 
     
    855974      !! The local variables for the correction term 
    856975      INTEGER  :: jk1, jis, jid, jjs, jjd 
     976      LOGICAL  :: ll_tmp1, ll_tmp2                  ! local logical variables 
    857977      REAL(wp) :: zuijk, zvijk, zpwes, zpwed, zpnss, zpnsd, zdeps 
    858978      REAL(wp) :: zrhdt1 
     
    861981      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp 
    862982      REAL(wp), POINTER, DIMENSION(:,:)   ::   zsshu_n, zsshv_n 
     983      REAL(wp), POINTER, DIMENSION(:,:)   ::  zcpx, zcpy    !W/D pressure filter 
    863984      !!---------------------------------------------------------------------- 
    864985      ! 
     
    866987      CALL wrk_alloc( jpi,jpj,jpk,   zdept, zrhh ) 
    867988      CALL wrk_alloc( jpi,jpj,       zsshu_n, zsshv_n ) 
     989      IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 
    868990      ! 
    869991      IF( kt == nit000 ) THEN 
     
    877999      znad = 1._wp 
    8781000      IF( ln_linssh )   znad = 0._wp 
     1001 
     1002      IF(ln_wd) THEN 
     1003        DO jj = 2, jpjm1 
     1004           DO ji = 2, jpim1  
     1005             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) & 
     1006                     & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) & 
     1007                     &  > rn_wdmin1 + rn_wdmin2 
     1008             ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) +& 
     1009                                                       & rn_wdmin1 + rn_wdmin2 
     1010 
     1011             IF(ll_tmp1) THEN 
     1012               zcpx(ji,jj) = 1.0_wp 
     1013             ELSE IF(ll_tmp2) THEN 
     1014               ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 
     1015               zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /& 
     1016                           &     (sshn(ji+1,jj) - sshn(ji,jj))) 
     1017             ELSE 
     1018               zcpx(ji,jj) = 0._wp 
     1019             END IF 
     1020       
     1021             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) & 
     1022                     & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) & 
     1023                     &  > rn_wdmin1 + rn_wdmin2 
     1024             ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) +& 
     1025                                                       & rn_wdmin1 + rn_wdmin2 
     1026 
     1027             IF(ll_tmp1.OR.ll_tmp2) THEN 
     1028               zcpy(ji,jj) = 1.0_wp 
     1029             ELSE IF(ll_tmp2) THEN 
     1030               ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 
     1031               zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) /& 
     1032                           &     (sshn(ji,jj+1) - sshn(ji,jj))) 
     1033             ELSE 
     1034               zcpy(ji,jj) = 0._wp 
     1035             END IF 
     1036           END DO 
     1037        END DO 
     1038        CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp ) 
     1039      ENDIF 
    8791040 
    8801041      ! Clean 3-D work arrays 
     
    9601121        END DO 
    9611122      END DO 
     1123 
     1124      CALL lbc_lnk (zsshu_n, 'U', 1.) 
     1125      CALL lbc_lnk (zsshv_n, 'V', 1.) 
    9621126 
    9631127      DO jj = 2, jpjm1 
     
    10571221                 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed) 
    10581222               ENDIF 
    1059 !!gm  Since umask(ji,:,:) = tmask(ji,:,:) * tmask(ji+1,:,:)  by definition 
    1060 !!gm      in the line below only * umask(ji,jj,jk)  is needed !! 
    1061                ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk) * tmask(ji,jj,jk) * tmask(ji+1,jj,jk) 
     1223               IF(ln_wd) THEN 
     1224                  zdpdx1 = zdpdx1 * zcpx(ji,jj) 
     1225                  zdpdx2 = zdpdx2 * zcpx(ji,jj) 
     1226               ENDIF 
     1227               ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk)  
    10621228            ENDIF 
    10631229 
     
    11141280                  zdpdy2 = zcoef0 * r1_e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd ) 
    11151281               ENDIF 
    1116 !!gm  Since vmask(:,jj,:) = tmask(:,jj,:) * tmask(:,jj+1,:)  by definition 
    1117 !!gm      in the line below only * vmask(ji,jj,jk)  is needed !! 
    1118                va(ji,jj,jk) = va(ji,jj,jk) + ( zdpdy1 + zdpdy2 ) * vmask(ji,jj,jk)*tmask(ji,jj,jk)*tmask(ji,jj+1,jk) 
     1282               IF(ln_wd) THEN 
     1283                  zdpdy1 = zdpdy1 * zcpy(ji,jj) 
     1284                  zdpdy2 = zdpdy2 * zcpy(ji,jj) 
     1285               ENDIF 
     1286 
     1287               va(ji,jj,jk) = va(ji,jj,jk) + (zdpdy1 + zdpdy2) * vmask(ji,jj,jk) 
    11191288            ENDIF 
    11201289               ! 
     
    11261295      CALL wrk_dealloc( jpi,jpj,jpk,   zdept, zrhh ) 
    11271296      CALL wrk_dealloc( jpi,jpj,       zsshu_n, zsshv_n ) 
     1297      IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 
    11281298      ! 
    11291299   END SUBROUTINE hpg_prj 
  • trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r6140 r6152  
    3232   USE phycst          ! physical constants 
    3333   USE dynvor          ! vorticity term 
     34   USE wet_dry         ! wetting/drying flux limter 
    3435   USE bdy_par         ! for lk_bdy 
    3536   USE bdytides        ! open boundary condition data 
     
    136137      LOGICAL  ::   ll_fw_start        ! if true, forward integration  
    137138      LOGICAL  ::   ll_init             ! if true, special startup of 2d equations 
     139      LOGICAL  ::   ll_tmp1, ll_tmp2            ! local logical variables used in W/D 
    138140      INTEGER  ::   ji, jj, jk, jn        ! dummy loop indices 
    139141      INTEGER  ::   ikbu, ikbv, noffset      ! local integers 
     
    153155      REAL(wp), POINTER, DIMENSION(:,:) :: zsshu_a, zsshv_a 
    154156      REAL(wp), POINTER, DIMENSION(:,:) :: zhf 
     157      REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy                 ! Wetting/Dying gravity filter coef. 
     158      REAL(wp), POINTER, DIMENSION(:,:) :: wduflt1, wdvflt1           ! Wetting/Dying velocity filter coef. 
    155159      !!---------------------------------------------------------------------- 
    156160      ! 
     
    162166      CALL wrk_alloc( jpi,jpj,   zwx, zwy, zssh_frc, zu_frc, zv_frc) 
    163167      CALL wrk_alloc( jpi,jpj,   zhup2_e, zhvp2_e, zhust_e, zhvst_e) 
    164       CALL wrk_alloc( jpi,jpj,   zsshu_a, zsshv_a                                   ) 
     168      CALL wrk_alloc( jpi,jpj,   zsshu_a, zsshv_a                  ) 
    165169      CALL wrk_alloc( jpi,jpj,   zhf ) 
     170      IF( ln_wd ) CALL wrk_alloc( jpi, jpj, zcpx, zcpy, wduflt1, wdvflt1 ) 
    166171      ! 
    167172      zmdi=1.e+20                               !  missing data indicator for masking 
     
    368373      !                                   ! ---------------------------------------------------- 
    369374      IF( .NOT.ln_linssh ) THEN                 ! Variable volume : remove surface pressure gradient 
    370          DO jj = 2, jpjm1  
    371             DO ji = fs_2, fs_jpim1   ! vector opt. 
    372                zu_trd(ji,jj) = zu_trd(ji,jj) - grav * (  sshn(ji+1,jj  ) - sshn(ji  ,jj  )  ) * r1_e1u(ji,jj) 
    373                zv_trd(ji,jj) = zv_trd(ji,jj) - grav * (  sshn(ji  ,jj+1) - sshn(ji  ,jj  )  ) * r1_e2v(ji,jj) 
    374             END DO 
    375          END DO 
     375        IF( ln_wd ) THEN                        ! Calculating and applying W/D gravity filters 
     376          wduflt1(:,:) = 1.0_wp 
     377          wdvflt1(:,:) = 1.0_wp 
     378          DO jj = 2, jpjm1 
     379             DO ji = 2, jpim1 
     380                ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj))   & 
     381                        & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj))   & 
     382                        &  > rn_wdmin1 + rn_wdmin2 
     383                ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj))   & 
     384                        &                                   + rn_wdmin1 + rn_wdmin2 
     385                IF(ll_tmp1) THEN 
     386                  zcpx(ji,jj)    = 1.0_wp 
     387                ELSEIF(ll_tmp2) THEN 
     388                   ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen here 
     389                  zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) & 
     390                        &          /(sshn(ji+1,jj) - sshn(ji,jj))) 
     391                ELSE 
     392                  zcpx(ji,jj)    = 0._wp 
     393                  wduflt1(ji,jj) = 0.0_wp 
     394                END IF 
     395 
     396                ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1))   & 
     397                        & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1))   & 
     398                        &  > rn_wdmin1 + rn_wdmin2 
     399                ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1))   & 
     400                        &                                   + rn_wdmin1 + rn_wdmin2 
     401                IF(ll_tmp1) THEN 
     402                   zcpy(ji,jj)    = 1.0_wp 
     403                ELSEIF(ll_tmp2) THEN 
     404                   ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen here 
     405                  zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) & 
     406                        &          /(sshn(ji,jj+1) - sshn(ji,jj))) 
     407                ELSE 
     408                  zcpy(ji,jj)    = 0._wp 
     409                  wdvflt1(ji,jj) = 0.0_wp 
     410                ENDIF 
     411 
     412             END DO 
     413           END DO 
     414 
     415           CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp ) 
     416 
     417           DO jj = 2, jpjm1 
     418              DO ji = 2, jpim1 
     419                 zu_trd(ji,jj) = ( zu_trd(ji,jj) - grav * ( sshn(ji+1,jj  ) - sshn(ji  ,jj ) )   & 
     420                        &                        * r1_e1u(ji,jj) ) * zcpx(ji,jj) * wduflt1(ji,jj) 
     421                 zv_trd(ji,jj) = ( zv_trd(ji,jj) - grav * ( sshn(ji  ,jj+1) - sshn(ji  ,jj ) )   & 
     422                        &                        * r1_e2v(ji,jj) ) * zcpy(ji,jj) * wdvflt1(ji,jj) 
     423              END DO 
     424           END DO 
     425 
     426         ELSE 
     427 
     428           DO jj = 2, jpjm1 
     429              DO ji = fs_2, fs_jpim1   ! vector opt. 
     430                 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * (  sshn(ji+1,jj  ) - sshn(ji  ,jj  )  ) * r1_e1u(ji,jj) 
     431                 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * (  sshn(ji  ,jj+1) - sshn(ji  ,jj  )  ) * r1_e2v(ji,jj)  
     432              END DO 
     433           END DO 
     434        ENDIF 
     435 
    376436      ENDIF 
    377437 
     
    405465      ! 
    406466      ! Note that the "unclipped" bottom friction parameter is used even with explicit drag 
    407       zu_frc(:,:) = zu_frc(:,:) + r1_hu_n(:,:) * bfrua(:,:) * zwx(:,:) 
    408       zv_frc(:,:) = zv_frc(:,:) + r1_hv_n(:,:) * bfrva(:,:) * zwy(:,:) 
    409       !        
     467      IF( ln_wd ) THEN 
     468        zu_frc(:,:) = zu_frc(:,:) + MAX(r1_hu_n(:,:) * bfrua(:,:),-1._wp / rdtbt) * zwx(:,:) 
     469        zv_frc(:,:) = zv_frc(:,:) + MAX(r1_hv_n(:,:) * bfrva(:,:),-1._wp / rdtbt) * zwy(:,:) 
     470      ELSE 
     471        zu_frc(:,:) = zu_frc(:,:) + r1_hu_n(:,:) * bfrua(:,:) * zwx(:,:) 
     472        zv_frc(:,:) = zv_frc(:,:) + r1_hv_n(:,:) * bfrva(:,:) * zwy(:,:) 
     473      END IF 
     474      ! 
    410475      !                                         ! Add top stress contribution from baroclinic velocities:       
    411476      IF (ln_bt_fw) THEN 
     
    500565         ub_e   (:,:) = 0._wp 
    501566         vb_e   (:,:) = 0._wp 
     567      ENDIF 
     568 
     569      IF( ln_wd ) THEN      !preserve the positivity of water depth 
     570                          !ssh[b,n,a] should have already been processed for this 
     571         sshbb_e(:,:) = MAX(sshbb_e(:,:), rn_wdmin1 - bathy(:,:)) 
     572         sshb_e(:,:)  = MAX(sshb_e(:,:) , rn_wdmin1 - bathy(:,:)) 
    502573      ENDIF 
    503574      ! 
     
    575646            zhup2_e (:,:) = hu_0(:,:) + zwx(:,:)                ! Ocean depth at U- and V-points 
    576647            zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:) 
     648            IF( ln_wd ) THEN 
     649              zhup2_e(:,:) = MAX(zhup2_e (:,:), rn_wdmin1) 
     650              zhvp2_e(:,:) = MAX(zhvp2_e (:,:), rn_wdmin1) 
     651            END IF 
    577652         ELSE 
    578653            zhup2_e (:,:) = hu_n(:,:) 
     
    612687         ENDIF 
    613688#endif 
     689         IF( ln_wd ) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt) 
    614690         ! 
    615691         ! Sum over sub-time-steps to compute advective velocities 
     
    626702         END DO 
    627703         ssha_e(:,:) = (  sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) )  ) * ssmask(:,:) 
     704         IF( ln_wd ) ssha_e(:,:) = MAX(ssha_e(:,:), rn_wdmin1 - bathy(:,:))  
    628705         CALL lbc_lnk( ssha_e, 'T',  1._wp ) 
    629706 
     
    672749         zsshp2_e(:,:) = za0 *  ssha_e(:,:) + za1 *  sshn_e (:,:) & 
    673750          &            + za2 *  sshb_e(:,:) + za3 *  sshbb_e(:,:) 
     751         IF( ln_wd ) THEN                   ! Calculating and applying W/D gravity filters 
     752           wduflt1(:,:) = 1._wp 
     753           wdvflt1(:,:) = 1._wp 
     754           DO jj = 2, jpjm1 
     755              DO ji = 2, jpim1 
     756                 ll_tmp1 = MIN( zsshp2_e(ji,jj), zsshp2_e(ji+1,jj) ) > MAX( -bathy(ji,jj), -bathy(ji+1,jj) ) & 
     757                        & .AND. MAX( zsshp2_e(ji,jj) + bathy(ji,jj), zsshp2_e(ji+1,jj) + bathy(ji+1,jj) )    & 
     758                        &                                  > rn_wdmin1 + rn_wdmin2 
     759                 ll_tmp2 = MAX( zsshp2_e(ji,jj), zsshp2_e(ji+1,jj) ) > MAX( -bathy(ji,jj), -bathy(ji+1,jj) ) & 
     760                        &                                  + rn_wdmin1 + rn_wdmin2 
     761                 IF(ll_tmp1) THEN 
     762                    zcpx(ji,jj) = 1._wp 
     763                 ELSE IF(ll_tmp2) THEN 
     764                    ! no worries about zsshp2_e(ji+1,jj)-zsshp2_e(ji,jj) = 0, it won't happen here 
     765                    zcpx(ji,jj) = ABS( (zsshp2_e(ji+1,jj) + bathy(ji+1,jj) - zsshp2_e(ji,jj) - bathy(ji,jj)) & 
     766                        &             / (zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj)) ) 
     767                 ELSE 
     768                    zcpx(ji,jj)    = 0._wp 
     769                    wduflt1(ji,jj) = 0._wp 
     770                 END IF 
     771 
     772                 ll_tmp1 = MIN( zsshp2_e(ji,jj), zsshp2_e(ji,jj+1) ) > MAX( -bathy(ji,jj), -bathy(ji,jj+1) ) & 
     773                        & .AND. MAX( zsshp2_e(ji,jj) + bathy(ji,jj), zsshp2_e(ji,jj+1) + bathy(ji,jj+1) )    & 
     774                        &                                  > rn_wdmin1 + rn_wdmin2 
     775                 ll_tmp2 = MAX( zsshp2_e(ji,jj), zsshp2_e(ji,jj+1) ) > MAX( -bathy(ji,jj), -bathy(ji,jj+1) ) & 
     776                        &                                  + rn_wdmin1 + rn_wdmin2 
     777                 IF(ll_tmp1) THEN 
     778                    zcpy(ji,jj) = 1._wp 
     779                 ELSE IF(ll_tmp2) THEN 
     780                    ! no worries about zsshp2_e(ji,jj+1)-zsshp2_e(ji,jj) = 0, it won't happen here 
     781                    zcpy(ji,jj) = ABS( (zsshp2_e(ji,jj+1) + bathy(ji,jj+1) - zsshp2_e(ji,jj) - bathy(ji,jj)) & 
     782                        &             / (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj)) ) 
     783                 ELSE 
     784                    zcpy(ji,jj)    = 0._wp 
     785                    wdvflt1(ji,jj) = 0._wp 
     786                 END IF 
     787              END DO 
     788            END DO 
     789            CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp ) 
     790         ENDIF 
    674791         ! 
    675792         ! Compute associated depths at U and V points: 
     
    688805               END DO 
    689806            END DO 
     807 
     808            IF( ln_wd ) THEN 
     809              zhust_e(:,:) = MAX(zhust_e (:,:), rn_wdmin1 ) 
     810              zhvst_e(:,:) = MAX(zhvst_e (:,:), rn_wdmin1 ) 
     811            END IF 
     812 
    690813         ENDIF 
    691814         ! 
     
    758881         ! 
    759882         ! Surface pressure trend: 
    760          DO jj = 2, jpjm1 
    761             DO ji = fs_2, fs_jpim1   ! vector opt. 
    762                ! Add surface pressure gradient 
    763                zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj) 
    764                zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj) 
    765                zwx(ji,jj) = zu_spg 
    766                zwy(ji,jj) = zv_spg 
    767             END DO 
    768          END DO 
     883 
     884         IF( ln_wd ) THEN 
     885           DO jj = 2, jpjm1 
     886              DO ji = 2, jpim1  
     887                 ! Add surface pressure gradient 
     888                 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj) 
     889                 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj) 
     890                 zwx(ji,jj) = zu_spg * zcpx(ji,jj)  
     891                 zwy(ji,jj) = zv_spg * zcpy(ji,jj) 
     892              END DO 
     893           END DO 
     894         ELSE 
     895           DO jj = 2, jpjm1 
     896              DO ji = fs_2, fs_jpim1   ! vector opt. 
     897                 ! Add surface pressure gradient 
     898                 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj) 
     899                 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj) 
     900                 zwx(ji,jj) = zu_spg 
     901                 zwy(ji,jj) = zv_spg 
     902              END DO 
     903           END DO 
     904         END IF 
     905 
    769906         ! 
    770907         ! Set next velocities: 
     
    789926            DO jj = 2, jpjm1 
    790927               DO ji = fs_2, fs_jpim1   ! vector opt. 
    791                   zhura = ssumask(ji,jj)/(hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - ssumask(ji,jj)) 
    792                   zhvra = ssvmask(ji,jj)/(hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - ssvmask(ji,jj)) 
     928 
     929                  IF( ln_wd ) THEN 
     930                    zhura = MAX(hu_0(ji,jj) + zsshu_a(ji,jj), rn_wdmin1) 
     931                    zhvra = MAX(hv_0(ji,jj) + zsshv_a(ji,jj), rn_wdmin1) 
     932                  ELSE 
     933                    zhura = hu_0(ji,jj) + zsshu_a(ji,jj) 
     934                    zhvra = hv_0(ji,jj) + zsshv_a(ji,jj) 
     935                  END IF 
     936                  zhura = ssumask(ji,jj)/(zhura + 1._wp - ssumask(ji,jj)) 
     937                  zhvra = ssvmask(ji,jj)/(zhvra + 1._wp - ssvmask(ji,jj)) 
    793938 
    794939                  ua_e(ji,jj) = (                hu_e(ji,jj)  *   un_e(ji,jj)   &  
     
    808953         ! 
    809954         IF( .NOT.ln_linssh ) THEN                     !* Update ocean depth (variable volume case only) 
    810             hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 
    811             hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 
     955            IF( ln_wd ) THEN 
     956              hu_e (:,:) = MAX(hu_0(:,:) + zsshu_a(:,:), rn_wdmin1) 
     957              hv_e (:,:) = MAX(hv_0(:,:) + zsshv_a(:,:), rn_wdmin1) 
     958            ELSE 
     959              hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 
     960              hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 
     961            END IF 
    812962            hur_e(:,:) = ssumask(:,:) / ( hu_e(:,:) + 1._wp - ssumask(:,:) ) 
    813963            hvr_e(:,:) = ssvmask(:,:) / ( hv_e(:,:) + 1._wp - ssvmask(:,:) ) 
     
    9361086      CALL wrk_dealloc( jpi,jpj,   zsshu_a, zsshv_a                                   ) 
    9371087      CALL wrk_dealloc( jpi,jpj,   zhf ) 
     1088      IF( ln_wd ) CALL wrk_dealloc( jpi, jpj, zcpx, zcpy, wduflt1, wdvflt1 ) 
    9381089      ! 
    9391090      IF ( ln_diatmb ) THEN 
  • trunk/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90

    r6140 r6152  
    3939   USE wrk_nemo       ! Memory Allocation 
    4040   USE timing         ! Timing 
     41   USE wet_dry         ! Wetting/Drying flux limting 
    4142 
    4243   IMPLICIT NONE 
     
    104105      !  
    105106      zcoef = 0.5_wp * r1_rau0 
     107 
     108      IF(ln_wd) CALL wad_lmt(sshb, zcoef * (emp_b(:,:) + emp(:,:)), z2dt) 
     109 
    106110      ssha(:,:) = (  sshb(:,:) - z2dt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * ssmask(:,:) 
    107111 
  • trunk/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90

    r6140 r6152  
    9292   USE diatmb          ! Top,middle,bottom output 
    9393   USE dia25h          ! 25h mean output 
     94   USE wet_dry         ! Wetting and drying setting   (wad_init routine) 
    9495 
    9596   IMPLICIT NONE 
     
    402403                            CALL     eos_init   ! Equation of state 
    403404      IF( lk_c1d        )   CALL     c1d_init   ! 1D column configuration 
     405                            CALL     wad_init   ! Wetting and drying options 
    404406                            CALL     dom_cfg    ! Domain configuration 
    405407                            CALL     dom_init   ! Domain 
  • trunk/NEMOGCM/SETTE/sette_rpt.sh

    r6140 r6152  
    2929    f2t=$vdir/$nam/$mach/$dorv/SHORT/tracer.stat 
    3030 
    31     if  [ ! -f $f1s ]; then  
    32       printf "%-20s %s\n" $nam " incomplete test"; 
    33       return;  
    34     fi 
    35     if  [ ! -f $f2s ]; then  
    36       printf "%-20s %s\n" $nam " incomplete test"; 
    37       return;  
    38     fi 
    39 # 
    40 # Check for tracer.stat files 
    41 # 
    42     ttest=-1 
    43     if  [  -f $f1t ]; then (( ttest++ )); fi 
    44     if  [  -f $f2t ]; then (( ttest++ )); fi 
     31    if  [ ! -f $f1s ] &&  [ ! -f $f1t ] ; then  
     32      printf "%-20s %s\n" $nam " incomplete test"; 
     33      return;  
     34    fi 
     35    if  [ ! -f $f2s ] &&  [ ! -f $f2t ] ; then  
     36      printf "%-20s %s\n" $nam " incomplete test"; 
     37      return;  
     38    fi 
     39# 
    4540    done_oce=0 
    4641 
    47     nl=(`wc -l $f2s`) 
    48     tail -${nl[0]} $f1s > f1.tmp$$ 
    49     cmp -s f1.tmp$$ $f2s 
    50     if [ $? == 0 ]; then 
    51       if [ $pass == 0 ]; then printf "%-20s %s\n" $nam  " solver.stat restartability  passed";fi 
    52     else 
    53       printf "%-20s %s\n" $nam  " solver.stat restartability  FAILED" 
    54 # 
    55 # Offer view of differences on the second pass 
    56 # 
    57       if [ $pass == 1 ]; then 
    58         echo "<return> to view solver.stat differences" 
    59         read y 
    60         sdiff f1.tmp$$ $f2s 
    61         echo "<return> to view ocean.output differences" 
    62         read y 
    63         sdiff $f1o $f2o | grep "|" 
    64         done_oce=1 
    65         echo "<return> to continue" 
    66         read y 
     42    if  [  -f $f1s ] && [  -f $f2s ]; then  
     43      nl=(`wc -l $f2s`) 
     44      tail -${nl[0]} $f1s > f1.tmp$$ 
     45      cmp -s f1.tmp$$ $f2s 
     46      if [ $? == 0 ]; then 
     47        if [ $pass == 0 ]; then  
     48          printf "%-20s %s %s\n" $nam  " solver.stat restartability  passed : " $dorv 
     49        fi 
     50      else 
     51        printf "%-20s %s %s\n" $nam  " solver.stat restartability  FAILED : " $dorv  
     52# 
     53# Offer view of differences on the second pass 
     54# 
     55        if [ $pass == 1 ]; then 
     56          echo "<return> to view solver.stat differences" 
     57          read y 
     58          sdiff f1.tmp$$ $f2s 
     59          echo "<return> to view ocean.output differences" 
     60          read y 
     61          sdiff $f1o $f2o | grep "|" 
     62          done_oce=1 
     63          echo "<return> to continue" 
     64          read y 
     65        fi 
    6766      fi 
    6867    fi 
     
    7069# Check tracer.stat files (if they exist) 
    7170# 
    72     if [ $ttest == 1 ]; then 
     71    if  [  -f $f1t ] && [  -f $f2t ]; then 
    7372      nl=(`wc -l $f2t`) 
    7473      tail -${nl[0]} $f1t > f1.tmp$$ 
    7574      cmp -s f1.tmp$$ $f2t 
    7675      if [ $? == 0 ]; then 
    77         if [ $pass == 0 ]; then printf "%-20s %s\n" $nam  " tracer.stat restartability  passed";fi 
    78       else 
    79         printf "%-20s %s\n" $nam  " solver.stat restartability  FAILED" 
     76        if [ $pass == 0 ]; then  
     77          printf "%-20s %s %s\n" $nam  " tracer.stat restartability  passed : " $dorv 
     78        fi 
     79      else 
     80        printf "%-20s %s %s\n" $nam  " tracer.stat restartability  FAILED : " $dorv  
    8081# 
    8182# Offer view of differences on the second pass 
     
    123124    f2t=$vdir/$nam/$mach/$dorv/$rep2/tracer.stat 
    124125 
    125     if  [ ! -f $f1s ]; then  
    126       printf "%-20s %s\n" $nam " incomplete test"; 
    127       return;  
    128     fi 
    129     if  [ ! -f $f2s ]; then  
    130       printf "%-20s %s\n" $nam " incomplete test"; 
    131       return;  
    132     fi 
    133 # 
    134 # Check for tracer.stat files 
    135 # 
    136     ttest=-1 
    137     if  [  -f $f1t ]; then (( ttest++ )); fi 
    138     if  [  -f $f2t ]; then (( ttest++ )); fi 
     126    if  [ ! -f $f1s ] && [ ! -f $f1t ] ; then  
     127      printf "%-20s %s\n" $nam " incomplete test"; 
     128      return;  
     129    fi 
     130    if  [ ! -f $f2s ] && [ ! -f $f2t ] ; then  
     131      printf "%-20s %s\n" $nam " incomplete test"; 
     132      return;  
     133    fi 
     134# 
    139135    done_oce=0 
    140136 
    141     cmp -s $f1s $f2s 
    142     if [ $? == 0 ]; then 
    143       if [ $pass == 0 ]; then printf "%-20s %s\n" $nam  " solver.stat reproducibility passed";fi 
    144     else 
    145       printf "%-20s %s\n" $nam  " solver.stat reproducibility FAILED" 
    146 # 
    147 # Offer view of differences on the second pass 
    148 # 
    149       if [ $pass == 1 ]; then 
    150         echo "<return> to view solver.stat differences" 
    151         read y 
    152         sdiff f1.tmp$$ $f2s 
    153         echo "<return> to view ocean.output differences" 
    154         read y 
    155         sdiff $f1o $f2o | grep "|" 
    156         done_oce=1 
    157         echo "<return> to continue" 
    158         read y 
     137    if  [ -f $f1s ] && [ -f $f2s ] ; then 
     138      cmp -s $f1s $f2s 
     139      if [ $? == 0 ]; then 
     140        if [ $pass == 0 ]; then  
     141          printf "%-20s %s %s\n" $nam  " solver.stat reproducibility passed : " $dorv 
     142        fi 
     143      else 
     144        printf "%-20s %s %s\n" $nam  " solver.stat reproducibility FAILED : " $dorv  
     145# 
     146# Offer view of differences on the second pass 
     147# 
     148        if [ $pass == 1 ]; then 
     149          echo "<return> to view solver.stat differences" 
     150          read y 
     151          sdiff f1.tmp$$ $f2s 
     152          echo "<return> to view ocean.output differences" 
     153          read y 
     154          sdiff $f1o $f2o | grep "|" 
     155          done_oce=1 
     156          echo "<return> to continue" 
     157          read y 
     158        fi 
    159159      fi 
    160160    fi 
     
    162162# Check tracer.stat files (if they exist) 
    163163# 
    164     if [ $ttest == 1 ]; then 
     164    if  [ -f $f1t ] && [ -f $f2t ] ; then 
    165165      cmp -s $f1t $f2t 
    166166      if [ $? == 0 ]; then 
    167         if [ $pass == 0 ]; then printf "%-20s %s\n" $nam  " tracer.stat reproducibility passed";fi 
    168       else 
    169         printf "%-20s %s\n" $nam  " solver.stat reproducibility FAILED" 
     167        if [ $pass == 0 ]; then           printf "%-20s %s %s\n" $nam  " tracer.stat reproducibility passed : " $dorv 
     168        fi 
     169      else 
     170        printf "%-20s %s %s\n" $nam  " tracer.stat reproducibility  FAILED : " $dorv 
    170171# 
    171172# Offer view of differences on the second pass 
Note: See TracChangeset for help on using the changeset viewer.