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.
wadlmt.F90 in branches/2013/dev_r4050_NOC_WaD/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2013/dev_r4050_NOC_WaD/NEMOGCM/NEMO/OPA_SRC/DYN/wadlmt.F90 @ 4389

Last change on this file since 4389 was 4389, checked in by hliu, 10 years ago

add a missed loop in wadlmt.F90

File size: 8.0 KB
Line 
1
2MODULE wadlmt
3   !!==============================================================================
4   !!                       ***  MODULE  wadlmt  ***
5   !! compute the flux limiter for water depth positivity in !wetting/drying case
6   !!==============================================================================
7   !! History :   
8   !!  NEMO      3.5  ! 2014-01  ((H.Liu)  Original code
9   !!----------------------------------------------------------------------
10
11   !!----------------------------------------------------------------------
12   !!   wad_lmt    : Compute the horizontal flux limiter and the limited velocity
13   !!                when wetting and drying happens
14   !!----------------------------------------------------------------------
15   USE oce             ! ocean dynamics and tracers
16   USE dom_oce         ! ocean space and time domain
17   USE sbc_oce, ONLY : ln_rnf   ! surface boundary condition: ocean
18   USE sbcrnf          ! river runoff
19   USE cla             ! cross land advection             (cla_div routine)
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   PUBLIC   wad_lmt    ! routine called by step.F90
31
32   !! * Substitutions
33#  include "domzgr_substitute.h90"
34#  include "vectopt_loop_substitute.h90"
35CONTAINS
36
37   SUBROUTINE wad_lmt( kt )
38      !!----------------------------------------------------------------------
39      !!                  ***  ROUTINE wad_lmt  ***
40      !!                   
41      !! ** Purpose :   generate flux limiters for wetting/drying
42      !!
43      !! ** Method  : - Prevent negative depth occurring (Not ready for Agrif)
44      !!
45      !! ** Action  : - calculate flux limiter and W/D flag
46      !!----------------------------------------------------------------------
47      INTEGER, INTENT(in) ::   kt       ! ocean time-step index
48      !
49      INTEGER  ::   ji, jj, jk, jk1     ! dummy loop indices
50      INTEGER  ::   z2dt, zflag         ! local scalar
51      REAL(wp) ::   zcoef, zdep1, zdep2 ! local scalars
52      REAL(wp) ::   zzflxp, zzflxn      ! local scalars
53      REAL(wp) ::   ztmp                ! local scalars
54      INTEGER,  POINTER,  DIMENSION(:,:) ::   iwdflg                   !: W/D cell flag (should become global)
55      REAL(wp), POINTER,  DIMENSION(:,:) ::   zwdlmtu, zwdlmtv         !: W/D flux limiters
56      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxp,  zflxn            ! specific 2D workspace
57      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxu,  zflxv            ! specific 2D workspace
58      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxu1, zflxv1           ! specific 2D workspace
59
60      !!----------------------------------------------------------------------
61      !
62
63      IF( nn_timing == 1 )  CALL timing_start('wad_lmt')
64
65      IF(ln_wad) THEN
66
67        CALL wrk_alloc( jpi, jpj, zflxp, zflxn, zflxu, zflxv, zflxu1, zflxv1 )
68        CALL wrk_alloc( jpi, jpj, zwdlmtu, zwdlmtv)
69        CALL wrk_alloc( jpi, jpj, iwdflg )
70        !
71       
72        IF(lwp) WRITE(numout,*)
73        IF(lwp) WRITE(numout,*) 'wad_lmt : wetting/drying limiters and velocity limiting'
74       
75        z2dt = 2. * rdt                                 ! Euler or leap-frog time step
76        IF( neuler == 0 .AND. kt == nit000 )  z2dt = rdt
77       
78        zflxp(:,:)   = 0._wp
79        zflxn(:,:)   = 0._wp
80        zflxu(:,:)   = 0._wp
81        zflxv(:,:)   = 0._wp
82
83        zwdlmtu(:,:)  = 1._wp
84        zwdlmtv(:,:)  = 1._wp
85       
86        ! Horizontal Flux in u and v direction
87        DO jk = 1, jpkm1 
88           DO jj = 1, jpjm1
89              DO ji = 1, jpim1
90                 zflxu(ji,jj) = zflxu(ji,jj) + fse3u(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)
91                 zflxv(ji,jj) = zflxv(ji,jj) + fse3v(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
92              END DO 
93           END DO 
94        END DO
95       
96        zflxu(:,:) = zflxu(:,:) * e2u(:,:)
97        zflxv(:,:) = zflxv(:,:) * e1v(:,:)
98       
99        DO jj = 2, jpjm1
100           DO ji = fs_2, fs_jpim1   ! vector opt.
101
102              IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE
103
104              zflxp(ji,jj) = max(zflxu(ji,jj), 0._wp) - min(zflxu(ji-1,jj),   0._wp) + &
105                           & max(zflxv(ji,jj), 0._wp) - min(zflxv(ji,  jj-1), 0._wp) 
106              zflxn(ji,jj) = min(zflxu(ji,jj), 0._wp) - max(zflxu(ji-1,jj),   0._wp) + &
107                           & min(zflxv(ji,jj), 0._wp) - max(zflxv(ji,  jj-1), 0._wp) 
108
109              zdep2 = bathy(ji,jj) + sshb(ji,jj) - rn_wadmin 
110              IF(zdep2 < 0._wp) THEN
111                zdep2 = 0._wp
112                sshb(ji,jj) = rn_wadmin - bathy(ji,jj)
113              END IF
114           ENDDO
115        END DO
116
117     
118        !! start limiter iterations
119        DO jk1 = 1, nn_waditr
120       
121           zflag = 0     ! flag indicating if any further iteration is needed?
122         
123           zflxu1(:,:) = zflxu(:,:) * zwdlmtu(:,:)
124           zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:)
125         
126           DO jj = 2, jpjm1
127              DO ji = fs_2, fs_jpim1   ! vector opt.
128       
129                 IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE
130       
131                 ztmp = e1t(ji,jj) * e2t(ji,jj)     !there must be an array ready for this
132
133                 zzflxp = max(zflxu1(ji,jj), 0._wp) - min(zflxu1(ji-1,jj),   0._wp) + &
134                        & max(zflxv1(ji,jj), 0._wp) - min(zflxv1(ji,  jj-1), 0._wp) 
135                 zzflxn = min(zflxu1(ji,jj), 0._wp) - max(zflxu1(ji-1,jj),   0._wp) + &
136                        & min(zflxv1(ji,jj), 0._wp) - max(zflxv1(ji,  jj-1), 0._wp) 
137         
138                 zdep1 = (zzflxp + zzflxn) * z2dt / ztmp
139                 zdep2 = bathy(ji,jj) + sshb(ji,jj) - rn_wadmin   ! this one can be moved out of the loops
140         
141                 IF(zdep1 > zdep2) THEN
142                   zflag = 1
143                   zcoef = ( ( zdep2 - rn_wadmine ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt )
144                   zcoef = max(zcoef, 0._wp)
145                   iwdflg(ji, jj) = 1
146                   IF(zflxu1(ji,  jj) > 0._wp) zwdlmtu(ji  ,jj) = zcoef
147                   IF(zflxu1(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = zcoef
148                   IF(zflxv1(ji,  jj) > 0._wp) zwdlmtv(ji  ,jj) = zcoef
149                   IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji-1,jj) = zcoef
150                 END IF
151              END DO ! ji loop
152           END DO  ! jj loop
153         
154           CALL lbc_lnk( zwdlmtu, 'U', 1. )
155           CALL lbc_lnk( zwdlmtv, 'V', 1. )
156         
157           IF(zflag == 0) EXIT
158         
159           IF(jk1 == nn_waditr) THEN
160              DO jj = 2, jpjm1
161                 DO ji = fs_2, fs_jpim1   ! vector opt.
162                 
163                 IF(zflxu1(ji,  jj  ) >= 0._wp) zwdlmtu(ji,  jj  ) = 0._wp
164                 IF(zflxu1(ji-1,jj  ) <  0._wp) zwdlmtu(ji-1,jj  ) = 0._wp
165                 IF(zflxv1(ji,  jj  ) >= 0._wp) zwdlmtv(ji,  jj  ) = 0._wp
166                 IF(zflxv1(ji,  jj-1) <  0._wp) zwdlmtv(ji,  jj-1) = 0._wp
167           
168                 CALL lbc_lnk( zwdlmtu, 'U', 1. )
169                 CALL lbc_lnk( zwdlmtv, 'V', 1. )
170                 END DO ! ji loop
171              END DO  ! jj loop
172           END IF
173
174        END DO  ! jk1 loop
175       
176        un(:,:,:) = un(:,:,:) * zwdlmtu(ji, jj) 
177        vn(:,:,:) = vn(:,:,:) * zwdlmtv(ji, jj) 
178
179        CALL lbc_lnk( un, 'U', -1. )
180        CALL lbc_lnk( vn, 'V', -1. )
181       
182        IF(zflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt!!!'
183       
184        !IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )          ! runoffs (update hdivn field)
185        !IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (update hdivn field)
186        !
187        !
188        CALL wrk_dealloc( jpi, jpj, zflxp, zflxn, zflxu, zflxv, zflxu1, zflxv1 )
189        CALL wrk_dealloc( jpi, jpj, zwdlmtu, zwdlmtv)
190        CALL wrk_dealloc( jpi, jpj, iwdflg )
191      !
192      END IF
193
194      IF( nn_timing == 1 )  CALL timing_stop('wad_lmt')
195   END SUBROUTINE wad_lmt
196   !!======================================================================
197END MODULE wadlmt
Note: See TracBrowser for help on using the repository browser.