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

source: branches/2016/dev_merge_2016/NEMOGCM/NEMO/NST_SRC/agrif_lim3_interp.F90 @ 7421

Last change on this file since 7421 was 7421, checked in by flavoni, 7 years ago

#1811 merge dev_CNRS_MERATOR_2016 with dev_merge_2016 branch
File size: 15.6 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      IF( PRESENT( kiter ) ) THEN  ! interpolation at the child sub-time step (for ice rheology)
57         zbeta = ( REAL(lim_nbstep) - REAL(kitermax - kiter) / REAL(kitermax) ) /  &
58            &    ( Agrif_Rhot() * REAL(Agrif_Parent(nn_fsbc)) / REAL(nn_fsbc) )
59      ELSE                         ! interpolation at the child time step
60         zbeta = REAL(lim_nbstep) / ( Agrif_Rhot() * REAL(Agrif_Parent(nn_fsbc)) / REAL(nn_fsbc) )
61      ENDIF
62      !
63      Agrif_SpecialValue=-9999.
64      Agrif_UseSpecialValue = .TRUE.
65      SELECT CASE(cd_type)
66      CASE('U')
67         CALL Agrif_Bc_variable( u_ice_id  , procname=interp_u_ice  , calledweight=zbeta )
68      CASE('V')
69         CALL Agrif_Bc_variable( v_ice_id  , procname=interp_v_ice  , calledweight=zbeta )
70      CASE('T')
71         CALL Agrif_Bc_variable( tra_ice_id, procname=interp_tra_ice, calledweight=zbeta )
72      END SELECT
73      Agrif_SpecialValue=0.
74      Agrif_UseSpecialValue = .FALSE.
75      !
76   END SUBROUTINE agrif_interp_lim3
77
78   !!------------------
79   !! Local subroutines
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 grid points,
87      !! put -9999 WHERE the parent grid is masked. The child solution will be found in the 9(?) points around
88      !!-----------------------------------------------------------------------
89      INTEGER , INTENT(in) :: i1, i2, j1, j2
90      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab
91      LOGICAL , INTENT(in) :: before
92      !!
93      REAL(wp) :: zrhoy
94      !!-----------------------------------------------------------------------
95      !
96      IF( before ) THEN  ! parent grid
97         ptab(:,:) = e2u(i1:i2,j1:j2) * u_ice_b(i1:i2,j1:j2)
98         WHERE( umask(i1:i2,j1:j2,1) == 0. )  ptab(:,:) = -9999.
99      ELSE               ! child grid
100         zrhoy = Agrif_Rhoy()
101         u_ice(i1:i2,j1:j2) = ptab(:,:) / ( e2u(i1:i2,j1:j2) * zrhoy ) * umask(i1:i2,j1:j2,1)
102      ENDIF
103      !
104   END SUBROUTINE interp_u_ice
105
106
107   SUBROUTINE interp_v_ice( ptab, i1, i2, j1, j2, before )
108      !!-----------------------------------------------------------------------
109      !!                    *** ROUTINE interp_v_ice ***
110      !!
111      !! i1 i2 j1 j2 are the index of the boundaries parent(when before) and child (when after)
112      !! To solve issues when parent grid is "land" masked but not all the corresponding child grid points,
113      !! put -9999 WHERE the parent grid is masked. The child solution will be found in the 9(?) points around
114      !!-----------------------------------------------------------------------     
115      INTEGER , INTENT(in) :: i1, i2, j1, j2
116      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab
117      LOGICAL , INTENT(in) :: before
118      !!
119      REAL(wp) :: zrhox
120      !!-----------------------------------------------------------------------
121      !
122      IF( before ) THEN  ! parent grid
123         ptab(:,:) = e1v(i1:i2,j1:j2) * v_ice_b(i1:i2,j1:j2)
124         WHERE( vmask(i1:i2,j1:j2,1) == 0. )  ptab(:,:) = -9999.
125      ELSE               ! child grid
126         zrhox = Agrif_Rhox()
127         v_ice(i1:i2,j1:j2) = ptab(:,:) / ( e1v(i1:i2,j1:j2) * zrhox ) * vmask(i1:i2,j1:j2,1)
128      ENDIF
129      !
130   END SUBROUTINE interp_v_ice
131
132
133   SUBROUTINE interp_tra_ice( ptab, i1, i2, j1, j2, k1, k2, before, nb, ndir )
134      !!-----------------------------------------------------------------------
135      !!                    *** ROUTINE interp_tra_ice ***                           
136      !!
137      !! i1 i2 j1 j2 are the index of the boundaries parent(when before) and child (when after)
138      !! To solve issues when parent grid is "land" masked but not all the corresponding child grid points,
139      !! put -9999 WHERE the parent grid is masked. The child solution will be found in the 9(?) points around
140      !!-----------------------------------------------------------------------
141      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
142      INTEGER , INTENT(in) :: i1, i2, j1, j2, k1, k2
143      LOGICAL , INTENT(in) :: before
144      INTEGER , INTENT(in) :: nb, ndir
145      !!
146      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztab
147      INTEGER  ::   ji, jj, jk, jl, jm
148      INTEGER  ::   imin, imax, jmin, jmax
149      REAL(wp) ::   zrhox, z1, z2, z3, z4, z5, z6, z7
150      LOGICAL  ::   western_side, eastern_side, northern_side, southern_side
151
152      !!-----------------------------------------------------------------------
153      ! clem: pkoi on n'utilise pas les quantités intégrées ici => before: * e1e2t ; after: * r1_e1e2t / rhox / rhoy
154      ! a priori c'est ok comme ca (cf ce qui est fait dans l'ocean). Je ne sais pas pkoi ceci dit
155      ALLOCATE( ztab(SIZE(a_i_b,1),SIZE(a_i_b,2),SIZE(ptab,3)) )
156     
157      IF( before ) THEN  ! parent grid
158         jm = 1
159         DO jl = 1, jpl
160            ptab(i1:i2,j1:j2,jm) = a_i_b  (i1:i2,j1:j2,jl) ; jm = jm + 1
161            ptab(i1:i2,j1:j2,jm) = v_i_b  (i1:i2,j1:j2,jl) ; jm = jm + 1
162            ptab(i1:i2,j1:j2,jm) = v_s_b  (i1:i2,j1:j2,jl) ; jm = jm + 1
163            ptab(i1:i2,j1:j2,jm) = smv_i_b(i1:i2,j1:j2,jl) ; jm = jm + 1
164            ptab(i1:i2,j1:j2,jm) = oa_i_b (i1:i2,j1:j2,jl) ; jm = jm + 1
165            DO jk = 1, nlay_s
166               ptab(i1:i2,j1:j2,jm) = e_s_b(i1:i2,j1:j2,jk,jl) ; jm = jm + 1
167            ENDDO
168            DO jk = 1, nlay_i
169               ptab(i1:i2,j1:j2,jm) = e_i_b(i1:i2,j1:j2,jk,jl) ; jm = jm + 1
170            ENDDO
171         ENDDO
172         
173         DO jk = k1, k2
174            WHERE( tmask(i1:i2,j1:j2,1) == 0. )  ptab(i1:i2,j1:j2,jk) = -9999.
175         ENDDO
176         
177      ELSE               ! child grid
178!! ==> The easiest interpolation is the following commented lines
179!!         jm = 1
180!!         DO jl = 1, jpl
181!!            a_i  (i1:i2,j1:j2,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
182!!            v_i  (i1:i2,j1:j2,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
183!!            v_s  (i1:i2,j1:j2,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
184!!            smv_i(i1:i2,j1:j2,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
185!!            oa_i (i1:i2,j1:j2,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
186!!            DO jk = 1, nlay_s
187!!               e_s(i1:i2,j1:j2,jk,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
188!!            ENDDO
189!!            DO jk = 1, nlay_i
190!!               e_i(i1:i2,j1:j2,jk,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
191!!            ENDDO
192!!         ENDDO
193
194!! ==> this is a more complex interpolation since we mix solutions over a couple of grid points
195!!     it is advised to use it for fields modified by high order schemes (e.g. advection UM5...)
196         ! record ztab
197         jm = 1
198         DO jl = 1, jpl
199            ztab(:,:,jm) = a_i_b  (:,:,jl) ; jm = jm + 1
200            ztab(:,:,jm) = v_i_b  (:,:,jl) ; jm = jm + 1
201            ztab(:,:,jm) = v_s_b  (:,:,jl) ; jm = jm + 1
202            ztab(:,:,jm) = smv_i_b(:,:,jl) ; jm = jm + 1
203            ztab(:,:,jm) = oa_i_b (:,:,jl) ; jm = jm + 1
204            DO jk = 1, nlay_s
205               ztab(:,:,jm) = e_s_b(:,:,jk,jl) ; jm = jm + 1
206            ENDDO
207            DO jk = 1, nlay_i
208               ztab(:,:,jm) = e_i_b(:,:,jk,jl) ; jm = jm + 1
209            ENDDO
210         ENDDO
211         !
212         ! borders of the domain
213         western_side  = (nb == 1).AND.(ndir == 1)  ;  eastern_side  = (nb == 1).AND.(ndir == 2)
214         southern_side = (nb == 2).AND.(ndir == 1)  ;  northern_side = (nb == 2).AND.(ndir == 2)
215         !
216         ! spatial smoothing
217         zrhox = Agrif_Rhox()
218         z1 =      ( zrhox - 1. ) * 0.5 
219         z3 =      ( zrhox - 1. ) / ( zrhox + 1. )
220         z6 = 2. * ( zrhox - 1. ) / ( zrhox + 1. )
221         z7 =    - ( zrhox - 1. ) / ( zrhox + 3. )
222         z2 = 1. - z1
223         z4 = 1. - z3
224         z5 = 1. - z6 - z7
225         !
226         ! Remove corners
227         imin = i1  ;  imax = i2  ;  jmin = j1  ;  jmax = j2
228         IF( (nbondj == -1) .OR. (nbondj == 2) )   jmin = 3
229         IF( (nbondj == +1) .OR. (nbondj == 2) )   jmax = nlcj-2
230         IF( (nbondi == -1) .OR. (nbondi == 2) )   imin = 3
231         IF( (nbondi == +1) .OR. (nbondi == 2) )   imax = nlci-2
232
233         ! smoothed fields
234         IF( eastern_side ) THEN
235            ztab(nlci,j1:j2,:) = z1 * ptab(nlci,j1:j2,:) + z2 * ptab(nlci-1,j1:j2,:)
236            DO jj = jmin, jmax
237               rswitch = 0.
238               IF( u_ice(nlci-2,jj) > 0._wp ) rswitch = 1.
239               ztab(nlci-1,jj,:) = ( 1. - umask(nlci-2,jj,1) ) * ztab(nlci,jj,:)  &
240                  &                +      umask(nlci-2,jj,1)   *  &
241                  &                ( ( 1. - rswitch ) * ( z4 * ztab(nlci,jj,:)   + z3 * ztab(nlci-2,jj,:) )  &
242                  &                  +      rswitch   * ( z6 * ztab(nlci-2,jj,:) + z5 * ztab(nlci,jj,:) + z7 * ztab(nlci-3,jj,:) ) )
243               ztab(nlci-1,jj,:) = ztab(nlci-1,jj,:) * tmask(nlci-1,jj,1)
244            END DO
245         ENDIF
246         !
247         IF( northern_side ) THEN
248            ztab(i1:i2,nlcj,:) = z1 * ptab(i1:i2,nlcj,:) + z2 * ptab(i1:i2,nlcj-1,:)
249            DO ji = imin, imax
250               rswitch = 0.
251               IF( v_ice(ji,nlcj-2) > 0._wp ) rswitch = 1.
252               ztab(ji,nlcj-1,:) = ( 1. - vmask(ji,nlcj-2,1) ) * ztab(ji,nlcj,:)  &
253                  &                +      vmask(ji,nlcj-2,1)   *  &
254                  &                ( ( 1. - rswitch ) * ( z4 * ztab(ji,nlcj,:)   + z3 * ztab(ji,nlcj-2,:) ) &
255                  &                  +      rswitch   * ( z6 * ztab(ji,nlcj-2,:) + z5 * ztab(ji,nlcj,:) + z7 * ztab(ji,nlcj-3,:) ) )
256               ztab(ji,nlcj-1,:) = ztab(ji,nlcj-1,:) * tmask(ji,nlcj-1,1)
257            END DO
258         END IF
259         !
260         IF( western_side) THEN
261            ztab(1,j1:j2,:) = z1 * ptab(1,j1:j2,:) + z2 * ptab(2,j1:j2,:)
262            DO jj = jmin, jmax
263               rswitch = 0.
264               IF( u_ice(2,jj) > 0._wp ) rswitch = 1.
265               ztab(2,jj,:) = ( 1. - umask(2,jj,1) ) * ztab(1,jj,:)  &
266                  &           +      umask(2,jj,1)   *   &
267                  &           ( ( 1. - rswitch ) * ( z4 * ztab(1,jj,:) + z3 * ztab(3,jj,:) ) &
268                  &             +      rswitch   * ( z6 * ztab(3,jj,:) + z5 * ztab(1,jj,:) + z7 * ztab(4,jj,:) ) )
269               ztab(2,jj,:) = ztab(2,jj,:) * tmask(2,jj,1)
270            END DO
271         ENDIF
272         !
273         IF( southern_side ) THEN
274            ztab(i1:i2,1,:) = z1 * ptab(i1:i2,1,:) + z2 * ptab(i1:i2,2,:)
275            DO ji = imin, imax
276               rswitch = 0.
277               IF( v_ice(ji,2) > 0._wp ) rswitch = 1.
278               ztab(ji,2,:) = ( 1. - vmask(ji,2,1) ) * ztab(ji,1,:)  &
279                  &           +      vmask(ji,2,1)   *  &
280                  &           ( ( 1. - rswitch ) * ( z4 * ztab(ji,1,:) + z3 * ztab(ji,3,:) ) &
281                  &             +      rswitch   * ( z6 * ztab(ji,3,:) + z5 * ztab(ji,1,:) + z7 * ztab(ji,4,:) ) )
282               ztab(ji,2,:) = ztab(ji,2,:) * tmask(ji,2,1)
283            END DO
284         END IF
285         !
286         ! Treatment of corners
287         IF( (eastern_side) .AND. ((nbondj == -1).OR.(nbondj == 2)) )  ztab(nlci-1,2,:)      = ptab(nlci-1,2,:)      ! East south
288         IF( (eastern_side) .AND. ((nbondj ==  1).OR.(nbondj == 2)) )  ztab(nlci-1,nlcj-1,:) = ptab(nlci-1,nlcj-1,:) ! East north
289         IF( (western_side) .AND. ((nbondj == -1).OR.(nbondj == 2)) )  ztab(2,2,:)           = ptab(2,2,:)           ! West south
290         IF( (western_side) .AND. ((nbondj ==  1).OR.(nbondj == 2)) )  ztab(2,nlcj-1,:)      = ptab(2,nlcj-1,:)      ! West north
291
292         ! retrieve ice tracers
293         jm = 1
294         DO jl = 1, jpl
295            a_i  (i1:i2,j1:j2,jl) = ztab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
296            v_i  (i1:i2,j1:j2,jl) = ztab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
297            v_s  (i1:i2,j1:j2,jl) = ztab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
298            smv_i(i1:i2,j1:j2,jl) = ztab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
299            oa_i (i1:i2,j1:j2,jl) = ztab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
300            DO jk = 1, nlay_s
301               e_s(i1:i2,j1:j2,jk,jl) = ztab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
302            ENDDO
303            DO jk = 1, nlay_i
304               e_i(i1:i2,j1:j2,jk,jl) = ztab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
305            ENDDO
306         ENDDO
307         
308      ENDIF
309     
310      DEALLOCATE( ztab )
311      !
312   END SUBROUTINE interp_tra_ice
313
314#else
315CONTAINS
316   SUBROUTINE agrif_lim3_interp_empty
317      !!---------------------------------------------
318      !!   *** ROUTINE agrif_lim3_interp_empty ***
319      !!---------------------------------------------
320      WRITE(*,*)  'agrif_lim3_interp : You should not have seen this print! error?'
321   END SUBROUTINE agrif_lim3_interp_empty
322#endif
323END MODULE agrif_lim3_interp
Note: See TracBrowser for help on using the repository browser.