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

Last change on this file since 13229 was 13229, checked in by francesca, 4 years ago

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