Changeset 5014


Ignore:
Timestamp:
2015-01-07T19:03:53+01:00 (6 years ago)
Author:
hliu
Message:

upload the modifications for W/D based on r:4826

Location:
branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC
Files:
9 edited

Legend:

Unmodified
Added
Removed
  • branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90

    r4671 r5014  
    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 
     
    257258   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: tmask, umask, vmask, fmask   !: land/ocean mask at T-, U-, V- and F-pts 
    258259 
    259    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   tpol, fpol          !: north fold mask (jperio= 3 or 4) 
     260   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   tpol, fpol           !: north fold mask (jperio= 3 or 4) 
     261 
    260262 
    261263#if defined key_noslip_accurate 
     
    263265   INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  :: nicoa, njcoa       !: ??? 
    264266#endif 
     267 
     268   !!---------------------------------------------------------------------- 
     269   !! critical depths,limiters,and masks for  Wetting and Drying 
     270   !! --------------------------------------------------------------------- 
     271 
     272   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wdmask             !: u- and v- limiter  
     273 
     274   LOGICAL,  PUBLIC, SAVE ::   ln_wd       !: key to turn on/off wetting/drying (T: on, F: off) 
     275   REAL(wp), PUBLIC, SAVE ::   rn_wdmin1   !: minimum water depth on dried cells 
     276   REAL(wp), PUBLIC, SAVE ::   rn_wdmin2   !: tolerrance of minimum water depth on dried cells 
     277   REAL(wp), PUBLIC, SAVE ::   rn_wdld     !: land elevation below which wetting/drying will be considered 
     278   INTEGER , PUBLIC, SAVE ::   nn_wdit     !: maximum number of iteration for W/D limiter 
     279 
     280   !LOGICAL,  PUBLIC, SAVE ::   ln_wd     =  .FALSE.  !: turn on wetting/drying (T: on, F: off) 
     281   !REAL(wp), PUBLIC, SAVE ::   rn_wdmin1 = 0.10_wp   !: minimum water depth on dried cells 
     282   !REAL(wp), PUBLIC, SAVE ::   rn_wdmin2 = 0.01_wp   !: tolerrance of minimum water depth on dried cells 
     283   !REAL(wp), PUBLIC, SAVE ::   rn_wdld   = 20.0_wp   !: land elevation below which wetting/drying will be considered 
     284   !INTEGER , PUBLIC, SAVE ::   nn_wdit   =   10      !: maximum number of iteration for W/D limiter 
    265285 
    266286   !!---------------------------------------------------------------------- 
     
    326346   INTEGER FUNCTION dom_oce_alloc() 
    327347      !!---------------------------------------------------------------------- 
    328       INTEGER, DIMENSION(11) :: ierr 
     348      INTEGER, DIMENSION(12) :: ierr 
    329349      !!---------------------------------------------------------------------- 
    330350      ierr(:) = 0 
     
    391411      ALLOCATE( npcoa(4,jpk), nicoa(2*(jpi+jpj),4,jpk), njcoa(2*(jpi+jpj),4,jpk), STAT=ierr(11) ) 
    392412#endif 
     413 
     414      IF(ln_wd) & 
     415      ALLOCATE( wdmask(jpi,jpj), STAT=ierr(12) ) 
    393416      ! 
    394417      dom_oce_alloc = MAXVAL(ierr) 
  • branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90

    r4624 r5014  
    8686                             CALL dom_zgr      ! Vertical mesh and bathymetry 
    8787                             CALL dom_msk      ! Masks 
    88       IF( ln_sco )           CALL dom_stiff    ! Maximum stiffness ratio/hydrostatic consistency 
     88      IF( ln_sco.AND.(.NOT.ln_wd) ) & 
     89                         &   CALL dom_stiff    ! Maximum stiffness ratio/hydrostatic consistency 
    8990      ! 
    9091      ht_0(:,:) = 0.0_wp                       ! Reference ocean depth at T-points 
  • branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90

    r4795 r5014  
    145145      fse3t_a(:,:,jpk) = e3t_0(:,:,jpk) 
    146146 
     147      IF(ln_wd) THEN 
     148        DO jj = 1, jpj 
     149          DO ji = 1, jpi 
     150            IF(e3t_0(ji,jj,1) < 0._wp) fse3t_a(ji,jj,1:2) = 0.5_wp * rn_wdmin1  
     151          ENDDO 
     152        ENDDO 
     153      END IF 
     154 
    147155      ! Reconstruction of all vertical scale factors at now and before time steps 
    148156      ! ============================================================================= 
     
    793801      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag 
    794802      !! * Local declarations 
    795       INTEGER ::   jk 
     803      INTEGER ::   ji, jj, jk 
    796804      INTEGER ::   id1, id2, id3, id4, id5     ! local integers 
    797805      !!---------------------------------------------------------------------- 
     
    867875            fse3t_n(:,:,:) = e3t_0(:,:,:) 
    868876            sshn(:,:) = 0.0_wp 
     877 
     878            IF(ln_wd) THEN 
     879              DO jj = 1, jpj 
     880                DO ji = 1, jpi 
     881                  IF(e3t_0(ji,jj,1) < 0._wp) THEN 
     882                    fse3t_b(ji,jj,1:2) = 0.5_wp * rn_wdmin1  
     883                    fse3t_n(ji,jj,1:2) = 0.5_wp * rn_wdmin1  
     884                    sshb(ji,jj) = rn_wdmin1 - bathy(ji,jj) 
     885                    sshn(ji,jj) = rn_wdmin1 - bathy(ji,jj) 
     886                  ENDIF 
     887                ENDDO 
     888              ENDDO 
     889            END IF 
     890 
    869891            IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN 
    870892               tilde_e3t_b(:,:,:) = 0.0_wp 
  • branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90

    r4687 r5014  
    1717   !!            3.4  ! 2012-08  (J. Siddorn) added Siddorn and Furner stretching function 
    1818   !!            3.4  ! 2012-12  (R. Bourdalle-Badie and G. Reffray)  modify C1D case   
     19   !!            3.6  ! 2014-09  (H. Liu) Modifications for Wetting/Drying 
    1920   !!---------------------------------------------------------------------- 
    2021 
     
    451452         END DO 
    452453         ! 
    453          !                                            ! ================ ! 
    454454      ELSEIF( ntopo == 1 ) THEN                       !   read in file   ! (over the local domain) 
    455455         !                                            ! ================ ! 
     
    11211121      REAL(wp) ::   zrmax, ztaper   ! temporary scalars 
    11221122      REAL(wp) ::   zrfact 
     1123      REAL(wp) ::   zsmth 
    11231124      ! 
    11241125      REAL(wp), POINTER, DIMENSION(:,:  ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2 
     
    11761177      bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) ) 
    11771178 
     1179      IF(.NOT.ln_wd) THEN                       
    11781180      DO jj = 1, jpj 
    11791181         DO ji = 1, jpi 
     
    11811183         END DO 
    11821184      END DO 
     1185      END IF 
    11831186      !                                        ! ============================= 
    11841187      !                                        ! Define the envelop bathymetry   (hbatt) 
     
    11871190      zenv(:,:) = bathy(:,:) 
    11881191      ! 
     1192      IF(.NOT.ln_wd) THEN     
    11891193      ! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing 
    11901194      DO jj = 1, jpj 
     
    12021206         END DO 
    12031207      END DO 
     1208      END IF 
     1209 
    12041210      ! apply lateral boundary condition   CAUTION: keep the value when the lbc field is zero 
    12051211      CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' ) 
     
    12411247         zri(:,:) = 0._wp 
    12421248         zrj(:,:) = 0._wp 
     1249          
     1250         !IF(ln_wd) THEN                     !extend the smoothed region to cover the W/D zones 
     1251         !  zsmth = -rn_wdld     
     1252         !ELSE 
     1253           zsmth = 0._wp                     ! The original form (only smooth ocean points) 
     1254         !ENDIF 
     1255 
    12431256         DO jj = 1, nlcj 
    12441257            DO ji = 1, nlci 
    12451258               iip1 = MIN( ji+1, nlci )      ! force zri = 0 on last line (ji=ncli+1 to jpi) 
    12461259               ijp1 = MIN( jj+1, nlcj )      ! force zrj = 0 on last raw  (jj=nclj+1 to jpj) 
    1247                IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(iip1,jj) > 0._wp)) THEN 
     1260               IF( (zenv(ji,jj) > zsmth) .AND. (zenv(iip1,jj) > zsmth)) THEN 
    12481261                  zri(ji,jj) = ( zenv(iip1,jj  ) - zenv(ji,jj) ) / ( zenv(iip1,jj  ) + zenv(ji,jj) ) 
    12491262               END IF 
    1250                IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(ji,ijp1) > 0._wp)) THEN 
     1263               IF( (zenv(ji,jj) > zsmth) .AND. (zenv(ji,ijp1) > zsmth)) THEN 
    12511264                  zrj(ji,jj) = ( zenv(ji  ,ijp1) - zenv(ji,jj) ) / ( zenv(ji  ,ijp1) + zenv(ji,jj) ) 
    12521265               END IF 
     
    12571270            END DO 
    12581271         END DO 
     1272 
    12591273         IF( lk_mpp )   CALL mpp_max( zrmax )   ! max over the global domain 
    12601274         ! 
     
    12711285      END DO                                                !     End loop     ! 
    12721286      !                                                     ! ================ ! 
    1273       DO jj = 1, jpj 
    1274          DO ji = 1, jpi 
    1275             zenv(ji,jj) = MAX( zenv(ji,jj), rn_sbot_min ) ! set all points to avoid undefined scale value warnings 
    1276          END DO 
    1277       END DO 
     1287      IF(ln_wd) THEN 
     1288        DO jj = 1, jpj 
     1289           DO ji = 1, jpi 
     1290              zenv(ji,jj) = MAX( zenv(ji,jj), -rn_wdld )    ! filt out land bathy data  
     1291           END DO 
     1292        END DO 
     1293      ELSE 
     1294        DO jj = 1, jpj 
     1295           DO ji = 1, jpi 
     1296              zenv(ji,jj) = MAX( zenv(ji,jj), rn_sbot_min ) ! set all points to avoid undefined scale value warnings 
     1297           END DO 
     1298        END DO 
     1299      ENDIF 
    12781300      ! 
    12791301      ! Envelope bathymetry saved in hbatt 
     
    13051327      IF(lwp) THEN 
    13061328         WRITE(numout,*) 
    1307          WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min 
     1329         IF(.NOT.ln_wd) THEN 
     1330           WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min 
     1331         ELSE 
     1332           WRITE(numout,*) ' zgr_sco: minimum positive depth of the envelop topography set to : ', rn_sbot_min 
     1333           WRITE(numout,*) ' zgr_sco: minimum negative depth of the envelop topography set to : ', -rn_wdld 
     1334         ENDIF 
    13081335      ENDIF 
    13091336      hbatu(:,:) = rn_sbot_min 
     
    13181345        END DO 
    13191346      END DO 
     1347 
     1348      IF(ln_wd) THEN               !avoid the zero depth on T- (U-,V-,F-) points 
     1349        DO jj = 1, jpj 
     1350          DO ji = 1, jpi 
     1351            IF(ABS(hbatt(ji,jj)) < rn_wdmin1) & 
     1352              & hbatt(ji,jj) = SIGN(1._wp, hbatt(ji,jj)) * rn_wdmin1 
     1353            IF(ABS(hbatu(ji,jj)) < rn_wdmin1) & 
     1354              & hbatu(ji,jj) = SIGN(1._wp, hbatu(ji,jj)) * rn_wdmin1 
     1355            IF(ABS(hbatv(ji,jj)) < rn_wdmin1) & 
     1356              & hbatv(ji,jj) = SIGN(1._wp, hbatv(ji,jj)) * rn_wdmin1 
     1357            IF(ABS(hbatf(ji,jj)) < rn_wdmin1) & 
     1358              & hbatf(ji,jj) = SIGN(1._wp, hbatf(ji,jj)) * rn_wdmin1 
     1359          END DO 
     1360        END DO 
     1361      END IF 
    13201362      !  
    13211363      ! Apply lateral boundary condition 
    13221364!!gm  ! CAUTION: retain non zero value in the initial file this should be OK for orca cfg, not for EEL 
     1365 
    13231366      zhbat(:,:) = hbatu(:,:)   ;   CALL lbc_lnk( hbatu, 'U', 1._wp ) 
    13241367      DO jj = 1, jpj 
    13251368         DO ji = 1, jpi 
    13261369            IF( hbatu(ji,jj) == 0._wp ) THEN 
     1370               !No worries about the following line when ln_wd == .true. 
    13271371               IF( zhbat(ji,jj) == 0._wp )   hbatu(ji,jj) = rn_sbot_min 
    13281372               IF( zhbat(ji,jj) /= 0._wp )   hbatu(ji,jj) = zhbat(ji,jj) 
     
    13301374         END DO 
    13311375      END DO 
     1376 
    13321377      zhbat(:,:) = hbatv(:,:)   ;   CALL lbc_lnk( hbatv, 'V', 1._wp ) 
    13331378      DO jj = 1, jpj 
     
    13391384         END DO 
    13401385      END DO 
     1386 
    13411387      zhbat(:,:) = hbatf(:,:)   ;   CALL lbc_lnk( hbatf, 'F', 1._wp ) 
    13421388      DO jj = 1, jpj 
     
    13501396 
    13511397!!bug:  key_helsinki a verifer 
    1352       hift(:,:) = MIN( hift(:,:), hbatt(:,:) ) 
    1353       hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) ) 
    1354       hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) ) 
    1355       hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) ) 
     1398      IF(.NOT.ln_wd) THEN 
     1399       hift(:,:) = MIN( hift(:,:), hbatt(:,:) ) 
     1400       hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) ) 
     1401       hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) ) 
     1402       hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) ) 
     1403      END IF 
    13561404 
    13571405      IF( nprint == 1 .AND. lwp )   THEN 
     
    13951443      CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 
    13961444 
    1397       fsdepw(:,:,:) = gdepw_0 (:,:,:) 
    1398       fsde3w(:,:,:) = gdep3w_0(:,:,:) 
    1399       ! 
    1400       where (e3t_0   (:,:,:).eq.0.0)  e3t_0(:,:,:) = 1.0 
    1401       where (e3u_0   (:,:,:).eq.0.0)  e3u_0(:,:,:) = 1.0 
    1402       where (e3v_0   (:,:,:).eq.0.0)  e3v_0(:,:,:) = 1.0 
    1403       where (e3f_0   (:,:,:).eq.0.0)  e3f_0(:,:,:) = 1.0 
    1404       where (e3w_0   (:,:,:).eq.0.0)  e3w_0(:,:,:) = 1.0 
    1405       where (e3uw_0  (:,:,:).eq.0.0)  e3uw_0(:,:,:) = 1.0 
    1406       where (e3vw_0  (:,:,:).eq.0.0)  e3vw_0(:,:,:) = 1.0 
     1445      !fsdepw(:,:,:) = gdepw_0 (:,:,:) 
     1446      !fsde3w(:,:,:) = gdep3w_0(:,:,:) 
     1447      ! 
     1448      IF(.NOT.ln_wd) THEN 
     1449        where (e3t_0   (:,:,:).eq.0.0)  e3t_0(:,:,:) = 1.0 
     1450        where (e3u_0   (:,:,:).eq.0.0)  e3u_0(:,:,:) = 1.0 
     1451        where (e3v_0   (:,:,:).eq.0.0)  e3v_0(:,:,:) = 1.0 
     1452        where (e3f_0   (:,:,:).eq.0.0)  e3f_0(:,:,:) = 1.0 
     1453        where (e3w_0   (:,:,:).eq.0.0)  e3w_0(:,:,:) = 1.0 
     1454        where (e3uw_0  (:,:,:).eq.0.0)  e3uw_0(:,:,:) = 1.0 
     1455        where (e3vw_0  (:,:,:).eq.0.0)  e3vw_0(:,:,:) = 1.0 
     1456      END IF 
    14071457 
    14081458#if defined key_agrif 
     
    14461496               IF( scobot(ji,jj) >= fsdept(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk ) 
    14471497            END DO 
    1448             IF( scobot(ji,jj) == 0._wp               )   mbathy(ji,jj) = 0 
    1449          END DO 
    1450       END DO 
     1498 
     1499            IF(ln_wd) THEN 
     1500              IF( scobot(ji,jj) <= -(rn_wdld - rn_wdmin2)) THEN 
     1501                mbathy(ji,jj) = 0 
     1502              ELSEIF(scobot(ji,jj) <= rn_wdmin1) THEN 
     1503                mbathy(ji,jj) = 2 
     1504              ENDIF 
     1505            ELSE 
     1506              IF( scobot(ji,jj) == 0._wp )   mbathy(ji,jj) = 0 
     1507            ENDIF 
     1508         END DO 
     1509      END DO 
     1510 
    14511511      IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ),   & 
    14521512         &                                                       ' MAX ', MAXVAL( mbathy(:,:) ) 
     
    15681628      INTEGER  ::   ji, jj, jk           ! dummy loop argument 
    15691629      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars 
     1630      REAL(wp) ::   ztmpu, ztmpv, ztmpf 
    15701631      ! 
    15711632      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 
     
    16211682      DO ji = 1, jpim1 
    16221683         DO jj = 1, jpjm1 
     1684            ! extended for Wetting/Drying case 
     1685            ztmpu = hbatt(ji,jj)+hbatt(ji+1,jj) 
     1686            ztmpv = hbatt(ji,jj)+hbatt(ji,jj+1) 
     1687            ztmpf = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) 
    16231688            DO jk = 1, jpk 
    1624                z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) )   & 
    1625                   &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
    1626                z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) )   & 
    1627                   &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
    1628                z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk)     & 
    1629                   &                + hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) )   & 
    1630                   &              / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 
    1631                z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) )   & 
    1632                   &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
    1633                z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) )   & 
    1634                   &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
     1689               IF(ABS(ztmpu) < rn_wdmin1.AND.ln_wd) THEN 
     1690                 z_esigtu3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) ) 
     1691                 z_esigwu3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji+1,jj,jk) ) 
     1692               ELSE 
     1693                 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) & 
     1694                        &              / ztmpu 
     1695                 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) & 
     1696                        &              / ztmpu 
     1697               END IF 
     1698 
     1699               IF(ABS(ztmpv) < rn_wdmin1.AND.ln_wd) THEN 
     1700                 z_esigtv3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji,jj+1,jk) ) 
     1701                 z_esigwv3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji,jj+1,jk) ) 
     1702               ELSE 
     1703                 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) & 
     1704                        &              / ztmpv 
     1705                 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) & 
     1706                        &              / ztmpv 
     1707               END IF 
     1708 
     1709               IF(ABS(ztmpf) < rn_wdmin1.AND.ln_wd) THEN 
     1710                 z_esigtf3(ji,jj,jk) = 0.25_wp * ( z_esigt3(ji,jj,jk)   + z_esigt3(ji+1,jj,jk)  & 
     1711                    &                            + z_esigt3(ji,jj+1,jk) + z_esigt3(ji+1,jj+1,jk) ) 
     1712               ELSE 
     1713                 z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk)  & 
     1714                    &                + hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) & 
     1715                    &              / ztmpf 
     1716               END IF 
     1717 
     1718               ! 
     1719               !z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) )   & 
     1720               !   &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
     1721               !z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) )   & 
     1722               !   &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
     1723               !z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk)     & 
     1724               !   &                + hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) )   & 
     1725               !   &              / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 
     1726               !z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) )   & 
     1727               !   &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
     1728               !z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) )   & 
     1729               !   &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
    16351730               ! 
    16361731               e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 
     
    16721767      REAL(wp) ::   zsmth               ! smoothing around critical depth 
    16731768      REAL(wp) ::   zzs, zzb           ! Surface and bottom cell thickness in sigma space 
     1769      REAL(wp) ::   ztmpu, ztmpv, ztmpf 
    16741770      ! 
    16751771      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 
     
    17501846        DO jj=1,jpj-1 
    17511847 
    1752           DO jk = 1, jpk 
    1753                 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) / & 
    1754                                     ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
    1755                 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) / & 
    1756                                     ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
    1757                 z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) +  & 
    1758                                       hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / & 
    1759                                     ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 
    1760                 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) / & 
    1761                                     ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
    1762                 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) / & 
    1763                                     ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
     1848           ! extend to suit for Wetting/Drying case 
     1849           ztmpu = hbatt(ji,jj)+hbatt(ji+1,jj) 
     1850           ztmpv = hbatt(ji,jj)+hbatt(ji,jj+1) 
     1851           ztmpf = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) 
     1852           DO jk = 1, jpk 
     1853              IF(ABS(ztmpu) < 1.e-5.AND.ln_wd) THEN 
     1854                z_esigtu3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) ) 
     1855                z_esigwu3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji+1,jj,jk) ) 
     1856              ELSE 
     1857                z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) & 
     1858                       &              / ztmpu 
     1859                z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) & 
     1860                       &              / ztmpu 
     1861              END IF 
     1862 
     1863              IF(ABS(ztmpv) < 1.e-5.AND.ln_wd) THEN 
     1864                z_esigtv3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji,jj+1,jk) ) 
     1865                z_esigwv3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji,jj+1,jk) ) 
     1866              ELSE 
     1867                z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) & 
     1868                       &              / ztmpv 
     1869                z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) & 
     1870                       &              / ztmpv 
     1871              END IF 
     1872 
     1873              IF(ABS(ztmpf) < 1.e-5.AND.ln_wd) THEN 
     1874                z_esigtf3(ji,jj,jk) = 0.25_wp * ( z_esigt3(ji,jj,jk)   + z_esigt3(ji+1,jj,jk)  & 
     1875                   &                            + z_esigt3(ji,jj+1,jk) + z_esigt3(ji+1,jj+1,jk) ) 
     1876              ELSE 
     1877                z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk)  & 
     1878                   &                + hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) & 
     1879                   &              / ztmpf 
     1880              END IF 
     1881 
     1882             !z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) / & 
     1883             !                    ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
     1884             !z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) / & 
     1885             !                    ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
     1886             !z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) +  & 
     1887             !                      hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / & 
     1888             !                    ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 
     1889             !z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) / & 
     1890             !                    ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
     1891             !z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) / & 
     1892             !                    ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
    17641893 
    17651894             e3t_0(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk) 
  • branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90

    r4624 r5014  
    1616   !!            3.4  !  2011-11  (H. Liu) hpg_prj: Original code for s-coordinates 
    1717   !!                 !           (A. Coward) suppression of hel, wdj and rot options 
     18   !!            3.6? !  2014-09  (H. Liu) add Wetting/Drying pressure filter 
    1819   !!---------------------------------------------------------------------- 
    1920 
     
    369370      !! 
    370371      INTEGER  ::   ji, jj, jk                 ! dummy loop indices 
    371       REAL(wp) ::   zcoef0, zuap, zvap, znad   ! temporary scalars 
     372      REAL(wp) ::   zcoef0, zuap, zvap, znad, ztmp   ! temporary scalars 
     373      LOGICAL  ::   ll_tmp1, ll_tmp2           ! local logical variables 
    372374      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zhpi, zhpj 
     375      REAL(wp), POINTER, DIMENSION(:,:)   ::  zcpx, zcpy    !W/D pressure filter 
    373376      !!---------------------------------------------------------------------- 
    374377      ! 
    375378      CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj ) 
     379      IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 
    376380      ! 
    377381      IF( kt == nit000 ) THEN 
     
    386390      IF ( lk_vvl ) THEN   ;     znad = 1._wp          ! Variable volume 
    387391      ELSE                 ;     znad = 0._wp         ! Fixed volume 
     392      ENDIF 
     393 
     394      IF(ln_wd) THEN 
     395        DO jj = 2, jpjm1 
     396           DO ji = 2, jpim1  
     397             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) 
     398             ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) +& 
     399                                                       & rn_wdmin1 + rn_wdmin2 
     400 
     401             IF(ll_tmp1) THEN 
     402               zcpx(ji,jj) = 1.0_wp 
     403             ELSE IF(ll_tmp2) THEN 
     404               ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 
     405               zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /& 
     406                           &     (sshn(ji+1,jj) - sshn(ji,jj))) 
     407             ELSE 
     408               zcpx(ji,jj) = 0._wp 
     409             END IF 
     410       
     411             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) 
     412             ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) +& 
     413                                                       & rn_wdmin1 + rn_wdmin2 
     414 
     415             IF(ll_tmp1) THEN 
     416               zcpy(ji,jj) = 1.0_wp 
     417             ELSE IF(ll_tmp2) THEN 
     418               ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 
     419               zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) /& 
     420                           &     (sshn(ji,jj+1) - sshn(ji,jj))) 
     421             ELSE 
     422               zcpy(ji,jj) = 0._wp 
     423             END IF 
     424           END DO 
     425        END DO 
     426        CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp ) 
    388427      ENDIF 
    389428 
     
    401440            zvap = -zcoef0 * ( rhd   (ji,jj+1,1) + rhd   (ji,jj,1) + 2._wp * znad )   & 
    402441               &           * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj) 
     442 
     443            IF(ln_wd) THEN 
     444              zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 
     445              zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj)  
     446              zuap = zuap * zcpx(ji,jj) 
     447              zvap = zvap * zcpy(ji,jj) 
     448            ENDIF 
     449 
    403450            ! add to the general momentum trend 
    404451            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap 
     
    423470               zvap = -zcoef0 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   & 
    424471                  &           * ( fsde3w(ji  ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) 
     472 
     473               IF(ln_wd) THEN 
     474                 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 
     475                 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj)  
     476                 zuap = zuap * zcpx(ji,jj) 
     477                 zvap = zvap * zcpy(ji,jj) 
     478               ENDIF 
     479 
    425480               ! add to the general momentum trend 
    426481               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap 
     
    431486      ! 
    432487      CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj ) 
     488      IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 
    433489      ! 
    434490   END SUBROUTINE hpg_sco 
     
    448504      REAL(wp) ::   z1_10, cffu, cffx   !    "         " 
    449505      REAL(wp) ::   z1_12, cffv, cffy   !    "         " 
     506      LOGICAL  ::   ll_tmp1, ll_tmp2    ! local logical variables 
    450507      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zhpi, zhpj 
    451508      REAL(wp), POINTER, DIMENSION(:,:,:) ::  dzx, dzy, dzz, dzu, dzv, dzw 
    452509      REAL(wp), POINTER, DIMENSION(:,:,:) ::  drhox, drhoy, drhoz, drhou, drhov, drhow 
    453510      REAL(wp), POINTER, DIMENSION(:,:,:) ::  rho_i, rho_j, rho_k 
     511      REAL(wp), POINTER, DIMENSION(:,:)   ::  zcpx, zcpy    !W/D pressure filter 
    454512      !!---------------------------------------------------------------------- 
    455513      ! 
     
    457515      CALL wrk_alloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow ) 
    458516      CALL wrk_alloc( jpi, jpj, jpk, rho_i, rho_j, rho_k,  zhpi,  zhpj        ) 
    459       ! 
     517      IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 
     518      ! 
     519      !!---------------------------------------------------------------------- 
     520      ! 
     521      ! 
     522      IF(ln_wd) THEN 
     523        DO jj = 2, jpjm1 
     524           DO ji = 2, jpim1  
     525             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) 
     526             ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) +& 
     527                                                       & rn_wdmin1 + rn_wdmin2 
     528 
     529             IF(ll_tmp1) THEN 
     530               zcpx(ji,jj) = 1.0_wp 
     531             ELSE IF(ll_tmp2) THEN 
     532               ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 
     533               zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /& 
     534                           &     (sshn(ji+1,jj) - sshn(ji,jj))) 
     535             ELSE 
     536               zcpx(ji,jj) = 0._wp 
     537             END IF 
     538       
     539             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) 
     540             ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) +& 
     541                                                       & rn_wdmin1 + rn_wdmin2 
     542 
     543             IF(ll_tmp1) THEN 
     544               zcpy(ji,jj) = 1.0_wp 
     545             ELSE IF(ll_tmp2) THEN 
     546               ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 
     547               zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) /& 
     548                           &     (sshn(ji,jj+1) - sshn(ji,jj))) 
     549             ELSE 
     550               zcpy(ji,jj) = 0._wp 
     551             END IF 
     552           END DO 
     553        END DO 
     554        CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp ) 
     555      ENDIF 
     556 
    460557 
    461558      IF( kt == nit000 ) THEN 
     
    628725            zhpi(ji,jj,1) = ( rho_k(ji+1,jj  ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) / e1u(ji,jj) 
    629726            zhpj(ji,jj,1) = ( rho_k(ji  ,jj+1,1) - rho_k(ji,jj,1) - rho_j(ji,jj,1) ) / e2v(ji,jj) 
     727            IF(ln_wd) THEN 
     728              zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 
     729              zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj)  
     730            ENDIF 
    630731            ! add to the general momentum trend 
    631732            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) 
     
    647748                  &           + (  ( rho_k(ji,jj+1,jk) - rho_k(ji,jj,jk  ) )    & 
    648749                  &               -( rho_j(ji,jj  ,jk) - rho_j(ji,jj,jk-1) )  ) / e2v(ji,jj) 
     750               IF(ln_wd) THEN 
     751                 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 
     752                 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj)  
     753               ENDIF 
    649754               ! add to the general momentum trend 
    650755               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) 
     
    657762      CALL wrk_dealloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow ) 
    658763      CALL wrk_dealloc( jpi, jpj, jpk, rho_i, rho_j, rho_k,  zhpi,  zhpj        ) 
     764      IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 
    659765      ! 
    660766   END SUBROUTINE hpg_djc 
     
    682788      !! The local variables for the correction term 
    683789      INTEGER  :: jk1, jis, jid, jjs, jjd 
     790      LOGICAL  :: ll_tmp1, ll_tmp2                  ! local logical variables 
    684791      REAL(wp) :: zuijk, zvijk, zpwes, zpwed, zpnss, zpnsd, zdeps 
    685792      REAL(wp) :: zrhdt1 
     
    687794      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zdept, zrhh 
    688795      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp 
    689       !!---------------------------------------------------------------------- 
     796      REAL(wp), POINTER, DIMENSION(:,:)   ::  zcpx, zcpy    !W/D pressure filter 
     797      !!---------------------------------------------------------------------- 
     798      ! 
    690799      ! 
    691800      CALL wrk_alloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp ) 
    692801      CALL wrk_alloc( jpi,jpj,jpk, zdept, zrhh ) 
     802      IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 
    693803      ! 
    694804      IF( kt == nit000 ) THEN 
     
    703813      znad = 0.0_wp 
    704814      IF( lk_vvl ) znad = 1._wp 
     815 
     816      IF(ln_wd) THEN 
     817        DO jj = 2, jpjm1 
     818           DO ji = 2, jpim1  
     819             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) 
     820             ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) +& 
     821                                                       & rn_wdmin1 + rn_wdmin2 
     822 
     823             IF(ll_tmp1) THEN 
     824               zcpx(ji,jj) = 1.0_wp 
     825             ELSE IF(ll_tmp2) THEN 
     826               ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 
     827               zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /& 
     828                           &     (sshn(ji+1,jj) - sshn(ji,jj))) 
     829             ELSE 
     830               zcpx(ji,jj) = 0._wp 
     831             END IF 
     832       
     833             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) 
     834             ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) +& 
     835                                                       & rn_wdmin1 + rn_wdmin2 
     836 
     837             IF(ll_tmp1.OR.ll_tmp2) THEN 
     838               zcpy(ji,jj) = 1.0_wp 
     839             ELSE IF(ll_tmp2) THEN 
     840               ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 
     841               zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) /& 
     842                           &     (sshn(ji,jj+1) - sshn(ji,jj))) 
     843             ELSE 
     844               zcpy(ji,jj) = 0._wp 
     845             END IF 
     846           END DO 
     847        END DO 
     848        CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp ) 
     849      ENDIF 
    705850 
    706851      ! Clean 3-D work arrays 
     
    8621007               ENDIF 
    8631008 
    864                ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * & 
    865                &           umask(ji,jj,jk) * tmask(ji,jj,jk) * tmask(ji+1,jj,jk) 
     1009               IF(ln_wd) THEN 
     1010                  zdpdx1 = zdpdx1 * zcpx(ji,jj) 
     1011                  zdpdx2 = zdpdx2 * zcpx(ji,jj) 
     1012                ENDIF 
     1013                ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk)  
    8661014            ENDIF 
    8671015 
     
    9191067               ENDIF 
    9201068 
    921                va(ji,jj,jk) = va(ji,jj,jk) + (zdpdy1 + zdpdy2)*& 
    922                &              vmask(ji,jj,jk)*tmask(ji,jj,jk)*tmask(ji,jj+1,jk) 
     1069               IF(ln_wd) THEN 
     1070                  zdpdy1 = zdpdy1 * zcpy(ji,jj) 
     1071                  zdpdy2 = zdpdy2 * zcpy(ji,jj) 
     1072                ENDIF 
     1073 
     1074               va(ji,jj,jk) = va(ji,jj,jk) + (zdpdy1 + zdpdy2) * vmask(ji,jj,jk) 
    9231075            ENDIF 
    9241076 
     
    9301082      CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp ) 
    9311083      CALL wrk_dealloc( jpi,jpj,jpk, zdept, zrhh ) 
     1084      IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 
    9321085      ! 
    9331086   END SUBROUTINE hpg_prj 
  • branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r4770 r5014  
    1111   !!             3.5  ! 2013-07  (J. Chanut) Switch to Forward-backward time stepping 
    1212   !!             3.6  ! 2013-11  (A. Coward) Update for z-tilde compatibility 
     13   !!             3.?  ! 2014-09  (H. Liu) Add Wetting/Drying components 
    1314   !!--------------------------------------------------------------------- 
    1415#if defined key_dynspg_ts   ||   defined key_esopa 
     
    1718   !!---------------------------------------------------------------------- 
    1819   !!   dyn_spg_ts  : compute surface pressure gradient trend using a time- 
    19    !!                 splitting scheme and add to the general trend  
     20   !!                 splitting scheme and add to the general trend 
    2021   !!---------------------------------------------------------------------- 
    2122   USE oce             ! ocean dynamics and tracers 
     
    2627   USE dynvor          ! vorticity term 
    2728   USE bdy_par         ! for lk_bdy 
    28    USE bdytides        ! open boundary condition data      
     29   USE bdytides        ! open boundary condition data 
    2930   USE bdydyn2d        ! open boundary conditions on barotropic variables 
    3031   USE sbctide         ! tides 
     
    3839   USE zdf_oce         ! Bottom friction coefts 
    3940   USE wrk_nemo        ! Memory Allocation 
    40    USE timing          ! Timing     
     41   USE timing          ! Timing 
    4142   USE sbcapr          ! surface boundary condition: atmospheric pressure 
     43   USE wadlmt          ! wetting/drying flux limter 
    4244   USE dynadv, ONLY: ln_dynadv_vec 
    4345#if defined key_agrif 
    4446   USE agrif_opa_interp ! agrif 
    4547#endif 
    46 #if defined key_asminc    
     48#if defined key_asminc 
    4749   USE asminc          ! Assimilation increment 
    4850#endif 
     
    5153   PRIVATE 
    5254 
    53    PUBLIC dyn_spg_ts        ! routine called in dynspg.F90  
     55   PUBLIC dyn_spg_ts        ! routine called in dynspg.F90 
    5456   PUBLIC dyn_spg_ts_alloc  !    "      "     "    " 
    5557   PUBLIC dyn_spg_ts_init   !    "      "     "    " 
     
    9799      ALLOCATE( wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj), STAT= ierr(2) ) 
    98100 
    99       IF( ln_dynvor_een ) ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , &  
     101      IF( ln_dynvor_een ) ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , & 
    100102                             &      ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(3) ) 
    101103 
     
    110112      !!---------------------------------------------------------------------- 
    111113      !! 
    112       !! ** Purpose :    
     114      !! ** Purpose : 
    113115      !!      -Compute the now trend due to the explicit time stepping 
    114       !!      of the quasi-linear barotropic system.  
     116      !!      of the quasi-linear barotropic system. 
    115117      !! 
    116       !! ** Method  :   
     118      !! ** Method  : 
    117119      !!      Barotropic variables are advanced from internal time steps 
    118120      !!      "n"   to "n+1" if ln_bt_fw=T 
    119       !!      or from  
     121      !!      or from 
    120122      !!      "n-1" to "n+1" if ln_bt_fw=F 
    121123      !!      thanks to a generalized forward-backward time stepping (see ref. below). 
     
    126128      !!      -Compute barotropic advective velocities at step "n" : un_adv, vn_adv 
    127129      !!      These are used to advect tracers and are compliant with discrete 
    128       !!      continuity equation taken at the baroclinic time steps. This  
     130      !!      continuity equation taken at the baroclinic time steps. This 
    129131      !!      ensures tracers conservation. 
    130132      !!      -Update 3d trend (ua, va) with barotropic component. 
    131133      !! 
    132       !! References : Shchepetkin, A.F. and J.C. McWilliams, 2005:  
    133       !!              The regional oceanic modeling system (ROMS):  
     134      !! References : Shchepetkin, A.F. and J.C. McWilliams, 2005: 
     135      !!              The regional oceanic modeling system (ROMS): 
    134136      !!              a split-explicit, free-surface, 
    135       !!              topography-following-coordinate oceanic model.  
    136       !!              Ocean Modelling, 9, 347-404.  
     137      !!              topography-following-coordinate oceanic model. 
     138      !!              Ocean Modelling, 9, 347-404. 
    137139      !!--------------------------------------------------------------------- 
    138140      ! 
    139141      INTEGER, INTENT(in)  ::   kt   ! ocean time-step index 
    140142      ! 
    141       LOGICAL  ::   ll_fw_start        ! if true, forward integration  
    142       LOGICAL  ::   ll_init             ! if true, special startup of 2d equations 
    143       INTEGER  ::   ji, jj, jk, jn        ! dummy loop indices 
    144       INTEGER  ::   ikbu, ikbv, noffset      ! local integers 
    145       REAL(wp) ::   zraur, z1_2dt_b, z2dt_bf    ! local scalars 
     143      LOGICAL  ::   ll_fw_start                 ! if true, forward integration 
     144      LOGICAL  ::   ll_init                         ! if true, special startup of 2d equations 
     145      LOGICAL  ::   ll_tmp1, ll_tmp2            ! local logical variables used in W/D 
     146      INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices 
     147      INTEGER  ::   ikbu, ikbv, noffset         ! local integers 
     148      REAL(wp) ::   zraur, z1_2dt_b, z2dt_bf    ! local scalars 
    146149      REAL(wp) ::   zx1, zy1, zx2, zy2         !   -      - 
    147       REAL(wp) ::   z1_12, z1_8, z1_4, z1_2    !   -      - 
     150      REAL(wp) ::   z1_12, z1_8, z1_4, z1_2       !   -      - 
    148151      REAL(wp) ::   zu_spg, zv_spg             !   -      - 
    149       REAL(wp) ::   zhura, zhvra               !   -      - 
    150       REAL(wp) ::   za0, za1, za2, za3           !   -      - 
     152      REAL(wp) ::   zhura, zhvra                     !   -      - 
     153      REAL(wp) ::   za0, za1, za2, za3                 !   -      - 
     154 
     155      REAL(wp) ::   ztmp                       ! temporary vaialbe 
    151156      ! 
    152157      REAL(wp), POINTER, DIMENSION(:,:) :: zun_e, zvn_e, zsshp2_e 
     
    156161      REAL(wp), POINTER, DIMENSION(:,:) :: zsshu_a, zsshv_a 
    157162      REAL(wp), POINTER, DIMENSION(:,:) :: zhf 
     163      REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy                 ! Wetting/Dying gravity filter coef. 
    158164      !!---------------------------------------------------------------------- 
    159165      ! 
     
    167173      CALL wrk_alloc( jpi, jpj, zsshu_a, zsshv_a                                   ) 
    168174      CALL wrk_alloc( jpi, jpj, zhf ) 
     175      IF(ln_wd) CALL wrk_alloc( jpi, jpj, zcpx, zcpy ) 
    169176      ! 
    170177      !                                         !* Local constant initialization 
    171       z1_12 = 1._wp / 12._wp  
    172       z1_8  = 0.125_wp                                    
     178      z1_12 = 1._wp / 12._wp 
     179      z1_8  = 0.125_wp 
    173180      z1_4  = 0.25_wp 
    174       z1_2  = 0.5_wp      
     181      z1_2  = 0.5_wp 
    175182      zraur = 1._wp / rau0 
    176183      ! 
    177       IF( kt == nit000 .AND. neuler == 0 ) THEN    ! reciprocal of baroclinic time step  
     184      IF( kt == nit000 .AND. neuler == 0 ) THEN         ! reciprocal of baroclinic time step 
    178185        z2dt_bf = rdt 
    179186      ELSE 
    180187        z2dt_bf = 2.0_wp * rdt 
    181188      ENDIF 
    182       z1_2dt_b = 1.0_wp / z2dt_bf  
    183       ! 
    184       ll_init = ln_bt_av                           ! if no time averaging, then no specific restart  
     189      z1_2dt_b = 1.0_wp / z2dt_bf 
     190      ! 
     191      ll_init = ln_bt_av                                ! if no time averaging, then no specific restart 
    185192      ll_fw_start = .FALSE. 
    186193      ! 
    187                                                        ! time offset in steps for bdy data update 
     194                                                        ! time offset in steps for bdy data update 
    188195      IF (.NOT.ln_bt_fw) THEN ; noffset=-2*nn_baro ; ELSE ;  noffset = 0 ; ENDIF 
    189196      ! 
    190       IF( kt == nit000 ) THEN                !* initialisation 
     197      IF( kt == nit000 ) THEN                   !* initialisation 
    191198         ! 
    192199         IF(lwp) WRITE(numout,*) 
     
    245252            IF ( .not. ln_sco ) THEN 
    246253! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case 
    247 !     Set it to zero for the time being  
     254!     Set it to zero for the time being 
    248255!              IF( rn_hmin < 0._wp ) THEN    ;   jk = - INT( rn_hmin )                                      ! from a nb of level 
    249256!              ELSE                          ;   jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 )  ! from a depth 
     
    264271            END DO 
    265272            CALL lbc_lnk( zhf, 'F', 1._wp ) 
    266             ! JC: TBC. hf should be greater than 0  
     273            ! JC: TBC. hf should be greater than 0 
    267274            DO jj = 1, jpj 
    268275               DO ji = 1, jpi 
     
    274281      ENDIF 
    275282      ! 
    276       ! If forward start at previous time step, and centered integration,  
     283      ! If forward start at previous time step, and centered integration, 
    277284      ! then update averaging weights: 
    278285      IF ((.NOT.ln_bt_fw).AND.((neuler==0).AND.(kt==nit000+1))) THEN 
     
    280287         CALL ts_wgt(ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2) 
    281288      ENDIF 
    282                            
     289 
    283290      ! ----------------------------------------------------------------------------- 
    284291      !  Phase 1 : Coupling between general trend and barotropic estimates (1st step) 
    285292      ! ----------------------------------------------------------------------------- 
    286       !       
     293      ! 
    287294      ! 
    288295      !                                   !* e3*d/dt(Ua) (Vertically integrated) 
     
    293300      DO jk = 1, jpkm1 
    294301         zu_frc(:,:) = zu_frc(:,:) + fse3u_n(:,:,jk) * ua(:,:,jk) * umask(:,:,jk) 
    295          zv_frc(:,:) = zv_frc(:,:) + fse3v_n(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)          
     302         zv_frc(:,:) = zv_frc(:,:) + fse3v_n(:,:,jk) * va(:,:,jk) * vmask(:,:,jk) 
    296303      END DO 
    297304      ! 
     
    311318      !                                   !* barotropic Coriolis trends (vorticity scheme dependent) 
    312319      !                                   ! -------------------------------------------------------- 
    313       zwx(:,:) = un_b(:,:) * hu(:,:) * e2u(:,:)        ! now fluxes  
     320      zwx(:,:) = un_b(:,:) * hu(:,:) * e2u(:,:)        ! now fluxes 
    314321      zwy(:,:) = vn_b(:,:) * hv(:,:) * e1v(:,:) 
    315322      ! 
     
    353360         END DO 
    354361         ! 
    355       ENDIF  
     362      ENDIF 
    356363      ! 
    357364      !                                   !* Right-Hand-Side of the barotropic momentum equation 
    358365      !                                   ! ---------------------------------------------------- 
    359366      IF( lk_vvl ) THEN                         ! Variable volume : remove surface pressure gradient 
    360          DO jj = 2, jpjm1  
    361             DO ji = fs_2, fs_jpim1   ! vector opt. 
    362                zu_trd(ji,jj) = zu_trd(ji,jj) - grav * (  sshn(ji+1,jj  ) - sshn(ji  ,jj  )  ) / e1u(ji,jj) 
    363                zv_trd(ji,jj) = zv_trd(ji,jj) - grav * (  sshn(ji  ,jj+1) - sshn(ji  ,jj  )  ) / e2v(ji,jj) 
    364             END DO 
    365          END DO 
     367        IF(ln_wd) THEN                    ! Calculating and applying W/D gravity filters 
     368          DO jj = 2, jpjm1 
     369             DO ji = 2, jpim1 
     370                ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) 
     371                ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) + & 
     372                                                          & rn_wdmin1 + rn_wdmin2 
     373                IF(ll_tmp1) THEN 
     374                  zcpx(ji,jj) = 1.0_wp 
     375                ELSE IF(ll_tmp2) THEN 
     376                   ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 
     377                  zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /& 
     378                              &     (sshn(ji+1,jj) - sshn(ji,jj))) 
     379                ELSE 
     380                   zcpx(ji,jj) = 0._wp 
     381                END IF 
     382 
     383                ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) 
     384                ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) + & 
     385                                                          & rn_wdmin1 + rn_wdmin2 
     386                IF(ll_tmp1) THEN 
     387                   zcpy(ji,jj) = 1.0_wp 
     388                ELSE IF(ll_tmp2) THEN 
     389                   ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 
     390                  zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) + bathy(ji,jj)) /& 
     391                              &     (sshn(ji,jj+1) - sshn(ji,jj))) 
     392                ELSE 
     393                  zcpy(ji,jj) = 0._wp 
     394                END IF 
     395             END DO 
     396           END DO 
     397           CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp ) 
     398        ENDIF 
     399 
     400        IF(ln_wd) THEN                    ! Calculating and applying W/D gravity filters 
     401           DO jj = 2, jpjm1 
     402              DO ji = 2, jpim1 
     403                 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj  ) - sshn(ji  ,jj ) ) / & 
     404                               &                 e1u(ji,jj) * zcpx(ji,jj) 
     405                 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji  ,jj+1) - sshn(ji  ,jj ) ) / & 
     406                               &                 e2v(ji,jj) * zcpy(ji,jj) 
     407              END DO 
     408           END DO 
     409         ELSE 
     410           DO jj = 2, jpjm1 
     411              DO ji = fs_2, fs_jpim1   ! vector opt. 
     412                 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * (  sshn(ji+1,jj  ) - sshn(ji  ,jj  )  ) / e1u(ji,jj) 
     413                 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * (  sshn(ji  ,jj+1) - sshn(ji  ,jj  )  ) / e2v(ji,jj) 
     414              END DO 
     415           END DO 
     416        END IF 
     417 
    366418      ENDIF 
    367419 
     
    371423             zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * vmask(ji,jj,1) 
    372424          END DO 
    373       END DO  
    374       ! 
    375       !                 ! Add bottom stress contribution from baroclinic velocities:       
     425      END DO 
     426      ! 
     427      !                                         ! Add bottom stress contribution from baroclinic velocities: 
    376428      IF (ln_bt_fw) THEN 
    377          DO jj = 2, jpjm1                           
     429         DO jj = 2, jpjm1 
    378430            DO ji = fs_2, fs_jpim1   ! vector opt. 
    379                ikbu = mbku(ji,jj)        
    380                ikbv = mbkv(ji,jj)     
     431               ikbu = mbku(ji,jj) 
     432               ikbv = mbkv(ji,jj) 
    381433               zwx(ji,jj) = un(ji,jj,ikbu) - un_b(ji,jj) ! NOW bottom baroclinic velocities 
    382434               zwy(ji,jj) = vn(ji,jj,ikbv) - vn_b(ji,jj) 
     
    386438         DO jj = 2, jpjm1 
    387439            DO ji = fs_2, fs_jpim1   ! vector opt. 
    388                ikbu = mbku(ji,jj)        
    389                ikbv = mbkv(ji,jj)     
     440               ikbu = mbku(ji,jj) 
     441               ikbv = mbkv(ji,jj) 
    390442               zwx(ji,jj) = ub(ji,jj,ikbu) - ub_b(ji,jj) ! BEFORE bottom baroclinic velocities 
    391443               zwy(ji,jj) = vb(ji,jj,ikbv) - vb_b(ji,jj) 
     
    397449      zu_frc(:,:) = zu_frc(:,:) + hur(:,:) * bfrua(:,:) * zwx(:,:) 
    398450      zv_frc(:,:) = zv_frc(:,:) + hvr(:,:) * bfrva(:,:) * zwy(:,:) 
    399       !        
     451      ! 
    400452      IF (ln_bt_fw) THEN                        ! Add wind forcing 
    401453         zu_frc(:,:) =  zu_frc(:,:) + zraur * utau(:,:) * hur(:,:) 
     
    404456         zu_frc(:,:) =  zu_frc(:,:) + zraur * z1_2 * ( utau_b(:,:) + utau(:,:) ) * hur(:,:) 
    405457         zv_frc(:,:) =  zv_frc(:,:) + zraur * z1_2 * ( vtau_b(:,:) + vtau(:,:) ) * hvr(:,:) 
    406       ENDIF   
     458      ENDIF 
    407459      ! 
    408460      IF ( ln_apr_dyn ) THEN                    ! Add atm pressure forcing 
    409461         IF (ln_bt_fw) THEN 
    410             DO jj = 2, jpjm1               
     462            DO jj = 2, jpjm1 
    411463               DO ji = fs_2, fs_jpim1   ! vector opt. 
    412464                  zu_spg =  grav * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj) ) /e1u(ji,jj) 
     
    417469            END DO 
    418470         ELSE 
    419             DO jj = 2, jpjm1               
     471            DO jj = 2, jpjm1 
    420472               DO ji = fs_2, fs_jpim1   ! vector opt. 
    421473                  zu_spg =  grav * z1_2 * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj)    & 
     
    427479               END DO 
    428480            END DO 
    429          ENDIF  
     481         ENDIF 
    430482      ENDIF 
    431483      !                                   !* Right-Hand-Side of the barotropic ssh equation 
     
    450502      ! 
    451503      ! ----------------------------------------------------------------------- 
    452       !  Phase 2 : Integration of the barotropic equations  
     504      !  Phase 2 : Integration of the barotropic equations 
    453505      ! ----------------------------------------------------------------------- 
    454506      ! 
    455507      !                                             ! ==================== ! 
    456508      !                                             !    Initialisations   ! 
    457       !                                             ! ==================== !   
    458       ! Initialize barotropic variables:       
     509      !                                             ! ==================== ! 
     510      ! Initialize barotropic variables: 
    459511      IF( ll_init )THEN 
    460512         sshbb_e(:,:) = 0._wp 
     
    465517         vb_e   (:,:) = 0._wp 
    466518      ENDIF 
    467       ! 
    468       IF (ln_bt_fw) THEN                  ! FORWARD integration: start from NOW fields                     
    469          sshn_e(:,:) = sshn (:,:)             
    470          zun_e (:,:) = un_b (:,:)             
     519 
     520      IF(ln_wd) THEN      !preserve the positivity of water depth 
     521                          !ssh[b,n,a] should have already been processed for this 
     522         sshbb_e(:,:) = MAX(sshbb_e(:,:), rn_wdmin1 - bathy(:,:)) 
     523         sshb_e(:,:)  = MAX(sshb_e(:,:) , rn_wdmin1 - bathy(:,:)) 
     524      ENDIF 
     525      ! 
     526      IF (ln_bt_fw) THEN                  ! FORWARD integration: start from NOW fields 
     527         sshn_e(:,:) = sshn (:,:) 
     528         zun_e (:,:) = un_b (:,:) 
    471529         zvn_e (:,:) = vn_b (:,:) 
    472530         ! 
    473          hu_e  (:,:) = hu   (:,:)        
    474          hv_e  (:,:) = hv   (:,:)  
    475          hur_e (:,:) = hur  (:,:)     
     531         hu_e  (:,:) = hu   (:,:) 
     532         hv_e  (:,:) = hv   (:,:) 
     533         hur_e (:,:) = hur  (:,:) 
    476534         hvr_e (:,:) = hvr  (:,:) 
    477535      ELSE                                ! CENTRED integration: start from BEFORE fields 
    478536         sshn_e(:,:) = sshb (:,:) 
    479          zun_e (:,:) = ub_b (:,:)          
     537         zun_e (:,:) = ub_b (:,:) 
    480538         zvn_e (:,:) = vb_b (:,:) 
    481539         ! 
    482          hu_e  (:,:) = hu_b (:,:)        
    483          hv_e  (:,:) = hv_b (:,:)  
    484          hur_e (:,:) = hur_b(:,:)     
     540         hu_e  (:,:) = hu_b (:,:) 
     541         hv_e  (:,:) = hv_b (:,:) 
     542         hur_e (:,:) = hur_b(:,:) 
    485543         hvr_e (:,:) = hvr_b(:,:) 
    486544      ENDIF 
     
    489547      ! 
    490548      ! Initialize sums: 
    491       ua_b  (:,:) = 0._wp       ! After barotropic velocities (or transport if flux form)           
     549      ua_b  (:,:) = 0._wp       ! After barotropic velocities (or transport if flux form) 
    492550      va_b  (:,:) = 0._wp 
    493551      ssha  (:,:) = 0._wp       ! Sum for after averaged sea level 
     
    506564         ! 
    507565         ! Set extrapolation coefficients for predictor step: 
    508          IF ((jn<3).AND.ll_init) THEN      ! Forward            
    509            za1 = 1._wp                                           
    510            za2 = 0._wp                         
    511            za3 = 0._wp                         
    512          ELSE                              ! AB3-AM4 Coefficients: bet=0.281105  
     566         IF ((jn<3).AND.ll_init) THEN      ! Forward 
     567           za1 = 1._wp 
     568           za2 = 0._wp 
     569           za3 = 0._wp 
     570         ELSE                              ! AB3-AM4 Coefficients: bet=0.281105 
    513571           za1 =  1.781105_wp              ! za1 =   3/2 +   bet 
    514572           za2 = -1.06221_wp               ! za2 = -(1/2 + 2*bet) 
     
    524582            ! Extrapolate Sea Level at step jit+0.5: 
    525583            zsshp2_e(:,:) = za1 * sshn_e(:,:)  + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:) 
     584 
     585 
    526586            ! 
    527587            DO jj = 2, jpjm1                                    ! Sea Surface Height at u- & v-points 
     
    535595               END DO 
    536596            END DO 
     597 
    537598            CALL lbc_lnk( zwx, 'U', 1._wp )    ;   CALL lbc_lnk( zwy, 'V', 1._wp ) 
    538599            ! 
    539600            zhup2_e (:,:) = hu_0(:,:) + zwx(:,:)                ! Ocean depth at U- and V-points 
    540601            zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:) 
     602            IF(ln_wd) THEN 
     603              zhup2_e(:,:) = MAX(zhup2_e (:,:), rn_wdmin1) 
     604              zhvp2_e(:,:) = MAX(zhvp2_e (:,:), rn_wdmin1) 
     605            END IF 
    541606         ELSE 
    542607            zhup2_e (:,:) = hu(:,:) 
     
    552617         ! 
    553618#if defined key_agrif 
    554          ! Set fluxes during predictor step to ensure  
     619         ! Set fluxes during predictor step to ensure 
    555620         ! volume conservation 
    556621         IF( (.NOT.Agrif_Root()).AND.ln_bt_fw ) THEN 
     
    577642         ENDIF 
    578643#endif 
     644         IF(ln_wd) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt) 
    579645         ! 
    580646         ! Sum over sub-time-steps to compute advective velocities 
     
    584650         ! 
    585651         ! Set next sea level: 
    586          DO jj = 2, jpjm1                                  
     652         DO jj = 2, jpjm1 
    587653            DO ji = fs_2, fs_jpim1   ! vector opt. 
    588654               zhdiv(ji,jj) = (   zwx(ji,jj) - zwx(ji-1,jj)   & 
     
    590656            END DO 
    591657         END DO 
     658         !IF(ln_wd) THEN 
     659 
    592660         ssha_e(:,:) = (  sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) )  ) * tmask(:,:,1) 
     661         IF(ln_wd) ssha_e(:,:) = MAX(ssha_e(:,:), rn_wdmin1 - bathy(:,:))  
    593662         CALL lbc_lnk( ssha_e, 'T',  1._wp ) 
    594663 
     
    600669         IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jn ) 
    601670#endif 
    602          !   
     671         ! 
    603672         ! Sea Surface Height at u-,v-points (vvl case only) 
    604          IF ( lk_vvl ) THEN                                 
     673         IF ( lk_vvl ) THEN 
    605674            DO jj = 2, jpjm1 
    606675               DO ji = 2, jpim1      ! NO Vector Opt. 
     
    613682               END DO 
    614683            END DO 
     684 
    615685            CALL lbc_lnk( zsshu_a, 'U', 1._wp )   ;   CALL lbc_lnk( zsshv_a, 'V', 1._wp ) 
    616          ENDIF    
    617          !                                  
     686         ENDIF 
     687         ! 
    618688         ! Half-step back interpolation of SSH for surface pressure computation: 
    619689         !---------------------------------------------------------------------- 
    620690         IF ((jn==1).AND.ll_init) THEN 
    621691           za0=1._wp                        ! Forward-backward 
    622            za1=0._wp                            
     692           za1=0._wp 
    623693           za2=0._wp 
    624694           za3=0._wp 
     
    627697           za1=-0.1666666666666_wp          ! za1 = gam 
    628698           za2= 0.0833333333333_wp          ! za2 = eps 
    629            za3= 0._wp               
    630          ELSE                               ! AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880  
    631            za0=0.614_wp                     ! za0 = 1/2 +   gam + 2*eps     
    632            za1=0.285_wp                     ! za1 = 1/2 - 2*gam - 3*eps  
     699           za3= 0._wp 
     700         ELSE                               ! AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880 
     701           za0=0.614_wp                     ! za0 = 1/2 +   gam + 2*eps 
     702           za1=0.285_wp                     ! za1 = 1/2 - 2*gam - 3*eps 
    633703           za2=0.088_wp                     ! za2 = gam 
    634704           za3=0.013_wp                     ! za3 = eps 
     
    640710         ! 
    641711         ! Compute associated depths at U and V points: 
    642          IF ( lk_vvl.AND.(.NOT.ln_dynadv_vec) ) THEN       
    643             !                                         
    644             DO jj = 2, jpjm1                             
     712         IF ( lk_vvl.AND.(.NOT.ln_dynadv_vec) ) THEN 
     713            ! 
     714            DO jj = 2, jpjm1 
    645715               DO ji = 2, jpim1 
    646716                  zx1 = z1_2 * umask(ji  ,jj,1) *  r1_e12u(ji  ,jj)    & 
     
    650720                     &       * ( e12t(ji ,jj  ) * zsshp2_e(ji  ,jj  )  & 
    651721                     &       +   e12t(ji ,jj+1) * zsshp2_e(ji  ,jj+1) ) 
    652                   zhust_e(ji,jj) = hu_0(ji,jj) + zx1  
     722                  zhust_e(ji,jj) = hu_0(ji,jj) + zx1 
    653723                  zhvst_e(ji,jj) = hv_0(ji,jj) + zy1 
    654724               END DO 
    655725            END DO 
     726 
     727            IF(ln_wd) THEN 
     728              zhust_e(:,:) = MAX(zhust_e (:,:), rn_wdmin1 ) 
     729              zhvst_e(:,:) = MAX(zhvst_e (:,:), rn_wdmin1 ) 
     730            END IF 
     731            !shall we call lbc_lnk for zhu(v)st_e() here? 
     732 
    656733         ENDIF 
    657734         ! 
     
    692769                  zu_trd(ji,jj) = + z1_12 / e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) & 
    693770                     &                                    + ftnw(ji+1,jj) * zwy(ji+1,jj  ) & 
    694                      &                                    + ftse(ji,jj  ) * zwy(ji  ,jj-1) &  
     771                     &                                    + ftse(ji,jj  ) * zwy(ji  ,jj-1) & 
    695772                     &                                    + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) 
    696                   zv_trd(ji,jj) = - z1_12 / e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) &  
     773                  zv_trd(ji,jj) = - z1_12 / e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) & 
    697774                     &                                    + ftse(ji,jj+1) * zwx(ji  ,jj+1) & 
    698                      &                                    + ftnw(ji,jj  ) * zwx(ji-1,jj  ) &  
     775                     &                                    + ftnw(ji,jj  ) * zwx(ji-1,jj  ) & 
    699776                     &                                    + ftne(ji,jj  ) * zwx(ji  ,jj  ) ) 
    700777               END DO 
    701778            END DO 
    702             !  
     779            ! 
    703780         ENDIF 
    704781         ! 
     
    720797         ! 
    721798         ! Surface pressure trend: 
    722          DO jj = 2, jpjm1 
    723             DO ji = fs_2, fs_jpim1   ! vector opt. 
    724                ! Add surface pressure gradient 
    725                zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj) 
    726                zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj) 
    727                zwx(ji,jj) = zu_spg 
    728                zwy(ji,jj) = zv_spg 
    729             END DO 
    730          END DO 
     799         IF(ln_wd) THEN                    ! Calculating and applying W/D gravity filters 
     800           DO jj = 2, jpjm1 
     801              DO ji = 2, jpim1 
     802                 ll_tmp1 = MIN(zsshp2_e(ji,jj), zsshp2_e(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) 
     803                 ll_tmp2 = MAX(zsshp2_e(ji,jj), zsshp2_e(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) + & 
     804                                                           & rn_wdmin1 + rn_wdmin2 
     805                 IF(ll_tmp1) THEN 
     806                   zcpx(ji,jj) = 1.0_wp 
     807                 ELSE IF(ll_tmp2) THEN 
     808                    ! no worries about zsshp2_e(ji+1,jj)-zsshp2_e(ji,jj) = 0, it won't happen ! here 
     809                   zcpx(ji,jj) = ABS((zsshp2_e(ji+1,jj) + bathy(ji+1,jj) - zsshp2_e(ji,jj) - bathy(ji,jj)) /& 
     810                               & (zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj))) 
     811                 ELSE 
     812                    zcpx(ji,jj) = 0._wp 
     813                 END IF 
     814 
     815                 ll_tmp1 = MIN(zsshp2_e(ji,jj), zsshp2_e(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) 
     816                 ll_tmp2 = MAX(zsshp2_e(ji,jj), zsshp2_e(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) + & 
     817                                                           & rn_wdmin1 + rn_wdmin2 
     818                 IF(ll_tmp1) THEN 
     819                    zcpy(ji,jj) = 1.0_wp 
     820                 ELSE IF(ll_tmp2) THEN 
     821                    ! no worries about zsshp2_e(ji,jj+1)-zsshp2_e(ji,jj) = 0, it won't happen ! here 
     822                   zcpy(ji,jj) = ABS((zsshp2_e(ji,jj+1) + bathy(ji,jj+1) - zsshp2_e(ji,jj) - bathy(ji,jj)) /& 
     823                               & (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj))) 
     824                 ELSE 
     825                   zcpy(ji,jj) = 0._wp 
     826                 END IF 
     827              END DO 
     828            END DO 
     829            CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp ) 
     830         ENDIF 
     831 
     832         IF(ln_wd) THEN 
     833           DO jj = 2, jpjm1 
     834              DO ji = 2, jpim1  
     835                 ! Add surface pressure gradient 
     836                 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj) 
     837                 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj) 
     838                 zwx(ji,jj) = zu_spg * zcpx(ji,jj) 
     839                 zwy(ji,jj) = zv_spg * zcpy(ji,jj) 
     840              END DO 
     841           END DO 
     842         ELSE 
     843           DO jj = 2, jpjm1 
     844              DO ji = fs_2, fs_jpim1   ! vector opt. 
     845                 ! Add surface pressure gradient 
     846                 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj) 
     847                 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj) 
     848                 zwx(ji,jj) = zu_spg 
     849                 zwy(ji,jj) = zv_spg 
     850              END DO 
     851           END DO 
     852         END IF 
     853 
    731854         ! 
    732855         ! Set next velocities: 
    733          IF( ln_dynadv_vec .OR. (.NOT. lk_vvl) ) THEN    ! Vector form 
     856         IF( ln_dynadv_vec .OR. (.NOT. lk_vvl) ) THEN   ! Vector form 
    734857            DO jj = 2, jpjm1 
    735858               DO ji = fs_2, fs_jpim1   ! vector opt. 
    736                   ua_e(ji,jj) = (                                zun_e(ji,jj)   &  
     859                  ua_e(ji,jj) = (                                zun_e(ji,jj)   & 
    737860                            &     + rdtbt * (                      zwx(ji,jj)   & 
    738861                            &                                 + zu_trd(ji,jj)   & 
    739                             &                                 + zu_frc(ji,jj) ) &  
     862                            &                                 + zu_frc(ji,jj) ) & 
    740863                            &   ) * umask(ji,jj,1) 
    741864 
     
    748871            END DO 
    749872 
    750          ELSE                 ! Flux form 
     873         ELSE                                           ! Flux form 
    751874            DO jj = 2, jpjm1 
    752875               DO ji = fs_2, fs_jpim1   ! vector opt. 
    753876 
    754                   zhura = umask(ji,jj,1)/(hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - umask(ji,jj,1)) 
    755                   zhvra = vmask(ji,jj,1)/(hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - vmask(ji,jj,1)) 
    756  
    757                   ua_e(ji,jj) = (                hu_e(ji,jj)  *  zun_e(ji,jj)   &  
    758                             &     + rdtbt * ( zhust_e(ji,jj)  *    zwx(ji,jj)   &  
     877                  IF(ln_wd) THEN 
     878                    zhura = MAX(hu_0(ji,jj) + zsshu_a(ji,jj), rn_wdmin1) 
     879                    zhvra = MAX(hv_0(ji,jj) + zsshv_a(ji,jj), rn_wdmin1) 
     880                  ELSE 
     881                    zhura = hu_0(ji,jj) + zsshu_a(ji,jj) 
     882                    zhvra = hv_0(ji,jj) + zsshv_a(ji,jj) 
     883                  END IF 
     884 
     885                  zhura = umask(ji,jj,1)/(zhura + 1._wp - umask(ji,jj,1)) 
     886                  zhvra = vmask(ji,jj,1)/(zhvra + 1._wp - vmask(ji,jj,1)) 
     887 
     888                  ua_e(ji,jj) = (                hu_e(ji,jj)  *  zun_e(ji,jj)   & 
     889                            &     + rdtbt * ( zhust_e(ji,jj)  *    zwx(ji,jj)   & 
    759890                            &               + zhup2_e(ji,jj)  * zu_trd(ji,jj)   & 
    760891                            &               +      hu(ji,jj)  * zu_frc(ji,jj) ) & 
     
    771902         ! 
    772903         IF( lk_vvl ) THEN                             !* Update ocean depth (variable volume case only) 
    773             !                                          !  ----------------------------------------------         
    774             hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 
    775             hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 
     904            !                                          !  ---------------------------------------------- 
     905            IF(ln_wd) THEN 
     906              hu_e (:,:) = MAX(hu_0(:,:) + zsshu_a(:,:), rn_wdmin1) 
     907              hv_e (:,:) = MAX(hv_0(:,:) + zsshv_a(:,:), rn_wdmin1) 
     908            ELSE 
     909              hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 
     910              hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 
     911            END IF 
     912 
    776913            hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1._wp - umask(:,:,1) ) 
    777914            hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1._wp - vmask(:,:,1) ) 
     
    781918         !                                                 !  ----------------------- 
    782919         ! 
    783          CALL lbc_lnk( ua_e  , 'U', -1._wp )               ! local domain boundaries  
     920         CALL lbc_lnk( ua_e  , 'U', -1._wp )               ! local domain boundaries 
    784921         CALL lbc_lnk( va_e  , 'V', -1._wp ) 
    785922 
    786 #if defined key_bdy   
     923#if defined key_bdy 
    787924                                                           ! open boundaries 
    788925         IF( lk_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, zun_e, zvn_e, hur_e, hvr_e, ssha_e ) 
    789926#endif 
    790 #if defined key_agrif                                                            
     927#if defined key_agrif 
    791928         IF( .NOT.Agrif_Root() )  CALL agrif_dyn_ts( jn )  ! Agrif 
    792929#endif 
     
    807944         !                                             !* Sum over whole bt loop 
    808945         !                                             !  ---------------------- 
    809          za1 = wgtbtp1(jn)                                     
     946         za1 = wgtbtp1(jn) 
    810947         IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN    ! Sum velocities 
    811             ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:)  
    812             va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:)  
    813          ELSE                                ! Sum transports 
     948            ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) 
     949            va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) 
     950         ELSE                                              ! Sum transports 
    814951            ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) * hu_e (:,:) 
    815952            va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) * hv_e (:,:) 
    816953         ENDIF 
    817          !                                   ! Sum sea level 
     954         !                                                 ! Sum sea level 
    818955         ssha(:,:) = ssha(:,:) + za1 * ssha_e(:,:) 
    819956         !                                                 ! ==================== ! 
     
    841978      ! 
    842979      ! Set advection velocity correction: 
    843       IF (((kt==nit000).AND.(neuler==0)).OR.(.NOT.ln_bt_fw)) THEN      
     980      IF (((kt==nit000).AND.(neuler==0)).OR.(.NOT.ln_bt_fw)) THEN 
    844981         un_adv(:,:) = zu_sum(:,:)*hur(:,:) 
    845982         vn_adv(:,:) = zv_sum(:,:)*hvr(:,:) 
     
    8941031      ! 
    8951032      ! 
    896 #endif       
     1033#endif 
    8971034      ! 
    8981035      !                                   !* write time-spliting arrays in the restart 
     
    9051042      CALL wrk_dealloc( jpi, jpj, zsshu_a, zsshv_a                                   ) 
    9061043      CALL wrk_dealloc( jpi, jpj, zhf ) 
     1044      IF(ln_wd) CALL wrk_alloc( jpi, jpj, zcpx, zcpy ) 
    9071045      ! 
    9081046      IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_ts') 
     
    9181056      LOGICAL, INTENT(in) ::   ll_av      ! temporal averaging=.true. 
    9191057      LOGICAL, INTENT(in) ::   ll_fw      ! forward time splitting =.true. 
    920       INTEGER, INTENT(inout) :: jpit      ! cycle length     
     1058      INTEGER, INTENT(inout) :: jpit      ! cycle length 
    9211059      REAL(wp), DIMENSION(3*nn_baro), INTENT(inout) ::   zwgt1, & ! Primary weights 
    9221060                                                         zwgt2    ! Secondary weights 
    923        
     1061 
    9241062      INTEGER ::  jic, jn, ji                      ! temporary integers 
    9251063      REAL(wp) :: za1, za2 
     
    9301068 
    9311069      ! Set time index when averaged value is requested 
    932       IF (ll_fw) THEN  
     1070      IF (ll_fw) THEN 
    9331071         jic = nn_baro 
    9341072      ELSE 
     
    9381076      ! Set primary weights: 
    9391077      IF (ll_av) THEN 
    940            ! Define simple boxcar window for primary weights  
    941            ! (width = nn_baro, centered around jic)      
     1078           ! Define simple boxcar window for primary weights 
     1079           ! (width = nn_baro, centered around jic) 
    9421080         SELECT CASE ( nn_bt_flt ) 
    9431081              CASE( 0 )  ! No averaging 
     
    9471085              CASE( 1 )  ! Boxcar, width = nn_baro 
    9481086                 DO jn = 1, 3*nn_baro 
    949                     za1 = ABS(float(jn-jic))/float(nn_baro)  
     1087                    za1 = ABS(float(jn-jic))/float(nn_baro) 
    9501088                    IF (za1 < 0.5_wp) THEN 
    9511089                      zwgt1(jn) = 1._wp 
     
    9561094              CASE( 2 )  ! Boxcar, width = 2 * nn_baro 
    9571095                 DO jn = 1, 3*nn_baro 
    958                     za1 = ABS(float(jn-jic))/float(nn_baro)  
     1096                    za1 = ABS(float(jn-jic))/float(nn_baro) 
    9591097                    IF (za1 < 1._wp) THEN 
    9601098                      zwgt1(jn) = 1._wp 
     
    9691107         jpit = jic 
    9701108      ENDIF 
    971      
     1109 
    9721110      ! Set secondary weights 
    9731111      DO jn = 1, jpit 
     
    9991137      ! 
    10001138      IF( TRIM(cdrw) == 'READ' ) THEN 
    1001          CALL iom_get( numror, jpdom_autoglo, 'ub2_b'  , ub2_b  (:,:) )    
    1002          CALL iom_get( numror, jpdom_autoglo, 'vb2_b'  , vb2_b  (:,:) )  
     1139         CALL iom_get( numror, jpdom_autoglo, 'ub2_b'  , ub2_b  (:,:) ) 
     1140         CALL iom_get( numror, jpdom_autoglo, 'vb2_b'  , vb2_b  (:,:) ) 
    10031141         IF( .NOT.ln_bt_av ) THEN 
    1004             CALL iom_get( numror, jpdom_autoglo, 'sshbb_e'  , sshbb_e(:,:) )    
    1005             CALL iom_get( numror, jpdom_autoglo, 'ubb_e'    ,   ubb_e(:,:) )    
     1142            CALL iom_get( numror, jpdom_autoglo, 'sshbb_e'  , sshbb_e(:,:) ) 
     1143            CALL iom_get( numror, jpdom_autoglo, 'ubb_e'    ,   ubb_e(:,:) ) 
    10061144            CALL iom_get( numror, jpdom_autoglo, 'vbb_e'    ,   vbb_e(:,:) ) 
    1007             CALL iom_get( numror, jpdom_autoglo, 'sshb_e'   ,  sshb_e(:,:) )  
    1008             CALL iom_get( numror, jpdom_autoglo, 'ub_e'     ,    ub_e(:,:) )    
     1145            CALL iom_get( numror, jpdom_autoglo, 'sshb_e'   ,  sshb_e(:,:) ) 
     1146            CALL iom_get( numror, jpdom_autoglo, 'ub_e'     ,    ub_e(:,:) ) 
    10091147            CALL iom_get( numror, jpdom_autoglo, 'vb_e'     ,    vb_e(:,:) ) 
    10101148         ENDIF 
     
    10121150         ! Read time integrated fluxes 
    10131151         IF ( .NOT.Agrif_Root() ) THEN 
    1014             CALL iom_get( numror, jpdom_autoglo, 'ub2_i_b'  , ub2_i_b(:,:) )    
     1152            CALL iom_get( numror, jpdom_autoglo, 'ub2_i_b'  , ub2_i_b(:,:) ) 
    10151153            CALL iom_get( numror, jpdom_autoglo, 'vb2_i_b'  , vb2_i_b(:,:) ) 
    10161154         ENDIF 
     
    10221160         ! 
    10231161         IF (.NOT.ln_bt_av) THEN 
    1024             CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e'  , sshbb_e(:,:) )  
     1162            CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e'  , sshbb_e(:,:) ) 
    10251163            CALL iom_rstput( kt, nitrst, numrow, 'ubb_e'    ,   ubb_e(:,:) ) 
    10261164            CALL iom_rstput( kt, nitrst, numrow, 'vbb_e'    ,   vbb_e(:,:) ) 
     
    10701208      CALL wrk_alloc( jpi, jpj, zcu ) 
    10711209      ! 
    1072       IF (lk_vvl) THEN  
     1210      IF (lk_vvl) THEN 
    10731211         DO jj = 1, jpj 
    10741212            DO ji =1, jpi 
     
    10931231      ! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax 
    10941232      IF (ln_bt_nn_auto) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax) 
    1095        
     1233 
    10961234      rdtbt = rdt / FLOAT(nn_baro) 
    10971235      zcmax = zcmax * rdtbt 
    1098                      ! Print results 
     1236                                                        ! Print results 
    10991237      IF(lwp) WRITE(numout,*) 
    11001238      IF(lwp) WRITE(numout,*) 'dyn_spg_ts : split-explicit free surface' 
     
    11291267           CASE( 0 )  ;   IF(lwp) WRITE(numout,*) '           Dirac' 
    11301268           CASE( 1 )  ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = nn_baro' 
    1131            CASE( 2 )  ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = 2*nn_baro'  
     1269           CASE( 2 )  ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = 2*nn_baro' 
    11321270           CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1,2' ) 
    11331271      END SELECT 
     
    11421280      ENDIF 
    11431281      IF ( zcmax>0.9_wp ) THEN 
    1144          CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_baro !' )           
     1282         CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_baro !' ) 
    11451283      ENDIF 
    11461284      ! 
     
    11651303      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag 
    11661304      WRITE(*,*) 'ts_rst    : You should not have seen this print! error?', kt, cdrw 
    1167    END SUBROUTINE ts_rst   
     1305   END SUBROUTINE ts_rst 
    11681306   SUBROUTINE dyn_spg_ts_init( kt )       ! Empty routine 
    11691307      INTEGER         , INTENT(in) ::   kt         ! ocean time-step 
     
    11711309   END SUBROUTINE dyn_spg_ts_init 
    11721310#endif 
    1173     
     1311 
    11741312   !!====================================================================== 
    11751313END MODULE dynspg_ts 
  • branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90

    r4486 r5014  
    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.?  !  2014-09  (H. Liu) add wetting and drying   
    1112   !!---------------------------------------------------------------------- 
    1213 
     
    4243   USE wrk_nemo        ! Memory Allocation 
    4344   USE timing          ! Timing 
     45   USE wadlmt          ! Wetting/Drying flux limting 
    4446 
    4547   IMPLICIT NONE 
     
    9496      ENDIF 
    9597      ! 
    96       CALL div_cur( kt )                              ! Horizontal divergence & Relative vorticity 
    97       ! 
    9898      z2dt = 2._wp * rdt                              ! set time step size (Euler/Leapfrog) 
    9999      IF( neuler == 0 .AND. kt == nit000 )   z2dt = rdt 
     100      ! 
     101      z1_rau0 = 0.5_wp * r1_rau0 
     102      ! 
     103      IF(ln_wd) CALL wad_lmt(sshb, z1_rau0 * (emp_b(:,:) + emp(:,:)), z2dt) 
     104 
     105      CALL div_cur( kt )                              ! Horizontal divergence & Relative vorticity 
     106      ! 
    100107 
    101108      !                                           !------------------------------! 
     
    110117      ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations. 
    111118      !  
    112       z1_rau0 = 0.5_wp * r1_rau0 
    113119      ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * tmask(:,:,1) 
    114120 
  • branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90

    r4723 r5014  
    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.?! 2014-09  (H. Liu) Add Wetting and Drying  
    3233   !!---------------------------------------------------------------------- 
    3334 
     
    233234      NAMELIST/namcfg/ cp_cfg, cp_cfz, jp_cfg, jpidta, jpjdta, jpkdta, jpiglo, jpjglo, & 
    234235         &             jpizoom, jpjzoom, jperio 
     236      NAMELIST/namwad/ ln_wd, rn_wdmin1, rn_wdmin2, rn_wdld, nn_wdit 
    235237      !!---------------------------------------------------------------------- 
    236238      ! 
     
    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', lwp ) 
     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', lwp ) 
     269 
    259270 
    260271! Force values for AGRIF zoom (cf. agrif_user.F90) 
  • branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/step_oce.F90

    r4328 r5014  
    116116#endif 
    117117#if defined key_top 
    118    USE trcstp           ! passive tracer time-stepping      (trc_stp routine) 
     118   USE trcstp           ! passive tracer time-stepping     (trc_stp routine) 
    119119#endif 
     120   USE wadlmt           ! wetting/drying limiter           (wad_lmt routine) 
    120121   !!---------------------------------------------------------------------- 
    121122   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
Note: See TracChangeset for help on using the changeset viewer.