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/r12377_ticket2386/src/NST – NEMO

source: NEMO/branches/2020/r12377_ticket2386/src/NST/agrif_ice_interp.F90 @ 13540

Last change on this file since 13540 was 13540, checked in by andmirek, 4 years ago

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