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.
agrif_top_sponge.F90 in NEMO/trunk/src/NST – NEMO

source: NEMO/trunk/src/NST/agrif_top_sponge.F90 @ 15437

Last change on this file since 15437 was 15437, checked in by jchanut, 11 months ago

agrif fixes: i) Add masking prior interpolation on "before" tracers: these seem to be strangly unmasked which induces extrapolation issues. ii) Add flux limiter on tracer sponge (recall sponge is applied on anomalies) to prevent from generating spurious extrema.

  • Property svn:keywords set to Id
File size: 11.4 KB
Line 
1#define SPONGE_TOP
2
3MODULE agrif_top_sponge
4   !!======================================================================
5   !!                ***  MODULE agrif_top_sponge  ***
6   !! AGRIF :   sponge layer pakage for passive tracers (TOP)
7   !!======================================================================
8   !! History :  2.0  ! 2006-08  (R. Benshila, L. Debreu)  Original code
9   !!----------------------------------------------------------------------
10#if defined key_agrif && defined key_top
11   !!----------------------------------------------------------------------
12   !!   Agrif_Sponge_trc :
13   !!   interptrn_sponge : 
14   !!----------------------------------------------------------------------
15   USE par_oce
16   USE par_trc
17   USE oce
18   USE trc
19   USE dom_oce
20   USE agrif_oce
21   USE agrif_oce_sponge
22   USE vremap
23   !
24   USE in_out_manager
25   USE lib_mpp
26
27   IMPLICIT NONE
28   PRIVATE
29
30   PUBLIC Agrif_Sponge_trc, interptrn_sponge
31
32   !! * Substitutions
33#  include "domzgr_substitute.h90"
34   !!----------------------------------------------------------------------
35   !! NEMO/NST 4.0 , NEMO Consortium (2018)
36   !! $Id$
37   !! Software governed by the CeCILL license (see ./LICENSE)
38   !!----------------------------------------------------------------------
39CONTAINS
40
41   SUBROUTINE Agrif_Sponge_trc
42      !!----------------------------------------------------------------------
43      !!                   *** ROUTINE Agrif_Sponge_Trc ***
44      !!----------------------------------------------------------------------
45      REAL(wp) ::   zcoef   ! local scalar
46      !!----------------------------------------------------------------------
47      !
48#if defined SPONGE_TOP
49      !! Assume persistence:
50      zcoef = REAL(Agrif_rhot()-1,wp)/REAL(Agrif_rhot())
51
52      Agrif_SpecialValue    = 0._wp
53      Agrif_UseSpecialValue = .TRUE.
54      l_vremap              = ln_vert_remap
55      tabspongedone_trn     = .FALSE.
56      !
57      CALL Agrif_Bc_Variable( trn_sponge_id, calledweight=zcoef, procname=interptrn_sponge )
58      !
59      Agrif_UseSpecialValue = .FALSE.
60      l_vremap              = .FALSE.
61#endif
62      !
63   END SUBROUTINE Agrif_Sponge_Trc
64
65
66   SUBROUTINE interptrn_sponge( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before) 
67      !!----------------------------------------------------------------------
68      !!                 *** ROUTINE interptrn_sponge ***
69      !!----------------------------------------------------------------------
70      INTEGER                                     , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2, n1, n2
71      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) ::   tabres
72      LOGICAL                                     , INTENT(in   ) ::   before
73      !
74      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
75      INTEGER  ::   iku, ikv
76      REAL(wp) :: ztra, zabe1, zabe2, zbtr, zhtot
77      REAL(wp), DIMENSION(i1-1:i2,j1-1:j2,jpk) :: ztu, ztv
78      REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::trbdiff
79      ! vertical interpolation:
80      REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tabres_child
81      REAL(wp), DIMENSION(k1:k2,n1:n2-1) :: tabin, tabin_i
82      REAL(wp), DIMENSION(k1:k2) :: z_in, z_in_i, h_in_i
83      REAL(wp), DIMENSION(1:jpk) :: h_out, z_out
84      INTEGER :: N_in, N_out
85      !!----------------------------------------------------------------------
86      !
87      IF( before ) THEN
88         DO jn = 1, jptra
89            DO jk=k1,k2
90               DO jj=j1,j2
91                  DO ji=i1,i2
92                     tabres(ji,jj,jk,jn) = tr(ji,jj,jk,jn,Kbb_a) * tmask(ji,jj,jk) 
93                  END DO
94               END DO
95            END DO
96         END DO
97
98         IF ( l_vremap.OR.ln_zps ) THEN
99
100            ! Fill cell depths (i.e. gdept) to be interpolated
101            ! Warning: these are masked, hence extrapolated prior interpolation.
102            DO jj=j1,j2
103               DO ji=i1,i2
104                  tabres(ji,jj,k1,jptra+1) = 0.5_wp * tmask(ji,jj,k1) * e3t(ji,jj,k1,Kbb_a)
105                  DO jk=k1+1,k2
106                     tabres(ji,jj,jk,jptra+1) = tmask(ji,jj,jk) * &
107                        & ( tabres(ji,jj,jk-1,jptra+1) + 0.5_wp * (e3t(ji,jj,jk-1,Kbb_a)+e3t(ji,jj,jk,Kbb_a)) )
108                  END DO
109               END DO
110            END DO
111
112            ! Save ssh at last level:
113            IF ( .NOT.ln_linssh ) THEN
114               tabres(i1:i2,j1:j2,k2,jptra+1) = ssh(i1:i2,j1:j2,Kbb_a)*tmask(i1:i2,j1:j2,1) 
115            END IF 
116   
117         END IF
118
119      ELSE 
120         !
121         IF ( l_vremap ) THEN
122
123            IF (ln_linssh) tabres(i1:i2,j1:j2,k2,n2) = 0._wp
124
125            DO jj=j1,j2
126               DO ji=i1,i2
127
128                  tabres_child(ji,jj,:,:) = 0._wp 
129                  ! Build vertical grids:
130                  N_in = mbkt_parent(ji,jj)
131                  ! Input grid (account for partial cells if any):
132                  IF ( N_in > 0 ) THEN
133                     DO jk=1,N_in
134                        z_in(jk) = tabres(ji,jj,jk,n2) - tabres(ji,jj,k2,n2)
135                        tabin(jk,1:jptra) = tabres(ji,jj,jk,1:jptra)
136                     END DO
137                 
138                     ! Intermediate grid:
139                     DO jk = 1, N_in
140                        h_in_i(jk) = e3t0_parent(ji,jj,jk) * & 
141                          &       (1._wp + tabres(ji,jj,k2,n2)/(ht0_parent(ji,jj)*ssmask(ji,jj) + 1._wp - ssmask(ji,jj)))
142                     END DO
143                     z_in_i(1) = 0.5_wp * h_in_i(1)
144                     DO jk=2,N_in
145                        z_in_i(jk) = z_in_i(jk-1) + 0.5_wp * ( h_in_i(jk) + h_in_i(jk-1) )
146                     END DO
147                     z_in_i(1:N_in) = z_in_i(1:N_in)  - tabres(ji,jj,k2,n2)
148                  END IF
149                  ! Output (Child) grid:
150                  N_out = mbkt(ji,jj)
151                  DO jk=1,N_out
152                     h_out(jk) = e3t(ji,jj,jk,Kbb_a)
153                  END DO
154                  z_out(1) = 0.5_wp * h_out(1)
155                  DO jk=2,N_out
156                     z_out(jk) = z_out(jk-1) + 0.5_wp * ( h_out(jk)+h_out(jk-1) )
157                  END DO
158                  IF (.NOT.ln_linssh) z_out(1:N_out) = z_out(1:N_out)  - ssh(ji,jj,Kbb_a)
159
160                  ! Account for small differences in the free-surface
161                  IF ( sum(h_out(1:N_out)) > sum(h_in_i(1:N_in) )) THEN
162                     h_out(1) = h_out(1)  - ( sum(h_out(1:N_out))-sum(h_in_i(1:N_in)) )
163                  ELSE
164                     h_in_i(1)= h_in_i(1) - ( sum(h_in_i(1:N_in))-sum(h_out(1:N_out)) )
165                  END IF
166                  IF (N_in*N_out > 0) THEN
167                     CALL remap_linear(tabin(1:N_in,1:jptra),z_in(1:N_in),tabin_i(1:N_in,1:jptra),z_in_i(1:N_in),N_in,N_in,jptra)
168                     CALL reconstructandremap(tabin_i(1:N_in,1:jptra),h_in_i(1:N_in),tabres_child(ji,jj,1:N_out,1:jptra),h_out(1:N_out),N_in,N_out,jptra)
169!                     CALL remap_linear(tabin(1:N_in,1:jptra),z_in(1:N_in),tabres_child(ji,jj,1:N_out,1:jptra),z_out(1:N_in),N_in,N_out,jptra) 
170                  ENDIF
171               END DO
172            END DO
173
174            DO jj=j1,j2
175               DO ji=i1,i2
176                  DO jk=1,jpkm1
177                     trbdiff(ji,jj,jk,1:jptra) = (tr(ji,jj,jk,1:jptra,Kbb_a) - tabres_child(ji,jj,jk,1:jptra)) * tmask(ji,jj,jk)
178                  END DO
179               END DO
180            END DO
181
182         ELSE
183
184            IF ( Agrif_Parent(ln_zps) ) THEN ! Account for partial cells
185
186               DO jj=j1,j2
187                  DO ji=i1,i2
188                     !
189                     N_in  = mbkt(ji,jj) 
190                     N_out = mbkt(ji,jj) 
191                     z_in(1) = tabres(ji,jj,1,n2)
192                     tabin(1,1:jptra) = tabres(ji,jj,1,1:jptra)
193                     DO jk=2, N_in
194                        z_in(jk) = tabres(ji,jj,jk,n2)
195                        tabin(jk,1:jptra) = tabres(ji,jj,jk,1:jptra)
196                     END DO
197                     IF (.NOT.ln_linssh) z_in(1:N_in) = z_in(1:N_in) - tabres(ji,jj,k2,n2)
198
199                     z_out(1) = 0.5_wp * e3t(ji,jj,1,Kbb_a)
200                     DO jk=2, N_out
201                        z_out(jk) = z_out(jk-1) + 0.5_wp * (e3t(ji,jj,jk-1,Kbb_a) + e3t(ji,jj,jk,Kbb_a)) 
202                     END DO
203                     IF (.NOT.ln_linssh) z_out(1:N_out) = z_out(1:N_out) - ssh(ji,jj,Kbb_a)
204
205                     CALL remap_linear(tabin(1:N_in,1:jptra), z_in(1:N_in), tabres(ji,jj,1:N_out,1:jptra), &
206                                         &   z_out(1:N_out), N_in, N_out, jptra)
207                  END DO
208               END DO
209            ENDIF
210
211            DO jj=j1,j2
212               DO ji=i1,i2
213                  DO jk=1,jpkm1
214                     trbdiff(ji,jj,jk,1:jptra) = (tr(ji,jj,jk,1:jptra,Kbb_a) - tabres(ji,jj,jk,1:jptra))*tmask(ji,jj,jk)
215                  END DO
216               END DO
217            END DO
218
219         END IF
220
221         DO jn = 1, jptra           
222            DO jk = 1, jpkm1
223               ztu(i1-1:i2,j1-1:j2,jk) = 0._wp
224               DO jj = j1,j2
225                  DO ji = i1,i2-1
226                     zabe1 = rn_sponge_tra * r1_Dt * umask(ji,jj,jk) * e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a)
227                     ztu(ji,jj,jk) = zabe1 * fspu(ji,jj) * ( trbdiff(ji+1,jj  ,jk,jn) - trbdiff(ji,jj,jk,jn) ) 
228                  END DO
229               END DO
230               ztv(i1-1:i2,j1-1:j2,jk) = 0._wp
231               DO ji = i1,i2
232                  DO jj = j1,j2-1
233                     zabe2 = rn_sponge_tra * r1_Dt * vmask(ji,jj,jk) * e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm_a)
234                     ztv(ji,jj,jk) = zabe2 * fspv(ji,jj) * ( trbdiff(ji  ,jj+1,jk,jn) - trbdiff(ji,jj,jk,jn) )
235                  END DO
236               END DO
237               !
238               IF( ln_zps ) THEN      ! set gradient at partial step level
239                  DO jj = j1,j2
240                     DO ji = i1,i2
241                        ! last level
242                        iku = mbku(ji,jj)
243                        ikv = mbkv(ji,jj)
244                        IF( iku == jk )   ztu(ji,jj,jk) = 0._wp
245                        IF( ikv == jk )   ztv(ji,jj,jk) = 0._wp
246                     END DO
247                  END DO
248               ENDIF
249            END DO
250            !
251! JC: there is something wrong with the Laplacian in corners
252            DO jk = 1, jpkm1
253               DO jj = j1,j2
254                  DO ji = i1,i2
255                     IF (.NOT. tabspongedone_trn(ji,jj)) THEN
256                        zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm_a)
257                        ! horizontal diffusive trends
258                        ztra = zbtr * ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk)   & 
259                          &           + ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) &
260                          &   - rn_trelax_tra * r1_Dt * fspt(ji,jj) * trbdiff(ji,jj,jk,jn)
261                        ! add it to the general tracer trends
262                        tr(ji,jj,jk,jn,Krhs_a) = tr(ji,jj,jk,jn,Krhs_a) + ztra
263                     ENDIF
264                  END DO
265               END DO
266
267            END DO
268            !
269         END DO
270         !
271         tabspongedone_trn(i1:i2,j1:j2) = .TRUE.
272         !
273      ENDIF
274      !
275   END SUBROUTINE interptrn_sponge
276
277#else
278   !!----------------------------------------------------------------------
279   !!   Empty module                                           no TOP AGRIF
280   !!----------------------------------------------------------------------
281CONTAINS
282   SUBROUTINE agrif_top_sponge_empty
283      WRITE(*,*)  'agrif_top_sponge : You should not have seen this print! error?'
284   END SUBROUTINE agrif_top_sponge_empty
285#endif
286
287   !!======================================================================
288END MODULE agrif_top_sponge
Note: See TracBrowser for help on using the repository browser.