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/UKMO/ROMS_WAD_7832/NEMOGCM/CONFIG/CS15mini/MY_SRC – NEMO

source: branches/UKMO/ROMS_WAD_7832/NEMOGCM/CONFIG/CS15mini/MY_SRC/wet_dry.F90 @ 8544

Last change on this file since 8544 was 8544, checked in by deazer, 7 years ago

Added reference level (above all potential wet points) to avoid negative depth bathymetry as a work around.
Require reference level to be emedded into the domain configuration file
uses ht_0 instead of ht_wd. This si still in the code but should be removed in time.
File size: 20.7 KB
Line 
1MODULE wet_dry
2
3   !! includes updates to namelist namwad for diagnostic outputs of ROMS wetting and drying
4
5   !!==============================================================================
6   !!                       ***  MODULE  wet_dry  ***
7   !! Wetting and drying includes initialisation routine and routines to
8   !! compute and apply flux limiters and preserve water depth positivity
9   !! only effects if wetting/drying is on (ln_wd == .true.)
10   !!==============================================================================
11   !! History :  3.6  ! 2014-09  ((H.Liu)  Original code
12   !!                 ! will add the runoff and periodic BC case later
13   !!----------------------------------------------------------------------
14
15   !!----------------------------------------------------------------------
16   !!   wad_lmt    : Compute the horizontal flux limiter and the limited velocity
17   !!                when wetting and drying happens
18   !!----------------------------------------------------------------------
19   USE oce             ! ocean dynamics and tracers
20   USE dom_oce         ! ocean space and time domain
21   USE sbc_oce, ONLY : ln_rnf   ! surface boundary condition: ocean
22   USE sbcrnf          ! river runoff
23   USE in_out_manager  ! I/O manager
24   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
25   USE lib_mpp         ! MPP library
26   USE wrk_nemo        ! Memory Allocation
27   USE timing          ! Timing
28
29   IMPLICIT NONE
30   PRIVATE
31
32   !!----------------------------------------------------------------------
33   !! critical depths,filters, limiters,and masks for  Wetting and Drying
34   !! ---------------------------------------------------------------------
35
36   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) ::   wdmask         !: u- and v- limiter
37   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) ::   ht_wd          !: wetting and drying t-pt depths
38                                                                     !  (can include negative depths)
39   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) ::   wdramp, wdrampu, wdrampv ! for hpg limiting
40
41   LOGICAL,  PUBLIC  ::   ln_wd       !: Wetting/drying activation switch (T:on,F:off)
42   REAL(wp), PUBLIC  ::   rn_wdmin0   !: depth at which wetting/drying starts
43   LOGICAL,  PUBLIC  ::   ln_rwd      !: ROMS Wetting/drying activation switch (T:on,F:off)
44   REAL(wp), PUBLIC  ::   rn_wdmin1   !: minimum water depth on dried cells
45   REAL(wp), PUBLIC  ::   rn_wdmin2   !: tolerance of minimum water depth on dried cells
46   REAL(wp), PUBLIC  ::   rn_wdld     !: land elevation below which wetting/drying
47                                           !: will be considered
48   INTEGER , PUBLIC  ::   nn_wdit     !: maximum number of iteration for W/D limiter
49   LOGICAL,  PUBLIC  ::   ln_rwd_bc   !: ROMS scheme: True implies 3D velocities are set to the barotropic values at points
50                                      !: where the flow is from wet points on less than half the barotropic sub-steps 
51   LOGICAL,  PUBLIC  ::   ln_rwd_rmp  !: use a ramp for the rwd flux limiter between 2 rn_wdmin1 and rn_wdmin1 (rather than a cut-off at rn_wdmin1)     
52   REAL(wp), PUBLIC  ::   rn_ssh_ref  !: height of z=0 with respect to the geoid;
53   REAL(wp), PUBLIC  ::   rn_ht_0     !: the height at which ht_0 = 0   
54
55   LOGICAL,  PUBLIC  ::   ln_wd_diag  ! True implies wad diagnostic at chosen point are printed out
56   INTEGER , PUBLIC  ::   jn_wd_i, jn_wd_j, jn_wd_k   ! indices at which diagnostic outputs are generated
57
58   PUBLIC   wad_init                  ! initialisation routine called by step.F90
59   PUBLIC   wad_lmt                   ! routine called by sshwzv.F90
60   PUBLIC   wad_lmt_bt                ! routine called by dynspg_ts.F90
61
62   !! * Substitutions
63#  include "vectopt_loop_substitute.h90"
64CONTAINS
65
66   SUBROUTINE wad_init
67      !!----------------------------------------------------------------------
68      !!                     ***  ROUTINE wad_init  ***
69      !!                   
70      !! ** Purpose :   read wetting and drying namelist and print the variables.
71      !!
72      !! ** input   : - namwad namelist
73      !!----------------------------------------------------------------------
74      NAMELIST/namwad/ ln_wd, ln_rwd, rn_wdmin0, ln_rwd, rn_wdmin1, rn_wdmin2, rn_wdld, nn_wdit, ln_wd_diag, ln_rwd_bc, &
75                     &  rn_ht_0, jn_wd_i, jn_wd_j, jn_wd_k,ln_rwd_rmp
76
77      INTEGER  ::   ios                 ! Local integer output status for namelist read
78      INTEGER  ::   ierr                ! Local integer status array allocation
79      !!----------------------------------------------------------------------
80
81      REWIND( numnam_ref )              ! Namelist namwad in reference namelist
82                                        ! : Parameters for Wetting/Drying
83      READ  ( numnam_ref, namwad, IOSTAT = ios, ERR = 905)
84905   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namwad in reference namelist', .TRUE.) 
85      REWIND( numnam_cfg )              ! Namelist namwad in configuration namelist
86                                        ! : Parameters for Wetting/Drying
87      READ  ( numnam_cfg, namwad, IOSTAT = ios, ERR = 906)
88906   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namwad in configuration namelist', .TRUE. )
89      IF(lwm) WRITE ( numond, namwad )
90
91      IF(lwp) THEN                  ! control print
92         WRITE(numout,*)
93         WRITE(numout,*) 'wad_init : Wetting and drying initialization through namelist read'
94         WRITE(numout,*) '~~~~~~~~'
95         WRITE(numout,*) '   Namelist namwad'
96         WRITE(numout,*) '      Logical for NOC  wd scheme         ln_wd      = ', ln_wd
97         WRITE(numout,*) '      Logical for ROMS wd scheme         ln_rwd     = ', ln_rwd
98         WRITE(numout,*) '      Depth at which wet/drying starts rn_wdmin0    = ', rn_wdmin0
99         WRITE(numout,*) '      Minimum wet depth on dried cells rn_wdmin1    = ', rn_wdmin1
100         WRITE(numout,*) '      Tolerance of min wet depth     rn_wdmin2      = ', rn_wdmin2
101         WRITE(numout,*) '      land elevation threshold         rn_wdld      = ', rn_wdld
102         WRITE(numout,*) '      Max iteration for W/D limiter    nn_wdit      = ', nn_wdit
103         WRITE(numout,*) '      Logical for WAD diagnostics      ln_wd_diag   = ', ln_wd_diag
104         WRITE(numout,*) '      T => baroclinic u,v = 0 at dry pts: ln_rwd_bc = ', ln_rwd_bc     
105         WRITE(numout,*) '      use a ramp for rwd limiter:      ln_rwd_rmp   = ', ln_rwd_rmp
106         WRITE(numout,*) '      the height (z) at which ht_0 = 0:rn_ht_0      = ', rn_ht_0 
107         WRITE(numout,*) '      i-index for diagnostic point     jn_wd_i      = ', jn_wd_i
108         WRITE(numout,*) '      j-index for diagnostic point     jn_wd_j      = ', jn_wd_j
109         WRITE(numout,*) '      k-index for diagnostic point     jn_wd_k      = ', jn_wd_k
110      ENDIF
111      !
112      IF(ln_wd .OR. ln_rwd) THEN
113         ALLOCATE( wdmask(jpi,jpj), ht_wd(jpi,jpj),  STAT=ierr )
114         ALLOCATE( wdramp(jpi,jpj), wdrampu(jpi,jpj), wdrampv(jpi,jpj), STAT=ierr ) 
115         IF( ierr /= 0 ) CALL ctl_stop('STOP', 'wad_init : Array allocation error')
116      ENDIF
117      !
118   END SUBROUTINE wad_init
119
120
121   SUBROUTINE wad_lmt( sshb1, sshemp, z2dt )
122      !!----------------------------------------------------------------------
123      !!                  ***  ROUTINE wad_lmt  ***
124      !!                   
125      !! ** Purpose :   generate flux limiters for wetting/drying
126      !!
127      !! ** Method  : - Prevent negative depth occurring (Not ready for Agrif)
128      !!
129      !! ** Action  : - calculate flux limiter and W/D flag
130      !!----------------------------------------------------------------------
131      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   sshb1
132      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   sshemp
133      REAL(wp), INTENT(in) :: z2dt
134      !
135      INTEGER  ::   ji, jj, jk, jk1     ! dummy loop indices
136      INTEGER  ::   jflag               ! local scalar
137      REAL(wp) ::   zcoef, zdep1, zdep2 ! local scalars
138      REAL(wp) ::   zzflxp, zzflxn      ! local scalars
139      REAL(wp) ::   zdepwd              ! local scalar, always wet cell depth
140      REAL(wp) ::   ztmp                ! local scalars
141      REAL(wp), POINTER,  DIMENSION(:,:) ::   zwdlmtu, zwdlmtv         !: W/D flux limiters
142      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxp,  zflxn            ! local 2D workspace
143      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxu,  zflxv            ! local 2D workspace
144      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxu1, zflxv1           ! local 2D workspace
145      !!----------------------------------------------------------------------
146      !
147
148      IF( nn_timing == 1 )  CALL timing_start('wad_lmt')
149
150      IF(ln_wd) THEN
151
152        CALL wrk_alloc( jpi, jpj, zflxp, zflxn, zflxu, zflxv, zflxu1, zflxv1 )
153        CALL wrk_alloc( jpi, jpj, zwdlmtu, zwdlmtv)
154        !
155       
156        !IF(lwp) WRITE(numout,*)
157        !IF(lwp) WRITE(numout,*) 'wad_lmt : wetting/drying limiters and velocity limiting'
158       
159        jflag  = 0
160        zdepwd = 50._wp   !maximum depth on which that W/D could possibly happen
161
162       
163        zflxp(:,:)   = 0._wp
164        zflxn(:,:)   = 0._wp
165        zflxu(:,:)   = 0._wp
166        zflxv(:,:)   = 0._wp
167
168        zwdlmtu(:,:)  = 1._wp
169        zwdlmtv(:,:)  = 1._wp
170       
171        ! Horizontal Flux in u and v direction
172        DO jk = 1, jpkm1 
173           DO jj = 1, jpj
174              DO ji = 1, jpi
175                 zflxu(ji,jj) = zflxu(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)
176                 zflxv(ji,jj) = zflxv(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
177              END DO 
178           END DO 
179        END DO
180       
181        zflxu(:,:) = zflxu(:,:) * e2u(:,:)
182        zflxv(:,:) = zflxv(:,:) * e1v(:,:)
183       
184        wdmask(:,:) = 1
185        DO jj = 2, jpj
186           DO ji = 2, jpi 
187
188             IF( tmask(ji, jj, 1) < 0.5_wp ) CYCLE   ! we don't care about land cells
189             IF( ht_0(ji,jj)-rn_ssh_ref > zdepwd )      CYCLE   ! and cells which are unlikely to dry
190
191              zflxp(ji,jj) = max(zflxu(ji,jj), 0._wp) - min(zflxu(ji-1,jj),   0._wp) + &
192                           & max(zflxv(ji,jj), 0._wp) - min(zflxv(ji,  jj-1), 0._wp) 
193              zflxn(ji,jj) = min(zflxu(ji,jj), 0._wp) - max(zflxu(ji-1,jj),   0._wp) + &
194                           & min(zflxv(ji,jj), 0._wp) - max(zflxv(ji,  jj-1), 0._wp) 
195
196              zdep2 = ht_0(ji,jj) + sshb1(ji,jj) - rn_wdmin1
197              IF(zdep2 .le. 0._wp) THEN  !add more safty, but not necessary
198                sshb1(ji,jj) = rn_wdmin1 - ht_0(ji,jj)
199                IF(zflxu(ji,  jj) > 0._wp) zwdlmtu(ji  ,jj) = 0._wp
200                IF(zflxu(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = 0._wp
201                IF(zflxv(ji,  jj) > 0._wp) zwdlmtv(ji  ,jj) = 0._wp
202                IF(zflxv(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = 0._wp 
203                wdmask(ji,jj) = 0._wp
204              END IF
205           ENDDO
206        END DO
207
208
209! slwa
210! HPG limiter from jholt
211      wdramp(:,:) = min((ht_0(:,:) + sshb1(:,:) - rn_wdmin1)/(rn_wdmin0 - rn_wdmin1),1.0_wp)
212!      write(6,*)'wdramp ',wdramp(10,10),wdramp(10,11)
213!jth assume don't need a lbc_lnk here
214       DO jj = 1, jpjm1
215         DO ji = 1, jpim1
216           wdrampu(ji,jj) = min(wdramp(ji,jj),wdramp(ji+1,jj))
217           wdrampv(ji,jj) = min(wdramp(ji,jj),wdramp(ji,jj+1))
218         ENDDO
219       ENDDO
220        !wdrampu(:,:)=1.0_wp
221        !wdrampv(:,:)=1.0_wp
222! end HPG limiter
223
224
225     
226        !! start limiter iterations
227        DO jk1 = 1, nn_wdit + 1
228       
229         
230           zflxu1(:,:) = zflxu(:,:) * zwdlmtu(:,:)
231           zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:)
232           jflag = 0     ! flag indicating if any further iterations are needed
233         
234           DO jj = 2, jpj
235              DO ji = 2, jpi 
236       
237                 IF( tmask(ji, jj, 1) < 0.5_wp ) CYCLE
238                 IF( ht_0(ji,jj) > zdepwd )      CYCLE
239       
240                 ztmp = e1e2t(ji,jj)
241
242                 zzflxp = max(zflxu1(ji,jj), 0._wp) - min(zflxu1(ji-1,jj),   0._wp) + &
243                        & max(zflxv1(ji,jj), 0._wp) - min(zflxv1(ji,  jj-1), 0._wp) 
244                 zzflxn = min(zflxu1(ji,jj), 0._wp) - max(zflxu1(ji-1,jj),   0._wp) + &
245                        & min(zflxv1(ji,jj), 0._wp) - max(zflxv1(ji,  jj-1), 0._wp) 
246         
247                 zdep1 = (zzflxp + zzflxn) * z2dt / ztmp
248                 zdep2 = ht_0(ji,jj) + sshb1(ji,jj) - rn_wdmin1 - z2dt * sshemp(ji,jj)
249         
250                 IF( zdep1 > zdep2 ) THEN
251                   wdmask(ji, jj) = 0
252                   zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt )
253                   !zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zzflxp * z2dt )
254                   ! flag if the limiter has been used but stop flagging if the only
255                   ! changes have zeroed the coefficient since further iterations will
256                   ! not change anything
257                   IF( zcoef > 0._wp ) THEN
258                      jflag = 1 
259                   ELSE
260                      zcoef = 0._wp
261                   ENDIF
262                   IF(jk1 > nn_wdit) zcoef = 0._wp
263                   IF(zflxu1(ji,  jj) > 0._wp) zwdlmtu(ji  ,jj) = zcoef
264                   IF(zflxu1(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = zcoef
265                   IF(zflxv1(ji,  jj) > 0._wp) zwdlmtv(ji  ,jj) = zcoef
266                   IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = zcoef
267                 END IF
268              END DO ! ji loop
269           END DO  ! jj loop
270
271           CALL lbc_lnk( zwdlmtu, 'U', 1. )
272           CALL lbc_lnk( zwdlmtv, 'V', 1. )
273
274           IF(lk_mpp) CALL mpp_max(jflag)   !max over the global domain
275
276           IF(jflag == 0) EXIT
277         
278        END DO  ! jk1 loop
279       
280        DO jk = 1, jpkm1
281          un(:,:,jk) = un(:,:,jk) * zwdlmtu(:, :) 
282          vn(:,:,jk) = vn(:,:,jk) * zwdlmtv(:, :) 
283        END DO
284
285        CALL lbc_lnk( un, 'U', -1. )
286        CALL lbc_lnk( vn, 'V', -1. )
287      !
288        un_b(:,:) = un_b(:,:) * zwdlmtu(:, :)
289        vn_b(:,:) = vn_b(:,:) * zwdlmtv(:, :)
290        CALL lbc_lnk( un_b, 'U', -1. )
291        CALL lbc_lnk( vn_b, 'V', -1. )
292       
293        IF(jflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt!!!'
294       
295        !IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )          ! runoffs (update hdivn field)
296        !IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (update hdivn field)
297        !
298        !
299        CALL wrk_dealloc( jpi, jpj, zflxp, zflxn, zflxu, zflxv, zflxu1, zflxv1 )
300        CALL wrk_dealloc( jpi, jpj, zwdlmtu, zwdlmtv)
301        !
302      ENDIF
303      !
304      IF( nn_timing == 1 )  CALL timing_stop('wad_lmt')
305      !
306   END SUBROUTINE wad_lmt
307
308
309   SUBROUTINE wad_lmt_bt( zflxu, zflxv, sshn_e, zssh_frc, rdtbt )
310      !!----------------------------------------------------------------------
311      !!                  ***  ROUTINE wad_lmt  ***
312      !!                   
313      !! ** Purpose :   limiting flux in the barotropic stepping (dynspg_ts)
314      !!
315      !! ** Method  : - Prevent negative depth occurring (Not ready for Agrif)
316      !!
317      !! ** Action  : - calculate flux limiter and W/D flag
318      !!----------------------------------------------------------------------
319      REAL(wp), INTENT(in)    ::   rdtbt    ! ocean time-step index
320      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zflxu,  zflxv, sshn_e, zssh_frc 
321      !
322      INTEGER  ::   ji, jj, jk, jk1     ! dummy loop indices
323      INTEGER  ::   jflag               ! local scalar
324      REAL(wp) ::   z2dt
325      REAL(wp) ::   zcoef, zdep1, zdep2 ! local scalars
326      REAL(wp) ::   zzflxp, zzflxn      ! local scalars
327      REAL(wp) ::   zdepwd              ! local scalar, always wet cell depth
328      REAL(wp) ::   ztmp                ! local scalars
329      REAL(wp), POINTER,  DIMENSION(:,:) ::   zwdlmtu, zwdlmtv         !: W/D flux limiters
330      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxp,  zflxn            ! local 2D workspace
331      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxu1, zflxv1           ! local 2D workspace
332      !!----------------------------------------------------------------------
333      !
334      IF( nn_timing == 1 )  CALL timing_start('wad_lmt_bt')
335
336      IF(ln_wd) THEN
337
338        CALL wrk_alloc( jpi, jpj, zflxp, zflxn, zflxu1, zflxv1 )
339        CALL wrk_alloc( jpi, jpj, zwdlmtu, zwdlmtv)
340        !
341       
342        !IF(lwp) WRITE(numout,*)
343        !IF(lwp) WRITE(numout,*) 'wad_lmt_bt : wetting/drying limiters and velocity limiting'
344       
345        jflag  = 0
346        zdepwd = 50._wp   !maximum depth that ocean cells can have W/D processes
347
348        z2dt = rdtbt   
349       
350        zflxp(:,:)   = 0._wp
351        zflxn(:,:)   = 0._wp
352
353        zwdlmtu(:,:)  = 1._wp
354        zwdlmtv(:,:)  = 1._wp
355       
356        ! Horizontal Flux in u and v direction
357       
358        DO jj = 2, jpj
359           DO ji = 2, jpi 
360
361             IF( tmask(ji, jj, 1 ) < 0.5_wp) CYCLE   ! we don't care about land cells
362             IF( ht_0(ji,jj) > zdepwd )      CYCLE   ! and cells which are unlikely to dry
363
364              zflxp(ji,jj) = max(zflxu(ji,jj), 0._wp) - min(zflxu(ji-1,jj),   0._wp) + &
365                           & max(zflxv(ji,jj), 0._wp) - min(zflxv(ji,  jj-1), 0._wp) 
366              zflxn(ji,jj) = min(zflxu(ji,jj), 0._wp) - max(zflxu(ji-1,jj),   0._wp) + &
367                           & min(zflxv(ji,jj), 0._wp) - max(zflxv(ji,  jj-1), 0._wp) 
368
369              zdep2 = ht_0(ji,jj) + sshn_e(ji,jj) - rn_wdmin1
370              IF(zdep2 .le. 0._wp) THEN  !add more safety, but not necessary
371                sshn_e(ji,jj) = rn_wdmin1 - ht_0(ji,jj)
372                IF(zflxu(ji,  jj) > 0._wp) zwdlmtu(ji  ,jj) = 0._wp
373                IF(zflxu(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = 0._wp
374                IF(zflxv(ji,  jj) > 0._wp) zwdlmtv(ji  ,jj) = 0._wp
375                IF(zflxv(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = 0._wp 
376              END IF
377           ENDDO
378        END DO
379
380     
381        !! start limiter iterations
382        DO jk1 = 1, nn_wdit + 1
383       
384         
385           zflxu1(:,:) = zflxu(:,:) * zwdlmtu(:,:)
386           zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:)
387           jflag = 0     ! flag indicating if any further iterations are needed
388         
389           DO jj = 2, jpj
390              DO ji = 2, jpi 
391       
392                 IF( tmask(ji, jj, 1 ) < 0.5_wp) CYCLE
393                 IF( ht_0(ji,jj) > zdepwd )      CYCLE
394       
395                 ztmp = e1e2t(ji,jj)
396
397                 zzflxp = max(zflxu1(ji,jj), 0._wp) - min(zflxu1(ji-1,jj),   0._wp) + &
398                        & max(zflxv1(ji,jj), 0._wp) - min(zflxv1(ji,  jj-1), 0._wp) 
399                 zzflxn = min(zflxu1(ji,jj), 0._wp) - max(zflxu1(ji-1,jj),   0._wp) + &
400                        & min(zflxv1(ji,jj), 0._wp) - max(zflxv1(ji,  jj-1), 0._wp) 
401         
402                 zdep1 = (zzflxp + zzflxn) * z2dt / ztmp
403                 zdep2 = ht_0(ji,jj) + sshn_e(ji,jj) - rn_wdmin1 - z2dt * zssh_frc(ji,jj)
404         
405                 IF(zdep1 > zdep2) THEN
406                   zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt )
407                   !zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zzflxp * z2dt )
408                   ! flag if the limiter has been used but stop flagging if the only
409                   ! changes have zeroed the coefficient since further iterations will
410                   ! not change anything
411                   IF( zcoef > 0._wp ) THEN
412                      jflag = 1 
413                   ELSE
414                      zcoef = 0._wp
415                   ENDIF
416                   IF(jk1 > nn_wdit) zcoef = 0._wp
417                   IF(zflxu1(ji,  jj) > 0._wp) zwdlmtu(ji  ,jj) = zcoef
418                   IF(zflxu1(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = zcoef
419                   IF(zflxv1(ji,  jj) > 0._wp) zwdlmtv(ji  ,jj) = zcoef
420                   IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = zcoef
421                 END IF
422              END DO ! ji loop
423           END DO  ! jj loop
424
425           CALL lbc_lnk( zwdlmtu, 'U', 1. )
426           CALL lbc_lnk( zwdlmtv, 'V', 1. )
427
428           IF(lk_mpp) CALL mpp_max(jflag)   !max over the global domain
429
430           IF(jflag == 0) EXIT
431         
432        END DO  ! jk1 loop
433       
434        zflxu(:,:) = zflxu(:,:) * zwdlmtu(:, :) 
435        zflxv(:,:) = zflxv(:,:) * zwdlmtv(:, :) 
436
437        CALL lbc_lnk( zflxu, 'U', -1. )
438        CALL lbc_lnk( zflxv, 'V', -1. )
439       
440        IF(jflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt_bt!!!'
441       
442        !IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )          ! runoffs (update hdivn field)
443        !IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (update hdivn field)
444        !
445        !
446        CALL wrk_dealloc( jpi, jpj, zflxp, zflxn, zflxu1, zflxv1 )
447        CALL wrk_dealloc( jpi, jpj, zwdlmtu, zwdlmtv)
448        !
449      END IF
450
451      IF( nn_timing == 1 )  CALL timing_stop('wad_lmt')
452   END SUBROUTINE wad_lmt_bt
453
454   !!==============================================================================
455END MODULE wet_dry
Note: See TracBrowser for help on using the repository browser.