Changeset 5066 for branches/UKMO


Ignore:
Timestamp:
2015-02-06T17:02:20+01:00 (6 years ago)
Author:
rfurner
Message:

added current state of wetting and drying code to test…note it does not work

Location:
branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC
Files:
1 added
9 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90

    r4488 r5066  
    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 
     
    262263   INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  :: nicoa, njcoa       !: ??? 
    263264#endif 
     265 
     266   !!---------------------------------------------------------------------- 
     267   !! critical depths,limiters,and masks for  Wetting and Drying 
     268   !! --------------------------------------------------------------------- 
     269 
     270   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wdmask             !: u- and v- limiter 
     271 
     272   LOGICAL,  PUBLIC, SAVE ::   ln_wd       !: key to turn on/off wetting/drying (T: on, F: off) 
     273   REAL(wp), PUBLIC, SAVE ::   rn_wdmin1   !: minimum water depth on dried cells 
     274   REAL(wp), PUBLIC, SAVE ::   rn_wdmin2   !: tolerrance of minimum water depth on dried cells 
     275   REAL(wp), PUBLIC, SAVE ::   rn_wdld     !: land elevation below which wetting/drying will be considered 
     276   INTEGER , PUBLIC, SAVE ::   nn_wdit     !: maximum number of iteration for W/D limiter 
    264277 
    265278   !!---------------------------------------------------------------------- 
     
    325338   INTEGER FUNCTION dom_oce_alloc() 
    326339      !!---------------------------------------------------------------------- 
    327       INTEGER, DIMENSION(11) :: ierr 
     340      INTEGER, DIMENSION(12) :: ierr 
    328341      !!---------------------------------------------------------------------- 
    329342      ierr(:) = 0 
     
    390403#endif 
    391404      ! 
     405      IF(ln_wd) & 
     406      ALLOCATE( wdmask(jpi,jpj), STAT=ierr(12) ) 
     407      ! 
    392408      dom_oce_alloc = MAXVAL(ierr) 
    393409      ! 
  • branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90

    r4624 r5066  
    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/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90

    r4624 r5066  
    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(mbathy(ji,jj) == 2 .AND. e3t_0(ji,jj,1) <= 0.5_wp * rn_wdmin1) THEN 
     151             fse3t_a(ji,jj,1:2) = 0.5_wp * rn_wdmin1 
     152            END IF 
     153          ENDDO 
     154        ENDDO 
     155      END IF 
     156 
    147157      ! Reconstruction of all vertical scale factors at now and before time steps 
    148158      ! ============================================================================= 
     
    793803      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag 
    794804      !! * Local declarations 
    795       INTEGER ::   jk 
     805      INTEGER ::   ji, jj, jk 
    796806      INTEGER ::   id1, id2, id3, id4, id5     ! local integers 
    797807      !!---------------------------------------------------------------------- 
     
    867877            fse3t_n(:,:,:) = e3t_0(:,:,:) 
    868878            sshn(:,:) = 0.0_wp 
     879 
     880            IF(ln_wd) THEN 
     881              DO jj = 1, jpj 
     882                DO ji = 1, jpi 
     883                  !IF(e3t_0(ji,jj,1) < 0._wp) THEN 
     884                  IF(mbathy(ji,jj) == 2 .AND. e3t_0(ji,jj,1) <= 0.5_wp * rn_wdmin1) THEN 
     885                    fse3t_b(ji,jj,1:2) = 0.5_wp * rn_wdmin1 
     886                    fse3t_n(ji,jj,1:2) = 0.5_wp * rn_wdmin1 
     887                    sshb(ji,jj) = rn_wdmin1 - bathy(ji,jj) 
     888                    sshn(ji,jj) = rn_wdmin1 - bathy(ji,jj) 
     889                  ENDIF 
     890                ENDDO 
     891              ENDDO 
     892            END IF 
     893 
    869894            IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN 
    870895               tilde_e3t_b(:,:,:) = 0.0_wp 
  • branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90

    r4754 r5066  
    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 
     
    11211122      REAL(wp) ::   zrmax, ztaper   ! temporary scalars 
    11221123      REAL(wp) ::   zrfact 
     1124      REAL(wp) ::   zsmth 
    11231125      ! 
    11241126      REAL(wp), POINTER, DIMENSION(:,:  ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2 
     
    11761178      bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) ) 
    11771179 
    1178       DO jj = 1, jpj 
    1179          DO ji = 1, jpi 
    1180            IF( bathy(ji,jj) > 0._wp )   bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) ) 
    1181          END DO 
    1182       END DO 
     1180      IF (.NOT. ln_wd) THEN 
     1181        DO jj = 1, jpj 
     1182           DO ji = 1, jpi 
     1183             IF( bathy(ji,jj) > 0._wp )   bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) ) 
     1184           END DO 
     1185        END DO 
     1186      ENDIF 
    11831187      !                                        ! ============================= 
    11841188      !                                        ! Define the envelop bathymetry   (hbatt) 
     
    11871191      zenv(:,:) = bathy(:,:) 
    11881192      ! 
    1189       ! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing 
    1190       DO jj = 1, jpj 
    1191          DO ji = 1, jpi 
    1192            IF( bathy(ji,jj) == 0._wp ) THEN 
    1193              iip1 = MIN( ji+1, jpi ) 
    1194              ijp1 = MIN( jj+1, jpj ) 
    1195              iim1 = MAX( ji-1, 1 ) 
    1196              ijm1 = MAX( jj-1, 1 ) 
    1197              IF( (bathy(iip1,jj) + bathy(iim1,jj) + bathy(ji,ijp1) + bathy(ji,ijm1) +              & 
    1198         &         bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0._wp ) THEN 
    1199                zenv(ji,jj) = rn_sbot_min 
     1193      IF (.NOT. ln_wd) THEN 
     1194        ! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing 
     1195        DO jj = 1, jpj 
     1196           DO ji = 1, jpi 
     1197             IF( bathy(ji,jj) == 0._wp ) THEN 
     1198               iip1 = MIN( ji+1, jpi ) 
     1199               ijp1 = MIN( jj+1, jpj ) 
     1200               iim1 = MAX( ji-1, 1 ) 
     1201               ijm1 = MAX( jj-1, 1 ) 
     1202               IF( (bathy(iip1,jj) + bathy(iim1,jj) + bathy(ji,ijp1) + bathy(ji,ijm1) +              & 
     1203          &         bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0._wp ) THEN 
     1204                 zenv(ji,jj) = rn_sbot_min 
     1205               ENDIF 
    12001206             ENDIF 
    1201            ENDIF 
    1202          END DO 
    1203       END DO 
     1207           END DO 
     1208        END DO 
     1209      END IF 
    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         zsmth = 0._wp 
    12431250         DO jj = 1, nlcj 
    12441251            DO ji = 1, nlci 
    12451252               iip1 = MIN( ji+1, nlci )      ! force zri = 0 on last line (ji=ncli+1 to jpi) 
    12461253               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 
     1254               IF( (zenv(ji,jj) > zsmth) .AND. (zenv(iip1,jj) > zsmth)) THEN 
    12481255                  zri(ji,jj) = ( zenv(iip1,jj  ) - zenv(ji,jj) ) / ( zenv(iip1,jj  ) + zenv(ji,jj) ) 
    12491256               END IF 
    1250                IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(ji,ijp1) > 0._wp)) THEN 
     1257               IF( (zenv(ji,jj) > zsmth) .AND. (zenv(ji,ijp1) > zsmth)) THEN 
    12511258                  zrj(ji,jj) = ( zenv(ji  ,ijp1) - zenv(ji,jj) ) / ( zenv(ji  ,ijp1) + zenv(ji,jj) ) 
    12521259               END IF 
     
    12711278      END DO                                                !     End loop     ! 
    12721279      !                                                     ! ================ ! 
    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 
     1280      IF(ln_wd) THEN 
     1281        DO jj = 1, jpj 
     1282           DO ji = 1, jpi 
     1283              zenv(ji,jj) = MAX( zenv(ji,jj), -rn_wdld )    ! fill out land bathy data 
     1284           END DO 
     1285        END DO 
     1286      ELSE 
     1287        DO jj = 1, jpj 
     1288           DO ji = 1, jpi 
     1289              zenv(ji,jj) = MAX( zenv(ji,jj), rn_sbot_min ) ! set all points to avoid undefined scale value warnings 
     1290           END DO 
     1291        END DO 
     1292      END IF 
    12781293      ! 
    12791294      ! Envelope bathymetry saved in hbatt 
     
    13051320      IF(lwp) THEN 
    13061321         WRITE(numout,*) 
    1307          WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min 
     1322         IF (.NOT. ln_wd ) THEN 
     1323           WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min 
     1324         ELSE 
     1325           WRITE(numout,*) ' zgr_sco: minimum positive depth of the envelop topography set to : ', rn_sbot_min 
     1326           WRITE(numout,*) ' zgr_sco: minimum negative depth of the envelop topography set to : ', -rn_wdld 
     1327         ENDIF 
    13081328      ENDIF 
    13091329      hbatu(:,:) = rn_sbot_min 
     
    13181338        END DO 
    13191339      END DO 
     1340 
     1341      IF(ln_wd) THEN               !avoid the zero depth on T- (U-,V-,F-) points 
     1342        DO jj = 1, jpj 
     1343          DO ji = 1, jpi 
     1344            IF(ABS(hbatt(ji,jj)) < rn_wdmin1) & 
     1345              & hbatt(ji,jj) = SIGN(1._wp, hbatt(ji,jj)) * rn_wdmin1 
     1346            IF(ABS(hbatu(ji,jj)) < rn_wdmin1) & 
     1347              & hbatu(ji,jj) = SIGN(1._wp, hbatu(ji,jj)) * rn_wdmin1 
     1348            IF(ABS(hbatv(ji,jj)) < rn_wdmin1) & 
     1349              & hbatv(ji,jj) = SIGN(1._wp, hbatv(ji,jj)) * rn_wdmin1 
     1350            IF(ABS(hbatf(ji,jj)) < rn_wdmin1) & 
     1351              & hbatf(ji,jj) = SIGN(1._wp, hbatf(ji,jj)) * rn_wdmin1 
     1352          END DO 
     1353        END DO 
     1354      END IF 
    13201355      !  
    13211356      ! Apply lateral boundary condition 
     
    13251360         DO ji = 1, jpi 
    13261361            IF( hbatu(ji,jj) == 0._wp ) THEN 
     1362               !No worries about the following line when ln_wd == .true. 
    13271363               IF( zhbat(ji,jj) == 0._wp )   hbatu(ji,jj) = rn_sbot_min 
    13281364               IF( zhbat(ji,jj) /= 0._wp )   hbatu(ji,jj) = zhbat(ji,jj) 
     
    13501386 
    13511387!!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(:,:) ) 
     1388 
     1389      IF (.NOT. ln_wd) THEN 
     1390        hift(:,:) = MIN( hift(:,:), hbatt(:,:) ) 
     1391        hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) ) 
     1392        hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) ) 
     1393        hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) ) 
     1394      END IF 
    13561395 
    13571396      IF( nprint == 1 .AND. lwp )   THEN 
     
    13951434      CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 
    13961435 
    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 
     1436      !fsdepw(:,:,:) = gdepw_0 (:,:,:) 
     1437      !fsde3w(:,:,:) = gdep3w_0(:,:,:) 
     1438      ! 
     1439      IF (.NOT. ln_wd) THEN 
     1440        where (e3t_0   (:,:,:).eq.0.0)  e3t_0(:,:,:) = 1.0 
     1441        where (e3u_0   (:,:,:).eq.0.0)  e3u_0(:,:,:) = 1.0 
     1442        where (e3v_0   (:,:,:).eq.0.0)  e3v_0(:,:,:) = 1.0 
     1443        where (e3f_0   (:,:,:).eq.0.0)  e3f_0(:,:,:) = 1.0 
     1444        where (e3w_0   (:,:,:).eq.0.0)  e3w_0(:,:,:) = 1.0 
     1445        where (e3uw_0  (:,:,:).eq.0.0)  e3uw_0(:,:,:) = 1.0 
     1446        where (e3vw_0  (:,:,:).eq.0.0)  e3vw_0(:,:,:) = 1.0 
     1447      END IF 
    14071448 
    14081449#if defined key_agrif 
     
    14491490               IF( scobot(ji,jj) >= fsdept(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk ) 
    14501491#endif 
    1451                IF( scobot(ji,jj) == 0._wp            )   mbathy(ji,jj) = 0 
    1452             END DO 
     1492            END DO 
     1493            IF(ln_wd) THEN 
     1494              IF( scobot(ji,jj) <= -(rn_wdld - rn_wdmin2)) THEN 
     1495                mbathy(ji,jj) = 0 
     1496              ELSEIF(scobot(ji,jj) <= rn_wdmin1) THEN 
     1497#if key_surge 
     1498                mbathy(ji,jj) = 1 
     1499#else 
     1500                mbathy(ji,jj) = 2 
     1501#endif 
     1502              ENDIF 
     1503            ELSE 
     1504              IF( scobot(ji,jj) == 0._wp            )   mbathy(ji,jj) = 0 
     1505            END IF 
    14531506         END DO 
    14541507      END DO 
     
    15081561         DO jj = 1, jpj 
    15091562 
    1510             IF( hbatt(ji,jj) > 0._wp) THEN 
     1563            IF( bathy(ji,jj) > 0._wp) THEN 
    15111564               DO jk = 1, mbathy(ji,jj) 
    15121565                 ! check coordinate is monotonically increasing 
    1513                  IF (fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN 
    1514                     WRITE(ctmp1,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk 
    1515                     WRITE(numout,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk 
     1566                 !RF test... 
     1567                 !IF (fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN 
     1568                 !   WRITE(ctmp1,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk 
     1569                 !   WRITE(numout,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk 
     1570                 IF (fse3w(ji,jj,jk) < 0._wp .OR. fse3t(ji,jj,jk) < 0._wp ) THEN 
     1571                    WRITE(ctmp1,*) 'ERROR zgr_sco :   e3w   or e3t   < 0  at point (i,j,k)= ', ji, jj, jk 
     1572                    WRITE(numout,*) 'ERROR zgr_sco :   e3w   or e3t   < 0  at point (i,j,k)= ', ji, jj, jk 
    15161573                    WRITE(numout,*) 'e3w',fse3w(ji,jj,:) 
    15171574                    WRITE(numout,*) 'e3t',fse3t(ji,jj,:) 
     
    15721629      INTEGER  ::   ji, jj, jk           ! dummy loop argument 
    15731630      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars 
     1631      REAL(wp) ::   ztmpu,  ztmpv,  ztmpf 
     1632      REAL(wp) ::   ztmpu1, ztmpv1, ztmpf1 
    15741633      ! 
    15751634      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 
     
    16251684      DO ji = 1, jpim1 
    16261685         DO jj = 1, jpjm1 
     1686            ! extended for Wetting/Drying case 
     1687            ztmpu = hbatt(ji,jj)+hbatt(ji+1,jj) 
     1688            ztmpv = hbatt(ji,jj)+hbatt(ji,jj+1) 
     1689            ztmpf = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) 
     1690            ztmpu1 = hbatt(ji,jj)*hbatt(ji+1,jj) 
     1691            ztmpv1 = hbatt(ji,jj)*hbatt(ji,jj+1) 
     1692            ztmpf1 = MIN(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) * & 
     1693                   & MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) 
    16271694            DO jk = 1, jpk 
    1628                z_esigtu3(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)+hbatt(ji+1,jj) ) 
    1630                z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) )   & 
    1631                   &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
    1632                z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk)     & 
    1633                   &                + hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) )   & 
    1634                   &              / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 
    1635                z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) )   & 
    1636                   &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
    1637                z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) )   & 
    1638                   &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
     1695               IF((ztmpu1 < 0._wp.OR.ABS(ztmpu) < rn_wdmin1).AND.ln_wd) THEN 
     1696                 z_esigtu3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) ) 
     1697                 z_esigwu3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji+1,jj,jk) ) 
     1698               ELSE 
     1699                 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) & 
     1700                        &              / ztmpu 
     1701                 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) & 
     1702                        &              / ztmpu 
     1703               END IF 
     1704 
     1705               IF((ztmpv1 < 0._wp.OR.ABS(ztmpv) < rn_wdmin1).AND.ln_wd) THEN 
     1706                 z_esigtv3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji,jj+1,jk) ) 
     1707                 z_esigwv3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji,jj+1,jk) ) 
     1708               ELSE 
     1709                 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) & 
     1710                        &              / ztmpv 
     1711                 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) & 
     1712                        &              / ztmpv 
     1713               END IF 
     1714 
     1715               IF((ztmpf1 < 0._wp.OR.ABS(ztmpf) < rn_wdmin1).AND.ln_wd) THEN 
     1716                 z_esigtf3(ji,jj,jk) = 0.25_wp * ( z_esigt3(ji,jj,jk)   + z_esigt3(ji+1,jj,jk)  & 
     1717                    &                            + z_esigt3(ji,jj+1,jk) + z_esigt3(ji+1,jj+1,jk) ) 
     1718               ELSE 
     1719                 z_esigtf3(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+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) & 
     1721                    &              / ztmpf 
     1722               END IF 
     1723 
    16391724               ! 
    16401725               e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 
     
    16761761      REAL(wp) ::   zsmth               ! smoothing around critical depth 
    16771762      REAL(wp) ::   zzs, zzb           ! Surface and bottom cell thickness in sigma space 
     1763      REAL(wp) ::   ztmpu, ztmpv, ztmpf 
     1764      REAL(wp) ::   ztmpu1, ztmpv1, ztmpf1 
    16781765      ! 
    16791766      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 
     
    17541841        DO jj=1,jpj-1 
    17551842 
    1756           DO jk = 1, jpk 
    1757                 z_esigtu3(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)+hbatt(ji+1,jj) ) 
    1759                 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) / & 
    1760                                     ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
    1761                 z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) +  & 
    1762                                       hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / & 
    1763                                     ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 
    1764                 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) / & 
    1765                                     ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
    1766                 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) / & 
    1767                                     ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
     1843           ! extend to suit for Wetting/Drying case 
     1844           ztmpu = hbatt(ji,jj)+hbatt(ji+1,jj) 
     1845           ztmpv = hbatt(ji,jj)+hbatt(ji,jj+1) 
     1846           ztmpf = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) 
     1847           ztmpu1 = hbatt(ji,jj)*hbatt(ji+1,jj) 
     1848           ztmpv1 = hbatt(ji,jj)*hbatt(ji,jj+1) 
     1849           ztmpf1 = MIN(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) * & 
     1850                  & MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) 
     1851           DO jk = 1, jpk 
     1852              IF((ztmpu1 < 0._wp.OR.ABS(ztmpu) < rn_wdmin1).AND.ln_wd) THEN 
     1853                z_esigtu3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) ) 
     1854                z_esigwu3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji+1,jj,jk) ) 
     1855              ELSE 
     1856                z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) & 
     1857                       &              / ztmpu 
     1858                z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) & 
     1859                       &              / ztmpu 
     1860              END IF 
     1861 
     1862              IF((ztmpv1 < 0._wp.OR.ABS(ztmpv) < rn_wdmin1).AND.ln_wd) THEN 
     1863                z_esigtv3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji,jj+1,jk) ) 
     1864                z_esigwv3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji,jj+1,jk) ) 
     1865              ELSE 
     1866                z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) & 
     1867                       &              / ztmpv 
     1868                z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) & 
     1869                       &              / ztmpv 
     1870              END IF 
     1871 
     1872              IF((ztmpf1 < 0._wp.OR.ABS(ztmpf) < rn_wdmin1).AND.ln_wd) THEN 
     1873                z_esigtf3(ji,jj,jk) = 0.25_wp * ( z_esigt3(ji,jj,jk)   + z_esigt3(ji+1,jj,jk)  & 
     1874                   &                            + z_esigt3(ji,jj+1,jk) + z_esigt3(ji+1,jj+1,jk) ) 
     1875              ELSE 
     1876                z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk)  & 
     1877                   &                + hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) & 
     1878                   &              / ztmpf 
     1879              END IF 
    17681880 
    17691881             e3t_0(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk) 
  • branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90

    r4794 r5066  
    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 
     
    402441            zvap = -zcoef0 * ( 2._wp * znad )   & 
    403442               &           * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj) 
     443 
     444            IF(ln_wd) THEN 
     445              zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 
     446              zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 
     447              zuap = zuap * zcpx(ji,jj) 
     448              zvap = zvap * zcpy(ji,jj) 
     449            ENDIF 
     450 
    404451            ! add to the general momentum trend 
    405452            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap 
     
    424471               zvap = -zcoef0 * ( 2._wp * znad )   & 
    425472                  &           * ( fsde3w(ji  ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) 
     473 
     474               IF(ln_wd) THEN 
     475                 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 
     476                 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) 
     477                 zuap = zuap * zcpx(ji,jj) 
     478                 zvap = zvap * zcpy(ji,jj) 
     479               ENDIF 
     480 
    426481               ! add to the general momentum trend 
    427482               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap 
     
    445500            zvap = -zcoef0 * ( rhd   (ji,jj+1,1) + rhd   (ji,jj,1) + 2._wp * znad )   & 
    446501               &           * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj) 
     502 
     503            IF(ln_wd) THEN 
     504              zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 
     505              zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 
     506              zuap = zuap * zcpx(ji,jj) 
     507              zvap = zvap * zcpy(ji,jj) 
     508            ENDIF 
     509 
    447510            ! add to the general momentum trend 
    448511            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap 
     
    467530               zvap = -zcoef0 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   & 
    468531                  &           * ( fsde3w(ji  ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) 
     532 
     533               IF(ln_wd) THEN 
     534                 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 
     535                 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) 
     536                 zuap = zuap * zcpx(ji,jj) 
     537                 zvap = zvap * zcpy(ji,jj) 
     538               ENDIF 
     539 
    469540               ! add to the general momentum trend 
    470541               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap 
     
    476547      ! 
    477548      CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj ) 
     549      IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 
    478550      ! 
    479551   END SUBROUTINE hpg_sco 
     
    493565      REAL(wp) ::   z1_10, cffu, cffx   !    "         " 
    494566      REAL(wp) ::   z1_12, cffv, cffy   !    "         " 
     567      LOGICAL  ::   ll_tmp1, ll_tmp2    ! local logical variables 
    495568      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zhpi, zhpj 
    496569      REAL(wp), POINTER, DIMENSION(:,:,:) ::  dzx, dzy, dzz, dzu, dzv, dzw 
    497570      REAL(wp), POINTER, DIMENSION(:,:,:) ::  drhox, drhoy, drhoz, drhou, drhov, drhow 
    498571      REAL(wp), POINTER, DIMENSION(:,:,:) ::  rho_i, rho_j, rho_k 
     572      REAL(wp), POINTER, DIMENSION(:,:)   ::  zcpx, zcpy    !W/D pressure filter 
    499573      !!---------------------------------------------------------------------- 
    500574      ! 
     
    502576      CALL wrk_alloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow ) 
    503577      CALL wrk_alloc( jpi, jpj, jpk, rho_i, rho_j, rho_k,  zhpi,  zhpj        ) 
     578      IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 
    504579      ! 
    505580 
     
    508583         IF(lwp) WRITE(numout,*) 'dyn:hpg_djc : hydrostatic pressure gradient trend' 
    509584         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate case, density Jacobian with cubic polynomial scheme' 
     585      ENDIF 
     586 
     587      IF(ln_wd) THEN 
     588        DO jj = 2, jpjm1 
     589           DO ji = 2, jpim1 
     590             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) 
     591             ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) +& 
     592                                                       & rn_wdmin1 + rn_wdmin2 
     593 
     594             IF(ll_tmp1) THEN 
     595               zcpx(ji,jj) = 1.0_wp 
     596             ELSE IF(ll_tmp2) THEN 
     597               ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 
     598               zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /& 
     599                           &     (sshn(ji+1,jj) - sshn(ji,jj))) 
     600             ELSE 
     601               zcpx(ji,jj) = 0._wp 
     602             END IF 
     603      
     604             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) 
     605             ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) +& 
     606                                                       & rn_wdmin1 + rn_wdmin2 
     607 
     608             IF(ll_tmp1) THEN 
     609               zcpy(ji,jj) = 1.0_wp 
     610             ELSE IF(ll_tmp2) THEN 
     611               ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 
     612               zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) /& 
     613                           &     (sshn(ji,jj+1) - sshn(ji,jj))) 
     614             ELSE 
     615               zcpy(ji,jj) = 0._wp 
     616             END IF 
     617           END DO 
     618        END DO 
     619        CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp ) 
    510620      ENDIF 
    511621 
     
    673783            zhpi(ji,jj,1) = ( rho_k(ji+1,jj  ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) / e1u(ji,jj) 
    674784            zhpj(ji,jj,1) = ( rho_k(ji  ,jj+1,1) - rho_k(ji,jj,1) - rho_j(ji,jj,1) ) / e2v(ji,jj) 
     785            IF(ln_wd) THEN 
     786              zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 
     787              zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 
     788            ENDIF 
    675789            ! add to the general momentum trend 
    676790            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) 
     
    692806                  &           + (  ( rho_k(ji,jj+1,jk) - rho_k(ji,jj,jk  ) )    & 
    693807                  &               -( rho_j(ji,jj  ,jk) - rho_j(ji,jj,jk-1) )  ) / e2v(ji,jj) 
     808               IF(ln_wd) THEN 
     809                 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 
     810                 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) 
     811               ENDIF 
    694812               ! add to the general momentum trend 
    695813               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) 
     
    702820      CALL wrk_dealloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow ) 
    703821      CALL wrk_dealloc( jpi, jpj, jpk, rho_i, rho_j, rho_k,  zhpi,  zhpj        ) 
     822      IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 
    704823      ! 
    705824   END SUBROUTINE hpg_djc 
     
    727846      !! The local variables for the correction term 
    728847      INTEGER  :: jk1, jis, jid, jjs, jjd 
     848      LOGICAL  :: ll_tmp1, ll_tmp2                  ! local logical variables 
    729849      REAL(wp) :: zuijk, zvijk, zpwes, zpwed, zpnss, zpnsd, zdeps 
    730850      REAL(wp) :: zrhdt1 
     
    732852      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zdept, zrhh 
    733853      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp 
     854      REAL(wp), POINTER, DIMENSION(:,:)   ::   zcpx, zcpy    !W/D pressure filter 
    734855      !!---------------------------------------------------------------------- 
    735856      ! 
    736857      CALL wrk_alloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp ) 
    737858      CALL wrk_alloc( jpi,jpj,jpk, zdept, zrhh ) 
     859      IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 
    738860      ! 
    739861      IF( kt == nit000 ) THEN 
     
    748870      znad = 0.0_wp 
    749871      IF( lk_vvl ) znad = 1._wp 
     872 
     873      IF(ln_wd) THEN 
     874        DO jj = 2, jpjm1 
     875           DO ji = 2, jpim1 
     876             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) 
     877             ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) +& 
     878                                                       & rn_wdmin1 + rn_wdmin2 
     879 
     880             IF(ll_tmp1) THEN 
     881               zcpx(ji,jj) = 1.0_wp 
     882             ELSE IF(ll_tmp2) THEN 
     883               ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 
     884               zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /& 
     885                           &     (sshn(ji+1,jj) - sshn(ji,jj))) 
     886             ELSE 
     887               zcpx(ji,jj) = 0._wp 
     888             END IF 
     889      
     890             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) 
     891             ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) +& 
     892                                                       & rn_wdmin1 + rn_wdmin2 
     893 
     894             IF(ll_tmp1.OR.ll_tmp2) THEN 
     895               zcpy(ji,jj) = 1.0_wp 
     896             ELSE IF(ll_tmp2) THEN 
     897               ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 
     898               zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) /& 
     899                           &     (sshn(ji,jj+1) - sshn(ji,jj))) 
     900             ELSE 
     901               zcpy(ji,jj) = 0._wp 
     902             END IF 
     903           END DO 
     904        END DO 
     905        CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp ) 
     906      ENDIF 
    750907 
    751908      ! Clean 3-D work arrays 
     
    9071064               ENDIF 
    9081065 
    909                ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * & 
    910                &           umask(ji,jj,jk) * tmask(ji,jj,jk) * tmask(ji+1,jj,jk) 
     1066               IF(ln_wd) THEN 
     1067                  zdpdx1 = zdpdx1 * zcpx(ji,jj) 
     1068                  zdpdx2 = zdpdx2 * zcpx(ji,jj) 
     1069                ENDIF 
     1070                ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk)  
    9111071            ENDIF 
    9121072 
     
    9641124               ENDIF 
    9651125 
    966                va(ji,jj,jk) = va(ji,jj,jk) + (zdpdy1 + zdpdy2)*& 
    967                &              vmask(ji,jj,jk)*tmask(ji,jj,jk)*tmask(ji,jj+1,jk) 
     1126               IF(ln_wd) THEN 
     1127                  zdpdy1 = zdpdy1 * zcpy(ji,jj) 
     1128                  zdpdy2 = zdpdy2 * zcpy(ji,jj) 
     1129               ENDIF 
     1130 
     1131               va(ji,jj,jk) = va(ji,jj,jk) + (zdpdy1 + zdpdy2) * vmask(ji,jj,jk) 
    9681132            ENDIF 
    9691133 
     
    9751139      CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp ) 
    9761140      CALL wrk_dealloc( jpi,jpj,jpk, zdept, zrhh ) 
     1141      IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 
    9771142      ! 
    9781143   END SUBROUTINE hpg_prj 
  • branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r4624 r5066  
    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 
     
    4041   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 
     
    137139      INTEGER, INTENT(in)  ::   kt   ! ocean time-step index 
    138140      ! 
    139       LOGICAL  ::   ll_fw_start        ! if true, forward integration  
    140       LOGICAL  ::   ll_init             ! if true, special startup of 2d equations 
    141       INTEGER  ::   ji, jj, jk, jn        ! dummy loop indices 
    142       INTEGER  ::   ikbu, ikbv, noffset      ! local integers 
    143       REAL(wp) ::   zraur, z1_2dt_b, z2dt_bf    ! local scalars 
    144       REAL(wp) ::   zx1, zy1, zx2, zy2         !   -      - 
    145       REAL(wp) ::   z1_12, z1_8, z1_4, z1_2    !   -      - 
    146       REAL(wp) ::   zu_spg, zv_spg             !   -      - 
    147       REAL(wp) ::   zhura, zhvra               !   -      - 
    148       REAL(wp) ::   za0, za1, za2, za3           !   -      - 
     141      LOGICAL  ::   ll_fw_start                 ! if true, forward integration 
     142      LOGICAL  ::   ll_init                         ! if true, special startup of 2d equations 
     143      LOGICAL  ::   ll_tmp1, ll_tmp2            ! local logical variables used in W/D 
     144      INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices 
     145      INTEGER  ::   ikbu, ikbv, noffset         ! local integers 
     146      REAL(wp) ::   zraur, z1_2dt_b, z2dt_bf    ! local scalars 
     147      REAL(wp) ::   zx1, zy1, zx2, zy2          !   -      - 
     148      REAL(wp) ::   z1_12, z1_8, z1_4, z1_2     !   -      - 
     149      REAL(wp) ::   zu_spg, zv_spg              !   -      - 
     150      REAL(wp) ::   zhura, zhvra                !   -      - 
     151      REAL(wp) ::   za0, za1, za2, za3          !   -      - 
     152 
     153      REAL(wp) ::   ztmp                       ! temporary vaialbe 
    149154      ! 
    150155      REAL(wp), POINTER, DIMENSION(:,:) :: zun_e, zvn_e, zsshp2_e 
     
    154159      REAL(wp), POINTER, DIMENSION(:,:) :: zsshu_a, zsshv_a 
    155160      REAL(wp), POINTER, DIMENSION(:,:) :: zhf 
     161      REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy                 ! Wetting/Dying gravity filter coef. 
    156162      !!---------------------------------------------------------------------- 
    157163      ! 
     
    165171      CALL wrk_alloc( jpi, jpj, zsshu_a, zsshv_a                                   ) 
    166172      CALL wrk_alloc( jpi, jpj, zhf ) 
     173      IF(ln_wd) CALL wrk_alloc( jpi, jpj, zcpx, zcpy ) 
    167174      ! 
    168175      !                                         !* Local constant initialization 
     
    173180      zraur = 1._wp / rau0 
    174181      ! 
    175       IF( kt == nit000 .AND. neuler == 0 ) THEN    ! reciprocal of baroclinic time step  
     182      IF( kt == nit000 .AND. neuler == 0 ) THEN  ! reciprocal of baroclinic time step  
    176183        z2dt_bf = rdt 
    177184      ELSE 
     
    180187      z1_2dt_b = 1.0_wp / z2dt_bf  
    181188      ! 
    182       ll_init = ln_bt_av                           ! if no time averaging, then no specific restart  
     189      ll_init = ln_bt_av                       ! if no time averaging, then no specific restart  
    183190      ll_fw_start = .FALSE. 
    184191      ! 
    185                                                        ! time offset in steps for bdy data update 
     192                                                     ! time offset in steps for bdy data update 
    186193      IF (.NOT.ln_bt_fw) THEN ; noffset=-2*nn_baro ; ELSE ;  noffset = 0 ; ENDIF 
    187194      ! 
    188       IF( kt == nit000 ) THEN                !* initialisation 
     195      IF( kt == nit000 ) THEN             !* initialisation 
    189196         ! 
    190197         IF(lwp) WRITE(numout,*) 
     
    365372      !                                   ! ---------------------------------------------------- 
    366373      IF( lk_vvl ) THEN                         ! Variable volume : remove surface pressure gradient 
    367          DO jj = 2, jpjm1  
    368             DO ji = fs_2, fs_jpim1   ! vector opt. 
    369                zu_trd(ji,jj) = zu_trd(ji,jj) - grav * (  sshn(ji+1,jj  ) - sshn(ji  ,jj  )  ) / e1u(ji,jj) 
    370                zv_trd(ji,jj) = zv_trd(ji,jj) - grav * (  sshn(ji  ,jj+1) - sshn(ji  ,jj  )  ) / e2v(ji,jj) 
    371             END DO 
    372          END DO 
     374        IF(ln_wd) THEN                    ! Calculating and applying W/D gravity filters 
     375          DO jj = 2, jpjm1 
     376             DO ji = 2, jpim1 
     377                ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) 
     378                ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) + & 
     379                                                          & rn_wdmin1 + rn_wdmin2 
     380                IF(ll_tmp1) THEN 
     381                  zcpx(ji,jj) = 1.0_wp 
     382                ELSE IF(ll_tmp2) THEN 
     383                   ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 
     384                  zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /& 
     385                              &     (sshn(ji+1,jj) - sshn(ji,jj))) 
     386                ELSE 
     387                   zcpx(ji,jj) = 0._wp 
     388                END IF 
     389 
     390                ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) 
     391                ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) + & 
     392                                                          & rn_wdmin1 + rn_wdmin2 
     393                IF(ll_tmp1) THEN 
     394                   zcpy(ji,jj) = 1.0_wp 
     395                ELSE IF(ll_tmp2) THEN 
     396                   ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 
     397                  zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) + bathy(ji,jj)) /& 
     398                              &     (sshn(ji,jj+1) - sshn(ji,jj))) 
     399                ELSE 
     400                  zcpy(ji,jj) = 0._wp 
     401                END IF 
     402             END DO 
     403           END DO 
     404           CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp ) 
     405        ENDIF 
     406 
     407        IF(ln_wd) THEN                    ! Calculating and applying W/D gravity filters 
     408           DO jj = 2, jpjm1 
     409              DO ji = 2, jpim1 
     410                 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj  ) - sshn(ji  ,jj ) ) / & 
     411                               &                 e1u(ji,jj) * zcpx(ji,jj) 
     412                 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji  ,jj+1) - sshn(ji  ,jj ) ) / & 
     413                               &                 e2v(ji,jj) * zcpy(ji,jj) 
     414              END DO 
     415           END DO 
     416         ELSE 
     417           DO jj = 2, jpjm1 
     418              DO ji = fs_2, fs_jpim1   ! vector opt. 
     419                 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * (  sshn(ji+1,jj  ) - sshn(ji  ,jj  )  ) / e1u(ji,jj) 
     420                 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * (  sshn(ji  ,jj+1) - sshn(ji  ,jj  )  ) / e2v(ji,jj) 
     421              END DO 
     422           END DO 
     423        END IF 
    373424      ENDIF 
    374425 
     
    464515      !                                             ! ==================== !   
    465516      ! Initialize barotropic variables:       
     517 
     518      IF(ln_wd) THEN      !preserve the positivity of water depth 
     519                          !ssh[b,n,a] should have already been processed for this 
     520         sshbb_e(:,:) = MAX(sshbb_e(:,:), rn_wdmin1 - bathy(:,:)) 
     521         sshb_e(:,:)  = MAX(sshb_e(:,:) , rn_wdmin1 - bathy(:,:)) 
     522      ENDIF 
     523 
    466524      IF (ln_bt_fw) THEN                  ! FORWARD integration: start from NOW fields                     
    467525         sshn_e(:,:) = sshn (:,:)             
     
    537595            zhup2_e (:,:) = hu_0(:,:) + zwx(:,:)                ! Ocean depth at U- and V-points 
    538596            zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:) 
     597            IF(ln_wd) THEN 
     598              zhup2_e(:,:) = MAX(zhup2_e (:,:), rn_wdmin1) 
     599              zhvp2_e(:,:) = MAX(zhvp2_e (:,:), rn_wdmin1) 
     600            END IF 
    539601         ELSE 
    540602            zhup2_e (:,:) = hu(:,:) 
     
    575637         ENDIF 
    576638#endif 
     639         IF(ln_wd) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt) 
    577640         ! 
    578641         ! Sum over sub-time-steps to compute advective velocities 
     
    589652         END DO 
    590653         ssha_e(:,:) = (  sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) )  ) * tmask(:,:,1) 
     654         IF(ln_wd) ssha_e(:,:) = MAX(ssha_e(:,:), rn_wdmin1 - bathy(:,:))  
    591655         CALL lbc_lnk( ssha_e, 'T',  1._wp ) 
    592656 
     
    652716               END DO 
    653717            END DO 
     718 
     719            IF(ln_wd) THEN 
     720              zhust_e(:,:) = MAX(zhust_e (:,:), rn_wdmin1 ) 
     721              zhvst_e(:,:) = MAX(zhvst_e (:,:), rn_wdmin1 ) 
     722            END IF 
     723            !shall we call lbc_lnk for zhu(v)st_e() here? 
     724 
    654725         ENDIF 
    655726         ! 
     
    717788         zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * zvn_e(:,:) * hvr_e(:,:) 
    718789         ! 
    719          ! Surface pressure trend: 
    720          DO jj = 2, jpjm1 
    721             DO ji = fs_2, fs_jpim1   ! vector opt. 
    722                ! Add surface pressure gradient 
    723                zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj) 
    724                zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj) 
    725                zwx(ji,jj) = zu_spg 
    726                zwy(ji,jj) = zv_spg 
    727             END DO 
    728          END DO 
     790         IF(ln_wd) THEN                    ! Calculating and applying W/D gravity filters 
     791           DO jj = 2, jpjm1 
     792              DO ji = 2, jpim1 
     793                 ll_tmp1 = MIN(zsshp2_e(ji,jj), zsshp2_e(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) 
     794                 ll_tmp2 = MAX(zsshp2_e(ji,jj), zsshp2_e(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) + & 
     795                                                           & rn_wdmin1 + rn_wdmin2 
     796                 IF(ll_tmp1) THEN 
     797                   zcpx(ji,jj) = 1.0_wp 
     798                 ELSE IF(ll_tmp2) THEN 
     799                    ! no worries about zsshp2_e(ji+1,jj)-zsshp2_e(ji,jj) = 0, it won't happen ! here 
     800                   zcpx(ji,jj) = ABS((zsshp2_e(ji+1,jj) + bathy(ji+1,jj) - zsshp2_e(ji,jj) - bathy(ji,jj)) /& 
     801                               & (zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj))) 
     802                 ELSE 
     803                    zcpx(ji,jj) = 0._wp 
     804                 END IF 
     805 
     806                 ll_tmp1 = MIN(zsshp2_e(ji,jj), zsshp2_e(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) 
     807                 ll_tmp2 = MAX(zsshp2_e(ji,jj), zsshp2_e(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) + & 
     808                                                           & rn_wdmin1 + rn_wdmin2 
     809                 IF(ll_tmp1) THEN 
     810                    zcpy(ji,jj) = 1.0_wp 
     811                 ELSE IF(ll_tmp2) THEN 
     812                    ! no worries about zsshp2_e(ji,jj+1)-zsshp2_e(ji,jj) = 0, it won't happen ! here 
     813                   zcpy(ji,jj) = ABS((zsshp2_e(ji,jj+1) + bathy(ji,jj+1) - zsshp2_e(ji,jj) - bathy(ji,jj)) /& 
     814                               & (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj))) 
     815                 ELSE 
     816                   zcpy(ji,jj) = 0._wp 
     817                 END IF 
     818              END DO 
     819            END DO 
     820            CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp ) 
     821         ENDIF 
     822 
     823         IF(ln_wd) THEN 
     824           DO jj = 2, jpjm1 
     825              DO ji = 2, jpim1 
     826                 ! Add surface pressure gradient 
     827                 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj) 
     828                 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj) 
     829                 zwx(ji,jj) = zu_spg * zcpx(ji,jj) 
     830                 zwy(ji,jj) = zv_spg * zcpy(ji,jj) 
     831              END DO 
     832           END DO 
     833         ELSE 
     834           DO jj = 2, jpjm1 
     835              DO ji = fs_2, fs_jpim1   ! vector opt. 
     836                 ! Add surface pressure gradient 
     837                 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj) 
     838                 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj) 
     839                 zwx(ji,jj) = zu_spg 
     840                 zwy(ji,jj) = zv_spg 
     841              END DO 
     842           END DO 
     843         END IF 
    729844         ! 
    730845         ! Set next velocities: 
     
    750865               DO ji = fs_2, fs_jpim1   ! vector opt. 
    751866 
    752                   zhura = umask(ji,jj,1)/(hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - umask(ji,jj,1)) 
    753                   zhvra = vmask(ji,jj,1)/(hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - vmask(ji,jj,1)) 
    754  
    755                   ua_e(ji,jj) = (                hu_e(ji,jj)  *  zun_e(ji,jj)   &  
    756                             &     + rdtbt * ( zhust_e(ji,jj)  *    zwx(ji,jj)   &  
     867                  IF(ln_wd) THEN 
     868                    zhura = MAX(hu_0(ji,jj) + zsshu_a(ji,jj), rn_wdmin1) 
     869                    zhvra = MAX(hv_0(ji,jj) + zsshv_a(ji,jj), rn_wdmin1) 
     870                  ELSE 
     871                    zhura = hu_0(ji,jj) + zsshu_a(ji,jj) 
     872                    zhvra = hv_0(ji,jj) + zsshv_a(ji,jj) 
     873                  END IF 
     874 
     875                  zhura = umask(ji,jj,1)/(zhura + 1._wp - umask(ji,jj,1)) 
     876                  zhvra = vmask(ji,jj,1)/(zhvra + 1._wp - vmask(ji,jj,1)) 
     877 
     878                  ua_e(ji,jj) = (                hu_e(ji,jj)  *  zun_e(ji,jj)   & 
     879                            &     + rdtbt * ( zhust_e(ji,jj)  *    zwx(ji,jj)   & 
    757880                            &               + zhup2_e(ji,jj)  * zu_trd(ji,jj)   & 
    758881                            &               +      hu(ji,jj)  * zu_frc(ji,jj) ) & 
     
    770893         IF( lk_vvl ) THEN                             !* Update ocean depth (variable volume case only) 
    771894            !                                          !  ----------------------------------------------         
    772             hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 
    773             hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 
     895            IF(ln_wd) THEN 
     896              hu_e (:,:) = MAX(hu_0(:,:) + zsshu_a(:,:), rn_wdmin1) 
     897              hv_e (:,:) = MAX(hv_0(:,:) + zsshv_a(:,:), rn_wdmin1) 
     898            ELSE 
     899              hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 
     900              hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 
     901            END IF 
     902 
    774903            hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1._wp - umask(:,:,1) ) 
    775904            hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1._wp - vmask(:,:,1) ) 
     
    9031032      CALL wrk_dealloc( jpi, jpj, zsshu_a, zsshv_a                                   ) 
    9041033      CALL wrk_dealloc( jpi, jpj, zhf ) 
     1034      IF(ln_wd) CALL wrk_alloc( jpi, jpj, zcpx, zcpy ) 
    9051035      ! 
    9061036      IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_ts') 
  • branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90

    r4486 r5066  
    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  
     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      ! 
    101107      !                                           !------------------------------! 
    102108      !                                           !   After Sea Surface Height   ! 
     
    110116      ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations. 
    111117      !  
    112       z1_rau0 = 0.5_wp * r1_rau0 
    113118      ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * tmask(:,:,1) 
    114119 
  • branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90

    r4794 r5066  
    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 ) 
    259269 
    260270! Force values for AGRIF zoom (cf. agrif_user.F90) 
     
    311321         WRITE( numond, namctl ) 
    312322         WRITE( numond, namcfg ) 
     323         WRITE( numond, namwad ) 
    313324      ENDIF 
    314325 
     
    527538         WRITE(numout,*) '      lateral cond. type (between 0 and 6) jperio = ', jperio    
    528539      ENDIF 
     540 
     541      IF(lwp) THEN                  ! control print 
     542         WRITE(numout,*) 
     543         WRITE(numout,*) 'namwad  : wetting and drying initialisation through namelist read' 
     544         WRITE(numout,*) '~~~~~~~ ' 
     545         WRITE(numout,*) '   Namelist namwad' 
     546         WRITE(numout,*) '   key to turn on/off wetting/drying              ln_wd      = ', ln_wd 
     547         WRITE(numout,*) '   minimum water depth on dried cells             rn_wdmin1  = ', rn_wdmin1 
     548         WRITE(numout,*) '   tolerrance of min water depth on dried cells   rn_wdmin2  = ', rn_wdmin2 
     549         WRITE(numout,*) '   land elevation below which wad considered      rn_wdld    = ', rn_wdld 
     550         WRITE(numout,*) '   maximum number of iteration for W/D limiter    nn_wdit    = ', nn_wdit 
     551      ENDIF  
     552 
    529553      !                             ! Parameter control 
    530554      ! 
  • branches/UKMO/2014_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/step_oce.F90

    r4328 r5066  
    118118   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.