source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/NST_SRC/agrif_lim3_interp.F90 @ 9167

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

fix issues but agrif + lim3 is still not restartable because of interpolation when ice thermodynamics is activated

File size: 18.5 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   ! 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) / ( 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(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(:,:,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(:,:,jm) * tmask(i1:i2,j1:j2,1)
211                  jm = jm + 1
212               END DO
213               !
214            END DO
215            !
216         ELSE                          ! ==> complex interpolation (only one ghost cell available)
217            !! Use a more complex interpolation since we mix solutions over a couple of grid points
218            !! it is advised to use it for fields modified by high order schemes (e.g. advection UM5...)
219            !! clem: for some reason (I don't know why), the following lines do not work
220            !!       with mpp (or in realistic configurations?). It makes the model crash
221            !        I think there is an issue with Agrif_SpecialValue here (not taken into account properly)
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 = nlcj-2
263            IF( (nbondi == -1) .OR. (nbondi == 2) )   imin = 3
264            IF( (nbondi == +1) .OR. (nbondi == 2) )   imax = nlci-2
265
266            ! smoothed fields
267            IF( eastern_side ) THEN
268               ztab(nlci,j1:j2,:) = z1 * ptab(nlci,j1:j2,:) + z2 * ptab(nlci-1,j1:j2,:)
269               DO jj = jmin, jmax
270                  rswitch = 0.
271                  IF( u_ice(nlci-2,jj) > 0._wp ) rswitch = 1.
272                  ztab(nlci-1,jj,:) = ( 1. - umask(nlci-2,jj,1) ) * ztab(nlci,jj,:)  &
273                     &                +      umask(nlci-2,jj,1)   *  &
274                     &                ( ( 1. - rswitch ) * ( z4 * ztab(nlci,jj,:)   + z3 * ztab(nlci-2,jj,:) )  &
275                     &                  +      rswitch   * ( z6 * ztab(nlci-2,jj,:) + z5 * ztab(nlci,jj,:) + z7 * ztab(nlci-3,jj,:) ) )
276                  ztab(nlci-1,jj,:) = ztab(nlci-1,jj,:) * tmask(nlci-1,jj,1)
277               END DO
278            ENDIF
279            !
280            IF( northern_side ) THEN
281               ztab(i1:i2,nlcj,:) = z1 * ptab(i1:i2,nlcj,:) + z2 * ptab(i1:i2,nlcj-1,:)
282               DO ji = imin, imax
283                  rswitch = 0.
284                  IF( v_ice(ji,nlcj-2) > 0._wp ) rswitch = 1.
285                  ztab(ji,nlcj-1,:) = ( 1. - vmask(ji,nlcj-2,1) ) * ztab(ji,nlcj,:)  &
286                     &                +      vmask(ji,nlcj-2,1)   *  &
287                     &                ( ( 1. - rswitch ) * ( z4 * ztab(ji,nlcj,:)   + z3 * ztab(ji,nlcj-2,:) ) &
288                     &                  +      rswitch   * ( z6 * ztab(ji,nlcj-2,:) + z5 * ztab(ji,nlcj,:) + z7 * ztab(ji,nlcj-3,:) ) )
289                  ztab(ji,nlcj-1,:) = ztab(ji,nlcj-1,:) * tmask(ji,nlcj-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(nlci-1,2,:)      = ptab(nlci-1,2,:)      ! East south
321            IF( (eastern_side) .AND. ((nbondj ==  1).OR.(nbondj == 2)) )  ztab(nlci-1,nlcj-1,:) = ptab(nlci-1,nlcj-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,nlcj-1,:)      = ptab(2,nlcj-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         ! integrated values
358         vt_i (i1:i2,j1:j2) = SUM(      v_i(i1:i2,j1:j2,:),            dim=3 )
359         vt_s (i1:i2,j1:j2) = SUM(      v_s(i1:i2,j1:j2,:),            dim=3 )
360         at_i (i1:i2,j1:j2) = SUM(      a_i(i1:i2,j1:j2,:),            dim=3 )
361         et_s(i1:i2,j1:j2)  = SUM( SUM( e_s(i1:i2,j1:j2,:,:), dim=4 ), dim=3 )
362         et_i(i1:i2,j1:j2)  = SUM( SUM( e_i(i1:i2,j1:j2,:,:), dim=4 ), dim=3 )
363         !
364      ENDIF
365     
366      DEALLOCATE( ztab )
367      !
368   END SUBROUTINE interp_tra_ice
369
370#else
371   !!----------------------------------------------------------------------
372   !!   Empty module                                             no sea-ice
373   !!----------------------------------------------------------------------
374CONTAINS
375   SUBROUTINE agrif_lim3_interp_empty
376      WRITE(*,*)  'agrif_lim3_interp : You should not have seen this print! error?'
377   END SUBROUTINE agrif_lim3_interp_empty
378#endif
379
380   !!======================================================================
381END MODULE agrif_lim3_interp
Note: See TracBrowser for help on using the repository browser.