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.
wetdry.F90 in branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/wetdry.F90 @ 5870

Last change on this file since 5870 was 5870, checked in by acc, 8 years ago

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