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_CNRS_2017/NEMOGCM/NEMO/NST_SRC – NEMO

source: branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/NST_SRC/agrif_lim3_interp.F90 @ 8882

Last change on this file since 8882 was 8882, checked in by flavoni, 6 years ago

dev_CNRS_2017 branch: merged dev_r7881_ENHANCE09_RK3 with trunk r8864
File size: 17.8 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   ! called by agrif_user.F90
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   ! local scalar
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')   ;   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_lim3
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_b(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_b(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_b (i1:i2,j1:j2,jl)
163            ptab(i1:i2,j1:j2,jm+1) = v_i_b (i1:i2,j1:j2,jl)
164            ptab(i1:i2,j1:j2,jm+2) = v_s_b (i1:i2,j1:j2,jl)
165            ptab(i1:i2,j1:j2,jm+3) = sv_i_b(i1:i2,j1:j2,jl)
166            ptab(i1:i2,j1:j2,jm+4) = oa_i_b(i1:i2,j1:j2,jl)
167            jm = jm + 5
168            DO jk = 1, nlay_s
169               ptab(i1:i2,j1:j2,jm) = e_s_b(i1:i2,j1:j2,jk,jl)   ;   jm = jm + 1
170            END DO
171            DO jk = 1, nlay_i
172               ptab(i1:i2,j1:j2,jm) = e_i_b(i1:i2,j1:j2,jk,jl)   ;   jm = jm + 1
173            END DO
174         END DO
175         
176         DO jk = k1, k2
177            WHERE( tmask(i1:i2,j1:j2,1) == 0._wp )   ptab(i1:i2,j1:j2,jk) = Agrif_SpecialValue
178         END DO
179         !
180      ELSE               ! child grid
181         !
182         IF( nbghostcells > 1 ) THEN   ! ==> The easiest interpolation is used
183            !
184            jm = 1
185            DO jl = 1, jpl
186               !
187               DO jj = j1, j2
188                  DO ji = i1, i2
189                     a_i (ji,jj,jl) = ptab(ji,jj,jm  ) * tmask(ji,jj,1)
190                     v_i (ji,jj,jl) = ptab(ji,jj,jm+1) * tmask(ji,jj,1)
191                     v_s (ji,jj,jl) = ptab(ji,jj,jm+2) * tmask(ji,jj,1)
192                     sv_i(ji,jj,jl) = ptab(ji,jj,jm+3) * tmask(ji,jj,1)
193                     oa_i(ji,jj,jl) = ptab(ji,jj,jm+4) * tmask(ji,jj,1)
194                  END DO
195               END DO
196               jm = jm + 5
197               !
198               DO jk = 1, nlay_s
199                  e_s(i1:i2,j1:j2,jk,jl) = ptab(:,:,jm) * tmask(i1:i2,j1:j2,1)
200                  jm = jm + 1
201               END DO
202               !
203               DO jk = 1, nlay_i
204                  e_i(i1:i2,j1:j2,jk,jl) = ptab(:,:,jm) * tmask(i1:i2,j1:j2,1)
205                  jm = jm + 1
206               END DO
207               !
208            END DO
209            !
210         ELSE                          ! ==> complex interpolation (only one ghost cell available)
211            !! Use a more complex interpolation since we mix solutions over a couple of grid points
212            !! it is advised to use it for fields modified by high order schemes (e.g. advection UM5...)
213            !! clem: for some reason (I don't know why), the following lines do not work
214            !!       with mpp (or in realistic configurations?). It makes the model crash
215            !        I think there is an issue with Agrif_SpecialValue here (not taken into account properly)
216            ! record ztab
217            jm = 1
218            DO jl = 1, jpl
219               ztab(:,:,jm  ) = a_i  (:,:,jl)
220               ztab(:,:,jm+1) = v_i  (:,:,jl)
221               ztab(:,:,jm+2) = v_s  (:,:,jl)
222               ztab(:,:,jm+3) = sv_i(:,:,jl)
223               ztab(:,:,jm+4) = oa_i(:,:,jl)
224               jm = jm + 5
225               DO jk = 1, nlay_s
226                  ztab(:,:,jm) = e_s(:,:,jk,jl)
227                  jm = jm + 1
228               END DO
229               DO jk = 1, nlay_i
230                  ztab(:,:,jm) = e_i(:,:,jk,jl)
231                  jm = jm + 1
232               END DO
233               !
234            END DO
235            !
236            ! borders of the domain
237            western_side  = (nb == 1).AND.(ndir == 1)  ;  eastern_side  = (nb == 1).AND.(ndir == 2)
238            southern_side = (nb == 2).AND.(ndir == 1)  ;  northern_side = (nb == 2).AND.(ndir == 2)
239            !
240            ! spatial smoothing
241            zrhox = Agrif_Rhox()
242            z1 =      ( zrhox - 1. ) * 0.5 
243            z3 =      ( zrhox - 1. ) / ( zrhox + 1. )
244            z6 = 2. * ( zrhox - 1. ) / ( zrhox + 1. )
245            z7 =    - ( zrhox - 1. ) / ( zrhox + 3. )
246            z2 = 1. - z1
247            z4 = 1. - z3
248            z5 = 1. - z6 - z7
249            !
250            ! Remove corners
251            imin = i1  ;  imax = i2  ;  jmin = j1  ;  jmax = j2
252            IF( (nbondj == -1) .OR. (nbondj == 2) )   jmin = 3
253            IF( (nbondj == +1) .OR. (nbondj == 2) )   jmax = nlcj-2
254            IF( (nbondi == -1) .OR. (nbondi == 2) )   imin = 3
255            IF( (nbondi == +1) .OR. (nbondi == 2) )   imax = nlci-2
256
257            ! smoothed fields
258            IF( eastern_side ) THEN
259               ztab(nlci,j1:j2,:) = z1 * ptab(nlci,j1:j2,:) + z2 * ptab(nlci-1,j1:j2,:)
260               DO jj = jmin, jmax
261                  rswitch = 0.
262                  IF( u_ice(nlci-2,jj) > 0._wp ) rswitch = 1.
263                  ztab(nlci-1,jj,:) = ( 1. - umask(nlci-2,jj,1) ) * ztab(nlci,jj,:)  &
264                     &                +      umask(nlci-2,jj,1)   *  &
265                     &                ( ( 1. - rswitch ) * ( z4 * ztab(nlci,jj,:)   + z3 * ztab(nlci-2,jj,:) )  &
266                     &                  +      rswitch   * ( z6 * ztab(nlci-2,jj,:) + z5 * ztab(nlci,jj,:) + z7 * ztab(nlci-3,jj,:) ) )
267                  ztab(nlci-1,jj,:) = ztab(nlci-1,jj,:) * tmask(nlci-1,jj,1)
268               END DO
269            ENDIF
270            !
271            IF( northern_side ) THEN
272               ztab(i1:i2,nlcj,:) = z1 * ptab(i1:i2,nlcj,:) + z2 * ptab(i1:i2,nlcj-1,:)
273               DO ji = imin, imax
274                  rswitch = 0.
275                  IF( v_ice(ji,nlcj-2) > 0._wp ) rswitch = 1.
276                  ztab(ji,nlcj-1,:) = ( 1. - vmask(ji,nlcj-2,1) ) * ztab(ji,nlcj,:)  &
277                     &                +      vmask(ji,nlcj-2,1)   *  &
278                     &                ( ( 1. - rswitch ) * ( z4 * ztab(ji,nlcj,:)   + z3 * ztab(ji,nlcj-2,:) ) &
279                     &                  +      rswitch   * ( z6 * ztab(ji,nlcj-2,:) + z5 * ztab(ji,nlcj,:) + z7 * ztab(ji,nlcj-3,:) ) )
280                  ztab(ji,nlcj-1,:) = ztab(ji,nlcj-1,:) * tmask(ji,nlcj-1,1)
281               END DO
282            END IF
283            !
284            IF( western_side) THEN
285               ztab(1,j1:j2,:) = z1 * ptab(1,j1:j2,:) + z2 * ptab(2,j1:j2,:)
286               DO jj = jmin, jmax
287                  rswitch = 0.
288                  IF( u_ice(2,jj) < 0._wp ) rswitch = 1.
289                  ztab(2,jj,:) = ( 1. - umask(2,jj,1) ) * ztab(1,jj,:)  &
290                     &           +      umask(2,jj,1)   *   &
291                     &           ( ( 1. - rswitch ) * ( z4 * ztab(1,jj,:) + z3 * ztab(3,jj,:) ) &
292                     &             +      rswitch   * ( z6 * ztab(3,jj,:) + z5 * ztab(1,jj,:) + z7 * ztab(4,jj,:) ) )
293                  ztab(2,jj,:) = ztab(2,jj,:) * tmask(2,jj,1)
294               END DO
295            ENDIF
296            !
297            IF( southern_side ) THEN
298               ztab(i1:i2,1,:) = z1 * ptab(i1:i2,1,:) + z2 * ptab(i1:i2,2,:)
299               DO ji = imin, imax
300                  rswitch = 0.
301                  IF( v_ice(ji,2) < 0._wp ) rswitch = 1.
302                  ztab(ji,2,:) = ( 1. - vmask(ji,2,1) ) * ztab(ji,1,:)  &
303                     &           +      vmask(ji,2,1)   *  &
304                     &           ( ( 1. - rswitch ) * ( z4 * ztab(ji,1,:) + z3 * ztab(ji,3,:) ) &
305                     &             +      rswitch   * ( z6 * ztab(ji,3,:) + z5 * ztab(ji,1,:) + z7 * ztab(ji,4,:) ) )
306                  ztab(ji,2,:) = ztab(ji,2,:) * tmask(ji,2,1)
307               END DO
308            END IF
309            !
310            ! Treatment of corners
311            IF( (eastern_side) .AND. ((nbondj == -1).OR.(nbondj == 2)) )  ztab(nlci-1,2,:)      = ptab(nlci-1,2,:)      ! East south
312            IF( (eastern_side) .AND. ((nbondj ==  1).OR.(nbondj == 2)) )  ztab(nlci-1,nlcj-1,:) = ptab(nlci-1,nlcj-1,:) ! East north
313            IF( (western_side) .AND. ((nbondj == -1).OR.(nbondj == 2)) )  ztab(2,2,:)           = ptab(2,2,:)           ! West south
314            IF( (western_side) .AND. ((nbondj ==  1).OR.(nbondj == 2)) )  ztab(2,nlcj-1,:)      = ptab(2,nlcj-1,:)      ! West north
315           
316            ! retrieve ice tracers
317            jm = 1
318            DO jl = 1, jpl
319               !
320               DO jj = j1, j2
321                  DO ji = i1, i2
322                     a_i (ji,jj,jl) = ztab(ji,jj,jm  ) * tmask(ji,jj,1)
323                     v_i (ji,jj,jl) = ztab(ji,jj,jm+1) * tmask(ji,jj,1)
324                     v_s (ji,jj,jl) = ztab(ji,jj,jm+2) * tmask(ji,jj,1)
325                     sv_i(ji,jj,jl) = ztab(ji,jj,jm+3) * tmask(ji,jj,1)
326                     oa_i (ji,jj,jl) = ztab(ji,jj,jm+4) * tmask(ji,jj,1)
327                  END DO
328               END DO
329               jm = jm + 5
330               !
331               DO jk = 1, nlay_s
332                  e_s(i1:i2,j1:j2,jk,jl) = ztab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1)
333                  jm = jm + 1
334               END DO
335               !
336               DO jk = 1, nlay_i
337                  e_i(i1:i2,j1:j2,jk,jl) = ztab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1)
338                  jm = jm + 1
339               END DO
340               !
341            END DO
342         
343         ENDIF  ! nbghostcells=1
344         
345         ! integrated values
346         vt_i (i1:i2,j1:j2) = SUM( v_i(i1:i2,j1:j2,:), dim=3 )
347         vt_s (i1:i2,j1:j2) = SUM( v_s(i1:i2,j1:j2,:), dim=3 )
348         at_i (i1:i2,j1:j2) = SUM( a_i(i1:i2,j1:j2,:), dim=3 )
349         et_s(i1:i2,j1:j2)  = SUM( SUM( e_s(i1:i2,j1:j2,:,:), dim=4 ), dim=3 )
350         et_i(i1:i2,j1:j2)  = SUM( SUM( e_i(i1:i2,j1:j2,:,:), dim=4 ), dim=3 )
351         !
352      ENDIF
353     
354      DEALLOCATE( ztab )
355      !
356   END SUBROUTINE interp_tra_ice
357
358#else
359   !!----------------------------------------------------------------------
360   !!   Empty module                                             no sea-ice
361   !!----------------------------------------------------------------------
362CONTAINS
363   SUBROUTINE agrif_lim3_interp_empty
364      WRITE(*,*)  'agrif_lim3_interp : You should not have seen this print! error?'
365   END SUBROUTINE agrif_lim3_interp_empty
366#endif
367
368   !!======================================================================
369END MODULE agrif_lim3_interp
Note: See TracBrowser for help on using the repository browser.