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/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/DYN/wet_dry.F90 @ 7433

Last change on this file since 7433 was 7433, checked in by timgraham, 8 years ago

#1811 Fix for wetting and drying
File size: 22.4 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 wrk_nemo        ! Memory Allocation
24   USE timing          ! Timing
25
26   IMPLICIT NONE
27   PRIVATE
28
29   !!----------------------------------------------------------------------
30   !! critical depths,filters, limiters,and masks for  Wetting and Drying
31   !! ---------------------------------------------------------------------
32
33   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) ::   wdmask         !: u- and v- limiter
34
35   LOGICAL,  PUBLIC  ::   ln_wd       !: Wetting/drying activation switch (T:on,F:off)
36   REAL(wp), PUBLIC  ::   rn_wdmin1   !: minimum water depth on dried cells
37   REAL(wp), PUBLIC  ::   rn_wdmin2   !: tolerrance of minimum water depth on dried cells
38   REAL(wp), PUBLIC  ::   rn_wdld     !: land elevation below which wetting/drying
39                                           !: will be considered
40   INTEGER , PUBLIC  ::   nn_wdit     !: maximum number of iteration for W/D limiter
41
42   PUBLIC   wad_init                  ! initialisation routine called by step.F90
43   PUBLIC   wad_lmt                   ! routine called by sshwzv.F90
44   PUBLIC   wad_lmt_bt                ! routine called by dynspg_ts.F90
45   PUBLIC   wad_istate                ! routine called by istate.F90 and domvvl.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), 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  ::   zflag               ! 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), POINTER,  DIMENSION(:,:) ::   zwdlmtu, zwdlmtv         !: W/D flux limiters
115      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxp,  zflxn            ! local 2D workspace
116      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxu,  zflxv            ! local 2D workspace
117      REAL(wp), POINTER,  DIMENSION(:,:) ::   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        CALL wrk_alloc( jpi,jpj,   zflxp, zflxn, zflxu, zflxv, zflxu1, zflxv1 )
126        CALL wrk_alloc( jpi,jpj,   zwdlmtu, zwdlmtv)
127        !
128       
129        !IF(lwp) WRITE(numout,*)
130        !IF(lwp) WRITE(numout,*) 'wad_lmt : wetting/drying limiters and velocity limiting'
131       
132        zflag  = 0
133        zdepwd = 50._wp   !maximum depth on which that W/D could possibly happen
134
135       
136        zflxp(:,:)   = 0._wp
137        zflxn(:,:)   = 0._wp
138        zflxu(:,:)   = 0._wp
139        zflxv(:,:)   = 0._wp
140
141        zwdlmtu(:,:)  = 1._wp
142        zwdlmtv(:,:)  = 1._wp
143       
144        ! Horizontal Flux in u and v direction
145        DO jk = 1, jpkm1 
146           DO jj = 1, jpj
147              DO ji = 1, jpi
148                 zflxu(ji,jj) = zflxu(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)
149                 zflxv(ji,jj) = zflxv(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
150              END DO 
151           END DO 
152        END DO
153       
154        zflxu(:,:) = zflxu(:,:) * e2u(:,:)
155        zflxv(:,:) = zflxv(:,:) * e1v(:,:)
156       
157        wdmask(:,:) = 1
158        DO jj = 2, jpj
159           DO ji = 2, jpi 
160
161             IF( tmask(ji,jj,1) == 0._wp  )   CYCLE   ! we don't care about land cells
162             IF( ht_0 (ji,jj)   >  zdepwd )   CYCLE   ! and cells which will unlikely go dried out
163
164              zflxp(ji,jj) = max(zflxu(ji,jj), 0._wp) - min(zflxu(ji-1,jj),   0._wp) + &
165                           & max(zflxv(ji,jj), 0._wp) - min(zflxv(ji,  jj-1), 0._wp) 
166              zflxn(ji,jj) = min(zflxu(ji,jj), 0._wp) - max(zflxu(ji-1,jj),   0._wp) + &
167                           & min(zflxv(ji,jj), 0._wp) - max(zflxv(ji,  jj-1), 0._wp) 
168
169              zdep2 = ht_0(ji,jj) + sshb1(ji,jj) - rn_wdmin1
170              IF(zdep2 < 0._wp) THEN  !add more safty, but not necessary
171                !zdep2 = 0._wp
172                sshb1(ji,jj) = rn_wdmin1 - ht_0(ji,jj)
173              END IF
174           ENDDO
175        END DO
176
177     
178        !! start limiter iterations
179        DO jk1 = 1, nn_wdit + 1
180       
181         
182           zflxu1(:,:) = zflxu(:,:) * zwdlmtu(:,:)
183           zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:)
184         
185           DO jj = 2, jpj
186              DO ji = 2, jpi 
187       
188                 wdmask(ji,jj) = 0
189                 IF( tmask(ji,jj,1) < 0.5_wp) CYCLE
190                 IF( ht_0(ji,jj) > zdepwd) CYCLE
191       
192                 ztmp = e1e2t(ji,jj)
193
194                 zzflxp = max(zflxu1(ji,jj), 0._wp) - min(zflxu1(ji-1,jj),   0._wp) + &
195                        & max(zflxv1(ji,jj), 0._wp) - min(zflxv1(ji,  jj-1), 0._wp) 
196                 zzflxn = min(zflxu1(ji,jj), 0._wp) - max(zflxu1(ji-1,jj),   0._wp) + &
197                        & min(zflxv1(ji,jj), 0._wp) - max(zflxv1(ji,  jj-1), 0._wp) 
198         
199                 zdep1 = (zzflxp + zzflxn) * z2dt / ztmp
200                 zdep2 = ht_0(ji,jj) + sshb1(ji,jj) - rn_wdmin1 - z2dt * sshemp(ji,jj)  ! this one can be moved out of the loop
201         
202                 IF(zdep1 > zdep2) THEN
203                   zflag = 1
204                   wdmask(ji, jj) = 0
205                   zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt )
206                   zcoef = max(zcoef, 0._wp)
207                   IF(jk1 > nn_wdit) zcoef = 0._wp
208                   IF(zflxu1(ji,  jj) > 0._wp) zwdlmtu(ji  ,jj) = zcoef
209                   IF(zflxu1(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = zcoef
210                   IF(zflxv1(ji,  jj) > 0._wp) zwdlmtv(ji  ,jj) = zcoef
211                   IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = zcoef
212                 END IF
213              END DO ! ji loop
214           END DO  ! jj loop
215
216           CALL lbc_lnk( zwdlmtu, 'U', 1. )
217           CALL lbc_lnk( zwdlmtv, 'V', 1. )
218
219           IF(lk_mpp) CALL mpp_max(zflag)   !max over the global domain
220
221           IF(zflag == 0) EXIT
222         
223           zflag = 0     ! flag indicating if any further iteration is needed?
224        END DO  ! jk1 loop
225       
226        DO jk = 1, jpkm1
227          un(:,:,jk) = un(:,:,jk) * zwdlmtu(:, :) 
228          vn(:,:,jk) = vn(:,:,jk) * zwdlmtv(:, :) 
229        END DO
230
231        CALL lbc_lnk( un, 'U', -1. )
232        CALL lbc_lnk( vn, 'V', -1. )
233      !
234        un_b(:,:) = un_b(:,:) * zwdlmtu(:, :)
235        vn_b(:,:) = vn_b(:,:) * zwdlmtv(:, :)
236        CALL lbc_lnk( un_b, 'U', -1. )
237        CALL lbc_lnk( vn_b, 'V', -1. )
238       
239        IF(zflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt!!!'
240       
241        !IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )          ! runoffs (update hdivn field)
242        !IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (update hdivn field)
243        !
244        !
245        CALL wrk_dealloc( jpi, jpj, zflxp, zflxn, zflxu, zflxv, zflxu1, zflxv1 )
246        CALL wrk_dealloc( jpi, jpj, zwdlmtu, zwdlmtv)
247        !
248      ENDIF
249      !
250      IF( nn_timing == 1 )  CALL timing_stop('wad_lmt')
251      !
252   END SUBROUTINE wad_lmt
253
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      IF( nn_timing == 1 )  CALL timing_start('wad_lmt_bt')
281
282      IF(ln_wd) THEN
283
284        CALL wrk_alloc( jpi, jpj, zflxp, zflxn, zflxu1, zflxv1 )
285        CALL wrk_alloc( jpi, jpj, zwdlmtu, zwdlmtv)
286        !
287       
288        !IF(lwp) WRITE(numout,*)
289        !IF(lwp) WRITE(numout,*) 'wad_lmt_bt : wetting/drying limiters and velocity limiting'
290       
291        zflag  = 0
292        zdepwd = 50._wp   !maximum depth that ocean cells can have W/D processes
293
294        z2dt = rdtbt   
295       
296        zflxp(:,:)   = 0._wp
297        zflxn(:,:)   = 0._wp
298        !zflxu(:,:)   = 0._wp
299        !zflxv(:,:)   = 0._wp
300
301        zwdlmtu(:,:)  = 1._wp
302        zwdlmtv(:,:)  = 1._wp
303       
304        ! Horizontal Flux in u and v direction
305       
306        DO jj = 2, jpj
307           DO ji = 2, jpi 
308
309             IF(tmask(ji,jj,1) < 0.5_wp) CYCLE   ! we don't care about land cells
310             IF(ht_0 (ji,jj)   > zdepwd) CYCLE       ! and cells which will unlikely go dried out
311
312              zflxp(ji,jj) = max(zflxu(ji,jj), 0._wp) - min(zflxu(ji-1,jj),   0._wp) + &
313                           & max(zflxv(ji,jj), 0._wp) - min(zflxv(ji,  jj-1), 0._wp) 
314              zflxn(ji,jj) = min(zflxu(ji,jj), 0._wp) - max(zflxu(ji-1,jj),   0._wp) + &
315                           & min(zflxv(ji,jj), 0._wp) - max(zflxv(ji,  jj-1), 0._wp) 
316
317              zdep2 = ht_0(ji,jj) + sshn_e(ji,jj) - rn_wdmin1
318           ENDDO
319        END DO
320
321     
322        !! start limiter iterations
323        DO jk1 = 1, nn_wdit + 1
324       
325         
326           zflxu1(:,:) = zflxu(:,:) * zwdlmtu(:,:)
327           zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:)
328         
329           DO jj = 2, jpj
330              DO ji = 2, jpi 
331       
332                 IF(tmask(ji,jj,1) < 0.5_wp) CYCLE
333                 IF(ht_0 (ji,jj)   > zdepwd) CYCLE
334       
335                 ztmp = e1e2t(ji,jj)
336
337                 zzflxp = max(zflxu1(ji,jj), 0._wp) - min(zflxu1(ji-1,jj),   0._wp) + &
338                        & max(zflxv1(ji,jj), 0._wp) - min(zflxv1(ji,  jj-1), 0._wp) 
339                 zzflxn = min(zflxu1(ji,jj), 0._wp) - max(zflxu1(ji-1,jj),   0._wp) + &
340                        & min(zflxv1(ji,jj), 0._wp) - max(zflxv1(ji,  jj-1), 0._wp) 
341         
342                 zdep1 = (zzflxp + zzflxn) * z2dt / ztmp
343                 zdep2 = ht_0(ji,jj) + sshn_e(ji,jj) - rn_wdmin1   ! this one can be moved out of the loop
344                 zdep2 = zdep2 - z2dt * zssh_frc(ji,jj)
345         
346                 IF(zdep1 > zdep2) THEN
347                   zflag = 1
348                   !wdmask(ji, jj) = 1
349                   zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt )
350                   zcoef = max(zcoef, 0._wp)
351                   IF(jk1 > nn_wdit) zcoef = 0._wp
352                   IF(zflxu1(ji,  jj) > 0._wp) zwdlmtu(ji  ,jj) = zcoef
353                   IF(zflxu1(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = zcoef
354                   IF(zflxv1(ji,  jj) > 0._wp) zwdlmtv(ji  ,jj) = zcoef
355                   IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = zcoef
356                 END IF
357              END DO ! ji loop
358           END DO  ! jj loop
359
360           CALL lbc_lnk( zwdlmtu, 'U', 1. )
361           CALL lbc_lnk( zwdlmtv, 'V', 1. )
362
363           IF(lk_mpp) CALL mpp_max(zflag)   !max over the global domain
364
365           IF(zflag == 0) EXIT
366         
367           zflag = 0     ! flag indicating if any further iteration is needed?
368        END DO  ! jk1 loop
369       
370        zflxu(:,:) = zflxu(:,:) * zwdlmtu(:, :) 
371        zflxv(:,:) = zflxv(:,:) * zwdlmtv(:, :) 
372
373        CALL lbc_lnk( zflxu, 'U', -1. )
374        CALL lbc_lnk( zflxv, 'V', -1. )
375       
376        IF(zflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt_bt!!!'
377       
378        !IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )          ! runoffs (update hdivn field)
379        !IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (update hdivn field)
380        !
381        !
382        CALL wrk_dealloc( jpi, jpj, zflxp, zflxn, zflxu1, zflxv1 )
383        CALL wrk_dealloc( jpi, jpj, zwdlmtu, zwdlmtv)
384        !
385      END IF
386
387      IF( nn_timing == 1 )  CALL timing_stop('wad_lmt')
388   END SUBROUTINE wad_lmt_bt
389
390   SUBROUTINE wad_istate
391      !!----------------------------------------------------------------------
392      !!                   ***  ROUTINE wad_istate  ***
393      !!
394      !! ** Purpose :   Initialization of the dynamics and tracers for WAD test
395      !!      configurations (channels or bowls with initial ssh gradients)
396      !!
397      !! ** Method  : - set temperature field
398      !!              - set salinity field
399      !!              - set ssh slope (needs to be repeated in domvvl_rst_init to
400      !!                               set vertical metrics )
401      !!----------------------------------------------------------------------
402      !
403      INTEGER  ::   ji, jj            ! dummy loop indices
404      REAL(wp) ::   zi, zj
405      !!----------------------------------------------------------------------
406      !
407      ! Uniform T & S in all test cases
408      tsn(:,:,:,jp_tem) = 10._wp
409      tsb(:,:,:,jp_tem) = 10._wp
410      tsn(:,:,:,jp_sal) = 35._wp
411      tsb(:,:,:,jp_sal) = 35._wp
412      SELECT CASE ( jp_cfg ) 
413         !                                        ! ====================
414         CASE ( 1 )                               ! WAD 1 configuration
415            !                                     ! ====================
416            !
417            IF(lwp) WRITE(numout,*)
418            IF(lwp) WRITE(numout,*) 'istate_wad : Closed box with EW linear bottom slope'
419            IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
420            !
421            do ji = 1,jpi
422             sshn(ji,:) = ( -5.5_wp + 5.5_wp*FLOAT(mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1)
423            end do
424            !                                     ! ====================
425         CASE ( 2 )                               ! WAD 2 configuration
426            !                                     ! ====================
427            !
428            IF(lwp) WRITE(numout,*)
429            IF(lwp) WRITE(numout,*) 'istate_wad : Parobolic EW channel, mid-range initial ssh slope'
430            IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
431            !
432            do ji = 1,jpi
433             sshn(ji,:) = ( -5.5_wp + 3.9_wp*FLOAT(jpidta - mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1)
434            end do
435            !                                     ! ====================
436         CASE ( 3 )                               ! WAD 3 configuration
437            !                                     ! ====================
438            !
439            IF(lwp) WRITE(numout,*)
440            IF(lwp) WRITE(numout,*) 'istate_wad : Parobolic EW channel, extreme initial ssh slope' 
441            IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
442            !
443            do ji = 1,jpi
444             sshn(ji,:) = ( -7.5_wp + 6.9_wp*FLOAT(jpidta - mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1)
445            end do
446
447            !
448            !                                     ! ====================
449         CASE ( 4 )                               ! WAD 4 configuration
450            !                                     ! ====================
451            !
452            IF(lwp) WRITE(numout,*)
453            IF(lwp) WRITE(numout,*) 'istate_wad : Parobolic bowl, mid-range initial ssh slope' 
454            IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
455            !
456            DO ji = 1, jpi
457               zi = MAX(1.0-FLOAT((mig(ji)-25)**2)/400.0, 0.0 )
458               DO jj = 1, jpj
459                  zj = MAX(1.0-FLOAT((mjg(jj)-17)**2)/144.0, 0.0 )
460                  sshn(ji,jj) = -8.5_wp + 8.5_wp*zi*zj
461               END DO
462            END DO
463
464            !
465            !                                    ! ===========================
466         CASE ( 5 )                              ! WAD 5 configuration
467            !                                    ! ====================
468            !
469            IF(lwp) WRITE(numout,*)
470            IF(lwp) WRITE(numout,*) 'istate_wad : Double slope with shelf'
471            IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
472            !
473            ! Needed rn_wdmin2 increased to 0.01 for this case?
474            do ji = 1,jpi
475             sshn(ji,:) = ( -5.5_wp + 9.0_wp*FLOAT(mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1)
476            end do
477
478            !
479            !                                     ! ===========================
480         CASE ( 6 )                               ! WAD 6 configuration
481            !                                     ! ====================
482            !
483            IF(lwp) WRITE(numout,*)
484            IF(lwp) WRITE(numout,*) 'istate_wad : Parobolic EW channel with gaussian ridge' 
485            IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
486            !
487            do ji = 1,jpi
488             !6a
489             sshn(ji,:) = ( -5.5_wp + 9.0_wp*FLOAT(jpidta - mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1)
490             !Some variations in initial slope that have been tested
491             !6b
492             !sshn(ji,:) = ( -5.5_wp + 6.5_wp*FLOAT(jpidta - mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1)
493             !6c
494             !sshn(ji,:) = ( -5.5_wp + 7.5_wp*FLOAT(jpidta - mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1)
495             !6d
496             !sshn(ji,:) = ( -4.5_wp + 8.0_wp*FLOAT(jpidta - mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1)
497            end do
498
499            !
500            !                                    ! ===========================
501         CASE DEFAULT                            ! NONE existing configuration
502            !                                    ! ===========================
503            WRITE(ctmp1,*) 'WAD test with a ', jp_cfg,' option is not coded'
504            !
505            CALL ctl_stop( ctmp1 )
506            !
507      END SELECT
508      !
509      ! Apply minimum wetdepth criterion
510      !
511      do jj = 1,jpj
512         do ji = 1,jpi
513            IF( bathy(ji,jj) + sshn(ji,jj) < rn_wdmin1 ) THEN
514               sshn(ji,jj) = tmask(ji,jj,1)*( rn_wdmin1 - bathy(ji,jj) )
515            ENDIF
516         end do
517      end do
518      sshb = sshn
519      ssha = sshn
520      !
521   END SUBROUTINE wad_istate
522
523   !!==============================================================================
524END MODULE wet_dry
Note: See TracBrowser for help on using the repository browser.