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.
wadlmt.F90 in branches/2013/dev_r4050_NOC_WaD/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2013/dev_r4050_NOC_WaD/NEMOGCM/NEMO/OPA_SRC/DYN/wadlmt.F90 @ 4355

Last change on this file since 4355 was 4355, checked in by hliu, 10 years ago

update Wetting/drying branch. Will keep commit in the next a few days.
File size: 6.9 KB
Line 
1
2MODULE wadlmt
3   !!==============================================================================
4   !!                       ***  MODULE  wadlmt  ***
5   !! compute the flux limiter for water depth positivity in !wetting/drying case
6   !!==============================================================================
7   !! History :   
8   !!  NEMO      3.5  ! 2014-01  ((H.Liu)  Original code
9   !!----------------------------------------------------------------------
10
11   !!----------------------------------------------------------------------
12   !!   wad_lmt    : Compute the horizontal flux limiter and the limited velocity
13   !!                when wetting and drying happens
14   !!----------------------------------------------------------------------
15   USE oce             ! ocean dynamics and tracers
16   USE dom_oce         ! ocean space and time domain
17   USE sbc_oce, ONLY : ln_rnf   ! surface boundary condition: ocean
18   USE sbcrnf          ! river runoff
19   USE cla             ! cross land advection             (cla_div routine)
20   USE in_out_manager  ! I/O manager
21   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
22   USE lib_mpp         ! MPP library
23   USE wrk_nemo        ! Memory Allocation
24   USE timing          ! Timing
25
26   IMPLICIT NONE
27   PRIVATE
28
29
30   PUBLIC   wad_lmt    ! routine called by step.F90
31
32   SUBROUTINE wad_lmt( kt )
33      !!----------------------------------------------------------------------
34      !!                  ***  ROUTINE wad_lmt  ***
35      !!                   
36      !! ** Purpose :   generate flux limiters for wetting/drying
37      !!
38      !! ** Method  : - Prevent negative depth occurring (Not ready for Agrif)
39      !!
40      !! ** Action  : - update: uwdlmt(:,:), vwdlmt(:,:)
41      !!----------------------------------------------------------------------
42      INTEGER, INTENT(in) ::   kt       ! ocean time-step index
43      !
44      INTEGER  ::   ji, jj, jk, jk1     ! dummy loop indices
45      INTEGER  ::   zflag, z2dt         ! local scalar
46      REAL(wp) ::   zflxp, zflxn        ! local scalars
47      REAL(wp) ::   zcoef, zdep1, zdep2 ! local scalars
48      REAL(wp) ::   ztmp                ! local scalars
49      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxu,  zflxv   ! specific 2D workspace
50      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxu1, zflxv1  ! specific 2D workspace
51      REAL(wp), POINTER,  DIMENSION(:,:) ::   uwdlmt, vwdlmt  !: W/D flux limiters
52      !!----------------------------------------------------------------------
53      !
54
55      IF( nn_timing == 1 )  CALL timing_start('wad_lmt')
56
57      IF(ln_wad) THEN
58
59        CALL wrk_alloc( jpi, jpj, zflxu, zflxv, zflxu1, zflxv1, uwdlmt, vwdlmt)
60        !
61       
62        IF(lwp) WRITE(numout,*)
63        IF(lwp) WRITE(numout,*) 'wad_lmt : wetting/drying limiters and velocity limiting'
64       
65        z2dt = 2. * rdt                                 ! Euler or leap-frog time step
66        IF( neuler == 0 .AND. kt == nit000 )  z2dt = rdt
67       
68       
69        zflxu(:,:)  = 0._wp
70        zflxv(:,:)  = 0._wp
71        uwdlmt(:,:) = 1._wp
72        vwdlmt(:,:) = 1._wp
73       
74       
75        ! Horizontal Flux in u and v direction
76        DO jk = 1, jpkm1 
77           DO jj = 1, jpjm1
78              DO ji = 1, jpim1
79                 zflxu(ji,jj) = zflxu(ji,jj) + fse3u(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)
80                 zflxv(ji,jj) = zflxv(ji,jj) + fse3v(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
81              END DO 
82           END DO 
83        END DO
84       
85        zflxu(:,:) = zflxu(:,:) * e2u(:,:)
86        zflxv(:,:) = zflxv(:,:) * e1v(:,:)
87       
88        DO jk1 = 1, 10
89       
90        zflag = 0     ! flag marking that if we need further iteration of W/D limiters?
91       
92        zflxu1(:,:) = zflxu(:,:) * uwdlmt(:,:)
93        zflxv1(:,:) = zflxv(:,:) * vwdlmt(:,:)
94       
95       
96        DO jj = 2, jpjm1
97           DO ji = fs_2, fs_jpim1   ! vector opt.
98              zflxp = max(zflxu1(ji,jj), 0._wp) - min(zflxu1(ji-1,jj), 0._wp) + &
99                    & max(zflxv1(ji,jj), 0._wp) - min(zflxv1(ji,jj-1), 0._wp) 
100              zflxn = min(zflxu1(ji,jj), 0._wp) - max(zflxu1(ji-1,jj), 0._wp) + &
101                    & min(zflxv1(ji,jj), 0._wp) - max(zflxv1(ji,jj-1), 0._wp) 
102       
103              IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE 
104       
105              !IF( .NOT. AGRIF_Root() ) THEN    !Wetting and Drying but not for AGRIF
106                 IF (((nbondi ==  1).OR.(nbondi == 2)).AND.(ji == nlci-1)) CYCLE ! east
107                 IF (((nbondi == -1).OR.(nbondi == 2)).AND.(ji == 2     )) CYCLE ! west
108                 IF (((nbondj ==  1).OR.(nbondj == 2)).AND.(jj == nlcj-1)) CYCLE ! north
109                 IF (((nbondj == -1).OR.(nbondj == 2)).AND.(jj == 2     )) CYCLE ! south
110              !END IF
111       
112             
113              ztmp = e1t(ji,jj) * e2t(ji,jj)
114       
115              zdep1 = (zflxp + zflxn) * z2dt / ztmp
116              zdep2 = bathy(ji,jj) + sshb(ji,jj) - rn_wadmin 
117       
118              IF(zdep2 < 0._wp) THEN
119                !WRITE(numout,*) 'depth less than minimum depth, cell(ji,jj):', ji,jj
120                zdep2 = 0._wp
121                sshb(ji,jj) = rn_wad_dep1 - bathy(ji,jj)
122              END IF
123       
124              IF(zdep1 > zdep2) THEN
125                zflag = 1
126                zcoef = ( ( zdep2 + rn_wadmine ) * ztmp + zflxn * z2dt ) / ( zflxp * z2dt )
127                IF(zflxu1(ji,  jj) >= 0._wp) uwdlmt(ji,jj  ) = zcoef
128                IF(zflxu1(ji-1,jj) <  0._wp) uwdlmt(ji-1,jj) = zcoef
129                IF(zflxv1(ji,  jj) >= 0._wp) vwdlmt(ji,jj  ) = zcoef
130                IF(zflxv1(ji,jj-1) <  0._wp) vwdlmt(ji,jj-1) = zcoef
131              END IF
132           END DO ! ji loop
133        END DO  ! jj loop
134       
135        CALL lbc_lnk( uwdlmt, 'U', 1. )   ;   CALL lbc_lnk( vwdlmt , 'V', 1. ) 
136       
137        IF(zflag == 0) EXIT
138       
139        END DO  ! jk1 loop
140       
141        DO jk = 1, jpkm1 
142           DO jj = 2, jpjm1
143              DO ji = fs_2, fs_jpim1   ! vector opt.
144               un(ji,jj,jk)   = 0.5_wp * ( sign(1._wp, un(ji,jj,jk)   + 1._wp ) * uwdlmt(ji,jj) 
145               vn(ji,jj,jk)   = 0.5_wp * ( sign(1._wp, vn(ji,jj,jk)   + 1._wp ) * vwdlmt(ji,jj) 
146               un(ji-1,jj,jk) = 0.5_wp * ( 1._wp - sign(1._wp, un(ji-1,jj,jk) ) * uwdlmt(ji-1,jj) 
147               vn(ji,jj-1,jk) = 0.5_wp * ( 1._wp - sign(1._wp, vn(ji,jj-1,jk) ) * vwdlmt(ji,jj-1) 
148              END DO
149           END DO
150        END DO
151       
152        IF(zflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt!!!'
153       
154        !IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )          ! runoffs (update hdivn field)
155        !IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (update hdivn field)
156        !
157        !
158        CALL wrk_dealloc( jpi, jpj, zflxu, zflxv, zflxu1, zflxv1, uwdlmt, vwdlmt)
159      !
160      END IF
161
162      IF( nn_timing == 1 )  CALL timing_stop('wad_lmt')
163   END SUBROUTINE wad_lmt
164   !!======================================================================
165END MODULE wadlmt
Note: See TracBrowser for help on using the repository browser.