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

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

Changeset 5954


Ignore:
Timestamp:
2015-11-30T15:04:08+01:00 (8 years ago)
Author:
rfurner
Message:

added surge code from 2014_Surge_Modelling branch

Location:
branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM
Files:
1 added
30 edited
4 copied

Legend:

Unmodified
Added
Removed
  • branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/CONFIG/SHARED/namelist_ref

    r5942 r5954  
    3131   nn_leapy    =       0   !  Leap year calendar (1) or not (0) 
    3232   ln_rstart   = .false.   !  start from rest (F) or from a restart file (T) 
     33   ln_rstdate  = .false.   !  Name restart dump by date instead of time step  
    3334   nn_euler    =       1   !  = 0 : start with forward time step if ln_rstart=T 
    3435   nn_rstctl   =       0   !  restart control ==> activated only if ln_rstart=T 
     
    336337   rn_vfac     = 0.        !  multiplicative factor for ocean/ice velocity 
    337338                           !  in the calculation of the wind stress (0.=absolute winds or 1.=relative winds) 
     339   ln_charnock = .false.   ! logical flag for charnock wind stress in surge model(true) or not(false) 
    338340/ 
    339341!----------------------------------------------------------------------- 
     
    12921294   rn_htrmax         =  200.0   ! max. depth of transition range 
    12931295/ 
     1296-----------------------------------------------------------------------  
     1297&namwad       !   Wetting and Drying namelist  
     1298!-----------------------------------------------------------------------  
     1299   ln_wd = .false.   !: key to turn on/off wetting/drying (T: on, F: off)  
     1300   rn_wdmin1=0.1     !: minimum water depth on dried cells  
     1301   rn_wdmin2 = 0.01  !: tolerrance of minimum water depth on dried cells  
     1302   rn_wdld   = 20.0  !: land elevation below which wetting/drying will be considered  
     1303   nn_wdit   =   10  !: maximum number of iteration for W/D limiter  
     1304/  
  • branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/CONFIG/cfg.txt

    r5416 r5954  
    1111ORCA2_LIM OPA_SRC LIM_SRC_2 NST_SRC 
    1212ORCA2_OFF_PISCES OPA_SRC OFF_SRC TOP_SRC 
     13AMM7_SURGE OPA_SRC 
  • branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90

    r4990 r5954  
    509509         ALLOCATE( dta_global(jpbdtau, 1, jpk) ) 
    510510         IF ( icount>0 ) ALLOCATE( dta_global2(jpbdtas, nrimmax, jpk) ) 
     511 
     512         nbidta(:,:,:)=0._wp 
     513         nbjdta(:,:,:)=0._wp 
     514         nbrdta(:,:,:)=0._wp 
    511515         !  
    512516      ENDIF 
  • branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90

    r5123 r5954  
    1010   !!            3.5  ! 2012     (S. Mocavero, I. Epicoco) Add arrays associated 
    1111   !!                             to the optimization of BDY communications 
     12   !!            3.6.?! 2014     (H. Liu) Add arrays associated with Wetting and Drying 
    1213   !!---------------------------------------------------------------------- 
    1314 
     
    270271   INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  :: nicoa, njcoa       !: ??? 
    271272#endif 
     273 
     274   !!---------------------------------------------------------------------- 
     275   !! critical depths,limiters,and masks for  Wetting and Drying 
     276   !! --------------------------------------------------------------------- 
     277 
     278   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wdmask             !: u- and v- limiter 
     279 
     280   LOGICAL,  PUBLIC, SAVE ::   ln_wd       !: key to turn on/off wetting/drying (T: on, F: off) 
     281   REAL(wp), PUBLIC, SAVE ::   rn_wdmin1   !: minimum water depth on dried cells 
     282   REAL(wp), PUBLIC, SAVE ::   rn_wdmin2   !: tolerrance of minimum water depth on dried cells 
     283   REAL(wp), PUBLIC, SAVE ::   rn_wdld     !: land elevation below which wetting/drying will be considered 
     284   INTEGER , PUBLIC, SAVE ::   nn_wdit     !: maximum number of iteration for W/D limiter 
    272285 
    273286   !!---------------------------------------------------------------------- 
     
    408421#endif 
    409422      ! 
     423      IF(ln_wd) & 
     424      ALLOCATE( wdmask(jpi,jpj), STAT=ierr(12) ) 
     425      ! 
    410426      dom_oce_alloc = MAXVAL(ierr) 
    411427      ! 
  • branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90

    r5363 r5954  
    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 
     
    137138      NAMELIST/namrun/ cn_ocerst_indir, cn_ocerst_outdir, nn_stocklist, ln_rst_list,               & 
    138139         &             nn_no   , cn_exp    , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl,   & 
     140         &             ln_rstdate                                                                  & 
    139141         &             nn_it000, nn_itend  , nn_date0    , nn_leapy     , nn_istate , nn_stock ,   & 
    140142         &             nn_write, ln_dimgnnn, ln_mskland  , ln_cfmeta    , ln_clobber, nn_chunksz, nn_euler 
     
    174176         WRITE(numout,*) '      restart output directory        cn_ocerst_outdir= ', cn_ocerst_outdir 
    175177         WRITE(numout,*) '      restart logical                 ln_rstart  = ', ln_rstart 
     178         WRITE(numout,*) '      datestamping of restarts        ln_rstdate = ', ln_rstdate 
     179         WRITE(numout,*) '      restart input directory         cn_ocerst_indir= ', cn_ocerst_indir 
     180         WRITE(numout,*) '      restart output directory        cn_ocerst_outdir= ', cn_ocerst_outdir 
    176181         WRITE(numout,*) '      start with forward time step    nn_euler   = ', nn_euler 
    177182         WRITE(numout,*) '      control of time step            nn_rstctl  = ', nn_rstctl 
  • branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90

    r5506 r5954  
    147147      fse3t_a(:,:,jpk) = e3t_0(:,:,jpk) 
    148148 
     149!      IF(ln_wd) THEN 
     150!        DO jj = 1, jpj 
     151!          DO ji = 1, jpi 
     152!            IF(mbathy(ji,jj) == 2 .AND. e3t_0(ji,jj,1) <= 0.5_wp * rn_wdmin1) THEN 
     153!             fse3t_a(ji,jj,1:2) = 0.5_wp * rn_wdmin1 
     154!             fse3t_n(ji,jj,1:2) = 0.5_wp * rn_wdmin1 
     155!             fse3t_b(ji,jj,1:2) = 0.5_wp * rn_wdmin1  
     156!            END IF 
     157!          ENDDO 
     158!        ENDDO 
     159!      END IF 
     160 
    149161      ! Reconstruction of all vertical scale factors at now and before time steps 
    150162      ! ============================================================================= 
     
    717729            DO jj = 1, jpjm1 
    718730               DO ji = 1, fs_jpim1   ! vector opt. 
    719                   pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * r1_e12u(ji,jj)                                   & 
     731!                  pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * r1_e12u(ji,jj)                                   & 
     732                  pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)   * r1_e12u(ji,jj)        & 
    720733                     &                       * (   e12t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     & 
    721734                     &                           + e12t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) ) 
     
    735748            DO jj = 1, jpjm1 
    736749               DO ji = 1, fs_jpim1   ! vector opt. 
    737                   pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) * r1_e12v(ji,jj)                                   & 
    738                      &                       * (   e12t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     & 
    739                      &                           + e12t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) ) 
     750!                  pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) * r1_e12v(ji,jj)                                   & 
     751                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)   * r1_e12v(ji,jj)        & 
     752                        &                       * (   e12t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     & 
     753                        &                           + e12t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) ) 
    740754               END DO 
    741755            END DO 
     
    753767            DO jj = 1, jpjm1 
    754768               DO ji = 1, fs_jpim1   ! vector opt. 
    755                   pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) * r1_e12f(ji,jj)               & 
     769!                  pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) * r1_e12f(ji,jj)               & 
     770                  pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)      & 
     771                     &                       * r1_e12f(ji,jj)                                                     & 
    756772                     &                       * (   e12u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     & 
    757773                     &                           + e12u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) ) 
     
    817833      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag 
    818834      !! * Local declarations 
    819       INTEGER ::   jk 
     835      INTEGER ::   ji, jj, jk 
    820836      INTEGER ::   id1, id2, id3, id4, id5     ! local integers 
    821837      !!---------------------------------------------------------------------- 
     
    905921            fse3t_n(:,:,:) = e3t_0(:,:,:) 
    906922            sshn(:,:) = 0.0_wp 
     923 
     924            IF(ln_wd) THEN 
     925              DO jj = 1, jpj 
     926                DO ji = 1, jpi 
     927                  !IF(e3t_0(ji,jj,1) < 0._wp) THEN 
     928                  !IF(mbathy(ji,jj) == 2 .AND. e3t_0(ji,jj,1) <= 0.5_wp * rn_wdmin1) THEN 
     929                  IF( e3t_0(ji,jj,1) <= 0.5_wp * rn_wdmin1) THEN 
     930                    fse3t_b(ji,jj,:) = 0.5_wp * rn_wdmin1 
     931                    fse3t_n(ji,jj,:) = 0.5_wp * rn_wdmin1 
     932                    fse3t_a(ji,jj,:) = 0.5_wp * rn_wdmin1  
     933                    sshb(ji,jj) = rn_wdmin1 - bathy(ji,jj) 
     934                    sshn(ji,jj) = rn_wdmin1 - bathy(ji,jj) 
     935                    ssha(ji,jj) = rn_wdmin1 - bathy(ji,jj) 
     936                  ENDIF 
     937                ENDDO 
     938              ENDDO 
     939            END IF 
     940 
    907941            IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN 
    908942               tilde_e3t_b(:,:,:) = 0.0_wp 
  • branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90

    r5506 r5954  
    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   !!            3.6  ! 2014-11  (P. Mathiot and C. Harris) add ice shelf capabilitye   
    2021   !!---------------------------------------------------------------------- 
     
    316317      zrefdep = 10._wp - 0.1_wp * MINVAL( e3w_1d )                   ! ref. depth with tolerance (10% of minimum layer thickness) 
    317318      nlb10 = MINLOC( gdepw_1d, mask = gdepw_1d > zrefdep, dim = 1 ) ! shallowest W level Below ~10m 
     319      nlb10 = MAX(nlb10, 2)                                          ! prevent nla10 = 0  
    318320      nla10 = nlb10 - 1                                              ! deepest    W level Above ~10m 
    319321!!gm end bug 
     
    17961798      REAL(wp) ::   zrmax, ztaper   ! temporary scalars 
    17971799      REAL(wp) ::   zrfact 
     1800      REAL(wp) ::   zsmth 
    17981801      ! 
    17991802      REAL(wp), POINTER, DIMENSION(:,:  ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2 
     
    18511854      bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) ) 
    18521855 
    1853       DO jj = 1, jpj 
    1854          DO ji = 1, jpi 
    1855            IF( bathy(ji,jj) > 0._wp )   bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) ) 
    1856          END DO 
    1857       END DO 
     1856      IF (.NOT. ln_wd) THEN 
     1857        DO jj = 1, jpj 
     1858           DO ji = 1, jpi 
     1859             IF( bathy(ji,jj) > 0._wp )   bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) ) 
     1860           END DO 
     1861        END DO 
     1862      ENDIF 
    18581863      !                                        ! ============================= 
    18591864      !                                        ! Define the envelop bathymetry   (hbatt) 
     
    18621867      zenv(:,:) = bathy(:,:) 
    18631868      ! 
    1864       ! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing 
    1865       DO jj = 1, jpj 
    1866          DO ji = 1, jpi 
    1867            IF( bathy(ji,jj) == 0._wp ) THEN 
    1868              iip1 = MIN( ji+1, jpi ) 
    1869              ijp1 = MIN( jj+1, jpj ) 
    1870              iim1 = MAX( ji-1, 1 ) 
    1871              ijm1 = MAX( jj-1, 1 ) 
    1872              IF( (bathy(iip1,jj) + bathy(iim1,jj) + bathy(ji,ijp1) + bathy(ji,ijm1) +              & 
    1873         &         bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0._wp ) THEN 
    1874                zenv(ji,jj) = rn_sbot_min 
     1869      IF (.NOT. ln_wd) THEN 
     1870        ! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing 
     1871        DO jj = 1, jpj 
     1872           DO ji = 1, jpi 
     1873             IF( bathy(ji,jj) == 0._wp ) THEN 
     1874               iip1 = MIN( ji+1, jpi ) 
     1875               ijp1 = MIN( jj+1, jpj ) 
     1876               iim1 = MAX( ji-1, 1 ) 
     1877               ijm1 = MAX( jj-1, 1 ) 
     1878               IF( (bathy(iip1,jj) + bathy(iim1,jj) + bathy(ji,ijp1) + bathy(ji,ijm1) +              & 
     1879          &         bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0._wp ) THEN 
     1880                 zenv(ji,jj) = rn_sbot_min 
     1881               ENDIF 
    18751882             ENDIF 
    1876            ENDIF 
    1877          END DO 
    1878       END DO 
     1883           END DO 
     1884        END DO 
     1885      END IF 
    18791886      ! apply lateral boundary condition   CAUTION: keep the value when the lbc field is zero 
    18801887      CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' ) 
     
    19161923         zri(:,:) = 0._wp 
    19171924         zrj(:,:) = 0._wp 
     1925         zsmth = 0._wp 
    19181926         DO jj = 1, nlcj 
    19191927            DO ji = 1, nlci 
    19201928               iip1 = MIN( ji+1, nlci )      ! force zri = 0 on last line (ji=ncli+1 to jpi) 
    19211929               ijp1 = MIN( jj+1, nlcj )      ! force zrj = 0 on last raw  (jj=nclj+1 to jpj) 
    1922                IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(iip1,jj) > 0._wp)) THEN 
     1930               IF( (zenv(ji,jj) > zsmth) .AND. (zenv(iip1,jj) > zsmth)) THEN 
    19231931                  zri(ji,jj) = ( zenv(iip1,jj  ) - zenv(ji,jj) ) / ( zenv(iip1,jj  ) + zenv(ji,jj) ) 
    19241932               END IF 
    1925                IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(ji,ijp1) > 0._wp)) THEN 
     1933               IF( (zenv(ji,jj) > zsmth) .AND. (zenv(ji,ijp1) > zsmth)) THEN 
    19261934                  zrj(ji,jj) = ( zenv(ji  ,ijp1) - zenv(ji,jj) ) / ( zenv(ji  ,ijp1) + zenv(ji,jj) ) 
    19271935               END IF 
     
    19461954      END DO                                                !     End loop     ! 
    19471955      !                                                     ! ================ ! 
    1948       DO jj = 1, jpj 
    1949          DO ji = 1, jpi 
    1950             zenv(ji,jj) = MAX( zenv(ji,jj), rn_sbot_min ) ! set all points to avoid undefined scale value warnings 
    1951          END DO 
    1952       END DO 
     1956      IF(ln_wd) THEN 
     1957        DO jj = 1, jpj 
     1958           DO ji = 1, jpi 
     1959              zenv(ji,jj) = MAX( zenv(ji,jj), -rn_wdld )    ! fill out land bathy data 
     1960           END DO 
     1961        END DO 
     1962      ELSE 
     1963        DO jj = 1, jpj 
     1964           DO ji = 1, jpi 
     1965              zenv(ji,jj) = MAX( zenv(ji,jj), rn_sbot_min ) ! set all points to avoid undefined scale value warnings 
     1966           END DO 
     1967        END DO 
     1968      END IF 
    19531969      ! 
    19541970      ! Envelope bathymetry saved in hbatt 
     
    19801996      IF(lwp) THEN 
    19811997         WRITE(numout,*) 
    1982          WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min 
     1998         IF (.NOT. ln_wd ) THEN 
     1999           WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min 
     2000         ELSE 
     2001           WRITE(numout,*) ' zgr_sco: minimum positive depth of the envelop topography set to : ', rn_sbot_min 
     2002           WRITE(numout,*) ' zgr_sco: minimum negative depth of the envelop topography set to : ', -rn_wdld 
     2003         ENDIF 
    19832004      ENDIF 
    19842005      hbatu(:,:) = rn_sbot_min 
     
    19932014        END DO 
    19942015      END DO 
     2016 
     2017      IF(ln_wd) THEN               !avoid the zero depth on T- (U-,V-,F-) points 
     2018        DO jj = 1, jpj 
     2019          DO ji = 1, jpi 
     2020            IF(ABS(hbatt(ji,jj)) < rn_wdmin1) & 
     2021              & hbatt(ji,jj) = SIGN(1._wp, hbatt(ji,jj)) * rn_wdmin1 
     2022            IF(ABS(hbatu(ji,jj)) < rn_wdmin1) & 
     2023              & hbatu(ji,jj) = SIGN(1._wp, hbatu(ji,jj)) * rn_wdmin1 
     2024            IF(ABS(hbatv(ji,jj)) < rn_wdmin1) & 
     2025              & hbatv(ji,jj) = SIGN(1._wp, hbatv(ji,jj)) * rn_wdmin1 
     2026            IF(ABS(hbatf(ji,jj)) < rn_wdmin1) & 
     2027              & hbatf(ji,jj) = SIGN(1._wp, hbatf(ji,jj)) * rn_wdmin1 
     2028          END DO 
     2029        END DO 
     2030      END IF 
    19952031      !  
    19962032      ! Apply lateral boundary condition 
     
    20002036         DO ji = 1, jpi 
    20012037            IF( hbatu(ji,jj) == 0._wp ) THEN 
     2038               !No worries about the following line when ln_wd == .true. 
    20022039               IF( zhbat(ji,jj) == 0._wp )   hbatu(ji,jj) = rn_sbot_min 
    20032040               IF( zhbat(ji,jj) /= 0._wp )   hbatu(ji,jj) = zhbat(ji,jj) 
     
    20252062 
    20262063!!bug:  key_helsinki a verifer 
    2027       hift(:,:) = MIN( hift(:,:), hbatt(:,:) ) 
    2028       hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) ) 
    2029       hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) ) 
    2030       hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) ) 
     2064 
     2065      IF (.NOT. ln_wd) THEN 
     2066        hift(:,:) = MIN( hift(:,:), hbatt(:,:) ) 
     2067        hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) ) 
     2068        hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) ) 
     2069        hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) ) 
     2070      END IF 
    20312071 
    20322072      IF( nprint == 1 .AND. lwp )   THEN 
     
    20702110      CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 
    20712111 
    2072       fsdepw(:,:,:) = gdepw_0 (:,:,:) 
    2073       fsde3w(:,:,:) = gdep3w_0(:,:,:) 
    2074       ! 
    2075       where (e3t_0   (:,:,:).eq.0.0)  e3t_0(:,:,:) = 1.0 
    2076       where (e3u_0   (:,:,:).eq.0.0)  e3u_0(:,:,:) = 1.0 
    2077       where (e3v_0   (:,:,:).eq.0.0)  e3v_0(:,:,:) = 1.0 
    2078       where (e3f_0   (:,:,:).eq.0.0)  e3f_0(:,:,:) = 1.0 
    2079       where (e3w_0   (:,:,:).eq.0.0)  e3w_0(:,:,:) = 1.0 
    2080       where (e3uw_0  (:,:,:).eq.0.0)  e3uw_0(:,:,:) = 1.0 
    2081       where (e3vw_0  (:,:,:).eq.0.0)  e3vw_0(:,:,:) = 1.0 
     2112      !fsdepw(:,:,:) = gdepw_0 (:,:,:) 
     2113      !fsde3w(:,:,:) = gdep3w_0(:,:,:) 
     2114      ! 
     2115      IF (.NOT. ln_wd) THEN 
     2116        where (e3t_0   (:,:,:).eq.0.0)  e3t_0(:,:,:) = 1.0 
     2117        where (e3u_0   (:,:,:).eq.0.0)  e3u_0(:,:,:) = 1.0 
     2118        where (e3v_0   (:,:,:).eq.0.0)  e3v_0(:,:,:) = 1.0 
     2119        where (e3f_0   (:,:,:).eq.0.0)  e3f_0(:,:,:) = 1.0 
     2120        where (e3w_0   (:,:,:).eq.0.0)  e3w_0(:,:,:) = 1.0 
     2121        where (e3uw_0  (:,:,:).eq.0.0)  e3uw_0(:,:,:) = 1.0 
     2122        where (e3vw_0  (:,:,:).eq.0.0)  e3vw_0(:,:,:) = 1.0 
     2123      END IF 
    20822124 
    20832125#if defined key_agrif 
     
    21192161         DO ji = 1, jpi 
    21202162            DO jk = 1, jpkm1 
     2163#if key_surge 
     2164               IF( scobot(ji,jj) >= fsdept(ji,jj,jk) )   mbathy(ji,jj) = MAX( 1, jk ) 
     2165#else 
    21212166               IF( scobot(ji,jj) >= fsdept(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk ) 
    2122             END DO 
    2123             IF( scobot(ji,jj) == 0._wp               )   mbathy(ji,jj) = 0 
     2167#endif 
     2168            END DO 
     2169            IF(ln_wd) THEN 
     2170              IF( scobot(ji,jj) <= -(rn_wdld - rn_wdmin2)) THEN 
     2171                mbathy(ji,jj) = 0 
     2172              ELSEIF(scobot(ji,jj) <= rn_wdmin1) THEN 
     2173#if key_surge 
     2174                mbathy(ji,jj) = 1 
     2175#else 
     2176                mbathy(ji,jj) = 2 
     2177#endif 
     2178              ENDIF 
     2179            ELSE 
     2180              IF( scobot(ji,jj) == 0._wp            )   mbathy(ji,jj) = 0 
     2181            END IF 
    21242182         END DO 
    21252183      END DO 
     
    21792237         DO jj = 1, jpj 
    21802238 
    2181             IF( hbatt(ji,jj) > 0._wp) THEN 
     2239            IF( bathy(ji,jj) > 0._wp) THEN 
    21822240               DO jk = 1, mbathy(ji,jj) 
    21832241                 ! check coordinate is monotonically increasing 
    2184                  IF (fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN 
    2185                     WRITE(ctmp1,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk 
    2186                     WRITE(numout,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk 
     2242                 !RF test... 
     2243                 !IF (fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN 
     2244                 !   WRITE(ctmp1,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk 
     2245                 !   WRITE(numout,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk 
     2246                 IF (fse3w(ji,jj,jk) < 0._wp .OR. fse3t(ji,jj,jk) < 0._wp ) THEN 
     2247                    WRITE(ctmp1,*) 'ERROR zgr_sco :   e3w   or e3t   < 0  at point (i,j,k)= ', ji, jj, jk 
     2248                    WRITE(numout,*) 'ERROR zgr_sco :   e3w   or e3t   < 0  at point (i,j,k)= ', ji, jj, jk 
    21872249                    WRITE(numout,*) 'e3w',fse3w(ji,jj,:) 
    21882250                    WRITE(numout,*) 'e3t',fse3t(ji,jj,:) 
     
    22432305      INTEGER  ::   ji, jj, jk           ! dummy loop argument 
    22442306      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars 
     2307      REAL(wp) ::   ztmpu,  ztmpv,  ztmpf 
     2308      REAL(wp) ::   ztmpu1, ztmpv1, ztmpf1 
    22452309      ! 
    22462310      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 
     
    22962360      DO ji = 1, jpim1 
    22972361         DO jj = 1, jpjm1 
     2362            ! extended for Wetting/Drying case 
     2363            ztmpu = hbatt(ji,jj)+hbatt(ji+1,jj) 
     2364            ztmpv = hbatt(ji,jj)+hbatt(ji,jj+1) 
     2365            ztmpf = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) 
     2366            ztmpu1 = hbatt(ji,jj)*hbatt(ji+1,jj) 
     2367            ztmpv1 = hbatt(ji,jj)*hbatt(ji,jj+1) 
     2368            ztmpf1 = MIN(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) * & 
     2369                   & MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) 
    22982370            DO jk = 1, jpk 
    2299                z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) )   & 
    2300                   &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
    2301                z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) )   & 
    2302                   &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
    2303                z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk)     & 
    2304                   &                + hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) )   & 
    2305                   &              / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 
    2306                z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) )   & 
    2307                   &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
    2308                z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) )   & 
    2309                   &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
     2371               IF((ztmpu1 < 0._wp.OR.ABS(ztmpu) < rn_wdmin1).AND.ln_wd) THEN 
     2372                 z_esigtu3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) ) 
     2373                 z_esigwu3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji+1,jj,jk) ) 
     2374               ELSE 
     2375                 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) & 
     2376                        &              / ztmpu 
     2377                 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) & 
     2378                        &              / ztmpu 
     2379               END IF 
     2380 
     2381               IF((ztmpv1 < 0._wp.OR.ABS(ztmpv) < rn_wdmin1).AND.ln_wd) THEN 
     2382                 z_esigtv3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji,jj+1,jk) ) 
     2383                 z_esigwv3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji,jj+1,jk) ) 
     2384               ELSE 
     2385                 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) & 
     2386                        &              / ztmpv 
     2387                 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) & 
     2388                        &              / ztmpv 
     2389               END IF 
     2390 
     2391               IF((ztmpf1 < 0._wp.OR.ABS(ztmpf) < rn_wdmin1).AND.ln_wd) THEN 
     2392                 z_esigtf3(ji,jj,jk) = 0.25_wp * ( z_esigt3(ji,jj,jk)   + z_esigt3(ji+1,jj,jk)  & 
     2393                    &                            + z_esigt3(ji,jj+1,jk) + z_esigt3(ji+1,jj+1,jk) ) 
     2394               ELSE 
     2395                 z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk)  & 
     2396                    &                + hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) & 
     2397                    &              / ztmpf 
     2398               END IF 
     2399 
    23102400               ! 
    23112401               e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 
     
    23472437      REAL(wp) ::   zsmth               ! smoothing around critical depth 
    23482438      REAL(wp) ::   zzs, zzb           ! Surface and bottom cell thickness in sigma space 
     2439      REAL(wp) ::   ztmpu, ztmpv, ztmpf 
     2440      REAL(wp) ::   ztmpu1, ztmpv1, ztmpf1 
    23492441      ! 
    23502442      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 
     
    24252517        DO jj=1,jpj-1 
    24262518 
    2427           DO jk = 1, jpk 
    2428                 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) / & 
    2429                                     ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
    2430                 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) / & 
    2431                                     ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
    2432                 z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) +  & 
    2433                                       hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / & 
    2434                                     ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 
    2435                 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) / & 
    2436                                     ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 
    2437                 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) / & 
    2438                                     ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 
     2519           ! extend to suit for Wetting/Drying case 
     2520           ztmpu = hbatt(ji,jj)+hbatt(ji+1,jj) 
     2521           ztmpv = hbatt(ji,jj)+hbatt(ji,jj+1) 
     2522           ztmpf = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) 
     2523           ztmpu1 = hbatt(ji,jj)*hbatt(ji+1,jj) 
     2524           ztmpv1 = hbatt(ji,jj)*hbatt(ji,jj+1) 
     2525           ztmpf1 = MIN(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) * & 
     2526                  & MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) 
     2527           DO jk = 1, jpk 
     2528              IF((ztmpu1 < 0._wp.OR.ABS(ztmpu) < rn_wdmin1).AND.ln_wd) THEN 
     2529                z_esigtu3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) ) 
     2530                z_esigwu3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji+1,jj,jk) ) 
     2531              ELSE 
     2532                z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) & 
     2533                       &              / ztmpu 
     2534                z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) & 
     2535                       &              / ztmpu 
     2536              END IF 
     2537 
     2538              IF((ztmpv1 < 0._wp.OR.ABS(ztmpv) < rn_wdmin1).AND.ln_wd) THEN 
     2539                z_esigtv3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji,jj+1,jk) ) 
     2540                z_esigwv3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji,jj+1,jk) ) 
     2541              ELSE 
     2542                z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) & 
     2543                       &              / ztmpv 
     2544                z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) & 
     2545                       &              / ztmpv 
     2546              END IF 
     2547 
     2548              IF((ztmpf1 < 0._wp.OR.ABS(ztmpf) < rn_wdmin1).AND.ln_wd) THEN 
     2549                z_esigtf3(ji,jj,jk) = 0.25_wp * ( z_esigt3(ji,jj,jk)   + z_esigt3(ji+1,jj,jk)  & 
     2550                   &                            + z_esigt3(ji,jj+1,jk) + z_esigt3(ji+1,jj+1,jk) ) 
     2551              ELSE 
     2552                z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk)  & 
     2553                   &                + hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) & 
     2554                   &              / ztmpf 
     2555              END IF 
    24392556 
    24402557             e3t_0(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk) 
  • branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90

    r5332 r5954  
    204204      INTEGER  :: ji, jj, jk 
    205205      REAL(wp) ::   zsal = 35.50 
     206#if defined key_surge 
     207      REAL(wp) ::   ztem = 10.0 
     208#endif 
    206209      !!---------------------------------------------------------------------- 
    207210      ! 
     
    210213      IF(lwp) WRITE(numout,*) '~~~~~~~~~~   and constant salinity (',zsal,' psu)' 
    211214      ! 
     215#if defined key_surge 
     216      tsn(:,:,:,jp_tem) = ztem * tmask(:,:,:) 
     217      tsb(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) 
     218#else 
    212219      DO jk = 1, jpk 
    213220         tsn(:,:,jk,jp_tem) = (  ( ( 7.5 - 0. * ABS( gphit(:,:) )/30. ) * ( 1.-TANH((fsdept(:,:,jk)-80.)/30.) )   & 
     
    215222         tsb(:,:,jk,jp_tem) = tsn(:,:,jk,jp_tem) 
    216223      END DO 
     224#endif 
    217225      tsn(:,:,:,jp_sal) = zsal * tmask(:,:,:) 
    218226      tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) 
  • branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90

    r5224 r5954  
    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   !!            3.6  !  2014-11  (P. Mathiot) hpg_isf: original code for ice shelf cavity 
    1920   !!---------------------------------------------------------------------- 
     
    379380      !! 
    380381      INTEGER  ::   ji, jj, jk                 ! dummy loop indices 
    381       REAL(wp) ::   zcoef0, zuap, zvap, znad   ! temporary scalars 
     382      REAL(wp) ::   zcoef0, zuap, zvap, znad, ztmp   ! temporary scalars 
     383      LOGICAL  ::   ll_tmp1, ll_tmp2           ! local logical variables 
    382384      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zhpi, zhpj 
     385      REAL(wp), POINTER, DIMENSION(:,:)   ::  zcpx, zcpy    !W/D pressure filter 
    383386      !!---------------------------------------------------------------------- 
    384387      ! 
    385388      CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj ) 
     389      IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 
    386390      ! 
    387391      IF( kt == nit000 ) THEN 
     
    398402      ENDIF 
    399403 
     404      IF(ln_wd) THEN 
     405        DO jj = 2, jpjm1 
     406           DO ji = 2, jpim1 
     407             IF ( tmask(ji+1,jj,1) == 0._wp) THEN 
     408                zcpx(ji,jj) = 1.0_wp 
     409             ELSE 
     410                ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) 
     411                ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) +& 
     412                                                          & rn_wdmin1 + rn_wdmin2 
     413    
     414                IF(ll_tmp1) THEN 
     415                  zcpx(ji,jj) = 1.0_wp 
     416                ELSE IF(ll_tmp2) THEN 
     417                   ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 
     418                   zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /& 
     419                               &     (sshn(ji+1,jj) - sshn(ji,jj))) 
     420                ELSE 
     421                  zcpx(ji,jj) = 0._wp 
     422                END IF 
     423             ENDIF 
     424         
     425             IF ( tmask(ji,jj+1,1) == 0._wp) THEN 
     426                zcpy(ji,jj) = 1.0_wp 
     427             ELSE 
     428                ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) 
     429                ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) +& 
     430                                                       & rn_wdmin1 + rn_wdmin2 
     431 
     432                IF(ll_tmp1) THEN 
     433                  zcpy(ji,jj) = 1.0_wp 
     434                ELSE IF(ll_tmp2) THEN 
     435                  ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 
     436                  zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) /& 
     437                              &     (sshn(ji,jj+1) - sshn(ji,jj))) 
     438                ELSE 
     439                  zcpy(ji,jj) = 0._wp 
     440                END IF 
     441             ENDIF 
     442           END DO 
     443        END DO 
     444        CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp ) 
     445      ENDIF 
     446 
     447#if defined key_surge 
     448      ! Surface value 
     449      DO jj = 2, jpjm1 
     450         DO ji = fs_2, fs_jpim1   ! vector opt. 
     451            ! hydrostatic pressure gradient along s-surfaces 
     452            zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( fse3w(ji+1,jj  ,1) * ( znad )   & 
     453               &                                  - fse3w(ji  ,jj  ,1) * ( znad ) ) 
     454            zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( fse3w(ji  ,jj+1,1) * ( znad )   & 
     455               &                                  - fse3w(ji  ,jj  ,1) * ( znad ) ) 
     456            ! s-coordinate pressure gradient correction 
     457            zuap = -zcoef0 * ( 2._wp * znad )   & 
     458               &           * ( fsde3w(ji+1,jj,1) - fsde3w(ji,jj,1) ) / e1u(ji,jj) 
     459            zvap = -zcoef0 * ( 2._wp * znad )   & 
     460               &           * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj) 
     461 
     462            IF(ln_wd) THEN 
     463              zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 
     464              zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 
     465              zuap = zuap * zcpx(ji,jj) 
     466              zvap = zvap * zcpy(ji,jj) 
     467            ENDIF 
     468 
     469            ! add to the general momentum trend 
     470            ua(ji,jj,1) = ua(ji,jj,1) + ( zhpi(ji,jj,1) + zuap ) * umask(ji,jj,1) 
     471            va(ji,jj,1) = va(ji,jj,1) + ( zhpj(ji,jj,1) + zvap ) * vmask(ji,jj,1) 
     472         END DO 
     473      END DO 
     474 
     475      ! interior value (2=<jk=<jpkm1) 
     476      DO jk = 2, jpkm1 
     477         DO jj = 2, jpjm1 
     478            DO ji = fs_2, fs_jpim1   ! vector opt. 
     479               ! hydrostatic pressure gradient along s-surfaces 
     480               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj)   & 
     481                  &           * (  fse3w(ji+1,jj,jk) * ( 2*znad )   & 
     482                  &              - fse3w(ji  ,jj,jk) * ( 2*znad )  ) 
     483               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj)   & 
     484                  &           * (  fse3w(ji,jj+1,jk) * ( 2*znad )   & 
     485                  &              - fse3w(ji,jj  ,jk) * ( 2*znad )  ) 
     486               ! s-coordinate pressure gradient correction 
     487               zuap = -zcoef0 * ( 2._wp * znad )   & 
     488                  &           * ( fsde3w(ji+1,jj  ,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj) 
     489               zvap = -zcoef0 * ( 2._wp * znad )   & 
     490                  &           * ( fsde3w(ji  ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) 
     491 
     492               IF(ln_wd) THEN 
     493                 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 
     494                 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) 
     495                 zuap = zuap * zcpx(ji,jj) 
     496                 zvap = zvap * zcpy(ji,jj) 
     497               ENDIF 
     498 
     499               ! add to the general momentum trend 
     500               ua(ji,jj,jk) = ua(ji,jj,jk) + ( zhpi(ji,jj,jk) + zuap ) * umask(ji,jj,jk) 
     501               va(ji,jj,jk) = va(ji,jj,jk) + ( zhpj(ji,jj,jk) + zvap ) * vmask(ji,jj,jk) 
     502            END DO 
     503         END DO 
     504      END DO 
     505      ! 
     506#else 
    400507      ! Surface value 
    401508      DO jj = 2, jpjm1 
     
    411518            zvap = -zcoef0 * ( rhd   (ji,jj+1,1) + rhd   (ji,jj,1) + 2._wp * znad )   & 
    412519               &           * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj) 
     520 
     521            IF(ln_wd) THEN 
     522              zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 
     523              zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 
     524              zuap = zuap * zcpx(ji,jj) 
     525              zvap = zvap * zcpy(ji,jj) 
     526            ENDIF 
     527 
    413528            ! add to the general momentum trend 
    414529            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap 
     
    433548               zvap = -zcoef0 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   & 
    434549                  &           * ( fsde3w(ji  ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) 
     550 
     551               IF(ln_wd) THEN 
     552                 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 
     553                 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) 
     554                 zuap = zuap * zcpx(ji,jj) 
     555                 zvap = zvap * zcpy(ji,jj) 
     556               ENDIF 
     557 
    435558               ! add to the general momentum trend 
    436559               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap 
     
    439562         END DO 
    440563      END DO 
    441       ! 
     564#endif 
     565      ! 
     566 
    442567      CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj ) 
     568      IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 
    443569      ! 
    444570   END SUBROUTINE hpg_sco 
     
    719845      REAL(wp) ::   z1_10, cffu, cffx   !    "         " 
    720846      REAL(wp) ::   z1_12, cffv, cffy   !    "         " 
     847      LOGICAL  ::   ll_tmp1, ll_tmp2    ! local logical variables 
    721848      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zhpi, zhpj 
    722849      REAL(wp), POINTER, DIMENSION(:,:,:) ::  dzx, dzy, dzz, dzu, dzv, dzw 
    723850      REAL(wp), POINTER, DIMENSION(:,:,:) ::  drhox, drhoy, drhoz, drhou, drhov, drhow 
    724851      REAL(wp), POINTER, DIMENSION(:,:,:) ::  rho_i, rho_j, rho_k 
     852      REAL(wp), POINTER, DIMENSION(:,:)   ::  zcpx, zcpy    !W/D pressure filter 
    725853      !!---------------------------------------------------------------------- 
    726854      ! 
     
    728856      CALL wrk_alloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow ) 
    729857      CALL wrk_alloc( jpi, jpj, jpk, rho_i, rho_j, rho_k,  zhpi,  zhpj        ) 
     858      IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 
    730859      ! 
    731860 
     
    734863         IF(lwp) WRITE(numout,*) 'dyn:hpg_djc : hydrostatic pressure gradient trend' 
    735864         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate case, density Jacobian with cubic polynomial scheme' 
     865      ENDIF 
     866 
     867      IF(ln_wd) THEN 
     868        DO jj = 2, jpjm1 
     869           DO ji = 2, jpim1 
     870             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) 
     871             ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) +& 
     872                                                       & rn_wdmin1 + rn_wdmin2 
     873 
     874             IF(ll_tmp1) THEN 
     875               zcpx(ji,jj) = 1.0_wp 
     876             ELSE IF(ll_tmp2) THEN 
     877               ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 
     878               zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /& 
     879                           &     (sshn(ji+1,jj) - sshn(ji,jj))) 
     880             ELSE 
     881               zcpx(ji,jj) = 0._wp 
     882             END IF 
     883      
     884             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) 
     885             ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) +& 
     886                                                       & rn_wdmin1 + rn_wdmin2 
     887 
     888             IF(ll_tmp1) THEN 
     889               zcpy(ji,jj) = 1.0_wp 
     890             ELSE IF(ll_tmp2) THEN 
     891               ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 
     892               zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) /& 
     893                           &     (sshn(ji,jj+1) - sshn(ji,jj))) 
     894             ELSE 
     895               zcpy(ji,jj) = 0._wp 
     896             END IF 
     897           END DO 
     898        END DO 
     899        CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp ) 
    736900      ENDIF 
    737901 
     
    8991063            zhpi(ji,jj,1) = ( rho_k(ji+1,jj  ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) / e1u(ji,jj) 
    9001064            zhpj(ji,jj,1) = ( rho_k(ji  ,jj+1,1) - rho_k(ji,jj,1) - rho_j(ji,jj,1) ) / e2v(ji,jj) 
     1065            IF(ln_wd) THEN 
     1066              zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 
     1067              zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 
     1068            ENDIF 
    9011069            ! add to the general momentum trend 
    9021070            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) 
     
    9181086                  &           + (  ( rho_k(ji,jj+1,jk) - rho_k(ji,jj,jk  ) )    & 
    9191087                  &               -( rho_j(ji,jj  ,jk) - rho_j(ji,jj,jk-1) )  ) / e2v(ji,jj) 
     1088               IF(ln_wd) THEN 
     1089                 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 
     1090                 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) 
     1091               ENDIF 
    9201092               ! add to the general momentum trend 
    9211093               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) 
     
    9281100      CALL wrk_dealloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow ) 
    9291101      CALL wrk_dealloc( jpi, jpj, jpk, rho_i, rho_j, rho_k,  zhpi,  zhpj        ) 
     1102      IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 
    9301103      ! 
    9311104   END SUBROUTINE hpg_djc 
     
    9511124      !! The local variables for the correction term 
    9521125      INTEGER  :: jk1, jis, jid, jjs, jjd 
     1126      LOGICAL  :: ll_tmp1, ll_tmp2                  ! local logical variables 
    9531127      REAL(wp) :: zuijk, zvijk, zpwes, zpwed, zpnss, zpnsd, zdeps 
    9541128      REAL(wp) :: zrhdt1 
     
    9571131      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp 
    9581132      REAL(wp), POINTER, DIMENSION(:,:)   ::   zsshu_n, zsshv_n 
     1133      REAL(wp), POINTER, DIMENSION(:,:)   ::   zcpx, zcpy    !W/D pressure filter 
    9591134      !!---------------------------------------------------------------------- 
    9601135      ! 
     
    9621137      CALL wrk_alloc( jpi,jpj,jpk, zdept, zrhh ) 
    9631138      CALL wrk_alloc( jpi,jpj, zsshu_n, zsshv_n ) 
     1139      IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy ) 
    9641140      ! 
    9651141      IF( kt == nit000 ) THEN 
     
    9741150      znad = 0.0_wp 
    9751151      IF( lk_vvl ) znad = 1._wp 
     1152 
     1153      IF(ln_wd) THEN 
     1154        DO jj = 2, jpjm1 
     1155           DO ji = 2, jpim1 
     1156             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) 
     1157             ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) +& 
     1158                                                       & rn_wdmin1 + rn_wdmin2 
     1159 
     1160             IF(ll_tmp1) THEN 
     1161               zcpx(ji,jj) = 1.0_wp 
     1162             ELSE IF(ll_tmp2) THEN 
     1163               ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 
     1164               zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /& 
     1165                           &     (sshn(ji+1,jj) - sshn(ji,jj))) 
     1166             ELSE 
     1167               zcpx(ji,jj) = 0._wp 
     1168             END IF 
     1169      
     1170             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) 
     1171             ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) +& 
     1172                                                       & rn_wdmin1 + rn_wdmin2 
     1173 
     1174             IF(ll_tmp1.OR.ll_tmp2) THEN 
     1175               zcpy(ji,jj) = 1.0_wp 
     1176             ELSE IF(ll_tmp2) THEN 
     1177               ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 
     1178               zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) /& 
     1179                           &     (sshn(ji,jj+1) - sshn(ji,jj))) 
     1180             ELSE 
     1181               zcpy(ji,jj) = 0._wp 
     1182             END IF 
     1183           END DO 
     1184        END DO 
     1185        CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp ) 
     1186      ENDIF 
    9761187 
    9771188      ! Clean 3-D work arrays 
     
    10521263        END DO 
    10531264      END DO 
     1265 
     1266      CALL lbc_lnk (zsshu_n, 'U', 1.) 
     1267      CALL lbc_lnk (zsshv_n, 'V', 1.) 
    10541268 
    10551269      DO jj = 2, jpjm1 
     
    11501364               ENDIF 
    11511365 
    1152                ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * & 
    1153                &           umask(ji,jj,jk) * tmask(ji,jj,jk) * tmask(ji+1,jj,jk) 
     1366               IF(ln_wd) THEN 
     1367                  zdpdx1 = zdpdx1 * zcpx(ji,jj) 
     1368                  zdpdx2 = zdpdx2 * zcpx(ji,jj) 
     1369                ENDIF 
     1370                ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk)  
    11541371            ENDIF 
    11551372 
     
    12071424               ENDIF 
    12081425 
    1209                va(ji,jj,jk) = va(ji,jj,jk) + (zdpdy1 + zdpdy2)*& 
    1210                &              vmask(ji,jj,jk)*tmask(ji,jj,jk)*tmask(ji,jj+1,jk) 
     1426               IF(ln_wd) THEN 
     1427                  zdpdy1 = zdpdy1 * zcpy(ji,jj) 
     1428                  zdpdy2 = zdpdy2 * zcpy(ji,jj) 
     1429               ENDIF 
     1430 
     1431               va(ji,jj,jk) = va(ji,jj,jk) + (zdpdy1 + zdpdy2) * vmask(ji,jj,jk) 
    12111432            ENDIF 
    12121433 
     
    12191440      CALL wrk_dealloc( jpi,jpj,jpk, zdept, zrhh ) 
    12201441      CALL wrk_dealloc( jpi,jpj, zsshu_n, zsshv_n ) 
     1442      IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy ) 
    12211443      ! 
    12221444   END SUBROUTINE hpg_prj 
  • branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90

    r5328 r5954  
    137137         DO jj = 2, jpjm1 
    138138            DO ji = fs_2, fs_jpim1   ! vector opt. 
    139                ua(ji,jj,jk) = ua(ji,jj,jk) - ( zhke(ji+1,jj  ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj) 
    140                va(ji,jj,jk) = va(ji,jj,jk) - ( zhke(ji  ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj) 
     139               ua(ji,jj,jk) = ua(ji,jj,jk) - ( zhke(ji+1,jj  ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj) * umask(ji,jj,jk) 
     140               va(ji,jj,jk) = va(ji,jj,jk) - ( zhke(ji  ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj) * vmask(ji,jj,jk) 
    141141            END DO  
    142142         END DO 
  • branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilap.F90

    r4990 r5954  
    140140          
    141141         ! Multiply by the eddy viscosity coef. (at u- and v-points) 
    142          zlu(:,:,jk) = zlu(:,:,jk) * ( fsahmu(:,:,jk) * (1-nkahm_smag) + nkahm_smag) 
    143  
    144          zlv(:,:,jk) = zlv(:,:,jk) * ( fsahmv(:,:,jk) * (1-nkahm_smag) + nkahm_smag) 
     142         zlu(:,:,jk) = zlu(:,:,jk) * ( fsahmu(:,:,jk) * (1-nkahm_smag) + nkahm_smag) * umask(:,:,jk) 
     143  
     144         zlv(:,:,jk) = zlv(:,:,jk) * ( fsahmv(:,:,jk) * (1-nkahm_smag) + nkahm_smag) * vmask(:,:,jk) 
    145145          
    146146         ! Contravariant "laplacian" 
     
    170170               zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) 
    171171               zut(ji,jj,jk) = (  zlu(ji,jj,jk) - zlu(ji-1,jj  ,jk)   & 
    172                   &             + zlv(ji,jj,jk) - zlv(ji  ,jj-1,jk) ) / zbt 
     172                  &             + zlv(ji,jj,jk) - zlv(ji  ,jj-1,jk) ) / zbt * tmask(ji,jj,jk) 
    173173            END DO 
    174174         END DO 
     
    197197                  &  + ( zut(ji,jj+1,jk) - zut(ji  ,jj,jk) ) / e2v(ji,jj) 
    198198               ! add it to the general momentum trends 
    199                ua(ji,jj,jk) = ua(ji,jj,jk) + zua * ( fsahmu(ji,jj,jk)*nkahm_smag +(1 -nkahm_smag )) 
    200                va(ji,jj,jk) = va(ji,jj,jk) + zva * ( fsahmv(ji,jj,jk)*nkahm_smag +(1 -nkahm_smag )) 
     199               ua(ji,jj,jk) = ua(ji,jj,jk) + zua * ( fsahmu(ji,jj,jk)*nkahm_smag +(1 -nkahm_smag )) * umask(ji,jj,jk) 
     200               va(ji,jj,jk) = va(ji,jj,jk) + zva * ( fsahmv(ji,jj,jk)*nkahm_smag +(1 -nkahm_smag )) * vmask(ji,jj,jk) 
    201201            END DO 
    202202         END DO 
  • branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90

    r5120 r5954  
    160160            DO jj = 2, jpjm1 
    161161               DO ji = fs_2, fs_jpim1   ! vector opt. 
    162                   ua(ji,jj,jk) = ua(ji,jj,jk) + spgu(ji,jj) 
    163                   va(ji,jj,jk) = va(ji,jj,jk) + spgv(ji,jj) 
     162                  ua(ji,jj,jk) = ua(ji,jj,jk) + spgu(ji,jj) * umask(ji,jj,jk) 
     163                  va(ji,jj,jk) = va(ji,jj,jk) + spgv(ji,jj) * vmask(ji,jj,jk) 
    164164               END DO 
    165165            END DO 
  • branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r5942 r5954  
    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 
     
    4142   USE timing          ! Timing     
    4243   USE sbcapr          ! surface boundary condition: atmospheric pressure 
     44   USE wadlmt          ! wetting/drying flux limter 
    4345   USE dynadv, ONLY: ln_dynadv_vec 
    4446#if defined key_agrif 
     
    140142      INTEGER, INTENT(in)  ::   kt   ! ocean time-step index 
    141143      ! 
    142       LOGICAL  ::   ll_fw_start        ! if true, forward integration  
    143       LOGICAL  ::   ll_init             ! if true, special startup of 2d equations 
    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           !   -      - 
     144      LOGICAL  ::   ll_fw_start                 ! if true, forward integration 
     145      LOGICAL  ::   ll_init                         ! if true, special startup of 2d equations 
     146      LOGICAL  ::   ll_tmp1, ll_tmp2            ! local logical variables used in W/D 
     147      INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices 
     148      INTEGER  ::   ikbu, ikbv, noffset         ! local integers 
     149      REAL(wp) ::   zraur, z1_2dt_b, z2dt_bf    ! local scalars 
     150      REAL(wp) ::   zx1, zy1, zx2, zy2          !   -      - 
     151      REAL(wp) ::   z1_12, z1_8, z1_4, z1_2     !   -      - 
     152      REAL(wp) ::   zu_spg, zv_spg              !   -      - 
     153      REAL(wp) ::   zhura, zhvra                !   -      - 
     154      REAL(wp) ::   za0, za1, za2, za3          !   -      - 
     155 
     156      REAL(wp) ::   ztmp                       ! temporary vaialbe 
    152157      ! 
    153158      REAL(wp), POINTER, DIMENSION(:,:) :: zun_e, zvn_e, zsshp2_e 
     
    157162      REAL(wp), POINTER, DIMENSION(:,:) :: zsshu_a, zsshv_a 
    158163      REAL(wp), POINTER, DIMENSION(:,:) :: zhf 
     164      REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy                 ! Wetting/Dying gravity filter coef. 
    159165      !!---------------------------------------------------------------------- 
    160166      ! 
     
    168174      CALL wrk_alloc( jpi, jpj, zsshu_a, zsshv_a                                   ) 
    169175      CALL wrk_alloc( jpi, jpj, zhf ) 
     176      IF(ln_wd) CALL wrk_alloc( jpi, jpj, zcpx, zcpy ) 
    170177      ! 
    171178      !                                         !* Local constant initialization 
     
    176183      zraur = 1._wp / rau0 
    177184      ! 
    178       IF( kt == nit000 .AND. neuler == 0 ) THEN    ! reciprocal of baroclinic time step  
     185      IF( kt == nit000 .AND. neuler == 0 ) THEN  ! reciprocal of baroclinic time step  
    179186        z2dt_bf = rdt 
    180187      ELSE 
     
    183190      z1_2dt_b = 1.0_wp / z2dt_bf  
    184191      ! 
    185       ll_init = ln_bt_av                           ! if no time averaging, then no specific restart  
     192      ll_init = ln_bt_av                       ! if no time averaging, then no specific restart  
    186193      ll_fw_start = .FALSE. 
    187194      ! 
    188                                                        ! time offset in steps for bdy data update 
    189       IF (.NOT.ln_bt_fw) THEN ; noffset=-nn_baro ; ELSE ;  noffset = 0 ; ENDIF 
    190       ! 
    191       IF( kt == nit000 ) THEN                !* initialisation 
     195                                                     ! time offset in steps for bdy data update 
     196      IF (.NOT.ln_bt_fw) THEN ; noffset=-nn_baro ; ELSE ;  noffset = 0 ; ENDIF  
     197      ! 
     198      IF( kt == nit000 ) THEN             !* initialisation 
    192199         ! 
    193200         IF(lwp) WRITE(numout,*) 
     
    378385      !                                   !* Right-Hand-Side of the barotropic momentum equation 
    379386      !                                   ! ---------------------------------------------------- 
     387 
    380388      IF( lk_vvl ) THEN                         ! Variable volume : remove surface pressure gradient 
    381          DO jj = 2, jpjm1  
    382             DO ji = fs_2, fs_jpim1   ! vector opt. 
    383                zu_trd(ji,jj) = zu_trd(ji,jj) - grav * (  sshn(ji+1,jj  ) - sshn(ji  ,jj  )  ) / e1u(ji,jj) 
    384                zv_trd(ji,jj) = zv_trd(ji,jj) - grav * (  sshn(ji  ,jj+1) - sshn(ji  ,jj  )  ) / e2v(ji,jj) 
    385             END DO 
    386          END DO 
     389        IF(ln_wd) THEN                    ! Calculating and applying W/D gravity filters 
     390          DO jj = 2, jpjm1 
     391             DO ji = 2, jpim1 
     392                IF ( tmask(ji+1,jj,1) == 0._wp ) THEN 
     393                   zcpx = 1._wp 
     394                ELSE 
     395                   ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj)*tmask(ji,jj,1), -bathy(ji+1,jj)*tmask(ji+1,jj,1)) & 
     396                           & .and. MAX(sshn(ji,jj) + bathy(ji,jj)*tmask(ji,jj,1), sshn(ji+1,jj) + bathy(ji+1,jj)*tmask(ji+1,jj,1)) & 
     397                           &  > rn_wdmin1 + rn_wdmin2 
     398                   ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj)*tmask(ji,jj,1), -bathy(ji+1,jj)*tmask(ji+1,jj,1)) + & 
     399                                                          & rn_wdmin1 + rn_wdmin2 
     400                   IF(ll_tmp1) THEN 
     401                      zcpx(ji,jj) = 1.0_wp 
     402                   ELSE IF(ll_tmp2) THEN 
     403                      ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 
     404                      zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj)*tmask(ji+1,jj,1) - sshn(ji,jj) - bathy(ji,jj)*tmask(ji,jj,1)) /& 
     405                                  &     (sshn(ji+1,jj) - sshn(ji,jj))) 
     406                   ELSE 
     407                      zcpx(ji,jj) = 0._wp 
     408                   END IF 
     409                ENDIF 
     410 
     411                IF ( tmask(ji,jj+1,1) == 0._wp ) THEN 
     412                   zcpy = 1._wp 
     413                ELSE 
     414                   ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj)*tmask(ji,jj,1), -bathy(ji,jj+1)*tmask(ji,jj+1,1)) & 
     415                           & .and. MAX(sshn(ji,jj) + bathy(ji,jj)*tmask(ji,jj,1), sshn(ji,jj+1) + bathy(ji,jj+1)*tmask(ji,jj+1,1)) & 
     416                           &  > rn_wdmin1 + rn_wdmin2 
     417                   ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj)*tmask(ji,jj,1), -bathy(ji,jj+1)*tmask(ji,jj+1,1)) + & 
     418                                                             & rn_wdmin1 + rn_wdmin2 
     419                   IF(ll_tmp1) THEN 
     420                      zcpy(ji,jj) = 1.0_wp 
     421                   ELSE IF(ll_tmp2) THEN 
     422                      ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 
     423                     zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1)*tmask(ji,jj+1,1) - sshn(ji,jj) - bathy(ji,jj)*tmask(ji,jj,1)) /& 
     424                                 &     (sshn(ji,jj+1) - sshn(ji,jj))) 
     425                   ELSE 
     426                     zcpy(ji,jj) = 0._wp 
     427                   END IF 
     428                ENDIF 
     429             END DO 
     430           END DO 
     431           CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp ) 
     432        ENDIF 
     433 
     434        IF(ln_wd) THEN                    ! Calculating and applying W/D gravity filters 
     435           DO jj = 2, jpjm1 
     436              DO ji = 2, jpim1 
     437                 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj  ) - sshn(ji  ,jj ) ) / & 
     438                               &                 e1u(ji,jj) * zcpx(ji,jj) 
     439                 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji  ,jj+1) - sshn(ji  ,jj ) ) / & 
     440                               &                 e2v(ji,jj) * zcpy(ji,jj) 
     441              END DO 
     442           END DO 
     443         ELSE 
     444           DO jj = 2, jpjm1 
     445              DO ji = fs_2, fs_jpim1   ! vector opt. 
     446                 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * (  sshn(ji+1,jj  ) - sshn(ji  ,jj  )  ) / e1u(ji,jj) 
     447                 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * (  sshn(ji  ,jj+1) - sshn(ji  ,jj  )  ) / e2v(ji,jj) 
     448              END DO 
     449           END DO 
     450        END IF 
    387451      ENDIF 
    388452 
     
    416480      ! 
    417481      ! Note that the "unclipped" bottom friction parameter is used even with explicit drag 
    418       zu_frc(:,:) = zu_frc(:,:) + hur(:,:) * bfrua(:,:) * zwx(:,:) 
    419       zv_frc(:,:) = zv_frc(:,:) + hvr(:,:) * bfrva(:,:) * zwy(:,:) 
     482      IF(ln_wd) THEN 
     483        zu_frc(:,:) = zu_frc(:,:) + MAX(hur(:,:) * bfrua(:,:),-1._wp / rdtbt) * zwx(:,:) 
     484        zv_frc(:,:) = zv_frc(:,:) + MAX(hvr(:,:) * bfrva(:,:),-1._wp / rdtbt) * zwy(:,:) 
     485      ELSE 
     486        zu_frc(:,:) = zu_frc(:,:) + hur(:,:) * bfrua(:,:) * zwx(:,:) 
     487        zv_frc(:,:) = zv_frc(:,:) + hvr(:,:) * bfrva(:,:) * zwy(:,:) 
     488      END IF 
    420489      !        
    421490      IF (ln_bt_fw) THEN                        ! Add wind forcing 
     
    488557      ENDIF 
    489558      ! 
     559      IF(ln_wd) THEN      !preserve the positivity of water depth 
     560                          !ssh[b,n,a] should have already been processed for this 
     561         sshbb_e(:,:jj) = MAX( sshbb_e(:,:), (rn_wdmin1 - bathy(:,:)) ) *tmask(:,:,1) 
     562         sshb_e(:,:jj)  = MAX( sshb_e(:,:) , (rn_wdmin1 - bathy(:,:)) ) *tmask(:,:,1) 
     563      ENDIF 
     564 
    490565      IF (ln_bt_fw) THEN                  ! FORWARD integration: start from NOW fields                     
    491566         sshn_e(:,:) = sshn (:,:)             
     
    523598         ! Update only tidal forcing at open boundaries 
    524599#if defined key_tide 
    525          IF ( lk_bdy .AND. lk_tide ) CALL bdy_dta_tides( kt, kit=jn, time_offset=(noffset+1) ) 
    526          IF ( ln_tide_pot .AND. lk_tide ) CALL upd_tide( kt, kit=jn, time_offset=noffset ) 
     600         IF ( lk_bdy .AND. lk_tide ) CALL bdy_dta_tides( kt, kit=jn, time_offset=(noffset+1) )  
     601         IF ( ln_tide_pot .AND. lk_tide ) CALL upd_tide( kt, kit=jn, time_offset=noffset )  
    527602#endif 
    528603         ! 
     
    561636            zhup2_e (:,:) = hu_0(:,:) + zwx(:,:)                ! Ocean depth at U- and V-points 
    562637            zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:) 
     638            IF(ln_wd) THEN 
     639              zhup2_e(:,:) = MAX(zhup2_e (:,:), rn_wdmin1) 
     640              zhvp2_e(:,:) = MAX(zhvp2_e (:,:), rn_wdmin1) 
     641            END IF 
    563642         ELSE 
    564643            zhup2_e (:,:) = hu(:,:) 
     
    599678         ENDIF 
    600679#endif 
     680 
     681         IF(ln_wd) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt) 
    601682         ! 
    602683         ! Sum over sub-time-steps to compute advective velocities 
     
    613694         END DO 
    614695         ssha_e(:,:) = (  sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) )  ) * tmask(:,:,1) 
     696         IF(ln_wd) ssha_e(:,:) = MAX(ssha_e(:,:), (rn_wdmin1 - bathy(:,:)) ) * tmask(:,:,1) 
    615697         CALL lbc_lnk( ssha_e, 'T',  1._wp ) 
    616698 
     
    676758               END DO 
    677759            END DO 
     760 
     761            IF(ln_wd) THEN 
     762              zhust_e(:,:) = MAX(zhust_e (:,:), rn_wdmin1 ) 
     763              zhvst_e(:,:) = MAX(zhvst_e (:,:), rn_wdmin1 ) 
     764            END IF 
     765            !shall we call lbc_lnk for zhu(v)st_e() here? 
     766 
    678767         ENDIF 
    679768         ! 
     
    738827         ! 
    739828         ! Add bottom stresses: 
    740          zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * zun_e(:,:) * hur_e(:,:) 
    741          zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * zvn_e(:,:) * hvr_e(:,:) 
    742          ! 
    743          ! Surface pressure trend: 
    744          DO jj = 2, jpjm1 
    745             DO ji = fs_2, fs_jpim1   ! vector opt. 
    746                ! Add surface pressure gradient 
    747                zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj) 
    748                zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj) 
    749                zwx(ji,jj) = zu_spg 
    750                zwy(ji,jj) = zv_spg 
    751             END DO 
    752          END DO 
     829         IF(ln_wd) THEN 
     830           zu_trd(:,:) = zu_trd(:,:) + MAX(bfrua(:,:) * hur_e(:,:), -1._wp / rdtbt) * zun_e(:,:) 
     831           zv_trd(:,:) = zv_trd(:,:) + MAX(bfrva(:,:) * hvr_e(:,:), -1._wp / rdtbt) * zvn_e(:,:) 
     832         ELSE 
     833           zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * zun_e(:,:) * hur_e(:,:) 
     834           zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * zvn_e(:,:) * hvr_e(:,:) 
     835         END IF 
     836          
     837         IF(ln_wd) THEN                    ! Calculating and applying W/D gravity filters 
     838           DO jj = 2, jpjm1 
     839              DO ji = 2, jpim1 
     840                 IF ( tmask(ji+1,jj,1) == 0._wp ) THEN 
     841                    zcpx = 1._wp 
     842                 ELSE 
     843                    ll_tmp1 = MIN(zsshp2_e(ji,jj), zsshp2_e(ji+1,jj)) > MAX(-bathy(ji,jj)*tmask(ji,jj,1), -bathy(ji+1,jj)*tmask(ji+1,jj,1)) & 
     844                           & .and. MAX(zsshp2_e(ji,jj) + bathy(ji,jj)*tmask(ji,jj,1), zsshp2_e(ji+1,jj) + bathy(ji+1,jj)*tmask(ji+1,jj,1))  & 
     845                           &  > rn_wdmin1 + rn_wdmin2 
     846                    ll_tmp2 = MAX(zsshp2_e(ji,jj), zsshp2_e(ji+1,jj)) > MAX(-bathy(ji,jj)*tmask(ji,jj,1), -bathy(ji+1,jj)*tmask(ji+1,jj,1)) + & 
     847                                                           & rn_wdmin1 + rn_wdmin2 
     848                    IF(ll_tmp1) THEN 
     849                      zcpx(ji,jj) = 1.0_wp 
     850                    ELSE IF(ll_tmp2) THEN 
     851                       ! no worries about zsshp2_e(ji+1,jj)-zsshp2_e(ji,jj) = 0, it won't happen ! here 
     852                      zcpx(ji,jj) = ABS((zsshp2_e(ji+1,jj) + bathy(ji+1,jj)*tmask(ji+1,jj,1) - zsshp2_e(ji,jj) - bathy(ji,jj)*tmask(ji,jj,1)) /& 
     853                                 & (zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj))) 
     854                    ELSE 
     855                       zcpx(ji,jj) = 0._wp 
     856                    END IF 
     857                 ENDIF 
     858 
     859                 IF ( tmask(ji,jj+1,1) == 0._wp ) THEN 
     860                    zcpy = 1._wp 
     861                 ELSE 
     862                    ll_tmp1 = MIN(zsshp2_e(ji,jj), zsshp2_e(ji,jj+1)) > MAX(-bathy(ji,jj)*tmask(ji,jj,1), -bathy(ji,jj+1)*tmask(ji,jj+1,1))& 
     863                           & .and. MAX(zsshp2_e(ji,jj) + bathy(ji,jj)*tmask(ji,jj,1), zsshp2_e(ji,jj+1) + bathy(ji,jj+1)*tmask(ji,jj+1,1)) & 
     864                           &  > rn_wdmin1 + rn_wdmin2 
     865                    ll_tmp2 = MAX(zsshp2_e(ji,jj), zsshp2_e(ji,jj+1)) > MAX(-bathy(ji,jj)*tmask(ji,jj,1), -bathy(ji,jj+1)*tmask(ji,jj+1,1)) + & 
     866                                                           & rn_wdmin1 + rn_wdmin2 
     867                    IF(ll_tmp1) THEN 
     868                       zcpy(ji,jj) = 1.0_wp 
     869                    ELSE IF(ll_tmp2) THEN 
     870                       ! no worries about zsshp2_e(ji,jj+1)-zsshp2_e(ji,jj) = 0, it won't happen ! here 
     871                      zcpy(ji,jj) = ABS((zsshp2_e(ji,jj+1) + bathy(ji,jj+1)*tmask(ji,jj+1,1) - zsshp2_e(ji,jj) - bathy(ji,jj)*tmask(ji,jj,1)) /& 
     872                                  & (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj))) 
     873                    ELSE 
     874                      zcpy(ji,jj) = 0._wp 
     875                    END IF 
     876                 ENDIF 
     877              END DO 
     878            END DO 
     879            CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp ) 
     880         ENDIF 
     881 
     882         IF(ln_wd) THEN 
     883           DO jj = 2, jpjm1 
     884              DO ji = 2, jpim1 
     885                 ! Add surface pressure gradient 
     886                 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj) 
     887                 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj) 
     888                 zwx(ji,jj) = zu_spg * zcpx(ji,jj) 
     889                 zwy(ji,jj) = zv_spg * zcpy(ji,jj) 
     890              END DO 
     891           END DO 
     892         ELSE 
     893           DO jj = 2, jpjm1 
     894              DO ji = fs_2, fs_jpim1   ! vector opt. 
     895                 ! Add surface pressure gradient 
     896                 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj) 
     897                 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj) 
     898                 zwx(ji,jj) = zu_spg 
     899                 zwy(ji,jj) = zv_spg 
     900              END DO 
     901           END DO 
     902         END IF 
    753903         ! 
    754904         ! Set next velocities: 
     
    774924               DO ji = fs_2, fs_jpim1   ! vector opt. 
    775925 
    776                   zhura = umask(ji,jj,1)/(hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - umask(ji,jj,1)) 
    777                   zhvra = vmask(ji,jj,1)/(hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - vmask(ji,jj,1)) 
    778  
    779                   ua_e(ji,jj) = (                hu_e(ji,jj)  *  zun_e(ji,jj)   &  
    780                             &     + rdtbt * ( zhust_e(ji,jj)  *    zwx(ji,jj)   &  
     926                  IF(ln_wd) THEN 
     927                    zhura = MAX(hu_0(ji,jj) + zsshu_a(ji,jj), rn_wdmin1) 
     928                    zhvra = MAX(hv_0(ji,jj) + zsshv_a(ji,jj), rn_wdmin1) 
     929                  ELSE 
     930                    zhura = hu_0(ji,jj) + zsshu_a(ji,jj) 
     931                    zhvra = hv_0(ji,jj) + zsshv_a(ji,jj) 
     932                  END IF 
     933 
     934                  zhura = umask(ji,jj,1)/(zhura + 1._wp - umask(ji,jj,1)) 
     935                  zhvra = vmask(ji,jj,1)/(zhvra + 1._wp - vmask(ji,jj,1)) 
     936 
     937                  ua_e(ji,jj) = (                hu_e(ji,jj)  *  zun_e(ji,jj)   & 
     938                            &     + rdtbt * ( zhust_e(ji,jj)  *    zwx(ji,jj)   & 
    781939                            &               + zhup2_e(ji,jj)  * zu_trd(ji,jj)   & 
    782940                            &               +      hu(ji,jj)  * zu_frc(ji,jj) ) & 
     
    794952         IF( lk_vvl ) THEN                             !* Update ocean depth (variable volume case only) 
    795953            !                                          !  ----------------------------------------------         
    796             hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 
    797             hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 
     954            IF(ln_wd) THEN 
     955              hu_e (:,:) = MAX(hu_0(:,:) + zsshu_a(:,:), rn_wdmin1 * umask(:,:,1) ) 
     956              hv_e (:,:) = MAX(hv_0(:,:) + zsshv_a(:,:), rn_wdmin1 * vmask(:,:,1) ) 
     957            ELSE 
     958              hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 
     959              hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 
     960            END IF 
     961 
    798962            hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1._wp - umask(:,:,1) ) 
    799963            hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1._wp - vmask(:,:,1) ) 
     
    9261090      CALL wrk_dealloc( jpi, jpj, zsshu_a, zsshv_a                                   ) 
    9271091      CALL wrk_dealloc( jpi, jpj, zhf ) 
     1092      IF(ln_wd) CALL wrk_dealloc( jpi, jpj, zcpx, zcpy ) 
    9281093      ! 
    9291094      IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_ts') 
  • branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90

    r5029 r5954  
    285285               zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1) 
    286286               zx2 = zwx(ji  ,jj) + zwx(ji  ,jj+1) 
    287                pua(ji,jj,jk) = pua(ji,jj,jk) + zfact2 / e1u(ji,jj) * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 
    288                pva(ji,jj,jk) = pva(ji,jj,jk) - zfact2 / e2v(ji,jj) * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )  
     287               pua(ji,jj,jk) = pua(ji,jj,jk) + zfact2 / e1u(ji,jj) * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) * umask(ji,jj,jk) 
     288               pva(ji,jj,jk) = pva(ji,jj,jk) - zfact2 / e2v(ji,jj) * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 ) * vmask(ji,jj,jk) 
    289289            END DO   
    290290         END DO   
     
    404404               zcva =-zfact2 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 ) 
    405405               ! mixed vorticity trend added to the momentum trends 
    406                ua(ji,jj,jk) = ua(ji,jj,jk) + zcua + zua 
    407                va(ji,jj,jk) = va(ji,jj,jk) + zcva + zva 
     406               ua(ji,jj,jk) = ua(ji,jj,jk) + ( zcua + zua ) * umask(ji,jj,jk) 
     407               va(ji,jj,jk) = va(ji,jj,jk) + ( zcva + zva ) * vmask(ji,jj,jk) 
    408408            END DO   
    409409         END DO   
     
    522522               zvau =-zfact1 / e2v(ji,jj) * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1)   & 
    523523                  &                         + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) 
    524                pua(ji,jj,jk) = pua(ji,jj,jk) + zuav * ( zwz(ji  ,jj-1) + zwz(ji,jj) ) 
    525                pva(ji,jj,jk) = pva(ji,jj,jk) + zvau * ( zwz(ji-1,jj  ) + zwz(ji,jj) ) 
     524               pua(ji,jj,jk) = pua(ji,jj,jk) + zuav * ( zwz(ji  ,jj-1) + zwz(ji,jj) ) * umask(ji,jj,jk) 
     525               pva(ji,jj,jk) = pva(ji,jj,jk) + zvau * ( zwz(ji-1,jj  ) + zwz(ji,jj) ) * vmask(ji,jj,jk) 
    526526            END DO   
    527527         END DO   
     
    691691               zva = - zfac12 / e2v(ji,jj) * (  ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji  ,jj+1)   & 
    692692                  &                           + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) ) 
    693                pua(ji,jj,jk) = pua(ji,jj,jk) + zua 
    694                pva(ji,jj,jk) = pva(ji,jj,jk) + zva 
     693               pua(ji,jj,jk) = pua(ji,jj,jk) + zua * umask(ji,jj,jk) 
     694               pva(ji,jj,jk) = pva(ji,jj,jk) + zva * vmask(ji,jj,jk) 
    695695            END DO   
    696696         END DO   
  • branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90

    r5120 r5954  
    124124               zva = - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) ) 
    125125               !                         ! add the trends to the general momentum trends 
    126                ua(ji,jj,jk) = ua(ji,jj,jk) + zua 
    127                va(ji,jj,jk) = va(ji,jj,jk) + zva 
     126               ua(ji,jj,jk) = ua(ji,jj,jk) + zua * umask(ji,jj,jk) 
     127               va(ji,jj,jk) = va(ji,jj,jk) + zva * vmask(ji,jj,jk) 
    128128            END DO   
    129129         END DO   
  • branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90

    r5942 r5954  
    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 
     
    3940   USE wrk_nemo        ! Memory Allocation 
    4041   USE timing          ! Timing 
     42   USE wadlmt          ! Wetting/Drying flux limting 
    4143 
    4244   IMPLICIT NONE 
     
    9193      ENDIF 
    9294      ! 
    93       CALL div_cur( kt )                              ! Horizontal divergence & Relative vorticity 
    94       ! 
    9595      z2dt = 2._wp * rdt                              ! set time step size (Euler/Leapfrog) 
    9696      IF( neuler == 0 .AND. kt == nit000 )   z2dt = rdt 
    97  
     97      ! 
     98      z1_rau0 = 0.5_wp * r1_rau0 
     99      ! 
     100      IF(ln_wd) CALL wad_lmt(sshb, z1_rau0 * (emp_b(:,:) + emp(:,:)), z2dt) 
     101 
     102      CALL div_cur( kt )                              ! Horizontal divergence & Relative vorticity 
     103      ! 
    98104      !                                           !------------------------------! 
    99105      !                                           !   After Sea Surface Height   ! 
     
    107113      ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations. 
    108114      !  
    109       z1_rau0 = 0.5_wp * r1_rau0 
    110115      ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * ssmask(:,:) 
    111116 
  • branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/ICB/icbrst.F90

    r5341 r5954  
    2424   USE icb_oce        ! define iceberg arrays 
    2525   USE icbutl         ! iceberg utility routines 
     26! added for ln_rstdate, should include separate branch at later date... 
     27   USE phycst         ! for rday 
     28   USE ioipsl, ONLY : ju2ymds    ! for calendar 
    2629 
    2730   IMPLICIT NONE 
     
    235238      TYPE(iceberg), POINTER :: this 
    236239      TYPE(point)  , POINTER :: pt 
    237       !!---------------------------------------------------------------------- 
     240! included for ln_rstdate.... 
     241      INTEGER             ::   iyear, imonth, iday 
     242      REAL (wp)           ::   zsec 
     243      CHARACTER(LEN=20)   ::   clkt     ! ocean time-step deine as a character 
     244      !!---------------------------------------------------------------------- 
     245 
     246      IF ( ln_rstdate ) THEN 
     247         CALL ju2ymds( fjulday + rdttra(1) / rday, iyear, imonth, iday, zsec )            
     248         WRITE(clkt, '(i4.4,2i2.2)') iyear, imonth, iday 
     249      ELSE 
     250         IF( kt > 999999999 ) THEN   ;   WRITE(clkt, *       ) kt 
     251         ELSE                        ;   WRITE(clkt, '(i8.8)') kt 
     252         ENDIF 
     253      ENDIF 
    238254 
    239255      ! Assume we write iceberg restarts to same directory as ocean restarts. 
    240256      cl_path = TRIM(cn_ocerst_outdir) 
    241257      IF( cl_path(LEN_TRIM(cl_path):) /= '/' ) cl_path = TRIM(cl_path) // '/' 
     258 
     259! changed for ln_rstdate.... 
    242260      IF( lk_mpp ) THEN 
    243          WRITE(cl_filename,'(A,"_icebergs_",I8.8,"_restart_",I4.4,".nc")') TRIM(cexper), kt, narea-1 
     261         WRITE(cl_filename,'(A,"_icebergs_",A,"_restart_",I4.4,".nc")') TRIM(cexper), TRIM(ADJUSTL(clkt)), narea-1 
    244262      ELSE 
    245          WRITE(cl_filename,'(A,"_icebergs_",I8.8,"_restart.nc")') TRIM(cexper), kt 
     263         WRITE(cl_filename,'(A,"_icebergs_",A,"_restart.nc")') TRIM(cexper), TRIM(ADJUSTL(clkt)) 
    246264      ENDIF 
    247265      IF (nn_verbose_level >= 0) WRITE(numout,'(2a)') 'icebergs, write_restart: creating ',TRIM(cl_path)//TRIM(cl_filename) 
  • branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/IOM/in_out_manager.F90

    r5518 r5954  
    4848   LOGICAL       ::   ln_clobber       !: clobber (overwrite) an existing file 
    4949   INTEGER       ::   nn_chunksz       !: chunksize (bytes) for NetCDF file (works only with iom_nf90 routines) 
     50! included for rstdate and rstdir, should be moved to seperate branch! 
     51   LOGICAL       ::   ln_rstdate       !: datestamping of restarts 
     52   CHARACTER(lc) ::   cn_ocerst_indir  !: restart input directory 
     53   CHARACTER(lc) ::   cn_ocerst_outdir !: restart output directory 
    5054#if defined key_netcdf4 
    5155   !!---------------------------------------------------------------------- 
  • branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90

    r5407 r5954  
    2424   USE trdmxl_oce      ! ocean active mixed layer tracers trends variables 
    2525   USE divcur          ! hor. divergence and curl      (div & cur routines) 
     26! for ln_rstdate.... 
     27   USE ioipsl, ONLY : ju2ymds    ! for calendar 
    2628 
    2729   IMPLICIT NONE 
     
    5860      CHARACTER(LEN=50)   ::   clname   ! ocean output restart file name 
    5961      CHARACTER(lc)       ::   clpath   ! full path to ocean output restart file 
     62! Included for ln_rstdate 
     63      INTEGER             ::   iyear, imonth, iday 
     64      REAL (wp)           ::   zsec 
     65 
    6066      !!---------------------------------------------------------------------- 
    6167      ! 
     
    8187      IF( kt == nitrst - 1 .OR. nstock == 1 .OR. ( kt == nitend .AND. .NOT. lrst_oce ) ) THEN 
    8288         IF( nitrst <= nitend .AND. nitrst > 0 ) THEN  
    83             ! beware of the format used to write kt (default is i8.8, that should be large enough...) 
    84             IF( nitrst > 999999999 ) THEN   ;   WRITE(clkt, *       ) nitrst 
    85             ELSE                            ;   WRITE(clkt, '(i8.8)') nitrst 
     89            IF ( ln_rstdate ) THEN 
     90               CALL ju2ymds( fjulday + rdttra(1) / rday, iyear, imonth, iday, zsec )            
     91               WRITE(clkt, '(i4.4,2i2.2)') iyear, imonth, iday 
     92            ELSE 
     93               ! beware of the format used to write kt (default is i8.8, that should be large enough...) 
     94               IF( nitrst > 999999999 ) THEN   ;   WRITE(clkt, *       ) nitrst 
     95               ELSE                            ;   WRITE(clkt, '(i8.8)') nitrst 
     96               ENDIF 
    8697            ENDIF 
    8798            ! create the file 
     
    89100            clpath = TRIM(cn_ocerst_outdir) 
    90101            IF( clpath(LEN_TRIM(clpath):) /= '/' ) clpath = TRIM(clpath) // '/' 
     102            ! 
    91103            IF(lwp) THEN 
    92104               WRITE(numout,*) 
     
    102114               ENDIF 
    103115            ENDIF 
    104             ! 
     116               ! 
    105117            CALL iom_open( TRIM(clpath)//TRIM(clname), numrow, ldwrt = .TRUE., kiolib = jprstlib ) 
    106118            lrst_oce = .TRUE. 
     
    172184      LOGICAL        ::   llok 
    173185      CHARACTER(lc)  ::   clpath   ! full path to ocean output restart file 
     186      CHARACTER(lc)  ::   clpath   ! full path to ocean output restart file 
    174187      !!---------------------------------------------------------------------- 
    175188      ! 
  • branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90

    r5942 r5954  
    9191   REAL(wp) ::   rn_zqt      ! z(q,t) : height of humidity and temperature measurements 
    9292   REAL(wp) ::   rn_zu       ! z(u)   : height of wind measurements 
     93   LOGICAL  ::   ln_charnock ! logical flag for charnock wind stress in surge model(true) or not(false) 
    9394 
    9495   !! * Substitutions 
     
    151152         &                  sn_wndi, sn_wndj, sn_humi  , sn_qsr ,           & 
    152153         &                  sn_qlw , sn_tair, sn_prec  , sn_snow,           & 
    153          &                  sn_tdif, rn_zqt,  rn_zu 
     154         &                  sn_tdif, rn_zqt , rn_zu, ln_charnock  
    154155      !!--------------------------------------------------------------------- 
    155156      ! 
     
    247248      INTEGER  ::   ji, jj               ! dummy loop indices 
    248249      REAL(wp) ::   zcoef_qsatw, zztmp   ! local variable 
     250      REAL(wp) ::   z_z0, z_Cd1          ! local variable 
     251      REAL(wp) ::   i                    ! local variable 
     252      REAL(wp) ::   charn_const=0.0275    ! local variable 
    249253      REAL(wp), DIMENSION(:,:), POINTER ::   zwnd_i, zwnd_j    ! wind speed components at T-point 
    250254      REAL(wp), DIMENSION(:,:), POINTER ::   zqsatw            ! specific humidity at pst 
     
    300304      !      I   Radiative FLUXES                                                     ! 
    301305      ! ----------------------------------------------------------------------------- ! 
    302  
     306       
    303307      ! ocean albedo assumed to be constant + modify now Qsr to include the diurnal cycle                    ! Short Wave 
    304308      zztmp = 1. - albo 
     309#if defined key_surge 
     310      qsr(:,:)=0._wp 
     311      zqlw(:,:) = 0._wp 
     312#else 
    305313      IF( ln_dm2dc ) THEN   ;   qsr(:,:) = zztmp * sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) * tmask(:,:,1) 
    306314      ELSE                  ;   qsr(:,:) = zztmp *          sf(jp_qsr)%fnow(:,:,1)   * tmask(:,:,1) 
     
    308316 
    309317      zqlw(:,:) = (  sf(jp_qlw)%fnow(:,:,1) - Stef * zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:)  ) * tmask(:,:,1)   ! Long  Wave 
     318#endif 
    310319      ! ----------------------------------------------------------------------------- ! 
    311320      !     II    Turbulent FLUXES                                                    ! 
    312321      ! ----------------------------------------------------------------------------- ! 
    313  
    314       ! ... specific humidity at SST and IST 
    315       zqsatw(:,:) = zcoef_qsatw * EXP( -5107.4 / zst(:,:) ) 
    316  
    317       ! ... NCAR Bulk formulae, computation of Cd, Ch, Ce at T-point : 
    318       CALL turb_core_2z( rn_zqt, rn_zu, zst, sf(jp_tair)%fnow, zqsatw, sf(jp_humi)%fnow, wndm,   & 
    319          &               Cd, Ch, Ce, zt_zu, zq_zu ) 
    320      
     322      IF (ln_charnock) THEN 
     323          Cd(:,:)=0.0001_wp 
     324          DO jj = 1,jpj 
     325             DO ji = 1,jpi 
     326                z_Cd1=0._wp 
     327                i=1 
     328                !Iterate 
     329                DO WHILE((abs(Cd(ji,jj)-z_Cd1))>1E-6) 
     330                z_Cd1=Cd(ji,jj) 
     331                z_z0=charn_const*z_Cd1*wndm(ji,jj)**2/grav 
     332                Cd(ji,jj)=(0.41_wp/log(10._wp/z_z0))**2 
     333                i=i+1 
     334                ENDDO 
     335             ENDDO 
     336          ENDDO 
     337      ELSE 
     338 
     339        ! ... specific humidity at SST and IST 
     340        zqsatw(:,:) = zcoef_qsatw * EXP( -5107.4 / zst(:,:) )  
     341 
     342        ! ... NCAR Bulk formulae, computation of Cd, Ch, Ce at T-point : 
     343        CALL turb_core_2z( rn_zqt, rn_zu, zst, sf(jp_tair)%fnow, zqsatw, sf(jp_humi)%fnow, wndm,   & 
     344           &               Cd, Ch, Ce, zt_zu, zq_zu ) 
     345 
     346      ENDIF 
     347 
    321348      ! ... tau module, i and j component 
    322349      DO jj = 1, jpj 
     
    352379      !  Turbulent fluxes over ocean 
    353380      ! ----------------------------- 
     381#if ! defined key_surge 
    354382      IF( ABS( rn_zu - rn_zqt) < 0.01_wp ) THEN 
    355383         !! q_air and t_air are (or "are almost") given at 10m (wind reference height) 
     
    363391      ENDIF 
    364392      zqla (:,:) = Lv * zevap(:,:)                                                              ! Latent Heat 
     393#endif 
    365394 
    366395      IF(ln_ctl) THEN 
     
    379408      ! ----------------------------------------------------------------------------- ! 
    380409      ! 
     410#if defined key_surge 
     411      emp (:,:) = 0._wp 
     412      qns(:,:)  = 0._wp 
     413#else 
    381414      emp (:,:) = (  zevap(:,:)                                          &   ! mass flux (evap. - precip.) 
    382415         &         - sf(jp_prec)%fnow(:,:,1) * rn_pfac  ) * tmask(:,:,1) 
     
    389422         &     + sf(jp_snow)%fnow(:,:,1) * rn_pfac                                &   ! add solid  precip heat content at min(Tair,Tsnow) 
    390423         &     * ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) 
     424#endif 
    391425      ! 
    392426#if defined key_lim3 
  • branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90

    r5942 r5954  
    250250         nstop = nstop + 1 
    251251      ENDIF 
     252      IF ( lk_surge .and. .not. ( ln_blk_core .or. ln_ana ) ) & 
     253         &   CALL ctl_stop( ' surge model only compatible with analytical fluxes or core formulae' )      
    252254      IF(lwp) THEN 
    253255         WRITE(numout,*) 
  • branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/SBC/tide.h90

    r4292 r5954  
    2828   Wave(18) = tide(  'L2'     , 0.006694 ,    2   ,  2 , -1 ,  2 , -1 ,  0  , +180  ,  2   , -2   ,  0   ,  0   , 0 ,  215    ) 
    2929   Wave(19) = tide(  'T2'     , 0.006614 ,    2   ,  2 ,  0 , -1 ,  0 ,  1  ,    0  ,  0   ,  0   ,  0   ,  0   , 0 ,    0    ) 
     30   ! 
     31   !             !! name_tide , equitide , nutide , nt , ns , nh , np , np1 , shift , nksi , nnu0 , nnu1 , nnu2 , R , formula !! 
     32   Wave(20) = tide(  'MNS2'   , 0.000000 ,    2   ,  2 , -5 ,  4 ,  1 ,  0  ,    0  ,  4   , -4   ,  0   ,  0   , 0 ,    6    ) 
     33   Wave(21) = tide(  'Lam2'   , 0.001760 ,    2   ,  2 , -1 ,  0 ,  1 ,  0  , +180  ,  2   , -2   ,  0   ,  0   , 0 ,   78    ) 
     34   Wave(22) = tide(  'MSN2'   , 0.000000 ,    2   ,  2 ,  1 ,  0 ,  1 ,  0  ,    0  ,  2   , -2   ,  0   ,  2   , 0 ,    6    ) 
     35   Wave(23) = tide(  '2SM2'   , 0.000000 ,    2   ,  2 ,  2 , -2 ,  0 ,  0  ,    0  , -2   ,  2   ,  0   ,  0   , 0 ,   16    ) 
     36   Wave(24) = tide(  'MO3'    , 0.000000 ,    3   ,  3 , -4 ,  1 ,  0 ,  0  ,  +90  ,  2   , -2   ,  0   ,  0   , 0 ,   13    ) 
     37   Wave(25) = tide(  'MK3'    , 0.000000 ,    3   ,  3 , -2 ,  3 ,  0 ,  0  ,  -90  ,  2   , -2   , -1   ,  0   , 0 ,   10    ) 
     38   Wave(26) = tide(  'MN4'    , 0.000000 ,    4   ,  4 , -5 ,  4 ,  1 ,  0  ,    0  ,  4   , -4   ,  0   ,  0   , 0 ,    1    ) 
     39   Wave(27) = tide(  'MS4'    , 0.000000 ,    4   ,  4 , -2 ,  2 ,  0 ,  0  ,    0  ,  2   , -2   ,  0   ,  0   , 0 ,    2    ) 
     40   Wave(28) = tide(  'M6'     , 0.000000 ,    6   ,  6 , -6 ,  6 ,  0 ,  0  ,    0  ,  6   , -6   ,  0   ,  0   , 0 ,    4    ) 
     41   Wave(29) = tide(  '2MS6'   , 0.000000 ,    6   ,  6 , -4 ,  4 ,  0 ,  0  ,    0  ,  4   , -4   ,  0   ,  0   , 0 ,    6    ) 
     42   Wave(30) = tide(  '2MK6'   , 0.000000 ,    6   ,  6 , -4 ,  6 ,  0 ,  0  ,    0  ,  4   , -4   ,  0   , -2   , 0 ,    5    ) 
     43   Wave(31) = tide(  '3M2S2'  , 0.000000 ,    2   , 2  , -6 ,  6 ,  0 ,  0  ,    0  ,  6   , -6   ,  0   ,  0   , 0 ,   12    ) 
  • branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/SBC/tide_mod.F90

    r5215 r5954  
    1616   PUBLIC   tide_init_Wave   ! called by tideini and diaharm modules 
    1717 
    18    INTEGER, PUBLIC, PARAMETER ::   jpmax_harmo = 19   !: maximum number of harmonic 
     18   INTEGER, PUBLIC, PARAMETER ::   jpmax_harmo = 31   !: maximum number of harmonic 
    1919 
    2020   TYPE, PUBLIC ::    tide 
    21       CHARACTER(LEN=4) ::   cname_tide 
     21      CHARACTER(LEN=5) ::   cname_tide 
    2222      REAL(wp)         ::   equitide 
    2323      INTEGER          ::   nutide 
  • branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/SBC/tideini.F90

    r5215 r5954  
    4848    INTEGER  :: ji, jk 
    4949    INTEGER, INTENT( in ) ::   kt     ! ocean time-step 
    50     CHARACTER(LEN=4), DIMENSION(jpmax_harmo) :: clname 
     50    CHARACTER(LEN=5), DIMENSION(jpmax_harmo) :: clname 
    5151    INTEGER  ::   ios                 ! Local integer output status for namelist read 
    5252    ! 
  • branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfini.F90

    r5386 r5954  
    115115         ioptio = ioptio+1 
    116116      ENDIF 
    117       IF( ioptio == 0 .OR. ioptio > 1 .AND. .NOT. lk_esopa )   & 
    118          &   CALL ctl_stop( ' one and only one vertical diffusion option has to be defined ' ) 
     117      IF( ioptio > 1 .AND. .NOT. lk_esopa )   & 
     118         &   CALL ctl_stop( 'only one vertical diffusion option has to be defined ' ) 
     119      IF( ioptio == 0 .AND. .NOT. ( lk_esopa .OR. lk_surge ) )   & 
     120         &   CALL ctl_stop( 'a vertical diffusion option has to be defined ' ) 
    119121      IF( ( lk_zdfric .OR. lk_zdfgls .OR. lk_zdfkpp ) .AND. ln_isfcav )   & 
    120122         &   CALL ctl_stop( ' only zdfcst and zdftke were tested with ice shelves cavities ' ) 
     
    151153      ENDIF 
    152154      IF ( ioptio > 1 .AND. .NOT. lk_esopa )   CALL ctl_stop( ' chose between ln_zdfnpc and ln_zdfevd' ) 
    153       IF( ioptio == 0 .AND. .NOT.( lk_zdftke .OR. lk_zdfgls .OR. lk_zdfkpp ) )           & 
    154          CALL ctl_stop( ' except for TKE, GLS or KPP physics, a convection scheme is',   & 
     155      IF( ioptio == 0 .AND. .NOT.( lk_zdftke .OR. lk_zdfgls .OR. lk_zdfkpp .OR. lk_surge ) )           & 
     156         CALL ctl_stop( ' except for TKE, GLS or KPP physics (or running in surge mode), a convection scheme is',   & 
    155157         &              ' required: ln_zdfevd or ln_zdfnpc logicals' ) 
    156158 
  • branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90

    r5942 r5954  
    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 
     
    228229      NAMELIST/namcfg/ cp_cfg, cp_cfz, jp_cfg, jpidta, jpjdta, jpkdta, jpiglo, jpjglo, & 
    229230         &             jpizoom, jpjzoom, jperio, ln_use_jattr 
     231      NAMELIST/namwad/ ln_wd, rn_wdmin1, rn_wdmin2, rn_wdld, nn_wdit 
    230232      !!---------------------------------------------------------------------- 
    231233      ! 
     
    253255      READ  ( numnam_cfg, namcfg, IOSTAT = ios, ERR = 904 ) 
    254256904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcfg in configuration namelist', .TRUE. )    
     257 
     258      REWIND( numnam_ref )              ! Namelist namwad in reference namelist : Parameters for Wetting/Drying 
     259      READ  ( numnam_ref, namwad, IOSTAT = ios, ERR = 905) 
     260905   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namwad in reference namelist', lwp ) 
     261 
     262      REWIND( numnam_cfg )              ! Namelist namwad in configuration namelist : Parameters for Wetting/Drying 
     263      READ  ( numnam_cfg, namwad, IOSTAT = ios, ERR = 906) 
     264906   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namwad in configuration namelist', lwp ) 
    255265 
    256266! Force values for AGRIF zoom (cf. agrif_user.F90) 
     
    311321         WRITE( numond, namctl ) 
    312322         WRITE( numond, namcfg ) 
     323         WRITE( numond, namwad ) 
    313324      ENDIF 
    314325 
     
    416427         &                  CALL zdf_ddm_init      ! double diffusive mixing 
    417428      !                                         ! Lateral physics 
     429#if ! defined key_surge 
    418430                            CALL ldf_tra_init      ! Lateral ocean tracer physics 
     431#endif 
    419432                            CALL ldf_dyn_init      ! Lateral ocean momentum physics 
    420433      IF( lk_ldfslp     )   CALL ldf_slp_init      ! slope of lateral mixing 
    421434 
     435#if ! defined key_surge 
    422436      !                                     ! Active tracers 
    423437                            CALL tra_qsr_init   ! penetrative solar radiation qsr 
     
    428442                            CALL tra_ldf_init   ! lateral mixing 
    429443                            CALL tra_zdf_init   ! vertical mixing and after tracer fields 
     444#endif 
    430445 
    431446      !                                     ! Dynamics 
     
    435450                            CALL dyn_ldf_init   ! lateral mixing 
    436451                            CALL dyn_hpg_init   ! horizontal gradient of Hydrostatic pressure 
     452#if ! defined key_surge 
    437453                            CALL dyn_zdf_init   ! vertical diffusion 
     454#endif 
    438455                            CALL dyn_spg_init   ! surface pressure gradient 
    439456 
     
    520537         WRITE(numout,*) '      use file attribute if exists as i/p j-start ln_use_jattr = ', ln_use_jattr 
    521538      ENDIF 
     539 
     540      IF(lwp) THEN                  ! control print 
     541         WRITE(numout,*) 
     542         WRITE(numout,*) 'namwad  : wetting and drying initialisation through namelist read' 
     543         WRITE(numout,*) '~~~~~~~ ' 
     544         WRITE(numout,*) '   Namelist namwad' 
     545         WRITE(numout,*) '   key to turn on/off wetting/drying              ln_wd      = ', ln_wd 
     546         WRITE(numout,*) '   minimum water depth on dried cells             rn_wdmin1  = ', rn_wdmin1 
     547         WRITE(numout,*) '   tolerrance of min water depth on dried cells   rn_wdmin2  = ', rn_wdmin2 
     548         WRITE(numout,*) '   land elevation below which wad considered      rn_wdld    = ', rn_wdld 
     549         WRITE(numout,*) '   maximum number of iteration for W/D limiter    nn_wdit    = ', nn_wdit 
     550      ENDIF  
     551 
    522552      !                             ! Parameter control 
    523553      ! 
  • branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/par_oce.F90

    r5118 r5954  
    101101#endif 
    102102 
     103#if defined key_vectopt_loop 
     104   LOGICAL, PUBLIC, PARAMETER ::   lk_vopt_loop = .TRUE.   !: vector optimization flag 
     105#else 
     106   LOGICAL, PUBLIC, PARAMETER ::   lk_vopt_loop = .FALSE.  !: vector optimization flag 
     107#endif 
     108#if defined key_surge 
     109   LOGICAL, PUBLIC, PARAMETER ::   lk_surge     = .TRUE.   ! flag for 2d baratropic modelling 
     110#else 
     111   LOGICAL, PUBLIC, PARAMETER ::   lk_surge     = .FALSE.  ! flag for 2d baratropic modelling 
     112#endif 
     113 
    103114   !!---------------------------------------------------------------------- 
    104115   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
  • branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/step.F90

    r5510 r5954  
    246246      ! Active tracers                              (ua, va used as workspace) 
    247247      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     248#if ! defined key_surge 
    248249                             tsa(:,:,:,:) = 0.e0            ! set tracer trends to zero 
    249250 
     
    294295                             CALL tra_nxt( kstp )                ! tracer fields at next time step 
    295296      ENDIF 
    296  
     297#endif 
    297298      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    298299      ! Dynamics                                    (tsa used as workspace) 
     
    309310 
    310311                               CALL dyn_bfr( kstp )         ! bottom friction 
     312#if ! defined key_surge 
    311313                               CALL dyn_zdf( kstp )         ! vertical diffusion 
     314#endif 
    312315      ELSE 
    313316                               ua(:,:,:) = 0.e0             ! set dynamics trends to zero 
     
    328331                               CALL dyn_hpg( kstp )         ! horizontal gradient of Hydrostatic pressure 
    329332                               CALL dyn_bfr( kstp )         ! bottom friction 
     333#if ! defined key_surge 
    330334                               CALL dyn_zdf( kstp )         ! vertical diffusion 
     335#endif 
    331336                               CALL dyn_spg( kstp, indic )  ! surface pressure gradient 
    332337      ENDIF 
  • branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/NEMO/OPA_SRC/step_oce.F90

    r5501 r5954  
    116116   USE trcstp           ! passive tracer time-stepping      (trc_stp routine) 
    117117#endif 
     118   USE wadlmt           ! wetting/drying limiter           (wad_lmt routine) 
    118119   !!---------------------------------------------------------------------- 
    119120   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
  • branches/UKMO/dev_r5518_Surge_Modelling/NEMOGCM/fcm-make/inc/common.cfg

    r4363 r5954  
    66extract.path-root[nemo] = NEMOGCM 
    77extract.path-excl[nemo] = / \ 
    8                         \ NEMO/OPA_SRC/TRD/trdmod_trc.F90 \ 
     8                        \ NEMO/OPA_SRC/TRD/trdtrc.F90 \ 
    99                        \ NEMO/LIM_SRC_2/limrhg.F90 \ 
    1010                        \ TOOLS/OBSTOOLS/src/obs_prof_io.F90 \ 
Note: See TracChangeset for help on using the changeset viewer.