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

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

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

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

Location:
trunk/NEMOGCM/NEMO/OPA_SRC/DOM
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90

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

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