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_ice_interp.F90 in NEMO/branches/2020/dev_r12558_HPC-08_epico_Extra_Halo/src/NST – NEMO

source: NEMO/branches/2020/dev_r12558_HPC-08_epico_Extra_Halo/src/NST/agrif_ice_interp.F90 @ 12807

Last change on this file since 12807 was 12807, checked in by smasson, 5 years ago

Extra_Halo: input file only over inner domain + new variables names, see #2366
  • Property svn:keywords set to Id
File size: 18.1 KB
Line 
1MODULE agrif_ice_interp
2   !!=====================================================================================
3   !!                       ***  MODULE agrif_ice_interp ***
4   !! Nesting module :  interp surface ice boundary condition from a parent grid
5   !!=====================================================================================
6   !! History :  2.0   !  04-2008  (F. Dupont)               initial version
7   !!            3.4   !  09-2012  (R. Benshila, C. Herbaut) update and EVP
8   !!            4.0   !  2018     (C. Rousset)              SI3 compatibility
9   !!----------------------------------------------------------------------
10#if defined key_agrif && defined key_si3 
11   !!----------------------------------------------------------------------
12   !!   'key_si3'                                         SI3 sea-ice model
13   !!   'key_agrif'                                       AGRIF library
14   !!----------------------------------------------------------------------
15   !!  agrif_interp_ice    : interpolation of ice at "after" sea-ice time step
16   !!  agrif_interp_u_ice   : atomic routine to interpolate u_ice
17   !!  agrif_interp_v_ice   : atomic routine to interpolate v_ice
18   !!  agrif_interp_tra_ice : atomic routine to interpolate ice properties
19   !!----------------------------------------------------------------------
20   USE par_oce
21   USE dom_oce
22   USE sbc_oce
23   USE ice
24   USE agrif_ice
25   USE phycst , ONLY: rt0
26   
27   IMPLICIT NONE
28   PRIVATE
29
30   PUBLIC   agrif_interp_ice   ! called by agrif_user.F90
31
32   !!----------------------------------------------------------------------
33   !! NEMO/NST 4.0 , NEMO Consortium (2018)
34   !! $Id$
35   !! Software governed by the CeCILL license (see ./LICENSE)
36   !!----------------------------------------------------------------------
37
38CONTAINS
39
40   SUBROUTINE agrif_interp_ice( cd_type, kiter, kitermax )
41      !!-----------------------------------------------------------------------
42      !!                 *** ROUTINE agrif_interp_ice  ***
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   ! local scalar
52      !!-----------------------------------------------------------------------
53      !
54      IF( Agrif_Root() .OR. nn_ice==0 )  RETURN   ! do not interpolate if inside Parent Grid 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 ice sub-time step (only for ice rheology)
59            zbeta = ( REAL(nbstep_ice) - REAL(kitermax - kiter) / REAL(kitermax) ) /  &
60               &    ( Agrif_Rhot() * REAL(Agrif_Parent(nn_fsbc)) / REAL(nn_fsbc) )
61         ELSE                         ! interpolation at the child ice time step
62            zbeta = REAL(nbstep_ice) / ( Agrif_Rhot() * REAL(Agrif_Parent(nn_fsbc)) / REAL(nn_fsbc) )
63         ENDIF
64      CASE('T')
65            zbeta = REAL(nbstep_ice) / ( 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')   ;   CALL Agrif_Bc_variable( u_ice_id  , procname=interp_u_ice  , calledweight=zbeta )
72      CASE('V')   ;   CALL Agrif_Bc_variable( v_ice_id  , procname=interp_v_ice  , calledweight=zbeta )
73      CASE('T')   ;   CALL Agrif_Bc_variable( tra_ice_id, procname=interp_tra_ice, calledweight=zbeta )
74      END SELECT
75      Agrif_SpecialValue    = 0._wp
76      Agrif_UseSpecialValue = .FALSE.
77      !
78   END SUBROUTINE agrif_interp_ice
79
80
81   SUBROUTINE interp_u_ice( ptab, i1, i2, j1, j2, before )
82      !!-----------------------------------------------------------------------
83      !!                     *** ROUTINE interp_u_ice ***
84      !!
85      !! i1 i2 j1 j2 are the index of the boundaries parent(when before) and child (when after)
86      !! To solve issues when parent grid is "land" masked but not all the corresponding child
87      !! grid points, put Agrif_SpecialValue WHERE the parent grid is masked.
88      !! The child solution will be found in the 9(?) points around
89      !!-----------------------------------------------------------------------
90      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
91      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
92      LOGICAL                         , INTENT(in   ) ::   before
93      !!
94      REAL(wp) ::   zrhoy   ! local scalar
95      !!-----------------------------------------------------------------------
96      !
97      IF( before ) THEN  ! parent grid
98         ptab(:,:) = e2u(i1:i2,j1:j2) * u_ice(i1:i2,j1:j2)
99         WHERE( umask(i1:i2,j1:j2,1) == 0. )   ptab(i1:i2,j1:j2) = Agrif_SpecialValue
100      ELSE               ! child grid
101         zrhoy = Agrif_Rhoy()
102         u_ice(i1:i2,j1:j2) = ptab(i1:i2,j1:j2) / ( e2u(i1:i2,j1:j2) * zrhoy ) * umask(i1:i2,j1:j2,1)
103      ENDIF
104      !
105   END SUBROUTINE interp_u_ice
106
107
108   SUBROUTINE interp_v_ice( ptab, i1, i2, j1, j2, before )
109      !!-----------------------------------------------------------------------
110      !!                    *** ROUTINE interp_v_ice ***
111      !!
112      !! i1 i2 j1 j2 are the index of the boundaries parent(when before) and child (when after)
113      !! To solve issues when parent grid is "land" masked but not all the corresponding child
114      !! grid points, put Agrif_SpecialValue WHERE the parent grid is masked.
115      !! The child solution will be found in the 9(?) points around
116      !!-----------------------------------------------------------------------     
117      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
118      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
119      LOGICAL                         , INTENT(in   ) ::   before
120      !!
121      REAL(wp) ::   zrhox   ! local scalar
122      !!-----------------------------------------------------------------------
123      !
124      IF( before ) THEN  ! parent grid
125         ptab(:,:) = e1v(i1:i2,j1:j2) * v_ice(i1:i2,j1:j2)
126         WHERE( vmask(i1:i2,j1:j2,1) == 0. )   ptab(i1:i2,j1:j2) = Agrif_SpecialValue
127      ELSE               ! child grid
128         zrhox = Agrif_Rhox()
129         v_ice(i1:i2,j1:j2) = ptab(i1:i2,j1:j2) / ( e1v(i1:i2,j1:j2) * zrhox ) * vmask(i1:i2,j1:j2,1)
130      ENDIF
131      !
132   END SUBROUTINE interp_v_ice
133
134
135   SUBROUTINE interp_tra_ice( ptab, i1, i2, j1, j2, k1, k2, before, nb, ndir )
136      !!-----------------------------------------------------------------------
137      !!                    *** ROUTINE interp_tra_ice ***                           
138      !!
139      !! i1 i2 j1 j2 are the index of the boundaries parent(when before) and child (when after)
140      !! To solve issues when parent grid is "land" masked but not all the corresponding child
141      !! grid points, put Agrif_SpecialValue WHERE the parent grid is masked.
142      !! The child solution will be found in the 9(?) points around
143      !!-----------------------------------------------------------------------
144      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) ::   ptab
145      INTEGER                               , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2
146      LOGICAL                               , INTENT(in   ) ::   before
147      INTEGER                               , INTENT(in   ) ::   nb, ndir
148      !!
149      INTEGER  ::   ji, jj, jk, jl, jm
150      INTEGER  ::   imin, imax, jmin, jmax
151      LOGICAL  ::   western_side, eastern_side, northern_side, southern_side
152      REAL(wp) ::   zrhox, z1, z2, z3, z4, z5, z6, z7
153      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztab
154      !!-----------------------------------------------------------------------
155      ! tracers are not multiplied by grid cell here => before: * e1e2t ; after: * r1_e1e2t / rhox / rhoy
156      ! and it is ok since we conserve tracers (same as in the ocean).
157      ALLOCATE( ztab(SIZE(a_i,1),SIZE(a_i,2),SIZE(ptab,3)) )
158     
159      IF( before ) THEN  ! parent grid
160         jm = 1
161         DO jl = 1, jpl
162            ptab(i1:i2,j1:j2,jm  ) = a_i (i1:i2,j1:j2,jl)
163            ptab(i1:i2,j1:j2,jm+1) = v_i (i1:i2,j1:j2,jl)
164            ptab(i1:i2,j1:j2,jm+2) = v_s (i1:i2,j1:j2,jl)
165            ptab(i1:i2,j1:j2,jm+3) = sv_i(i1:i2,j1:j2,jl)
166            ptab(i1:i2,j1:j2,jm+4) = oa_i(i1:i2,j1:j2,jl)
167            ptab(i1:i2,j1:j2,jm+5) = a_ip(i1:i2,j1:j2,jl)
168            ptab(i1:i2,j1:j2,jm+6) = v_ip(i1:i2,j1:j2,jl)
169            ptab(i1:i2,j1:j2,jm+7) = t_su(i1:i2,j1:j2,jl)
170            jm = jm + 8
171            DO jk = 1, nlay_s
172               ptab(i1:i2,j1:j2,jm) = e_s(i1:i2,j1:j2,jk,jl)   ;   jm = jm + 1
173            END DO
174            DO jk = 1, nlay_i
175               ptab(i1:i2,j1:j2,jm) = e_i(i1:i2,j1:j2,jk,jl)   ;   jm = jm + 1
176            END DO
177         END DO
178         
179         DO jk = k1, k2
180            WHERE( tmask(i1:i2,j1:j2,1) == 0._wp )   ptab(i1:i2,j1:j2,jk) = Agrif_SpecialValue
181         END DO
182         !
183      ELSE               ! child grid
184         !
185!         IF( nbghostcells > 1 ) THEN   ! ==> The easiest interpolation is used
186            !
187            jm = 1
188            DO jl = 1, jpl
189               !
190               DO jj = j1, j2
191                  DO ji = i1, i2
192                     a_i (ji,jj,jl) = ptab(ji,jj,jm  ) * tmask(ji,jj,1)
193                     v_i (ji,jj,jl) = ptab(ji,jj,jm+1) * tmask(ji,jj,1)
194                     v_s (ji,jj,jl) = ptab(ji,jj,jm+2) * tmask(ji,jj,1)
195                     sv_i(ji,jj,jl) = ptab(ji,jj,jm+3) * tmask(ji,jj,1)
196                     oa_i(ji,jj,jl) = ptab(ji,jj,jm+4) * tmask(ji,jj,1)
197                     a_ip(ji,jj,jl) = ptab(ji,jj,jm+5) * tmask(ji,jj,1)
198                     v_ip(ji,jj,jl) = ptab(ji,jj,jm+6) * tmask(ji,jj,1)
199                     t_su(ji,jj,jl) = ptab(ji,jj,jm+7) * tmask(ji,jj,1)
200                  END DO
201               END DO
202               jm = jm + 8
203               !
204               DO jk = 1, nlay_s
205                  e_s(i1:i2,j1:j2,jk,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1)
206                  jm = jm + 1
207               END DO
208               !
209               DO jk = 1, nlay_i
210                  e_i(i1:i2,j1:j2,jk,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1)
211                  jm = jm + 1
212               END DO
213               !
214            END DO
215            !
216!!==> clem: this interpolation does not work because it creates negative values, due
217!!          to negative coefficients when mixing points (for ex. z7)
218!!
219!         ELSE                          ! ==> complex interpolation (only one ghost cell available)
220!            !! Use a more complex interpolation since we mix solutions over a couple of grid points
221!            !! it is advised to use it for fields modified by high order schemes (e.g. advection UM5...)
222!            ! record ztab
223!            jm = 1
224!            DO jl = 1, jpl
225!               ztab(:,:,jm  ) = a_i (:,:,jl)
226!               ztab(:,:,jm+1) = v_i (:,:,jl)
227!               ztab(:,:,jm+2) = v_s (:,:,jl)
228!               ztab(:,:,jm+3) = sv_i(:,:,jl)
229!               ztab(:,:,jm+4) = oa_i(:,:,jl)
230!               ztab(:,:,jm+5) = a_ip(:,:,jl)
231!               ztab(:,:,jm+6) = v_ip(:,:,jl)
232!               ztab(:,:,jm+7) = t_su(:,:,jl)
233!               jm = jm + 8
234!               DO jk = 1, nlay_s
235!                  ztab(:,:,jm) = e_s(:,:,jk,jl)
236!                  jm = jm + 1
237!               END DO
238!               DO jk = 1, nlay_i
239!                  ztab(:,:,jm) = e_i(:,:,jk,jl)
240!                  jm = jm + 1
241!               END DO
242!               !
243!            END DO
244!            !
245!            ! borders of the domain
246!            western_side  = (nb == 1).AND.(ndir == 1)  ;  eastern_side  = (nb == 1).AND.(ndir == 2)
247!            southern_side = (nb == 2).AND.(ndir == 1)  ;  northern_side = (nb == 2).AND.(ndir == 2)
248!            !
249!            ! spatial smoothing
250!            zrhox = Agrif_Rhox()
251!            z1 =      ( zrhox - 1. ) * 0.5
252!            z3 =      ( zrhox - 1. ) / ( zrhox + 1. )
253!            z6 = 2. * ( zrhox - 1. ) / ( zrhox + 1. )
254!            z7 =    - ( zrhox - 1. ) / ( zrhox + 3. )
255!            z2 = 1. - z1
256!            z4 = 1. - z3
257!            z5 = 1. - z6 - z7
258!            !
259!            ! Remove corners
260!            imin = i1  ;  imax = i2  ;  jmin = j1  ;  jmax = j2
261!            IF( (nbondj == -1) .OR. (nbondj == 2) )   jmin = 3
262!            IF( (nbondj == +1) .OR. (nbondj == 2) )   jmax = jpj-2
263!            IF( (nbondi == -1) .OR. (nbondi == 2) )   imin = 3
264!            IF( (nbondi == +1) .OR. (nbondi == 2) )   imax = jpi-2
265!
266!            ! smoothed fields
267!            IF( eastern_side ) THEN
268!               ztab(jpi,j1:j2,:) = z1 * ptab(jpi,j1:j2,:) + z2 * ptab(jpi-1,j1:j2,:)
269!               DO jj = jmin, jmax
270!                  rswitch = 0.
271!                  IF( u_ice(jpi-2,jj) > 0._wp ) rswitch = 1.
272!                  ztab(jpi-1,jj,:) = ( 1. - umask(jpi-2,jj,1) ) * ztab(jpi,jj,:)  &
273!                     &               +      umask(jpi-2,jj,1)   *  &
274!                     &               ( (1. - rswitch) * ( z4 * ztab(jpi  ,jj,:) + z3 * ztab(jpi-2,jj,:) )  &
275!                     &                 +     rswitch  * ( z6 * ztab(jpi-2,jj,:) + z5 * ztab(jpi  ,jj,:) + z7 * ztab(jpi-3,jj,:) ) )
276!                  ztab(jpi-1,jj,:) = ztab(jpi-1,jj,:) * tmask(jpi-1,jj,1)
277!               END DO
278!            ENDIF
279!            !
280!            IF( northern_side ) THEN
281!               ztab(i1:i2,jpj,:) = z1 * ptab(i1:i2,jpj,:) + z2 * ptab(i1:i2,jpj-1,:)
282!               DO ji = imin, imax
283!                  rswitch = 0.
284!                  IF( v_ice(ji,jpj-2) > 0._wp ) rswitch = 1.
285!                  ztab(ji,jpj-1,:) = ( 1. - vmask(ji,jpj-2,1) ) * ztab(ji,jpj,:)  &
286!                     &               +      vmask(ji,jpj-2,1)   *  &
287!                     &               ( (1. - rswitch) * ( z4 * ztab(ji,jpj  ,:) + z3 * ztab(ji,jpj-2,:) ) &
288!                     &                 +     rswitch  * ( z6 * ztab(ji,jpj-2,:) + z5 * ztab(ji,jpj  ,:) + z7 * ztab(ji,jpj-3,:) ) )
289!                  ztab(ji,jpj-1,:) = ztab(ji,jpj-1,:) * tmask(ji,jpj-1,1)
290!               END DO
291!            END IF
292!            !
293!            IF( western_side) THEN
294!               ztab(1,j1:j2,:) = z1 * ptab(1,j1:j2,:) + z2 * ptab(2,j1:j2,:)
295!               DO jj = jmin, jmax
296!                  rswitch = 0.
297!                  IF( u_ice(2,jj) < 0._wp ) rswitch = 1.
298!                  ztab(2,jj,:) = ( 1. - umask(2,jj,1) ) * ztab(1,jj,:)  &
299!                     &           +      umask(2,jj,1)   *   &
300!                     &           ( ( 1. - rswitch ) * ( z4 * ztab(1,jj,:) + z3 * ztab(3,jj,:) ) &
301!                     &             +      rswitch   * ( z6 * ztab(3,jj,:) + z5 * ztab(1,jj,:) + z7 * ztab(4,jj,:) ) )
302!                  ztab(2,jj,:) = ztab(2,jj,:) * tmask(2,jj,1)
303!               END DO
304!            ENDIF
305!            !
306!            IF( southern_side ) THEN
307!               ztab(i1:i2,1,:) = z1 * ptab(i1:i2,1,:) + z2 * ptab(i1:i2,2,:)
308!               DO ji = imin, imax
309!                  rswitch = 0.
310!                  IF( v_ice(ji,2) < 0._wp ) rswitch = 1.
311!                  ztab(ji,2,:) = ( 1. - vmask(ji,2,1) ) * ztab(ji,1,:)  &
312!                     &           +      vmask(ji,2,1)   *  &
313!                     &           ( ( 1. - rswitch ) * ( z4 * ztab(ji,1,:) + z3 * ztab(ji,3,:) ) &
314!                     &             +      rswitch   * ( z6 * ztab(ji,3,:) + z5 * ztab(ji,1,:) + z7 * ztab(ji,4,:) ) )
315!                  ztab(ji,2,:) = ztab(ji,2,:) * tmask(ji,2,1)
316!               END DO
317!            END IF
318!            !
319!            ! Treatment of corners
320!            IF( (eastern_side) .AND. ((nbondj == -1).OR.(nbondj == 2)) )  ztab(jpi-1,2    ,:) = ptab(jpi-1,    2,:)   ! East south
321!            IF( (eastern_side) .AND. ((nbondj ==  1).OR.(nbondj == 2)) )  ztab(jpi-1,jpj-1,:) = ptab(jpi-1,jpj-1,:)   ! East north
322!            IF( (western_side) .AND. ((nbondj == -1).OR.(nbondj == 2)) )  ztab(    2,    2,:) = ptab(    2,    2,:)   ! West south
323!            IF( (western_side) .AND. ((nbondj ==  1).OR.(nbondj == 2)) )  ztab(    2,jpj-1,:) = ptab(    2,jpj-1,:)   ! West north
324!           
325!            ! retrieve ice tracers
326!            jm = 1
327!            DO jl = 1, jpl
328!               !
329!               DO jj = j1, j2
330!                  DO ji = i1, i2
331!                     a_i (ji,jj,jl) = ztab(ji,jj,jm  ) * tmask(ji,jj,1)
332!                     v_i (ji,jj,jl) = ztab(ji,jj,jm+1) * tmask(ji,jj,1)
333!                     v_s (ji,jj,jl) = ztab(ji,jj,jm+2) * tmask(ji,jj,1)
334!                     sv_i(ji,jj,jl) = ztab(ji,jj,jm+3) * tmask(ji,jj,1)
335!                     oa_i(ji,jj,jl) = ztab(ji,jj,jm+4) * tmask(ji,jj,1)
336!                     a_ip(ji,jj,jl) = ztab(ji,jj,jm+5) * tmask(ji,jj,1)
337!                     v_ip(ji,jj,jl) = ztab(ji,jj,jm+6) * tmask(ji,jj,1)
338!                     t_su(ji,jj,jl) = ztab(ji,jj,jm+7) * tmask(ji,jj,1)
339!                  END DO
340!               END DO
341!               jm = jm + 8
342!               !
343!               DO jk = 1, nlay_s
344!                  e_s(i1:i2,j1:j2,jk,jl) = ztab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1)
345!                  jm = jm + 1
346!               END DO
347!               !
348!               DO jk = 1, nlay_i
349!                  e_i(i1:i2,j1:j2,jk,jl) = ztab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1)
350!                  jm = jm + 1
351!               END DO
352!               !
353!            END DO
354!         
355!         ENDIF  ! nbghostcells=1
356         
357         DO jl = 1, jpl
358            WHERE( tmask(i1:i2,j1:j2,1) == 0._wp )   t_su(i1:i2,j1:j2,jl) = rt0   ! to avoid a division by 0 in sbcblk.F90
359         END DO
360         !
361      ENDIF
362     
363      DEALLOCATE( ztab )
364      !
365   END SUBROUTINE interp_tra_ice
366
367#else
368   !!----------------------------------------------------------------------
369   !!   Empty module                                             no sea-ice
370   !!----------------------------------------------------------------------
371CONTAINS
372   SUBROUTINE agrif_ice_interp_empty
373      WRITE(*,*)  'agrif_ice_interp : You should not have seen this print! error?'
374   END SUBROUTINE agrif_ice_interp_empty
375#endif
376
377   !!======================================================================
378END MODULE agrif_ice_interp
Note: See TracBrowser for help on using the repository browser.