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/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2016/dev_r6711_SIMPLIF_6_aerobulk/NEMOGCM/NEMO/OPA_SRC/DYN/wet_dry.F90 @ 7163

Last change on this file since 7163 was 7163, checked in by gm, 7 years ago

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