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 NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/OCE/DYN – NEMO

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/OCE/DYN/wet_dry.F90 @ 11053

Last change on this file since 11053 was 11053, checked in by davestorkey, 5 years ago

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap : Merge in latest changes from main branch and finish conversion of "h" variables. NB. This version still doesn't work!

  • Property svn:keywords set to Id
File size: 20.0 KB
RevLine 
[6152]1MODULE wet_dry
[9023]2
3   !! includes updates to namelist namwad for diagnostic outputs of ROMS wetting and drying
4
[6152]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
[9023]9   !! only effects if wetting/drying is on (ln_wd_il == .true. or ln_wd_dl==.true. )
[6152]10   !!==============================================================================
[7646]11   !! History :  3.6  ! 2014-09  ((H.Liu)  Original code
[6152]12   !!                 ! will add the runoff and periodic BC case later
13   !!----------------------------------------------------------------------
14
15   !!----------------------------------------------------------------------
[9019]16   !!   wad_init      : initialisation of wetting and drying
17   !!   wad_lmt       : horizontal flux limiter and limited velocity when wetting and drying happens
18   !!   wad_lmt_bt    : same as wad_lmt for the barotropic stepping (dynspg_ts)
[6152]19   !!----------------------------------------------------------------------
[9019]20   USE oce            ! ocean dynamics and tracers
21   USE dom_oce        ! ocean space and time domain
22   USE sbc_oce  , ONLY: ln_rnf   ! surface boundary condition: ocean
23   USE sbcrnf         ! river runoff
24   !
25   USE in_out_manager ! I/O manager
26   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
27   USE lib_mpp        ! MPP library
[9124]28   USE timing         ! timing of the main modules
[6152]29
30   IMPLICIT NONE
31   PRIVATE
32
33   !!----------------------------------------------------------------------
34   !! critical depths,filters, limiters,and masks for  Wetting and Drying
35   !! ---------------------------------------------------------------------
36
[9019]37   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) ::   wdmask   !: u- and v- limiter
38   !                                                           !  (can include negative depths)
[9023]39   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) ::   wdramp, wdrampu, wdrampv !: for hpg limiting
[6152]40
[9023]41   LOGICAL,  PUBLIC  ::   ln_wd_il    !: Wetting/drying il activation switch (T:on,F:off)
42   LOGICAL,  PUBLIC  ::   ln_wd_dl    !: Wetting/drying dl activation switch (T:on,F:off)
43   REAL(wp), PUBLIC  ::   rn_wdmin0   !: depth at which wetting/drying starts
[6152]44   REAL(wp), PUBLIC  ::   rn_wdmin1   !: minimum water depth on dried cells
[9023]45   REAL(wp), PUBLIC  ::   r_rn_wdmin1 !: 1/minimum water depth on dried cells
46   REAL(wp), PUBLIC  ::   rn_wdmin2   !: tolerance of minimum water depth on dried cells
[10499]47   REAL(wp), PUBLIC  ::   rn_wd_sbcdep   !: Depth at which to taper sbc fluxes
48   REAL(wp), PUBLIC  ::   rn_wd_sbcfra   !: Fraction of SBC at taper depth
[9019]49   REAL(wp), PUBLIC  ::   rn_wdld     !: land elevation below which wetting/drying will be considered
[6152]50   INTEGER , PUBLIC  ::   nn_wdit     !: maximum number of iteration for W/D limiter
[9023]51   LOGICAL,  PUBLIC  ::   ln_wd_dl_bc !: DL scheme: True implies 3D velocities are set to the barotropic values at points
52                                      !: where the flow is from wet points on less than half the barotropic sub-steps 
53   LOGICAL,  PUBLIC  ::  ln_wd_dl_rmp !: use a ramp for the dl flux limiter between 2 rn_wdmin1 and rn_wdmin1 (rather than a cut-off at rn_wdmin1)     
54   REAL(wp), PUBLIC  ::   ssh_ref     !: height of z=0 with respect to the geoid;
[6152]55
[9023]56   LOGICAL,  PUBLIC  ::   ll_wd       !: Wetting/drying activation switch if either ln_wd_il or ln_wd_dl
57
[6152]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"
[9019]64   !!----------------------------------------------------------------------
[6152]65CONTAINS
66
67   SUBROUTINE wad_init
68      !!----------------------------------------------------------------------
69      !!                     ***  ROUTINE wad_init  ***
70      !!                   
71      !! ** Purpose :   read wetting and drying namelist and print the variables.
72      !!
73      !! ** input   : - namwad namelist
74      !!----------------------------------------------------------------------
[9124]75      INTEGER  ::   ios, ierr   ! Local integer
[9019]76      !!
[9124]77      NAMELIST/namwad/ ln_wd_il, ln_wd_dl   , rn_wdmin0, rn_wdmin1, rn_wdmin2, rn_wdld,   &
[10499]78         &             nn_wdit , ln_wd_dl_bc, ln_wd_dl_rmp, rn_wd_sbcdep,rn_wd_sbcfra
[6152]79      !!----------------------------------------------------------------------
[9019]80      !
81      REWIND( numnam_ref )              ! Namelist namwad in reference namelist : Parameters for Wetting/Drying
[6152]82      READ  ( numnam_ref, namwad, IOSTAT = ios, ERR = 905)
[9168]83905   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namwad in reference namelist', .TRUE.) 
[9019]84      REWIND( numnam_cfg )              ! Namelist namwad in configuration namelist : Parameters for Wetting/Drying
[6152]85      READ  ( numnam_cfg, namwad, IOSTAT = ios, ERR = 906)
[9168]86906   IF( ios >  0 )   CALL ctl_nam ( ios , 'namwad in configuration namelist', .TRUE. )
[6152]87      IF(lwm) WRITE ( numond, namwad )
[9019]88      !
[10499]89      IF( rn_wd_sbcfra>=1 )   CALL ctl_stop( 'STOP', 'rn_wd_sbcfra >=1 : must be < 1' )
[6152]90      IF(lwp) THEN                  ! control print
91         WRITE(numout,*)
92         WRITE(numout,*) 'wad_init : Wetting and drying initialization through namelist read'
[7646]93         WRITE(numout,*) '~~~~~~~~'
[6152]94         WRITE(numout,*) '   Namelist namwad'
[9023]95         WRITE(numout,*) '      Logical for Iter Lim wd option   ln_wd_il     = ', ln_wd_il
96         WRITE(numout,*) '      Logical for Dir. Lim wd option   ln_wd_dl     = ', ln_wd_dl
97         WRITE(numout,*) '      Depth at which wet/drying starts rn_wdmin0    = ', rn_wdmin0
[6152]98         WRITE(numout,*) '      Minimum wet depth on dried cells rn_wdmin1    = ', rn_wdmin1
[9023]99         WRITE(numout,*) '      Tolerance of min wet depth       rn_wdmin2    = ', rn_wdmin2
[6152]100         WRITE(numout,*) '      land elevation threshold         rn_wdld      = ', rn_wdld
101         WRITE(numout,*) '      Max iteration for W/D limiter    nn_wdit      = ', nn_wdit
[9023]102         WRITE(numout,*) '      T => baroclinic u,v=0 at dry pts: ln_wd_dl_bc = ', ln_wd_dl_bc     
103         WRITE(numout,*) '      use a ramp for rwd limiter:  ln_wd_dl_rwd_rmp = ', ln_wd_dl_rmp
[10499]104         WRITE(numout,*) '      cut off depth sbc for wd   rn_wd_sbcdep       = ', rn_wd_sbcdep
105         WRITE(numout,*) '      fraction to start sbc wgt rn_wd_sbcfra        = ', rn_wd_sbcfra
[7646]106      ENDIF
[9023]107      IF( .NOT. ln_read_cfg ) THEN
[9042]108         IF(lwp) WRITE(numout,*) '      No configuration file so seting ssh_ref to zero  '
[9124]109         ssh_ref=0._wp
[9023]110      ENDIF
111
[9124]112      r_rn_wdmin1 = 1 / rn_wdmin1
[9023]113      ll_wd = .FALSE.
[9124]114      IF( ln_wd_il .OR. ln_wd_dl ) THEN
[9023]115         ll_wd = .TRUE.
116         ALLOCATE( wdmask(jpi,jpj),   STAT=ierr )
117         ALLOCATE( wdramp(jpi,jpj), wdrampu(jpi,jpj), wdrampv(jpi,jpj), STAT=ierr ) 
[6152]118         IF( ierr /= 0 ) CALL ctl_stop('STOP', 'wad_init : Array allocation error')
119      ENDIF
[7646]120      !
[6152]121   END SUBROUTINE wad_init
122
[7646]123
[11053]124   SUBROUTINE wad_lmt( psshb1, psshemp, z2dt, Kmm, puu, pvv )
[6152]125      !!----------------------------------------------------------------------
126      !!                  ***  ROUTINE wad_lmt  ***
127      !!                   
128      !! ** Purpose :   generate flux limiters for wetting/drying
129      !!
130      !! ** Method  : - Prevent negative depth occurring (Not ready for Agrif)
131      !!
132      !! ** Action  : - calculate flux limiter and W/D flag
133      !!----------------------------------------------------------------------
[11053]134      REAL(wp), DIMENSION(:,:)            , INTENT(inout) ::   psshb1
135      REAL(wp), DIMENSION(:,:)            , INTENT(in   ) ::   psshemp
136      REAL(wp)                            , INTENT(in   ) ::   z2dt
137      INTEGER                             , INTENT(in   ) ::   Kmm       ! time level index
138      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::   puu, pvv  ! velocity arrays
[6152]139      !
140      INTEGER  ::   ji, jj, jk, jk1     ! dummy loop indices
[7646]141      INTEGER  ::   jflag               ! local scalar
[6152]142      REAL(wp) ::   zcoef, zdep1, zdep2 ! local scalars
143      REAL(wp) ::   zzflxp, zzflxn      ! local scalars
144      REAL(wp) ::   zdepwd              ! local scalar, always wet cell depth
145      REAL(wp) ::   ztmp                ! local scalars
[9019]146      REAL(wp),  DIMENSION(jpi,jpj) ::   zwdlmtu, zwdlmtv   ! W/D flux limiters
147      REAL(wp),  DIMENSION(jpi,jpj) ::   zflxp  ,  zflxn    ! local 2D workspace
148      REAL(wp),  DIMENSION(jpi,jpj) ::   zflxu  ,  zflxv    ! local 2D workspace
149      REAL(wp),  DIMENSION(jpi,jpj) ::   zflxu1 , zflxv1    ! local 2D workspace
[6152]150      !!----------------------------------------------------------------------
[9124]151      IF( ln_timing )   CALL timing_start('wad_lmt')      !
[6152]152      !
[9023]153      DO jk = 1, jpkm1
[11053]154         puu(:,:,jk,Kmm) = puu(:,:,jk,Kmm) * zwdlmtu(:,:) 
155         pvv(:,:,jk,Kmm) = pvv(:,:,jk,Kmm) * zwdlmtv(:,:) 
[9023]156      END DO
[9019]157      jflag  = 0
[9124]158      zdepwd = 50._wp      ! maximum depth on which that W/D could possibly happen
159      !
[9019]160      zflxp(:,:)   = 0._wp
161      zflxn(:,:)   = 0._wp
162      zflxu(:,:)   = 0._wp
163      zflxv(:,:)   = 0._wp
[9124]164      !
165      zwdlmtu(:,:) = 1._wp
166      zwdlmtv(:,:) = 1._wp
167      !
168      DO jk = 1, jpkm1     ! Horizontal Flux in u and v direction
[11053]169         zflxu(:,:) = zflxu(:,:) + e3u(:,:,jk,Kmm) * puu(:,:,jk,Kmm) * umask(:,:,jk)
170         zflxv(:,:) = zflxv(:,:) + e3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm) * vmask(:,:,jk)
[9019]171      END DO
172      zflxu(:,:) = zflxu(:,:) * e2u(:,:)
173      zflxv(:,:) = zflxv(:,:) * e1v(:,:)
[9124]174      !
175      wdmask(:,:) = 1._wp
[9019]176      DO jj = 2, jpj
177         DO ji = 2, jpi 
[9124]178            !
179            IF( tmask(ji,jj,1)        < 0.5_wp )   CYCLE    ! we don't care about land cells
180            IF( ht_0(ji,jj) - ssh_ref > zdepwd )   CYCLE    ! and cells which are unlikely to dry
181            !
182            zflxp(ji,jj) = MAX( zflxu(ji,jj) , 0._wp ) - MIN( zflxu(ji-1,jj  ) , 0._wp )   &
183               &         + MAX( zflxv(ji,jj) , 0._wp ) - MIN( zflxv(ji,  jj-1) , 0._wp ) 
184            zflxn(ji,jj) = MIN( zflxu(ji,jj) , 0._wp ) - MAX( zflxu(ji-1,jj  ) , 0._wp )   &
185               &         + MIN( zflxv(ji,jj) , 0._wp ) - MAX( zflxv(ji,  jj-1) , 0._wp ) 
186            !
[11053]187            zdep2 = ht_0(ji,jj) + psshb1(ji,jj) - rn_wdmin1
[9124]188            IF( zdep2 <= 0._wp ) THEN     ! add more safty, but not necessary
[11053]189               psshb1(ji,jj) = rn_wdmin1 - ht_0(ji,jj)
[9023]190               IF(zflxu(ji,  jj) > 0._wp) zwdlmtu(ji  ,jj) = 0._wp
191               IF(zflxu(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = 0._wp
192               IF(zflxv(ji,  jj) > 0._wp) zwdlmtv(ji  ,jj) = 0._wp
193               IF(zflxv(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = 0._wp 
[9019]194               wdmask(ji,jj) = 0._wp
[9023]195            END IF
[9124]196         END DO
[9023]197      END DO
[9124]198      !
199      !           ! HPG limiter from jholt
[11053]200      wdramp(:,:) = min((ht_0(:,:) + psshb1(:,:) - rn_wdmin1)/(rn_wdmin0 - rn_wdmin1),1.0_wp)
[9124]201      !jth assume don't need a lbc_lnk here
[9023]202      DO jj = 1, jpjm1
203         DO ji = 1, jpim1
[9124]204            wdrampu(ji,jj) = MIN( wdramp(ji,jj) , wdramp(ji+1,jj) )
205            wdrampv(ji,jj) = MIN( wdramp(ji,jj) , wdramp(ji,jj+1) )
[9019]206         END DO
207      END DO
[9124]208      !           ! end HPG limiter
209      !
210      !
211      DO jk1 = 1, nn_wdit + 1      !==  start limiter iterations  ==!
212         !
[9019]213         zflxu1(:,:) = zflxu(:,:) * zwdlmtu(:,:)
214         zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:)
215         jflag = 0     ! flag indicating if any further iterations are needed
[9124]216         !
[9019]217         DO jj = 2, jpj
218            DO ji = 2, jpi 
[9124]219               IF( tmask(ji, jj, 1) < 0.5_wp )   CYCLE
220               IF( ht_0(ji,jj)      > zdepwd )   CYCLE
221               !
[9019]222               ztmp = e1e2t(ji,jj)
[9124]223               !
224               zzflxp = MAX( zflxu1(ji,jj) , 0._wp ) - MIN( zflxu1(ji-1,jj  ) , 0._wp)   &
225                  &   + MAX( zflxv1(ji,jj) , 0._wp ) - MIN( zflxv1(ji,  jj-1) , 0._wp) 
226               zzflxn = MIN( zflxu1(ji,jj) , 0._wp ) - MAX( zflxu1(ji-1,jj  ) , 0._wp)   &
227                  &   + MIN( zflxv1(ji,jj) , 0._wp ) - MAX( zflxv1(ji,  jj-1) , 0._wp) 
228               !
[9019]229               zdep1 = (zzflxp + zzflxn) * z2dt / ztmp
[11053]230               zdep2 = ht_0(ji,jj) + psshb1(ji,jj) - rn_wdmin1 - z2dt * psshemp(ji,jj)
[9124]231               !
[9019]232               IF( zdep1 > zdep2 ) THEN
[9124]233                  wdmask(ji, jj) = 0._wp
[9019]234                  zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt )
235                  !zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zzflxp * z2dt )
236                  ! flag if the limiter has been used but stop flagging if the only
237                  ! changes have zeroed the coefficient since further iterations will
238                  ! not change anything
[9124]239                  IF( zcoef > 0._wp ) THEN   ;   jflag = 1 
240                  ELSE                       ;   zcoef = 0._wp
[9019]241                  ENDIF
[9124]242                  IF( jk1 > nn_wdit )   zcoef = 0._wp
243                  IF( zflxu1(ji  ,jj  ) > 0._wp )   zwdlmtu(ji  ,jj  ) = zcoef
244                  IF( zflxu1(ji-1,jj  ) < 0._wp )   zwdlmtu(ji-1,jj  ) = zcoef
245                  IF( zflxv1(ji  ,jj  ) > 0._wp )   zwdlmtv(ji  ,jj  ) = zcoef
246                  IF( zflxv1(ji  ,jj-1) < 0._wp )   zwdlmtv(ji  ,jj-1) = zcoef
247               ENDIF
248            END DO
249         END DO
[10425]250         CALL lbc_lnk_multi( 'wet_dry', zwdlmtu, 'U', 1., zwdlmtv, 'V', 1. )
[9124]251         !
[10425]252         CALL mpp_max('wet_dry', jflag)   !max over the global domain
[9124]253         !
254         IF( jflag == 0 )   EXIT
255         !
[9019]256      END DO  ! jk1 loop
[9124]257      !
[9019]258      DO jk = 1, jpkm1
[11053]259         puu(:,:,jk,Kmm) = puu(:,:,jk,Kmm) * zwdlmtu(:,:) 
260         pvv(:,:,jk,Kmm) = pvv(:,:,jk,Kmm) * zwdlmtv(:,:) 
[9019]261      END DO
[11053]262      uu_b(:,:,Kmm) = uu_b(:,:,Kmm) * zwdlmtu(:, :)
263      vv_b(:,:,Kmm) = vv_b(:,:,Kmm) * zwdlmtv(:, :)
[9124]264      !
265!!gm TO BE SUPPRESSED ?  these lbc_lnk are useless since zwdlmtu and zwdlmtv are defined everywhere !
[11053]266      CALL lbc_lnk_multi( 'wet_dry', puu(:,:,:,Kmm)  , 'U', -1., pvv(:,:,:,Kmm)  , 'V', -1. )
267      CALL lbc_lnk_multi( 'wet_dry', uu_b(:,:,Kmm), 'U', -1., vv_b(:,:,Kmm), 'V', -1. )
[9124]268!!gm
[6152]269      !
[9124]270      IF(jflag == 1 .AND. lwp)   WRITE(numout,*) 'Need more iterations in wad_lmt!!!'
[7646]271      !
[11053]272      !IF( ln_rnf      )   CALL sbc_rnf_div( hdiv )          ! runoffs (update hdiv field)
[9023]273      !
[9124]274      IF( ln_timing )   CALL timing_stop('wad_lmt')      !
[9023]275      !
[6152]276   END SUBROUTINE wad_lmt
277
[7646]278
[6152]279   SUBROUTINE wad_lmt_bt( zflxu, zflxv, sshn_e, zssh_frc, rdtbt )
280      !!----------------------------------------------------------------------
281      !!                  ***  ROUTINE wad_lmt  ***
282      !!                   
283      !! ** Purpose :   limiting flux in the barotropic stepping (dynspg_ts)
284      !!
285      !! ** Method  : - Prevent negative depth occurring (Not ready for Agrif)
286      !!
287      !! ** Action  : - calculate flux limiter and W/D flag
288      !!----------------------------------------------------------------------
[9124]289      REAL(wp)                , INTENT(in   ) ::   rdtbt    ! ocean time-step index
[6152]290      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zflxu,  zflxv, sshn_e, zssh_frc 
291      !
292      INTEGER  ::   ji, jj, jk, jk1     ! dummy loop indices
[9124]293      INTEGER  ::   jflag               ! local integer
[6152]294      REAL(wp) ::   z2dt
295      REAL(wp) ::   zcoef, zdep1, zdep2 ! local scalars
296      REAL(wp) ::   zzflxp, zzflxn      ! local scalars
297      REAL(wp) ::   zdepwd              ! local scalar, always wet cell depth
298      REAL(wp) ::   ztmp                ! local scalars
[9019]299      REAL(wp), DIMENSION(jpi,jpj) ::   zwdlmtu, zwdlmtv         !: W/D flux limiters
300      REAL(wp), DIMENSION(jpi,jpj) ::   zflxp,  zflxn            ! local 2D workspace
301      REAL(wp), DIMENSION(jpi,jpj) ::   zflxu1, zflxv1           ! local 2D workspace
[6152]302      !!----------------------------------------------------------------------
[9124]303      IF( ln_timing )   CALL timing_start('wad_lmt_bt')      !
[6152]304      !
[9124]305      jflag  = 0
306      zdepwd = 50._wp   ! maximum depth that ocean cells can have W/D processes
[9023]307      !
[9019]308      z2dt = rdtbt   
[9124]309      !
[9019]310      zflxp(:,:)   = 0._wp
311      zflxn(:,:)   = 0._wp
312      zwdlmtu(:,:) = 1._wp
313      zwdlmtv(:,:) = 1._wp
[9124]314      !
315      DO jj = 2, jpj      ! Horizontal Flux in u and v direction
[9019]316         DO ji = 2, jpi 
[9124]317            !
318            IF( tmask(ji, jj, 1 ) < 0.5_wp) CYCLE   ! we don't care about land cells
319            IF( ht_0(ji,jj) > zdepwd )      CYCLE   ! and cells which are unlikely to dry
320            !
321            zflxp(ji,jj) = MAX( zflxu(ji,jj) , 0._wp ) - MIN( zflxu(ji-1,jj  ) , 0._wp )   &
322               &         + MAX( zflxv(ji,jj) , 0._wp ) - MIN( zflxv(ji,  jj-1) , 0._wp ) 
323            zflxn(ji,jj) = MIN( zflxu(ji,jj) , 0._wp ) - MAX( zflxu(ji-1,jj  ) , 0._wp )   &
324               &         + MIN( zflxv(ji,jj) , 0._wp ) - MAX( zflxv(ji,  jj-1) , 0._wp ) 
325            !
[9023]326            zdep2 = ht_0(ji,jj) + sshn_e(ji,jj) - rn_wdmin1
[9124]327            IF( zdep2 <= 0._wp ) THEN  !add more safety, but not necessary
[9023]328              sshn_e(ji,jj) = rn_wdmin1 - ht_0(ji,jj)
[9124]329              IF( zflxu(ji  ,jj  ) > 0._wp)   zwdlmtu(ji  ,jj  ) = 0._wp
330              IF( zflxu(ji-1,jj  ) < 0._wp)   zwdlmtu(ji-1,jj  ) = 0._wp
331              IF( zflxv(ji  ,jj  ) > 0._wp)   zwdlmtv(ji  ,jj  ) = 0._wp
332              IF( zflxv(ji  ,jj-1) < 0._wp)   zwdlmtv(ji  ,jj-1) = 0._wp 
333            ENDIF
334         END DO
[9019]335      END DO
[9124]336      !
337      DO jk1 = 1, nn_wdit + 1      !! start limiter iterations
338         !
[9019]339         zflxu1(:,:) = zflxu(:,:) * zwdlmtu(:,:)
340         zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:)
341         jflag = 0     ! flag indicating if any further iterations are needed
[9124]342         !
[9019]343         DO jj = 2, jpj
344            DO ji = 2, jpi 
[9124]345               !
346               IF( tmask(ji, jj, 1 ) < 0.5_wp )   CYCLE
347               IF( ht_0(ji,jj)       > zdepwd )   CYCLE
348               !
[9019]349               ztmp = e1e2t(ji,jj)
[9124]350               !
351               zzflxp = max(zflxu1(ji,jj), 0._wp) - min(zflxu1(ji-1,jj),   0._wp)   &
352                  &   + max(zflxv1(ji,jj), 0._wp) - min(zflxv1(ji,  jj-1), 0._wp) 
353               zzflxn = min(zflxu1(ji,jj), 0._wp) - max(zflxu1(ji-1,jj),   0._wp)   &
354                  &   + min(zflxv1(ji,jj), 0._wp) - max(zflxv1(ji,  jj-1), 0._wp) 
[6152]355         
[9019]356               zdep1 = (zzflxp + zzflxn) * z2dt / ztmp
[9023]357               zdep2 = ht_0(ji,jj) + sshn_e(ji,jj) - rn_wdmin1 - z2dt * zssh_frc(ji,jj)
[6152]358         
[9019]359               IF(zdep1 > zdep2) THEN
[9023]360                 zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt )
361                 !zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zzflxp * z2dt )
362                 ! flag if the limiter has been used but stop flagging if the only
363                 ! changes have zeroed the coefficient since further iterations will
364                 ! not change anything
365                 IF( zcoef > 0._wp ) THEN
366                    jflag = 1 
367                 ELSE
368                    zcoef = 0._wp
369                 ENDIF
370                 IF(jk1 > nn_wdit) zcoef = 0._wp
371                 IF(zflxu1(ji,  jj) > 0._wp) zwdlmtu(ji  ,jj) = zcoef
372                 IF(zflxu1(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = zcoef
373                 IF(zflxv1(ji,  jj) > 0._wp) zwdlmtv(ji  ,jj) = zcoef
374                 IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = zcoef
375               END IF
[9019]376            END DO ! ji loop
377         END DO  ! jj loop
[9124]378         !
[10425]379         CALL lbc_lnk_multi( 'wet_dry', zwdlmtu, 'U', 1., zwdlmtv, 'V', 1. )
[9124]380         !
[10425]381         CALL mpp_max('wet_dry', jflag)   !max over the global domain
[9124]382         !
383         IF(jflag == 0)   EXIT
384         !
[9019]385      END DO  ! jk1 loop
[9124]386      !
[9019]387      zflxu(:,:) = zflxu(:,:) * zwdlmtu(:, :) 
388      zflxv(:,:) = zflxv(:,:) * zwdlmtv(:, :) 
[9124]389      !
390!!gm THIS lbc_lnk is useless since it is already done at the end of the jk1-loop
[10425]391      CALL lbc_lnk_multi( 'wet_dry', zflxu, 'U', -1., zflxv, 'V', -1. )
[9124]392!!gm end
393      !
394      IF( jflag == 1 .AND. lwp )   WRITE(numout,*) 'Need more iterations in wad_lmt_bt!!!'
395      !
[11053]396      !IF( ln_rnf      )   CALL sbc_rnf_div( hdiv )          ! runoffs (update hdiv field)
[9019]397      !
[9124]398      IF( ln_timing )   CALL timing_stop('wad_lmt_bt')      !
[9019]399      !
[6152]400   END SUBROUTINE wad_lmt_bt
[7646]401
402   !!==============================================================================
[6152]403END MODULE wet_dry
Note: See TracBrowser for help on using the repository browser.