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/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/wadlmt.F90 @ 5520

Last change on this file since 5520 was 5520, checked in by hliu, 9 years ago

removed two bugs, and some redundant priniting out info for Wetting/Drying?
File size: 13.7 KB
Line 
1
2MODULE wadlmt
3   !!==============================================================================
4   !!                       ***  MODULE  wadlmt  ***
5   !! compute and apply flux limiters and preserve water depth positivity
6   !! only effects if wetting/drying is on (ln_wd == .true.
7   !!==============================================================================
8   !! History :   
9   !!  NEMO      3.6  ! 2014-09  ((H.Liu)  Original code
10   !!                 ! will add the runoff and periodic BC case later
11   !!----------------------------------------------------------------------
12
13   !!----------------------------------------------------------------------
14   !!   wad_lmt    : Compute the horizontal flux limiter and the limited velocity
15   !!                when wetting and drying happens
16   !!----------------------------------------------------------------------
17   USE oce             ! ocean dynamics and tracers
18   USE dom_oce         ! ocean space and time domain
19   USE sbc_oce, ONLY : ln_rnf   ! surface boundary condition: ocean
20   USE sbcrnf          ! river runoff
21   USE cla             ! cross land advection             (cla_div routine)
22   USE in_out_manager  ! I/O manager
23   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
24   USE lib_mpp         ! MPP library
25   USE wrk_nemo        ! Memory Allocation
26   USE timing          ! Timing
27
28   IMPLICIT NONE
29   PRIVATE
30
31
32   PUBLIC   wad_lmt, wad_lmt_bt    ! routine called by step.F90
33
34   !! * Substitutions
35#  include "domzgr_substitute.h90"
36#  include "vectopt_loop_substitute.h90"
37CONTAINS
38
39   SUBROUTINE wad_lmt( sshb1, sshemp, z2dt )
40      !!----------------------------------------------------------------------
41      !!                  ***  ROUTINE wad_lmt  ***
42      !!                   
43      !! ** Purpose :   generate flux limiters for wetting/drying
44      !!
45      !! ** Method  : - Prevent negative depth occurring (Not ready for Agrif)
46      !!
47      !! ** Action  : - calculate flux limiter and W/D flag
48      !!----------------------------------------------------------------------
49      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   sshb1
50      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   sshemp
51      REAL(wp), INTENT(in) :: z2dt
52      !
53      INTEGER  ::   ji, jj, jk, jk1     ! dummy loop indices
54      INTEGER  ::   zflag               ! local scalar
55      REAL(wp) ::   zcoef, zdep1, zdep2 ! local scalars
56      REAL(wp) ::   zzflxp, zzflxn      ! local scalars
57      REAL(wp) ::   zdepwd              ! local scalar, always wet cell depth
58      REAL(wp) ::   ztmp                ! local scalars
59      REAL(wp), POINTER,  DIMENSION(:,:) ::   zwdlmtu, zwdlmtv         !: W/D flux limiters
60      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxp,  zflxn            ! local 2D workspace
61      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxu,  zflxv            ! local 2D workspace
62      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxu1, zflxv1           ! local 2D workspace
63
64      !!----------------------------------------------------------------------
65      !
66
67      IF( nn_timing == 1 )  CALL timing_start('wad_lmt')
68
69      IF(ln_wd) THEN
70
71        CALL wrk_alloc( jpi, jpj, zflxp, zflxn, zflxu, zflxv, zflxu1, zflxv1 )
72        CALL wrk_alloc( jpi, jpj, zwdlmtu, zwdlmtv)
73        !
74       
75        !IF(lwp) WRITE(numout,*)
76        !IF(lwp) WRITE(numout,*) 'wad_lmt : wetting/drying limiters and velocity limiting'
77       
78        zflag  = 0
79        zdepwd = 50._wp   !maximum depth on which that W/D could possibly happen
80
81       
82        zflxp(:,:)   = 0._wp
83        zflxn(:,:)   = 0._wp
84        zflxu(:,:)   = 0._wp
85        zflxv(:,:)   = 0._wp
86
87        zwdlmtu(:,:)  = 1._wp
88        zwdlmtv(:,:)  = 1._wp
89       
90        ! Horizontal Flux in u and v direction
91        DO jk = 1, jpkm1 
92           DO jj = 1, jpjm1
93              DO ji = 1, jpim1
94                 zflxu(ji,jj) = zflxu(ji,jj) + fse3u(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)
95                 zflxv(ji,jj) = zflxv(ji,jj) + fse3v(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
96              END DO 
97           END DO 
98        END DO
99       
100        zflxu(:,:) = zflxu(:,:) * e2u(:,:)
101        zflxv(:,:) = zflxv(:,:) * e1v(:,:)
102       
103        DO jj = 2, jpjm1
104           DO ji = 2, jpim1 
105
106             IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE   ! we don't care about land cells
107             IF(bathy(ji,jj) > zdepwd) CYCLE       ! and cells which will unlikely go dried out
108
109              zflxp(ji,jj) = max(zflxu(ji,jj), 0._wp) - min(zflxu(ji-1,jj),   0._wp) + &
110                           & max(zflxv(ji,jj), 0._wp) - min(zflxv(ji,  jj-1), 0._wp) 
111              zflxn(ji,jj) = min(zflxu(ji,jj), 0._wp) - max(zflxu(ji-1,jj),   0._wp) + &
112                           & min(zflxv(ji,jj), 0._wp) - max(zflxv(ji,  jj-1), 0._wp) 
113
114              zdep2 = bathy(ji,jj) + sshb1(ji,jj) - rn_wdmin1
115              IF(zdep2 < 0._wp) THEN  !add more safty, but not necessary
116                !zdep2 = 0._wp
117                sshb1(ji,jj) = rn_wdmin1 - bathy(ji,jj)
118              END IF
119           ENDDO
120        END DO
121
122     
123        !! start limiter iterations
124        DO jk1 = 1, nn_wdit + 1
125       
126         
127           zflxu1(:,:) = zflxu(:,:) * zwdlmtu(:,:)
128           zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:)
129         
130           DO jj = 2, jpjm1
131              DO ji = 2, jpim1 
132       
133                 wdmask(ji,jj) = 0
134                 IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE
135                 IF(bathy(ji,jj) > zdepwd) CYCLE
136       
137                 !ztmp = e1t(ji,jj) * e2t(ji,jj)     !there must be an array ready for this
138                 ztmp = e12t(ji,jj)
139
140                 zzflxp = max(zflxu1(ji,jj), 0._wp) - min(zflxu1(ji-1,jj),   0._wp) + &
141                        & max(zflxv1(ji,jj), 0._wp) - min(zflxv1(ji,  jj-1), 0._wp) 
142                 zzflxn = min(zflxu1(ji,jj), 0._wp) - max(zflxu1(ji-1,jj),   0._wp) + &
143                        & min(zflxv1(ji,jj), 0._wp) - max(zflxv1(ji,  jj-1), 0._wp) 
144         
145                 zdep1 = (zzflxp + zzflxn) * z2dt / ztmp
146                 zdep2 = bathy(ji,jj) + sshb1(ji,jj) - rn_wdmin1 - z2dt * sshemp(ji,jj)  ! this one can be moved out of the loop
147         
148                 IF(zdep1 > zdep2) THEN
149                   zflag = 1
150                   wdmask(ji, jj) = 1
151                   zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt )
152                   zcoef = max(zcoef, 0._wp)
153                   IF(jk1 > nn_wdit) zcoef = 0._wp
154                   IF(zflxu1(ji,  jj) > 0._wp) zwdlmtu(ji  ,jj) = zcoef
155                   IF(zflxu1(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = zcoef
156                   IF(zflxv1(ji,  jj) > 0._wp) zwdlmtv(ji  ,jj) = zcoef
157                   IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji-1,jj) = zcoef
158                 END IF
159              END DO ! ji loop
160           END DO  ! jj loop
161
162           CALL lbc_lnk( zwdlmtu, 'U', 1. )
163           CALL lbc_lnk( zwdlmtv, 'V', 1. )
164
165           IF(lk_mpp) CALL mpp_max(zflag)   !max over the global domain
166
167           IF(zflag == 0) EXIT
168         
169           zflag = 0     ! flag indicating if any further iteration is needed?
170        END DO  ! jk1 loop
171       
172        DO jk = 1, jpkm1
173          un(:,:,jk) = un(:,:,jk) * zwdlmtu(:, :) 
174          vn(:,:,jk) = vn(:,:,jk) * zwdlmtv(:, :) 
175        END DO
176
177        CALL lbc_lnk( un, 'U', -1. )
178        CALL lbc_lnk( vn, 'V', -1. )
179       
180        IF(zflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt!!!'
181       
182        !IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )          ! runoffs (update hdivn field)
183        !IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (update hdivn field)
184        !
185        !
186        CALL wrk_dealloc( jpi, jpj, zflxp, zflxn, zflxu, zflxv, zflxu1, zflxv1 )
187        CALL wrk_dealloc( jpi, jpj, zwdlmtu, zwdlmtv)
188      !
189      END IF
190
191      IF( nn_timing == 1 )  CALL timing_stop('wad_lmt')
192   END SUBROUTINE wad_lmt
193
194   SUBROUTINE wad_lmt_bt( zflxu, zflxv, sshn_e, zssh_frc, rdtbt )
195      !!----------------------------------------------------------------------
196      !!                  ***  ROUTINE wad_lmt  ***
197      !!                   
198      !! ** Purpose :   limiting flux in the barotropic stepping (dynspg_ts)
199      !!
200      !! ** Method  : - Prevent negative depth occurring (Not ready for Agrif)
201      !!
202      !! ** Action  : - calculate flux limiter and W/D flag
203      !!----------------------------------------------------------------------
204      REAL(wp), INTENT(in)    ::   rdtbt    ! ocean time-step index
205      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zflxu,  zflxv, sshn_e, zssh_frc 
206      !
207      INTEGER  ::   ji, jj, jk, jk1     ! dummy loop indices
208      INTEGER  ::   zflag         ! local scalar
209      REAL(wp) ::   z2dt
210      REAL(wp) ::   zcoef, zdep1, zdep2 ! local scalars
211      REAL(wp) ::   zzflxp, zzflxn      ! local scalars
212      REAL(wp) ::   zdepwd              ! local scalar, always wet cell depth
213      REAL(wp) ::   ztmp                ! local scalars
214      REAL(wp), POINTER,  DIMENSION(:,:) ::   zwdlmtu, zwdlmtv         !: W/D flux limiters
215      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxp,  zflxn            ! local 2D workspace
216      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxu1, zflxv1           ! local 2D workspace
217
218      !!----------------------------------------------------------------------
219      !
220
221      IF( nn_timing == 1 )  CALL timing_start('wad_lmt_bt')
222
223      IF(ln_wd) THEN
224
225        CALL wrk_alloc( jpi, jpj, zflxp, zflxn, zflxu1, zflxv1 )
226        CALL wrk_alloc( jpi, jpj, zwdlmtu, zwdlmtv)
227        !
228       
229        !IF(lwp) WRITE(numout,*)
230        !IF(lwp) WRITE(numout,*) 'wad_lmt_bt : wetting/drying limiters and velocity limiting'
231       
232        zflag  = 0
233        zdepwd = 50._wp   !maximum depth that ocean cells can have W/D processes
234
235        z2dt = rdtbt   
236       
237        zflxp(:,:)   = 0._wp
238        zflxn(:,:)   = 0._wp
239        !zflxu(:,:)   = 0._wp
240        !zflxv(:,:)   = 0._wp
241
242        zwdlmtu(:,:)  = 1._wp
243        zwdlmtv(:,:)  = 1._wp
244       
245        ! Horizontal Flux in u and v direction
246       
247        !zflxu(:,:) = zflxu(:,:) * e2u(:,:)
248        !zflxv(:,:) = zflxv(:,:) * e1v(:,:)
249       
250        DO jj = 2, jpjm1
251           DO ji = 2, jpim1 
252
253             IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE   ! we don't care about land cells
254             IF(bathy(ji,jj) > zdepwd) CYCLE       ! and cells which will unlikely go dried out
255
256              zflxp(ji,jj) = max(zflxu(ji,jj), 0._wp) - min(zflxu(ji-1,jj),   0._wp) + &
257                           & max(zflxv(ji,jj), 0._wp) - min(zflxv(ji,  jj-1), 0._wp) 
258              zflxn(ji,jj) = min(zflxu(ji,jj), 0._wp) - max(zflxu(ji-1,jj),   0._wp) + &
259                           & min(zflxv(ji,jj), 0._wp) - max(zflxv(ji,  jj-1), 0._wp) 
260
261              zdep2 = bathy(ji,jj) + sshn_e(ji,jj) - rn_wdmin1
262              IF(zdep2 < 0._wp) THEN  !add more safty, but not necessary
263                !zdep2 = 0._wp
264               sshn_e(ji,jj) = rn_wdmin1 - bathy(ji,jj)
265              END IF
266           ENDDO
267        END DO
268
269     
270        !! start limiter iterations
271        DO jk1 = 1, nn_wdit + 1
272       
273         
274           zflxu1(:,:) = zflxu(:,:) * zwdlmtu(:,:)
275           zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:)
276         
277           DO jj = 2, jpjm1
278              DO ji = 2, jpim1 
279       
280                 wdmask(ji,jj) = 0
281                 IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE
282                 IF(bathy(ji,jj) > zdepwd) CYCLE
283       
284                 !ztmp = e1t(ji,jj) * e2t(ji,jj)     !there must be an array ready for this
285                 ztmp = e12t(ji,jj)
286
287                 zzflxp = max(zflxu1(ji,jj), 0._wp) - min(zflxu1(ji-1,jj),   0._wp) + &
288                        & max(zflxv1(ji,jj), 0._wp) - min(zflxv1(ji,  jj-1), 0._wp) 
289                 zzflxn = min(zflxu1(ji,jj), 0._wp) - max(zflxu1(ji-1,jj),   0._wp) + &
290                        & min(zflxv1(ji,jj), 0._wp) - max(zflxv1(ji,  jj-1), 0._wp) 
291         
292                 zdep1 = (zzflxp + zzflxn) * z2dt / ztmp
293                 zdep2 = bathy(ji,jj) + sshn_e(ji,jj) - rn_wdmin1   ! this one can be moved out of the loop
294                 zdep2 = zdep2 - z2dt * zssh_frc(ji,jj)
295         
296                 IF(zdep1 > zdep2) THEN
297                   zflag = 1
298                   !wdmask(ji, jj) = 1
299                   zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt )
300                   zcoef = max(zcoef, 0._wp)
301                   IF(jk1 > nn_wdit) zcoef = 0._wp
302                   IF(zflxu1(ji,  jj) > 0._wp) zwdlmtu(ji  ,jj) = zcoef
303                   IF(zflxu1(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = zcoef
304                   IF(zflxv1(ji,  jj) > 0._wp) zwdlmtv(ji  ,jj) = zcoef
305                   IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji-1,jj) = zcoef
306                 END IF
307              END DO ! ji loop
308           END DO  ! jj loop
309
310           CALL lbc_lnk( zwdlmtu, 'U', 1. )
311           CALL lbc_lnk( zwdlmtv, 'V', 1. )
312
313           IF(lk_mpp) CALL mpp_max(zflag)   !max over the global domain
314
315           IF(zflag == 0) EXIT
316         
317           zflag = 0     ! flag indicating if any further iteration is needed?
318        END DO  ! jk1 loop
319       
320        zflxu(:,:) = zflxu(:,:) * zwdlmtu(:, :) 
321        zflxv(:,:) = zflxv(:,:) * zwdlmtv(:, :) 
322
323        CALL lbc_lnk( zflxu, 'U', -1. )
324        CALL lbc_lnk( zflxv, 'V', -1. )
325       
326        IF(zflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt_bt!!!'
327       
328        !IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )          ! runoffs (update hdivn field)
329        !IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (update hdivn field)
330        !
331        !
332        CALL wrk_dealloc( jpi, jpj, zflxp, zflxn, zflxu1, zflxv1 )
333        CALL wrk_dealloc( jpi, jpj, zwdlmtu, zwdlmtv)
334      !
335      END IF
336
337      IF( nn_timing == 1 )  CALL timing_stop('wad_lmt')
338   END SUBROUTINE wad_lmt_bt
339   !!======================================================================
340END MODULE wadlmt
Note: See TracBrowser for help on using the repository browser.