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/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/wet_dry.F90 @ 9019

Last change on this file since 9019 was 9019, checked in by timgraham, 6 years ago

Merge of dev_CNRS_2017 into branch
File size: 17.6 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_init      : initialisation of wetting and drying
14   !!   wad_lmt       : horizontal flux limiter and limited velocity when wetting and drying happens
15   !!   wad_lmt_bt    : same as wad_lmt for the barotropic stepping (dynspg_ts)
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   !
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 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(:,:) ::   wdmask   !: u- and v- limiter
35   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) ::   ht_wd    !: wetting and drying t-pt depths
36   !                                                           !  (can include negative depths)
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 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"
50   !!----------------------------------------------------------------------
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      INTEGER  ::   ios, ierr   ! Local integer
62      !!
63      NAMELIST/namwad/ ln_wd, rn_wdmin1, rn_wdmin2, rn_wdld, nn_wdit
64      !!----------------------------------------------------------------------
65      !
66      REWIND( numnam_ref )              ! Namelist namwad in reference namelist : 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 : 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        !!gm DOCTOR names: should start with p !
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      IF( ln_timing )   CALL timing_start('wad_lmt')
121      !
122      !IF(lwp) WRITE(numout,*)
123      !IF(lwp) WRITE(numout,*) 'wad_lmt : wetting/drying limiters and velocity limiting'
124      !
125      jflag  = 0
126      zdepwd = 50._wp   !maximum depth on which that W/D could possibly happen
127      !
128      zflxp(:,:)   = 0._wp
129      zflxn(:,:)   = 0._wp
130      zflxu(:,:)   = 0._wp
131      zflxv(:,:)   = 0._wp
132      !
133      zwdlmtu(:,:) = 1._wp
134      zwdlmtv(:,:) = 1._wp
135      !
136      ! Horizontal Flux in u and v direction
137      DO jk = 1, jpkm1 
138         DO jj = 1, jpj
139            DO ji = 1, jpi
140               zflxu(ji,jj) = zflxu(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)
141               zflxv(ji,jj) = zflxv(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
142            END DO 
143         END DO 
144      END DO
145      !
146      zflxu(:,:) = zflxu(:,:) * e2u(:,:)
147      zflxv(:,:) = zflxv(:,:) * e1v(:,:)
148      !
149      wdmask(:,:) = 1
150      DO jj = 2, jpj
151         DO ji = 2, jpi 
152            !
153            IF( tmask(ji, jj, 1) < 0.5_wp )   CYCLE   ! we don't care about land cells
154            IF( ht_wd(ji,jj)     > zdepwd )   CYCLE   ! and cells which are unlikely to dry
155            !
156            zflxp(ji,jj) = MAX( zflxu(ji,jj) , 0._wp ) - MIN( zflxu(ji-1,jj) , 0._wp )   &
157               &         + MAX( zflxv(ji,jj) , 0._wp ) - MIN( zflxv(ji,jj-1) , 0._wp ) 
158            zflxn(ji,jj) = MIN( zflxu(ji,jj) , 0._wp ) - MAX( zflxu(ji-1,jj) , 0._wp )   &
159               &         + MIN( zflxv(ji,jj) , 0._wp ) - MAX( zflxv(ji,jj-1) , 0._wp ) 
160            !
161            zdep2 = ht_wd(ji,jj) + sshb1(ji,jj) - rn_wdmin1
162            IF( zdep2 .le. 0._wp) THEN  !add more safty, but not necessary
163               sshb1(ji,jj) = rn_wdmin1 - ht_wd(ji,jj)
164               IF( zflxu(ji,  jj) > 0._wp )   zwdlmtu(ji  ,jj) = 0._wp
165               IF( zflxu(ji-1,jj) < 0._wp )   zwdlmtu(ji-1,jj) = 0._wp
166               IF( zflxv(ji,  jj) > 0._wp )   zwdlmtv(ji  ,jj) = 0._wp
167               IF( zflxv(ji,jj-1) < 0._wp )   zwdlmtv(ji,jj-1) = 0._wp 
168               wdmask(ji,jj) = 0._wp
169            ENDIF
170         END DO
171      END DO
172      !!
173      !! start limiter iterations
174      DO jk1 = 1, nn_wdit + 1
175         !
176         zflxu1(:,:) = zflxu(:,:) * zwdlmtu(:,:)
177         zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:)
178         jflag = 0     ! flag indicating if any further iterations are needed
179         !
180         DO jj = 2, jpj
181            DO ji = 2, jpi 
182               !
183               IF( tmask(ji,jj,1) < 0.5_wp )   CYCLE
184               IF( ht_wd(ji,jj)   > zdepwd )   CYCLE
185               !
186               ztmp = e1e2t(ji,jj)
187               !
188               zzflxp = MAX( zflxu1(ji,jj) , 0._wp ) - MIN( zflxu1(ji-1,jj) , 0._wp )   &
189                  &   + MAX( zflxv1(ji,jj) , 0._wp ) - MIN( zflxv1(ji,jj-1) , 0._wp ) 
190               zzflxn = MIN( zflxu1(ji,jj) , 0._wp ) - MAX( zflxu1(ji-1,jj) , 0._wp )   &
191                  &   + MIN( zflxv1(ji,jj) , 0._wp ) - MAX( zflxv1(ji,jj-1) , 0._wp ) 
192               !
193               zdep1 = (zzflxp + zzflxn) * z2dt / ztmp
194               zdep2 = ht_wd(ji,jj) + sshb1(ji,jj) - rn_wdmin1 - z2dt * sshemp(ji,jj)
195               !
196               IF( zdep1 > zdep2 ) THEN
197                  wdmask(ji, jj) = 0
198                  zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt )
199                  !zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zzflxp * z2dt )
200                  ! flag if the limiter has been used but stop flagging if the only
201                  ! changes have zeroed the coefficient since further iterations will
202                  ! not change anything
203                  IF( zcoef > 0._wp ) THEN   ;   jflag = 1 
204                  ELSE                       ;   zcoef = 0._wp
205                  ENDIF
206                  IF( jk1 > nn_wdit )   zcoef = 0._wp
207                  IF( zflxu1(ji,  jj) > 0._wp )   zwdlmtu(ji  ,jj) = zcoef
208                  IF( zflxu1(ji-1,jj) < 0._wp )   zwdlmtu(ji-1,jj) = zcoef
209                  IF( zflxv1(ji,  jj) > 0._wp )   zwdlmtv(ji  ,jj) = zcoef
210                  IF( zflxv1(ji,jj-1) < 0._wp )   zwdlmtv(ji,jj-1) = zcoef
211               ENDIF
212            END DO ! ji loop
213         END DO  ! jj loop
214         !
215         CALL lbc_lnk( zwdlmtu, 'U', 1. )
216         CALL lbc_lnk( zwdlmtv, 'V', 1. )
217         !
218         IF(lk_mpp)   CALL mpp_max(jflag)   !max over the global domain
219         !
220         IF(jflag == 0)   EXIT
221         !
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!!gm ==> Andrew  : the lbclnk below is useless since above lbclnk is applied on zwdlmtu/v
230!!                                             and un, vn always with lbclnk
231      CALL lbc_lnk( un, 'U', -1. )
232      CALL lbc_lnk( vn, 'V', -1. )
233!!gm end
234      !
235      un_b(:,:) = un_b(:,:) * zwdlmtu(:,:)
236      vn_b(:,:) = vn_b(:,:) * zwdlmtv(:,:)
237!!gm ==> Andrew   : probably same as above
238      CALL lbc_lnk( un_b, 'U', -1. )
239      CALL lbc_lnk( vn_b, 'V', -1. )
240!!gm end
241       
242      IF(jflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt!!!'
243       
244      !IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )          ! runoffs (update hdivn field)
245      !IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (update hdivn field)
246      !
247      !
248      IF( ln_timing )   CALL timing_stop('wad_lmt')
249      !
250   END SUBROUTINE wad_lmt
251
252
253   SUBROUTINE wad_lmt_bt( zflxu, zflxv, sshn_e, zssh_frc, rdtbt )
254      !!----------------------------------------------------------------------
255      !!                  ***  ROUTINE wad_lmt  ***
256      !!                   
257      !! ** Purpose :   limiting flux in the barotropic stepping (dynspg_ts)
258      !!
259      !! ** Method  : - Prevent negative depth occurring (Not ready for Agrif)
260      !!
261      !! ** Action  : - calculate flux limiter and W/D flag
262      !!----------------------------------------------------------------------
263      REAL(wp), INTENT(in)    ::   rdtbt    ! ocean time-step index
264      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zflxu,  zflxv, sshn_e, zssh_frc 
265      !
266      INTEGER  ::   ji, jj, jk, jk1     ! dummy loop indices
267      INTEGER  ::   jflag               ! local scalar
268      REAL(wp) ::   z2dt
269      REAL(wp) ::   zcoef, zdep1, zdep2 ! local scalars
270      REAL(wp) ::   zzflxp, zzflxn      ! local scalars
271      REAL(wp) ::   zdepwd              ! local scalar, always wet cell depth
272      REAL(wp) ::   ztmp                ! local scalars
273      REAL(wp), DIMENSION(jpi,jpj) ::   zwdlmtu, zwdlmtv         !: W/D flux limiters
274      REAL(wp), DIMENSION(jpi,jpj) ::   zflxp,  zflxn            ! local 2D workspace
275      REAL(wp), DIMENSION(jpi,jpj) ::   zflxu1, zflxv1           ! local 2D workspace
276      !!----------------------------------------------------------------------
277      !
278      IF( ln_timing )  CALL timing_start('wad_lmt_bt')
279      !       
280      !IF(lwp) WRITE(numout,*)
281      !IF(lwp) WRITE(numout,*) 'wad_lmt_bt : wetting/drying limiters and velocity limiting'
282       
283      jflag  = 0
284      zdepwd = 50._wp   !maximum depth that ocean cells can have W/D processes
285
286      z2dt = rdtbt   
287       
288      zflxp(:,:)   = 0._wp
289      zflxn(:,:)   = 0._wp
290
291      zwdlmtu(:,:) = 1._wp
292      zwdlmtv(:,:) = 1._wp
293       
294      ! Horizontal Flux in u and v direction
295       
296      DO jj = 2, jpj
297         DO ji = 2, jpi 
298            !
299            IF( tmask(ji,jj,1) < 0.5_wp )   CYCLE   ! we don't care about land cells
300            IF( ht_wd(ji,jj)   > zdepwd )   CYCLE   ! and cells which are unlikely to dry
301            !
302            zflxp(ji,jj) = MAX( zflxu(ji,jj) , 0._wp ) - MIN( zflxu(ji-1,jj) , 0._wp )   &
303               &         + MAX( zflxv(ji,jj) , 0._wp ) - MIN( zflxv(ji,jj-1) , 0._wp ) 
304            zflxn(ji,jj) = MIN( zflxu(ji,jj) , 0._wp ) - MAX( zflxu(ji-1,jj) , 0._wp )   &
305               &         + MIN( zflxv(ji,jj) , 0._wp ) - MAX( zflxv(ji,jj-1) , 0._wp ) 
306
307            zdep2 = ht_wd(ji,jj) + sshn_e(ji,jj) - rn_wdmin1
308            IF(zdep2 .le. 0._wp) THEN  !add more safety, but not necessary
309               sshn_e(ji,jj) = rn_wdmin1 - ht_wd(ji,jj)
310               IF( zflxu(ji,  jj) > 0._wp )   zwdlmtu(ji  ,jj) = 0._wp
311               IF( zflxu(ji-1,jj) < 0._wp )   zwdlmtu(ji-1,jj) = 0._wp
312               IF( zflxv(ji,  jj) > 0._wp )   zwdlmtv(ji  ,jj) = 0._wp
313               IF( zflxv(ji,jj-1) < 0._wp )   zwdlmtv(ji,jj-1) = 0._wp 
314            ENDIF
315         END DO
316      END DO
317
318     
319      !! start limiter iterations
320      DO jk1 = 1, nn_wdit + 1
321         !
322         zflxu1(:,:) = zflxu(:,:) * zwdlmtu(:,:)
323         zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:)
324         jflag = 0     ! flag indicating if any further iterations are needed
325         !
326         DO jj = 2, jpj
327            DO ji = 2, jpi 
328               !
329               IF( tmask(ji,jj, 1 ) < 0.5_wp  )   CYCLE
330               IF( ht_wd(ji,jj)      > zdepwd )   CYCLE
331               !
332               ztmp = e1e2t(ji,jj)
333               !
334               zzflxp = MAX( zflxu1(ji,jj) , 0._wp ) - MIN( zflxu1(ji-1,jj) , 0._wp )   &
335                  &   + MAX( zflxv1(ji,jj) , 0._wp ) - MIN( zflxv1(ji,jj-1) , 0._wp ) 
336               zzflxn = MIN( zflxu1(ji,jj) , 0._wp ) - MAX( zflxu1(ji-1,jj) , 0._wp )   &
337                  &   + MIN( zflxv1(ji,jj) , 0._wp ) - MAX( zflxv1(ji,jj-1) , 0._wp ) 
338         
339               zdep1 = (zzflxp + zzflxn) * z2dt / ztmp
340               zdep2 = ht_wd(ji,jj) + sshn_e(ji,jj) - rn_wdmin1 - z2dt * zssh_frc(ji,jj)
341         
342               IF(zdep1 > zdep2) THEN
343                  zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt )
344                  !zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zzflxp * z2dt )
345                  ! flag if the limiter has been used but stop flagging if the only
346                  ! changes have zeroed the coefficient since further iterations will
347                  ! not change anything
348                  IF( zcoef > 0._wp ) THEN
349                      jflag = 1 
350                  ELSE
351                      zcoef = 0._wp
352                  ENDIF
353                  IF( jk1 > nn_wdit )   zcoef = 0._wp
354                  IF( zflxu1(ji,  jj) > 0._wp )   zwdlmtu(ji  ,jj) = zcoef
355                  IF( zflxu1(ji-1,jj) < 0._wp )   zwdlmtu(ji-1,jj) = zcoef
356                  IF( zflxv1(ji,  jj) > 0._wp )   zwdlmtv(ji  ,jj) = zcoef
357                  IF( zflxv1(ji,jj-1) < 0._wp )   zwdlmtv(ji,jj-1) = zcoef
358               ENDIF
359            END DO ! ji loop
360         END DO  ! jj loop
361         !
362         CALL lbc_lnk( zwdlmtu, 'U', 1. )
363         CALL lbc_lnk( zwdlmtv, 'V', 1. )
364         !
365         IF(lk_mpp)   CALL mpp_max(jflag)   !max over the global domain
366         !
367         IF( jflag == 0 )   EXIT
368         !   
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(jflag == 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      IF( ln_timing )  CALL timing_stop('wad_lmt_bt')
383      !
384   END SUBROUTINE wad_lmt_bt
385
386   !!==============================================================================
387END MODULE wet_dry
Note: See TracBrowser for help on using the repository browser.