source: branches/2016/dev_r6393_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/wet_dry.F90 @ 6986

Last change on this file since 6986 was 6986, checked in by acc, 5 years ago

Branch dev_r6393_NOC_WAD. Introduced some WAD test cases, tidied and corrected WAD code

File size: 21.2 KB
Line 
1
2MODULE wet_dry
3   !!==============================================================================
4   !!                       ***  MODULE  wet_dry  ***
5   !! Wetting and drying includes initialisation routine and routines to
6   !! compute and apply flux limiters and preserve water depth positivity
7   !! only effects if wetting/drying is on (ln_wd == .true.)
8   !!==============================================================================
9   !! History :   
10   !!  NEMO      3.6  ! 2014-09  ((H.Liu)  Original code
11   !!                 ! will add the runoff and periodic BC case later
12   !!----------------------------------------------------------------------
13
14   !!----------------------------------------------------------------------
15   !!   wad_lmt    : Compute the horizontal flux limiter and the limited velocity
16   !!                when wetting and drying happens
17   !!----------------------------------------------------------------------
18   USE oce             ! ocean dynamics and tracers
19   USE dom_oce         ! ocean space and time domain
20   USE sbc_oce, ONLY : ln_rnf   ! surface boundary condition: ocean
21   USE sbcrnf          ! river runoff
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 wrk_nemo        ! Memory Allocation
26   USE timing          ! Timing
27
28   IMPLICIT NONE
29   PRIVATE
30
31   !!----------------------------------------------------------------------
32   !! critical depths,filters, limiters,and masks for  Wetting and Drying
33   !! ---------------------------------------------------------------------
34
35   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) ::   wduflt, wdvflt !: u- and v- filter
36   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) ::   wdmask         !: u- and v- limiter
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
42                                           !: will be considered
43   INTEGER , PUBLIC  ::   nn_wdit     !: maximum number of iteration for W/D limiter
44
45   PUBLIC   wad_init                  ! initialisation routine called by step.F90
46   PUBLIC   wad_lmt                   ! routine called by sshwzv.F90
47   PUBLIC   wad_lmt_bt                ! routine called by dynspg_ts.F90
48   PUBLIC   wad_istate                ! routine called by istate.F90 and domvvl.F90
49
50   !! * Substitutions
51#  include "vectopt_loop_substitute.h90"
52CONTAINS
53
54   SUBROUTINE wad_init
55      !!----------------------------------------------------------------------
56      !!                     ***  ROUTINE wad_init  ***
57      !!                   
58      !! ** Purpose :   read wetting and drying namelist and print the variables.
59      !!
60      !! ** input   : - namwad namelist
61      !!----------------------------------------------------------------------
62      NAMELIST/namwad/ ln_wd, rn_wdmin1, rn_wdmin2, rn_wdld, nn_wdit
63      INTEGER  ::   ios                 ! Local integer output status for namelist read
64      INTEGER  ::   ierr                ! Local integer status array allocation
65      !!----------------------------------------------------------------------
66
67      REWIND( numnam_ref )              ! Namelist namwad in reference namelist
68                                        ! : Parameters for Wetting/Drying
69      READ  ( numnam_ref, namwad, IOSTAT = ios, ERR = 905)
70905   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namwad in reference namelist', .TRUE.) 
71      REWIND( numnam_cfg )              ! Namelist namwad in configuration namelist
72                                        ! : Parameters for Wetting/Drying
73      READ  ( numnam_cfg, namwad, IOSTAT = ios, ERR = 906)
74906   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namwad in configuration namelist', .TRUE. )
75      IF(lwm) WRITE ( numond, namwad )
76
77      IF(lwp) THEN                  ! control print
78         WRITE(numout,*)
79         WRITE(numout,*) 'wad_init : Wetting and drying initialization through namelist read'
80         WRITE(numout,*) '~~~~~~~ '
81         WRITE(numout,*) '   Namelist namwad'
82         WRITE(numout,*) '      Logical activation                 ln_wd      = ', ln_wd
83         WRITE(numout,*) '      Minimum wet depth on dried cells rn_wdmin1    = ', rn_wdmin1
84         WRITE(numout,*) '      Tolerance of min wet depth     rn_wdmin2      = ', rn_wdmin2
85         WRITE(numout,*) '      land elevation threshold         rn_wdld      = ', rn_wdld
86         WRITE(numout,*) '      Max iteration for W/D limiter    nn_wdit      = ', nn_wdit
87       ENDIF
88
89      IF(ln_wd) THEN
90         ALLOCATE( wduflt(jpi,jpj), wdvflt(jpi,jpj), wdmask(jpi,jpj), STAT=ierr )
91         IF( ierr /= 0 ) CALL ctl_stop('STOP', 'wad_init : Array allocation error')
92      ENDIF
93   END SUBROUTINE wad_init
94
95   SUBROUTINE wad_lmt( sshb1, sshemp, z2dt )
96      !!----------------------------------------------------------------------
97      !!                  ***  ROUTINE wad_lmt  ***
98      !!                   
99      !! ** Purpose :   generate flux limiters for wetting/drying
100      !!
101      !! ** Method  : - Prevent negative depth occurring (Not ready for Agrif)
102      !!
103      !! ** Action  : - calculate flux limiter and W/D flag
104      !!----------------------------------------------------------------------
105      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   sshb1
106      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   sshemp
107      REAL(wp), INTENT(in) :: z2dt
108      !
109      INTEGER  ::   ji, jj, jk, jk1     ! dummy loop indices
110      INTEGER  ::   zflag               ! local scalar
111      REAL(wp) ::   zcoef, zdep1, zdep2 ! local scalars
112      REAL(wp) ::   zzflxp, zzflxn      ! local scalars
113      REAL(wp) ::   zdepwd              ! local scalar, always wet cell depth
114      REAL(wp) ::   ztmp                ! local scalars
115      REAL(wp), POINTER,  DIMENSION(:,:) ::   zwdlmtu, zwdlmtv         !: W/D flux limiters
116      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxp,  zflxn            ! local 2D workspace
117      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxu,  zflxv            ! local 2D workspace
118      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxu1, zflxv1           ! local 2D workspace
119
120      !!----------------------------------------------------------------------
121      !
122
123      IF( nn_timing == 1 )  CALL timing_start('wad_lmt')
124
125      IF(ln_wd) THEN
126
127        CALL wrk_alloc( jpi, jpj, zflxp, zflxn, zflxu, zflxv, zflxu1, zflxv1 )
128        CALL wrk_alloc( jpi, jpj, zwdlmtu, zwdlmtv)
129        !
130       
131        !IF(lwp) WRITE(numout,*)
132        !IF(lwp) WRITE(numout,*) 'wad_lmt : wetting/drying limiters and velocity limiting'
133       
134        zflag  = 0
135        zdepwd = 50._wp   !maximum depth on which that W/D could possibly happen
136
137       
138        zflxp(:,:)   = 0._wp
139        zflxn(:,:)   = 0._wp
140        zflxu(:,:)   = 0._wp
141        zflxv(:,:)   = 0._wp
142
143        zwdlmtu(:,:)  = 1._wp
144        zwdlmtv(:,:)  = 1._wp
145       
146        ! Horizontal Flux in u and v direction
147        DO jk = 1, jpkm1 
148           DO jj = 1, jpjm1
149              DO ji = 1, jpim1
150                 zflxu(ji,jj) = zflxu(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)
151                 zflxv(ji,jj) = zflxv(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
152              END DO 
153           END DO 
154        END DO
155       
156        zflxu(:,:) = zflxu(:,:) * e2u(:,:)
157        zflxv(:,:) = zflxv(:,:) * e1v(:,:)
158       
159        wdmask(:,:) = 1
160        DO jj = 2, jpjm1
161           DO ji = 2, jpim1 
162
163             IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE   ! we don't care about land cells
164             IF(bathy(ji,jj) > zdepwd) CYCLE       ! and cells which will unlikely go dried out
165
166              zflxp(ji,jj) = max(zflxu(ji,jj), 0._wp) - min(zflxu(ji-1,jj),   0._wp) + &
167                           & max(zflxv(ji,jj), 0._wp) - min(zflxv(ji,  jj-1), 0._wp) 
168              zflxn(ji,jj) = min(zflxu(ji,jj), 0._wp) - max(zflxu(ji-1,jj),   0._wp) + &
169                           & min(zflxv(ji,jj), 0._wp) - max(zflxv(ji,  jj-1), 0._wp) 
170
171              zdep2 = bathy(ji,jj) + sshb1(ji,jj) - rn_wdmin1
172              IF(zdep2 .le. 0._wp) THEN  !add more safty, but not necessary
173                !zdep2 = 0._wp
174                sshb1(ji,jj) = rn_wdmin1 - bathy(ji,jj)
175                wdmask(ji,jj) = 0._wp
176              END IF
177           ENDDO
178        END DO
179
180     
181        !! start limiter iterations
182        DO jk1 = 1, nn_wdit + 1
183       
184         
185           zflxu1(:,:) = zflxu(:,:) * zwdlmtu(:,:)
186           zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:)
187         
188           DO jj = 2, jpjm1
189              DO ji = 2, jpim1 
190       
191                 IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE
192                 IF(bathy(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 = bathy(ji,jj) + sshb1(ji,jj) - rn_wdmin1 - z2dt * sshemp(ji,jj)  ! this one can be moved out of the loop
203         
204                 IF(zdep1 > zdep2) THEN
205                   zflag = 1
206                   wdmask(ji, jj) = 0
207                   zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt )
208                   zcoef = max(zcoef, 0._wp)
209                   IF(jk1 > nn_wdit) zcoef = 0._wp
210                   IF(zflxu1(ji,  jj) > 0._wp) zwdlmtu(ji  ,jj) = zcoef
211                   IF(zflxu1(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = zcoef
212                   IF(zflxv1(ji,  jj) > 0._wp) zwdlmtv(ji  ,jj) = zcoef
213                   IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = zcoef
214                 END IF
215              END DO ! ji loop
216           END DO  ! jj loop
217
218           CALL lbc_lnk( zwdlmtu, 'U', 1. )
219           CALL lbc_lnk( zwdlmtv, 'V', 1. )
220
221           IF(lk_mpp) CALL mpp_max(zflag)   !max over the global domain
222
223           IF(zflag == 0) EXIT
224         
225           zflag = 0     ! flag indicating if any further iteration is needed?
226        END DO  ! jk1 loop
227       
228        DO jk = 1, jpkm1
229          un(:,:,jk) = un(:,:,jk) * zwdlmtu(:, :) 
230          vn(:,:,jk) = vn(:,:,jk) * zwdlmtv(:, :) 
231        END DO
232
233        CALL lbc_lnk( un, 'U', -1. )
234        CALL lbc_lnk( vn, 'V', -1. )
235      !
236        un_b(:,:) = un_b(:,:) * zwdlmtu(:, :)
237        vn_b(:,:) = vn_b(:,:) * zwdlmtv(:, :)
238        CALL lbc_lnk( un_b, 'U', -1. )
239        CALL lbc_lnk( vn_b, 'V', -1. )
240       
241        IF(zflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt!!!'
242       
243        !IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )          ! runoffs (update hdivn field)
244        !IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (update hdivn field)
245        !
246        !
247        CALL wrk_dealloc( jpi, jpj, zflxp, zflxn, zflxu, zflxv, zflxu1, zflxv1 )
248        CALL wrk_dealloc( jpi, jpj, zwdlmtu, zwdlmtv)
249      !
250      END IF
251
252      IF( nn_timing == 1 )  CALL timing_stop('wad_lmt')
253   END SUBROUTINE wad_lmt
254
255   SUBROUTINE wad_lmt_bt( zflxu, zflxv, sshn_e, zssh_frc, rdtbt )
256      !!----------------------------------------------------------------------
257      !!                  ***  ROUTINE wad_lmt  ***
258      !!                   
259      !! ** Purpose :   limiting flux in the barotropic stepping (dynspg_ts)
260      !!
261      !! ** Method  : - Prevent negative depth occurring (Not ready for Agrif)
262      !!
263      !! ** Action  : - calculate flux limiter and W/D flag
264      !!----------------------------------------------------------------------
265      REAL(wp), INTENT(in)    ::   rdtbt    ! ocean time-step index
266      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zflxu,  zflxv, sshn_e, zssh_frc 
267      !
268      INTEGER  ::   ji, jj, jk, jk1     ! dummy loop indices
269      INTEGER  ::   zflag         ! local scalar
270      REAL(wp) ::   z2dt
271      REAL(wp) ::   zcoef, zdep1, zdep2 ! local scalars
272      REAL(wp) ::   zzflxp, zzflxn      ! local scalars
273      REAL(wp) ::   zdepwd              ! local scalar, always wet cell depth
274      REAL(wp) ::   ztmp                ! local scalars
275      REAL(wp), POINTER,  DIMENSION(:,:) ::   zwdlmtu, zwdlmtv         !: W/D flux limiters
276      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxp,  zflxn            ! local 2D workspace
277      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxu1, zflxv1           ! local 2D workspace
278
279      !!----------------------------------------------------------------------
280      !
281
282      IF( nn_timing == 1 )  CALL timing_start('wad_lmt_bt')
283
284      IF(ln_wd) THEN
285
286        CALL wrk_alloc( jpi, jpj, zflxp, zflxn, zflxu1, zflxv1 )
287        CALL wrk_alloc( jpi, jpj, zwdlmtu, zwdlmtv)
288        !
289       
290        !IF(lwp) WRITE(numout,*)
291        !IF(lwp) WRITE(numout,*) 'wad_lmt_bt : wetting/drying limiters and velocity limiting'
292       
293        zflag  = 0
294        zdepwd = 50._wp   !maximum depth that ocean cells can have W/D processes
295
296        z2dt = rdtbt   
297       
298        zflxp(:,:)   = 0._wp
299        zflxn(:,:)   = 0._wp
300        !zflxu(:,:)   = 0._wp
301        !zflxv(:,:)   = 0._wp
302
303        zwdlmtu(:,:)  = 1._wp
304        zwdlmtv(:,:)  = 1._wp
305       
306        ! Horizontal Flux in u and v direction
307       
308        !zflxu(:,:) = zflxu(:,:) * e2u(:,:)
309        !zflxv(:,:) = zflxv(:,:) * e1v(:,:)
310       
311        DO jj = 2, jpjm1
312           DO ji = 2, jpim1 
313
314             IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE   ! we don't care about land cells
315             IF(bathy(ji,jj) > zdepwd) CYCLE       ! and cells which will unlikely go dried out
316
317              zflxp(ji,jj) = max(zflxu(ji,jj), 0._wp) - min(zflxu(ji-1,jj),   0._wp) + &
318                           & max(zflxv(ji,jj), 0._wp) - min(zflxv(ji,  jj-1), 0._wp) 
319              zflxn(ji,jj) = min(zflxu(ji,jj), 0._wp) - max(zflxu(ji-1,jj),   0._wp) + &
320                           & min(zflxv(ji,jj), 0._wp) - max(zflxv(ji,  jj-1), 0._wp) 
321
322              zdep2 = bathy(ji,jj) + sshn_e(ji,jj) - rn_wdmin1
323              IF(zdep2 < 0._wp) THEN  !add more safty, but not necessary
324                !zdep2 = 0._wp
325               sshn_e(ji,jj) = rn_wdmin1 - bathy(ji,jj)
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         
338           DO jj = 2, jpjm1
339              DO ji = 2, jpim1 
340       
341                 !wdmask(ji,jj) = 0
342                 IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE
343                 IF(bathy(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 = bathy(ji,jj) + sshn_e(ji,jj) - rn_wdmin1   ! this one can be moved out of the loop
354                 zdep2 = zdep2 - z2dt * zssh_frc(ji,jj)
355         
356                 IF(zdep1 > zdep2) THEN
357                   zflag = 1
358                   !wdmask(ji, jj) = 1
359                   zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt )
360                   zcoef = max(zcoef, 0._wp)
361                   IF(jk1 > nn_wdit) zcoef = 0._wp
362                   IF(zflxu1(ji,  jj) > 0._wp) zwdlmtu(ji  ,jj) = zcoef
363                   IF(zflxu1(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = zcoef
364                   IF(zflxv1(ji,  jj) > 0._wp) zwdlmtv(ji  ,jj) = zcoef
365                   IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = zcoef
366                 END IF
367              END DO ! ji loop
368           END DO  ! jj loop
369
370           CALL lbc_lnk( zwdlmtu, 'U', 1. )
371           CALL lbc_lnk( zwdlmtv, 'V', 1. )
372
373           IF(lk_mpp) CALL mpp_max(zflag)   !max over the global domain
374
375           IF(zflag == 0) EXIT
376         
377           zflag = 0     ! flag indicating if any further iteration is needed?
378        END DO  ! jk1 loop
379       
380        zflxu(:,:) = zflxu(:,:) * zwdlmtu(:, :) 
381        zflxv(:,:) = zflxv(:,:) * zwdlmtv(:, :) 
382
383        CALL lbc_lnk( zflxu, 'U', -1. )
384        CALL lbc_lnk( zflxv, 'V', -1. )
385       
386        IF(zflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt_bt!!!'
387       
388        !IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )          ! runoffs (update hdivn field)
389        !IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (update hdivn field)
390        !
391        !
392        CALL wrk_dealloc( jpi, jpj, zflxp, zflxn, zflxu1, zflxv1 )
393        CALL wrk_dealloc( jpi, jpj, zwdlmtu, zwdlmtv)
394      !
395      END IF
396
397      IF( nn_timing == 1 )  CALL timing_stop('wad_lmt')
398   END SUBROUTINE wad_lmt_bt
399
400   SUBROUTINE wad_istate
401      !!----------------------------------------------------------------------
402      !!                   ***  ROUTINE wad_istate  ***
403      !!
404      !! ** Purpose :   Initialization of the dynamics and tracers for WAD test
405      !!      configurations (channels or bowls with initial ssh gradients)
406      !!
407      !! ** Method  : - set temperature field
408      !!              - set salinity field
409      !!              - set ssh slope (needs to be repeated in domvvl_rst_init to
410      !!                               set vertical metrics )
411      !!----------------------------------------------------------------------
412      !
413      INTEGER  ::   ji, jj            ! dummy loop indices
414      !!----------------------------------------------------------------------
415      !
416      ! Uniform T & S in all test cases
417      tsn(:,:,:,jp_tem) = 10._wp
418      tsb(:,:,:,jp_tem) = 10._wp
419      tsn(:,:,:,jp_sal) = 35._wp
420      tsb(:,:,:,jp_sal) = 35._wp
421      SELECT CASE ( jp_cfg ) 
422         !                                        ! ====================
423         CASE ( 1 )                               ! WAD 1 configuration
424            !                                     ! ====================
425            !
426            IF(lwp) WRITE(numout,*)
427            IF(lwp) WRITE(numout,*) 'istate_wad : Closed box with EW linear bottom slope'
428            IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
429            !
430            do ji = 1,jpi
431             sshn(ji,:) = ( -5.5_wp + 5.5_wp*FLOAT(mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1)
432            end do
433            !                                     ! ====================
434         CASE ( 2 )                               ! WAD 2 configuration
435            !                                     ! ====================
436            !
437            IF(lwp) WRITE(numout,*)
438            IF(lwp) WRITE(numout,*) 'istate_wad : Parobolic EW channel, mid-range initial ssh slope'
439            IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
440            !
441            do ji = 1,jpi
442             sshn(ji,:) = ( -5.5_wp + 3.9_wp*FLOAT(jpidta - mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1)
443            end do
444            !                                     ! ====================
445         CASE ( 3 )                               ! WAD 3 configuration
446            !                                     ! ====================
447            !
448            IF(lwp) WRITE(numout,*)
449            IF(lwp) WRITE(numout,*) 'istate_wad : Parobolic EW channel, extreme initial ssh slope' 
450            IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
451            !
452            do ji = 1,jpi
453             sshn(ji,:) = ( -7.5_wp + 6.9_wp*FLOAT(jpidta - mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1)
454            end do
455
456            !
457            !                                     ! ====================
458         CASE ( 4 )                               ! WAD 4 configuration
459            !                                     ! ====================
460            !
461            IF(lwp) WRITE(numout,*)
462            IF(lwp) WRITE(numout,*) 'istate_wad : Parobolic bowl, mid-range initial ssh slope' 
463            IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
464            !
465            do ji = 1,jpi
466             sshn(ji,:) = ( -5.5_wp + 3.9_wp*FLOAT(jpidta - mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1)
467             ! very small displacement test:
468             !sshn(ji,:) = ( -0.05_wp + 0.05_wp*FLOAT(jpidta - mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1)
469            end do
470
471            !
472            !                                    ! ===========================
473         CASE DEFAULT                            ! NONE existing configuration
474            !                                    ! ===========================
475            WRITE(ctmp1,*) 'WAD test with a ', jp_cfg,' option is not coded'
476            !
477            CALL ctl_stop( ctmp1 )
478            !
479      END SELECT
480      !
481      ! Apply minimum wetdepth criterion
482      !
483      do jj = 1,jpj
484         do ji = 1,jpi
485            IF( bathy(ji,jj) + sshn(ji,jj) < 0.4_wp) THEN
486               sshn(ji,jj) = tmask(ji,jj,1)*(0.4_wp - bathy(ji,jj))
487            ENDIF
488         end do
489      end do
490      sshb = sshn
491      ssha = sshn
492      !
493   END SUBROUTINE wad_istate
494
495   !!=====================================================================
496END MODULE wet_dry
Note: See TracBrowser for help on using the repository browser.