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_ice_interp.F90 in NEMO/trunk/src/NST – NEMO

source: NEMO/trunk/src/NST/agrif_ice_interp.F90 @ 9598

Last change on this file since 9598 was 9598, checked in by nicolasmartin, 6 years ago

Reorganisation plan for NEMO repository: changes to make compilation succeed with new structure
Juste one issue left with AGRIF_NORDIC with AGRIF preprocessing
Standardisation of routines header with version 4.0 and year 2018
Fix for some broken symlinks
File size: 18.3 KB
Line 
1MODULE agrif_ice_interp
2   !!=====================================================================================
3   !!                       ***  MODULE agrif_ice_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_si3 
12   !!----------------------------------------------------------------------
13   !!   'key_si3'  :                                 LIM 3.6 sea-ice model
14   !!   'key_agrif' :                                 AGRIF library
15   !!----------------------------------------------------------------------
16   !!  agrif_interp_si3    : 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   USE phycst , ONLY: rt0
27   
28   IMPLICIT NONE
29   PRIVATE
30
31   PUBLIC   agrif_interp_si3   ! called by agrif_user.F90
32
33   !!----------------------------------------------------------------------
34   !! NEMO/NST 3.6 , NEMO Consortium (2018)
35   !! $Id: agrif_ice_interp.F90 6204 2016-01-04 13:47:06Z cetlod $
36   !! Software governed by the CeCILL licence     (./LICENSE)
37   !!----------------------------------------------------------------------
38
39CONTAINS
40
41   SUBROUTINE agrif_interp_si3( cd_type, kiter, kitermax )
42      !!-----------------------------------------------------------------------
43      !!                 *** ROUTINE agrif_rhg_lim3  ***
44      !!
45      !!  ** Method  : simple call to atomic routines using stored values to
46      !!  fill the boundaries depending of the position of the point and
47      !!  computing factor for time interpolation
48      !!-----------------------------------------------------------------------
49      CHARACTER(len=1), INTENT(in   )           ::   cd_type
50      INTEGER         , INTENT(in   ), OPTIONAL ::   kiter, kitermax
51      !!
52      REAL(wp) ::   zbeta   ! local scalar
53      !!-----------------------------------------------------------------------
54      !
55      IF( Agrif_Root() .OR. nn_ice==0 )  RETURN   ! do not interpolate if inside Parent Grid or if child domain does not have ice
56      !
57      SELECT CASE( cd_type )
58      CASE('U','V')
59         IF( PRESENT( kiter ) ) THEN  ! interpolation at the child ice sub-time step (only for ice rheology)
60            zbeta = ( REAL(lim_nbstep) - REAL(kitermax - kiter) / REAL(kitermax) ) /  &
61               &    ( Agrif_Rhot() * REAL(Agrif_Parent(nn_fsbc)) / REAL(nn_fsbc) )
62         ELSE                         ! interpolation at the child ice time step
63            zbeta = REAL(lim_nbstep) / ( Agrif_Rhot() * REAL(Agrif_Parent(nn_fsbc)) / REAL(nn_fsbc) )
64         ENDIF
65      CASE('T')
66            zbeta = REAL(lim_nbstep) / ( Agrif_Rhot() * REAL(Agrif_Parent(nn_fsbc)) / REAL(nn_fsbc) )
67      END SELECT
68      !
69      Agrif_SpecialValue    = -9999.
70      Agrif_UseSpecialValue = .TRUE.
71      SELECT CASE( cd_type )
72      CASE('U')   ;   CALL Agrif_Bc_variable( u_ice_id  , procname=interp_u_ice  , calledweight=zbeta )
73      CASE('V')   ;   CALL Agrif_Bc_variable( v_ice_id  , procname=interp_v_ice  , calledweight=zbeta )
74      CASE('T')   ;   CALL Agrif_Bc_variable( tra_ice_id, procname=interp_tra_ice, calledweight=zbeta )
75      END SELECT
76      Agrif_SpecialValue    = 0._wp
77      Agrif_UseSpecialValue = .FALSE.
78      !
79   END SUBROUTINE agrif_interp_si3
80
81
82   SUBROUTINE interp_u_ice( ptab, i1, i2, j1, j2, before )
83      !!-----------------------------------------------------------------------
84      !!                     *** ROUTINE interp_u_ice ***
85      !!
86      !! i1 i2 j1 j2 are the index of the boundaries parent(when before) and child (when after)
87      !! To solve issues when parent grid is "land" masked but not all the corresponding child
88      !! grid points, put Agrif_SpecialValue WHERE the parent grid is masked.
89      !! The child solution will be found in the 9(?) points around
90      !!-----------------------------------------------------------------------
91      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
92      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
93      LOGICAL                         , INTENT(in   ) ::   before
94      !!
95      REAL(wp) ::   zrhoy   ! local scalar
96      !!-----------------------------------------------------------------------
97      !
98      IF( before ) THEN  ! parent grid
99         ptab(:,:) = e2u(i1:i2,j1:j2) * u_ice(i1:i2,j1:j2)
100         WHERE( umask(i1:i2,j1:j2,1) == 0. )   ptab(i1:i2,j1:j2) = Agrif_SpecialValue
101      ELSE               ! child grid
102         zrhoy = Agrif_Rhoy()
103         u_ice(i1:i2,j1:j2) = ptab(i1:i2,j1:j2) / ( e2u(i1:i2,j1:j2) * zrhoy ) * umask(i1:i2,j1:j2,1)
104      ENDIF
105      !
106   END SUBROUTINE interp_u_ice
107
108
109   SUBROUTINE interp_v_ice( ptab, i1, i2, j1, j2, before )
110      !!-----------------------------------------------------------------------
111      !!                    *** ROUTINE interp_v_ice ***
112      !!
113      !! i1 i2 j1 j2 are the index of the boundaries parent(when before) and child (when after)
114      !! To solve issues when parent grid is "land" masked but not all the corresponding child
115      !! grid points, put Agrif_SpecialValue WHERE the parent grid is masked.
116      !! The child solution will be found in the 9(?) points around
117      !!-----------------------------------------------------------------------     
118      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
119      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
120      LOGICAL                         , INTENT(in   ) ::   before
121      !!
122      REAL(wp) ::   zrhox   ! local scalar
123      !!-----------------------------------------------------------------------
124      !
125      IF( before ) THEN  ! parent grid
126         ptab(:,:) = e1v(i1:i2,j1:j2) * v_ice(i1:i2,j1:j2)
127         WHERE( vmask(i1:i2,j1:j2,1) == 0. )   ptab(i1:i2,j1:j2) = Agrif_SpecialValue
128      ELSE               ! child grid
129         zrhox = Agrif_Rhox()
130         v_ice(i1:i2,j1:j2) = ptab(i1:i2,j1:j2) / ( e1v(i1:i2,j1:j2) * zrhox ) * vmask(i1:i2,j1:j2,1)
131      ENDIF
132      !
133   END SUBROUTINE interp_v_ice
134
135
136   SUBROUTINE interp_tra_ice( ptab, i1, i2, j1, j2, k1, k2, before, nb, ndir )
137      !!-----------------------------------------------------------------------
138      !!                    *** ROUTINE interp_tra_ice ***                           
139      !!
140      !! i1 i2 j1 j2 are the index of the boundaries parent(when before) and child (when after)
141      !! To solve issues when parent grid is "land" masked but not all the corresponding child
142      !! grid points, put Agrif_SpecialValue WHERE the parent grid is masked.
143      !! The child solution will be found in the 9(?) points around
144      !!-----------------------------------------------------------------------
145      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) ::   ptab
146      INTEGER                               , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2
147      LOGICAL                               , INTENT(in   ) ::   before
148      INTEGER                               , INTENT(in   ) ::   nb, ndir
149      !!
150      INTEGER  ::   ji, jj, jk, jl, jm
151      INTEGER  ::   imin, imax, jmin, jmax
152      LOGICAL  ::   western_side, eastern_side, northern_side, southern_side
153      REAL(wp) ::   zrhox, z1, z2, z3, z4, z5, z6, z7
154      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztab
155      !!-----------------------------------------------------------------------
156      ! tracers are not multiplied by grid cell here => before: * e1e2t ; after: * r1_e1e2t / rhox / rhoy
157      ! and it is ok since we conserve tracers (same as in the ocean).
158      ALLOCATE( ztab(SIZE(a_i,1),SIZE(a_i,2),SIZE(ptab,3)) )
159     
160      IF( before ) THEN  ! parent grid
161         jm = 1
162         DO jl = 1, jpl
163            ptab(i1:i2,j1:j2,jm  ) = a_i (i1:i2,j1:j2,jl)
164            ptab(i1:i2,j1:j2,jm+1) = v_i (i1:i2,j1:j2,jl)
165            ptab(i1:i2,j1:j2,jm+2) = v_s (i1:i2,j1:j2,jl)
166            ptab(i1:i2,j1:j2,jm+3) = sv_i(i1:i2,j1:j2,jl)
167            ptab(i1:i2,j1:j2,jm+4) = oa_i(i1:i2,j1:j2,jl)
168            ptab(i1:i2,j1:j2,jm+5) = a_ip(i1:i2,j1:j2,jl)
169            ptab(i1:i2,j1:j2,jm+6) = v_ip(i1:i2,j1:j2,jl)
170            ptab(i1:i2,j1:j2,jm+7) = t_su(i1:i2,j1:j2,jl)
171            jm = jm + 8
172            DO jk = 1, nlay_s
173               ptab(i1:i2,j1:j2,jm) = e_s(i1:i2,j1:j2,jk,jl)   ;   jm = jm + 1
174            END DO
175            DO jk = 1, nlay_i
176               ptab(i1:i2,j1:j2,jm) = e_i(i1:i2,j1:j2,jk,jl)   ;   jm = jm + 1
177            END DO
178         END DO
179         
180         DO jk = k1, k2
181            WHERE( tmask(i1:i2,j1:j2,1) == 0._wp )   ptab(i1:i2,j1:j2,jk) = Agrif_SpecialValue
182         END DO
183         !
184      ELSE               ! child grid
185         !
186!         IF( nbghostcells > 1 ) THEN   ! ==> The easiest interpolation is used
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                     a_ip(ji,jj,jl) = ptab(ji,jj,jm+5) * tmask(ji,jj,1)
199                     v_ip(ji,jj,jl) = ptab(ji,jj,jm+6) * tmask(ji,jj,1)
200                     t_su(ji,jj,jl) = ptab(ji,jj,jm+7) * tmask(ji,jj,1)
201                  END DO
202               END DO
203               jm = jm + 8
204               !
205               DO jk = 1, nlay_s
206                  e_s(i1:i2,j1:j2,jk,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1)
207                  jm = jm + 1
208               END DO
209               !
210               DO jk = 1, nlay_i
211                  e_i(i1:i2,j1:j2,jk,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1)
212                  jm = jm + 1
213               END DO
214               !
215            END DO
216            !
217!!==> clem: this interpolation does not work because it creates negative values, due
218!!          to negative coefficients when mixing points (for ex. z7)
219!!
220!         ELSE                          ! ==> complex interpolation (only one ghost cell available)
221!            !! Use a more complex interpolation since we mix solutions over a couple of grid points
222!            !! it is advised to use it for fields modified by high order schemes (e.g. advection UM5...)
223!            ! record ztab
224!            jm = 1
225!            DO jl = 1, jpl
226!               ztab(:,:,jm  ) = a_i (:,:,jl)
227!               ztab(:,:,jm+1) = v_i (:,:,jl)
228!               ztab(:,:,jm+2) = v_s (:,:,jl)
229!               ztab(:,:,jm+3) = sv_i(:,:,jl)
230!               ztab(:,:,jm+4) = oa_i(:,:,jl)
231!               ztab(:,:,jm+5) = a_ip(:,:,jl)
232!               ztab(:,:,jm+6) = v_ip(:,:,jl)
233!               ztab(:,:,jm+7) = t_su(:,:,jl)
234!               jm = jm + 8
235!               DO jk = 1, nlay_s
236!                  ztab(:,:,jm) = e_s(:,:,jk,jl)
237!                  jm = jm + 1
238!               END DO
239!               DO jk = 1, nlay_i
240!                  ztab(:,:,jm) = e_i(:,:,jk,jl)
241!                  jm = jm + 1
242!               END DO
243!               !
244!            END DO
245!            !
246!            ! borders of the domain
247!            western_side  = (nb == 1).AND.(ndir == 1)  ;  eastern_side  = (nb == 1).AND.(ndir == 2)
248!            southern_side = (nb == 2).AND.(ndir == 1)  ;  northern_side = (nb == 2).AND.(ndir == 2)
249!            !
250!            ! spatial smoothing
251!            zrhox = Agrif_Rhox()
252!            z1 =      ( zrhox - 1. ) * 0.5
253!            z3 =      ( zrhox - 1. ) / ( zrhox + 1. )
254!            z6 = 2. * ( zrhox - 1. ) / ( zrhox + 1. )
255!            z7 =    - ( zrhox - 1. ) / ( zrhox + 3. )
256!            z2 = 1. - z1
257!            z4 = 1. - z3
258!            z5 = 1. - z6 - z7
259!            !
260!            ! Remove corners
261!            imin = i1  ;  imax = i2  ;  jmin = j1  ;  jmax = j2
262!            IF( (nbondj == -1) .OR. (nbondj == 2) )   jmin = 3
263!            IF( (nbondj == +1) .OR. (nbondj == 2) )   jmax = nlcj-2
264!            IF( (nbondi == -1) .OR. (nbondi == 2) )   imin = 3
265!            IF( (nbondi == +1) .OR. (nbondi == 2) )   imax = nlci-2
266!
267!            ! smoothed fields
268!            IF( eastern_side ) THEN
269!               ztab(nlci,j1:j2,:) = z1 * ptab(nlci,j1:j2,:) + z2 * ptab(nlci-1,j1:j2,:)
270!               DO jj = jmin, jmax
271!                  rswitch = 0.
272!                  IF( u_ice(nlci-2,jj) > 0._wp ) rswitch = 1.
273!                  ztab(nlci-1,jj,:) = ( 1. - umask(nlci-2,jj,1) ) * ztab(nlci,jj,:)  &
274!                     &                +      umask(nlci-2,jj,1)   *  &
275!                     &                ( ( 1. - rswitch ) * ( z4 * ztab(nlci,jj,:)   + z3 * ztab(nlci-2,jj,:) )  &
276!                     &                  +      rswitch   * ( z6 * ztab(nlci-2,jj,:) + z5 * ztab(nlci,jj,:) + z7 * ztab(nlci-3,jj,:) ) )
277!                  ztab(nlci-1,jj,:) = ztab(nlci-1,jj,:) * tmask(nlci-1,jj,1)
278!               END DO
279!            ENDIF
280!            !
281!            IF( northern_side ) THEN
282!               ztab(i1:i2,nlcj,:) = z1 * ptab(i1:i2,nlcj,:) + z2 * ptab(i1:i2,nlcj-1,:)
283!               DO ji = imin, imax
284!                  rswitch = 0.
285!                  IF( v_ice(ji,nlcj-2) > 0._wp ) rswitch = 1.
286!                  ztab(ji,nlcj-1,:) = ( 1. - vmask(ji,nlcj-2,1) ) * ztab(ji,nlcj,:)  &
287!                     &                +      vmask(ji,nlcj-2,1)   *  &
288!                     &                ( ( 1. - rswitch ) * ( z4 * ztab(ji,nlcj,:)   + z3 * ztab(ji,nlcj-2,:) ) &
289!                     &                  +      rswitch   * ( z6 * ztab(ji,nlcj-2,:) + z5 * ztab(ji,nlcj,:) + z7 * ztab(ji,nlcj-3,:) ) )
290!                  ztab(ji,nlcj-1,:) = ztab(ji,nlcj-1,:) * tmask(ji,nlcj-1,1)
291!               END DO
292!            END IF
293!            !
294!            IF( western_side) THEN
295!               ztab(1,j1:j2,:) = z1 * ptab(1,j1:j2,:) + z2 * ptab(2,j1:j2,:)
296!               DO jj = jmin, jmax
297!                  rswitch = 0.
298!                  IF( u_ice(2,jj) < 0._wp ) rswitch = 1.
299!                  ztab(2,jj,:) = ( 1. - umask(2,jj,1) ) * ztab(1,jj,:)  &
300!                     &           +      umask(2,jj,1)   *   &
301!                     &           ( ( 1. - rswitch ) * ( z4 * ztab(1,jj,:) + z3 * ztab(3,jj,:) ) &
302!                     &             +      rswitch   * ( z6 * ztab(3,jj,:) + z5 * ztab(1,jj,:) + z7 * ztab(4,jj,:) ) )
303!                  ztab(2,jj,:) = ztab(2,jj,:) * tmask(2,jj,1)
304!               END DO
305!            ENDIF
306!            !
307!            IF( southern_side ) THEN
308!               ztab(i1:i2,1,:) = z1 * ptab(i1:i2,1,:) + z2 * ptab(i1:i2,2,:)
309!               DO ji = imin, imax
310!                  rswitch = 0.
311!                  IF( v_ice(ji,2) < 0._wp ) rswitch = 1.
312!                  ztab(ji,2,:) = ( 1. - vmask(ji,2,1) ) * ztab(ji,1,:)  &
313!                     &           +      vmask(ji,2,1)   *  &
314!                     &           ( ( 1. - rswitch ) * ( z4 * ztab(ji,1,:) + z3 * ztab(ji,3,:) ) &
315!                     &             +      rswitch   * ( z6 * ztab(ji,3,:) + z5 * ztab(ji,1,:) + z7 * ztab(ji,4,:) ) )
316!                  ztab(ji,2,:) = ztab(ji,2,:) * tmask(ji,2,1)
317!               END DO
318!            END IF
319!            !
320!            ! Treatment of corners
321!            IF( (eastern_side) .AND. ((nbondj == -1).OR.(nbondj == 2)) )  ztab(nlci-1,2,:)      = ptab(nlci-1,2,:)      ! East south
322!            IF( (eastern_side) .AND. ((nbondj ==  1).OR.(nbondj == 2)) )  ztab(nlci-1,nlcj-1,:) = ptab(nlci-1,nlcj-1,:) ! East north
323!            IF( (western_side) .AND. ((nbondj == -1).OR.(nbondj == 2)) )  ztab(2,2,:)           = ptab(2,2,:)           ! West south
324!            IF( (western_side) .AND. ((nbondj ==  1).OR.(nbondj == 2)) )  ztab(2,nlcj-1,:)      = ptab(2,nlcj-1,:)      ! West north
325!           
326!            ! retrieve ice tracers
327!            jm = 1
328!            DO jl = 1, jpl
329!               !
330!               DO jj = j1, j2
331!                  DO ji = i1, i2
332!                     a_i (ji,jj,jl) = ztab(ji,jj,jm  ) * tmask(ji,jj,1)
333!                     v_i (ji,jj,jl) = ztab(ji,jj,jm+1) * tmask(ji,jj,1)
334!                     v_s (ji,jj,jl) = ztab(ji,jj,jm+2) * tmask(ji,jj,1)
335!                     sv_i(ji,jj,jl) = ztab(ji,jj,jm+3) * tmask(ji,jj,1)
336!                     oa_i(ji,jj,jl) = ztab(ji,jj,jm+4) * tmask(ji,jj,1)
337!                     a_ip(ji,jj,jl) = ztab(ji,jj,jm+5) * tmask(ji,jj,1)
338!                     v_ip(ji,jj,jl) = ztab(ji,jj,jm+6) * tmask(ji,jj,1)
339!                     t_su(ji,jj,jl) = ztab(ji,jj,jm+7) * tmask(ji,jj,1)
340!                  END DO
341!               END DO
342!               jm = jm + 8
343!               !
344!               DO jk = 1, nlay_s
345!                  e_s(i1:i2,j1:j2,jk,jl) = ztab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1)
346!                  jm = jm + 1
347!               END DO
348!               !
349!               DO jk = 1, nlay_i
350!                  e_i(i1:i2,j1:j2,jk,jl) = ztab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1)
351!                  jm = jm + 1
352!               END DO
353!               !
354!            END DO
355!         
356!         ENDIF  ! nbghostcells=1
357         
358         DO jl = 1, jpl
359            WHERE( tmask(i1:i2,j1:j2,1) == 0._wp )   t_su(i1:i2,j1:j2,jl) = rt0   ! to avoid a division by 0 in sbcblk.F90
360         END DO
361         !
362      ENDIF
363     
364      DEALLOCATE( ztab )
365      !
366   END SUBROUTINE interp_tra_ice
367
368#else
369   !!----------------------------------------------------------------------
370   !!   Empty module                                             no sea-ice
371   !!----------------------------------------------------------------------
372CONTAINS
373   SUBROUTINE agrif_ice_interp_empty
374      WRITE(*,*)  'agrif_ice_interp : You should not have seen this print! error?'
375   END SUBROUTINE agrif_ice_interp_empty
376#endif
377
378   !!======================================================================
379END MODULE agrif_ice_interp
Note: See TracBrowser for help on using the repository browser.