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 @ 14476

Last change on this file since 14476 was 14170, checked in by jchanut, 3 years ago

#2222, 2129: 1) Corrected ssh initialization from parent in line with what has been introduced by Sibylle 2) Fixed bug in dyn interp with expliciit free surface 3) Added check on number of levels in child grid without vertical remapping (must be < jpk_parent) 4) Removed the constrain on initialization from parent only when starting from climatology (requires Euler first step though).
  • 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)
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.