Changeset 5842


Ignore:
Timestamp:
2015-10-30T16:34:24+01:00 (5 years ago)
Author:
hliu
Message:

Wetting and Drying update based on r:5803

Location:
branches/2015/dev_r5803_NOC_WAD/NEMOGCM
Files:
1 added
9 edited

Legend:

Unmodified
Added
Removed
  • branches/2015/dev_r5803_NOC_WAD/NEMOGCM/CONFIG/AMM12/cpp_AMM12.fcm

    r4245 r5842  
    1  bld::tool::fppkeys  key_bdy key_tide key_dynspg_ts key_ldfslp  key_zdfgls  key_vvl key_diainstant key_mpp_mpi key_iomput 
     1 bld::tool::fppkeys  key_bdy key_tide key_dynspg_ts key_ldfslp  key_zdfgls  key_vvl key_diainstant key_mpp_mpi  
  • branches/2015/dev_r5803_NOC_WAD/NEMOGCM/CONFIG/cfg.txt

    r5656 r5842  
    55C1D_PAPA OPA_SRC 
    66GYRE_BFM OPA_SRC TOP_SRC 
    7 AMM12 OPA_SRC 
    87ORCA2_LIM_PISCES OPA_SRC LIM_SRC_2 NST_SRC TOP_SRC 
    98GYRE OPA_SRC 
     
    1110ORCA2_LIM OPA_SRC LIM_SRC_2 NST_SRC 
    1211ORCA2_OFF_PISCES OPA_SRC OFF_SRC TOP_SRC 
     12AMM12 OPA_SRC 
     13IRS18WD OPA_SRC 
  • branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90

    r5123 r5842  
    1010   !!            3.5  ! 2012     (S. Mocavero, I. Epicoco) Add arrays associated 
    1111   !!                             to the optimization of BDY communications 
     12   !!            3.6.?! 2014     (H. Liu) Add arrays associated with Wetting and Drying 
    1213   !!---------------------------------------------------------------------- 
    1314 
     
    270271   INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  :: nicoa, njcoa       !: ??? 
    271272#endif 
     273 
     274   !!---------------------------------------------------------------------- 
     275   !! critical depths,filters, limiters,and masks for  Wetting and Drying 
     276   !! --------------------------------------------------------------------- 
     277 
     278   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wduflt, wdvflt     !: u- and v- filter 
     279   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wdmask             !: u- and v- limiter  
     280 
     281   LOGICAL,  PUBLIC, SAVE ::   ln_wd       !: key to turn on/off wetting/drying (T: on, F: off) 
     282   REAL(wp), PUBLIC, SAVE ::   rn_wdmin1   !: minimum water depth on dried cells 
     283   REAL(wp), PUBLIC, SAVE ::   rn_wdmin2   !: tolerrance of minimum water depth on dried cells 
     284   REAL(wp), PUBLIC, SAVE ::   rn_wdld     !: land elevation below which wetting/drying will be considered 
     285   INTEGER , PUBLIC, SAVE ::   nn_wdit     !: maximum number of iteration for W/D limiter 
     286 
     287   !LOGICAL,  PUBLIC, SAVE ::   ln_wd     =  .FALSE.  !: turn on wetting/drying (T: on, F: off) 
     288   !REAL(wp), PUBLIC, SAVE ::   rn_wdmin1 = 0.10_wp   !: minimum water depth on dried cells 
     289   !REAL(wp), PUBLIC, SAVE ::   rn_wdmin2 = 0.01_wp   !: tolerrance of minimum water depth on dried cells 
     290   !REAL(wp), PUBLIC, SAVE ::   rn_wdld   = 20.0_wp   !: land elevation below which wetting/drying will be considered 
     291   !INTEGER , PUBLIC, SAVE ::   nn_wdit   =   10      !: maximum number of iteration for W/D limiter 
    272292 
    273293   !!---------------------------------------------------------------------- 
     
    407427      ALLOCATE( npcoa(4,jpk), nicoa(2*(jpi+jpj),4,jpk), njcoa(2*(jpi+jpj),4,jpk), STAT=ierr(12) ) 
    408428#endif 
     429 
     430      IF(ln_wd) & 
     431      ALLOCATE( wduflt(jpi,jpj), wdvflt(jpi,jpj), wdmask(jpi,jpj), STAT=ierr(12) ) 
    409432      ! 
    410433      dom_oce_alloc = MAXVAL(ierr) 
  • branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90

    r5506 r5842  
    147147      fse3t_a(:,:,jpk) = e3t_0(:,:,jpk) 
    148148 
     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 
    149161      ! Reconstruction of all vertical scale factors at now and before time steps 
    150162      ! ============================================================================= 
     
    701713      !                                                                !   =  'U', 'V', 'W, 'F', 'UW' or 'VW' 
    702714      !! * Local declarations 
     715      REAL(wp) ::  zwad                                                ! = 1.0 when ln_wd = .true. 
     716                                                                       ! = 0.0 when ln_wd = .false. 
     717                                                                       !  
    703718      INTEGER ::   ji, jj, jk                                          ! dummy loop indices 
    704719      LOGICAL ::   l_is_orca                                           ! local logical 
     
    706721      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_interpol') 
    707722         ! 
     723      IF(ln_wd) THEN 
     724        zwad = 1.0_wp 
     725      ELSE 
     726        zwad = 0.0_wp 
     727      END IF 
     728 
    708729      l_is_orca = .FALSE. 
    709730      IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) l_is_orca = .TRUE.      ! ORCA R2 configuration - will need to correct some locations 
     
    717738            DO jj = 1, jpjm1 
    718739               DO ji = 1, fs_jpim1   ! vector opt. 
    719                   pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * r1_e12u(ji,jj)                                   & 
     740                  !pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * r1_e12u(ji,jj)                                   & 
     741                  pe3_out(ji,jj,jk) = 0.5_wp * (umask(ji,jj,jk) * (1.0_wp - zwad) + zwad) * r1_e12u(ji,jj)        & 
    720742                     &                       * (   e12t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     & 
    721743                     &                           + e12t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) ) 
     
    735757            DO jj = 1, jpjm1 
    736758               DO ji = 1, fs_jpim1   ! vector opt. 
    737                   pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) * r1_e12v(ji,jj)                                   & 
     759                  !pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) * r1_e12v(ji,jj)                                   & 
     760                  pe3_out(ji,jj,jk) = 0.5_wp * (vmask(ji,jj,jk) * (1.0_wp - zwad) + zwad) * r1_e12v(ji,jj)        & 
    738761                     &                       * (   e12t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     & 
    739762                     &                           + e12t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) ) 
     
    753776            DO jj = 1, jpjm1 
    754777               DO ji = 1, fs_jpim1   ! vector opt. 
    755                   pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) * r1_e12f(ji,jj)               & 
     778                  !pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) * r1_e12f(ji,jj)               & 
     779                  pe3_out(ji,jj,jk) = 0.5_wp * (umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zwad) + zwad)     & 
     780                     &                       * r1_e12f(ji,jj)                                                     & 
    756781                     &                       * (   e12u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     & 
    757782                     &                           + e12u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) ) 
     
    771796         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing 
    772797         DO jk = 2, jpk 
    773             pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * tmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) )   & 
    774                &                            +            0.5_wp * tmask(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3t_0(:,:,jk  ) ) 
     798            !pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * tmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) )   & 
     799            !   &                            +            0.5_wp * tmask(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3t_0(:,:,jk  ) ) 
     800            pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * (tmask(:,:,jk) * (1.0_wp - zwad) + zwad) ) & 
     801               &                            * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) )                         & 
     802               &                            +            0.5_wp * (tmask(:,:,jk) * (1.0_wp - zwad) + zwad)   & 
     803               &                            * ( pe3_in(:,:,jk  ) - e3t_0(:,:,jk  ) ) 
    775804         END DO 
    776805         !               ! -------------------------------------- ! 
     
    781810         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing 
    782811         DO jk = 2, jpk 
    783             pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * umask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) )   & 
    784                &                             +            0.5_wp * umask(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3u_0(:,:,jk  ) ) 
     812            !pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * umask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) )   & 
     813            !   &                             +            0.5_wp * umask(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3u_0(:,:,jk  ) ) 
     814            pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * (umask(:,:,jk) * (1.0_wp - zwad) + zwad) ) & 
     815               &                             * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) )                         & 
     816               &                             +            0.5_wp * (umask(:,:,jk) * (1.0_wp - zwad) + zwad)   & 
     817               &                             * ( pe3_in(:,:,jk  ) - e3u_0(:,:,jk  ) ) 
    785818         END DO 
    786819         !               ! -------------------------------------- ! 
     
    791824         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing 
    792825         DO jk = 2, jpk 
    793             pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * vmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) )   & 
    794                &                             +            0.5_wp * vmask(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3v_0(:,:,jk  ) ) 
     826            !pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * vmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) )   & 
     827            !   &                             +            0.5_wp * vmask(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3v_0(:,:,jk  ) ) 
     828            pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * (vmask(:,:,jk) * (1.0_wp - zwad) + zwad) ) & 
     829               &                             * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) )                         & 
     830               &                             +            0.5_wp * (vmask(:,:,jk) * (1.0_wp - zwad) + zwad)   & 
     831               &                             * ( pe3_in(:,:,jk  ) - e3v_0(:,:,jk  ) ) 
    795832         END DO 
    796833      END SELECT 
     
    817854      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag 
    818855      !! * Local declarations 
    819       INTEGER ::   jk 
     856      INTEGER ::   ji, jj, jk 
    820857      INTEGER ::   id1, id2, id3, id4, id5     ! local integers 
    821858      !!---------------------------------------------------------------------- 
     
    905942            fse3t_n(:,:,:) = e3t_0(:,:,:) 
    906943            sshn(:,:) = 0.0_wp 
     944 
     945            IF(ln_wd) THEN 
     946              DO jj = 1, jpj 
     947                DO ji = 1, jpi 
     948                  !IF(e3t_0(ji,jj,1) < 0._wp) THEN 
     949                  !IF(mbathy(ji,jj) == 2 .AND. e3t_0(ji,jj,1) <= 0.5_wp * rn_wdmin1) THEN 
     950                  IF( e3t_0(ji,jj,1) <= 0.5_wp * rn_wdmin1) THEN 
     951                    fse3t_b(ji,jj,:) = 0.5_wp * rn_wdmin1  
     952                    fse3t_n(ji,jj,:) = 0.5_wp * rn_wdmin1  
     953                    fse3t_a(ji,jj,:) = 0.5_wp * rn_wdmin1  
     954                    sshb(ji,jj) = rn_wdmin1 - bathy(ji,jj) 
     955                    sshn(ji,jj) = rn_wdmin1 - bathy(ji,jj) 
     956                    ssha(ji,jj) = rn_wdmin1 - bathy(ji,jj) 
     957                  ENDIF 
     958                ENDDO 
     959              ENDDO 
     960            END IF 
     961 
    907962            IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN 
    908963               tilde_e3t_b(:,:,:) = 0.0_wp 
  • branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90

    r5656 r5842  
    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 
     
    18101811      REAL(wp) ::   zrmax, ztaper   ! temporary scalars 
    18111812      REAL(wp) ::   zrfact 
     1813      REAL(wp) ::   zsmth 
    18121814      ! 
    18131815      REAL(wp), POINTER, DIMENSION(:,:  ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2 
     
    18651867      bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) ) 
    18661868 
     1869      IF(.NOT.ln_wd) THEN                       
    18671870      DO jj = 1, jpj 
    18681871         DO ji = 1, jpi 
     
    18701873         END DO 
    18711874      END DO 
     1875      END IF 
    18721876      !                                        ! ============================= 
    18731877      !                                        ! Define the envelop bathymetry   (hbatt) 
     
    18761880      zenv(:,:) = bathy(:,:) 
    18771881      ! 
     1882      IF(.NOT.ln_wd) THEN     
    18781883      ! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing 
    18791884      DO jj = 1, jpj 
     
    18911896         END DO 
    18921897      END DO 
     1898      END IF 
     1899 
    18931900      ! apply lateral boundary condition   CAUTION: keep the value when the lbc field is zero 
    18941901      CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' ) 
     
    19301937         zri(:,:) = 0._wp 
    19311938         zrj(:,:) = 0._wp 
     1939          
     1940         !IF(ln_wd) THEN                     !extend the smoothed region to cover the W/D zones 
     1941         !  zsmth = -rn_wdld     
     1942         !ELSE 
     1943           zsmth = 0._wp                     ! The original form (only smooth ocean points) 
     1944         !ENDIF 
     1945 
    19321946         DO jj = 1, nlcj 
    19331947            DO ji = 1, nlci 
    19341948               iip1 = MIN( ji+1, nlci )      ! force zri = 0 on last line (ji=ncli+1 to jpi) 
    19351949               ijp1 = MIN( jj+1, nlcj )      ! force zrj = 0 on last raw  (jj=nclj+1 to jpj) 
    1936                IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(iip1,jj) > 0._wp)) THEN 
     1950               IF( (zenv(ji,jj) > zsmth) .AND. (zenv(iip1,jj) > zsmth)) THEN 
    19371951                  zri(ji,jj) = ( zenv(iip1,jj  ) - zenv(ji,jj) ) / ( zenv(iip1,jj  ) + zenv(ji,jj) ) 
    19381952               END IF 
    1939                IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(ji,ijp1) > 0._wp)) THEN 
     1953               IF( (zenv(ji,jj) > zsmth) .AND. (zenv(ji,ijp1) > zsmth)) THEN 
    19401954                  zrj(ji,jj) = ( zenv(ji  ,ijp1) - zenv(ji,jj) ) / ( zenv(ji  ,ijp1) + zenv(ji,jj) ) 
    19411955               END IF 
     
    19601974      END DO                                                !     End loop     ! 
    19611975      !                                                     ! ================ ! 
    1962       DO jj = 1, jpj 
    1963          DO ji = 1, jpi 
    1964             zenv(ji,jj) = MAX( zenv(ji,jj), rn_sbot_min ) ! set all points to avoid undefined scale value warnings 
    1965          END DO 
    1966       END DO 
     1976      IF(ln_wd) THEN 
     1977        DO jj = 1, jpj 
     1978           DO ji = 1, jpi 
     1979              !zenv(ji,jj) = MAX( zenv(ji,jj), -rn_wdld )    ! filt out land bathy data  
     1980           END DO 
     1981        END DO 
     1982      ELSE 
     1983        DO jj = 1, jpj 
     1984           DO ji = 1, jpi 
     1985              zenv(ji,jj) = MAX( zenv(ji,jj), rn_sbot_min ) ! set all points to avoid undefined scale value warnings 
     1986           END DO 
     1987        END DO 
     1988      ENDIF 
    19671989      ! 
    19681990      ! Envelope bathymetry saved in hbatt 
     
    19942016      IF(lwp) THEN 
    19952017         WRITE(numout,*) 
    1996          WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min 
     2018         IF(.NOT.ln_wd) THEN 
     2019           WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min 
     2020         ELSE 
     2021           WRITE(numout,*) ' zgr_sco: minimum positive depth of the envelop topography set to : ', rn_sbot_min 
     2022           WRITE(numout,*) ' zgr_sco: minimum negative depth of the envelop topography set to : ', -rn_wdld 
     2023         ENDIF 
    19972024      ENDIF 
    19982025      hbatu(:,:) = rn_sbot_min 
     
    20072034        END DO 
    20082035      END DO 
     2036 
     2037      IF(ln_wd) THEN               !avoid the zero depth on T- (U-,V-,F-) points 
     2038        DO jj = 1, jpj 
     2039          DO ji = 1, jpi 
     2040            IF(ABS(hbatt(ji,jj)) < rn_wdmin1) & 
     2041              & hbatt(ji,jj) = SIGN(1._wp, hbatt(ji,jj)) * rn_wdmin1 
     2042            IF(ABS(hbatu(ji,jj)) < rn_wdmin1) & 
     2043              & hbatu(ji,jj) = SIGN(1._wp, hbatu(ji,jj)) * rn_wdmin1 
     2044            IF(ABS(hbatv(ji,jj)) < rn_wdmin1) & 
     2045              & hbatv(ji,jj) = SIGN(1._wp, hbatv(ji,jj)) * rn_wdmin1 
     2046            IF(ABS(hbatf(ji,jj)) < rn_wdmin1) & 
     2047              & hbatf(ji,jj) = SIGN(1._wp, hbatf(ji,jj)) * rn_wdmin1 
     2048          END DO 
     2049        END DO 
     2050      END IF 
    20092051      !  
    20102052      ! Apply lateral boundary condition 
     
    20142056         DO ji = 1, jpi 
    20152057            IF( hbatu(ji,jj) == 0._wp ) THEN 
     2058               !No worries about the following line when ln_wd == .true. 
    20162059               IF( zhbat(ji,jj) == 0._wp )   hbatu(ji,jj) = rn_sbot_min 
    20172060               IF( zhbat(ji,jj) /= 0._wp )   hbatu(ji,jj) = zhbat(ji,jj) 
     
    20392082 
    20402083!!bug:  key_helsinki a verifer 
    2041       hift(:,:) = MIN( hift(:,:), hbatt(:,:) ) 
    2042       hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) ) 
    2043       hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) ) 
    2044       hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) ) 
     2084      IF(.NOT.ln_wd) THEN 
     2085       hift(:,:) = MIN( hift(:,:), hbatt(:,:) ) 
     2086       hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) ) 
     2087       hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) ) 
     2088       hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) ) 
     2089      END IF 
    20452090 
    20462091      IF( nprint == 1 .AND. lwp )   THEN 
     
    20872132      fsde3w(:,:,:) = gdep3w_0(:,:,:) 
    20882133      ! 
    2089       where (e3t_0   (:,:,:).eq.0.0)  e3t_0(:,:,:) = 1.0 
    2090       where (e3u_0   (:,:,:).eq.0.0)  e3u_0(:,:,:) = 1.0 
    2091       where (e3v_0   (:,:,:).eq.0.0)  e3v_0(:,:,:) = 1.0 
    2092       where (e3f_0   (:,:,:).eq.0.0)  e3f_0(:,:,:) = 1.0 
    2093       where (e3w_0   (:,:,:).eq.0.0)  e3w_0(:,:,:) = 1.0 
    2094       where (e3uw_0  (:,:,:).eq.0.0)  e3uw_0(:,:,:) = 1.0 
    2095       where (e3vw_0  (:,:,:).eq.0.0)  e3vw_0(:,:,:) = 1.0 
     2134      IF(.NOT.ln_wd) THEN 
     2135        where (e3t_0   (:,:,:).eq.0.0)  e3t_0(:,:,:) = 1.0 
     2136        where (e3u_0   (:,:,:).eq.0.0)  e3u_0(:,:,:) = 1.0 
     2137        where (e3v_0   (:,:,:).eq.0.0)  e3v_0(:,:,:) = 1.0 
     2138        where (e3f_0   (:,:,:).eq.0.0)  e3f_0(:,:,:) = 1.0 
     2139        where (e3w_0   (:,:,:).eq.0.0)  e3w_0(:,:,:) = 1.0 
     2140        where (e3uw_0  (:,:,:).eq.0.0)  e3uw_0(:,:,:) = 1.0 
     2141        where (e3vw_0  (:,:,:).eq.0.0)  e3vw_0(:,:,:) = 1.0 
     2142      END IF 
    20962143 
    20972144#if defined key_agrif 
     
    21352182               IF( scobot(ji,jj) >= fsdept(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk ) 
    21362183            END DO 
    2137             IF( scobot(ji,jj) == 0._wp               )   mbathy(ji,jj) = 0 
    2138          END DO 
    2139       END DO 
     2184 
     2185            IF(ln_wd) THEN 
     2186              IF( scobot(ji,jj) <= -(rn_wdld - rn_wdmin2)) THEN 
     2187                mbathy(ji,jj) = 0 
     2188              ELSEIF(scobot(ji,jj) <= rn_wdmin1) THEN 
     2189                mbathy(ji,jj) = 2 
     2190              ENDIF 
     2191            ELSE 
     2192              IF( scobot(ji,jj) == 0._wp )   mbathy(ji,jj) = 0 
     2193            ENDIF 
     2194         END DO 
     2195      END DO 
     2196 
    21402197      IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ),   & 
    21412198         &                                                       ' MAX ', MAXVAL( mbathy(:,:) ) 
     
    22572314      INTEGER  ::   ji, jj, jk           ! dummy loop argument 
    22582315      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars 
     2316      REAL(wp) ::   ztmpu,  ztmpv,  ztmpf 
     2317      REAL(wp) ::   ztmpu1, ztmpv1, ztmpf1 
    22592318      ! 
    22602319      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 
     
    23102369      DO ji = 1, jpim1 
    23112370         DO jj = 1, jpjm1 
     2371            ! extended for Wetting/Drying case 
     2372            ztmpu = hbatt(ji,jj)+hbatt(ji+1,jj) 
     2373            ztmpv = hbatt(ji,jj)+hbatt(ji,jj+1) 
     2374            ztmpf = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) 
     2375            ztmpu1 = hbatt(ji,jj)*hbatt(ji+1,jj) 
     2376            ztmpv1 = hbatt(ji,jj)*hbatt(ji,jj+1) 
     2377            ztmpf1 = MIN(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) * & 
     2378                   & MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) 
    23122379            DO jk = 1, jpk 
    2313                z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) )   & 
    2314                   &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
    2315                z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) )   & 
    2316                   &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
    2317                z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk)     & 
    2318                   &                + hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) )   & 
    2319                   &              / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 
    2320                z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) )   & 
    2321                   &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
    2322                z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) )   & 
    2323                   &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
     2380               IF((ztmpu1 < 0._wp.OR.ABS(ztmpu) < rn_wdmin1).AND.ln_wd) THEN 
     2381                 z_esigtu3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) ) 
     2382                 z_esigwu3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji+1,jj,jk) ) 
     2383               ELSE 
     2384                 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) & 
     2385                        &              / ztmpu 
     2386                 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) & 
     2387                        &              / ztmpu 
     2388               END IF 
     2389 
     2390               IF((ztmpv1 < 0._wp.OR.ABS(ztmpv) < rn_wdmin1).AND.ln_wd) THEN 
     2391                 z_esigtv3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji,jj+1,jk) ) 
     2392                 z_esigwv3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji,jj+1,jk) ) 
     2393               ELSE 
     2394                 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) & 
     2395                        &              / ztmpv 
     2396                 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) & 
     2397                        &              / ztmpv 
     2398               END IF 
     2399 
     2400               IF((ztmpf1 < 0._wp.OR.ABS(ztmpf) < rn_wdmin1).AND.ln_wd) THEN 
     2401                 z_esigtf3(ji,jj,jk) = 0.25_wp * ( z_esigt3(ji,jj,jk)   + z_esigt3(ji+1,jj,jk)  & 
     2402                    &                            + z_esigt3(ji,jj+1,jk) + z_esigt3(ji+1,jj+1,jk) ) 
     2403               ELSE 
     2404                 z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk)  & 
     2405                    &                + hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) & 
     2406                    &              / ztmpf 
     2407               END IF 
     2408 
     2409               ! 
     2410               !z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) )   & 
     2411               !   &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
     2412               !z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) )   & 
     2413               !   &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
     2414               !z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk)     & 
     2415               !   &                + hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) )   & 
     2416               !   &              / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 
     2417               !z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) )   & 
     2418               !   &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
     2419               !z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) )   & 
     2420               !   &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
    23242421               ! 
    23252422               e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 
     
    23612458      REAL(wp) ::   zsmth               ! smoothing around critical depth 
    23622459      REAL(wp) ::   zzs, zzb           ! Surface and bottom cell thickness in sigma space 
     2460      REAL(wp) ::   ztmpu, ztmpv, ztmpf 
     2461      REAL(wp) ::   ztmpu1, ztmpv1, ztmpf1 
    23632462      ! 
    23642463      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 
     
    24392538        DO jj=1,jpj-1 
    24402539 
    2441           DO jk = 1, jpk 
    2442                 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) / & 
    2443                                     ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
    2444                 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) / & 
    2445                                     ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
    2446                 z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) +  & 
    2447                                       hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / & 
    2448                                     ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 
    2449                 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) / & 
    2450                                     ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
    2451                 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) / & 
    2452                                     ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
     2540           ! extend to suit for Wetting/Drying case 
     2541           ztmpu = hbatt(ji,jj)+hbatt(ji+1,jj) 
     2542           ztmpv = hbatt(ji,jj)+hbatt(ji,jj+1) 
     2543           ztmpf = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) 
     2544           ztmpu1 = hbatt(ji,jj)*hbatt(ji+1,jj) 
     2545           ztmpv1 = hbatt(ji,jj)*hbatt(ji,jj+1) 
     2546           ztmpf1 = MIN(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) * & 
     2547                  & MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) 
     2548           DO jk = 1, jpk 
     2549              IF((ztmpu1 < 0._wp.OR.ABS(ztmpu) < rn_wdmin1).AND.ln_wd) THEN 
     2550                z_esigtu3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) ) 
     2551                z_esigwu3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji+1,jj,jk) ) 
     2552              ELSE 
     2553                z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) & 
     2554                       &              / ztmpu 
     2555                z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) & 
     2556                       &              / ztmpu 
     2557              END IF 
     2558 
     2559              IF((ztmpv1 < 0._wp.OR.ABS(ztmpv) < rn_wdmin1).AND.ln_wd) THEN 
     2560                z_esigtv3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji,jj+1,jk) ) 
     2561                z_esigwv3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji,jj+1,jk) ) 
     2562              ELSE 
     2563                z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) & 
     2564                       &              / ztmpv 
     2565                z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) & 
     2566                       &              / ztmpv 
     2567              END IF 
     2568 
     2569              IF((ztmpf1 < 0._wp.OR.ABS(ztmpf) < rn_wdmin1).AND.ln_wd) THEN 
     2570                z_esigtf3(ji,jj,jk) = 0.25_wp * ( z_esigt3(ji,jj,jk)   + z_esigt3(ji+1,jj,jk)  & 
     2571                   &                            + z_esigt3(ji,jj+1,jk) + z_esigt3(ji+1,jj+1,jk) ) 
     2572              ELSE 
     2573                z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk)  & 
     2574                   &                + hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) & 
     2575                   &              / ztmpf 
     2576              END IF 
     2577 
     2578             !z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) / & 
     2579             !                    ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
     2580             !z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) / & 
     2581             !                    ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
     2582             !z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) +  & 
     2583             !                      hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / & 
     2584             !                    ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 
     2585             !z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) / & 
     2586             !                    ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
     2587             !z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) / & 
     2588             !                    ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
    24532589 
    24542590             e3t_0(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk) 
  • branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90

    r5224 r5842  
    1717   !!                 !           (A. Coward) suppression of hel, wdj and rot options 
    1818   !!            3.6  !  2014-11  (P. Mathiot) hpg_isf: original code for ice shelf cavity 
     19   !!            3.6? !  2015-11  (H. Liu) add Wetting/Drying pressure filter 
    1920   !!---------------------------------------------------------------------- 
    2021 
     
    378379      INTEGER, INTENT(in) ::   kt    ! ocean time-step index 
    379380      !! 
    380       INTEGER  ::   ji, jj, jk                 ! dummy loop indices 
    381       REAL(wp) ::   zcoef0, zuap, zvap, znad   ! temporary scalars 
     381      INTEGER  ::   ji, jj, jk, jii, jjj                 ! dummy loop indices 
     382      REAL(wp) ::   zcoef0, zuap, zvap, znad, ztmp       ! temporary scalars 
     383      LOGICAL  ::   ll_tmp1, ll_tmp2, ll_tmp3            ! local logical variables 
    382384      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zhpi, zhpj 
     385      REAL(wp), POINTER, DIMENSION(:,:)   ::  zcpx, zcpy !W/D pressure filter 
    383386      !!---------------------------------------------------------------------- 
    384387      ! 
    385388      CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj ) 
     389      IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 
    386390      ! 
    387391      IF( kt == nit000 ) THEN 
     
    397401      ELSE                 ;     znad = 0._wp         ! Fixed volume 
    398402      ENDIF 
     403 
     404      IF(ln_wd) THEN 
     405        DO jj = 2, jpjm1 
     406           DO ji = 2, jpim1  
     407             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj))  
     408             ll_tmp2 = MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) > rn_wdmin1 + rn_wdmin2 
     409             ll_tmp3 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) + & 
     410                                                       & rn_wdmin1 + rn_wdmin2 
     411 
     412             IF(ll_tmp1.AND.ll_tmp2) THEN 
     413               zcpx(ji,jj) = 1.0_wp 
     414               wduflt(ji,jj) = 1.0_wp 
     415             ELSE IF(ll_tmp3) THEN 
     416               ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 
     417               zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) / & 
     418                           &     (sshn(ji+1,jj) - sshn(ji,jj))) 
     419               wduflt(ji,jj) = 1.0_wp 
     420             ELSE 
     421               zcpx(ji,jj) = 0._wp 
     422               wduflt(ji,jj) = 0.0_wp 
     423             END IF 
     424       
     425             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1))  
     426             ll_tmp2 = MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) > rn_wdmin1 + rn_wdmin2 
     427             ll_tmp3 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) + & 
     428                                                       & rn_wdmin1 + rn_wdmin2 
     429 
     430             IF(ll_tmp1.AND.ll_tmp2) THEN 
     431               zcpy(ji,jj) = 1.0_wp 
     432               wdvflt(ji,jj) = 1.0_wp 
     433             ELSE IF(ll_tmp3) THEN 
     434               ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 
     435               zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) / & 
     436                           &     (sshn(ji,jj+1) - sshn(ji,jj))) 
     437               wdvflt(ji,jj) = 1.0_wp 
     438             ELSE 
     439               zcpy(ji,jj) = 0._wp 
     440               wdvflt(ji,jj) = 0.0_wp 
     441             END IF 
     442           END DO 
     443        END DO 
     444        CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp ) 
     445      ENDIF 
     446 
     447 
     448      !jii=jidbg-nimpp+1;jjj=jjdbg-njmpp+1 
    399449 
    400450      ! Surface value 
     
    411461            zvap = -zcoef0 * ( rhd   (ji,jj+1,1) + rhd   (ji,jj,1) + 2._wp * znad )   & 
    412462               &           * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj) 
     463 
     464 
     465            IF(ln_wd) THEN 
     466 
     467              zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 
     468              zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj)  
     469              zuap = zuap * zcpx(ji,jj) 
     470              zvap = zvap * zcpy(ji,jj) 
     471            ENDIF 
     472 
    413473            ! add to the general momentum trend 
    414474            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap 
     
    433493               zvap = -zcoef0 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   & 
    434494                  &           * ( fsde3w(ji  ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) 
     495 
     496               IF(ln_wd) THEN 
     497                 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 
     498                 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj)  
     499                 zuap = zuap * zcpx(ji,jj) 
     500                 zvap = zvap * zcpy(ji,jj) 
     501               ENDIF 
     502 
    435503               ! add to the general momentum trend 
    436504               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap 
     
    441509      ! 
    442510      CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj ) 
     511      IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 
    443512      ! 
    444513   END SUBROUTINE hpg_sco 
     
    719788      REAL(wp) ::   z1_10, cffu, cffx   !    "         " 
    720789      REAL(wp) ::   z1_12, cffv, cffy   !    "         " 
     790      LOGICAL  ::   ll_tmp1, ll_tmp2    ! local logical variables 
    721791      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zhpi, zhpj 
    722792      REAL(wp), POINTER, DIMENSION(:,:,:) ::  dzx, dzy, dzz, dzu, dzv, dzw 
    723793      REAL(wp), POINTER, DIMENSION(:,:,:) ::  drhox, drhoy, drhoz, drhou, drhov, drhow 
    724794      REAL(wp), POINTER, DIMENSION(:,:,:) ::  rho_i, rho_j, rho_k 
     795      REAL(wp), POINTER, DIMENSION(:,:)   ::  zcpx, zcpy    !W/D pressure filter 
    725796      !!---------------------------------------------------------------------- 
    726797      ! 
     
    728799      CALL wrk_alloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow ) 
    729800      CALL wrk_alloc( jpi, jpj, jpk, rho_i, rho_j, rho_k,  zhpi,  zhpj        ) 
    730       ! 
     801      IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 
     802      ! 
     803      !!---------------------------------------------------------------------- 
     804      ! 
     805      ! 
     806      IF(ln_wd) THEN 
     807        DO jj = 2, jpjm1 
     808           DO ji = 2, jpim1  
     809             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) & 
     810                     & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) & 
     811                     &  > rn_wdmin1 + rn_wdmin2 
     812             ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) +& 
     813                                                       & rn_wdmin1 + rn_wdmin2 
     814 
     815             IF(ll_tmp1) THEN 
     816               zcpx(ji,jj) = 1.0_wp 
     817             ELSE IF(ll_tmp2) THEN 
     818               ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 
     819               zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /& 
     820                           &     (sshn(ji+1,jj) - sshn(ji,jj))) 
     821             ELSE 
     822               zcpx(ji,jj) = 0._wp 
     823             END IF 
     824       
     825             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) & 
     826                     & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) & 
     827                     &  > rn_wdmin1 + rn_wdmin2 
     828             ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) +& 
     829                                                       & rn_wdmin1 + rn_wdmin2 
     830 
     831             IF(ll_tmp1) THEN 
     832               zcpy(ji,jj) = 1.0_wp 
     833             ELSE IF(ll_tmp2) THEN 
     834               ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 
     835               zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) /& 
     836                           &     (sshn(ji,jj+1) - sshn(ji,jj))) 
     837             ELSE 
     838               zcpy(ji,jj) = 0._wp 
     839             END IF 
     840           END DO 
     841        END DO 
     842        CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp ) 
     843      ENDIF 
     844 
    731845 
    732846      IF( kt == nit000 ) THEN 
     
    8991013            zhpi(ji,jj,1) = ( rho_k(ji+1,jj  ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) / e1u(ji,jj) 
    9001014            zhpj(ji,jj,1) = ( rho_k(ji  ,jj+1,1) - rho_k(ji,jj,1) - rho_j(ji,jj,1) ) / e2v(ji,jj) 
     1015            IF(ln_wd) THEN 
     1016              zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 
     1017              zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj)  
     1018            ENDIF 
    9011019            ! add to the general momentum trend 
    9021020            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) 
     
    9181036                  &           + (  ( rho_k(ji,jj+1,jk) - rho_k(ji,jj,jk  ) )    & 
    9191037                  &               -( rho_j(ji,jj  ,jk) - rho_j(ji,jj,jk-1) )  ) / e2v(ji,jj) 
     1038               IF(ln_wd) THEN 
     1039                 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 
     1040                 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj)  
     1041               ENDIF 
    9201042               ! add to the general momentum trend 
    9211043               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) 
     
    9281050      CALL wrk_dealloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow ) 
    9291051      CALL wrk_dealloc( jpi, jpj, jpk, rho_i, rho_j, rho_k,  zhpi,  zhpj        ) 
     1052      IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 
    9301053      ! 
    9311054   END SUBROUTINE hpg_djc 
     
    9511074      !! The local variables for the correction term 
    9521075      INTEGER  :: jk1, jis, jid, jjs, jjd 
     1076      LOGICAL  :: ll_tmp1, ll_tmp2                  ! local logical variables 
    9531077      REAL(wp) :: zuijk, zvijk, zpwes, zpwed, zpnss, zpnsd, zdeps 
    9541078      REAL(wp) :: zrhdt1 
     
    9571081      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp 
    9581082      REAL(wp), POINTER, DIMENSION(:,:)   ::   zsshu_n, zsshv_n 
     1083      REAL(wp), POINTER, DIMENSION(:,:)   ::  zcpx, zcpy    !W/D pressure filter 
    9591084      !!---------------------------------------------------------------------- 
    9601085      ! 
     
    9621087      CALL wrk_alloc( jpi,jpj,jpk, zdept, zrhh ) 
    9631088      CALL wrk_alloc( jpi,jpj, zsshu_n, zsshv_n ) 
     1089      IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 
    9641090      ! 
    9651091      IF( kt == nit000 ) THEN 
     
    9741100      znad = 0.0_wp 
    9751101      IF( lk_vvl ) znad = 1._wp 
     1102 
     1103      IF(ln_wd) THEN 
     1104        DO jj = 2, jpjm1 
     1105           DO ji = 2, jpim1  
     1106             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) & 
     1107                     & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) & 
     1108                     &  > rn_wdmin1 + rn_wdmin2 
     1109             ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) +& 
     1110                                                       & rn_wdmin1 + rn_wdmin2 
     1111 
     1112             IF(ll_tmp1) THEN 
     1113               zcpx(ji,jj) = 1.0_wp 
     1114             ELSE IF(ll_tmp2) THEN 
     1115               ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 
     1116               zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /& 
     1117                           &     (sshn(ji+1,jj) - sshn(ji,jj))) 
     1118             ELSE 
     1119               zcpx(ji,jj) = 0._wp 
     1120             END IF 
     1121       
     1122             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) & 
     1123                     & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) & 
     1124                     &  > rn_wdmin1 + rn_wdmin2 
     1125             ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) +& 
     1126                                                       & rn_wdmin1 + rn_wdmin2 
     1127 
     1128             IF(ll_tmp1.OR.ll_tmp2) THEN 
     1129               zcpy(ji,jj) = 1.0_wp 
     1130             ELSE IF(ll_tmp2) THEN 
     1131               ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 
     1132               zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) /& 
     1133                           &     (sshn(ji,jj+1) - sshn(ji,jj))) 
     1134             ELSE 
     1135               zcpy(ji,jj) = 0._wp 
     1136             END IF 
     1137           END DO 
     1138        END DO 
     1139        CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp ) 
     1140      ENDIF 
    9761141 
    9771142      ! Clean 3-D work arrays 
     
    10521217        END DO 
    10531218      END DO 
     1219 
     1220      CALL lbc_lnk (zsshu_n, 'U', 1.) 
     1221      CALL lbc_lnk (zsshv_n, 'V', 1.) 
    10541222 
    10551223      DO jj = 2, jpjm1 
     
    11501318               ENDIF 
    11511319 
    1152                ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * & 
    1153                &           umask(ji,jj,jk) * tmask(ji,jj,jk) * tmask(ji+1,jj,jk) 
     1320               IF(ln_wd) THEN 
     1321                  zdpdx1 = zdpdx1 * zcpx(ji,jj) 
     1322                  zdpdx2 = zdpdx2 * zcpx(ji,jj) 
     1323                ENDIF 
     1324                ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk)  
    11541325            ENDIF 
    11551326 
     
    12071378               ENDIF 
    12081379 
    1209                va(ji,jj,jk) = va(ji,jj,jk) + (zdpdy1 + zdpdy2)*& 
    1210                &              vmask(ji,jj,jk)*tmask(ji,jj,jk)*tmask(ji,jj+1,jk) 
     1380               IF(ln_wd) THEN 
     1381                  zdpdy1 = zdpdy1 * zcpy(ji,jj) 
     1382                  zdpdy2 = zdpdy2 * zcpy(ji,jj) 
     1383                ENDIF 
     1384 
     1385               va(ji,jj,jk) = va(ji,jj,jk) + (zdpdy1 + zdpdy2) * vmask(ji,jj,jk) 
    12111386            ENDIF 
    12121387 
     
    12191394      CALL wrk_dealloc( jpi,jpj,jpk, zdept, zrhh ) 
    12201395      CALL wrk_dealloc( jpi,jpj, zsshu_n, zsshv_n ) 
     1396      IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 
    12211397      ! 
    12221398   END SUBROUTINE hpg_prj 
  • branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r5656 r5842  
    4141   USE timing          ! Timing     
    4242   USE sbcapr          ! surface boundary condition: atmospheric pressure 
     43   USE wadlmt          ! wetting/drying flux limter 
    4344   USE dynadv, ONLY: ln_dynadv_vec 
    4445#if defined key_agrif 
     
    142143      LOGICAL  ::   ll_fw_start        ! if true, forward integration  
    143144      LOGICAL  ::   ll_init             ! if true, special startup of 2d equations 
     145      LOGICAL  ::   ll_tmp1, ll_tmp2            ! local logical variables used in W/D 
    144146      INTEGER  ::   ji, jj, jk, jn        ! dummy loop indices 
    145147      INTEGER  ::   ikbu, ikbv, noffset      ! local integers 
     
    157159      REAL(wp), POINTER, DIMENSION(:,:) :: zsshu_a, zsshv_a 
    158160      REAL(wp), POINTER, DIMENSION(:,:) :: zhf 
     161      REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy                 ! Wetting/Dying gravity filter coef. 
     162      REAL(wp), POINTER, DIMENSION(:,:) :: wduflt1, wdvflt1           ! Wetting/Dying velocity filter coef. 
    159163      !!---------------------------------------------------------------------- 
     164 
     165      ! tempoary debugging integer 
     166      INTEGER :: jidbg, jjdbg 
     167      jidbg = 163; jjdbg = 179 
    160168      ! 
    161169      IF( nn_timing == 1 )  CALL timing_start('dyn_spg_ts') 
     
    168176      CALL wrk_alloc( jpi, jpj, zsshu_a, zsshv_a                                   ) 
    169177      CALL wrk_alloc( jpi, jpj, zhf ) 
     178      IF(ln_wd) CALL wrk_alloc( jpi, jpj, zcpx, zcpy, wduflt1, wdvflt1 ) 
    170179      ! 
    171180      !                                         !* Local constant initialization 
     
    379388      !                                   ! ---------------------------------------------------- 
    380389      IF( lk_vvl ) THEN                         ! Variable volume : remove surface pressure gradient 
    381          DO jj = 2, jpjm1  
    382             DO ji = fs_2, fs_jpim1   ! vector opt. 
    383                zu_trd(ji,jj) = zu_trd(ji,jj) - grav * (  sshn(ji+1,jj  ) - sshn(ji  ,jj  )  ) / e1u(ji,jj) 
    384                zv_trd(ji,jj) = zv_trd(ji,jj) - grav * (  sshn(ji  ,jj+1) - sshn(ji  ,jj  )  ) / e2v(ji,jj) 
    385             END DO 
    386          END DO 
    387       ENDIF 
    388  
     390        IF(ln_wd) THEN                    ! Calculating and applying W/D gravity filters 
     391          DO jj = 2, jpjm1 
     392             DO ji = 2, jpim1 
     393                ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) & 
     394                        & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) & 
     395                        &  > rn_wdmin1 + rn_wdmin2 
     396                ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) + & 
     397                                                          & rn_wdmin1 + rn_wdmin2 
     398                IF(ll_tmp1) THEN 
     399                  zcpx(ji,jj) = 1.0_wp 
     400                  wduflt1(ji,jj) = 1.0_wp 
     401                ELSE IF(ll_tmp2) THEN 
     402                   ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 
     403                  zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /& 
     404                              &     (sshn(ji+1,jj) - sshn(ji,jj))) 
     405                  wduflt1(ji,jj) = 1.0_wp 
     406                ELSE 
     407                  zcpx(ji,jj) = 0._wp 
     408                  wduflt1(ji,jj) = 0.0_wp 
     409                END IF 
     410 
     411                !for w/d debugging, delete it when finished 
     412                !   zcpx(ji,jj) = 0._wp 
     413                ! end debugging part 
     414                    
     415                ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) & 
     416                        & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) & 
     417                        &  > rn_wdmin1 + rn_wdmin2 
     418                ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) + & 
     419                                                          & rn_wdmin1 + rn_wdmin2 
     420                IF(ll_tmp1) THEN 
     421                   zcpy(ji,jj) = 1.0_wp 
     422                   wdvflt1(ji,jj) = 1.0_wp 
     423                ELSE IF(ll_tmp2) THEN 
     424                   ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 
     425                  zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) /& 
     426                              &     (sshn(ji,jj+1) - sshn(ji,jj))) 
     427                   wdvflt1(ji,jj) = 1.0_wp 
     428                ELSE 
     429                  zcpy(ji,jj) = 0._wp 
     430                  wdvflt1(ji,jj) = 0.0_wp 
     431                END IF 
     432 
     433                !for w/d debugging, delete it when finished 
     434                !   zcpy(ji,jj) = 0._wp 
     435                ! end debugging part 
     436             END DO 
     437           END DO 
     438           CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp ) 
     439        ENDIF 
     440 
     441        IF(ln_wd) THEN                    ! Calculating and applying W/D gravity filters 
     442           DO jj = 2, jpjm1 
     443              DO ji = 2, jpim1 
     444                 zu_trd(ji,jj) = (zu_trd(ji,jj) - grav * ( sshn(ji+1,jj  ) - sshn(ji  ,jj ) ) / & 
     445                               &                 e1u(ji,jj)) * zcpx(ji,jj) 
     446                 zv_trd(ji,jj) = (zv_trd(ji,jj) - grav * ( sshn(ji  ,jj+1) - sshn(ji  ,jj ) ) / & 
     447                               &                 e2v(ji,jj)) * zcpy(ji,jj) 
     448                 IF(ji+nimpp-1 == jidbg .and. jj+njmpp-1 == jjdbg) THEN 
     449                   write(333,*) 'in spt_ts:  zutrdx, zutrdy, hpgx, hpgy, zcpx, zcpy' 
     450                   write(333,*) zu_trd(ji,jj), zv_trd(ji,jj), & 
     451                             &  -grav*(sshn(ji+1,jj)-sshn(ji,jj))/e1u(ji,jj)*zcpx(ji,jj), & 
     452                             &  -grav*(sshn(ji,jj+1)-sshn(ji,jj))/e2v(ji,jj)*zcpy(ji,jj), & 
     453                             &  zcpx(ji,jj), zcpy(ji,jj)  
     454                   write(334,*) 'in spt_ts:  sshn_ij, sshn_ip1j, sshn_ijp1, bathy_ij, bathy_ip1j, bathy_ijp1' 
     455                   write(334,*) sshn(ji,jj), sshn(ji+1,jj), sshn(ji,jj+1), bathy(ji,jj), bathy(ji+1,jj), bathy(ji,jj+1)  
     456                 END IF 
     457              END DO 
     458           END DO 
     459         ELSE 
     460           DO jj = 2, jpjm1 
     461              DO ji = fs_2, fs_jpim1   ! vector opt. 
     462                 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * (  sshn(ji+1,jj  ) - sshn(ji  ,jj  )  ) / e1u(ji,jj) 
     463                 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * (  sshn(ji  ,jj+1) - sshn(ji  ,jj  )  ) / e2v(ji,jj)  
     464              END DO 
     465           END DO 
     466        END IF 
     467 
     468      ENDIF 
     469 
     470      IF(ln_wd) THEN 
     471      DO jj = 2, jpjm1                          ! Remove coriolis term (and possibly spg) from barotropic trend 
     472         DO ji = fs_2, fs_jpim1 
     473             zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * umask(ji,jj,1) * wduflt1(ji,jj) 
     474             zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * vmask(ji,jj,1) * wdvflt1(ji,jj) 
     475          END DO 
     476      END DO 
     477      ELSE 
    389478      DO jj = 2, jpjm1                          ! Remove coriolis term (and possibly spg) from barotropic trend 
    390479         DO ji = fs_2, fs_jpim1 
     
    392481             zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * vmask(ji,jj,1) 
    393482          END DO 
    394       END DO  
     483      END DO 
     484      END IF 
    395485      ! 
    396486      !                 ! Add bottom stress contribution from baroclinic velocities:       
     
    416506      ! 
    417507      ! Note that the "unclipped" bottom friction parameter is used even with explicit drag 
    418       zu_frc(:,:) = zu_frc(:,:) + hur(:,:) * bfrua(:,:) * zwx(:,:) 
    419       zv_frc(:,:) = zv_frc(:,:) + hvr(:,:) * bfrva(:,:) * zwy(:,:) 
    420       !        
     508      IF(ln_wd) THEN 
     509        zu_frc(:,:) = zu_frc(:,:) + MAX(hur(:,:) * bfrua(:,:),-1._wp / rdtbt) * zwx(:,:) 
     510        zv_frc(:,:) = zv_frc(:,:) + MAX(hvr(:,:) * bfrva(:,:),-1._wp / rdtbt) * zwy(:,:) 
     511      ELSE 
     512        zu_frc(:,:) = zu_frc(:,:) + hur(:,:) * bfrua(:,:) * zwx(:,:) 
     513        zv_frc(:,:) = zv_frc(:,:) + hvr(:,:) * bfrva(:,:) * zwy(:,:) 
     514      END IF 
     515      ! 
    421516      IF (ln_bt_fw) THEN                        ! Add wind forcing 
    422517         zu_frc(:,:) =  zu_frc(:,:) + zraur * utau(:,:) * hur(:,:) 
     
    487582         vb_e   (:,:) = 0._wp 
    488583      ENDIF 
    489       ! 
    490       IF (ln_bt_fw) THEN                  ! FORWARD integration: start from NOW fields                     
    491          sshn_e(:,:) = sshn (:,:)             
    492          zun_e (:,:) = un_b (:,:)             
     584 
     585      IF(ln_wd) THEN      !preserve the positivity of water depth 
     586                          !ssh[b,n,a] should have already been processed for this 
     587         sshbb_e(:,:) = MAX(sshbb_e(:,:), rn_wdmin1 - bathy(:,:)) 
     588         sshb_e(:,:)  = MAX(sshb_e(:,:) , rn_wdmin1 - bathy(:,:)) 
     589      ENDIF 
     590      ! 
     591      IF (ln_bt_fw) THEN                  ! FORWARD integration: start from NOW fields 
     592         sshn_e(:,:) = sshn (:,:) 
     593         zun_e (:,:) = un_b (:,:) 
    493594         zvn_e (:,:) = vn_b (:,:) 
    494595         ! 
     
    561662            zhup2_e (:,:) = hu_0(:,:) + zwx(:,:)                ! Ocean depth at U- and V-points 
    562663            zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:) 
     664            IF(ln_wd) THEN 
     665              zhup2_e(:,:) = MAX(zhup2_e (:,:), rn_wdmin1) 
     666              zhvp2_e(:,:) = MAX(zhvp2_e (:,:), rn_wdmin1) 
     667            END IF 
    563668         ELSE 
    564669            zhup2_e (:,:) = hu(:,:) 
     
    599704         ENDIF 
    600705#endif 
     706         IF(ln_wd) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt) 
    601707         ! 
    602708         ! Sum over sub-time-steps to compute advective velocities 
     
    613719         END DO 
    614720         ssha_e(:,:) = (  sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) )  ) * tmask(:,:,1) 
     721         IF(ln_wd) ssha_e(:,:) = MAX(ssha_e(:,:), rn_wdmin1 - bathy(:,:))  
    615722         CALL lbc_lnk( ssha_e, 'T',  1._wp ) 
    616723 
     
    660767          &            + za2 *  sshb_e(:,:) + za3 *  sshbb_e(:,:) 
    661768 
     769         IF(ln_wd) THEN                    ! Calculating and applying W/D gravity filters 
     770           DO jj = 2, jpjm1 
     771              DO ji = 2, jpim1 
     772                 ll_tmp1 = MIN(zsshp2_e(ji,jj), zsshp2_e(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj))& 
     773                        & .and. MAX(zsshp2_e(ji,jj) + bathy(ji,jj), zsshp2_e(ji+1,jj) + bathy(ji+1,jj)) & 
     774                        &  > rn_wdmin1 + rn_wdmin2 
     775                 ll_tmp2 = MAX(zsshp2_e(ji,jj), zsshp2_e(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) + & 
     776                                                           & rn_wdmin1 + rn_wdmin2 
     777                 IF(ll_tmp1) THEN 
     778                   zcpx(ji,jj) = 1.0_wp 
     779                   wduflt1(ji,jj) = 1.0_wp 
     780                 ELSE IF(ll_tmp2) THEN 
     781                    ! no worries about zsshp2_e(ji+1,jj)-zsshp2_e(ji,jj) = 0, it won't happen ! here 
     782                   zcpx(ji,jj) = ABS((zsshp2_e(ji+1,jj) + bathy(ji+1,jj) - zsshp2_e(ji,jj) - bathy(ji,jj)) /& 
     783                               & (zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj))) 
     784                   wduflt1(ji,jj) = 1.0_wp 
     785                 ELSE 
     786                    zcpx(ji,jj) = 0._wp 
     787                    wduflt1(ji,jj) = 0.0_wp 
     788                 END IF 
     789 
     790                !for w/d debugging, delete it when finished 
     791                !   zcpx(ji,jj) = 0._wp 
     792                ! end debugging part 
     793 
     794                 ll_tmp1 = MIN(zsshp2_e(ji,jj), zsshp2_e(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1))& 
     795                        & .and. MAX(zsshp2_e(ji,jj) + bathy(ji,jj), zsshp2_e(ji,jj+1) + bathy(ji,jj+1)) & 
     796                        &  > rn_wdmin1 + rn_wdmin2 
     797                 ll_tmp2 = MAX(zsshp2_e(ji,jj), zsshp2_e(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) + & 
     798                                                           & rn_wdmin1 + rn_wdmin2 
     799                 IF(ll_tmp1) THEN 
     800                    zcpy(ji,jj) = 1.0_wp 
     801                    wdvflt1(ji,jj) = 1.0_wp 
     802                 ELSE IF(ll_tmp2) THEN 
     803                    ! no worries about zsshp2_e(ji,jj+1)-zsshp2_e(ji,jj) = 0, it won't happen ! here 
     804                   zcpy(ji,jj) = ABS((zsshp2_e(ji,jj+1) + bathy(ji,jj+1) - zsshp2_e(ji,jj) - bathy(ji,jj)) /& 
     805                               & (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj))) 
     806                    wdvflt1(ji,jj) = 1.0_wp 
     807                 ELSE 
     808                   zcpy(ji,jj) = 0._wp 
     809                    wdvflt1(ji,jj) = 0.0_wp 
     810                 END IF 
     811 
     812                !for w/d debugging, delete it when finished 
     813                !   zcpy(ji,jj) = 0._wp 
     814                ! end debugging part 
     815              END DO 
     816            END DO 
     817            CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp ) 
     818         ENDIF 
    662819         ! 
    663820         ! Compute associated depths at U and V points: 
     
    676833               END DO 
    677834            END DO 
     835 
     836            IF(ln_wd) THEN 
     837              zhust_e(:,:) = MAX(zhust_e (:,:), rn_wdmin1 ) 
     838              zhvst_e(:,:) = MAX(zhvst_e (:,:), rn_wdmin1 ) 
     839            END IF 
     840            !shall we call lbc_lnk for zhu(v)st_e() here? 
     841 
    678842         ENDIF 
    679843         ! 
     
    742906         ! 
    743907         ! Surface pressure trend: 
    744          DO jj = 2, jpjm1 
    745             DO ji = fs_2, fs_jpim1   ! vector opt. 
    746                ! Add surface pressure gradient 
    747                zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj) 
    748                zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj) 
    749                zwx(ji,jj) = zu_spg 
    750                zwy(ji,jj) = zv_spg 
    751             END DO 
    752          END DO 
     908 
     909         IF(ln_wd) THEN 
     910           DO jj = 2, jpjm1 
     911              DO ji = 2, jpim1  
     912                 ! Add surface pressure gradient 
     913                 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj) 
     914                 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj) 
     915                 ! for W/D debugging 
     916                 !zwx(ji,jj) = zu_spg * zcpx(ji,jj) * 0._wp 
     917                 !zwy(ji,jj) = zv_spg * zcpy(ji,jj) * 0._wp 
     918                 IF(ji+nimpp-1 == jidbg .and. jj+njmpp-1 == jjdbg) THEN 
     919                   write(444,*) 'in spt_ts_intg:  zu_spg, zv_spg, zcpx, zcpy' 
     920                   write(444,*) zu_spg, zv_spg, zcpx(ji,jj), zcpy(ji,jj) 
     921                   write(445,*) 'in spt_ts_intg:  sshn_ij, sshn_ip1j, sshn_ijp1, bathy_ij, bathy_ip1j, bathy_ijp1' 
     922                   write(445,*) zsshp2_e(ji,jj), zsshp2_e(ji+1,jj), zsshp2_e(ji,jj+1), & 
     923                    &           bathy(ji,jj), bathy(ji+1,jj), bathy(ji,jj+1)  
     924                 END IF 
     925                 zwx(ji,jj) = zu_spg * zcpx(ji,jj)  
     926                 zwy(ji,jj) = zv_spg * zcpy(ji,jj) 
     927              END DO 
     928           END DO 
     929         ELSE 
     930           DO jj = 2, jpjm1 
     931              DO ji = fs_2, fs_jpim1   ! vector opt. 
     932                 ! Add surface pressure gradient 
     933                 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj) 
     934                 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj) 
     935                 zwx(ji,jj) = zu_spg 
     936                 zwy(ji,jj) = zv_spg 
     937              END DO 
     938           END DO 
     939         END IF 
     940 
    753941         ! 
    754942         ! Set next velocities: 
     
    774962               DO ji = fs_2, fs_jpim1   ! vector opt. 
    775963 
    776                   zhura = umask(ji,jj,1)/(hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - umask(ji,jj,1)) 
    777                   zhvra = vmask(ji,jj,1)/(hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - vmask(ji,jj,1)) 
    778  
    779                   ua_e(ji,jj) = (                hu_e(ji,jj)  *  zun_e(ji,jj)   &  
    780                             &     + rdtbt * ( zhust_e(ji,jj)  *    zwx(ji,jj)   &  
     964                  IF(ln_wd) THEN 
     965                    zhura = MAX(hu_0(ji,jj) + zsshu_a(ji,jj), rn_wdmin1) 
     966                    zhvra = MAX(hv_0(ji,jj) + zsshv_a(ji,jj), rn_wdmin1) 
     967                  ELSE 
     968                    zhura = hu_0(ji,jj) + zsshu_a(ji,jj) 
     969                    zhvra = hv_0(ji,jj) + zsshv_a(ji,jj) 
     970                  END IF 
     971 
     972                  zhura = umask(ji,jj,1)/(zhura + 1._wp - umask(ji,jj,1)) 
     973                  zhvra = vmask(ji,jj,1)/(zhvra + 1._wp - vmask(ji,jj,1)) 
     974 
     975                  ua_e(ji,jj) = (                hu_e(ji,jj)  *  zun_e(ji,jj)   & 
     976                            &     + rdtbt * ( zhust_e(ji,jj)  *    zwx(ji,jj)   & 
    781977                            &               + zhup2_e(ji,jj)  * zu_trd(ji,jj)   & 
    782978                            &               +      hu(ji,jj)  * zu_frc(ji,jj) ) & 
     
    793989         ! 
    794990         IF( lk_vvl ) THEN                             !* Update ocean depth (variable volume case only) 
    795             !                                          !  ----------------------------------------------         
    796             hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 
    797             hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 
     991            !                                          !  ---------------------------------------------- 
     992            IF(ln_wd) THEN 
     993              hu_e (:,:) = MAX(hu_0(:,:) + zsshu_a(:,:), rn_wdmin1) 
     994              hv_e (:,:) = MAX(hv_0(:,:) + zsshv_a(:,:), rn_wdmin1) 
     995            ELSE 
     996              hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 
     997              hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 
     998            END IF 
     999 
    7981000            hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1._wp - umask(:,:,1) ) 
    7991001            hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1._wp - vmask(:,:,1) ) 
     
    9251127      CALL wrk_dealloc( jpi, jpj, zsshu_a, zsshv_a                                   ) 
    9261128      CALL wrk_dealloc( jpi, jpj, zhf ) 
     1129      IF(ln_wd) CALL wrk_dealloc( jpi, jpj, zcpx, zcpy, wduflt1, wdvflt1 ) 
    9271130      ! 
    9281131      IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_ts') 
  • branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90

    r5656 r5842  
    99   !!             -   !  2010-09  (D.Storkey and E.O'Dea) bug fixes for BDY module 
    1010   !!            3.3  !  2011-10  (M. Leclair) split former ssh_wzv routine and remove all vvl related work 
     11   !!            3.?  !  2015-11  (H. Liu) add wetting and drying   
    1112   !!---------------------------------------------------------------------- 
    1213 
     
    3839   USE wrk_nemo        ! Memory Allocation 
    3940   USE timing          ! Timing 
     41   USE wadlmt          ! Wetting/Drying flux limting 
    4042 
    4143   IMPLICIT NONE 
     
    9092      ENDIF 
    9193      ! 
    92       CALL div_cur( kt )                              ! Horizontal divergence & Relative vorticity 
    93       ! 
    9494      z2dt = 2._wp * rdt                              ! set time step size (Euler/Leapfrog) 
    9595      IF( neuler == 0 .AND. kt == nit000 )   z2dt = rdt 
     96      ! 
     97      z1_rau0 = 0.5_wp * r1_rau0 
     98      ! 
     99      IF(ln_wd) CALL wad_lmt(sshb, z1_rau0 * (emp_b(:,:) + emp(:,:)), z2dt) 
     100 
     101      CALL div_cur( kt )                              ! Horizontal divergence & Relative vorticity 
     102      ! 
    96103 
    97104      !                                           !------------------------------! 
     
    106113      ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations. 
    107114      !  
    108       z1_rau0 = 0.5_wp * r1_rau0 
    109115      ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * ssmask(:,:) 
    110116 
  • branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90

    r5656 r5842  
    3030   !!            3.4  ! 2011-11  (C. Harris) decomposition changes for running with CICE 
    3131   !!                 ! 2012-05  (C. Calone, J. Simeon, G. Madec, C. Ethe) Add grid coarsening  
     32   !!            3.6.?! 2015-11  (H. Liu) Add Wetting and Drying  
    3233   !!---------------------------------------------------------------------- 
    3334 
     
    232233      NAMELIST/namcfg/ cp_cfg, cp_cfz, jp_cfg, jpidta, jpjdta, jpkdta, jpiglo, jpjglo, & 
    233234         &             jpizoom, jpjzoom, jperio, ln_use_jattr 
     235      NAMELIST/namwad/ ln_wd, rn_wdmin1, rn_wdmin2, rn_wdld, nn_wdit 
    234236      !!---------------------------------------------------------------------- 
    235237      ! 
     
    257259      READ  ( numnam_cfg, namcfg, IOSTAT = ios, ERR = 904 ) 
    258260904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcfg in configuration namelist', .TRUE. )    
     261 
     262      REWIND( numnam_ref )              ! Namelist namwad in reference namelist : Parameters for Wetting/Drying 
     263      READ  ( numnam_ref, namwad, IOSTAT = ios, ERR = 905) 
     264905   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namwad in reference namelist', .TRUE. ) 
     265 
     266      REWIND( numnam_cfg )              ! Namelist namwad in configuration namelist : Parameters for Wetting/Drying 
     267      READ  ( numnam_cfg, namwad, IOSTAT = ios, ERR = 906) 
     268906   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namwad in configuration namelist', .TRUE. ) 
     269 
     270 
    259271 
    260272! Force values for AGRIF zoom (cf. agrif_user.F90) 
     
    315327         WRITE( numond, namctl ) 
    316328         WRITE( numond, namcfg ) 
     329         WRITE( numond, namwad ) 
    317330      ENDIF 
    318331 
Note: See TracChangeset for help on using the changeset viewer.