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/2019/dev_r11842_SI3-10_EAP/src/NST – NEMO

source: NEMO/branches/2019/dev_r11842_SI3-10_EAP/src/NST/agrif_ice_interp.F90 @ 13662

Last change on this file since 13662 was 13662, checked in by clem, 3 years ago

update to almost r4.0.4
  • 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   !!  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) = v_il(i1:i2,j1:j2,jl)
170            ptab(i1:i2,j1:j2,jm+8) = t_su(i1:i2,j1:j2,jl)
171            jm = jm + 9
172            DO jk = 1, nlay_s
173               ptab(i1:i2,j1:j2,jm) = e_s(i1:i2,j1:j2,jk,jl)   ;   jm = jm + 1
174            END DO
175            DO jk = 1, nlay_i
176               ptab(i1:i2,j1:j2,jm) = e_i(i1:i2,j1:j2,jk,jl)   ;   jm = jm + 1
177            END DO
178         END DO
179         
180         DO jk = k1, k2
181            WHERE( tmask(i1:i2,j1:j2,1) == 0._wp )   ptab(i1:i2,j1:j2,jk) = Agrif_SpecialValue
182         END DO
183         !
184      ELSE               ! child grid
185         !
186!         IF( nbghostcells > 1 ) THEN   ! ==> The easiest interpolation is used
187            !
188            jm = 1
189            DO jl = 1, jpl
190               !
191               DO jj = j1, j2
192                  DO ji = i1, i2
193                     a_i (ji,jj,jl) = ptab(ji,jj,jm  ) * tmask(ji,jj,1)
194                     v_i (ji,jj,jl) = ptab(ji,jj,jm+1) * tmask(ji,jj,1)
195                     v_s (ji,jj,jl) = ptab(ji,jj,jm+2) * tmask(ji,jj,1)
196                     sv_i(ji,jj,jl) = ptab(ji,jj,jm+3) * tmask(ji,jj,1)
197                     oa_i(ji,jj,jl) = ptab(ji,jj,jm+4) * tmask(ji,jj,1)
198                     a_ip(ji,jj,jl) = ptab(ji,jj,jm+5) * tmask(ji,jj,1)
199                     v_ip(ji,jj,jl) = ptab(ji,jj,jm+6) * tmask(ji,jj,1)
200                     v_il(ji,jj,jl) = ptab(ji,jj,jm+7) * tmask(ji,jj,1)
201                     t_su(ji,jj,jl) = ptab(ji,jj,jm+8) * tmask(ji,jj,1)
202                  END DO
203               END DO
204               jm = jm + 9
205               !
206               DO jk = 1, nlay_s
207                  e_s(i1:i2,j1:j2,jk,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1)
208                  jm = jm + 1
209               END DO
210               !
211               DO jk = 1, nlay_i
212                  e_i(i1:i2,j1:j2,jk,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1)
213                  jm = jm + 1
214               END DO
215               !
216            END DO
217            !
218!!==> clem: this interpolation does not work because it creates negative values, due
219!!          to negative coefficients when mixing points (for ex. z7)
220!!
221!         ELSE                          ! ==> complex interpolation (only one ghost cell available)
222!            !! Use a more complex interpolation since we mix solutions over a couple of grid points
223!            !! it is advised to use it for fields modified by high order schemes (e.g. advection UM5...)
224!            ! record ztab
225!            jm = 1
226!            DO jl = 1, jpl
227!               ztab(:,:,jm  ) = a_i (:,:,jl)
228!               ztab(:,:,jm+1) = v_i (:,:,jl)
229!               ztab(:,:,jm+2) = v_s (:,:,jl)
230!               ztab(:,:,jm+3) = sv_i(:,:,jl)
231!               ztab(:,:,jm+4) = oa_i(:,:,jl)
232!               ztab(:,:,jm+5) = a_ip(:,:,jl)
233!               ztab(:,:,jm+6) = v_ip(:,:,jl)
234!               ztab(:,:,jm+7) = v_il(:,:,jl)
235!               ztab(:,:,jm+8) = t_su(:,:,jl)
236!               jm = jm + 9
237!               DO jk = 1, nlay_s
238!                  ztab(:,:,jm) = e_s(:,:,jk,jl)
239!                  jm = jm + 1
240!               END DO
241!               DO jk = 1, nlay_i
242!                  ztab(:,:,jm) = e_i(:,:,jk,jl)
243!                  jm = jm + 1
244!               END DO
245!               !
246!            END DO
247!            !
248!            ! borders of the domain
249!            western_side  = (nb == 1).AND.(ndir == 1)  ;  eastern_side  = (nb == 1).AND.(ndir == 2)
250!            southern_side = (nb == 2).AND.(ndir == 1)  ;  northern_side = (nb == 2).AND.(ndir == 2)
251!            !
252!            ! spatial smoothing
253!            zrhox = Agrif_Rhox()
254!            z1 =      ( zrhox - 1. ) * 0.5
255!            z3 =      ( zrhox - 1. ) / ( zrhox + 1. )
256!            z6 = 2. * ( zrhox - 1. ) / ( zrhox + 1. )
257!            z7 =    - ( zrhox - 1. ) / ( zrhox + 3. )
258!            z2 = 1. - z1
259!            z4 = 1. - z3
260!            z5 = 1. - z6 - z7
261!            !
262!            ! Remove corners
263!            imin = i1  ;  imax = i2  ;  jmin = j1  ;  jmax = j2
264!            IF( (nbondj == -1) .OR. (nbondj == 2) )   jmin = 3
265!            IF( (nbondj == +1) .OR. (nbondj == 2) )   jmax = nlcj-2
266!            IF( (nbondi == -1) .OR. (nbondi == 2) )   imin = 3
267!            IF( (nbondi == +1) .OR. (nbondi == 2) )   imax = nlci-2
268!
269!            ! smoothed fields
270!            IF( eastern_side ) THEN
271!               ztab(nlci,j1:j2,:) = z1 * ptab(nlci,j1:j2,:) + z2 * ptab(nlci-1,j1:j2,:)
272!               DO jj = jmin, jmax
273!                  rswitch = 0.
274!                  IF( u_ice(nlci-2,jj) > 0._wp ) rswitch = 1.
275!                  ztab(nlci-1,jj,:) = ( 1. - umask(nlci-2,jj,1) ) * ztab(nlci,jj,:)  &
276!                     &                +      umask(nlci-2,jj,1)   *  &
277!                     &                ( ( 1. - rswitch ) * ( z4 * ztab(nlci,jj,:)   + z3 * ztab(nlci-2,jj,:) )  &
278!                     &                  +      rswitch   * ( z6 * ztab(nlci-2,jj,:) + z5 * ztab(nlci,jj,:) + z7 * ztab(nlci-3,jj,:) ) )
279!                  ztab(nlci-1,jj,:) = ztab(nlci-1,jj,:) * tmask(nlci-1,jj,1)
280!               END DO
281!            ENDIF
282!            !
283!            IF( northern_side ) THEN
284!               ztab(i1:i2,nlcj,:) = z1 * ptab(i1:i2,nlcj,:) + z2 * ptab(i1:i2,nlcj-1,:)
285!               DO ji = imin, imax
286!                  rswitch = 0.
287!                  IF( v_ice(ji,nlcj-2) > 0._wp ) rswitch = 1.
288!                  ztab(ji,nlcj-1,:) = ( 1. - vmask(ji,nlcj-2,1) ) * ztab(ji,nlcj,:)  &
289!                     &                +      vmask(ji,nlcj-2,1)   *  &
290!                     &                ( ( 1. - rswitch ) * ( z4 * ztab(ji,nlcj,:)   + z3 * ztab(ji,nlcj-2,:) ) &
291!                     &                  +      rswitch   * ( z6 * ztab(ji,nlcj-2,:) + z5 * ztab(ji,nlcj,:) + z7 * ztab(ji,nlcj-3,:) ) )
292!                  ztab(ji,nlcj-1,:) = ztab(ji,nlcj-1,:) * tmask(ji,nlcj-1,1)
293!               END DO
294!            END IF
295!            !
296!            IF( western_side) THEN
297!               ztab(1,j1:j2,:) = z1 * ptab(1,j1:j2,:) + z2 * ptab(2,j1:j2,:)
298!               DO jj = jmin, jmax
299!                  rswitch = 0.
300!                  IF( u_ice(2,jj) < 0._wp ) rswitch = 1.
301!                  ztab(2,jj,:) = ( 1. - umask(2,jj,1) ) * ztab(1,jj,:)  &
302!                     &           +      umask(2,jj,1)   *   &
303!                     &           ( ( 1. - rswitch ) * ( z4 * ztab(1,jj,:) + z3 * ztab(3,jj,:) ) &
304!                     &             +      rswitch   * ( z6 * ztab(3,jj,:) + z5 * ztab(1,jj,:) + z7 * ztab(4,jj,:) ) )
305!                  ztab(2,jj,:) = ztab(2,jj,:) * tmask(2,jj,1)
306!               END DO
307!            ENDIF
308!            !
309!            IF( southern_side ) THEN
310!               ztab(i1:i2,1,:) = z1 * ptab(i1:i2,1,:) + z2 * ptab(i1:i2,2,:)
311!               DO ji = imin, imax
312!                  rswitch = 0.
313!                  IF( v_ice(ji,2) < 0._wp ) rswitch = 1.
314!                  ztab(ji,2,:) = ( 1. - vmask(ji,2,1) ) * ztab(ji,1,:)  &
315!                     &           +      vmask(ji,2,1)   *  &
316!                     &           ( ( 1. - rswitch ) * ( z4 * ztab(ji,1,:) + z3 * ztab(ji,3,:) ) &
317!                     &             +      rswitch   * ( z6 * ztab(ji,3,:) + z5 * ztab(ji,1,:) + z7 * ztab(ji,4,:) ) )
318!                  ztab(ji,2,:) = ztab(ji,2,:) * tmask(ji,2,1)
319!               END DO
320!            END IF
321!            !
322!            ! Treatment of corners
323!            IF( (eastern_side) .AND. ((nbondj == -1).OR.(nbondj == 2)) )  ztab(nlci-1,2,:)      = ptab(nlci-1,2,:)      ! East south
324!            IF( (eastern_side) .AND. ((nbondj ==  1).OR.(nbondj == 2)) )  ztab(nlci-1,nlcj-1,:) = ptab(nlci-1,nlcj-1,:) ! East north
325!            IF( (western_side) .AND. ((nbondj == -1).OR.(nbondj == 2)) )  ztab(2,2,:)           = ptab(2,2,:)           ! West south
326!            IF( (western_side) .AND. ((nbondj ==  1).OR.(nbondj == 2)) )  ztab(2,nlcj-1,:)      = ptab(2,nlcj-1,:)      ! West north
327!           
328!            ! retrieve ice tracers
329!            jm = 1
330!            DO jl = 1, jpl
331!               !
332!               DO jj = j1, j2
333!                  DO ji = i1, i2
334!                     a_i (ji,jj,jl) = ztab(ji,jj,jm  ) * tmask(ji,jj,1)
335!                     v_i (ji,jj,jl) = ztab(ji,jj,jm+1) * tmask(ji,jj,1)
336!                     v_s (ji,jj,jl) = ztab(ji,jj,jm+2) * tmask(ji,jj,1)
337!                     sv_i(ji,jj,jl) = ztab(ji,jj,jm+3) * tmask(ji,jj,1)
338!                     oa_i(ji,jj,jl) = ztab(ji,jj,jm+4) * tmask(ji,jj,1)
339!                     a_ip(ji,jj,jl) = ztab(ji,jj,jm+5) * tmask(ji,jj,1)
340!                     v_ip(ji,jj,jl) = ztab(ji,jj,jm+6) * tmask(ji,jj,1)
341!                     v_il(ji,jj,jl) = ztab(ji,jj,jm+7) * tmask(ji,jj,1)
342!                     t_su(ji,jj,jl) = ztab(ji,jj,jm+8) * tmask(ji,jj,1)
343!                  END DO
344!               END DO
345!               jm = jm + 9
346!               !
347!               DO jk = 1, nlay_s
348!                  e_s(i1:i2,j1:j2,jk,jl) = ztab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1)
349!                  jm = jm + 1
350!               END DO
351!               !
352!               DO jk = 1, nlay_i
353!                  e_i(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!            END DO
358!         
359!         ENDIF  ! nbghostcells=1
360         
361         DO jl = 1, jpl
362            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
363         END DO
364         !
365      ENDIF
366     
367      DEALLOCATE( ztab )
368      !
369   END SUBROUTINE interp_tra_ice
370
371#else
372   !!----------------------------------------------------------------------
373   !!   Empty module                                             no sea-ice
374   !!----------------------------------------------------------------------
375CONTAINS
376   SUBROUTINE agrif_ice_interp_empty
377      WRITE(*,*)  'agrif_ice_interp : You should not have seen this print! error?'
378   END SUBROUTINE agrif_ice_interp_empty
379#endif
380
381   !!======================================================================
382END MODULE agrif_ice_interp
Note: See TracBrowser for help on using the repository browser.