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 @ 4375

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

updated gravity filters
File size: 7.4 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   !! * Substitutions
33#  include "domzgr_substitute.h90"
34#  include "vectopt_loop_substitute.h90"
35CONTAINS
36
37   SUBROUTINE wad_lmt( kt )
38      !!----------------------------------------------------------------------
39      !!                  ***  ROUTINE wad_lmt  ***
40      !!                   
41      !! ** Purpose :   generate flux limiters for wetting/drying
42      !!
43      !! ** Method  : - Prevent negative depth occurring (Not ready for Agrif)
44      !!
45      !! ** Action  : - update: uwdlmt(:,:), vwdlmt(:,:)
46      !!----------------------------------------------------------------------
47      INTEGER, INTENT(in) ::   kt       ! ocean time-step index
48      !
49      INTEGER  ::   ji, jj, jk, jk1     ! dummy loop indices
50      INTEGER  ::   zflag, z2dt         ! local scalar
51      REAL(wp) ::   zcoef, zdep1, zdep2 ! local scalars
52      REAL(wp) ::   ztmp                ! local scalars
53      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxp,  zflxn   ! specific 2D workspace
54      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxu,  zflxv   ! specific 2D workspace
55      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxu1, zflxv1  ! specific 2D workspace
56      REAL(wp), POINTER,  DIMENSION(:,:) ::   wdlmt           !: W/D flux limiters
57
58      !!----------------------------------------------------------------------
59      !
60
61      IF( nn_timing == 1 )  CALL timing_start('wad_lmt')
62
63      IF(ln_wad) THEN
64
65        CALL wrk_alloc( jpi, jpj, zflxp, zflxn, zflxu, zflxv, zflxu1, zflxv1, wdlmt )
66        !
67       
68        IF(lwp) WRITE(numout,*)
69        IF(lwp) WRITE(numout,*) 'wad_lmt : wetting/drying limiters and velocity limiting'
70       
71        z2dt = 2. * rdt                                 ! Euler or leap-frog time step
72        IF( neuler == 0 .AND. kt == nit000 )  z2dt = rdt
73       
74        zflxp(:,:)  = 0._wp
75        zflxn(:,:)  = 0._wp
76        zflxu(:,:)  = 0._wp
77        zflxv(:,:)  = 0._wp
78        wdlmt(:,:)  = 1._wp
79       
80        ! Horizontal Flux in u and v direction
81        DO jk = 1, jpkm1 
82           DO jj = 1, jpjm1
83              DO ji = 1, jpim1
84                 zflxu(ji,jj) = zflxu(ji,jj) + fse3u(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)
85                 zflxv(ji,jj) = zflxv(ji,jj) + fse3v(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
86              END DO 
87           END DO 
88        END DO
89       
90        zflxu(:,:) = zflxu(:,:) * e2u(:,:)
91        zflxv(:,:) = zflxv(:,:) * e1v(:,:)
92       
93        DO jj = 2, jpjm1
94           DO ji = fs_2, fs_jpim1   ! vector opt.
95
96              IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE
97
98              zflxp(ji,jj) = max(zflxu(ji,jj), 0._wp) - min(zflxu(ji-1,jj),   0._wp) + &
99                           & max(zflxv(ji,jj), 0._wp) - min(zflxv(ji,  jj-1), 0._wp) 
100              zflxn(ji,jj) = min(zflxu(ji,jj), 0._wp) - max(zflxu(ji-1,jj),   0._wp) + &
101                           & min(zflxv(ji,jj), 0._wp) - max(zflxv(ji,  jj-1), 0._wp) 
102
103              zdep2 = bathy(ji,jj) + sshb(ji,jj) - rn_wadmin 
104              IF(zdep2 < 0._wp) THEN
105                zdep2 = 0._wp
106                sshb(ji,jj) = rn_wadmin - bathy(ji,jj)
107              END IF
108           ENDDO
109        END DO
110
111     
112        !! start limiter iteration
113        DO jk1 = 1, nn_waditr
114       
115           zflag = 0     ! flag indicating if any further iteration is needed?
116         
117           zflxu1(:,:) = zflxu(:,:) * wdlmt(:,:)
118           zflxv1(:,:) = zflxv(:,:) * wdlmt(:,:)
119         
120           DO jj = 2, jpjm1
121              DO ji = fs_2, fs_jpim1   ! vector opt.
122       
123                 IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE
124       
125                 ztmp = e1t(ji,jj) * e2t(ji,jj)
126
127                 zflxn(ji,jj) = min(zflxu1(ji,jj), 0._wp) - max(zflxu1(ji-1,jj),   0._wp) + &
128                              & min(zflxv1(ji,jj), 0._wp) - max(zflxv1(ji,  jj-1), 0._wp) 
129         
130                 zdep1 = (zflxp(ji,jj) * wdlmt(ji,jj) + zflxn(ji,jj)) * z2dt / ztmp
131                 zdep2 = bathy(ji,jj) + sshb(ji,jj) - rn_wadmin 
132         
133                 IF(zdep1 > zdep2) THEN
134                   zflag = 1
135                   zcoef = ( ( zdep2 - rn_wadmine ) * ztmp - zflxn(ji,jj) * z2dt ) / ( zflxp(ji,jj) * z2dt )
136                   zcoef = max(zcoef, 0._wp)
137                   IF(zflxu1(ji,  jj  ) >= 0._wp) wdlmt(ji,  jj  ) = zcoef
138                   IF(zflxu1(ji-1,jj  ) <  0._wp) wdlmt(ji-1,jj  ) = zcoef
139                   IF(zflxv1(ji,  jj  ) >= 0._wp) wdlmt(ji,  jj  ) = zcoef
140                   IF(zflxv1(ji,  jj-1) <  0._wp) wdlmt(ji,  jj-1) = zcoef
141                 END IF
142              END DO ! ji loop
143           END DO  ! jj loop
144         
145           CALL lbc_lnk( wdlmt, 'T', 1. )
146         
147           IF(zflag == 0) EXIT
148         
149           IF(jk1 == nn_waditr) THEN
150              IF(zflxu1(ji,  jj  ) >= 0._wp) wdlmt(ji,  jj  ) = 0._wp
151              IF(zflxu1(ji-1,jj  ) <  0._wp) wdlmt(ji-1,jj  ) = 0._wp
152              IF(zflxv1(ji,  jj  ) >= 0._wp) wdlmt(ji,  jj  ) = 0._wp
153              IF(zflxv1(ji,  jj-1) <  0._wp) wdlmt(ji,  jj-1) = 0._wp
154           END IF
155
156        END DO  ! jk1 loop
157       
158        DO jk = 1, jpkm1 
159           DO jj = 2, jpjm1
160              DO ji = fs_2, fs_jpim1   ! vector opt.
161               un(ji,  jj,  jk) = ( sign(0.5_wp, un(ji,  jj,  jk))  + 0.5_wp ) * wdlmt(ji  ,jj) 
162               vn(ji,  jj,  jk) = ( sign(0.5_wp, vn(ji,  jj,  jk))  + 0.5_wp ) * wdlmt(ji  ,jj) 
163               un(ji-1,jj,  jk) = (-sign(0.5_wp, un(ji-1,jj,  jk))  + 0.5_wp ) * wdlmt(ji-1,jj) 
164               vn(ji,  jj-1,jk) = (-sign(0.5_wp, vn(ji,  jj-1,jk))  + 0.5_wp ) * wdlmt(ji,  jj-1) 
165              END DO
166           END DO
167        END DO
168       
169        IF(zflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt!!!'
170       
171        !IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )          ! runoffs (update hdivn field)
172        !IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (update hdivn field)
173        !
174        !
175        CALL wrk_dealloc( jpi, jpj, zflxp, zflxn, zflxu, zflxv, zflxu1, zflxv1, wdlmt )
176      !
177      END IF
178
179      IF( nn_timing == 1 )  CALL timing_stop('wad_lmt')
180   END SUBROUTINE wad_lmt
181   !!======================================================================
182END MODULE wadlmt
Note: See TracBrowser for help on using the repository browser.