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

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

Last change on this file since 7030 was 7030, checked in by acc, 8 years ago

Branch dev_r6393_NOC_WAD. Added over-topping test case to WAD_TEST_CASES (jp_cfg=6). Also removed last references to redundant wduflt and wdvflt arrays
File size: 22.1 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(:,:) ::   wdmask         !: u- and v- limiter
36
37   LOGICAL,  PUBLIC  ::   ln_wd       !: Wetting/drying activation switch (T:on,F:off)
38   REAL(wp), PUBLIC  ::   rn_wdmin1   !: minimum water depth on dried cells
39   REAL(wp), PUBLIC  ::   rn_wdmin2   !: tolerrance of minimum water depth on dried cells
40   REAL(wp), PUBLIC  ::   rn_wdld     !: land elevation below which wetting/drying
41                                           !: will be considered
42   INTEGER , PUBLIC  ::   nn_wdit     !: maximum number of iteration for W/D limiter
43
44   PUBLIC   wad_init                  ! initialisation routine called by step.F90
45   PUBLIC   wad_lmt                   ! routine called by sshwzv.F90
46   PUBLIC   wad_lmt_bt                ! routine called by dynspg_ts.F90
47   PUBLIC   wad_istate                ! routine called by istate.F90 and domvvl.F90
48
49   !! * Substitutions
50#  include "vectopt_loop_substitute.h90"
51CONTAINS
52
53   SUBROUTINE wad_init
54      !!----------------------------------------------------------------------
55      !!                     ***  ROUTINE wad_init  ***
56      !!                   
57      !! ** Purpose :   read wetting and drying namelist and print the variables.
58      !!
59      !! ** input   : - namwad namelist
60      !!----------------------------------------------------------------------
61      NAMELIST/namwad/ ln_wd, rn_wdmin1, rn_wdmin2, rn_wdld, nn_wdit
62      INTEGER  ::   ios                 ! Local integer output status for namelist read
63      INTEGER  ::   ierr                ! Local integer status array allocation
64      !!----------------------------------------------------------------------
65
66      REWIND( numnam_ref )              ! Namelist namwad in reference namelist
67                                        ! : Parameters for Wetting/Drying
68      READ  ( numnam_ref, namwad, IOSTAT = ios, ERR = 905)
69905   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namwad in reference namelist', .TRUE.) 
70      REWIND( numnam_cfg )              ! Namelist namwad in configuration namelist
71                                        ! : Parameters for Wetting/Drying
72      READ  ( numnam_cfg, namwad, IOSTAT = ios, ERR = 906)
73906   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namwad in configuration namelist', .TRUE. )
74      IF(lwm) WRITE ( numond, namwad )
75
76      IF(lwp) THEN                  ! control print
77         WRITE(numout,*)
78         WRITE(numout,*) 'wad_init : Wetting and drying initialization through namelist read'
79         WRITE(numout,*) '~~~~~~~ '
80         WRITE(numout,*) '   Namelist namwad'
81         WRITE(numout,*) '      Logical activation                 ln_wd      = ', ln_wd
82         WRITE(numout,*) '      Minimum wet depth on dried cells rn_wdmin1    = ', rn_wdmin1
83         WRITE(numout,*) '      Tolerance of min wet depth     rn_wdmin2      = ', rn_wdmin2
84         WRITE(numout,*) '      land elevation threshold         rn_wdld      = ', rn_wdld
85         WRITE(numout,*) '      Max iteration for W/D limiter    nn_wdit      = ', nn_wdit
86       ENDIF
87
88      IF(ln_wd) THEN
89         ALLOCATE( wdmask(jpi,jpj), STAT=ierr )
90         IF( ierr /= 0 ) CALL ctl_stop('STOP', 'wad_init : Array allocation error')
91      ENDIF
92   END SUBROUTINE wad_init
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
122      IF( nn_timing == 1 )  CALL timing_start('wad_lmt')
123
124      IF(ln_wd) THEN
125
126        CALL wrk_alloc( jpi, jpj, zflxp, zflxn, zflxu, zflxv, zflxu1, zflxv1 )
127        CALL wrk_alloc( jpi, jpj, zwdlmtu, zwdlmtv)
128        !
129       
130        !IF(lwp) WRITE(numout,*)
131        !IF(lwp) WRITE(numout,*) 'wad_lmt : wetting/drying limiters and velocity limiting'
132       
133        zflag  = 0
134        zdepwd = 50._wp   !maximum depth on which that W/D could possibly happen
135
136       
137        zflxp(:,:)   = 0._wp
138        zflxn(:,:)   = 0._wp
139        zflxu(:,:)   = 0._wp
140        zflxv(:,:)   = 0._wp
141
142        zwdlmtu(:,:)  = 1._wp
143        zwdlmtv(:,:)  = 1._wp
144       
145        ! Horizontal Flux in u and v direction
146        DO jk = 1, jpkm1 
147           DO jj = 1, jpj
148              DO ji = 1, jpi
149                 zflxu(ji,jj) = zflxu(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)
150                 zflxv(ji,jj) = zflxv(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
151              END DO 
152           END DO 
153        END DO
154       
155        zflxu(:,:) = zflxu(:,:) * e2u(:,:)
156        zflxv(:,:) = zflxv(:,:) * e1v(:,:)
157       
158        wdmask(:,:) = 1
159        DO jj = 2, jpj
160           DO ji = 2, jpi 
161
162             IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE   ! we don't care about land cells
163             IF(bathy(ji,jj) > zdepwd) CYCLE       ! and cells which will unlikely go dried out
164
165              zflxp(ji,jj) = max(zflxu(ji,jj), 0._wp) - min(zflxu(ji-1,jj),   0._wp) + &
166                           & max(zflxv(ji,jj), 0._wp) - min(zflxv(ji,  jj-1), 0._wp) 
167              zflxn(ji,jj) = min(zflxu(ji,jj), 0._wp) - max(zflxu(ji-1,jj),   0._wp) + &
168                           & min(zflxv(ji,jj), 0._wp) - max(zflxv(ji,  jj-1), 0._wp) 
169
170              zdep2 = bathy(ji,jj) + sshb1(ji,jj) - rn_wdmin1
171              IF(zdep2 .le. 0._wp) THEN  !add more safty, but not necessary
172                !zdep2 = 0._wp
173                sshb1(ji,jj) = rn_wdmin1 - bathy(ji,jj)
174                wdmask(ji,jj) = 0._wp
175              END IF
176           ENDDO
177        END DO
178
179     
180        !! start limiter iterations
181        DO jk1 = 1, nn_wdit + 1
182       
183         
184           zflxu1(:,:) = zflxu(:,:) * zwdlmtu(:,:)
185           zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:)
186         
187           DO jj = 2, jpj
188              DO ji = 2, jpi 
189       
190                 IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE
191                 IF(bathy(ji,jj) > zdepwd) CYCLE
192       
193                 ztmp = e1e2t(ji,jj)
194
195                 zzflxp = max(zflxu1(ji,jj), 0._wp) - min(zflxu1(ji-1,jj),   0._wp) + &
196                        & max(zflxv1(ji,jj), 0._wp) - min(zflxv1(ji,  jj-1), 0._wp) 
197                 zzflxn = min(zflxu1(ji,jj), 0._wp) - max(zflxu1(ji-1,jj),   0._wp) + &
198                        & min(zflxv1(ji,jj), 0._wp) - max(zflxv1(ji,  jj-1), 0._wp) 
199         
200                 zdep1 = (zzflxp + zzflxn) * z2dt / ztmp
201                 zdep2 = bathy(ji,jj) + sshb1(ji,jj) - rn_wdmin1 - z2dt * sshemp(ji,jj)  ! this one can be moved out of the loop
202         
203                 IF(zdep1 > zdep2) THEN
204                   zflag = 1
205                   wdmask(ji, jj) = 0
206                   zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt )
207                   zcoef = max(zcoef, 0._wp)
208                   IF(jk1 > nn_wdit) zcoef = 0._wp
209                   IF(zflxu1(ji,  jj) > 0._wp) zwdlmtu(ji  ,jj) = zcoef
210                   IF(zflxu1(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = zcoef
211                   IF(zflxv1(ji,  jj) > 0._wp) zwdlmtv(ji  ,jj) = zcoef
212                   IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = zcoef
213                 END IF
214              END DO ! ji loop
215           END DO  ! jj loop
216
217           CALL lbc_lnk( zwdlmtu, 'U', 1. )
218           CALL lbc_lnk( zwdlmtv, 'V', 1. )
219
220           IF(lk_mpp) CALL mpp_max(zflag)   !max over the global domain
221
222           IF(zflag == 0) EXIT
223         
224           zflag = 0     ! flag indicating if any further iteration is needed?
225        END DO  ! jk1 loop
226       
227        DO jk = 1, jpkm1
228          un(:,:,jk) = un(:,:,jk) * zwdlmtu(:, :) 
229          vn(:,:,jk) = vn(:,:,jk) * zwdlmtv(:, :) 
230        END DO
231
232        CALL lbc_lnk( un, 'U', -1. )
233        CALL lbc_lnk( vn, 'V', -1. )
234      !
235        un_b(:,:) = un_b(:,:) * zwdlmtu(:, :)
236        vn_b(:,:) = vn_b(:,:) * zwdlmtv(:, :)
237        CALL lbc_lnk( un_b, 'U', -1. )
238        CALL lbc_lnk( vn_b, 'V', -1. )
239       
240        IF(zflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt!!!'
241       
242        !IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )          ! runoffs (update hdivn field)
243        !IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (update hdivn field)
244        !
245        !
246        CALL wrk_dealloc( jpi, jpj, zflxp, zflxn, zflxu, zflxv, zflxu1, zflxv1 )
247        CALL wrk_dealloc( jpi, jpj, zwdlmtu, zwdlmtv)
248      !
249      END IF
250
251      IF( nn_timing == 1 )  CALL timing_stop('wad_lmt')
252   END SUBROUTINE wad_lmt
253
254   SUBROUTINE wad_lmt_bt( zflxu, zflxv, sshn_e, zssh_frc, rdtbt )
255      !!----------------------------------------------------------------------
256      !!                  ***  ROUTINE wad_lmt  ***
257      !!                   
258      !! ** Purpose :   limiting flux in the barotropic stepping (dynspg_ts)
259      !!
260      !! ** Method  : - Prevent negative depth occurring (Not ready for Agrif)
261      !!
262      !! ** Action  : - calculate flux limiter and W/D flag
263      !!----------------------------------------------------------------------
264      REAL(wp), INTENT(in)    ::   rdtbt    ! ocean time-step index
265      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zflxu,  zflxv, sshn_e, zssh_frc 
266      !
267      INTEGER  ::   ji, jj, jk, jk1     ! dummy loop indices
268      INTEGER  ::   zflag         ! local scalar
269      REAL(wp) ::   z2dt
270      REAL(wp) ::   zcoef, zdep1, zdep2 ! local scalars
271      REAL(wp) ::   zzflxp, zzflxn      ! local scalars
272      REAL(wp) ::   zdepwd              ! local scalar, always wet cell depth
273      REAL(wp) ::   ztmp                ! local scalars
274      REAL(wp), POINTER,  DIMENSION(:,:) ::   zwdlmtu, zwdlmtv         !: W/D flux limiters
275      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxp,  zflxn            ! local 2D workspace
276      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxu1, zflxv1           ! local 2D workspace
277
278      !!----------------------------------------------------------------------
279      !
280
281      IF( nn_timing == 1 )  CALL timing_start('wad_lmt_bt')
282
283      IF(ln_wd) THEN
284
285        CALL wrk_alloc( jpi, jpj, zflxp, zflxn, zflxu1, zflxv1 )
286        CALL wrk_alloc( jpi, jpj, zwdlmtu, zwdlmtv)
287        !
288       
289        !IF(lwp) WRITE(numout,*)
290        !IF(lwp) WRITE(numout,*) 'wad_lmt_bt : wetting/drying limiters and velocity limiting'
291       
292        zflag  = 0
293        zdepwd = 50._wp   !maximum depth that ocean cells can have W/D processes
294
295        z2dt = rdtbt   
296       
297        zflxp(:,:)   = 0._wp
298        zflxn(:,:)   = 0._wp
299
300        zwdlmtu(:,:)  = 1._wp
301        zwdlmtv(:,:)  = 1._wp
302       
303        ! Horizontal Flux in u and v direction
304       
305        DO jj = 2, jpj
306           DO ji = 2, jpi 
307
308             IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE   ! we don't care about land cells
309             IF(bathy(ji,jj) > zdepwd) CYCLE       ! and cells which will unlikely go dried out
310
311              zflxp(ji,jj) = max(zflxu(ji,jj), 0._wp) - min(zflxu(ji-1,jj),   0._wp) + &
312                           & max(zflxv(ji,jj), 0._wp) - min(zflxv(ji,  jj-1), 0._wp) 
313              zflxn(ji,jj) = min(zflxu(ji,jj), 0._wp) - max(zflxu(ji-1,jj),   0._wp) + &
314                           & min(zflxv(ji,jj), 0._wp) - max(zflxv(ji,  jj-1), 0._wp) 
315
316              zdep2 = bathy(ji,jj) + sshn_e(ji,jj) - rn_wdmin1
317           ENDDO
318        END DO
319
320     
321        !! start limiter iterations
322        DO jk1 = 1, nn_wdit + 1
323       
324         
325           zflxu1(:,:) = zflxu(:,:) * zwdlmtu(:,:)
326           zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:)
327         
328           DO jj = 2, jpj
329              DO ji = 2, jpi 
330       
331                 IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE
332                 IF(bathy(ji,jj) > zdepwd) CYCLE
333       
334                 ztmp = e1e2t(ji,jj)
335
336                 zzflxp = max(zflxu1(ji,jj), 0._wp) - min(zflxu1(ji-1,jj),   0._wp) + &
337                        & max(zflxv1(ji,jj), 0._wp) - min(zflxv1(ji,  jj-1), 0._wp) 
338                 zzflxn = min(zflxu1(ji,jj), 0._wp) - max(zflxu1(ji-1,jj),   0._wp) + &
339                        & min(zflxv1(ji,jj), 0._wp) - max(zflxv1(ji,  jj-1), 0._wp) 
340         
341                 zdep1 = (zzflxp + zzflxn) * z2dt / ztmp
342                 zdep2 = bathy(ji,jj) + sshn_e(ji,jj) - rn_wdmin1   ! this one can be moved out of the loop
343                 zdep2 = zdep2 - z2dt * zssh_frc(ji,jj)
344         
345                 IF(zdep1 > zdep2) THEN
346                   zflag = 1
347                   zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt )
348                   zcoef = max(zcoef, 0._wp)
349                   IF(jk1 > nn_wdit) zcoef = 0._wp
350                   IF(zflxu1(ji,  jj) > 0._wp) zwdlmtu(ji  ,jj) = zcoef
351                   IF(zflxu1(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = zcoef
352                   IF(zflxv1(ji,  jj) > 0._wp) zwdlmtv(ji  ,jj) = zcoef
353                   IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = zcoef
354                 END IF
355              END DO ! ji loop
356           END DO  ! jj loop
357
358           CALL lbc_lnk( zwdlmtu, 'U', 1. )
359           CALL lbc_lnk( zwdlmtv, 'V', 1. )
360
361           IF(lk_mpp) CALL mpp_max(zflag)   !max over the global domain
362
363           IF(zflag == 0) EXIT
364         
365           zflag = 0     ! flag indicating if any further iteration is needed?
366        END DO  ! jk1 loop
367       
368        zflxu(:,:) = zflxu(:,:) * zwdlmtu(:, :) 
369        zflxv(:,:) = zflxv(:,:) * zwdlmtv(:, :) 
370
371        CALL lbc_lnk( zflxu, 'U', -1. )
372        CALL lbc_lnk( zflxv, 'V', -1. )
373       
374        IF(zflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt_bt!!!'
375       
376        !
377        !
378        CALL wrk_dealloc( jpi, jpj, zflxp, zflxn, zflxu1, zflxv1 )
379        CALL wrk_dealloc( jpi, jpj, zwdlmtu, zwdlmtv)
380      !
381      END IF
382
383      IF( nn_timing == 1 )  CALL timing_stop('wad_lmt')
384   END SUBROUTINE wad_lmt_bt
385
386   SUBROUTINE wad_istate
387      !!----------------------------------------------------------------------
388      !!                   ***  ROUTINE wad_istate  ***
389      !!
390      !! ** Purpose :   Initialization of the dynamics and tracers for WAD test
391      !!      configurations (channels or bowls with initial ssh gradients)
392      !!
393      !! ** Method  : - set temperature field
394      !!              - set salinity field
395      !!              - set ssh slope (needs to be repeated in domvvl_rst_init to
396      !!                               set vertical metrics )
397      !!----------------------------------------------------------------------
398      !
399      INTEGER  ::   ji, jj            ! dummy loop indices
400      REAL(wp) ::   zi, zj
401      !!----------------------------------------------------------------------
402      !
403      ! Uniform T & S in all test cases
404      tsn(:,:,:,jp_tem) = 10._wp
405      tsb(:,:,:,jp_tem) = 10._wp
406      tsn(:,:,:,jp_sal) = 35._wp
407      tsb(:,:,:,jp_sal) = 35._wp
408      SELECT CASE ( jp_cfg ) 
409         !                                        ! ====================
410         CASE ( 1 )                               ! WAD 1 configuration
411            !                                     ! ====================
412            !
413            IF(lwp) WRITE(numout,*)
414            IF(lwp) WRITE(numout,*) 'istate_wad : Closed box with EW linear bottom slope'
415            IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
416            !
417            do ji = 1,jpi
418             sshn(ji,:) = ( -5.5_wp + 5.5_wp*FLOAT(mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1)
419            end do
420            !                                     ! ====================
421         CASE ( 2 )                               ! WAD 2 configuration
422            !                                     ! ====================
423            !
424            IF(lwp) WRITE(numout,*)
425            IF(lwp) WRITE(numout,*) 'istate_wad : Parobolic EW channel, mid-range initial ssh slope'
426            IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
427            !
428            do ji = 1,jpi
429             sshn(ji,:) = ( -5.5_wp + 3.9_wp*FLOAT(jpidta - mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1)
430            end do
431            !                                     ! ====================
432         CASE ( 3 )                               ! WAD 3 configuration
433            !                                     ! ====================
434            !
435            IF(lwp) WRITE(numout,*)
436            IF(lwp) WRITE(numout,*) 'istate_wad : Parobolic EW channel, extreme initial ssh slope' 
437            IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
438            !
439            do ji = 1,jpi
440             sshn(ji,:) = ( -7.5_wp + 6.9_wp*FLOAT(jpidta - mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1)
441            end do
442
443            !
444            !                                     ! ====================
445         CASE ( 4 )                               ! WAD 4 configuration
446            !                                     ! ====================
447            !
448            IF(lwp) WRITE(numout,*)
449            IF(lwp) WRITE(numout,*) 'istate_wad : Parobolic bowl, mid-range initial ssh slope' 
450            IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
451            !
452            DO ji = 1, jpi
453               zi = MAX(1.0-FLOAT((mig(ji)-25)**2)/400.0, 0.0 )
454               DO jj = 1, jpj
455                  zj = MAX(1.0-FLOAT((mjg(jj)-17)**2)/144.0, 0.0 )
456                  sshn(ji,jj) = -8.5_wp + 8.5_wp*zi*zj
457               END DO
458            END DO
459
460            !
461            !                                    ! ===========================
462         CASE ( 5 )                              ! WAD 5 configuration
463            !                                    ! ====================
464            !
465            IF(lwp) WRITE(numout,*)
466            IF(lwp) WRITE(numout,*) 'istate_wad : Double slope with shelf'
467            IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
468            !
469            ! Needed rn_wdmin2 increased to 0.01 for this case?
470            do ji = 1,jpi
471             sshn(ji,:) = ( -5.5_wp + 9.0_wp*FLOAT(mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1)
472            end do
473
474            !
475            !                                     ! ===========================
476         CASE ( 6 )                               ! WAD 6 configuration
477            !                                     ! ====================
478            !
479            IF(lwp) WRITE(numout,*)
480            IF(lwp) WRITE(numout,*) 'istate_wad : Parobolic EW channel with gaussian ridge' 
481            IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
482            !
483            do ji = 1,jpi
484             !6a
485             sshn(ji,:) = ( -5.5_wp + 9.0_wp*FLOAT(jpidta - mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1)
486             !Some variations in initial slope that have been tested
487             !6b
488             !sshn(ji,:) = ( -5.5_wp + 6.5_wp*FLOAT(jpidta - mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1)
489             !6c
490             !sshn(ji,:) = ( -5.5_wp + 7.5_wp*FLOAT(jpidta - mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1)
491             !6d
492             !sshn(ji,:) = ( -4.5_wp + 8.0_wp*FLOAT(jpidta - mig(ji))/FLOAT(jpidta-1))*tmask(ji,:,1)
493            end do
494
495            !
496            !                                    ! ===========================
497         CASE DEFAULT                            ! NONE existing configuration
498            !                                    ! ===========================
499            WRITE(ctmp1,*) 'WAD test with a ', jp_cfg,' option is not coded'
500            !
501            CALL ctl_stop( ctmp1 )
502            !
503      END SELECT
504      !
505      ! Apply minimum wetdepth criterion
506      !
507      do jj = 1,jpj
508         do ji = 1,jpi
509            IF( bathy(ji,jj) + sshn(ji,jj) < rn_wdmin1 ) THEN
510               sshn(ji,jj) = tmask(ji,jj,1)*( rn_wdmin1 - bathy(ji,jj) )
511            ENDIF
512         end do
513      end do
514      sshb = sshn
515      ssha = sshn
516      !
517   END SUBROUTINE wad_istate
518
519   !!=====================================================================
520END MODULE wet_dry
Note: See TracBrowser for help on using the repository browser.