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

source: branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/NST_SRC/agrif_lim3_interp.F90 @ 8738

Last change on this file since 8738 was 8738, checked in by dancopsey, 6 years ago

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