source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/OPA_SRC/DYN/wet_dry.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 3 years ago

All wrk_alloc removed

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