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_lim3_interp.F90 in branches/2017/dev_r8127_AGRIF_LIM3_GHOST/NEMOGCM/NEMO/NST_SRC – NEMO

source: branches/2017/dev_r8127_AGRIF_LIM3_GHOST/NEMOGCM/NEMO/NST_SRC/agrif_lim3_interp.F90 @ 8129

Last change on this file since 8129 was 8129, checked in by clem, 8 years ago

make those things work: ghostcells>1 + nn_ice(child)=0 + fix timing
File size: 16.5 KB
Line 
1MODULE agrif_lim3_interp
2   !!=====================================================================================
3   !!                       ***  MODULE agrif_lim3_interp ***
4   !! Nesting module :  interp surface ice boundary condition from a parent grid
5   !! Sea-Ice model  :  LIM 3.6 Sea ice model time-stepping
6   !!=====================================================================================
7   !! History :  2.0   !  04-2008  (F. Dupont)  initial version
8   !!            3.4   !  09-2012  (R. Benshila, C. Herbaut) update and EVP
9   !!            3.6   !  05-2016  (C. Rousset)  Add LIM3 compatibility
10   !!----------------------------------------------------------------------
11#if defined key_agrif && defined key_lim3 
12   !!----------------------------------------------------------------------
13   !!   'key_lim3'  :                                 LIM 3.6 sea-ice model
14   !!   'key_agrif' :                                 AGRIF library
15   !!----------------------------------------------------------------------
16   !!  agrif_interp_lim3    : interpolation of ice at "after" sea-ice time step
17   !!  agrif_interp_u_ice   : atomic routine to interpolate u_ice
18   !!  agrif_interp_v_ice   : atomic routine to interpolate v_ice
19   !!  agrif_interp_tra_ice : atomic routine to interpolate ice properties
20   !!----------------------------------------------------------------------
21   USE par_oce
22   USE dom_oce
23   USE sbc_oce
24   USE ice
25   USE agrif_ice
26   
27   IMPLICIT NONE
28   PRIVATE
29
30   PUBLIC agrif_interp_lim3
31
32   !!----------------------------------------------------------------------
33   !! NEMO/NST 3.6 , NEMO Consortium (2016)
34   !! $Id: agrif_lim3_interp.F90 6204 2016-01-04 13:47:06Z cetlod $
35   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
36   !!----------------------------------------------------------------------
37
38CONTAINS
39
40   SUBROUTINE agrif_interp_lim3( cd_type, kiter, kitermax )
41      !!-----------------------------------------------------------------------
42      !!                 *** ROUTINE agrif_rhg_lim3  ***
43      !!
44      !!  ** Method  : simple call to atomic routines using stored values to
45      !!  fill the boundaries depending of the position of the point and
46      !!  computing factor for time interpolation
47      !!-----------------------------------------------------------------------
48      CHARACTER(len=1), INTENT( in )           :: cd_type
49      INTEGER         , INTENT( in ), OPTIONAL :: kiter, kitermax
50      !!
51      REAL(wp) :: zbeta
52      !!-----------------------------------------------------------------------
53      !
54      IF( Agrif_Root() .OR. nn_ice==0 )  RETURN   ! clem2017: do not interpolate if inside Parent domain or if child domain does not have ice
55      !
56      SELECT CASE(cd_type)
57      CASE('U','V')
58         IF( PRESENT( kiter ) ) THEN  ! interpolation at the child sub-time step (only for ice rheology)
59            zbeta = ( REAL(lim_nbstep) - REAL(kitermax - kiter) / REAL(kitermax) ) /  &
60               &    ( Agrif_Rhot() * REAL(Agrif_Parent(nn_fsbc)) / REAL(nn_fsbc) )
61         ELSE                         ! interpolation at the child time step
62            zbeta = REAL(lim_nbstep) / ( Agrif_Rhot() * REAL(Agrif_Parent(nn_fsbc)) / REAL(nn_fsbc) )
63         ENDIF
64      CASE('T')
65            zbeta = REAL(lim_nbstep-1) / ( Agrif_Rhot() * REAL(Agrif_Parent(nn_fsbc)) / REAL(nn_fsbc) )
66      END SELECT
67      !
68      Agrif_SpecialValue=-9999.
69      Agrif_UseSpecialValue = .TRUE.
70      SELECT CASE(cd_type)
71      CASE('U')
72         CALL Agrif_Bc_variable( u_ice_id  , procname=interp_u_ice  , calledweight=zbeta )
73      CASE('V')
74         CALL Agrif_Bc_variable( v_ice_id  , procname=interp_v_ice  , calledweight=zbeta )
75      CASE('T')
76         CALL Agrif_Bc_variable( tra_ice_id, procname=interp_tra_ice, calledweight=zbeta )
77      END SELECT
78      Agrif_SpecialValue=0.
79      Agrif_UseSpecialValue = .FALSE.
80      !
81   END SUBROUTINE agrif_interp_lim3
82
83   !!------------------
84   !! Local subroutines
85   !!------------------
86   SUBROUTINE interp_u_ice( ptab, i1, i2, j1, j2, before )
87      !!-----------------------------------------------------------------------
88      !!                     *** ROUTINE interp_u_ice ***
89      !!
90      !! i1 i2 j1 j2 are the index of the boundaries parent(when before) and child (when after)
91      !! To solve issues when parent grid is "land" masked but not all the corresponding child grid points,
92      !! put -9999 WHERE the parent grid is masked. The child solution will be found in the 9(?) points around
93      !!-----------------------------------------------------------------------
94      INTEGER , INTENT(in) :: i1, i2, j1, j2
95      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab
96      LOGICAL , INTENT(in) :: before
97      !!
98      REAL(wp) :: zrhoy
99      !!-----------------------------------------------------------------------
100      !
101      IF( before ) THEN  ! parent grid
102         ptab(:,:) = e2u(i1:i2,j1:j2) * u_ice_b(i1:i2,j1:j2)
103         WHERE( umask(i1:i2,j1:j2,1) == 0. )  ptab(:,:) = -9999.
104      ELSE               ! child grid
105         zrhoy = Agrif_Rhoy()
106         u_ice(i1:i2,j1:j2) = ptab(:,:) / ( e2u(i1:i2,j1:j2) * zrhoy ) * umask(i1:i2,j1:j2,1)
107      ENDIF
108      !
109   END SUBROUTINE interp_u_ice
110
111
112   SUBROUTINE interp_v_ice( ptab, i1, i2, j1, j2, before )
113      !!-----------------------------------------------------------------------
114      !!                    *** ROUTINE interp_v_ice ***
115      !!
116      !! i1 i2 j1 j2 are the index of the boundaries parent(when before) and child (when after)
117      !! To solve issues when parent grid is "land" masked but not all the corresponding child grid points,
118      !! put -9999 WHERE the parent grid is masked. The child solution will be found in the 9(?) points around
119      !!-----------------------------------------------------------------------     
120      INTEGER , INTENT(in) :: i1, i2, j1, j2
121      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab
122      LOGICAL , INTENT(in) :: before
123      !!
124      REAL(wp) :: zrhox
125      !!-----------------------------------------------------------------------
126      !
127      IF( before ) THEN  ! parent grid
128         ptab(:,:) = e1v(i1:i2,j1:j2) * v_ice_b(i1:i2,j1:j2)
129         WHERE( vmask(i1:i2,j1:j2,1) == 0. )  ptab(:,:) = -9999.
130      ELSE               ! child grid
131         zrhox = Agrif_Rhox()
132         v_ice(i1:i2,j1:j2) = ptab(:,:) / ( e1v(i1:i2,j1:j2) * zrhox ) * vmask(i1:i2,j1:j2,1)
133      ENDIF
134      !
135   END SUBROUTINE interp_v_ice
136
137
138   SUBROUTINE interp_tra_ice( ptab, i1, i2, j1, j2, k1, k2, before, nb, ndir )
139      !!-----------------------------------------------------------------------
140      !!                    *** ROUTINE interp_tra_ice ***                           
141      !!
142      !! i1 i2 j1 j2 are the index of the boundaries parent(when before) and child (when after)
143      !! To solve issues when parent grid is "land" masked but not all the corresponding child grid points,
144      !! put -9999 WHERE the parent grid is masked. The child solution will be found in the 9(?) points around
145      !!-----------------------------------------------------------------------
146      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
147      INTEGER , INTENT(in) :: i1, i2, j1, j2, k1, k2
148      LOGICAL , INTENT(in) :: before
149      INTEGER , INTENT(in) :: nb, ndir
150      !!
151      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztab
152      INTEGER  ::   ji, jj, jk, jl, jm
153      INTEGER  ::   imin, imax, jmin, jmax
154      REAL(wp) ::   zrhox, z1, z2, z3, z4, z5, z6, z7
155      LOGICAL  ::   western_side, eastern_side, northern_side, southern_side
156      INTEGER  ::   ind1, ind2, ind3
157
158      !!-----------------------------------------------------------------------
159      ! tracers are not multiplied by grid cell here => before: * e12t ; after: * r1_e12t / rhox / rhoy
160      ! and it is ok since we conserve tracers (same as in the ocean).
161      ALLOCATE( ztab(SIZE(a_i_b,1),SIZE(a_i_b,2),SIZE(ptab,3)) )
162     
163      IF( before ) THEN  ! parent grid
164         jm = 1
165         DO jl = 1, jpl
166            ptab(i1:i2,j1:j2,jm) = a_i_b  (i1:i2,j1:j2,jl) ; jm = jm + 1
167            ptab(i1:i2,j1:j2,jm) = v_i_b  (i1:i2,j1:j2,jl) ; jm = jm + 1
168            ptab(i1:i2,j1:j2,jm) = v_s_b  (i1:i2,j1:j2,jl) ; jm = jm + 1
169            ptab(i1:i2,j1:j2,jm) = smv_i_b(i1:i2,j1:j2,jl) ; jm = jm + 1
170            ptab(i1:i2,j1:j2,jm) = oa_i_b (i1:i2,j1:j2,jl) ; jm = jm + 1
171            DO jk = 1, nlay_s
172               ptab(i1:i2,j1:j2,jm) = e_s_b(i1:i2,j1:j2,jk,jl) ; jm = jm + 1
173            ENDDO
174            DO jk = 1, nlay_i
175               ptab(i1:i2,j1:j2,jm) = e_i_b(i1:i2,j1:j2,jk,jl) ; jm = jm + 1
176            ENDDO
177         ENDDO
178         
179         DO jk = k1, k2
180            WHERE( tmask(i1:i2,j1:j2,1) == 0. )  ptab(i1:i2,j1:j2,jk) = -9999.
181         ENDDO
182         
183      ELSE               ! child grid
184!! ==> The easiest interpolation is the following commented lines
185         jm = 1
186         DO jl = 1, jpl
187            a_i  (i1:i2,j1:j2,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
188            v_i  (i1:i2,j1:j2,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
189            v_s  (i1:i2,j1:j2,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
190            smv_i(i1:i2,j1:j2,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
191            oa_i (i1:i2,j1:j2,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
192            DO jk = 1, nlay_s
193               e_s(i1:i2,j1:j2,jk,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
194            ENDDO
195            DO jk = 1, nlay_i
196               e_i(i1:i2,j1:j2,jk,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
197            ENDDO
198         ENDDO
199
200!! ==> this is a more complex interpolation since we mix solutions over a couple of grid points
201!!     it is advised to use it for fields modified by high order schemes (e.g. advection UM5...)
202!!        clem: for some reason (I don't know why), the following lines do not work
203!!              with mpp (or in realistic configurations?). It makes the model crash
204!         ! record ztab
205!         jm = 1
206!         DO jl = 1, jpl
207!            ztab(:,:,jm) = a_i  (:,:,jl) ; jm = jm + 1
208!            ztab(:,:,jm) = v_i  (:,:,jl) ; jm = jm + 1
209!            ztab(:,:,jm) = v_s  (:,:,jl) ; jm = jm + 1
210!            ztab(:,:,jm) = smv_i(:,:,jl) ; jm = jm + 1
211!            ztab(:,:,jm) = oa_i (:,:,jl) ; jm = jm + 1
212!            DO jk = 1, nlay_s
213!               ztab(:,:,jm) = e_s(:,:,jk,jl) ; jm = jm + 1
214!            ENDDO
215!            DO jk = 1, nlay_i
216!               ztab(:,:,jm) = e_i(:,:,jk,jl) ; jm = jm + 1
217!            ENDDO
218!         ENDDO
219!         !
220!         ! borders of the domain
221!         western_side  = (nb == 1).AND.(ndir == 1)  ;  eastern_side  = (nb == 1).AND.(ndir == 2)
222!         southern_side = (nb == 2).AND.(ndir == 1)  ;  northern_side = (nb == 2).AND.(ndir == 2)
223!         !
224!         ! spatial smoothing
225!         zrhox = Agrif_Rhox()
226!         z1 =      ( zrhox - 1. ) * 0.5
227!         z3 =      ( zrhox - 1. ) / ( zrhox + 1. )
228!         z6 = 2. * ( zrhox - 1. ) / ( zrhox + 1. )
229!         z7 =    - ( zrhox - 1. ) / ( zrhox + 3. )
230!         z2 = 1. - z1
231!         z4 = 1. - z3
232!         z5 = 1. - z6 - z7
233!         !
234!         ! Remove corners
235!         imin = i1  ;  imax = i2  ;  jmin = j1  ;  jmax = j2
236!         !!clem2017 ghost
237!         ind1 =     nbghostcells
238!         ind2 = 1 + nbghostcells
239!         ind3 = 2 + nbghostcells
240!         IF( (nbondj == -1) .OR. (nbondj == 2) )   jmin = ind3
241!         IF( (nbondj == +1) .OR. (nbondj == 2) )   jmax = nlcj-ind2
242!         IF( (nbondi == -1) .OR. (nbondi == 2) )   imin = ind3
243!         IF( (nbondi == +1) .OR. (nbondi == 2) )   imax = nlci-ind2
244!         !!clem2017 ghost
245!
246!         ! smoothed fields
247!         IF( eastern_side ) THEN
248!            ztab(nlci,j1:j2,:) = z1 * ptab(nlci,j1:j2,:) + z2 * ptab(nlci-1,j1:j2,:)
249!            DO jj = jmin, jmax
250!               rswitch = 0.
251!               IF( u_ice(nlci-2,jj) > 0._wp ) rswitch = 1.
252!               ztab(nlci-1,jj,:) = ( 1. - umask(nlci-2,jj,1) ) * ztab(nlci,jj,:)  &
253!                  &                +      umask(nlci-2,jj,1)   *  &
254!                  &                ( ( 1. - rswitch ) * ( z4 * ztab(nlci,jj,:)   + z3 * ztab(nlci-2,jj,:) )  &
255!                  &                  +      rswitch   * ( z6 * ztab(nlci-2,jj,:) + z5 * ztab(nlci,jj,:) + z7 * ztab(nlci-3,jj,:) ) )
256!               ztab(nlci-1,jj,:) = ztab(nlci-1,jj,:) * tmask(nlci-1,jj,1)
257!            END DO
258!         ENDIF
259!         !
260!         IF( northern_side ) THEN
261!            ztab(i1:i2,nlcj,:) = z1 * ptab(i1:i2,nlcj,:) + z2 * ptab(i1:i2,nlcj-1,:)
262!            DO ji = imin, imax
263!               rswitch = 0.
264!               IF( v_ice(ji,nlcj-2) > 0._wp ) rswitch = 1.
265!               ztab(ji,nlcj-1,:) = ( 1. - vmask(ji,nlcj-2,1) ) * ztab(ji,nlcj,:)  &
266!                  &                +      vmask(ji,nlcj-2,1)   *  &
267!                  &                ( ( 1. - rswitch ) * ( z4 * ztab(ji,nlcj,:)   + z3 * ztab(ji,nlcj-2,:) ) &
268!                  &                  +      rswitch   * ( z6 * ztab(ji,nlcj-2,:) + z5 * ztab(ji,nlcj,:) + z7 * ztab(ji,nlcj-3,:) ) )
269!               ztab(ji,nlcj-1,:) = ztab(ji,nlcj-1,:) * tmask(ji,nlcj-1,1)
270!            END DO
271!         END IF
272!         !
273!         IF( western_side) THEN
274!            ztab(1,j1:j2,:) = z1 * ptab(1,j1:j2,:) + z2 * ptab(2,j1:j2,:)
275!            DO jj = jmin, jmax
276!               rswitch = 0.
277!               IF( u_ice(2,jj) < 0._wp ) rswitch = 1.
278!               ztab(2,jj,:) = ( 1. - umask(2,jj,1) ) * ztab(1,jj,:)  &
279!                  &           +      umask(2,jj,1)   *   &
280!                  &           ( ( 1. - rswitch ) * ( z4 * ztab(1,jj,:) + z3 * ztab(3,jj,:) ) &
281!                  &             +      rswitch   * ( z6 * ztab(3,jj,:) + z5 * ztab(1,jj,:) + z7 * ztab(4,jj,:) ) )
282!               ztab(2,jj,:) = ztab(2,jj,:) * tmask(2,jj,1)
283!            END DO
284!         ENDIF
285!         !
286!         IF( southern_side ) THEN
287!            ztab(i1:i2,1,:) = z1 * ptab(i1:i2,1,:) + z2 * ptab(i1:i2,2,:)
288!            DO ji = imin, imax
289!               rswitch = 0.
290!               IF( v_ice(ji,2) < 0._wp ) rswitch = 1.
291!               ztab(ji,2,:) = ( 1. - vmask(ji,2,1) ) * ztab(ji,1,:)  &
292!                  &           +      vmask(ji,2,1)   *  &
293!                  &           ( ( 1. - rswitch ) * ( z4 * ztab(ji,1,:) + z3 * ztab(ji,3,:) ) &
294!                  &             +      rswitch   * ( z6 * ztab(ji,3,:) + z5 * ztab(ji,1,:) + z7 * ztab(ji,4,:) ) )
295!               ztab(ji,2,:) = ztab(ji,2,:) * tmask(ji,2,1)
296!            END DO
297!         END IF
298!         !
299!         ! Treatment of corners
300!         IF( (eastern_side) .AND. ((nbondj == -1).OR.(nbondj == 2)) )  ztab(nlci-1,2,:)      = ptab(nlci-1,2,:)      ! East south
301!         IF( (eastern_side) .AND. ((nbondj ==  1).OR.(nbondj == 2)) )  ztab(nlci-1,nlcj-1,:) = ptab(nlci-1,nlcj-1,:) ! East north
302!         IF( (western_side) .AND. ((nbondj == -1).OR.(nbondj == 2)) )  ztab(2,2,:)           = ptab(2,2,:)           ! West south
303!         IF( (western_side) .AND. ((nbondj ==  1).OR.(nbondj == 2)) )  ztab(2,nlcj-1,:)      = ptab(2,nlcj-1,:)      ! West north
304!
305!         ! retrieve ice tracers
306!         jm = 1
307!         DO jl = 1, jpl
308!            a_i  (i1:i2,j1:j2,jl) = ztab(i1:i2,j1:j2,jm) ; jm = jm + 1
309!            v_i  (i1:i2,j1:j2,jl) = ztab(i1:i2,j1:j2,jm) ; jm = jm + 1
310!            v_s  (i1:i2,j1:j2,jl) = ztab(i1:i2,j1:j2,jm) ; jm = jm + 1
311!            smv_i(i1:i2,j1:j2,jl) = ztab(i1:i2,j1:j2,jm) ; jm = jm + 1
312!            oa_i (i1:i2,j1:j2,jl) = ztab(i1:i2,j1:j2,jm) ; jm = jm + 1
313!            DO jk = 1, nlay_s
314!               e_s(i1:i2,j1:j2,jk,jl) = ztab(i1:i2,j1:j2,jm) ; jm = jm + 1
315!            ENDDO
316!            DO jk = 1, nlay_i
317!               e_i(i1:i2,j1:j2,jk,jl) = ztab(i1:i2,j1:j2,jm) ; jm = jm + 1
318!            ENDDO
319!         ENDDO
320       
321         ! integrated values
322         vt_i (i1:i2,j1:j2) = SUM( v_i(i1:i2,j1:j2,:), dim=3 )
323         vt_s (i1:i2,j1:j2) = SUM( v_s(i1:i2,j1:j2,:), dim=3 )
324         at_i (i1:i2,j1:j2) = SUM( a_i(i1:i2,j1:j2,:), dim=3 )
325         et_s(i1:i2,j1:j2)  = SUM( SUM( e_s(i1:i2,j1:j2,:,:), dim=4 ), dim=3 )
326         et_i(i1:i2,j1:j2)  = SUM( SUM( e_i(i1:i2,j1:j2,:,:), dim=4 ), dim=3 )
327
328      ENDIF
329     
330      DEALLOCATE( ztab )
331      !
332   END SUBROUTINE interp_tra_ice
333
334#else
335CONTAINS
336   SUBROUTINE agrif_lim3_interp_empty
337      !!---------------------------------------------
338      !!   *** ROUTINE agrif_lim3_interp_empty ***
339      !!---------------------------------------------
340      WRITE(*,*)  'agrif_lim3_interp : You should not have seen this print! error?'
341   END SUBROUTINE agrif_lim3_interp_empty
342#endif
343END MODULE agrif_lim3_interp
Note: See TracBrowser for help on using the repository browser.