source: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/wet_dry.F90 @ 7646

Last change on this file since 7646 was 7646, checked in by timgraham, 4 years ago

Merge of dev_merge_2016 into trunk. UPDATE TO ARCHFILES NEEDED for XIOS2.
LIM_SRC_s/limrhg.F90 to follow in next commit due to change of kind (I'm unable to do it in this commit).
Merged using the following steps:

1) svn merge —reintegrate svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk .
2) Resolve minor conflicts in sette.sh and namelist_cfg for ORCA2LIM3 (due to a change in trunk after branch was created)
3) svn commit
4) svn switch svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk
5) svn merge svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/2016/dev_merge_2016 .
6) At this stage I checked out a clean copy of the branch to compare against what is about to be committed to the trunk.
6) svn commit #Commit code to the trunk

In this commit I have also reverted a change to Fcheck_archfile.sh which was causing problems on the Paris machine.

File size: 18.1 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   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) ::   ht_wd          !: wetting and drying t-pt depths
35                                                                     !  (can include negative depths)
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
48   !! * Substitutions
49#  include "vectopt_loop_substitute.h90"
50CONTAINS
51
52   SUBROUTINE wad_init
53      !!----------------------------------------------------------------------
54      !!                     ***  ROUTINE wad_init  ***
55      !!                   
56      !! ** Purpose :   read wetting and drying namelist and print the variables.
57      !!
58      !! ** input   : - namwad namelist
59      !!----------------------------------------------------------------------
60      NAMELIST/namwad/ ln_wd, rn_wdmin1, rn_wdmin2, rn_wdld, nn_wdit
61      INTEGER  ::   ios                 ! Local integer output status for namelist read
62      INTEGER  ::   ierr                ! Local integer status array allocation
63      !!----------------------------------------------------------------------
64
65      REWIND( numnam_ref )              ! Namelist namwad in reference namelist
66                                        ! : Parameters for Wetting/Drying
67      READ  ( numnam_ref, namwad, IOSTAT = ios, ERR = 905)
68905   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namwad in reference namelist', .TRUE.) 
69      REWIND( numnam_cfg )              ! Namelist namwad in configuration namelist
70                                        ! : Parameters for Wetting/Drying
71      READ  ( numnam_cfg, namwad, IOSTAT = ios, ERR = 906)
72906   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namwad in configuration namelist', .TRUE. )
73      IF(lwm) WRITE ( numond, namwad )
74
75      IF(lwp) THEN                  ! control print
76         WRITE(numout,*)
77         WRITE(numout,*) 'wad_init : Wetting and drying initialization through namelist read'
78         WRITE(numout,*) '~~~~~~~~'
79         WRITE(numout,*) '   Namelist namwad'
80         WRITE(numout,*) '      Logical activation                 ln_wd      = ', ln_wd
81         WRITE(numout,*) '      Minimum wet depth on dried cells rn_wdmin1    = ', rn_wdmin1
82         WRITE(numout,*) '      Tolerance of min wet depth     rn_wdmin2      = ', rn_wdmin2
83         WRITE(numout,*) '      land elevation threshold         rn_wdld      = ', rn_wdld
84         WRITE(numout,*) '      Max iteration for W/D limiter    nn_wdit      = ', nn_wdit
85      ENDIF
86      !
87      IF(ln_wd) THEN
88         ALLOCATE( wdmask(jpi,jpj), ht_wd(jpi,jpj),  STAT=ierr )
89         IF( ierr /= 0 ) CALL ctl_stop('STOP', 'wad_init : Array allocation error')
90      ENDIF
91      !
92   END SUBROUTINE wad_init
93
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  ::   jflag               ! 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      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        jflag  = 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( ht_wd(ji,jj) > zdepwd )      CYCLE   ! and cells which are unlikely to dry
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 = ht_wd(ji,jj) + sshb1(ji,jj) - rn_wdmin1
171              IF(zdep2 .le. 0._wp) THEN  !add more safty, but not necessary
172                sshb1(ji,jj) = rn_wdmin1 - ht_wd(ji,jj)
173                IF(zflxu(ji,  jj) > 0._wp) zwdlmtu(ji  ,jj) = 0._wp
174                IF(zflxu(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = 0._wp
175                IF(zflxv(ji,  jj) > 0._wp) zwdlmtv(ji  ,jj) = 0._wp
176                IF(zflxv(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = 0._wp 
177                wdmask(ji,jj) = 0._wp
178              END IF
179           ENDDO
180        END DO
181
182     
183        !! start limiter iterations
184        DO jk1 = 1, nn_wdit + 1
185       
186         
187           zflxu1(:,:) = zflxu(:,:) * zwdlmtu(:,:)
188           zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:)
189           jflag = 0     ! flag indicating if any further iterations are needed
190         
191           DO jj = 2, jpj
192              DO ji = 2, jpi 
193       
194                 IF( tmask(ji, jj, 1) < 0.5_wp ) CYCLE
195                 IF( ht_wd(ji,jj) > zdepwd )      CYCLE
196       
197                 ztmp = e1e2t(ji,jj)
198
199                 zzflxp = max(zflxu1(ji,jj), 0._wp) - min(zflxu1(ji-1,jj),   0._wp) + &
200                        & max(zflxv1(ji,jj), 0._wp) - min(zflxv1(ji,  jj-1), 0._wp) 
201                 zzflxn = min(zflxu1(ji,jj), 0._wp) - max(zflxu1(ji-1,jj),   0._wp) + &
202                        & min(zflxv1(ji,jj), 0._wp) - max(zflxv1(ji,  jj-1), 0._wp) 
203         
204                 zdep1 = (zzflxp + zzflxn) * z2dt / ztmp
205                 zdep2 = ht_wd(ji,jj) + sshb1(ji,jj) - rn_wdmin1 - z2dt * sshemp(ji,jj)
206         
207                 IF( zdep1 > zdep2 ) THEN
208                   wdmask(ji, jj) = 0
209                   zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt )
210                   !zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zzflxp * z2dt )
211                   ! flag if the limiter has been used but stop flagging if the only
212                   ! changes have zeroed the coefficient since further iterations will
213                   ! not change anything
214                   IF( zcoef > 0._wp ) THEN
215                      jflag = 1 
216                   ELSE
217                      zcoef = 0._wp
218                   ENDIF
219                   IF(jk1 > nn_wdit) zcoef = 0._wp
220                   IF(zflxu1(ji,  jj) > 0._wp) zwdlmtu(ji  ,jj) = zcoef
221                   IF(zflxu1(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = zcoef
222                   IF(zflxv1(ji,  jj) > 0._wp) zwdlmtv(ji  ,jj) = zcoef
223                   IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = zcoef
224                 END IF
225              END DO ! ji loop
226           END DO  ! jj loop
227
228           CALL lbc_lnk( zwdlmtu, 'U', 1. )
229           CALL lbc_lnk( zwdlmtv, 'V', 1. )
230
231           IF(lk_mpp) CALL mpp_max(jflag)   !max over the global domain
232
233           IF(jflag == 0) EXIT
234         
235        END DO  ! jk1 loop
236       
237        DO jk = 1, jpkm1
238          un(:,:,jk) = un(:,:,jk) * zwdlmtu(:, :) 
239          vn(:,:,jk) = vn(:,:,jk) * zwdlmtv(:, :) 
240        END DO
241
242        CALL lbc_lnk( un, 'U', -1. )
243        CALL lbc_lnk( vn, 'V', -1. )
244      !
245        un_b(:,:) = un_b(:,:) * zwdlmtu(:, :)
246        vn_b(:,:) = vn_b(:,:) * zwdlmtv(:, :)
247        CALL lbc_lnk( un_b, 'U', -1. )
248        CALL lbc_lnk( vn_b, 'V', -1. )
249       
250        IF(jflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt!!!'
251       
252        !IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )          ! runoffs (update hdivn field)
253        !IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (update hdivn field)
254        !
255        !
256        CALL wrk_dealloc( jpi, jpj, zflxp, zflxn, zflxu, zflxv, zflxu1, zflxv1 )
257        CALL wrk_dealloc( jpi, jpj, zwdlmtu, zwdlmtv)
258        !
259      ENDIF
260      !
261      IF( nn_timing == 1 )  CALL timing_stop('wad_lmt')
262      !
263   END SUBROUTINE wad_lmt
264
265
266   SUBROUTINE wad_lmt_bt( zflxu, zflxv, sshn_e, zssh_frc, rdtbt )
267      !!----------------------------------------------------------------------
268      !!                  ***  ROUTINE wad_lmt  ***
269      !!                   
270      !! ** Purpose :   limiting flux in the barotropic stepping (dynspg_ts)
271      !!
272      !! ** Method  : - Prevent negative depth occurring (Not ready for Agrif)
273      !!
274      !! ** Action  : - calculate flux limiter and W/D flag
275      !!----------------------------------------------------------------------
276      REAL(wp), INTENT(in)    ::   rdtbt    ! ocean time-step index
277      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zflxu,  zflxv, sshn_e, zssh_frc 
278      !
279      INTEGER  ::   ji, jj, jk, jk1     ! dummy loop indices
280      INTEGER  ::   jflag               ! local scalar
281      REAL(wp) ::   z2dt
282      REAL(wp) ::   zcoef, zdep1, zdep2 ! local scalars
283      REAL(wp) ::   zzflxp, zzflxn      ! local scalars
284      REAL(wp) ::   zdepwd              ! local scalar, always wet cell depth
285      REAL(wp) ::   ztmp                ! local scalars
286      REAL(wp), POINTER,  DIMENSION(:,:) ::   zwdlmtu, zwdlmtv         !: W/D flux limiters
287      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxp,  zflxn            ! local 2D workspace
288      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxu1, zflxv1           ! local 2D workspace
289      !!----------------------------------------------------------------------
290      !
291      IF( nn_timing == 1 )  CALL timing_start('wad_lmt_bt')
292
293      IF(ln_wd) THEN
294
295        CALL wrk_alloc( jpi, jpj, zflxp, zflxn, zflxu1, zflxv1 )
296        CALL wrk_alloc( jpi, jpj, zwdlmtu, zwdlmtv)
297        !
298       
299        !IF(lwp) WRITE(numout,*)
300        !IF(lwp) WRITE(numout,*) 'wad_lmt_bt : wetting/drying limiters and velocity limiting'
301       
302        jflag  = 0
303        zdepwd = 50._wp   !maximum depth that ocean cells can have W/D processes
304
305        z2dt = rdtbt   
306       
307        zflxp(:,:)   = 0._wp
308        zflxn(:,:)   = 0._wp
309
310        zwdlmtu(:,:)  = 1._wp
311        zwdlmtv(:,:)  = 1._wp
312       
313        ! Horizontal Flux in u and v direction
314       
315        DO jj = 2, jpj
316           DO ji = 2, jpi 
317
318             IF( tmask(ji, jj, 1 ) < 0.5_wp) CYCLE   ! we don't care about land cells
319             IF( ht_wd(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
326              zdep2 = ht_wd(ji,jj) + sshn_e(ji,jj) - rn_wdmin1
327              IF(zdep2 .le. 0._wp) THEN  !add more safety, but not necessary
328                sshn_e(ji,jj) = rn_wdmin1 - ht_wd(ji,jj)
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              END IF
334           ENDDO
335        END DO
336
337     
338        !! start limiter iterations
339        DO jk1 = 1, nn_wdit + 1
340       
341         
342           zflxu1(:,:) = zflxu(:,:) * zwdlmtu(:,:)
343           zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:)
344           jflag = 0     ! flag indicating if any further iterations are needed
345         
346           DO jj = 2, jpj
347              DO ji = 2, jpi 
348       
349                 IF( tmask(ji, jj, 1 ) < 0.5_wp) CYCLE
350                 IF( ht_wd(ji,jj) > zdepwd )      CYCLE
351       
352                 ztmp = e1e2t(ji,jj)
353
354                 zzflxp = max(zflxu1(ji,jj), 0._wp) - min(zflxu1(ji-1,jj),   0._wp) + &
355                        & max(zflxv1(ji,jj), 0._wp) - min(zflxv1(ji,  jj-1), 0._wp) 
356                 zzflxn = min(zflxu1(ji,jj), 0._wp) - max(zflxu1(ji-1,jj),   0._wp) + &
357                        & min(zflxv1(ji,jj), 0._wp) - max(zflxv1(ji,  jj-1), 0._wp) 
358         
359                 zdep1 = (zzflxp + zzflxn) * z2dt / ztmp
360                 zdep2 = ht_wd(ji,jj) + sshn_e(ji,jj) - rn_wdmin1 - z2dt * zssh_frc(ji,jj)
361         
362                 IF(zdep1 > zdep2) THEN
363                   zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt )
364                   !zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zzflxp * z2dt )
365                   ! flag if the limiter has been used but stop flagging if the only
366                   ! changes have zeroed the coefficient since further iterations will
367                   ! not change anything
368                   IF( zcoef > 0._wp ) THEN
369                      jflag = 1 
370                   ELSE
371                      zcoef = 0._wp
372                   ENDIF
373                   IF(jk1 > nn_wdit) zcoef = 0._wp
374                   IF(zflxu1(ji,  jj) > 0._wp) zwdlmtu(ji  ,jj) = zcoef
375                   IF(zflxu1(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = zcoef
376                   IF(zflxv1(ji,  jj) > 0._wp) zwdlmtv(ji  ,jj) = zcoef
377                   IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji,jj-1) = zcoef
378                 END IF
379              END DO ! ji loop
380           END DO  ! jj loop
381
382           CALL lbc_lnk( zwdlmtu, 'U', 1. )
383           CALL lbc_lnk( zwdlmtv, 'V', 1. )
384
385           IF(lk_mpp) CALL mpp_max(jflag)   !max over the global domain
386
387           IF(jflag == 0) EXIT
388         
389        END DO  ! jk1 loop
390       
391        zflxu(:,:) = zflxu(:,:) * zwdlmtu(:, :) 
392        zflxv(:,:) = zflxv(:,:) * zwdlmtv(:, :) 
393
394        CALL lbc_lnk( zflxu, 'U', -1. )
395        CALL lbc_lnk( zflxv, 'V', -1. )
396       
397        IF(jflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt_bt!!!'
398       
399        !IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )          ! runoffs (update hdivn field)
400        !IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (update hdivn field)
401        !
402        !
403        CALL wrk_dealloc( jpi, jpj, zflxp, zflxn, zflxu1, zflxv1 )
404        CALL wrk_dealloc( jpi, jpj, zwdlmtu, zwdlmtv)
405        !
406      END IF
407
408      IF( nn_timing == 1 )  CALL timing_stop('wad_lmt')
409   END SUBROUTINE wad_lmt_bt
410
411   !!==============================================================================
412END MODULE wet_dry
Note: See TracBrowser for help on using the repository browser.