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.
wet_dry.F90 in branches/NERC/dev_r5549_BDY_ZEROGRAD/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/NERC/dev_r5549_BDY_ZEROGRAD/NEMOGCM/NEMO/OPA_SRC/DYN/wet_dry.F90 @ 6808

Last change on this file since 6808 was 6808, checked in by jamesharle, 8 years ago

merge with trunk@6232 for consistency with SSB code
File size: 16.7 KB
Line 
1
2MODULE wet_dry
3   !!==============================================================================
4   !!                       ***  MODULE  wet_dry  ***
5   !! Wetting and drying includes initialisation routine and routines to
6   !! compute and apply flux limiters and preserve water depth positivity
7   !! only effects if wetting/drying is on (ln_wd == .true.)
8   !!==============================================================================
9   !! History :   
10   !!  NEMO      3.6  ! 2014-09  ((H.Liu)  Original code
11   !!                 ! will add the runoff and periodic BC case later
12   !!----------------------------------------------------------------------
13
14   !!----------------------------------------------------------------------
15   !!   wad_lmt    : Compute the horizontal flux limiter and the limited velocity
16   !!                when wetting and drying happens
17   !!----------------------------------------------------------------------
18   USE oce             ! ocean dynamics and tracers
19   USE dom_oce         ! ocean space and time domain
20   USE sbc_oce, ONLY : ln_rnf   ! surface boundary condition: ocean
21   USE sbcrnf          ! river runoff
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   !! critical depths,filters, limiters,and masks for  Wetting and Drying
33   !! ---------------------------------------------------------------------
34
35   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) ::   wduflt, wdvflt !: u- and v- filter
36   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) ::   wdmask         !: u- and v- limiter
37
38   LOGICAL,  PUBLIC  ::   ln_wd       !: Wetting/drying activation switch (T:on,F:off)
39   REAL(wp), PUBLIC  ::   rn_wdmin1   !: minimum water depth on dried cells
40   REAL(wp), PUBLIC  ::   rn_wdmin2   !: tolerrance of minimum water depth on dried cells
41   REAL(wp), PUBLIC  ::   rn_wdld     !: land elevation below which wetting/drying
42                                           !: will be considered
43   INTEGER , PUBLIC  ::   nn_wdit     !: maximum number of iteration for W/D limiter
44
45   PUBLIC   wad_init                  ! initialisation routine called by step.F90
46   PUBLIC   wad_lmt                   ! routine called by sshwzv.F90
47   PUBLIC   wad_lmt_bt                ! routine called by dynspg_ts.F90
48
49   !! * Substitutions
50#  include "vectopt_loop_substitute.h90"
51CONTAINS
52
53   SUBROUTINE wad_init
54      !!----------------------------------------------------------------------
55      !!                     ***  ROUTINE wad_init  ***
56      !!                   
57      !! ** Purpose :   read wetting and drying namelist and print the variables.
58      !!
59      !! ** input   : - namwad namelist
60      !!----------------------------------------------------------------------
61      NAMELIST/namwad/ ln_wd, rn_wdmin1, rn_wdmin2, rn_wdld, nn_wdit
62      INTEGER  ::   ios                 ! Local integer output status for namelist read
63      INTEGER  ::   ierr                ! Local integer status array allocation
64      !!----------------------------------------------------------------------
65
66      REWIND( numnam_ref )              ! Namelist namwad in reference namelist
67                                        ! : Parameters for Wetting/Drying
68      READ  ( numnam_ref, namwad, IOSTAT = ios, ERR = 905)
69905   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namwad in reference namelist', .TRUE.) 
70      REWIND( numnam_cfg )              ! Namelist namwad in configuration namelist
71                                        ! : Parameters for Wetting/Drying
72      READ  ( numnam_cfg, namwad, IOSTAT = ios, ERR = 906)
73906   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namwad in configuration namelist', .TRUE. )
74      IF(lwm) WRITE ( numond, namwad )
75
76      IF(lwp) THEN                  ! control print
77         WRITE(numout,*)
78         WRITE(numout,*) 'wad_init : Wetting and drying initialization through namelist read'
79         WRITE(numout,*) '~~~~~~~ '
80         WRITE(numout,*) '   Namelist namwad'
81         WRITE(numout,*) '      Logical activation                 ln_wd      = ', ln_wd
82         WRITE(numout,*) '      Minimum wet depth on dried cells rn_wdmin1    = ', rn_wdmin1
83         WRITE(numout,*) '      Tolerance of min wet depth     rn_wdmin2      = ', rn_wdmin2
84         WRITE(numout,*) '      land elevation threshold         rn_wdld      = ', rn_wdld
85         WRITE(numout,*) '      Max iteration for W/D limiter    nn_wdit      = ', nn_wdit
86       ENDIF
87
88      IF(ln_wd) THEN
89         ALLOCATE( wduflt(jpi,jpj), wdvflt(jpi,jpj), wdmask(jpi,jpj), STAT=ierr )
90         IF( ierr /= 0 ) CALL ctl_stop('STOP', 'wad_init : Array allocation error')
91      ENDIF
92   END SUBROUTINE wad_init
93
94   SUBROUTINE wad_lmt( sshb1, sshemp, z2dt )
95      !!----------------------------------------------------------------------
96      !!                  ***  ROUTINE wad_lmt  ***
97      !!                   
98      !! ** Purpose :   generate flux limiters for wetting/drying
99      !!
100      !! ** Method  : - Prevent negative depth occurring (Not ready for Agrif)
101      !!
102      !! ** Action  : - calculate flux limiter and W/D flag
103      !!----------------------------------------------------------------------
104      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   sshb1
105      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   sshemp
106      REAL(wp), INTENT(in) :: z2dt
107      !
108      INTEGER  ::   ji, jj, jk, jk1     ! dummy loop indices
109      INTEGER  ::   zflag               ! local scalar
110      REAL(wp) ::   zcoef, zdep1, zdep2 ! local scalars
111      REAL(wp) ::   zzflxp, zzflxn      ! local scalars
112      REAL(wp) ::   zdepwd              ! local scalar, always wet cell depth
113      REAL(wp) ::   ztmp                ! local scalars
114      REAL(wp), POINTER,  DIMENSION(:,:) ::   zwdlmtu, zwdlmtv         !: W/D flux limiters
115      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxp,  zflxn            ! local 2D workspace
116      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxu,  zflxv            ! local 2D workspace
117      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxu1, zflxv1           ! local 2D workspace
118
119      !!----------------------------------------------------------------------
120      !
121
122      IF( nn_timing == 1 )  CALL timing_start('wad_lmt')
123
124      IF(ln_wd) THEN
125
126        CALL wrk_alloc( jpi, jpj, zflxp, zflxn, zflxu, zflxv, zflxu1, zflxv1 )
127        CALL wrk_alloc( jpi, jpj, zwdlmtu, zwdlmtv)
128        !
129       
130        !IF(lwp) WRITE(numout,*)
131        !IF(lwp) WRITE(numout,*) 'wad_lmt : wetting/drying limiters and velocity limiting'
132       
133        zflag  = 0
134        zdepwd = 50._wp   !maximum depth on which that W/D could possibly happen
135
136       
137        zflxp(:,:)   = 0._wp
138        zflxn(:,:)   = 0._wp
139        zflxu(:,:)   = 0._wp
140        zflxv(:,:)   = 0._wp
141
142        zwdlmtu(:,:)  = 1._wp
143        zwdlmtv(:,:)  = 1._wp
144       
145        ! Horizontal Flux in u and v direction
146        DO jk = 1, jpkm1 
147           DO jj = 1, jpjm1
148              DO ji = 1, jpim1
149                 zflxu(ji,jj) = zflxu(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)
150                 zflxv(ji,jj) = zflxv(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
151              END DO 
152           END DO 
153        END DO
154       
155        zflxu(:,:) = zflxu(:,:) * e2u(:,:)
156        zflxv(:,:) = zflxv(:,:) * e1v(:,:)
157       
158        DO jj = 2, jpjm1
159           DO ji = 2, jpim1 
160
161             IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE   ! we don't care about land cells
162             IF(bathy(ji,jj) > zdepwd) CYCLE       ! and cells which will unlikely go dried out
163
164              zflxp(ji,jj) = max(zflxu(ji,jj), 0._wp) - min(zflxu(ji-1,jj),   0._wp) + &
165                           & max(zflxv(ji,jj), 0._wp) - min(zflxv(ji,  jj-1), 0._wp) 
166              zflxn(ji,jj) = min(zflxu(ji,jj), 0._wp) - max(zflxu(ji-1,jj),   0._wp) + &
167                           & min(zflxv(ji,jj), 0._wp) - max(zflxv(ji,  jj-1), 0._wp) 
168
169              zdep2 = bathy(ji,jj) + sshb1(ji,jj) - rn_wdmin1
170              IF(zdep2 < 0._wp) THEN  !add more safty, but not necessary
171                !zdep2 = 0._wp
172                sshb1(ji,jj) = rn_wdmin1 - bathy(ji,jj)
173              END IF
174           ENDDO
175        END DO
176
177     
178        !! start limiter iterations
179        DO jk1 = 1, nn_wdit + 1
180       
181         
182           zflxu1(:,:) = zflxu(:,:) * zwdlmtu(:,:)
183           zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:)
184         
185           DO jj = 2, jpjm1
186              DO ji = 2, jpim1 
187       
188                 wdmask(ji,jj) = 0
189                 IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE
190                 IF(bathy(ji,jj) > zdepwd) CYCLE
191       
192                 ztmp = e1e2t(ji,jj)
193
194                 zzflxp = max(zflxu1(ji,jj), 0._wp) - min(zflxu1(ji-1,jj),   0._wp) + &
195                        & max(zflxv1(ji,jj), 0._wp) - min(zflxv1(ji,  jj-1), 0._wp) 
196                 zzflxn = min(zflxu1(ji,jj), 0._wp) - max(zflxu1(ji-1,jj),   0._wp) + &
197                        & min(zflxv1(ji,jj), 0._wp) - max(zflxv1(ji,  jj-1), 0._wp) 
198         
199                 zdep1 = (zzflxp + zzflxn) * z2dt / ztmp
200                 zdep2 = bathy(ji,jj) + sshb1(ji,jj) - rn_wdmin1 - z2dt * sshemp(ji,jj)  ! this one can be moved out of the loop
201         
202                 IF(zdep1 > zdep2) THEN
203                   zflag = 1
204                   wdmask(ji, jj) = 1
205                   zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt )
206                   zcoef = max(zcoef, 0._wp)
207                   IF(jk1 > nn_wdit) zcoef = 0._wp
208                   IF(zflxu1(ji,  jj) > 0._wp) zwdlmtu(ji  ,jj) = zcoef
209                   IF(zflxu1(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = zcoef
210                   IF(zflxv1(ji,  jj) > 0._wp) zwdlmtv(ji  ,jj) = zcoef
211                   IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji-1,jj) = zcoef
212                 END IF
213              END DO ! ji loop
214           END DO  ! jj loop
215
216           CALL lbc_lnk( zwdlmtu, 'U', 1. )
217           CALL lbc_lnk( zwdlmtv, 'V', 1. )
218
219           IF(lk_mpp) CALL mpp_max(zflag)   !max over the global domain
220
221           IF(zflag == 0) EXIT
222         
223           zflag = 0     ! flag indicating if any further iteration is needed?
224        END DO  ! jk1 loop
225       
226        DO jk = 1, jpkm1
227          un(:,:,jk) = un(:,:,jk) * zwdlmtu(:, :) 
228          vn(:,:,jk) = vn(:,:,jk) * zwdlmtv(:, :) 
229        END DO
230
231        CALL lbc_lnk( un, 'U', -1. )
232        CALL lbc_lnk( vn, 'V', -1. )
233       
234        IF(zflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt!!!'
235       
236        !IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )          ! runoffs (update hdivn field)
237        !IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (update hdivn field)
238        !
239        !
240        CALL wrk_dealloc( jpi, jpj, zflxp, zflxn, zflxu, zflxv, zflxu1, zflxv1 )
241        CALL wrk_dealloc( jpi, jpj, zwdlmtu, zwdlmtv)
242      !
243      END IF
244
245      IF( nn_timing == 1 )  CALL timing_stop('wad_lmt')
246   END SUBROUTINE wad_lmt
247
248   SUBROUTINE wad_lmt_bt( zflxu, zflxv, sshn_e, zssh_frc, rdtbt )
249      !!----------------------------------------------------------------------
250      !!                  ***  ROUTINE wad_lmt  ***
251      !!                   
252      !! ** Purpose :   limiting flux in the barotropic stepping (dynspg_ts)
253      !!
254      !! ** Method  : - Prevent negative depth occurring (Not ready for Agrif)
255      !!
256      !! ** Action  : - calculate flux limiter and W/D flag
257      !!----------------------------------------------------------------------
258      REAL(wp), INTENT(in)    ::   rdtbt    ! ocean time-step index
259      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zflxu,  zflxv, sshn_e, zssh_frc 
260      !
261      INTEGER  ::   ji, jj, jk, jk1     ! dummy loop indices
262      INTEGER  ::   zflag         ! local scalar
263      REAL(wp) ::   z2dt
264      REAL(wp) ::   zcoef, zdep1, zdep2 ! local scalars
265      REAL(wp) ::   zzflxp, zzflxn      ! local scalars
266      REAL(wp) ::   zdepwd              ! local scalar, always wet cell depth
267      REAL(wp) ::   ztmp                ! local scalars
268      REAL(wp), POINTER,  DIMENSION(:,:) ::   zwdlmtu, zwdlmtv         !: W/D flux limiters
269      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxp,  zflxn            ! local 2D workspace
270      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxu1, zflxv1           ! local 2D workspace
271
272      !!----------------------------------------------------------------------
273      !
274
275      IF( nn_timing == 1 )  CALL timing_start('wad_lmt_bt')
276
277      IF(ln_wd) THEN
278
279        CALL wrk_alloc( jpi, jpj, zflxp, zflxn, zflxu1, zflxv1 )
280        CALL wrk_alloc( jpi, jpj, zwdlmtu, zwdlmtv)
281        !
282       
283        !IF(lwp) WRITE(numout,*)
284        !IF(lwp) WRITE(numout,*) 'wad_lmt_bt : wetting/drying limiters and velocity limiting'
285       
286        zflag  = 0
287        zdepwd = 50._wp   !maximum depth that ocean cells can have W/D processes
288
289        z2dt = rdtbt   
290       
291        zflxp(:,:)   = 0._wp
292        zflxn(:,:)   = 0._wp
293        !zflxu(:,:)   = 0._wp
294        !zflxv(:,:)   = 0._wp
295
296        zwdlmtu(:,:)  = 1._wp
297        zwdlmtv(:,:)  = 1._wp
298       
299        ! Horizontal Flux in u and v direction
300       
301        !zflxu(:,:) = zflxu(:,:) * e2u(:,:)
302        !zflxv(:,:) = zflxv(:,:) * e1v(:,:)
303       
304        DO jj = 2, jpjm1
305           DO ji = 2, jpim1 
306
307             IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE   ! we don't care about land cells
308             IF(bathy(ji,jj) > zdepwd) CYCLE       ! and cells which will unlikely go dried out
309
310              zflxp(ji,jj) = max(zflxu(ji,jj), 0._wp) - min(zflxu(ji-1,jj),   0._wp) + &
311                           & max(zflxv(ji,jj), 0._wp) - min(zflxv(ji,  jj-1), 0._wp) 
312              zflxn(ji,jj) = min(zflxu(ji,jj), 0._wp) - max(zflxu(ji-1,jj),   0._wp) + &
313                           & min(zflxv(ji,jj), 0._wp) - max(zflxv(ji,  jj-1), 0._wp) 
314
315              zdep2 = bathy(ji,jj) + sshn_e(ji,jj) - rn_wdmin1
316              IF(zdep2 < 0._wp) THEN  !add more safty, but not necessary
317                !zdep2 = 0._wp
318               sshn_e(ji,jj) = rn_wdmin1 - bathy(ji,jj)
319              END IF
320           ENDDO
321        END DO
322
323     
324        !! start limiter iterations
325        DO jk1 = 1, nn_wdit + 1
326       
327         
328           zflxu1(:,:) = zflxu(:,:) * zwdlmtu(:,:)
329           zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:)
330         
331           DO jj = 2, jpjm1
332              DO ji = 2, jpim1 
333       
334                 wdmask(ji,jj) = 0
335                 IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE
336                 IF(bathy(ji,jj) > zdepwd) CYCLE
337       
338                 ztmp = e1e2t(ji,jj)
339
340                 zzflxp = max(zflxu1(ji,jj), 0._wp) - min(zflxu1(ji-1,jj),   0._wp) + &
341                        & max(zflxv1(ji,jj), 0._wp) - min(zflxv1(ji,  jj-1), 0._wp) 
342                 zzflxn = min(zflxu1(ji,jj), 0._wp) - max(zflxu1(ji-1,jj),   0._wp) + &
343                        & min(zflxv1(ji,jj), 0._wp) - max(zflxv1(ji,  jj-1), 0._wp) 
344         
345                 zdep1 = (zzflxp + zzflxn) * z2dt / ztmp
346                 zdep2 = bathy(ji,jj) + sshn_e(ji,jj) - rn_wdmin1   ! this one can be moved out of the loop
347                 zdep2 = zdep2 - z2dt * zssh_frc(ji,jj)
348         
349                 IF(zdep1 > zdep2) THEN
350                   zflag = 1
351                   !wdmask(ji, jj) = 1
352                   zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt )
353                   zcoef = max(zcoef, 0._wp)
354                   IF(jk1 > nn_wdit) zcoef = 0._wp
355                   IF(zflxu1(ji,  jj) > 0._wp) zwdlmtu(ji  ,jj) = zcoef
356                   IF(zflxu1(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = zcoef
357                   IF(zflxv1(ji,  jj) > 0._wp) zwdlmtv(ji  ,jj) = zcoef
358                   IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji-1,jj) = zcoef
359                 END IF
360              END DO ! ji loop
361           END DO  ! jj loop
362
363           CALL lbc_lnk( zwdlmtu, 'U', 1. )
364           CALL lbc_lnk( zwdlmtv, 'V', 1. )
365
366           IF(lk_mpp) CALL mpp_max(zflag)   !max over the global domain
367
368           IF(zflag == 0) EXIT
369         
370           zflag = 0     ! flag indicating if any further iteration is needed?
371        END DO  ! jk1 loop
372       
373        zflxu(:,:) = zflxu(:,:) * zwdlmtu(:, :) 
374        zflxv(:,:) = zflxv(:,:) * zwdlmtv(:, :) 
375
376        CALL lbc_lnk( zflxu, 'U', -1. )
377        CALL lbc_lnk( zflxv, 'V', -1. )
378       
379        IF(zflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt_bt!!!'
380       
381        !IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )          ! runoffs (update hdivn field)
382        !IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (update hdivn field)
383        !
384        !
385        CALL wrk_dealloc( jpi, jpj, zflxp, zflxn, zflxu1, zflxv1 )
386        CALL wrk_dealloc( jpi, jpj, zwdlmtu, zwdlmtv)
387      !
388      END IF
389
390      IF( nn_timing == 1 )  CALL timing_stop('wad_lmt')
391   END SUBROUTINE wad_lmt_bt
392END MODULE wet_dry
Note: See TracBrowser for help on using the repository browser.