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

source: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/NST_SRC/agrif_lim3_interp.F90 @ 7953

Last change on this file since 7953 was 7953, checked in by gm, 7 years ago

#1880 (HPC-09): add zdfphy (the ZDF manager) + remove all key_...
File size: 16.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   ! called by ???
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() )  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')   ;   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 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   ! local scalar
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   ! local scalar
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      INTEGER  ::   ji, jj, jk, jl, jm
147      INTEGER  ::   imin, imax, jmin, jmax
148      LOGICAL  ::   western_side, eastern_side, northern_side, southern_side
149      REAL(wp) ::   zrhox, z1, z2, z3, z4, z5, z6, z7
150      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztab
151      !!-----------------------------------------------------------------------
152      ! tracers are not multiplied by grid cell here => before: * e1e2t ; after: * r1_e1e2t / rhox / rhoy
153      ! and it is ok since we conserve tracers (same as in the ocean).
154      ALLOCATE( ztab(SIZE(a_i_b,1),SIZE(a_i_b,2),SIZE(ptab,3)) )
155     
156      IF( before ) THEN  ! parent grid
157         jm = 1
158         DO jl = 1, jpl
159            ptab(i1:i2,j1:j2,jm) = a_i_b  (i1:i2,j1:j2,jl)   ;   jm = jm + 1
160            ptab(i1:i2,j1:j2,jm) = v_i_b  (i1:i2,j1:j2,jl)   ;   jm = jm + 1
161            ptab(i1:i2,j1:j2,jm) = v_s_b  (i1:i2,j1:j2,jl)   ;   jm = jm + 1
162            ptab(i1:i2,j1:j2,jm) = smv_i_b(i1:i2,j1:j2,jl)   ;   jm = jm + 1
163            ptab(i1:i2,j1:j2,jm) = oa_i_b (i1:i2,j1:j2,jl)   ;   jm = jm + 1
164            DO jk = 1, nlay_s
165               ptab(i1:i2,j1:j2,jm) = e_s_b(i1:i2,j1:j2,jk,jl)   ;   jm = jm + 1
166            END DO
167            DO jk = 1, nlay_i
168               ptab(i1:i2,j1:j2,jm) = e_i_b(i1:i2,j1:j2,jk,jl)   ;   jm = jm + 1
169            END DO
170         END DO
171         
172         DO jk = k1, k2
173            WHERE( tmask(i1:i2,j1:j2,1) == 0._wp )   ptab(i1:i2,j1:j2,jk) = -9999.
174         END DO
175         
176      ELSE               ! child grid
177!! ==> The easiest interpolation is the following commented lines
178         jm = 1
179         DO jl = 1, jpl
180            a_i  (i1:i2,j1:j2,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1)   ;   jm = jm + 1
181            v_i  (i1:i2,j1:j2,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1)   ;   jm = jm + 1
182            v_s  (i1:i2,j1:j2,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1)   ;   jm = jm + 1
183            smv_i(i1:i2,j1:j2,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1)   ;   jm = jm + 1
184            oa_i (i1:i2,j1:j2,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1)   ;   jm = jm + 1
185            DO jk = 1, nlay_s
186               e_s(i1:i2,j1:j2,jk,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1)   ;   jm = jm + 1
187            END DO
188            DO jk = 1, nlay_i
189               e_i(i1:i2,j1:j2,jk,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1)   ;   jm = jm + 1
190            END DO
191         END DO
192
193!! ==> this is a more complex interpolation since we mix solutions over a couple of grid points
194!!     it is advised to use it for fields modified by high order schemes (e.g. advection UM5...)
195!!        clem: for some reason (I don't know why), the following lines do not work
196!!              with mpp (or in realistic configurations?). It makes the model crash
197!         ! record ztab
198!         jm = 1
199!         DO jl = 1, jpl
200!            ztab(:,:,jm) = a_i  (:,:,jl) ; jm = jm + 1
201!            ztab(:,:,jm) = v_i  (:,:,jl) ; jm = jm + 1
202!            ztab(:,:,jm) = v_s  (:,:,jl) ; jm = jm + 1
203!            ztab(:,:,jm) = smv_i(:,:,jl) ; jm = jm + 1
204!            ztab(:,:,jm) = oa_i (:,:,jl) ; jm = jm + 1
205!            DO jk = 1, nlay_s
206!               ztab(:,:,jm) = e_s(:,:,jk,jl) ; jm = jm + 1
207!            ENDDO
208!            DO jk = 1, nlay_i
209!               ztab(:,:,jm) = e_i(:,:,jk,jl) ; jm = jm + 1
210!            ENDDO
211!         ENDDO
212!         !
213!         ! borders of the domain
214!         western_side  = (nb == 1).AND.(ndir == 1)  ;  eastern_side  = (nb == 1).AND.(ndir == 2)
215!         southern_side = (nb == 2).AND.(ndir == 1)  ;  northern_side = (nb == 2).AND.(ndir == 2)
216!         !
217!         ! spatial smoothing
218!         zrhox = Agrif_Rhox()
219!         z1 =      ( zrhox - 1. ) * 0.5
220!         z3 =      ( zrhox - 1. ) / ( zrhox + 1. )
221!         z6 = 2. * ( zrhox - 1. ) / ( zrhox + 1. )
222!         z7 =    - ( zrhox - 1. ) / ( zrhox + 3. )
223!         z2 = 1. - z1
224!         z4 = 1. - z3
225!         z5 = 1. - z6 - z7
226!         !
227!         ! Remove corners
228!         imin = i1  ;  imax = i2  ;  jmin = j1  ;  jmax = j2
229!         IF( (nbondj == -1) .OR. (nbondj == 2) )   jmin = 3
230!         IF( (nbondj == +1) .OR. (nbondj == 2) )   jmax = nlcj-2
231!         IF( (nbondi == -1) .OR. (nbondi == 2) )   imin = 3
232!         IF( (nbondi == +1) .OR. (nbondi == 2) )   imax = nlci-2
233!
234!         ! smoothed fields
235!         IF( eastern_side ) THEN
236!            ztab(nlci,j1:j2,:) = z1 * ptab(nlci,j1:j2,:) + z2 * ptab(nlci-1,j1:j2,:)
237!            DO jj = jmin, jmax
238!               rswitch = 0.
239!               IF( u_ice(nlci-2,jj) > 0._wp ) rswitch = 1.
240!               ztab(nlci-1,jj,:) = ( 1. - umask(nlci-2,jj,1) ) * ztab(nlci,jj,:)  &
241!                  &                +      umask(nlci-2,jj,1)   *  &
242!                  &                ( ( 1. - rswitch ) * ( z4 * ztab(nlci,jj,:)   + z3 * ztab(nlci-2,jj,:) )  &
243!                  &                  +      rswitch   * ( z6 * ztab(nlci-2,jj,:) + z5 * ztab(nlci,jj,:) + z7 * ztab(nlci-3,jj,:) ) )
244!               ztab(nlci-1,jj,:) = ztab(nlci-1,jj,:) * tmask(nlci-1,jj,1)
245!            END DO
246!         ENDIF
247!         !
248!         IF( northern_side ) THEN
249!            ztab(i1:i2,nlcj,:) = z1 * ptab(i1:i2,nlcj,:) + z2 * ptab(i1:i2,nlcj-1,:)
250!            DO ji = imin, imax
251!               rswitch = 0.
252!               IF( v_ice(ji,nlcj-2) > 0._wp ) rswitch = 1.
253!               ztab(ji,nlcj-1,:) = ( 1. - vmask(ji,nlcj-2,1) ) * ztab(ji,nlcj,:)  &
254!                  &                +      vmask(ji,nlcj-2,1)   *  &
255!                  &                ( ( 1. - rswitch ) * ( z4 * ztab(ji,nlcj,:)   + z3 * ztab(ji,nlcj-2,:) ) &
256!                  &                  +      rswitch   * ( z6 * ztab(ji,nlcj-2,:) + z5 * ztab(ji,nlcj,:) + z7 * ztab(ji,nlcj-3,:) ) )
257!               ztab(ji,nlcj-1,:) = ztab(ji,nlcj-1,:) * tmask(ji,nlcj-1,1)
258!            END DO
259!         END IF
260!         !
261!         IF( western_side) THEN
262!            ztab(1,j1:j2,:) = z1 * ptab(1,j1:j2,:) + z2 * ptab(2,j1:j2,:)
263!            DO jj = jmin, jmax
264!               rswitch = 0.
265!               IF( u_ice(2,jj) < 0._wp ) rswitch = 1.
266!               ztab(2,jj,:) = ( 1. - umask(2,jj,1) ) * ztab(1,jj,:)  &
267!                  &           +      umask(2,jj,1)   *   &
268!                  &           ( ( 1. - rswitch ) * ( z4 * ztab(1,jj,:) + z3 * ztab(3,jj,:) ) &
269!                  &             +      rswitch   * ( z6 * ztab(3,jj,:) + z5 * ztab(1,jj,:) + z7 * ztab(4,jj,:) ) )
270!               ztab(2,jj,:) = ztab(2,jj,:) * tmask(2,jj,1)
271!            END DO
272!         ENDIF
273!         !
274!         IF( southern_side ) THEN
275!            ztab(i1:i2,1,:) = z1 * ptab(i1:i2,1,:) + z2 * ptab(i1:i2,2,:)
276!            DO ji = imin, imax
277!               rswitch = 0.
278!               IF( v_ice(ji,2) < 0._wp ) rswitch = 1.
279!               ztab(ji,2,:) = ( 1. - vmask(ji,2,1) ) * ztab(ji,1,:)  &
280!                  &           +      vmask(ji,2,1)   *  &
281!                  &           ( ( 1. - rswitch ) * ( z4 * ztab(ji,1,:) + z3 * ztab(ji,3,:) ) &
282!                  &             +      rswitch   * ( z6 * ztab(ji,3,:) + z5 * ztab(ji,1,:) + z7 * ztab(ji,4,:) ) )
283!               ztab(ji,2,:) = ztab(ji,2,:) * tmask(ji,2,1)
284!            END DO
285!         END IF
286!         !
287!         ! Treatment of corners
288!         IF( (eastern_side) .AND. ((nbondj == -1).OR.(nbondj == 2)) )  ztab(nlci-1,2,:)      = ptab(nlci-1,2,:)      ! East south
289!         IF( (eastern_side) .AND. ((nbondj ==  1).OR.(nbondj == 2)) )  ztab(nlci-1,nlcj-1,:) = ptab(nlci-1,nlcj-1,:) ! East north
290!         IF( (western_side) .AND. ((nbondj == -1).OR.(nbondj == 2)) )  ztab(2,2,:)           = ptab(2,2,:)           ! West south
291!         IF( (western_side) .AND. ((nbondj ==  1).OR.(nbondj == 2)) )  ztab(2,nlcj-1,:)      = ptab(2,nlcj-1,:)      ! West north
292!
293!         ! retrieve ice tracers
294!         jm = 1
295!         DO jl = 1, jpl
296!            a_i  (i1:i2,j1:j2,jl) = ztab(i1:i2,j1:j2,jm) ; jm = jm + 1
297!            v_i  (i1:i2,j1:j2,jl) = ztab(i1:i2,j1:j2,jm) ; jm = jm + 1
298!            v_s  (i1:i2,j1:j2,jl) = ztab(i1:i2,j1:j2,jm) ; jm = jm + 1
299!            smv_i(i1:i2,j1:j2,jl) = ztab(i1:i2,j1:j2,jm) ; jm = jm + 1
300!            oa_i (i1:i2,j1:j2,jl) = ztab(i1:i2,j1:j2,jm) ; jm = jm + 1
301!            DO jk = 1, nlay_s
302!               e_s(i1:i2,j1:j2,jk,jl) = ztab(i1:i2,j1:j2,jm) ; jm = jm + 1
303!            ENDDO
304!            DO jk = 1, nlay_i
305!               e_i(i1:i2,j1:j2,jk,jl) = ztab(i1:i2,j1:j2,jm) ; jm = jm + 1
306!            ENDDO
307!         ENDDO
308       
309         ! integrated values
310         vt_i (i1:i2,j1:j2) = SUM( v_i(i1:i2,j1:j2,:), dim=3 )
311         vt_s (i1:i2,j1:j2) = SUM( v_s(i1:i2,j1:j2,:), dim=3 )
312         at_i (i1:i2,j1:j2) = SUM( a_i(i1:i2,j1:j2,:), dim=3 )
313         et_s(i1:i2,j1:j2)  = SUM( SUM( e_s(i1:i2,j1:j2,:,:), dim=4 ), dim=3 )
314         et_i(i1:i2,j1:j2)  = SUM( SUM( e_i(i1:i2,j1:j2,:,:), dim=4 ), dim=3 )
315         !
316      ENDIF
317     
318      DEALLOCATE( ztab )
319      !
320   END SUBROUTINE interp_tra_ice
321
322#else
323   !!----------------------------------------------------------------------
324   !!   Empty module                                             no sea-ice
325   !!----------------------------------------------------------------------
326CONTAINS
327   SUBROUTINE agrif_lim3_interp_empty
328      WRITE(*,*)  'agrif_lim3_interp : You should not have seen this print! error?'
329   END SUBROUTINE agrif_lim3_interp_empty
330#endif
331
332   !!======================================================================
333END MODULE agrif_lim3_interp
Note: See TracBrowser for help on using the repository browser.