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

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

add wadlmt.F90 which seems not uploaded correctly last time
File size: 13.7 KB
RevLine 
[5017]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 that ocean cells can have W/D processes
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        un(:,:,:) = un(:,:,:) * zwdlmtu(ji, jj) 
173        vn(:,:,:) = vn(:,:,:) * zwdlmtv(ji, jj) 
174
175        CALL lbc_lnk( un, 'U', -1. )
176        CALL lbc_lnk( vn, 'V', -1. )
177       
178        IF(zflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt!!!'
179       
180        !IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )          ! runoffs (update hdivn field)
181        !IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (update hdivn field)
182        !
183        !
184        CALL wrk_dealloc( jpi, jpj, zflxp, zflxn, zflxu, zflxv, zflxu1, zflxv1 )
185        CALL wrk_dealloc( jpi, jpj, zwdlmtu, zwdlmtv)
186      !
187      END IF
188
189      IF( nn_timing == 1 )  CALL timing_stop('wad_lmt')
190   END SUBROUTINE wad_lmt
191
192   SUBROUTINE wad_lmt_bt( zflxu, zflxv, sshn_e, zssh_frc, rdtbt )
193      !!----------------------------------------------------------------------
194      !!                  ***  ROUTINE wad_lmt  ***
195      !!                   
196      !! ** Purpose :   limiting flux in the barotropic stepping (dynspg_ts)
197      !!
198      !! ** Method  : - Prevent negative depth occurring (Not ready for Agrif)
199      !!
200      !! ** Action  : - calculate flux limiter and W/D flag
201      !!----------------------------------------------------------------------
202      REAL(wp), INTENT(in)    ::   rdtbt    ! ocean time-step index
203      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zflxu,  zflxv, sshn_e, zssh_frc 
204      !
205      INTEGER  ::   ji, jj, jk, jk1     ! dummy loop indices
206      INTEGER  ::   zflag         ! local scalar
207      REAL(wp) ::   z2dt
208      REAL(wp) ::   zcoef, zdep1, zdep2 ! local scalars
209      REAL(wp) ::   zzflxp, zzflxn      ! local scalars
210      REAL(wp) ::   zdepwd              ! local scalar, always wet cell depth
211      REAL(wp) ::   ztmp                ! local scalars
212      REAL(wp), POINTER,  DIMENSION(:,:) ::   zwdlmtu, zwdlmtv         !: W/D flux limiters
213      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxp,  zflxn            ! local 2D workspace
214      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxu1, zflxv1           ! local 2D workspace
215
216      !!----------------------------------------------------------------------
217      !
218
219      IF( nn_timing == 1 )  CALL timing_start('wad_lmt_bt')
220
221      IF(ln_wd) THEN
222
223        CALL wrk_alloc( jpi, jpj, zflxp, zflxn, zflxu1, zflxv1 )
224        CALL wrk_alloc( jpi, jpj, zwdlmtu, zwdlmtv)
225        !
226       
227        IF(lwp) WRITE(numout,*)
228        IF(lwp) WRITE(numout,*) 'wad_lmt_bt : wetting/drying limiters and velocity limiting'
229       
230        zflag  = 0
231        zdepwd = 50._wp   !maximum depth that ocean cells can have W/D processes
232
233        z2dt = rdtbt   
234       
235        zflxp(:,:)   = 0._wp
236        zflxn(:,:)   = 0._wp
237        zflxu(:,:)   = 0._wp
238        zflxv(:,:)   = 0._wp
239
240        zwdlmtu(:,:)  = 1._wp
241        zwdlmtv(:,:)  = 1._wp
242       
243        ! Horizontal Flux in u and v direction
244       
245        zflxu(:,:) = zflxu(:,:) * e2u(:,:)
246        zflxv(:,:) = zflxv(:,:) * e1v(:,:)
247       
248        DO jj = 2, jpjm1
249           DO ji = 2, jpim1 
250
251             IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE   ! we don't care about land cells
252             IF(bathy(ji,jj) > zdepwd) CYCLE       ! and cells which will unlikely go dried out
253
254              zflxp(ji,jj) = max(zflxu(ji,jj), 0._wp) - min(zflxu(ji-1,jj),   0._wp) + &
255                           & max(zflxv(ji,jj), 0._wp) - min(zflxv(ji,  jj-1), 0._wp) 
256              zflxn(ji,jj) = min(zflxu(ji,jj), 0._wp) - max(zflxu(ji-1,jj),   0._wp) + &
257                           & min(zflxv(ji,jj), 0._wp) - max(zflxv(ji,  jj-1), 0._wp) 
258
259              zdep2 = bathy(ji,jj) + sshn_e(ji,jj) - rn_wdmin1
260              IF(zdep2 < 0._wp) THEN  !add more safty, but not necessary
261                !zdep2 = 0._wp
262               sshn_e(ji,jj) = rn_wdmin1 - bathy(ji,jj)
263              END IF
264           ENDDO
265        END DO
266
267     
268        !! start limiter iterations
269        DO jk1 = 1, nn_wdit + 1
270       
271         
272           zflxu1(:,:) = zflxu(:,:) * zwdlmtu(:,:)
273           zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:)
274         
275           DO jj = 2, jpjm1
276              DO ji = 2, jpim1 
277       
278                 wdmask(ji,jj) = 0
279                 IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE
280                 IF(bathy(ji,jj) > zdepwd) CYCLE
281       
282                 !ztmp = e1t(ji,jj) * e2t(ji,jj)     !there must be an array ready for this
283                 ztmp = e12t(ji,jj)
284
285                 zzflxp = max(zflxu1(ji,jj), 0._wp) - min(zflxu1(ji-1,jj),   0._wp) + &
286                        & max(zflxv1(ji,jj), 0._wp) - min(zflxv1(ji,  jj-1), 0._wp) 
287                 zzflxn = min(zflxu1(ji,jj), 0._wp) - max(zflxu1(ji-1,jj),   0._wp) + &
288                        & min(zflxv1(ji,jj), 0._wp) - max(zflxv1(ji,  jj-1), 0._wp) 
289         
290                 zdep1 = (zzflxp + zzflxn) * z2dt / ztmp
291                 zdep2 = bathy(ji,jj) + sshn_e(ji,jj) - rn_wdmin1   ! this one can be moved out of the loop
292                 zdep2 = zdep2 - z2dt * zssh_frc(ji,jj)
293         
294                 IF(zdep1 > zdep2) THEN
295                   zflag = 1
296                   !wdmask(ji, jj) = 1
297                   zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt )
298                   zcoef = max(zcoef, 0._wp)
299                   IF(jk1 > nn_wdit) zcoef = 0._wp
300                   IF(zflxu1(ji,  jj) > 0._wp) zwdlmtu(ji  ,jj) = zcoef
301                   IF(zflxu1(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = zcoef
302                   IF(zflxv1(ji,  jj) > 0._wp) zwdlmtv(ji  ,jj) = zcoef
303                   IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji-1,jj) = zcoef
304                 END IF
305              END DO ! ji loop
306           END DO  ! jj loop
307
308           CALL lbc_lnk( zwdlmtu, 'U', 1. )
309           CALL lbc_lnk( zwdlmtv, 'V', 1. )
310
311           IF(lk_mpp) CALL mpp_max(zflag)   !max over the global domain
312
313           IF(zflag == 0) EXIT
314         
315           zflag = 0     ! flag indicating if any further iteration is needed?
316        END DO  ! jk1 loop
317       
318        zflxu(:,:) = zflxu(:,:) * zwdlmtu(ji, jj) 
319        zflxv(:,:) = zflxv(:,:) * zwdlmtv(ji, jj) 
320
321        CALL lbc_lnk( zflxu, 'U', -1. )
322        CALL lbc_lnk( zflxv, 'V', -1. )
323       
324        IF(zflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt_bt!!!'
325       
326        !IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )          ! runoffs (update hdivn field)
327        !IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (update hdivn field)
328        !
329        !
330        CALL wrk_dealloc( jpi, jpj, zflxp, zflxn, zflxu1, zflxv1 )
331        CALL wrk_dealloc( jpi, jpj, zwdlmtu, zwdlmtv)
332      !
333      END IF
334
335      IF( nn_timing == 1 )  CALL timing_stop('wad_lmt')
336   END SUBROUTINE wad_lmt_bt
337   !!======================================================================
338END MODULE wadlmt
Note: See TracBrowser for help on using the repository browser.