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
RevLine 
[1271]1#define SPONGE_TOP
2
[5656]3MODULE agrif_top_sponge
[6140]4   !!======================================================================
5   !!                ***  MODULE agrif_top_sponge  ***
[9019]6   !! AGRIF :   sponge layer pakage for passive tracers (TOP)
[6140]7   !!======================================================================
8   !! History :  2.0  ! 2006-08  (R. Benshila, L. Debreu)  Original code
9   !!----------------------------------------------------------------------
[9019]10#if defined key_agrif && defined key_top
[6140]11   !!----------------------------------------------------------------------
12   !!   Agrif_Sponge_trc :
13   !!   interptrn_sponge : 
14   !!----------------------------------------------------------------------
[1271]15   USE par_oce
[5656]16   USE par_trc
[1271]17   USE oce
[6140]18   USE trc
[1271]19   USE dom_oce
20   USE agrif_oce
[9570]21   USE agrif_oce_sponge
[12377]22   USE vremap
[6140]23   !
24   USE in_out_manager
[2715]25   USE lib_mpp
[1271]26
27   IMPLICIT NONE
28   PRIVATE
29
[5656]30   PUBLIC Agrif_Sponge_trc, interptrn_sponge
[1271]31
[14148]32   !! * Substitutions
33#  include "domzgr_substitute.h90"
[1271]34   !!----------------------------------------------------------------------
[9598]35   !! NEMO/NST 4.0 , NEMO Consortium (2018)
[2528]36   !! $Id$
[10068]37   !! Software governed by the CeCILL license (see ./LICENSE)
[1271]38   !!----------------------------------------------------------------------
[5656]39CONTAINS
[1271]40
[5656]41   SUBROUTINE Agrif_Sponge_trc
[6140]42      !!----------------------------------------------------------------------
43      !!                   *** ROUTINE Agrif_Sponge_Trc ***
44      !!----------------------------------------------------------------------
[9019]45      REAL(wp) ::   zcoef   ! local scalar
[6140]46      !!----------------------------------------------------------------------
47      !
[1271]48#if defined SPONGE_TOP
[14086]49      !! Assume persistence:
[9056]50      zcoef = REAL(Agrif_rhot()-1,wp)/REAL(Agrif_rhot())
[14086]51
[6140]52      Agrif_SpecialValue    = 0._wp
[1271]53      Agrif_UseSpecialValue = .TRUE.
[14086]54      l_vremap              = ln_vert_remap
[6140]55      tabspongedone_trn     = .FALSE.
[14086]56      !
[9019]57      CALL Agrif_Bc_Variable( trn_sponge_id, calledweight=zcoef, procname=interptrn_sponge )
[14086]58      !
[1271]59      Agrif_UseSpecialValue = .FALSE.
[14086]60      l_vremap              = .FALSE.
[1271]61#endif
[6140]62      !
[1271]63   END SUBROUTINE Agrif_Sponge_Trc
64
[5656]65
[14086]66   SUBROUTINE interptrn_sponge( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before) 
[6140]67      !!----------------------------------------------------------------------
[14086]68      !!                 *** ROUTINE interptrn_sponge ***
[6140]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      !
[5656]74      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
[14086]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
[9031]79      ! vertical interpolation:
[14086]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
[9031]84      INTEGER :: N_in, N_out
[6140]85      !!----------------------------------------------------------------------
[3680]86      !
[6140]87      IF( before ) THEN
[9031]88         DO jn = 1, jptra
89            DO jk=k1,k2
90               DO jj=j1,j2
91                  DO ji=i1,i2
[12377]92                     tabres(ji,jj,jk,jn) = tr(ji,jj,jk,jn,Kbb_a)
[9031]93                  END DO
94               END DO
95            END DO
96         END DO
97
[14086]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.
[9031]102            DO jj=j1,j2
103               DO ji=i1,i2
[14086]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
[9031]109               END DO
110            END DO
[14086]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):
[14170]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
[14086]137                 
[14170]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
[14086]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
[9031]171               END DO
[14086]172            END DO
[9031]173
[14086]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
[9031]181
[14086]182         ELSE
[12377]183
[14086]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           
[5656]222            DO jk = 1, jpkm1
[14086]223               ztu(i1-1:i2,j1-1:j2,jk) = 0._wp
224               DO jj = j1,j2
[5656]225                  DO ji = i1,i2-1
[14086]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) ) 
[6140]228                  END DO
229               END DO
[14086]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
[6140]237               !
[14086]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
[5656]263                     ENDIF
[6140]264                  END DO
265               END DO
[14086]266
[6140]267            END DO
268            !
269         END DO
270         !
[14086]271         tabspongedone_trn(i1:i2,j1:j2) = .TRUE.
272         !
[5656]273      ENDIF
[14086]274      !
[5656]275   END SUBROUTINE interptrn_sponge
276
[1271]277#else
[9019]278   !!----------------------------------------------------------------------
279   !!   Empty module                                           no TOP AGRIF
280   !!----------------------------------------------------------------------
[1271]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
[6140]287   !!======================================================================
[1271]288END MODULE agrif_top_sponge
Note: See TracBrowser for help on using the repository browser.