source: trunk/NEMOGCM/NEMO/NST_SRC/agrif_lim3_interp.F90 @ 7795

Last change on this file since 7795 was 7761, checked in by clem, 4 years ago

make AGRIF and LIM3 fully compatible

File size: 16.2 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
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
52      !!-----------------------------------------------------------------------
53      !
54      IF( Agrif_Root() )  RETURN
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')
72         CALL Agrif_Bc_variable( u_ice_id  , procname=interp_u_ice  , calledweight=zbeta )
73      CASE('V')
74         CALL Agrif_Bc_variable( v_ice_id  , procname=interp_v_ice  , calledweight=zbeta )
75      CASE('T')
76         CALL Agrif_Bc_variable( tra_ice_id, procname=interp_tra_ice, calledweight=zbeta )
77      END SELECT
78      Agrif_SpecialValue=0.
79      Agrif_UseSpecialValue = .FALSE.
80      !
81   END SUBROUTINE agrif_interp_lim3
82
83   !!------------------
84   !! Local subroutines
85   !!------------------
86   SUBROUTINE interp_u_ice( ptab, i1, i2, j1, j2, before )
87      !!-----------------------------------------------------------------------
88      !!                     *** ROUTINE interp_u_ice ***
89      !!
90      !! i1 i2 j1 j2 are the index of the boundaries parent(when before) and child (when after)
91      !! To solve issues when parent grid is "land" masked but not all the corresponding child grid points,
92      !! put -9999 WHERE the parent grid is masked. The child solution will be found in the 9(?) points around
93      !!-----------------------------------------------------------------------
94      INTEGER , INTENT(in) :: i1, i2, j1, j2
95      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab
96      LOGICAL , INTENT(in) :: before
97      !!
98      REAL(wp) :: zrhoy
99      !!-----------------------------------------------------------------------
100      !
101      IF( before ) THEN  ! parent grid
102         ptab(:,:) = e2u(i1:i2,j1:j2) * u_ice_b(i1:i2,j1:j2)
103         WHERE( umask(i1:i2,j1:j2,1) == 0. )  ptab(:,:) = -9999.
104      ELSE               ! child grid
105         zrhoy = Agrif_Rhoy()
106         u_ice(i1:i2,j1:j2) = ptab(:,:) / ( e2u(i1:i2,j1:j2) * zrhoy ) * umask(i1:i2,j1:j2,1)
107      ENDIF
108      !
109   END SUBROUTINE interp_u_ice
110
111
112   SUBROUTINE interp_v_ice( ptab, i1, i2, j1, j2, before )
113      !!-----------------------------------------------------------------------
114      !!                    *** ROUTINE interp_v_ice ***
115      !!
116      !! i1 i2 j1 j2 are the index of the boundaries parent(when before) and child (when after)
117      !! To solve issues when parent grid is "land" masked but not all the corresponding child grid points,
118      !! put -9999 WHERE the parent grid is masked. The child solution will be found in the 9(?) points around
119      !!-----------------------------------------------------------------------     
120      INTEGER , INTENT(in) :: i1, i2, j1, j2
121      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab
122      LOGICAL , INTENT(in) :: before
123      !!
124      REAL(wp) :: zrhox
125      !!-----------------------------------------------------------------------
126      !
127      IF( before ) THEN  ! parent grid
128         ptab(:,:) = e1v(i1:i2,j1:j2) * v_ice_b(i1:i2,j1:j2)
129         WHERE( vmask(i1:i2,j1:j2,1) == 0. )  ptab(:,:) = -9999.
130      ELSE               ! child grid
131         zrhox = Agrif_Rhox()
132         v_ice(i1:i2,j1:j2) = ptab(:,:) / ( e1v(i1:i2,j1:j2) * zrhox ) * vmask(i1:i2,j1:j2,1)
133      ENDIF
134      !
135   END SUBROUTINE interp_v_ice
136
137
138   SUBROUTINE interp_tra_ice( ptab, i1, i2, j1, j2, k1, k2, before, nb, ndir )
139      !!-----------------------------------------------------------------------
140      !!                    *** ROUTINE interp_tra_ice ***                           
141      !!
142      !! i1 i2 j1 j2 are the index of the boundaries parent(when before) and child (when after)
143      !! To solve issues when parent grid is "land" masked but not all the corresponding child grid points,
144      !! put -9999 WHERE the parent grid is masked. The child solution will be found in the 9(?) points around
145      !!-----------------------------------------------------------------------
146      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
147      INTEGER , INTENT(in) :: i1, i2, j1, j2, k1, k2
148      LOGICAL , INTENT(in) :: before
149      INTEGER , INTENT(in) :: nb, ndir
150      !!
151      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztab
152      INTEGER  ::   ji, jj, jk, jl, jm
153      INTEGER  ::   imin, imax, jmin, jmax
154      REAL(wp) ::   zrhox, z1, z2, z3, z4, z5, z6, z7
155      LOGICAL  ::   western_side, eastern_side, northern_side, southern_side
156
157      !!-----------------------------------------------------------------------
158      ! tracers are not multiplied by grid cell here => before: * e12t ; after: * r1_e12t / rhox / rhoy
159      ! and it is ok since we conserve tracers (same as in the ocean).
160      ALLOCATE( ztab(SIZE(a_i_b,1),SIZE(a_i_b,2),SIZE(ptab,3)) )
161     
162      IF( before ) THEN  ! parent grid
163         jm = 1
164         DO jl = 1, jpl
165            ptab(i1:i2,j1:j2,jm) = a_i_b  (i1:i2,j1:j2,jl) ; jm = jm + 1
166            ptab(i1:i2,j1:j2,jm) = v_i_b  (i1:i2,j1:j2,jl) ; jm = jm + 1
167            ptab(i1:i2,j1:j2,jm) = v_s_b  (i1:i2,j1:j2,jl) ; jm = jm + 1
168            ptab(i1:i2,j1:j2,jm) = smv_i_b(i1:i2,j1:j2,jl) ; jm = jm + 1
169            ptab(i1:i2,j1:j2,jm) = oa_i_b (i1:i2,j1:j2,jl) ; jm = jm + 1
170            DO jk = 1, nlay_s
171               ptab(i1:i2,j1:j2,jm) = e_s_b(i1:i2,j1:j2,jk,jl) ; jm = jm + 1
172            ENDDO
173            DO jk = 1, nlay_i
174               ptab(i1:i2,j1:j2,jm) = e_i_b(i1:i2,j1:j2,jk,jl) ; jm = jm + 1
175            ENDDO
176         ENDDO
177         
178         DO jk = k1, k2
179            WHERE( tmask(i1:i2,j1:j2,1) == 0. )  ptab(i1:i2,j1:j2,jk) = -9999.
180         ENDDO
181         
182      ELSE               ! child grid
183!! ==> The easiest interpolation is the following commented lines
184         jm = 1
185         DO jl = 1, jpl
186            a_i  (i1:i2,j1:j2,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
187            v_i  (i1:i2,j1:j2,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
188            v_s  (i1:i2,j1:j2,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
189            smv_i(i1:i2,j1:j2,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
190            oa_i (i1:i2,j1:j2,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
191            DO jk = 1, nlay_s
192               e_s(i1:i2,j1:j2,jk,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
193            ENDDO
194            DO jk = 1, nlay_i
195               e_i(i1:i2,j1:j2,jk,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
196            ENDDO
197         ENDDO
198
199!! ==> this is a more complex interpolation since we mix solutions over a couple of grid points
200!!     it is advised to use it for fields modified by high order schemes (e.g. advection UM5...)
201!!        clem: for some reason (I don't know why), the following lines do not work
202!!              with mpp (or in realistic configurations?). It makes the model crash
203!         ! record ztab
204!         jm = 1
205!         DO jl = 1, jpl
206!            ztab(:,:,jm) = a_i  (:,:,jl) ; jm = jm + 1
207!            ztab(:,:,jm) = v_i  (:,:,jl) ; jm = jm + 1
208!            ztab(:,:,jm) = v_s  (:,:,jl) ; jm = jm + 1
209!            ztab(:,:,jm) = smv_i(:,:,jl) ; jm = jm + 1
210!            ztab(:,:,jm) = oa_i (:,:,jl) ; jm = jm + 1
211!            DO jk = 1, nlay_s
212!               ztab(:,:,jm) = e_s(:,:,jk,jl) ; jm = jm + 1
213!            ENDDO
214!            DO jk = 1, nlay_i
215!               ztab(:,:,jm) = e_i(:,:,jk,jl) ; jm = jm + 1
216!            ENDDO
217!         ENDDO
218!         !
219!         ! borders of the domain
220!         western_side  = (nb == 1).AND.(ndir == 1)  ;  eastern_side  = (nb == 1).AND.(ndir == 2)
221!         southern_side = (nb == 2).AND.(ndir == 1)  ;  northern_side = (nb == 2).AND.(ndir == 2)
222!         !
223!         ! spatial smoothing
224!         zrhox = Agrif_Rhox()
225!         z1 =      ( zrhox - 1. ) * 0.5
226!         z3 =      ( zrhox - 1. ) / ( zrhox + 1. )
227!         z6 = 2. * ( zrhox - 1. ) / ( zrhox + 1. )
228!         z7 =    - ( zrhox - 1. ) / ( zrhox + 3. )
229!         z2 = 1. - z1
230!         z4 = 1. - z3
231!         z5 = 1. - z6 - z7
232!         !
233!         ! Remove corners
234!         imin = i1  ;  imax = i2  ;  jmin = j1  ;  jmax = j2
235!         IF( (nbondj == -1) .OR. (nbondj == 2) )   jmin = 3
236!         IF( (nbondj == +1) .OR. (nbondj == 2) )   jmax = nlcj-2
237!         IF( (nbondi == -1) .OR. (nbondi == 2) )   imin = 3
238!         IF( (nbondi == +1) .OR. (nbondi == 2) )   imax = nlci-2
239!
240!         ! smoothed fields
241!         IF( eastern_side ) THEN
242!            ztab(nlci,j1:j2,:) = z1 * ptab(nlci,j1:j2,:) + z2 * ptab(nlci-1,j1:j2,:)
243!            DO jj = jmin, jmax
244!               rswitch = 0.
245!               IF( u_ice(nlci-2,jj) > 0._wp ) rswitch = 1.
246!               ztab(nlci-1,jj,:) = ( 1. - umask(nlci-2,jj,1) ) * ztab(nlci,jj,:)  &
247!                  &                +      umask(nlci-2,jj,1)   *  &
248!                  &                ( ( 1. - rswitch ) * ( z4 * ztab(nlci,jj,:)   + z3 * ztab(nlci-2,jj,:) )  &
249!                  &                  +      rswitch   * ( z6 * ztab(nlci-2,jj,:) + z5 * ztab(nlci,jj,:) + z7 * ztab(nlci-3,jj,:) ) )
250!               ztab(nlci-1,jj,:) = ztab(nlci-1,jj,:) * tmask(nlci-1,jj,1)
251!            END DO
252!         ENDIF
253!         !
254!         IF( northern_side ) THEN
255!            ztab(i1:i2,nlcj,:) = z1 * ptab(i1:i2,nlcj,:) + z2 * ptab(i1:i2,nlcj-1,:)
256!            DO ji = imin, imax
257!               rswitch = 0.
258!               IF( v_ice(ji,nlcj-2) > 0._wp ) rswitch = 1.
259!               ztab(ji,nlcj-1,:) = ( 1. - vmask(ji,nlcj-2,1) ) * ztab(ji,nlcj,:)  &
260!                  &                +      vmask(ji,nlcj-2,1)   *  &
261!                  &                ( ( 1. - rswitch ) * ( z4 * ztab(ji,nlcj,:)   + z3 * ztab(ji,nlcj-2,:) ) &
262!                  &                  +      rswitch   * ( z6 * ztab(ji,nlcj-2,:) + z5 * ztab(ji,nlcj,:) + z7 * ztab(ji,nlcj-3,:) ) )
263!               ztab(ji,nlcj-1,:) = ztab(ji,nlcj-1,:) * tmask(ji,nlcj-1,1)
264!            END DO
265!         END IF
266!         !
267!         IF( western_side) THEN
268!            ztab(1,j1:j2,:) = z1 * ptab(1,j1:j2,:) + z2 * ptab(2,j1:j2,:)
269!            DO jj = jmin, jmax
270!               rswitch = 0.
271!               IF( u_ice(2,jj) < 0._wp ) rswitch = 1.
272!               ztab(2,jj,:) = ( 1. - umask(2,jj,1) ) * ztab(1,jj,:)  &
273!                  &           +      umask(2,jj,1)   *   &
274!                  &           ( ( 1. - rswitch ) * ( z4 * ztab(1,jj,:) + z3 * ztab(3,jj,:) ) &
275!                  &             +      rswitch   * ( z6 * ztab(3,jj,:) + z5 * ztab(1,jj,:) + z7 * ztab(4,jj,:) ) )
276!               ztab(2,jj,:) = ztab(2,jj,:) * tmask(2,jj,1)
277!            END DO
278!         ENDIF
279!         !
280!         IF( southern_side ) THEN
281!            ztab(i1:i2,1,:) = z1 * ptab(i1:i2,1,:) + z2 * ptab(i1:i2,2,:)
282!            DO ji = imin, imax
283!               rswitch = 0.
284!               IF( v_ice(ji,2) < 0._wp ) rswitch = 1.
285!               ztab(ji,2,:) = ( 1. - vmask(ji,2,1) ) * ztab(ji,1,:)  &
286!                  &           +      vmask(ji,2,1)   *  &
287!                  &           ( ( 1. - rswitch ) * ( z4 * ztab(ji,1,:) + z3 * ztab(ji,3,:) ) &
288!                  &             +      rswitch   * ( z6 * ztab(ji,3,:) + z5 * ztab(ji,1,:) + z7 * ztab(ji,4,:) ) )
289!               ztab(ji,2,:) = ztab(ji,2,:) * tmask(ji,2,1)
290!            END DO
291!         END IF
292!         !
293!         ! Treatment of corners
294!         IF( (eastern_side) .AND. ((nbondj == -1).OR.(nbondj == 2)) )  ztab(nlci-1,2,:)      = ptab(nlci-1,2,:)      ! East south
295!         IF( (eastern_side) .AND. ((nbondj ==  1).OR.(nbondj == 2)) )  ztab(nlci-1,nlcj-1,:) = ptab(nlci-1,nlcj-1,:) ! East north
296!         IF( (western_side) .AND. ((nbondj == -1).OR.(nbondj == 2)) )  ztab(2,2,:)           = ptab(2,2,:)           ! West south
297!         IF( (western_side) .AND. ((nbondj ==  1).OR.(nbondj == 2)) )  ztab(2,nlcj-1,:)      = ptab(2,nlcj-1,:)      ! West north
298!
299!         ! retrieve ice tracers
300!         jm = 1
301!         DO jl = 1, jpl
302!            a_i  (i1:i2,j1:j2,jl) = ztab(i1:i2,j1:j2,jm) ; jm = jm + 1
303!            v_i  (i1:i2,j1:j2,jl) = ztab(i1:i2,j1:j2,jm) ; jm = jm + 1
304!            v_s  (i1:i2,j1:j2,jl) = ztab(i1:i2,j1:j2,jm) ; jm = jm + 1
305!            smv_i(i1:i2,j1:j2,jl) = ztab(i1:i2,j1:j2,jm) ; jm = jm + 1
306!            oa_i (i1:i2,j1:j2,jl) = ztab(i1:i2,j1:j2,jm) ; jm = jm + 1
307!            DO jk = 1, nlay_s
308!               e_s(i1:i2,j1:j2,jk,jl) = ztab(i1:i2,j1:j2,jm) ; jm = jm + 1
309!            ENDDO
310!            DO jk = 1, nlay_i
311!               e_i(i1:i2,j1:j2,jk,jl) = ztab(i1:i2,j1:j2,jm) ; jm = jm + 1
312!            ENDDO
313!         ENDDO
314       
315         ! integrated values
316         vt_i (i1:i2,j1:j2) = SUM( v_i(i1:i2,j1:j2,:), dim=3 )
317         vt_s (i1:i2,j1:j2) = SUM( v_s(i1:i2,j1:j2,:), dim=3 )
318         at_i (i1:i2,j1:j2) = SUM( a_i(i1:i2,j1:j2,:), dim=3 )
319         et_s(i1:i2,j1:j2)  = SUM( SUM( e_s(i1:i2,j1:j2,:,:), dim=4 ), dim=3 )
320         et_i(i1:i2,j1:j2)  = SUM( SUM( e_i(i1:i2,j1:j2,:,:), dim=4 ), dim=3 )
321
322      ENDIF
323     
324      DEALLOCATE( ztab )
325      !
326   END SUBROUTINE interp_tra_ice
327
328#else
329CONTAINS
330   SUBROUTINE agrif_lim3_interp_empty
331      !!---------------------------------------------
332      !!   *** ROUTINE agrif_lim3_interp_empty ***
333      !!---------------------------------------------
334      WRITE(*,*)  'agrif_lim3_interp : You should not have seen this print! error?'
335   END SUBROUTINE agrif_lim3_interp_empty
336#endif
337END MODULE agrif_lim3_interp
Note: See TracBrowser for help on using the repository browser.