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 5842 for branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

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

Wetting and Drying update based on r:5803

Location:
branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • 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) 
Note: See TracChangeset for help on using the changeset viewer.